Numerical study of the stability of solutions for the half-space Ginzburg–Landau model

Based on a shooting alternative that allows one to numerically solve the one-dimensional system of Ginzburg–Landau in an unbounded domain, a numerical study of the stability of solutions of this system is performed here. This stability notion, from a physical point of view, means that each solution of the system is identified as stable when it minimizes the corresponding Ginzburg–Landau functional. As opposed to a previous paper, the present one is concerned with a more general study since the weak and large regimes of the Ginzburg–Landau parameter are considered and the initial data are no longer subject to the de Gennes condition. Certain conjectures regarding the superheating field are also investigated numerically.


Introduction
The Ginzburg-Landau theory allows one to describe the states of a superconducting material in an exterior magnetic field H e . When the material is a film delimited by two planes (x = −d, x = d), with the exterior field parallel to the surface, the determination of states from this theory involves a functional depending on the inner magnetic potential A and on a real-valued wave function f (cf. e.g. [1][2][3][4][5]). In fact, the states of the film are then characterized from global or local minima of this functional. When we are concerned with pairs (f, A), with f even and A odd on [−d, d], also called symmetric solutions, we thus consider the reduced Ginzburg-Landau functional, where d is proportional to the thickness of the film and h is proportional to |H e |. In the definition of ε d , it is assumed that the edge of the film is at x = 0. Also, (f, A) ∈ (H 1 ((0, d))) 2 with A(d) = 0. The so-called normal phase corresponds to the case where f ≡ 0 with A ′ ≡ h, whereas the superconducting phase corresponds to the context where f ≡ 0. The Ginzburg-Landau parameter, κ > 0, characterizes the properties of the material. Its values determine the type of superconductor according to the type * LAMFA, CNRS, Université de Picardie, 33 r. Saint-Leu, 80039 Amiens Cedex 1, France. Email: Pierre.Delcastillo@upicardie.fr † IECN, CNRS, Université de Nancy-I, 54506 Vandoeuvre-lès-Nancy Cedex, France. Email: Seraphin.Mefire@iecn.unancy.fr of phase transition which takes place between the normal phase and the superconducting phase: κ small describes what is known as a type I superconductor and κ large as a type II superconductor.
The critical value of κ that is usually given to separate type I and type II superconductors is κ = 1 √ 2 . For a type I superconductor, there is a critical value h ⋆ c such that if h < h ⋆ c , the material is entirely superconducting and the magnetic field is expelled from the sample away from a boundary layer. This is called the Meissner effect. If h > h ⋆ c , superconductivity is destroyed and the material is in the normal state. Regarding applications, there have been proposals to use superheated type I superconductors as detectors for elementary particles, so that the sample acts as a superconducting "bubble chamber" (cf. e.g. [6]). The passage of a sufficiently energetic particle through the sample would initiate the transition to the normal state. For a type II superconductor, the phase transition is different and there are two critical values: h c 1 and h c 2 . For h < h c 1 , the exterior magnetic field is expelled from the sample and there is a Meissner effect as for type I superconductors. When h increases above h c 1 , superconductivity is not destroyed straight away, since the superconducting and the normal phases coexist in the form of filaments or vortices. As h increases further, the vortices become more numerous until the critical value h c 2 is reached at which superconductivity is destroyed. For h > h c 2 , the material is in the normal state. The way superconductivity is nucleated is highly dependent on d and κ (see e.g. [7,Chap. 3]). Further discussions of the development of the above model as well as for the non reduced model, and physical significances, are given in [2], [8][9][10][11][12].
When the width d of the film is large (in the sense that κd is large), we are concerned with a functional, ε ∞ , formally introduced by setting d = +∞ in the definition of ε d after a renormalization obtained by adding the term ( 1 2 − h 2 )d. In fact, (1 − f ) ∈ H 1 ((0, +∞)), A ∈ H 1 ((0, +∞))}. The corresponding Ginzburg-Landau equations, expressing the necessary conditions for the existence of minima, are then with the initial boundary conditions and (1 − f, A) ∈ (H 2 ((0, +∞))) 2 . The problem (1) -(2) is called the half-space model and was studied in [3,4]. Following [13], any pair (f, A) satisfying also (1) - (2) is such that f > 0 and is thus called a positive solution. On the other hand, f ≤ 1. With the choice of the variational space H ∞ , the normal solutions are eliminated and a particular role is attributed to the positive solutions since it is then enforced that f tends to 1, as x tends to +∞. For any fixed κ, we know (see [2]) that the set of h, such that there exists a solution of (1) - (2), is a bounded interval [0, h + ]. We then introduce the superheating field, h sh (κ), defined as the supremum of this interval. This critical field is very important for many applications. For instance, measuring the superheating field provides one of the few methods for experimentally determining the characteristic of the superconductor (see [14,Chap. 16] for other properties and developments). The Ginzburg-Landau equations corresponding to ε d constitute a system, denoted by (GL) s d , introduced similarly as (GL) ∞ but on (0, d) and with in addition to (2): f ′ (d) = 0, A(d) = 0. Let us recall that Bolley, Foucher and Helffer have shown in [1,15] that the system (1) -(2) is close to (GL) s d for κ small, with κd large. More precisely, they have proved the existence of a local superheating field for (GL) s d which is near to the global superheating field h sh (κ) for (1) -(2) when κd is large.
In [16], Bolley and Del Castillo have established the following existence result, regarding (GL) ∞ . Theorem 1.1 For all κ > 0 and for all f 0 ∈ (0, 1], there exists a unique solution (f, A) of (1) such that: We can then introduce (see also [17]), for κ > 0, called the response map. Hereafter, we denote a solution of (1) -(2) by (f, A; h); the same notation will represent also a solution of (1) and (3), according to the response map.
In [18], we numerically studied the stability of solutions of (1) -(2), in the weak−κ regime and in the context of a localization of these solutions. Precisely, we were concerned with pairs (f, A) satisfying (1) - (2) and in accordance with the de Gennes condition (see [19]), . The case f (0) = 1 √ 2 was taken into account in the study of the stability when the intensity h, of the exterior magnetic field, was considered near the superheating field.
Based on a formal approach, namely the method of matched asymptotic expansions, Dolgert, Di Bartolo and Dorsey [20] introduce in the weak−κ regime an expansion of h sh (κ) in powers of κ 1 2 . It results from this paper that the numerical values of h sh (κ), obtained from a computation of solutions by a Newton-type method, agree quite well with the formal values given by the introduced asymptotic expansion. These authors have studied moreover the stability of solutions of (1) -(2) in the weak−κ regime, with the same formal approach.
As opposed to [18], here we numerically study the stability of solutions of (1) -(2): in the weak and large−κ regimes, in the context where the pairs (f (0), h) are not subject to the de Gennes condition, and for f (0) ∈ (0, 1).
In contrast with the approach of [20], we are here concerned with a numerical framework that leads us both to a determination of the superheating field and to a study of the stability of solutions in each regime of κ. This framework requires the combination of a shooting approach (introduced rigorously in [17]) with a semi-implicit method (which is A−stable) for the numerical computation of solutions, as well as Hermite element approximations in the stage of the stability study. The efficiency of this framework is distinguished by the fact that small values of the initial datum f 0 = f (0) can be here considered, contrary with an approach using the method of matched asymptotic expansions as in [20], where the formal pair, from which the expansion of h sh (κ) is determined, cannot be an approximate solution of (1) -(3) for such an initial value. The robustness of the present framework is noticed since it also allows us to study the stability of solutions between the two regimes of κ.
This paper is subdivided into five sections. In Section 2, we are concerned with some properties of the solutions of (GL) ∞ , and of σ κ . After describing in Section 3 the numerical approximation of solutions of (1) -(2), based on a shooting alternative, we numerically inspect σ κ and other maps associated with (GL) ∞ , as well as certain conjectures regarding the superheating field. In Section 4, we numerically study the stability of the solutions of (1) -(2) in the weak and large−κ regimes as well as between these regimes. Our numerical investigations lead us mainly to conjecture that the particular initial value f 0 , from which h sh (κ) is determined, corresponds to the initial datum for which the change of stability of the associated solution occurs -independently of the values considered here for κ. Our conclusions are given in Section 5.

Preliminaries
Let us start by recalling some general properties (see [2,21]) of the solutions of (GL) ∞ . Proposition 2.1 Let (f, A) be a solution of (1) - (2). The pair (f, A) is such that: • f is increasing on [0, +∞) and we have, for all x ∈ [0, +∞), • A is strictly negative and increasing on [0, +∞); we have • the following energy conservation law, is satisfied.

Proposition 2.2
We have: In the weak−κ regime, we have In the large−κ regime, we have Remark 2.1 Considering (5) at the point x = 0, and according to the definition of σ κ , it results that From Theorem 1.1, we can introduce the following maps where f (x; f 0 ) and A(x; f 0 ) represent the values of the solution (f, A; h) at the point x, and for f 0 ∈ (0, 1] given. Based on (6), we can show that these maps have the following properties.

Proposition 2.3 We have
, as follows: where the superscript "T " denotes the transpose.

Numerical inspections of certain maps and conjectures
In this section, we numerically study the behavior of the maps F, A introduced in (9), and inspect certain conjectures regarding the superheating field. We are thus concerned at first with the numerical approximation of the solutions of (GL) ∞ associated with the initial conditions in (2).

Shooting alternative and numerical approximations
A way to overcome the numerical difficulties due to the conditions at infinity consists of approximating the solutions of (1) -(2) from computations involving an associated initial value problem. Such a problem is obtained by replacing the conditions at +∞ involving f and A, by suitable initial conditions for f and A. It consists of finding f and A, for (f 0 , A 0 ; h) given, such that: where I max := [0, T max ) is the maximal interval of existence. As any pair (f, A), solution of (1) -(2), satisfies the energy conservation law (5), we assume that As indicated also in Proposition 2.1, the solutions of (1) -(2) are defined on (0, +∞) and bounded.
This is not the case for most of the solutions of (GL) ivp . In order to select the trajectories of (GL) ivp that correspond to those of (1) - (2) and (13), we make use of a shooting alternative, suitable in all the regimes of κ. This alternative is introduced from the definition of two disjoint open sets , and is formulated as follows (see Bolley and Helffer [17]): (13) is such that f is increasing and becomes larger than 1. Moreover, the function A is strictly negative on I max and becomes decreasing. (13) has the property that there exist

Moreover, A is increasing on
I max and becomes greater than 0.
Based on this alternative, we numerically determine the particular value σ κ (f 0 ) of h which leads to a solution of (1) -(2), for f 0 and A 0 considered as in (13). More precisely, after fixing f 0 ∈ (0, 1) at first, for each h, we take A 0 in accordance with (13). We start at x = 0 with the initial data of (GL) ivp , and compute the functions f , f ′ , A and A ′ by iteration on an increasing sequence (x k ) k , Runge-Kutta like method (with a variable step of size ζ i ). By testing on such a sequence the criteria on f , f ′ , A and A ′ provided from the contexts i) and ii) of Theorem 3.1, we can deduce whether the considered h is above or below the value σ κ (f 0 ). We mention here that our simulations reveal that only criteria on f , f ′ can be involved in the large−κ regime and only criteria on A, A ′ can be enforced in the weak−κ regime. Concretely, for . This dichotomy process leads to two values h n and h n+1 , h n < h n+1 , for some n, such that: Before presenting the numerical method used for solving (GL) ivp , let us reformulate this system as follows: where Each solution (f, A) of (1) -(2) will be determined from a solution Y of (14), (13), with lim x→+∞ Y 1 (x) = 1 and lim x→+∞ Y 3 (x) = 0. As in [25], in the context of a numerical study of the superheating field in the weak−κ regime, we use the semi-implicit Runge-Kutta method of order 3, with two intermediate steps, described below.
This is an A-stable method (see [26,Chap. 5]), which appears efficient for solving stiff problems. We are here concerned of course with an initial value problem which is very sensitive to small variations of the initial data, since we are looking for trajectories tending to the unstable stationary solution The numerical approximation of the solution Y is achieved on an ordered sequence of points x i , and by H i the variable steplength of the mentioned sequence, we determine the approximation Y i+1 of Y (x i+1 ) by solving the following system: where The parameter H i is determined at each step from comparative calculations involving two different steplengths and a local a posteriori estimate (see e.g. [27, Chap. 11] for such calculations). On the other hand, we use a Newton-type method for solving (15).
The present approach for approximating a solution of (1) -(2) is then based on the combination of the shooting alternative with this numerical method.

Numerical inspections
From the computation of each pair (f, A) described as above, we numerically inspect here at first the behavior of F(x, f 0 ) and of A(x, f 0 ), with respect to f 0 ∈ (0, 1) as well as to x ∈ IR + , in the weak and large−κ regimes.     Let us now inspect the behavior of the response map σ κ . The results in the weak−κ regime are represented in Figure 5 whereas those obtained in the large−κ regime are shown by Figure 6.
In the simulations, the sequence of values considered for f 0 is enriched, around its element providing a coarse evaluation of the superheating field, in order to determine with more accuracy the particular initial value that is associated with the superheating field.  Before inspecting the conjectures regarding h sh , let us examine how our numerical computations fit with the evaluation of h sh based on the Parr formula [28], which is recalled in (7) and rigorously proved in [23,24]. In this way, we reproduce in the second column of Table 1 the values provided by the first two terms of the asymptotic expansion in (7) (7) indicates that results fit well with the Parr formula, and validates in a certain sense our computations. For a more detailed numerical inspection of the Parr formula, we refer to [25].
We are now interested in the numerical inspection of two conjectures: an asymptotic expansion of h sh proposed by Chapman [30], and an estimate (cf. [22]) that involves also h sh .
The expansion of the superheating field conjectured in [30] concerns the large−κ regime and is expressed as follows: where C ≃ 0.33.
Let us mention that the existence of a family of solutions of the second Painlevé equation which arose in the approach of [30] is investigated in [31], as well as some properties associated with the  Table 2: Values of f sh 0 and of f sh 0,num in the weak−κ regime (f sh 0,num being associated with h sh num presented in Table 1).
existence of the superheating field.
In order to numerically inspect formula (17), we reproduce in the second column of Table 3  The numerical results presented in Tables 1 and 3, as well as those obtained with the intermediary values, κ = 0.6, 0.7, 0.707, 0.8, 0.9, that are respectively h sh num = 1.023937, 0.982080, 0.979540, 0.949568, 0.923575, suggest that the following estimate involving the superheating field, and conjectured in [22], appears also to be numerically satisfied.
By distinguishing also by f sh 0,num , in Table 4 and for the large−κ regime, the particular values of f 0 deriving from numerical computations and associated with h sh num presented in Table 3, it appears that these values become (very slowly) smaller when κ increases. Also, it results (see e.g. Figure 6) that σ κ is increasing on (0, f sh 0,num ] and decreasing on [f sh 0,num , 1). The results of Tables 3 -4 and of Figure 6 lead us to conjecture in the large−κ regime that there exists a mapθ, defined on (0, +∞), withθ(κ) decreasing when κ becomes large, such that    Table 4: Values of f sh 0,num in the large−κ regime (associated with h sh num presented in Table 3).

Numerical study of the stability of solutions of (GL) ∞
In [13], Bolley and Helffer proved that the functional ε ∞ is semi-bounded when h < 1 Our investigations, based on the numerical approach of Subsection 3.1 and on a framework involving a generalized eigenvalue problem (see also [18]), will allow us to provide numerical answers to this question.

Generalized eigenvalue problem
Let us start by introducing some notation. For f 0 ∈ (0, 1], let (f, A; h) be the unique solution of (1) -(3) obtained from Theorem 1.1. Let us set: We have studied in [18] the operator L f,A , defined on To check if the pair (f, A) is a local minimum of ε ∞ , we have just to analyze the Hessian of ε ∞ at (f, A; h). This Hessian is defined on This particular value of f 0 will thus be excluded in our numerical investigations.
Let us denote by . , . the usual scalar product on (L 2 (IR + )) 2 and by . the associated norm.
Our approach is here reduced to the numerical inspection of this spectrum, σ dis (L f,A ), for studying then the stability of the pair (f, A). This approach involves, as in the framework developed in [18], is defined on In the present approach, we shall numerically inspect then the bottom of the discrete spectrum of L f,A by considering the weak eigenvalue problem associated with Lf ,Ã , and that consists of finding Hereafter, the usual scalar product on (L 2 ([0, d])) 2 is represented by ( . , . ) (L 2 ([0,d])) 2 . Also, the (usual) norm associated with any functional space X is denoted by . X below. the injection ofH in H is compact. Since there exist constants γ > 0 and α > 0 such that a(U, U ) + the problem that consists of finding a pair (λ, U ) ∈ IR ×H, U = (u 1 , u 2 ) T , such that: admits at least one solution. In fact, from [33, p. 137], there exists (λ n , U n ) ∈ IR×H, U n = (u n 1 , u n 2 ) T , with (λ n ) n an increasing sequence, such that a(U n , V ) = λ n (U n , V ) (L 2 ([0,d])) 2 ∀ V ∈H. Taking On the other hand, it results from standard arguments that U n ∈D and satisfies: Lf ,Ã (U n ) = λ n U n .
Following Remark 4.1, the weak problem (21) admits at least one solution.
Remark 4.2 From (18) - (19), there exists a real constant C such that: The operator L f,A being self-adjoint and its essential spectrum being equal to [1, +∞), from [32, p. 120], the lowest eigenvalue of L f,A , denoted by α ∞ (f 0 ), can then be expressed as follows: Representing similarly the lowest eigenvalue of Lf ,Ã by α d (f 0 ), it follows that for any ǫ > 0, there As in [18], we also consider here Hermite finite elements (see e.g. [34,Chap. 3]) that appear well suited to an accurate discretization of the space D. Let us denote by {x j } 1≤j≤N +1 the set of where P 3 is the space of polynomials of degree less than or equal to 3. The discrete space associated with H 2 00 ([0, d]) is defined as follows: The discrete problem associated with (21) consists of finding (α k , We are thus concerned with (22) in numerical simulations, where A and f are respectively replaced by their numerical approximations A app and f app , obtained from the numerical approach of Subsection 3.1. As in [18], we consider the usual Jacobi method (see e.g. [34,Chap. 6]) for the computation of eigenvalues, based on the linear matrix system deriving from (22), where an adaptive Simpson formula is applied for calculating the integral terms.

Numerical results
The numerical study of the stability of the solutions of (GL) ∞ is here performed through the numerical inspection of the map φ k,d : is the lowest eigenvalue resulting from (22).
We first investigate the weak−κ regime. In the experiments associated with the results of Figure   7, we consider κ = 0.1, and d = 1000.   Table 5: Values of f eig 0,num in the weak−κ regime.
A study of the influence of the mesh size k on the computation of the minimal eigenvalue of the discrete spectrum has also been performed. By considering now β k 1 ,k 2 Namely, for each considered value of κ, it appears that there exists a critical value f eig 0,num such that The values obtained from simulations, with respect to κ, and corresponding to f eig 0,num , are presented in Table 5.
The numerical results regarding the stability of the solutions of (GL) ∞ are then here in accordance with those obtained in a formal way by Dolgert et al. [20] in the weak−κ regime.
We are now concerned with the experiments in the large−κ regime. Figure 9 presents the results obtained with κ = 1.3 and d = 300. It appears for these considerations that there exists a critical a value f eig 0,num such that for f 0 > f eig 0,num , we have φ k,d (f 0 ) > 0. A study, in the large−κ regime, on the influence of the bound of the truncated domain on the computation of the minimal eigenvalue of the discrete spectrum, has been performed. We observe that the agreement is quite good with d = 300; for instance, β k 1000,300 is smaller than 10 −7 , when κ = 1.3. We thus fix d = 300 for the following experiments. Those associated with Figure 10 Table 6 reproduces the values obtained from simulations, with respect to κ, and corresponding to  Table 6: Values of f eig 0,num in the large−κ regime.  Table 7: Values of f sh 0,num and f eig 0,num in the intermediary−κ regime.
Since the same remark was observed in the weak−κ regime, we are then led to conjecture that: in each regime of κ, σ κ (f 0 ) = h sh (κ) if and only if φ k,d (f 0 ) = 0.
The results presented in Table 7 suggest that this remark appears also to be satisfied in the context of intermediary values of κ.
The values of f eig 0,num become, very slowly, smaller when κ increases; the consideration of κ = 21 for example provides f eig 0,num = 0.297935 which is not significantly smaller than the value that was obtained for κ = 15. Our results suggest in the large−κ regime that there existκ 0 > 0 and a map θ decreasing on [κ 0 , +∞), such that for all f 0 ≤θ(κ) the solution (f, A; h) is unstable, and for all f 0 >θ(κ) the solution (f, A; h) is stable.

Conclusions
A numerical inspection of certain maps and conjectures associated with the solutions of the half-space Ginzburg-Landau system has been performed in the weak and large−κ regimes. These conjectures involve the superheating field h sh ; namely, the asymptotic expansion of h sh proposed by Chapman [30], in the large−κ regime, as well as the estimate (cf. [22]) involving h sh in all the regimes of κ appear numerically satisfied. On the other hand, we have numerically studied the stability of solutions of this system in all the regimes of κ. In the weak−κ regime, our numerical results suggest that there exist κ 0 > 0 and a map θ defined on [0, κ 0 ], with θ(κ) ∈ (0, 1), θ(0) = 1 √ 2 , such that for all f 0 ≤ θ(κ) the solution (f, A; h) is unstable, and for all f 0 > θ(κ) the solution (f, A; h) is stable.
Typically, these results extend those that we had obtained in [18] from the framework based on a localization of solutions of the system (GL) ∞ satisfying the de Gennes condition. In the large−κ regime, our numerical results suggest that there existκ 0 > 0 and a mapθ decreasing on [κ 0 , +∞), withθ(κ) ∈ (0, 1), such that for all f 0 ≤θ(κ) the solution (f, A; h) is unstable, and for all f 0 >θ(κ) the solution (f, A; h) is stable.
Our numerical investigations lead us to conjecture that σ κ (f 0,c ) = h sh (κ) if and only if 0 is the lowest eigenvalue of L f,A , where f (0) = f 0,c . In other words, such a value f 0,c corresponds to the initial datum for which the change of stability occurs.
The results obtained here provide the numerical answers to several questions that remain open from a theoretical point of view and concern therefore the stability of solutions of the half-space Ginzburg-Landau system as well as the study of the two conjectures involving (in the large−κ regime and in all the regimes of κ respectively) the field h sh . The theoretical proof or a numerical inspection of the asymptotic expansion of h sh in powers of κ 1 2 , introduced in [29] formally, at any order, and extending thus in the weak−κ regime the Parr formula (see (7)), remains also to be performed. Another numerical work is the one that would investigate, in the large−κ regime, the extension of the asymptotic expansion of h sh proposed by Chapman (see (17)), by describing then the associated third term.