Introduction

The numerical simulation of superconductivity has attracted increasing attention during the last years, specially due to the appearance of high-temperature superconductors. At the mesoscopic level, the governing equations are provided by the Ginzburg-Landau theory, and are frequently referred to as TDGL equations (for Time Dependent Ginzburg-Landau equations). These are coupled nonlinear partial differential equations for the (complex) order parameter $\psi$ and for the electromagnetic vector potential ${\bf A}$ (the scalar potential is usually eliminated through an appropriate choice of gauge).

Numerical approximations to the TDGL equations have been derived using both finite element [1] and finite difference [2,3,4,5,6,7] methods. Most physical applications use a specific finite difference method that we will refer to as $\psi U$-method. The unknowns in the $\psi U$-method in two spatial dimensions are the order parameter $\psi$ and two auxiliary fields, ${{\cal{U}}}^x$ and ${{\cal{U}}}^y$, that are related to ${\bf A}$ by

$${\cal{U}}^x (x,y,t) \exp \left ( - \imath \int_{x_o}^x A_x(\xi, y,t) ~d\xi \right )$$ = (1)

$${\cal{U}}^y (x,y,t) \exp \left ( - \imath \int_{y_o}^y A_y(x, \eta,t) ~d\eta \right )$$ = (2)

The point (x_o, y_o) is arbitrary, $\imath = \sqrt{-1}$. Such variables were first introduced in lattice gauge theories [8]. To our knowledge they were first applied to the TDGL equations by Liu et al [3]. The$% \psi U$-method has since proved useful in the numerical simulation of many superconductivity phenomena [3,5,6,7].

The TDGL equations coupled with Maxwell equations with the zero-scalar potential gauge choice lead to the following mathematical problem:

$${\partial_t \psi } -\frac 1\eta \left[ \left( -\imath\nabla-% {\bf A}\right) ^2\psi +\left( 1-T\right) \left( \vert\psi \vert^2-1\right) \psi \right] +\tilde{f}$$ = (3)

$${\partial_t {\bf A}} \left( 1-T\right) {\rm Re} \left[ {\overline \psi} \left( -\imath\nabla -{\bf A}\right) \psi \right] -\kappa ^2\nabla \times \nabla \times {\bf A}$$ = (4)

where lengths have been scaled in units of $\xi (0)$, time in units of $% t_0=\pi \hbar /(96k_BT_c)$, ${\bf A}$ in units of $H_{c2}(0)\xi (0)$ and temperatures in units of T_c. It has been assumed that the coherence length $\xi$ obeys $\xi (T)=\xi (0)(1-T)^{-1/2}$, where T is the temperature in units of the critical temperature T_c, and that the Ginzburg-Landau parameter $\kappa$ is independent of temperature. $\eta$ is a positive constant (a ratio of characteristic times for $\psi$ and ${\bf A}$),$% \tilde{f}$ a random force simulating thermal fluctuations, k_B the Boltzmann constant and H_c2 the upper critical magnetic field for type-II superconductors (see [9]). ${\rm Re}$ stands for ``real part of''.

Eqs. (3)-(4) are to be solved in a bounded domain$% \Omega$, complemented with initial conditions for $\psi$ and ${\bf A}$, together with the following boundary conditions:

- Boundary condition for ${\bf A}$:

A given applied magnetic field H_e in the z-direction, possibly time-dependent but spatially uniform, is assumed. Continuity of the field thus implies

$$B_z := \hat{\bf e}_z \cdot \nabla \times {\bf A} = H_e$$ (5)

- Boundary condition for $\psi$:

Zero supercurrent density perpendicular to the boundary is imposed, namely,

$$\hat{\nu} \cdot \left( -\imath \nabla -{\bf A}\right) \psi =0$$ (6)

where $\hat{\nu}$ denotes the unit normal to the superconductor-vacuum interface. This automatically implies that the normal current perpendicular to the boundary also vanishes, since the total current across the superconductor-vacuum interface is zero. To see this, let ${\bf J}_s=(1-T) {\rm Re} [{\overline \psi} \left( -\imath\nabla -{\bf A}\right) \psi]$ denote the supercurrent density and ${\bf J}_n = -\partial_t {\bf A}$ the normal current density. Rewrite Eq. 4 as

$${\bf J}_n + {\bf J}_s = \kappa ^2\nabla \times \nabla \times {\bf A}$$

Projection of this equation along the normal $\hat{\nu}=(\nu_x,\nu_y,0)$ leads to

$$\hat{\nu} \cdot {\bf J}_s +\hat{\nu} \cdot {\bf J}_n = \kappa^2 (\nu_x \partial_y - \nu_y \partial_x) B_z$$ (7)

Since $(\nu_x \partial_y - \nu_y \partial_x)$ is nothing but the tangential derivative, and considering that the applied field H_e is uniform, the right-hand side of (7) identically vanishes showing that the total current across the boundary is zero.

For later use, we define the magnetization M_z as (see, e.g., [9])

$$M_z(t)=\frac{\int \left( B_z(x,y,t)-H_e\right) dx~dy}{4\pi \int dx~dy}$$