Electromagnetic superconductivity of vacuum induced by strong magnetic field: numerical evidence in lattice gauge theory

Using numerical simulations of quenched SU(2) gauge theory we demonstrate that an external magnetic field leads to spontaneous generation of quark condensates with quantum numbers of electrically charged rho mesons if the strength of the magnetic field exceeds the critical value eBc = 0.927(77) GeV^2 or Bc =(1.56 \pm 0.13) 10^{16} Tesla. The condensation of the charged rho mesons in strong magnetic field is a key feature of the magnetic-field-induced electromagnetic superconductivity of the vacuum.

The strong magnetic field causes exotic effects in hot quark-gluon matter, the well-known example is the chiral magnetic effect [3]. In the absence of matter the magnetic field background leads also to unusual effects like the magnetic catalysis [4], shift of finite-temperature transitions in QCD [5] and anisotropic conductivity [6].
It was recently suggested that a sufficiently strong external magnetic field should turn the vacuum into an electromagnetic superconductor [7,8]. The superconductivity emerges due to spontaneous condensation of electrically charged vector particles, ρ ± mesons, if the magnetic field exceeds the critical strength where m ρ = 775.5 MeV is the mass of the ρ meson and e is the elementary electric charge. In terms of the up-quark (u) and down-quark (d) fields the suggested ρ-meson condensate should have the following form [8]: where ρ = ρ(x ⊥ ) is a certain complex-valued periodic function of the coordinates x ⊥ = (x 1 , x 2 ) of a plane which is perpendicular to the magnetic field B = (0, 0, B).
The superconducting vacuum should have many unusual features. Firstly, no matter is required to create the superconductor so that the electromagnetic superconductivity appears literally "from nothing". Secondly, the superconductivity is anisotropic so that the vacuum acts as a superconductor along the magnetic-field axis only. Thirdly, the superconductivity is inhomogeneous because the ρ-meson condensate is not uniform in the B-transverse directions due to the presence a new type of topological defects, the ρ vortices. Fourthly, the net electric charge of the superconducting vacuum is zero despite of the presence of the charged condensates (2) [7,8].
The spontaneous generation of the ρ-meson condensate (2) -which plays a role of the Cooper pair condensate in the conventional superconductivity -is the key feature of the vacuum superconductor mechanism [7,8].
The appearance of the ρ-meson condensate was found analytically in a phenomenological model based on the vector meson dominance (VMD) [7], in the Nambu-Jona-Lasinio (NJL) model [8], and in holographical approaches based on gauge/gravity duality [9]. We use numerical simulations of the lattice gauge theory to demonstrate that strong magnetic field indeed leads to emergence of the superconducting condensate of the charged ρ mesons.
In QCD the charged ρ meson field is identified with ρ µ =ūγ µ d. The condensation pattern (2) corresponds to the condensate of the ρ mesons with the spins aligned along the axis of the magnetic field, with the s z = +1 projection of the spin onto the z axis. It is convenient to introduce two combinations of the negatively-charged ρ-meson fields 2 , which correspond to the spin projections s z = ±1, respectively. Indeed, according to the simplified arguments of Ref. [7], the invariant masses M n,sz of the ρ mesons states in the magnetic field B should behave as follows, with the nonnegative integer n and the spin projection onto the z axis s z = −1, 0, +1. The ground state is identified with the quantum numbers n = 0 and s z = +1, and the charged ρ mesons should get condensed, M 2 0,+1 < 0, if the magnetic field exceeds the critical value (1).
The simplest way to check numerically the possible appearance of the superconducting condensate (2) is to calculate the equal-time correlation function along the direction of the magnetic field, where the separation in the transverse coordinates of the two probes is set to zero, x ⊥ = (0, 0). The long-distance behavior of the s z = +1 correlation function (5) should expose the expected emergence of the condensate (2) due to the factorization property: 2 One can equivalently work with the positively-charged ρ-meson fields,d(γ 1 ± iγ 2 )u. The magnitudes of the positive and negative condensate are equivalent because the vacuum state is electrically neutral. Our results on vacuum condensation are the same for positive and negatively charged operators.
while the ρ mesons with the opposite orientation of the spins, s z = −1, should not be condensed 3 : We calculate the correlation functions (5) numerically, using lattice Monte-Carlo simulations of quenched SU (2) lattice gauge theory following numerical setup of Ref. [6]. The quark fields are introduced by the overlap lattice Dirac operator D with exact chiral symmetry [10]. The correlation function (5) is a linear combination of the current-current correlators in the vector meson channel. The vector correlator is represented in terms of Dirac propagators in fixed background of Abelian and non-Abelian gauge fields and is then averaged over an equilibrium ensemble of non-Abelian gauge fields A µ : where S Y M [A µ ] is the lattice action for gluons A µ . A uniform time-independent magnetic field B is introduced into the Dirac operator D f for the flavor f = u, d in a standard way by substituting the su(2)-valued vector potential A µ with the u(2)-valued potential A µ ij → A µ ij + δ ij q f F µν x ν /2, where q f is the electric charge of the corresponding quark, q u = +2e/3, q d = −e/3, and i, j are color indices. We also introduce an additional twist for fermions in order to account for periodic boundary conditions in spatial directions [11,12]. For technical reasons the bare quark mass m 0 is fixed at a small value am 0 = 0.01, where a is the lattice spacing. The vector correlation functions depend very weakly on the bare quark mass if they are calculated with the help of the overlap Dirac operator [13].
Our numerical approach is done in two complimentary ways. Firstly, we study in details the superconducting condensate for eleven values of the magnetic field B at a relatively small symmetric 14 4 lattice at a single lattice gauge coupling β = 3.281. These parameters correspond to the physical volume of the lattice is L 4 = (1.44 fm) 4 and the lattice spacing a = 0.103 fm [14].
Then we use an heuristic fitting method to find the superconducting condensate without taking a long-range limit (6) because the factorization (6) does not work well in too small volume. Secondly, we consider a set of lattices of various physical volumes L 4 and four values of magnetic field strengths B, and then utilize a conventional fitting procedure to extract the condensate η = η(L). The extrapolation to the infinite volume, L ≡ la → ∞, shows that these two methods give the same results. In both approaches the ultraviolet lattice artifacts are reduced with the help of the tadpole-improved Symanzik gluon action [15].
In our first approach we use 30 configurations of the gluon gauge field for each value of the background magnetic field. The periodicity of the lattice leads to the quantization (k ∈ Z) of the magnetic field, because of the requirement In Eq. (9) the integer k = 0, 1, . . . , L 2 s /2 determines the number of elementary magnetic fluxes which pass through the boundary of the lattice in the (x 1 , x 2 ) plane.
The maximal possible value of the fluxes k = l 2 /2 = 98 corresponds to an extremely large magnetic field with the magnetic length L B ∼ (eB) −1/2 being of the order of the lattice spacing, L B ∼ a. In order to avoid associated ultraviolet artifacts, in our simulations we limit the maximal value of the fluxes by k max = 10 l 2 /2, so that our maximal magnetic field, eB max = 3.54 GeV 2 is much larger than the expected critical magnetic field (1).
In Fig. 1 we show correlator (5) for a few relatively small values of the magnetic field. As we have anticipated, the s z = ± 1 correlators split in the external magnetic field [in Fig. 1 we show both spin orientations for B = B min ]. The observed splitting reflects the change in the relevant lowest energies (4), M 2 ± = m 2 ρ ∓ |eB|, so that the expected hierarchy of the masses, M + <m ρ <M − , is encoded in the slopes of the correlators The splitting of the s z = ±1 masses was also found in SU (3) lattice gauge theory at weak magnetic fields [16]. We have checked that the no-condensation property (7) for the s z = − 1 correlator G − is valid for all studied values of the magnetic field. In the weak field region, B < B c , the long distance behavior of the correlator G + is expected to be proportional to the function e −µ|z| , or cosh[(|z| − L/2)µ] in a finite volume with periodic boundary conditions. Here µ is a massive parameter and L = la is the physical lattice size.
We have found, however, that our numerical data for the correlators (5) are consistent with the cosh-like behavior only in a very narrow interval of the coordinate z. Therefore, in at B < B c we used a heuristic fit function which reduces to the exponential ("no-condensate") behavior in the thermodynamic limit L → ∞. In Eq. (11) µ > 0, γ > 0 and C weak > 0 are the fitting parameters. The function (11) which works surprisingly well at weak values of the magnetic field (9) with k = 0, 1, 2. The values of χ 2 /d.o.f. are shown in Table 1. The heuristic function fits nicely our numerical data for the G +  Fig. 1.
We have found that in the high-strength region, eB > 1 GeV 2 , the numerical data for the correlator G + (z) can well be described by another heuristic function for all separations z excluding the ultraviolet region with z a (and its periodic mirror). In Eqs. (12) and (13) µ > 0, η > 0 and C high < 0 are the fitting parameters. The G + correlator in the strong magnetic field region and the corresponding best fits (12) are shown in Fig. 2.
The parameter η in the fitting function (12) plays a role of the charged ρ-meson condensate, η ≡ | ρ |, because in the thermodynamic limit, L → ∞, the function V (z) reduces to an exponential e −µ|z| so that lim z→∞ lim ls→∞ G fit,strong Thus, the fits of the G + correlators provide us with the values of the condensate of the charged ρ mesons, Fig. 3. The corresponding values of χ 2 /d.o.f. are shown in Table 1. In the weak field domain the fitting does not converge properly due to presence of flat directions in the fitting parameter space.
In the weak field region the ρ meson condensate η = | ρ | vanishes, while at higher values of the magnetic field the condensate deviates spontaneously from zero signaling the presence of the superconducting phase. We find that the dependence of the ρ-meson condensate on the magnetic field can be described by the linear function:  The best linear fit (shown by the dashed line in Fig. 3) of the condensate allows us to determine the critical magnetic field of the insulator-superconductor transition, or B c = (1.56 ± 0.13) · 10 16 T, in satisfactory agreement with the theoretical relation (1), eB c = m 2 ρ , for the quenched mass of the ρ meson in SU (2) lattice gauge theory, m ρ ∼ 1.1 GeV 2 [17]. The prefactor in Eq. (15) is C ρ = 7.5(5) MeV. Notice that in our quenched model the exponent in Eq. (15) is ν = 1 while the mean field methods both in the bosonic VMD model [7] and in the fermionic NJL model [8] predict ν = 1/2 (so that theoretically η ∼ √ B − B c for B B c ). The behaviour of the fitting parameter µ in the fitting functions at both sides of the critical phase transition (11) and (12), (13), are shown in Fig. 4.
The unusual forms of the fit functions (11) and (12) is used to absorb the finite volume effects. In order to support our small-volume results we have performed an infinite-volume extrapolation of the condensate η obtained at larger lattices (l = 17 . . . 21 with L ≈ 1.65 . . . 2.2 fm). The condensate was    Table 2: The parameters of the lattices used in the thermodynamic extrapolation for nonzero values of B and the corresponding best fit parameters obtained with the help of Eq. (17). The data is visualised in Fig. 5.
obtained by fitting of the numerical data by the standard function, where A, m and η are the fitting parameters. The fitting parameters and parameters of the lattice at B = 0 are shown in Table 2. It turns out that in the insulator phase, B < B c , the data for η(L) can very well be fitted by the exponentially decaying (in the L → ∞ limit) function where L 0 and C are fitting parameters. The condensate in the infinite volume tends zero at B < B c , as expected. The fits are shown in Fig. 5 by the solid lines. The corresponding slopes are L 0 = 0.42(2) fm and L 0 = 0.90(8) fm for eB = 0 and eB = 1.07 GeV 2 , respectively.
At higher values of B, the condensate shows plateaux as L increases. We get the L → ∞ extrapolation for the condensate by averaging the data for η(L) at two largest values of the lattice size L (the horizontal dotted lines in Fig. 5). The extrapolated data -shown by the blue squares in Figure 3 -agrees quantitatively well with our small-volume analysis. The nonzero values of the extrapolated condensates are η = 0.0046(4) GeV 2 and η = 0.0095(1) GeV 2 for eB = 1.28 GeV 2 and eB = 2.14 GeV 2 , respectively. Thus, our numerical results support the theoretical prediction of Refs. [7,8] that the superconducting charged condensate of the ρ mesons forms spontaneously at strong magnetic field. Using the simulations of the quenched QCD vacuum, we determined the critical magnetic field (16) which is remarkably close to the theoretical prediction (1).
Finally, we would like to stress that theoretical calculations show that the condensate in the vacuum ground state should be an inhomogeneous function of spatial coordinates [7,8]. The ground state can be represented as a coherent static lattice of the topological (vortex-like) defects in the ρmeson condensates, the ρ vortices, which are directed along the magnetic field axis. Qualitatively, the ρ-vortex state is very similar to the Abrikosov vortex lattices observed in the type-II superconductors in a background of a strong magnetic field [18].
It turns out, however, that in the QCD vacuum the energy gap between the lowest vortex energy state (given by a triangular vortex lattice) and excited vortex lattice states is parametrically very small [19], implying that the spatial lattice order of the vortex state may be destroyed by quantum (or thermal) fluctuations. The latter fact indicates that the actual vortex structures in the superconducting phase may resemble a much less ordered but persistent "spaghetti state", where the correlation functions, given by Eq. (5) and/or Eq. (8), get additional suppression factors due to almost random vortex motion. The investigation of the detailed features of the superconducting ground state is currently underway [20].