Mathematical and numerical analysis of time-dependent Ginzburg–Landau equations in nonconvex polygons based on Hodge decomposition This work is supported in part by the National Natural Science Foundation of China (NSFC) under grants No. 11301262, No. 11471031, No. 91430216, and the US National Science Foundation (NSF) through grants DMS-1115530 and DMS-1419040.

# Mathematical and numerical analysis of time-dependent Ginzburg–Landau equations in nonconvex polygons based on Hodge decomposition ††thanks: This work is supported in part by the National Natural Science Foundation of China (NSFC) under grants No. 11301262, No. 11471031, No. 91430216, and the US National Science Foundation (NSF) through grants DMS-1115530 and DMS-1419040.

Buyang Li Department of Mathematics, Nanjing University, Nanjing, 210093, Jiangsu, P.R. China. (buyangli@nju.edu.cn)    Zhimin Zhang Beijing Computational Science Research Center, Beijing, 100084, P.R. China.
Department of Mathematics, Wayne State University, Detroit, MI 48202, USA. (zzhang@math.wayne.edu)
###### Abstract

We prove well-posedness of time-dependent Ginzburg–Landau system in a nonconvex polygonal domain, and decompose the solution as a regular part plus a singular part. We see that the magnetic potential is not in in general, and the finite element method (FEM) may give incorrect solutions. To remedy this situation, we reformulate the equations into an equivalent system of elliptic and parabolic equations based on the Hodge decomposition, which avoids direct calculation of the magnetic potential. The essential unknowns of the reformulated system admit solutions and can be solved correctly by the FEMs. We then propose a decoupled and linearized FEM to solve the reformulated equations and present error estimates based on proved regularity of the solution. Numerical examples are provided to support our theoretical analysis and show the efficiency of the method.

Key words. superconductivity, reentrant corner, singularity, well-posedness, finite element method, convergence, Hodge decomposition

AMS subject classifications. 35Q56, 35K61, 65M12, 65M60

## 1 Introduction

The Ginzburg–Landau theory, initially introduced by Ginzburg and Landau [16] and subsequently extended to the time-dependent case by Gor’kov and Eliashberg [18], are widely used to describe the phenomena of superconductivity in both low and high temperatures [11, 22]. In a two-dimensional domain , the time-dependent Ginzburg–Landau model (TDGL) is governed by two equations (with the Lorentz gauge),

 η∂ψ∂t+(iκ∇+A)2ψ+(|ψ|2−1)ψ−iηκψ∇⋅A=0, \hb@xt@.01(1.1) ∂A∂t+∇×(∇×A)−∇(∇⋅A)+Re[ψ∗(iκ∇+A)ψ]=∇×f, \hb@xt@.01(1.2)

where and are given positive constants, the order parameter is an unknown complex scalar function and denotes the complex conjugate of , the real-vector valued function denotes the unknown magnetic potential, and the scalar function denotes the external magnetic field, and we have used the notations

 ∇×A=∂A2∂x1−∂A1∂x2,∇⋅A=∂A1∂x1+∂A2∂x2, ∇×f=(∂f∂x2,−∂f∂x1),∇ψ=(∂ψ∂x1,∂ψ∂x2).

The natural boundary and initial conditions for this problem are

 ∇ψ⋅n=0,A⋅n=0,∇×A=f,on∂Ω×(0,T), \hb@xt@.01(1.3) ψ(x,0)=ψ0(x),A(x,0)=A0(x),inΩ, \hb@xt@.01(1.4)

where denotes the unit outward normal vector on the boundary .

The TDGL has been widely studied both theoretically and numerically. Existence and uniqueness of the solution for (LABEL:PDE1)-(LABEL:PDE2) in a smooth domain were proved by Chen et al. [8], where equivalence of (LABEL:PDE1)-(LABEL:PDE2) to the Ginzburg–Landau equations under the temporal gauge was proved. Various numerical methods for solving the TDGL were reviewed in [12, 14]. In contrast with the many numerical approximation schemes, numerical analysis of the model seems very limited so far. Error analysis of a Galerkin finite element method (FEM) with an implicit backward Euler time-stepping scheme was presented in [7, 13], where optimal-order convergence rate of the numerical solution was proved for sufficiently regular solution. A linearized Crank–Nicolson scheme was proposed in [24] for a regularized TDGL under the temporal gauge without error analysis. An alternating Crank–Nicolson scheme was proposed in [25] and error estimates were presented for a regularized TDGL under the grid-ratio restriction , where and are the time-step size and spatial mesh size. Although convergence of the numerical solutions has been proved in [7, 13, 25] in smooth domains, these error estimates may not hold in a domain with corners, where the regularity of the solution may not satisfy the conditions required in the analysis. It has been reported in [15, 24] that the numerical solution of the magnetic potential by the FEM often exhibits undesired singularities around a corner. To resolve this problem, a mixed FEM was proposed in [6] to approximate the triple , in a finite element subspace of , which requires less regularity of intuitively, and error estimates of the finite element solution were presented under the assumption that is in . Recently, an optimal-order error estimate of the FEM with a linearized Crank–Nicolson scheme was presented in [15] without restriction on the grid ratio, but the analysis requires stronger regularity of the solution and the domain. On one hand, existing theoretical and numerical analysis of the model all require the magnetic potential to be in . In a domain with reentrant corners, however, the magnetic potential may not be in and well-posedness of the TDGL remains open. On the other hand, numerical approximations of the TDGL in domains with reentrant corners are important for physicists to study the effects of surface defects in superconductivity [2, 26], which are often accomplished by solving (LABEL:PDE1)-(LABEL:PDE2) directly with the finite element or finite difference methods, without being aware of the danger of these numerical methods.

In this paper, we study the TDGL in a nonconvex polygon, possibly with reentrant corners. We shall prove that the system (LABEL:PDE1)-(LABEL:init) is well-posed, with for some which depends on the interior angles of the reentrant corners. As shown in the numerical examples, with such low-regularity, the FEM may give an incorrect solution for the magnetic potential , which further pollutes the numerical solution of due to the coupling of equations. We are interested in reformulating (LABEL:PDE1)-(LABEL:init) into an equivalent form which can be solved correctly by the FEMs, as they are preferred when using software packages and when other equations are coupled with the Ginzburg–Landau equations. Our idea is to apply the Hodge decomposition , and consider the projection of (LABEL:PDE2) onto the divergence-free and curl-free subspaces, respectively. Then (LABEL:PDE1)-(LABEL:init) is reformulated as

 η∂ψ∂t+(iκ∇+A)2ψ+(|ψ|2−1)ψ−iηκψ∇⋅A=0, \hb@xt@.01(1.5) \hb@xt@.01(1.6) Δq=∇⋅(Re[ψ∗(iκ∇+A)ψ]) \hb@xt@.01(1.7) ∂u∂t−Δu=f−p, \hb@xt@.01(1.8) ∂v∂t−Δv=−q, \hb@xt@.01(1.9)

with the boundary and initial conditions

 ∇ψ⋅n=0,p=0,∇q⋅n=0,u=0,∇v⋅n=0,on ∂Ω×(0,T], \hb@xt@.01(1.10) ψ(x,0)=ψ0(x),u(x,0)=u0(x),v(x,0)=v0(x),in Ω, \hb@xt@.01(1.11)

where and are just the divergence-free and curl-free parts of , respectively, and are defined by

 {−Δu0=∇×A0in  Ω,u0=0on  ∂Ω,and{Δv0=∇⋅A0in  Ω,∂nv0=0on  ∂Ω,

with . We shall prove that the solution of the projected TDGL (LABEL:RFPDE1)-(LABEL:RFinit) coincides with the solution of (LABEL:PDE1)-(LABEL:init). Then we propose a decoupled and linearized FEM to solve (LABEL:RFPDE1)-(LABEL:RFinit), and establish error estimates based on proved regularity of the solution. Our main results are presented in Section 2, and we prove these results in Section 3-5. In Section 6, we present numerical examples to support our theoretical analysis. Due to limitations on pages, derivations of the system (LABEL:RFPDE1)-(LABEL:RFinit) are presented in a separate paper [21], where the efficiency of the method is shown via numerical simulations in comparison with the traditional approaches of solving the TDGL directly under the temporal gauge and the Lorentz gauge.

## 2 Main results

For any nonnegative integer , we let , and denote the the conventional Sobolev spaces of real-valued and complex-valued functions defined in , respectively, with , , and ; see [1]. For a positive real number , with , we define via the complex interpolation; see [3]. We denote , , , , and let denote the subspace of consisting of functions whose traces are zero on . For any two functions we define

 (f,g)=∫Ωf(x)g(x)∗dx,

where denotes the complex conjugate of , and define

 Lp=Lp×Lp,Hs=Hs×Hs,H1n(Ω):={a∈H1×H1:a⋅n=0 on ∂Ω}, Hn(curl,div)={a∈L2:∇×a∈L2, ∇⋅a∈L2 and a⋅n=0 on ∂Ω}, H(curl)={g∈L2:∇×g∈L2}.
###### Definition 2.1

(Weak solutions of (LABEL:PDE1)-(LABEL:init))   Let denote the maximal interior angle of the nonconvex polygon . The pair is called a weak solution of (LABEL:PDE1)-(LABEL:init) if

 ∂tψ,Δψ∈L2((0,T);L2),|ψ|≤1  a.e.~{}in~{}\,Ω×(0,T), A∈C([0,T];L2)∩L∞((0,T);Hn(curl,div)), ∂tA∈L2((0,T);L2),∇×A,∇⋅A∈L2((0,T);H1),

for any , with , , and the variational equations

 ∫T0[(η∂ψ∂t,φ)+((iκ∇+A)ψ,(iκ∇+A)φ)]dt +∫T0[((|ψ|2−1)ψ−iηκψ∇⋅A,φ)]dt=0, \hb@xt@.01(2.1) ∫T0[(∂A∂t,a)+(∇×A,∇×a)+(∇⋅A,∇⋅a)]dt =∫T0[(f,∇×a)−(Re[ψ∗(iκ∇+A)ψ],a)]dt, \hb@xt@.01(2.2)

hold for all and .

###### Definition 2.2

(Weak solutions of (LABEL:RFPDE1)-(LABEL:RFinit))   Let denote the maximal interior angle of the nonconvex polygon . The quintuple is called a weak solution of (LABEL:RFPDE1)-(LABEL:RFinit) if

 ∂tψ,Δψ∈L2((0,T);L2),|ψ|≤1  a.e.~{}in~{}\,Ω×(0,T), p∈L∞((0,T);˚H1),q∈L∞((0,T);H1),u∈C([0,T];˚H1),v∈C([0,T];H1), ∂tu,∂tv,Δu,Δv∈L∞((0,T);L2)∩L2((0,T);H1),

for any , with , , , and the variational equations

 ∫T0[(η∂ψ∂t,φ)+((iκ∇+A)ψ,(iκ∇+A)φ)]dt +∫T0[((|ψ|2−1)ψ−iηκψ∇⋅A,φ)]dt=0, \hb@xt@.01(2.3) ∫T0(∇p,∇ξ)dt=∫T0(Re[ψ∗(iκ∇+A)ψ],∇×ξ)dt \hb@xt@.01(2.4) ∫T0(∇q,∇ζ)dt=∫T0(Re[ψ∗(iκ∇+A)ψ],∇ζ)dt \hb@xt@.01(2.5) ∫T0[(∂u∂t,θ)+(∇u,∇θ)]dt=∫T0(f−p,θ)dt, \hb@xt@.01(2.6) ∫T0[(∂v∂t,ϑ)+(∇v,∇ϑ)]dt=−∫T0(q,ϑ)dt, \hb@xt@.01(2.7)

hold for all , and .

Our first result is the well-posedness and equivalence of the systems (LABEL:PDE1)-(LABEL:init) and (LABEL:RFPDE1)-(LABEL:RFinit), which are presented in the following theorem.

###### Theorem 2.1

(Well-posedness and equivalence of the two systems)
If , , and a.e. in , then the system (LABEL:PDE1)-(LABEL:init) admits a unique weak solution in the sense of Definition LABEL:DefWSol, and the system (LABEL:RFPDE1)-(LABEL:RFinit) admits a unique solution which coincides with the solution of (LABEL:PDE1)-(LABEL:init).

Moreover, if we let , , be the reentrant corners of the domain , then the solution has the decomposition

 ψ(x,t)=Ψ(x,t)+m∑j=1αj(t)Φ(|x−xj|)|x−xj|π/ωjcos(πΘj(x)/ωj), A=∇×u+∇v

with

 u(x,t)=˜u(x,t)+m∑j=1βj(t)Φ(|x−xj|)|x−xj|π/ωjsin(πΘj(x)/ωj), v(x,t)=˜v(x,t)+m∑j=1γj(t)Φ(|x−xj|)|x−xj|π/ωjcos(πΘj(x)/ωj),

where , , is a given smooth cut-off function which equals in a neighborhood of , is the angle shown in Figure LABEL:LshapeD, and .

Further regularity of the solution is presented below, which is needed in the analysis of the convergence of the numerical solution.

###### Theorem 2.2

(Further regularity)
If , , , , , , , a.e. in , and the compatibility conditions

 ∂nψ0=0and∇×A0=f(⋅,0)on∂Ω

are satisfied, then the solution of (LABEL:RFPDE1)-(LABEL:RFinit) possesses the regularity

 ψ∈C([0,T];H1+s),∂tψ∈L2((0,T);H1+s),∂ttψ∈L2((0,T);L2), p,q∈L∞((0,T);H1),u,v∈C([0,T];H1+s), ∂tu,∂tv∈L2((0,T);H1+s),∂ttu,∂ttv∈L2((0,T);L2)

for any .

To solve the reformulated system (LABEL:RFPDE1)-(LABEL:RFinit), we propose a decoupled and linearized Galerkin FEM. For this purpose, we let be a quasi-uniform triangulation of the domain and denote the mesh size by . Let denote the space of complex-valued piecewise linear functions subject to the triangulation, let denote the space of real-valued piecewise linear functions, and set . Clearly, , and are finite dimensional subspaces of , and , respectively. Let be the commonly used Lagrange interpolation operator onto the finite element spaces. For any positive integer , we let be a uniform partition of the time interval and set . For any sequence of functions , we define , and we define a cut-off function by

 χ(z)=z/max(|z|,1),∀ z∈C,

which is Lipschitz continuous and satisfies that , .

We look for , and satisfying the equations

 (Dτψn+1h,φ)+((iκ−1∇+Anh)ψn+1h,(iκ−1∇+Anh)φ) \hb@xt@.01(2.8) (∇pn+1h,∇ξ)=(Re[χ(ψnh)∗(iκ−1∇ψn+1h+Anhψn+1h)],∇×ξ) \hb@xt@.01(2.9) (∇qn+1h,∇ζ)=(Re[χ(ψnh)∗(iκ−1∇ψn+1h+Anhψn+1h)],∇ζ) \hb@xt@.01(2.10) (Dτun+1h,θ)+(∇un+1h,∇θ)=(fn+1−pn+1h,θ) \hb@xt@.01(2.11) (Dτvn+1h,ϑ)+(∇vn+1h,∇ϑ)=(−qn+1h,ϑ), \hb@xt@.01(2.12)

for all , and , with , where and are solved from

 (∇u0h,∇ξ)=(A,∇×ξ),∀ ξ∈˚V1h, \hb@xt@.01(2.13) (∇v0h,∇ζ)=(A,∇⋅ζ), ∀ ζ∈V1h, \hb@xt@.01(2.14)

and is the Lagrange interpolation of .

For the proposed scheme, we have the following theorem concerning the convergence of the numerical solution.

###### Theorem 2.3

(Convergence of the finite element solution)
The finite element system
(LABEL:FEMEq1)-(LABEL:FEMEq5) admits a unique solution when and, under the assumptions of Theorem LABEL:MainTHM2,

 max1≤n≤N(∥un−unh∥H1+∥vn−vnh∥H1+∥An−Anh∥L2+∥ψn−ψnh∥L2)≤C(τ+hs),

where is a positive constant independent of and .

In the rest part of this paper, we prove Theorem LABEL:MainTHM1LABEL:MainTHM3. To simplify the notations, we denote by a generic positive constant which may be different at each occurrence but is independent of , and .

## 3 Proof of Theorem LABEL:MainTHM1

In this section, we prove well-posedness of the Ginzburg–Landau equations in a nonconvex polygon and equivalence of the two formulations (LABEL:PDE1)-(LABEL:init) and (LABEL:RFPDE1)-(LABEL:RFinit). Compared with smooth domains, in a nonconvex polygon, the space is not equivalent to and is not embedded into for large . Convergence of the nonlinear terms of the approximating solutions needs to be proved based on the weaker embedding in the compactness argument, and uniqueness of solution needs to be proved based on weaker regularity of the solution.

### 3.1 Preliminaries

Firstly, we cite a lemma concerning the regularity of Poisson’s equations in a nonconvex polygon [10, 17].

###### Lemma 3.1

The solution of the Poisson equations

 {Δw=gin  Ω,w=0on  ∂Ω,and{Δw=gin  Ω,∂nw=0on  ∂Ω,

satisfies that the Neumann problem requires

 ∥w∥W1,ps+∥w∥H1+s≤Cs∥g∥L2,∀ s∈(1/2,π/ω),

where when .

Secondly, we introduce a lemma concerning the embedding of into .

###### Lemma 3.2

for any .

Proof. From [6] we know that has the decomposition , where and are the solutions of

 {−Δu=∇×Ain  Ω,u=0on  ∂Ω,and{Δv=∇⋅Ain  Ω,∂nv=0on  ∂Ω,

respectively, with . For the two Poisson’s equations, Lemma LABEL:HsPoissEq implies that

 ∥u∥H1+s+∥v∥H1+s≤Cs(∥∇×A∥L2+∥∇⋅A∥L2),∀ s∈(1/2,π/ω).

Thirdly, we introduce a lemma concerning the embedding of discrete Sobolev spaces.

###### Lemma 3.3

Let , with . If we define and by

 (Δhθh,φ)=−(∇θh,∇φ),∀φ∈˚V1h, (Δhϑh,φ)=−(∇ϑh,∇φ),∀φ∈V1h,

then

 ∥∇θh∥L4≤C∥Δhθh∥L2,and∥∇ϑh∥L4≤C∥Δhϑh∥L2.

Proof. Let be the solution of the Poisson’s equation

 Δθ=Δhθh

with the Dirichlet boundary condition on . Then for any , which implies that, via the standard -norm error estimate and Lemma LABEL:HsPoissEq,

 ∥∇(θ−θh)∥L2≤C∥θ∥H1+shs≤C∥Δhθh∥L2hs.

Since , by applying the inverse inequality we obtain that

 ∥∇(Ihθ−θh)∥L4≤Ch−1/2∥∇(Ihθ−θh)∥L2≤C∥Δhθh∥L2hs−1/2≤C∥Δhθh∥L2.

Thus . The proof for is similar.

### 3.2 Existence of weak solutions for (Label:pde1)-(LABEL:init)

In this subsection, we prove existence of weak solutions for the system (LABEL:PDE1)-(LABEL:init) by constructing approximating solutions in finite dimensional spaces and then applying a compactness argument. Firstly, we need the following lemma to control the order parameter pointwisely.

###### Lemma 3.4

For any given , the equation (LABEL:PDE1) has at most one weak solution in the sense of (LABEL:VPDE1). If the solution exists then it satisfies that a.e. in .

Proof. From Lemma LABEL:LemCDReg we see that and so . Uniqueness of the solution can be proved easily based on the regularity assumption of . To prove a.e. in , we integrate (LABEL:PDE1) against and consider the real part, where denotes the positive part of . For any we have

 ∫Ω(η4(|ψ(x,t′)|2−1)2+)dx+∫t′0∫Ω(|ψ|2−1)2+|ψ|2dxdt =−∫t′0Re∫Ω(iκ∇ψ+Aψ)(−iκ∇+A)[ψ∗(|ψ|2−1)+]dxdt =−∫t′0∫Ω∣∣∣iκ∇ψ+Aψ∣∣∣2(|ψ|2−1)+dxdt +∫t′0Re∫{|ψ|2>1}(iκ∇ψ+Aψ)ψ∗(iκψ∇ψ∗+iκψ∗∇ψ)dxdt =−∫t′0∫Ω∣∣∣iκ∇ψ+Aψ∣∣∣2(|ψ|2−1)+dxdt −∫t′0Re∫{|ψ|2>1}(|ψ|2|∇ψ|2+(ψ∗)2∇ψ⋅∇ψ)dxdt ≤0,

which implies that . Thus a.e. in .

Secondly, we construct approximating solutions in finite dimensional spaces. For this purpose, we let be the eigenfunctions of the Neumann Laplacian, which form a basis of . Let be defined by

 (Mu,v)=(∇×u,∇×