2309.10025v2

# Emergent Chern-Simons Interactions in 3+1 Dimensions **Authors**: M. Stålhammar, D. Rudneva, T.H. Hansson, F. Wilczek1,2,3,4,512345,, 5 > 33{}^{3}start_FLOATSUPERSCRIPT 3 end_FLOATSUPERSCRIPTCenter for Theoretical Physics, MIT, Cambridge, Massachusetts 02139 USA44{}^{4}start_FLOATSUPERSCRIPT 4 end_FLOATSUPERSCRIPTT. D. Lee Institute and Wilczek Quantum Center,Shanghai Jiao Tong University, Shanghai, China55{}^{5}start_FLOATSUPERSCRIPT 5 end_FLOATSUPERSCRIPTArizona State University, Tempe, Arizona 85287, USA Abstract Parity-violating superconductors can support a low-dimension local interaction that becomes, upon condensation, a purely spatial Chern-Simons term. Solutions to the resulting generalized London equations can be obtained from solutions of the ordinary London equations with a complex penetration depth, and suggest several remarkable physical phenomena. The problem of flux exclusion by a sphere brings in an anapole moment, the problem of current-carrying wires brings in an azimuthal magnetic field, and the problem of vortices brings in currents along the vortices. We demonstrate that interactions of this kind, together with a conceptually related dimensionally reduced Chern-Simons interaction, can arise from physically plausible microscopic interactions. pacs: 74.20.De, 74.20.Rp, 03.65.Vf I Introduction The principles of symmetry and locality allow us to survey interactions that are likely to emerge in the description of materials at low energy in a systematic way. This Ginzburg-Landau or effective field theory approach has proved to be a fruitful guide to low-energy dynamics, phase transitions, and response to external fields in many applications. In this approach, the focus is on possible interactions represented by local operators of low mass dimension. In this paper we shall study terms in the effective Lagrangian, specifically $\vec{\beta}·\vec{A}×\vec{E}-A_{0}\vec{\beta}·\vec{B}$ and especially $\vec{A}·\vec{B}$ , that break discrete symmetries. For general reasons related to gauge symmetry it appears that the second of these terms cannot appear in a normal material in its thermodynamic ground state Nazarov (1986). However, as demonstrated below, it is permitted in superconducting states that can support persistent current. We will provide examples of microscopic models that illustrate both terms. These terms have noteworthy phenomenological implications, which we will explore further and exemplify below. Interactions mediated by Lagrangian densities of the Chern-Simons form $$ {\cal L}\leavevmode\nobreak\ \propto\leavevmode\nobreak\ \epsilon^{\alpha\beta% \gamma}A_{\alpha}\partial_{\beta}A_{\gamma}, \tag{1} $$ involving a gauge field $A_{\alpha}$ , and their multi-field and non-Abelian generalizations, have attracted much attention in recent years, mostly in the context of (2+1)-dimensional systems, where they can have a topological character. In its most straightforward application, the gauge field $\vec{A}$ is the gauge field of electromagnetism. Then these terms directly induce, and parametrize, interesting aspects of electromagnetic response that, at a heuristic level, derive from current-field mixing. Here we will examine a different appearance of interactions of this kind, in 3+1 dimensions, where we take all the indices to be spatial: $$ {\cal L}_{\rm CSt}=-\frac{\beta}{2}\vec{A}\cdot\vec{B}. \tag{2} $$ Terms of this kind are not relativistically invariant, and they also violate parity (but not time reversal). That does not forbid their appearance, since many materials, such as those based on crystals lacking an inversion center, and others described below, violate those symmetries. More seriously, such terms embody only a limited form of gauge symmetry. Under a local gauge transformation $A_{\alpha}→ A_{\alpha}+∂_{\alpha}\Lambda$ , we have $$ \epsilon^{\alpha\beta\gamma}A_{\alpha}\partial_{\beta}A_{\gamma}\leavevmode% \nobreak\ \rightarrow\leavevmode\nobreak\ \epsilon^{\alpha\beta\gamma}A_{% \alpha}\partial_{\beta}A_{\gamma}\,+\,\partial_{\alpha}(\epsilon^{\alpha\beta% \gamma}\Lambda\partial_{\beta}A_{\gamma}), \tag{3} $$ so the change in the bulk interaction can be cast into a surface term; but in possible applications the surface term requires careful consideration. Notably, in the context of the quantum Hall effect it is connected to the existence of edge modes and is canceled through an anomalous surface theory Stone (1991). In addition to the spatial Chern-Simons term Eq. (2), we shall also consider the term $$ \displaystyle\mbox{$$}{\cal L}_{\rm CSs}=\frac{1}{2}\,\bigl{(}\vec{\beta}\cdot% \vec{A}\times\vec{E}-A_{0}\vec{\beta}\cdot\vec{B}\bigr{)}. \tag{4} $$ In a superconductor we can generate a spatial Chern-Simons interaction from a conventional, manifestly gauge invariant interaction by condensation, viz.: $$ {\rm Re}\,\phi^{\dagger}i\vec{D}\phi\cdot\vec{B}\leavevmode\nobreak\ % \rightarrow\leavevmode\nobreak\ qv^{2}\vec{A}\cdot\vec{B}, \tag{5} $$ where $\vec{D}\equiv\vec{∇}-iq\vec{A}$ is the covariant derivative, $\phi→\langle\phi\rangle\equiv v$ through condensation, and where we used the London gauge. This is similar to how condensation generates a photon mass term $\propto A^{2}$ from the kinetic energy $\propto\phi^{\dagger}\vec{D}^{2}\phi$ . We can expect such terms to arise even in $s$ -wave superconductors that violate parity symmetry, for example those based on chiral crystals, on organic superconductors subject to chiral selection of the base molecules, generic $s$ -wave superconductors incorporating chiral dopants, or non-centrosymmetric superconductors Samokhin (2004); Agterberg (2003); Bauer and Sigrist (2012); Samoilenka and Babaev (2020). Spatial Chern-Simons terms are also known to be directly connected to the chiral magnetic effect (CME), and arise naturally in various Weyl systems as a consequence of driving or domain wall motion Hannukainen et al. (2020, 2021); Zhong et al. (2016); Landsteiner (2016). The presence of these terms have, however, been debated in the past years, for instance, in Refs Vazifeh and Franz (2013); Başar et al. (2014) (see, also Refs. Hosur and Qi ; Landsteiner (2016) for a summary). Care must be taken to distinguish the equilibrium ground state from driven states, and also to take the correct order of the limits $\vec{q}→ 0$ and $\omega→ 0$ Chen et al. (2013); Chang and Yang (2015). In Sec. V below, we will further investigate the presence of both $\vec{A}·\vec{B}$ , and $\vec{\beta}·\vec{A}×\vec{E}$ terms in equilibrium Weyl systems. Heuristically, one identifies quantities $\vec{\tilde{j}}$ that appear in Lagrangian densities of the form $\vec{A}·\vec{\tilde{j}}$ as effective currents, since they will appear as such in the Maxwell equations. (If $\vec{\tilde{j}}$ depends explicitly on $\vec{A}$ , slight complications ensue.) Famously, the London diamagnetic current ${\vec{j}}_{d}\propto\vec{A}$ is characteristic of the superconducting photon mass. Following this heuristic, the spatial Chern-Simons term Eq. (2) gives us a current ${\vec{j}}_{\text{CS}}$ proportional to the magnetic field: $$ \vec{j}_{\text{CS}}\propto\vec{B}. \tag{6} $$ This yields unusual, interesting, and potentially important phenomenological consequences, which is the subject of Sect. II From a broader theoretical perspective, a natural term descending from a Lorentz invariant effective action is $\epsilon^{\mu\nu\sigma\omega}\beta_{\mu}A_{\nu}F_{\sigma\omega}$ , where $\beta_{\mu}$ is an axial (i.e., unnatural) four-vector. Constant values of $\beta$ violate Lorentz invariance, and a constant $\beta_{0}$ can be powerfully constrained phenomenologically using astronomical data Carroll et al. (1990). But interactions of the form Eqs. (2) and (4) arise naturally from the canonical axion coupling to electromagnetic fields $$ {\cal L}\leavevmode\nobreak\ \propto\leavevmode\nobreak\ a\epsilon^{\alpha% \beta\gamma\delta}\partial_{\alpha}A_{\beta}\partial_{\gamma}A_{\delta} \tag{7} $$ for the simplest space-time variations of $a$ , corresponding to the axion background $a=\beta t$ and $a=\vec{\beta}·\vec{x}$ , respectively, where the latter describes an “axion wind” background, i.e., one that is constant in time but varies linearly in a spatial direction. Of course, both can occur together. The expression $\vec{A}·\vec{B}$ appears in many places in the literature on magnetohydrodynamics Moffatt (1969), where it is used to characterize the magnetic field configurations in the plasmas. In particular, Eq. (6) describes a force-free field since the Lorentz force on a current parallel with the magnetic field vanishes Woltjer (1958). We stress, however, that here we are interested in the response of materials where $\vec{A}·\vec{B}$ is part of the effective action, and thus determines the response to external fields. II Phenomena in Actively Chiral Superconductors For ease of reference, and in view of their connection with chirality and optical activity, we shall refer to superconductors that incorporate a purely spatial Chern-Simons term ${\cal L}_{\text{CS}}$ as actively chiral superconductors. We will work with the Lagrangian density $$ {\cal L}\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{1}{2}E^{2}-\frac{1}{2% }B^{2}-\frac{\beta}{2}\vec{A}\cdot\vec{B}-\frac{\gamma}{2}A^{2}. \tag{8} $$ II.1 Plane waves and stability From Eq. (8) we derive, after fixing the gauge $A_{0}=0$ and adopting the plane-wave ansatz $$ \vec{A}\leavevmode\nobreak\ =\leavevmode\nobreak\ \vec{\varepsilon}\exp i(\vec% {k}\cdot\vec{x}-\omega t) \tag{9} $$ the equations of motion $$ \displaystyle\vec{k}\cdot\vec{\varepsilon}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ 0, \displaystyle(\omega^{2}-k^{2}-\gamma)\vec{\varepsilon}\,\mp\,\beta\vec{k}% \times\vec{\varepsilon}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ 0. \tag{10} $$ The eigen-polarizations are transverse and circular. Indeed, with $$ \displaystyle\vec{k}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ (0,0,k), \displaystyle\vec{\varepsilon}\leavevmode\nobreak\ \displaystyle\propto \displaystyle\leavevmode\nobreak\ (1,\pm i,0), \tag{12} $$ we find the dispersion relations $$ (\omega^{2}-k^{2}-\gamma)\mp\beta k\leavevmode\nobreak\ =\leavevmode\nobreak\ 0. \tag{13} $$ The two circular polarizations propagate with different velocities. This gives rise to optical activity, i.e., rotation of the plane of linear polarization as the (transverse) wave propagates. For stability in time we require that for real $k$ the $\omega$ that solve the dispersion relation are real. This gives us the stability condition $$ 4\gamma\geq\beta^{2}. \tag{14} $$ The same condition also ensures the positivity of the energy. Indeed, since the electric field contribution is manifestly positive, at issue is only the positivity of the magnetic energy $$ {\cal E}\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{1}{2}\,\int B^{2}+% \beta\vec{A}\cdot\vec{B}+\gamma A^{2}. \tag{15} $$ We can write this as $$ \displaystyle{\cal E}=\frac{1}{4}\,\int \displaystyle\left(1+\frac{\beta}{2\sqrt{\gamma}}\right)\left(\vec{B}+\sqrt{% \gamma}\vec{A}\right)^{2}\,+ \displaystyle\left(1-\frac{\beta}{2\sqrt{\gamma}}\right)\left(\vec{B}-\sqrt{% \gamma}\vec{A}\right)^{2}. \tag{16} $$ When Eq. (14) is satisfied the coefficients of these two manifestly positive terms will both be non-negative. The stability condition Eq. (14) requires, for $\beta≠ 0$ , that $\gamma>0$ . Thus, it requires a non-zero effective photon mass, such as we have in superconductivity. Note that if the condition Eq. (14) is relaxed, the modes at very low $k$ will still be stable since the $A^{2}$ term dominates, and the same will be true for large $k$ modes where the $B^{2}$ term dominates. There will however be a region of intermediate $k$ where the $\vec{A}·\vec{B}$ term will give an instability. In a more complete theory this instability could be cured by higher order terms and this would open for a non-trivial magnetic structure in the ground state. We will revisit this subject in a slightly different context in the Sec. III. II.2 Connection to Optical Activity The optical activity of a material is usually described in terms of a frequency and momentum dependent dielectric constant and/or magnetic permeability. We shall consider the latter case and write the magnetic energy density as $$ \displaystyle\mbox{$$}{\mathcal{E}}_{m} \displaystyle= \displaystyle\frac{1}{2}B^{i}\mu^{-1}_{ij}(\omega,\vec{\nabla})B^{j} \displaystyle\rightarrow \displaystyle\frac{1}{2\mu}B^{2}+\frac{\alpha(\omega)}{2}\vec{B}\cdot\vec{% \nabla}\times\vec{B}\,, \tag{17} $$ where in the second line we put $\mu_{ij}(\omega,\vec{∇})=\mu(\delta_{ij}-\alpha(\omega)\epsilon_{ikj}% ∇^{k})$ and expanded to leading order in $\vec{∇}$ . The energy corresponding to the second term can be rewritten as | | $\displaystyle\mbox{$$}E_{c}$ | $\displaystyle=$ | $\displaystyle\frac{\alpha(\omega)}{2}∈t_{V}d^{3}x\,(\vec{∇}×\vec{A% })·(\vec{∇}×\vec{B})$ | | | --- | --- | --- | --- | --- | | $\displaystyle=$ | $\displaystyle-\frac{\alpha(\omega)}{2}∈t_{V}d^{3}x\,\vec{A}·∇^{2}% \vec{B}$ | | | | | $\displaystyle+$ | $\displaystyle\frac{\alpha(\omega)}{2}∈t_{\delta V}dS_{i}\,A_{j}(∂_{i}% B_{j}-∂_{j}B_{i}).$ | | | | We already mentioned that the first term on the second line is not gauge invariant, but in this case it is easy to show that the surface term in the last line, as expected, restores gauge invariance, so there is no need for an additional surface theory. Let us now assume that we have a superconductor with randomly implanted optically active impurities, that will add a term $E_{c}$ to the free energy functional of the superconductor. To leading order in $\alpha$ , we can then just substitute the London relation $∇^{2}\vec{B}=\lambda_{L}^{-2}\vec{B}$ in Eq. (17), and assuming that $\alpha(\omega)$ can be approximated by a constant $\alpha$ at low frequencies we obtain the low-energy Lagrangian in Eq. (8) if we identify $\beta=\alpha(0)/\lambda_{L}^{2}$ . II.3 Solution schema We are interested in solving the equation $$ \vec{\nabla}\times\vec{\nabla}\times\vec{B}+\beta\vec{\nabla}\times\vec{B}+% \gamma\vec{B}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0. \tag{19} $$ Equation (19) is a generalization of the famous London equation for superconducting magnetostatics, which is the special case $\beta=0$ . In the London equation, $\gamma$ represents the inverse square of the penetration depth. As we now demonstrate, one can generate solutions to Eq. (19) out of solutions to the London equation with a complex coefficient. Indeed, inserting the superposition ansatz $$ \vec{B}\leavevmode\nobreak\ =\leavevmode\nobreak\ \vec{B}_{a}\,+\,\kappa\vec{% \nabla}\times\vec{B}_{a} \tag{20} $$ into Eq. (19) leads to $$ \displaystyle\left(1+\beta\kappa\right)\vec{\nabla}\times\vec{\nabla}\times% \vec{B}_{a}\,+\,\gamma\vec{B}_{a} \displaystyle+\kappa\vec{\nabla}\times\vec{\nabla}\times\vec{\nabla}\times\vec% {B}_{a}\,+\,\left(\beta+\gamma\kappa\right)\vec{\nabla}\times\vec{B}_{a}\, \displaystyle\leavevmode\nobreak\ =\leavevmode\nobreak\ 0, \tag{21} $$ and therefore when $$ \vec{\nabla}\times\vec{\nabla}\times\vec{B}_{a}+\alpha\vec{B}_{a}\leavevmode% \nobreak\ =\leavevmode\nobreak\ 0, \tag{22} $$ to $$ \displaystyle\left[-\alpha\,\left(1+\beta\kappa\right)\,+\,\gamma\,\right]\,% \vec{B}_{a} \displaystyle\ +\left(\,-\alpha\kappa+\beta+\gamma\kappa\,\right)\,\vec{\nabla% }\times\vec{B}_{a}\leavevmode\nobreak\ =\leavevmode\nobreak\ 0. \tag{23} $$ Thus, if we enforce the algebraic relations $$ \displaystyle-\alpha\,(1+\beta\kappa)\,+\,\gamma\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ 0, \displaystyle-\alpha\kappa+\beta+\gamma\kappa\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ 0, \tag{24} $$ then $\vec{B}$ will satisfy the generalized London equation Eq. (19). We are given $\beta,\gamma$ and seek to solve for $\alpha,\kappa$ . From Eqs. (II.3) we derive a quadratic equation for $\kappa$ , that is solved by $$ \kappa\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{-\beta\,\pm\,i\sqrt{4% \gamma-\beta^{2}}}{2\gamma}. \tag{25} $$ Here we see that the realistic situation $4\gamma-\beta^{2}>0$ brings in complex numbers. Having gotten $\kappa$ in terms of $\beta,\gamma$ it is straight forward to further arrive at $$ \alpha\leavevmode\nobreak\ =\leavevmode\nobreak\ \gamma\,+\,\frac{\beta}{2}% \left(-\beta\,\mp\,i\sqrt{4\gamma-\beta^{2}}\,\right). \tag{26} $$ Thus, we have two complex conjugate solutions for our auxiliary inverse square penetration depth. An immediate physical implication is that we can expect oscillations to accompany the exponential damping of fields (and currents) we usually encounter as we penetrate a superconductor. Ultimately we want real solutions of our field equations. Since our auxiliary equations are linear, we can simply use the real and imaginary parts of their solutions. Note that since the auxiliary equations are complex conjugates of one another, they both lead us to the same real fields. This solution scheme embodies in a precise form the concept of field-current mixing that we anticipated heuristically. Indeed, since $\vec{∇}×\vec{B}_{a}$ is the London diamagnetic current associated to $\vec{B}_{a}$ , the solution $\vec{B}$ defined in Eq. (20) is a linear combination of its field and current. Let us further note the curious fact that our construction in Eq. (20) leads to (complex-valued) fields $\vec{B}$ that, like $\vec{B}_{a}$ , satisfy Eq. (22) and thus, in view of Eq. (19), $$ \beta\,\vec{\nabla}\times\vec{B}\leavevmode\nobreak\ =\leavevmode\nobreak\ (% \alpha-\gamma)\,\vec{B}, \tag{27} $$ or $$ \vec{\nabla}\times\vec{B}\leavevmode\nobreak\ =\frac{1}{2}\left(-\beta\,\mp\,i% \sqrt{4\gamma-\beta^{2}}\right)\,\vec{B}. \tag{28} $$ In the critical case $4\gamma-\beta^{2}=0$ , we get force-free fields. The solutions to the generalized London equation govern the magnetostatics for parity-violating superconductors. Examples are the noncentrosymmetric superconductors Samokhin (2004); Agterberg (2003); Bauer and Sigrist (2012); Samoilenka and Babaev (2020), which allows for Lifshitz invariants in the Ginzburg-Landau free energy. These terms are linear in covariant derivatives and magnetic fields, and can give rise to to helical phases and effects similar to the CME, see Ref. Bauer and Sigrist (2012) for an extensive summary on the topic. We now turn to the explicit solutions of the generalized London equations in various geometries with the purpose to demonstrate how pertinent effects, such as as induced currents and magnetic field profiles, are affected by the geometry of the superconductor. II.4 Slab geometry A relatively simple, yet physically significant and mathematically transparent situation to analyze is the half space or slab geometry. Thus, we imagine our superconductor to fill the half space $x>0$ while in the remaining half-space we have a constant magnetic field $$ \vec{B}^{\rm ext.}=\hat{z}B_{0},\ \ \ \ \ (x<0). \tag{29} $$ Here we match onto the solution of the ordinary London equation Eq. (22) proportional to $$ \vec{B}_{a}\leavevmode\nobreak\ =\leavevmode\nobreak\ \hat{z}B_{0}e^{-\sqrt{% \alpha}x}, \tag{30} $$ whose curl $$ \vec{\nabla}\times\vec{B}_{a}\leavevmode\nobreak\ =\leavevmode\nobreak\ \hat{y% }B_{0}\sqrt{\alpha}e^{-\sqrt{\alpha}x} \tag{31} $$ is the diamagnetic screening current. We will build our solution using this auxiliary form with $B_{0}=1$ . We invoke Eq. (26) to choose $$ \sqrt{\alpha}\leavevmode\nobreak\ =\leavevmode\nobreak\ \sqrt{\gamma-\frac{% \beta^{2}}{4}}+i\frac{\beta}{2}\leavevmode\nobreak\ \equiv\leavevmode\nobreak% \ p+iq, \tag{32} $$ where positive square roots are understood throughout. Note that in order to get solutions that fall off as $x→∞$ , we must take roots with a positive real part; the remaining choice associated with the $\mp$ in $\alpha$ has no effect on our final result, and we have chosen the lower sign. With that preparation, we can use our solution scheme to solve the generalized London equation Eq. (19). After some algebra, we arrive at $$ \displaystyle B_{z}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ e^{-px}\cos qx, \displaystyle B_{y}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -e^{-px}\sin qx, \tag{33} $$ $$ \displaystyle B_{z}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -e^{-px}\sin qx, \displaystyle B_{y}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -e^{-px}\cos qx, \tag{34} $$ for the real and imaginary parts. Finally, to ensure continuity of the magnetic field at the boundary we take the linear combination of these two solutions that has $B_{y}(0)=0$ and $B_{z}(0)=B_{0}$ . This gives us the magnetic field $$ \displaystyle B_{z}(x)\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ B_{0}e^{-px}\cos qx, \displaystyle B_{y}(x)\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -B_{0}e^{-px}\sin qx, \tag{35} $$ and the current $$ \displaystyle j_{z}(x) \displaystyle=\leavevmode\nobreak\ \frac{\partial B_{y}}{\partial x}=-B_{0}e^{% -px}(-p\sin qx+q\cos qx), \displaystyle j_{y}(x) \displaystyle=\leavevmode\nobreak\ -\frac{\partial B_{z}}{\partial x}=B_{0}e^{% -px}(p\cos qx+q\sin qx), \tag{36} $$ inside the superconductor. As anticipated, this solution displays three qualitatively new features relative to the usual London ( $\beta=0$ ) case. Most profoundly, there is a current running parallel to the external field direction. Secondly, there is an induced perpendicular magnetic field in the interior. Thirdly, the interior fields and currents have an oscillatory character. We can also consider a slab, occupying the region $0≤ x≤ a$ . We can use the same solution inside the superconductor, matched to the constant field $$ \displaystyle B_{z}(x\geq a)\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -B_{0}e^{-pa}\sin qa, \displaystyle B_{y}(x\geq a)\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -B_{0}e^{-pa}\cos qa. \tag{37} $$ This represents the result of applying $B=B_{0}\hat{z}$ at $x≤ 0$ and a rotated (and damped) field at $x≥ a$ . Here we see a close analogy, in magnetostatics, to optical activity (accompanied by absorption). II.5 Intrinsic solenoid (trapped flux and model vortex) We can notionally insert a solenoid into our superconductor, and ask that it be generated self-consistently by screening currents within the superconductor. This is of interest in itself, and also allows us to anticipate and model, within the relatively simple and parameter-sparse context of the (modified) London equations, properties of trapped flux and of quantized magnetic vortices. <details> <summary>2309.10025v2/extracted/5431462/Current_Magnetic_Field_Cylinder_small.jpg Details</summary>  ### Visual Description # Technical Document Extraction: Oscillatory Behavior Analysis ## Image Overview The image contains four subplots (a)-(d) depicting oscillatory behavior of physical quantities as functions of radial distance \( r \). All subplots share identical parameters: \( \beta = 15 \), \( \gamma = 57 \). Each plot uses a logarithmic y-axis scale, with distinct variables and amplitude ranges. --- ### Subplot (a): \( B_\phi \) Oscillation **Title**: \( B_\phi \) Oscillation (\( \beta = 15, \gamma = 57 \)) **Axes**: - \( x \)-axis: \( r \) (radial distance), range [1, 5] - \( y \)-axis: \( B_\phi \times 10^1 \), range [-8, 0] **Key Observations**: 1. **Trend**: Oscillatory waveform with decreasing amplitude over \( r \). 2. **Data Points**: - \( r = 1 \): \( B_\phi = 0 \times 10^1 \) - \( r = 2 \): \( B_\phi \approx -8 \times 10^1 \) (minimum) - \( r = 3 \): \( B_\phi \approx 0 \times 10^1 \) (zero-crossing) - \( r = 4 \): \( B_\phi \approx 8 \times 10^1 \) (peak) - \( r = 5 \): \( B_\phi \approx 0 \times 10^1 \) (zero-crossing) **Legend**: None present. --- ### Subplot (b): \( B_z \) Oscillation **Title**: \( B_z \) Oscillation (\( \beta = 15, \gamma = 57 \)) **Axes**: - \( x \)-axis: \( r \), range [1, 5] - \( y \)-axis: \( B_z \times 10^1 \), range [-6, 1] **Key Observations**: 1. **Trend**: Damped oscillations with amplitude reduction. 2. **Data Points**: - \( r = 1 \): \( B_z = 1 \times 10^1 \) (peak) - \( r = 2 \): \( B_z \approx -6 \times 10^1 \) (minimum) - \( r = 3 \): \( B_z \approx 0 \times 10^1 \) (zero-crossing) - \( r = 4 \): \( B_z \approx 6 \times 10^1 \) (peak) - \( r = 5 \): \( B_z \approx 0 \times 10^1 \) (zero-crossing) **Legend**: None present. --- ### Subplot (c): \( j_\phi \) Oscillation **Title**: \( j_\phi \) Oscillation (\( \beta = 15, \gamma = 57 \)) **Axes**: - \( x \)-axis: \( r \), range [1, 5] - \( y \)-axis: \( j_\phi \times 10^2 \), range [-6, 0] **Key Observations**: 1. **Trend**: Sustained oscillations with minimal damping. 2. **Data Points**: - \( r = 1 \): \( j_\phi = 0 \times 10^2 \) - \( r = 2 \): \( j_\phi \approx -6 \times 10^2 \) (minimum) - \( r = 3 \): \( j_\phi \approx 0 \times 10^2 \) (zero-crossing) - \( r = 4 \): \( j_\phi \approx 6 \times 10^2 \) (peak) - \( r = 5 \): \( j_\phi \approx 0 \times 10^2 \) (zero-crossing) **Legend**: None present. --- ### Subplot (d): \( j_z \) Oscillation **Title**: \( j_z \) Oscillation (\( \beta = 15, \gamma = 57 \)) **Axes**: - \( x \)-axis: \( r \), range [1, 5] - \( y \)-axis: \( j_z \times 10^2 \), range [-4, 4] **Key Observations**: 1. **Trend**: Persistent oscillations with consistent amplitude. 2. **Data Points**: - \( r = 1 \): \( j_z = 4 \times 10^2 \) (peak) - \( r = 2 \): \( j_z \approx -4 \times 10^2 \) (minimum) - \( r = 3 \): \( j_z \approx 0 \times 10^2 \) (zero-crossing) - \( r = 4 \): \( j_z \approx 4 \times 10^2 \) (peak) - \( r = 5 \): \( j_z \approx 0 \times 10^2 \) (zero-crossing) **Legend**: None present. --- ### Cross-Subplot Analysis - **Parameter Consistency**: All subplots share \( \beta = 15 \), \( \gamma = 57 \), suggesting these parameters govern the oscillatory dynamics. - **Amplitude Scaling**: - \( B_\phi \) and \( B_z \): \( \times 10^1 \) scale (subplots a, b). - \( j_\phi \) and \( j_z \): \( \times 10^2 \) scale (subplots c, d). - **Waveform Similarity**: All subplots exhibit sinusoidal-like oscillations with phase shifts. --- ### Missing Elements - **Legend**: No legend is present to correlate line colors (blue) with data series. - **Units**: Explicit units for \( B_\phi, B_z, j_\phi, j_z \) are not provided (e.g., Tesla, Amperes). --- ### Conclusion The subplots collectively illustrate oscillatory behavior of magnetic (\( B_\phi, B_z \)) and current density (\( j_\phi, j_z \)) components as functions of radial distance \( r \). Damping effects are evident in \( B_\phi \) and \( B_z \), while \( j_\phi \) and \( j_z \) maintain sustained oscillations. Further analysis requires unit specifications and legend clarification. </details> FIG. 1: Azimuthal and longitudinal components of the magnetic field and the currents inside the superconductor with a solenoid for $\beta=15$ , $\gamma=57$ , $B_{0}=100$ , and $R=1$ , displaying a (spatially damped) longitudinal current. Notably, the applied, constant and longitudinal magnetic field in the solenoid gives rise to an azimuthal field component inside the superconductor that further oscillates in a damped fashion, just as the longitudinal component, as a function of $r$ . In cylindrical coordinates, our solenoid is defined by $$ \vec{B}(r,z,\phi)\leavevmode\nobreak\ =\leavevmode\nobreak\ B_{0}\hat{z},\ \ % \ r\leq R, \tag{38} $$ and it joins on to a solution of Eq. (19) for $r≥ R$ . The self-consistency condition is that there are no singular surface currents, which we enforce by demanding continuity of the tangential magnetic fields at $r=R$ . <details> <summary>2309.10025v2/extracted/5431462/Spherical_Currents_2.jpg Details</summary>  ### Visual Description # Technical Document Extraction ## Panel (a) - **Graph Title**: Not explicitly labeled - **Parameters**: β = 15, γ = 65 - **Axes**: - **X-axis**: r (ranging from 0 to 1) - **Y-axis**: Bφ (ranging from -20 to 20) - **Key Trends**: - Sinusoidal curve with a single peak at r ≈ 0.5 - Bφ oscillates between ~-20 and ~20 - Baseline at Bφ = 0 ## Panel (d) - **Graph Title**: Not explicitly labeled - **Parameters**: β = 20, γ = 200 - **Axes**: - **X-axis**: r (ranging from 0 to 1) - **Y-axis**: Bφ (scaled ×10³, ranging from 0 to 4×10³) - **Key Trends**: - Single sharp peak at r ≈ 0.8 - Bφ rises from 0 to ~4×10³, then drops sharply - Minimal baseline variation ## Panel (b) - **Streamline Plot**: - **Parameters**: β = 15, γ = 65, x = 0 - **Axes**: - **X-axis**: y (ranging from -2 to 2) - **Y-axis**: z (ranging from -2 to 2) - **Flow Characteristics**: - Circular flow pattern around a central object - Streamlines form concentric loops with inward/outward arrows - Central object appears as a shaded gray region ## Panel (c) - **Streamline Plot**: - **Parameters**: β = 15, γ = 65, x = 0.5 - **Axes**: - **X-axis**: y (ranging from -2 to 2) - **Y-axis**: z (ranging from -2 to 2) - **Flow Characteristics**: - Distorted flow pattern compared to x = 0 - Streamlines converge toward the central object - Increased density near the object's boundary ## Panel (e) - **Streamline Plot**: - **Parameters**: β = 20, γ = 200, x = 0 - **Axes**: - **X-axis**: y (ranging from -2 to 2) - **Y-axis**: z (ranging from -2 to 2) - **Flow Characteristics**: - Similar circular pattern to panel (b) - Streamlines are denser and more tightly packed - Central object retains shaded gray appearance ## Panel (f) - **Streamline Plot**: - **Parameters**: β = 20, γ = 200, x = 0.5 - **Axes**: - **X-axis**: y (ranging from -2 to 2) - **Y-axis**: z (ranging from -2 to 2) - **Flow Characteristics**: - Highly distorted flow with elongated streamlines - Central object shows asymmetric flow distortion - Increased turbulence near the object's edges ## Cross-Reference Summary | Panel | β | γ | x | Flow Pattern | Bφ Behavior | |-------|-----|-----|-------|-----------------------|---------------------------| | (a) | 15 | 65 | N/A | Sinusoidal | Single peak at r ≈ 0.5 | | (d) | 20 | 200 | N/A | Sharp peak | Single sharp peak at r ≈ 0.8 | | (b) | 15 | 65 | 0 | Circular | N/A | | (c) | 15 | 65 | 0.5 | Distorted circular | N/A | | (e) | 20 | 200 | 0 | Circular (denser) | N/A | | (f) | 20 | 200 | 0.5 | Distorted elongated | N/A | ## Notes - All streamline plots use the same coordinate system (y, z) with identical axis ranges (-2 to 2). - Central object geometry is consistent across panels but flow dynamics vary with β, γ, and x. - No legends explicitly link colors to variables; shading of central object is uniform gray. - No textual annotations or legends present in the diagrams. </details> FIG. 2: Expulsion of the magnetic field by a superconducting unit sphere for representative parameter values. The magnetic field forms closed loops inside the superconducting sphere. Panels (a) and (d) display how the parity violating azimuthal component oscillates inside the sphere, for polar angle $\theta=\frac{\pi}{3}$ . The field configuration gives rise to an anapole moment, as discussed in the text. The auxiliary solution of the ordinary London equation brings in the Bessel function $K_{1}$ , which dies exponentially at infinity: $$ B_{a}(r,z,\phi)\leavevmode\nobreak\ =\leavevmode\nobreak\ B_{0}\hat{z}\,\frac{% K_{1}({\sqrt{\alpha}}r)}{K_{1}({\sqrt{\alpha}}R)},\ \ \ r\geq R. \tag{39} $$ From the curl of this field, we infer the azimuthal diamagnetic screening current. With this starting point, we can invoke the machinery of our solution schema to generate solutions of our generalized London equation in the exterior (superconducting) region. Details are spelled out in Appendix B. Let us mention how the qualitative novelties we observed above get manifested here: within the superconductor we find longitudinal current flows $j_{z}$ , azimuthal magnetic fields $B_{\phi}$ , and oscillatory behavior (possibly damped) of all the fields and currents as functions of $r$ . These general findings are displayed in Fig. 1 for some exemplary values of $\beta$ and $\gamma$ . It is possible to consider cylindrical shells, fields imposed from the outside, and so forth, both analytically and numerically (and, presumably, experimentally), based on the same ideas. II.6 Sphere geometry Another accessible problem, often considered to be the paradigmatic Meissner effect, is the superconducting sphere exposed to a constant external magnetic field. The auxiliary reference problem here was solved and presented by London himself in his classic book London (1950). One finds, for spheres much larger than the penetration depth, the field canceled or expelled by azimuthal diamagnetic screening currents near the surface of the sphere. Magnetic field lines with the superconductor get routed to within that penetration region, as displayed by Fig. 2. In additional to the imposed field, one finds a calculable magnetic dipole arising from the circulating currents. Since the auxiliary solution is expressed in terms of exponentials, we can use our solution schema to generate completely explicit solutions of the modified equations in terms of exponentials and trigonometric functions. The behavior of the magnetic field is illustrated in Fig. 2 for various penetration depths, and calculational details along with the full solution are spelled out in Appendix B. In Fig. 3, we see that the currents that run along the surface of the sphere and return in a (squashed) toroidal fashion give no external moment, but represent a form of what are called anapole moments in the literature. An anapole moment, or a magnetic toroidal moment, is a term in the multipole expansion of the electromagnetic field that violates both $\mathcal{P}$ and $\mathcal{T}$ symmetry. The anapole moment is given by $$ T_{i}=\frac{1}{10}\int\left[r_{i}\left(\vec{r}\cdot\vec{j}\right)-2r^{2}J_{i}% \right]d^{3}x, \tag{40} $$ where $r_{i}$ are the Cartesian coordinates, and $\vec{j}$ the current. Using the explicit solutions of the magnetic field inside the sphere (see Appendix B for their explicit appearance), the current is given by $\vec{j}_{\text{sphere}}=\vec{∇}×\vec{B}^{\text{in}}_{\text{sphere}}$ , it can be shown that $T_{x}$ and $T_{y}$ are identically zero. However, $T_{z}$ is finite and for the two cases illustrated in Fig. 2 given by, $$ \displaystyle T_{z}\left(\beta=15,\gamma=65\right) \displaystyle=-15.5804, \displaystyle T_{z}\left(\beta=20,\gamma=200\right) \displaystyle=2182.52. \tag{41} $$ Appendix B contains a general expression for the $z$ component of the anapole moment. <details> <summary>2309.10025v2/extracted/5431462/Spherical_Currents_new.jpg Details</summary>  ### Visual Description # Technical Document Extraction: Subplot Analysis ## Subplot (a) - **Title**: β = 15, γ = 65 - **Y-axis**: \( j_r \times 10^8 \) - **X-axis**: \( r \) (0 to 1) - **Data**: Smooth curve decreasing from ~8×10⁸ at \( r=0 \) to ~2×10⁸ at \( r=1 \), with a concave-up inflection near \( r=0.5 \). ## Subplot (b) - **Title**: β = 15, γ = 65 - **Y-axis**: \( j_\theta \times 10^2 \) - **X-axis**: \( r \) (0 to 1) - **Data**: Smooth curve starting at ~-2×10² at \( r=0 \), peaking at ~2×10² near \( r=0.5 \), then declining to ~-2×10² at \( r=1 \). ## Subplot (c) - **Title**: β = 15, γ = 65 - **Y-axis**: \( j_\phi \times 10^2 \) - **X-axis**: \( r \) (0 to 1) - **Data**: Smooth curve starting at 0 at \( r=0 \), dipping to ~-2×10² near \( r=0.5 \), then rising to ~2×10² at \( r=1 \). ## Subplot (d) - **Title**: β = 20, γ = 200 - **Y-axis**: \( j_r \times 10^3 \) - **X-axis**: \( r \) (0 to 1) - **Data**: Nearly flat curve with a shallow concave-up shape, rising from ~0.5×10³ at \( r=0 \) to ~5×10³ at \( r=1 \). ## Subplot (e) - **Title**: β = 20, γ = 200 - **Y-axis**: \( j_\theta \times 10^4 \) - **X-axis**: \( r \) (0 to 1) - **Data**: Nearly flat curve with a shallow concave-down shape, dipping from ~0.5×10⁴ at \( r=0 \) to ~-2×10⁴ at \( r=1 \). ## Subplot (f) - **Title**: β = 20, γ = 200 - **Y-axis**: \( j_\phi \times 10^4 \) - **X-axis**: \( r \) (0 to 1) - **Data**: Nearly flat curve with a sharp concave-down drop, starting at ~0 at \( r=0 \), rising slightly to ~2×10⁴ near \( r=0.5 \), then plunging to ~-6×10⁴ at \( r=1 \). ## Key Observations 1. **Parameter Impact**: - Higher β (20 vs. 15) and γ (200 vs. 65) correlate with increased magnitude in \( j_r \), \( j_\theta \), and \( j_\phi \). - \( j_\phi \) exhibits the most dramatic change under higher γ (e.g., -6×10⁴ vs. -2×10²). 2. **Trend Consistency**: - All \( j_r \) curves show monotonic or near-monotonic behavior. - \( j_\theta \) and \( j_\phi \) exhibit oscillatory patterns, with amplitude scaling with γ. 3. **Scale Sensitivity**: - \( j_r \) and \( j_\theta \) magnitudes increase by ~100× when β and γ double. - \( j_\phi \) scales by ~30× under the same parameter changes. </details> FIG. 3: Currents in spherical coordinates inside a superconducting unit sphere generated by the external and constant magnetic, for polar angle $\theta=\frac{\pi}{3}$ . Note that the radial component of the current vanishes at the boundary. No current escapes from the sphere, but there is an anapole moment. III Actively Chiral Magnetism Chiral materials generally, and not only superconductors, support a $\vec{B}·(\vec{∇}×\vec{B})$ term in the effective Lagrangian density. At optical frequencies, it gives optical activity. One should also consider its effect in magneto-statics. In that context, the most relevant terms are $$ -{\cal L}\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{1}{2\mu}B^{2}+\kappa% \vec{B}\cdot\left(\vec{\nabla}\times\vec{B}\right)+\frac{\lambda}{2}\left(\vec% {\nabla}\times\vec{B}\right)^{2}. \tag{43} $$ This Lagrangian density bears a close family resemblance to the effective Lagrangian we used in our analysis of chirally active superconductivity.. Indeed, the substitution $\vec{∇}×\vec{B}→\vec{A}$ brings it into the form we analyzed above. A simple heuristic consideration suggests the common occurrence of this term in organic (chiral) diamagnetism. Imagine a helical molecule segment along which diamagnetic currents can flow. If we apply magnetic flux along the axis of the helix, the diamagnetic current that along the helix, regarded vectorially, will have a component along the applied magnetic field direction. But the current sources $\vec{∇}×\vec{B}$ , so this correlation represents a $\vec{B}·\vec{∇}×\vec{B}$ . If the magnetic field is off-axis only its component along the axis will be operative, but the same logic applies. Helices of the same chirality will all contribute with the same sign. We can write the energy density as $$ \displaystyle\frac{1}{4}\left[\left(\frac{1}{\mu}+\frac{\kappa}{\sqrt{\lambda% \mu}}\right)\left(\vec{B}+\sqrt{\lambda\mu}\vec{\nabla}\times\vec{B}\right)^{2% }\right. \displaystyle\ \ \ \ \ \left.+\left(\frac{1}{\mu}-\frac{\kappa}{\sqrt{\lambda% \mu}}\right)\left(\vec{B}-\sqrt{\lambda\mu}\vec{\nabla}\times\vec{B}\right)^{2% }\right], \tag{44} $$ from which we see that we have stability for $$ \frac{\lambda}{\mu}\geq\kappa^{2}, \tag{45} $$ and of course $\mu,\lambda≥ 0$ . Alternatively we can consider plane waves $$ \vec{A}\leavevmode\nobreak\ =\leavevmode\nobreak\ \left(\begin{array}[]{c}1\\ \pm i\\ 0\end{array}\right)e^{ikz}, \tag{46} $$ with energy density proportional to $$ \frac{1}{\mu}\pm 2\kappa k+\lambda k^{2}. \tag{47} $$ Here, positivity of the energy density (for real $k$ ) leads again to Eq. (45). As long as $\lambda>0$ , we can stabilize the model by adding a $(B^{2})^{2}$ term. Taking $\lambda,\kappa,k>0$ , the minimum energy density taking into account only quadratic terms occurs, according to Eq. (47), at $$ k_{c}=\frac{\kappa}{\lambda}, \tag{48} $$ with the lower choice of sign, where it has the value $$ \varepsilon\leavevmode\nobreak\ \equiv\leavevmode\nobreak\ \frac{1}{\mu}-\frac% {\kappa^{2}}{\lambda}. \tag{49} $$ When $\varepsilon<0$ , we can lower the energy by bringing in fields of the form. $$ \displaystyle B_{x}\leavevmode\nobreak\ \displaystyle\propto \displaystyle\leavevmode\nobreak\ \cos k_{c}z, \displaystyle B_{y}\leavevmode\nobreak\ \displaystyle\propto \displaystyle\leavevmode\nobreak\ \sin k_{c}z, \displaystyle B_{z}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ 0. \tag{50} $$ Note that this instability does not require $\mu<0$ , i.e., instability toward ordinary (i.e., $k=0$ ) ferromagnetism, though of course it includes that possibility. To describe a stable system, we must bring in a $(B^{2})^{2}$ penalty term that limits the amplitude of the spontaneously developed structure. Equations (III) represent fields that are constant within $z={\rm const.}$ planes whose direction rotates periodically within the $x-y$ plane as $z$ varies. In other words, we see here magnetic fields characteristic of optical activity frozen in time. A point of interest is that because the coupling $\frac{\lambda}{2}(\vec{∇}×\vec{B})^{2}$ required to stabilize the $\kappa\vec{B}·(\vec{∇}×\vec{B})$ optical activity term contains a larger number of derivatives than the minimal Maxwell $\frac{1}{2\mu}B^{2}$ term it cannot be regarded as a uniformly small perturbation, even when $\lambda$ is small. Indeed, it changes the nature of the boundary value problem. If we fix the gauge $\vec{∇}·\vec{A}=0$ and (for simplicity) set $\kappa=0$ , varying ${\cal L}$ leads to the equation $$ \left[\frac{1}{\mu}+\lambda\left(\nabla^{2}\right)\,\right]\,\nabla^{2}\vec{A}% \leavevmode\nobreak\ =\leavevmode\nobreak\ 0. \tag{51} $$ Any harmonic vector field will solve this equation, and very naively one might expect that for small $\lambda$ such fields provide excellent approximate solutions in general. But in regions where $\vec{A}$ varies rapidly the second term comes in strongly, and other solutions may be physically appropriate. In particular, the boundary conditions one must apply at surfaces where the value of $\lambda$ changes (notably, at boundaries between our chiral magnetic materials and conventional materials, or empty space) one must enforce additional continuity of normal derivatives, beyond what is usually required for harmonic fields, and this may require substantial adjustments of candidate solutions near the boundary. Closely related mathematical issues arise in hydrodynamics, where they have stimulated the development of boundary layer theory. The solution schema we used in the $\vec{A}·\vec{B}$ problem continues to work in this new context, so we can leverage known solutions of Eq. (51) to get solutions of the full ( $\kappa≠ 0$ ) equations in various geometries. IV Phenomenology of Axion Wind Materials The axion wind term $$ {\cal L}_{w}\ \propto\vec{\beta}_{i}\cdot\epsilon^{i\alpha\beta\gamma}A_{% \alpha}F_{\beta\gamma}\ \propto\ \vec{\beta}\cdot(\vec{A}\times\vec{E})-A_{0}% \vec{\beta}\cdot\vec{B} \tag{52} $$ breaks rotation symmetry (and time-reversal symmetry), so its phenomenology is more complicated. Here we confine ourselves to a simple but important general observation and a calculation of its effect on wave propagation. For simplicity, we will take $\vec{\beta}=\beta\hat{z}$ . Whereas the term $\propto\vec{A}·\vec{B}$ brings in spatial derivatives in all directions, and thereby involves the material as a whole, the term $\vec{\beta}·(\vec{A}×\vec{E})$ does not bring in derivatives in the $\hat{z}$ direction, and in that sense reduces to a stack of planar terms. (Of course, other terms in the Lagrangian will link the planes.) Within each plane, we have, in effect, a (2+1) dimensional Chern-Simons theory. Indeed, through the alternative formulation $\epsilon^{3\alpha\beta\gamma}A_{\alpha}∂_{\beta}A_{\gamma}$ of our term, we see that in isolation it literally represents a stack of independent (2+1)-dimensional Chern-Simons theories. This interpretation indicates that in a bounded sample there will be massless surface modes, as in the quantum Hall effect, whose anomalies cancel the surface terms that otherwise obstruct full gauge invariance of the bulk theory. From the Lagrangian $$ \displaystyle\mbox{$$}{\mathcal{L}}=\frac{1}{2}E^{2}-\frac{1}{2}B^{2}+\frac{% \beta}{2}\left[\hat{z}\cdot\left(\vec{A}\times\vec{E}\right)-A_{0}B_{z}\right], \tag{53} $$ we derive the equations of motion $$ \displaystyle\vec{\nabla}\cdot\vec{B}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ 0, \displaystyle\vec{\nabla}\times\vec{E}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -\frac{\partial\vec{B}}{\partial t}, \displaystyle\vec{\nabla}\cdot\vec{E}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \beta B_{z}, \displaystyle\vec{\nabla}\times\vec{B}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\partial\vec{E}}{\partial t}-\beta\hat% {z}\times\vec{E}. \tag{54} $$ Thus we have effective charge and current densities $$ \displaystyle\rho_{e}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \beta B_{z}, \displaystyle\vec{j}_{e}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ -\beta\hat{z}\times\vec{E}, \tag{55} $$ that automatically satisfy the conservation equation. Having imposed the corresponding variational equation, we can set the non-dynamical field $A_{0}=0$ . For plane waves propagating in the $\hat{z}$ direction, we use the ansatz $A=\epsilon\,e^{i(kz-\omega t)}$ . We find that transverse circular polarizations lead to uncoupled dispersion relations, in the forms $$ \displaystyle\vec{A} \displaystyle=\left(\begin{array}[]{c}1\\ \pm i\\ 0\end{array}\right)e^{i(kx-\omega t)}, \displaystyle 0 \displaystyle=\omega^{2}-k^{2}\mp\beta\omega. \tag{60} $$ Here again we find that the different circular polarizations travel at different velocities, so there is optical activity. Unlike before, however, here we have no zone of instability. For plane waves propagating in the $\hat{x}$ direction, we use the ansatz $A=\epsilon e^{i(kx-\omega t)}$ . One eigenmode does not feel the new term at all: $$ \displaystyle\vec{A} \displaystyle=\left(\begin{array}[]{c}0\\ 0\\ 1\end{array}\right)e^{i(kx-\omega t)}, \displaystyle 0 \displaystyle=\omega^{2}-k^{2}. \tag{001} $$ The other eigenmode is more unusual. It is $$ \vec{A}\leavevmode\nobreak\ =\leavevmode\nobreak\ \left(\begin{array}[]{c}i% \frac{\beta}{\omega}\\ 1\\ 0\end{array}\right)e^{i(kx-\omega t)}, \tag{67} $$ with the dispersion relation $$ 0\leavevmode\nobreak\ =\leavevmode\nobreak\ \omega^{2}-k^{2}-\beta^{2}. \tag{68} $$ This dispersion relation is characteristic of a massive excitation. The polarization, which is never transverse, becomes increasingly longitudinal at low frequencies. It is straightforward, though lengthy, to calculate the general case $\vec{k}=k(\sin\theta\hat{x}+\cos\theta\hat{z})$ . Here we only record the dispersion relation: $$ \omega^{2}\leavevmode\nobreak\ =\leavevmode\nobreak\ k^{2}\,+\frac{\beta^{2}}{% 2}\pm\beta\sqrt{k^{2}\cos^{2}\theta+\frac{\beta^{2}}{4}}. \tag{69} $$ <details> <summary>2309.10025v2/extracted/5431462/Diagrams_Traingle.jpg Details</summary>  ### Visual Description # Technical Diagram Analysis ## Diagram Type Feynman diagram depicting particle interactions. ## Key Components 1. **Vertices**: - **Left Vertex**: Labeled $ A_\mu^5 $ (muon field with superscript 5). - **Top Right Vertex**: Labeled $ A_\nu $ (nu field). - **Bottom Right Vertex**: Labeled $ A_\rho $ (rho field). 2. **Lines**: - **Wavy Lines**: Represent gauge boson propagators (e.g., photons, W/Z bosons). - One wavy line connects to $ A_\mu^5 $. - Two wavy lines connect to $ A_\nu $. - Two wavy lines connect to $ A_\rho $. 3. **Central Vertex**: - A single vertex connects all three labeled vertices ($ A_\mu^5 $, $ A_\nu $, $ A_\rho $). ## Structural Details - **Flow**: - $ A_\mu^5 $ emits a wavy line to the central vertex. - The central vertex splits into two wavy lines: one to $ A_\nu $ and one to $ A_\rho $. - **Symmetry**: - $ A_\nu $ and $ A_\rho $ receive identical wavy line connections from the central vertex. ## Absent Elements - No legends, axis titles, or data tables present. - No numerical data or categorical labels beyond vertex/line identifiers. ## Interpretation The diagram illustrates a theoretical interaction where a muon field ($ A_\mu^5 $) couples to a central vertex, which mediates interactions between $ A_\nu $ and $ A_\rho $ fields via gauge bosons. The superscript "5" on $ A_\mu $ may indicate a specific representation or state in the underlying theory. </details> FIG. 4: Given an effective axial vector potential, two-photon response includes the classic VVA triangle anomaly graph. V Microscopic Models Based on Semimetals Previously, we discussed, in general terms several situations where we can expect emergent Chern-Simons terms to arise. In this section, we will describe in detail a specific construction, inspired by the appearance of anomalies in quantum field theories, that gives rise to them. The $\mathcal{P}$ or $\mathcal{T}$ -breaking terms we seek are quadratic in gauge potentials and involve cross products. A natural way for those to arise, in Feynman graphs, is through the structure $\operatorname{Tr}[(\vec{a}·\vec{\sigma})(\vec{b}·\vec{\sigma})(\vec{c}% ·\vec{\sigma})]$ when integrating over fermions in vacuum polarization loops. This is similar to the structure that gives rise to chiral anomalies in relativistic theories in even dimensions. Inspired by these thoughts, we are led to consider Dirac and Weyl materials. V.1 Axion wind terms in Weyl semimetals In WSMs the different chiralities can be split in momentum space by imposing stress. The vector separating the nodes is axial, and thus provides a candidate for the axial vector $\vec{k}$ . We now show, by an explicit calculation of the vacuum current, that such a splitting does lead to a $\sim\vec{k}·(\vec{A}×\vec{E})$ term in $S_{\text{eff}}$ Hosur and Qi . Consider an inversion symmetric Weyl semimetal with Hamiltonian $$ \displaystyle\mathcal{H} \displaystyle=\sum_{a=\pm}\sum_{\vec{q}}c_{a}^{\dagger}(\vec{q}_{a})h_{a}(\vec% {q}_{a})c_{a} \displaystyle=\sum_{\vec{q}}\tilde{\Psi}^{\dagger}(\vec{q})\begin{pmatrix}h_{-% }(\vec{q}_{-})&0\\ 0&h_{+}(\vec{q}_{+})\end{pmatrix}\tilde{\Psi}(\vec{q}), \tag{72} $$ where $$ \displaystyle h_{\pm}(\vec{q}^{\pm}) \displaystyle=q^{\pm}_{x}\sigma^{x}+q^{\pm}_{y}\sigma^{y}+q^{\pm}_{z}\sigma^{z% }-b_{0}^{\pm}\sigma^{0}, \displaystyle\vec{q}_{\pm} \displaystyle=\left[k_{x}\pm b_{x},k_{y}\pm b_{y},\mp\left(k_{z}\pm b_{z}% \right)\right], \tag{73} $$ and the wave functions are defined as $$ \tilde{\Psi}(\vec{q})=\begin{pmatrix}c_{-}(\vec{q})\\ c_{+}(\vec{q})\end{pmatrix}. \tag{75} $$ The Weyl semimetal has a pair of nodes located at $\vec{k}=±\vec{b}$ at energies $b_{0}^{±}$ . Here, we have set the product of Planck’s constant and the Fermi energy, $\hbar v_{F}=1$ . In the following, to simplify notation the sum over momenta will be left implicit. The Hamiltonian can be re-written as $$ \displaystyle\mathcal{H} \displaystyle=\tilde{\Psi}^{\dagger}\begin{pmatrix}\left(\vec{k}+\vec{b}\right% )\cdot\vec{\sigma}+b_{0}^{-}&0\\ 0&\sigma^{z}\left[-\left(\vec{k}-\vec{b}\right)\cdot\vec{\sigma}+b_{0}^{+}% \right]\sigma^{z}\end{pmatrix}\tilde{\Psi} \displaystyle=\Psi^{\dagger}\begin{pmatrix}\left(\vec{k}+\vec{b}\right)\cdot% \vec{\sigma}+b_{0}^{-}&0\\ 0&\ -\left(\vec{k}-\vec{b}\right)\cdot\vec{\sigma}+b_{0}^{+}\end{pmatrix}\Psi, \tag{78} $$ with $\Psi(\vec{q})=\left(c_{-}(\vec{q}),\sigma^{z}c_{+}(\vec{q})\right)^{\text{T}}$ . Using the chiral representation of the Weyl matrices, $$ \gamma^{0}=\begin{pmatrix}0&\sigma^{0}\\ \sigma^{0}&0\end{pmatrix},\quad\vec{\gamma}=\begin{pmatrix}0&\vec{\sigma}\\ -\vec{\sigma}&0\end{pmatrix},\quad\gamma^{5}=\begin{pmatrix}-\sigma^{0}&0\\ 0&\sigma^{0}\end{pmatrix}, \tag{82} $$ the Hamiltonian becomes $$ \displaystyle\mathcal{H} \displaystyle=\bar{\Psi}\left[-\vec{\gamma}\cdot\vec{k}+\vec{\gamma}\cdot\vec{% b}\gamma^{5}+\gamma^{0}\left(P_{-}b_{0}^{-}+P_{+}b_{0}^{+}\right)\right]\Psi, \tag{83} $$ with $P_{±}=\frac{1}{2}\left({\bf 1}±\gamma^{5}\right)$ and $\bar{\Psi}=\Psi^{\dagger}\gamma^{0}$ . Heisenberg equations of motion then yield $$ \displaystyle\left[\not{\partial}-\vec{\gamma}\cdot\vec{b}\gamma^{5}-\gamma^{0% }\left(P_{-}b_{0}^{-}+P_{+}b_{0}^{+}\right)\right]\Psi \displaystyle=0\,, \tag{84} $$ from which we can extract the Lagrangian: $$ \mathcal{L}=\bar{\Psi}\left[i\not{\partial}-\vec{\gamma}\cdot\vec{b}\gamma^{5}% -\gamma^{0}\left(P_{-}b_{0}^{-}+P_{+}b_{0}^{+}\right)\right]\,. \tag{85} $$ Coupling to the electromagnetic field via minimal coupling gives $$ \mathcal{L}=\bar{\Psi}\left[i\not{\partial}+\not{A}-\vec{\gamma}\cdot\vec{b}% \gamma^{5}-\gamma^{0}\left(P_{-}b_{0}^{-}+P_{+}b_{0}^{+}\right)\right]\Psi\,. \tag{86} $$ By shifting the zeroth component of the gauge field $A$ such that $A_{0}→ A_{0}-\frac{1}{2}\left(b_{0}^{+}+b_{0}^{-}\right)$ , and defining $b_{0}=\frac{1}{2}\left(b_{0}^{-}-b_{0}^{+}\right)$ , it can be written as $$ \mathcal{L}=\bar{\Psi}\left(i\not{\partial}+\not{A}+\not{b}\gamma^{5}\right)\Psi. \tag{87} $$ Thus, the shifts in momentum and energy of the Weyl nodes can be recast into an effective axial gauge field (of a very special form) in the Lagrangian. A constant $b_{\mu}$ can be written $b_{\mu}=∂_{\mu}\xi$ , with $\xi=b_{\mu}x^{\mu}$ . Naively, such a $b_{\mu}$ can be eliminated by the chiral rotation $\psi→ e^{i\xi\gamma^{5}}\psi$ . However, a nonzero contribution to the effective action arises from the triangle anomaly, Fig. 4. After a partial integration, it becomes (see, e.g., Ref. Nair (2006)) $$ S_{\text{top}}=\frac{1}{8\pi^{2}}\int d^{4}x\epsilon^{\mu\nu\rho\sigma}% \operatorname{Tr}\left(b_{\mu}A_{\nu}\partial_{\rho}A_{\sigma}\right). \tag{88} $$ Here, $\epsilon^{\mu\nu\rho\sigma}$ is the Levi-Civita symbol with convention $\epsilon^{0123}=1$ . We will now break down the individual components of this variation and write them out explicitly to understand their physical significance: $$ \displaystyle S_{\text{top}} \displaystyle=\frac{1}{8\pi^{2}}\int d^{4}x\left(\epsilon^{0ijk}b_{0}A_{i}% \partial_{j}A_{k}+\epsilon^{i0jk}b_{i}A_{0}\partial_{j}A_{k}+\epsilon^{ij0k}b_% {i}A_{j}\partial_{0}A_{k}+\epsilon^{ijk0}b_{i}A_{j}\partial_{k}A_{0}\right) \displaystyle=\frac{1}{8\pi^{2}}\int d^{4}x\left(b_{0}\vec{A}\cdot\vec{B}-A_{0% }\vec{b}\cdot\vec{B}+\vec{b}\cdot\vec{A}\times\vec{E}\right). \tag{89} $$ Here, we again used that $b_{\mu}$ are constants, in order to perform integration by parts. At first glance it seems that an $\vec{A}·\vec{B}$ term appears already here, and its presence has been debated in recent years Landsteiner (2016); Hosur and Qi ; Hannukainen et al. (2020, 2021); Vazifeh and Franz (2013); Başar et al. (2014); Chen et al. (2013); Chang and Yang (2015). Here, however, we have calculated the current in the vacuum, meaning that we have assumed that each of the Weyl cones are filled up exactly to the nodes, despite them being at different energies. Having finite density amounts to filling particles at the nodes to levels set by the chemical potentials $\mu_{+}$ and $\mu_{-}$ , and the total current becomes $\vec{j}=\frac{\mu_{5}-b_{0}}{2\pi^{2}}\vec{B}$ , where $\mu_{5}=\mu_{+}-\mu_{-}$ , and the energy shift $b_{0}$ plays the role of an $A_{0}^{5}$ -term Landsteiner (2016). In realistic systems, this is not a stable situation, since there are scattering processes that transfer electrons between the cones, so the (generalized) grand canonical equilibrium amounts to having both Weyl cones filled up to the same Fermi energy rather than to the separate ones. This amounts to $\mu_{5}=A_{0}^{5}$ and thus no CME, so we are left with $$ S_{\text{top}}=\frac{1}{8\pi^{2}}\int d^{4}x\left(\vec{b}\cdot\vec{A}\times% \vec{E}-A_{0}\vec{b}\cdot\vec{B}\right), \tag{90} $$ where $\vec{b}$ , as promised, is the (constant) axial vector that gives rise to the $\vec{b}·\vec{A}×\vec{E}$ term in $S_{\text{eff}}$ . Note that there is an additional term $\sim\vec{b}·\vec{B}$ , whose strength is determined by $A_{0}$ . V.2 $k_{0}\vec{A}·\vec{B}$ term from a flux biased Weyl supercondutor It is more difficult to generate an $S_{\text{top}}[\vec{E},\vec{B}]$ with an $\vec{A}·\vec{B}$ term by this mechanism. Formally, in the static limit, such a term corresponds to an imbalance between Weyl nodes of positive and negative chirality, which is disallowed in a system with bounded energy bands according to fermion doubling theorems. Here, we further want to clarify an important difference to the work by Zhong et. al. in Ref. Zhong et al. (2016), where the current associated to the $\vec{A}·\vec{B}$ term remains finite. This is due to the presence of an electromagnetic wave, resulting in a non-trivial electric field and this a driving in the system. Here, we consider the static case where we assume a constant magnetic field and no electric fields. Imbalances between chiral nodes are, however, known to arise in several contexts including Floquet systems Sun et al. (2018); Higashikawa et al. (2019); Zhu et al. (2020), and situations where effects that only gap out nodes of one particular chirality are present. Examples of the latter include chirality locking charge density waves Li et al. and certain Weyl superconductors O’Brien et al. (2017); Pacholski et al. (2020); Madsen et al. (2021). To illustrate that terms on the form $\vec{A}·\vec{B}$ indeed do arise in physically realizable systems, we will show it for a single Weyl node Hamiltonian originating from a flux-biased Weyl superconductor. For the sake of completeness, we first outline how the considered system is set up, referring the reader to the very insightful original work in Ref. O’Brien et al. (2017) for further details. The parent Hamiltonian is taken as, $$ \displaystyle\mathcal{H} \displaystyle=\sum_{\vec{k}}\Psi^{\dagger}_{\vec{k}}H(\vec{k})\Psi_{\vec{k}},% \quad\Psi_{\vec{k}}=\left(\psi_{\vec{k}},\sigma^{y}\psi^{\dagger}_{-\vec{k}}% \right), \displaystyle H(\vec{k}) \displaystyle=\begin{pmatrix}H_{0}(\vec{k}-e\vec{A})&\Delta_{0}\\ \Delta_{0}^{*}&-\sigma^{y}H_{0}^{*}(-\vec{k}-e\vec{A})\sigma^{y}\end{pmatrix}, \displaystyle H_{0}(\vec{k}) \displaystyle=\sum_{i}\tau^{z}\sigma^{i}\sin k_{i}+\tau^{0}\left(\beta\sigma^{% z}-\mu\sigma^{0}\right)+m_{\vec{k}}\tau^{x}\sigma^{0}, \displaystyle m_{\vec{k}} \displaystyle=m_{0}+\sum_{i}(1-\cos k_{i}). \tag{91} $$ $\tau_{i}$ and $\sigma_{i}$ , $i=x,y,z$ , are orbital and spin Pauli matrices, respectively, $\beta$ a magnetization, $\mu$ a chemical potential, $\vec{A}$ the electromagnetic vector potential, and $\Delta_{0}$ the BCS-pairing potential. A system with this Hamiltonian can be obtained by stacking alternate layers of topological insulators and conventional BCS superconductors, which introduces a coupling between the Weyl nodes centered at $(0,0,±\sqrt{\beta^{2}-m_{0}^{2}})$ of $H_{0}$ and their corresponding particle-hole conjugates. The number of ungapped particle-hole conjugate Weyl cones determine the topological phase of the Weyl superconductor, phases that can be accessed in an externally controllable way, as explained in Ref. O’Brien et al. (2017). This is done by coupling a flux-bias circuit to the material slab. This will alter the Hamiltonian, as the flux bias will be taken into account as a constant shift in the vector potential $\vec{A}$ . In this particular setup of, the flux bias gives a contribution $\Lambda/e$ to $A_{z}$ . As a result, the Weyl nodes appear at $(0,0,b_{±})$ and $(0,0,-b_{±})$ , with $$ b_{\pm}^{2}=\left(\sqrt{\beta^{2}-m_{0}^{2}}\pm\Lambda\right)^{2}-\Delta_{0}^{% 2}, \tag{97} $$ meaning that when $$ \left|\sqrt{\beta^{2}-m_{0}^{2}}-\Lambda\right|<\Delta_{0}<\sqrt{\beta^{2}-m_{% 0}^{2}}+\Lambda, \tag{98} $$ one of the two pairs of particle-hole conjugate Weyl nodes are gapped out, leaving only the nodes with positive chirality. In the regime where Weyl nodes of only one chirality are gapped out, we are left with two nodes of the same chirality. The contribution from the current from these respective nodes will then add up instead of canceling one another (as for nodes of opposite chirality). The two nodes are described by essentially the same Hamiltonian, $$ \tilde{\mathcal{H}}_{\alpha}=\sum_{\vec{k}}\tilde{\psi}^{\dagger}_{\vec{k}}% \left[\sum_{i}\nu_{i}(\delta k_{i}-Q_{i}A_{i})\sigma^{i}-Q_{0}\mu\sigma^{0}% \right]\tilde{\psi}_{\vec{k}}, \tag{99} $$ where, $\vec{k}=(0,0,b_{\alpha})+\delta\vec{k}$ , $\vec{\nu}=(1,1,-\kappa)$ , $Q_{0}=\kappa$ , ${\vec{Q}}=e(\kappa,\kappa,1/\kappa)$ , and $$ \displaystyle\mbox{$$}\kappa\approx\sqrt{1-\frac{\Delta_{0}^{2}}{(\beta+% \Lambda)^{2}}}. \tag{100} $$ Using standard procedure, this can be recast as a Lagrangian, $$ \displaystyle\mbox{$$}\mathcal{L}_{\alpha}=\bar{\psi}_{\alpha}(i\tilde{\not{% \partial}}+\tilde{\not{A}}_{\alpha})\psi_{\alpha}, \tag{101} $$ where $\tilde{\not{∂}}=\gamma^{0}∂_{0}-\nu_{i}\gamma^{i}∂_{i}$ and the left-handed chiral gauge field is $\tilde{A}_{\alpha}=(A_{0}-Q_{0}\mu,\,\nu_{i}Q_{i}A_{i}-\nu_{i}b_{\alpha;i})$ . We now set $A_{0}=0$ , and take $\vec{b}_{\alpha}$ to be constant, but allow for an $\vec{x}$ -dependent chemical potential $\mu(\vec{x})$ . The left-handed chiral anomaly is $\tilde{∂}_{\mu}J^{\mu}=-\frac{e^{2}}{32\pi^{2}}\epsilon^{\mu\nu\sigma% \lambda}\tilde{F}_{\mu\nu}\tilde{F}_{\sigma\lambda}$ , where $\tilde{F}$ is the field strength related to $\tilde{A}$ . From this, we can extract the topological part of the chiral current: $$ \displaystyle\mbox{$$}J^{i}=-\frac{e^{2}}{16\pi^{2}}\epsilon^{i0jk}\tilde{A}_{% 0}\tilde{\partial}_{i}\tilde{A}_{k}\,. \tag{102} $$ where $\tilde{A}_{0}=-Q_{0}\mu$ . This expression is well-defined and finite also for constant $\mu$ , but a direct calculation of $J^{i}$ will give a logarithmically divergent result, as shown in Appendix A. This is no contradiction with (102), since for constant $\mu$ $∂_{i}J^{i}=0$ , there can be an extra contribution that is not determined by the anomaly. We believe that the limiting procedure $∂_{i}\mu→ 0$ , which parallels the derivation in Ref. O’Brien et al. (2017), gives the correct result. It is moreover consistent with a physically motivated subtraction procedure, which is explained in Appendix A. We now express the components of the current in terms of the original fields $\vec{A}$ and $\vec{B}$ . Recalling that we have two nodes, and that $J\equiv J_{em}=2J_{L}$ , we get (for details, see Appendix A) $$ \displaystyle J^{x} \displaystyle=\frac{\kappa e^{2}\mu}{4\pi^{2}}\left[B^{x}+\left(1-\kappa^{2}% \right)\partial_{z}A_{y}\right], \displaystyle J^{y} \displaystyle=\frac{\kappa e^{2}\mu}{4\pi^{2}}\left[B^{y}-\left(1-\kappa^{2}% \right)\partial_{z}A_{x}\right], \displaystyle J^{z} \displaystyle=\frac{\kappa e^{2}\mu}{4\pi^{2}}B^{z}\,. \tag{103} $$ Integrating Eqs. (103) we get the topological action, $$ S_{\text{top}}[A]=-\frac{\kappa e^{2}\mu}{8\pi^{2}}\int d^{3}x\,\left[\vec{A}% \cdot\vec{B}+2\left(1-\kappa^{2}\right)A_{x}\partial_{z}A_{y}\right], \tag{104} $$ which is the central result of this section. The electromagnetic response in this regime can be expected to support a CME, i.e., a current in the direction on the externally applied magnetic field, which is often read off directly from the anomaly equation. In our approach the dynamics is determined by considering the full action, as in Sect. II, where we do find qualitative effects of that kind. In prospective experiments it would be natural to take the field perpendicular to the stacked planes, and to require that the slabs should be thinner than the penetration depth, so averaging over layers is justified. The term $\sim A_{x}∂_{z}A_{y}$ is not gauge invariant, and contrary to $\sim\vec{A}·\vec{B}$ , not even in the bulk. Since we started from a gauge invariant theory, this means that in a full calculation the gauge invariance would be restored by extra terms that would effectively amount to the substitution $2e\vec{A}→ 2e\vec{A}+\vec{∇}\phi$ , where $\phi$ is the phase of the superconducting order parameter. Using Eq. (100), the amplitude of the term $\sim A_{x}∂_{z}A_{y}$ is $\sim\frac{\Delta_{0}^{2}}{(\beta+\Lambda)}$ and could be made small in certain parameter ranges. In Appendix A, we shall give an alternative derivation of Eq. (104), which is more in line with the derivation in this paper and does not rely on the chiral anomaly. We note that the expression for $J^{z}$ is the same as derived in Ref. O’Brien et al. (2017) and that in Ref. O’Brien et al. (2017), it is shown that this current component is accompanied by a surface counterflow sourced by Fermi arcs. This is to ensure that the system does not have a thermal current in the ground state, as required by thermodynamics. Although we have not done an explicit calculation, we believe that the same is true for the other components in Eqs. (103). We should also note that there is a hidden assumption in the preceding derivation, in that we neglected the possibility of adding a Wess-Zumino counter term, but took the perturbative result at face value. Usually, this ambiguity is fixed by requiring gauge invariance, but in our superconducting context it is less clear. Still, it is reassuring that the simple argument based on the anomaly gives the same result as the direct calculation in Appendix A if we there make a physically motivated subtraction inspired by the treatment in Ref. O’Brien et al. (2017). Importantly, it was shown there that the expression for $J^{z}$ (which was the only one considered in that paper) agrees numerically with a direct calculation in the parent eight-band theory, Eq. (91). Although we believe that the connection to the anomaly cannot be a coincidence, we presently lack a sound theoretical argument excluding any Wess-Zumino term. VI Summary and Outlook We have motivated the consideration of emergent Chern-Simons interactions in 3+1 dimensions, displayed some of their striking phenomenological consequences, and indicated how they might be realized in plausible material systems. These interactions arise in two forms: $\vec{A}·\vec{B}$ and $\hat{n}·\vec{A}×\vec{E}-A_{0}\hat{n}·\vec{B}$ . Highlights for the first type include a precise form of current-field mixing, non-dissipative complex penetration depths, and anapole moments. This type generically arises in $s$ -wave superconductors that break parity symmetry. It is closely related, at a mathematical level, to optical activity. We suggest that it can be achieved in chirally purified organic superconductors or, more generally, by chiral doping or through parity-violating crystalline structures. We also calculated its appearance in a microscopic model based on superconducting Weyl semimetals, where it arises through a mechanism closely related to the chiral anomaly of quantum field theory. Highlights for the second type include massless boundary excitations and unusual bulk effects in electromagnetic wave propagation. This type appears to be comparatively easy to achieve in the Weyl semimetal context. We also extended the optical activity analogy in a slightly different direction—formally, towards higher rather than lower orders of gradient—to define actively chiral magnets. These do not bring in Chern-Simons terms, but physical intuition and mathematical techniques carry over. This extension frees us of the constraints of superconducitivity (notably, cryogenic temperatures and magnetic screening) and opens up many possibilities for realization in organic magnetism and metamaterials, as well as naturally occurring materials. The next, crucial development for this work will be to bring its mathematical paradise down to earth in concrete material realizations. Acknowledgement: M.S. and T.H.H. thanks J. Hannukainen and J.H. Bardarson for insightful discussions. MS acknowledges fruitful discussions with E.J. Bergholtz at an early stage of this project. F.W. is supported by the U.S. Department of Energy under Grant Contract No. DE-SC0012567, by the European Research Council under Grant No. 742104, and by the Swedish Research Council under Contract No. 335-2014-7424. Appendix A ALTERNATIVE DERIVATION OF EQS. (103) AND (104) <details> <summary>2309.10025v2/extracted/5431462/Bubble_Diagram_Sum.jpg Details</summary>  ### Visual Description The image depicts a Feynman diagram illustrating two possible photon emission processes involving a muon (μ) and a photon (A_μ). Below is a detailed technical breakdown: --- ### **Components and Flow** 1. **Left Diagram**: - A muon (μ) is represented by a solid black circle. - A photon (A_μ) is emitted from the muon, depicted as a wavy line labeled with a momentum vector **q** (→q). - The photon propagates to the right, interacting with another photon (A_ν) via a wavy line labeled **q** (→q). - The muon’s initial momentum is labeled **k** (→k), and its final momentum after photon emission is labeled **k - q** (→k - →q). 2. **Right Diagram**: - A photon (A_μ) is emitted from the muon, depicted as a wavy line labeled **q** (→q). - The photon propagates to the right, interacting with another photon (A_ν) via a wavy line labeled **q** (→q). - The muon’s initial momentum is labeled **k** (→k), and its final momentum after photon emission is labeled **k - q** (→k - →q). 3. **Summation**: - The two diagrams are connected by a "+" symbol, indicating the total amplitude is the sum of both processes. --- ### **Key Labels and Notation** - **μ**: Muon particle (solid black circle). - **A_μ, A_ν**: Photon fields (wavy lines). - **→k**: Initial momentum of the muon. - **→q**: Momentum of the emitted photon. - **→k - →q**: Final momentum of the muon after photon emission. --- ### **Interpretation** The diagram represents quantum electrodynamics (QED) interactions where a muon emits a photon (A_μ) and scatters off another photon (A_ν). The two diagrams account for different photon emission paths: 1. Photon emission before scattering (left diagram). 2. Photon emission after scattering (right diagram). The momentum conservation is enforced by the labels **k**, **q**, and **k - q**, ensuring the total momentum is preserved in each process. --- This description captures all textual and symbolic information from the diagram, omitting no labels or structural details. </details> FIG. 5: Feynman diagrams for the polarization tensor to leading order in $\mu$ . In this appendix we recalculate Eq. (104) using straight-forward diagrammatic perturbation theory. To get the pertinent static electromagnetic response function to quadratic order, we evaluate the Feynman diagram in Fig. 5 to linear order in $\mu$ and $\vec{q}$ . The relevant integrand is, $$ \displaystyle\mbox{$$}\operatorname{Tr}\left[\sigma^{i}\frac{1}{\not{k}}\frac{% 1}{\not{k}}\sigma^{j}\frac{1}{\not{k}-\not{q}}+\sigma^{i}\frac{1}{\not{k}}% \sigma^{j}\frac{1}{\not{k}-\not{q}}\frac{1}{\not{k}-\not{q}}\right]=\frac{1}{k% ^{2}}\operatorname{Tr}\left[\sigma^{i}\sigma^{j}(\not{k}-\not{q})+\sigma^{i}% \not{k}\,\sigma^{j}\right]\frac{1}{(k-q)^{2}}=-\frac{2i}{k^{4}}\epsilon^{ijk}q% _{k}+O(q), \tag{105} $$ where we used the four-vector notation $k_{\mu}=(\omega,\vec{k})$ and the cyclic property of the trace. Restoring $\mu$ , the energy and momentum integrals, and the minus sign due to the fermion loop, we get the polarization tensor: $$ \Pi^{ij}=2i\mu q_{k}\epsilon^{ijk}\int\frac{d^{3}k}{(2\pi)^{3}}\int\frac{d% \omega}{2\pi}\frac{1}{(\omega^{2}+k^{2})^{2}}. \tag{106} $$ The appearance of $\Pi^{ij}$ makes it clear that the contribution from two nodes of positive chirality symmetrically shifted from the origin of momentum space will indeed add up instead of cancel. This can be seen by shifting $k→ k± b$ , and make an expansion for small $b$ . The integral in Eq. (106) is logarithmically divergent both in the infrared and the ultraviolet. The infrared divergence is clearly a result of expanding to order $\mu$ , and is regulated if the full $\mu$ -dependence is kept. However, to follow as closely as possible to Ref. O’Brien et al. (2017), we shall instead regulate the infrared by a finite temperature $T=1/k_{B}\beta$ . The ultraviolet divergence is a consequence of that, naively, there is a contribution to the current from the whole Dirac sea, as will be discussed below. Using the standard Euclidean formulation of finite temperature QFT, the polarization tensor Eq. (106) at temperature $T$ , becomes, $$ \Pi^{ij}=2i\mu q_{k}\epsilon^{ijk}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{\beta}% \sum_{n}\frac{1}{(\omega^{2}_{n}+k^{2})^{2}}, \tag{107} $$ where $T=\frac{1}{k_{\text{B}}\beta}$ , $k_{\text{B}}$ is the Boltzmann constant, and $\omega_{n}=\frac{2\pi}{\beta}\left(n+\frac{1}{2}\right)$ are the fermionic Matsubara poles. Rewriting the sum in Eq. (107) as $$ \frac{1}{\beta}\sum_{n}\frac{1}{(\omega_{n}^{2}+k^{2})^{2}}=\frac{1}{\beta}% \left(-\frac{1}{2k}\right)\frac{\partial}{\partial k}\sum_{n}\frac{1}{\omega^{% 2}_{n}+k^{2}}=\frac{1}{\beta}\left(-\frac{1}{2k}\right)\frac{\partial}{% \partial k}\frac{\beta^{2}}{(2\pi)^{2}}\sum_{n}\frac{1}{(n+\frac{1}{2})^{2}+% \left(\frac{\beta k}{2\pi}\right)^{2}}. \tag{108} $$ and using $\sum_{n=-∞}^{∞}\frac{1}{(n+\frac{1}{2})^{2}+(Ax)^{2}}=\frac{\pi}{Ax}% \tanh(A\pi x)$ , the finite $T$ polarization tensor becomes, $$ \frac{1}{\beta}\sum_{n}\frac{1}{(\omega_{n}^{2}+k^{2})^{2}}=\left(-\frac{1}{2k% }\right)\frac{\partial}{\partial k}\frac{1}{2k}\tanh\left(\frac{\beta k}{2}% \right)=\frac{\text{sech}^{2}\left(\frac{\beta k}{2}\right)\left[\sinh\left(% \beta k\right)-\beta k\right]}{8k^{3}}. \tag{109} $$ Since the integrand is isotropic, we use spherical coordinates and $∈t\frac{d^{3}k}{(2\pi)^{3}}=∈t\frac{dk}{2\pi^{2}}k^{2}$ , and after inserting the Jacobian factor relating $\tilde{A}$ to $A$ in the integral measure, we get $$ \Pi^{ij}=2iQ_{0}\mu q_{k}\epsilon^{ijk}\frac{1}{8\pi^{2}}\int_{0}^{\infty}% \frac{dk}{|\nu_{x}\nu_{y}\nu_{z}|}\left[\frac{2\tanh\left(\frac{\beta k}{2}% \right)}{k}-\frac{\beta}{2}\text{sech}^{2}\left(\frac{\beta k}{2}\right)\right]. \tag{110} $$ The first term Eq. (110) is divergent, while the second is convergent. Subtracting the divergent piece, we arrive at the final result for the polarization tensor, $$ \Pi^{ij}=\frac{2iQ_{0}\mu\nu_{k}q_{k}\epsilon^{ijk}}{|\nu_{x}\nu_{y}\nu_{z}|}% \left(-\frac{1}{8\pi^{2}}+\text{Div}\right). \tag{111} $$ This subtraction will be discussed below. The finite part of the effective action, after substituting $q_{k}→-i∂_{k}$ , becomes $$ S_{eff}[A]=\frac{1}{2}\nu_{i}Q_{i}A_{i}\Pi^{ij}\nu_{j}Q_{j}A_{j}=-\frac{Q_{0}% \mu\nu_{i}Q_{i}\nu_{j}Q_{j}\nu_{k}}{8\pi^{2}|\nu_{x}\nu_{y}\nu_{z}|}A_{i}% \partial_{k}A_{j}\epsilon^{ijk}=-\frac{1}{8\pi^{2}}\operatorname{sgn}{(\nu_{x}% \nu_{y}\nu_{z})}Q_{0}\mu Q_{i}Q_{j}A_{i}\partial_{k}A_{j}\epsilon^{ijk}\,, \tag{112} $$ and, finally, the current $$ \displaystyle J^{l} \displaystyle=\frac{\delta S_{eff}[A]}{\delta A_{l}} \displaystyle=-\operatorname{sgn}{(\nu_{x}\nu_{y}\nu_{z})}\frac{\mu Q_{i}Q_{j}% }{8\pi^{2}}\partial_{k}A_{j}\epsilon^{ijk}\frac{\delta A_{i}}{\delta A_{l}} \displaystyle=-\operatorname{sgn}{(\nu_{x}\nu_{y}\nu_{z})}\frac{\mu Q_{l}Q_{j}% }{8\pi^{2}}\partial_{k}A_{j}\epsilon^{ljk}, \tag{113} $$ where the index $l$ in the right hand side is not summed over. Inserting the system parameters from the Hamiltonian in Eq. (99) gives $$ \displaystyle S_{\text{top}}[A] \displaystyle=-\frac{\mu\kappa}{8\pi^{2}}Q_{i}Q_{k}A_{i}\partial_{j}A_{k}% \epsilon^{ijk} \displaystyle=-\frac{\mu\kappa}{8\pi^{2}}\left[Q_{x}A_{x}\left(Q_{z}\partial_{% y}A_{z}-Q_{y}\partial_{z}A_{y}\right)+Q_{y}A_{y}\left(Q_{x}\partial_{z}A_{x}-Q% _{z}\partial_{x}A_{z}\right)+Q_{z}A_{z}\left(Q_{y}\partial_{x}A_{y}-Q_{z}% \partial_{y}A_{x}\right)\right] \displaystyle=-\frac{\mu e^{2}\kappa}{8\pi^{2}}\left[A_{x}\left(\partial_{y}A_% {z}-\kappa^{2}\partial_{z}A_{y}\right)+A_{y}\left(\kappa^{2}\partial_{z}A_{x}-% \partial_{x}A_{z}\right)+A_{z}\left(\partial_{x}A_{y}-\partial_{y}A_{x}\right)\right] \displaystyle=-\frac{\mu e^{2}\kappa}{8\pi^{2}}\left\{A_{x}\left[\partial_{y}A% _{z}-\partial_{z}A_{y}+\left(1-\kappa^{2}\right)\partial_{z}A_{y}\right]+A_{y}% \left[\partial_{z}A_{x}-\partial_{x}A_{z}-\left(1-\kappa^{2}\right)\partial_{z% }A_{x}\right]+A_{z}\left(\partial_{x}A_{y}-\partial_{y}A_{x}\right)\right\} \displaystyle=-\frac{\mu e^{2}\kappa}{8\pi^{2}}\left[-\vec{A}\cdot\vec{B}+% \left(1-\kappa^{2}\right)\left(A_{y}\partial_{z}A_{x}-A_{x}\partial_{z}A_{y}% \right)\right] \displaystyle=\frac{\mu e^{2}\kappa}{8\pi^{2}}\left[\vec{A}\cdot\vec{B}+\left(% 1-\kappa^{2}\right)\left(A_{y}\partial_{z}A_{x}-A_{x}\partial_{z}A_{y}\right)% \right]\,, \tag{114} $$ and finally restoring $\hbar$ , using the notation $e^{*}=\kappa e$ , the current becomes $$ J^{l}=\frac{\delta S_{\text{top}}[A]}{\delta A_{l}}=\frac{\mu ee^{*}}{h^{2}}% \left[B^{l}+\left(1-\kappa^{2}\right)\left(\delta^{l}_{y}\partial_{z}A_{x}-% \delta^{l}_{x}\partial_{z}A_{y}\right)\right]. \tag{115} $$ Using the gauge where $\vec{A}=\left(0,Bx,\Lambda/e\right)$ , the only surviving component of the current reads, $$ J^{z}=\frac{ee^{*}\mu}{h^{2}}B_{z}. \tag{116} $$ which agrees with the result of Ref. O’Brien et al. (2017). We now return to discuss the subtraction of the logarithmic UV divergence in the integral Eq. (110). Note that the second, convergent, term has support only close to the Fermi surface (i.e., at $k≈ 0$ ) while the first, UV divergent part gets contributions from the full Fermi sea. Since anomalies are expressed in IR phenomena, it is plausible to subtract the first term and keep the second. This was done in Ref. O’Brien et al. (2017), where it was shown that it gives results consistent with a numerical simulation of the full Hamiltonian Eq. (91). It is reassuring that these various methods give the same result for the current, and satisfying that our derivations make a close connection to the chiral anomaly. Appendix B DETAILS OF SPHERE AND CYLINDER SOLUTIONS B.1 Sphere Following Eqs. (26) and (32), we define $\alpha$ and $\sqrt{\alpha}$ : $$ \displaystyle\alpha_{1}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \gamma\,+\,\frac{\beta}{2}\bigl{(}-\beta\,+% \,i\sqrt{4\gamma-\beta^{2}}\,\bigr{)},\quad \displaystyle\sqrt{\alpha}_{1}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \sqrt{\gamma-\frac{\beta^{2}}{4}}+i\frac{% \beta}{2}=p+iq \tag{117} $$ The reference solution has the form $$ \displaystyle B_{r} \displaystyle= \displaystyle\frac{2}{\alpha r^{3}}\left[\sinh\left(\sqrt{\alpha}r\right)-% \sqrt{\alpha}r\cosh\left(\sqrt{\alpha}r\right)\right]\cos\theta, \displaystyle B_{\theta} \displaystyle= \displaystyle\frac{1}{\alpha r^{3}}\left[\left(1+\alpha r^{2}\right)\sinh\left% (\sqrt{\alpha}r\right)-\sqrt{\alpha}\ r\cosh\left(\sqrt{\alpha}r\right)\right]% \sin\theta, \displaystyle B_{\phi} \displaystyle= \displaystyle\frac{1}{r^{2}}\left[\sinh\left(\sqrt{\alpha}r\right)-\sqrt{% \alpha}r\cosh\left(\sqrt{\alpha}r\right)\right]\sin\theta \tag{119} $$ and since we must take $\alpha$ and $\sqrt{\alpha}$ complex, we get complex fields, whose real and imaginary parts satisfy our Eq. (19) separately, that we must combine to satisfy the boundary conditions. The real part $B_{r}$ , after considerable algebra, reads $$ \displaystyle B^{r}_{r}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\cos\theta}{\gamma^{2}r^{3}} \displaystyle\times \displaystyle\sqrt{4\gamma-\beta^{2}}\cosh[\frac{1}{2}r\sqrt{4\gamma-\beta^{2}% }]\left(\beta\sin\frac{\beta r}{2}-\gamma r\cos\frac{\beta r}{2}\right)-\sinh[% \frac{1}{2}r\sqrt{4\gamma-\beta^{2}}]\left(\left(\beta^{2}-2\gamma\right)\cos% \frac{\beta r}{2}+\beta\gamma r\sin\frac{\beta r}{2}\right) \displaystyle B^{r}_{\theta}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\sin\theta}{2\gamma^{2}r^{3}} \displaystyle\times \displaystyle\sqrt{4\gamma-\beta^{2}}\cosh[\frac{1}{2}r\sqrt{4\gamma-\beta^{2}% }]\left(\beta\sin\frac{\beta r}{2}-\gamma r\cos\frac{\beta r}{2}\right)-\sinh[% \frac{1}{2}r\sqrt{4\gamma-\beta^{2}}]\left(\left(\beta^{2}-2\gamma\left(\gamma r% ^{2}+1\right)\right)\cos\frac{\beta r}{2}+\beta\gamma r\sin\frac{\beta r}{2}\right) \displaystyle B^{r}_{\phi}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\sin\theta}{2\gamma r^{2}} \displaystyle\times \displaystyle\sinh[\frac{1}{2}r\sqrt{4\gamma-\beta^{2}}]\left(2\gamma r\sin% \frac{\beta r}{2}+\beta\cos\frac{\beta r}{2}\right)-\sqrt{4\gamma-\beta^{2}}% \sin\frac{\beta r}{2}\cosh[\frac{1}{2}r\sqrt{4\gamma-\beta^{2}}] \tag{122} $$ and the imaginary part $B^{i}$ reads $$ \displaystyle B^{i}_{r}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\cos\theta}{\gamma^{2}r^{3}} \displaystyle\times \displaystyle\cosh\left[\frac{1}{2}r\sqrt{4\gamma-\beta^{2}}\right]\left(\beta% \gamma r\cos\frac{\beta r}{2}-\left(\beta^{2}-2\gamma\right)\sin\frac{\beta r}% {2}\right)-\sqrt{4\gamma-\beta^{2}}\sinh\left[\frac{1}{2}r\sqrt{4\gamma-\beta^% {2}}\right]\left(\gamma r\sin\frac{\beta r}{2}+\beta\cos\frac{\beta r}{2}\right) \displaystyle B^{i}_{\theta}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\sin\theta}{2\gamma^{2}r^{3}} \displaystyle\times \displaystyle\cosh\left[\frac{1}{2}r\sqrt{4\gamma-\beta^{2}}\right]\left(2% \gamma\bigl{(}\left(\gamma r^{2}+1\right)-\beta^{2}\right)\sin\frac{\beta r}{2% }+\beta\gamma r\cos\frac{\beta r}{2}\bigr{)} \displaystyle- \displaystyle\sqrt{4\gamma-\beta^{2}}\sinh\left[\frac{1}{2}r\sqrt{4\gamma-% \beta^{2}}\right]\left(\gamma r\sin\frac{\beta r}{2}+\beta\cos\frac{\beta r}{2% }\right) \displaystyle B^{i}_{\phi}\leavevmode\nobreak\ \displaystyle= \displaystyle\leavevmode\nobreak\ \frac{\sin\theta}{2\gamma r^{2}} \displaystyle\times \displaystyle\sqrt{4\gamma-\beta^{2}}\cos\frac{\beta r}{2}\sinh\left[\frac{1}{% 2}r\sqrt{4\gamma-\beta^{2}}\right]+\cosh\left[\frac{1}{2}r\sqrt{4\gamma-\beta^% {2}}\right]\left(\beta\sin\frac{\beta r}{2}-2\gamma r\cos\frac{\beta r}{2}\right) \tag{123} $$ We note that this is the solution inside the superconducting sphere, and that it is thus valid only if $r≤ R$ . To satisfy the boundary condition $j_{r}(R)=0$ , we must take a linear superposition $B^{r}+\eta B^{i}$ to get the full solution inside the superconducting sphere. One finds $$ \displaystyle\eta=\frac{\sinh[\frac{1}{2}R\sqrt{4\gamma-\beta^{2}}]\left(2% \gamma R\sin\frac{\beta R}{2}+\beta\cos\frac{\beta R}{2}\right)-\sqrt{4\gamma-% \beta^{2}}\sin\frac{\beta R}{2}\cosh[\frac{1}{2}R\sqrt{4\gamma-\beta^{2}}]}{% \cosh[\frac{1}{2}R\sqrt{4\gamma-\beta^{2}}]\left(2\gamma R\cos\frac{\beta R}{2% }-\beta\sin\frac{\beta R}{2}\right)-\sqrt{4\gamma-\beta^{2}}\cos\frac{\beta R}% {2}\sinh[\frac{1}{2}R\sqrt{4\gamma-\beta^{2}}]}. \tag{124} $$ Thus, the final solution for the magnetic field inside the superconducting sphere reads $$ \vec{B}^{\text{in}}_{\text{sphere}}=\left(B_{r}^{r},B^{r}_{\theta},B^{r}_{\phi% }\right)+\eta\left(B^{i}_{r},B^{i}_{\theta},B^{i}_{\phi}\right) \tag{125} $$ Outside the sphere, i.e., for $r>R$ , the magnetic field takes the usual dipole expansion form according to the solution of London, and reads London (1950), $$ \displaystyle B^{\text{out}}_{r}(r,\theta) \displaystyle=\left(H_{0}+\frac{2M}{r^{3}}\right)\cos{\theta}, \displaystyle B^{\text{out}}_{\theta}(r,\theta) \displaystyle=\left(-H_{0}+\frac{M}{r^{3}}\right)\sin{\theta}, \tag{126} $$ with $$ \displaystyle M \displaystyle=\frac{R^{3}}{3}\left[\frac{B^{\text{in}}_{r}(R,\theta)}{\cos{% \theta}}+\frac{B^{\text{in}}_{\theta}(R,\theta)}{\sin{\theta}}\right], \displaystyle H_{0} \displaystyle=\frac{1}{3}\left[\frac{B^{\text{in}}_{r}(R,\theta)}{\cos{\theta}% }-2\frac{B^{\text{in}}_{\theta}(R,\theta)}{\sin{\theta}}\right]. \tag{128} $$ From the internal solutions, one can calculate the corresponding currents $\vec{j}_{\text{sphere}}=\vec{∇}×\vec{B}^{\text{in}}_{\text{sphere}}$ , which can be used to explicitly calculate the magnetic anapole moment, which is given by $$ T_{i}=\frac{1}{10c}\int\left[r_{i}\left(\vec{r}\cdot\vec{j}\right)-2r^{2}J_{i}% \right]d^{3}x. \tag{130} $$ In the present situation, $T_{x}$ and $T_{y}$ are both zero, but $T_{z}$ takes a finite value. On the unit sphere, and in units where $c=1$ , it explicitly reads, $$ T_{z}=-2\pi\frac{-6\beta\left(1+\gamma\right)\cosh{\left(2\delta^{-1}\right)}+% 2\left[-3\beta\left(\gamma-1\right)\cos{\beta}+\left(3\beta^{2}-\gamma^{2}% \right)\sin{\beta}\right]+\beta\delta\left[-3\beta^{2}+\gamma\left(12+\gamma% \right)\right]\sinh{\left(2\delta^{-1}\right)}}{3\gamma^{2}\left[\delta\cosh{% \left(\delta^{-1}\right)}\left(2\gamma\cos{\frac{\beta}{2}}-\beta\sin{\frac{% \beta}{2}}\right)-2\cos{\frac{\beta}{2}}\sinh{\left(\delta^{-1}\right)}\right]}, \tag{131} $$ where $\delta=\frac{2}{\sqrt{4\gamma-\beta^{2}}}$ is the penetration depth. B.2 Cylinder In the case of a cylinder, we take real and imaginary parts of an ansatz consisting of modified Bessel functions of the second kind: $$ B_{r}=0,\quad B_{\phi}=\sqrt{\alpha}\kappa K_{1}(r\sqrt{\alpha}),\quad B_{z}=K% _{0}(r\sqrt{\alpha}) $$ . As in the spherical case, we take the real and imaginary parts independently and find coefficients to satisfy boundary conditions $B_{z}=B_{0},\quad B_{\phi}=0$ , fixing $R=1$ . Since this brings in Bessel functions, it must be done numerically for specific values $\beta$ and $\gamma$ . We looked at a few cases to get numerical solutions and here are the answers that we get for the azimuthal components of the magnetic field. For $\beta=15,\gamma=57$ , we find $$ \displaystyle(-4.77814B_{0})\,{\rm Re}\left\{K_{1}\left[\left(\frac{15i}{2}+% \frac{\sqrt{3}}{2}\right)r\right]\right\} \displaystyle+ \displaystyle(-2.14645B_{0})\,{\rm Im\/}\left\{K_{1}\left[\left(\frac{15i}{2}+% \frac{\sqrt{3}}{2}\right)r\right]\right\}. \tag{132} $$ For $\beta=10,\gamma=4$ , we find $$ \displaystyle(-69.1343B_{0})\,{\rm Re}\left\{K_{1}\left[\left(5i+\sqrt{15}% \right)r\right]\right\} \displaystyle+ \displaystyle(69.2327B_{0})\,{\rm Im\/}\left\{K_{1}\left[\left(5i+\sqrt{15}% \right)r\right]\right\}. \tag{133} $$ For $\beta=1,\gamma=200$ , we find $$ \displaystyle(2.06099\times 10^{6}B_{0})\,{\rm Re}\left\{K_{1}\left[\left(% \frac{i}{2}+\frac{\sqrt{799}}{2}\right)r\right]\right\} \displaystyle+ \displaystyle(3.61166\times 10^{6}B_{0})\,{\rm Im\/}\left\{K_{1}\left[\left(% \frac{i}{2}+\frac{\sqrt{799}}{2}\right)r\right]\right\}. \tag{134} $$ References - Nazarov (1986) Y. V. Nazarov, Zh. Eksp. Teor. Fiz, 91, 332 (1986). - Stone (1991) M. Stone, Annals of Physics 207, 38 (1991). - Samokhin (2004) K. V. Samokhin, Phys. Rev. B 70, 104521 (2004). - Agterberg (2003) D. Agterberg, Physica C: Superconductivity 387, 13 (2003), proceedings of the 3rd Polish-US Workshop on Superconductivity and Magnetism of Advanced Materials. - Bauer and Sigrist (2012) E. Bauer and M. Sigrist, Non-Centrosymmetric Superconductors, Lecture Notes in Physics (Springer Berlin, Heidelberg, 2012). - Samoilenka and Babaev (2020) A. Samoilenka and E. Babaev, Phys. Rev. B 102, 184517 (2020). - Hannukainen et al. (2020) J. D. Hannukainen, Y. Ferreiros, A. Cortijo, and J. H. Bardarson, Phys. Rev. B 102, 241401(R) (2020). - Hannukainen et al. (2021) J. D. Hannukainen, A. Cortijo, J. H. Bardarson, and Y. Ferreiros, SciPost Phys. 10, 102 (2021). - Zhong et al. (2016) S. Zhong, J. E. Moore, and I. Souza, Phys. Rev. Lett. 116, 077201 (2016). - Landsteiner (2016) K. Landsteiner, Acta Phys. Pol. B 47 (2016), 10.5506/APhysPolB.47.2617. - Vazifeh and Franz (2013) M. M. Vazifeh and M. Franz, Phys. Rev. Lett. 111, 027201 (2013). - Başar et al. (2014) G. Başar, D. E. Kharzeev, and H.-U. Yee, Phys. Rev. B 89, 035142 (2014). - (13) P. Hosur and X. Qi, Comptes Rendus Physique https://doi.org/10.1016/j.crhy.2013.10.010. - Chen et al. (2013) Y. Chen, S. Wu, and A. A. Burkov, Phys. Rev. B 88, 125105 (2013). - Chang and Yang (2015) M.-C. Chang and M.-F. Yang, Phys. Rev. B 91, 115203 (2015). - Carroll et al. (1990) S. M. Carroll, G. B. Field, and R. Jackiw, Phys. Rev. D 41, 1231 (1990). - Moffatt (1969) H. K. Moffatt, J . Fluid Mech. 35, 117 (1969). - Woltjer (1958) L. Woltjer, Proceedings of the National Academy of Sciences 44, 489 (1958). - London (1950) F. London, Superfluids, Vol 1: Macroscopic Theory of Superconductivity (Dover Publications Inc., Mineola, New York, 1950). - Nair (2006) V. Nair, Quantum Field Theory: A Modern Perspective, Graduate Texts in Contemporary Physics (Springer New York, NY, 2006). - Sun et al. (2018) X.-Q. Sun, M. Xiao, T. c. v. Bzdušek, S.-C. Zhang, and S. Fan, Phys. Rev. Lett. 121, 196401 (2018). - Higashikawa et al. (2019) S. Higashikawa, M. Nakagawa, and M. Ueda, Phys. Rev. Lett. 123, 066403 (2019). - Zhu et al. (2020) Y. Zhu, T. Qin, X. Yang, G. Xianlong, and Z. Liang, Phys. Rev. Res. 2, 033045 (2020). - (24) G. Li, H. Yang, P. Jiang, C. Wang, Q. Cheng, S. Tian, G. Han, C. Shen, X. Lin, H. Lei, W. Ji, Z. Wang, and H.-J. Gao, Nature Communications , 2914. - O’Brien et al. (2017) T. E. O’Brien, C. W. J. Beenakker, and i. d. I. Adagideli, Phys. Rev. Lett. 118, 207701 (2017). - Pacholski et al. (2020) M. Pacholski, C. Beenakker, and İ. Adagideli, Annals of Physics 417, 168103 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond. - Madsen et al. (2021) K. A. Madsen, P. W. Brouwer, and M. Breitkreiz, Phys. Rev. B 104, 035109 (2021).