Magnetic structure of an imbalanced Fermi gas in an optical lattice

We analyze the repulsive fermionic Hubbard model on square and cubic lattices with spin imbalance and in the presence of a parabolic confinement. We analyze the magnetic structure as a function of the repulsive interaction strength and polarization. In the first part of the paper we perform unrestricted Hartree-Fock calculations for the 2D case and find that above a critical interaction strength $U_c$ the system turns ferromagnetic at the edge of the trap, in agreement with the ferromagnetic Stoner instability of a homogeneous system away from half-filling. For $U<U_c$ we find a canted antiferromagnetic structure in the Mott region in the center and a partially polarized compressible edge. The antiferromagnetic order in the Mott plateau is perpendicular to the direction of the imbalance. In this regime the same qualitative behavior is expected for 2D and 3D systems. In the second part of the paper we give a general discussion of magnetic structures above $U_c$. We argue that spin conservation leads to nontrivial textures, both in the ferromagnetic polarization at the edge and for the Neel order in the Mott plateau. We discuss differences in magnetic structures for 2D and 3D cases.


I. INTRODUCTION
Cold atoms constitute a promising route to simulate model Hamiltonians of strongly correlated many-body physics with accurate control of system parameters [1,2]. After major experimental breakthroughs with ultracold bosonic atoms like the Bose-Einstein condensation (BEC) of alkali metal gases [3,4] or the observation of the superfluid-Mott-insulator transition in a bosonic Hubbard model [5], the field of ultracold atoms is currently addressing problems of strongly correlated fermionic systems [6][7][8][9]. Arguably, the most prominent goal is the understanding of the phase diagram of the fermionic Hubbard model, which is believed to be of major importance for high-temperature superconductivity [6,[10][11][12][13][14][15]. A twocomponent Fermi gas in an optical lattice is well described by the single-band Hubbard model whenever the energy gap to higher bands is much larger than on-site interaction, temperature, and chemical potential [1,2,11]. Only recently the fermionic Mott transition has been realized experimentally [16,17]. The major challenge for studying the magnetism of the fermionic Hubbard model is to reach temperatures below the Néel temperature [18,19]. In addition to the preparation of the antiferromagnetic state, characterization tools have to be developed to allow a clear identification of the magnetic structure. Possible experimental techniques include Bragg spectroscopy [20,21], local measurements of the magnetization [22,23], noise correlations [24,25], or the recently realized quantum-gas microscope [26].
The experimental control of spin imbalance in Fermi gases offered a unique way to study pairing phenomena beyond the standard BCS picture for attractive interactions [27,28]. Motivated by these results, we address in this work the effect of spin imbalance on the repulsive fermionic Hubbard model [29,30]. While we study strong optical lattices, where a singleband Hubbard model is realized, the magnetic structure of * bwunsch@physics.harvard.edu weak-to-intermediate lattice strength including multiple bands has also been discussed [31]. We find rich physics arising from the interplay between antiferromagnetic and Stoner ferromagnetic instabilities and spin imbalance.
The magnetic order of the two-dimensional (2D) repulsive Hubbard model has been extensively studied in the past (see Ref. [14]). Cold atoms in optical lattices differ in several ways from typical condensed-matter systems. First, there is a superposed external confinement potential, which divides the system in an incompressible Mott state in the center of the trap and in a compressible region at the edge. Second, the total spin is conserved, which means that we need to minimize the energy of the system given a global magnetization rather than a finite Zeeman field. One interesting problem concerns the spatial distribution of the imbalance between Mott plateau and edge, and it turns out that the solution strongly depends on the interaction strength. The constraint of spin conservation affects the ferromagnetic instability at the edge by enforcing nontrivial spin textures [32,33] which also affects the Néel order in the Mott plateau in the center, as we will discuss in Sec. IV.
In this work we study the repulsive fermionic Hubbard model including a parabolic confinement potential. In the first part of this work, we perform unrestricted Hartree-Fock calculations for the 2D case. Relevant physics for this system can be identified based on the mean-field phase diagram for the repulsive 2D homogeneous Hubbard model [34]. Up to a critical interaction strength U c , it predicts antiferromagnetic order close to half -filling and paramagnetic order elsewhere. In the spirit of a local density approximation, one might then expect that cold fermionic atoms in an optical lattice have antiferromagnetic correlations in spatial regions with one atom per site and are paramagnetic elsewhere. In order to account for a finite imbalance, the system has to change its magnetic structure. Using an unrestricted Hartree-Fock approach for the 2D system, we find a canted antiferromagnet in the Mott plateau in the trap center and a partially polarized edge. We note that canted antiferromagnetic order close to half-filling has been reported previously in Ref. [35]. With spin polarization along the z direction, the canted antiferromagnet accommodates the imbalance forming a constant z component of the local magnetization, and simultaneously it benefits from the superexchange interaction by building up an alternating magnetic order perpendicular to the z direction. Fixing the global imbalance and increasing the interaction strength results in more imbalance flowing to the edge.
Above a critical interaction strength U c , the unrestricted Hartree-Fock calculation predicts that the system turns ferromagnetic at the edge of the trap, in agreement with the ferromagnetic Stoner instability of a homogeneous system away from half-filling. Furthermore, the orientation of the antiferromagnetic order in the Mott plateau is perpendicular to the direction of the ferromagnet in the edge. Spin conservation has again a strong impact on the magnetic structure of the system since a uniformly polarized ferromagnetic edge together with an antiferromagnetic Mott plateau are generally not allowed. We will discuss spin textures in 2D and threedimensional (3D) lattices for U > U c , which fulfill spin conservation and which show the two prominent features predicted by the mean-field calculation, namely a) magnetic instabilities toward ferromagnetism in the compressible edge and antiferromagnetism in the Mott plateau and b) at the interface between Mott plateau and compressible edge, the orientation of the antiferromagnet and the ferromagnet are perpendicular to each other.
We are aware that the chosen mean-field approach generally overestimates symmetry breaking, and therefore the critical on-site interaction strength, U c , corresponding to the appearance of an intrinsic ferromagnetic edge, will presumably be higher than the one predicted here. However, intrinsic ferromagnetism away from half-filling is expected for sufficiently large interaction strength [32,36], and in fact, experimental indications for itinerant ferromagnetism in a Fermi gas of ultracold atoms have been reported recently in Ref. [37]. Given the tunability of the ratio between on-site interaction and nearest-neighbor hopping, U/t, the interaction strength required for the presented phase separation should be accessible in experiment (U/t = 150 have been reported in Ref. [16]).
This article is organized as follows. In Sec. II, we introduce the model, and in Sec. III, we calculate the magnetic structure for U < U c within an unrestricted Hartree-Fock approach. The topology of the intrinsically ferromagnetic edge arising for U > U c is addressed in Sec IV and in the Appendix. Finally in Sec. V, we summarize our findings and comment on the experimental significance of our results.

II. MODEL
We consider the fermionic single-band Hubbard model on a 2D and 3D cubic lattice with an external parabolic confining potential. The Hamiltonian is where σ ∈ {↑, ↓} labels the two fermionic components, which are the eigenstates of the z component of a spin algebra. These two components can either be the hyperfine state of the trapped fermions or even correspond to different atomic species. c iσ denotes the annihilation operator for a particle with spin σ at site i, whereas n iσ = c † iσ c iσ and n i = σ n iσ are the spin resolved and total occupation of site i. U is the on-site interaction, and t is the nearest-neighbor hopping. Finally r i denotes the distance of site i from the trap center measured in units of the lattice spacing a and α = mω 2 a 2 /2 characterizes the strength of the external confinement. The associated energy scale is the confinement strength at the edge of the atom cloud with one atom per site, denoted by V t . In 2D, V t = N α/π, where N is the particle number.

III. UNRESTRICTED HARTREE-FOCK APPROACH IN 2D
We now apply a Hartree-Fock mean-field decoupling in the spin and the density channel. Since the trap breaks translational invariance, the mean-field parameters will be site-dependent. Allowing for arbitrary spin and density at each site, we obtain the following mean-field Hamiltonian [38]: iα σ α,β c iβ )/2 denotes the spin operator at site i ( σ is the vector of Pauli matrices) and M i = S i is the local magnetization. Magnetization and density are determined self-consistently for fixed total particle number N . In the following, we assume zero temperature. The energy of the self-consistent solution is given by the sum over the lowest N single-particle energies of the Hamiltonian (2) plus the constant energy E 0 = U i ( M 2 i − n i 2 /4). An important subclass of self-consistent solutions are the ones with collinear magnetization where M y (i) = M x (i) = 0 on all sites. In particular, the generic phases of the homogeneous Hubbard model [34] have a collinear magnetization; either ferromagnetic M z (i) = M, antiferromagnetic M z (i) = (−1) i M, or paramagnetic M z (i) = 0. However, we will show that generally the combination of trapping potential and imbalance will lead to a non-collinear-magnetization profile.
We are interested in the ground state for a given imbalance, characterized by the polarization P = (N ↑ − N ↓ )/(N ↑ + N ↓ ), which is an experimentally controllable parameter [27]. The imbalance is conserved since the two components correspond to different internal states of the atoms (typically different hyperfine states) and transitions between these states are energetically forbidden unless they are driven by additional lasers. The single-particle eigenstates of the Hamiltonian in Eq. (2) only have well-defined spin if the magnetization is collinear. Generally, an expectation S z = 0 can be tuned by spin-dependent chemical potentials or equivalently by a fictitious magnetic field in z direction H z = −BS z .
The parabolic confinement will decrease the density away from the trap center. In a local density approximation, a cross section through the trap corresponds to a cut through the (n, U ) phase diagram at constant interaction U . Polarization can most easily be accommodated by ferromagnetism, but also antiferromagnetic and paramagnetic regions can account for finite imbalance. As discussed in the Introduction, in a canted antiferromaget, a spatially constant component aligned with the field is added to the alternating component perpendicular to the imbalance. The paramagnetic region can be partially polarized in the spirit of Pauli paramagnetism, where the polarization is proportional to the applied field. In the following, we show that canted antiferromagnetic order is realized at half-filling, and we study how the imbalance is distributed between Mott plateau and edge as a function of interaction and imbalance. Self-consistent solutions of the Hubbard model on the two-dimensional square lattice (2) have either a collinear or coplanar magnetization [35,38], and we can set M y = 0 without loss of generality. However, we note that enforcing vanishing global in-plane magnetization can lead to nontrivial three-dimensional topologies for the intrinsic ferromagnet [32,33], which will be discussed in Sec. IV.
A. The homogeneous system at half-filling Figure 1 shows the mean-field energies of canted and collinear solutions as a function of increasing imbalance for the homogeneous system at half-filling. A rough explanation of why the canted antiferromagnetic order is favored can be given within the mean-field Heisenberg model. Here the energy increases only quadratically with polarization for the canted order but linearly with polarization for collinear magnetization. Since the solutions are the same at the extreme values P = 0 and P = 1, the ground state is always a canted antiferromagnet.  particle number, N = 540, correspond to V t = 3.4t which is smaller than the on-site interaction so that double occupancies are absent.

B. Magnetization profile in the trap
Within the Mott plateau, we find canted antiferromagnetic order, as expected from the analysis of the homogeneous system. The cross sections of the spin resolved densities and the local magnetization in panels (b) and (c) of Fig. 2 show that the edge is partially polarized and does not have antiferromagnetic order, although the x component of the magnetization extends into the edge.
We now consider the distribution of a fixed imbalance for various on-site repulsions. Figure 3 illustrates that increasing interaction moves the imbalance to the edge. (We define the Mott plateau through |n i − 1| < 0.05.) Above a critical interaction strength (of order U c ≈ 10t), the edge is fully polarized and the Mott plateau is a pure antiferromagnet. The maximum in the majority density at the border of the Mott plateau can be understood by recalling that in the homogeneous system for strong interactions, there is a first-order phase transition between an antiferromagnet close to half-filling and a ferromagnet at finite doping [34]. By decreasing interactions below U c , the canting in the Mott plateau increases and the polarization at the edge decreases.
Next we describe the magnetic structure as a function of the global polarization, P , keeping the other parameters fixed. For U = 5t, the upper panel of Fig. 4 shows that both the polarization in the center with canted antiferromagnetic order and in the partially polarized edge increases linearly with the global polarization. The polarization at the edge is always larger than in the center until the Mott plateau disappears close to full polarization.
We now discuss the case of strong interaction (i.e., U > U c ). Here the edge is intrinsically ferromagnetic. As shown in Fig. 3, at U = 12t, the edge is already fully ferromagnetic in the absence of any fictitious magnetic field that is otherwise used to fix a certain global imbalance. Given the total number of atoms in the trap, N , and the number of atoms in the edge, N 01 , this defines a critical polarization P c = N 01 /N, which is P c ≈ 0.5 in Figs. 3 and 4. Our mean-field approach predicts for P < P c and U > U c a spatially uniform ferromagnetic edge with a direction other than the z direction. This implies a finite global in-plane magnetization. However, as we discuss in detail in the next section, such a solution which is forbidden by spin conservation and the preferred ferromagnetic order in the edge will have nontrivial spin textures for P < P c and U > U c . For now we restrict our discussion to P > P c and U > U c . Then the ferromagnetic order at the edge points in the z direction and the antiferromagnet in the Mott plateau is canted as shown in the lower panel of Fig 4. We now increase the number of particles so that the center of the trap is more than half-filled. In agreement with the symmetry of the homogeneous Hubbard model around half-filling, we find that the edge between the Mott plateau and double-occupied sites shows similar features as the outer edge discussed above. Figure 5 shows the magnetization profile and the spin-resolved densities. Here V t = 15.7 which is larger than the chosen on-site interaction. The Mott plateau is formed on a ring and has canted antiferromagnetic order. Moving away from the Mott ring, the antiferromagnetic order rapidly vanishes and the edge is strongly polarized. In fact, for this rather large value of U , we see a small maximum of the majority component at the outer edge and a minimum in the minority component at the inner edge.

IV. NONTRIVIAL SPIN TEXTURES FOR U > U c
The Hartree-Fock calculation predicts that above a critical interaction strength U c , the edge of the atom cloud turns ferromagnetic, even in absence of any imbalance or fictitious magnetic field. In the previous section, we defined a critical polarization, P c , corresponding to a fully polarized ferromagnetic edge along the z direction and an antiferromagnetic Mott plateau. In this section, we discuss qualitatively the magnetic structure for U > U c and P < P c .
A cold atom experiment is prepared from a paramagnetic state with no optical lattice. Controlling the imbalance between the two fermion species, the initial state is characterized by where P is the polarization and N the number of atoms. Since there is no coupling between the effective spin degree of freedom and the rest of the experimental system, the same constraints apply in the presence of an optical lattice and with strong on-site interaction U [32,33]. This additional constraint is always fulfilled in our mean-field treatment except for U > U c and P < P c , where a spatially uniform ferromagnetic edge is predicted with a direction other than the z direction. However, such a solution leads to a finite global in-plane magnetization, which is forbidden by the boundary condition. In order to fulfill Eq. (3), itinerant ferromagnetism in cold atom systems can have nontrivial topology as shown recently for balanced systems with filling factor less than unity everywhere [32,33].
In the following, we discuss the magnetic structure for U > U c and P < P c , both for 2D and 3D systems. We look for magnetic structures that fulfill spin conservation (3) and which show the two prominent features predicted by the mean-field calculation, namely a) magnetic instabilities toward ferromagnetism in the compressible edge and antiferromagnetism in the Mott plateau and b) at the interface between Mott plateau and compressible edge, the orientation of the antiferromagnet and the ferromagnet should be perpendicular to each other. Our qualitative analysis is based on the Ginzburg-Landau-type free energy functional (see Refs. [32] and [33]): where ρ is the positive stiffness constant, M 0 is the magnitude of the favored magnetization, and β > 0 determines the cost of amplitude fluctuations. The favored spin texture for strong interactions, U > U c , is determined by minimizing the total energy under the constraint of Eq. (3). In our qualitative analysis, we neglect that at the edge, the system parameters in Eq. (4) depend on the radius. This allows us to write the total energy of a spin structure, as a sum of three contributions: the energies of the spin structures at the edge, inside the Mott plateau, and at the interface of both regions. We note that the energy scale related with the spin structure of the ferromagnetic edge is of the order t and thus much bigger than the small superexchange t 2 /U that determines the spin structure in the Mott plateau. Therefore, we first minimize the free energy of the intrinsically ferromagnetic edge. The remaining two energy terms describe the interface between ferromagnetic and antiferromagnetic order at the edge of the Mott plateau and the energy of the spin structure in the Mott plateau. Based on the different scaling with the system size, we argue that the interface term dominates for large systems. While the interface term scales with r D−1 M , where r M is the radius of the Mott plateau and D denotes the dimension, the antiferromagnet scales like ln r M in 2D and like r M ln r M in 3D, as we will show. We minimize the interface term by choosing the orientation of the ferromagnetic and the antiferromagnetic order to be perpendicular to each other at the interface between Mott plateau and compressible edge. In the following, we discuss solutions, where the Mott plateau has no net imbalance. In fact, in the limit of large interactions U → ∞, the superexchange t 2 /U vanishes, so that one could allow for a strong polarization of the edge, by polarizing the Mott plateau in the opposite way. As estimated in Appendix B, such a solution is however higher in energy for realistic interaction strengths.

A. 2D lattice
We argue that (i) in presence of the Mott plateau a vortex structure for the ferromagnetic edge should be energetically 013616-5 favored as depicted in Fig. 6(a), and (ii) a finite imbalance should result in a vortex structure of the ferromagnetic order parameter in the xy plane together with a small z component (see Fig. 7). An important experimental consequence is a strong z component of the antiferromagnetic order in the center, which is aligned perpendicular to the ferromagnetic order in the edge. In Appendix A, we derive the energy of the different topological orders of the ferromagnetic edge: vortex, domain wall, or Skyrmion. These structures are illustrated in Fig. 6. It turns out that for realistic parameters, the vortex is lowest in energy. For finite imbalance, the edge will then be described by a ferromagnetic vortex in the xy plane and a constant z component. The energetically preferred direction of the antiferromagnetic order in the Mott plateau is perpendicular to that of the ferromagnet at the edge. The antiferromagnet in the Mott plateau will therefore have a small in-plane magnetization forming a vortex, which grows with increasing imbalance, and a strong z component as illustrated in Fig. 7.

B. 3D lattice
Similar arguments can be applied to a 3D system. Taking into account the boundary condition of vanishing global magnetization in balanced systems and by applying Eq. (4), one finds that the preferred structure of the ferromagnetic edge in a balanced system is a hedgehog [32,33]. As shown in Appendix A, the energetically preferred antiferromagnetic order in the center should then be either a planar vortex structure with M AF = ±M 0 e φ or a 3D spherical vortex M AF = M 0 e θ , where e φ = (− sin φ, cos φ, 0) and e θ = (cos θ cos φ, cos θ sin φ, − sin θ ) are spherical unit vectors.
Both solutions are illustrated in Fig. 8. They guarantee that at the edge of the Mott plateau, where the antiferromagnetic order of the center of the trap has an interface with the ferromagnet order at the edge, the orientations of the antiferromagnet and the ferromagnet are perpendicular to each other. A violation of this requirement would cost an energy that scales with the area of the interface r 2 M . Deformations of the perfect Néel order in the center of the trap, either in amplitude or phase, are minimized and the corresponding energy scales as r M ln(r M /a). For perfectly balanced systems, the vortex within the Mott plateau could lie in any plane. Imbalance will deform the hedgehog leading to a net z component (see Fig. 9). While this does not affect the energy of a vortex in the xy plane, While the ferromagnetic edge always has a hedgehog structure, the Mott plateau has either a planar vortex structure (a) or a 3D "spherical" vortex structure (b). While both structures have the same energy for the balanced system, the planar vortex is favored by finite imbalance (see Fig. 9).
FIG. 9. (Color online) Illustration of magnetic structures in 3D for an imbalanced system with U > U c and 0 < P < P c . Outer arrows indicate the magnetization at the ferromagnetic edge. Inner arrows illustrate the staggered magnetization in the antiferromagnetic Mott plateau. Note that finite imbalance only deforms the hedgehog structure of the ferromagnet edge, while the antiferromagnetic order in the Mott plateau is unchanged [see Fig. 8(a)].
it increases the energy for the vortices in other planes or for the spherical vortex. Therefore, we expect that for imbalanced systems in 3D with U > U c , the antiferromagnetic order in the Mott plateau will form a planar vortex structure in the xy plane as in Fig. 9. In contrast to the 2D case where we expect a strong z component of the antiferromagnetic order for U > U c and P < P c , we expect a vanishing z component in 3D.

V. DISCUSSION
In this work, we studied an interacting two-component Fermi gas on a 2D and 3D cubic lattice subject to a parabolic external confinement. We analyzed the magnetic structure as a function of the repulsive interaction strength and spin imbalance. Applying an unrestricted Hartree-Fock calculation for a 2D system, we identified the critical interaction strength U c where the edge turns ferromagnetic and analyzed the spatial distribution of a finite imbalance between the two Fermi components for U < U c . We found that the system has canted antiferromagnetic structure at half-filling with antiferromagnetic ordering in the plane perpendicular to the imbalance and is partially polarized elsewhere. Fixing the global imbalance and increasing the interaction strength results in more imbalance flowing to the edge. We expect the same qualitative behavior for 3D in that regime.
In the second part of the work, we gave a general discussion of the magnetic structure above U c both for 2D and 3D. We showed that spin conservation generally leads to nontrivial spin textures, both in the Mott plateau and at the edge. We predict that the edge has non-vanishing in-plane magnetization with a vortex structure in 2D and a hedgehog structure in 3D. We furthermore expect that for U > U c and small imbalance, the antiferromagnetic order in the Mott plateau has a finite z component in 2D, while in 3D a vanishing z component of the antiferromagnetic order in the Mott plateau is predicted.
We expect our findings to have clear experimental signatures if temperatures below the Néel temperature can be reached. A phase-contrast image [27] showing the density of each component separately can test our predictions of a Mott plateau with ferromagnetic borders. Detection of a canted antiferromagnet in the Mott plateau requires direct access to the order parameter. This can be achieved for instance through noise correlations [24] or by measuring the local magnetization [22,23,26]. Additionally, one can use Bragg spectroscopy [20,21] where the double-unit cell of the antiferromagnet results in additional Bragg peaks. Furthermore, the intensity of the additional Bragg peaks can then be used to measure the strength of the z component of the antiferromagnet.

APPENDIX A: GINZBURG-LANDAU THEORY
Following Ref. [32], we apply a Ginzburg-Landau-type description of the magnetism based on Eq. (4) to analyze the magnetic structure for U > U c , where the edge is intrinsically ferromagnetic. By enforcing a vanishing global in-plane magnetization, the ferromagnetic edge acquires nontrivial topology. For the energy estimate, we consider three energy contributions. The most relevant contribution is the free energy of the intrinsically ferromagnetic edge. Thereafter the contribution of the interface between ferromagnetic and antiferromagnetic order at the edge of the Mott plateau has to be taken into account, which is minimized by choosing the orientation of the ferromagnet and the antiferromagnetic to be perpendicular to each other. Finally the free energy of the antiferromagnetic Mott plateau has to be minimized.
We simplify our calculation by assuming constant parameters ρ, β, and M 0 in Eq. (4), thus neglecting a radial dependence of these parameters due to the trapping potential [33]. We denote the radius of the atom cloud by R c and the radius of the Mott plateau by r M .

A. 2D lattice
In a 2D system we expect the magnetization at the edge to form a vortex-like structure. Furthermore, we claim that for a small imbalance, the vortex will lie in the xy plane with a uniform magnetization component pointing in the z direction. The energetically preferred direction of the antiferromagnetic order in the Mott plateau is perpendicular to that of the ferromagnet at the edge. At the interface, the antiferromagnet in the Mott plateau will have an in-plane magnetization forming a vortex and a z component. The lowest energy corresponds to the maximally allowed z component of the antiferromagnetic order parameter, thus minimizing the in-plane vortex.
We now give quantitative arguments for the physics described in the preceding discussion based on a comparison of the energies of a ferromagnetic edge with different topologies: either a vortex, a domain wall, or a Skyrmion as depicted in Figs. 6 and 10. First we discuss the balanced system. For the vortex, the direction of magnetization is independent of radius, but it rotates by 2π on each circumference. A particular realization of a vortex is M V = M 0 e r . However, for the balanced system, there is global rotation invariance and the plane of the vortex is arbitrary. Using Eq. (4), the energy cost of a vortex is given by E V = πρM 2 0 ln(R c /r M ). Even in the absence of a Mott plateau, the lattice spacing, a 0 , gives a natural cutoff for the core energy leading to E V < πρM 2 0 ln(R c /a 0 ). A vortex naturally fulfills the requirement of vanishing global magnetization in all three spatial directions.
Another possibility is the formation of a domain wall. In the inner ring r M < r < r 0 , there is a uniform polarization (e.g., M = M 0 e z ), and within a finite region, r 0 < r < r 0 + L, the sign of the magnetization is inverted [e.g., M = M 0 (1 − 2(r − r 0 )/L) e z ]. In the outer ring, r 0 + L < r < R c , the magnetization points in opposite direction (e.g., M = −M 0 e z ). While the inner and outer rings have a perfect uniform ferromagnetic order, the domain wall is energetically costly due to the suppression of the amplitude of the order parameter. The energy cost is given by E D = πρM 2 0 (r 0 /L + 1/2)[4 + 4L 2 /(15ξ 2 )], with ξ = ρ/(βM 2 0 ) denoting the coherence length. r 0 and L are not independent of each other but related by the condition of vanishing global magnetization. In absence of any Mott plateau, r M = 0, the smallest allowed value is r 0 /L ≈ 0.6 which increases with r M . Neglecting the term containing the coherence length, we therefore obtain a lower bound for the energy of the domain wall: E D > πρM 2 0 4. Finally, we estimate the energy of a Skyrmion. The magnetization is uniform (e.g., M = M 0 e z ) in the inner ring, r M < r < r 0 , and then it rotates by an angle aπ around a local axis in a ring of width L, r 0 < r < r 0 + L [e.g., M = M 0 cos( r−r 0 L aπ) e z + M 0 sin( r−r 0 L aπ) e r ]. For a = 1, the magnetization in the outer ring is inverted, while for other angles, it has a vortex structure [e.g., M = M 0 cos(aπ) e z + M 0 sin(aπ) e r ]. The Skyrmion interpolates between the inner and outer rings by tilting the order parameter, keeping the amplitude of the magnetization fixed in constrast to the domain wall where the amplitude is suppressed. In the region, r 0 < r < r 0 + L, the magnetization of the Skyrmion changes in the radial direction and along the circumference. The radial dependence of the magnetization gives rise to an energy contribution given by E S = πρM 2 0 (r 0 /L + 1/2)(aπ 2 ). Again the variables r 0 , L, and a are not independent of each other but related by the condition of vanishing global magnetization. By minimizing this energy only, and neglecting the energy cost of the change of magnetization along the circumference, we get a lower bound for the Skyrmion energy: E S > πρM 2 0 [r M /(R c − r M ) + 1/2]π 2 . According to these estimates, the lower bound of the energy for the domain wall is larger than the total energy of the vortex for r M > exp(−4)R c ≈ R c /50, and the lower bound for the Skyrmion is larger than the vortex energy for r M > exp(−5)R c ≈ R c /150. In fact, the real minima for both Skyrmion and domain wall will be larger. Since the radius of the whole atomic cloud is about 50 lattice sites [16,17], our conservative estimate shows that the vortex should be favored for practically any size of the Mott plateau. We note that for the results shown in the main part of this article, r M /R c ≈ 1/2. For the balanced system there is global rotation invariance and the plane of the vortex structure is arbitrary. However, in the presence of a finite imbalance, the energetically preferred magnetization will be a vortex in the xy plane with a uniform ferromagnetic z component. Assuming that the directions of the ferromagnetic edge and the antiferromagnet in the Mott plateau are perpendicular, we expect the direction of antiferromagnet to have a large z component.

B. 3D lattice
Minimizing the free energy in Eq. (4), one finds that the preferred structure of the ferromagnetic edge in a balanced 3D system is a hedgehog [32,33]. We now explain why we expect the antiferromagnetic order in the center to have a planar vortex structure for U > U c and P < P c . First, at the edge of the Mott plateau, the preferred direction of the antiferromagnetic order is perpendicular to the orientation of the ferromagnet order in the edge. A violation of this requirement will cost an energy that scales with the area of the interface r 2 M . In order to fulfill the boundary condition at the edge of the Mott plateau, the antiferromagnetic order in the trap center can neither have perfect Néel order nor a hedgehog configuration, since the latter needs to be oriented in the radial direction. One possibility is to build up the 3D magnetization from the preferred 2D solution for each plane z, which is given by where e z = (0, 0, 1) and e ρ = [cos φ, sin φ, 0] are cylindrical unit vectors. However, this solution is not realized in 3D, since the change in the z component of the magnetization between different planes costs a large energy that scales with the volume of the Mott plateau E AF ∝ r 3 M /a 2 . In fact, the preferred magnetic orders in the Mott plateau are either planar vortex structures like M AF = M 0 e φ or 3D solutions like M AF = M 0 e θ , where e φ = (− sin φ, cos φ, 0) and e θ = (cos θ cos φ, cos θ sin φ, − sin θ ) are spherical unit vectors. For the balanced system, these solutions have the same energy given by E AF ≈ 4πM 2 0 ρr M ln(r M /a). However, imbalance will deform the hedgehog at the edge of the trap leading to a net z component. Such a deformation increases the energy of these solutions except for a vortex in xy plane. Therefore, we expect the magnetization profile in 3D for U > U c and small imbalance to be given by a (slightly deformed) hedgehog ferromagnet at the edge of the trap and an antiferromagnetic order with a vortex structure in the xy plane in the center of the trap. In contrast with the 2D case, where we expect a strong z component in the antiferromagnetic order for U > U c and P < P c , we expect a vanishing z component of the antiferromagnetic order in the Mott plateau in 3D.

APPENDIX B: POLARIZING THE MOTT PLATEAU
In Sec. IV, we propose spin structures for U > U c and P < P c that minimize the total energy while fulfilling the spin conservation (3). The constraints (3) prohibit the formation of a uniform ferromagnet at the edge of the trap if simultaneously the Mott plateau has an antiferromagnetic structure with zero net imbalance. However, since the constraints (3) apply to the whole system, one could imagine a system consisting of a fully polarized ferromagnetic edge and a Mott plateau strongly polarized in the opposite direction, such that the global imbalance is small or even zero. We now justify why such solutions are energetically more costly than the ones proposed in Sec. IV.
We therefore discuss the magnetic structure of a balanced Fermi gas in a 3D trap. The system can be divided into a Mott plateau for radius r < r M and an edge for radius r M < r < R c . In Sec. IV, we claimed that the ferromagnetic structure at the edge forms a hedgehog. Applying the Ginzburg-Landau-type free energy in Eq. (4), the energy of a hedgehog can be estimated. We strongly simplify our calculation by assuming a constant density at the edge [32]. The magnitude M 0 of the ferromagnetic magnetization is therefore constant along with the stiffness, ρ, in Eq. (4). We now estimate the energy of the ferromagnetic hedgehog as E F = 8πρM 2 0 (R c − r M ). The stiffness of a homogeneous Fermi gas is given by ρ = 1/(12k 2 F χ 0 ) =h 2 /(36m n), where k F is the Fermi wave vector, χ 0 the magnetic susceptibility, m the mass of the fermions, and n is the density (see Ref. [33]). For sufficiently small densities, the mass of a particle hopping between nearest neighbors in a 3D cubic lattice is given by m =h 2 /(ta 2 ), where a is the lattice constant and t the hopping matrix element between nearest-neighbor sites. The stiffness on a 3D cubic lattice is therefore given by ρ ≈ ta 2 /(36n), and the energy of the balanced hedgehog becomes E ≈ (2π/9)t(R c − r M )/(a 4 n) tN 1/3 E , where N E is the number of atoms at the edge of the trap. This energy could be gained by uniformly polarizing the edge. However, due to the conservation of the total imbalance, the Mott plateau would then also be polarized by P M = N E /N M , where N M denotes the number of atoms in the Mott plateau. The corresponding cost in energy can be estimated as E AF ≈ P 2 M 2N M 4t 2 /U = 8t 2 N 2 E /(UN M ), where 2N M is the number of nearest neighbors in the Mott plateau and 4t 2 /U is the superexchange. The energy cost of polarizing the Mott plateau is smaller than the energy gain of forming a uniform ferromagnetic edge if U/t > 36N 5/3 E /(πN M ). This is not satisfied for realistic particle numbers N E , N M > 10 3 and N E /N M 1. We thus conclude that in 3D for U > U C and P < P c , the Mott plateau is not significantly polarized. The magnetic structures that minimize the total energy and fulfill Eq. (3) are therefore the ones presented in Sec. IV. We expect similar arguments to hold in 2D.