Quantum Simulation/Analog quantum simulation

Resource type: this resource contains a lecture or lecture notes.

Ultracold atoms

We have seen in the previous chapters that classical simulation methods do not work for many interesting many-body problems, most notably involving fermions or frustrated interactions in two or three spatial dimensions. The root of this difficulty arises from the Hilbert space dimension growing exponentially with the size of the system. Therefore, it looks like we need to use a simulator in which the available resources scale in a similar fashion, meaning the simulator has to be a quantum system itself! This was first recognized by Richard Feynman in 1982 and is widely thought to mark the birth of quantum simulation. A good quantum simulator is a physical system which can be well controlled so it can mimic basically any other quantum system that might be of interest. Here, ultracold atoms have been proven to be very attractive:

1. Individual atoms are well isolated from the environment,
2. Their properties can be widely tuned with electric fields, magnetic fields, microwaves, and lasers.
3. Theoretical predictions can be derived from first principles.
4. Ultracold temperatures in the nanokelvin range a accessible using laser cooling methods.

Ultracold atoms can be either bosonic or fermionic, which depends on the total spin of all constituents. As neutral atoms contain as many protons as electrons, their contribution to the spin will always be an integer number. Consequently, the statistical properties are set by the number of neutrons, i.e., we have bosons for an even number of neutrons and fermions for an odd number of neutrons. A prominent example of a bosonic atom is $^{87}{\textrm {Rb}}$ (50 neutrons), while an important fermionic isotope is $^{40}{\textrm {K}}$ (21 neutrons). Note that in the case of potassium, we can actually have both a fermionic and a bosonic ($^{39}{\textrm {K}}$) isotope.

A many-body system of ultracold atoms is most coveniently described in the formalism of the second quantization. Let us first consider spinless bosons, for which we choose a Fock space for each possible mode as a suitable basis, i.e., eigenstates of the particle number operator for each mode, $\{|n\rangle _{k}\}$. The total Hilbert space is then formed as the tensor product of all these single-mode Fock spaces, ${\mathcal {H}}=\otimes _{k}{\mathcal {H}}_{k}$. The creation and annihilation operators for the Fock states satisfy the commutation relations known from the harmonic oscillator:

$[a_{k},a_{l}]=[a_{k}^{\dagger },a_{l}^{\dagger }]=0$

$[a_{k},a_{l}^{\dagger }]=\delta _{kl}.$

Additionally, the annihilation operator vanishes when applied to the vacuum state of the mode, $a_{k}|0\rangle _{k}=0$, and the number operator for each mode is simply $n=a_{k}^{\dagger }a_{k}$.

For fermions, the commutation relations are replaced by anticommutation relations:

$\{a_{k},a_{l}\}=\{a_{k}^{\dagger },a_{l}^{\dagger }\}=0$

$\{a_{k},a_{l}^{\dagger }\}=\delta _{kl}.$

A consequence of these relations is the Pauli exclusion principle, which states that two fermions cannot occupy the same state, i.e., $a_{k}^{\dagger }a_{k}^{\dagger }|0\rangle _{k}=0$. As operators between different modes no longer commute (they anti-commute instead), we have to define an ordering of the modes. Usually, a multi-mode Fock state is defined as

$|n_{1},n_{2},\ldots n_{i},\ldots \rangle =(a_{i}^{\dagger })_{i}^{n}(a_{2}^{\dagger })_{2}^{n}(a_{1}^{\dagger })_{1}^{n}|0,0,\ldots ,0,\ldots \rangle ,$

i.e., the first mode appears on the right. To understand the subtle point about the importance of ordering, consider this two-mode example:

$a_{1}^{\dagger }|0,1\rangle =a_{1}^{\dagger }a_{2}^{\dagger }|0,0\rangle =-a_{2}^{\dagger }a_{1}^{\dagger }|0,0\rangle =-|1,1\rangle .$

Accounting for spin is straightforward in this formalism and only requires to replace the operators $a_{k}$ (and its hermitian conjugate) by $a_{k,\sigma }$, where $\sigma$ is the spin variable. We can also define field operators, which will allow us to express observables and the Hamiltonian in a convenient way:

$\psi (r)={\frac {1}{\sqrt {V}}}\sum _{k}e^{ikr}a_{k}$

$\psi ^{\dagger }(r)={\frac {1}{\sqrt {V}}}\sum _{k}e^{-ikr}a_{k}^{\dagger },$

where $V$ is the quantization volume. Depending on the statistics of the creation and annihilation operators, the field operators will obey bosonic or fermionic statistics as well. Since the field creation operator is the Fourier transform of the mode creation operator $a_{k}^{\dagger }$, $\psi ^{\dagger }(r)$ acting on the vacuum will simply create a particle at position $r$. We can use this to define the density operator

$n(r)=\psi ^{\dagger }(r)\psi (r),$

which by integration over the full space gives the total particle number, as expected:

$N=\int \mathrm {d} r\psi ^{\dagger }(r)\psi (r)=\int \mathrm {d} r{\frac {1}{V}}\sum _{k,k'}e^{i(k-k')r}a_{k}^{\dagger }a_{k}=\sum _{k}a_{k}^{\dagger }a_{k}.$

We are now in the position to formulate the Hamiltonian for the many-body system of ultracold atoms. The kinetic energy is diagonal when using creation and annihilation operators,

$H_{\mathrm {kin} }=\sum _{k}{\frac {\hbar ^{2}k^{2}}{2m}}a_{k}^{\dagger }a_{k}=\int \mathrm {d} r{\frac {\hbar ^{2}}{2m}}[\nabla \psi ^{\dagger }(r)][\nabla \psi (r)].$

Likewise, the energy associated to an external potential (e.g., describing a trapping potential) is given by

$H_{\mathrm {pot} }=\int \mathrm {d} rn(r)V(r)=\int \mathrm {d} r\psi ^{\dagger }(r)\psi (r)V(r)$.

Finally, the only thing that is missing is the interaction term. In most cases the interaction will only depend on the spatial separation of two particles and hence can be expressed as

$H_{\mathrm {int} }={\frac {1}{2}}\int \mathrm {d} r\mathrm {d} r'\psi ^{\dagger }(r)\psi ^{\dagger }(r')V_{\mathrm {int} }(r-r')\psi (r')\psi (r)={\frac {1}{V}}\sum _{k.k',q}{\tilde {V}}_{q}a_{k+q}^{\dagger }a_{k'-q}^{\dagger }a_{k'}a_{k},$

where ${\tilde {V}}_{q}$ is the Fourier transform of the interaction potential. Note that for fermions, the ordering of the creation and annihilation operators is important, while for bosons it only matters that creation operators appear to the left of annihilation operators.

To develop an understanding of the interactions arising between ultracold atoms, we first have to look into their electronic structure in some detail. We focus on alkali atoms like Rb as they are the most commonly used ones in ultracold atoms experiments. In their electronic ground states, two Rb atoms have their single valence electron in the $|5s\rangle$, i.e., $|5s,5s\rangle$ in the two-atom basis. The atoms are neutral and their magnetic dipole moments of $1\,\mu _{B}$ are negligible, so the dominant contribution comes from a van der Waals interaction behaving as $-C_{6}/r^{6}$, which is a second-order off-resonant dipole-dipole interaction. Here, the dominant term comes from the dipole transitions to the first electronic excited state, resulting in

$C_{6}\sim {\frac {\langle 5_{s},5_{s}|x_{1}x_{2}|5_{p},5p\rangle \langle 5p,5p|x_{1}x_{2}|5s,5s\rangle }{2\Delta }},$

where $x_{1}$ and $x_{2}$ are the coordinates for the valence electrons and $\Delta$ is the energyh gap to the first excited state. Using the known values for Rb , we obtain a value for the van der Waals coefficient of $C_{6}=3100\mathrm {a.u.}$, which is not too far off from the experimentally measure value of $C_{6}=4707\mathrm {a.u.}$ . This large value of the van der Waals coefficient can be understood as the single valence electron in alkali atoms behaving almost hydrogenic (in hydrogen, the excitation gap $\Delta$ vanishes for $ns-np$ transitions).

While the van der Waals interaction gives a correct description at large separation, at short ranges core repulsion will kick in, leading to a molecular potential supporting many bound states. However, in many cases the short range details are not important. We consider the two-atom Schrödinger equation,

$\left[-{\frac {1}{2\mu }}\Delta +V({\textbf {r}})\right]\psi ({\textbf {r}})=E\psi ({\textbf {r}}),$

with $\mu$ being the reduced mass. At long distances, we can look at its asymptotic behavior ,

$\psi ({\textbf {r}})=e^{ikx}+f(\theta ,\phi ){\frac {e^{ikr}}{r}}.$

This first term in this expression describes an incoming plane wave, while the second term describes an outgoing spherical wave after the scattering event between the two atoms. For radially symmetric interaction potentials, the dependence on the azimuthal angle drops out and the scattering amplitude $f(\theta )$ can be expressed with the help of Legendre polynomials as

$f(\theta )=\sum _{l=0}^{\infty }f_{l}P_{l}(\cos \theta ).$

In this form, we have an expression in terms of different partial waves with angular momentum $l$. In analogy to atomic angular momentum, the $l=0$ channel is called $s$--wave scattering, the $l=1$ channel is called $p$-wave scattering and so on. The partial wave amplitudes $f_{l}$ take the form

$f_{l}={\frac {2l+1}{2ik}}(e^{2i\delta _{l}}-1),$

where $\delta _{l}$ denotes the scattering phase shift. At low temperatures, we are interested in how the scattering amplitude behaves in the limit $k\to 0$. In particular, for the van der Waals interaction we have

$\delta _{0}=-ka_{s}$

$\delta _{1}\sim k^{3}$

$\delta _{l}\sim k^{4}\;\;\;l>1$

Here, we have introduced the $s$-wave scattering length $a_{s}$, in which the details of the scattering process have been absorbed. The typical scale is given by the effective range $r_{e}$

$a_{s}\sim r_{e}=(2\mu C_{6})^{1/4},$

however, close to resonances in the molecular channel (so-called Feshbach resonances), the scattering length will substantially deviate from this typical scale . Consequently, at low $k$, we can expect $s$-wave scattering to be the most relevant term. Then, we can use the analytic expression for the scattering amplitude

$f(k)={\frac {1}{-1/a_{s}+r_{e}k^{2}/2-ik}}.$

As long as the energy scales are lower than the characteristic scale $1/(2\mu r_{e})^{2}$, we can neglect the contribution from the effective range and are in the regime where the scattering is completely described by the scattering length (typical temperatures: $<1\,\mathrm {mK}$). This is the regime of ultracold atoms.

We are now looking for an interaction potential in position space that reproduces this scattering behavior. This is achieved by the pseudopotential

$V({\textbf {r}})={\frac {4\pi a_{s}}{2\mu }}\delta ({\textbf {r}}){\frac {\partial }{\partial r}}r.$

However, the term involving the partial derivative is only important when the wave function is singular at the origin and can be neglected in most cases. Then, the interaction term in the many-body Hamiltonian simplifies to

$H_{\mathrm {int} }={\frac {1}{2}}\int \mathrm {d} r\mathrm {d} r'\psi ^{\dagger }(r)\psi ^{\dagger }(r')V_{\mathrm {int} }(r-r')\psi (r')\psi (r)={\frac {g}{2}}\int \mathrm {d} r\psi ^{\dagger }(r)\psi ^{\dagger }(r)\psi (r)\psi (r).$

Note that for fermions this term vanishes due to the Pauli exclusion principle. In this case, the dominant scattering occurs in the much weaker $p$-wave channel.

Optical lattices

The ground state of the bosonic many-body Hamiltonian with a contact interaction is a Bose-Einstein condensate, i.e., even in the presence of strong interactions, the $k=0$ mode gets macroscopically occupied, i.e.,

$\lim _{N\to \infty }\langle a_{0}^{\dagger }a_{0}\rangle /N=n_{c}>0,$

where $n_{c}$ is the condensate fraction. Closely related is the concept of off-diagonal long-range order (ODLRO),

$\lim _{r\to \infty }\langle \psi ^{\dagger }(r)\psi (0)\rangle =n_{c},$

from which we see that a Bose-Einstein condensate is characterized by phase coherence.

If we want to use an ultracold Bose gas as a quantum simulator for typical many-body problems arising in condensed matter physics, we should put the atoms into a periodic potentials similarly how it is done for electrons in a solid state lattice provided by the atomic nuclei. Such periodic structures can be efficiently realized for ultracold atoms using laser potentials. In particular, we consider near-resonant light acting on Rb atoms in their electronic ground state. Then, the laser will create a perturbation of the form

$V=dE_{0}\cos(\omega t-k_{L}x)+\mathrm {h.c.} .$,

where $dis$ the dipole operator, $E_{0}$ is the strength of the electric field created by the laser, $\omega$ is the laser frequency, and $k_{L}$ its wavevector. As previously, we assume that only a single transition to one of the $5p$ is relevant. Then, the potential can be cast onto a two-level problem of the form

$V=\Delta |5p\rangle \langle 5p|+[\Omega \cos(\omega t-k_{L}x)|5s\rangle \langle 5p|+\mathrm {h.c.} ],$

where we have introduced the Rabi frequency $\Omega =d_{5s-5p}E_{0}$. We now go into the rotating frame of the laser field, i.e., we make the transformation $|5p\rangle \to |5_{p}\rangle \exp(i\omega t)$. Inserting this into the Schrödinger equation shifts the $5p$ level by the frequency $\omega$ and leads to the detuning $\delta =\Delta -\omega$. In the rotating-wave approximation, we neglect fast rotating terms on the order of $2\omega$ and obtain the effective Hamiltonian

$H=\delta |5p\rangle \langle 5p|+\left[{\frac {\Omega }{2}}\cos(k_{L}x)|5s\rangle \langle 5p|+\mathrm {h.c.} \right].$

If the detuning is much larger than the Rabi frequency, we can integrate out the excited state in second-order perturbation theory and obtain the effective potential for the ground state atoms,

$V(r)={\frac {\Omega ^{2}}{2\Delta }}\cos ^{2}(k_{L}x).$

We are now interested in the dynamics of ultracold atoms in such an optical lattice potential. We first look at the behavior of a single atom. The laser potential is periodic, hence the solution will be given by Bloch waves of the form

$\phi _{nk}(r)=e^{ikr}u_{nk}(r),$

where $u_{nk}(r)$ are periodic functions. The Bloch wavefunctions are characterized by a quasi-momentum $k$, restricted to the first Brillouin zone of the reciprocal lattice, and a band index $n$. Dealing with Bloch waves is often inconvenient, however, therefore it is usually better to work with Wannier functions $w_{n}(r-R_{i})$, which are well localized in space around a particular lattice minimum $R_{i}$. They are given by the Fourier transform of the Bloch waves over the first Brillouin zone,

$w_{n}(r-R_{i})=\int {\frac {\mathrm {d} k}{v}}\phi _{nk}(r),$

where $v$ is the volume of the first Brillouin zone. Note that the definition of the Wannier functions is not unique, as the Bloch waves are defined only up to an arbitrary phase. We can, however, require the orthonormality relation,

$\int \mathrm {d} rw_{n}^{*}(r-R_{i})w_{m}(r-R_{j})=\delta _{mn}\delta _{ij}.$

We can also express the field operators in terms of the Wannier functions,

$\psi (r)=\sum _{R_{i},n}w_{n}(r-R_{i})a_{R,n}.$

Here, the operator $a_{R,n}$ annihilates are particle in a Wannier state corresponding to the lattice site $R_{i}$ and the Bloch band $n$. For reasonably deep lattices, the bands are well separated, and we may restrict our analysis to the lowest Bloch band.

Assuming separable lattice potentials, the problem reduces to the one-dimensional Schrödinger equation

$\left(-{\frac {1}{2m}}{\frac {\mathrm {d} ^{2}}{\mathrm {d} x^{2}}}-{\frac {V_{0}}{2}}\cos 2k_{L}x\right)\psi (x)=E\psi (x).$

Using the substitutions $v=k_{L}x$ and $q=V_{0}/(4E_{r})$, where $E_{r}=k_{L}^{2}/2m$ is the recoil energy of the laser, this can be mapped onto a Mathieu equation,

${\frac {\mathrm {d} ^{2}}{\mathrm {d} x^{2}}}y(v)+(E-2q\cos 2x)y(v)=0.$

The Mathieu equation has symmetric solutions with eigenvalues $a_{r}$ and antisymmetric solutions with eigenvalues $b_{r}$. Here, the lowest eigenvalue in the lowest band corresponds to the eigenvalue $a_{0}$, while $b_{1}$ is the eigenvalue of the highest state in the first band. Then, in the limit of deep optical lattices ($q\gg 1$) we can use an (unproven) relation to obtain the width $W$ of the first Bloch band ,

$W=b_{1}-a_{0}={\frac {16}{\sqrt {\pi }}}\left({\frac {V_{0}}{E_{r}}}\right)^{3/4}e^{-4{\sqrt {V/E_{r}}}}.$

In the lowest band, the dispersion relation is given by

$E_{k}=-2J\cos ka$

where $a=\pi /k_{L}$ is the lattice spacing and the characteristic energy scale $J$ is related to the bandwidth by $W=4J$. We can also express the Hamiltonian using the Wannier basis,

$H=\sum _{R_{i},R_{j}}J(R_{i}-R_{j})a_{R_{i}}^{\dagger }a_{R_{j}},$

where the function $J(R_{i}-R_{j})$ describes how the tunneling matrix elements depend on the distance between the lattice sites. It is important to note that the Wannier functions are well localized and decay exponentially on larger distances , so the dominant contribution is given by the hopping between adjacent lattice sites for sufficiently deep lattices. Then, the strength of this nearest-neighbor hopping is simply given by $J$.

The Bose-Hubbard model

Let us now look at the effect of interactions. Here, we assume that the interaction is weaker than the separation to the first excited Bloch band, i.e., we can still reduce our analysis to the lowest band. Additionally, we restrict the effect the effect of interactions to occur within the same lattice site, this is justified if the lattice is deep enough the Wannier functions decay rapidly. Then, we can approximate the optical lattice potential by a harmonic confinement and the Wannier function reduces to the ground state wave function of the harmonic oscillator,

$w(x)\approx {\frac {1}{\sqrt {{\sqrt {\pi }}a_{0}}}}e^{-x^{2}/(2a_{0}^{2})},$

with the oscillator length being given by

$a_{0}=\left({\frac {1}{4m^{2}V_{0}E_{r}}}\right)^{1/4}.$

The matrix element for the contact interaction on a lattice site is then given by <ref="Bloch2008" />

$U=g\int \mathrm {d} r|w(r)|^{4}={\sqrt {\frac {8}{\pi }}}k_{L}a_{0}E_{r}\left({\frac {V_{0}}{E_{r}}}\right)^{3/4}.$

Note that this interaction occurs for each pairwise combination of two-particles, so for $n_{i}$ particles on lattice site $i$, this interaction has to be weighted by a factor of $n_{i}(n_{i}-1)/2$.

Finally, we work in the grand-canonical ensemble, so there will be a chemical potential $\mu$ associated with the creation of a particle. This chemical potential can be site-dependent, in particular if we take the shape of the (typically harmonic) trapping potential into account. Using the localized Wannier states in the lowest Bloch band, the many-body Hamiltonian then takes the form

$H=-J\sum _{\langle ij\rangle }\left(a_{i}a_{j}^{\dagger }+\mathrm {h.c.} \right)+{\frac {U}{2}}\sum _{i}n_{i}(n_{i}-1)-\mu \sum _{i}n_{i}.$

This model is called the Hubbard model and is one of the most studied models of condensed matter physics. Here, we have found a way to quantum simulate its bosonic variant with ultracold atoms.

Let us now turn to the ground state phase diagram for the Bose-Hubbard model. For $J=0$, the Hamiltonian is classical and simplifies to a sum of on-site problems. Depending on the parameter $\mu /U$, there are different possible integer fillings per site, with phases having filling factor $n$ satisfying the relation

$n-1\leq {\frac {\mu }{U}}\leq n.$

If $J$ is nonzero but small, we can do perturbation theory in $J$. Then, in second order, we will create virtual particle or hole excitations. These excitations can move around, but essentially remain confined to the site where they originated. Consequently, this regime of the Bose-Hubbard model is an insulating phase (a Mott insulator).

On the other hand, in the limit of large $J$, the bosons will simply accumulate at the bottom of the single-particle band and form a Bose-Einstein condensate there. Therefore, there has to be a quantum phase transition in between separating these two regions.

To understand the behavior close to this insulator-superfluid transition, we make use of a variational ansatz that accounts for small fluctuations around the Mott insulator with a lattice filling $n^{*}$, taking the form of a product wave function ,

$|\psi \rangle =\prod _{i}|\psi _{0}\rangle _{i}.$

In terms of the localized states, the on-site wave function takes the form

$|\psi _{0}\rangle ={\sqrt {1-\varepsilon }}|n^{*}\rangle +{\sqrt {\varepsilon }}|n^{*}\pm 1\rangle ,$

where $\varepsilon$ is a variational parameter. As the wavefunction treats fluctuations increasing or decreasing the particle number equally, this ansatz is only correct when the average density satisfies $\langle n\rangle =n^{*}$, i.e., at the center of the Mott phases given by $\mu /U=n^{*}-1/2$. Then, the expectation value of the hopping reduces to

$\langle \psi |a_{i}a_{j}^{\dagger }|\psi \rangle =(\langle \psi _{0}|a|\psi _{0}\rangle )^{2}=(\langle \psi _{0}|{\sqrt {1-2\varepsilon }}{\sqrt {n^{*}}}|n^{*}-1\rangle +{\sqrt {\varepsilon }}{\sqrt {n^{*}+1}}|n^{*}\rangle )^{2}=\varepsilon (1-2\varepsilon )({\sqrt {n^{*}}}+{\sqrt {n^{*}+1}})^{2}.$

For a lattice with $z$ nearest neighbors the total variational energy per lattice site is given by

$\langle H\rangle /N=-Jz\varepsilon (1-2\varepsilon )({\sqrt {n^{*}}}+{\sqrt {n^{*}+1}})^{2}+U/2[2\varepsilon +n^{*}(n^{*}-1)]-\mu n^{*}.$

After minimization with respect to $\varepsilon$, we find a critical value $U_{c}$, at which $\varepsilon$ vanishes (and remains exactly zero for larger values of $U$), given by

${\frac {U_{c}}{zJ}}=({\sqrt {n^{*}}}+{\sqrt {n^{*}+1}})^{2}.$

This marks the quantum phase transition between the Mott insulator and the superfluid. For $n^{*}=1$, we find $U_{c}/zJ=3+{\sqrt {2}}=5.828$. For a three-dimensional square lattice ($z=6$), quantum Monte-Carlo results predict $U_{c}/zJ=4.89$ , i.e., our simple variational ansatz predicts a value that is only 20% too large. The overestimation of the superfluid phase is also quite natural: our variational wave function neglects any correlations between the sites and thus corresponds to a mean-field approximation . Mean-field treatments are well-known to overestimate ordered phases such as superfluid as the order can be destroyed by the quantum fluctuation around the mean field that are being neglected. In lower dimensions, quantum fluctuations are even more important and the performance of mean-field theories becomes worse.

To map out the shape of the phase diagram in more detail, let us investigate what happens at small $J$ near the boundary between two Mott phases. There, only the two states $|n^{*}\rangle$ and $|n^{*}+1\rangle$ will be relevant. We use an ansatz of the form

$|\psi _{0}\rangle ={\sqrt {1-\alpha }}|n^{*}\rangle +{\sqrt {\alpha }}|n^{*}+1\rangle .$

Depending on our choice of $\mu /U$, either $\alpha$ or $1-\alpha$ will be a small quantity for small enough values of $zJ/U$. In the latter case, we find the lower branch of the Mott lobe as

${\frac {\mu _{c}}{U}}=n-1+{\frac {2nzJ}{U}},$

while the other case results in the upper branch given by

${\frac {\mu _{c}}{U}}=n-{\frac {2(n-1)zJ}{U}}.$

The entire phase diagram can then be constructed by interpolation between these known limits.

Time-of-flight images of a two-dimensional Bose gas. For shallow lattices, the pronounced diffraction peaks indicate that the gas is in a superfluid phase. Increasing the lattice depth leads to stronger interactions and eventually results in the formation of a Mott insulator (taken from ).

Experimentally, the Mott transition can be probed in a time-of-flight experiment. For this, the trapping potentials are switched off and the atomic cloud is allowed to expand in free space. Then, an absorption image of the expanded cloud will correspond to the initial momentum distribution of the atoms. In the superfluid, we expect a pronounced peak at $k=0$, which gets repeated when going to a higher-order Brillouin zones. In the Mott insulator, however, there is no preferred momentum so we expect a very broad distribution without any distinct features, see Fig. 1. The first experimental quantum simulation of the Mott insulator to superfluid transition were carried out in 2002 by Greiner et al., where they found a critical value of $U/zJ\approx 6$ .

Fermi-Hubbard model

While the bosonic variant of the Hubbard model is well understood, this is not the case for the fermionic version, especially since the sign problem makes quantum Monte-Carlo calculations prohibitively hard. Using ultracold atoms, quantum simulation of the Fermi-Hubbard model can be achieved in a very analogous manner. However, since fermions cannot interact via an on-site interaction due to the Pauli principle, it is necessary to include a hyperfine spin degree as well so that fermions in different spin states can undergo a contact interaction, which is also characterized by its $s$-wave scattering length. Including spin, the Fermi-Hubbard model reads

$H=-J\sum _{\langle ij\rangle \sigma }\left(a_{i,\sigma }a_{j,\sigma }^{\dagger }+\mathrm {h.c.} \right)+U\sum _{i}n_{i\uparrow }n_{i\downarrow }-\mu \sum _{i\sigma }n_{i\sigma }.$

Despite its hardness, several things can be observed about the Fermi-Hubbard model. In the case of half filling and in the limit of $J\ll U$, we can employ perturbation theory in $J$. At $J=0$, all the sites are filled with a single fermion, and there are no double occupancies anywhere in the system. However, for finite $J$, we can virtually create holes $|0\rangle$ and doublons $|D\rangle$. Due to the Pauli principle, this can only happen when the spins are opposite on the neighboring sites. In lowest order, the following terms are relevant:

$\langle \uparrow \downarrow |a_{i}a_{i+1}^{\dagger }|D0\rangle \langle D0|a_{i}^{\dagger }a_{i+1}|\uparrow \downarrow \rangle =-1$

$\langle \uparrow \downarrow |a_{i}a_{i+1}^{\dagger }|D0\rangle \langle D0|a_{i}^{\dagger }a_{i+1}|\downarrow \uparrow \rangle =1.$

Then, the model can be represented by local spin-1/2 degrees of freedom and reduces to the antiferromagnetic Heisenberg model,

$H={\frac {zJ^{2}}{U}}\sum _{\langle ij\rangle }{\textbf {S}}_{i}{\textbf {S}}_{j}.$

On the two-dimensional square lattice, there is no sign problem, and quantum Monte-Carlo results have found the existence of an antiferromagnetically ordered ground state . Away from half-filling, the situation is much more subtle, however, there is the striking possibility of the appearance of an unconventional $d$-wave superconducting phase, which might explain the mystery of high-temperature superconductivity in copper oxides .

The status of experiments with ultracold fermions in optical lattices is not as far developed as for bosons, as the cooling process is more challenging. In addition, the energy scales associated with antiferromagnetic order or even $d$-wave superconductivity are a fraction of $J^{2}/U$, which requires to cool down to extremely low temperatures. In 2008, two groups achieved the realization of a Mott insulating phase with fermions .