Axion strings are superconducting

We explore the cosmological consequences of the superconductivity of QCD axion strings. Axion strings can support a sizeable chiral electric current and charge, which alters their early universe dynamics. Shrinking axion string loops can become effectively stable remnants called vortons, supported by the electromagnetic force of the string current. Generically, vortons produced by axion strings overclose the universe, unless there are efficient current leakage processes. Furthermore, if a primordial magnetic field (PMF) exists in the early universe, a large current is induced on axion strings, creating a significant drag force from interactions with the surrounding plasma. As a result, the strings are slowed down, which leads to an orders of magnitude enhancement in the number of strings per Hubble volume. Finally, we study the implications for the QCD axion relic abundance. The QCD axion window is shifted by orders of magnitude in some parts of our parameter space.


Introduction
One of the long standing puzzles in theoretical physics is the tiny value of θ < 10 −10 1, usually referred to as the strong CP problem. The strong CP problem cannot be solved anthropically [1]. A possible solution is by the Peccei-Quinn (PQ) mechanism [2,3], in which a broken anomalous chiral U(1) symmetry dynamically sets the θ angle to zero [4]. PQ symmetry breaking produces a light pseudo-Goldstone boson called the axion [5,6]. The axion is also a dark matter (DM) candidate.
If the PQ symmetry is broken during inflation, the corresponding axion field is frozen during the inflationary era. Eventually, when the Hubble expansion H becomes smaller than the axion mass H m a , at which point the axion starts behaving like a massive particle, and the initial axion field value determines the relic axion density [7][8][9]. Quantum fluctuations of the axion are generated, and the isocurvature fluctuation is generally large. To avoid observational constraints [10], the inflationary Hubble scale H inf must be as small as H inf 10 7 GeV [11,12]. In the case that PQ symmetry breaking occurs after inflation, the Kibble-Zurek mechanism [13][14][15][16][17] leads to the formation of cosmic axion strings, topological defects associated with U(1) PQ symmetry breaking. In particular, the generated axion strings are vortices around which the phase of the scalar field rotates by 2π. Axion strings are topologically stable because π 1 (U (1) PQ ) = Z, and the U (1) PQ symmetry is restored in their core. The dynamics of axion strings is crucial for understanding the relic density of axion DM. In this paper we focus on post inflationary PQ breaking, and explore novel dynamics that governs the evolution of the resulting axion string network.
In standard cosmic string cosmology, an intially large number of axion strings appear at the PQ phase transition [15,17]. As soon as they are created, they begin to reconnect with each other and the total length of string per Hubble horizon diminishes. Due to friction caused by interactions with the thermal plasma, the number of strings in the Hubble horizon, ξ, may be larger than O(1) at first [18]. As the universe cools, plasma friction weakens and the string abundance eventually enters a scaling regime, where ξ ∼ O(1) of strings exist per Hubble volume. The scaling regime persists over many orders of magnitude in temperature until T ∼ > T QCD . The scaling of strings in the early universe has been extensively studied by numerical simulations [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. The scaling behavior is the result of a balance between string length entering the horizon, and string length which is lost by the emission of unstable string loops due to efficient string reconnection.
At T ∼ > T QCD , the axion gets a mass from QCD non-perturbative effects and consequently domain walls are formed between the axion strings. At this point, each axion string becomes the boundary of N DW domain walls, where N DW is the coefficient of the PQ-QCD anomaly. Generically, when N DW > 1 the domain walls cannot decay 1 , and so they can overclose the universe. For this reason, most post-inflationary QCD axion scenarios are KSVZ-like [38,39] with N DW = 1. In this case, the string-domain wall systems are unstable, and they annihilate into cold axions which form the bulk of the DM relic density.
In this paper, we focus on the KSVZ scenario as an example with N DW = 1 and point out an important effect -the axion string is (EM and color) superconducting [40]. 2 While our discussion is completely generic and can be applied to any axion scenario, we use the KSVZ example for a quantitative discussion. The superconductivity of axion strings is widely known in the context of Chern-Simons theory, but has been overlooked in phenomenological cosmological studies. We find that axion string superconductivity changes the density of axion strings in the early universe, and hence the relic axion abundance. In particular, the presence of a primordial magnetic field (PMF), the axion relic abundance could be considerably larger than previous expecta-tions. We also study if superconductivity allows for metastable remnants of the axion string, called vortons, as discussed previously in Ref. [41][42][43][44][45][46] for local strings supporting vector-like currents [47]. We find that these are generically unstable in our scenario, as a result of zero mode decay induced by the curvature of the vortons.

Fermion Zero-Modes and Axion String Superconductivity
As a first step, we provide a quick review of the superconductivity of KSVZ axion strings. The argument is very general, and applies to other axion models as well. This section also serves to define our notation, and to introduce features of the model that will play a role in our analysis.
In addition to the Standard Model (SM) fields, the minimal KSVZ model [38,39] has two chiral fermions ψ L and ψ R which transform as 3 under color and are weak SU (2) singlets, and a color-singlet complex scalar Φ with the interaction Lagrangian (2.1) The phase of y Φ can be absorbed into Φ so that y Φ is real. The chiral transformation is anomalous under QCD, θ → θ − α, and is the PQ symmetry in this theory. This example has N DW = 1. We assume Φ acquires a vacuum expectation value (VEV), Φ = f a / √ 2 which breaks the PQ symmetry and gives ψ a mass m ψ = y Φ f a / √ 2. The phase of Φ is the axion a, i.e. neglecting radial fluctuations

3)
f a is the axion decay constant, and a/f a has periodicity 2π. As we are considering post inflationary PQ breaking, the KSVZ fermions are abundant in the early universe, and so they must decay to SM particles before big-bang nucleosynthesis (BBN) (see [37] and references within). To allow for decays, there have to be interactions between ψ and SM particles. For example, KSVZ fermions with the quantum numbers of vertorlike quarks (VLQ) have Yukawa couplings via the Higgs field H to the SM quark doublet where y is the Yukawa coupling of the KSVZ fermion to the SM Higgs. The Yukawa term can be of either form, depending on whether ψ has weak hypercharge Y = −1/3 or Y = 2/3, i.e. the quantum numbers of the SM fields b R or t R . The interactions Eq. (2.4) preserve the PQ symmetry Eq. (2.2). Note that because the gauge group of the standard model particles is actually limited to G SM ≡ SU(3) × SU(2) × U(1) / Z 6 , ψ must have Y = 2/3 mod 1, so that it can decay to the standard model particles. The above two examples are the simplest choices that allow for direct Yukawa couplings without need for additional particles. The axion string is a topological defect around which the axion (the phase of Φ) winds by 2π. An infinitely straight axion string in the z-direction has the form where r ⊥ = x 2 + y 2 is the transverse distance from the string, and θ = tan −1 (y/x) is the azimuthal angle. The radial solution h(r ⊥ ) depends on the axion potential V (Φ), with h(r ⊥ ) → 1 as r ⊥ → ∞, and h(0) = 0 so that Φ is non-singular at the origin.
In the string background Eq. (2.5) for Φ, the Dirac operator for the KSVZ fermions in the xy-plane has a zero eigenvalue, as guaranteed by an index theorem [48]. As a result, the 3 + 1 dimensional theory has a chiral massless fermion moving in the zdirection, which is a KSVZ fermion bound to the string [49]. The transverse size of the massless fermion wavefunction is of order 1/m ψ . Because the KSVZ fermion interaction with Φ in Eq. (2.1) is chiral, the axion string has a trapped chiral fermion left (downward) moving fermion and its CPT conjugate anti-particle, originating from the KSVZ fields ψ L /ψ R . Replacing Φ → Φ * in Eq. (2.1) would result in a right (upward) moving trapped mode. Since ψ carries color and hypercharge, the zero-mode carries color and U (1) current along the string.
The 1 + 1 dimensional theory of the chiral fermion has a U (1) gauge anomaly. Below the electroweak phase transition, we simply consider the 1+1 dimensional gauge anomaly of U (1) EM , where N c = 3, e ψ is the electric charge of ψ multiplied by e, and E z is the electric field at the string core in the z direction, leading to a violation of current conservation on the string. 3 On the other hand the full 3 + 1 dimensional theory is anomaly free, and has a conserved electromagnetic current. This "paradox" was resolved by Callan and Harvey [40] by noting that there is a Goldstone-Wilczek current [50] in the bulk, leading to charge inflow onto the string exactly equal to Eq. (2.6). Here F µν is the U (1) field strength tensor, and we use the sign convention 0123 = +1. 4 The above current is deposited on the string in the form of massless fermion zero-modes. This mechanism is referred to as anomaly inflow. We have started from an explicit UV model to see the chiral anomaly on the string. However, the anomaly is determined solely by the IR Lagrangian. Let us derive the anomaly inflow from the coupling between the axion a and gauge fields. Away from the string core, the fermion ψ is massive and can be integrated out to give a low-energy Lagrangian for the interaction of the axion with the electromagnetic field Here, we omit the color gauge field for simplicity. The effective Lagrangian preserves electromagnetic gauge invariance, since the original 3 + 1 dimensional theory had no electromagnetic anomaly, and was gauge invariant. From Eq. (2.8), the electromagnetic current in the bulk is which is the same as the Goldstone-Wilczek current Eq. (2.7). The divergence of J µ a is which naively vanishes since F µν is antisymmetric. However, the axion field has a winding number around the origin, leading to [40,52,53] ( 13) and is exactly canceled by the zero-mode current divergence Eq. (2.6). In other words, the theory on the string must produce the same anomaly with the opposite sign to the inflow current to cancel the divergence of this current regardless of the UV theory, as long as the UV theory preserves electromagnetic gauge invariance. In the KSVZ model, we showed that the current is carried by a chiral fermion on the string so the trapped current density is chiral, 5 ρ = −I, (2.14) where ρ is the charge per unit length and I is the current. The relative minus sign is because the zero-mode is left-moving. Taking the divergence If the system is placed in a background electric field that is z-independent, I does not depend on z, and i.e. the current increases linearly with the applied electric field. Thus, there is no resistivity and the axion string is superconducting. A similar argument holds for the color current. 6 A microscopic description of the current is useful for phenomenological applications later in this paper. The current is proportional to the number density of zero mode fermions traveling along the string. If the zero mode fermion states are occupy energy all levels up to the Fermi energy ε F , the current on the string is since the 1-dimensional density of states is dk/(2π) = dE/(2π) for each color. The Fermi energy picture allows us to understand the evolution of the current, Eq. (2.16) 5 See Ref. [52,53] for a discussion on regularization. 6 Since the KSVZ fermion ψ is a color triple and the interaction with Φ does not depend on color, there is a zero-mode ψ 0,α for each color component α = 1, 2, 3 of ψ. ψ zero-modes can be combined to form color-singlet "baryons" αβγ ψ 0,α ψ 0,β ψ 0,γ . A background electric field is color-singlet, and will populate the string with such color singlet states. from another viewpoint. In this language, the current evolution is due to the Fermi energy changing as a result of the applied electric field, and so new fermions "appear from the vacuum" as negative energy states become positive energy states [54]. For Fermi energy ε F , the energy per unit length in bound fermions ε I = (ε F /2)(N c ε F )/(2π), so that neglecting fermion-fermion interactions. In addition to the energy of the zero modes, axion strings also have an intrinsic tension µ ln(Lf a )f 2 a , where L is a typical distance between strings. 7 For the fermion energy is negligible compared with the intrinsic string tension. Unless otherwise stated, the magnitude of the current is assumed to satisfy this condition, so that previous studies on axion string interactions, which do not include Eq. (2.18), can be used. Nevertheless, interactions between axion strings and charged particles are modified because of the trapped currents, as studied in later sections. Superconducting strings with zero-mode currents were first discussed in detail in Ref. [47]. This paper considers a vector-like UV theory with zero-modes propagating in both directions along the string. An important difference from this earlier analysis is the chiral nature of axion strings. One consequence of chirality is that electromagnetic fields can induce currents around the strings via Goldstone-Wilczek current inflow, Eq. (2.7). In a large enough box (such as the entire universe), the axion winding is trivial and the Goldstone-Wilczek current inflow from infinity vanishes. However, the Goldstone-Wilczek current can still transport charge between strings, so that string loops develop charge even though the whole system remains neutral. Charge inflow and electrodynamics in the presence of chiral axion strings were studied in Refs. [40,51,52,55], and the rotation of the polarization of electromagnetic waves propagating in the string background was studied in Refs. [52,[55][56][57].

Current Leakage
Besides the generation of zero mode currents on the string by the ambient EM field, we also have to take into account current leakage processes, which play an important role in the dynamics. In this section, we study mechanisms responsible for current leakage. These have mostly been analyzed for non-anomalous superconducting strings and are irrelevant in our case. For example, ref. [47,58] considered a vector-like theory with zero-modes moving in both directions along the string. In this case, the zero-modes can collide and produce fermion pairs that escape to infinity if the Fermi energy of the zero-modes exceeds m ψ , the mass of the KSVZ fermion at infinity. This pair-production process limits the Fermi energy on the string [47] by ε F ≤ m ψ , which translates into a maximum current I ≤ 2e ψ m ψ /(2π), where 1/(2π) is the density of states in 1 + 1 dimension. The overall factor of 2 is because a vector-like theory can have positive charges moving in one direction, and negative charges moving in the opposite direction.
In contrast, the axion string has chiral zero-modes, so zero mode particles and antiparticles are massless and move in the same direction. For this reason, they cannot annihilate among themselves and produce massive fermions which escape to infinity; the process is forbidden by energy-momentum conservation and, there is no maximum current from zero mode collisions. For straight strings this can also be seen by making a Lorentz transformation [59]. The collision between a zero mode particle and its antiparticle on the long string is boost-equivalent to a massless particle decaying into massive particles, and so is kinematically forbidden.
Another potential source of leakage is the Goldstone-Wilczek current itself. This current, induced by the ambient EM field, can change the number of the chiral zeromode on strings. However, for the long strings we the GW current and the zero mode current on the strings are in equilibrium. In particular, in the case of string loops, the Goldstone-Wilczek ccan never change the net charge of the loop but only redistribute it. This is because the overall axion winding around the loop (as seen from far away) is trivial.

Leakage from plasma scattering
Though the current on axion strings cannot dissipate by internal collisions among zero modes, it can be reduced by scattering processes involving a plasma particle incident on the string.
As in the calculation of the photo-ionization of hydrogen, here the rate for the the interaction of plasma particles with the current depends on the geometric overlap of the zero mode wavefuntion and the plane wave of the incoming plasma particle. In the hydrogen case, the ground state photo-ionization rate depends on the Fourier transform of the ground-state wavefunction, and is suppressed by (ωa 0 ) 3 at low frequencies, where ω is the incident photon frequency and a 0 is the Bohr radius [60, §37]. In our problem, the zero-modes is localized in two transverse dimensions, so the cross section for ionization is suppressed by p 2 ⊥ /m 2 ψ , the ratio of the fermion zero-mode size and the wavelength of the incident particle.
At temperatures T ∼ m ψ ∼ y Φ f a , thermal plasma particles which scatter off the string have enough kinetic energy to knock the zero mode ψ into the bulk. The cross section for such scattering process is geometric, σ ∼ f −2 a , but the process is Boltzmann suppressed for T < m ψ . Consequently, it decouples below T ∼ m ψ /b, b ∼ 20 (for f a ∼ m ψ ∼ 10 8-12 GeV), when its rate drops below the Hubble expansion rate.
At temperatures much smaller than the fermion mass, T m ψ , the current can no longer leak into ψ particles, but it can still leak by producing SM particles via the interaction Eq. (2.4), the rate of which depends on the Yukawa coupling y. Prototypical scattering processes leading to current decrease arē via t-and s-channel Q exchange, where ψ 0 string is the string zero-mode. The cross-section for these processes, which are related by crossing, is of order where g s is the QCD coupling constant. The effective coupling of the zero-mode is suppressed by the transverse wavelength of the incoming fermion relative to the size of the string core, p ⊥ /m ψ with p ⊥ ∼ T , times the coupling y from Eq. (2.4). The 1/(16πs) factor is a dimensional estimate of the phase space where s ∼ T ε F is the center-of-mass energy squared, ε F is the zero-mode Fermi energy and m ψ = y Φ f a , with y Φ from Eq. (2.1), and we have used Eq. (2.17) to rewrite ε F in I. The resulting rate to ionize a zero-mode is then where we take v ⊥ ∼ 1 because the zero mode is relativistic, and h = g B + (3/4)g F is the effective number of degrees of freedom which can interact with the zero-mode. Comparing this rate to the Hubble expansion rate, we see that the rate Eq. (3.4) where g = g B + (7/8)g F for all particle species. For temperatures lower than T scat we can neglect current dissipation by scattering off the thermal background. The last factor in Eq. (3.4) is order unity, and will be neglected.

Leakage from string oscillations
Another potential leakage mechanism is the interaction of fermion zero-modes with the oscillation modes of the string. We can compute the decay rate by quantizing the oscillations, which are then treated as real scalars ("Nambu-Goto" bosons) on the string world-sheet. 8 The scalar is an elementary degrees of freedom in the Nambu-Goto action, labeled as X µ . The interaction in the language of the world-sheet and bulk theories is written as On the world sheet, this process can be described by an effective operator which respects the translational invariance of X µ , where we have taken the interaction coefficient to be unity. Here, α is a coordinate on the world sheet. Switching to the canonical normalization for X µ = f −1 a ϕ µ , O is a dimension-three operator with an overall 1/f a suppression. Therefore, the coupling X µ + ψ 0 string → ψ * is suppressed by 1/f a , while the off-shell ψ * propagator gives an additional 1/(y Φ f a ) suppression, so the overall cross section is where s 2 /(16π) is a dimensional estimate of the phase space and the squared amplitude, and we have used Eq. (2.17). σ X is a cross section in 1 + 1 dimension for X µ , and hence is dimensionless. Multiplying by T /π for the average number of X µ modes gives the dissipation rate and falls below the Hubble rate Γ X < H for T < T X , T X = 3.8 × 10 2 GeV 1 y 2 f a 10 10 GeV The last factor in Eq. (3.9) is order unity, and will be neglected.

Overall leakage rate
The overall leakage rate is given by which in practice means whichever rate is the larger one, and so T leak , the temperature at which current destruction is relevant is given by We summarize various scales appearing in this section in Fig. 1 for typical values of couplings. The string network appears at T ∼ f a and the Fermion zero-modes on the string escapes as the massive bulk Fermions by the collision between bulk particles until their temperature drops to m ψ /b. The scattering between standard model particles and X µ decouples at T scat and T x , respectively. The smaller one is the temperature when the leakage due to the scattering stops. The decay of current on strings with finite curvature is not shown here, but is discussed in following sections.

Additional leakage processes for loops
For closed string loops, boost invariance is broken by 1/R, where R is the radius of curvature of the string, and the zero-mode can annihilate/decay into massive fermions. In particular, the trapped KSVZ fermions can decay directly into SM particles via processes such as ψ → Q + h. The asymptotic expansion of the corresponding decay rate for large R is computed in Appendix A. If the final state particles are massless, the decay rate is using the suppression factor in Eq. (A.14). Here k is the momentum of the zero-mode k ∼ 2πI/(N c e ψ ) from Eq. (2.17). The decay rate vanishes as R −2/3 as R → ∞. In general, particles get a mass through the interaction and we need to take it into account. It changes the asymptotic form of the decay rate to using the suppression factor in Eq. (A.16), where m 1 , m 2 are the masses of the two particles. The decay rate of the current depends on the curvature and environment of the strings. We will discuss individual situations in later sections. Another potential leakage mechanism on loops is the collision of zero-modes moving past each other in opposite directions. This interaction is exponentially small in the radius of the loop ∼ e −m ψ R , from the asymptotic form of the zero-mode solution for large r ⊥ .

Vortons
An interesting consequence of the charge and current on the string is that it impedes the shrinkage of string loops. Suppose we start with a string loop with length L and current I. We assume, for the moment, that as long as Eq. (2.19) is satisfied, the loop shrinks and loses its energy as axion radiation [62]. As discussed in Sec. 3, when the temperature of the universe is lower than T leak , current leakage by scattering with the thermal plasma is suppressed. Below we will address the curvature induced leakage mechanism for loops, but since it becomes active only for small loop radii, we first analyse loop stability in the absence of this leakage mechanism. In the absence of (curvature induced) leakage, the charge on a string loop, is conserved and constant. As the loop shrinks, the current grows and eventually saturates. This is seen from the total energy of the loop, which is the sum of the tension and the current energy for L f −1 a and I f a . The latter is given The total energy is where µ Cf 2 a is the tension of the string and we define C as As discussed above, this formula is valid only for I e ψ C/2πf a . However, we assume that this is even valid for I ∼ e ψ C/2πf a for the purpose of estimating the stabilization radius of string loops, as was done in previous studies [42][43][44][45][46]. This was shown for the case of vector-like superconducting strings in [63], and we expect it to hold in our case as well. 9 Setting aside the logarithmic L dependence of C, the minimum of the energy of the loop is reached when (4.5) or equivalently Consequently, the mass of the loop at the minimum is given as Qf a e ψ . (4.7) The loop stops shrinking once the current reaches this value. Such static loops are called vortons, first predicted for the superconducting cosmic string in Ref. [41] and extensively studied in the '90s [42][43][44][45][46]. Later, it was analytically shown that the static loops exist if the chiral current reaches a certain value [64,65]. In addition to the trapped current, the vorton also has a cloud of charge carried by the Goldstone-Wilczek current of the gauge field around the vorton [51]. However, this contribution is suppressed by N c e 2 ψ /4π 2 , and we ignore it for simplicity.
To estimate whether our loops are indeed stable or decay as a result of current leakage mechanisms, we consider three different ways in which they can decay. First, we consider vorton decay by tunneling of the current across the vorton. The tunneling decay rate is of order ∼ e −m ψ Lv and is negligible for Next, we consider the leakage from scattering against plasma particles /string oscillations, analysed in detail in section 3. These lead to the dissipation ratė with Γ leak given in Eq. (3.10). This is because the rate Eq. (3.10) is the rate to ionize a single zero-mode state. The total ionization rate is given by multiplying by N 0 , the zero-mode occupation number, and the current loss rate is thus given by multiplying by e ψ N 0 . On the other hand, I = e ψ N 0 , leading to Eq. (4.9). Consequently, for a string loop with large enough charge, the loops that form below T leak can potentially be stabilized by their current and become vortons. In contrast, loops formed at higher temperatures lose their current to plasma/oscillation leakage processes, and subsequently decay to axions. Finally, we consider the effect of curvature-mediated leakage, discussed in section 3.4. To estimate this rate, we have to first know the typical stabilization radius of our potential vortons. If the curvature-mediated leakage rate is too large at this radius, our vortons will decay.
To estimate the typical stabilization radius in the absence of leakage, we can make use of Eq. 4.5, provided that we know the value of the overall charge Q. We now show how to estimate it. The charge Q is induced on the long strings/loops by thermal fluctuations in the early universe, which lead to an average electric field of E z ≈ T 2 , with a coherence time of 1/T . Eq. (2.16) then gives at a temperature T . I fluctuates over a thermal coherence scale T −1 , and can have either sign. While fluctuations result in non-zero charges on the strings, global charge conservation is maintained via the Goldstone-Wilczek current. The local charge density is ρ = −I. Thus the total charge on a string loop of initial length L is taking the r.m.s. average over L T independent segments. Thus, for L T −1 , Eq. (4.8) is satisfied and all loops can potentially end up as vortons below T leak . We estimate the initial loop length at temperature T as the string correlation length where ξ denotes the number of the strings per Hubble volume. The vorton charge is then (4.12) Here, since the charge on the string moves in the same direction, the Goldstone-Wilczek current is assumed to average the charge to the r.m.s. value. By Eq. 4.5-4.7, the vorton size and mass are then given by 14) The stabilization length Eq. (4.13) was derived while neglecting the curvature mediated leakage process of section 3.4. To estimate the rate for this process, we can use the stabilization length as a proxy for the radius of curvature in Eq. (3.13).
If we assume the daughter particles, the quark and Higgs, are massless, the decay rate of the current, Eq. (3.12), is to be compared with the Hubble expansion rate for k ∼ 2π Nce ψ I v and R ∼ L v . If we take the KSVZ Yukawa to be y 10 −15 , so that the Ψ decay well before T BBN ∼ MeV, Eq. (3.12) then tells us that the current on the would-be vortons quickly dissipates, and they subsequently decay. If, on the other hand y ∼ > 10 −15 so that the bulk ψ decays at relatively late times before BBN, T 10 T BBN , a dominant fraction of the vortons can still remain post BBN, which is in serious tension with observation [66]. In appendix B, we further estimate the resulting vorton density for this scenario. Our estimate is also relevant for other UV models for the KSVZ fermion which could potentially lead to the decay of the bulk Ψ before BBN, while suppressing the leakage Eq. (3.12).
We stress that our estimate for the curvature-mediated leakage is sensitive to the effective mass of the decay products. For example, we should take into account the plasma mass for the outgoing quark and Higgs. According to Eq. (3.13), the decay rate depends on these masses exponentially, i.e.
with m q(H) the quark (Higgs) effective mass. (Un)fortunately, this dependence is not large enough to suppress the decay rate. The situation could be different if for some reason the decay products get an even bigger effective mass in the vicinity of the string -for example by filling the Landau levels from its induced magnetic field. We leave the careful estimation of effective masses for future work, bearing in mind that they could have important implications for the dynamics. In particular, if the effective masses of the decay products end up suppressing the curvature-induced leakage, the resulting metastable "vortons" could be in tension with BBN bounds.

Axion String Evolution in a Primordial Magnetic Field
It is widely believed that the intergalactic magnetic field played a role in the formation of the ∼ 10 µG galactic magnetic fields observed today, via the galactic dynamo effect [67]. The intergalactic magnetic field (IGMF) is characterized by two parameters, the magnitude of the magnetic field B and the typical coherence length λ. An observational lower bound on the IGMF today is set by the non-observation of secondary photons from the emission of highly collimated gamma rays by blazars. Blazar-produced gamma rays collide with extra-galacting beckground (EBL) photons, and produce highly boosted electron-positron pairs. These, in turn, undergo inverse Compton scattering on the CMB, producing secondary photons which could in principle be detected, if they are collinear enough with the original gamma rays. However, in the presence of a large enough IGMF, the electrons and positrons bend and no longer produce collinear secondary photons. Recent non-observations of secondary photon put a lower limit of B 3 × 10 −16 G [68] on the IGMF today for λ 0.1 Mpc. For shorter coherence lengths, the lower bound becomes even stronger by the factor 0.1 Mpc/λ. Though the exact origin of today's IGMF isn't known, one very attractive possibility is that it originated from a primordial magnetic field (PMF). There is a vast literature exploring the different ways a PMF could have been generated in the early universe, a process referred to as magnetogenesis [69,70]. The strongest upper bound on the PMF scenario is B 10 −11 G, obtained from CMB perturbations [71].
In this section, we point out that if a PMF is generated in the early universe, it induces very large currents on axion strings, significantly altering their evolution. As a result, the relic abundance of axions may change by many orders of magnitude. Our estimates are based on a simplified analytic model [72], and so should be taken with a grain of salt until a full simulation is conducted. However, we believe that they capture the right physics, and should be explored further using numerical tools.
For later convenience, we define the comoving magnitude of the magnetic field B ≡ BR 2 and the comoving coherent scaleλ ≡ λR −1 . The dynamics of the magnetic field in the early universe is governed by its interaction with the highly conducting primordial plasma, and is accurately described by magnetohydrodynamics (MHD). The cosmological evolution of plasma-coupled PMF is a complicated problem that was studied by many authors [69,73], and we will not describe it in detail. Here we only include a simplified discussion, and refer the readers to exhaustive reviews [69,70] for more details.
The key quantity in our simplified analysis is the Alfvén velocity v A ∼ ρ B /ρ tot ∝ B [73], where ρ B is the energy density of the magnetic field and ρ tot is the total energy density of the plasma in the early universe. As the universe evolves, the primordial plasma becomes turbulent at scales smaller than the sound horizon v A τ , where τ is conformal time. Consequently, B fluctuations with a correlation length smaller than the sound horizon get damped (processed), while those with a longer correlation length remain unprocessed. Assuming the typical coherence scale is initially small enough, the coherence scale at conformal time τ is given by λ = v A τ ∝ Bτ . This simple linear relation between B and λ is maintained as the universe expands.
The evolution of B is determined by the the unprocessed component of the magnetic field and the MHD plasma. Let n denote the spectral index of the unprocessed part of the magnetic field and the plasma;P (k) ∼k −n for the wavenumber smaller thanλ −1 , whereP is the comoving power spectrum of the magnetic field and of the MHD plasma and is related to the comoving energy density of the MHD flowρ byρ ∼P (λ −1 )λ −3 . Sinceρ ∼B 2 atλ, their conformal time dependence is given as [69] The spectral index n can take several different values depending on the dynamics of magnetogenesis. For simplicity, we focus on a PMF coupled to a compressible plasma, where n = 0, and the helical PMF, whose evolution is equivalent to the case n = −2 [69]: • Compressible plasma: Another case we consider, is a PMF generated during inflation -and so is acausal today. The spectrum for this PMF is expected to be scale invariant, 10 and so the typical coherence scale cannot be defined. However, by a slight abuse notation we defineλ = v a τ . The evolution of the magnetic field andλ is then described by We stress that in reality the magnetic field in the inflationary scenario can be correlated on super horizon scales. In the following, whenever we mention the scale invariant magnetic field, we simply take it to be comovingly constant (B=const.) and coherent over the entire Hubble horizon.
When analyzing the effect of the PMF on axion strings in the early universe, we will make use of the average comoving magnetic fieldB(T ) as a function of temperature. It is useful to relate this quantity to the intergalactic magnetic field B 0 observed today. Naively, the magnetic field decouples from matter at recombination, and the comoving amplitude of the magnetic field remains roughly constant afterwards. From Eq. (5.1), we haveB where T 0 is the current temperature of the photon and T rec(eq) is the temperature recombination (the matter-radiation equality). Indeed, a detailed study [73], including viscosity effects and MHD dynamics post matter-radiation equality results in a similar conclusion for T 10 MeV; (5.2) In this paper, we use the above formula for the PMF as a function of temperature. Note that the coherence length at each time is automatically determined from the Alfvén velocity onceB is fixed. For T > T rec [73], where r(T ) is the ratio between the magnetic energy and a power of the entropy density [75], The ratio r would have been a constant if not for the coupling of the PMF to the primordial plasma. 11 After recombination, the coherence length is frozen as well. In terms ofB 0 = B 0 , the coherent length of the PMF is [68,73,76] λ 0 11 pc B 0 10 −13 G . (5.5) Combining this with the lower bound on the IGMF introduced in the beginning of this section, we obtain Thus, we take 10 −13 G as the reference value of the causal (n = −3) PMF today. For acausal (inflationary) PMF, we take the direct lower bound 3 × 10 −16 G, as the PMF is assumed to be coherent on a 0.1 Mpc scale.
Having described in some detail the dynamics of a PMF in the early universe, we turn to consider its effect on the cosmic evolution of QCD axion strings. For each one of the PMF scenarios we consider, we vary the initial magnetic field 12 B mag so that its value today is consistent with observations. Since axion strings are superconducting, their interaction with the PMF yields sizeable electric currents, which, as we shall see, back-react on their evolution.
By Lorentz symmetry, a boosted magnetic field becomes an electric field. When the axion string moving with the velocity v s is embedded in the PMF, an electric current develops on the string generated by the effective electric field E ∼ Bv s . The evolution of the current is governed by the electric field, Hubble expansion and current dissipation, we obtain the current evolution equation [77], with the addition of a leakage termİ whereİ denotes the derivative with respect to the physical time and I does the derivative with respect to the spacial coördinate on the string. The first term on the right hand side is a direct consequence of Eq. (2.7), while the second term accounts for the redshift of the current due to Hubble expansion. The factor β encodes the effect of the small scale structure of the string; if the string is infinitely long, the current redshift is as usual, β = 1 whereas if the string is very wiggly, Hubble expansion straightens the wiggles first so that β ∼ 0. We assume β ∼ 1 in this paper. The current destruction rate Γ leak is given by Eq. (3.8) and Eq. (3.3), and we conservatively assume y = y Φ = 1 in the following calculations, corresponding to an unpressed leakage rate. Here, we may ignore the curvature-induced current decay, Eq. (3.13), since the radius of curvature of the long strings is sufficiently small. We elaborate on this point below. For simplicity, before solving Eq. (5.7) numerically, we first take its root-mean-square average over a Hubble time. Since the string velocity is coherent over the typical inter-string distance, L ≡ t/ √ ξ, the root-mean-square current I RMS (t) includes a suppression factor of √ coh τ current with coh = min(λ, L) and τ current ≡ min(H −1 , Γ −1 leak ). In writing Eq. (5.7), we implicitly assumed that the main mode of energy loss by the strings is string reconnections, which results is a loss of both string length and charge. However, if some of the energy of the strings is directly radiated into axions without the loss of charge, this could induce an additional blue-shift factor, which enhances the current. This effect has been previously considered in [78,79]. However, we use the above formula as a conservative estimate.
According to Eq. (5.7), the current at a given Hubble time is roughly where B eff ≡ B √ coh τ current . If it were not for friction and coh ∼ τ current , v s would be O(1) and I ∼ e 2 ψ 2π r(T )M p , which would be much larger than the thermal fluctuation given in Eq. (4.10). Thus, potentially, the current on the axion string by the PMF may significantly affect the evolution of the string. Of course, the current on the string induces friction, as discussed on Sec. 4. The string velocity is decelerated by the friction force as Eq. (B.11). Thus, in order to estimate the current, we need to follow the evolution of the string including friction and evaluate the velocity.

String Evolution with the PMF
The early universe evolution of cosmic string networks has been the subject of a vast number of analytic and numerical analyses. These have mostly been performed for nonsuperconducting local and global strings, as well as for local superconducting strings. We currently do not know of any dedicated simulations of the properties of anomalous superconducting global strings, and so we encourage the study of their reconnection and emission properties. In the meantime, we extrapolate some of these properties from those of global/local-superconducting strings. In particular, we assume that: • Their reconnection properties are similar to those of local superconducting strings and global strings up to O(1) corrections.
• Their axion emission rate is similar to that of non-superconducting global strings -since it is orthogonal to the dynamics of the zero mode currents on the strings.
• The plasma friction on the superconducting strings is similar to that of nonanomalous global superconducting strings [80]. This should not be different for axion strings, since it is not linked to the PQ symmetry being anomalous.
Under these assumptions, we can analyse the early universe evolution of superconducting axion strings using the analytic velocity dependent one-scale model (VOS) [62,72], which averages the time evolution of the string energy and velocity using the Nambu-Goto action. In terms of the average inter-string distance L and the average string velocity v s , the equations are In the first equation, the first term on the right hand side is from Hubble expansion. The second term encodes the effect of string loop emission by string reconnection. c L is a O(1) factor of the loop-chopping efficiency, taken as c L ∼ 0.66 [81], and we assume it remains constant. The third term is from friction due to the current. The friction relaxation timescale, τ fric , is defined in Eq. (B.9) for the thermal current, but the formula also holds for any current, regardless of its source. As for the second equation, here the first term is from acceleration of the string by its curvature, which depends on the small scale structure of the string, encoded in the factor k s = 0.25 [81]. The second term is deceleration by Hubble expansion and friction. Finally, we note that the VOS model for global strings has been discussed in [81], and it involves two modifications to Eq. 5.9. The first is the logarithmic dependence of the string tension on the correlation length L. We have explicitly incorporated this effect into all of our equations. The second is an extra term sv 6 /C on the right hand side of the L equation, representing the energy loss axion radiation as in [81]. We checked that it leads to a negligible modification of our dynamics (it is subleading with respect to c L v), and so we omit it for simplicity. Lastly, we checked that O(1) changes in the values of c L and k s lead to O(1) changes in the resulting number of strings. Below we sill see that the current leads to at least a 4 order of magnitude effect in the number of strings, and so we simply fix our values for k s and c L to those of [81].
Since τ fric is determined by the current on the string as calculated in Eq. (B.9), we combine the Hubble-time average of Eq. (5.7) and Eq. (5.9) to determine the evolution of our current carrying string network embedded in a PMF. For the PMF time evolution, we simply assume that the magnetic field instantaneously appears at T = T init and evolves according to Eq. (5.2) for the causal PMF, n = −3. In the following, we vary T init from around 10 GeV to 10 3 GeV as an example. In this setup, the string evolution experiences four distinct stages: (1) the evolution before the PMF appears (2) a transient stage in which the current grows and the string decelerates (3) the stretching phase (4) the friction dominated regime. We find the string never enters the scaling regime after a PMF appears. Below we introduce these stages, and the key processes that characterize them.
In Eq. (5.7), we ignored the curvature-induced current decay. Typical curvature of strings are of order the string correlation length, L. Since the string network is in the thermal bath, the standard model quark and Higgs have the thermal masses of order g 2 s T and g W T , respectively, where g s(W ) is the strong (weak) gauge coupling. Thus, we may ignore Eq. (3.13) if the momentum of the current k satisfies k √ LT 3 . In the calculation presented below, we confirmed that this condition is satisfied.

Evolution before the PMF appears
Before the PMF appears, the evolution of the axion string follows the standard evolution [13,15,17,62,72,82], with the extra effect of the thermally induced current on the string. According to the standard picture, the string evolution consists of two phases, the friction dominated regime and the scaling regime.
Just after the phase transition, numerous strings appear [15,17]. They strongly interact with each other and quickly decrease forming loops. The number of strings immediately stops decreasing after the phase transition , as friction due to the thermal current becomes effective, and so the strings stop moving. As a result, the number of the string per Hubble horizon, ξ, obeys Eq. (B.12). As the temperature of the universe drops, friction becomes weaker and eventually negligible. The Kibble temperature, the temperature when friction becomes negligible, is when τ fric ∼ H −1 . To be explicit, it is T Kibble = π 2 g 90 µ M p 140 C g 106.75 1/2 f a 10 10 GeV 2 GeV . (5.10) • Scaling regime. Below T Kibble , the universe is cold enough and friction is negligible. The string velocity becomes O(1) constant and only ξ = O(1) = constant strings exist in one Hubble volume due to efficient reconnection. As a result, the ratio between the string and radiation energy is constant. This is the scaling regime of the string. Using the VOS model, one may calculate the typical distance and velocity of the string assuming L ∝ t, v = constant as Note that in recent years there has been a growing interest in a potential logarithmic violation of this scaling law, ξ > 1, backed by dedicated simulations [30][31][32][33][34][35][36]81, 83] though we do not discuss this further.

Transient stage: current growth and string deceleration
At time T init , the current growth is "turned on" according to Eq. (5.7). Then the current grows rapidly, and there is large plasma friction on the string. The friction relaxation time τ fric is much shorter than the Hubble time leading to string deceleration, which in turn feeds back on the current growth equation via its velocity dependence. We call this the transient stage. The whole process ends in a much shorter period than the Hubble time. Since the VOS equation, Eq. (5.9) assumes time-averaging over scales ∼ H −1 , such shorter time scale physics cannot be discussed using the VOS equation. For simplicity, we assume L does not change during the transient stage and solve Eq. (5.7) taking spatial averages for I and v s . In Fig. 2 we present the numerical solution of this initial stage. As we can see in the plot, the current I grows rapidly until it reaches a plateau, and the velocity decreases accordingly. Indeed, the transient stage ends very quickly, which justifies our assumption. We can estimate the plateau values of I and v by taking R, H, B, λ and L as constants equal to their values at t init . The plateau values are then obtained by demanding We also show these plateau values for v s and I in Fig. 2. Numerically, we use these values as initial conditions of the next phase.

Stretching phase
The initial transient phase terminates when v s reaches a temporary plateau at v s 1. At this point reconnections are suppressed, because the mean inter-string distance is larger than v s H −1 , and so the number of the string per comoving Hubble volume is constant. This phase is called the stretching regime, since the string is stretched by the Hubble expansion. Let us see how this dynamics emerges from in VOS model. Neglecting the velocity dependent terms in the L(t) Eq. 5.9, we obtain dL dt = LH .

Plasma friction domination
The stretching regime ends when friction becomes weaker, v s increases, and v s H −1 becomes comparable to the typical string distance. At this point reconnections commence, but the plasma friction is still strong and the whole network evolution again follows the friction dominated regime. The typical string distance and velocity asymptotically follow Eq. (B.12) with L = t/ √ ξ. Let us make a rough estimate of ξ. From Eq. (B.12), we obtain where the subscript 0 denotes the asymptotic value. On the other hand, the current is Thus, we expect as the asymptotic value. For temperatures lower than T leak , current dissipation is ineffective and we expect Note that for all regions of interest, ξ 1 and the axion string never enters the scaling regime if the PMF exists, as we numerically show later.

Numerical result
In Fig. 5.1.5 we plot the number of strings per Hubble volume ξ = (t/L) 2 and the velocity v s as a function of temperature. We take f a = 10 10 GeV, T init = 10 3 GeV and B 0 = 3 · 10 −13 G as an example and show the result for α = 1/3 and 3/5. For comparison, we also show ξ and v s without the PMF.
The figure demonstrates the dynamics discussed above. First, before T init , the axion string is at the onset of the scaling regime. Without the PMF, the thermal current decays as the temperature drops and the string network enters the scaling regime, although there is an O(1) change in ξ, Eq. (B.12). However, if the PMF appears, it produces a large current on the string and the string is immediately decelerated by plasma friction. For several orders of magnitude in temperature after T init , the string velocity is too small, and string reconnections are suppressed. Thus, the number of strings per comoving Hubble horizon is conserved and ξ grows as t ∼ T 2 . This regime of string evolution is called the "stretching phase." As the universe gets cold, the PMF and friction force become weaker, allowing the string to accelerates because of its curvature. Eventually, the string network starts to reconnect again and the network goes over to another friction dominated regime. The temperature at which the stretching phase ends depends on the magnitude of the current on the strings, which in turn, depends on the PMF. The larger α is, the faster the PMF decays. Since we assume the present magnitude of the PMF is the same for both α = 1/3 and 3/5, the PMF at early times is stronger for larger values of α. Thus, axion strings embedded in a PMF with α = 3/5 spend more time in the stretching regime.
For causal magnetic fields, the temperature of magnetogenesis is bounded from above by the requirement that the PMF does not dominate the expansion of the universe. For that reason, we take T init = 10 3 GeV, which is the relevant upper bound for α = 3/5. For this value of T init , the α = 1/3 case enters the friction dominated regime at around T ∼ 10 GeV, while in the α = 3/5 case, the network is approximately still in the stretching regime as late as 1 GeV. In comparison, if T init is lower than 10 3 GeV, the network is still in the stretching regime at the time of domain wall formation. Eventually, The tension of the domain wall wins over the friction force and the string-domain wall network shrinks and fragments.

Axion Relic Abundance
Here, we estimate the relic density of axion radiated from axion strings embedded in a PMF. Naively, an increased number of strings per Hubble in leads to a proportional increase in the number density of radiated axions. However, as reference [36] points out, the actual increase in the axion relic density is roughly proportional to √ C ξ. The  reason for that is as follows: The contribution to the relic abundance is biggest when the axion mass starts affecting the dynamics. Since the inter-string distance is H −1 / √ ξ, each emitted axion carries a momentum of order √ ξH and the number of emitted axions is then Assuming the number of axions is conserved, the final relic abundance is proportional to √ ξ. In reality, three corrections should be added to the above estimates. First, since the tension of the axion string is not constant but proportional to C, we have to include this factor in our estimate. Secondly, we need to take into account how the axion mass varies with temperature. As a result, Ref. [36] concludes where Ω ax is the relic abundance of the axion today, γ is the temperature dependence of the axion mass at around 1 GeV, m a ∝ T −γ/2 and we take γ 8. The estimate Eq. (5.18) has also been validated by numerical simulations in [36]. Finally, we have to account for the different emission properties of current-carrying axion strings with respect to current-less axion strings. We do this in two stages. First, we need to account for the different reconnection/loop emission rate in our scenario. This can be directly accounted for by comparing the energy loss rate of our long strings to that of standard axion strings. To account for this difference, we multiply Eq. 5.18 by the suppression factor S ≡ Γ a /(ρ s /t) , (5.19) which could be as small as 10 −2 in the stretching regime. 13 To calculate Γ a in our scenario, we evaluateρ s for our network evolution, as compared toρ free s , an equivalent network of non-interacting strings with the same instantaneous string density at time t, ρ s | t = ρ free s | t . The difference between the two terms is the instantaneous energy radiated into loops, and eventually into the axion field [30]: In Fig. 4 we present the rate Γ a calculated in this manner, for the network evolution described in Fig. 5.1.5. As can be seen in the plot, this rate is suppressed in the stretching regime, but grows to 1 in the friction dominated regime. Note that even in the presence of plasma friction, the dominant source of energy loss for our strings is reconnections, which go like c L v s in Eq. 5.9, and not plasma friction, which enters as v 2 s in the equation for L (equivalently, the string energy). This is directly evident from Fig. 4.
Finally, we need to account for the different emission properties of loops in our scenario with respect to the standard case, where all of the loop energy is effectively radiated to axions. In our case, however, the energy of the shrinking loops can also be transferred to the plasma by the friction. If this effect dominates, it leads to a large suppression in the density of radiated axion DM. This is an important point to explore in future work, but it is beyond the scope of the current paper. In particular, the resulting axion relic density strongly depends of loop emission properties at later time, T ∼ 1 GeV. In this case the effect of domain walls must be taken into account. Here we do not delve further into this analysis, but simply assume that axion emission prevails as the dominant energy loss for loops at T ∼ 1 GeV.
In Fig. 5 we show contour plots of ξ at 1 GeV for the causal PMF, n = −3, as a function of the temperature T init at which current starts to accumulate, and the magnetic field B 0 today. f a is taken to be 10 10 GeV here. Note that throughout this section, we assume current leakage is efficient enough, y = 1 in Eq. (3.11). The shaded region is excluded by the magnetic field dominating the universe at magetogenesis [69]. We also confirm that the total energy of the current does not exceed the energy of the magnetic field so that its backreaction on the PMF is negligible.
We present two different contour plots for α = 1/3 and α = 3/5. For lower T init , ξ at 1 GeV does not depend on the magnitude of the PMF. This is because the string network is still marginally in the stretching regime at 1 GeV. Indeed, ξ for α = 1/3 and α = 3/5 are almost the same. On the other hand, in particular for α = 1/3, if  . Energy emission rate from axion strings into the axion field, normalized by the approximation ρ s /t. In the friction dominated regime, this ratio becomes unity. In the stretching regime, the actual rate is suppressed by up to 10 −2 due to the low reconnection rate.
T init is high enough, ξ becomes independent of T init . There, the string network enters the friction dominated regime at 1 GeV and ξ gets closes to the asymptotic value, ξ 0 . In Fig. 6, we show ξ for the inflationary PMF in terms of the magnetic field today. Note that since inflation ends before the Peccei-Quinn phase transition, we may assume T init = T PQ = f a . Since the inflationary PMF exists before the PQ symmetry breaking, the network enters the friction dominated regime well above T = 1 GeV and ξ takes its asymptotic value ξ ∝ B 4/7 0 (see Eq. (5.16) with coh ∝ ξ −1/2 ). Finally, we show the axion abundance. In Fig. 7, we show Ω ax determined at 1 GeV for the causal PMF, as a function of the temperature T init at which current starts to accumulate, and the magnetic field B 0 today. We refer to Ref. [30] for the value of Ω ax when B 0 = 0 and scale it by S(Cξ) 1/2+1/(4+γ) . The parameter region where vorton abundance could exceed the dark matter abundance is indicated by the green shaded region. Note that this part of the parameter space is likely still viable, as long as y 10 −15 and so the vortons are unstable. Compared with Fig. 5, Ω ax is suppressed by the suppression factor. In Fig. 5, we find for lower T init , ξ is independent of B 0 , since the string network is in the stretching phase and the ratio of the comoving Hubble horizon solely determines ξ. However, the suppression factor is smaller for larger B 0 , since the friction becomes larger and the reconnection occurs less, resulting in smaller S. As a result, for fixed T init 10 2 GeV, where the network is in the stretching phase, the abundance is larger for the smaller B 0 (> B min ). Also, we see ξ does not depend on α if the string network is in the stretching regime. However, the PMF is stronger for α = 3/5 scenario and S becomes smaller, as we discuss. Consequently, for fixed T init 10 2 GeV, the axion abundance for α = 1/3 is larger. In Fig. 8, we show Ω ax determined at 1 GeV and Ω v for the inflationary PMF in terms of the magnetic field B 0 today (scaled by the observed lower bound B min ≡ 3 · 10 −16 G), for f a = 10 10 GeV. As we have discussed, for the inflationary PMF, the string network enters the friction dominated regime at 1 GeV. Thus, the suppression factor is no longer effective and Ω ax scales as ∼ √ ξ. To derive the implications for the cosmological bound on f a , we now present our results for the axion relic density, while scanning over f a . For lower values of f a , the string tension is lower and the effect of friction is larger. This makes ξ larger for these small values of f a . On the other hand, Ω ax (ξ = 1) increases with f a , as is well known. The axion abundance is determined by the balancing of these two factors.
In the case where the network is deep in the friction dominated regime by T ∼ . ξ at T = 1 GeV as a function of B 0 today, scaled by the observed lower bound B inf min ≡ 3 · 10 −16 G, for the inflationary PMF, whose coherent length is assumed to be larger than the current Hubble horizon. f a is taken to be 10 10 GeV.
1 GeV, we can use Eq. (5.16) to estimate ξ ∼ µ −2/3 with ξ dependence of coh ignored, and so Ω ax ∼ f −2/3 a . On the other hand, Ω ax (ξ = 1) ∼ f a . Thus, in total, we expect Ω ax ∼ f 1/3 a in this limit. In practice the network is either in the stretching regime or on the cusp between the two regimes, and so there is an extra f a dependent contribution from the suppression factor S.
Our numerical results are shown in Fig. 9. For α = 1/3, we show Ω ax /Ω DM vs f a for 3 different values of B 0 . The string network is not deep in the stretching regime and the suppression factor S| 1 GeV is O(1) in the entire parameter space. Thus, as we expect, Ω ax is a mildly increasing function of f a , although the exponent is different from 1/3, since the string network is only at the cusp of the friction dominated phase. We find, for example, if B 0 = B min , the correct relic abundance is given for f a ∼ 10 7 GeV. For α = 3/5, we show Ω ax /Ω DM in terms of f a with changing T init . In this case, since the magnetic field in the early universe is larger, the whole network may be in the moderate stretching phase and the suppression factor, S, may be small.
For smaller f a , the eventual axion relic density is affected not only by the PMF, but also by friction from the thermal current on the strings, prior to magnetogenesis. This sets a different initial condition with a higher value of ξ at the onset of the stretching regime. As a result, the number of strings in the Hubble volume at T init is larger for smaller f a and the transition from the stretching to the friction dominated phase occurs earlier. Thus, Ω ax /Ω DM is no longer monotonically decreasing in f a , as can be seen in the right panel of Fig. 9 and in Fig. 10.
Finally, in Fig. 10 we present a 2D scan of Ω ax /Ω DM as a function of f a and T init for the minimal value B 0 = 10 −13 G of the magnetic field today. As mentioned before, for T init 10 2 GeV, the string network for both α = 1/3 and α = 3/5 is in the stretching regime, and so their value of ξ is simlar. Even so, the value of Ω ax for α = 3/5 is smaller than the one for α = 1/3, because of its greater friction that leads to a slower string velocity and a smaller value of S.
The axion abundance is also presented for the inflationary scenario in Fig. 11. In this case, we may assume the asymptotic value of ξ. As a result, Ω ax /Ω DM scales like a B 2/7 as one can check in the figure. As can be seen from the plot, the inflationary PMF scenario with axion strings is in tension with the astrophysical limit of f a 10 8 GeV from supernova cooling [84].   Figure 9. Ω ax /Ω DM as a function of f a for a fixed T init and B 0 today, scaled by the observed lower bound B min ≡ 10 −13 G, for the causal PMF. Left: α = 1/3 and we change B 0 and see how Ω ax behaves for T init = 10 3 GeV. Right: α = 3/5 and we change T init and see how Ω ax behaves for B 0 = B min .

Conclusion and Discussion
In this paper, we point out the comsmological implications of the chiral superconductivity of axion strings. We examined the possibility that the axion string network leaves stable remnants, called the vortons, and found that they largely decay away as a result of curvature induced current leakage. If this leakage is suppressed by a small coupling y ∼ 10 −15 , then vortons could potentially lead to tension with BBN bounds, as well as to overclosing the universe. In addition, we analyzed the evolution of the axion string network in the presence of a PMF within the analytic VOS model of string evolution. We found that the induced current on the string is sizeable, and so it generates a dominant friction force on the string which suppresses string reconnections. This increases the number of strings per Hubble volume by many orders of magnitude, and could greatly modify the relic abundance of axion DM. We stress that our results have been obtained within the simplified VOS model for cosmic string evolution, and so need to be validated by future numerical studies.
where (r, φ, z) are cylindrical coordinates, and L = 2πR is the length of the loop. kR 1 must be an integer. ψ 0 is constant in theẑ direction between [−∆/2, ∆/2] and for radii r between R − ∆/2 and R + ∆/2. Eq. (A.1) models the transverse form of the zero-mode wavefunction as a step function. The details of the functional form of ψ 0 in the transverse direction are not important; the consequences of energy-momentum conservation arise from the t and φ dependence of ψ 0 . We are interested in the decay of zero-modes with momentum k∆ 1, so that their momentum is smaller than the inverse size of the zero-mode. The energy of the zero-mode is ω = k.
For free-particle decays k → q 1 + q 2 + . . ., the matrix element involves a space-time integral The time integral gives an energy-conserving δ-function δ(q 0 1 +. . .+q 0 n −k 0 ) and the space integral gives a momentum-conserving δ-function δ(q 1 + . . . + q n − k). We consider a matrix element q 1 , . . . q n |ψ 0 where the final states are plane-waves, but the intial state is the zero-mode trapped on a loop. The time integral still gives an energy-conserving δ-function δ(q 0 1 + . . . + q 0 n − ω). However, the space integral In the regime of interest, |q| ∆ 1 so the Bessel function is approximately constant, over the integration domain, and (using L = 2πR) The square of the matrix element has the factor which is the analog of where V is the volume of space. We want the limit of Eq. (A.8) for q R 1, i.e. for a large loop. From Ref. [86, §8.4], Using this asymptotic expansion in Eq. (A.8) with gives For the two-body decay ψ 0 → q 1 + q 2 , the decay rate using the asymptotic expansion Eq. (A.12) is We can find the decay rate for large R by an asymptotic expansion of Eq. (A.13). If the final state particles are massless, then q , with q = q 1 + q 2 can be as large as q = k, if the two massless particles are emitted in the same direction. The asymptotic form of the integral for large R is dominated by this region, and is where Γ 0 is the free-particle decay rate (i.e. for a free particle not trapped on the string). 14 Eq. (A.14) has (k∆) 2 and 1/(kR) 2/3 suppression factors relative to the free particle decay rate. The final result vanishes as R → ∞, as it must by the boost invariance argument for a straight string.
If the final particles are massive, with masses m 1 and m 2 , then the maximum value of q is when both particles are emitted in the same direction with the same velocity, with boost factor γ given by and the smallest value of k 2 − q 2 is (m 1 + m 2 ) 2 . The asymptotic form of the integral is dominated by this point, and is provided the exponent is much larger than unity.

B Vorton Abundance for Small y
In this appendix we estimate the vorton relic density in the absence of leakage processes, and without a PMF. While in the case of KSVZ fermions with a Yukawa coupling we have seen that vorton are generically unstable, here we provide the complementary analysis for the case that the leakage processes are suppressed. Though we do not have a particular UV model in mind that allows for the decay of KSVZ fermions while preventing current leakage, we are interested in the projected vorton abundance in this case.
In order to estimate the vorton energy density, we first estimate the number of vortons formed at a given Hubble time. From the evolution of the string network, we can calculate the total length of strings disappering from the network at that time. The network reduces its length by radiating axions or forming loops. Analytic estimates show the string shrinkage by axion emissions is at most the same order contribution as one by loop formations [62]. For simplicity, we assume all the extra length forms loops and eventually becomes vortons.
The comoving number density of strings (per unit area) is where R(t) is the scale factor of the universe. n str would have been conserved if not for the string reconnections. Thus, the total length of the string which forms loops, L tot (t), is proportional to the time derivative of n str ; Within a Hubble time, where we take R(t) = 1 without any loss of generality. The loops of string eventually ends up as vortons. Thus, the number of the vortons at a given Hubble time is given by the ratio between L loop, tot and the typical size of each loop, L loop . Here, for simplicity, we assume that the sizes of all loops at a Hubble time are the same. For L loop = L 0 = ξ −1/2 H −1 , the number of the loops is denotes the efficiency of string reconnections. The vorton energy density at a given Hubble time is ρ v E v N loop H 3 . In order to compare it with the dark matter number density, it is convenient to normalize it by the entropy density s: which is a comoving quantity. In the above estimate, assumed that the vorton does not annihilates after it forms. This may be a reasonable assumption for particles, but since the size of the loop is macroscopic at the beginning, there may be some annihilation with other vortons. The opposite limit from no annihilation is to assume that the vorton annihilates efficiently, and only the net current asymmetry of the initial vorton distribution survives. In this case, we can simply replace N loop → N loop , obtaining the same equation as Eq. (B.6) with k ξ 5/4 → √ k ξ 1/2 . Note that we have focused on the electromagnetic current, but the same dynamics holds for the color current on the axion string as well. A similar analysis of the color current gives the same order-of-magnitude estimate of the vorton abundance.
So far, we have assumed that all loops formed at a Hubble time are initially the same size, L 0 = ξ −1/2 H −1 . In reality, the sizes of loops distributes over the wide range of scales [24,87]. If we take smaller loop sizesL, the vorton charge and mass decreases by L /L 0 whereas the number density increases by L 0 /L for the fixed energy of loops. On the other hand, the estimation of the total charge asymmetry inside a horizon does not change. Thus, if we included more detailed spectrum of string loops, the naïve estimation would increase but the conservative one would not change.
We now estimate the value of ξ, the number of strings per Hubble. In the early universe, the friction force due to scattering with plasma particle decelerates the string and ξ > 1 in general [61,88]. The scattering cross section between the string and thermal particles with a non-zero Peccei-Quinn charge can be as large as allowed by unitarity, λ ∼ 1/T [18]. In the KSVZ scenario, such PQ charged particles are absent in the thermal bath. The axion may exist in the universe, but the scattering cross section between the axion and the string is expected to be geometrical. This is an analogue of the scattering between the magnetic monopole and the photon as the magnetic dual of the axion couples to the string. However, since the string supports the (color) current, a (color) magnetic fields winding the string play a role of an effective thickness of the string. This results a scattering rate as large as the unitarity bound. Ref. [79,80,89] estimate the cross section as where ρ is the energy density of the thermal plasma. This roughly corresponds to the radius where the pressure between the thermal plasma and the magnetic field induced by the current balances. According to the standard picture [82], ξ and the string velocity v s are determined in the following way. First, the frictional force per unit length of string is where is the thermal average for the scattering particles, v rel is the relative velocity between the particles and the string, n is the number density of the scattering particles and q is the momentum transfer. The typical time scale due to friction is the relaxation time of the string momentum, For a typical string distance scale L ∼ ξ −1/2 H −1 , the force on the string per unit length is from the string curvature F str µ L . (B.10) The string accelerates for a time ∼ τ fric and the string velocity is thus limited to v s F str µ τ fric µ LT 2 I . (B.11) On the other hand, ξ is also given as L v s H −1 if we assume the reconnection process occurs, so thatṅ str = 0. Then, v s ξ − 1 2 , ξ T 2 I µH , (B.12) which justify the assumption. Note that since v s < 1, ξ stops decreasing when it reaches a value of order one. For simplicity, we thus assume the number of the string per Hubble horizon is max (ξ, 1). In the above estimation, we assume the number of the string just after the PQ phase transition at T PQ , where we assume T PQ = f a , is more than ξ(T PQ ) so that the number of the string is monotonically decreasing after the phase transition, which is a reasonable assumption [17]. Using Eq. (B.12) with Eq. (B.6), we can calculate ρ v /s for a given Hubble time and obtain ρ v /s ∼ T 15/4 . Thus, the temperature T v when most vortons form is T v = T leak and the relic abundance of vortons today is Ω v = ρv(T leak ) s(T leak ) s 0 ρ 0 , where s 0 (ρ 0 ) is the entropy (energy) density of the universe today.
We plot the vorton abundance normalized by the dark matter density in terms of T v in Fig. 12. As we have discussed, we show two estimations, assuming no annihilation or efficient annihilation of vortons, as a solid and dashed line, respectively. As the universe gets cold and ξ decreases to ξ ∼ O(1), the two estimates are expected to coincide since the number vortons in the Hubble horizon is ξ ∼ O(1). Indeed, one can see that after the kink in the plot, at which temperature ξ reaches O(1), the solid and dashed line go together.
Note that our analysis of vortons is different from previous work [42][43][44][45][46] in several respects. As far as we can tell, our work is the first to show that vortons are inevitably generated by axion strings. Secondly, our estimate for the overall charge per vorton differs from previous analyses [46]. The underlying reason for our different estimate is that the Goldstone-Wilczek in the bulk serves to average the current and the charge over the string loop, leading to an overall charge of Q ∝ √ T L instead of Q ∝ T L. Consequently, our analysis results in smaller vorton abundances. Finally, though previous analyses concluded that the chiral current on the vortons is completely conserved, we explore several scattering processes that lead to current leakage from the string. When current leakage is efficient enough, the vorton abundance is greatly suppressed, and there is no risk of vortons overclosing the universe. In particular, since T leak < f a /b < f a , the vorton abundance is insensitive to physics at the PQ scale. After the QCD phase transition, each string loop becomes the boundary of a domain wall. As the whole system shrinks, the edge current on the string increases. When the current reaches Eq. (4.6), where the energy of the string is at the minimum in terms of the length, the shrinking of the loop tends to be halted since it is no more energetically favored. On the other hand, the wall is energetically likely to continue to shrink. In order that the wall shrinks with the perimeter constant, it is expected that the domain wall-string system is twisted and eventually breaks into many vortons whose size is around the domain wall thickness. Vortons produced by this process are negligible since the energy is only a tiny fraction of the total energy of the string-domain wall network giving the axion dark matter. Note that if the size of the vorton is larger than the axion compton wavelength (i.e. the domain wall thickness), the same process happens and the vorton breaks into many small vortons. However, the total energy of the vorton is conserved, since it is proportional to the total charge, which is conserved during this breaking process.
It is worth asking whether the superconductivity alters the domain wall dynamics. A domain wall supports some topological field theory to compensate the anomaly on the string around the wall [90]. For example, the electromagnetic anomaly on the string infers that the Chern-Simons theory for the electromagnetism lives on the domain wall, just like the (anomalous) Hall effect [91]. The wall thickness is of order of the axion mass and much thicker than T −1 , the typical scale of the thermal bath particles. Thus, we expect such IR structure of the wall does not qualitatively change the picture of the time evolution of the wall, although detailed studies may be needed to determine the precise abundance of axions from domain walls.
Finally, let us comment on the potential detection of vortons from the early universe. The charge of the vorton is gigantic 10 3 as Eq. (4.8) requires. The situation is similar to an atomic nucleus with a very large atomic number. Once the charge exceeds 1/α ∼ 137, Schwinger pair production [92] occurs and the charge is screened to be Q 137, although the magnetic flux cannot be screened. Similarly, the quark and gluon cloud neutralize the color charge of the vorton after QCD confinement. The color magnetic flux is Higgsed due to the dual Meissner effect [93,94]. Thus, the vorton looks like a charged exotic hadron with a large magnetic dipole moment. The detection of such an object is not easy [95,96]. The vorton is screened by the Earth's crust and underground direct detection experiments cannot detect it. In addition, the vorton is much heavier than f a and the number density is low. Thus, searches for exotic hadrons do not place any constraints. Consequently, we require Ω v < Ω DM based on Ref. [46] although a more detailed discussion of the constraints on vorton abundance remains to be studied.