Incompatibility of time-dependent Bogoliubov--de-Gennes and Ginzburg--Landau equations

We study the time-dependent Bogoliubov--de-Gennes equations for generic translation-invariant fermionic many-body systems. For initial states that are close to thermal equilibrium states at temperatures near the critical temperature, we show that the magnitude of the order parameter stays approximately constant in time and, in particular, does not follow a time-dependent Ginzburg--Landau equation, which is often employed as a phenomenological description and predicts a decay of the order parameter in time. The full non-linear structure of the equations is necessary to understand this behavior.


Introduction
The Ginzburg-Landau (GL) model [1] is a paradigm for the phenomenological description of phase transitions in physical systems. In the static case, its relation to the microscopic BCS theory [2] was clarified by Gor'kov [3] (see also [4] and [5] for alternative approaches), and a mathematically rigorous derivation of the GL model from BCS theory was recently given in [6,7]. This work is concerned with the analogous question in the time-dependent case. Arguments for the validity of a time-dependent GL equation can be found in the literature, starting with [8][9][10] in the case of superconductors. These attempts were critically analyzed in [11] where it was argued that the equation can only hold in a gapless regime, however; we refer to [12] for a thorough discussion. By applying similar arguments in the study of superfluid cold gases, it was stated in [13,14] that the non-linear time-dependent GL equation applies to such systems for temperatures T slightly above the critical temperature.
The main message of our present work is that a time-dependent GL equation cannot be derived from the time-dependent Bogoliubov-de-Gennes (BdG) equations near the critical temperature, in contrast to the static case. We shall consider the BdG equations for a translation-invariant system, which can be derived in a standard way applying the BCS approximation, either to the Heisenberg equations of motion for the fermion field operators, or the time-dependent Green's function [12,15]. We consider a system which is initially close to a thermal equilibrium state near the critical temperature, with non-vanishing order parameter. We then show that, in contrast to what would be expected from GL theory, the order parameter does not decay in time. Interestingly, this is a purely non-linear effect. If, instead, we consider the corresponding linear equation, then the solution indeed decays exponentially in time, on a timescale that can be calculated via the imaginary part of the corresponding resonance pole, which turns out to be proportional to the inverse of the distance to the critical temperature. The non-linear terms formally contribute a term cubic in the order parameter, which indeed resembles a GL type equation; however, small denominators close to the Fermi sphere invalidate such formal arguments, and prevent the decay on all timescales.
Our claims are mathematically rigorous and are not based on any ad hoc approximations; they are confirmed numerically in [16] in the case of a onedimensional system with contact interactions.
The present work can be viewed as a continuation of a recent series of studies of the mathematical properties of the BCS theory of superconductivity and superfluidity [17][18][19][20]. It is motivated by the current interest concerning the applicability of the theory to cold gases, in particular concerning the BCS-BEC crossover [21]. In the BCS regime [22,23], a rigorous proof of the emergence of a static GL equation close to the critical temperature was given in [6,7]. In the BEC regime, the prediction [13,14,24,25] of the emergence of the Gross-Pitaevskii equation in the low-density limit was rigorously established, both in the static [26] and the dynamical case [27].

Model and Main Result
The starting point of our analysis is the BCS model [2,17,22]. The state of a (translation-invariant, three-dimensional 1 ) system is described in terms of two quantities, the momentum distribution γ (k) = a † k a k and the pair densityα(k) = a k a −k , determining the Cooper pair wave-function via Fourier transform as α(x − y) = (2π) −3/2 α(k)e ik·(x−y) d 3 k. They can be conveniently combined to the 2 × 2 matrix which satisfies 0 ≤ (k) ≤ I C 2 for every k. We suppress spin in our notation; the pair densityα is assumed, for simplicity, to be a spin singlet. The equilibrium state at temperature T ≥ 0 and chemical potential μ is determined by minimizing the pressure functional with entropy We find it convenient to choose units such that = 2m = 1, with m denoting the particle mass. In (1), V (x) denotes a local two-body interaction potential, as appropriate for the description of superfluid gases. Our results can easily be generalized to include non-local interactions as well, which effectively arise via interactions through phonons in solids, for instance. It was shown in [17] that the critical temperature T c for the model (1) is the unique T for which the operator K T (−i∇) + V (x) has 0 as its lowest eigenvalue, where That is, for T ≥ T c the pressure functional (1) is minimized by the normal state n (k) with γ n (k) = (1 + e (k 2 −μ)/T ) −1 andα n (k) ≡ 0, while for T < T c the pairing densityα(k) does not vanish identically, showing the transition to a superfluid (or superconducting) phase. Strictly speaking, this characterization only applies if T c is strictly positive, which we shall assume henceforth. We shall also assume that K T c (−i∇) + V (x) has a unique normalized eigenvector α * corresponding to the eigenvalue 0; for radial V (x) this corresponds to the assumption that α * is an s-wave. 2 Smoothness and decay properties of α * are analyzed in detail in [6, with ψ ∈ C of order √ T c − T , and ξ(k) √ T c − T for small T c − T . The order parameter ψ is, in fact, determined by minimizing the GL type expression for a suitable constant C GL > 0. This follows from the analysis in [6], which is actually more general and not restricted to the translation-invariant case considered here.
The time-dependent BCS equation, also known as the BdG equation, has the form with effective Hamiltonian It can be derived from the Heisenberg equations of motion for the fermion field operators, applying the same BCStype approximations as in the static case [12,28,29]. Alternatively, it also follows from the time-dependent Green's function method introduced in [15] (see also [30]). While only certain interaction terms are retained in the BCS approximation, the equation (4) allows for pair creation and annihilation and hence the superfluid density is in general not a conserved quantity. Note that H t depends on t through t , hence this equation is non-linear. It is also important to note that the pressure functional (1) is preserved by the time evolution (4).
In this paper, we study the time evolution (4) for initial states close to the normal state, for temperatures close to T c . Closeness is measured in terms of a small parameter, which we call h. Let us assume that |T − T c | ≤ h 2 , and also that the initial state 0 satisfies for some constant C > 0 (independent of h). 3 This is satisfied, for instance, by thermal equilibrium states in the temperature range considered, but also by states where the order parameter ψ in (2) is modified by a factor of order one. In addition to the assumption T c > 0, we assume that μ > 0 and that α * does not vanish identically on the Fermi sphere, i.e., which is satisfied generically. With α t denoting the pairing density of t we define, for every time t, the complex number ψ t = h −1 α * |α t . 3 Throughout the paper, we denote by C generic constants, even if they take different values at different places.
For convenience, we multiply by h −1 for ψ to be of order one. Our main result states that for small h, |ψ t | remains approximately constant, uniformly in time. THEOREM 1. Let t be the solution of (4) with initial state 0 satisfying (5), and |T − T c | ≤ h 2 . Then there exists a constant C > 0 such that, for small h, for all times t.
We remark that the conditions of the theorem allow |ψ 0 | to take any value of order one. The result states that this value remains constant to leading order in h. This holds true even for T > T c , in which case the normal state n minimizes the pressure functional (1). In other words, the order parameter does not tend towards the minimum of the GL energy (3), as would be expected on the basis of a timedependent GL equation of the form [13,14] id∂ t ψ = aψ + b|ψ| 2 ψ with a ∈ R, b > 0 and d ∈ C with d > 0. In fact, for T > T c one has a > 0, in which case the solution to (8) goes to zero as t → ∞, in contrast to our main result (7). Moreover, at T − T c < 0 (but small), one could for instance start with a state with the "wrong" ψ, i.e., with α(k) the equilibrium pairing density multiplied by a complex number of modulus not equal to one, and our theorem states that this structure will be preserved at all times. We emphasize that our results do not rule out the validity of the time-dependent Ginzburg-Landau equation, in general, which has been successfully employed over several decades. What they show, however, is that such an equation cannot be derived from the Bogoliubov-de-Gennes equations, which also appear prominently in the physics literature. From the point of view of mathematical physics, it thus remains a challenging open problem to unveil the relevant additional physical effects which are responsible for the possible emergence of a time-dependent GL equation.
Proof. The proof of Theorem 1 is divided into three steps.
Step 1. The first step is to show that the energy bound F( ) − F( n ) ≤ Ch 4 implies a decomposition of the form for an appropriate constant C > 0. These bounds follow from the analysis of [6]; we shall sketch the main ideas here. From the bound [6, Lemma 1] on the relative entropy of with respect to n we deduce that F( ) − F( n ) can be bounded from below by For T ≥ T c , we can bound K T ≥ K T c in the first term. Since α * is, by definition, the non-degenerate zero eigenvector of K T c + V , with a spectral gap κ > 0, we conclude that the first term in (10) is bounded from below by ξ |(K T c + V )ξ ≥ κ ξ 2 2 in this case. Moreover, the second term is bounded from below by 2T γ − γ n 2 2 since K T ≥ 2T . Hence, we obtain the second and third bound in (9). Moreover, with the aid of the last term in (10) it is not difficult to show that |ψ| ≤ C, concluding the proof of (9) for T ≥ T c . The case T < T c is very similar, using the fact that K T ≥ K T c + 2(T − T c ) ≥ K T c − 2h 2 instead, and we refer to [6] for the details.
Step 2. Since the dynamics (4) is unitary, the eigenvalues of the 2 × 2 matrix t (k) are conserved, for all k ∈ R 3 . A simple computation shows that they are of the form 1/2 ± s t (k), with that is, also s t (k) is independent of t.
After integration over δ , the right side of (13) is thus bounded from below by Ch 2 δ |ψ t | 2 − |ψ 0 | 2 − Ch 3 δ 1/2 − Ch 4 using again the Cauchy-Schwarz inequality and (9). Together with the bound on the left side above, this implies that The choice δ = h leads to the claim (7).

Comparison with Linear Case
It is interesting to observe that our main result, namely the fact that |ψ t | remains approximately constant in time, crucially depends on the non-linear terms in the time-dependent BCS equation. Let us explain this point in more detail. The equation for α t in (4) is given by Writing again γ t (k) = γ n (k) + η t (k) this can be rewritten abstractly as where S denotes the operator K T (−i∇) + V (x), L is multiplication by 2 − 4γ n (k) = 2 tanh((k 2 − μ)/2T ), V is multiplication by V (x) in x-space and η t is multiplication by η t (k) in k-space. Consider, for simplicity, the case T > T c ; in this case, the operator S is positive and the solution to (14) satisfies for t > 0, where is the unitary evolution generated by S 1/2 L S 1/2 . In the second term on the right side of (15), we can use (12) to express η s in terms of α s ; this leads to a non-linear equation for α t . Let us neglect for a moment this second term, and let us focus on the linear dynamics S −1/2 U (t)S 1/2 α 0 . The spectrum of S 1/2 L S 1/2 coincides with the one of L S. Its continuous spectrum can easily be seen to cover the halfline [−2μ; ∞), since L S − 2(k 2 − μ) is relatively compact with respect to k 2 . Moreover, for T = T c , L S has an eigenvalue 0 associated with the eigenvector α * , which is embedded in the continuous spectrum. Perturbation theory predicts that for T > T c the zero eigenvalue turns into a complex resonance λ, with real and imaginary parts of the order T − T c .
It is particularly simple to find the resonance λ of S 1/2 L S 1/2 if the potential V is rank one, i.e., of the form V = −|ϕ ϕ| for a ϕ ∈ L 2 (R 3 ) which we assume to be radial for simplicity. In this case, T c is determined by ϕ, K −1 T c ϕ = 1 and α * is proportional to K −1 T c ϕ. To compute λ, we use complex dilation. For θ ∈ R, we define the unitary operator u(θ ) by Alternatively, u(θ ) = e iθ A with A = x · k + k · x. Assuming ϕ to be an analytic vector for A, we can extend ϕ θ = u(θ )ϕ to a strip −β < Im θ ≤ 0 below the real axis. In this way, we can also define the operators L θ and S θ for all −β < Im θ ≤ 0. The resonance λ then satisfies the eigenvalue equation L θ S θ χ θ = λχ θ , which is equivalent to Multiplying from the left with φ θ |, we obtain Note that λ vanishes at T = T c . We can use implicit differentiation with respect to T to expand around T = T c , letting θ → 0 afterwards. This yields to leading order in T − T c , where the integral in the last factor is understood in the principal value (p.v.) sense. The fact that Im λ < 0 suggests that the corresponding state decays exponentially in time. In particular, one would expect from the linear evolution (16) that the order parameter satisfies |ψ t | ≈ |ψ 0 |e t Im λ to leading order in T − T c , i.e., it decays on a timescale of the order (T − T c ) −1 . Such a decay was in fact predicted in [12][13][14]. The meaning of the ≈ sign here can be made precise following the analysis of [31], but the details are not relevant here. A comparison with the statement of Theorem 1 shows the importance of the second, non-linear term on the right side of (15) for understanding the behavior of α t . Let us examine it closer. The function η t (k) is determined by Equation (12). Away from the Fermi sphere, the second term on the left side of (12) dominates, and we have to leading order, which leads to a cubic equation for the evolution of α t . In terms of ψ t this would even resemble the cubic GL term. On the Fermi sphere the denominator on the right side of (17) vanishes, however. As a consequence, η t is much larger, of the order ||α t (k)| 2 − |α 0 (k)| 2 | 1/2 according to (12). The latter expression equals h|α * (k)|||ψ t | 2 − |ψ 0 | 2 | 1/2 to leading order. Sinceα * does not vanish on the Fermi sphere, we conclude that |ψ t | |ψ 0 |; otherwise, η t (k) would be too large to satisfy (9). This is exactly the mechanism used in the proof of Theorem 1, explaining why the non-linear term in (15) plays such an important role in determining the behavior of ψ t .