The high controllability of analogue quantum simulators using ultracold atoms in optical lattices pushes forward the frontiers in the experimental investigation of the fermionic Hubbard model. The bilayer Hubbard model is a step beyond the two-dimensional Hubbard model that extends the latter by incorporating a coupling between two two-dimensional Hubbard systems. This is also a step forward in the idea of analogue quantum simulate real materials such as the copper oxide high temperature superconductors which possess a coupled layer structure. This thesis is dedicated to the experimental implementation of an analogue quantum simulator for a bilayer Hubbard system with cold atoms in optical mono-and bi-chromatic lattices. The measurement of competing magnetic order in the bilayer Hubbard system is a result of the work during the course of this thesis. This measurement requires a high controllability of the system which goes along with a precise calibration and fundamental characterisation of the implemented bilayer Hubbard system. Here, the calibration of the interaction strength by means of a comparison between data and theoretical predictions becomes possible through one major outcome of this thesis. This is a method to compute interacting Wannier functions in an optical superlattice. A further major outcome of this thesis is the measurement of thermodynamics, density fluctuations and entropy in the bilayer system considered as a monolayer Hubbard system with reservoir. Further outcomes of this thesis for a fundamental characterisation of the bilayer Hubbard system are:1. The characterisation of the Hubbard band insulator.2. The computation of the potential map for an optical bi-chromatic superlattice. The calibration of the optical bi-chromatic superlattice by comparing experimentaldata to theoretical predictions. Here, the band projection position operator method to compute non-interacting Wannier functions in a superlattice was successfully implemented in this thesis. The latter are the starting point for the newly developed method to compute interacting Wannier functions.V Furthermore, I thank Prof. Dr. Reinhard Kienberger and Prof. Dr. Florian Schreck for their support during the Master's thesis. VII Bibliography 133List of Figures 1451 This is in contrast to the scattering force which light exerts on atoms when its frequency is close to the resonance of the atomic transition, i.e. ∆ ≈ 0. This scattering force saturates for strong intensity and decreases with 1/∆ 2 while the dipole force does not saturate and only decreases with 1/∆. 2 Here, ΩR is the Rabi frequency in the two level picture.