Quantum Monte Carlo simulations of antiferromagnetism in ultracold fermionson optical lattices within real-space dynamical mean-field theory

Quantum Monte Carlo simulations of antiferromagnetism in ultracold fermions
on optical lattices within real-space dynamical mean-field theory

N. Blümer E. V. Gorelik Institute of Physics, Johannes Gutenberg University, 55099 Mainz, Germany

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.

dynamical mean-field theory, quantum Monte Carlo, ultracold fermions, optical lattices, antiferromagnetism
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 ],


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.

Figure 1: RDMFT-QMC results for square lattice (, ). First row: as evidenced by the local magnetization, a large AF core is strongly polarized (up to 90%) for (left column); with increasing , both extent and magnitude of the AF order decay until a paramagnetic phase is reached for . Second row: the double occupancy is strongly enhanced in AF regions; in contrast, the density (third row) shows little dependence.

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):


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.

Figure 2: Main panel: RDMFT-QMC estimates of entropy for square lattice (, ); QMC data for the average entropy per particle (circles, for discretization ) and for the entropy of the central site (squares, also for ) are hardly distinguishable from corresponding estimates for (crosses). Inset: temperature derivative of the central density for representative values of the chemical potential .

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.

Figure 3: RDMFT-QMC estimates of the AF order parameter, i.e., the staggered magnetization, for square lattice (, , half filling at central site) for various trap stengths (symbols) as a function of the effective chemical potential . Even for nearly realistically weak trap strength (), the transition is much smoother, by proximity effects, than estimated from LDA (using single-site DMFT and QMC; dash-dotted line).

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 .

Figure 4: RDMFT-QMC estimates for cubic lattice (, ) using periodic boundary conditions, either for a full 3D lattice or using the slab approximation (see text). For the density (top panel), the slab results (dashed lines) are hardly distinguishable from the full calculations (solid and dotted lines); finite-size effects are very small. The agreement is also very good for the order parameter (bottom panel), even close to the Néel temperature , with the exception of spikes corresponding to edges of the 3D systems. Insets: central density and AF ordering pattern from slab calculation.

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.


  • (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.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description