Harmonic oscillators with Neumann condition on the half-line

This paper is devoted to the computation of the minimum Θ0 of the first eigenvalues for the Neumann realization of harmonic oscillators on the half-line. We propose an algorithm to determine this minimum and we estimate the accuracy of these computations. We also give numerical computations of constants appearing in superconductivity theory.

IRMAR, ENS Cachan Bretagne, Univ. Rennes 1, CNRS, UEB av Robert Schuman, F-35170 Bruz, France (Communicated by the associate editor name) Abstract. We consider the spectrum of the family of one-dimensional selfadjoint operators −d 2 /dt 2 + (t − ζ) 2 , ζ ∈ R on the half-line with Neumann boundary condition. It is well known that the first eigenvalue µ(ζ) of this family of harmonic oscillators has a unique minimum when ζ ∈ R. This paper is devoted to the accurate computations of this minimum Θ 0 and Φ(0) where Φ is the associated positive normalized eigenfunction. We propose an algorithm based on finite element method to determine this minimum and we give a sharp estimate of the numerical accuracy. We compare these results with a finite element method.

Motivation .
Even if an accurate computation of the spectral quantities of the family of harmonic oscillators H(ζ), ζ ∈ R, can be interessant by itself, the parameters Θ 0 and Φ(0) appear frequently in the analysis of the superconductivity modelled by Ginzburg-Landau theory. Let us recall now some results in this way.
Let Ω ⊂ R 2 be a bounded simply-connected domain with Lipschitz boundary. The Ginzburg-Landau functional involves a wave function ψ and a magnetic potential A and reads for any (ψ, A) ∈ W 1,2 (Ω; C)×{A = A 0 +Ã withÃ ∈Ḣ 1 (R 2 , R 2 ), divÃ = 0}, where A 0 (x) = 1/2(−x 2 , x 1 ). We use the notationḢ 1 (R 2 ) for the homogeneous Sobolev spaces. The function ψ is called the order parameter and its modulus informs us about the localization of the superconductivity. The parameter κ is a characteristic of the material and we only consider superconductors of type II for which κ is large. The physical interessant states are the minimizers of the functional E κ,H and they satisfy the Euler-Lagrange equation. When H is very large, the unique minimizer is the trivial state (0, A 0 ), called the normal state, and there is no superconduting property for the sample in this case. Decreasing the applied magnetic field, the sample becomes superconducting. Thus we define the critical field H C3 as the value of H where the transition between the normal and superconducting state takes place: The first rigorous definition of the critical field H C3 appeared in [28]. The calculation of this critical field H C3 for large values of κ has been the focus of much activity (see [22,2,27,28,29,25,14,15,16]). In the works [14,15,16,17], the definition of H C3 in the case of samples with smooth section has been clarified and the asymptotic is given by: Proposition 2 (see [16]). Suppose Ω is a bounded simply-connected domain in R 2 with smooth boundary. Let κ max be the maximal curvature of ∂Ω. Then It was realized that the asymptotics of the critical field is completely determined by the linear eigenvalue problem. Indeed, if we denote by µ (n) (h) the n-th eigenvalue of the magnetic Neumann operator P h = (ih∇ + A 0 ) 2 defined on D(P h ) = {u ∈ H 2 (Ω)|ν · (ih∇ + A 0 )u |∂Ω = 0}, then the asymptotics of µ (n) (h) was established by Fournais-Helffer in [15]: Proposition 3 (see [15]). Suppose that Ω is a smooth bounded and simply connected domain of R 2 , that the curvature ∂Ω s → κ(s) at the boundary has a unique maximum κ max reached at s = s 0 and that the maximum is non-degenerate, i. e. k 2 := −κ (s 0 ) = 0. Then for all n ∈ N, there exists a sequence {ξ (n) j } ∞ j=1 ⊂ R such that µ (n) (h) admits the following asymptotic expansion (for h → 0): HARMONIC OSCILLATORS ON THE HALF-LINE   3 To carry through an analysis of the critical field H C3 in the case of domains with corners, a linear spectral problem, studied in depth in [5,6,7,8], is usefull. Let us first give estimates for the Schrödinger operator in a model geometry: the infinite sector.
Proposition 4 (see [6]). Let G α be the sector in R 2 with opening α and Q α be the Neumann realization of the Schrödinger operator (i∇ + A 0 ) 2 on G α . We denote by µ k (α) the k-th smallest element of the spectrum given by the max-min principle. Then: 1. The infimum of the essential spectrum of Q α is equal to Θ 0 . 2. For all α ∈ (0, π/2], µ 1 (α) < Θ 0 and µ 1 (π) = Θ 0 . 3. Let α ∈ (0, 2π), k ≥ 1 be such that µ k (α) < Θ 0 and Ψ α k an associated normalized eigenfunction. Then Ψ α k satisfies the following exponential decay estimate: Thanks to the model situation given by the analysis of the angular sector, we are able to determine the asymptotic expansion of the low-lying eigenmodes of the Schrödinger operator on curvilinear polygons: Proposition 5 (see [7]). Let Ω be a bounded curvilinear polygon, Σ be the set of its vertices, α s be the angle at the vortex s. We denote by Λ n the n-th eigenvalue of the model operator ⊕ s∈Σ Q αs , and µ (n) (h) the n-th smallest eigenvalue of P h . Let n be such that Λ n < Θ 0 . There exists (m (n) j ) j≥1 such that µ (n) (h) admits the following asymptotic expansion (for h → 0): If Ω is a bounded convex polygon, there exists r n > 0 and for any ε > 0, C ε > 0 such that For non constant magnetic field, the low-lying eigenvalues admit an asymptotic expansion in power of √ h. These results highlight the importance of comparing µ k (α) with Θ 0 and then of computing precisely Θ 0 . It is also natural to wonder for which angle α we have µ k (α) < Θ 0 . It was conjectured in [1,8] that µ 1 is strictly increasing from (0, π) onto (0, Θ 0 ) and is equal to Θ 0 on [π, 2π). This conjecture is based on numerical computations and could be strengthened with an accurate estimate of Θ 0 . As in the case of smooth domains, spectral informations produce results about the minimizers of the Ginzburg-Landau functional for domains with corners. We obtain in particular a complete asymptotics of H C3 for large values of κ in terms of linear spectral data and precise estimates on the location of nucleation of superconductivity for magnetic field strengths just below the critical field: Proposition 6 (see [10]). Let Ω be a curvilinear polygon and Λ 1 = min s∈Σ µ 1 (α s ). There exists a real-valued sequence {η j } ∞ j=1 such that Let µ ∈ (Λ 1 , Θ 0 ) and define Σ = {s ∈ Σ|µ 1 (α) ≤ µ}. There exist constants κ 0 , M , This Agmon type estimate describes how superconductivity can nucleate successively in the corners, ordered according to their spectral parameter µ 1 (α s ) seeing that µ 1 (α s ) < Θ 0 . This reinforces the interest to compare precisely µ 1 (α) and Θ 0 .
When we consider the Schrödinger operator in dimension 3, see [23,24,30], we have to analyze some new operators: the Neumann realization of is the angle that makes the magnetic field with the boundary at each point (approximated by the tangent plane). We first make a Fourier transform in r. When θ = 0, we are led to the so-called de Gennes operator H(ζ) on the half-line (see [12] and this present paper). If θ = 0, we perform a translation in s and a rescaling. Thus we are reduced to a Schrödinger operator with an electric potential on the half-plane This operator is deeply studied in [9], both theoretically and numerically. The authors prove an isotropic estimate and anisotropic estimate for the eigenfunctions. They also analyze the asymptotics when θ → 0. In particular, they prove the following result: Proposition 7. We have the following upper-bound for the n-th eigenvalue σ n (θ) of L θ : σ n (θ) ≤ Θ 0 cos θ + (2n − 1) sin θ, ∀n ≥ 1.
For all n ≥ 1, there exists a sequence (β (n) j ) j≥0 such that σ n (θ) is an eigenvalue for θ small enough and admits the following asymptotic expansion (for θ → 0): If we denote by n(θ) the number of eigenvalues of L θ below the essential spectrum, we have with (1): If we bound from below Θ 0 by 1, 0.6, 0.591, we lower-bound shows that n(π/2000) is greater than 0, 127 and 130 respectively. A greater approximation of Θ 0 we have, a greater lower-bound of n(θ) we deduce.
1.3. Main results. An estimate of Θ 0 by 0.59010 was already given in [13], using the Weber functions but there is no mention of the accuracy of this estimate. Using an integral representation [11], Chapman approximates Θ 0 by 0.59 without any estimate of the error. In the literature, we can find some estimates of Θ 0 but there is no mention of the accuracy of the computations. To our knowledge, we do not find any computation of Φ(0). The aim of this article is to give accurate estimates of Θ 0 and Φ(0) and of the error between exact values and numerical computations. The numerical method implemented here is very standard since we use finite difference and finite element methods.
To estimate the numerical accuracy, we first establish in Section 2 error estimates on eigenmodes: Theorem 2.1 quantifies the gap between the eigenvalue Θ 0 and the energy associated with a quasi-mode for the operator H(ζ). In Theorem 2.2, we prove H 1 -estimate between the normalized eigenfunction Φ associated with Θ 0 for the operator H(ζ 0 ) and a normalized quasi-mode for H(ζ). We deduce in Theorem 2.3 an estimate of Φ(0). In Section 3, we construct an adequate quasimode combining the finite difference method and analysis of the ODE theory for the differential equations depending on parameters. We implement this method in Subsection 3.5 and obtain an accurate approximation of Θ 0 and Φ(0): Section 4 presents computations with the finite element method. From a numerical point of view, we also mention papers [4,3] which deal with the numerical computations for the bottom of the spectrum of −d 2 /dt 2 + (t − ζ) 2 on a symmetric interval using a finite difference method.
2. Error estimates on eigenmodes. This section concerns the analysis of the operator H(ζ) and error estimates between Θ 0 and the energy associated with a quasi-mode for H(ζ).
Theorem 2.2. Let ζ ∈ R and ϕ ζ be a normalized and positive function of D. With Notation 2, we assume η ζ ≤ µ 2 (ζ 0 ). Then To prove this result, we use an estimate of quasi-modes established in [21, Proposition 4.1.1, p. 30] : where λ min S is the smallest eigenvalues of S = ( Ψ j , Ψ k H ) and d the non-symmetric distance defined by d( Proof. Theorem 2.2 We apply Proposition 8 with N = 1, A = H(ζ 0 ), Ψ 1 = ϕ ζ , E the space spanned by ϕ ζ and F the space spanned by Φ. We first connect the distance d with the norm ϕ ζ − Φ L 2 (R + ) by noticing that we estimate r ζ L 2 (R + ) using the orthogonality relation r ζ , ϕ ζ = 0: Relations (9), (10) and Let us now estimate the L 2 -norm of (ϕ ζ − Φ ). An integration by parts gives: On the other hand, We deduce from (11), (12) and Theorem 2.1 a upper-bound for the L 2 -norm of Φ − ϕ ζ : We deduce now an estimate for ϕ ζ − Φ at point t = 0.
Theorem 2.3. Using the same notation and assumptions as Theorem 2.2, we have Proof. As Φ − ϕ ∈ H 1 (R + ), it suffices to write We conclude with the Cauchy-Schwarz inequality.
3. Construction of a quasi-mode by a finite difference method. Theorem 2.1 gives bounds for Θ 0 as soon as we get quasi-modes for the operator H(ζ).
Of course, the closer ζ is from ζ 0 , the better the bounds. A heuristic approach based on finite difference method and the ODE theory gives a sequence of approximated values for ϕ ζ . Then we use this sequence to construct a test-function with energy as small as possible and thus try and give a good approximation of Θ 0 . We organize this approximation in several steps: For L large enough, the function µ N,N (ζ, ·) is increasing on (L, +∞) and Proof with u ζ,L a normalized eigenfunction associated with µ N,N (ζ, L). The positivity of the first derivative is directly deduced for L large enough.

3.2.
Finite difference scheme. Instead of looking for a normalized eigenfunction, we impose the value of Φ at t = 0. Therefore, we try to determine (ζ 0 , Φ) ∈ R + × D such that: Varying parameter ζ 0 , it is natural to look for a function ϕ ζ and satisfying: The system (18) is numerically solved by a finite difference scheme. Let h be step of discretization. We determine recursively an approximationφ ζ j of ϕ ζ (jh) for any integer j ≥ 0. For this, ϕ ζ (jh) and ϕ ζ (0) are classically approximated respectively by (φ ζ j+1 − 2φ ζ j +φ ζ j−1 )/h 2 and (φ ζ 1 −φ ζ 0 )/h. The boundary condition at t = 0 determines completely the sequence (φ ζ j ) j≥0 : 3.3. Dependence on ζ of the sequence (φ ζ j ). The change of variables x = t − ζ in the eigenmode equation leads to the second order differential equation: The Sturm-Liouville equation (cf [18,20,31,12]) admits a basis of fundamental solutions u ± ζ with u − ζ = O(exp(−x 2 /2)) and u + ζ = O(x −(1+ζ 2 )/2 exp(x 2 /2)) at infinity. By a change of variable, we deduce that the solution ϕ ζ of problem (18) is a linear combination of an exponentially increasing function denoting by f + ζ and an exponentially decreasing function f − ζ . Moreover f + ζ → +∞ and f − ζ → 0 as t → +∞. Thus, there exist constants a ζ and b ζ which depend continously on ζ such that: We now use this dependence on ζ to determine Θ 0 . Indeed, for ζ = ζ 0 , ϕ ζ0 = Φ is integrable and then b ζ0 = 0. To determine Θ 0 , it is then enough to find the smallest ζ such that the solution ϕ ζ is bounded. Furthermore, we know that the eigenfunction Φ associated with the first eigenvalue Θ 0 and normalized with Φ(0) = 1, holds strictly positive. The positivity of Φ gives a criterion to select functions which constitute a good quasi-modes. Indeed, if for some ζ, the sequence (φ ζ j ) has positive and strictly negative coefficients, then the coefficient b ζ in the decomposition (21) of the associated interpolated functionφ ζ is negative and consequently ζ > ζ 0 . At the opposite, the parameter b ζ is positive for ζ < ζ 0 .

3.4.
Construction of quasi-modes. Discretization (19) gives two behaviors for (φ ζ j ) j (see Figures 1 and 2) and we modify coefficients of (φ ζ j ) j consequently:  • The sequence (φ ζ j ) j remains positive (see Figure 1). We determine n the smallest integer where the sequence (φ ζ j ) j reaches its minimum and we denote L = nh. The restriction ofφ ζ on (0, L) makes a better quasi-mode and we have µ N,N (ζ, L) ≤μ(ζ, L) withμ(ζ, L) the energy of (φ ζ j ) j computed on [0, L]. Nevertheless, as we can not compare µ N,N (ζ, L) and Θ 0 for any L, we modify the sequence by translation so that the minimum equals to 0 and dilation to keep the normalizationφ ζ 1 = 1. We then define the new sequence: The energy associated with a regular interpolation of (ϕ ζ j ) j gives a upperbound of Θ 0 according to Lemma 3.1. The initial sequence (see Figure 1) corresponds to b ζ > 0 in the decomposition (21).
• The sequence (φ ζ j ) j has positive and negative terms (see Figure 2). Let n be the smallest integer such thatφ ζ j0 < 0. We set Lemma 3.1 bounds from above Θ 0 by the energy of the function constructed from (ϕ ζ j ) j . For the initial sequence, b ζ < 0 in the decomposition (21). We would like to find a sequence such that b ζ = 0 and our algorithm could then be seen as a shooting method.
3.6. Estimates of the second eigenvalue. To apply Theorem 2.1, we need an estimate of the second eigenvalue µ 2 (ζ 0 ) of H(ζ 0 ). For this point, we do not need to be very accurate and so we consider the matrix A ζ defined by the discretization of H(ζ) for ζ ∈ [0.76818, 076819]. If we denote by A ζ i,j the coefficients of the matrix A ζ , we have:

3.7.
Accurate estimate for Θ 0 and Φ(0). Lemma 3.3. We have this first coarse bound: The upper-bound was proved in [12] and recalled in Proposition 1. Let us prove the lower-bound. For any ζ ∈ R, we write Choosing ζ = 0 and using Proposition 1, we deduce the lower-bound.
We apply Algorithm 3.2 for h such that 1/h ∈ {100 × k, k = 10, . . . , 40}. For each value, we obtain characteristic values as in Table 1 and we complete this table by computing the lower-bound of Θ 0 given by Theorem 2.1, a lower-bound and a upper-bound for Φ(0) given in Theorem 2.3. To make these computations, we need a lower-bound of |ζ − ζ 0 |. We start with the coarse estimate of Lemma 3.3 and we improve this estimate at each step of the algorithm with the new bounds of Θ 0 . Using the upper-bound µ 2 (ζ 0 ) ≥ 3.315, we obtain 4. Finite element method. In this section, we use a finite element method to analyze the dependence of µ k (ζ) with ζ. We compute the eigenvalues of the operator −d 2 /dt 2 +(t−ζ) 2 on [0, L] with Dirichlet condition on t = L and Neumann condition on t = 0. Let V a disrete variational space. We denote by (μ k (ζ),φ k,ζ ) the k-th discrete eigenpair of the operator in V: The computed eigenvaluesμ k (ζ) give a upper-bound of µ k (ζ). We omit the subscript k when k = 1. Figure 3 illustrates the fact that the minimum of ζ → µ k (ζ) is achieved on the curve ζ → ζ 2 . We observe also the convergence of ζ → µ k (ζ) to 2k − 1 as ζ → +∞. For these computations, we use a finite element method with 10 elements of degree Q 10 on [0, 10]. Let us now use the finite element method to approximate Θ 0 and Φ 0 . With this method, we do not have exact estimate of the error but only a upper-bound for Θ 0 . To determine accurately ζ 0 , we use a finite element method of degree Q 8 or Q 10 and nbel elements. The computational domain is [0, L] and we impose Dirichlet condition on t = L. We compute the first eigenvalueμ(ζ) and compare it with ζ 2 . These computations give also an accurate value for Φ(0) and a 1 . Letφ ζ be the computed normalized positive eigenfunction associated withμ(ζ). Then, we  Tables 2 and 3 give the results of these computations. In particular we obtain approximation for Θ 0 , Φ(0) and a 1 : Θ 0 = 0.590106125,Φ(0) = 0.873043139,ȃ 1 = 0.765188147.
In Table 2, we study the influence of L. As soon L is larger than 7, the cutoff do not interfere the computations which are then satisfactory. In Table 3, we present the computations at a fixed number of degrees of freedom. We observe that the numerical simulations are comparable for degrees larger than 4. Notice that computed valuesμ(ζ) in Tables 2 and 3 provide better upper-bounds for Θ 0 than in Proposition 9.   Table 3. Computation with the finite element method with a constant number of degrees of freedom, L = 10, ζ a = 0.768183653140 and ζ b = 0.768183653141.