Comparing the full time-dependent Bogoliubov–de-Gennes equations to their linear approximation: A numerical investigation
In this paper we report on the results of a numerical study of the nonlinear time-dependent Bardeen–Cooper–Schrieffer (BCS) equations, often also denoted as Bogoliubov–de–Gennes (BdG) equations, for a one-dimensional system of fermions with contact interaction. We show that, even above the critical temperature, the full equations and their linear approximation give rise to completely different evolutions. In contrast to its linearization, the full nonlinear equation does not show any diffusive behavior in the order parameter. This means that the order parameter does not follow a Ginzburg–Landau-type of equation, in accordance with a recent theoretical result in FHSchS (). We include a full description on the numerical implementation of the partial differential BCS/BdG equations.
When Bardeen, Cooper and Schrieffer, shortly BCS, published one of the most famous papers in physics in 1957 BCS (), giving the first microscopic explanation for superconductivity, a phenomenological theory for the phenomenon had already been around. Systems close to the critical temperature were described with the help of a macroscopic phase-transition parameter introduced by Ginzburg and Landau in 1950 Ginzburg-Landau (). Their theory was the first one to allow for the description of the spatial dependence of the superconductivity inside superconducting alloys and the first with which to explain type-II superconductors and the hexogonally shaped penetrations by magnetic flux.
As the Ginzburg–Landau theory yields reliable results on the large scale, soon the question arose as to whether this model can be understood as a macroscopic limit of BCS theory for systems close to the critical temperature. Gorkov gave a positive answer to this question for the stationary case shortly after the publication of BCS, see Gorkov (). A rigorous mathematical proof of the convergence was achieved some years ago proveGLst ().
But what remains unclear and controversial up to this day, in particular in terms of a rigorous derivation, is the question whether the time evolution of superconducting systems close to the critical temperature are governed by a Ginzburg–Landau type of equation. After first attempts for a derivation of the macroscopic limit had been presented SteSu (); Schmid (); AbrTs (), Gorkov and Eliasberg pointed out that a nonlinear equation could only be valid in a gapless regime GorEli (). Still, in Cyrot (); Randeriaetal (); Randeria () the authors made a case for a time-dependent Ginzburg–Landau equation for superfluid gases at temperatures slightly above the critical one. The argument is based on the assumption that the nonlinear terms in the BCS/BdG equations only lead to small perturbations but do not quantitatively change the system’s behavior. In more detail, this would mean that the projection of the Cooper pair density onto the center of mass direction is governed by a nonlinear dispersive equation. However, it has been argued recently in FHSchS () that for a translation invariant homogeneous system close to equilibrium, the full BCS/BdG equations and their linearization do not yield the same behavior at temperatures close to the critical one. In particular, dissipative behavior can only be expected for the linear approximation of the BCS equations but not for the full equations.
With our work, we demonstrate this result by means of a thorough numerical study of the long-term evolution of the BCS equations and their linearization for spatially homogeneous systems close to equilibrium at temperatures slightly above the critical one. For decreasing values of the parameter , defined via , we evolve the full and the linear system over a long time span and track the behavior of the Cooper pair density and the order parameter. For each values of the small parameter , we find clear differences between the full equation and its linearization. Additionally, we see that the full BCS /BdG equations yield oscillations in the order parameter about a constant value. Such a behavior has long been predicted for and already been observed in out-of-equilibrium systems, see, e.g., VolkovKogan (); Barankov1 (); Yuzbashyan (). Although the focus of our study is not on oscillations in particular but rather on the long-term behavior of the equations in general, it is interesting that we can replicate such oscillations for systems close to equilibrium.
In the realm of numerical analysis, the treatment of quantum dynamical systems has been of huge interest for many decades (see bluebook () for an extensive overview). Various evolution schemes for the linear Schrödinger equation in varying settings have been proposed, see, e.g., FeitFleckSteiger (); GrayManolopoulos (); BlanesCasasMurua (); TalEzerKosloff (); ParkLight (). Nonlinear Schrödinger equations such as the Gross–Pitaevskii equation and equations arising from the Hartree and Hartree–Fock approximation of the quantum state have also been devoted attention to, see, e.g., Bao (); Thalhammer (); TVZPG (); GaucklerLubich () and LubichMCTDH (); Lubichvar (). Regarding the BCS regime, the stationary equations have been treated numerically in LewinPaul () and the time-dependent BCS/BdG equations have been considered from an analytical perspective in HLLSch (). But, the above-mentioned studies of the out-of-equilibrium dynamics of the BCS equations (VolkovKogan (); Barankov1 (); Yuzbashyan ()) notwithstanding, to the best of our knowledge the coupled nonlinear time-dependent BCS equations have not been paid much attention to from a numerical point of view. Therefore, we come up with a reliable integration algorithm for the evolution of the system.
The paper is organized as follows: First, we introduce the system we are considering and the physical background in Section 2. This is followed by a brief summary of the theoretical results of FHSchS () in Section 3. Then, we present our numerical results for the linear and the full equation in Section 4. Finally, we summarize our main results in Section 5. A detailed discussion of the initial setup of the system and the numerical implementation is provided in the Appendix Sections A and B, respectively.
2 Physical and mathematical background
2.1 Energy functional and BCS equations
In mathematical terms, BCS theory is a special case of a generalized Hartree–Fock variational principle, itself described by Bogoliubov-theory, for the density operators and acting on the considered Hilbert space . Those matrices are put together to form the two-by-two operator-valued matrix
see, e.g., BachLiebSolovej () for an introduction. The entries of the matrix can be represented by means of their momentum distribution and the pair density , determining the Cooper pair wave-function via Fourier transform as . We suppress spin in our notation; the pair density is assumed, for simplicity, to be a spin singlet. For a one-dimensional translation invariant system of fermions at temperature interacting via a potential , the BCS pressure functional per unit volume is given by
where the entropy is defined as
The evolutions of and are given by the time-dependent BCS equations which are also known as Bogoliubov–de–Gennes equations deGennes (). In momentum space they can be written conveniently in the self-consistent form
Here, the subscript indicates the time-dependence and the Hamiltonian is defined as
with denoting the convolution. Calculating the upper-left and upper-right entries of the matrix-valued equation (4), we arrive at the system of coupled nonlinear equations
2.2 Contact interactions
In this paper we concentrate on attractive contact interactions, i.e., potentials of the form
which lead to exactly solvable systems in the stationary case. Not only is such a potential the most interesting one from a physical model point-of-view but also does it allow us to implement the terms including a convolution in the equations of motion conveniently as we will illustrate in the numerics Section B.
2.3 Initial values
In this work we consider initial data which, in the stationary case, could be described by the Ginzburg-Landau energy functional for temperatures close to the critical temperature , i.e., for a small parameter . For temperatures above , the free energy is minimized by the so-called normal state for which , . For initial data to be within the range of Ginzburg–Landau, they have to satisfy
This condition can be complied with by choosing
where and take the special form
In our simulations we choose a temperature which is slightly above the critical temperature for the setting under consideration and set the initial value for the gap parameter to a non-vanishing value. We explain how to obtain the critical temperature for our setting and how to find physically reasonable initial values for in the Appendix Section A.
2.4 Ginzburg–Landau and macroscopic parameter
For the stationary case it is well known that the Ginzburg–Landau theory emerges as the macroscopic limit of the BCS theory. To be more specific, define as the translation invariant minimizer of the BCS functional which, in case of the contact interaction (8), can be calculated via
Then, for the Cooper pair density corresponding to the non-translation invariant minimizer of , the quantity
is an approximate solution of the stationary Ginzburg–Landau equation, see, e.g., FHSS3 (). This told, if there were an analogous relation between the time-dependent BCS and the GL equations, the order parameter
should, close to , approximately satisfy a conventional time-dependent Ginzburg–Landau (TDGL) equation. In the spatially homogeneous case we are studying in this work, the conventional TDGL equation takes the form
with some appropriate parameters and , see, e.g., Cyrot () and (Randeria, , Eq. (18)). The parameter depends on . Crucially, has the same sign as . Thus, the TDGL equation is dissipative for temperatures above by definition. This implies that if could be described by the TDGL for small it should decay over time. However, we will demonstrate in Section 4 that this is not the case, at least for the full non-linear equation. The same conclusion has been reached by an analytical investigation recently as we will outline in Section 3.
2.5 The linear approximation
Let us decompose the particle density as
However, close to the Fermi-surface the quantity is not small but the dominant part in the non-linear evolution. Consequently, the full BCS equations (6)-(7) and the linearization (20) give rise to very different evolutions. Namely, Eq. (20) yields a dissipative behavior in whereas the full equations do not as is shown formally in FHSchS () and as we confirm by our numerical experiments below. Let us briefly summarize the results of FHSchS () in the next Section.
3 Recent mathematical results
The BCS time-evolution (4) is studied analytically in FHSchS (). Based on the work proveGLst () the authors prove in (FHSchS, , Theorem 1) that does not vanish for any times. More precisely, it is shown in a very general setting that, if the initial state is close to the energy of the normal state, i.e., , then the corresponding satisfies
for an appropriate constant independent of .
On the other hand, it is shown in FHSchS () that the solution of the linearized equation (20) tends to exponentially fast compared to the system’s time scale of . In detail, using strategies from perturbation theory, it can be derived that
holds, where is a resonance of order which emerges from the zero-eigenvalue at of the linear operator .
Furthermore, using the methods of FHSchS (), it is straightforward to derive the following bound on the derivative
In other words, although the solution tends to the constant in the limit , its derivate might well oscillate more and more –in line with according predictions for systems which are suddenly perturbed out of equilibrium (VolkovKogan ()). These findings are well reproduced in our numerical experiments as we show now.
In this work we are interested in a qualitative study of the differences between the full BCS/BdG equations and their linearization. Thus, without loss of generality, we can work in dimensionless units and set the constant of the contact interaction and the chemical potential to and , respectively. The initial data for the simulations are obtained as outlined in the Appendix Section A. For this, we approximate the integrals in Eqs. (25) and (26) by the sum over the discrete momenta we take into account. For the sake of reproducibility we add the thus-obtained values for and to the results of our simulations. For more details on the discretization of the equations under consideration we refer the interested reader to the Appendix Section B.
4.1 Gap as a function of
In order not to have to calculate the initial value of the gap parameter which depends on the crucial parameter at the start of each evoluation again, we calculate with the procedure outlined in the Appendix Section A for various once. The interesting result is illustrated in Fig. 1 where we can see that depends more or less linearly on the crucial parameter.
Finally, with both and at hand, we are able to present the results of the simulations. Doing so, we take into account that at temperatures , physically interesting dynamics are expected to occur on a time-scale of . Therefore, we always set or in the following.
4.2 Results for
We plot the scaled -norm of , which in the discrete setting is given by the sum over the discrete momenta as
For both quantities, the linear equation leads to exponential decay. The full equation, in contrast, coincides with the linear approximation only for a short period after which both and grow again.
4.3 Results for
Again, the linear evolution equation is clearly diffusive while the full equation yields a similar behavior only for very small times. After a short decline in the beginning of the simulation, and seem to oscillate. Similar oscillations have been predicted by VolkovKogan () and observed for suddenly perturbed non-equilibrium systems in Yuzbashyan (). Although, as compared to these studies, we work on systems close to equilibrium and slightly above the critical temperature, it is interesting to see that our long-term evolutions show oscilllations which resemble the ones predicted for the out-of-equilibrium case.
4.4 Results for
The conclusions we can draw from these two plots are the same as for , the only difference being the faster oscilllations in line with the bound (23). Most importantly, even for this small value of , we only observe diffusion for the linear approximation which, belying its name, does not approximate the BCS equation for reasonably long time intervals. Let us summarize our results in the concluding Section.
We have introduced a reliable integration scheme for the time-dependent BCS equation and its linear approximation in spatially homogeneous settings. With the help of these algorithms, we could perform numerical long-term studies for systems close to equilibrium in order to investigate the time-evolution of the order parameter at the limit close to the critical temperature. The study shows very clearly that, opposed to the linear case, the full BCS equation does not yield any decay over time in the order parameter . Since the conventional time dependent Ginzburg–Landau equation is dissipative above the critical temperature by definition, it cannot give a valid macroscopic limit of the full time-dependent BCS/BdG equations. It can only be seen as the limit of the linearization of the full equations but the effects of this linearization could clearly be shown not to be negligible in the considered regime. We thus confirm the analysis provided in FHSchS ().
In addition, when evolving the system as described by the non-linear
BCS/BdG equations, we observed oscillations in the Cooper pair density and in the order parameter about a finite value which are similar
to oscillations which have been observed for out-of-equilibrium systems in various works.
Acknowledgements: We want to thank Christian Lubich for useful remarks and discussions. This work was partially funded by the DFG grant GRK 1838.
Declaration of contribution Contribution to the original manuscript was: C. H. 30%, J. S. 70%.
Contribution to the update of the manuscript was: C. H. 50%, J. S. 50%.
Appendix A Criticial temperature and initial energy gap
For translation invariant systems with contact interaction, the cricital temperature is well-known to be given implicitly by
see, e.g. Leggett (); NRS (); Randeria (). The energy gap between the superconducting state and the normal state at temperatures beneath the cricital temperature, in turn, can be obtained from the relation
In order to calculate the critical temperature and a realistic initial value for the gap parameter we thus proceed as follows: For a given value of the small crucial parameter , we first determine the critical temperature and set . For this temperature we then search the corresponding gap following the above definition (26) and set as its initial state. Finally, as we are interested in simulations for temperatures slightly above the critical temperature, we put and insert this into Eq. (12) together with the just-determined . This yields physically realistic conditions which satisfy the energy constraint (9).
Appendix 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 space but with a large enough period.
b.1 Finite extension and discrete system
In the Ginzburg-Landau regime, one often takes into account external potentials that vary on a scale of and, consequently, lead to variations of the system which occur over intervals of that very scale. Thus, a valid model system should have an extension no smaller than those physical variations. But, in order to avoid artificial effects due to the periodicity, it is necessary to enlarge this extensions by some multiples of . For convenience we furthermore include a factor of , wherefore we consider systems with period , . The kernels of the density operators are now functions on . In order to simplify the notation, we introduce macroscopic variables via . We end up in a -periodic setting for which the inner product of two functions and is just
The self-consistent BCS equations are now given by
with the Hamiltonian
and the Fourier transform of . Please note that in the present discrete case the convolution of two summable series and has to be understood as
b.2 The equations for a delta potential
For systems on a large torus with a contact interaction (8), we can easily see that
where is the state given by for all integers . With this, the equations of motion take the convenient form
for the nonlinear case and
for the linear case.
Up to now we are still left with an infinite-dimensional system of equations. In order to solve these numerically, we have to introduce a suitable finite-dimensional subspace.
b.3 Space discretization
As the BCS equations are given in their momentum space representation, it is most convenient to use the so-called Fourier collocation. This means that a -periodic function is approximated by
where the coefficients are obtained by the discrete Fourier transform of the values , . Mathematically speaking we work on the subspace spanned by the first eigenfunctions of the Laplacian on . As a consequence, the evolution of the system is given by the -dimensional system of ordinary differential equations (ODE)
and accordingly for the linear case. From numerical analysis it is well known that (35) yields a very good approximation to -periodic functions with the discretization error decreasing rapidly as a function of , see, e.g. bluebook (), Chapter III.1.3.
For practical reasons we set to be an integer power of so that for a given , , the corresponding distribution at the discrete points can be computed efficiently with the well-known fast Fourier transform (FFT). As we want to resolve phenomena happening on the microscopic scale , we choose
b.4 Solving the system of ordinary differential equations
We first notice that the Hamiltonian is self-adjoint. Thus, the time-evolution of is a unitary transformation and, hence, its eigenvalues are preserved. With regard to definition (1), the eigenvalues can be readily computed as
and we see that the equality
holds. Solving this for , we get
where we have defined the auxiliary function
b.5 Time discretization
Putting it in a formal way, the system we have to integrate is given by
The right hand side of our initial problem can be written as the sum of two terms,
where represents the linear part which resembles the kinetic part in the Schrödinger-equation and is the nonlinear part. Let denote a time step and the smooth map between and . Given the special form (46) of the differential equation, one can approximate numerically by
This is the well-known Strang splitting. Applying it successively yields an approximation to the exact solution at times , , the error of which decreases quadratically as a function of the step size , see, e.g. hairerlubichwanner (), Chapter II.5.
The advantage of the Strang splitting is that can be calculated exactly as
As for , it has to be approximated due to the nonlinearity. For this, we choose a simple Runge-Kutta scheme as proposed by nr () whose numerical error is small compared to the error expacted from the splitting remark (). Before starting the simulations, we still need to fix the mentioned discretization parameters and . In our case, itself depends on three parameters, cf. Eq. (38). As is the semiclassical parameter we want to vary throughout the study, we have to choose reasonable values for the remaining quantities , and . We first consider .
b.6 Fixing the time discretization parameter
The step size has to be chosen small enough for both the numerical approximation of and the Strang splitting to give accurate results. For our simulations it turned out that reliable results can only be expacted for a step size inversely proportional to . Playing safe we include a small factor and set . As a measure for the time integrator’s accuracy, we consider the discrete analogon of the free energy introduced in Eq. (2) above, which is given by
A short calculation yields that this quantity is conserved under the exact flow of the corresponding initial value problem. Therefore, the reliability of a numerical integration scheme can be checked by tracking the relative error , defined by
along the numerical evolution. Recurring to a constant of motion as a criterion of accurateness is a much applied procedure in various computational fields, see, e.g. Seyrich (); Lukes (). Following this line of reasoning, we have verified the accuracy of our time integrator for every simulation presented below. As an example, we show the plot of corresponding 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 time step has to be inversely proportional to the dimension of the subspace we are approximating our system on. Furthermore, every time step requires a computational effort which grows linearly with . Consequently, the complete CPU time for a simulation over a given time interval is quadratic in . So the dimension of the subspace and, thus, the related and should be the smallest possible. In order to check how small a we can choose without any significant loss of accuracy, we fix and and calculate the cricital temperature via the discretized version of Eq. (25),
for different values of . The result can be seen in Fig. 9.
For different values of and we get the same plot. We see that for the critical temperature is still slightly too small. However, when comparing the evolutions obtained with to the according ones for , the relevant figures are indistinguighable from each another. For the sake of efficiency, we thus fix for the rest of this work.
As for the extension of our interval, , we have to choose it large enough so that the solution cannot reach the boundaries during the simulation. As, by construction, we work with a periodic setting, a solution reaching one end of the interval would enter again at the other end, thus leading to unphysical interference. As an example of this numerical artifact, we consider the case , and plot the modulus of in Fig. 10.
We observe oscillations for larger which should not show up in reality, cf. Subsection 4.3. Whenever we encountered such an artifact, we successively increased by factors of until the artifact vanished.
- (1) R. L. Frank, C. Hainzl, B. Schlein, R. Seiringer, preprint, arXiv:1504.05885
- (2) J. Bardeen, L. N. Cooper and J. R. Schrieffer, Physical Review 108, 1175 (1957)
- (3) E. P. Gross, Il Nuovo Cimento 10, 454 (1961)
- (4) L. P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961)
- (5) L. D. Landau and V. L. Ginzburg, Zh. Eksp. Teor. Fiz. 20, 1064 (1950)
- (6) L. P. Gorkov, Sov. Phys. JETP 9, 1364 (1959)
- (7) R. L. Frank, C. Hainzl, R. Seiringer and J. P. Solovej, J. Amer. Math. Soc. 25, 667 (2012)
- (8) A. J. Leggett, Lecture Notes in Physics, vol. 115 (Springer, Berlin, 1980) pp. 13-27
- (9) P. Noziéres and S. Schmitt-Rink. Journal of Low Temperature Physics, 59:195–211, 1985.
- (10) M. Drechsler and W. Zwerger, Ann. Phys. 504, 15 (1992).
- (11) C. Hainzl and R. Seiringer, Lett. Math. Phys. 100, 119 (2012)
- (12) C. Hainzl and B. Schlein, Journal of Functional Analysis 265, 399 (2013)
- (13) P. Pieri and G.C. Strinati, Phys. Rev. Lett. 91, 030401 (2003)
- (14) M. J. Stephen and H. Suhl, Phys. Rev. Lett. 13, 797 (1964)
- (15) A. Schmid, Z. Phys. 216, 336 (1966)
- (16) E. Abrahams and T. Tsuneto, Phys. Rev. 152, 416 (1966)
- (17) L. P. Gorkov and G. M. Eliasberg, Sov. Phys. JETP 27, 328 (1968)
- (18) M. Cyrot, Rep. Prog. Phys. 36, 103 (1973)
- (19) C. A. R. Sá de Melo, M. Randeria and J. R. Engelbrecht, Phys. Rev. Lett. 71, 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 Dynamics: Reduced Models and Numerical Analysis (Europ. Math. Soc., Zürich, 2008)
- (22) M. D. Feit, J. A. Fleck Jr. and A. Steiger, J. Comput. Phys. 47, 412 (1982)
- (23) S. K. Gray and D. E. Manolopoulos, J. Comput. Phys. 104, 7099 (1996)
- (24) S. Blanes, F. Casas and A. Murua, J. Comput. Phys. 124, 234105 (2006)
- (25) H. Tal-Ezer and R. Kosloff, J. Chem. Phys. 81, 3967 (1984)
- (26) T. J. Park, J. Chem. Phys. 85, 5870 (1986)
- (27) W. Bao, D. Jaksch and P. A. Markowich, J. Comput. Phys. 187, 318 (2003)
- (28) M. Caliari, C. Neuhauser and M. Thalhammer, J. Comput. Phys. 228, 822 (2009)
- (29) L. Gauckler and C. Lubich, Found. Comput. Math. 10, 275 (2010)
- (30) Y.-F. Tang, L. Vázquez, F. Zhang and V. M. Pérez-García, Computers Math. Applic. 32, 73 (1996)
- (31) C. Lubich, Appl. Numer. Math. 48, 355 (2004)
- (32) C. Lubich, Math. Comp. 74, 765 (2005)
- (33) M. Léwin and S. Paul, ESAIM: Mathematical Modelling and Numerical Analysis 48, 53 (2014)
- (34) C. Hainzl, E. Lenzmann, M. Léwin and B. Schlein, Ann. H. Poincaré 11, 1023 (2010)
- (35) A. . Volkov and S. M. Kogan, Sov. Phys. JETP 38, 1018 (1974)
- (36) R. A. Barankov, L. S. Levitov and B. Z. Spivak, Phys. Rev. Lett. 93, 160401 (2004)
- (37) E. A. Yuzbashyan, O. Tsyplyatyev and B. L. Altshuler, Phys. Rev. Lett. 96, 097005 (2006)
- (38) V. Bach, E. H. Lieb and J. P. Solovej, Journal of statistical physics 76, 3 (1994)
- (39) P. G. de Gennes and P. A. Pincus Superconductivity of metals and alloys (WA Benjamin, New York, 1966), Vol. 86.
- (40) C. Hainzl, E. Hamza, R. Seiringer and J. P. Solovej, Comm. Math. Phys 281, 349 (2008)
- (41) R.L. Frank, C. Hainzl, R. Seiringer, J.P. Solovej, Operator Theory: Advances and Applications Volume 227, 57–88 (2013),
- (42) E. Hairer, C. Lubich and G. Wanner, Geometric numerical integration. Structure-preserving algorithms for ordinary differential 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 87, 084064 (2013)
- (45) G. Lukes-Gerakopoulos, Phys. Rev. D 89, 043002 (2014)
- (46) Note that in the linear case can also be calculated exactly.