Examples

We detail in the following a few numerical examples. These depict some of the typical phenomena that occur in superconducting systems and how they are modeled by the described method. The reader is referred to [7] for an application of the model to the study of vortex arrays in superconducting thin films. It is also possible to extend the formulation to consider d-wave superconductors, as has been done in [10,11].

Consider first the case of a square sample with no hole (case I, see Fig. 3), with dimensions $32 \, \xi(0) \times 32 \,\xi(0)$. We use, with the units defined in Section 1, a_x = a_y = 0.5, $\kappa=2$, $\eta=1$ and T=0.5, with a noise constant of E₀=10^-5. With these definitions H_c1=0.04 and H_c2=0.5 for the bulk material, while H_c3=0.85 for a semi-infinite domain. Numerical limitations arise in the choice of the time step due to the forward-Euler treatment of the equations. A practical rule for time step selection is

$$\Delta t \leq \min \left \{ \frac{h^2 \eta}{4},\frac{h^2}{4 \kappa^2} \right \}$$ (33)

where

$$h^2 = \frac{2}{\frac{1}{a_x^2}+\frac{1}{a_y^2}}$$

In this case we choose $\Delta t = 0.015$, which satisfies stability. External field begins at H_e = 0 and is linearly incremented to H_e=1 along the 10⁶ simulated time steps. Variables are homogeneously initialized to a perfect Meissner state, $\psi(t=0) = 1$, ${\bf A}(t=0)=0$.

Figure 3: Scheme of the cases considered.

This case is particularly simple and runs at about 2500 steps per minute on a personal computer. The magnetization curve that is obtained can be seen in Fig. 4 (a). For applied fields smaller than H_e =0.203 the sample remains in a Meissner state, but at this field an instability develops that leads to the entrance of four vortices with the consequent jump in the magnetization. Similar vortex-entrance events occur at H_e = 0.216 (8 vortices), 0.252 (4 vortices), 0.270 (4 vortices), 0.288 (4 vortices), 0.312 (4 vortices), etc. In Figs. 5 we show the distribution of the modulus of the order parameter on the sample for several values of the applied field. It can be observed that the arrangement of the vortices is strongly affected by the finite size of the sample and its square symmetry. Larger samples allow for the obtention of hexagonal vortex lattices (an example can be seen in Fig. 6). Also notice that for H_e > 0.4 the superconductivity is strongly depressed in the sample's interior, but recovers near the boundary. This gradually leads, for $H_e \sim 0.5$ or greater, to surface superconductivity in a layer a few coherence lengths thick. The corners always remain the points where the order parameter is maximal.

Figure 4: Magnetization curves for the different cases considered.

(a)

(b)

The series of minima in the magnetization curves deserve special attention. Similar extrema were measured by Guimpel et al [12] and by Brongersma et al [13]. A related phenomenon was reported by Hünnekes et al [14]. In [7] it was argued that the minima reflect the magnetization behavior of the superconducting sheet at the sample surface and not rearrangements in the vortex lattice (which do occur). The simple case reported here confirms this argument, since minima extend far beyond H_c2 (H_e=0.5) and must thus come from a surface effect.

Figure 5: Colour graph of the modulus of the order parameter for different applied fields. Case I. H_e = 0.160 (a), 0.200 (b), 0.205 (c), 0.215 (d), 0.220 (e), 0.250 (f), 0.260 (g), 0.270 (h), 0.280 (i), 0.290 (j), 0.320 (k), 0.330 (l), 0.350 (m), 0.400 (n), 0.500 (o), 0.700 (p). Colours run from yellow ($\vert\psi\vert\sim 1$) to blue ($\vert\psi\vert\sim 0$).

(a)	(b)	(c)	(d)
(e)	(f)	(g)	(h)
(i)	(j)	(k)	(l)
(m)	(n)	(o)	(p)

Figure 6: Vortex arrangement at H_e=0.4 in a square sample with dimensions $96 \xi(0) \times 96 \xi(0)$.

Figure 7: Colour graph of the modulus of the order parameter for different applied fields. Case I, return (applied field descending from 1 to 0). H_e = 0.900 (a), 0.600 (b), 0.460 (c), 0.450 (d), 0.445 (e), 0.430 (f), 0.400 (g), 0.350 (h), 0.300 (i), 0.200 (j), 0.100 (k), 0.000 (l).Colours run from yellow ($\vert\psi\vert\sim 1$) to blue ($\vert\psi\vert\sim 0$).

(a)	(b)	(c)	(d)
(e)	(f)	(g)	(h)
(i)	(j)	(k)	(l)

In Fig. 4 (a) the hysteretic behavior of the system is clearly observed. Beginning at the calculated solution for H_e=1, negative increments in the applied field were imposed to show the hysteresis. Some snapshots of $\vert\psi\vert$ can be seen in Fig 7. For any given field, there are many more vortices in the sample when the applied field is decreasing than increasing. The lower (absolute) values in the magnetization suggest a smaller barrier for vortices to leave the system than the barrier for vortices to enter the system. The results are however not conclusive since in the simulation the rate of change of the applied field is rather high and time-dependent terms play a role. The three vortices in Fig. 7 (l), for example, leave the sample if the field is kept constant at H_e=0 during another 35000 time steps.

Consider now a hollow sample. Let the hole be a centered square, with dimensions $8 \xi(0) \times 8 \xi (0)$ (case II), or $16 \xi(0) \times 16 \xi(0)$ (case III). The same process of increasing and decreasing the applied field as before is conducted. The magnetization curves can be observed in Fig. 4 (b). Essentially the same structure as before arises, but at zero field in case II and III there remain fluxoids ``trapped'' in the sample. These fluxoids (5 for case II and 13 for case III) do not leave the sample if the simulation is continued keeping H_e=0 for as many as 10⁶ additional time steps.

Vortex arrangements are also different, as shown in Figs. 8-9. It is interesting to remark how the first four vortices that enter the system are ``captured'' by the hole in case II, and similarly for the first sixteen ones in case III. In Fig. 8 (b) the instant before the capture has been depicted, and similarly in Fig. 9 (f). The dynamics can be better grasped looking at animated sequences. Some are available at the same webpage mentioned above, or can be reproduced with the program in a few hours of CPU time.

We end up here the example section, since it is not the purpose of this article to discuss particular applications of the proposed method but rather to illustrate some physically interesting cases that can be simulated with the program.

Figure 8: Colour graph of the modulus of the order parameter for different applied fields. Case II. H_e = 0.202 [0] (a), 0.206 [4] (b), 0.210 [4](c), 0.230 [12](d), 0.250 [12] (e), 0.270 [20] (f), 0.300 [24] (g), 0.310 [24] (h), 0.320 [28] (i), 0.330 [28] (j), 0.340 [32] (k), 0.350 [32] (l), 0.37 [36] (m), 0.400 [40] (n), 0.450 [52] (o), 0.550 (p).Colours run from yellow ($\vert\psi\vert\sim 1$) to blue ($\vert\psi\vert\sim 0$). Between brackets the winding number along the external boundary is reported (above H_c2 the winding number is not evaluated because round-off pollutes the results).

(a)	(b)	(c)	(d)
(e)	(f)	(g)	(h)
(i)	(j)	(k)	(l)
(m)	(n)	(o)	(p)

Figure 9: Colour graph of the modulus of the order parameter for different applied fields. Case III. H_e = 0.200 [0] (a), 0.202 [0] (b), 0.218 [4](c), 0.220 [8](d), 0.230 [8] (e), 0.234 [16] (f), 0.250 [16] (g), 0.270 [24] (h), 0.300 [28] (i), 0.330 [28] (j), 0.340 [32] (k), 0.350 [32] (l), 0.400 [40] (m), 0.420 [44] (n), 0.450 [52] (o), 0.550 (p).Colours run from yellow ($\vert\psi\vert\sim 1$) to blue ($\vert\psi\vert\sim 0$). Between brackets the winding number along the external boundary is reported.

(a)	(b)	(c)	(d)
(e)	(f)	(g)	(h)
(i)	(j)	(k)	(l)
(m)	(n)	(o)	(p)