Ring solitons and soliton sacks in imbalanced fermionic systems

We show that in superfluids with fermionic imbalance and uniform ground state, there are solitons representing multiple local minima of the free energy landscape. These solitons have nontrivial soliton-soliton and soliton-vortex interactions and can form complicated bound states in the form of"soliton sacks".

Solitons have long been understood to have profound consequences for the physical properties of fermionic systems [1][2][3][4].
Recently new methods were developed to create and observe solitons in superfluid ultracold atoms [4], opening up a route to explore new regimes and properties [5][6][7][8][9]. We will focus on the existence of solitons in so-called imbalanced fermionic systems. Such superfluids exhibit pairing between fermions with different magnitude of Fermi momenta. For example such pairing has been considered in the context of dense quark matter, [10], mixtures of different ultracold atoms [11][12][13][14][15], and superconductors [16][17][18][19][20][21][22][23][24][25][26][27][28]. When the effects of imbalance are strong, the ground state of such a system can spontaneously break translation symmetry, by inducing periodic modulation in the complex order parameter. Two commonly considered inhomogeneous ground states are the Fulde-Ferrell (FF) state [29], which exhibits purely phase modulation, and the Larkin-Ovchinnikov (LO) state [30], which consists of purely density modulation. These two are jointly referred to as FFLO states. Situations where both phase and density modulate, have also been found [13,28]. Conversely if the imbalance is weak, the ground state remains uniform. In this paper we will show that fermionic imbalance can nonetheless change the properties of the system, even when the ground state is uniform, through the existence of energetically stable solitonic excitations.
The Ginburg-Landau (GL) free energy has been derived from the microscopic theory for various imbalanced systems [12,28,31,32], where fermionic population imbalance diminishes the coefficient of the second-order gradient term. This leads to the required inclusion of higher order gradient terms. We will use the GL model, originally derived in [31], which has been shown to be sufficient when describing effects far from the boundary [33,34], where the free energy density reads written in dimensionless units, where the field = | | i is the complex order parameter. The two parameters 1 and 2 are positive constants and, for systems with a two-dimensional Fermi-surface are given to be 1 = 8∕3 and 2 = 16∕3 (see supplementary material for details of rescaling). Therefore the model can be described by a single parameter , which depends on both the temperature and the population imbalance of the system.
Since the second order gradient term is strictly negative, there exists the possibility of non-uniform ground states. This can be seen in the two-dimensional model, where the configuration that minimises the free energy transitions from a uniform state , to a density-modulating LO state at LO ≃ 0.857. This inhomogeneous state minimizes the free energy until the transition to the normal state at = 4∕3. We note that the soliton solutions we find, exist for the parameter region 0.56 ≲ < LO , which is comparable to the size of the LO regime.
One feature to note regarding Eq. (1), is that it is mathematically related to the Swift-Hohenberg equation, which is also described by a fourth order partial differential equation. This model, initially formulated in the context of thermal convection [35], is commonly used to model pattern formation. The main differences between the GL model described by Eq. (1) and the standard Swift-Hohenberg equation are that the order parameter is complex-valued and there is an additional coupling between the field strength and the gradient in the terms proportional to 2 in Eq. (1). A more detailed comparison between the models can be found in the supplementary material.
We turn now to the numerical solutions of the free energy in Eq. (1). We used the nonlinear conjugate gradient flow method both in finite element (FreeFEM [36]) and finite difference schemes, which produced consistent results.
We find numerically that the two-dimensional GL free energy defined by the density in Eq. (1) has a number of local minima, in the form of solitons. The simplest subset of these soliton solutions retain the rotational spatial symmetry of the model, which we coin "ring solitons". These radial solutions take the form = ( ) i where is a constant and ( ) is a real profile function that modulates continuously between being positive and negative, before decaying to it's ground state value ±| U | as → ∞. This leads us to characterise the solutions by the number of radial nodes ( ( ) = 0) they exhibit ( ). In the cases where is constant, we assume without loss of generality that = 0 and thus is a real field. The = 1 solutions are displayed in FIG. 1 where ( ) changes sign once. It is important to note that the energy density deviation from the uniform state, plotted for = 0.7 in FIG. 1, decays such that the total soliton energy is finite. This means that, since entropy scales with system size, solitons will be thermally induced in the thermodynamic limit.
We found stable solutions for ≥ 1 ≃ 0.56, which suggests that a system with sufficiently weak imbalance will not support these solitonic excitations. As increases, the size of 7 . The lower panel shows cross-sections of the order parameter, for multiple values of , of = 1 solutions, which were found to be stable for ≥ 1 ≃ 0.56. The nodal radius increases with , such that for significantly large , the order parameter interpolates between the two ground state values ±| U |.
the soliton also increases, and the order parameter approaches | U | at the center of the soliton.
For higher values of we find that > 1 solutions become stable, first at = 2 ≃ 0.733 for = 2, followed by solutions with three and four nodal rings ( = 3, 4) at 3 ≃ 0.784 and 4 ≃ 0.806 respectively. The = 1 to 4 solutions are plotted in figure FIG. 2 along with their energies and increasing nodal radii (i.e. ( ) = 0). The = 1 solution has the lowest excitation energy above the constant ground state. As the LO transition is approached ( → LO ) the energies decrease, becoming zero relative to the constant ground state at the transition. At = LO , a state that modulates between the vacua values indefinitely become stable, similar to the LO modulating ground state. Therefore, despite only presenting solutions with four or less nodal rings, we expect that as the FFLO transition is approached, solutions with any number of concentric nodal rings become stable. Namely, for all natural numbers , there exists an such that for ∈ ( , LO ), a solitonic excitation with concentric nodal rings is stable.
Surprisingly the radial configurations of solitons that we found numerically can be approximated by a simple logistic function to reasonable accuracy. We first note that the LO state can be approximated as successive kink-like modulations. We can then approximate the transition value LO to the LO state, by calculating when it is energetically favourable for a single kink-like modulation to appear. Namely when the total en- concentric nodal circles, evaluated at the corresponding value of they become stable . The associated excitation energy − U (left) and nodal radii (right) of these solutions, as a function of , is shown in the lower panels. As increases (the LO transition is approached), the excitation energy of the solition approaches zero and its nodal radii diverges.
ergy deviation from the constant ground state, of an infinite system, minimised with respect to , of = U tanh( ) becomes zero. This was calculated to be ≃ 0.858, which is remarkably close to the numerically computed value LO ≃ 0.857. This in turn leads us to approximate the = 1 soliton in a similar way, by a radial kink-like profile = U (tanh ( + ) − tanh ( − ) − 1). If we then substitute this approximation into Eq. (1), the total energy and nodal radius, when numerically minimised with respect to and , is within an average of 1% of the true numerical solution. Nonetheless, this approximation does not capture the asymptotics of the solution well.
We can understand the stability and nature of these solitons by considering a crude approximation. Consider an = 1 soliton in the vicinity of the LO phase transition ( ≲ LO ). Inspecting the numerical solutions in FIG. 1 suggests the approximation that ≃ ±| U | everywhere except for a small finite region centred on the nodal radius . We then approximate the various terms of the energy density as being independent of , except for ∇ 2 ≃ 2 + . We reiterate that this crude approximation is valid only when the nodal radius is large due to being close to the LO transition. This gives the total excitation energy − U ∝ ∫ + + , where , and depend on . These terms then have the following physical interpretation: tot = ∫ corresponds to the energy per unit length of a straight nodal line, which in a uniform ground state is positive. Hence its contribution to the energy decreases as the radius becomes smaller. This shrinking is balanced by the term tot = ∫ (where ∝ ( ) 2 ) which represents the energy cost associated with increasing the curvature of the nodal ring ∝ 1 . These competing contributions lead to stable ring solitons with ≃ √ tot tot . Note, that → ∞ as → LO , due to tot → 0 in this limit. The argument above can be extended to > 1, demonstrating that solitons with any number of nodal rings are expected to be stable if is sufficiently close to LO .
We can understand the long range nature of the solutions by considering the linearised theory. We consider the field far from the soliton centre, such that we can write it as a small perturbation about it's ground state = U + . By assuming that any terms of order ( 2 ) or higher are negligible, we acquire the tractable linearised equation, described in detail in the supplementary material. The solution to this linearised equation gives the asymptotic form of the field as, where (1) 0 is the zeroth order Hankel function of the first kind and is a complex constant. In our model = . For our parameters we have that and are positive. Note, that → ∞ − √ cos + ∞ for → ∞, where ∞ and ∞ are some real constants. Hence the tails of the solitons decay exponentially over the length scale 1∕ , while their amplitude oscillates with period 2 ∕ . This effect is present in the states shown in FIG. 1, but becomes visible only if the axes scales are changed. One can interpret this oscillatory decay as a complex coherence length, as discussed in the context of other superconducting models [37].
As we understand the asymptotic form of the solitons we can approximate the long-range inter-soliton forces. We do this by considering two point sources that replicate the asymptotic fields of the interacting solitons and calculate the interaction energy between them [38]. The total deviation from the ground state is given by the superposition of the two asymptotic fields (given in Eq. (2)) with constants (1) and (2) , located at 1 and 2 respectively, where the distance | 1 − 2 | is large. Details of the calculations are given in the supplementary material. This gives the interaction energy which is an oscillating function of separation distance. This predicts that there will be weakly bound states with period 2 ∕ for solitons at large distances. It is interesting to compare the radial solitons reported here with the related, yet distinct, radial solutions to the stationary Swift-Hohenberg (SH) equation [39,40]. The similarity is that both solutions exhibit stable radial oscillations. Apart from the previously mentioned differences (complex order parameter and the additional terms proportional to 2 in Eq. (1)), there is a key difference in that the ground state in [39,40] is the homogeneous solution = 0, around which the solutions oscillate. In contrast the nonlinear part of our solution modulate between two antipodal points on the (1) ground state manifold. At large distances, both the SH and our solutions decay exponentially, exhibiting oscillatory tails. We find howewer that radially symmetric solutions, represent only a small fraction of the solitonic solutions in the model described by Eq. (1) . Namely, as the LO transition is approached, we find more structurally complicated solutions that break rotational symmetry, examples of which are plotted in FIG. 3. These can be interpreted as soliton sacks, where a larger soliton confines a group of smaller solitons. This confinement is a completely nonlinear effect and cannot be explained by the asymptotic inter-soliton forces. Such solutions are reminiscent of the ostensibly unrelated Skyrmion sack/bag solutions, which attract substantial interest in chiral superconductors [41,42] and magnets [43,44].
While we have demonstrated a rich spectrum of new solutions in imbalanced systems, the natural question is how they interact with the familiar soliton excitations, namely vortices. We consider a regime away from the LO instability (exemplified by the choice = 0.7) where ordinary vortex solutions exist (namely away from the regime where a vortex core induced FFLO state reported in [45]). We find that the solitons and vortices form bound states, shown in FIG. 4. The energy of this bound state is lower than the combined energy of a separate single vortex and soliton, but the energy of the bound state is larger than the energy of a single ordinary vortex.
Finally our results prompt the question of whether or not these solitons exist over the background of another imbalanced state with uniform density: namely the Fulde-Ferrell (FF) state, where the background phase modulates. An example of a microscopically derived Ginzburg-Landau model for the FF state can be found in [32]. To model the FF state, without fine-tuning, we chose parameters = 0.5, 1 = 2 and In conclusion, we have shown that fermionic imbalance leads to solitonic excitations. These solitons constitute a number of local minima of the free energy landscape and should be observable in a fluctuating, quenched or driven system. Some of the solutions we find are related, but distinctly different, to solutions of the Swift-Hohenberg equation. The solitons have nontrivial interactions, leading to multi-solitonic bound states, soliton sacks and composite vortex-soliton excitations. In ultracold atoms, such states could be created and observed by imprinting methods and standard density-sensitive techniques. In a superconductor, such solitons could be observed via scanning tunneling microscopy in superconducting films subjected to in-plane magnetic fields. the UK Engineering and Physical Sciences Research Council through grant EP/P024688/1. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Center in Linkoping, Sweden.  Supplementary Material: Ring solitons and soliton sacks in imbalanced fermionic systems

RESCALING OF MICROSCOPICALLY DERIVED GINZBURG-LANDAU ENERGY FUNCTIONAL
In this study we use the Ginzburg-Landau free energy expansion derived in [31], starting from a microscopic model for a spin-imbalanced superfluid. The resulting free energy functional reads where Δ is the superfluid order parameter and the coefficients , , , , are functions of the temperature and the Zeeman splitting . The coefficients , , and are given by where (0) is the electron density of states at the Fermi surface, is the critical temperature at = 0 and where = 1 2 − 2 and Ψ ( ) is the polygamma function of order . The remaining coefficients are given in terms of and as =̂ 2 , =̂ 4 , and =̂ 2 , where is the Fermi velocity and̂ ,̂ ,̂ are positive constants that depend on the dimensionality of the Fermi surface. The numerical values of̂ ,̂ ,̂ in one, two and three dimensions are given in TABLE S1. The possibility of inhomogeneous ground states arise in the parameter regime in which the gradient coeffient is negative. Since shares sign with the quartic coefficient , the inclusion of positive higher order terms, both in density and momentum, is necessary.
For convenience we perform the following rescaling where |Δ 0 | 2 = − 2 , 0 = 2 4 , 2 0 = − 2 and the rescaled free energỹ reads wherẽ denotes the gradient with respect tõ and 1 = 2̂ 2 and 2 =̂ ̂ ̂ . The values of 1 and 2 in one, two and three dimensions are listed in TABLE S1. In the main text, we drop the tilde notation, but still work in the rescaled model. . Note, that Eq. (S.7) becomes the static cubic-quintic Swift-Hohenberg equation, in the limit = 0.

LINEARISATION AND SOLITON INTERACTIONS
In this section we will show the technical details how the linearised solutions presented in the paper for both the long-range behaviour of the field and the inter-soliton interaction were found. From our numerical solutions it follows that there is a class of solutions which can be described by a real field, thus in the analysis below we can restrict the field to be real. The resulting equation of motion reads where ( ) = 2 − 2 4 + 6 is the potential density. However, as we are interested in the behaviour of the soliton far from it's center, we write the field as, where we assume the deviation from the uniform ground state U has only radial dependence, is real and is small. which can be written, by introducing = ∇ 2 , as a system of coupled differential equations By linear transformation to the eigenbasis of the matrix , the two equations decouple into where 2 = ± √ − 2 are the two eigenvalues of the matrix . We identify Eq. (S.14) as the two-dimensional Bessel equation, where for each eigenvalue of , the solution in general is given as a superposition of four Bessel functions. Two of the Bessel functions can be discarded directly by considering the asymptotic behaviour at → ∞ and we obtain which means that the inverse of the imaginary part of sets the decay length scale while the inverse of the real part of sets the oscillatory length scale. The small deviation is some superposition of + and − . However, since is real valued, we can use that Re (1,2) 0 (− * ) = −Re (1,2) 0 ( ) and thus reduce to a superposition of (1) 0 ( ) and (2) 0 ( * ). Lastly we use that where is some complex valued constant. Next we follow the point source approach [38] to determine the intersoliton forces. To that end let us determine the point source that replicates the asymptotic field given by Eq. (S.17). The point source is defined by the equation 1 ∇ 4 + 2 ∇ 2 + = (S. 18) and can be found by considering the limit | | = → 0, where (1) 0 ( | |) → 2i ln | |. Using that ∇ 2 ln | | = 2 ( ) we find where = + i . Having calculated the appropriate point source, we now consider two solitons, centred around 1 and 2 respectively, where the distance | 1 − 2 | is large. To estimate the interaction between the two solitons, we assume that the total field is given by the superposition ( ) = U + 1 ( ) + 2 ( ), where is the deviation from the ground state U for one soliton, centred around . That is ( ) approaches ( ) far from its center. The main assumption here is that field of the first soliton at the position of the second is the same as it would have been if there were no soliton there, and vice versa. This is quite a crude approximation since the centers contributes significantly to the interaction energy. The interaction energy int = 12 − 1 − 2 is derived by expanding to first order in the value of 1 near 2 and vice versa. After several integrations by parts we can get rid of the exact solitons 1,2 in favor of their asymptotics 1,2 , resulting in the interaction energy (1) (2) (1) 0 ( | 1 − 2 |) , (S. 23) where (1,2) are the constants associated with the two soliton asymptotics respectively. The oscillatory nature of (1) 0 implies that the solitons will be weakly bound at large distances.