\documentclass[twoside]{article}
\usepackage{amsmath, graphicx}
\pagestyle{myheadings}

\markboth{\hfil Computational simulation of vortex phenomena \hfil EJDE/Conf/10}
{EJDE/Conf/10 \hfil John W. Neuberger \& Robert J. Renka \hfil}

\begin{document}
\setcounter{page}{241}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
Fifth Mississippi State Conference on Differential Equations and
Computational Simulations, \newline
Electronic Journal of Differential Equations,
Conference 10, 2003, pp 241--250. \newline
http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp  ejde.math.swt.edu (login: ftp)}
 \vspace{\bigskipamount} \\
%
  Computational simulation of vortex phenomena in superconductors
%
\thanks{ {\em Mathematics Subject Classifications:} 35J50, 64K10, 81V99.
\hfil\break\indent
{\em Key words:} Ginzburg-Landau functional, superconductivity, Sobolev gradient.
\hfil\break\indent
\copyright 2003 Southwest Texas State University. \hfil\break\indent
Published February 28, 2003. } }

\date{}
\author{John W. Neuberger \& Robert J. Renka}
\maketitle

\begin{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.
\end{abstract}

\numberwithin{equation}{section}

\section{Introduction}

    We have developed numerical methods and codes for finding critical
points of the Ginzburg-Landau energy functional \cite{NR}.  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.

\section{Ginzburg-Landau Functional}

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

\begin{equation}\label{e1}
\begin{aligned} E(u,A)  = & E_0 + \frac{1}{2}\int_{\Omega\setminus\Omega_0}
             \Big( \| \nabla u - iuA \|^2 + (\nabla \times A - H_0)^2 +
             \frac{\kappa^2}{2} (|u|^2 - 1)^2 \Big)  \\
   &  + \frac{1}{2}\int_{\Omega_0} (\nabla \times A - H_0)^2 ,
\end{aligned}
\end{equation}
where $\Omega$ is a bounded region in $R^2$, possibly including holes
whose union is denoted by $\Omega_0$, $u \in H^{1,2}(\Omega,C)$ is a
complex wave function (the order parameter) with $|u|^2$ representing
the (probability) density of superconducting electrons (Cooper pairs),
and $A \in H^{1,2}(\Omega,R^2)$ is the vector potential so that
$H = \nabla \times A$ (third component) is the induced magnetic field.
The applied magnetic field (with direction perpendicular to the plane)
is denoted by $H_0$, and $\kappa$ is the
Ginzburg-Landau parameter.  The above formulation of $E$, due to
Jacob Rubinstein \cite{Rub}, can be converted to the more commonly
used formulation of \cite{DGP} by scaling $A$ and $H_0$ by $1/\kappa$
and scaling $E$ by $2/\kappa^2$.

    The critical points of $E$ are solutions to the stationary
Ginzburg-Landau equation (the Euler equation):
\begin{equation}\label{e2}
  \left(\nabla - iA \right)^2 u + \kappa^2 \left(1-|u|^2 \right)u = 0
       \quad \mbox{on } \Omega \setminus \Omega_0 ,
\end{equation}
and
\begin{equation}\label{e3}
  \nabla \times (H-H_0) = \mathop{\rm Im}[u^*(\nabla-iA)u]
       \quad \mbox{on } \Omega ,
\end{equation}
with the corresponding natural boundary conditions
\begin{equation}\label{e4}
  (\nabla - iA)u \cdot n = 0
       \quad \mbox{on } \partial(\Omega \setminus \Omega_0) ,
\end{equation}
\begin{equation}\label{e5}
  H = H_0 \quad \mbox{on } \partial (\Omega) ,
\end{equation}
where $n$ denotes the outward-directed unit normal vector to the
boundary $\partial(\Omega \setminus \Omega_0)$.  Note that, for
$u = \rho e^{i\phi}$, the supercurrent is
\begin{equation}\label{e6}
  J_s  =  \nabla \times (H-H_0)
       =  -\frac{i}{2}(u^* \nabla u-u \nabla u^*) - |u|^2 A
       =  \rho^2 (\nabla \phi - A) ,
\end{equation}
and
  \[ (\nabla - iA)u \cdot n = [\nabla \rho + i \rho(\nabla \phi -A)]
     e^{i \phi} \cdot n = 0  \]
implying that
\begin{equation}\label{e7}
  \nabla \rho \cdot n = 0 \;\; \mbox{and} \;\; J_s \cdot n = 0
       \quad \mbox{on} \quad \partial (\Omega \setminus \Omega_0) .
\end{equation}
Also, $u=0, J_s=0$, and $H-H_0$ is constant in $\Omega_0$.

\section{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 \cite{DGP2, Mu}.  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 \ref{e1} 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 \cite{N, NR}.  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
$\Omega$.  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 $\Omega = [0,\Delta x] \times [0,\Delta y]$ into a
$k_x$-cell by $k_y$-cell rectangular grid with mesh widths
$h_x = \Delta x/k_x$ and $h_y = \Delta y/k_y$.  Any of the cells
may be designated as part of $\Omega_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 $\kappa$ 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(\Omega)$.

\section{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 $\kappa < \sqrt{2}/2$,
or repel each other if $\kappa > \sqrt{2}/2$ \cite{J}.  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.

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig1.eps}
  \end{center}
  \caption{Magnetic field at iteration 200}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig2.eps}
  \end{center}
  \caption{Magnetic field at iteration 1000}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig3.eps}
  \end{center}
   \caption{Magnetic field at iteration 1700}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig4.eps}
  \end{center}
  \caption{Magnetic field at iteration 2000}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig5.eps}
  \end{center}
  \caption{Energy versus iteration count}
\end{figure}

    A typical test result is depicted in Figures 4.1--4.5.  For
$\kappa = \sqrt{2}/2$ and a 30 by 30 rectangle with a 64-cell by
64-cell mesh, we obtained a sequence of near-critical 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$.

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig6.eps}
  \end{center}
  \caption{Vortex and antivortex on a domain with holes}
\end{figure}

    We next constructed an initial estimate with vortices at the
same locations and reduced $\kappa$ 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 $\kappa$ to $1.0$,
and constructed an 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.

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig7.eps}
  \end{center}
    \caption{Two vortices on a domain with holes, $\kappa = .5$}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig8.eps}
  \end{center}
  \caption{Degree-2 vortex on a domain with no holes, $\kappa = .5$}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig9.eps}
  \end{center}
  \caption{Magnetic field on a domain with holes, $\kappa = 1.0$}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
    \includegraphics[width=0.7\textwidth]{fig10.eps}
    \end{center}
    \caption{Magnetic field on a domain with no holes, $\kappa = 1.0$}
\end{figure}

    For our final tests we duplicated the configuration used in
\cite{DGP2}:  $\Delta x = \Delta y = 1, \kappa = 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.

\begin{figure}[tbh]
  \begin{center}
  \includegraphics[width=0.7\textwidth]{fig11.eps}
    \end{center}
    \caption{Electron density on a domain with no holes}
\end{figure}

\begin{figure}[tbh]
  \begin{center}
   \includegraphics[width=0.7\textwidth]{fig12.eps}
   \end{center}
   \caption{Electron density on a domain with 10 holes}
\end{figure}


\begin{thebibliography}{99} \frenchspacing

\bibitem{DGP}  Qiang Du, M. Gunzburger, and Janet S. Peterson,
             {\em Analysis and Approximation of the Ginzburg- Landau Model
             of Superconductivity}, SIAM Review 34 (1992), 54--81.

\bibitem{DGP2}  Qiang Du, M. Gunzburger, and Janet S. Peterson, {\em
             Computational simulation of type-II superconductivity
             including pinning phenomena}, Physical Review B 51,
             22 (June 1995), 16194--16203.

\bibitem{J}  A. Jaffe and C. Taubes, {\em Vortices and Monopoles},
             Birk\"{a}user, 1980.

\bibitem{Mu}  Mo Mu, {\em A linearized Crank-Nicolson-Galerkin
             method for the Ginzburg-Landau model}, SIAM J. Sci. Comput.
             18 (1997), 1028--1039.

\bibitem{N}  J. W. Neuberger, {\em Sobolev gradients and
             differential equations}, Springer Lecture Notes in Mathematics
             \#1670, 1997.

\bibitem{Ne}  J. W. Neuberger and R. J. Renka, {\em Numerical
             Calculation of Singularities for Ginzburg-Landau Functionals},
             EJDE 1997, 10 (1997), 1--4.

\bibitem{NR}  J. W. Neuberger and R. J. Renka, {\em Sobolev
             Gradients and the Ginzburg-Landau Functional}, SIAM J. Sci.
             Comput. 20 (1998), 582--590.

\bibitem{Rub}  J. Rubinstein,
             {\em Six lectures on superconductivity.  Interfaces and
             Transitions.}, CRM Proc. Lec. Notes 13, Amer. Math. Soc., 1998.

\end{thebibliography}

\noindent\textsc{John W. Neuberger}\\
        Department of Mathematics,
        University of North Texas,\\
        Denton, TX 76203-1430, USA\\
        e-mail: jwn@unt.edu \smallskip

\noindent\textsc{Robert J. Renka}\\
        Department of Computer Sciences,
        University of North Texas,\\
        Denton, TX 76203-1366, USA\\
        e-mail: renka@cs.unt.edu

\end{document}