Disordered structures in ultracold spin-imbalanced Fermi gas

We investigate properties of spin-imbalanced ultracold Fermi gas in a large range of spin polarizations at low temperatures. We present results of microscopic calculations based on mean-field and density functional theory approaches, with no symmetry constraints. At low polarization values we predict the structure of the system as consisting of several spin-polarized droplets. As the polarization increases, the system self-organizes into a disordered structures similar to liquid crystals, and energetically they can compete with ordered structures such as grid-like domain walls. At higher polarizations the system starts to develop regularities that, in principle, can be called supersolid, where periodic density modulation and pairing correlations coexist. The robustness of the results has been checked with respect to temperature effects, dimensionality, and the presence of a trapping potential. Dynamical stability has also been investigated.


I. INTRODUCTION
The crossover between Bardeen-Cooper-Schrieffer (BCS) type superfluid and Bose-Einstein condensate (BEC) of diatomic pairs has been, in the last years, a subject of considerable theoretical and experimental efforts [1][2][3][4]. The research is driven by developments in the field of ultracold Fermi gas. In the case of spin-symmetric gas, the phase diagram as a function of temperature and interaction strength is well-established, apart from the presence of a pseudogap regime which is still the subject of ongoing discussion [5,6]. In contrast, spin-imbalanced systems lack an undisputed phase diagram [7][8][9]. There are various proposed phases. For high values of the spin polarization, but still below the critical value above which the system transit into a normal state, a homogeneous gapless superfluidity may emerge [10][11][12]. An alternative scenario was proposed by Larkin and Ovchinnikov (LO) [13] and simultaneously by Fulde and Ferrel (FF) [14]. It is described by the spatial modulation of the pairing field (order parameter) ∆ due to the mismatch between the Fermi surfaces of the paired species. The FF-type of modulation is described as ∆(r) ∼ |∆|e iq·r where q is the momentum of the Cooper pair due to the mismatch. It breaks the time-reversal symmetry while preserving the translational symmetry. In contrast, LO-type of modulation breaks the translational symmetry while preserving the time-reversal symmetry and has the form of ∆(r) ∼ |∆| cos(q · r). This so-called LOFF phase remains one of the long-sought phenomena in ultracold Fermi gases. To date, the predictions (both theoretically and experimentally) point out that the LOFF phase occupies a small portion of the phase diagram [7][8][9]15] and it is still an elusive phenomenon. The presence of the LOFF phase and its stability has been a subject of intensive research. In particular, the interplay between LOFF and phase separation as a main competing scenario has been discussed [16,17].
In the case of the LOFF phase, it is implicitly assumed that the long-range order is set in the system producing highly symmetric configurations. Their appearance has consequences for the spectrum of low-energy excitation modes. This issue has been raised and discussed in Refs. [18,19] where low-energy effective Hamiltonian has been proposed FIG. 1. Spatial distribution of absolute value of the order parameter |∆(r)| for several values of global spin polarization P. They correspond to self-consistent solutions of energies close to the ground state energies. The size of the simulation box is 192k −1 F × 192k −1 F , with Fermi momentum kF ≈ 1. The coupling constant g was adjusted to mimic strength of the pairing correlation as for the unitary Fermi gas. At nodal lines/regions ∆ ≈ 0 (blue color) the order parameter changes its sign. assuming a certain type of ordering or the deviations from regular structures. In most of these studies, Ginzburg-Landau model is employed based on its relation to microscopic approach [18,[20][21][22]. However, little is known about the possible structures that lacks regularities, in particular whether they can compete with the ordered structures. This brings us to the question concerning energetics and stability of configuration with broken symmetries. Although the issue of stability has been raised by various authors, it is usually considered in terms of expansions around the symmetric configuration [18,19,23], whereas highly nonordered structures have never been studied. There are two aspects that require to shed light upon. The first one concerns the energy differences between ordered and highly nonordered structures and their susceptibility with respect to thermal fluctuations. The second one concerns the typical time scales and the character of the evolution of these structures. Insight into both aspects can be provided by microscopic approach which explicitly takes into account fermionic degrees of freedom to account for possible shell effects coming from single-quasiparticle states. At the technical level, applying the framework is very challenging as one needs to execute calculations with no symmetry constraints. Nevertheless, such studies are particularly timely in the light of experimental efforts to detect LOFF-like structures [24][25][26][27][28][29][30][31]. They also provide information about the characteristics of low-energy modes that may appear in such systems.
In this paper, we focus on the solutions of a spin-imbalanced ultracold Fermi gas as a function of a wide range of global polarization P = (N ↑ − N ↓ ) / (N ↑ + N ↓ ), where N ↑ , and N ↓ are the numbers of spin-up and spin-down particles, respectively. In order to explore the richness of possible configurations, we break the translational symmetry of the pairing field ∆(r) by perturbing it with a spatial, complex random noise and subsequently perform unconstrained energy minimization. This symmetry breaking procedure can be thought of as mimicking thermal fluctuations which occur during the cooling process.
The sample of obtained solution is presented in Fig. 1. For low polarization values (e.g., P = 5%), the randomly perturbed system produces several distinct spin-polarized droplets called ferrons [32][33][34][35]. The increase in the spinimbalance induces elongation of nodal lines, which eventually form more complex disordered structures extended over the whole system. Here, for medium polarizations such as P = 20 − 25%, the nodal structure of the order parameter resembles rather a liquid crystal. There is a vast number of different nodal structures which are practically energetically degenerate. As the total length of nodal lines becomes comparable with the system size, the number of different ways to arrange the nodal lines decreases; the system starts to develop an order in configuration. For high polarization values (i.e., P 30%), the system arranges itself in different regions where each region shows a unique lattice-like character. For extremely high polarization (P 40%), one may expect to find a single lattice structure that resembles a supersolid. However, as the polarization increases, the mean distance between the nodal lines becomes comparable to the coherence length ξ. This degrades the superfluid phase, and eventually, the system transits to the normal state. In the following, we argue that the structures, as shown in Fig. 1, are more likely to be observed experimentally compared to typically considered scenarios like LOFF state.

II. THEORETICAL FRAMEWORK
We consider Fermi gas on the spatial lattice with periodic boundary conditions at finite temperatures. The stationary configuration has been found through minimization of the grand canonical potential Ω = E − µ ↑ N ↑ − µ ↓ N ↓ − T S, where N ↑↓ = N ↑↓ , and ... denotes the thermal average. Entropy and temperature are denoted by S and T , respectively. The energy density functional has been chosen as (we set units m = = k B = 1): τ is defined as the kinetic energy density and ν is the anomalous density which can be associated with condensate wave function of Cooper pairs. These densities read (σ = {↑, ↓}): For completeness we define normal density which is main experimental observable. The Fermi-Dirac function, f (E) = 1 + e E/T −1 is introduced to study finite temperature effects. The single quasiparticle energies E n and wave functions ϕ n (r) = [u n↑ (r), u n↓ (r), v n↑ (r), v n↓ (r)] T are determined from the eigenequation Hϕ n (r) = E n ϕ n (r), where H is the Bogoliubov-de-Gennes (BdG) Hamiltonian: In this case the mean-field hamiltonian is spin independent and reads h ↑ (r) = h ↓ (r) = − 1 2 ∇ 2 and µ σ is the chemical potential for a given spin component (σ =↑, ↓). The order parameter is related to the anomalous density ∆(r) = −gν(r), and g is the s-wave coupling constant. In order to check stability of the predictions with respect to the choice of the energy density functional, we employ also more advanced form which is particularly suited for spin asymmetric systems at the unitary regime [36][37][38][39]. This is so-called Asymmetric Superfluid Local Density Approximation (ASLDA) and reads As compared to the pure BdG approach, the above functional comprises of polarization-dependent effective mass α σ , the mean-field term D, which scales with a density as a Fermi gas and is polarization dependent too. The paring coupling function g also acquires density dependence. The functions α σ , D and g were chosen in such a way to assure compatibility with quantum Monte Carlo calculations for both spin-symmetric and spin-imbalanced unitary Fermi gas, for details see [37].
Most of the results presented here were obtained within BdG functional for 2D geometries. It was motivated from one side by geometries of traps which can be nowadays realized experimentally [40,41], but also by numerical complexity, which was reduced in this case to reasonable limits. In order to investigate the robustness of the results, we have also performed selected 3D studies by applying full ASLDA functional. The coupling constant g can be related to the scattering length. However, explicit relation is not needed in the context of these studies. Instead, we fix the value of g in such a way to obtain ∆/ε F ≈ 0.5 when considering the spin-symmetric system (P = 0), which corresponds to the so-called unitary (strongly interacting) regime [42]. The Fermi energy ε F = k 2 F /2 is related to the (2D) Fermi momentum k F = 2π(n ↑ + n ↓ ). The BdG equations were solved numerically on a spatial lattice of size N x × N y with lattice spacing dx = 1 (definition of the length unit). The particle numbers N ↑ + N ↓ were adjusted to satisfy k F ≈ 1. The minimization procedure starts with the uniform solution, which has been perturbed initially by the random external potential ∆ ext (r) = ∆ r (r) · exp{iφ r (r)}. The amplitude ∆ r (r) has a positive real value constrained not to exceed 1% of the pairing strength in symmetric system. Both, ∆ r and phase φ r were randomly selected at each lattice point. The perturbation was switched off after 50 self-consistent iterations, which amounts to less than 5% of the total simulation effort. Eventually, a converged solution was obtained in an unconstrained manner. The calculations were executed by means of the W-SLDA Toolkit [38,43,44]. In appendix E we provide reproducibility packs, allowing for numerical reproductions of the presented below data.

III. SELF-CONSISTENT SOLUTIONS AND THEIR ENERGIES
We have examined self-consistent solutions obtained according to two methodologies: using a random starting point (see above) and performing entirely unconstrained minimization or imposing a certain constraint on the order parameter which allows us to obtain a self-consistent solution with required symmetry. Both approaches have limitations, and none guarantee that the lowest energy state is obtained. The unconstrained minimization process can get locked in a local minimum, while the constrained approach scans only a selected subspace of possible solutions. Examples of obtained solutions are presented in Fig. 2 for two selected spin imbalances P = 5% (top) and 20% (bottom) and with imposed constraints (a-c) or for the unconstrained process (d-e). The details on imprinting the nodal lines are presented in Appendix A.
In the regime of low spin imbalances (e.g., P = 5%), the randomly seeded minimization process converges to states consisting of many ferrons randomly distributed in space, see Fig. 2 (d-e, top). Within this class of solutions, we find that states with a smaller number of ferrons (but larger in size) have lower energies. This can be understood as an additional energy cost related to the curvature of nodal lines favoring large ferrons [35]. Guided by this observation, Size of error bars denotes the numerical accuracy of energy determination in self-consistent calculations.
we have constrained the solver to search for a solution consisting of only one ferron in the simulation domain (a). We find that it has the lowest energy among all considered cases. Other identified cases, with comparable energy, are obtained if: the nodal lines preserve rotational symmetry, which can be linked to radial variant of LO state ∆(r) ∼ |∆(r)| cos(qr) (b) [45]; the transitional symmetry is imposed, like in standard LO case (c). Note that the energy differences for these configurations are small. The ferron state has already been identified as the object of enhanced stability [32]. Ginzburg-Landau approach, valid close to the critical temperature, also points to the state consisting of ferrons (called by authors ring solitons) as a candidate for the ground state [35]. Here we confirm that the state consisting of relatively large ferrons (radius by order of magnitude larger than the coherence length) is the predicted ground state in the low spin-polarization regime.
In Ref. [33], it has been shown that in 2D, the ferron size has a linear dependence on the spin imbalance. Thus, as they grow with increasing the polarization, they will start to overlap, which we define as transition to medium polarization regime. Example of such case, P = 20%, is shown in bottom panels (d-e) of Fig. 2. The obtained self-consistent solution depends on the initial random perturbation. They have comparable energies, which results in a vast possibility of different structures. We confront the energetics of disorder solutions with ordered ones. Here, the ordered structures are obtained by breaking the translational symmetry, e.g., imprinting a series of parallel nodal lines (a,b) or a nodal grid (c). Again, these ordered states can be interpreted as LO-like scenarios where the pairing field modulates as ∆(r) ∼ cos (q · r) in one or both directions [36]. The main outcome of our calculations is that one cannot discriminate which one has lower energy: LO state or disorder state? Up to numerical precision, both of them have the same energies; compare panels c and d. This suggests that the ground state in this polarization regime is either highly degenerate or there exists a large number of states with energy almost equal to the ground state energy.
Nodal lines separate regions where the pairing field changes sign. Let us define nodal area A nod = δ(arg[∆(r)] − π) dr, so it measures an area where the phase of the pairing field is shifted by π. When expressing it as a fraction of the total area A tot = (N x × N y )dx 2 the relevant values span interval [0 − 0.5]. Fig. 3 shows that for P 30%, the fraction A nod /A tot ≈ 0.5, which we define as the high polarization regime. In this regime the effective repulsion between nodal lines increases and is directly related to the gradient term in Ginzburg-Landau model. Therefore the system favors ordered lattice-like structures with equidistant distribution of nodal lines. Indeed the solutions of unconstrained calculations converge to structures where order parameter behaves like ∼ cos (q · r), at least locally, see Fig. 1. It has to be noted that rapid cooling of highly polarized system may lead to domains where the order is different due to Kibble-Zurek mechanism [46]. At the border of these domains the order parameter changes rapidly.
With increasing polarization, the size of these ordered clusters increases and eventually becomes comparable with the system size. Not surprisingly, in the high polarization regime, the energy difference between LO state and states where LO order is present only locally is indistinguishable at the numerical level. Moreover, as the structures get denser, the average distance between two nodal points becomes comparable with the coherence length, with insufficient space for the pairing to recover fully within the structures. The maximal value of the pairing strength starts to decrease with the spin-polarization, eventually reaching a critical value, and the system transits to the normal state, see Fig. 3.

IV. ROBUSTNESS OF THE RESULTS
Before we conclude, let us demonstrate the stability of the results. First, we checked the solutions' stability with respect to the temperature, see Fig. 4. The state consisting of ferrons or disordered structures survive to temperatures of the order of T c /4 − T c /3 and can compete energetically with the ordered solutions. As the temperature increases, the BCS coherence length increases as well, and structures whose sizes become comparable to this length scale melt. Except for this quantitative change, we do not find at finite temperatures any qualitative modifications as compared to the T → 0 limit.
Next, we confirmed that the similar structures are predicted if we consider 3D unitary Fermi gas and apply more elaborate energy density functional based on ASLDA. Imposing translational symmetry of the 3D system along one direction (by using ansatz ϕ n (r) → ϕ n (x, y)e ikzz ), we find that the disordered structures can be found as selfconsistent solutions of ASLDA theory, see Fig. 4. The results of 3D calculations are represented with 2D cuts in the plane perpendicular to the z axis. For example, a single ferron represented as a ring in 2D would correspond to a cylindrical tube in 3D. The quantitative differences between BdG and ASLDA approaches are related to different values of P, which separate low-medium-high polarization regimes.
Finally, we considered the box-like trap [26,40,[47][48][49][50], to verify feasibility of an experimental detection. In Fig. 5, we present the density contrast C(r) = [n ↑ (r) − n ↓ (r)]/[n is a density far away from the nodal line. For experimentally accessible temperature, of about 1/3 of the critical temperature, we find that the contrast C provides strong enough signature to identify a disordered state, assuming imaging resolution to be of order of the coherence length. These estimates do not consider imaging integration along the third direction present in experimental realization, which may decrease the contrast. For sufficiently small heights of the trapping potential that is comparable to the coherence length (as in the experiment [40]), the impact of averaging over the third dimension is expected to be negligible.

V. RELAXATION TIME SCALES
The configurations that have been investigated do not necessarily form stable minima. It is likely that the plethora of configurations can evolve into one another, and the natural question arises concerning the time scale of this evolution. Moreover, as mentioned before, the order parameter in the high polarization regime may have a different structure within certain domains. This may appear due to fast cooling (Kibble-Zurek mechanism) creating rapidly varying order parameter at the boundaries between domains. These regions store additional energy, which can be released as the system evolves towards equilibrium. In order to quantify these characteristic time scales, we have evolved the specific configurations corresponding to two regimes of medium and high polarization, as established in the previous sections. The results showing the change in the order parameter's structure within the evolution time tε F ≈ 10 3 are presented in Fig. 6. As expected, in the medium polarization regime, due to the weak interaction between nodal lines, we observe that the order parameter's structure is altered, but there is no qualitative change; the new configuration has a similar liquid crystal-like structure as the initial one. The situation is different for the high polarization regime. The initial system possessing discontinuity of the order parameter relaxes within the same time scale into a fully ordered system. Clearly, the strong repulsion between nodal lines leads to this effect. It should be noted that the final configuration in the high polarization regime shows a periodic modulation in the density and the order parameter with a global phase coherence which can be identified as a supersolid [36,51].

VI. CONCLUSIONS
We have demonstrated that the energy landscape of the spin-imbalanced system is very complex, with many configurations of comparable energies. We have explicitly demonstrated that the LOFF-type solutions represent only a tiny fraction of possible self-consistent solutions. The other identified solutions are characterized by disorder states, which cannot be parametrized in terms of simple ansatz for the pairing field, like ∆ ∼ cos (q · r). The number of possibilities for arranging the nodal lines for disorder states is vast as compared to the ordered cases. During a cooling process, where we do not have control over the order parameter distribution, the system will likely settle in the minima corresponding to the disordered state. This scenario becomes even more likely when increasing the cooling rate. Then the Kibble-Zurek mechanism will naturally favor disordered structures. For these reasons, we expect lack of perfect ordering in experimentally produced spin-imbalanced systems (especially in the low and medium polarization regime), contrary to standard predictions concerning the LOFF state detection, see [31] for a recent review. Disordered structures are moreover not stable in a sense they will evolve in time. In the low polarization limit the time evolution does not change qualitatively the order parameter distribution. The situation is different in the high polarization limit where the system evolves towards a globally regular lattice-like structure. All these issues challenge the experimental verification of the proposed scenario here. For example, averaging over many realizations will wash out the signal from the density distribution, as in each case, the system can relax towards different local minimums.
Recent studies indicate that the LOFF phase is expected to be unstable against low energy fluctuations, like paring fluctuations [18,19,52,53]. In these works, the LOFF state is defined as a state where all Cooper pairs have welldefined momentum ±q. Thus, our results are compatible with expected fragility of the LOFF, as we expect the disordered states, with no well-defined wave vector for the order parameter fluctuations, are the ones of experimental importance.
The constraint above forces the paring field to have phase shifted by π (far away from the center r > R 2 ) in comparison to the value close to the center (r < R 1 ). Note that the prescription for the ∆ is incomplete as its form is undefined for r ∈ [R 1 ; R 2 ]. In this way, we let the solver decide the optimal radius of the ferron at a given spin-polarization. After τ iterations (it), the imprinting scheme is turned off, and the code finds the ground state solution. The value of τ differs for various imprinting procedures, but it is usually approximately ten times lower than the duration of the self-consistent process. The similar methodology was applied when searching for solutions with other types of regularities.

Appendix B: Energy corrections
Our self-consistent algorithm is designed to converge to a solution with the requested number of particles N (req.) σ by independently adjusting the chemical potentials µ σ . Practically, the constraint for the particle number is satisfied only with some accuracy Although N,σ is very small, the differences in desired and obtained numbers of particles should be considered when comparing energies between states of the same spin-polarization. Hence, we introduce a correction to the energy: and uncertainty generated by this correction This uncertainty sets the size of error bars presented in the figures in the main text.

Appendix C: Simulated annealing
To diminish the probability that the self-consistent solver terminates in a local energy minimum, we applied a technique similar to well known Simulated Annealing (SA) method. It aims to take the system out of equilibrium, allowing for a significant position change in the energy landscape. The computation process is presented schematically in Fig.(8). We executed calculations for T = 0 until the algorithm found a self-consistent solution. Next, we increased the temperature to a new value and started the iteration process again. After a substantial number of self-consistent iterations, the temperature was brought back to zero. Then, the solver searched again for the ground state solution, starting from a modified finite-temperature state. The value of the finite temperature was chosen to make the smallest ferrons unstable, simultaneously leaving larger structures relatively unchanged. i ∂ ∂t u n,↑ (r, t) v n,↓ (r, t) = h ↑ (r, t) ∆(r, t) ∆ * (r, t) −h * ↓ (r, t) u n,↑ (r, t) v n,↓ (r, t) . (D1) Single-particle hamiltonian h and pairing field are defined according to formulas provided in the main text, with the main difference that now all densities also acquire time dependence. During the computation, the total energy and particle numbers' conservation laws were monitored. They were satisfied up to high precision of 2 · 10 −7 for N σ (t)/N σ (0) and 5 · 10 −6 for E(t)/E(0). A detailed description of the numerical framework for TDBdG can be found in [38].

Appendix E: Reproducibility packs
We include a set of self-explanatory examples for readers to be able to redo calculations presented in a publication on their own. We prepared supplementary files comprising all details required to reproduce three representative types of our calculations.
• BdG_P25_disordered.tar -an example of a disordered state's creation with energy minimization according to the BdG energy density functional.
• BdG_P20_grid.tar -the one where grid-like pairing potential is enforced.
Each of them requires access to considerable computing resources with W-SLDA toolkit [38,43,44] installed.