Extended Ginzburg-Landau formalism: systematic expansion in small deviation from the critical temperature

# Extended Ginzburg-Landau formalism: systematic expansion in small deviation from the critical temperature

## Abstract

Based on the Gor’kov formalism for a clean -wave superconductor, we develop an extended version of the single-band Ginzburg-Landau (GL) theory by means of a systematic expansion in the deviation from the critical temperature , i.e., . We calculate different contributions to the order parameter and the magnetic field: the leading contributions ( in the order parameter and in the magnetic field) are controlled by the standard Ginzburg-Landau (GL) theory, while the next-to-leading terms ( in the gap and in the magnetic field) constitute the extended GL (EGL) approach. We derive the free-energy functional for the extended formalism and the corresponding expression for the current density. To illustrate the usefulness of our formalism, we calculate, in a semi-analytical form, the temperature-dependent correction to the GL parameter at which the surface energy becomes zero, and analytically, the temperature dependence of the thermodynamic critical field. We demonstrate that the EGL formalism is not just a mathematical extension to the theory - variations of both the gap and the thermodynamic critical field with temperature calculated within the EGL theory are found in very good agreement with the full BCS results down to low temperatures, which dramatically improves the applicability of the formalism compared to its standard predecessor.

###### pacs:
74.20.De, 74.20.Fg, 74.25.Bt, 74.25.Ha

## I Introduction

In 1950 Ginzburg and Landau proposed a phenomenological theory of superconductivity (the GL theory) based on a specific form of the free energy functional constructed in the vicinity of the critical temperature from Landau’s theory of second-order transitions. GL () Minimization of this functional generates the system of two GL equations that give the spatial distribution of the superconducting order parameter (the condensate wave function) and of the magnetic field in a superconductor. Over the years, the GL approach has been enormously successful in describing various properties of superconductors (see e.g. Ref. degen, ), and has been extensively used particularly in the last decade in the domain of mesoscopic superconductivity, see, e.g., Ref. mesoGL, .

In 1957 Bardeen, Cooper and Schrieffer (BCS) developed the well-known microscopic theory of the conventional superconductivity BCS () and then, two years later, Gor’kov showed that the GL equations can be obtained from the BCS formalism. gor () However, in spite of the availability of the microscopic theory, the GL approach remains still an optimal choice for many practical calculations when the spatial distribution of the pair condensate and, thus, of the magnetic field are nontrivial, e.g., for multiple-vortex configurations. The simple differential structure of the local GL equations in comparison with the microscopic theory offers clear advantages, including the possibility of analytical derivations in many important cases.

A desire to develop an extension to the GL theory, with the idea of improving the formalism while retaining at least some advantages of its original formulation, stimulated significant efforts by many theorists. Several GL-like theories of different complexity were proposed. In the earliest developments, tewordt (); werthammer () the so-called “local superconductor” formalism was attempted, being a complicated synthesis of the BCS and GL approaches. The GL theory with nonlocal corrections (i.e., including higher powers of the gradients of the order parameter) was used in studies of the anisotropy of the upper critical field (see Ref. takanaka, ) and the vortex structure (see Ref. ichioka, ) in -wave superconductors. Recently, various extensions to the GL theory were introduced in the context of studying the Fulde-Ferrell-Larkin-Ovchinnikov state. buzdin (); adachi (); mineev () In all the above examples, extending the GL theory was based on the expansion of the self-consistent gap equation by including higher powers of the order parameter and its spatial gradients phenomenologically. note1 () However, accounting for such higher-power terms is not as straightforward as it may seem. The fundamental problem here is to select properly all relevant contributions of the same order of magnitude (accuracy). This (serious) issue does not arise in the derivation of the original GL theory for a single-band superconductor, where only the first non-linear term and the second-order (leading) gradient of the condensate wave function are included. However, as recently shown in Ref. kogan, , a similar selection performed for the GL theory of a two-band superconductor leads to the appearance of incomplete higher-order contributions. Such incomplete terms may cause misleading conclusions and should be avoided.

To tackle this problem, one needs to work with a single small parameter in the expansion. In the present case, such a small parameter is the proximity to the critical temperature, i.e., . Indeed, the standard GL approach can be seen as the lowest-order theory in the -expansion of the self-consistent gap equation, see, e.g., Refs. kogan, and shanenko, . However, next orders in are also of great importance, for example, to capture the physics of different healing lengths of different condensates in multi-band superconductors. shanenko (); lucia () In the present paper, we show that next orders in are also important in the single-band case, surprisingly improving the GL theory. We obtain a systematic series expansion of the self-consistent gap equation for a single-band, -wave, clean superconductor, using as the governing small parameter. In the derivation we employ a technique in the spirit of the asymptotic expansion methods used extensively in the applied mathematics. babich () Similarly to the asymptotic analysis in other models, we obtain a hierarchical system of the so-called transport equations which need to be solved recursively starting from the lowest order. Using this method, corrections to the standard GL theory can in principle be calculated with an arbitrary accuracy. However, unlike asymptotic expansions for linear models, here the complexity of the higher-order equations increases rapidly and their solution cannot be obtained in the general form. Based on the Gor’kov Green function formalism, we derive and investigate the first three orders of the -expansion of the gap equation, i.e., with . To the order , we find the equation for the critical temperature. Collecting the terms proportional to , we obtain the standard GL theory giving the lowest-order (in ) contributions to the superconducting condensate, i.e., , and to the magnetic field, i.e., . Then, by matching the terms of the order , we derive equations for the next-to-leading corrections to the superconducting order parameter and magnetic field, and , respectively. The equations controlling the order parameter and the magnetic field up to the next-to-leading order in constitute the extended GL formalism (EGL). To illustrate the power of our extension to the GL theory, we investigate the energy associated with a surface between the superconducting (S) and normal (N) phases. In particular, we calculate the temperature dependence of , the value of the GL parameter at which the surface energy becomes zero. It is important to stress here that contrary to other available extensions of the GL theory discussed above, our formalism is not much more complicated than the standard GL theory. As is shown in the calculation of the S/N surface energy, plenty of important information can be obtained from the EGL formalism analytically.

Note that as is known, the GL theory is heavily used also in the study of thermal fluctuations. Here we do not address this issue but deal with an extension to the GL formalism in the mean-field level. Our aim is to expand the validity domain of the GL theory down to lower temperatures, which will be useful for the problems with a strongly nonuniform distribution of the pair condensate, e.g., for multiple vortex solution in the presence of stripes and clusters of vortices. In this case any fully microscopic approach will be an very complicated and rather time consuming task.

The paper is organized as follows. In Sec. II we introduce a general approach to construct the asymptotic expansion in for the self-consistent gap equation. In order to illustrate the main ideas behind the method, a simpler case of zero magnetic field is considered here. Sec. III presents the generalization to a nonzero magnetic field. Such a generalization requires the normal-metal Green function beyond the traditional Peierls (phase) approximation, and the corresponding expression is presented and discussed. The series expansion of the free-energy functional up to the next-to-leading order in is given in Sec. IV. The complete set of equations for the order parameter and the magnetic field in the EGL approach is derived by finding the stationary point of the functional in Sec. V. In Sec. VI we estimate the accuracy of the EGL formalism by comparing its results for the uniform order parameter and the critical magnetic field with those of the standard GL approach and the BCS theory. Here we demonstrate that the temperature in which the GL theory is valid, is dramatically increased due to the extension. In Sec. VII we investigate the S/N surface energy and calculate the temperature-dependent correction to . Finally, Sec. VIII presents the summary of our results and our conclusions.

The main text of the paper contains only the key formulas, the general ideas and the main steps of the derivation. Readers interested in details are referred to Appendices. In particular, Appendix A shows how to calculate the coefficients appearing in the -expansion of the gap equation in the absence of a magnetic field. Appendix B presents details of our calculations for the normal-metal Green function beyond the Peierls phase approximation. Appendix C generalizes the calculations given in Appendix A to the case of a nonzero magnetic field.

## Ii Series expansion in τ of the gap equation at zero magnetic field

As is known since the classical work by Gor’kov, gor (); degen (); fett () the GL equations can be derived from the microscopic BCS theory in the most elegant way via the Green function formalism. For the sake of clarity of presentation of our main ideas, the case of zero magnetic field is considered first, while the generalization to a nonzero magnetic field is given in the next section. The goal of our work is to construct the extended GL formalism through the expansion of the Gor’kov equations in , with the critical temperature (). We start by writing the Gor’kov equations as the Dyson equation for the Green function in the Gor’kov-Nambu matrix representation, which reads as (see, e.g., Ref. fett, ; zagoskin, ; kopnin, )

 ˇGω=ˇG(0)ω+ˇG(0)ωˇΔˇGω, (1)

with

 ˇGω=(GωFω˜Fω˜Gω),ˇG(0)ω=(G(0)ω00˜G(0)ω), (2)

where is the fermionic Matsubara frequency ( is an integer and is set to unity) and the matrix gap operator in Eq. (1) is defined by

 ˇΔ=(0^Δ^Δ∗0),⟨r|^Δ|r′⟩=δ(r−r′)Δ(r′), (3)

where the superconducting order parameter obeys the self-consistent gap equation, i.e.,

 Δ(r)=−gT∑ωFω(r,r), (4)

with the (Gor’kov) coupling constant. As usual, the sum in the r.h.s. of Eq. (4) is assumed to be properly restricted to avoid the ultraviolet divergence. Equations (1) and (2) further give

 Fω(r,r′)=∫d3yG(0)ω(r,y)Δ(y)˜Gω(y,r′), (5a) ˜Gω(r,r′)=˜G(0)ω(r,r′) +∫d3y˜G(0)ω(r,y)^Δ∗(y)Fω(y,r′). (5b)

For the normal-state temperature Green function we have (at zero magnetic field)

 G(0)ω(r,y)=∫d3k(2π)3eik(r−y)iℏω−ξk, (6)

with the single-particle energy measured from the chemical potential , and .

The Gor’kov equations (5) supplemented by Eq. (6) make it possible to express the anomalous (Gor’kov) Green function in terms of and . Then, inserting this expression into Eq. (4), one obtains the self-consistent gap equation. Solution to the gap equation can be represented in the form of a perturbation series over powers of  (which is small in the vicinity of the critical temperature ):

 Δ(r) =∫d3yKa(r,y)Δ(y)+∫3∏j=1d3yjKb(r,{y}3) ×Δ(y1)Δ∗(y2)Δ(y3)+∫5∏j=1d3yjKc(r,{y}5) ×Δ(y1)Δ∗(y2)Δ(y3)Δ∗(y4)Δ(y5)+…, (7)

where and the integral kernels are given by

 Ka(r,y)=−gT∑ωG(0)ω(r,y)˜G(0)ω(y,r), Kb(r,{y}3)=−gT∑ωG(0)ω(r,y1)˜G(0)ω(y1,y2) ×G(0)ω(y2,y3)˜G(0)ω(y3,r), (8) Kc(r,{y}5)=−gT∑ωG(0)ω(r,y1)˜G(0)ω(y1,y2) ×G(0)ω(y2,y3)˜G(0)ω(y3,y4)G(0)ω(y4,y5)˜G(0)ω(y5,r).

Equation (7) can be truncated to a desired order, which yields a non-linear integral equation. The latter is further converted into a non-linear partial differential equation by using the gradient expansion

 (9)

with the scalar product. In particular, the GL equation is obtained when keeping only the first two terms, including and , in Eq. (7) and the second-order spatial derivatives in the gradient expansion (9).

While Eq. (7) is a regular series expansion of the gap equation (4), the partial differential equation mentioned above is not. The gradient expansion introduces a second small parameter together with the corresponding truncation approximation, and the relation between the order parameter and its gradients is not known a priory. As a result one cannot compare the accuracy of the relevant terms in the expansion and the truncation procedure becomes ill-defined. This problem does not appear in the derivation of the GL equation where one keeps only the lowest second order gradient term. In order to extend the GL formalism, one has to deal with a single small parameter in the system that can be used to compare the relevant contributions. This small parameter follows from the solution of the GL equation. When the order parameter decays as . Also the solution reveals the scaling length , i.e., the GL coherence length, which determines the spatial variations of the order parameter in the vicinity of and dictates that , or, using the short-hand notation, . Thus, the parameter controls both relevant quantities and can be used to produce a single-small-parameter series expansion of the gap equation (4).

A systematic expansion of the gap equation in can be facilitated by introducing the scaling transformation for the order parameter, the coordinates and the spatial derivatives of the order parameter in the following form:

 Δ=τ1/2¯Δ,r=τ−1/2¯r,∇r=τ1/2∇¯r. (10)

Note that in terms of the scaled coordinates, the typical spatial variation of the order parameter occurs on a scale that is independent of to the leading order. After the transformation given by Eq. (10), the parameter enters the expansion in Eqs. (7) and (9) as follows. In Eq. (9), only coordinate changes. The scaling of the difference does not change the expressions as it is an integration variable, and thus the scaling will not be applied to it. As is now a function of , each derivative in the expansion (9) adds a factor . Inserting Eq. (9) with these factors into the expansion given by Eq. (7) and taking into account the factor for the order parameter in Eq. (10) we arrive at a simple mnemonic rule to count the minimal order of each term in the expansion of the gap equation: each occurrence of or in the formulas adds the factor . The final form of the expansion is obtained by calculating the relevant coefficients through the evaluation of the remaining integrals. As those coefficients depend on temperature, they can also represented as series in . The formulated procedure allows one to calculate the -expansion for the gap equation to arbitrary order. However, in practice, calculations of higher orders become more and more complicated. In this work we limit ourselves to the analysis of the self-consistent gap equation in the first three orders, i.e., up to the order . Collecting terms of the order , we obtain the equation for . Working in the order , we recover the standard GL theory producing the leading contribution to , i.e., . The order yields the equation that controls the next-to-leading contribution to , i.e.,  (this is what we call the EGL formalism). Details of the selection of all the necessary terms in Eq. (7) that contribute to one of the three orders mentioned above, are given in Appendix A. The final result reads ()

 τ1/2g ¯Δ=a1τ1/2¯Δ+a2τ3/2¯∇2¯Δ+a3τ5/2¯∇2(¯∇2¯Δ) − b1τ3/2|¯Δ|2¯Δ−b2τ5/2[2¯Δ|¯∇¯Δ|2+3¯Δ∗(¯∇¯Δ)2 + ¯Δ2¯∇2¯Δ∗+4|¯Δ|2¯∇2¯Δ]+c1τ5/2|¯Δ|4¯Δ, (11)

where the coefficients , and are obtained from the integrals with the kernels , and , respectively, and they are given by

 a1=AT−a(τ+τ22+O(τ3)),ATN(0)=ln(2eγℏωDπTc), b1=b(1+2τ+O(τ2)),b=N(0)7ζ(3)8π2T2c, c1=c(1+O(τ)),c=N(0)93ζ(5)128π4T4c, a2=K(1+2τ+O(τ2)),K=b6ℏ2v2F, a3=Q(1+O(τ)),Q=c30ℏ4v4F, b2=L(1+O(τ)),L=c9ℏ2v2F, (12)

where and is the DOS at the Fermi energy, with the Fermi velocity; denotes the Debye (cut-off) frequency, is the Euler constant and is the Riemann zeta-function. It is of importance to note that Eq. (11) contains only half-integer powers of . The reason for this is two-fold. First, due to the structure of Eq. (7), there appear only odd integer powers of the order parameter. Second, the spherical symmetry dictates that the integrals with an odd number of operators are equal to zero.

The solution to the gap equation (11) must also be sought in the form of a series expansion in . Based on Eqs. (11) and (12), we are able to introduce

 ¯Δ(r)=¯Δ0(r)+τ¯Δ1(r)+…. (13)

Substituting this into Eq. (11) and collecting terms of the same order we obtain a set of equations for each .

Collecting terms of the order , we obtain

 (g−1−AT)¯Δ0=0. (14)

The solution to this equation, i.e., gives the ordinary BCS expression for the critical temperature, i.e., .

In the order one recovers the standard GL equation for the leading contribution to the order parameter :

 a¯Δ0+b|¯Δ0|2¯Δ0−K¯∇2¯Δ0=0. (15)

Note that its standard form with the temperature dependent -coefficient is obtained by multiplying all terms by the factor and returning to the unscaled quantities, i.e.,

 aτΔ0+b|Δ0|2Δ0−K∇2Δ0=0.

Finally, collecting all terms of the order , we arrive at the equation for , i.e., the next-to-leading term in the order parameter:

 a¯Δ1+b(2|¯Δ0|2¯Δ1+¯Δ20¯Δ∗1)−K¯∇2¯Δ1=F, (16)

where is given by

 F= −a2¯Δ0+2K¯∇2¯Δ0+Q¯∇2(¯∇2Δ0) −2b|¯Δ0|2¯Δ0−L[2¯Δ0|¯∇¯Δ0|2+3¯Δ∗0(¯∇¯Δ0)2 +¯Δ20¯∇2¯Δ∗0+4|¯Δ0|2¯∇2¯Δ0]+c|¯Δ0|4¯Δ0. (17)

This is a linear differential inhomogeneous equation to be solved after is found from Eq. (15). Note that similar features in the -expansion of the gap equation appear for a two-band superconductor, as well. shanenko ()

We note again that in principle, one can continue the procedure up to arbitrary order in , obtaining corrections to the standard GL theory with desired accuracy. While the equation for is nonlinear, the higher order contributions to will be controlled by inhomogeneous linear differential equations. Such a system of equations is solved recursively, starting from the lowest order, since solutions for previous orders will appear in the higher order equations, but not vice versa. The solution to the system will thus be uniquely defined (when the relevant boundary conditions are specified), ensuring consistency of the developed expansion.

We also remark that the structure of Eq. (16) makes it possible to conclude that the next-to-leading term is not trivially proportional to . For that reason, the spatial profile of is different compared to . This means that the characteristic length for the spatial variations of in EGL differs from the standard GL coherence length. However, both lengths have the same dependence on , i.e., .

## Iii Series expansion in τ of the gap equation for nonzero magnetic field

The magnetic field enters the formalism developed in the previous section via changes in the normal-metal Green function. In the Gor’kov derivation the field-induced corrections are accounted for through the field-dependent Peierls phase factor as

 G(0)Gor(rt,r′t′)=eieℏcr∫r′A⋅dqG(0)B=0(rt,r′t′), (18)

where is the vector potential, i.e., . It is of importance to note that the integral in the exponent is evaluated along the straight line connecting and . This approximation leads to Eq. (15), where the gauge-invariant derivative replaces . Obtaining corrections to the GL equation requires the Green function to be calculated with an accuracy sufficient to produce the complete set of terms contributing up to the order in the -expansion of the gap equation. Taking into account Eqs. (11) and (12), one concludes that the normal-metal Green function must be calculated with the accuracy .

Accounting for the magnetic field in the expansion can be done by noting that the critical magnetic field in a superconducting system is proportional to . Similarly a solution of the GL equations for the field also changes as . Thus the derivation of the expansion for the field corrections can be conveniently done by introducing the following scaling for the magnetic field:

 A=τ1/2¯A, B=rotA=τ¯B, (19)

so that the critical field quantities become independent on in the first order. We note that the spatial dependence of the magnetic field can also be represented in the form of the gradient expansion as it is done for the order parameter in Eq. (9). As the characteristic scale for the spatial variations of the magnetic field is defined by the magnetic penetration depth , which has the same dependence as the GL coherence length , we can again employ the scaled spatial coordinates as those in Eq. (10). So, the gradient expansion for the magnetic field follows the same rule as for the order parameter: each introduces an additional factor .

The -expansion of the Green function can now be done using the path integral method, by accounting of the quantum fluctuations around the classical trajectory. For details of the calculations we refer the reader to Appendix B. The final expression for the Green function that contains all contributions to the order reads as

 G(0)ω(r,r′)=eieℏcr∫r′A⋅dq{1+e2τ2¯B2(r)24m2c2 ×[∂2ω+iℏm(r−r′)2⊥∂ω]+O(τ5/2)}G(0)ω,B=0(r,r′), (20)

where the phase factor in the exponent is the same as in Eq. (18). We remark that the Peierls phase also contains terms of higher orders than . However, they do not enter the final equations for the next-to-leading contribution to the order parameter and magnetic field. It is simply convenient to represent the Green function in the form with the Peierls factor because it naturally leads to the appearance of the gauge invariant spatial derivatives of the order parameter. One also notes that this factor is written using the unscaled quantities. In fact, the scaling does not affect it. In addition, we do not scale the difference of the coordinates as it is an integration variable. The new field-dependent term in Eq. (20) is proportional to . It does not follow from the Peierls phase and, so, is not present in the approximation given by (18). It is interesting that the field gradients do not contribute to the correction to the Peierls approximation.

The field-modified expansion of the gap in powers of the order parameter and its spatial derivatives are obtained by substituting Eq. (20) into the kernels in Eq. (8) and proceeding in a manner similar to that discussed in Appendix A for Eq. (11). It is convenient to remove the phase factor from the Green functions and introduce the “two-point” auxiliary order parameter, i.e.,

 ¯Δ(¯r,¯r′)=e−2ieℏc¯r∫¯r′¯A⋅d¯q¯Δ(¯r). (21)

Then, as shown in Appendix C, the modification of Eq. (11) due to the phase factor is that is replaced by , and the limit is implied after all relevant calculations (we remark that such a limit is not permutable with the differentiating with respect to ). In particular, for the first three terms in the modified Eq. (11) we have

 a1τ1/2limr′→r¯Δ(¯r,¯r′)=a1τ1/2¯Δ, (22a) a2τ3/2limr′→r¯∇2r¯Δ(¯r,¯r′)=a2τ3/2¯D2¯Δ, (22b) a3τ5/2limr′→r¯∇2r(¯∇2r¯Δ(¯r,¯r′))= =a3τ5/2[¯D4−4ie3ℏc¯rot¯B⋅¯D+4e2ℏ2c2¯B2]¯Δ, (22c)

with , the gauge invariant gradient, and . The terms related to the kernel in the field-modified Eq. (11) read

 − b1τ3/2limr′→r|¯Δ(¯r,¯r′)|2¯Δ(¯r,¯r′)=−b1τ3/2|¯Δ|2¯Δ, (23a) − b2τ5/2limr′→r[2¯Δ(¯r,¯r′)|¯∇r¯Δ(¯r,¯r′)|2+3¯Δ∗(¯r,¯r′) ×(¯∇r¯Δ(¯r,¯r′))2+¯Δ2(¯r,¯r′)¯∇2r¯Δ∗(¯r,¯r′) +4|¯Δ(¯r,¯r′)|2¯∇2r¯Δ(¯r,¯r′)]=−b2τ5/2[2¯Δ|¯D¯Δ|2 +3¯Δ∗(¯D¯Δ)2+¯Δ2(¯D2¯Δ)∗+4|¯Δ|2¯D2¯Δ], (23b)

whereas the term coming from the integral with the kernel is of the form

 c1τ5/2limr′→r|¯Δ(¯r,¯r′)|4¯Δ(¯r,¯r′)=c1τ5/2|¯Δ|4¯Δ. (24)

In addition to the contributions that appear due to the Peierls factor in the Green function, we also obtain an extra contribution to the r.h.s. of Eq. (11) due to the terms proportional to in Eq. (20). This contribution comes only from the integral involving the kernel  (when keeping terms up to the order ) and reads

 −a4τ5/2¯B2(¯r)limr′→r¯Δ(¯r,¯r′)=−a4τ5/2¯B2¯Δ, (25)

where , with given by Eq. (12).

Using the modified Eq. (11) together with Eqs. (22)-(24) and collecting terms of the same order in , one can generalize Eqs. (15) and (16) to the case of a nonzero magnetic field. However, since the magnetic field is also a variable in the GL theory, it needs to be found self-consistently from a complementary set of equations. The most elegant way to derive the complete set of equations for the magnetic field and for the order parameter is based on the free-energy functional that accounts for the energy associated with the presence of the magnetic field and the superconducting pairing. This functional must be also represented as a series expansion in , which is the subject of the next section.

## Iv Free-energy functional

The free-energy functional can be obtained, e.g., by using the path integral methods developed for the BCS theory. popov () Its expansion in reads

 Fs =Fn,B=0+∫d3rB2(r)8π +1g∫d3rd3y[δ(r−y)−Ka(r,y)]Δ∗(r)Δ(y) −12g∫d3r3∏j=1d3yjKb(r,{y}3)Δ∗(r)Δ(y1) ×Δ∗(y2)Δ(y3)−13g∫d3r5∏j=1d3yjKc(r,{y}5) ×Δ∗(r)Δ(y1)Δ∗(y2)Δ(y3)Δ∗(y4)Δ(y5)−…, (26)

where denotes the free energy of the normal state ( stays for the zero magnetic field). Here it is worth noting that, generally,

 K∗a(r,y) =Ka(y,r), K∗b(c)(r,{y}3(5)) =Kb(c)({y}3(5),r),

which means that the r.h.s. of Eq. (26) is a real quantity (as necessary for the free energy). In addition, we have

 Kb(r,y1,y2,y3) =Kb(y2,y3,r,y1), Kc(r,y1,y2,y3,y4,y5) =Kc(y2,y3,r,y1,y4,y5) =Kc(y2,y3,y4,y5,r,y1),

which makes it possible to immediately find that the extremum condition of this functional with respect to leads to Eq. (7).

Series expansion of the free-energy functional given by Eq. (26) is obtained by following essentially a similar approach as is used in previous sections. The scaling transformation for the order parameter, coordinates and magnetic field is introduced, and the gradient expansion for the order parameter (9) is substituted into Eq. (26). As before, the coefficients of the series expansion appears as the integrals with kernels . As we are interested in the terms that are used to derive the GL equations and the equations for the next-to-leading contributions to the order parameter and the magnetic field, we must expand the functional up to the order . This follows from the rules for counting the powers of introduced in Sec. II. A calculation of the series expansion in of the free-energy functional is greatly simplified by noting that all terms in Eq. (26) can be derived from the corresponding terms in Eq. (11) by multiplying the latter by , where is the number of -factors of the corresponding integral term in Eq. (26). As a result, we obtain the functional, which can be represented in the following symmetric (real) form (here and below we omit bars over the scaled quantities unless it causes confusion):

 fs =fn,B=0+B28π+1τ(g−1−a1)|Δ|2+a2|DΔ|2 −τa3(|D2Δ|2+13rotB⋅i+4e2ℏ2c2B2|Δ|2) +τa4B2|Δ|2+b12|Δ|4−τb22[8|Δ|2|DΔ|2 +(Δ∗)2(DΔ)2+Δ2(D∗Δ∗)2]−τc13|Δ|6, (27)

where  (with the scaling ) and is given by

 i=i2eℏc(ΔD∗Δ∗−Δ∗DΔ). (28)

We stress that is not the current density . However, when replacing and in Eq. (28), we find being proportional to , the leading contribution to the current density , see the next section. We also note that the representation of the functional in a real symmetric form as in Eq. (27) implicitly relies on the disappearance of the corresponding surface integrals. It is easy to check that this is ensured by the boundary conditions discussed in Sec. V.

Obtaining the final series expansion in for the free-energy density requires the -expansion for the coefficients given in Eqs. (12) and (25), for the order parameter given by Eq. (13), and for the magnetic field. The latter is expressed in the form

 A =A0+τA1+…, B=B0+τB1+…. (29)

After collecting the relevant terms, the free-energy density is written in the form

 fs−fn,B=0=f0+τf1+O(τ2), (30)

where the leading-order term (the standard GL functional) is specified by

 f0 =B208π+a|Δ0|2+b2|Δ0|4+K|D0Δ0|2, (31a) and the next-to-leading contribution f1 can be written as a sum of two terms, i.e., f1=f(0)1+f(1)1 with f(0)1=a2|Δ0|2+2K|D0Δ0|2−Q(|D02Δ0|2 +13rotB0⋅i0+4e2ℏ2c2B20|Δ0|2)+b36e2ℏ2m2c2B20|Δ0|2 +b|Δ0|4−L2[8|Δ0|2|D0Δ0|2+(Δ∗0)2(D0Δ0)2 +Δ20(D∗0Δ∗0)2]−c3|Δ0|6 (31b) and f(1)1 =B0⋅B14π+(a+b|Δ0|2)(Δ∗0Δ1+Δ0Δ∗1) +K[(D0Δ0⋅D∗0Δ∗1+c.c.)−A1⋅i0]. (31c)

Here is equal to with the substitution , and is obtained from Eq. (28) by and .

## V Complete set of equations for B≠0

A system of coupled equations for both and are obtained from the extremum condition of the free energy given by Eq. (27). The subsequent expansion of the obtained expressions yields the GL equations as well as the equations for the next-to-leading contributions to the order parameter and the magnetic field. Alternatively, this full set of equations can also be obtained by finding the extremum of the free energy with respect to and in such a way that all terms appearing in Eq. (30) are taken separately. Note that finding the extremum of the functional with respect to and yields the same set of equations. This property follows from Eqs. (13) and (29).

By searching for the minimum of the functional with the density , we reproduce the standard GL equations:

 aΔ0+b|Δ0|2Δ0−KD20Δ0=0, (32a) rotB0=4πcj0, (32b)

where . Likewise, the same equations are obtained by finding the extremum of with respect to and .

The minimum of the functional with the density with respect to and yields the following equations:

 aΔ1+b(2|Δ0|2Δ1+Δ20Δ∗1)−KD20Δ1=F, (33a) rotB1=4πcj1. (33b)

The r.h.s. of Eq. (33a) is given by

 F =−a2Δ0+2KD20Δ0+Q[D20(D20Δ0) −i4e3ℏcrotB0⋅D0Δ0+4e2ℏ2c2B20Δ0] −b36e2ℏ2m2c2B20Δ0−2b|Δ0|2Δ0−L[2Δ0|D0Δ0|2 +3Δ∗0(D0Δ0)2+Δ20(D20Δ0)∗+4|Δ0|2D20Δ0