Phase Transitions of Repulsive Two-Component Fermi Gases in Two Dimensions

We predict the phase separations of two-dimensional Fermi gases with repulsive contact-type interactions between two spin components. Using density-potential functional theory with systematic semiclassical approximations, we address the long-standing problem of itinerant ferromagnetism in realistic settings. We reveal a universal transition from the paramagnetic state at small repulsive interactions towards ferromagnetic density profiles at large interaction strengths, with intricate particle-number dependent phases in between. Building on quantum Monte Carlo results for uniform systems, we benchmark our simulations against Hartree-Fock calculations for a small number of trapped fermions. We thereby demonstrate that our employed corrections to the mean-field interaction energy and especially to the Thomas-Fermi kinetic energy functional are necessary for reliably predicting properties of trapped mesoscopic Fermi gases. The density patterns of the ground state survive at low finite temperatures and confirm the Stoner-type polarization behavior across a universal interaction parameter, albeit with substantial quantitative differences that originate in the trapping potential and the quantum-corrected kinetic energy. We also uncover a zoo of metastable configurations that are energetically comparable to the ground-state density profiles and are thus likely to be observed in experiments. We argue that our density-functional approach can be easily applied to interacting multi-component Fermi gases in general.


I. INTRODUCTION
For almost a century the interacting many-body problem of quantum mechanics has been proven highly demanding both conceptually and practically. Despite decades of intense efforts, Kohn-Sham density functional theory (DFT) [1], the first-principles orbital-based workhorse of computational chemistry and materials science [2,3], remains inapplicable to large systems that are relevant in technological applications and at the forefront of fundamental research. Especially, ab initio descriptions of quantum gases demand a new angle of investigation.
Orbital-free DFT is the only available method for routinely and reliably computing quantum systems that harbor thousands to millions of interacting particles in nonperiodic confinement [4][5][6]. Among the various flavors of orbital-free DFT, density potential functional theory (DPFT) is uniquely qualified for reliably extracting the intricate phases of interacting Fermi gases [7][8][9][10][11][12][13][14][15]. DPFT reduces the many-body problem to two self-consistent equations for the single-particle density and an effective potential that includes the interaction effects. Its capacity in simulating trapped quantum gases, especially for two-dimensional (2D) setups, extends beyond the capabilities of conventional DFT methods, which are either limited to small particle numbers [16,17], periodic confinement [18], or rely on ad-hoc parameterizations of the kinetic energy [19,20], although systematic gradient corrections in 2D are available for electronic systems [21]. DPFT is a scalable approach that enables systematic semiclassical expansions beyond the Thomas-Fermi (TF) approximation across one-, two-, and three-dimensional geometries. The conceptual, theoretical, and numerical work on DPFT over the past years has identified DPFT as an efficient, accurate, and versatile approach for targeting large-scale many-body quantum systems with arbitrary constituents, interactions, and geometries. It has been applied to (i) noninteracting systems for benchmarking purposes [12][13][14]22], (ii) systems in one [23], two [12,13,24], and three [14,22] dimensions, (iii) small [14,22,23] and large [13,[22][23][24] particle numbers, (iv) layered graphene materials [24], (v) atomic physics [7][8][9][10][11]22], (vi) chemistry [22], and (vii) interacting Fermi gases [12]. The overarching feature of all these studies is the systematic methodology of DPFT, whose approximations are universally applicable to a large class of quantum systems. Our DPFT approach naturally accounts for inhomogeneities of large, trapped systems beyond the common local density approximation (LDA) and reliably yields candidates for the ground-state densities.
DPFT thus provides a natural platform to study fermion systems whose components can undergo spatial segregation due to repulsive interactions [25,26]. One seminal example of such a behavior is itinerant ferromagnetism in metals such as iron or nickel [27,28], where valence electrons spontaneously form spin-polarized domains. A quantum-mechanical description of this phenomenon has been proposed by Stoner in his meanfield model, which favors a ferromagnetic state thanks to a short-range screened Coulomb interaction that overcomes the Fermi pressure [29]. That is, same-spin electrons congregate to form regions with nonzero net magnetization at the expense of increased kinetic energy. This simple model has fostered qualitative analyses of many-electron systems, but in other fermionic systems mechanisms beyond the short-range repulsion may suppress phase separation [30,31].
Experimental efforts of preparing a ferromagnetic state in an ultracold atomic system of a balanced mixture of the two lowest hyperfine states of lithium-6 date back to the late 2000s. The initial attempts of observing a para-to ferromagnetic transition proved inconclusive, though some signatures, such as an increase of the kinetic energy, supported its existence [53,55]. The ambiguity came from an alternative explanation of the rapid molecule formation that could produce similar results. To circumvent this problem, the system was prepared in an artificial domain structure, where each of the components initially reside in their respective half of the harmonic trap [54,57]. Such a setup showed stability for a finite time, which was later confirmed by more advanced time-resolved studies of the competition between pairing and ferromagnetic instabilities [58,59].
The theoretical treatment, on the other hand, has been continuously refined in recent history. The analysis of the purely repulsive ground state of a balanced system has been studied with many different approaches, based on second-order perturbation theory [34], Landau's Fermi liquid theory with state-of-the-art Quantum Monte Carlo simulations [36,38,39], lowest-order constraint variational calculation [61], nonperturbative ladder approximation [62], large-N expansion, dimensional -expansion [63], and pseudo-Schrödinger evolution [44,47]. In a three-dimensional (3D) geometry, each of these methods suggested the existence of the ferromagnetic transition, though with a varying critical value of interaction strength, depending on the approach that was used.
The 2D setup considered here potentially offers an escape from the stability problem, as three-body recombination processes are less important in lower dimensions [27]. However, the pairing mechanisms in 2D differ from their 3D counterparts and may still preclude a stable phase-separated state [43]. Moreover, experimental data suggests that in an impurity limit, relaxation to the bound state plays a crucial role [64]. Also theoretical approaches, including mean-field, perturbative and diagrammatic expansions [65,66], polaronic approach [67,68], and quantum Monte Carlo methods [69][70][71][72], have not unambiguously predicted the Stoner transition, not even for purely repulsive mixtures. The subtle interplay between interaction and kinetic energies, an essence of the Stoner ferromagnetism, is greatly affected by quantum correlations and as such, evaluation of beyond-mean-field effects is crucial for such investigations. Also other types of interactions such as Rabi coupling [73] and dipolar forces [74] have been analyzed recently in search of stable ferromagnetic phases in 2D ultracold gases.
In this work, we aim at an unambiguous and quantitatively reliable picture of the phase transitions of a repulsive balanced mixture in realistic settings, that is, for large particle numbers and an inhomogeneous trapping potential. We base our predictions on two different approaches of kinetic energy evaluation, which are discussed in section II. One is based on DPFT, which we introduce in section II A and systematically approximate in section II B; see also appendix A. The other is the multi-component Hartree-Fock (HF) method, see section II C and appendix B. One option for assessing the robustness of a phase transition in realistic settings is to vary the interaction energy functional. section III specifies how we utilize two different interactions-a bare contact, viz., mean-field repulsion and a quantum Monte Carlo energy functional in local density approximation, which we term 'renormalized' contact interaction henceforth. We perform a thorough analysis of emerging partially separated density profiles and their dependence on either of the two energy functionals. We present our main results in section IV. By comparing with DPFT simulations for the mean-field contact interaction in section IV A, we elucidate the inadequacy of the TF model for describing multi-component Fermi gases with contact-type interactions. In section IV B, we benchmark the consequently required quantum corrections beyond the TF approximation against HF results for both interaction energy functionals. This enables us in sections IV C and IV D to (i) reliably predict semiclassical DPFT density profiles of both ground-and metastable states for large particle numbers across interaction strengths and (ii) compare the polarization curve of the resulting phase transition against the Stoner-type QMC prediction for uniform systems. Section V summarizes our findings and points at further, potentially fruitful continuations of our work.

II. KINETIC ENERGY FUNCTIONALS BEYOND THE THOMAS-FERMI APPROXIMATION
A. Multi-component density-potential functional theory A prerequisite for accurate densities and energies is a sufficiently accurate kinetic energy. Ideally, one employs the computationally most efficient TF approximation of the kinetic energy, which is adequate for selected systems. However, the TF density n TF often cannot even qualitatively describe the physics, for example, of a contact-interacting two-component Fermi gas [44]; see also figure 7 in appendix A. In any case, quantum-corrected density formulae have to either validate or replace n TF . While von-Weizsäcker-type gradient corrections to the TF approximation of the kinetic energy density functional are successfully used for three-dimensional geometries [75,76], also in the context of ultracold Fermi gases [44,77], attempts of systematically deriving its 2D analog have produced ambiguous results at best [78][79][80][81][82][83][84]. Density potential functional theory presents an unambiguous solution to this dilemma [12,13]. What sets DPFT apart from other orbital-free DFT approaches are systematic quantum corrections to the TF approximation, which are not available for the commonly used density functional E kin [n] of the kinetic energy in 2D. We are equipped with two such approximation schemes that rely on the explicitly available expression of the Legendre transform of E kin [n], the potential functional E 1 [V ], which is expressed as a single-particle trace that can be systematically approximated using semiclassical techniques. One scheme delivers nonlocal density formulae from a splitoperator approximation of the quantum-mechanical propagator [14,[22][23][24]. The other scheme is based on the Wigner function formalism and Airy-averaging techniques [10,12,13,22,23]. Both schemes are wellestablished with a track-record of excellent accuracy and computational efficiency for a large variety of systems from harmonium to Fermi gases to electron-hole distributions in layered materials [10,[12][13][14][22][23][24]. However, in this work we shall focus on the first scheme, since the semilocal 'Airy-averaged' densities inherit some of the shortcomings of the inadequate TF model when applied to multi-component contact-interacting systems; see appendix A for further details.
In this section, we present the straightforward multi-species extension of the DPFT formalism. In any orbital-free DFT approach, the stationary points of the constrained density functional of the total energy deliver the ground-state densities n(r) = {n 1 (r), . . . , n S (r)} that integrate to the chosen particle numbers N = {N 1 , . . . , N S } for S species of particles (enforced through the chemical potentials µ = {µ 1 , . . . , µ S }), thereby producing the proper trade-off between kinetic (E kin ), external (E ext ), and interaction energy (E int ). In DPFT, we introduce an auxiliary variable, the effective potential energy for species s ∈ {1, . . . , S}, such that the Legendre transform of the kinetic energy functional transforms equation (1) into Strictly, equation (4) holds only for independent particles, but it can be made exact by transferring the interacting part of the kinetic energy into the interaction energy functional E int [n]. For each species s, the V s -and n s -variations at the stationary points of E[V , n, µ] obey and δn s : respectively. The µ s -variation, combined with equation (6), reproduces the particle-number constraint (dr) n s (r) = N s .
Equation (6) states that the particle density is a functional of the effective potential and immediately yields the particle density in the noninteracting case (V s = V ext s ) for any given µ s . Conversely, equation (7) declares V s as a functional of all densities n, such that a self-consistent solution of equations (6)-(8) for any interaction functional E int [n] produces the ground-state density, much like in the Kohn-Sham scheme, but without resorting to orbitals. We initialize the self-consistent loop with n until all densities n have converged with the help of mixing parameters θ = {θ 1 , . . . , θ S }. The chemical potentials µ are adjusted in each iteration to enforce the particle number constraints of equations (8). Further details of the numerical implementation of equation (9) and implications of equations (5)- (8) are well documented in the literature, see references [10][11][12][13]24] and references therein.
The key element of DPFT is the potential functional E 1 [V − µ], which captures the effects of the kinetic energy in place of the density functional E kin [n]. The explicit form of E kin [n] is unknown even for noninteracting systems and its approximations can compete with Kohn-Sham computations only for selected systems. In contrast, E 1 [V − µ] is explicitly available for independent particles [11] in terms of single-particle traces, which can be approximated systematically: for temperatures T ≥ 0, with the function of the single-particle Hamiltonian The single-particle position and momentum operators for D Cartesian dimensions are R and P , respectively.
Here and in the following we omit arguments of functions for brevity wherever the command of clarity permits. Although the (unknown) interacting part ∆E kin of the kinetic energy can formally be transferred into the interaction energy such that equation (10) becomes exact as part of the total energy in equation (5), we neglect ∆E kin altogether-a procedure that comes with an excellent track-record also for DPFT [10,[12][13][14][22][23][24] and reiterates the fact that ∆E kin is often of secondary importance. In the limit of zero temperature we obtain the ground-state version of equation (11); the step function Θ( ) is the zero-temperature limit of Θ T (µ s − H s ) = 1 + e (Hs−µs)/kBT −1 . Any approximation of the single-particle trace in equation (10) yields an according approximation for the particle density in equation (6). We can benchmark semiclassical approximations of E 1 unambiguously if the interaction functional is known exactly or for any noninteracting system, as done in references [12-14, 22, 23].

B. Densities and energies from Suzuki-Trotter-factorized time-evolution operator
Our approximation schemes for DPFT are based on semiclassical expansions of the trace in equation (10), which includes a degeneracy factor g (e.g., spin-multiplicity). We reiterate some of the results in [14] and begin with realizing that equations (6) and (10) which invites tailored Suzuki-Trotter (ST) factorizations of the unitary time-evolution operator U (t) = e − it H . In equation (14), we make use of the Fourier transform of the step function Θ( ), and the integration path from t = −∞ to t = ∞ crosses the imaginary t axis in the lower half-plane. For notational brevity, we drop the subscript s in all species-dependent variables (the here developed formulae hold for all s individually).
We obtain a hierarchy of approximations of n(r) from appropriate coefficients a k and b k of the ansatz where the exponential factors are multiplied from left to right in order of increasing k-values. Equation (15) creates a series of increasingly accurate semiclassical approximations beyond the TF approximation without a gradient expansion. Reference [14] reports particle densities n ν based on up to ν = 7 exponential factors. The (generically) most accurate approximation U 7 in [14] has been proven highly accurate for a variety of systems [14,[85][86][87][88]. However, the computational cost of n ν (r) for ν > 3 in the currently available spatially explicit formulations can be prohibitive at high spatial resolutions and in the case of slow convergence of equation (9). We therefore focus on a three-factor approximation, which has been successfully used in layered 2D materials [24]. Choosing a 1 = 0, a 2 = 1, b 1 = b 2 = 1/2, we obtain with the Bessel function J D ( ) of order D and the effective Fermi wavenumber where [z] + = z Θ(z). The approximate density formula n 3 (r) is the quantum-corrected successor of the TF density with solid angle Ω D in D dimensions, which is obtained from the two-factor splitting a 1 = b 1 = 1 that neglects the noncommutativity of R and P . In all formulae of this work that exhibit the dimensionless constant the quantities of energy are given in units of E and those of length in units of L. For example, µ in equation (18) is implicit for µ/E, and n TF (r) comes in units of L −D , which are not explicitly exhibited in equation (18). For the concrete examples in the sections below we use harmonic oscillator units E = ω and L = /(m ω), which imply U = 1 and V ext (r) = E 2 (r/L) 2 . In contrast to the local TF density, whose computational cost scales with size G of the numerical grid, n 3 (r) is a fully nonlocal expression, which samples the effective potential V in a neighborhood of the position r. This effective averaging over whole regions of position space provides a quantum correction to the TF density and allows for density tails in the classically forbidden region of V , but also lets the computational cost of n 3 scale like G 2 . In this work we make extensive use of the zero-temperature expression n 3 from equation (16), which is efficient enough for isotropic calculations. For anisotropic calculations at T = 0, it is expedient to rephrase n 3 in terms of Fast Fourier Transforms. We derive the according expression n F 3 as follows. Upon approximating the time-evolution operator by U 3 and inverting the Fourier transform of Θ( ) in equation (14), we arrive at where for the Fourier transform of n 3 (r). Defining r 2 = r + r 1 , we express the last integral in equation (21) as such that where for n 3 (r) to distinguish equation (25) from the (numerically identical) expression of equation (16). The computational cost of n F 3 still scales like G 2 , but is reduced by a factor of ∼10-40 since only exponentials (not Bessel functions) have to be evaluated.
The kinetic energy E kin is a functional of n, not of V , but its U 3 -approximated stationary value can be calculated from the ground-state V , see appendix A. Incidentally, equation (26) coincides with the TF kinetic energy E TF kin , whose ground-state values (evaluated with the ground-state effective potential) equal those of the density functional (evaluated with the ground-state density) upon translating 2 between n and V via equation (18). Quantitative differences between E at the stationary point of the total energy are generally expected since n 3 = n TF and V 3 = V TF at the ground state.

C. Hartree-Fock approach to trapped multi-component fermion gases
In this work, we benchmark n 3 -based DPFT densities against Hartree-Fock (HF) results. The derivation of the HF equations for the two-component spin mixture i (r) and ϕ (2) i (r), with i = 1, ..., N/2, are spin-orbitals of the first and the second spin component, respectively. We deploy the interaction terms δEint[n] δn 1/2 (r) as defined below through equations (30) and (35) for the mean-field and the renormalized contact interaction, respectively. The one-particle densities associated with the spin components s sum to the total one-particle density n(r) = n 1 (r) + n 2 (r).

III. INTERACTION ENERGY OF A TWO-COMPONENT REPULSIVE FERMI GAS
Many studies have shown that the simple mean field interaction functional E MF int [n = (n 1 , n 2 )] = (dr) α n 1 (r) n 2 (r) (30) needs to be renormalized in order to reproduce experimental data [43]. Different interaction regimes of homogeneous two-dimensional Fermi gases are commonly defined through the dimensionless gas parameter where k F = √ 2πn is the Fermi wavenumber, and a 2D is the scattering length in 2D. Out of the two common definitions of a 2D that are based on the two-body scattering problem, we use the one in which the energy of the two-body bound state equals b = −4 2 / m a 2 2D e 2γ , where γ ≈ 0.577 is the Euler-Mascheroni constant. The many-body ground state of a two-component ultracold Fermi gas is called the lower (attractive) branch of the energy spectrum. For this lower branch, the crossover from the BEC regime of tightly bound dimers (η > 0 and η 1) to the BCS superfluid regime (η < 0 and η −1) has been thoroughly investigated. The upper (repulsive) branch refers to the excited state of the many-body spectrum that exhibits repulsive behavior. It is usually associated with the para-to ferromagnetic transition of the Fermi gas, above which the polarized mixture is favored energetically over the spin-balanced one. Such a 'Stoner' phase transition is named after Edmund Stoner who devised an early mean-field description of ferromagnetism in the 1930s. Stoner ferromagnetism has been studied extensively, both theoretically and experimentally, mostly in 3D settings. However, the decay into the energetically more favorable lower branch is the main obstacle towards the experimental realization of phase-separated states and invites the precise study of competition between pairing and anticorrelating dynamics. However, this issue can be resolved by introducing an artificial domain structure or through a fast interaction quench.
To date, theoretical studies have been focusing mainly on the homogeneous mixture. The usual approach for describing the homogeneous system is based on the Jastrow-Slater ansatz for the many-body wave function that includes only two-body correlations. With this ansatz, the energy can be minimized, for example, by a Quantum Monte Carlo (QMC) scheme or by introducing constraints as in the lowest-order constrained variational (LOCV) approximation. While QMC is believed to be the most accurate method available, popular alternatives include perturbative methods, either through a diagrammatic expansion, a polaron approach or a large-N expansion.
In this work, we parameterize the interaction energy functional using QMC results. For a homogeneous system, it is expedient to introduce the ratio between the interaction energy and the kinetic energy E 0 kin of the noninteracting system: which is expressed in terms of the dimensional gas parameter η. To compute β, we use the CASINO package with a smooth pseudopotential devised by Whitehead et al. [72] that captures the effective interactions between fermions from different species. In choosing the closed-shell structure of 49 + 49 fermions, we reduce finite-size effects. We approximate the resulting parameterization by a polynomial which is depicted in figure 1 (left), together with its derivative β (η). Note that we restrict ourselves to positive values of η, which corresponds to k F a 2D < 1 and ensures that the mixture stays on the repulsive branch.
Since we address inhomogeneous systems using a local density approximation, we introduce local Fermi wavenumbers k s F = √ 4πn s (here, s ∈ {1, 2}) and local gas parameters η s = −1/ log (k s F a 2D ) as direct extensions of their homogeneous versions. In order to consider locally spin-imbalanced systems, we perform a simple symmetrization to obtain the total energy, consisting of kinetic and interaction terms with local Fermi energies s F = 2 2 πn s /m. In the spirit of the local density approximation, equation (34) invites us to address heterogeneous setups using the approximate density functional of the interaction energy, with and Before studying inhomogeneous systems, let us analyze a homogeneous one with a fixed density n = n 1 + n 2 . Defining relative densities x = n 1 /n and 1 − x = n 2 /n, we rewrite the total energy as The original idea of Stoner is to extract the phase diagram of the homogeneous mixture by minimizing equation (38) at fixed n, resulting in the sample's polarization P = |n 1 − n 2 |/(n 1 + n 2 ) = |2x − 1| as a function of the gas parameter η of equation (31), see figure 1 (right), which identifies the onset of polarization at η ≈ 1.15. A very similar value (η 0 ≈ 1.22) has been obtained recently with fixed-node QMC calculations utilizing hard-and soft-disk potentials [90].
FIG. 1. Left: The ratio β between the interaction energy and the kinetic energy of the noninteracting system and its derivative as a function of the dimensionless parameter η. We evaluate β through a quantum Monte Carlo approach with the help of the CASINO package, using a smooth ultratransferable potential [72] as an inter-particle interaction. Right: The ground-state polarization of the uniform two-component Fermi mixture. The sample becomes partially polarized at a critical interaction strength of η ≈ 1.15.
To calculate a phase diagram similar to figure 1 (right) for a trapped mixture, we cannot rely on a universal parameter η since the local Fermi energy and the polarization change alongside the local density at fixed fermion number. Analogous to studies of mixtures in 3D, we therefore introduce a universal gas parameter for a trapped mixture in 2D, based on the Fermi wavenumber k 0 F = (8N ) 1/4 /a ho of the noninteracting system at the center of the trap and the geometric average a ho = /m √ ω x ω y of harmonic oscillator lengths. We compute the total polarization (viz., net magnetization) P of the trapped mixture as a spatial integral of the local polarization: Our primary objective is the reliable prediction of experimentally relevant density distributions for ultracold two-component Fermi gases with contact-type interactions. While DPFT as our principal tool of investigation can be applied to virtually all such settings and beyond, here we focus on harmonically trapped gases in 2D, subjected to (i) the mean-field contact interaction of equation (30) and (ii) the renormalized contact interaction of equation (35). Its simplicity makes the former a popular approximation, but we find that the (more realistic) renormalized interaction implies quite different density profiles for intermediate interaction strengths. For both choices of the interaction functional and for strong interactions, however, we predict a ferromagnetic state that separates both components into two semi-disks with minimal interface. We also find that by increasing the particle number in the ferromagnetic phase, we diminish the overlap of the two components across the interface, which likely results in decreased dimer formation. In light of our successful benchmarking of DPFT predictions against HF results, we argue that itinerant ferromagnetism on the repulsive branch of the many-body spectrum is a real and robust phenomenon of 2D Fermi gases.
We will begin by demonstrating that the quasi-classical TF approximation, which is supposed to become accurate for large particle numbers, is in fact inapplicable to any particle number for strong interactions. It is therefore necessary to go beyond the TF approximation. Here, we use the systematically quantum-corrected density formula n 3 , which also proves sufficient as it captures the essential features of Hartree-Fock results for up to N = N 1 + N 2 = 110 particles (which is close to the practical limit of high-throughput HF calculations): First, we obtain an excellent quantitative agreement of the density profiles n 3 and n HF , especially for strong interactions. Second, both methods reveal a transition from the paramagnetic phase at small interaction strengths to intricate partially separated profiles to an almost complete segregation into two semi-disks. Third, we predict similar polarization curves (across particle numbers and with both DPFT and HF) that measure this transition according to equation (40). Here, we exhibit a striking disparity to the QMC results for uniform systems, which highlights the crucial influence of the trapping potential.
Our simulations for up to 10000 particles thus present a reliable picture of the phase transitions of the repulsive two-component Fermi gas and inform experimenters about which real-space density profiles to expect across interaction strengths. Of particular importance in this respect are the many qualitatively different configurations of metastable states, which we encounter at intermediate interaction strengths. These metastable density profiles often have energy differences, both relative to each other and to the ground state, of the order of 10 −6 -10 −3 and are therefore likely to be observed in experiments. We detail our findings in the following sections.

A. Lessons from the TF model
Using the TF-approximated kinetic energy E TF kin = c 2 (dr) n 2 1 + n 2 2 , where c = 2π 2 /(mg), we can analytically solve the two resulting variational equations and for the mean-field contact interaction of equation (30). Setting α = c, we circumvent the fine-tuning problem of equation (41), which then yields N 1 = N 2 from µ 1 = µ 2 . That is, a balanced mixture N 1 = N 2 , implying µ 1 = µ 2 and, hence, n 1 = n 2 = n/2 with total energy in 2D does not separate in the TF model-with the proviso that equations (41) and (42) follow from the TF energy functional only if both n 1 and n 2 are strictly positive. Indeed, at the boundary (n i = 0 and n j =i = n) of the support of E TF [n 1 , n 2 ], we find for arbitrary domains D 1 and D 2 = R 2 \D 1 (i.e., E MF int = 0) that yield N i = N/2 = Di (dr) n i (r). Because ofẼ TF < E TF ⇔ α > c = 2π (in harmonic oscillator units), we find a universal (viz., N -independent) phase transition at α = 2π, beyond which any fully polarized phase (characterized by domains D i ) is energetically favored over the paramagnetic phase that features n 1 = n 2 everywhere; figure 7 in appendix A illustrates the analogous situation for the renormalized contact interaction. Therefore, as shown in the following, it is no coincidence that the transition towards the complete split into two semi-disks occurs at α ≈ 2π for all particle numbers and is heralded by slow convergence of the intricate phase patterns (with both DPFT and HF), which we encounter for α 2π. Clearly, the TF approximation cannot adequately describe such a two-component system with repulsive contact interaction, since we could increase the actual kinetic energy not E TF kin of the TF profiles at will by fragmenting the domains D i at fine scales.
Since a nonlocal treatment of the kinetic energy is necessary, we will deploy the semiclassical DPFT framework with the density formulae n 3 of equations (16) and (25). Indeed, n 3 delivers a particle-numberindependent transition from the paramagnetic phase at α 6.2 across α = 2π towards a ferromagnetic phase for α 6.3, see figure 2. The separation into two domains comes with minimal interface, as one would intuitively expect for strong enough repulsion. In between the para-and ferromagnetic states, we observe a zoo of particle-number-dependent phases; see reference [91] for a comprehensive set of plots beyond those shown in figure 2. As the particle number increases, the density profiles at α < 2π become increasingly fragmented into partially spin-polarized domains, with many different domain configurations that are close in energy. Consequently, the convergence of the self-consistent DPFT loop in equation (9) requires a large number of iterations at high spatial resolution and high numerical accuracy. In fact, for N s ≥ 55 at α = 6.23 we find the domains in a state of perpetual transformation: The three according plots in figure 2 are snapshots after ∼ 10 5 iterations.
FIG. 2. The ground-state (n 3 -)approximated local polarizations n1 − n2 for the case of the mean-field contactinteraction functional of equation (30). Up to α 6.2 (< 2π) both density profiles are identical (paramagnetic phase). Strong repulsive interactions, that is, α 6.3 (> 2π), separate both species into two semi-disks (ferromagnetic phase), encircled by a ring of almost identical profiles, akin to the outer region of the cross-sections shown for the renormalized contact interaction in figure 3 below.
These phase transitions of the contact-interacting two-component Fermi gas in 2D are clearly different from (and more diverse than) their 3D counterparts, which merely evolve from a symmetric phase to a splitting into two semi-spheres through an intermediate isotropic separation [44]. Also the transition window (6.2 α 6.3) contrasts with the 3D situation, where the phase transition sets in (and is completed) at smaller α for larger N 1/2 . Furthermore, when constraining the densities to isotropic profiles, we find no separation at all for α < 2π, with energies slightly above those of the anisotropic ground-state separations shown in figure 2.
In summary, DPFT yields a transition into the ferromagnetic phase at α = 2π for both n TF and n 3 . Unlike n TF , however, which does not predict any segregations for α < 2π, n 3 segregates the two fermion components into intricate patterns in a small window below α = 2π. We demonstrate in section IV C below that the (more realistic) renormalized rather than the mean-field contact interaction should be deployed for reliable simulations since the respective density profiles of this intermediate phase differ markedly.

B. Benchmarking DPFT against Hartree-Fock
In order to gain confidence in the DPFT predictions for the mesoscopic particle numbers realized in experiments on contact-interacting ultracold Fermi gases, we now benchmark our n 3 -based density profiles against HF calculations for particle numbers up to N 1 + N 2 = 110. Figure 3 demonstrates that n 3 captures all essential features of the HF densities, both qualitatively and quantitatively. This includes overall structure, cloud sizes, evanescent tails, as well as extent and gradient of the interface between the two fermionic species. The latter property is of particular importance for reliable estimates of dimer formation rates [57]. Although n 3 even predicts nontrivial HF features like the marginally asymmetric separation depicted in the top row of figure 3, the similar feature of n HF in the middle row is not captured. There, however, n 3 accurately exhibits the partial polarization magnitudes of n HF as well as the revival of the minority component in the outer cloud ring. That is, n 3 presents a reliable density expression at least for strong interactions that clearly yield a ferromagnetic state in both the HF and DPFT framework. This close match between both methods is less pronounced for intermediate interaction strengths, where the complex structures in the outer cloud regions are somewhat disparate and only HF predicts an asymmetric separation (bottom row of figure 3). As we will discuss in the next section, the regime of intermediate interaction strengths is the richest in terms of possible metastable configurations that are energetically comparable and poses, therefore, numerically challenging for both HF and DPFT.

C. Phase transitions and total polarizations of Fermi gases with renormalized contact interaction
In section (IV B), we established the high quality of the quantum-corrected DPFT formula n 3 , see equations (16) and (25), by benchmarking against HF results for mean-field and renormalized contact interactions for moderate particle numbers. We now present the main results of this work for balanced mixtures with N = N 1 + N 2 up to 10000 as required for describing realistic experimental setups [27]. Figure 4 shows two universal features of the local polarizations n 1 − n 2 in n 3 -approximation across various balanced mixtures (N 1 = N 2 ) and interaction strengths (viz., values of the universal gas parameter η 0 ). First, η 0 2 implies a paramagnetic phase, where both fermion components have the same density profile. Second, when exceeding a critical interaction strength η c 0 , the two components transit into a bipartite splitting towards an almost complete separation into two semi-disks, whose interface becomes more pronounced with increasing particle number. Both features are in line with our results on the mean-field-interacting gas, see figure 2 and reference [91]. The depicted sequences (columns in figure 4) of plots for fixed N 1/2 are snapshots of the parato ferromagnetic phase transition.
The noise imposed on the densities at the start of the self-consistent loop of equation (9) breaks the spherical symmetry and lets the two fermion clouds equilibrate with a random orientation that differs from run to run, if η 0 is large enough to induce anisotropic ground-state separations. The smooth density profiles presented here cannot be obtained using the TF approximation (see also figure 7) and are distinctly different from the profiles for the mean-field contact interaction, see figure 2. Figure 4 shows a rich zoo of particle-numberdependent phases at intermediate interaction strengths, from isotropic as well as symmetry-broken partial separations to intricate domain wall structures for larger particle numbers. That is, realistic descriptions of contact-interacting two-component Fermi gases in 2D require interaction density functionals that supersede the mean-field kernel α n 1 n 2 of equation (30).
As we increase the particle number beyond ∼ 100, we find more and more qualitatively disparate phases with minute energy differences and increasingly intricate separation patterns in the window 2 η 0 3. For N 1000, the unambiguous identification of the ground state densities is thus both more cumbersome and less relevant since the ground state will less likely be encountered in experiments. For example, the two metastable configurations for N 1/2 = 500 at η 0 = 2.19 have energies of E = 28733.2 and E = 28733.4, close to the E = 28727.7 for the ground state configuration, which is an isotropic separation that emerges smoothly from the configuration at η 0 = 2.06. These two pictures are snapshots of two separate runs of the selfconsistent DPFT loop taken after 30000 and 100000 iterations, respectively. Evidently, the intermediate interaction regime easily promotes transformations between energetically comparable, though qualitatively disparate, density profiles during equilibration. We have made similar observations with the HF calculations for smaller particle numbers and expect that such transient states can also be seen in the laboratory. When employing limited numerical accuracy, one cannot determine whether configurations that are energetically above the ground-state energy are (i) transient or (ii) present a metastable state of a local minimum, even if the self-consistent DPFT loop has converged. In section IV D, we will illuminate how to identify metastability within our DPFT approach and showcase further examples that complement the global picture of the phase transition displayed in figure 4.
Finally, we reach the regime of mesoscopic particle numbers (N 1/2 = 5000, shown in the last column of figure 4) that are commonly studied in ultracold-gas experiments with (quasi-)2D geometries. The parato ferromagnetic transition parallels what we find for smaller particle numbers, but we expose important aspects that are difficult to extrapolate from the results on small N 100. Most importantly, the sharp and seemingly random interfaces that emerge during the self-consistent DPFT equilibration are reminiscent of the 'fragmented' TF regime of degenerate domain structures. Indeed, even a substantial restructuring of the domains has only marginal impact on the energy. It is thus not surprising that the convergence towards the ground state profiles demands considerable computational effort. We thus present, as for the case of N 1/2 = 500, metastable configurations whose energies are slightly above those of the isotropic ground state for η 0 η c 0 . Even in the ferromagnetic regime, the DPFT equilibration can get stuck despite various measures that aid the convergence (see section IV D for details). Therefore, and for the sole purpose of obtaining a clean ferromagnetic separation, we converge with an inter-specific Coulomb interaction of strength γ superimposed onto the renormalized contact interaction, followed by an equilibration with gradual reduction of γ, see appendix A for details. The result is a symmetrically split configuration with slightly lower energy than configurations whose domain walls are marginally shifted relative to the ground-state configuration. As an example, we present such a situation in figure 4 for η 0 = 5.77. We verified for small particle numbers that the intricate separation patterns at intermediate interaction strengths are recovered when this artificial interaction is gradually switched off. We recognize in retrospect that the TF model proves qualitatively correct regarding key properties for large particle numbers: The TF densities of mesoscopic mixtures subjected to the mean-field contact interaction share three characteristic features with the quantum-corrected densities for both the mean-field and the renormalized contact interaction. First, the para-to ferromagnetic transition becomes relatively sharp as the particle number increases, both in terms of the 2D scattering length and in terms of η 0 . Second, the ferromagnetic TF solutions, which are degenerate for any domain fragmentation, are structurally similar to the metastable n 3 configurations that we encounter for N s = 5000- figure 4 shows two such examples at η 0 = 2.52 and η 0 = 5.77, respectively. Third, the density gradients at the domain interfaces of the quantum-corrected mixtures increase with particle number. Judging from our simulations for up to 10000 particles, we expect this trend to continue toward the quasi-classical limit embodied by the TF model, which features vertical interfaces.
The universal character of the phase transitions for the harmonically confined fermion mixture illustrated in figure 4 is summarized with figure 5, which shows the polarization P of equation (40) as a function of the universal gas parameter η 0 of equation (39). Both the onset of phase separation and a dampened increase towards full polarization are universal across particle numbers and follow qualitatively the homogeneous case shown in figure 1. The crossing of metastable-and ground-state polarizations shown in the inset of figure 5 parallels the crossing of the respective energies near η c 0 (see also figure 6), where the mixtures transit into the ferromagnetic phase. For example, for N 1/2 = 15 around η 0 = 3, the polarization curves of the metastable states, which exhibit anisotropic (isotropic) separations for η 0 3 (η 0 3), cross the polarization curves of the ground states, which exhibit isotropic (anisotropic) separations for η 0 3 (η 0 3).
FIG. 5. The polarization (P-)curves for different particle numbers collapse onto a unique function of the universal gas parameter η0. In particular, we find a universal onset of phase separation across particle numbers. Accepting the disparities between the P-curves of different Ns as error estimates of our DPFT approach, we find the polarizations of the metastable states (marked by squares with color code of the respective P-curves; displayed in figure 6) to approximately equal those of the ground states, but they are located near kinks of the P-curves that indicate the transition into the ferromagnetic state at η c 0 . The inset magnifies the according region of the phase diagram.

D. Metastable density profiles
For the contact-type interactions studied in this work, the self-consistent DPFT loop of equation (9) requires ∼ 10 2 -10 6 iterations with admixing parameters θ s ∼ 0.2-0.6 for the densities to converge, in the sense of the correlation measure of the vectorized densities n = {n(r 1 ), . . . , n(r G )} on the numerical grid of size G to fall below a predefined threshold. We declare convergence once χ < 10 −12 . The required number of iterations tends to be higher in the vicinity of η c 0 , where the density profiles transition between qualitatively different phases within a narrow window of η 0 . Depending on the interaction strength, we have encountered qualitatively different density profiles with relative energy differences of ∼ 10 −6 -10 −3 . Metastable profiles like those shown in figure 6 (panels 1a-5a) are therefore likely to be encountered in experiments alongside the actual ground-state profiles (panels 1b-5b).
The self-consistent loop of equation (9) can converge into different local minima. The following recipe increases the probability of converging to the actual ground-state densities-in our case, the global minimizers of the (3 -approximated) total energy functional-rather than to a metastable state. We always superimpose white noise 3 with relative magnitude of 10% on the initial densities n (0) to ensure that anisotropic density profiles are probed during the early stages of the self-consistent loop. Together with the sampling of the whole range of possible interaction strengths, this allows us to identify phase transitions between isotropic and anisotropic ground-state densities. In order to sample a sufficiently large space of densities, we impose noise on n (i) throughout the self-consistent loop akin to simulated annealing, with the noise magnitude incrementally decaying towards machine precision at the end of the loop. We converge to nearly machine precision in order to phase out intermediate metastable states on the way to the global energy minimum. But even with these measures in place, equation (9) occasionally gets trapped in local energy minima, such that repeated runs are called for, especially near transitions between qualitatively different phases. Although the sequence of ground-state separation patterns for N 1/2 = 15 in figure 4 suggests a sharp transition (around η 0 = 3) from isotropic separations towards the complete (symmetric) split into two semi-disks, panel 3b of figure 6 emerges smoothly from the metastable domain pattern encountered at lower η 0 (panel 2a). We make an analogous observation for N 1/2 = 55, where the pattern with η 0 = 2.26 in figure 4 emerges smoothly from the metastable branch, exemplified by the panels 4a and 5a. The contour density plots identify candidates for metastable profiles (panels 1a-5a; energies denoted by E * ) and ground-state densities (panels 1b-5b) with energies E, marginally below E * .
We note that labeling density profiles as metastable is a challenging enterprise if based on (an approximate) energy evaluation alone. The minute energy differences in concert with numerical uncertainty are one reason: Differing density patterns and associated energies may emerge from different accuracy criteria for the numerics that have the potential to confuse the energy-based identification of metastability. Therefore, a high accuracy of internal computations and a high spatial resolution is required. But, more importantly, an unambiguous statement on metastability based on energy alone is possible only if the employed method (e.g., the energy functional) is exact and if the errors of the numerical procedures that target the ground state are negligible. Then, a density profile with higher energy is not the ground state density. Rigorously identifying the groundstate density is off the table in all other settings-in particular, if the energy functional is approximate. That said, making sure that numerical errors are insignificant and taking the approximate energy functional as a given, we can and do use energy as a criterion that decides on metastability if the corresponding DPFT calculations are initialized with identical input (except for the initial noise, which phases out after many loop iterations).
The multi-particle ground state for a rotationally symmetric trap is isotropic. We thus argue that the symmetry-broken density profiles, as obtained from our approximate HF and DFT schemes, represent some of the information encoded in the correlation functions (for example, density-density-correlations), which can be anisotropic also for a rotationally invariant many-body Hamiltonian with interactions [92]. It is then conceivable that the actual isotropic ground-state can be constructed as an appropriate superposition of anisotropic states, each of which giving rise to 'single-shot' density profiles akin to the approximate HF and DFT outcomes. However, such a superposition presents a fine-tuning problem that is irrelevant for both experiments and simulations-unless the isotropic ground-state density is robust against anisotropic perturbations that inevitably emerge both in the laboratory and in simulations. We have the latter situation, for example, if η 0 is well below η c 0 . In our DPFT scheme, different types and magnitudes of perturbations on nonequilibrated densities during the self-consistent loop of equation (9) determine the minimizers (equilibrated density patterns) at the (local or global) energy minima. For instance, initializing equation (9) with an isotropic density and omitting noise, we observe convergence to an isotropic local minimizer and miss the anisotropic global minimizer because the local minimum is robust enough against numerical rounding-off errors.

V. CONCLUSIONS
We mapped the para-to ferromagnetic phase transition of repulsive two-component Fermi gases in two dimensions beyond the local density approximation. By recovering the essential features of Hartree-Fock (HF) densities from density-potential functional theory (DPFT), supplied with systematically quantum-corrected semiclassical expressions for the particle densities that become more accurate for larger particle numbers, we gained quantitatively reliable density profiles that experimenters can expect to observe in realistic settings. Mapping both ground-state and metastable density configurations across interaction strengths for up to 10000 fermions, we predicted that strong contact-type interactions segregate the two fermion species into two semidisks. We also found that the overlap of both species in this ferromagnetic phase can be reduced by increasing the particle number, which will likely suppress dimer formation.
We revealed several universal features of this phase transition across system sizes. All curves obtained from integrating the local polarizations |n 1 − n 2 | essentially collapse onto a single graph and summarize the deviation from the paramagnetic state as a function of the contact-interaction strength. This confirms the long-standing prediction of a Stoner-type polarization behavior across particle numbers in terms of a universal gas parameter η 0 , but we found the (nonuniform) harmonic trap responsible for stark quantitative differences to Quantum-Monte Carlo results for a uniform setting. Apart from the successful benchmarking against HF results, we gained confidence in the reliability of the DPFT results by observing that the density profiles transit smoothly from the para-to the ferromagnetic phase when increasing η 0 incrementally. In a sense, this observation also holds as η 0 crosses a critical value 2.0 η c 0 3.1. Then, the anisotropic phase separations, which are metastable for η 0 η c 0 , smoothly transform into the anisotropic ground-state phase, which shows isotropic separations for η 0 η c 0 . With increasing particle number, we observe more and more qualitatively different intricate patterns of partially separated profiles, which are robust against small finite temperature, see appendix A. They typically have relative energy differences of the order of 10 −4 or less, which makes them likely to be observed in experiment in lieu of the actual ground state. In this regime of intermediate interaction strength, our markedly differing results from employing (i) the mean-field and (ii) the renormalized contact interaction, point in favor of the latter.
We also demonstrated the need to improve upon the Thomas-Fermi (TF) approximation for the kinetic energy, which fails to provide even a qualitatively reliable picture of the phase transition. We suspect the TF model's poor performance to originate in the inter-species contact-type interactions considered in this work. It is conceivable, however, that large systems with nonlocal interactions may be adequately addressed with the TF model. In any case, systematic corrections to the TF approximations have to either validate or replace the TF density expression. The semiclassical DPFT framework presented here is uniquely qualified to execute this task. Our DPFT code is part of the C++ software package 'mpDPFT', available at https://doi.org/10.5281/zenodo.4774448. Imbalanced mixtures can be addressed without any modifications of the code, and implementations of additional density-dependent interaction functionals are straightforward. DPFT is a multi-purpose tool for addressing quantum-many-body problems in one-, two-, and threedimensional geometries, such as multi-component dipolar Fermi gases. They are recently coming into focus with first experimental realizations [93][94][95] and are candidates for showing even richer phase transitions than the ones predicted here for contact-type interactions.
Benchmarking of DPFT densities against exact results. In figure 7 we benchmark the semiclassical density formulae n 3 and n T Ai against exact densities of a harmonically trapped Fermi gas in 2D. For n T Ai , we use temperatures T i such that k B T 1 = 10 −6 ω, k B T 2 = 10 −6 E ex (with the exact energy E ex = 770 ω of the two-component system of N = N 1 + N 2 = 55 + 55 particles), and k B T 3 = ω. Both n T1 Ai and n T2 Ai are close to n ex and showcase the capacity of semiclassical DPFT for describing the region around the quantumclassical boundary. We also see that T 2 is small enough for extracting ground-state properties from the finite-temperature formula n T Ai , see equation (A1) in appendix A. However, as argued in the following, n T Ai is not suitable to address multi-component systems with contact-type interactions due to the bulk properties of n T Ai inherited from the TF density-despite n T Ai being generically more accurate than n 3 . Ai ≈ E T 2 Ai ≈ 770.61 ω, E T 3 Ai ≈ 770.36 ω) are close to the exact quantities, as expected for over one hundred particles. The Friedel oscillations in the bulk of nex are not reproduced by our semiclassical densities, which rather give an approximate average account of the exact oscillations. Center: The densities n T 1/2 Ai approximate nex well over ∼ 20 orders of magnitude, while n 3 exhibits unphysical Bessel-functioninduced oscillations far into the classically forbidden region. This imperfection of n 3 is of minor concern for the present work since the phase separations between the two fermion species develop in the bulk. Right: For the two-component mixture at renormalized contact interaction strength η0 = 12.99, the TF approximation of n1 − n2 simply amplifies the noise imposed at the start of the self-consistent loop of equation (9) because different positions are decoupled for contact-type interactions (see section IV A, and [44] for the 3D case of the mean-field contact interaction). This shortcoming of TF is inherited by n T Ai and contrasts with the performance of the nonlocal n 3 approximation, whose converged smooth density profiles (for η0 = 12.99) roughly envelop the 'TF noise' depicted here, cf. figure 4.
The leading term in equation (A1) with the natural logarithm log( ) recovers the (finite-temperature) TF density for uniform effective potentials (i.e., ∇V = 0 everywhere). Equation (A1) is exact up to the leading gradient correction O(∇ 2 ) , and thus presents a systematic correction to the TF approximation O(∇ 0 ) . The 'Airy-average' in equation (A1) also contains higher-order gradient corrections that provide an accurate density tail across the boundary of classically allowed and forbidden regions of V , where the TF approximation can fail epically, even if supplemented with the leading gradient correction [13]. However, in this work we focus on the bulk, not on the evanescent tails. The two fermion components separate in the bulk, where n T TF is the dominant component of (and transfers its inadequacy to) n T Ai . Indeed, we find that the semilocal nature of n T Ai , which stems from the derivatives of V , does not prevent the convergence into random domains of partial polarization, very similar to figure 7 (right). In contrast, n 3 is fully nonlocal and retains less of the TF characteristics in the bulk. Both these features promote smooth density profiles when using n 3 even for large particle numbers.
Numerical details on n 3 . The approximate nature of n 3 can lead to locally negative densities (exhibited by the 'gaps' in the evanescent region of the graph in the central panel of figure 7), although this effect can be regarded negligible for 10 particles. Since the renormalized contact interaction requires strictly positive densities everywhere, cf. equation (31), we replace n 3 by [n 3 ] + + 10 −16 when evaluating equations (36) and (37). We determine the chemical potentials µ s in each iteration of the DPFT loop via an adaptive bisection algorithm, enforcing a relative accuracy of at least 10 −6 for the particle numbers N s . n 3 at finite temperature. The derivation of the finite-temperature version of n 3 (r) will be given elsewhere [22]. Here, Γ ( ) denotes the Gamma function, with τ = µ/(k B T ), and is easily tabulated for all required values of κ = ( k) 2 /(8m k B T ) = U k 2 /(k B T ), where k is the magnitude of the wave vector k. The computational cost of n T 3 (r) scales like G log G. This contrasts with equation (16), where the density n 3 (r) at each of the G grid points, indexed by r, requires a summation over the whole grid. Naturally, there is a trade-off between grid size and accurate enough evaluation of the y-integral-as a rule of thumb, n T 3 outperforms n 3 for G 50 3 ≈ 350 2 and n F 3 for G 100 3 = 1000 2 . That is, the computational efficiency of n T 3 exceeds that of n F 3 for most 3D applications. In 2D settings, however, it is expedient to use n F 3 instead of n T 3 , unless very high spatial resolution is required or if finite-temperature observables are targeted.
Our results in figure 8 demonstrate that n T 3 transforms smoothly into the ground-state density n 3 as T tends to zero. We thus expect that the intricate phase separations at intermediate interaction strengths can also be observed in experiments performed at small but finite temperatures. Regularization with a transient Coulomb interaction. For N 1/2 = 5000 in the ferromagnetic phase, we arrive at the ground-state configuration that splits the Fermi components symmetrically by initializing the self-consistent DPFT loop with a symmetrically split configuration. We aim at a smooth DPFT equilibration into the actual ground-state profiles by counteracting the tendency of large particle numbers at strong contacttype interactions to simply enhance initially imposed noise. While this could in principle be achieved with the finite-temperature formula n T 3 , which washes out any noise and oscillatory structures if T is high enough, see figure 8, n T 3 incurs high computational cost as T comes close to zero. We therefore regularize by adding an artificial mutually repulsive long-range (viz., nonlocal) interaction to the effective potential V s , specifically, the Coulomb interaction potential in 2D. We express V Coul s in terms of Fourier transforms for an efficient numerical implementation and regularize the Coulomb singularity in reciprocal space that is introduced through the numerical discretization. In the case of N 1/2 = 5000, we employ γ = 1/4 in harmonic oscillator units.
Kinetic energy in U 3 -approximation. We obtain approximations of the (ground-state) kinetic energy (dr) ∇ 2 r n (1) (r; r ) r =r (A6) by deriving approximate one-body reduced density matrices n (1) (r; r ) in terms of the effective potential V .
With H = T + V , where T = P 2 /(2m) is the single-particle kinetic-energy operator, equation (3)  In the spirit of equation (14), we thus identify the approximate one-body reduced density matrix n e it µ r | U 3 (t) |r = g (dr ) which is consistent with the Suzuki-Trotter approximation inherent to n 3 . Equation (A8) follows the structure of equation (16) with b = r 2 + (r + r − r ) 2 and yields equation (26): which can also be used in lieu of the finite-temperature kinetic energy E In the single-particle approach each fermion is described by a spinor φ i (j) built from spin-orbitals ϕ i (r j , s j ), where r and s are spatial and spin coordinates, respectively. The indices i and j take values from 1 to N , where N is the total number of fermions. The spin-orbitals obey S s=1 (dr) ϕ * i (r, s) ϕ j (r, s) = δ ij , where S is the number of spin components. For noninteracting fermions the solution of the multi-particle Schrödinger equation is where P is a permutation of N elements and each φ i (j) is the solution of the single-particle Schrödinger equation. The lowest-lying set of solutions has to be taken for the ground state. The total one-particle density is given by n(r) = S s=1 n s (r), where n s (r) = N i=1 |ϕ i (r, s)| 2 is the one-particle density of the spin component s.
When the interaction is turned on, Equation (B2) can be taken as a variational ansatz. The ground state of a system is found by minimizing the total energy functional with respect to {ϕ * i } and {ϕ j }. This functional is given by is the Hamiltonian of our system. Equation (B3) can be rewritten as where h i is the average one-particle energy, K ij is the average interaction energy between the states ϕ i and ϕ j , and J ij is the interchange energy: (dr)(dr ) ϕ * i (r, s) ϕ * j (r , s ) V ss int (r − r ) ϕ i (r , s ) ϕ j (r, s).