Brief numerical analysis of (3+1) Ginzburg-Landau equations

In this contribution, we show the implementation of the Link-variable method for solve the complete set of acopled non-linear time dependent Ginzburg-Landau differential equations in a three-dimensional homogeneous and isotropic mesoscopic superconducting system. In this case, the sample is immersed in an external magnetic d at zero applied current. The effects of demagnetization are taken in count and we show the order parameter and its phase in zero field cooling and field cooling process. This numerical analysis shows good results when this solution is applied to a superconducting cubic sample.


Introduction
The superconducting systems has been one of the most studied topics in the last years, since its discovery. This is due to the interesting electronic, thermodynamic and magnetic properties that systems exhibit, this added to the possible technological applications that they could have in the following generations. The main properties are the Meissner-Oschenfeld effect, which describes the shielding of the applied external field, in addition to the conduction of electric currents, without Ohminc loss, for a certain critical temperature characteristic of the sample T c [1][2][3][4]. However, in the last years, this branch of condensed matter has been revived quite a bit, due to the discovery of high-temperature superconductors critical in 1986 and more recently the study of multibands, multicomponents, topological, and type 1.5 superconductors, which they are part of the so-called unconventional superconductors, because they cannot be studied using the Bardeen-Cooper-Schrieffer (BCS) theory [5][6][7].
With this, one of the theories that have had extensive applications in conventional and unconventional superconductors, is the time dependence Ginzburg-Landau (TDGL) theory, which in general are the result of the application of the second-order phase transition theory by the Gibbs functional, obtaining a system of coupled nonlinear equations for the order parameter and the potential vector is due to non-linearity and the complexity of the differential equations obtained, with which in general (2 + 1) systems are studied, and even in these systems, the computational and storage cost is too great. The study of superconducting systems has been revived in recent years, thanks to the new phases discovered, such as topology [8], multibands [9][10][11], multicomponents [12,13], Mott and other effects in which the principles of the BCS theory are not applicable. With this, one of the main theories that have been used in the study  [14][15][16]. Through the application of TDGL, the dynamics of vortices and the main thermodynamic properties are generally studied [14]. However, the computational cost in mesoscopic systems is very high [17], and in (2 + 1) systems it involves long computation and storage times, in such two-dimensional systems, the real effects on the sample are not studied, such as the demagnetization effects that occur in three-dimensional samples, which of course has profound effects on the dynamics of the vortices and other properties [2,17], i.e., the real effect that the geometry of the sample has on the main properties of the superconductor. With this, in the present work, we show the computational implementation for a three-dimensional superconducting system, the work is written in such a way that its implementation in Matlab, is almost immediate because we will use the notation of the colon operator commonly used in that program. Finally, we will present the order electronic superconducting density (order parameter and its phase) for both field cooling (FC) and zero field cooling (ZFC) process [18].

Physical model
In the TDGL model, the superconducting order parameter is a complex pseudo-function ψ(r, t) = √ n s e iΦ , where n s , is the density of the super-electrons and Φ its phase. So, the time-dependent Ginzburg-Landau equations in dimensionless units for the ψ and the vector potential A, in the gauge of zero electric potenctial are given by [2,3,15,16].
In the Equation (1) and Equation (2), Re is the real part, ψ and ψ * are presented in units of (α/β) 1/2 , distance in units of coherence length ξ, time in units of π /96K B T c , A in units of ξ times the second critical field H c2 , temperature T in units of the critical temperature T c . κ is the Ginzburg-Landau parameter. We use the superconducting-dielectric boundary condition for the order parameter and magnetic field, they are given byn · (−i∇ + A)ψ| n = 0, and ∇ × A = B at the surface of the sample.n is the normal vector to the surface [15,16].

Mathematical model
The geometry of the problem that we investigate is illustrated in Figure 1(a). The superconducting domain covers the parallelepiped of high a × a × c. Due to the demagnetization effects, we need to consider a larger domain of dimensions L × L × z. (For more details see Ref. [17]).  Them first we will start the initial conditions, in the matrices ψ(1 : N x + 1, 1 : N y + 1, 1 : N z + 1) = 1.0 (Meissner-Oschenfeld state), A x = zeros(N x, N y + 1, N z + 1), A y = zeros(N x + 1, N y, N z + 1) and A z = zeros(N x + 1, N y + 1, N z) (Figure 1(b)), and the the link variables (Equation (3)).

Link variables (3+1) dimensions
Now we will use, the Link variable method [2,3], for establish the discrete form in this base. For notation we will use ψ n (i, j, k) is the n the time dependence and (i, j, k) the spatial domain (Equation (4) to Equation (14)); with this.
Now for the numerical implementation, we will use the M atlab notation for the internal box (i, j, k) = (x1 + 1 : x2, y1 + 1 : y2, z1 + 1 : z2). Now for the second TDGL equation, we will use the following change of variable.
The loops for the super-currents J x (i, j, k) = (1 : N x, 2 : N y, 2 : N z), J y (i, j, k) = (2 : N x, 1 : N y, 2 : N z) and J z (i, j, k) = (2 : N x, 2 : N y, 1 : N z) and is imaginary part. Now, for the last term of the second TDGL equation we take the usual definition for the rotor , and using the forward numerical derivade.
And finally for the second TDGL.

Results
Now, for numerical simulations, we take L = 40ξ, z = 20ξ, a = 10ξ and c = 5ξ, the size of the grid ∆x = ∆y = ∆z = 0.25, the Ginzburg-Landau parameter is κ = 1.0, the time step is ∆t = 10 −5 and the field step ∆H = 10 −3 . With this, the Figure 2(a) show the phase of the order parameter ∆Φ and the the Figure 2(b) show the superconducting electronic density |ψ| 2 for T = 0 and indicated magnetic fields. We observe the entry of vortices symmetrically into the system. As is well know, dark and bright regions represent values of the modulus of the order parameter (as well as ∆Φ/2π) from 0 to 1. The phase allows to determine the number of vortices in a given region, by counting the phase variation in a closed path around this region. If the vorticity in this region is N , then the phase changes by 2πN [19,20]. For low magnetic fields there is an even symmetry in the entry of the vortices, in accordance with the symmetry of the sample, however as the field increases more vortices enter the sample, breaking this symmetry. Now we see in Figure 3, we present the results for the vortex configuration in the zero field cooling process at H = 1.1 and indicates temperatures. In this case we appreciate non-stationary states with N = 11 vortices for all studied temperatures. This is a physically coherent result since the number of vortices depends on the external magnetic field. Therefore, we can observe that the symmetry of the superconducting nano-structures influences the magnetic properties of the sample. By improving the experimental process of cooling of these superconductors, the perspectives of their practical applications can be improved.