Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 10, 2003, pp 241–250. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Computational simulation of vortex phenomena in superconductors ^{∗}

John W. Neuberger & Robert J. Renka

Abstract We have developed a numerical method for approximating critical points of the Ginzburg-Landau functional. We briefly describe the capabilities of our codes and present some test results for the two-dimensional case in the form of plots of electron densities and magnetic fields. Our results demonstrate that vortices can be pinned to small holes corresponding to normal (non-superconducting) regions.

1 Introduction

We have developed numerical methods and codes for finding critical points of the Ginzburg-Landau energy functional [7]. The efficiency and versatility of our code enables us to very rapidly simulate the effects of changes in domain geometry, external magnetic field, and parameter values. We report the results of several tests here. Among our findings is the fact that, with no external magnetic field, vortices and antivortices can be pinned to arbitrary locations by placing small holes in the domain. These results have implications for the design of superconducting devices.

We describe the problem in Section 2, discuss our methods in Section 3, and present test results in Section 4.

2 Ginzburg-Landau Functional

The two-dimensional Ginzburg-Landau energy functional in nondimensionalized units is

\begin{matrix}E(u, A) = & E_{0} + \frac{1}{2} ∫_{Ω\Ω_{0}} (∥∇u − iuA∥^{2} + (∇ × A − H_{0})^{2} + \frac{κ^{2}}{2} (|u|^{2} − 1)^{2}) \\ + \frac{1}{2} ∫_{Ω_{0}}(∇ × A − H_{0})^{2}, \end{matrix}

(2.1) where Ω is a bounded region in R^{2}, possibly including holes whose union is denoted by Ω_{0}, u ∈ H^{1,2}(Ω, C) is a complex wave function (the order parameter) with |u|^{2} representing the (probability) density of superconducting electrons (Cooper pairs), and A ∈ H^{1,2}(Ω, R^{2}) is the vector potential so that H = ∇ × A (third component) is the induced magnetic field. The applied magnetic field (with direction perpendicular to the plane) is denoted by H_{0}, and κ is the Ginzburg-Landau parameter. The above formulation of E, due to Jacob Rubinstein [8], can be converted to the more commonly used formulation of [1] by scaling A and H_{0} by 1/κ and scaling E by 2/κ^{2}.

∗Mathematics Subject Classifications: 35J50, 64K10, 81V99.

Key words: Ginzburg-Landau functional, superconductivity, Sobolev gradient.

\begin{matrix}c \\ ◯ \end{matrix}2003 Southwest Texas State University.

Published February 28, 2003.

The critical points of E are solutions to the stationary Ginzburg-Landau equation (the Euler equation):

(∇ − iA)^{2} u + κ^{2} (1 − |u|^{2}) u = 0 on Ω \ Ω_{0},  (2.2)

and

∇ × (H − H_{0}) = Im[u^{∗}(∇ − iA)u] on Ω,  (2.3)

with the corresponding natural boundary conditions

(∇ − iA)u · n = 0 on ∂(Ω \ Ω_{0}),  (2.4)
H = H_{0} on ∂(Ω),  (2.5)

where n denotes the outward-directed unit normal vector to the boundary ∂(Ω\ Ω_{0}). Note that, for u = ρe^{iϕ}, the supercurrent is

J_{s} = ∇ × (H − H_{0}) = −\frac{i}{2}(u^{∗}∇u − u∇u^{∗}) − |u|^{2}A = ρ^{2}(∇ϕ − A), (2.6)

and

(∇ − iA)u · n = [∇ρ + iρ(∇ϕ − A)]e^{iϕ} · n = 0

implying that

∇ρ · n = 0 and J_{s} · n = 0 on ∂(Ω \ Ω_{0}).  (2.7)

Also, u = 0, J_{s} = 0, and H − H_{0} is constant in Ω_{0}.

3 Method and Code

The conventional approach to numerical solution of the Ginzburg-Landau equations is to evolve the time-dependent equations until a steady state is achieved [2, 4]. While no benchmark results are available, we believe our approach has an advantage in terms of simplicity and computational efficiency. Rather than dealing with the coupled nonlinear PDE’s and nonlinear boundary conditions, we treat the minimization problem directly; i.e., we minimize a finite difference discretization of 2.1 by a steepest descent method. However, in place of the discretized L_{2}-gradient, which leads to the notoriously inefficient standard method of steepest descent, we use a discretization of a Sobolev gradient [5, 7].

This may be viewed as a very effective preconditioning method in which the preconditioner arises in a natural manner, as contrasted with the more common experimental approach. The method applies to a much larger class of functionals, and while we have implemented it for a three-dimensional domain, we treat only the two-dimensional problem here.

By treating the least squares minimization problem directly, we not only avoid dealing with the nonlinear boundary conditions, but we also need not remove the freedom induced by gauge invariance. Furthermore, we are able to obtain different critical points for each tested configuration by simply varying the initial estimate. We have two methods for generating initial estimates: one random and one which allows an interactive user to place vortices and/or antivortices (of arbitrary degrees) at any locations in the domain Ω. The latter method produces good approximations to critical points in the case of no external magnetic field — H_{0} = 0.

Our finite difference discretization is based on partitioning the rectangle Ω = [0, Δx] × [0, Δy] into a k_{x}-cell by k_{y}-cell rectangular grid with mesh widths h_{x} = Δx/k_{x} and h_{y} = Δy/k_{y}. Any of the cells may be designated as part of Ω_{0} which is thus constrained to have boundary edges parallel to the axes (although any shape can be approximated). The unknowns are the grid-point values of the real and imaginary parts of u and the two components of A, and the approximation to E consists of a weighted sum of cell-center values obtained by averaging or differentiating grid-point values. The Ginzburg-Landau parameter κ need not be constant, thus allowing for different materials (as well as holes) in different portions of the domain. Similarly, the external field H_{0} may be specified as an arbitrary function in L_{2}(Ω).

4 Test Results

For a rectangular domain with no external magnetic field, a critical point of the energy functional cannot exhibit both vortices (spin up) and antivortices (spin down). More specifically, in the solution of the time-dependent Ginzburg-Landau equation, an initial estimate with a vortex/antivortex pair evolves in such a way that the vortex and antivortex are attracted to each other and eventually annihilate each other. A pair of vortices, on the other hand either attract each other, eventually merging into a degree-2 vortex if κ < \sqrt{2}/2, or repel each other if κ > \sqrt{2}/2 [3]. To some extent our descent procedure mimics the evolution of the time-dependent Ginzburg-Landau equations, and we sought to numerically test this theoretical characterization of vortex interaction.

We used rectangular domains of various sizes and interpolated random initial estimates on a coarse (two-cell by two-cell) mesh to a fine mesh (up to 128 cells in each direction), resulting in smooth initial estimates and fast convergence. We found that the interaction between widely separated vortices/antivortices is extremely weak. In particular, a (u, A) pair with both vortices and antivortices can be arbitrarily close to (and thus numerically indistinguishable from) a saddle point of the Ginzburg-Landau functional. Such a near-critical point (‘metastable’ state) must have physical significance, and the numerics thus reveal some interesting phenomena that is not apparent from the analysis.

Figure 1: Magnetic field at iteration 200

A typical test result is depicted in Figures 4.1–4.5. For κ = \sqrt{2}/2 and a 30 by 30 rectangle with a 64-cell by 64-cell mesh, we obtained a sequence of nearcritical points (saddle points) with three vortices and one antivortex (Figures 4.1–4.3), and a critical point (local minimum) with two vortices (Figure 4.4). Figures 4.1–4.4 display the magnetic field associated with iteration counts 200, 1000, 1700, and 2000, respectively. A plot of energy as a function of iteration count is depicted in Figure 4.5. It is clear from this plot that the decrease in energy is extremely slow as a vortex approaches the antivortex but then drops rapidly once there is overlap.

We now consider the effect of adding small holes (two cells by two cells in all cases) to the domain. We centered a hole at (23h, 31h) and another at (39h, 31h) for h = h_{x} = h_{y} = 30/64, and we constructed an initial estimate with a vortex centered at (23h, 31h) and an antivortex centered at (39h, 31h). The method rapidly converged to a critical point with a vortex and antivortex centered on the holes; i.e., the initial vortex locations did not change. The magnetic field is depicted in Figure 4.6. For the same initial estimate with no holes in the domain, the vortex and antivortex quickly annihilated each other resulting in the zero-energy solution associated with constant u, |u| = 1, and A = 0.

We next constructed an initial estimate with vortices at the same locations and reduced κ to 0.5. Again, the vortices remained centered on the holes. See Figure 4.7. Starting with the same initial estimate without holes, the vortices merged into a single degree-2 vortex — Figure 4.8. Finally, we placed the hole centers at (29h, 31h) and (33h, 31h), increased κ to 1.0, and constructed an

Figure 2: Magnetic field at iteration 1000

initial estimate with vortices centered on the holes. The critical point again retained the initial vortex locations (Figure 4.9), but when run without the holes, the vortices separated as predicted by theory (Figure 4.10). All of the above tests ran very quickly, requiring at most 1000 descent steps and only a few minutes execution time on a 550-Mhz Pentium III.

For our final tests we duplicated the configuration used in [2]: Δx = Δy = 1, κ = 20, and H_{0} = 200. (Note that our spatial variables have units equal to penetration depth while those used by Du, et al have units of coherence length). We took k_{x} = k_{y} = 64. With no holes in the domain and starting with u = 1, A = 0, we obtained the highly symmetric 20-vortex solution whose electron density |u|^{2} is depicted in Figure 4.11. Note, however, that this critical point is not unique. With different initial estimates we obtained different patterns of vortices. We then added 10 holes with roughly the same distribution as those used by Du, et al, and obtained the pattern shown in Figure 4.12. With one exception, every hole is within a vortex. The failure of the lower left hole to pin a vortex may be because it is too close to the boundary. Due to the large number of vortices these tests were much more expensive than those with H_{0} = 0 — about 20000 descent steps requiring 1.5 hours on the Pentium III.

References

[1] Qiang Du, M. Gunzburger, and Janet S. Peterson, Analysis and Approximation of the Ginzburg- Landau Model of Superconductivity, SIAM Review 34 (1992), 54–81.

Figure 3: Magnetic field at iteration 1700

[2] Qiang Du, M. Gunzburger, and Janet S. Peterson, Computational simulation of type-II superconductivity including pinning phenomena, Physical Review B 51, 22 (June 1995), 16194–16203.

[3] A. Jaffe and C. Taubes, Vortices and Monopoles, Birkäuser, 1980.

[4] Mo Mu, A linearized Crank-Nicolson-Galerkin method for the Ginzburg-Landau model, SIAM J. Sci. Comput. 18 (1997), 1028–1039.

[5] J. W. Neuberger, Sobolev gradients and differential equations, Springer Lecture Notes in Mathematics #1670, 1997.

[6] J. W. Neuberger and R. J. Renka, Numerical Calculation of Singularities for Ginzburg-Landau Functionals, EJDE 1997, 10 (1997), 1–4.

[7] J. W. Neuberger and R. J. Renka, Sobolev Gradients and the Ginzburg-Landau Functional, SIAM J. Sci. Comput. 20 (1998), 582–590.

[8] J. Rubinstein, Six lectures on superconductivity. Interfaces and Transitions., CRM Proc. Lec. Notes 13, Amer. Math. Soc., 1998.

John W. Neuberger Department of Mathematics, University of North Texas, Denton, TX 76203-1430, USA e-mail: jwn@unt.edu Robert J. Renka Department of Computer Sciences, University of North Texas, Denton, TX 76203-1366, USA e-mail: renka@cs.unt.edu

EJDE/Conf/10  John W. Neuberger & Robert J. Renka  247

Figure 4: Magnetic field at iteration 2000 Figure 5: Energy versus iteration count

248  Computational simulation of vortex phenomena  EJDE/Conf/10

Figure 6: Vortex and antivortex on a domain with holes Figure 7: Two vortices on a domain with holes, κ = .5 Figure 8: Degree-2 vortex on a domain with no holes, κ = .5

EJDE/Conf/10  John W. Neuberger & Robert J. Renka  249

Figure 9: Magnetic field on a domain with holes, κ = 1.0 Figure 10: Magnetic field on a domain with no holes, κ = 1.0

250  Computational simulation of vortex phenomena  EJDE/Conf/10

Figure 11: Electron density on a domain with no holes Figure 12: Electron density on a domain with 10 holes