Giant paramagnetic Meissner effect in multiband superconductors

Superconductors, ideally diamagnetic when in the Meissner state, can also exhibit paramagnetic behavior due to trapped magnetic flux. In the absence of pinning such paramagnetic response is weak, and ceases with increasing sample thickness. Here we show that in multiband superconductors paramagnetic response can be observed even in slab geometries, and can be far larger than any previous estimate - even multiply larger than the diamagnetic Meissner response for the same applied magnetic field. We link the appearance of this giant paramagnetic response to the broad crossover between conventional Type-I and Type-II superconductors, where Abrikosov vortices interact non-monotonically and multibody effects become important, causing unique flux configurations and their locking in the presence of surfaces.

The diamagnetic Meissner effect is one of the hallmarks of superconductivity, where applied magnetic field is ideally screened out of the superconductor when cooled below the critical temperature T c . However, many field-cooled experiments on various materials over the past two decades detected a paramagnetic response, i.e. enhanced magnetic field inside the sample, usually referred to as paramagnetic Meissner effect (PME) or Wohlleben effect. The materials in question range from elementary ones such as Nb 1,2 , to much more complex high-T c cuprates [3][4][5][6][7][8][9] . One proposed explanation for the enigmatic origin of PME in cuprates is based on the d-wave symmetry of the order parameter and the idea that π junctions formed due to Josephson coupling between grain boundaries can result in spontaneous current loops with significant magnetic moments [10][11][12][13][14][15] . A much simpler and more general explanation is the compression and trapping of magnetic flux on cooling. Using this picture, Koshelev and Larkin 16 calculated the magnitude of PME in thin stripes of conventional superconductors, and concluded that its theoretical maximum is ~27% of the full Meissner response for the given magnetic field. Kostić et al. 17 performed a supporting experiment on bulk Nb, and further concluded that polishing the sample surfaces strongly alters the PME, i.e. that surface barriers for flux entry and exit play an important role.
The appearance of PME due to flux compression is easiest to understand in the case of mesoscopic samples, where the influence of confining boundaries is crucial. There, in analogy to surface superconductivity, during field-cooling the superconducting order parameter nucleates at the sample surface, and traps a multiquanta (giant) vortex inside the sample. Such large and compressed flux may lead to paramagnetic response, as first predicted by Moshchalkov et al. using self-consistent Ginzburg-Landau simulations 18 , and subsequently verified experimentally by Geim et al. 19 . A simple consideration shows that the paramagnetic moment strongly depends on the sample thickness, so that in very thick samples it appears only at very large fields and scales with penetration depth λ over lateral size of the sample, while in thin plates it can be significant and scales with λ over thickness 20,21 . Therefore, the enigmatic PME in high-temperature superconductors becomes intrinsic for thin mesoscopic conventional superconductors.

Results
We consider a larger than mesoscopic two-band superconducting slab of width w (w/λ ranges from 15 to 50 for the considered parameters) in a parallel magnetic field (see Fig. 1), and report particular behavior of the sample magnetization as a function of the applied field [M(H) loops]. We primarily focus on two-band materials, but our findings can be qualitatively extrapolated to systems with more than two bands. The calculations are performed within the two-component Ginzburg-Landau (TCGL) theory (see Methods), where we have cautiously set a sufficiently high temperature T to ensure the qualitative and quantitative validity of our predictions in the context of recent debates [40][41][42][43][44][45] , and we used full microscopic expressions of all coefficients in the theory [44][45][46] . TCGL theory then comprises eight independent parameters, namely, the Fermi velocities of the bands v 1 and v 2 , the elements of the coupling matrix λ 11 , λ 22 and λ 12 = λ 21 , the total density of states N(0) as well as the partial density of states of the first band n 1 (note n 1 + n 2 = 1), and finally T c , which sets the energy scale W T 8 7 3 By fixing the unit of length ζ 1 and normalizing the order parameters by W, the parameters v 1 and T c are fixed, and we are left with six parameters in the model: λ 11 , λ 22 , λ 12 , v 1 /v 2 , n 1 and N(0). Instead of choosing N(0), we opt to show the GL parameter of the first (stronger) band-condensate cW hev n N 1 , which is an indicator for the expected magnetic behavior of the sample.
In what follows, we consider an infinitely thick slab of width w = 80ζ 1 , and use periodic boundary conditions in the longitudinal direction (with size of the unit cell l = 120ζ 1 , see Fig. 1). Without loss of generality, we take for the remaining microscopic parameters of the sample: κ 1 = 1.5, λ 11 = 1.55, λ 22 = 1.3, λ 12 = 0.09 and n 1 = 0.48. Note that such choice of parameters does not correspond to any particular material, and is actually by no means unique -since our main study will concern the dependence of the magnetic properties on the ratio of the Fermi velocities v 1 /v 2 and temperature. We demonstrate these properties via calculated magnetization [M(H)] loops while adiabatically sweeping the magnetic field up and down.
In Fig. 2   ζ = / (since ζ 1 is fixed as the unit of distance) and we are thereby driving the system into the Type-II magnetic behavior (since   In increasing field, all calculated magnetization loops exhibit a superheated Meissner state above the thermodynamic critical field H c , where the superheating field H sh agrees very well with the seminal calculations of Matricon and Saint-James for H sh (κ) of single-band materials 49 . At H = H sh , superconductivity is either destroyed (for v 1 /v 2 < 0.34) or a jump to the mixed state occurs (for v 1 /v 2 > 0.34). The delimiting value of v 1 /v 2 = 0.34 exactly satisfies the condition H c = H c2 . In decreasing magnetic field, the superconductivity nucleates at the surface superconductivity field H c3 46 . Indeed, the nucleated states were only superconducting at the surfaces of the slab, with a large normal domain in the interior of the slab. For further lowered field and v 1 /v 2 < 0.34 the normal domain remains trapped until abruptly expelled from the sample at the expulsion field H e . This analysis confirms that magnetic response of the system for v 1 /v 2 < 0.34 is the one of Type-I superconductors, since typical superheating-supercooling picture holds there, H c2 is smaller than H c , and no vortices are found in the paramagnetic branch where flux was trapped upon nucleation of surface superconductivity. However, while decreasing field for v 1 /v 2 > 0.34, where consequently H c2 > H c , the normal domain becomes unstable at field H d but is not expelled; instead, it spreads into a vortex configuration, stable down to persistently lower expulsion field H e as v 1 /v 2 is increased. Simultaneously, flux trapping becomes notably more efficient, so that the vortex exit is hampered in decreasing field and paramagnetic response increases to its maximum at H e . This tendency continues up to v 1 /v 2 ≈ 0.53, for which paramagnetic response is almost an order of magnitude larger than the Meissner response at H = H e , and approximately 30 times larger than the largest theoretical estimate of paramagnetic response to date (scaled to the diamagnetic response at a given field, see Ref. [16]). We therefore refer to this property as giant paramagnetic response (GPR). For v 1 /v 2 > 0.53, the cumulative paramagnetic response is still very large but gradually decreases, and magnetization curves in decreasing field connect to zero without any abrupt flux expulsion. In other words, we approach the Type-II limit, in which magnetization is expected to hover around zero for descending field in the presence of surface barriers 48 . In Fig. 3, we summarize the observed maximal amplitude, Max(4πM/H) in the entire field range, and the total cumulative paramagnetic response, , as a function of v 1 /v 2 , extracted from Fig. 2.
Based on Fig. 3, we argue that the giant paramagnetic response is characteristic for superconductors between conventional Type-I and Type-II 39 . Namely, this pronounced paramagnetic response is exactly found for sample parameters between the line H c (T) = H c2 (T) and the line where long-range vortex interaction changes sign (determined by effective GL parameter κ * calculated after Ref. [28]), with a maximum found close to the parametric line where surface energy (σ SN ) of the superconductor-normal metal (S-N) interface changes sign (determining the change in the polarity of the short-range vortex interaction 28 ). For the microscopic parameters considered here, we show this domain in Fig. 4(a), as a function of v 1 /v 2 and temperature. To test our hypothesis further, we calculated an additional set of M(H) loops, shown in Fig. 4(b), for fixed v 1 /v 2 = 0.55 and varied temperature indicated by yellow arrow in Fig. 4(a). From Fig. 4(b), we confirmed exactly the same behavior of the loops and relationship of the giant paramagnetic response (GPR) with the delimiting lines of the critical domain: for T ≥ 0.98T c the expected response of a Type-II slab is found, for 0.98T c > T > 0.91T c the paramagnetic response in decreasing field increases when crossing the long-range vortex attraction line, and finally a pronounced paramagnetic response followed by a jump to the Meissner state is observed when crossing the σ SN = 0 line. Besides being useful for reaffirming our conclusions, this temperature dependence of the GPR can also be directly verifiable experimentally. Here, the considered samples are ideally clean, but even in realistic samples where flux trapping is present even at zero field, the rise and fall of GPR as a function of temperature will be easily observable in the above discussed scenario. Note that in general, changing any of the parameters can drive the in silico material across the crossover between the types of superconductivity, and thereby change the paramagnetic response. GPR is only sensitive on the regime of superconductivity the material is in, i.e., where the taken parameter set lies in the reconstructed Fig. 4(a) -Type-I, Type-II superconductivity, or in between.

Discussion
What is the underlying mechanism for the giant paramagnetic response? In simple terms, it is the facilitated trapping of magnetic flux in the crossover domain between Type-I and Type-II superconductivity, since vortices attract in the entire range of parameters where GPR is observed. However, GPR is found to be particularly large for σ SN > 0, where vortex-vortex interaction is purely attractive and vortices should coalesce into larger normal domains. On the contrary, we observe that in decreasing field separate vortex cores are still visible, though strongly overlapping (see inset in Fig. 5). To clarify the dense vortex packing observed in Fig. 5, we calculate the multibody vortex-vortex interaction shown for several vortex clusters in Fig. 6. As a major surprise, we found that in this regime multibody vortex interactions become short-range repulsive and cause the formation of a vortex lattice. This is illustrated in Fig. 6(a) (for v 1 /v 2 = 0.47 and T = 0.94T c , i.e. σ SN > 0), where we show the calculated vortex-vortex interaction as a function of the distance between vortices (labelled d), for a vortex pair, a vortex trimer, a vortex diamond-like cluster and a hexagonal vortex cluster. The pairwise vortex interaction is purely attractive, as expected, but in the other cases the short-range repulsion arises so that energetically favorable vortex-vortex distance arises in mid-range (note that this favorable distance closely corresponds to the average vortex distance observed in Fig. 5(b)). An insight into the physics of this short-range repulsive interaction can be achieved by analysing the superconducting state inside, for example, the hexagonal vortex cluster shown in Fig. 6. With this aim, we computed the maximum of the Cooper-pair density, n max , inside that cluster for each band-condensate separately, shown as a function of vortex distance d in Fig. 6(b). We reveal that the Cooper-pair density in the second condensate vanishes inside the vortex cluster at the vortex distance where short-range repulsion arises. Hence we can conclude that inside the vortex cluster the physics is driven by the other condensate, which has Type-II character, hence the repulsive interaction of vortices prevails at short distances. It is known that multibody vortex interactions are more complex than a simple superposition of pairwise interactions (see Refs. [50][51][52][53]), but it has never been found before that multibody interactions can change the polarity of the vortex-vortex interaction. This is a key feature of the found mixed state for parameters of the system between σ SN = 0 and H c (T) = H c2 (T) lines in Fig. 4. In addition, we have plotted in Fig. 5(a) the number of vortices in the sample N v as a function of H in the downward branch of M(H) in Fig. 2 for v 1 /v 2 = 0.65 (in the Type-II limit) and v 1 /v 2 = 0.47 (inside the crossover region). The high retention of flux is clearly seen as a nonlinear behavior for v 1 /v 2 = 0.47 which contrasts the Type-II case in which N v is linearly decreasing towards the origin. We find that although the number of vortices in the states between σ SN = 0 and H c (T) = H c2 (T) lines slowly decreases with decreasing magnetic field, their favorable distance is approximately independent of field [see Fig. 5(b)]. This unconventional vortex state allows the penetration of the magnetic field in larger portions of the sample (inhomogeneous penetration, within but also between vortices), and clearly traps more flux than an ordinary vortex lattice, down to very low field -resulting in a more pronounced GPR. Due to the interlocking of vortices in this regime, the barrier for the expulsion of the entire vortex cluster in decreasing field corresponds to the Bean-Livingston barrier for a single vortex, which we confirmed by an independent calculation. Notice that as soon as the S-N surface energy changes sign, the barrier for single-vortex expulsion becomes nonzero at all fields. However, we have the simultaneous appearance of short-range vortex repulsion, which in effect diminishes the Bean-Livingston barrier and vortices are gradually expelled from the sample depending on their density and applied magnetic field. This manifests in the magnetization curves as a gradual decrease of the paramagnetic effect in decreasing field, down to zero for zero field. As the v 1 /v 2 ratio or temperature are further increased, vortices become increasingly repulsive and the paramagnetic response decreases to its conventional behavior for Type-II superconductors.
In summary, we revealed a possibility of giant paramagnetic response in slabs of multiband superconductors (to which many recently discovered metal-borides, iron-chalcogenides, iron-pnictides, belong), with magnitude similar or multiply larger than the Meissner response for the same applied magnetic field. We showed that such unique magnetic response occurs in the crossover region between conventional types of superconductivity, and is not captured by the standard textbook descriptions. On technological end, our findings open a new class of desirable materials which can be switched to either strongly enhance or fully remove the applied magnetic field while having low power consumption. Further work is needed to characterize the behavior of these materials under e.g. applied electric current and nanostructuring or downscaling.

Methods
In this work we used the two-component Ginzburg-Landau (TCGL) theory. In the TCGL framework, as given in Ref. 28, eight independent material parameters are needed for a system with both interband and magnetic coupling, namely, the Fermi velocity of the first band v 1 , the square of the ratio of the Fermi d Figure 6. (a) The vortex-vortex interaction energy, as a function of the distance between vortices, for parameters leading to pairwise vortex attraction (σ SN > 0, see open dots). Nevertheless, the short-range repulsion between vortices arises for clusters comprising more than two vortices (insets depict the cumulative Cooper-pair density distribution for the considered clusters). (b) Maximum of the Cooper-pair density, n max , inside the hexagonal vortex cluster for each band-condensate separately, shown as a function of vortex distance d between vortices. The shaded area delimits the short-range repulsion found for the hexagonal vortex cluster. , the elements of the coupling matrix λ 11 , λ 22 and λ 12 = λ 21 , the total density of states N(0) as well as the partial density of states of the first band n 1 (n 2 = 1 − n 1 ), and finally T c , which sets the energy scale W T 8 7 3 where j = 1, 2 is the band index, ( ) = / ( ) , and Γ = (N(0)λ 12 )/δ, with δ being the determinant of the coupling matrix, and S, S 1 and S 2 defined as in Ref. 41. The local magnetic field in the sample is denoted by  h and the external applied field by H.
Minimization of the free energy in Eq. (1) with respect to ψ j and A yields the Ginzburg-Landau equations: two for the order parameters ψ 1 and ψ 2 , and the equation for the vector potential (calculated from the supercurrent of the coupled condensate). Introducing the normalization for the order parameters by W, for the vector potential by A 0 = hc/4eπζ 1 , and for the lengths by v W 6 , the dimensionless TCGL equations are written as: where  denotes the real part of the expression. After the made choice of normalization units, we are left with six parameters: λ 11 , λ 22 , λ 12 , v 1 /v 2 , n 1 , and N(0). In our numerical experiment, the TCGL equations (2)-(4) were integrated self-consistently on a two dimensional grid with grid spacing a x = a y = ζ 1 , much smaller than any characteristic length scale at the considered temperature. The discretization was implemented by the link variable method which preserves the gauge invariance of these equations 54 . For the iterative solver, we combined a relaxation method with a stable and accurate semi-implicit algorithm 55 . Periodic boundary conditions were applied in the x direction whereas for the y direction we imposed Neumann boundary conditions at the superconductor-vacuum interface (for details of the numerical implementation, please see Ref. 54). Note that due to the infinite slab geometry, the surface magnetic field equals the applied one (the demagnetizing effects are negligible), and the simulation is effectively two-dimensional (in the (x,y) plane). ). The magnetic field is scaled to the thermodynamic critical field H c , for easier comprehension of the related physics. For calculation of H c for two-band superconductors, we refer to Ref. [28].