Comparing the full time-dependent Bogoliubov--de-Gennes equations to their linear approximation: A numerical investigation
aa r X i v : . [ m a t h - ph ] A p r EPJ manuscript No. (will be inserted by the editor)
Comparing the full time-dependent Bogoliubov–de-Gennesequations to their linear approximation: A numerical investigation
Christian Hainzl and Jonathan Seyrich Mathematisches Institut, Universit¨at T¨ubingen, Auf der Morgenstelle, 72076 T¨ubingen, Germanythe date of receipt and acceptance should be inserted later
Abstract.
In this paper we report on the results of a numerical study of the nonlinear time-dependentBardeen–Cooper–Schrieffer (BCS) equations, often also denoted as Bogoliubov–de–Gennes (BdG) equa-tions, for a one-dimensional system of fermions with contact interaction. We show that, even above thecritical temperature, the full equations and their linear approximation give rise to completely differentevolutions. In contrast to its linearization, the full nonlinear equation does not show any diffusive behaviorin the order parameter. This means that the order parameter does not follow a Ginzburg–Landau-typeof equation, in accordance with a recent theoretical result in [1]. We include a full description on thenumerical implementation of the partial differential BCS/BdG equations.
PACS.
XX.XX.XX No PACS code given
When Bardeen, Cooper and Schrieffer, shortly BCS, pub-lished one of the most famous papers in physics in 1957 [2],giving the first microscopic explanation for superconduc-tivity, a phenomenological theory for the phenomenon hadalready been around. Systems close to the critical tem-perature were described with the help of a macroscopicphase-transition parameter introduced by Ginzburg andLandau in 1950 [5]. Their theory was the first one to al-low for the description of the spatial dependence of thesuperconductivity inside superconducting alloys and thefirst with which to explain type-II superconductors andthe hexogonally shaped penetrations by magnetic flux.As the Ginzburg–Landau theory yields reliable resultson the large scale, soon the question arose as to whetherthis model can be understood as a macroscopic limit ofBCS theory for systems close to the critical temperature.Gorkov gave a positive answer to this question for the sta-tionary case shortly after the publication of BCS, see [6].A rigorous mathematical proof of the convergence wasachieved some years ago [7].But what remains unclear and controversial up to thisday, in particular in terms of a rigorous derivation, is thequestion whether the time evolution of superconductingsystems close to the critical temperature are governedby a Ginzburg–Landau type of equation. After first at-tempts for a derivation of the macroscopic limit had beenpresented [14,15,16], Gorkov and Eliasberg pointed outthat a nonlinear equation could only be valid in a gaplessregime [17]. Still, in [18,19,20] the authors made a casefor a time-dependent Ginzburg–Landau equation for su-perfluid gases at temperatures slightly above the critical one. The argument is based on the assumption that thenonlinear terms in the BCS/BdG equations only lead tosmall perturbations but do not quantitatively change thesystem’s behavior. In more detail, this would mean thatthe projection of the Cooper pair density onto the cen-ter of mass direction is governed by a nonlinear dispersiveequation. However, it has been argued recently in [1] thatfor a translation invariant homogeneous system close toequilibrium, the full BCS/BdG equations and their lin-earization do not yield the same behavior at temperaturesclose to the critical one. In particular, dissipative behaviorcan only be expected for the linear approximation of theBCS equations but not for the full equations.With our work, we demonstrate this result by means ofa thorough numerical study of the long-term evolution ofthe BCS equations and their linearization for spatially ho-mogeneous systems close to equilibrium at temperaturesslightly above the critical one. For decreasing values ofthe parameter h , defined via T = (1 + h ) T c , we evolvethe full and the linear system over a long time span andtrack the behavior of the Cooper pair density and the or-der parameter. For each values of the small parameter h ,we find clear differences between the full equation and itslinearization. Additionally, we see that the full BCS /BdGequations yield oscillations in the order parameter abouta constant value. Such a behavior has long been predictedfor and already been observed in out-of-equilibrium sys-tems, see, e.g., [35,36,37]. Although the focus of our studyis not on oscillations in particular but rather on the long-term behavior of the equations in general, it is interestingthat we can replicate such oscillations for systems close toequilibrium. Christian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length
In the realm of numerical analysis, the treatment ofquantum dynamical systems has been of huge interest formany decades (see [21] for an extensive overview). Vari-ous evolution schemes for the linear Schr¨odinger equationin varying settings have been proposed, see, e.g., [22,23,24,25,26]. Nonlinear Schr¨odinger equations such as theGross–Pitaevskii equation and equations arising from theHartree and Hartree–Fock approximation of the quantumstate have also been devoted attention to, see, e.g., [27,28,30,29] and [31,32]. Regarding the BCS regime, the sta-tionary equations have been treated numerically in [33]and the time-dependent BCS/BdG equations have beenconsidered from an analytical perspective in [34]. But,the above-mentioned studies of the out-of-equilibrium dy-namics of the BCS equations ([35,36,37]) notwithstand-ing, to the best of our knowledge the coupled nonlineartime-dependent BCS equations have not been paid muchattention to from a numerical point of view. Therefore,we come up with a reliable integration algorithm for theevolution of the system.The paper is organized as follows: First, we introducethe system we are considering and the physical backgroundin Section 2. This is followed by a brief summary of thetheoretical results of [1] in Section 3. Then, we presentour numerical results for the linear and the full equationin Section 4. Finally, we summarize our main results inSection 5. A detailed discussion of the initial setup of thesystem and the numerical implementation is provided inthe Appendix Sections A and B, respectively.
In mathematical terms, BCS theory is a special case ofa generalized Hartree–Fock variational principle, itself de-scribed by Bogoliubov-theory, for the density operators γ : H 7→ H and α : H 7→ H acting on the consideredHilbert space H . Those matrices are put together to formthe two-by-two operator-valued matrix Γ := (cid:18) γ αα − γ (cid:19) , (1)see, e.g., [38] for an introduction. The entries of the matrixcan be represented by means of their momentum distribu-tion ˆ γ ( k ) = h a † k a k i and the pair density ˆ α ( k ) = h a k a − k i ,determining the Cooper pair wave-function via Fouriertransform as α ( x − y ) = (2 π ) − / R ˆ α ( k ) e ik · ( x − y ) d k . Wesuppress spin in our notation; the pair density ˆ α is as-sumed, for simplicity, to be a spin singlet. For a one-dimensional translation invariant system of fermions attemperature T interacting via a potential V , the BCSpressure functional per unit volume is given by F T ( Γ ) = Z R ( p − µ )ˆ γ ( p )d p + Z R | α ( x ) | V ( x )d x − T S ( Γ ) , (2) where the entropy S is defined as S ( Γ ) := − Z R Tr C ( Γ ( p ) log Γ ( p )) d p. (3)The evolutions of α and γ are given by the time-dependentBCS equations which are also known as Bogoliubov–de–Gennes equations [39]. In momentum space they can bewritten conveniently in the self-consistent form ı ˙ Γ t ( p ) = [ H Γ t ( p ) , Γ t ( p )] . (4)Here, the subscript t indicates the time-dependence andthe Hamiltonian H Γ t ( p ) is defined as H Γ t ( p ) = (cid:18) p − µ
2[ ˆ V ∗ ˆ α t ]( p )2[ ˆ V ∗ ˆ α t ]( p ) µ − p (cid:19) , (5)with ∗ denoting the convolution. Calculating the upper-left and upper-right entries of the matrix-valued equa-tion (4), we arrive at the system of coupled nonlinearequations ı ˙ˆ γ t ( p ) = 2 h ( ˆ V ∗ ˆ α t )( p )ˆ α t ( p ) − ( ˆ V ∗ ˆ α t )( p )ˆ α t ( p ) i , (6) ı ˙ˆ α t ( p ) = 2( p − µ )ˆ α t ( p ) + 2( ˆ V ∗ ˆ α t )( p ) −
4( ˆ V ∗ ˆ α t )( p )ˆ γ t ( p ) . (7) In this paper we concentrate on attractive contact inter-actions, i.e., potentials of the form V ( x ) = − aδ ( x ) , a > , (8)which lead to exactly solvable systems in the stationarycase. Not only is such a potential the most interesting onefrom a physical model point-of-view but also does it allowus to implement the terms including a convolution in theequations of motion conveniently as we will illustrate inthe numerics Section B. In this work we consider initial data which, in the sta-tionary case, could be described by the Ginzburg-Landauenergy functional for temperatures T close to the criticaltemperature T c , i.e., T = T c + h for a small parameter h ∈ R . For temperatures above T c , the free energy is mini-mized by the so-called normal state Γ N for which α N = 0, γ N = / ( p − µ ) / T ) . For initial data Γ to be within therange of Ginzburg–Landau, they have to satisfy F T ( Γ ) − F T ( Γ N ) ≤ O ( h ) . (9)This condition can be complied with by choosing Γ = 11 + e H∆ / T (10) hristian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length 3 with H ∆ = (cid:18) p − µ − ∆ − ∆ µ − p (cid:19) , (11)where ∆ is a small parameter of the order of h , see, e.g.,[7]. Calculating the right-hand side of the matrix equa-tion (10) gives Γ = (cid:18) ˆ γ ˆ α ˆ α − ˆ γ (cid:19) (12)where ˆ γ and ˆ α take the special formˆ γ = 12 − p − µ (cid:18) √ ( p − µ ) + | ∆ | T (cid:19)p ( p − µ ) + | ∆ | (13)ˆ α = ∆ (cid:18) √ ( p − µ ) + | ∆ | T (cid:19)p ( p − µ ) + | ∆ | . (14)In our simulations we choose a temperature which is slightlyabove the critical temperature for the setting under con-sideration and set the initial value for the gap parameter ∆ to a non-vanishing value. We explain how to obtainthe critical temperature for our setting and how to findphysically reasonable initial values for ∆ in the AppendixSection A. For the stationary case it is well known that the Ginzburg–Landau theory emerges as the macroscopic limit of theBCS theory. To be more specific, define | α ∗ i as the trans-lation invariant minimizer of the BCS functional which, incase of the contact interaction (8), can be calculated viaˆ α ∗ ( p ) = ∆ (cid:18) √ ( p − µ ) + | ∆ | T (cid:19)p ( p − µ ) + | ∆ | . (15)Then, for the Cooper pair density | α i corresponding to thenon-translation invariant minimizer of F T , the quantity ψ := 1 h h α ∗ | α i (16)is an approximate solution of the stationary Ginzburg–Landau equation, see, e.g., [41]. This told, if there werean analogous relation between the time-dependent BCSand the GL equations, the order parameter ψ t := 1 h h α ∗ | α t i (17)should, close to T c , approximately satisfy a conventionaltime-dependent Ginzburg–Landau (TDGL) equation. Inthe spatially homogeneous case we are studying in thiswork, the conventional TDGL equation takes the form˙ ψ t = − c GL,1 ψ t − c GL,2 | ψ t | ψ t , (18) with some appropriate parameters c GL,1 and c GL,2 , see,e.g., [18] and [20, Eq. (18)]. The parameter c GL,1 dependson ( T − T c ) / ( h T c ) . Crucially, c GL,1 has the same sign as( T − T c ). Thus, the TDGL equation is dissipative for tem-peratures above T c by definition. This implies that if ψ t could be described by the TDGL for small h it should de-cay over time. However, we will demonstrate in Section 4that this is not the case, at least for the full non-linearequation. The same conclusion has been reached by ananalytical investigation recently as we will outline in Sec-tion 3. Let us decompose the particle density as γ t = γ + η t . (19)For states satisfying (9) η t appears to depend quadrati-cally on ˆ α t , see, e.g., [1, Eq. 11], and it seems legitimateto approximate the full equation by its linearization ı ˙ˆ α t ( p ) = 2( p − µ )ˆ α t ( p ) + 2( ˆ V ∗ ˆ α t )( p )(1 − γ ( p )) . (20)However, close to the Fermi-surface the quantity η t isnot small but the dominant part in the non-linear evo-lution. Consequently, the full BCS equations (6)-(7) andthe linearization (20) give rise to very different evolu-tions. Namely, Eq. (20) yields a dissipative behavior in ψ t whereas the full equations do not as is shown formallyin [1] and as we confirm by our numerical experiments be-low. Let us briefly summarize the results of [1] in the nextSection. The BCS time-evolution (4) is studied analytically in [1].Based on the work [7] the authors prove in [1, Theorem1] that | ψ t | does not vanish for any times. More precisely,it is shown in a very general setting that, if the initialstate Γ is close to the energy of the normal state, i.e., F T ( Γ ) − F T ( Γ N ) ≤ O ( h ), then the corresponding ψ t satisfies || ψ t | − | ψ || ≤ Ch / , (21)for an appropriate constant C independent of h .On the other hand, it is shown in [1] that the solution ofthe linearized equation (20) tends to 0 exponentially fastcompared to the system’s time scale of / h . In detail, usingstrategies from perturbation theory, it can be derived that | ψ t | ≈ | ψ | e t Im λ (22)holds, where λ is a resonance of order / h which emergesfrom the zero-eigenvalue at T = T c of the linear operator O = (cid:0) k − µ (cid:1) tanh − (cid:16) k − µ T (cid:17) + V . Christian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length
The combination of the bounds (21) and (22) showsclearly that the non-dissipative behavior of ψ t is a purelynon-linear effect which takes place solely in a tiny neigh-borhood of the Fermi surface.Furthermore, using the methods of [1], it is straight-forward to derive the following bound on the derivative | ˙ ψ t | = O ( / h ) . (23)In other words, although the solution | ψ t | tends to the con-stant | ψ | in the limit h →
0, its derivate might well oscil-late more and more –in line with according predictions forsystems which are suddenly perturbed out of equilibrium([35]). These findings are well reproduced in our numericalexperiments as we show now.
In this work we are interested in a qualitative study of thedifferences between the full BCS/BdG equations and theirlinearization. Thus, without loss of generality, we can workin dimensionless units and set the constant a of the contactinteraction and the chemical potential µ to a = 1 and µ = 1, respectively. The initial data for the simulationsare obtained as outlined in the Appendix Section A. Forthis, we approximate the integrals in Eqs. (25) and (26) bythe sum over the discrete momenta we take into account.For the sake of reproducibility we add the thus-obtainedvalues for T c and ∆ to the results of our simulations.For more details on the discretization of the equationsunder consideration we refer the interested reader to theAppendix Section B. h In order not to have to calculate the initial value of thegap parameter which depends on the crucial parameter h at the start of each evoluation again, we calculate ∆ with the procedure outlined in the Appendix Section Afor various h once. The interesting result is illustrated inFig. 1 where we can see that ∆ depends more or lesslinearly on the crucial parameter.Finally, with both T c and ∆ at hand, we are able topresent the results of the simulations. Doing so, we takeinto account that at temperatures T = T c + h , physicallyinteresting dynamics are expected to occur on a time-scaleof O ( / h ). Therefore, we always set t end = / h or t end = / h in the following. h = / We plot the scaled L -norm of α , which in the discretesetting is given by the sum over the K discrete momentaas 1 h k α t k = 1 h K / − X k = − K / | α Kt ( k ) | , (24) ∆ h Fig. 1.
The gap ∆ as a function of the semiclassical parameter h in semilogarithmic scale. / h || α t || t full BCSlinear BCS Fig. 2. / h k α t k as a function of integration time t for h = / .The physical parameters are T c = 0 .
19 and ∆ = 0 . as well as the modulus of the interesting macroscopic pa-rameter ψ t introduced in Eq. (17). We plot the resultsfor both the BCS equation (7) and its linear approxima-tion (20), see Fig. 2 and Fig. 3. For both quantities, thelinear equation leads to exponential decay. The full equa-tion, in contrast, coincides with the linear approximationonly for a short period after which both k α t k and | ψ t | grow again. h = / Here, too, we consider the scaled norm of the Cooper pairdensity and the modulus of the parameter ψ t . The resultsare shown in Fig 4 and Fig. 5. Again, the linear evolu-tion equation is clearly diffusive while the full equationyields a similar behavior only for very small times. Aftera short decline in the beginning of the simulation, k α t k and | ψ t | seem to oscillate. Similar oscillations have beenpredicted by [35] and observed for suddenly perturbednon-equilibrium systems in [37]. Although, as comparedto these studies, we work on systems close to equilibrium hristian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length 5 | ψ t | t full BCSlinear BCS Fig. 3. ψ t as a function of integration time t for h = / . Thephysical parameters are T c = 0 .
19 and ∆ = 0 . / h || α t || t full BCSlinear BCS Fig. 4. / h k α t k as a function of integration time t for h = / .The physical parameters are T c = 0 .
19 and ∆ = 0 . | ψ t | t full BCSlinear BCS Fig. 5. ψ t as a function of integration time t for h = / . Thephysical parameters are T c = 0 .
19 and ∆ = 0 . / h || α t || t full BCSlinear BCS Fig. 6. / h k α t k as a function of integration time t for h = / . The physical parameters are T c = 0 .
19 and ∆ = 0 . | ψ t | t full BCSlinear BCS Fig. 7. ψ t as a function of integration time t for h = / . Thephysical parameters are T c = 0 .
19 and ∆ = 0 . and slightly above the critical temperature, it is inter-esting to see that our long-term evolutions show oscill-lations which resemble the ones predicted for the out-of-equilibrium case. h = / Once more, we depict the time evolution of k α t k and | ψ t | ,cf. Figs. 6 and 7. The conclusions we can draw from thesetwo plots are the same as for h = / , the only differencebeing the faster oscilllations in line with the bound (23).Most importantly, even for this small value of h , we onlyobserve diffusion for the linear approximation which, be-lying its name, does not approximate the BCS equationfor reasonably long time intervals. Let us summarize ourresults in the concluding Section. We have introduced a reliable integration scheme for thetime-dependent BCS equation and its linear approxima-
Christian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length tion in spatially homogeneous settings. With the help ofthese algorithms, we could perform numerical long-termstudies for systems close to equilibrium in order to inves-tigate the time-evolution of the order parameter at thelimit close to the critical temperature. The study showsvery clearly that, opposed to the linear case, the full BCSequation does not yield any decay over time in the or-der parameter ψ t . Since the conventional time dependentGinzburg–Landau equation is dissipative above the criti-cal temperature by definition, it cannot give a valid macro-scopic limit of the full time-dependent BCS/BdG equa-tions. It can only be seen as the limit of the linearizationof the full equations but the effects of this linearizationcould clearly be shown not to be negligible in the consid-ered regime. We thus confirm the analysis provided in [1].In addition, when evolving the system as described bythe non-linear BCS/BdG equations, we observed oscilla-tions in the Cooper pair density and in the order param-eter about a finite value which are similar to oscillationswhich have been observed for out-of-equilibrium systemsin various works. Acknowledgements:
We want to thank Christian Lu-bich for useful remarks and discussions. This work waspartially funded by the DFG grant GRK 1838.
Declaration of contribution
Contribution to the orig-inal manuscript was: C. H. 30%, J. S. 70%.Contribution to the update of the manuscript was: C. H.50%, J. S. 50%.
A Criticial temperature and initial energy gap
For translation invariant systems with contact interaction,the cricital temperature T c is well-known to be given im-plicitly by 2 πa = Z R tanh (cid:16) p − µ T c (cid:17) p − µ d p, (25)see, e.g. [8,9,20]. The energy gap ∆ between the super-conducting state and the normal state at temperaturesbeneath the cricital temperature, in turn, can be obtainedfrom the relation2 πa = Z R tanh (cid:18) √ ( p − µ ) + | ∆ | T (cid:19)p ( p − µ ) + | ∆ | d p. (26)In order to calculate the critical temperature and a real-istic initial value for the gap parameter we thus proceedas follows: For a given value of the small crucial parame-ter h , we first determine the critical temperature T c andset T = T c − h . For this temperature we then search thecorresponding gap ∆ following the above definition (26)and set ∆ = ∆ as its initial state. Finally, as we areinterested in simulations for temperatures slightly abovethe critical temperature, we put T = T c + h and insert this into Eq. (12) together with the just-determined ∆ .This yields physically realistic conditions which satisfy theenergy constraint (9). B Numerical treatment of the equations
We want to model a system of infinite spacial extension,which, of course, is not possible to achieve on a machine.Therefore, we pretend our system to be periodic in spacebut with a large enough period.
B.1 Finite extension and discrete system
In the Ginzburg-Landau regime, one often takes into ac-count external potentials that vary on a scale of O ( / h )and, consequently, lead to variations of the system whichoccur over intervals of that very scale. Thus, a valid modelsystem should have an extension no smaller than thosephysical variations. But, in order to avoid artificial effectsdue to the periodicity, it is necessary to enlarge this exten-sions by some multiples of / h . For convenience we further-more include a factor of 2 π , wherefore we consider systemswith period πN / h , 1 < N ∈ N . The kernels of the den-sity operators are now functions on L ([0 , πN / h ] R ). Inorder to simplify the notation, we introduce macroscopicvariables via x mac := h / N x . We end up in a 2 π -periodicsetting for which the inner product of two functions f and g is just h f | g i = X k ∈ Z ˆ f ˆ g. (27)The self-consistent BCS equations are now given by ı ˙ Γ t ( k ) = [ H Γ t ( K ) , Γ t ( k )] , k ∈ Z , (28)with the Hamiltonian H Γ t = (cid:16) h N k − µ (cid:17)
2[ ˆ V Nh ∗ ˆ α t ] k
2[ ˆ V Nh ∗ ˆ α t ] k (cid:16) µ − h N k (cid:17) (29)and the Fourier transform ˆ V Nh of V Nh ( · ) := V ( N / h · ).Please note that in the present discrete case the convo-lution of two summable series a k and b k has to be under-stood as ( a ∗ b ) k = X j ∈ Z a k − j b j . (30) B.2 The equations for a delta potential
For systems on a large torus with a contact interaction (8),we can easily see thatˆ V Nh ∗ ˆ α t = − a h φ | α i , (31) hristian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length 7 where φ is the state given by φ ( k ) = 1 for all integers k . With this, the equations of motion take the convenientform ı ˙ˆ γ t ( k ) = 2 a h h φ | α i ˆ α t ( k ) − h φ | α i ˆ α t ( k ) i , (32) ı ˙ˆ α t ( k ) = 2 (cid:18) h N k − µ (cid:19) ˆ α t ( k ) + 2 a h φ | α i (2ˆ γ t ( k ) − ı ˙ˆ α t ( k ) = 2 (cid:18) h N k − µ (cid:19) ˆ α t ( k ) + 2 a h φ | α i (2ˆ γ ( k ) − B.3 Space discretization
As the BCS equations are given in their momentum spacerepresentation, it is most convenient to use the so-called
Fourier collocation . This means that a 2 π -periodic func-tion f ( x ) = P j ∈ Z ˆ f ( j ) e ıkx is approximated by f K ( x ) = K − X j = − K ˆ f K ( j ) e ıkx , (35)where the coefficients ˆ f K ( j ) are obtained by the discreteFourier transform of the values f j = f ( π / K j ), j = − K / , ..., K / −
1. Mathematically speaking we work on the subspace spannedby the first K eigenfunctions of the Laplacian on [0 , π ].As a consequence, the evolution of the system is given bythe K -dimensional system of ordinary differential equa-tions (ODE) ı ˙ˆ γ Kt ( k ) = 2 a h h φ K | α K i ˆ α Kt ( k ) − (cid:10) φ K | α K (cid:11) ˆ α Kt ( k ) i , (36) − K ≤ k ≤ K − ,ı ˙ˆ α Kt ( k ) = 2 (cid:18) h N k − µ (cid:19) ˆ α Kt ( k )+ 2 a (cid:10) φ K | α K (cid:11) (cid:0) γ Kt ( k ) − (cid:1) , (37) − K ≤ k ≤ K − , and accordingly for the linear case. From numerical anal-ysis it is well known that (35) yields a very good approx-imation to 2 π -periodic functions with the discretizationerror decreasing rapidly as a function of K , see, e.g. [21],Chapter III.1.3.For practical reasons we set K to be an integer powerof 2 so that for a given ˆ α Kt ( j ), j = − K / , ..., K / −
1, the corresponding distribution α Kt ( x ) at the discrete points x j = πK j can be computed efficiently with the well-knownfast Fourier transform (FFT). As we want to resolve phe-nomena happening on the microscopic scale O ( h / N ), wechoose K = M Nh (38)for a large enough integer M . Let us now explain how wesolve the system of ODE (36)–(37). B.4 Solving the system of ordinary differentialequations
We first notice that the Hamiltonian H Γ t is self-adjoint.Thus, the time-evolution of Γ t is a unitary transformationand, hence, its eigenvalues are preserved. With regard todefinition (1), the eigenvalues can be readily computed as λ , ( p ) = 12 ± s(cid:18) ˆ γ t ( p ) − (cid:19) + | ˆ α t ( p ) | , (39)and we see that the equality (cid:18) ˆ γ t ( p ) − (cid:19) + | ˆ α t ( p ) | = (cid:18) ˆ γ ( p ) − (cid:19) + | ˆ α ( p ) | (40)holds. Solving this for ˆ γ t , we getˆ γ t ( p ) = ( + p h ( p ) − | ˆ α t ( p ) | for p < µ , − p h ( p ) − | ˆ α t ( p ) | for p ≥ µ , (41)where we have defined the auxiliary function h ( p ) := (cid:18) ˆ γ ( p ) − (cid:19) + | ˆ α ( p ) | . (42)The signs in Eq. (41) can be inferred from the initial valueswe use in this work, cf. (12). They are such that ˆ γ ( p ) isgreater than / for µ > p and less than or equal to / for µ ≤ p .Inserting the discrete analogon of Eq. (41) into therelevant equation of motion (37), we get the nonlinearcoupled system of equations ı ˙ˆ α Kt ( k ) = 2 (cid:18) h N k − µ (cid:19) ˆ α Kt ( k ) ± a (cid:10) φ K | α K (cid:11) q h ( k ) − | ˆ α Kt ( k ) | , − K ≤ k ≤ K − . (43)This said, we now present our time integration algorithm. Christian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length
B.5 Time discretization
Putting it in a formal way, the system we have to integrateis given by ( d y ( t )d t = f ( y ( t )) , y (0) = y , (44)with y = (cid:0) ˆ α K ( − K / ) . . . ˆ α K ( K / − (cid:1) T ∈ R K . (45)The right hand side of our initial problem can be writtenas the sum of two terms, f ( y ) = f ( y ) + f ( y ) , (46)where f represents the linear part which resembles thekinetic part in the Schr¨odinger-equation and f is the non-linear part. Let τ denote a time step and Φ τ,f the smoothmap between y (0) and y ( τ ). Given the special form (46)of the differential equation, one can approximate Φ τ,f nu-merically by Φ num τ,f ( y ) = (cid:0) Φ τ / ,f ◦ Φ τ,f ◦ Φ τ / ,f (cid:1) ( y ) . (47)This is the well-known Strang splitting . Applying it suc-cessively yields an approximation to the exact solution attimes t = nτ , n = 1 , , ... , the error of which decreasesquadratically as a function of the step size τ , see, e.g. [42],Chapter II.5.The advantage of the Strang splitting is that Φ τ,h canbe calculated exactly as Φ τ,f ( · ) = e − ı (cid:16) h N k − µ (cid:17) τ · . (48)As for Φ τ,f , it has to be approximated due to the nonlin-earity. For this, we choose a simple Runge-Kutta schemeas proposed by [43] whose numerical error is small com-pared to the error expacted from the splitting [46]. Beforestarting the simulations, we still need to fix the mentioneddiscretization parameters τ and K . In our case, K itselfdepends on three parameters, cf. Eq. (38). As h is thesemiclassical parameter we want to vary throughout thestudy, we have to choose reasonable values for the remain-ing quantities M , N and τ . We first consider τ . B.6 Fixing the time discretization parameter
The step size has to be chosen small enough for both thenumerical approximation of Φ τ,f and the Strang splittingto give accurate results. For our simulations it turned outthat reliable results can only be expacted for a step size in-versely proportional to K . Playing safe we include a smallfactor and set τ = . / K . As a measure for the time inte-grator’s accuracy, we consider the discrete analogon of thefree energy introduced in Eq. (2) above, which is given by F KT ( Γ K ) = K / − X k = − K / (cid:18) h N k − µ (cid:19) ˆ γ K ( k )+ 12 π Z π V ( x mac ) | α K ( x mac ) | d x mac − T S ( Γ K ) . (49) ∆ F T ( t ) t Fig. 8.
Relative error ∆F T of the discretized free energyagainst integration time t in semilogarithmic scale for h = / .The physical parameters are T c = 0 .
19 and ∆ = 0 . A short calculation yields that this quantity is conservedunder the exact flow of the corresponding initial valueproblem. Therefore, the reliability of a numerical integra-tion scheme can be checked by tracking the relative error ∆F T , defined by ∆F T ( t ) = (cid:12)(cid:12)(cid:12)(cid:12) F T ( Γ Kt ) − F T ( Γ K ) F T ( Γ K ) (cid:12)(cid:12)(cid:12)(cid:12) , (50)along the numerical evolution. Recurring to a constantof motion as a criterion of accurateness is a much appliedprocedure in various computational fields, see, e.g. [44,45].Following this line of reasoning, we have verified the accu-racy of our time integrator for every simulation presentedbelow. As an example, we show the plot of ∆F T corre-sponding to the simulations of Subsection 4.3 in Fig. 8. B.7 Fixing the space discretization parameters
We have seen in the previous Subsection that the timestep has to be inversely proportional to the dimension ofthe subspace we are approximating our system on. Fur-thermore, every time step requires a computational ef-fort which grows linearly with K . Consequently, the com-plete CPU time for a simulation over a given time inter-val [0 , t end ] is quadratic in K . So the dimension of thesubspace and, thus, the related N and M should be thesmallest possible. In order to check how small a M wecan choose without any significant loss of accuracy, we fix N = 8 and h = / and calculate the cricital temperaturevia the discretized version of Eq. (25),2 πa = K / − X k = − K / tanh (cid:18) h N k − µ T c (cid:19) h N k − µ , (51)for different values of M . The result can be seen in Fig. 9.For different values of N and h we get the same plot. hristian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length 9 T c M Fig. 9.
Critical temperature T c as a function of the number ofmomenta per unit volume M for N = 8 and h = / . | ψ t | t Fig. 10. ψ t as a function of integration time t for N = 4 and h = / . The physical parameters are T c = 0 .
19 and ∆ = 0 . We see that for M = 256 the critical temperature isstill slightly too small. However, when comparing the evo-lutions obtained with M = 256 to the according onesfor M = 512, the relevant figures are indistinguighablefrom each another. For the sake of efficiency, we thus fix M = 256 for the rest of this work.As for the extension of our interval, N , we have tochoose it large enough so that the solution cannot reachthe boundaries during the simulation. As, by construction,we work with a periodic setting, a solution reaching oneend of the interval would enter again at the other end, thusleading to unphysical interference. As an example of thisnumerical artifact, we consider the case h = / , N = 4and plot the modulus of ψ t in Fig. 10. We observe oscil-lations for larger t which should not show up in reality,cf. Subsection 4.3. Whenever we encountered such an ar-tifact, we successively increased N by factors of 2 untilthe artifact vanished. References
1. R. L. Frank, C. Hainzl, B. Schlein, R. Seiringer, preprint,arXiv:1504.058852. J. Bardeen, L. N. Cooper and J. R. Schrieffer,
Physical Re-view , 1175 (1957)3. E. P. Gross,
Il Nuovo Cimento , 454 (1961)4. L. P. Pitaevskii, Sov. Phys. JETP , 451 (1961)5. L. D. Landau and V. L. Ginzburg, Zh. Eksp. Teor. Fiz. ,1064 (1950)6. L. P. Gorkov, Sov. Phys. JETP , 1364 (1959)7. R. L. Frank, C. Hainzl, R. Seiringer and J. P. Solovej, J.Amer. Math. Soc. , 667 (2012)8. A. J. Leggett, Lecture Notes in Physics, vol. 115 (Springer,Berlin, 1980) pp. 13-279. P. Nozi´eres and S. Schmitt-Rink.
Journal of Low Tempera-ture Physics , 59:195–211, 1985.10. M. Drechsler and W. Zwerger,
Ann. Phys. , 15 (1992).11. C. Hainzl and R. Seiringer,
Lett. Math. Phys. , 119(2012)12. C. Hainzl and B. Schlein,
Journal of Functional Analysis , 399 (2013)13. P. Pieri and G.C. Strinati,
Phys. Rev. Lett. , 030401(2003)14. M. J. Stephen and H. Suhl, Phys. Rev. Lett. , 797 (1964)15. A. Schmid, Z. Phys. , 336 (1966)16. E. Abrahams and T. Tsuneto,
Phys. Rev. , 416 (1966)17. L. P. Gorkov and G. M. Eliasberg,
Sov. Phys. JETP ,328 (1968)18. M. Cyrot, Rep. Prog. Phys. , 103 (1973)19. C. A. R. S´a de Melo, M. Randeria and J. R. Engelbrecht, Phys. Rev. Lett. , 3202 (1993)20. M. Randeria, Bose-Einstein Condensation , eds. A. Griffin,D. W. Snoke and S. Stringari (Cambridge University Press,July 1996), pp. 355–392.21. C. Lubich,
From Quantum to Classical Molecular Dynam-ics: Reduced Models and Numerical Analysis (Europ. Math.Soc., Z¨urich, 2008)22. M. D. Feit, J. A. Fleck Jr. and A. Steiger,
J. Comput.Phys. , 412 (1982)23. S. K. Gray and D. E. Manolopoulos, J. Comput. Phys. , 7099 (1996)24. S. Blanes, F. Casas and A. Murua,
J. Comput. Phys. ,234105 (2006)25. H. Tal-Ezer and R. Kosloff,
J. Chem. Phys. , 3967 (1984)26. T. J. Park, J. Chem. Phys. , 5870 (1986)27. W. Bao, D. Jaksch and P. A. Markowich, J. Comput. Phys. , 318 (2003)28. M. Caliari, C. Neuhauser and M. Thalhammer,
J. Comput.Phys. , 822 (2009)29. L. Gauckler and C. Lubich,
Found. Comput. Math. , 275(2010)30. Y.-F. Tang, L. V´azquez, F. Zhang and V. M. P´erez-Garc´ıa, Computers Math. Applic. , 73 (1996)31. C. Lubich, Appl. Numer. Math. , 355 (2004)32. C. Lubich, Math. Comp. , 765 (2005)33. M. L´ewin and S. Paul, ESAIM: Mathematical Modellingand Numerical Analysis , 53 (2014)34. C. Hainzl, E. Lenzmann, M. L´ewin and B. Schlein, Ann.H. Poincar´e , 1023 (2010)35. A. . Volkov and S. M. Kogan, Sov. Phys. JETP , 1018(1974)0 Christian Hainzl, Jonathan Seyrich: Title Suppressed Due to Excessive Length36. R. A. Barankov, L. S. Levitov and B. Z. Spivak, Phys. Rev.Lett. , 160401 (2004)37. E. A. Yuzbashyan, O. Tsyplyatyev and B. L. Altshuler, Phys. Rev. Lett. , 097005 (2006)38. V. Bach, E. H. Lieb and J. P. Solovej, Journal of statisticalphysics , 3 (1994)39. P. G. de Gennes and P. A. Pincus Superconductivity ofmetals and alloys (WA Benjamin, New York, 1966), Vol. 86.40. C. Hainzl, E. Hamza, R. Seiringer and J. P. Solovej,
Comm.Math. Phys , 349 (2008)41. R.L. Frank, C. Hainzl, R. Seiringer, J.P. Solovej,
Opera-tor Theory: Advances and Applications
Volume , 57–88(2013),42. E. Hairer, C. Lubich and G. Wanner,
Geometric numeri-cal integration. Structure-preserving algorithms for ordinarydifferential equations (Springer, Berlin, 2006), 2nd ed.43. W. Press, S. Teukolsky, W. Vetterling and B. Flannery,
Numerical Recipes in C. The art of scientific computing (Cambridge University Press, 1992), 2nd ed.44. J. Seyrich,
Phys. Rev. D , 084064 (2013)45. G. Lukes-Gerakopoulos, Phys. Rev. D , 043002 (2014)46. Note that in the linear case Φ τ,f2