{'fraction_non_alphanumeric': 0.02686553952376737, 'fraction_numerical': 0.006904487917146145, 'mean_word_length': 4.81296629791613, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-Dependent Ginsburg-Landau Simulations of Superconducting Vortices in Three Dimensions

Here we describe a development of computer algorithm to simulate the Time Dependent Ginzburg-Landau equation (TDGL) and its application to understand superconducting vortex dynamics in confined geometries. Our initial motivation to get involved in this task was trying to understand better our experimental measurements on dynamics of superconductors with vortices at high frequencies leading to microwave stimulated superconductivity due to presence of vortex (Lara, et al., Scientific Reports, 59187 (2015)).


Time-dependent Ginzburg-Landau simulations
High frequency dynamics in superconducting and magnetic vortex states has been subject of recent activities both from experimental and numerical points of view. While Landau Lifshitz Gilbert equation is simulated numerically to obtain static magnetization and spin wave modes [1], the Time Dependent Ginzburg-Landau equation is usually resolved numerically to understand the static and dynamic properties in the superconductong vortex state [2].
The initial attempts to simulate the vortex behavior as a function of time using the TDGL equation date from the early 1990s. The most popular method to tackle this problem is the link variable method [3,4,5,6,7,8], which uses finite differences (i.e. regular grids and derivatives are approximated by subtractions between neighbors). Some attempts have been made using finite elements [9], but the extra complications that this method has makes it much harder to implement, although it allows to better represent curved shapes.
The work presented below is based on the numerical approach exposed in [10], with important extensions to include arbitrary shapes and calculations in 3D. First, some general concepts about numerical solution of differential equations are explained, followed by a more detailed description of the particular method for solving TDGL.
2 Numerical method Numerical solution of differential equations There are three types of differential equations to represent physical phenomena: • Time dependent equations. To solve these equations, one integrates numerically some variable in time. An example of this is the classical three body problem, where one integrates in time the position of a body in the presence of two other bodies, taking into account only their graviational interaction. Time integration uses methods like the simple but effective Euler method, or the more complicated Runge-Kutta.
• Space dependent problems. This type of problems require integrating a differential equation independent of time, such as Laplace's equation for the electrostatic potential. When there is a spatial dependence, one resorts to methods like finite differences or finite elements. Typically an iterative method is used to find a convergent solution.
• Space and time dependent methods. This is the present case of TDGL equation, and it is a combination of the former two. One typically finds by spatial integration the time derivative of the variable of interest (in our case, dΨ(t) dt , in turn dependent on Ψ(t)) and integrates it in time. With the updated value of Ψ(t), dΨ(t)/dt is calculated again. The mesh has a constant spacing (not necessarily the same spacing in different directions).

Link variable method
We have used the usual approach to solving the TDGL equations known as link variable method. Normally in finite difference methods the interaction between neighboring cells is expressed (partially or completely) via the spatial derivatives. The interesting feature of the link variable method is that it relates the interaction between neighboring cells with the magnetic flux enclosed by each set of cells forming a closed path, which is also dependent on an external applied field. This is done by introducing the auxiliary "link variables" in all directions, of the form: Az(x,y,ζ,t)dζ being x, y, z the spatial components, and ξ, η, ζ are auxiliary integration variables in the same directions as x, y, z (not to be confused with the superconducting coherence length, ξ). x i are the positions where Ψ(x, y, z) is calculated, and x 0 , y 0 , z 0 are arbitrary positions that in the end up cancelling out after doing line integrals along closed paths.
Given a plane, XY for example, the rectangle formed by the cells in (i, j), (i + 1, j), (i, j + 1) and (i + 1, j + 1) will hold a circulation of the vector potential that can be related to magnetic flux at that position inside the superconductor via link variables: From such a rectangle the discretized distribution of magnetic field can be found in the x = i, y = j, z = k coordinates: ,k U y i,j,k = e −iH(i,j,k)∆x∆y and a similar calculation for the other two planes Y Z and ZX.

Discretization of the TDGL equations
We start with the TDGL equations in the following dimensionless , for zero scalar potential gauge: Here, lengths are scaled by the the coherence length at zero temperature The vector potential A is expressed in units of H c2 · ξ 0 . The coherence length dependence on temperature is assumed to be well described by ξ(T ) = ξ 0 / 1 − T /T c . κ is the Ginzburg-Landau ratio which decides if a superconductor is type I or II. η is a positive constant with value η = To reproduce the configuration of our experimental setup, we need to be able to have magnetic fields both parallel and perpendicular to the plane. With the more usual 2D configuration, one can only apply fields perpendicular to the plane, since to consider the effect of magnetic flux, it has to actually go through some closed path. If one tries to apply a field parallel to the plane (XY , for example) in a 2D simulation, because there is only one cell in the Z direction, no flux can enter. To have nonzero magnetic flux in-plane requires having more than one cell in the perpendicular to the plane direction. However, the more complicated geometry we want to reproduce cannot be achieved in 2D, since magnetic flux (used for boundary conditions) needs at least two layers to be accomodated in the simulated domain (flux is calculated via the circulation of link variables. Therefore, only flux perpendicular to a direction where the film has more than one layer of cells can be considered). Figure 1 shows more graphically the circulation of the link variables in every direction. Also, we need boundary conditions. The first condition is for Ψ at the boundaries. We will consider the case of samples isolated from their surroundings. This means that no supercurrent flows through the boundary, and is expressed as: A second boundary condition, for the magnetic field is used. It is through it that the magnetic field actually enters the sample. It has the form: Where n i is the unit vector in the direction of the component of the magnetic field H that we are interested in, parallel to a given surface, and H e,i is the value of the same component of the external magnetic field.
The equations integrated in time are four, one for the order parameter (first Ginzburg Landau equation), and the other three for the auxiliar variables known as link variables (see [10] for more details), that can be obtained by rearranging the three components of the second Ginzburg Landau equation for the vector potential. The boundary conditions, which depend on A, can also be written in terms of the link variables. Therefore, the vector potential A is not explicitely solved, and it is not needed either. Its value (more precisely, its circulation) is included implicitly in the link variables, which are used to recalculate the order parameter in each step of the simulation, making it unnecesary to explicitly solve A.
An important detail encountered when transitioning from 2D to 3D is described next. To solve the equations for the link variables in the boundary in 2D one uses the circulation (see figure 2 for more clarity): This set of four cells is what we will refer to as "square loop". By multiplying both sides by de appropriate link variable (conjugated or not) one can isolate the link variable at the edge of the sample for each of the smallest possible square loops to calculate its value from the value of the other three and the applied field.

Passing from 2D to 3D
If one tries now to do the same thing in 3D, there is an inconsistency, because at edges the unknown boundary link variable can take two different values, depending on whether it is calculated from the one or the other loop in which it takes part (figure 3). The workaround that we came up with, and works well, is to forget about calculating boundary link variables, and just force the flux through the boundary square loops. The boundary link variables are given arbitrary values. This works because these link variables really aren't used to calculate anything.

Hx
Hz L X L Z Figure 3: Link variables at the edges (in this figure, between the two green cells) have an uncertain value.
Having overcome this problem, everything works well. Then, one would expect to see tilted vortices under an inclined field, as indeed happens. Figure  4 shows two isosurfaces of early stages of the development of a vortex state under an inclined field in a cube in a simulation with κ = 2 and 50x50x50 cells. After some simulation time the system stabilizes and tilted vortices are observed (figure 5). In these scalar fields |Ψ(x, y, z)| one needs to select Isosurfaces at low values, for example, plot only the points with |Ψ(x, y, z)| = 0.1 to clearly see the vortices. In this case, a magnetic field from one corner to the oposite ((1, 1, 1) direction) is applied. Even without considering demagnetizing effects, vortices start at one border of the sample and end at another one which is not the opposite. The boundary conditions force the direction of magnetic field at the boundaries, and far from them the vortex curves to adapt. Thus, shorter vortices (closer to corners) are also more curved. Calculated isocurrents curl around vortices. Instead of plotting all (or a subset) of the arrows at each cell, as in figure 5 a), the supercurrent direction can be easier to understand if we only plot some isolines, as in figure 5 b). These isolines clearly show the supercurrent screening taking place in the direction perpendicular to the applied field ( figure  4 a), i.e., the streamlines fall in planes perpendicular to the (1, 1, 1) direction.

Results
Increasing field part of hysteresis loops Figure 6 shows field ramps and the global magnetization response. As expected from a superconductor, magnetization is negative for positive fields, creating a field opposing the external field, so that the internal field is zero in the superconductor. In these time dependent calculations, the speed at which the applied fields change is relevant. In the figure it is easy to see how ramping the field more slowly (more iterations) gives a curve with more jumps. This is so because vortices have more time to enter the sample discretely, overcoming at certain moments, many at once, boundary energy barriers. Faster ramps are less controlled, and vortices enter in a more disorganized manner, leading to a smoother curve, where jumps corresponding to entering vortices are smoothed out. The presence of jumps is also very dependent on the size of the sample. Their relative importance is larger in smaller samples. On the other hand, large samples give the more typical smooth "field penetration regions" found in magnetization measurements.

Steady flow of vortices between permanent magnets
In this numerical experiment, we verify the idea of converting magnetostatic energy into a steady flow of vortices in one direction. The setup is simple: a superconducting strip, with the appropriate width to hold only one row of vortices across its length, and two permanent magnets, each one in an end of the strip, with magnetization pointing in opposite directions, perpendicular to the plane of the strip ( figure 7).
The magnets are modeled by boundary conditions: a non uniform external field, nonzero at the strip ends (only in a few cells near the center, not the whole lateral side), and zero everywhere else.
Of course, this is an approximation, but the fast decay of 1/r 3 typical of dipolar magnetic field suggests that is still valid. If both magnets are equally strong, vortices and antivortices are created in the same amount, and they meet at the center of the strip, as long as they don't need to travel a large distance (if they have, they will exit the superconductor throught the long sides). Since they have different signs, they annihilate. What can be done now is tuning the values of the strength of the magnets, making one higher than the other. In this case, the point where the vortices meet shifts laterally, closer to the weaker magnet. This is so because the stronger magnet manages to get the vortices it creates sooner into the strip than the weaker magnet. We can even get to the point where there is only a steady flow of vortices from one magnet to the other (when one of them is not strong enough to create vortices, but still stronger than the zero field at the long sides of the strip). Practically, this would mean to tune the values of very small magnets, which does not seem very feasible or practical. Instead, we can apply a perpendicular to the plane magnetic field, which will oppose the field created by one of the permanent magnets, and reinforce that of the other magnet. One needs to be careful to apply low values to this external field, since if it is too large it can lead to the creation of vortices itself. Instead, by keeping its values low, we are creating a "slope" which will determine where vortices of different signs will annihilate. Figure 8 shows how the external field can influence this.
This idea could be implemented for different purposes, such as very precise magnetic field detectors, since moving vortices induce electric fields, which could be tracked by some array of sensors over the strip, to relate the annihilation position of vortices to the external field. Also, electrically charged nanoparticles could be transported by these moving vortices, until they are deposited in the annihilation place, where no longer electric field is induced after annihilation. Figure 7: A superconducting strip with two magnets at its ends.

Including holes
Another feature that could be added to the program, which fell out of the scope of this thesis due to lack of time to develop it completely, is being able to include holes in the calculations, also of arbitrary cross section but, at least to start with, "drilled" completely perpendicular to the boundary. Again, the idea shown in [10] could be adapted to 3 dimensions. An extra step would be necessary for including multiple holes: since for each hole, the magnetic field in that hole needs to be tracked for imposing boundary conditions, a function could scan the whole sample, to detect which cells belong to what hole. A more or less complicated procedure should be thought of to identify what is a hole, and what cells form it, based on the neighboring cells which have already been scanned and are empty. After that, the cells forming the boundary of each hole should be associated a "type" of cell, depending on where the neighboring empty cell lays with respect to them. With that, and knowing the field already trapped inside each hole, the boundary conditions could be calculated. Whenever a vortex falls inside a hole, the magnetic field in it will experience abrupt jumps. Including holes in the calculations would allow to simulate effects like periodic pinning, to calculate matching fields, as well as the interaction of vortices with defects of different shapes, for example, as is usually done in the context of ratchet effects with triangular defects [11,12].

Rectification effects
Rectification effects of the vortex lattice have attracted a lot of attention, specially in superconducting films with periodic pinning arrays, typically of triangular elements [11,12,13]. In general, nonlinear dynamics in the vortices in asymmetric potentials, under high frequency microwave or combined DC and AC (microwave) drives has recently returned into the focus of interest of the reseachers in vortex physics and applications. [14,15,16] The idea of rectification is to create a potential landscape that favors the motion of vortices in one direction with respect to the opposite direction. These systems are of interest both from the applied and fundamental point of view, since they allow to design mechanisms to control the vortex flow artificially. Some other works focus on motion of vortices in structures with ratchet-like shape [17,18]. Also, rectification effects have been found in systems that do not posses obvious ratchet potentials [19]. Those AC current induced DC voltages appear at high enough frequencies due to vortex rectification by each of the oppositely situated surface barriers [20]. Most of the work so far considers the motion of vortices already present in the sample, whether it is a uniform film with pinning centers, or patterned superconducting elements. In this case we have applied our program to the case of vortices entering and exiting the sample. Rectification effects in this case occur due to the very geometry of the sample and the energy barriers that vortices need to overcome in order to enter inside. In this set of simulations (figure 9 we have put our program to the test to verify two facts: First, adding rectangular notches to the boundary facilitates the entrance of vortices through them. This had been already found for single triangular notches in circular samples [21] and one side situated triangular notches in rectangular samples [22]. Second, a rectangular shape with the notches in the long sides, in the presence of a perpendicular to the plane field (which slowly oscillates in value), a larger number of vortices enter through these sides than they exit during an oscillation period. This means that the combination of this geometry with the notches is enough to create a net vortex flux for an external field of zero average value in time. In practice, if such mechanism proved to be true, a voltage could be measured between the sides of the sample.

Motion of vortices under a stationary current
This simulation reproduces the vortex flow under the influece of a current flowing through a superconducting strip. A simple approximation to simulate a current is to, instead of using any current at all, modify the magnetic fields at the boundaries parallel to the flow of current. In one of them, an extra field should be added, and in the other, the same value should be subtracted, since a current flowing uniformly will create a magnetic field around it. This simplification only works for 2D simulations, due to symmetry reasons, since for 3D simulations, the curvature of the magnetic field lines should be taken into account, but in 2D they are perfectly perpendicular to the plane, since this case is equivalent to a sample infinitely long in the z axis. This method has been used before [23]. Figure 10 shows snapshots of a simulation of the vortex motion under a current parallel to the long sides of the strip, every 20 timesteps.

Spontaneous nucleation of vortices
Usually it is of interest to create vortices in a controlled manner with an applied field, but that is not the only way. Spontaneous formation of vortex antivortex pairs has been observed due to the stray fields of magnetic elements placed on top of superconducting films [24]. A superconductor cooled from above the critical temperature can also develop pairs of vortex antivortex inside [25,26], still having a total zero magnetic flux. This configuration is unstable, and vortices and antivortices will find it favorable to approach and annihilate each other to reduce the magnetic energy. Next we explore the behavior of the program to Figure 11: Spontaneous vortex antivortex pairs formation starting from a random order parameter.
calculate the spontaneous formation of vortex antivortex pairs in a superconductor in the absence of an external magnetic field, starting from a random distribution of the order parameter centered around |Ψ| = 0.1. This is shown in figure 11 for a simulation with 1000 x 1000 cells (0.5ξ per cell) in the plane. This replicates the situation of abruptly quenching superconductivity by cooling down a sample in the normal state. After the vortices and antivortices are formed, they start annihilating each other if they happen to be close together. After several annihilations, there will be some vortices and antivortices left, far apart enough to not feel each other's attraction, and the final state will have vortices inside, but still a zero (or close to zero, due to vortices exiting through the boundary) total magnetic field inside.

Perspectives
The numerical solution of the time dependent Ginzburg Landau equation has been an active topic since the early 1990s, until the present moment. Currently there are a number of researchers using these simulations to study theoretically the behavior of increasingly complex superconducting systems. We have focused on a rather simple system, a thin uniform film, under conditions reproducing our measurements (crossed alternating and constant magnetic fields), as well as to predict what would happen under certain conditions, paying attention to the motion of vortices. Interesting results on static 3D simulations of vortices in parallel magnetic field can be found in [27]. Application of the developed method to high frequency response of the vortex system related with microwave stimulated superconductivity can be found in [2,28]. The use of these simulations is far from being exhausted. Currently an important effort is being made toward understanding multiband superconductors, such as MgB 2 , with several order parameters using numerical methods more complicated that the one presented here, where several order parameters are present and coupled among themselves [29]. Also, as has already happened in the field of electronic circuits design, especially in circuits working at high frequencies, simulations can become a must-have tool to design better performing circuits, without the need to build lots of circuits by trial and error before reaching an optimum configuration. We expect that the same situation will be reached soon with superconducting circuits, which are gaining more and more attention for applications such as superconducting antennas and quantum computers [30].