# Multi-particle instability in a spin-imbalanced Fermi gas

## Abstract

Weak attractive interactions in a spin-imbalanced Fermi gas induce a multi-particle instability, binding multiple fermions together. The maximum binding energy per particle is achieved when the ratio of the number of up- and down-spin particles in the instability is equal to the ratio of the up- and down-spin densities of states in momentum at the Fermi surfaces, to utilize the variational freedom of all available momentum states. We derive this result using an analytical approach, and verify it using exact diagonalization. The multi-particle instability extends the Cooper pairing instability of balanced Fermi gases to the imbalanced case, and could form the basis of a many-body state, analogously to the construction of the Bardeen–Cooper–Schrieffer theory of superconductivity out of Cooper pairs.

## I Introduction

Attractive interactions have a long and noteworthy history as the progenitors of strongly correlated states. One of the earliest yet most profound insights was that attractive interactions between up- and down-spin electrons may induce a pairing instability, resulting in the formation of Cooper pairs Cooper56 (). These Cooper pairs then form the basis of the many-body Bardeen–Cooper–Schrieffer (BCS) theory of superconductivity Bardeen57 (); Bardeen57a (). Furthermore, even when there are unequal numbers of up- and down-spin particles in a system, Fulde and Ferrell Fulde64 () and, separately, Larkin and Ovchinnikov Larkin65 () (FFLO) showed that it is still energetically favorable for up- and down-spin particles from their respective Fermi surfaces to form Cooper pairs, leading to a strongly correlated superconducting phase in spin-imbalanced Fermi gases Casalbuoni04 (); Bowers02 (). However, the density of states in momentum at the Fermi surface of the majority-spin particles is greater than that of the minority-spin species, so the number of bound pairs that can exist is limited by the number of minority-spin particles, leaving many of the majority-spin particles at their Fermi surface unpaired and so uncorrelated. We propose a multi-particle instability that involves multiple majority-spin particles for each minority-spin particle, allowing us to utilize all of the potential of the majority-spin particles for contributing correlation energy. We find that the number of particles involved in the instability per species is proportional to the density of states in momentum at their respective Fermi surfaces. The multi-particle instability has more binding energy per particle than a Cooper pair, so could replace the Cooper pair as the building block of a superconducting state in spin-imbalanced Fermi gases.

The prototypical experimental realization of an imbalanced Fermi gas is electrons in an external magnetic field. Most superconductors are destroyed by an external magnetic field, reverting to the normal phase. However some materials, including CeCoIn Bianchi03 () and -(BEDT-TTF)Cu(NCS) Mayaffre14 (), which are superconducting at zero magnetic field, with increasing field undergo a phase transition into an exotic second superconducting state, before a further transition into the normal phase. Other materials, including ErRhB Prozorov08 () and ErNiBC Canfield96 (), display overlap of ferromagnetism and superconductivity at zero applied field, and it has recently been suggested that BiSrCaCuO exhibits some characteristics of an FFLO-like phase in the pseudogap regime Hamidian16 (). Further possible realizations of FFLO superconductivity in spin-imbalanced Fermi gases include an ultracold atomic gas of fermions trapped in one dimension that displays a transition between superconducting phases Liao10 (), or a spin-orbit coupled superconductor with imbalanced Fermi surfaces Zhang13 (); Lo14 (). However, the exotic superconducting state has not been fully characterized in any of these systems, leaving the true nature of the ground state an open question.

We follow the prescription of Cooper Cooper56 () to study a multi-particle instability on top of the Fermi surfaces. Working in second quantization notation, we construct a trial wavefunction for a multi-particle instability of several majority-spin particles binding to a (potentially smaller) number of minority-spin particles to make the binding energy per particle larger than for a Cooper pair. The optimal ratio for the number of majority- to minority-spin particles is found to be the ratio of the densities of states in momentum at their respective Fermi surfaces.

To verify our multi-particle instability we analyze the system with exact diagonalization. We confirm that our second quantized wave function captures the crucial correlations of the exact solution, expose additional insights into the structure of the wavefunction, and verify our conclusion that the optimal number of particles in the instability is set by the ratio of the densities of states in momentum.

## Ii Theory

To explore the possibility of the multi-particle instability we study a two-spin fermionic system with an attractive contact interaction at zero temperature. The BCS Hamiltonian takes the form

(1) |

where is the spin index, is the single-particle dispersion for spin species and momentum , () is a creation (annihilation) operator for a fermionic particle, and is the strength of the attractive contact interaction. In the absence of interactions the ground state of the Hamiltonian in Eq. (1) is a filled Fermi sea,

(2) |

with species-dependent Fermi energies (corresponding to Fermi momenta ) and being the vacuum state. Without loss of generality we fix the number of particles in the Fermi sea of the up-spin species to be greater than or equal to that of the down-spin species. We follow the prescription of Cooper Cooper56 () and assume that the non-interacting ground state remains undisturbed for , and focus only on a few-particle instability at . We work in a general number of dimensions.

### ii.1 Fermi surface arcs

The idealized conceptual situation where we expect a multi-particle instability to be present consists of two Fermi surfaces for the different species that are different sizes, but otherwise geometrically similar. In Fig. 1 we show example Fermi surfaces for the up- and down-spin particles in dimensions. Above the Fermi surfaces are the unoccupied momentum states that can host the multi-particle instability, which for typical phonon-mediated interactions extend over a species-independent Debye momentum . We assume that , as for many conventional superconductors Kok37 (); Reynolds51 (); Horowitz52 ().

To construct the trial wavefunction for the multi-particle instability, we start by developing multi-particle basis states. To capture all possible correlations in the system, we require that the interaction term in the Hamiltonian can couple the different basis states. As the interaction term conserves momentum, all basis states must have the same total momentum. To construct these basis states, first consider the Cooper pair situation with only one up-spin and one down-spin particle in the instability. We start with a basis state that has both particles on their respective Fermi surfaces on opposite sides of the Fermi seas (at the momenta labeled and in Fig. (a)a). In systems with anisotropic Fermi surfaces, like many of the candidate systems for FFLO Aebi94 (); Shishido03 (); Lortz07 (), the Cooper pair (and, later, the multi-particle instability) will be dominated by the lowest-curvature parts of the Fermi surface, and so in a general dispersion we place the initial basis state at the points on the Fermi surfaces with the lowest curvature.

If we move away from these starting momenta, tangentially to the Fermi surfaces by equal and opposite momenta for the different species to conserve momentum, we eventually reach the Debye momentum above the Fermi surfaces where there are no more momentum states accessible via the interaction term (reach the outer edge of the shaded regions in Fig. (a)a). The tighter curvature of the down-spin species means we will first run out of allowed momentum states for the down-spin species (at the point in Fig. (a)a). The angular width of the allowed down-spin momentum states thus sets the angular width of the up-spin momentum states for Cooper pairs.

We refer to the allowed momentum states for the particles as forming ‘arcs’ on the Fermi surfaces. An idealized version of the available momentum states for the down-spin species is indicated in Fig. (a)a by the arc above the down-spin Fermi surface bounded by blue lines, with angular width . The corresponding up-spin species arc is shown bounded by red lines.

Because it was the down-spin species that exhausted its available momentum states first in the Cooper pair situation in Fig. (a)a, we wasted the opportunity for some up-spin species momentum states to become involved in the instability and so lower the energy of the system. We can make use of twice as many up-spin momentum states by duplicating the arc of up-spin momentum states that were available in the Cooper pair situation, offsetting the arcs so they do not intersect as required by Pauli exclusion, and placing a particle in each arc. If we allow either up-spin particle to interact with the down-spin particle, we have increased the variational freedom in the system and would generically expect the binding energy to become larger. Two such Fermi surface arcs for the up-spin species are shown in Fig. (b)b, bounded by red lines.

We can generalize the above argument to include more than two up-spin and more than one down-spin particles: in general we may have up-spin arcs and particles, and down-spin arcs and particles. However, if we include too many particles, the gradients of the Fermi surfaces of the different species will differ radically at the extremal Fermi surface arcs, and it will not be possible to move around one species’ arc without immediately pushing the other species out of their allowed momentum states. We bound the maximum extent of the Fermi surface arcs by noting that when the tangents to the species’ Fermi surfaces are parallel it is possible to move particles of both species simultaneously without either being forced from their allowed momentum states. For dispersions with inversion symmetry this is achieved when the total angular widths of the two species’ arcs are equal, shown in Fig. (b)b with the total width of the arcs of both species taking the value .

The densities of states in momentum of the occupied arcs, , describe the availability of momentum states throughout all arcs for each species, in our two-dimensional example being proportional to . The density of states in momentum per particle is then . The species with the smaller value of this ratio limits the angular size available for each particle to move in, and we refer to this as the ‘critical’ species, with particles and density of states in momentum .

To show that an instability with multiple particles in separate Fermi surface arcs is the energetically favorable solution for a broad class of spin-imbalanced systems, we follow the approach of Cooper Cooper56 () to construct a variational wavefunction for the multi-particle instability. We demonstrate that in spin-imbalanced systems the multi-particle instability gives an improved binding energy over traditional Cooper pairs.

### ii.2 Basis states

To formalize the above description of the Fermi surface arcs, we label the angular center of each arc, on the Fermi surface, by a -vector . These -vectors therefore satisfy and . For and so small the can be taken to be parallel, . All the momenta within a particular arc are described by , where the vectors indicate the positions of the particles within the Fermi surface arcs, and for small we have . This guarantees and so that the particle momenta lie near their corresponding Fermi surfaces. Examples of this labeling procedure are shown in Fig. 1.

The proposed multi-particle instability is an excitation of (,) correlated particles on top of the undisturbed Fermi seas, with each particle existing in a unique arc. This can be constructed out of basis states

(3) |

where is an matrix of particle momenta in spatial dimensions.

### ii.3 Trial wavefunction

The trial wavefunction for a system with a given set of vectors is a sum over basis states with optimizable coefficients ,

(4) |

where the sum is over all momentum components of each matrix , with the prime on the sum indicating that we only sum over such that , ensuring momentum conservation. We take the coefficients to be non-zero only if all of the momenta lie within their respective arcs of the Fermi surfaces. With Eq. (4) collapses to the trial wavefunction for a Cooper pair.

### ii.4 Kinetic energy

To find an analytic expression for the energy expectation value , we first focus on the kinetic energy term and linearize the dispersions near the Fermi surfaces, . Here is the momentum corresponding to the Fermi energy, which for small enough can be considered constant over the Fermi surface arcs, and is the derivative of the single-particle energy at the Fermi surface. For the dispersions involved, and , recall that and , and so and where is the projection of along (or equivalently, the radial component of ).

This linearity simplifies the full expression for the kinetic energy of our trial wavefunction. With kinetic energy operator , we find

(5) |

which may be simplified further by using the conservation of total momentum to define , giving

(6) |

with .

### ii.5 Potential energy

To evaluate the total energy of the wavefunction we also need to evaluate the effect of the potential energy operator . The interaction operator removes two particles, one of each spin species, from a basis state and then replaces them, having transferred momentum between them. For a general basis state there are ways of choosing the pairs of particles that are involved. We can formalize this procedure by defining

(7) |

and hence

(8) |

where the vectors form a set of standard basis vectors in particle-number space: each has one element that takes the value , with the remaining elements having value . These label the particles in the different arcs in Fig. (b)b. The effect of an outer product of a vector with a scattering vector is to construct the matrix , where the column containing is determined by the particular vector. We sum over all possible pairs of up- and down-spin particles.

### ii.6 Multi-particle instability

We are now ready to combine the effect of the kinetic and potential energies by projecting the full Schrödinger equation onto the state to calculate the energy expectation value . We find that

(9) |

which, following the approach of Cooper Cooper56 (), we divide by and sum over all to obtain

(10) |

Shifting the dummy momentum variables on the right hand side by to remove the and from the arguments of , we bring the implicit expression for the energy to the form

(11) |

where is the radial projection of . We can now separate the angular and radial parts of the sum over , and carry out the angular summation. The angular summation is limited by the critical species, giving a contribution of the density of available states , meaning that the whole sum over should be considered as over the critical species.

We can also make the substitution , which has the effect of restraining the dependence of the right hand side of Eq. (11) entirely to the parameters and the limits of the sums over . However, the momentum accounts only for single momentum-transfer events, which following the prescription of Cooper theory have a maximum radial width in momentum of the Debye momentum .

The maximum kinetic energy of a basis state is obtained when each particle is at the upper end of its Fermi surface arc, giving a total kinetic energy , and the minimum kinetic energy is obtained when each particle is at the bottom of its arc, for . The summation over may be extended to cover this range, giving an implicit expression for the energy of

(12) |

The only dependence on the in the implicit expression for the energy is in the coefficients , so we can factorize out from both sides of Eq. (12). We have also removed all dependence on from the expression, and so can explicitly carry out those summations to give a factor of . This leaves us with

(13) |

analogous to Eq. (4) of Cooper’s original paper Cooper56 ().

We have reduced the complexity of the multi-particle instability to a single summation with a multiplicative constant. In the same manner as Cooper’s original analysis we may now convert this summation to an integral and solve, finding the binding energy

(14) |

In the weakly interacting limit this binding energy simplifies to

(15) |

similar to the familiar form of the binding energy of a Cooper pair.

We wish to identify the number of particles in the energetically optimal multi-particle instability. The strongest dependence of the binding energy in Eq. (15) on and is in the exponential, with the binding energy being maximized when the argument of the exponential function is least negative. The values of which achieve this, and are therefore the optimal solutions for the system, can be deduced by symmetry to satisfy the relation

(16) |

i.e. the number of particles involved in the wavefunction per spin species is proportional to the density of states in momentum. This means that all of the available momentum states are involved in the instability, and so contributing the maximum possible binding energy. Eq. (16) suggests that in dimensions a multi-particle instability is energetically favorable over conventional pair instabilities in a spin-imbalanced system.

In the next subsection we shall analyze our expression for the binding energy in the light of Eq. (16), which gives a definite prediction for the energetically optimal instability in different systems. We shall then return to the trial wavefunction given by Eq. (4), and look further at its properties and limits.

### ii.7 Binding energy analysis

To build our intuition for the expression for the binding energy of the multi-particle instability found in Eq. (14), we now examine the binding energy as a function of the ratio of the number of particles . We render the binding energy dimensionless by normalizing by , the interaction strength; , the maximum interaction momentum; , in order to account for different system sizes; and , the number of critical species particles per density of states in momentum. Normalizing by this final ratio looks forward to the eventual creation of a many-body strongly-correlated state from multi-particle instabilities, with the number of instabilities merged being limited by the availability of critical species particles. We note, however, that at low interaction strengths the dominant term in the binding energy in Eq. (15) is the exponential, so the normalization could be chosen to be by the total number of particles without affecting the results below. This results in a measure of the binding energy per critical species particle of

(17) |

To further justify this measure of the binding energy per critical species particle, we first examine the strongly interacting limit of Eq. (14). Here, in terms of the normalized ratio of number of particles per species,

(18) |

the binding energy per critical species particle takes the simple form

(19) |

This expression is maximized at , that is when , in agreement with the expression in Eq. (16) for the weakly-interacting limit.

Away from the strongly- and weakly-interacting limits the optimal binding energy remains at . In Fig. 2 we show the binding energy per critical species particle from Eq. (14) as a function of imbalance for ratios of densities of states in momentum at an intermediate interaction strength . We take as an example system a free dispersion with in dimensions, so that , although similar results hold in other systems. The balanced system is shown by the gray line, with the conventional Cooper state, having , being the energetically optimal instability. This line is symmetric about on the log-log scale, which reflects the symmetry between spin species when . For the spin-imbalanced systems where , the energetically optimal instability is still found at , as predicted by Eq. (16). To the right of this there are too many up-spin particles in the instability, and to the left there are too many down-spin particles in the instability; this leads to the lines not being symmetric about their maxima, as in imbalanced systems including the wrong number of up-spin particles is not equivalent to including the wrong number of down-spin particles.

Having examined the result of Eq. (16) that the optimum ratio of number of particles is given by the ratio of densities of states in momentum, we now discuss the difference between instabilities with different numbers of particles, but the same ratio .

### ii.8 Instabilities with same ratio

The prediction given in Eq. (16) that the energetically optimal numbers of particles involved in the instability are related by only sets the ratio between and , but does not predict the absolute numbers of particles. To probe the effect of changing the absolute numbers of particles, we need to examine in more detail the effect of Pauli blocking.

The effect of Pauli blocking has been carefully analyzed Pogosov10 (); Pogosov11 (); Fan10 () for the product of two instabilities, and found to give only a small correction to the binding energy of two separate pairs (the correction going as the inverse of the number of available momentum states). This agreement with our result for a instability, up to small Pauli blocking corrections that vanish in the thermodynamic limit, supports our finding that the binding energy per critical species particle is independent of the total number of particles involved in the instability. We shall present further numerical evidence that captures the Pauli blocking corrections in Subsection III.4.

However, the effect of Pauli blocking will become more acute in a many-body state constructed from multi-particle instabilities. This suggests that in the limit of a large number of multi-particle instabilities in a system, instabilities with fewer total particles will be energetically favorable over instabilities with more particles but the same value of in a given system.

Having investigated the structure and binding energy of the proposed multi-particle instability, we now turn to some of its limits. We examine the conventional Cooper system, with balanced Fermi seas, and identify the predictions made for one-dimensional systems, recovering in both cases agreement with well-known results from the literature. We also briefly examine the strongly-interacting limit of the proposed multi-particle instability.

### ii.9 Cooper limit

The system studied originally by Cooper Cooper56 () is a balanced Fermi gas, and so has , , which we predict should have the optimal ratio in agreement with Cooper’s findings. Moreover, with our trial wavefunction Eq. (4) reproduces the conventional Cooper trial wavefunction Cooper56 (). Therefore, with a free dispersion the weakly-interacting binding energy given by Eq. (15) reduces to the familiar Cooper expression Cooper56 ()

(20) |

where the Debye energy and the density of states in energy .

### ii.10 One-dimensional limit

Although the discussion in previous subsections has focused on dimensions, our main prediction of also holds in dimension. Here the density of states in momentum is independent of the Fermi momentum, and so for both balanced and imbalanced systems. This suggests that a Cooper pair instability with should be energetically optimal for both balanced and imbalanced systems in dimension. This is in agreement with both analytical predictions Yang01 (); Hu07 (); Parish07 () and numerically exact calculations Feiguin07 (); Batrouni08 () that show an FFLO phase constructed from Cooper pairs is the ground state throughout a large part of the phase diagram of one-dimensional Fermi gases.

### ii.11 Strongly-interacting limit

In the limit of strong attractive interactions we expect the system to promote particles to the energy of the up-spin Fermi surface to reconstruct full rotational symmetry, similar to a breached superconductor Gubankova03 (); Liu03 (); Parish07a (). This turns the system effectively into one with balanced reconstructed Fermi surfaces, and so supporting conventional Cooper pair instabilities. In the strongly-interacting limit of a many-body theory built from Cooper pairs, the pair coherence length becomes small on the scale of the separation between pairs, and so the pairs can be considered tightly-bound dimers Massignan14 (); Nozieres85 ().

We have shown that the proposed multi-particle instability reduces to the well-studied Cooper problem in the balanced limit, and collapses to a pair instability in one dimension, which both link with previous results, and also reproduces a known result in the strongly-interacting limit. This gives us confidence that the multi-particle construction is also valid away from these limits. Having shown the strength of the formalism in reproducing these known limits, we now provide numerical evidence for the multi-particle instability being energetically optimal in a range of spin-imbalanced systems.

## Iii Exact diagonalization

### iii.1 Method

In order to provide further insights into our conclusion that the optimal ratio of number of particles in an instability is given by we turn to a numerical evaluation of the wavefunction and energy expectation value . To gain computational traction, we examine a reduced Hilbert space, taking only finitely many momentum states from the Fermi surfaces. We indicate this reduction in Hilbert space size in Fig. 3, where instead of considering all momentum states (gray points) or even all momentum states on the up- and down-spin Fermi surfaces (red and blue curves), we use just linear subsets from opposite sides of the Fermi surfaces. This allows us to focus on the angular extent of the Fermi arcs, the driving force behind the emergence of the multi-particle instability. We work in the strongly interacting limit, to minimize the effect of neglecting the radial component of the sum over momentum. We use systems with momentum states for up-spin particles, and momentum states for down-spin particles: the ratio then mimics the ratio of densities of states in momentum . Fig. 3 shows an example system with .

To numerically identify the ground state of the system of particles in a system with momentum states, we explicitly construct the combinations of particle momenta for each species, for a total of basis states. Note that we do not explicitly include the additional constraint of the separation into Fermi surface arcs used in the wavefunction Eq. (4). We then directly evaluate and diagonalize the matrix of interactions between these states, with the optimal instability being that with the most negative eigenvalue.

### iii.2 Binding energy

We investigate the dependence of the optimal binding energy on the ratio of number of particles in the instability in Fig. 4, where we plot the normalized binding energy per critical species particle against the ratio of number of particles per species, normalized by the inverse ratio of number of momentum states. This rescaling of ensures that our predicted optimal binding energies are located at , as in Fig. 2. We examine systems with different ratios of numbers of momentum states , with the lines in Fig. 4 coming from systems containing , , , and momentum states respectively.

We observe that, as predicted by Eq. (16), the optimal binding energy per critical species particle for each ratio of number of momentum states is obtained with a ratio of number of particles of . This is the principal result of our exact diagonalization investigation: our numerical study reproduces the result of our approximate analytical method.

To highlight that Cooper pairs are suboptimal in spin-imbalanced systems, we indicate the Cooper pair instability for each system in Fig. 4 with orange circles, from left to right for the , , , and systems. We note that for these Cooper pair states have lower binding energy per critical species particle than the proposed multi-particle instability, whilst for the optimal multi-particle instability is simply a Cooper pair, as predicted by Cooper Cooper56 ().

In Fig. 5 we confirm the convergence of our exact diagonalization results with respect to system size for an example ratio . The different blue lines in Fig. 5 correspond to exact diagonalization calculations of the binding energy using different numbers of up-spin particles, , and fixed . We observe a rapid convergence to the infinite size limit, with the system shown in Fig. 3 giving results within % of the infinite size limit for the ratios of numbers of particles. A slice through Fig. 5 at gives the line for in Fig. 4.

### iii.3 Fermi surface arcs

It is also illuminating to examine the wavefunctions of the energetically optimal instabilities. In Fig. 6 we show the basis states involved in the energetically optimal instability of the system. Each down-spin momentum state is part of a basis state with the two up-spin momentum states joined to it by lines of the same thickness and color. Thicker lines indicate higher weighting (larger ) of the basis states, and colours represent the separation in momentum between the up-spin species particles in the instability. The wavefunction comprises basis states that have spontaneously organized arcs of the up-spin Fermi surface: each plotted basis state has one up-spin particle in the left-hand half of the up-spin Fermi surface, and one particle in the right-hand half. This is in agreement with the use of arcs in the analytical wavefunction given by Eq. (4). In addition, the highest-weighted basis states are those at the angular center of the arcs, which are the momenta labeled in Section II.

The separation of the wavefunction into Fermi surface arcs is also indicated by the integrated weights of the basis states at each momentum state, which are shown by the small points above the up-spin momentum states and below the down-spin momentum states in Fig. 6. The integrated weights for the up-spin particles show a bimodal distribution, indicating a separation into arcs. The black lines are symmetric fits to the data points, showing the arcs to contain identical distributions of integrated weights. As expected for an system, the down-spin particle inhabits a single Fermi surface arc.

### iii.4 Instabilities with same ratio

Exact diagonalization may also be used to confirm the conclusion of Subsection II.8 that instabilities with fewer total particles are marginally energetically favorable over instabilities with the same ratio , but more particles. By examining the binding energy per particle of the simple and instabilities in balanced systems with , we observe that the instability does indeed have higher binding energy per particle at all finite interaction strengths. As predicted analytically Pogosov10 (), the difference scales as in the weakly-interacting limit, confirming the conclusion that instabilities with fewer total particles are energetically favorable in finite systems.

Our exact diagonalization results on this simplified system have supported the main claims and conclusions of the analytical arguments in Section II. The energetically optimal instability in a range of different spin-imbalanced systems has been shown to satisfy the relationship predicted in Eq. (16). The separation of the Fermi surface into arcs in the analytical wavefunction has also been justified by the emergence of such arcs in the numerical calculations, and we have provided evidence for which instabilities with the same ratio are most energetically favorable.

## Iv Discussion

We have shown that spin-imbalanced Fermi gases with attractive interactions support a multi-particle instability. The most energetically favorable instability contains up- and down-spin particles in the ratio , set by the ratio of the densities of states in momentum at the Fermi surfaces.

The proposed trial wavefunction for the multi-particle instability interpolates between the well-known Cooper wavefunction Cooper56 () in the limit of balanced Fermi surfaces and theoretical predictions Yang01 (); Hu07 (); Parish07 () of the FFLO phase in one dimension. This lends support to the contention that our trial wavefunction is also valid away from these limits.

We note that the physics presented here can be explored in few-body systems. Cold atoms in an harmonic trap Serwane11 (); Zurn12 () are an ideal system to explore few-particle physics, as the exact energy and expectation values such as the wavefunction symmetry may be directly measured Bugnion13 (); Bugnion13a (). Cold atom experiments may therefore be able to observe the scaling of the binding energy and spatial structure of the trial wavefunction proposed here, of which hints may previously have been seen numerically Bugnion13a ().

In real experiments the interaction between fermions will never be exactly the contact interaction from Eq. (1). In cold atom systems the interaction may be expanded as , where is the scattering length and is the effective range Keyserlingk13 (). Positive reduces the effective interaction strength, making the multi-particle instability less energetically favorable, whilst negative makes it more energetically favorable; however, is typically small on the scale of , and so the effect of the finite range interaction is also small. The screened Coulomb interaction, where is the Thomas-Fermi screening length, relevant for example to electron-hole systems Lozovik75 (); Perali13 (), has a similar effect, with the screening length taking on the same role as the effective range for cold-atom interactions, and so weakening the multi-particle instability relative to the purely-contact case. This weakening is also found in standard Cooper pairs, however, and so is unlikely to qualitatively change the conclusions in the manuscript. The next order term in the effective range expansion would go as : this term will discriminate between multi-particle instabilities and Cooper pairs, being a function of how many fermions near the Fermi surfaces are involved in the instability, but is not expected to have a large effect, as in our formalism both and are small.

In the same way that Cooper pairs form the conceptual basis of the Bardeen–Cooper–Schrieffer theory of superconductivity, it is expected that a many-body state may be constructed using the multi-particle instabilities presented here with even values of . By analogy to the relationship between the traditional Cooper result and the BCS order parameter, we expect that the order parameter of the future many-body superconducting theory should have a form that is reminiscent of Eq. (15). The many-body theory should not be limited to including a single type of multi-particle instability, and similarly to predictions made for the FFLO phase Bowers02 () may be constructed from multiple superposed multi-particle instabilities, forming a crystalline structure.

A natural tool to use to search for this exotic superconducting state is a spin-imbalanced ultracold fermionic gas Zwierlein06 (); Partridge06 (). This system allows fine control over the populations and interactions of the fermions, allowing experiments to focus on the potential of new physics. Previous spin-imbalanced ultracold fermionic gas experiments have used inhomogeneous optical trapping potentials, in which the region of space where multi-particle instability-based superconductivity is likely to be observable is very small. However, recent experimental developments have allowed the creation of homogeneous ultracold fermionic gases Mukherjee17 (), where the delicate novel superconducting state is likely to exist over larger regions of space, and so be easier to observe and characterize.

Such a strongly correlated state would present novel superconducting properties, including unusual Andreev reflection Andreev64 (), Josephson tunneling Josephson62 (), and SQUID Jaklevic64 () or other superconducting loop Kirsanskas15 () properties, due to the underlying multi-particle structure. With the underlying instabilities involving fermions, magnetic flux is likely to be quantized in units of , rather than for BCS superconductivity based on Cooper pairs. The superconducting order parameter would also exhibit unusual behavior, being necessarily complex due to the presence of non-antipodal -vectors, and oscillating with wavevectors , with interference due to similar -vectors giving rise to beats in the order parameter amplitude. The existence of a superconducting state constructed from multi-particle instabilities may also explain the lack of definitive observations of the conjectured FFLO state.

Data used for this paper are available online data ().

###### Acknowledgements.

The authors thank Adam Nahum, Johannes Hofmann, Johannes Knolle, Jens Paaske, Robin Reuvers, and Darryl Foo for useful discussions, and acknowledge the financial support of the EPSRC and the Royal Society.### References

- L.N. Cooper, Phys. Rev. 104, 1189 (1956).
- J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phys. Rev. 106, 162 (1957).
- J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phys. Rev. 108, 1175 (1957).
- P. Fulde and R.A. Ferrell, Phys. Rev. 135, A550 (1964).
- A.I. Larkin, Y.N. Ovchinnikov, Sov. Phys. JETP 20, 762 (1965).
- R. Casalbuoni and G. Nardulli, Rev. Mod. Phys. 76, 263 (2004).
- J.A. Bowers and K. Rajagopal, Phys. Rev. D 66, 065002 (2002).
- A. Bianchi, R. Movshovich, C. Capan, P.G. Pagliuso, and J.L. Sarrao, Phys. Rev. Lett. 91, 187004 (2003).
- H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, K. Miyagawa, K. Kanoda, and V.F. Mitrović, Nature Physics 10, 928 (2014).
- R. Prozorov, M.D. Vannette, S.A. Law, S.L. Bud’ko, and P.C. Canfield, Phys. Rev. B 77, 100503(R) (2008).
- P.C. Canfield, S.L. Bud’ko, and B.K. Cho, Physica C 262, 249 (1996).
- M.H. Hamidian, S.D. Edkins, S.H. Joo, A. Kostin, H. Eisaki, S. Uchida, M.J. Lawler, E.-A. Kim, A.P. Mackenzie, K. Fujita, J. Lee, and J.C. S. Davis, Nature 532, 343 (2016).
- Y. Liao, A.S.C. Rittner, T. Paprotta, W. Li, G.B. Partridge, R.G. Hulet, S.K. Baur, and E.J. Mueller, Nature 467, 567 (2010).
- W. Zhang, and W. Yi, Nature Communications 4, 2711 (2013).
- S.-T. Lo, S.-W. Lin, Y.-T. Wang, S.-D. Lin, and C.-T. Liang, Sci. Rep. 4, 5438 (2014).
- J.A. Kok and W.H. Keesom, Physica 4, 835 (1937).
- C.A. Reynolds, B. Serin, and L.B. Nesbitt, Phys. Rev. 84, 691 (1951).
- M. Horowitz, A.A. Silvidi, S.F. Malaker, and J.G. Daunt, Phys. Rev. 88, 1182 (1952).
- P. Aebi, J. Osterwalder, P. Schwaller, L. Schlapbach, M. Shimoda, T. Mochiku, and K. Kadowaki, Phys. Rev. Lett. 72, 2757 (1994).
- H. Shishido, T. Ueda, S. Hashimoto, T. Kubo, R. Settai, H. Harima, and Y. Ōnuki, J. Phys.: Condens. Matter 15, L499 (2003).
- R. Lortz, Y. Wang, A. Demuer, P.H.M. Böttger, B. Bergk, G. Zwicknagl, Y. Nakazawa, and J. Wosnitza, Phys. Rev. Lett. 99, 187002 (2007).
- W.V. Pogosov, M. Combescot, and M. Crouzeix, Phys. Rev. B 81, 174514 (2010).
- W.V. Pogosov and M. Combescot, Physica C 471, 566 (2011).
- J.D. Fan and Y.M. Malozovsky, J. Supercond. Nov. Magn. 23, 655 (2010).
- K. Yang, Phys. Rev. B 63, 140511(R), (2001).
- H. Hu, X.-J. Liu, and P.D. Drummond, Phys. Rev. Lett. 98, 070403 (2007).
- M.M. Parish, S.K. Baur, E.J. Mueller, and D.A. Huse, Phys. Rev. Lett. 99, 250403 (2007).
- A.E. Feiguin and F. Heidrich-Meisner, Phys. Rev. B 76, 220508(R) (2007).
- G.G. Batrouni, M.H. Huntley, V.G. Rousseau, and R.T. Scalettar, Phys. Rev. Lett. 100, 116405 (2008).
- E. Gubankova, W.V. Liu, and F. Wilczek, Phys. Rev. Lett. 91, 032001 (2003).
- W.V. Liu and F. Wilczek, Phys. Rev. Lett. 90, 047002 (2003).
- M.M. Parish, F.M. Marchetti, A. Lamacraft, and B.D. Simons, Phys. Rev. Lett. 98, 160402 (2007).
- P. Massignan, M. Zaccanti, and G.M. Bruun, Rep. Prog. Phys. 77, 034401 (2014).
- P. Nozières and S. Schmitt-Rink, J. Low Temp. Phys. 59, 195 (1985).
- F. Serwane, G. Zürn, T. Lompe, T.B. Ottenstein, A.N. Wenz, and S. Jochim, Science 332, 336 (2011).
- G. Zürn, F. Serwane, T. Lompe, A.N. Wenz, M.G. Ries, J.E. Bohn, and S. Jochim, Phys. Rev. Lett. 108, 075303 (2012).
- P.O. Bugnion and G.J. Conduit, Phys. Rev. A 87, 060502(R) (2013).
- P.O. Bugnion, J.A. Lofthouse, and G.J. Conduit, Phys. Rev. Lett. 111, 045301 (2013).
- C.W. von Keyserlingk and G.J. Conduit, Phys. Rev. B 87, 184424 (2013).
- Yu.E. Lozovik and V.I. Yudson, JETP Lett. 22, 274 (1975).
- A. Perali, D. Neilson, and A.R. Hamilton, Phys. Rev. Lett. 110, 146803 (2013).
- M.W. Zwierlein, A. Schirotzek, C.H. Schunck, and W. Ketterle, Science 311, 492 (2006).
- G.B. Partridge, W. Li, R.I. Kamar, Y. Liao, and R.G. Hulet, Science 311, 503 (2006).
- B. Mukherjee, Z. Yan, P.B. Patel, Z. Hadzibabic, T. Yefsah, J. Struck, and M.W. Zwierlein, Phys. Rev. Lett. 118, 123401 (2017).
- A.F. Andreev, Sov. Phys. JETP 19, 1228 (1964).
- B.D. Josephson, Phys. Lett. 1, 251 (1962).
- R.C. Jaklevic, J. Lambe, A.H. Silver, and J.E. Mercereau, Phys. Rev. Lett. 12, 159 (1964).
- G. Kiršanskas, M. Goldstein, K. Flensberg, L.I. Glazman, and J. Paaske, Phys. Rev. B 92, 235422 (2015).
- T.M. Whitehead and G.J. Conduit, Cambridge University DSpace repository, https://doi.org/10.17863/CAM.15832.