Accurate Symbolic Solution of Ginzburg-Landau Equations in the Circular Cell Approximation by Variational Method : Magnetization of Ideal Type II Superconductor

In this paper I access the degree of approximation of known symbolic approach to solving of Ginzburg-Landau (GL) equations using variational method and a concept of vortex lattice with circular unit cells, refine it in a clear and concise way, identify and eliminate the errors. Also, I will improve its accuracy by providing for the first time precise dependencies of the variational parameters; correct and calculate magnetisation, compare it with the one calculated numerically and conclude they agree within 98.5% or better for any value of the GL parameter k and at magnetic field 2 0.01 1 c B B ≤ ≤ , which is good basis for many engineering applications. As a result, a theoretical tool is developed using known symbolic solutions of GL equations with accuracy surpassing that of any other known symbolic solution and approaching that of numerical one.


Introduction
Much of basic superconductor behavior in magnetic field can be understood from the phenomenological model expressed by two Ginzburg-Landau (GL) equations [1] [2].Known numerical methods produce excellent and reliable results along with "…difficulty of a numerical solution of the complex-valued GL equations" [3] [4], the demands of calculation effort and time and being less transparent.Calculated magnetisation curves for vortex lattices with rectangular, hexagonal or circular unit cells coincide within line thickness [3].This fact mo-tivates solving GL equations symbolically that currently lags behind due to the complexity of the problem [5]- [10].The lower critical field B c1 can be accurately calculated symbolically, see Section 6.1.In strong magnetic fields the symbolic approximation [7] holds, see Section 4.2.At medium magnetic fields a symbolic approximation considering vortex lattice of circular unit cells is developed [5] [6] [8] [10].The drawbacks of [5] [6] are dealt with in [8] [9] [10] and will not be discussed here.As stated in [8] [9] [10], their method approximately describes magnetic properties of type II superconductors with periodic lattice over the entire range of magnetic fields and for any value of the GL parameter with λ, ξ London penetration depth and the coherence length respectively.Moreover, the magnetic field dependence of (reversible) magnetisation can be calculated in seconds.However, a simple check shows that magnetisation calculated with this method ( [8] [9] [10], as published) when compared to that calculated numerically [3], has relative differences exceeding 25%, which raises concerns about the accuracy.Clearly, the accuracy of the obtained solution of GL equations is vital: a symbolic one only has added value when its accuracy is comparable to that of the numerics.In this paper I will access the degree of approximation of the method [8] [9] [10], refine it in a clear, concise and rational way, identify and eliminate the errors.Furthermore, I will improve accuracy of the method by providing for the first time accurate dependencies of the variational parameters; correct and calculate magnetisation, compare it with the one calculated numerically and conclude they agree for any k within 98.5% or better at magnetic field 2 0.01 1 c B B ≤ ≤ , which is sufficient for many engineering applications.

Normalisation
In this paper I omit the time-dependent terms in the GL equations written in SI units [11] [12] and use almost the same normalisation as in [3] [11] [12], except that I normalise all length-related quantities by the penetration depth k λ ξ = so that both B c2 and λ are defined through ξ and Annex provides more details.
Below there are six mean (averaged over the unit cell area) magnetic field magnitudes (magnetic flux densities) in use: thermodynamic (equilibrium, applied) field B a , magnetisation M (from the induced currents), the resulting (total) field: with the local field b r defined by GL equations in Table 1; the lower, the upper and thermodynamic critical magnetic fields: B c1 , B c2 and     q r ϕ , Table 1.Two GL equations written in forms 1, 2 for a unit cell with axial symmetry.

Variational Method
The variational method [5] [6] [8] [9] [10] assumes for the dimensionless modulus ψ of the complex-valued order parameter ( ) ( ) where the notation variational parameters representing respectively the depression of the order parameter due to overlapping vortices and the effective core radius of a vortex; r, φ are the radial coordinate of cylindrical coordinate system and the phase angle respectively (a vortex line is centered on the z-axis, so that ( ) arctan y x ϕ = ; x, y, z are rectangular system coordinates, etc. [8]).Both variational parameters f ∞ and v ξ depend only on magnetic field b and on the GL parameter k (Figures 1-4), so that e.g., d d r r r r ρ ρ = .The order parameter ψ is interpreted as the local density of the Cooper pairs and its phase ϑ is de- termined by the electric potential [12].

Free Energy Density of Superconductor
Averaged over the cell area Helmholtz dimensionless free energy density F of a circular unit cell (with radius R) has two contributions [3] ,    ; F em -is related to super-current and to magnetic field; the first term of F core -to the local density of the Cooper pairs, and the second term: is the area of the circular (Wigner-Seitz) unit cell carrying one magnetic flux quantum, thus and Table A5 lists the scaling factors.

Ginzburg-Landau Equations
The GL equations are obtained by minimising the free energy of superconductor F with respect to e.g., f r and to q r (form 1, Table 1) or to f r and b kr (form 2).In the circular cell approximation both the order parameter and magnetic flux density have axial symmetry and for this case two stationary GL equations in two forms are listed in Table 1, where the five local quantities: r f (with are respectively moduli of the order parameter, of vector potential (satisfying the Coulomb gauge and having only φ-component [5]), of local magnetic flux density vector (having only z-component), of the dimensionless super-velocity (having only φ-component), and of the vortex current density vector (having only φ-component).The two unknowns are: r f and one of the following: r a or kr b respectively for the forms 1 or 2.
Complementing the Equations ( 1 q q b j q r r r r r r r r r r

Solving Ginsburg-Landau Equation(s)
From Equations ( 4) and ( 5) one gets: Since 0 v ξ ≠ (Section 3.4), also 0 r ρ ≠ everywhere, see Equation (3).On the other hand, Equation ( 11) is identical to Equation (4) everywhere except at 0 f ∞ = , which is the case at Therefore, I conclude that at 1 b → the approach [7] must complement [8] [9] [10] as further elaborated in Section 4.2.At 0 f ∞ ≠ (and with f ∞ independent on r) Equation ( 11) is the modified Bessel equation with the solution: where is magnetic field strength, 0 I and 0 K are the modified Bessel's functions respectively of the first and second kind and order 0; c 1 and c 2 are the integration constants set by the boundary conditions.So far the solution is the same as obtained in [8] [9] [10].Moreover, from Equation (12) one obtains: where 1 I and 1 K are the modified Bessel's functions respectively of the first and second kind and order 1.Note that the derivation of the solution (Equations ( 10)-( 13)) was skipped in [8] [9] [10] and as a result, the important restriction 0 f ∞ ≠ of the solution was hidden so far.
In order to stay focused and limit the size of the paper, below I simply accept Equations ( 14), (15) and will study namely this case in more detail.With the constants c 1 and c 2 defined, the solution allows calculating the local quantities as well as the mean quantities: equilibrium magnetic field, magnetisation, etc.In this paper I focus on magnetisation.

Variational Parameter f∞
At any 5 200 k ≤ ≤ all obtained data points of the dependence ( ) lapse at one curve fitted here by a cubic spline, see Figure 1 and Table A1.For this reason no data for the individual splines is given here for this range.This simplification causes an estimated error of 0.5% as explained in Section 5 and the error can be reduced by using more accurate data from the minimising of F. Furthermore, the obtained data are in reasonable agreement with that from [ [8], Equation ( 14)] (numbered as Equation (18) here), except at lower ( 0.25 b < ) and at higher ( 0.75 b > ) fields, see Figure 1: ( ) ( ) At 12 5 k k ≤ < the obtained dots deviate from this curve and the error when using Equation (18) can reach 5% as Figure 1 exemplifies.Therefore in this range instead of using Figure 1 or Equation (18), I recommend using Figure 1 or more accurately the data from Annex e.g., by constructing the individual splines.Importance of this correction becomes clear in Section 5.For the above reasons I am convinced that Equation (18) does not approximate the dependence of ( ) with "an accuracy of about 0.5% for arbitrary k and b" as stated in [8] [9] [10].In addition, I find unsatisfactory the agreement of the obtained data with the Equation (24) of [6]:

Variational Parameter ξv
Minimisation of F with respect to v ξ at 0 b → shows that Equation (15) of [8] is less accurate than the original Equation ( 16) of [5], especially at k < 10 (and that the error exceeds 1% up to k = 50), see Figure 2. Therefore, the correct equation to calculate the dependence ( ) The obtained (by minimizing F) dependence ( ) is similar and different from those described in [3] [5].On one hand at 75 k ≥ all data points collapse practically at one curve, see Figure 3 and Table A2 in Annex.For this reason only one cubic spline fit (namely, k = 100) is shown, the error of this simplification is below 0.5%.
The agreement of the obtained data with [[8], Equation ( 13)] is unsatisfactory, see Figure 3 and therefore I recommend calculating ( ) and any value of magnetic field b from the spline fit, Figure 3 and Table A2 in Annex.
On the other hand at 50 k < the obtained dependence ( ) [10] and from the single curve 75 k ≥ in Figure 3.For instance, the relative difference for k = 50 and 0.75 reaches 50% as Figure 4 shows.

k k ≥
and over the entire range of magnetic fields 0 1 b ≤ ≤ .In- stead I recommend using more accurate data from this paper, see Annex.The most accurate results are obtained through minimising the free energy density F of superconductor with respect to the variational parameters f ∞ and v ξ as described in this paper.This step is a must when aiming at the agreement better than 1% between magnetisation calculated symbolically and numerically, see Section 5.

The Error and the Correction
The error (present in Equations ( 10)-( 26) at 0 f ∞ → , see.Equations ( 10)-( 11)) is evident from the Abrikosov solution of the linearized GL equations [7] stating: for any k and 1 b → with A β equal to 1.1803 (1.1596) for the vortex lattice with square (hexagonal) unit cells [4] and to 1.1576 for the lattice with circular cells [9];  (29) From [ [7], Equation (19)] it is clear that at 3 1 b b ≤ ≤ the free energy density (see Equations ( 6)-( 8)) is (at minimum when A β is as close to 1 as possible): ( ) ( ) From Equations ( 30) and ( 16) (17) it is easy to check that for 1.1576A smooth transition from Equation (26) to Equation ( 27) can be achieved in several ways.In this paper I use the following approach.So far calculated from Equation ( 26 , so that the derivative of the magnetisation m 1 becomes equal to that set by Equation (28) and the transition from m 1 to m 4 is smooth since the higher derivatives are preserved.Furthermore, used for the comparison (in Section 5) magnetisation m 2 from [[3], Figure 7] is calculated numerically for the triangular vortex lattice, while the obtained here results are for the vortex lattice of circular cells, therefore accounting the respective ratio of the Abrikosov parameters (Equation ( 27)), more accurately: ( ) in this case.Clearly, more sophisticated methods of achieving the same result can be used (all required math for combining the solutions is present e.g., in [13]), but they are beyond the scope of this paper.
In conclusion, the correction makes the obtained solution compliant with [7] [13] (the same point m c2 and the same direction of the magnetisation curve at this point).It should be noted that this correction of the magnetisation uses the symbolic form of the theoretical result [7] at 1 b → and therefore the solution obtained here remains self-sufficient (even though it now uses two solutions of GL equations: [8] [9] [10] and [7]).So far I did not use any numerical results (the fact that calculated numerically magnetisation m 2 also agree with the conditions of Equations ( 27), (28) only means that these conditions are just).Cor-rected this way magnetisation m 1 is in excellent agreement with the conditions of Equations ( 27), (28) and is further compared to that calculated numerically (and with other symbolic methods) in Figures 5-8.

Comparison and Discussion
In Figures 5-8 the magnetisation m 1 is compared to that calculated numerically [3] and to those calculated symbolically [7]  × .An overview can be found elsewhere [3] [4] of the general features of reversible magnetisation, of its dependence on k and on magnetic field, etc. as follows from numerical solving of GL equations and these are not discussed here.Instead, I focus on the validation of the obtained analytical solution (Equations ( 20)-( 26) complemented by Equations ( 27)-( 30) and with the data in Figures 1-4 and Annex).This is achieved by comparing magnetization m 1 to magnetisation calculated using other methods.Namely, in Figures 5-8  where it is below 4%, see further elaborated in Section 6.1).where it is below 4%, see further elaborated in Section 6.1).As expected [3], the interpolation fit (m 3 ) totally fails to describe the data (m 2 ) quantitatively, since k < 3 and thus      where they are below 4% and 3.2% respectively, as further elaborated in Section 6.1).As expected [3], the interpolation fit (m 3 ) fails to de-scribe the data (m 2 ) quantitatively, since k < 3 and thus where it is below 2% respectively, as further elaborated in section 6.1).As expected [3] the interpolation fit (m 3 ) describes the data (m 2 ) quantitatively, since k > 3 and thus

Ginzburg-Landau Parameter k 200 =
The fragment of the magnetisation curve at k = 200 exemplifies that the error in the 1 st derivative of the magnetization m 5 is present at highest values of k and results in the noticeable difference where it is below 4% as elaborated in section 6.1).This result validates the advanced symbolic approach of this paper.On the other hand, for the first time numerical results [3] for hexagonal and circular unit cells are accurately validated over the entire ranges of k and m b in a trans- parent way with the (essentially independent) symbolic method of solving GL equations.Moreover, I find that the obtained close agreement of 1 m and 2 m makes some of the interpolation fits [ [3], Equations ( 15)-( 23)] obsolete.Clearly, the remaining discrepancy of 1 m and 2 m can be reduced by expanding the data sets in Annex (the simplification error is now about 0.5%), making the direct comparison of the underlying data, revising the symbolic solution and the boundary conditions, etc.This however is beyond the scope of this paper.It must be noted, that magnetisation calculated from Equations (20)-( 26) is rather sensitive to the errors in calculating the variational parameters.Therefore, in order to avoid the interpolation errors in calculating f ∞ and v ξ (that can be caused e.g., by splines), values of f ∞ , v ξ and m are best calculated at exactly the same values of b and k.
Presented in Annex data for the dependencies of f ∞ and v ξ on k and b contain the simplification error of up to 0.5% (caused by the limit in the size of the paper).I recommend that when aiming at the agreement better than 99% between magnetisation calculated symbolically and numerically, one has to use more accurate data obtained through the minimising the free energy density F of superconductor with respect to the variational parameters f ∞ and v ξ .Thus obviously "In order to achieve self-consistency in the theory the dependencies ( ) with 0.49693 The red dashed line is calculated from the symbolic expression for isolated flux line [5] [6]:

The Upper Critical Field bc2
As noted [8]  ≤ ≤ collapse practical- ly at the same curve.At lower magnetic fields there is a stratification depending on the k value, besides m 3 gives errors as expected [3] that are higher at lower k.
Note that in order to simplify comparisons with [8] [9] [10] in Tables A1-A4 the values of b are not corrected (by 0.985).Table A1.Spline data for the calculated variational parameter f ∞ at 0.75           Submit or recommend next manuscript to SCIRP and we will provide best service for you: Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.A wide selection of journals (inclusive of 9 subjects, more than 200 journals) Providing 24-hour high-quality service User-friendly online submission system Fair and swift peer-review system Efficient typesetting and proofreading procedure Display of the result of downloads and visits, as well as the number of cited articles Maximum dissemination of your research work Submit your manuscript at: http://papersubmission.scirp.org/Or contact jmp@scirp.org or B, defined externally (through ξ and BCS theory[2]) and at the end used as a scaling factor for all magnetic fields; 0 Φ is magnetic flux quantum.Note the local magnetic flux density ( ) ( ) 2 zk z c b r k B r B = = r b introduced (along with the four other local quantities) in

4 .
Here and below I use the dimensionless units: distance r, modulus of the magnetic vector potential ( ) a r ϕ , modulus of the super-velocity ( ) density (μ 0 is magnetic permeability of vacuum) and by the GL coefficients 0

Figure 1 .
Figure 1.Calculated variational parameter f ∞ vs. magnetic field b m , all dots from Annex for k ≥ 5 collapse at one curve approximated by the splines (visible here as the solid black line marked "k ≥ 5", the dashed black line is from Equation (18) and the solid red line is the cubic spline fit at k = 0.75, indicating the range of the deviation from the curve k ≥ 5).

2 f
∞ is independent on r.Replacing the variables: (25) after correcting the error (see Section 4.2) becomes the magnetization ( ) 1 m b − compared in Section 5 to those calculated numerically and with other symbolic methods.
examples in Section 5 will illustrate, at 1 b → the line 5 ([8] [9] [10] representing the magnetisation calculated from Equations (25) (26) without corrections) crosses the line 4 representing the m 4 and moreover it crosses the horizontal axis at 0.985 (instead of at 1 as Equation (27) implies).The error in ( ) 0 m b is visible at least at 0.5 1 a b ≤ ≤ .In this paper the error is eliminated in the following way.The Equations (21)-(26) are only valid at 0 f ∞ ≠ , see Equation (11).Otherwise, Equation (27) is the correct symbolic solution of GL equation(s) providing the missing additional conditions to Equations (21)-(26) and simply stating that in the plane ( ) 1 a m b − the point m c2 on the true magnetisation curve 1 m has the following coordinates: derivative through this point of the true magnetisation curve should be constant set by Equation (28) and that only depends on k and on the Abrikosov parameter 1 A β > [7].Moreover, using Equation (30) and min F (that follows from the minimising the free energy in Equations (16) (17)), we can now define more accurately the vague condition " 1 b → " (or "1 1 a b −  " [3]) as: ≤ ≤ with 3 b being defined by the condition: min 4 F F = .
) magnetisation 0 0.985 0 a b m = = .Thus Equation (26) produces the erroneous value of the upper critical field * 2 c b (see Section 6.2) due to ignoring the restriction 0 f ∞ ≠ , see the condition by Equations (10) and (11).Replacing b by 0.985 b (Equation (33)) when plotting ( ) m b satisfies the condition set by Equation (27) at 1 b → : 0 1 0 a b m = = .This step shifts the entire magnetisation curve (Equation (26)) parallel to itself and slightly to the right as exemplified in Section 5.Moreover, at 1 b → calculated from Equation (26) magneti- sation should have the same first derivative as set by Equation (28).This is obviously not the case as one can see in Section 5 (not only magnetisation curve calculated from [8] [9] [10] crosses the horizontal axis at 0.985 m b = , but it has the slope different from that set by Equation (27)).In this paper a compliance with this condition (Equation (28)) is achieved by introducing the correction: from Equation (26).This correction slightly rotates the entire magnetisation curve around the point ( ) [8] [9][10] through the entire range of magnetic fields 0 1 b ≤ ≤ and representative range of the GL parameter defined in Section 6.1).Moreover, in this range of k the 1 c b changes by more than 4 orders in magnitude: between 1 and 5 7.2 10 − I compare m 1 to m 2 -m 5 being respectively magnetisation calculated from: this work (m 1 ); [3] [4], numerically for hexagonal unit cells (m 2 ); [[3], Equation (19)] from the interpolation fit for m 2 with limited validity (m 3 ); [7], Equation (27) here (m 4 ) and [8] [9] [10] as published (m 5 ).As clear from Figures5-8, for any value of the GL parameter

Figure 1 ]
Figure 1]) set of the magnetisation curves at k = 0.85 is shown in Figure5(a).The relative difference

−
is too high.The data represented by magnetisation m 1 (and m 2 ) are in good agreement with these represented by m 4 at 0.8 1 b ≤ ≤ (so that both

Figure 6 .
Figure 6.(a) Calculated for k = 2 magnetisation m as function of the field

Figure 7 .
Figure 7. (a) Calculated for k = 10 magnetisation m as function of the field

Figure 8 .
Figure 8.(a) Calculated for k = 100 magnetisation m as function of the field and 1.3% respectively in the entire range 0 1 m b b ≤ = ≤ (except in the range: 0 0.01 m b ≤ ≤

Figure 1 ]−
Figure 1]) set of the magnetisation curves 2 k = and 5 k = is shown in Figure 6(a) and Figure 6(b) respectively.In Figure 6(b) the relative difference 1 2 2 m m m − is below 0.5% in the entire range of magnetic fields 0 1 m b b ≤ = ≤

−Figure 6 − 4 at−
Figure 6(b) and Figure 7(a) respectively.In Figure 7(a) the relative difference 1 2 2 m m m − is below 1% in the entire range 0 1 m b b ≤ = ≤ (except in the range: of the magnetisation m 2 caused by the digitalisation of the data [[3], Figure 7] is also visible in the figure.To summarise, over the entire ranges of the GL parameter 12

6 . The Critical Fields 6 . 1 .
respect to v ξ and f ∞ "[8] [9][10], to which I add: the minimisation procedure is straightforward, use Equations (6) (7) (8), (16) (17) and a standard optimisation program, such as commonly available Solver in Excel.I find the quote above missing the "how" part (almost as to write: in order to get solution of GL equations, solve them).The Lower Critical Field bc1 In this paper I calculate values of the field b c1 from Equations (25) (26) typically at 0 b = the Equation (25) diverges.In Figure 9 the boxes show calculated this way dependence ( ) 11 c b k of the lower critical field b c1 .The black solid line is calculated from the fit to the numerical results [3] [4]:

[ 9 ]Figure 10 .
Figure 10.Minimum value F min of the free energy density F (Equation (16) (17)) as function of the resulting magnetic field b m for the selected values of k (shown next to the corresponding lines).

Figure 11 .∞
Figure 11.Calculated for k = 34, 50 and 75 magnetisation as function of the field

Table A2 .
Spline data for the calculated variational parameter

Table A3 .
Spline data for the calculated variational parameter

Table A5 .
The dimensional and dimensionless quantities and scaling factors.