Normal and anomalous solitons in the theory of dynamical Cooper pairing

# Normal and anomalous solitons in the theory of dynamical Cooper pairing

Emil A. Yuzbashyan Center for Materials Theory, Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854, USA
###### Abstract

We obtain multi-soliton solutions of the time-dependent Bogoliubov-de Gennes equations or, equivalently, Gorkov equations that describe the dynamics of a fermionic condensate in the dissipationless regime. There are two kinds of solitons – normal and anomalous. At large times, normal multi-solitons asymptote to unstable stationary states of the BCS Hamiltonian with zero order parameter (normal states), while the anomalous ones tend to eigenstates characterized by a nonzero anomalous average. Under certain circumstances, multi-soliton solutions break up into sums of single solitons. In the linear analysis near the stationary states, solitons correspond to unstable modes. Generally, they are nonlinear extensions of these modes, so that a stationary state with unstable modes gives rise to a -soliton solution. We relate parameters of the multi-solitons to those of the asymptotic stationary state, which determines the conditions necessary for exciting solitons. We further argue that the dynamics in many physical situations is multi-soliton.

## I Introduction and summary of the results

Recent years have witnessed renewed experimental and theoretical interest in far from equilibrium phenomena in strongly interacting many-body systems at low temperatures. Examples include non-stationary Kondo and other impurity models[Goldhaber-Gordon, Mehta, ], quenched Luttinger liquids[Schumm, Gritsev, ], electron spin dynamics induced by hyperfine interactions[Ono, Al-Hassanieh, ] etc. On the theory side, there is a considerable effort to develop new approaches to nonequilibrium many-body physics. This presents a significant challenge as conventional techniques are often inadequate for the description of these phenomena. In particular, there have been major advances in the theory of dynamical fermionic pairing in the collisionless regime[galaiko, spectroscopy, ]. This problem is long known to be accurately described by the time-dependent Bogoliubov-de Gennes equations, which in this case are a set of coupled nonlinear integro-differential equationsvolkov (); galperin (); single (); gensol1 (). Nevertheless, it was not until recently that this nonlinear system was realized to be exactly solvablegensol (); gensol1 (); dicke (). The exact solution proved to be a unique approach to the problem of dynamical pairing and has been extensively exploited to obtain analytical information about its key physical properties. For example, a nonequilibrium “phase diagram” of a homogeneous Bardeen-Cooper-Schriffer (BCS) superfluid with a number of novel phases, as well as their responses to existing experimental probes were predicted analyticallyemilts (); Levitov2006 (); Emil2006 (); spectroscopy ().

However, while much attention was focused on the asymptotic states of the condensate at large times, the transient dynamics has not been fully explored. Most importantly, nonlinear integrable systems are known to exhibit a remarkable class of multi-soliton solutions that play a central role in understanding and predicting their properties. Physical solutions can often be represented as a superposition of solitons making a quantitative analysis possible. For instance, one can show that the dynamics giving rise to nonequilibrium “phases” mentioned above is multi-soliton, see below. The existence and properties of solutions of this type are also of a general interests from the point of view of nonlinear physics due to the nonlocal nature of the BCS problem distinguishing it from familiar integrable systems, such as nonlinear Shrödinger, Korteweg-de Vries, sine-Gordon etc.

In this paper we construct multi-soliton solutions to the dynamical fermionic pairing problem, see Figs. 15. We establish a direct correspondence between solitons and the stationary states of the mean-field BCS Hamiltonian, such that each soliton solution asymptotes to an eigenstate at times . This also identifies conditions necessary for exciting solitons. There are two distinct types of solitons – normal and anomalous. Normal solitons asymptote to stationary states that are simultaneous eigenstates of a Fermi gas and the mean-field BCS Hamiltonian and are characterized by a zero anomalous average. For anomalous solitons the asymptotic value of this average is finite, . As the separation between solitons is increased, the multi-soliton solution breaks up into a simple sum of single solitons (see e.g. Figs. 2 and 3)– one of the defining properties of solitons.

In the rest of this section, we briefly formulate the problem and then summarize the main results. The collisionless dynamics of a fermionic superfluid can be described by the Bogoliubov-de Gennes equationssingle (); gensol1 ()

 iddt(UmVm)=(ϵmΔ(t)Δ∗(t)−ϵm)(UmVm) (1)

where is the anomalous average, are the single fermion energies relative to the Fermi level , and is the coupling constant. These nonlinear equations are known to be integrable for any number of Bogoliubov amplitudes gensol (); gensol1 (); dicke (). In the continuum limit the summation in the expression for is replaced by integration and Eqs. (1) become integrable nonlinear integro-differential equations. Each solution of Eq. (1) yields a many-body wave function

 |Ψ(t)⟩=∏m[U∗m(t)+V∗m(t)^c†m↑^c†m↓]|0⟩, (2)

where the product is taken only over unblocked levels – levels that are either unoccupied or doubly occupied.

As we demonstrate in subsequent sections, solitons tend to unstable stationary states of the mean-field BCS Hamiltonian at . The mean-field Hamiltonian has two types of eigenstates – normal and anomalous. Anomalous stationary states have a nonzero constant value of the anomalous average and are solutions of Eq. (1) of the formBCS (); Tinkham () , where and are the eigenvectors and eigenvalues of the matrix on the right hand side of Eq. (1). There are two states for each . The BCS ground state has for all . A state where while describes a single excited pairBCS (); ander1 () and has energy above the ground state. We note also that an excited pair introduces a discontinuity in the average fermion occupation number , since changes sign at . Normal eigenstates have and amplitudes equal to either or . Both types of eigenstates are exact stationary states of the mean-field BCS Hamiltonian

 ^H=∑j;σ=↓,↑ϵj^c†jσ^cjσ−∑j,k(Δ^c†j↑^c†j↓+h.c.), (3)

where , and and are the creation and annihilation operators for the two fermion species. Normal states are also eigenstates of the Fermi gas – the first term in Eq. (3). For example, the Fermi ground state is a normal eigenstate with for and for , i.e. all single particle states below the Fermi level occupied and states above it empty.

A linear analysis of Eq. (1) around stationary states shows that some of them are unstableElihu (); single (). These states give rise to solitons, as is typical in integrable nonlinear dynamics. For example, a simple pendulum displays a soliton solution when started in its unstable equilibrium with zero velocity. In the phase space the soliton is the separatrix connecting the unstable equilibrium to itself. The same is true for example for the single soliton solution of the Korteweg-de Vries equationArnold (). Similarly in the present case the unstable modes start to grow exponentially and become solitons due to nonlinear effects. In a certain sense, solitons can be viewed as nonlinear extensions of the unstable modes.

Now let us summarize soliton solutions of the Bogoliubov-de Gennes equations (1). Detailed derivation and discussion of these as well as some other solutions can be found in subsequent sections. Here we present only the results for the amplitude of the order parameter .

### i.1 Normal solitons

First, we present normal multi-solitons derived in Sec. IV. The 1-soliton solution of a normal type has been previously obtained in Ref single, . The amplitude of the order parameter has the following form:

 |Δ(t)|=2γcosh(2γt+α), (4)

where is a real parameter. At the corresponding wave function (2) asymptotes to the Fermi ground state. The latter is unstable in the presence of the pairing interaction. Indeed, a linear analysis of equations of motion (1) around the Fermi ground state shows a single unstable normal mode, which grows as Elihu (); single () as can be seen from Eq. (4) in the limit. The plot of Eq. (4) is a single peak centered at , see Fig. 1. Its height (the amplitude of the soliton) and width are controlled by the parameter . In the present case , where is the ground state BCS gap. This parameter can be interpreted as an imaginary “frequency” of the unstable normal mode, while Eq. (4) as an extension of this mode to the nonlinear regime. Below we will see that the number of unstable modes and consequently the number of solitons corresponding to a given stationary state is related to the number of discontinuities in the average fermion occupation number in this state. Specifically, discontinuities (jumps) in lead to up to coupled normal solitons. The Fermi ground state has a single discontinuity at the Fermi level, giving rise to the 1-soliton solution.

The 2-normal-solitons have considerably richer structure. Let us give two examples. Both are described by

 |Δ(t)|=A∣∣ ∣∣h(t)h(t)¨h(t)−˙h2(t)∣∣ ∣∣, (5)

with different choices for the amplitude and function . One choice is and

 h(t)=eiϕ1cosh(2γ1t+α1)2γ1+eiϕ2cosh(2γ2t+α2)2γ2. (6)

Examples of for this case are plotted in Fig. 2a. The second option is and

 h(t)=e−2iμt+iϕ1cosh(2γt+α1−iβ)2γ+e2iμt+iϕ2cosh(2γt+α2+iβ)2γ, (7)

where is the phase of , i.e. . Fig. 2b shows graphs of obtained using Eq. (7). In both cases the wave function asymptotes to a normal eigenstate at . This eigenstate is obtained from the Fermi ground state by moving all fermions in the energy interval below the Fermi level to the symmetric interval above it. Eqs. (6) and (7) correspond to different values of . Note that this normal eigenstate has discontinuities at , 0, and in (see the insets in Fig. 2) thus leading to solitons. A linear analysis around this state yields two unstable modes that exponentially grow with rates in the first case and in the second. Which of the two cases is realized depends on the ratio . The values of or are also fixed by and , while and are arbitrary real parameters. They control the separation between solitons and their relative phase in the 2-soliton solution, respectively. For sufficiently large Eq. (5) always breaks up into a sum of two single solitons of the form (4). The graph of in this case displays two well separated peaks, each representing a single soliton, see Fig. 2.

The general -soliton solution of the normal type has the form (see Sec. IV.1)

 |Δ(t)|=22k−1∣∣∣Dk−1Dk∣∣∣, (8)

where is the following determinant:

 Dr=∣∣ ∣ ∣∣f…f(r−1)⋮⋮f(r−1)…f2(r−1)∣∣ ∣ ∣∣, (9)

is the th derivative of the function with respect to , and

 f(t)=2k∑j=1A(cj)e−2icjt∏m≠j(cj−cm), (10)

The set of complex parameters (“frequencies” of the unstable modes) is complex conjugate to itself. Let us order this set so that and for . The constants and are related as follows

 A(cl)=eαl+iϕl,A(ck+l)=−e−αl+iϕl, (11)

where and are arbitrary real parameters. The single soliton (4) is obtained from Eq. (8) by setting and . The 2-soliton corresponds to and or to get Eq. (6) or (7), respectively. See also Fig. 3 for examples of 3-normal-solitons.

At the -normal-soliton tends to a normal eigenstate that has at least discontinuities in the distribution function . Linearizing the Bogoliubov-de Gennes equations around this state at large negative , one obtains unstable normal modes that grow as for . When the differences between parameters are large the -soliton solution (8) breaks up into a sum of single solitons,

 |Δ(k)(t,{ci,αi,ϕi})|≈k∑i=1|Δ(1)(t,Im(ci),αi)|, (12)

where denotes the -normal-soliton (8) and stands for the single soliton (4) with and . In this case, the plot of shows well separated peaks (individual solitons) as illustrated in Figs. 2 and 3. set the amplitudes and widths of individual solitons, are the frequencies with which they “rotate” with respect to one another as in Eq. (7), where , and and determine the separation between the solitons and their relative phases, respectively.

### i.2 Anomalous solitons

Next, let us summarize the results of Sec. V for anomalous solitons. A single anomalous soliton has the following form (see also Fig. 4):

 Δ(t)−Δa=2(γ2−Δ2a)Δa±γcosh(λt+α), (13)

where and is an arbitrary real parameter as before. As the state of the system tends to an anomalous eigenstate with the value of the BCS order parameter equal to . In this state, all pairs in a certain energy interval around the Fermi level are excited, i.e. for . As a result, the distribution function has two jumps at (inset in Fig. 4). A linear analysis around this anomalous eigenstate shows a single unstable mode that grows as . In general, jumps in the distribution function of a stationary anomalous state give rise to up to anomalous solitons. Note also that the anomalous soliton (13) generalizes the normal one (4) and turns into it when .

A general -anomalous-soliton can also be constructed within our approach. However, here we present only an example of a 2-anomalous-soliton solution

 Δ(t)−Δa=h(t)2[h(t)¨h(t)−˙h2(t)], (14)

where

 h(t)=Δaλ21λ22±γ1cosh(λ1t+α1)λ21(λ21−λ22)±γ2cosh(λ2t+α2)λ22(λ22−λ21), (15)

The signs can be chosen independently of each other. As in the case of the 1-anomalous-soliton, at large times the wave function asymptotes to an anomalous stationary state with order parameter . This state has two unstable modes that grow exponentially with rates . The graph of the 2-anomalous-soliton solution (14) consists of two peaks or dips depending on the choice of signs in Eq. (15), see Fig. 5. The parameters take arbitrary real values. Their difference, , determines the separation between the peaks in time. Similarly to Eq. (12), at large separations the -anomalous-soliton turns into a sum of individual solitons of the form (13)

 Δ(k)(t,{λi,αi})−Δa≈k∑i=1[Δ(1)(t,λi,αi)−Δa], (16)

where is the single anomalous soliton (13) with and .

The rest of the paper is organized as follows: in the next section, we review the basic setup of the problem and the tools (Lax vector and separation variables) necessary for deriving solitons. In Sec. IV, we perform linear analysis of equations of motion around normal and anomalous stationary states. This section also provides examples of normal and anomalous eigenstates that give rise to one and two normal and anomalous solitons. Sections IV and V are devoted to a detailed derivation of soliton solutions and a discussion of their main properties.

## Ii Review of the basic setup and relevant previous results

In this section we discuss the basic setup of the problem and introduce our notation (see Refs. single, ; gensol1, for more details). We also review the properties of the exact solutiongensol (); gensol1 (); dicke (); emilts () of the equations of motion needed for obtaining and analyzing the multi-soliton solutions summarized in the previous section.

### ii.1 Notations and basic equations

Here we review the model Hamiltonian and its mean-field equations of motion (1). The latter can be reformulated as equations of motion for classical spins (angular momenta) – this is the form we will be primarily using. We also describe normal and anomalous stationary states in terms of the spin variables.

The dissipationless dynamics of a fermionic superfluid can be modeled by the reduced BCS Hamiltonianvolkov (); galperin (); single (); BCS ()

 ^H=∑jϵj^nj−g∑j,k^c†j↑^c†j↓^ck↓^ck↑, (17)

where . This description is valid in the weak coupling regime at times shorter than the energy relaxation time and for a system of size smaller than the BCS coherence length . Under these conditions, the BCS order parameter is uniform in space, the interaction matrix elements can be evaluated at the Fermi energy yielding a single coupling constant that is independent of and . The summations in Eq. (17) over and are restricted to single particle energies and , where is an ultraviolet cutoff for the pairing interaction. In metallic superconductors , where is the Debye frequency. For atomic fermions . The offdiagonal interactions – terms of the form with or can be neglected, since they are relevant only at times .

The validity of the mean-field approach is rooted in the fact that each pair creation operator in Eq. (17) interacts with the collective pairing field , which is expected to deviate little from its quantum mechanical average . For example, the mean-field is known to be exact for the description of the low-energy properties of the Hamiltonian (17) in the limit ander1 (); rich (); ander (), where and are the mean spacing between the single particle levels and the ground state gap, respectively. Note that the conditions and are compatible in the weak coupling regime.

We are interested in solving the Heisenberg equations of motion for Hamiltonian (17) to determine the evolution of various correlators, e.g. , , and . In mean-field approach, we replace the operator in the Heisenberg equations with its quantum mechanical average

 Δ(t)=g∑k⟨^ck↓(t)^ck↑(t)⟩. (18)

Further, introducing

 2szm=⟨^nm⟩−1,s−m≡sxm−isym=⟨^cm↓^cm↑⟩, (19)

we obtainander1 ()

 ˙sm=bm×sm,bm=(−2Δx,−2Δy,2ϵm), (20)

where and are the real and imaginary parts of . In terms of and , Green’s functions at coinciding times, Eqs. (20) are well-known Gorkov equationsGorkov (); volkov (). The above procedure leading to Gorkov equations is essentially equivalent to taking the time-dependent wave function of the system to have the product form (2) at all times. Then, the Schrödinger equation takes the form of the Bogoliubov-de Gennes equations (1), which are in turn equivalent to Eq. (20) with

 2szm=|Vm|2−|Um|2,s−m=UmV∗m. (21)

Equations of motion (20) are Hamilton’s equations for the following classical spin (interacting angular momentum) model:

 Hcl=n∑j=12ϵjszj−gn∑j,k=1s+js−k, (22)

with the usual angular momentum Poisson brackets etc. The summations in Eq. (22) are restricted to the subspace of unblocked (with occupation numbers ) single fermion levels . The Hamiltonian (17) does not have matrix elements connecting the unblocked levels to blocked ones (). The latter are decoupled and their occupation numbers are conserved by the evolution. Note also that Eq. (1) conserves the norm and therefore the length of the spins is fixed .

The normal and anomalous eigenstates of the BCS Hamiltonian discussed in the Introduction correspond to equilibrium spin configurations where each spin is either parallel or antiparallel to the field ander1 (). According to Eq. (21), every such arrangement of spins uniquely determines an eigenstate of the form (2) and vice versa. The anomalous eigenstates yield

 2szm=−emϵm√ϵ2m+Δ2a,2sxm=−emΔa√ϵ2m+Δ2a,em=±1, (23)

where the axis has been chosen so that the stationary value of the order parameter, , is real. The factor if the spin is parallel to the field and otherwise. The self-consistency condition reads

 ∑mem√ϵ2m+Δ2a=2g. (24)

This is the BCS gap equation, which determines the value of in the anomalous state. The configuration of spins with all is equivalent to the BCS ground state. In this case – the ground state gap – and Eq. (24) becomes in the continuum limit

 ∫D0dϵ√ϵ2+Δ20=2λλ=gνFV≡gδ, (25)

where is the density of states at the Fermi level and is the volume of the system. In Eq. (25) and throughout this paper we assume the weak coupling regime and a constant density of , , in the continuum limit. A non-constant density of states modifies the value of determined from Eq. (25) but will not affect any other equations derived in the rest of the paper. As we will see, these equations are confined to energies of order , while the density of states varies on an energy scale of order or larger. Using , we obtain from Eq. (25)

 Δ0=2De−1/λ. (26)

The configuration with only one flipped spin, and , corresponds to an excited state – it contains an excited pair and has energy relative to the ground state. Similarly, having two spins parallel to the field is equivalent to an eigenstate with two excited pairs etc.

Normal eigenstates are spin arrangements where each spin is along axis, i.e.

 2szm=±1≡lm,s−m=0. (27)

They are also equilibria of the classical Hamiltonian (22) according to Eq. (20). The Fermi ground state corresponds to , while other normal eigenstates can be obtained from this state by moving pairs of fermions from levels below the Fermi energy to levels above it and flipping the corresponding spins.

### ii.2 General properties of the dynamics

In this subsection we introduce the Lax vector constructiongensol1 (); dicke (), which plays a central role in analyzing the dynamics of the BCS Hamiltonian. We also define the separation variables and describe the general features of the dynamics.

The dynamics of the classical Hamiltonian (22) or, equivalently, Eqs. (20) and (1) turn out to be integrable. A convenient tool for their analysis is the Lax vector defined as

 L(u)=−^zg+n∑m=1smu−ϵm, (28)

where is a complex parameter, is a unit vector along axis, and is the total number of spins. The length of this vector is conserved by Eqs. (20) for any , i.e.

 dL2(u)dt=0. (29)

For this reason can be viewed as the generator of the integrals of motiongensol1 () for Eqs. (20). For example, its zeroes are conserved and constitute a set of independent integrals. Another possible choice for the integrals is e.g. the set of the residues of at the poles at . Note that

 L2(u)=Q2n(u)g2∏j(u−ϵj)2, (30)

where is a (spectral) polynomial of order . We also have

 L2(u)=L2z(u)+L−(u)L+(u), (31)

where and are the components of the Lax vector .

To obtain solitons, we need to introduce new dynamical variables sklyanin (); vadim () in which Eqs. (20) separate and can be integrated. The separation variables are defined in terms of the “old” dynamical variables as solutions of the following equation:

 L−(um)=n∑j=1s−jum−ϵj=0. (32)

This equation has solutions since, when is brought to a common denominator, its numerator is a polynomial of order . Consequently, there are separation variables . Eq. (32) can be inverted to obtain the spins in terms of the separation variables as follows

 s−j=J−∏k(ϵj−uk)∏k≠j(ϵj−ϵk), (33)

where as usual and is the total classical spin.

In terms of the separation variables Eqs. (20) read

 ˙uj=2i√Q2n(uj)∏m≠j(uj−um),j=1,…,n−1, (34)
 ˙J−=−2iJ−(n∑j=1ϵj+gJz2−n−1∑m=1um). (35)

An important observationemilts () is that main properties of the dynamics can be effectively discerned by analyzing the zeros of . According to Eq. (30), these are the roots of the spectral polynomial and we will often refer to their configuration in the complex plane as to the root diagram of . Since is positively defined, it has pairs of complex conjugate roots. For generic initial conditions all roots are distinct. In this case the dynamics of the system is quasiperiodic with incommensurate frequencies and any dynamical quantity, e.g. the order parameter , typically contains all frequencies.

Significant simplifications occur when some roots are degenerategensol (); dicke (). It is important to distinguish between real and complex double roots. Note that any real root of is automatically a double root (zero) because is positively defined. A real zero of must also be a zero of all three components of note1 (). Further, note from Eq. (32) that one of the separation variables must coincide with . In other words, it must be time-independent as it is “frozen” into the real root . Eq. (34) shows that this is an allowed solution of the equations of motion for the separation variables. This freezing of a separation variable can be translated into a genuine reduction of the number of degrees of freedom by one so that the dynamics of the Hamiltonian (22) with spins reduces to that of the same Hamiltonian but with spins. In general, real zeros (or equivalently complex zeros) mean a reduction of the dynamics to that of effective spins, see Ref. gensol, and dicke, for details. Below we will often encounter a situation when has a number of real zeros and consequently a number of separation variables are frozen. The remaining variables we call unfrozen.

Let be the separation variable frozen into the double zero of . Consider Eq. (34) for . Both the numerator and the denominator of the right hand side contain a factor , which cancels lowering the order of the polynomial under the square root by two. Suppose has double real zeros. Then, there are unfrozen separation variables . For these variables Eq. (34) can be brought to the following formgensol () with the help of Eq. (30):

 2k−1∑j=1uljduj∏m(uj−ϵm)√˜L2(uj)=2igdtδl,2k−2,l=0,…,2k−2, (36)

where is obtained from by removing all real zeros , i.e.

 ˜L2(u)=L2(u)∏m(u−cm)2.

## Iii Linear analysis around stationary states

In this section, we analyze equations of motion linearized in the vicinity of normal and anomalous stationary states. We show that the separation variables are the normal modes of the linearized problem. Some of the stationary states are unstable. As we will see in the next section, the corresponding normal modes become solitons in the nonlinear regime.

### iii.1 Frequencies of oscillations around stationary states

Here we show that the frequencies of small oscillations around normal and anomalous states are determined by the zeros of , see also Ref. emilts, . When one of the frequencies is complex, the state is unstable and the corresponding mode grows exponentially.

The linear analysis of equations of motion (20) around stationary states greatly simplifies in terms of separation variables. According to Eq. (34), stationary positions of are the roots of the polynomial (or equivalently the zeros of , see Eq. (30)). Let us determine the form of in the stationary states. Consider first the anomalous states (23). Using Eqs. (23), (28), and (24), we obtain

 L(u)=(Δa^x−u^z)Ls(u),L2(u)=(u2+Δ2a)L2s(u), (37)

where is a unit vector along axis and

 Ls(u)=n∑m=1em2(u−ϵm)√ϵ2m+Δ2a,em=±1. (38)

Note that when the right hand side of Eq. (38) is brought to a common denominator, the numerator is a polynomial of order . Therefore, and consequently have double zeros – the solutions of the equation . In addition, we see from Eq. (37) that there are two roots , i.e.

 Q2n(u)=(u2+Δ2a)n−1∏r=1(u−cr)2. (39)

When the spins deviate from their equilibrium positions (23), the roots of polynomial shift to note2 (). Linearizing Eq. (34) in deviations and around the stationary positions and using Eq. (39), we obtain

 δ˙ur=2i√c2r+Δ2a√(δur)2+b2r (40)

with a solution , where . Linearizing Eq. (33), one derives the spin variables in terms of . At this point we are interested only in the frequencies .

We conclude that the separation variables are indeed the normal modes of the linearized problem (since they contain a single frequency). The frequencies of small oscillations around anomalous stationary states are related to the double zeros of as . If any of has an imaginary part, the stationary state is unstable.

Next, consider linear analysis around normal eigenstates (27). In this case all spins are along axis. It follows from Eqs. (27) and (28) that , where

 Ln(u)=−1g+n∑j=1lj2(u−ϵj),lj=±1. (41)

We see that all zeros of are double zeros. There are of them as the numerator of is a polynomial of order , i.e. . As before, the stationary positions of separation variables are . Note however that there are only separation variables, so one of the zeros must remain vacant.

Next, we show that the frequencies of small oscillations around a normal stationary state are . When one of the zeros is complex, the oscillatory behavior is replaced with an exponential growth, i.e. the stationary state is unstable. Note that for small deviations from a normal state the components of the total spin are small. Therefore, in linear approximation we can set the separation variables to their equilibrium values in Eqs. (33) and (35), , i.e. only is time-dependent. As mentioned above, has a vacant zero (say ) which does not correspond to any separation variable. Eq. (35) yields

 −d(lnJ−)2idt=[n∑j=1ϵj+gJz2−n∑m=1cm]+cr=cr, (42)

where we used the fact that the contribution in square brackets vanishes. This can be seen by observing that, since are the zeros of , Eq. (41) can be written as

 Ln(u)=−1g∏m(u−cm)∏j(u−ϵj). (43)

Expanding the right hand sides of Eqs. (43) and (41) in , matching the coefficients at , and using (this follows from Eq. (27)), we see that the sum of terms in square brackets in Eq. (42) is indeed zero. It follows that and from Eq. (33) we also derive . Thus, the frequencies of oscillations around normal stationary states are .

### iii.2 Examples of root diagrams of L2(u)

Here we provide examples of root diagrams – configurations of solutions of the equation in the plane of complex , see Ref. emilts, for more examples. We saw that the zeros of evaluated in stationary states determine the frequencies of oscillations around them. Moreover, the most important features of the dynamics, e.g. the behavior of at large times, can be predicted by inspecting the root diagramemilts (), see also the discussion below Eq. (35). Similarly, we will see that the root diagram determines the number and properties of solitons corresponding to a given stationary state.

For simplicity, we assume particle-hole symmetry, i.e. the single fermion energies are symmetric with respect to zero (Fermi level). According to Eq. (23), this means

 sx(ϵm)=sx(−ϵm),sz(ϵm)=−sz(−ϵm),sy(ϵm)=−sy(−ϵm), (44)

where . These relations can also be derived from Eq. (19) using particle-hole transformation for fermion creation and annihilation operators . Note that relations (44) are preserved by equations of motion (20) and also imply .

#### iii.2.1 BCS ground state

As discussed below Eq. (24), the BCS ground state corresponds to . We note from Eq. (23) that spin components and in this state are continuous functions of single particle energy . It follows from Eqs. (37) and (38) that has two single roots at , see Fig. 6, and double roots that are the solutions of the following equation:

 Ls(u)=n∑m=112(u−ϵm)√ϵ2m+Δ20=0. (45)

All solutions are real. This can be seen by noting that changes sign between consecutive . Indeed, let be ordered so that . Since as and as , there is a point in the interval where . According to the previous subsection, this yields a frequency of small oscillations around the BCS ground state. In the continuum limit, when level spacings tend to zero, we have and . In this limit, densely fill the interval from to , i.e. has a line of double roots as shown in Fig. 6. Frequencies are also the energies of excited pairs – excitations obtained by flipping the spin from its ground state position antiparallel to the field to an equilibrium position parallel to the field , see Ref. ander1, for a discussion of this relationship between the frequencies and the excitation spectrum.

#### iii.2.2 Excited anomalous states

Next, consider two examples of anomalous stationary states obtained from the BCS ground state by flipping spins in certain energy intervals.

#### Example 1

First, let spins in the interval be flipped, i.e. . This means that the Cooper pairs in this energy interval are excitedBCS (); ander1 (). Eq. (23) implies that spin components and are discontinuous at (see the inset in Fig 7). As before has two single zeros at and double zeros that are the solutions of , where is given by Eq. (38). The difference is that in this case two of can be imaginary. Suppose and . Then, while and similarly for and we are no longer guaranteed real zeros of in intervals and as in the BCS ground state. Instead, can acquire two complex conjugate zeros. In the particle-hole symmetric cas