# Quantum Monte Carlo simulations of antiferromagnetism in ultracold fermions

on optical lattices within real-space dynamical mean-field theory

###### Abstract

We present a massively parallel quantum Monte Carlo based implementation of real-space dynamical mean-field theory for general inhomogeneous correlated fermionic lattice systems. As a first application, we study magnetic order in a binary mixture of repulsively interacting fermionic atoms harmonically trapped in an optical lattice. We explore temperature effects and establish signatures of the Néel transition in observables directly accessible in cold-atom experiments; entropy estimates are also provided. We demonstrate that the local density approximation (LDA) fails for ordered phases. In contrast, a “slab” approximation allows us to reach experimental system sizes with atoms without significant loss of accuracy.

###### keywords:

dynamical mean-field theory, quantum Monte Carlo, ultracold fermions, optical lattices, antiferromagnetism###### Pacs:

67.85.-d, 03.75.Ss, 71.10.Fd, 71.30.+h^{†}

^{†}journal: Computer Physics Communications

Starting with the achievement of Bose Einstein condensation, ultracold atomic gases have led to fascinating insights in quantum many-body phenomena Bloch08_RMP (). With the recent experimental successes in realizing quantum degenerate Koehl05 () and strongly interacting Fermi gases Esslinger08_Nature (); Bloch08_ferm () in optical lattices, such systems are considered as highly tunable quantum simulators of condensed matter Bloch08_Science (). The recent observation of the fermionic Mott transition in binary mixtures of on cubic lattices Esslinger08_Nature (); Bloch08_ferm () marks important progress in this respect. The next experimental challenge, taken up by many of the leading groups, is to reach and identify an ordered antiferromagnetic (AF) Néel phase.

Ultracold quantum gases on optical lattices have several advantages in comparison to solid state systems; however, the detection methods established for solids do not always find correspondence in cold-atom experiments. In particular, the detection of the AF order parameter is not straightforward. Therefore, it is important to identify fingerprints of AF phases that are directly accessible experimentally.

Due to the intrinsic inhomogeneity of trapped atomic clouds, even the qualitative interpretation of experimental data may depend on corresponding quantitative simulations. The dynamical mean-field theory (DMFT) is well established as a powerful, nonperturbative approach to interacting Fermi systems Georges96 (); Kotliar_Vollhardt (). Spatial inhomogeneities of the optical lattice can be captured within a real-space extension of the method (RDMFT) Snoek_NJP08 (); Helmes_PRL08 () or, approximately, by applying DMFT within a local density approximation (LDA). In either case, the numerical accuracy depends on the method chosen as DMFT impurity solver. Quantum Monte Carlo (QMC) based solvers are precise and even numerically cheap at or above the temperature () ranges relevant for AF ordering and accessible in cold-atom experiments.

In this paper, we introduce our Hirsch-Fye QMC Hirsch86 () based RDMFT implementation, with focus on computational aspects. Then, we turn to a first application, establishing the essential physics related to AF order in a relatively small system involving atoms. In order to minimize finite-size effects, we choose a square lattice geometry for which finite-temperature antiferromagnetism is not really physical; however, due to the mean-field character of the DMFT, the results are representative of three-dimensional cubic systems, up to a change in energy scales. Using selective simulations with atoms, we establish that the properties of large three-dimensional clouds can be accurately extracted from simulations for central two-dimensional slabs and that simulation boxes can be chosen smaller than one would naively expect.

Model and Methods – Binary mixtures of equivalent fermions in an optical lattice are well described by the Hubbard model [with (isotropic) trapping potential ],

(1) |

Here, , () are annihilation (creation) operators for a fermion with (pseudo) spin at site (with coordinates ), is the hopping amplitude between nearest-neighbor sites , is the on-site interaction, and is the chemical potential. We choose and use as the energy unit.

The RDMFT approach Snoek_NJP08 (); Helmes_PRL08 () is based on the following expression for the Green function matrix for Matsubara frequency (with spin indices suppressed):

(2) |

here the only approximation is the DMFT assumption of a local (i.e. site-diagonal) self-energy , which corresponds to a momentum independence of in translation-invariant systems Kotliar_Vollhardt (). The standard DMFT impurity problem, one for each inequivalent lattice site, is solved in this work using the Hirsch-Fye QMC algorithm Hirsch86 (). The discretization of the imaginary-time interval introduces a systematic error which can be eliminated by extrapolations Bluemer05 (); Bluemer07 (). In this work, we choose unless indicated otherwise.

Note that (2) only couples and at the same Matsubara frequency; thus, the matrix inversions can be performed in parallel for the 500 Matsubara frequencies explicitly taken into account. In contrast, the impurity equations couple different frequencies (and imaginary times), but are local in the site indices. Therefore, this part of the self-consistency problem can be performed in parallel for all inequivalent impurities. In addition, the importance sampling for each impurity may be distributed (via MPI) over some processes. In total, the program can saturate 500 CPU cores for large enough problems.

Results for square lattice – For the square lattice, the noninteracting dispersion extends from to , leading to a bandwidth of . The chosen interaction is already in the strongly correlated regime: at half filling, the fermions would always be localized. However, in the inhomogeneous trapped system, an outer shell has to remain delocalized due to the outward decay of the filling to zero. Consequently, staggered magnetic order can only be expected in the core region.

This is exactly what is seen in Fig. 1 (first row): at low (left column), the core shows a nearly perfect staggered magnetization. With increasing (from left to right), both the extent of the ordered region and its polarization decrease until the order is lost at the Néel temperature . Unfortunately, this most obvious signature of AF order is not directly accessible experimentally, primarily due to lack of single-site resolution. A measurement of particle density profiles (third row) also wouldn’t help in this case as they are practically unchanged (at this scale) through the transition. In contrast, the double occupancy, the probability of two particles occupying the same lattice site, provides a distinct signal: at high , it is featureless in the center, with a maximum value of about . Only at low temperatures, it is enhanced, by up to , in the developing central antiferromagnetic core. For measurements of the double occupancy, accuracies of about have already been established Esslinger08_Nature (); thus, only a minor refinement is needed to resolve this AF signal (when low enough can be reached).

Experimentally, the fermionic systems are prepared without an optical lattice and with negligible interactions; the initial temperature can, therefore, be determined by fitting the density profile with a Fermi distribution function. Then, the lattice is switched on; due to the localization, this also induces a contact interaction. Clearly, the temperature of a final equilibrium state will differ greatly from the initial value (in an unknown way); however, the process is essentially adiabatic. Therefore, only the entropy density is a reliable temperature scale.

Monte Carlo importance sampling cannot directly access entropy (or free energy) information. Thus, we revert to the thermodynamic relation and obtain the entropy density from integrating over temperature derivatives of particle densities (using the fact that corresponds to an empty system with zero entropy). The result is shown in Fig. 2: after a steep initial increase, the central entropy density (diamonds) develops a kink around (arrow) and then reaches a plateau with a value slightly above the mean field value of for a spin system. The average entropy density (circles) is always larger due to the low-density shell; in particular, the kink is much weaker and a significant slope remains above . Corresponding entropy data for smaller discretization (crosses) show excellent agreement, which establishes that systematic errors are insignificant. We should point out that fine grids in and are essential for this study since strong peaks in the integrand (see inset of Fig. 2) have to be resolved.

One might ask whether expensive RDMFT calculations as described above are really necessary for properly resolving the ordering phenomena. After all, in the case of the paramagnetic Mott metal-insulator transition, the LDA, which approximates the properties of each site by those of a homogeneous system with the same local potential, yields static observables in nearly perfect agreement with RDMFT Helmes_PRL08 (). However, as seen in Fig. 3, the order parameter profiles show enormous deviations from LDA even for the largest systems, which fully justifies the larger numerical effort of the RDMFT approach.

Efficient simulations for three dimensions – Let us, finally, leave the “toy” case of small two-dimensional systems and discuss strategies for realistic simulations of current experimental setups involving cubic lattices and particles. While such system sizes do not pose fundamental difficulties on the impurity part of the self-consistency problem, due to linear scaling and perfect parallelism, a full matrix inversion (with cubic scaling in the number of sites) as required in (2) is out of question, at least for dense-matrix algorithms. In our current implementation and for typical numerical accuracies, the cost of the matrix inversion dominates beyond some 1000 lattice sites.

Fortunately, the system sizes that have to be simulated explicitly can be reduced quite drastically. First of all, periodic boundary conditions should be used even for the inhomogeneous trapped systems. Then, the simulation box can be chosen much smaller than the atomic cloud provided that the AF core fn:isotropy () does not touch the boundaries, as shown in Fig. 4: while spikes are visible (at radial distances of about ) in for a system (lower panel, short-dashed lines) at , the data (long-dashed lines) are already smooth at this scale, indicating that finite-size effects become negligible. At , closer to the Néel temperature , the accuracy in improves further due to the smaller core.

Even larger savings are possible by restricting the simulation to a central slab, with periodic boundary conditions in the perpendicular direction. In practice, we neglect the (small) perpendicular component of the trapping potential in the slab, i.e., replace the spherical potential by a cylindrical one, which is an excellent approximation for large enough systems. Then, a thermodynamic limit is reached, within numerical accuracy, already at layer thickness of 4. Corresponding results, indicated by solid and dotted lines in Fig. 4, show negligible size dependence and agree well with the direct 3D results, independent of the temperature (density profiles for , nearly indistinguishable from those for , are omitted). In contrast, the LDA (dash-dotted lines) is completely off for .

In conclusion, we have implemented the RDMFT approach for inhomogeneous correlated Fermi systems on a QMC basis. In this work, we have shown, for a relatively small system, that the onset of antiferromagnetic order at low is signaled by an enhanced double occupancy. This effect and the impact of the interaction strenght are discussed on a more quantitative level for large cubic systems elsewhere Gorelik10 (), using the slab approximation established in this paper.

We thank W. Hofstetter, M. Snoek, and I. Titvinidze for valuable discussions and help in implementing RDMFT. Support by the DFG within the TR 49 and by the John von Neumann Institute for Computing is gratefully acknowledged.

## References

- (1) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
- (2) M. Köhl, H. Moritz, Th. Stöferle, K. Günter, T. Esslinger, Phys. Rev. Lett. 94, 080403 (2005).
- (3) R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature 455, 204 (2008).
- (4) U. Schneider et al., Science 322, 1520 (2008).
- (5) I. Bloch, Science 319, 1202 (2008).
- (6) A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- (7) G. Kotliar and D. Vollhardt, Physics Today 57, 53 (2004).
- (8) M. Snoek, I. Titvinidze, C. Töke, K. Byczuk, and W. Hofstetter, New J. Phys. 10, 093008 (2008).
- (9) R. W. Helmes, T. A. Costi, and A. Rosch, Phys. Rev. Lett. 100, 056403 (2008).
- (10) J. Hirsch and R. Fye, Phys. Rev. Lett. 56, 2521 (1986)
- (11) N. Blümer and E. Kalinowski, Physica B 359-361, 648 (2005).
- (12) N. Blümer, Phys. Rev. B 76, 205120 (2007).
- (13) Due to the uniform hopping, the AF ordering vector is .
- (14) E. V. Gorelik, I. Titvinidze, M. Snoek, W. Hofstetter, and N. Blümer, arXiv:1004.4857.