Antiferromagnetic order of strongly interacting fermions in a trap: real-space dynamical mean-field analysis

We apply dynamical mean-field theory to strongly interacting fermions in an inhomogeneous environment. With the help of this real-space dynamical mean-field theory (R-DMFT) we investigate antiferromagnetic states of repulsively interacting fermions with spin- in a harmonic potential. Within R-DMFT, antiferromagnetic order is found to be stable in spatial regions with total particle density close to one, but persists also in parts of the system where the local density significantly deviates from half filling. In systems with spin imbalance, we find that antiferromagnetism is gradually suppressed and phase separation emerges beyond a critical value of the spin imbalance.


Introduction
Ultracold atoms in optical lattices provide a new laboratory for interacting quantum many body systems [1].Bosonic atoms in optical lattices realize the Bose-Hubbard model [2] and undergo a superfluid-Mott insulating transition when the potential depth of the optical lattice is increased [3].Recent experiments directly observed correlated particle tunneling [4] and superexchange [5], which are the basic mechanisms underlying quantum antiferromagnetism.Also the Fermi surface of fermionic atoms in optical lattices [6] and fermionic superfluidity of attractively interacting lattice fermions [7] were recently observed, bringing the realization of strongly correlated many-fermion states nearby [8].A two-component mixture of repulsively interacting fermions (e.g. 6 Li or 40 K) at half filling is predicted to form a correlated paramagnetic Mott insulator state above the critical (Néel) temperature and antiferromagnetic order below this temperature [9].Reaching this temperature is predicted to be within today's experimental capabilities [10,11].
In contrast to solid state systems, lattice defects, impurities and phonons are absent in optical lattices.However, the spatial inhomogeneity due to the harmonic confinement potential is always present, leading to a spatially varying local density.Therefore, the concept of long range order is questionable and ordered phases are expected to develop on finite length scales.A Hartree-Fock static mean-field theory predicts that antiferromagnetism, with staggered magnetization on a finite length scale, coexists with paramagnetic states in various spatial patterns, e.g.antiferromagnetism in the center of the trap or antiferromagnetism in a ring surrounded by a particle-or a hole-doped atomic liquid [12].On the other hand, both commensurate and incommensurate spin-densitywaves have been predicted for the hole-doped Hubbard model [13,14,15].However, the existence and properties of any ordered state on a finite length scale are strongly sensitive to quantum and thermal fluctuations.Therefore a theoretical description that captures effects of strong correlations and spatial inhomogeneity in a unified framework is needed.In this paper we apply a Real-Space Dynamical Mean-Field Theory (R-DMFT), which is a comprehensive, thermodynamically consistent and conserving mean-field theory for correlated lattice fermions in the presence of an external inhomogeneous potential.The R-DMFT takes into account local correlations exactly [16,17,18,19].
We prove that for spin-1 2 lattice fermions with local repulsive interaction antiferromagnetic order exists at zero temperature when the harmonic potential is present.We find that antiferromagnetic order is stable in spatial regions with total particle density close to one, but persists also in parts of the system where the local density significantly deviates from half filling.We also show that for strong repulsion phase separation occurs in imbalanced mixtures, when the difference in the particlenumber of the spin components is large.For weaker repulsion a strong imbalance destroys the antiferromagnetic order, but does not lead to phase separation.These results are especially intriguing with respect to recent experiments on attractively interacting fermions with spin imbalance, which have led to still unresolved questions regarding the nature of the observed phase separation [20,21,22].

Model
Repulsively interacting fermions in an optical lattice almost perfectly implement the Hubbard Hamiltonian where n iσ = c † iσ c iσ , and c iσ (c † iσ ) are fermionic annihilation (creation) operators for an atom with spin σ at site i, J is the hopping amplitude between nearest neighbor sites ij , U > 0 is the on-site interaction, µ σ is the (spin-dependent) chemical potential and i is the harmonic confinement potential.Moreover we define μ ≡ 1 2 (µ ↑ +µ ↓ ) and ∆µ ≡ µ ↑ − µ ↓ .The parameters of this model are tunable in experiments by a change of the lattice amplitude and via Feshbach resonances [1].In the following, J = 1 sets the energy unit and we take the lattice constant to be a = 1.

Method
To obtain the ground state properties of this system, we apply a real-space extension of dynamical mean-field theory (DMFT) [16,17,18,19,23,24].Within R-DMFT the self-energy is taken to be local, which is exact in the infinite dimensional limit [23,24].However, in an inhomogeneous system it depends on the lattice site, i.e.Σ ijσ = Σ (i) σ δ ij , where δ ij is a Kronecker delta.Formerly, a similar scheme has been developped for systems with inhomogeneity in one direction [25].Only recently, models with full inhomogeneity have been investigated, in particular the Falicov-Kimball model [26,27], disorderded systems [28] and paramagnetic states of cold fermionic atoms [29].
In the R-DMFT method, the Hamiltonian is mapped onto a set of single site problems.The physics of lattice site i is described by the local effective action [18] which explicitly depend on the site index i.Here τ is imaginary time.The function ) is a local non-interacting propagator interpreted as a dynamical Weiss mean-field which simulates the effect of all other sites [18] and is determined selfconsistently as follows: firstly, given the local self-energies Σ (i) (σ, iω n ) obtained from solving the action (2), the interacting lattice Green's function follows from the Dyson equation in real-space representation where a boldface notation indicates that the quantities are matrices labeled by site indices i and j and ω n are the Matsubara frequencies.The non-interacting lattice Green's function is given by G where 1 is the unity matrix.The matrix elements J ij are hopping amplitudes for a given lattice structure and V ij = δ ij V i represents a spatially varying potential.Secondly, the diagonal elements of the lattice Green's function are identified with the interacting local Green's functions, i.e.
Finally, the Weiss mean-field is obtained from the local Dyson equation which closes the set of the self-consistency equations.
The most difficult step in this procedure is the solution of the local action (2).This step is, however, similar to the solution of the local action in a homogeneous DMFT calculation.The difference is that in the present case the Weiss field G (i) 0 (σ, τ ) is obtained via the Real-Space Dyson equation (3), which incorparates the effect of the spatial inhomogeneity.This implies that for the numerical solution of the local action we can use standard techniques, which have proven to be reliable and efficient.In the present manuscript we use the Numerical Renormalization Group (NRG) at T = 0 [33,34,35,36,37] to solve the single site problems.
In practice the self-consistent solution is obtained iteratively from the initial Weiss mean-fields G (i) 0 (σ, iω n ) which are chosen differently for different spin σ and lattice sites i.Then the solutions with staggered magnetization or phase separation are obtained naturally in contrast to the standard DMFT, where additional sublattice structure has to be added [18].
Within R-DMFT significantly larger systems can be investigated than those studied by quantum Monte Carlo [30,31,32], for which in two and three dimensions only homogeneous data are available.The computational effort scales polynomially with the number of lattice sites N within R-DMFT.The application of the real-space Dyson equation requires a sparse matrix inversion for each frequency, which scales as N 3/2 .The number of NRG calculations per R-DMFT-run is linear in N, but can be kept small due to symmetries.Moreover, the solution of the real-space Dyson equation can be parallelized over the frequencies and the NRG-calculations can be parallelized over the lattice sites.

Results
We apply this method to spin-1 2 fermions in a two-dimensional square lattice with harmonic confinement.In the context of cold atoms a two dimensional system can be realized by applying a highly anisotropic optical lattice, which divides the system into two-dimensional slices.Although not exact, R-DMFT is expected to be a good approximation for the two-dimensional situation at zero temperature, since the derivation of the DMFT equations is controlled by the small parameter 1/z = 1/4 on the square lattice.

Balanced mixture
First we consider the case of an equal mixture of spin-up and down atoms: N ↑ = N ↓ .We find that antiferromagnetic order is stable in the presence of the inhomogeneous harmonic potential.In Fig. 1 we present examples for the spatial dependence of the magnetization at different strengths of the confining potential and the chemical potentials.In the case that the lattice at the center of the trap is half-filled, antiferromagnetism appears in the center of the system (Fig. 1a).When the particle density in the center of the trap is higher, antiferromagnetic order forms in a ring enclosing a paramagnetic region (Fig. 1b).These results are particularly important for ongoing attempts to realize antiferromagnetic states in optical lattices.Namely, we predict that the observation of antiferromagnetic order does not critically depend on the number of atoms in the system.For sufficiently strong repulsion between the particles, the necessary condition to find antiferromagnetic order is to prepare the system such that the local filling factor approximates or exceeds one in at least part of the system.We find no evidence for phase separation or a paramagnetic insulating boundary layer for the N ↑ = N ↓ case.
The antiferromagnetic ground state of homogeneous fermions described by the Hubbard Hamiltonian (1) without trap is stable when the density of particles varies from n ≈ 0.8 to 1.2, depending on the interaction value U [38].On the contrary, in the presence of the external harmonic potential, antiferromagnetic order appears for much lower or higher local total densities.Indeed, in Fig. 2 we present examples of the local density n i and the local magnetization m i = 1 2 (n i↑ − n i↓ ) as a function of distance from the center (main panel) and along a cut through the system (inset) which proves that antiferromagnetic order extends from the center of the trap and disappears only when n i ≈ 0.5 in Fig. 2a).Similarly, Fig. 2b) shows that antiferromagnetic order is stable on a ring when the local density extends between 0.5 n i 1.5.
We also determine the local density and the local magnetization within the Thomas-Fermi approximation (TFA) to R-DMFT, where the external potential is only included by a spatially varying chemical potential [39].The agreement between the full R-DMFT and TFA results is very good in regions well within or outside the antiferromagnetic domain.Encouraged by this, Fig. 3 shows additional TFA+R-DMFT profiles that can be used to compare R-DMFT with experiments for realistic systems in two and three dimensions.However, the staggered magnetization decays abruptly within TFA as compared to the full R-DMFT solution, i.e. the Thomas-Fermi approximation to R-DMFT essentially reproduces results from the standard homogeneous DMFT, cf.Fig. 3.The wider stability regime of the antiferromagnetic order found within full R-DMFT is caused by a proximity effect; antiferromagnetic order is induced into parts of the systems where the local densities are too low to stabilize antiferromagnetism in the homogeneous case.On the other hand, the proximity of the paramagnetic state reduces the staggered magnetization when the local density is close to one.

Imbalanced mixture
We now proceed by investigating the imbalanced case [40], i.e.N ↑ = N ↓ .Imbalance between the two spin-components is induced by a nonzero chemical potential difference ∆µ = µ ↑ − µ ↓ , which corresponds to a magnetic field.In the experimental situation, the density imbalance can be highly tuned and is stable due to the suppression of spinflip scattering proccesses in cold-atomic gases.Representative results are presented in Fig. 4, where we plot the up-and down-component of the density along a cut through the system.Upon increasing the imbalance parameter ∆µ, we find suppression of the antiferromagnetic order and emergence of phase separation between the minority and majority species.The phase separation region starts to develop far away from the center of the trap at small ∆µ and gradually spreads toward the center.We thus find that the border of the anti-ferromagnetic domain is most sensitive to phase separation.This is indeed reasonable: the energy cost to polarize the antiferromagnetic state is the energy difference between an antiferromagnetic state and a ferromagnetic state.This is of the order J 2 /U, which small for the large interaction U considered here.The antiferromagnetic order is thereby more unstable for larger distances to the trap-center, because of the vicinity to the paramagnetic region.The energy cost to polarize the paramagnetic regime is higher, because in this case kinetic energy has to be paid, whereas in the anti-ferromagnetic domain the kinetic energy is already quenched because the particles are almost localized.Due to the proximity effect we find that the paramagnetic regime close to the insulating domain also gets phase-separated, which leads to a ringlike structure of the minority species.
At strong interaction, U = 10 in the case shown in Fig. 4, atoms with different spins ultimately tend to occupy different spatial regions to avoid the mutual interaction and the minority species is completely expelled from the trap center.At weaker interaction, however, we found that the imbalanced system still contains interpenetrating atoms with different spins and phase separation does not occur.This is shown in Fig. 5, where for ∆µ = 0.8 the antiferromagnetic order has completely disappeared, but the two spin components are still interpenetrating.The small oscillations in the component densities can be understood as Friedel oscillation due to the small size of the system.We note that in the case of imbalanced spin-mixtures the agreement between the TFA and the exact solution to R-DMFT is far less good than in the balanced case presented above.

Conclusions
In conclusion we used the R-DMFT to establish stability of antiferromagnetism for balanced fermionic spin- 1  2 systems in a trap and the appearance of phase separation for imbalanced mixtures.The antiferromagnetic order predicted here can be observed at low enough temperatures by Fourier-sampling of time-of-flight images via Raman pulses [41], by measuring spin correlation functions via local probes [42], probing noise correlations [43,44], polarization spectroscopy [45], and Bragg scattering [46].The effect of spatial inhomogeneity on these probes will be investigated within R-DMFT in future studies.Moreover, the R-DMFT scheme presented here opens up the possibility to study a variety of other strongly correlated systems in inhomogeneous environments.

Figure 2 .
Figure 2. Particle and spin density profiles determined within the exact solution of R-DMFT and within the Thomas-Fermi approximation (TFA) to R-DMFT.The main panel shows the local total density n i and the local magnetization m i as a function of a distance from the trap center.The inset shows n i and m i along the y = 1 2 line.The parameters in the a) and b) panel are the same as in the respective panels of Fig. 1.

Figure 3 .
Figure 3.Total density and staggered magnetization as a function of effective chemical potential obtained within the Thomas-Fermi approximation to R-DMFT for two dimensional square (main panel) and three dimensional cubic lattices (inset).Main panel: U = 10 (diamonds) and U = 20 (circles); inset: U = 30.