# Generation of electric current by magnetic field at the boundary: quantum scale anomaly vs. semiclassical Meissner current out-side of the conformal limit

M. N. Chernodub, V. A. Goy, A. V. Molochkov

► To cite this version:

M. N. Chernodub, V. A. Goy, A. V. Molochkov. Generation of electric current by magnetic field at the boundary: quantum scale anomaly vs. semiclassical Meissner current out-side of the conformal limit. Advanced Physics Research, 2023, 2300058, 10.1002/apxr.202300058 . hal-04104539v2

HAL Id: hal-04104539

https://hal.science/hal-04104539v2

Submitted on 28 Aug 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Generation of electric current by magnetic field at the boundary: quantum scale anomaly vs. semiclassical Meissner current outside of the conformal limit

Maxim Chernodub* Vladimir Goy Alexander Molochkov

Prof. M. N. Chernodub

Institut Denis Poisson UMR 7013, Université de Tours, 37200 France

Email Address: maxim.chernodub@univ-tours.fr

Dr. V. A. Goy, Prof. A. V. Molochkov

Pacific Quantum Center, Far Eastern Federal University, 690950 Vladivostok, Russia

Keywords: scalar QED, scale anomaly, magnetic field, superconducting, boundaries, lattice simulations.

The scale (conformal) anomaly can generate an electric current near the boundary of a system in the presence of a static magnetic field. The magnitude of this magnetization current, produced at zero temperature and in the absence of matter, is proportional to a beta function associated with the renormalization of the electric charge. Using first-principle lattice simulations, we investigate how the breaking of the scale symmetry affects this “scale magnetic effect” near a Dirichlet boundary in scalar QED (Abelian Higgs model). We demonstrate the interplay of the generated current with vortex excitations both in symmetric (normal) and broken (superconducting) phases and compare the results with the anomalous current produced in the conformal, scale-invariant regime. Possible experimental signatures of the effect in Dirac semimetals are discussed.

## 1 Introduction

Quantum anomalies lead to distinctive transport phenomena, such as the chiral magnetic effect and its generalizations, which produce various currents of electric and axial charges in diverse physical environments, including systems at finite density and finite temperature, at strong electromagnetic fields, and in vortical backgrounds [1, 2].

The appearance of anomalous currents is directly related to breaking a classical symmetry by quantum fluctuations, with the most vivid examples given by the axial anomaly and the mixed axial-gravitational anomaly [3, 4]. These phenomena involve the axial (chiral) degrees of freedom that are expected to apply to various physical environments, from quark-gluon plasma created in ultrarelativistic heavy-ion collisions [5] to condensed matter environments [6].

The anomalous breaking of the classical conformal (scale) invariance [7, 8, 9] can also be responsible for anomalous transport effects [10]. For example, the scale electromagnetic effects generate an electric current and produce electric charge accumulation in electromagnetic backgrounds in curved spacetime [11]. In condensed matter context, the curved spacetime can be modeled as a temperature gradient that drives a system out of equilibrium [12]. The anomalous phenomena are then seen as the Nernst effect [13] in which the magnitude of the anomalous electric current is controlled by the beta function related to the renormalization of the electric charge. A purely gravitational conformal anomaly can reveal itself in various solid-state systems by affecting, for example, thermodynamic characteristics such as pressure gradients [14] and thermal transport [15, 16].

In a bounded classically conformal system, the scale anomaly can generate another type of anomalous electric current which reveal itself in the presence of the background magnetic fields [17, 18] (see also Ref. [19] for an earlier study). The anomalous current flows along the boundary with a magnitude proportional to the beta function of the electric charge, similar to the currents generated by the scale electromagnetic effects in bulk. In addition, the conformal anomaly leads to the electric charge accumulation, which can have an observable imprint in Dirac semimetals at charge neutrality point [20]. Anomaly-induced boundary phenomena also include a large variety of other effects associated with different symmetries or curved boundaries [21, 22, 23, 24, 25].

The appearance of electric boundary currents in the presence of the background magnetic field is in no way a new phenomenon. For example, the Meissner current can also be generated by the motion of the
Cooper pairs in a thin layer of a superconductor in response to the applied magnetic field [26]. As the magnetic field penetrates the superconductor, it induces a circulating current in the material that flows in such a way as to generate an opposing magnetic field. This current is carried by the Cooper pairs, which can flow freely without any resistance due to their superconducting nature.

There are several conceptual differences between anomalous and superconducting boundary currents.

The anomalous current is generated, as the name says for it, by the quantum anomaly in the vacuum.

On the other hand, the Meissner current has a semiclassical non-anomalous nature and is produced in matter (in the presence of a superconducting condensate). Our paper investigates quantitative differences between these currents on the unified footing within the same system.

We numerically study, from first principles, the generation of an electric boundary current in the lattice formulation of the (3+1) dimensional Abelian Higgs model (AHM) following the analysis of Ref. [27].

Our choice of this model is governed by its distinct features. First of all, the AHM possesses two phases: a phase with a spontaneously broken U(1) symmetry (the condensed or superconducting phase) and an unbroken (Coulomb) phase with a massive scalar field. In the latter phase, the theory is known as “scalar QED”. The AHM also possesses a region near a second-order phase transition where the mass gap vanishes, so the model effectively approaches a conformal limit. Therefore, all three interesting cases – represented by the conformal (scale-invariant), symmetric (massive fields with explicitly broken gauge symmetry), and broken (gapped phase with spontaneously broken gauge symmetry) regions – can be studied within the same model.

In the next section, we briefly discuss a qualitative physical picture behind the generation of the boundary current. We show that one can generally expect the existence of both conformal (quantum) and non-conformal (semiclassical) components to the boundary current. As it follows from their names, the conformal electric current is essential in a region of parameter space close to the conformal point. At the same time, the non-conformal contribution plays a dominant role in a part of the phase diagram with a mass gap in the spectrum of the matter field(s).

In Section 3, we describe the lattice model, observables, and the parameters of our lattice simulations.

Then, Section 4 is devoted to a detailed description of our numerical results. Finally, our conclusions are summarized in the last section.

# 2 Boundary current

## 2.1 Anomalous generation of the boundary current and conformal symmetry

The scale (conformal) anomaly produces an electric current near a flat boundary in a general class of spatially bounded quantum field theories with $U(1)$ gauge symmetry [19, 28]. The current takes the following general form [28, 17]:

$$j(\boldsymbol{x}) = -\nu(\boldsymbol{x})e\boldsymbol{n} \times \boldsymbol{B}(\boldsymbol{x}), \quad (1)$$

where $\boldsymbol{B}$ is the strength of the magnetic field, $e = |e|$ is the elementary electric charge, and $\boldsymbol{n}$ is the outward spatial vector normal to the boundary surface. The anomalous factor

$$\nu(\boldsymbol{x}) \equiv \nu(x_{\perp}) = \frac{2\beta(e)}{e^3} \frac{1}{x_{\perp}}, \quad (2)$$

is a function of the distance $x_{\perp}$ from the boundary to the point $\boldsymbol{x} = (x_{\perp}, \boldsymbol{x}_{\parallel})$ at which the current is generated. In our notations, $x_{\perp} = x$ and $\boldsymbol{x}_{\parallel} = (y, z)$.

The anomalous coefficient (2) depends on the beta function $\beta$ associated with the renormalization of the electric charge. Equation (1) implies that the induced electric current is tangential to the boundary and, simultaneously, normal to the axis of the magnetic field.

The electric current at the boundary (1) is proportional to the beta function $\beta$, which enters the proportionality coefficient (2). The beta function is associated with the renormalization of electric charge:

$$\beta(e) = \mu \frac{\partial e(\mu)}{\partial \mu}, \quad (3)$$
where $\mu$ denotes the energy scale at which the electric charge is measured: $e = e(\mu)$. Thus, the interactions break the conformal symmetry of the system via renormalization of the electric charge: the properties of particles interacting at different energies are no more related to each other by a simple rescaling transformation. The presence of the beta function (3) highlights the scale-anomalous nature of the produced electric current (1).

The conformal anomaly is also called the scale anomaly because the conformal and scale properties of reasonably defined field theories are inherently connected to each other [29]. The same type of anomaly is also often called the trace anomaly because in the presence of the background classical electromagnetic field $F^{\mu\nu}$, the trace of energy-momentum tensor $T^{\mu\nu}$ is nonzero: $T_{\mu}^{\mu} = \beta(e)/(2e)F^{\mu\nu}F_{\mu\nu}$. The last relation is helpful for first-principle calculations of the equation of state in non-perturbative gauge theories [30]. The loss of tracelessness originates from radiative corrections, as the strength-energy tensor should otherwise be traceless in classical conformal field theories in 3+1 dimensions. More details on conformal anomalies and transport can be found in Ref. [10].

Our paper concentrates on scalar QED (sQED), which is much easier to simulate numerically than the usual QED with massless or light fermions. The one-loop $\beta$ function of sQED with one bosonic species coupled to a single Abelian gauge field is four times smaller compared to the conventional QED:

$$ \beta_{\text{sQED}}^{1-\text{loop}} = \frac{e^3}{48\pi^2}, \qquad (4) $$

The anomalous coefficient (2) is therefore given by the following simple expression:

$$ \nu_{\text{cQED}}^{1-\text{loop}}(x_{\perp}) = \frac{1}{24\pi^2} \frac{1}{x_{\perp}}. \qquad (5) $$

The direction of the anomalous electric current (1) is tangential to the boundary and normal to the external magnetic field $B$. Therefore, it is natural to characterize the system by calculating the total current density per unit area of the boundary. In the conformal limit, the total electric current density induced at the boundary (1),

$$ J_{\text{tot}} = \int_{0}^{\infty} j_{\parallel}(x_{\perp}) dx_{\perp}, \qquad (6) $$

is a linear function of the background magnetic field:

$$ J_{\text{tot}} = \gamma e B, \qquad (7) $$

where the proportionality coefficient $\gamma$, given by the integration of Eq. (1) along the normal direction $x_{\perp}$, diverges both in infrared and ultraviolet limits:

$$ \gamma_{\text{conf}}^{\text{th}} = \frac{1}{24\pi^2} \ln \frac{\lambda_{\text{IR}}}{\lambda_{\text{UV}}}. \qquad (8) $$

In the conformal limit, the infrared cutoff $\lambda_{\text{IR}}$ in Eq. (8) should be of the order of the size of the system, while the ultraviolet cutoff $\lambda_{\text{UV}}$ should be determined by a typical shortest scale of the material at which the continuum description of the particle's motion is no more applicable (for example, in a crystal, the value of the cutoff $\lambda_{\text{UV}}$ is of the order of an interatomic distance). In our simulations, the lattice spacing naturally gives the ultraviolet scale, $\lambda_{\text{UV}} = a$. For a typical lattice size $L = 32$ (given in the units of lattice spacing $a$, with more details given below), the most extended physical size along any axis is $L/2 = 16$, implying that in our simulations, the logarithmic factor in Eq. (8) gives a modest multiplicative correction, $\ln 16 \simeq 2.8$. Despite the ultraviolet divergence both in infrared and ultraviolet limits, the total anomalous current (8) gives a reasonable value even in the realistic systems where the ultraviolet scale, given by an interatomic distance with a typical value of a few Ångströms (we take $\lambda_{\text{UV}} = 1 \AA \equiv 10^{-10} \text{m}$) is many orders of magnitude smaller than the typical crystal size (we take $\lambda_{\text{IR}} = 1 \text{cm} \equiv 10^{-2} \text{m}$). Then, the logarithmic factor is $\ln(\lambda_{\text{UV}}/\lambda_{\text{IR}}) \simeq 18.4 \simeq 6\pi$.
### 2.2 A single model for three regimes

The anomalous boundary current (1) has been found in the conformal limit where the model has no dimensionful parameters that could fix a scale. However, the conformal symmetry can be broken, and this breaking can proceed via two mechanisms: either perturbatively, via radiative corrections, or non-perturbatively, via the condensate. In our paper, we ask a natural question: can the boundary current be generated in a non-conformal phase that possesses a mass scale so that the conformal symmetry is explicitly broken? Furthermore, how does the breaking of the U(1) gauge symmetry by the condensate affects the generation of the current? In the condensate phase, both the conformal symmetry and the continuous gauge symmetry are broken simultaneously, which allows us to ask yet another question: how does the interplay between these two symmetries contribute to the generation of the anomalous current? In order to address all three regimes, we consider the Abelian Higgs model with the Lagrangian:

$$ \mathcal{L} = -\frac{1}{4} F_{\mu\nu} F^{\mu\nu} + (D_{\mu}\phi)^{*} D^{\mu}\phi - V(\phi), \qquad (9) $$

where $F_{\mu\nu} = \partial_{\mu}A_{\nu} - \partial_{\nu}A_{\mu}$ is the field strength tensor of the gauge field $A_{\mu}$, $D_{\mu} = \partial_{\mu} - ieA_{\mu}$ is the covariant derivative acting on the complex scalar field $\phi$, and

$$ V(\phi) = m^2|\phi|^2 + \lambda|\phi|^4, \qquad (10) $$

is the potential of the scalar field.

The Abelian Higgs model can host two types of conformally broken regimes corresponding to the symmetric and broken phases.<fnref n="1" /> In the symmetric phase, the Abelian U(1) symmetry is unbroken so that the physical content of the phase is the massless vector field $A_{\mu}$ and the massive scalar field $\phi$ with a vanishing condensate $\langle\phi\rangle = 0$. In the broken phase, the U(1) symmetry is spontaneously broken by the nonzero condensate $\langle\phi\rangle \neq 0$, and both the vector field $A_{\mu}$ and a scalar excitation over the condensate, $\delta\phi \equiv \phi - \langle\phi\rangle$, acquire masses.

<figure><img src="image_2.png" /><figcaption>Figure 1: The illustration of the mechanisms of the electric current generation near the boundary in (a) the conformal regime and (b) the symmetry-broken phase. (a) Vacuum fluctuations produce particle-anti-particle pairs that follow skipping orbits without mutual annihilation, thus creating a uniform electric current along the reflective boundary. (b) Depletion of the scalar condensate acts as a well that attracts vortices to the boundary, thus enhancing the circular currents of the vortices and generating a coherent magnetization (Meissner) current along the boundary.</figcaption></figure>

### 2.3 Electric current generation

#### 2.3.1 Conformal current in the conformal regime

In the conformal limit of the zero-temperature model, the quantum fluctuations produce the boundary current [19, 28, 18]. The mechanism of the current generation is visualized in Fig. 1(a). Initially, vacuum fluctuations create particle-anti-particle pairs which, placed far from a boundary, annihilate within

<footnote n="1">¹In a lattice version of the Abelian Higgs model (9) the gauge field may be a compact field, which means the presence of the third, electrically confining phase. Our paper considers a weakly coupled region ($e \ll 1$) where the confining phase is not realized.</footnote>
the time interval stipulated by the uncertainty principle. The presence of the magnetic field background does not change the fate of these virtual particles, which will, in this case, follow close circular orbits. However, a virtual particle-antiparticle pair created near the reflective boundary does not close its shared orbit, provided the magnetic field is directed tangentially to the boundary. As a result, the virtual particles follow skipping orbits near the wall without mutual annihilation, thus creating an electric current along the boundary (notice that particles with opposite charges move in opposite directions along the wall, thus doubling the total current). The current profile near the boundary has an infinite thickness in the normal direction to the wall because – in the absence of any scales – the current decays polynomially (not exponentially) at large distances from the boundary.

The produced current has an entirely quantum origin since it originates from the vacuum at zero temperature without matter. The vanishing mass of particles facilitates the creation of their pairs from the vacuum, thus enhancing the boundary current and enforcing its long-ranged nature. The existence of the boundary current in the conformal limit has been confirmed numerically from the first principles in Ref. [27] to which we refer an interested reader. Below, we complete the results in the conformal point by calculations at the symmetry broken phase where the current generation is governed by an entirely different mechanism and symmetric phase, where the current generation has similar features to the conformal symmetry point.

### 2.3.2 Meissner current in the broken phase

Outside the conformal region, in the broken phase with a nonzero scalar condensate, the induced boundary current is generated by Abrikosov vortices via a simple semi-classical mechanism visualized in Fig. 1(b). An external magnetic flux parallel to the boundary creates several Abrikosov vortices equal to the number of elementary fluxes in the total magnetic flux. At the center of each vortex, the scalar condensate vanishes, and the magnetic field takes its maximum value. The radial profile of the magnetic field is shaped by a circular electric current that circumvents the vortex core in the plane normal to the core. The current is zero at the center of the vortex, takes its maximal value at the distance of the order of the penetration depth $\lambda$, and decays exponentially at longer distances.

The scalar condensate should vanish at the boundary due to the Dirichlet condition on the scalar field. The condensate is then restored to its vacuum expectation value at a certain distance from the boundary, which is of the order of the coherence length $\xi$ (the correlation length of the scalar field). The Dirichlet boundary should then attract the Abrikosov vortices because the energy of the vortex is a monotonically rising function of the value of the scalar condensate, while at the boundary, the condensate is smaller than in bulk. Therefore the vortices tend to be “pinned” to the Dirichlet boundary so that the two-dimensional density of the vortices is expected to be increased in the vicinity of the boundary. As a result, the circular electric currents of individual vortices cancel each other in the normal direction to the boundary while simultaneously generating a net electric current along the boundary.<fnref n="2" /> Thus, we get a boundary current that should decay exponentially along the normal to the boundary. The boundary current has a non-anomalous, semi-classical origin in the broken phase. In this semi-classical scenario, the thickness of the boundary current is set by the penetration depth of the superconductor.

This mechanism should be contrasted to the flux barrier at the surface of realistic superconductors with a vacuum (insulator) where the vortices experience tunneling due to the energy barrier at the boundary [26]. Likewise, to our example above, the background magnetic field generates the Meissner current along the boundary of the superconducting samples, although the boundary conditions in these cases are different. The current is generated despite the difference in the boundary conditions (the Dirichlet boundary condition in our case vs. the Neumann condition on the charged field in the Ginzburg-Landau formulation of a superconducting condensate [26]).

<footnote n="2">²One may expect that the result of our considerations is valid in type-I and type-II superconductivity regimes where vortices, respectively, attract and repel each other.</footnote>
### 2.3.3 Mixed regime in the symmetric phase

In the symmetric phase, the classical Abrikosov vortex solutions do not exist, as the penetration length is infinite (the photon is a massless particle), and no topological mechanism exists that could stabilize the classical solutions. However, the described mechanism will nevertheless work since the external magnetic field induces light (quantum) vortices which may be thought of as clumps of magnetic field shaped by the vacuum currents of a fluctuating scalar field. Such quantum vortices are often observed in the numerical simulations of scalar models [31].

The correlation length of the scalar field sets the size of the effective scalar core of the emergent vortices. Despite their non-classical nature, the fluctuating vortex lines which would tend to concentrate at the vicinity of the boundary where the Dirichlet condition naturally suppresses the value of the scalar field, thus making it close to the vanishing value of the scalar field in the center of the vortex. The elevated density of vortex lines will induce the (again, non-anomalous) boundary current, which is normal to the magnetic field. The thickness of the boundary current is set by the correlation length of the scalar field, which is the only length scale in the theory.

## 3 Lattice model

### 3.1 Action of the model and dynamical fields

The Euclidean lattice action of the Abelian Higgs model (9) is given by the following formula:

$$S = \beta_{\text{latt}} \sum_{x} \sum_{\mu<\nu=1}^{4} (1 - \cos \theta_{x,\mu\nu}) + \sum_{x} \sum_{\mu=1}^{4} \left| \phi_x - e^{i(\theta_{x\mu} + \theta_{x\mu}^B)} \phi_{x+\hat{\mu}} \right|^2 + \sum_{x} \left( -\kappa |\phi_x|^2 + \lambda |\phi_x|^4 \right), \quad (11)$$

where $x \equiv (x_1, x_2, x_3, x_4)$ is the Euclidean space-time coordinate, $\phi_x$ is the complex scalar field and $\theta_{x\mu}$ is the dynamical (quantum) vector gauge field. The vector gauge field $\theta_{\mu}^B$ is the background (classical) magnetic field.

The model is characterized by the lattice spacing $a$, which corresponds to the physical length of an elementary lattice link. In a naive continuum limit, $a \to 0$, the dimensionless gauge $\theta_l$ and scalar $\phi$ fields are related to their continuum counterparts (both of the dimension of mass) via the relations $A_\mu(x) = \theta_{x,\mu}/a$ and $\phi(x) = \phi_x/a$, respectively. The physical value of the lattice spacing $a$ is usually determined by matching dimensionless lattice results to known dimensionful quantities (for example, to the mass of a particular physical excitation).

In the continuum limit ($a \to 0$), the lattice gauge coupling $\beta_{\text{latt}}$ in Eq. (11) is related to the bare electric charge $e$ as follows: $\beta_{\text{latt}} = 1/e^2$ (here, we introduce the subscript “latt” in the lattice coupling “$\beta_{\text{latt}}$” in order to discriminate it from the beta function $\beta$). The bare couplings $\kappa$ and $\lambda$ are associated, respectively, with quadratic and quartic terms of the scalar field $\phi$ in the interaction potential (10) corresponding to the last line in Lagrangian (11).

Action (11) is invariant under the Abelian gauge transformations for the gauge field, $\theta_{x\mu} \to \theta_{x\mu} + \omega_x - \omega_{x+\hat{\mu}}$, and the scalar field, $\phi_x \to e^{i\omega_x} \phi_x$, where $\omega_x$ is an arbitrary real-valued scalar function defined at the sites of the lattice. The background gauge field $\theta_{x\mu}^B$ is insensitive to the gauge transformations.

In the lattice model (11), the gauge field $\theta_l$ is a compact field variable because the action is invariant under the discrete shifts $\theta_l \to \theta_l + 2\pi n_l$, where $n_l \in \mathbb{Z}$ is an arbitrary integer. The compactness of the gauge field automatically implies the existence of Abelian monopoles in the vacuum of the theory. The monopoles condense in the strong coupling region with $e \gtrsim 1$ (small $\beta_{\text{latt}}$), and the vacuum of the model becomes confining. As we are not interested in the confining effects, we keep the lattice gauge coupling sufficiently large, thus ensuring that the monopoles rarely appear in our numerical calculations. In other words, the compactness of the lattice gauge field $\theta_l$ does not influence our results.

The numerical simulations of the present article are similar to the ones of Ref. [27], albeit now we work in a different region of parameters. We use the symmetric lattice $L^4$ and the elongated lattice $L^3 \times L_x$
with $L_x \ge L$. Both lattices correspond to a zero-temperature theory. We use $L = 32$ and $L_x = 32, 48$ with periodic boundary conditions imposed at each direction.

In this article we simulated the model (11) at the fixed parameters $\beta_{\text{latt}} = 4$ and $\lambda = 10$. We generated field configurations using a Hybrid Monte Carlo algorithm [32, 33] and performed simulations on Nvidia GPU cards. To achieve acceptable statistics, we used about $10^6 \dots 10^7$ trajectories per each value of the background magnetic field. To accelerate calculations and reduce write operations, we accumulated mean values per every 100 trajectories only, as even in this case, our simulations generated about 1.5 TB of data. We also used binning for correct error estimations for our observables.

## 3.2 Background: magnetic field and boundary wall

We insert the reflective Dirichlet boundary for the scalar field at the middle of the lattice at $x = L_x/2$:

$$ \phi_x \Big|_{x_2=L_x/2} = 0. \qquad (12) $$

The static uniform magnetic background field is introduced along the $z \equiv x_3$ axis:

$$ \theta_{x,12}^{B} = \frac{2\pi k}{L^2}, \qquad (13) $$

where other components of the field-strength tensor are taken to be zero. We use the following parametrization of the background gauge field [34]:

$$ \theta_{x,2}^{B} = \frac{2\pi k}{L_x L_y} x, \qquad \theta_{x,1}^{B}|_{x=L_x-1} = -\frac{2\pi k}{L_y} y, \qquad (14) $$

with $\theta_{x,3}^B = \theta_{x,4}^B = 0$ and $k = 0, 1, \dots, L_x L_y/2$.

The integer number $k$ in Eq. (14) determines the strength (13) of the magnetic field. The quantization of the electromagnetic gauge field (14) is an essential infrared property needed to satisfy the periodic boundary conditions imposed on the gauge field:

$$ e^{i\theta_{x+L,\mu}^{B}} = e^{i\theta_{x,\mu}^{B}}. \qquad (15) $$

On the contrary, the bound on $k$ from above is an ultraviolet feature related to the presence of the shortest physical length $a$ as the value $k$ corresponds to the number of elementary fluxes introduced on the lattice in the $xy$ plane. The maximal value of $k$ in Eq. (14) corresponds to the lattice half-filled with Abrikosov vortices (so that the physical vortex density is given by the ultraviolet lattice cutoff $\rho_{\text{vort}} \sim a^{-2}$). We will work at weak magnetic fields $\theta^B$, far from the artificially large values of the integer number $k$. The lattice magnetic field (13) is related to the continuum field $B$ via the relation $\theta_{x,12}^B = eBa^2$. All dimensionful quantities in the rest of the paper are presented in lattice units, with the lattice spacing $a$ set to unity.

## 3.3 Observables: vortices and currents

The lattice density of the Abrikosov vortices is given by the integer number [35]

$$ \rho_{x,\mu\nu} = \frac{1}{2\pi} (l_{x,\mu} + l_{x+\hat{\mu},\nu} - l_{x+\hat{\nu},\mu} - l_{x,\nu}), \qquad (16) $$

defined at every plaquette $P = P_{x,\mu\nu}$ of the lattice. The vortex density (16) is the lattice curl of the gauge-invariant link variable

$$ l_{x,\mu} = \arg \left[ \phi_x^* e^{i(\theta_{x\mu} + \theta_{x\mu}^B)} \phi_{x+\hat{\mu}} \right]. \qquad (17) $$

The variable $\rho_P$ is equal to zero, $\rho_P = 0$, if no vortex pierces plaquette $P$, and it gives $\rho_P = +1$ ($\rho_P = -1$) if a vortex (an antivortex) pierces the plaquette $P$.

The electric current is given by the variation of the action (11) with respect to the gauge field $\theta_{x\mu}$:

$$ j_{x\mu} = -2 \left\langle \text{Im} \left[ \phi_x^* e^{i(\theta_{x\mu} + \theta_{x\mu}^B)} \phi_{x+\hat{\mu}} \right] \right\rangle. \qquad (18) $$
# 4 Numerical results

This section reports the results of the first-principle numerical simulations on the boundary current in both the symmetric and broken phases of the Abelian Higgs model defined in Section 3 and compares them with our theoretical expectations described earlier in Section 2.3.

## 4.1 Phase transition

First, we examine the symmetric phase far from the conformal regime. To eliminate the presence of Abelian monopoles and make the gauge field behave similarly to the continuum theory (9), we set the large gauge coupling to a fixed value of $\beta = 4$. In addition, we choose a small value for the quartic constant, namely $\lambda = 0.01$, which guarantees the absence of a second-order phase transition at any value of the two-point coupling $\kappa$. We work at the lattices $32^4$ and $32^3 \times 48$. In the middle of the lattice, at $y = 16$, we put the Dirichlet boundary (12) to simulate a reflective wall.

<figure><img src="image_3.png" /><figcaption>Figure 2: The (unrenormalized) quadratic fluctuations of the scalar field as the function of the lattice hopping parameter $\kappa$ in the vicinity of the first-order phase transition at fixed quartic coupling $\lambda = 0.01$ and the gauge coupling $\beta = 4$ at $32^4$ lattice at a vanishing magnetic field, $B = 0$. The discontinuity in the quadratic fluctuations marks the first-order phase transition.</figcaption></figure>

The behavior of the quadratic fluctuations of the scalar field $\langle|\phi|^2\rangle$ as the function of the hopping coupling $\kappa$, Fig. 2, demonstrates the existence of a first-order phase transition at $\kappa_c \simeq 0.121$ with the symmetric (broken) phase at smaller $\kappa < \kappa_c$ (larger $\kappa > \kappa_v$) values of the coupling, separated by a discontinuity in $\langle|\phi|^2\rangle$.

Therefore in the symmetric and broken phases, we choose the following set of lattice parameters, respectively:

$$ \beta = 4, \quad \lambda = 0.01, \quad \kappa = 0.12 \quad [\text{symmetric}], \qquad (19) $$

$$ \beta = 4, \quad \lambda = 0.01, \quad \kappa = 0.2 \quad [\text{broken}]. \qquad (20) $$

The presence of the background magnetic field can also influence the position of the phase transition point by restoring the symmetric (normal) phase from the broken (superconducting) phase. At the chosen point of parameters (20), this transition happens at a relatively large value of the magnetic field,

$$ eB_c \simeq 0.87, \qquad (21) $$

which is seen as a discontinuity in the expectation value $\langle|\phi|^2\rangle$, Fig. 3. Therefore, we consider the effect of the current generation with the fields lower than the critical value, $B < B_c$. On the contrary, the model expectedly remains in the symmetric phase at all considered values of the background magnetic field at our set of parameters (19). Below, we study these phases in more detail.
<figure><img src="image_4.png" /><figcaption>Figure 3: The quadratic fluctuations of the scalar field as the function of the magnetic field strength $eB$ at fixed couplings (20) corresponding to the broken phase at $32^3 \times 48$ lattice. The thin red vertical line approximately marks the position of the critical magnetic field, which corresponds to the first-order phase transition from the broken phase to the symmetric phase.</figcaption></figure>

### 4.2 Symmetric phase

First, we consider the symmetric phase represented by the set of parameters (19). In Fig. 4, we show the expectation value of the scalar field fluctuations $\langle|\phi|^2\rangle$ at a zero magnetic field. As expected, the fluctuations vanish at the boundary and recover quickly to the bulk value in a few lattice spacings from the boundary wall. This picture is practically independent of the strength of the background magnetic field.

<figure><img src="image_5.png" /><figcaption>Figure 4: Symmetric phase: The mean value of the scalar field. The light magenta line marks the position of the boundary. Simulations are performed at the parameter set (19).</figcaption></figure>

Figure 5 shows the expectation value of the electric current (18), revealing the presence of the electric current perpendicular to the external magnetic field. We calculated by averaging the local current over a moderate number of configurations to highlight the magnitude of quantum fluctuations. When the magnetic field is zero, the current appears as a small random vector (top panel in Fig. 5), with a vanishing expectation value if averaged over a sufficiently large number of configurations. However, as a nonzero magnetic field is applied tangentially to the boundary, the presence of the boundary electric currents becomes evident. The currents flow opposite directions at opposite sides of the boundary (two lower panels in Fig. 5). Increasing the strength of the magnetic field reduces the relative randomness of the generated boundary current while amplifying its intensity.

In Fig. 6, we show the local density of the boundary current $j_y$ as the function of the distance from the
<figure><img src="image_6.png" /><figcaption>Figure 5: Symmetric phase (19): The (normalized) electric current in the *xy* plane in the vicinity of the wall (positioned at *x* = 16 and directed along the *y* axis) at three values of increasing magnetic field. The current is calculated with a moderate number of lattice field configurations which allow us to highlight the magnitude of its fluctuations. The color of the arrows corresponds to the sign of the *y*-component of the current.</figcaption></figure>

wall $x_{\perp}$. We consider only positive values of $x_{\perp}$ since the current is perfectly anti-symmetric with respect to the inversion $x_{\perp} \rightarrow -x_{\perp}$. It turns out that the dependence of the boundary current on the distance from the wall $x_{\perp}$ can be excellently (with $\chi^2/d.o.f. \lesssim 1$) described by the following fitting function:

$$j_y^{\text{fit}}(x_{\perp}) = \frac{2MJ_{\text{tot}}}{\pi \cosh(Mx_{\perp})} \text{sign } x_{\perp}, \qquad (22)$$

with two fitting parameters: the mass parameter *M* and the total current *J*<sub>tot</sub>.

The mass parameter *M* controls the thickness of the current in the normal direction, $\lambda_{\perp} = M^{-1}$, an essential feature of the current profile in the broken phase. For example, an attempt to describe the data with a conformal-like fit $j_y^{\text{fit}} \sim 1/|x_{\perp}|$ with either infrared or ultraviolet cutoffs (or even with both) gives us unreliable results with poor qualities of the fits. The last fact is a natural consequence of the absence of conformal symmetry in a deep symmetric phase (19) which calls us to introduce the mass parameter *M*, which breaks the conformal symmetry. The fitting parameter *J*<sub>tot</sub> corresponds to the total current, obtained by integration of the current (22) along the whole transverse direction (6). The total current is a finite quantity, both in infrared and ultraviolet limits.

The localization of the electric current in the vicinity of the boundary is determined by the value of the mass of the fit (22), which is shown in Fig. 7(a) as the function of the magnetic field strength $eB$. The localization mass *M* increases with the background magnetic field, which is not an unusual effect since masses of elementary electrically charged scalar excitations should be rising, in a relativistic free field theory, as $M(eB) = \sqrt{m^2 + |eB|}$. The last property suggests us the form of the fit of the numerical data:

$$M^{\text{fit}}(eB) = \sqrt{M_0^2 + g|eB|}, \qquad (23)$$
<figure><img src="image_7.png" /><figcaption>Figure 6: Symmetric phase (19): The local electric current $j_y$ at the function of the distance $x_⊥$ from the boundary at various values of the magnetic field $eB$, Eq. (13). The dashed lines correspond to the best fits by function (22). The inset shows a zoom-in on the region close to the wall.</figcaption></figure>

where the factor $g$ takes into account the fact that the particles are interacting, and the increase of magnetic field drives the system deeper into the symmetric phase. The best fit, shown in Fig. 7(a) by the dashed line, gives us the parameters $M_0 = 0.53(1)$ and $g = 2.2(1)$. Thus, the stronger the magnetic field, the thinner the current density distribution at the wall.

<figure><img src="image_8.png" /><figcaption>Figure 7: Symmetric phase (19): (a) The screening mass factor $M$ of the current density near the boundary (22) vs. the magnetic field strength. The dashed line represents the best fit (23). (b) The total electric current generated at the boundary vs. the external magnetic field. The line represents the best fit by the linear function (24).</figcaption></figure>

In Fig. 7(b), we show the total boundary current as the function of the magnetic field. In a qualitative agreement with the general linear-response arguments (1), we find that the total current, similarly to the local current density, is a linear function of the background magnetic field:

$$J_{\text{tot}}^{\text{fit}} = \gamma eB. \qquad (24)$$

Function (24) matches our numerical data very well ($\chi^2/d.o.f. \sim 1$), as it seen in Fig. 7(b).

The best fit gives us for the proportionality coefficient $\gamma = 0.0074(1)$ which matches well the order of the theoretical coefficient (8) as, theoretically, we expect to obtain $\gamma^{\text{th}} = 1/(24\pi^2) \approx 0.0042$. Formally, if it is the conformal anomaly that was responsible for the generation of the boundary current, then the relation between the infrared and ultraviolet cutoffs in Eq. (8) would be quite reasonable as well: $\lambda_{\text{IR}} \approx 6\lambda_{\text{UV}}$.
<figure><img src="image_9.png" /><figcaption>Figure 8: Symmetric phase (19): The (normalized) local vortex density as the function of the distance to the wall at the $32^3 \times 48$ lattice with an extended spatial direction $L_x = 48$.</figcaption></figure>

The conformal anomaly cannot contribute much to the boundary current in a deep symmetric phase due to a large mass gap. In order to check the alternative, vortex-based semi-classical mechanism discussed in the introduction, in Fig. 8, we plot the local vortex density (16) as the function of the distance to the boundary. We normalize the local vortex density to the mean vortex density on the whole lattice. The numerical data is consistent with the vortex mechanism of the current generation represented in Fig. 1(b): the density of the vortices is elevated close to the boundary, implying that the circular electric currents that circumvent the vortex cores are generating the boundary current. Figure 8 also demonstrates the existence of a repulsion force between the vortices as the boundary layer of the vortices repels the vortices from the bulk, causing a noticeable drop in the vortex density deeper in bulk.

In other words, vortices are attracted to the wall causing a small, of the order of one percent, excess of the vortex density at the wall. In the symmetric phase, the quantum vortices repel each other similarly to the type-II superconductor. Therefore, some of the vortices are concentrated precisely at the position of the wall, making the vortex density elevated exactly at the boundary. Those boundary-pinned vertices are repelling the other vortices that emerge from the bulk, thus causing the drop in the vortex density near, but not strictly at, the position of the wall. This mechanism, highlighted qualitatively in Fig. 1(b), is responsible for the physical picture observed in Fig. 8.

We finish the study of the current in the symmetric phase by addressing the question whether the near-boundary behaviour of the current, determined numerically by the fitting function (22), related to the conformal behaviour, predicted analytically by Eqs. (1) and (2). The naive zero mass limit, $M \to 0$, does not recover the conformal $1/x_\perp$ behavior from the functional form $1/\cosh(Mx_\perp)$. Our attempts to fit the data by the function (22) with $\sinh(Mx_\perp)$ instead of $\cosh(Mx_\perp)$ lead to inconsistent results with a poor quality of fit. In detail, we found d.o.f $\simeq 10...20$ for the sinh-like behaviour in comparison with good d.o.f $\simeq 0.5...1.5$ for the cosh-fit function (22). These results imply that the point of the coupling space considered in this article (19) lies deeply in the symmetric phase which does not allow for this simple continuation to the conformal point (1).

## 4.3 Broken phase

Now we repeat all our calculations in the broken phase given by the initial set of parameters (20) at a larger lattice volume $32^3 \times 48$, where the longer lattice size extends one of the directions of the boundary, $L_x = 48$. This theory corresponds to zero temperature.

First, we visualize in Fig. 9 the electric current generated near the wall at increasing values of the background magnetic field. While the produced electric current has the same direction as in the symmetric phase, Fig. 5, one immediately notices the difference between the broken phase and the symmetric phase in the form of the current density profiles. In the symmetric phase, the current density reaches its maximum at the wall, but in the broken phase, the maximum takes place at a certain distance from the wall.
<figure><img src="image_10.png" /><figcaption>Figure 9: Broken phase (20): The (normalized) electric currents generated close to the wall as in Fig. 5 at various values of the background magnetic field B. Notice that the lowest panel corresponds to the symmetry restored phase at B > B_c with the critical magnetic field B_c given in Eq. (21).</figcaption></figure>

Figure 9 also shows that the current maximum occurs at longer distances from the boundary as the magnetic strength increases. However, at the highest magnetic field shown in the lowest panel of the same figure, the electric current changes its profile again by shifting its maximum to the boundary in close similarity with the symmetric phase. This sudden change in the behavior is not unexpected since the system experiences the transition to the symmetric phase as the magnetic strength achieves the critical value (21), and at B > B_c, the gauge symmetry is restored. Therefore, the lowest panel of Fig. 9 mimics the results in the symmetric phase qualitatively, Fig. 5. Notice that these figures show the normalized electric current, thus highlighting the shape of the current rather than its magnitude.

<figure><img src="image_11.png" /><figcaption>Figure 10: Broken phase (20): The electric density current $j_y$ at the function of the distance $x_⊥$ from the boundary at various values of the magnetic field $eB$. The dashed lines correspond to the best fits by function (25).</figcaption></figure>

The electric density current $j_y$ as a function of the distance from the boundary is shown in Fig. 10 at various magnetic field values below the critical value (21). This figure shows that the current density profile, which takes its maximum at a certain distance from the wall, can be described very well by the function

$$j_y^{\text{fit}}(x_{\perp}) = \frac{M^{1+\nu}|x|^{\nu}J_{\text{tot}}}{\Gamma(1+\nu)\exp(Mx_{\perp})}\text{sign } x_{\perp}, \qquad (25)$$
where the total current $J_{\text{tot}}$, the effective mass *M* and the power *ν* are the fitting parameters. The pref-actors of Eq. (25) are chosen in such a way that the quantity $J_{\text{tot}}$ corresponds to the total current ac-cording to the standard normalization (6).

Figure 10 shows that the generated electric current increases in magnitude and takes its maximum fur-ther from the wall as the magnetic field increases. This picture works while the strength of the mag-netic field still resides below the critical value (21). The total current, shown in Fig. 11(a), confirms this observation by showing a linear growth and the saturation of the current below the critical magnetic field (21), featuring, at the same time, the sudden drop of the current as the magnetic field exceeds the critical value above which the symmetry gets restored.

<figure><img src="image_12.png" /><figcaption>Figure 11: Broken phase (20): (a) The total current density $J_{\text{tot}}$, (b) effective mass parameter *M* and (c) the exponent *ν* as functions of the magnetic field *eB*. The inset in the upper plot shows the zoom-in on the symmetric phase above the critical magnetic field (21).</figcaption></figure>

Remarkably, Fig. 11 shows that in the symmetry-restored phase, the electric current is three (!) orders of magnitude smaller than in the symmetry-broken phase. In other words, the semi-classical mechanism of the current generation based on the semi-classical vortex solutions drastically exceeds the mechanism based on the quantum vortices, which emerge in the symmetric phase as unstable excitations. Both these mechanisms, illustrated qualitatively in Fig. 1(b), differ from the anomalous current generation catalyzed by the conformal anomaly [28].

Notice that the magnitude of the electric current produced in the symmetric phase, which is, in turn, generated by the strong magnetic field – corresponding to high values of the magnetic field in Fig. 11(a) – and the magnitude of the electric current in the genuine symmetric phase show in Fig. 7(b), are the same. This fact highlights some universality of our results, distinguishing clearly the mechanism of the current generation in the symmetric phase from the one in the broken phase.
Moreover, the total current that emerges in the conformal region – according to the earlier results of Ref. [27] – appears to be of the same order as the current generated in the symmetric phase. Therefore, while the mechanisms of the current generation in the conformal region and the symmetry restored phase are different – given by Fig. 1(a) and Fig. 1(b), respectively – they produce the current of the same magnitude which is, in turn, much smaller than the Meissner current generated in the symmetry broken phase. The effective mass obtained from the current distribution with the help of the fit (25) raises with the increase of magnetic field in agreement with a general tendency (23) of the mass of the charged scalar particle in the magnetic field background. The mass, shown in Fig. 11(b), experiences a sudden drop at the critical magnetic field (21) due to the disappearance of the condensate and then again raises with $eB$ following the same qualitative behaviour (23).

The power factor $\nu$ of the current density profile (25) increases with the increase of magnetic field as shown in Fig. 11(c), raising from a modest $\nu \simeq 3$ to a very high value $\nu \simeq 15$ close to the transition point (21).

<figure><img src="image_13.png" /><figcaption>Figure 12: Broken phase (20): (a) the scalar field fluctuations $\langle |\phi|^2 \rangle$ and (b) the normalized vortex density at a few values of the background magnetic field.</figcaption></figure>

Finally, in Figs. 12(a) and (b), we show the behavior of the scalar field and the normalized vortex density, respectively, as functions of the normal distance to the wall. The increasing magnetic field widens the hollow in the condensate near the boundary and creates a double layer of vortices near the wall. This nontrivial structure appears due to the specific geometry of our problem and the mutual repulsion of the vortices interacting through the wall. In addition, the vortices create a strong electric Meissner current due to the mechanism depicted in Fig. 1(b).

# 5 Summary and Conclusions

Our paper is concentrated on first-principle numerical results for electric current induced by a background magnetic field close to a boundary in the vacuum of the Abelian Higgs model. This model represents a scalar counterpart of quantum electrodynamics with charged scalar fields instead of electrically charged spinors (therefore, it is also called “scalar QED”). Our investigation is carried out at zero temperature. In the conformal limit of the zero-temperature model, the boundary current arises from quantum fluctuations [19, 28, 36] due to a simple mechanism [28]. Initially, vacuum fluctuations generate particle-anti-particle pairs which move along closed circular paths in a magnetic field background and then annihilate. However, their mutual orbit remains open when a virtual particle-anti-particle pair is formed near a reflective boundary with the magnetic field oriented tangentially to the boundary. Consequently, the virtual particles follow skipping orbits without fast annihilation, thus establishing the generation of electric current along the boundary. The masslessness of particles facilitates the creation of particle pairs from the vacuum, amplifying the boundary current and enforcing its long-range nature. The strength of the
current is proportional to the beta function responsible for the renormalization of the electric charge. Outside the conformal limit, the scalar electrodynamics can reside either in the symmetry-restored phase featuring a massive scalar field or in the regime with spontaneous symmetry breaking in which the charged scalar field is condensed. We found that in the symmetry-unbroken regime, the current-generation effect gets suppressed in its amplitude far from the boundary because of the short-range nature of fluctuations of the massive virtual particles in agreement with the analytical expectations [24]. However, the total current that emerges in the conformal region appears to be of the same order as the current generated in the symmetric phase.

Our numerical calculations show that the quantum boundary current gets drastically enhanced in the symmetry-broken phase due to the condensation of charged particles. We consider the reflective boundary condition of the Dirichlet type, which enforces the vanishing of the scalar field at the boundary and is qualitatively different from the superconducting (Neumann-type) boundary conditions. We consider a superconducting regime with a transparent semi-classical vortex-based mechanism to generate the non-conformal regime's boundary current. We find that the vortices are attracted to the boundary, and their circular currents at up coherently, thus generating the global magnetization current along the boundary of the material. The mechanism is similar to the standard Meissner effect in which the vortices are forced to flow around the surface, creating a circulating current that shields the superconductor from the external magnetic field [26]. As the superconducting condensate is forced to vanish at the boundary, the vortices near the boundary are more energetically favorable than in the bulk of the material.

The boundary current in the conformal regime possesses a much smaller amplitude than the current in the symmetry-broken phase since the former effect is a purely quantum phenomenon while the latter has a semi-classical nature. Despite the fact that the origin and magnitude of these effects differ, they lead to the same phenomenon, the dissipationless magnetization current along the boundary.

Finally, we would like to discuss experimentally observable consequences of the induced magnetization current due to the conformal anomaly in Dirac semimetals. In one-component massless QED, the magnetization current (1) leads to the enhanced edge conductivity which can be estimated as

$$ \delta\sigma_{xy} = \frac{2\beta}{e\hbar} \ln \frac{L_{IR}}{L_{UV}}, \qquad (26) $$

where $L_{IR}$ and $L_{UV}$ are infrared and ultraviolet distance scales, respectively. Outside the conformal limit, which is typical for most of the Dirac semimetals, the former length scale is inversely proportional to the mass which breaks the scale invariance. The ultraviolet length scale in Eq. (26) can be taken to be the inter-atomic distance $L_{UV} \simeq a$, which can be taken a few Angstroms. Taking, for the sake of estimation, the mass scale to be equal about 10 K (in temperature units) and $a$ of the order of a few Angstrom, one can estimate the logarithmic factor in Eq. (26) as $\ln(L_{IR}/L_{UV}) \simeq 4\pi$, which gives us

$$ \delta\sigma_{xy} = \frac{4c}{3\varepsilon v_F} \sigma^{\text{Hall}}, \qquad (27) $$

where $\sigma^{\text{Hall}} = e^2/h$ is the Hall conductivity. Due to the nature of logarithm function, a order-of-magnitude change in $L_{IR}$ does not affect the prefactor of Eq. (27) significantly. Notice that the QED beta function is four times bigger than the one in the scalar QED (3).

In Eq. (27), which takes peculiarities of the running of electric charge in Dirac semimetals [37], $\varepsilon \simeq 10$ is a dielectric constant and $v_F \simeq c/300$ is the Fermi velocity. The slowness of electrons enhances the beta function and makes significant contribution to the magnetization current resulting in the large value of conductivity associated with the magnetization current:

$$ \delta\sigma_{xy} \simeq 40\sigma^{\text{Hall}}. \qquad (28) $$

Notice that this estimation (28) does not refer to the Hall effect but to the magnetization current, which is confined to the boundaries to material and, therefore, does not generate any net transverse current.

The persistent boundary current creates accompanying magnetic field which can be probed, for example, using scanning probe based on a nitrogen-vacancy (NV) defect in a diamond as an atomic-size quantum
magnetometer [38]. The NV-based measurement is very sensitive to the magnetic field through the Zeeman effect, which allowed, for example, to probe the hydrodynamic transport of electrons in bulk crystals $WTe_2$ [39] at low temperatures. Another, perhaps, less complicated option is to measure the enhanced local conductivity at the edge of the crystal using microwave impedance microscopy which has been employed to observe the quantum spin Hall edges in monolayer $WTe_2$ in Ref. [40].

**Acknowledgements**

The work of VAG and AVM was supported by Grant No. 0657-2020-0015 of the Ministry of Science and Higher Education of Russia. The numerical simulations were performed at the computing cluster Vostok-1 of Far Eastern Federal University.

**References**

[1] D. E. Kharzeev, K. Landsteiner, A. Schmitt, H.-U. Yee, *Lect. Notes Phys.* **2013**, *871* 1.

[2] D. E. Kharzeev, *Prog. Part. Nucl. Phys.* **2014**, *75* 133.

[3] K. Landsteiner, E. Megias, F. Pena-Benitez, *Phys. Rev. Lett.* **2011**, *107* 021601.

[4] K. Landsteiner, E. Megias, F. Pena-Benitez, *Lect. Notes Phys.* **2013**, *871* 433.

[5] X.-G. Huang, *Rept. Prog. Phys.* **2016**, *79*, 7 076302.

[6] K. Landsteiner, *Acta Phys. Polon. B* **2016**, *47* 2617.

[7] D. M. Capper, M. J. Duff, *Nucl. Phys. B* **1974**, *82* 147.

[8] D. M. Capper, M. J. Duff, L. Halpern, *Phys. Rev. D* **1974**, *10* 461.

[9] S. Deser, M. J. Duff, C. J. Isham, *Nucl. Phys. B* **1976**, *111* 45.

[10] M. N. Chernodub, Y. Ferreiros, A. G. Grushin, K. Landsteiner, M. A. H. Vozmediano, *Phys. Rept.* **2022**, *977* 1.

[11] M. N. Chernodub, *Phys. Rev. Lett.* **2016**, *117*, 14 141601.

[12] J. M. Luttinger, *Phys. Rev.* **1964**, *135* A1505.

[13] M. N. Chernodub, A. Cortijo, M. A. H. Vozmediano, *Phys. Rev. Lett.* **2018**, *120*, 20 206601.

[14] M. N. Chernodub, C. Corianò, M. M. Maglio, *Phys. Lett. B* **2020**, *802* 135236.

[15] B. Bermond, M. Chernodub, A. G. Grushin, D. Carpentier **2022**.

[16] C. Northe, C. Zhang, R. Wawrzyńczak, J. Gooth, S. Galeski, E. M. Hankiewicz **2022**.

[17] C.-S. Chu, R.-X. Miao, *Phys. Rev. Lett.* **2018**, *121*, 25 251602.

[18] C.-S. Chu, R.-X. Miao, *JHEP* **2018**, *07* 005.

[19] D. M. McAvity, H. Osborn, *Class. Quant. Grav.* **1991**, *8* 603.

[20] M. N. Chernodub, M. A. H. Vozmediano, *Phys. Rev. Research.* **2019**, *1* 032002.

[21] C. P. Herzog, K.-W. Huang, *JHEP* **2017**, *10* 189.

[22] N. Andrei, et al., *J. Phys. A* **2020**, *53*, 45 453002.

[23] C.-S. Chu, R.-X. Miao, *Phys. Rev. D* **2020**, *102*, 4 046011.

[24] P.-J. Hu, Q.-L. Hu, R.-X. Miao, *Phys. Rev. D* **2020**, *101*, 12 125010.
[25] R.-X. Miao, Y.-Q. Zeng, *Phys. Lett. B* **2023**, *838* 137700.

[26] M. Tinkham, *Introduction to Superconductivity*, Dover Books on Physics Series. Dover Publications, **2004**.

[27] M. N. Chernodub, V. A. Goy, A. V. Molochkov, *Phys. Lett. B* **2019**, *789* 556.

[28] C.-S. Chu, R.-X. Miao, *Phys. Rev. Lett.* **2018**, *121*, 252502.

[29] Y. Nakayama, *Phys. Rept.* **2015**, *569* 1.

[30] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lutgemeier, B. Petersson, *Nucl. Phys. B* **1996**, *469* 419.

[31] J. Ranft, J. Kripfganz, G. Ranft, *Phys. Rev. D* **1983**, *28* 360.

[32] C. Gattringer, C. Lang, *Quantum chromodynamics on the lattice: an introductory presentation*, volume 788, Springer Science & Business Media, **2009**.

[33] I. Omelyan, I. Mryglod, R. Folk, *Physical Review E* **2002**, *65*, 056706.

[34] G. S. Bali, F. Bruckmann, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg, A. Schafer, K. K. Szabo, *JHEP* **2012**, *02* 044.

[35] M. B. Einhorn, R. Savit, *Phys. Rev. D* **1978**, *17* 2583.

[36] C.-S. Chu, R.-X. Miao, *JHEP* **2018**, *07* 005.

[37] H. Isobe, N. Nagaosa, *Physical Review B* **2012**, *86*, 16.

[38] P. Maletinsky, S. Hong, M. S. Grinolds, B. Hausmann, M. D. Lukin, R. L. Walsworth, M. Loncar, A. Yacoby, *Nature Nanotechnology* **2012**, *7*, 5320.

[39] U. Vool, A. Hamo, G. Varnavides, Y. Wang, T. X. Zhou, N. Kumar, Y. Dovzhenko, Z. Qiu, C. A. C. Garcia, A. T. Pierce, J. Gooth, P. Anikeeva, C. Felser, P. Narang, A. Yacoby, *Nature Physics* **2021**, *17*, 11216.

[40] Y. Shi, J. Kahn, B. Niu, Z. Fei, B. Sun, X. Cai, B. A. Francisco, D. Wu, Z.-X. Shen, X. Xu, D. H. Cobden, Y.-T. Cui, *Science Advances* **2019**, *5*, 2.