Signatures of Primordial Energy Injection from Axion Strings

Axion strings are horizon-size topological defects that may be produced in the early Universe. Ultra-light axion-like particles may form strings that persist to temperatures below that of big bang nucleosynthesis. Such strings have been considered previously as sources of gravitational waves and cosmic microwave background (CMB) polarization rotation. In this work we show, through analytic arguments and dedicated adaptive mesh refinement cosmological simulations, that axion strings deposit a sub-dominant fraction of their energy into high-energy Standard Model (SM) final states, for example, by the direct production of heavy radial modes that subsequently decay to SM particles. This high-energy SM radiation is absorbed by the primordial plasma, leading to novel signatures in precision big bang nucleosynthesis, the CMB power spectrum, and gamma-ray surveys. In particular, we show that CMB power spectrum data constrains axion strings with decay constants $f_a \lesssim 10^{12}$ GeV, up to model dependence on the ultraviolet completion, for axion masses $m_a \lesssim 10^{-29}$ eV; future CMB surveys could find striking evidence of axion strings with lower decay constants.

Axion strings are extended topological defects stretching over cosmological distances that may develop for high post-inflationary reheat temperatures.For example, if the axion arises as the pseudo-Goldstone boson of a spontaneously broken global U (1) Peccei-Quinn (PQ) symmetry, then axion strings develop so long as the reheat temperature is above the temperature of U (1) symmetry restoration.The axion strings are characterized by the property that the axion field, which is a periodic field, undergoes a full field excursion when traversing a circle encompassing the string core; to resolve the singularity at the string core, the heavy radial mode of the PQ complex scalar deviates from its vacuum expectation value (VEV) and sends the full complex scalar field to zero at the core center.Apart from at the string cores, the radial mode is otherwise frozen at its VEV for temperatures well below that of PQ symmetry breaking.The axion-string network evolves to maintain an approximate scaling solution, where there is roughly one string per Hubble patch at any time (see [1] for a review).Fig. 1 illustrates a snapshot of an axion-string network in the context of a cosmological adaptive mesh refinement (AMR) simulation performed in this work.
Axion strings were originally discussed in the context of the quantum chromodynamics (QCD) axion, which was introduced to solve the strong-CP problem [2][3][4][5].It Figure 1.A zoom-in of the axion-string network as realized in a PQ-Higgs field simulation, which is discussed in Sec.V. We show the energy density in axion radiation in a 3D volume rendering enclosing approximately 1.9 Hubble volumes at log(ms/H) ≈ 7.2; the string network evolves to maintain the scaling solution by emitting energy into relativistic axion modes.As we discuss in this work, however, axion strings can also efficiently produce high-energy SM radiation through the production and subsequent decay of heavy radial modes and through the direct production of SM Higgs bosons.Inset in this figure is the outline of the AMR grid structure, showing the different refinement level locations at this snapshot.The finest refinement level (yellow boxes) is nested within the intermediate level (red boxes) and is mostly localized around the string cores.The intermediate level covers a larger area and ensures proper numerical convergence of the outgoing radiation wherever necessary.The blue lines correspond to the string cores where the width of the line is matched to the string core width.Note that due to the periodic nature of the simulation box, some strings seem to end suddenly but are in reality continuing on the other side of the volume.Animations available here.
Axion strings may also develop in axion-like particle models; for ultralight axion masses m a below that of the QCD axion, the resulting string networks can persist near or below the epoch of big bang nucleosynthesis (BBN), leading to a number of observable signatures.These include, for example, gravitational wave production, contribution to the effective number of neutrino degrees of freedom N eff , and the polarization rotation of cosmic microwave background (CMB) photons.The string network persists until H ∼ m a , with H the Hubble parameter.If the axion domain wall number N dw is larger than unity then stable domain walls develop for H ≲ m a , which could themselves lead to novel cosmological signatures.If N dw = 1 then the string-domain-wall network collapses around H ∼ m a .In this work, we restrict to axion masses low enough at a given cosmological epoch that we do not have to consider domain walls or specifiy N dw , for simplicity.
It has been shown that axion strings with axion masses m a roughly less than 10 −18 eV, with decay constants f a ≳ 10 14 GeV, could produce detectable gravitational wave signatures at next-generation observatories [26] (see also [27][28][29][30]).The gravitational waves are sourced by the energy density in the evolving string network.If the network persists until after the CMB decouples (m a ≲ 10 −29 eV), then axion strings may rotate the polarization of CMB photons by an amount that is proportional to the electromagnetic anomaly coefficient; polarization observations of the CMB already constrain this scenario [31][32][33][34][35][36].Interestingly, the polarization signatures are independent of f a and only depend on the electromagnetic anomaly coefficient of the axion.
In this work, we point out that axion-like particle strings, which we refer to as axion strings for simplicity, also leave novel energy injection signatures in the primordial plasma.It is well established that the axion-string network evolves by emitting axions to maintain the scaling solution.We show, however, that the string network releases a subdominant fraction of its energy into heavy radial modes (see also [37,38]) that promptly decay to Standard Model (SM) final states and, under certain circumstances that we enumerate, high-energy SM Higgs particles.This high-energy SM radiation is absorbed by the primordial plasma and can undo the success of BBN, modify the CMB power spectrum, or even lead to observable gamma-ray signatures today.Understanding how to translate cosmological measurements of these observables to constraints on or evidence for axion strings requires a detailed understanding of how axion strings emit heavy radial modes.
While this work focuses on global strings, local string networks are closely related and result from gauging the U (1) PQ symmetry that gives rise to the axion as a Goldstone mode.The axion is then "eaten" by the abelian gauge field, which acquires a mass of order the scale of symmetry breaking.Thus, unlike for global strings local strings do not emit massless radiation, except gravitational wave emission and perhaps other massless radiation that they may couple to indirectly.High-energy SM radiation from cosmic string networks has been studied before in the context of local strings (see, e.g., [39][40][41][42][43][44][45][46]).It is debated whether or not local strings can directly emit heavy modes with masses of order the symmetry breaking scale (see, e.g., [39,40,47,48]), but if they can then they would be constrained through analogous probes to those studied in this work.Higgs condensates surrounding local strings have also been studied previously [42,43,46], in analogy with the Higgs configurations around global strings studied in this work.
The remainder of this article is organized as follows.In Sec.II we review constraints on and future probes of axion-like particle strings from the axion contribution to N eff , gravitational waves, CMB distortions, and CMB polarization rotation.In Sec.III, through analytic arguments and AMR simulations, we compute the amount of heavy radial mode emission from an axion string network.The emitted radial modes can, in turn, decay into SM degrees of freedom and inject energy into the primordial plasma.To compute this effect, in Sec.IV, we derive the branching ratios of the radial mode into SM particles.Among such SM final states, the radial mode can also decay into a pair of SM Higgs via a PQ scalar-Higgs quartic coupling.In particular, such a coupling would give rise to non-winding classical Higgs configurations surrounding axion strings.In Sec.V we compute the properties of such Higgs 'sheaths' analytically and using AMR simulations.Using these results, in Sec.VI we obtain constraints on axion strings from the energy injection signatures they would leave in the primordial plasma at the epochs of BBN and CMB decoupling, in addition to constraints arising from present-day gammaray surveys.We conclude in Sec.VII with a discussion of how some of our results may apply to local strings.

II. EXISTING PROBES OF AXION STRINGS
We describe some of the current and future probes of axion strings that have been previously discussed in order of decreasing axion mass m a , as the constraints typically become stronger as one allows the string network to persist to later times.Note that throughout this work we consider only the so-called "field theory axions," which are those that emerge as the pseudo-Goldstone bosons of spontaneously broken U (1) PQ symmetries.Axionlike particles are also motivated by the framework of the string axiverse [49,50], where the axions arise not from global symmetry breaking but rather as the zero modes of higher-dimensional gauge fields integrated over the compact manifolds in string theory compactifications [49,[51][52][53][54][55][56].Axion string production for string theory axions is more subtle than for field theory axions and is not considered here; we discuss the formation and signatures of string theory axion strings in [57] (see also [58]).
To produce axion strings in field theory realizations, the Universe must reheat after inflation to a temperature T RH > f a , large enough to restore the PQ symmetry.This is because finite-temperature corrections to the effective potential for the PQ field restore the PQ symmetry for temperatures T ≳ √ 3f a , where the numerical pre-factor is somewhat model dependent (see, e.g., [59,60]).The maximum reheat temperature, given a Hubble parameter during inflation H I , is determined by hypothesizing that the inflaton decays promptly at the end of inflation to the SM, such that the energy density after reheating ρ RH ∼ T 4  RH is equal to the energy density directly before reheating, which is roughly 3M 2 pl H 2 I , with M pl the reduced Planck mass.If the inflaton does not decay promptly then the reheat temperature may be lower.The scale of Hubble during inflation is constrained by the tensor-to-scalar ratio as measured by CMB anisotropy measurements.A combination of Planck data and data from the BICEP2/Keck Array constrain H I ≲ 5 × 10 13 GeV at 95% confidence level [61], which implies T RH ≲ 10 16 GeV.Thus, we conclude that f a ≲ 6 × 10 15 GeV, regardless of the axion mass m a , in order to produce axion strings.
The string network is relatively unconstrained at present until BBN.However, this may change in the future with the next-generation gravitational wave observatories.Axion strings evolve primarily through axion emission.As we discuss in this work, a sub-dominant fraction of the energy density is also emitted in the form of heavy radial modes, though the energy density in this fraction is suppressed relative to that in axions by an amount ∼log(f a /H) ∼ 10 2 , where H is evaluated at late times, such as during BBN or CMB decoupling.On the other hand, the axion string network also sources gravitational waves but with a rate heavily suppressed relative to axion emission by an amount ∝ (f a /M pl ) 2 .Gravitational wave emission thus has a negligible effect on the dynamics of axion strings for the decay constants of interest (f a ≲ 10 15 GeV).However, the spectrum of gravitational waves emitted by the string network may be detectable through next-generation lowfrequency gravitational wave observatories such as LISA and the square kilometer array (SKA) with pulsar timing; in particular, Ref. [26] concluded that string networks with f a ≳ few × 10 14 GeV for m a ≲ 10 −18 eV may be detectable.Note that the recently-detected gravitational wave signal at pulsar timing array experiments, including EPTA and NANOGrav, is not compatible with axion string emission because of the N eff bound discussed below [62][63][64][65][66].
At BBN, existing constraints on axion strings arise from the contribution of relativistic axions to N eff (see, e.g., [26]).The string network evolves by emitting relativistic axions; the energy density in axions is [20] , where log * ≡ log(m s /H), with H the Hubble parameter at the epoch of interest and where c 1 is an O(1) coefficient that we discuss later in this article.We may divide this energy density by the energy density in one species of neutrino to compute the contribution to N eff during radiation domination, in particular, at 1 MeV: with g * the effective number of degrees of freedom.At BBN, ∆N eff is constrained at 95% to be less than 0.46, which implies f a ≲ 7×10 14 GeV for m a ≲ 10 −18 eV.Note that it is more appropriate to account for the change in g * around the epoch of BBN, and a more detailed calculation incorporating this effect can be found in [26], which finds an approximately similar upper bound f a ≲ 9 × 10 14 GeV.Note that the BBN constraints leave a narrow range of decay constants that may be detectable with gravitational wave observations.In this work, we further constrain and narrow the parameter space that may be detectable in gravitational waves through radial mode emission.
If the string network persists to the epoch of CMB decoupling then additional probes arise from the gravitational effects of the string network imprinted on CMB anisotropies [67][68][69].No analyses of these anisotropies have been performed to-date that account for the scaling violation of the string network.However, extrapolating from existing results from global strings [70], to account for the larger values of ξ expected at the CMB decoupling epoch, suggest that f a ≲ 2 × 10 14 GeV are likely in tension with CMB anisotropy data [26].Axions emitted from strings contribute to ∆N eff during the CMB epoch as well, and the bound ∆N eff < 0.34 [71] leads to f a ≲ 3 × 10 14 GeV, with log * evaluated at recombination.Note that the CMB limit is stronger than that from BBN in part because of the larger log * value at the CMB epoch and in part because for the CMB N eff is measured in matter domination, at z ∼ 1100, while matter-radiation equality is at z ∼ 3400.
Axion strings persisting until after recombination (m a ≪ 10 −29 eV) may also rotate the polarization of CMB photons [31][32][33][34]36].This effect involves the axionelectromagnetic interaction, which may be parameterized by where F is the electromagnetic field strength (with F its dual), α EM the fine-structure constant, and A the mixed PQ-electromagnetic anomaly coefficient.If a photon propagates along a trajectory over which the background axion field changes by an amount ∆a, then the polarization angle of the photon will change by ∆Φ = Aα EM /(2πf a )∆a.Given that complete loops enclosing strings cores are characterized by ∆a = 2πf a field excursions, we expect that photons propagating through a background of cosmic axion strings will undergo polarization angle rotations on the order of ∆Φ ∼ Aα EM .Ref. [33] performed a dedicated search for this effect in CMB polarization data as measured by Planck and set the upper limit Aξ 0 ≲ 0.93, with ξ 0 ∼ 10 the expected number of strings per Hubble patch at recombination, though in detail this upper limit likely depends on the morphology of the string network assumed in [33].Interestingly, the polarization probes are independent of f a and close to probing theoretically motivated parameter space for which A ≳ 0.1 [72].In contrast, the signatures that we develop in this work are independent of A but directly probe f a .

III. PROMPT RADIAL MODE EMISSION FROM AXION STRINGS
We now turn to our calculation of the SM radiation produced by evolving axion strings.We start by considering SM particles generated by the decay of heavy radial modes that are generated by the evolving strings.First, we review the standard picture for axion string evolution in the context of a field theory UV completion with a PQ complex scalar field that undergoes spontaneous symmetry breaking.

A. PQ axion strings
Let us first recall the dynamics of axion strings in the scaling regime, at temperatures T ≪ f a but H(t) ≫ m a (see [1,20] for modern reviews).Note that the dynamics we describe here are valid both for axion-like particle strings and for the QCD axion string network at temperatures above the QCD phase transition.This is because the QCD axion mass is temperature dependent in the early Universe; it only becomes relevant at the QCD epoch, where it rises rapidly and quickly exceeds the Hubble parameter.Thus, at temperatures well above a GeV the dynamics of the QCD axion string are the same as those of axion-like particle strings with H(t) ≫ m a .Furthermore, in this limit, the domain wall number also does not affect the dynamics, since the domain wall number only plays a role in the evolution when the axion potential becomes important.
In the standard field theory UV completion the axion a arises as the (pseudo-) Goldstone boson of global PQ symmetry breaking of a complex scalar field Φ, which for T ≪ f a we represent as with s the radial mode. 1 In vacuum, the field Φ is subject to the Lagrangian In the thermal Universe the field Φ also has a thermal mass term V therm (Φ) = m 2 therm |Φ| 2 , with m 2 therm ≈ λ Φ T 2 /3 (see, e.g., [60]).The thermal potential restores the PQ symmetry at T ≫ f a .After PQ symmetry breaking the radial mode acquires a mass m s = √ 2λ Φ f a while the axion is massless.
After PQ symmetry breaking the radial mode abundance redshifts like matter until it decays to lighter states everywhere except at the location of the axion string cores.Axion strings are topologically protected solutions to the PQ equations of motion for which the axion acquires a 2πf a phase shift when traversing a path that encloses a string.The radial mode, which takes values s ≈ 0 far from the strings, reaches the value s = −f a , so that Φ = 0, at the string core.As we discuss more below, semi-analytic solutions are available for infinitely straight strings (see, e.g., [1] for a review), though in the cosmological context, the strings move, bend, combine, and disappear dynamically, thus requiring numerical simulations (see [20,23,25] for the current cutting-edge simulations).The string network primarily loses energy by radiating axions, though -importantly for this work - 1 For simplicity we assume a KSVZ-type UV completion [73,74].
a small fraction of the dissipated energy goes into heavy radial modes.
The axion strings have tension µ eff ≈ πf 2 a log(m s /H), to leading order in large log(m s /H).The log arises because the axion configuration that surrounds the string has an energy density that falls off slowly with distance, leading to a logarithmic divergence in the string tension; the large-distance cut-off is ∼H −1 , which is approximately the distance to the nearest string in the scaling regime.The energy density in the string network is ρ s ≈ 4ξH 2 µ eff , where ξ is the average number of strings per Hubble patch.Note that ξ is formally defined by ξ ≡ ℓt 2 /V at time t, where ℓ is the string length within a large volume V.In the scaling regime, ξ is approximately constant, regardless of whether the Universe is, e.g., matter or radiation dominated.Logarithmic derivations to the scaling solution are now understood to arise [23,25], with ξ ≈ c 1 log(m s /H) at large log, with c 1 ≈ 0.25 in the radiation dominated epoch [25].In App.D we show that c 1 ≈ 0.06 in matter domination.
At large log(m s /H) the rate of axion production is (e.g., [20]) Recent simulations suggest that the momentum-space distribution of radiated modes is nearly conformal: ∂Γ a /∂k ∝ (H/k) q , with q = 1.02 ± 0.03, for 1 ≪ k/H ≪ m s /H [25] with k being the momentum of radiated axion.

B. Radial mode emission: general expectations
In addition to radiating axions, the strings may also radiate radial modes.In the next subsection, we compute the radial mode emission rate by performing dedicated simulations, but in this section, we discuss our general expectations for this emission rate.At momenta k ≳ m s we do not expect that the PQ theory differentiates between axion emission and radial mode emission, since at these high energies the radial mode emission is mildly relativistic.On the other hand, radial mode emission is disallowed at frequencies less than m s .Since dΓ a /dk ∝ 1/k, as observed in simulations for H ≪ k ≪ m s , we conjecture that Γ is that with k ≲ m s .Since for highk the theory should not differentiate axion versus radial mode emission, we then conjecture the radial-mode emission rate Γ s ≈ Γ k≳m S a : to leading order in large log(m s /H), where c is an undetermined constant expected to be of order unity.A deeper understanding of the relation between ( 6) and ( 7) is found through the distribution of string loops.As discussed in [25], we may understand the spectrum dΓ a /dk ∼ 1/k observed in simulations through the observation, in the same simulations, that ℓdn ℓ /dℓ ≈ const.Here, n ℓ represents the number density of string loops with lengths less than ℓ at any given time.A loop of length ℓ radiates axions at a characteristic wavelength k ∼ 1/ℓ. 2 Moreover, it has been shown that string loops (and also kinks in long strings) radiate energy at a constant rate dE/dt regardless of the loop (kink) size [11,15,75,76].This implies a conformal spectrum of axion emission dΓ a /dk ∼ 1/k, with the high-k modes being emitted by loops and kinks with large curvature; if ℓ ≲ 1/m s , then these loops and kinks also emit radial modes with similar efficiency to axions.Indeed, in the numerical simulations described below, we find that radial modes are dominantly produced in regions of large string curvature.

Analytic estimate for radial mode emission rate from string tension
In the following subsections, we simulate the stringnetwork evolution to measure the constant c.First, however, we present a rough but insightful analytic argument that suggests a value c ∼ 0.1, which we later verify and refine with the numerical simulations.The general idea behind our approach below is that as the string network evolves it radiates energy into axions and radial modes, but that energy must come from the stored energy in the string tension.Thus, it is plausible that the fraction of radiated energy with k ≳ m s will be proportional to the fraction of the string tension that, in Fourier space, also has k ≳ m s .For k ≳ m s the axion and the radial mode are not qualitatively different, so we hypothesize that the modes emitted with k ≳ m s are split democratically between axions and radial modes.Below, we make this argument precise.
We parameterize the PQ profile describing a long string oriented along the z axis as where g(m s r) is a dimensionless function and {r, θ, z} describe a 3D cylindrical coordinate system.For a string solution, g ∼ m s r for small m s r, and g ∼ 1−1/(m s r) 2 for large m s r.To estimate the emission of radial modes, we can Fourier transform the gradient of the position space profile to see which momentum modes are supported with k ≳ m s .The position space expression for the energy 2 More precisely, axion string loops of size ℓ appear to radiate axions with a conformal instantaneous spectrum ∝ 1/k for 1 ℓ ≲ k ≲ ms [38]; in App.E we perform numerical simulations of collapsing string loops to verify this scaling and discuss in more detail how it leads to a conformal emission spectrum for the network as a whole.density is given by, where x = m s r and ′ denotes a derivative with respect to x.The string tension µ is then given by µ = dθdr rρ str .
Restricting to a given plane orthogonal to the infinite string, the 2D Fourier transform of the energy density is where J 0 is the 0 th Bessel function of the first kind.At small k we have ρ(k) ∝ log(k), which gives rise to the IR divergence in the string tension µ = ρ(0).Physically this divergence is cut off by Hubble providing the largest, relevant length scale at kIR ∼ H/m s .(Note that we define µ eff to be the effective tension computed with the IR cut-off k IR .)At large k, ρ(k) falls exponentially, as in Fig. 2. We estimate the part of the string tension relevant for axion emission as that with k > m s ; we compute this contribution as µ UV ≡ ρ(m s ) = 2πc UV f 2 a , where the constant c UV is defined with the specific normalization because we show in the following paragraph that it is related to the constant c in (7).Numerically, we find c UV ≈ 0.16.
Let us now discuss the relation between µ UV and Γ s .By comparing the evolution of the energy density of the string network in the scaling solution to that of the free-string network one may infer that the string network must emit energy with rate Γ tot = 8ξµH 3 [20].Let us assume that the modes with k > m s emitted from the strings are split equally between axions and radial modes; while this is almost certainly not completely true, it allows us to make an O(1) estimate for the radial mode emission rate.Then, we estimate that Γ s ≈ 8c UV H 3 ξπf 2 a .Comparing with (7) we thus estimate that c ≈ c UV ≈ 0.16.As we show in the following subsections, this estimate for c is similar to that we find in dedicated numerical simulations of the string network.

C. PQ simulations for radial mode radiation: setup
We simulate the evolution of axion strings with the Lagrangian as in (5) along with the thermal mass term for Φ.We fix λ Φ = 1 throughout this work for definiteness, in all of our simulations, such that the radial-mode mass is m s = √ 2f a .We follow the basic procedure outlined in [25] and the equations of motion, along with common technical details pertaining to the base code, can be found therein.Ref. [25] performed simulations of the axion string network to measure the axion radiation and infer the axion mass that gives rise to the observed DM abundance for the QCD axion.For our purposes, however, we are interested in the ratio of emission rates of radial-mode emission relative to axion emission.
Our code is based on the AMReX framework [77] and is capable of AMR.AMR allows for a dynamical grid with multiple refinement levels that track given spatial locations at higher spatial and temporal resolution (see Fig. 1 for an illustration of an AMR grid).While our setup is largely identical to that of [25], our different objective forces two major changes.First, because the radial mode spectrum is expected to peak at short wavelengths which propagate throughout the entire simulation volume, we do not employ any of the AMR capabilities (i.e., we use a single refinement level).This is because the AMR setup would miss some of the short-wavelength radiation of the radial mode since its amplitude is often too small compared to that of axion radiation to trigger our refinement criteria.(Note, however, that we do use multiple AMR refinement levels in our simulations including the Higgs field in Sec.V.) Instead, the simulation is performed on a static lattice containing initially 512 3 grid sites to ensure that the radial mode is resolved as well.As the co-moving width of axion strings decreases over time, the number of grid sites is increased by a factor of 2 3 through quartic interpolation every time the string width would be resolved by less than four grid sites, up to a final size of 4096 3 grid sites.
The second change we make relative to the simulation setup in [25] involves our initial conditions.Ref. [25] started with thermal initial conditions at η i for the complex scalar field Φ.Here, since we are unable to evolve to as large of log(m s /H) values as in [25] -since we are not using multiple refinement levels -we pre-evolve the state following the procedure in [20].The pre-evolution mitigates the effects of transient oscillations to the string cores that are relevant at small values of log(m s /H) and which would otherwise contaminate our emission measurement.The pre-evolution procedure is described in more detail in App.B and involves modified equations of motion that support a constant string width and a moderate amount of Hubble friction.The pre-evolved state is designed to be close to the attractor solution at our starting time log(m s /H) = 2.To avoid the reintroduction of transient radial-mode excitations when starting from a pre-evolved state due to the sudden change in the underlying physics we introduce a short adiabatic regime between log(m s /H) = 2 and log(m s /H) ∼ 3.In this regime, we smoothly interpolate between the two sets of equations of motion using a logistic function.(See App.B for more details.) The simulation volume is a periodic box with comoving side length L = 33/(R 1 H 1 ), where R 1 is the scale factor of the Friedmann-Robertson-Walker (FRW) metric at a reference time t 1 such that the Hubble parameter is H(t 1 ) = f a ≡ H 1 .Our simulation is evolved in conformal time η defined by R(t)/R(t 1 ) ≡ η.The simulation begins at η i ≈ 2.3 (log(m s /H) = 2), and we evolve until η f ≈ 21.7 (log(m s /H) ≈ 6.5).We use a low-storage strong-stability preserving third-order Runge-Kutta algorithm to advance the field with a time step size given by the Courant-Friedrichs-Lewy condition through ∆η = 0.33∆x.The Laplacian operator is computed through a finite-difference 7-point stencil.At the end of the simulation, our box contains approximately 3.5 Hubble volumes.
The simulations are performed at the NERSC Perlmutter supercomputer.Each simulation runs for approximately 30 minutes on 256 AMD EPYC 7763 CPUs (i.e., 128 total nodes, 256 total CPUs, and 16,384 total CPU cores).

D. PQ simulations for radial mode radiation: results
We are interested in measuring the amount of axion and radial mode radiation that is emitted from the string network over time.In Fig. 3 we illustrate the radial mode and axion radiation from a snapshot near the end of the simulation, at log(m s /H) ≈ 6.5.In the top panel, we show the time-derivative of the radial mode squared ( ṡ2 ), in logarithmic units, which is a proxy for the radial mode energy density.The string network is clearly visible in the top right panel; the bright regions away from the strings are regions of significant axion radiation.The zoom-in on the top left panel shows a string region producing large amounts of radial mode emission.That region of the string is characterized by its high curvature, which suggests that radial modes are predominantly produced from regions of the strings with high curvature of order  the radial mode mass itself (see also App.E).In contrast, the lower panel shows the axion time derivative squared ( ȧ2 ) for the same state as in the left panel.The axion radiation has support at longer wavelengths relative to radial mode radiation.Thus while the high-curvature region also produces significant axion radiation, the contrast versus the rest of the string regions is not as large.
To compute the energy densities more precisely we use the fact that away from the string cores both the axions and radial modes are free fields.At a given point x the energy density of a real, free scalar field X, which solves its classical equations of motion, is where m X is the field's mass and where we have applied the equation of motion to arrive at the second line.This implies that we can compute the average energy density over the simulation box, ρ X ≡ 1 4.0 4.5 5.0 5.5 6.0 6.5 log(m s /H) We may take X(x) to be either the axion a(x) or the radial mode s(x).X(k) is the Fourier transform of Ẋ(x) as extracted from the simulation, where Ẋ(x) is screened to avoid contributions from the strings itself [20,25].Explicitly, we use Ẋ(x) → Ẋ(x)(s(x)/f a + 1) 2 , since (s/f a + 1) 2 ≈ 1 away from strings but ≈ 0 near the core.The exact form of this screening has little effect on the result [25].In Fig. 4 we show the energy density as a function of time, displayed as log(m s /H), for both the axion and the radial mode.The energy density emitted by the string network is dominated by axion radiation; to compute the emission rates we need to take the appropriate time derivatives: Here, z characterizes how the average energy density of the field X red-shifts, with z = 3 for non-relativistic modes and z = 4 for relativistic modes.The axions have z = 4 since they are massless.The redshift factor z for the radial mode can be computed by [20,23,26] with and In Fig. 5 we show our results for z computed over time from the simulation output for the radial mode.In general, we find z ≈ 3.7, though it appears that z slightly .The quantity z computed for the radial mode, defined in ( 14) and ( 15), describes how the instantaneous emission red-shifts at production.Completely non-relativistic (relativistic) radiation has z = 3 (z = 4), with the free particles scaling with the scale factor like ρs ∝ R −z .The emitted spectrum of radial modes is semi-relativistic, giving the intermediate z shown.
decreases over time.Note that z can be no smaller than z = 3 and physically we expect z > 3 since the radial modes are produced semi-relativistically.The small decrease in z over time may be expected since as we go to larger log(m s /H) values we have a larger dynamical range between Hubble, which provides an IR cut-off, and the UV cut-off provided by m s ; this may account for an increased weight for the low-k part of the spectrum at later times.In our fiducial analysis below we take z as measured at each log(m s /H) step, though we show that our results are robust to changes in z.For example.we consider z as small as 3.3 and as large as 3.8.
In the left panel of Fig. 6 we show the instantaneous axion spectrum, divided by 8H 3 ξπf 2 a ; referring to (6) we expect Γ a /(8H 3 ξπf 2 a ) ≈ log(c a m s /H), where c a is a constant of order unity that accounts for finite contributions to the string tension and the precise form of the IR cutoff to the tension.We fit the model expectation to the data (for log(m s /H) ≥ 5) to get the best-fit curve shown in dashed black, which has c a ≈ 0.063 ± 0.002.(The grey band illustrated the 1σ band on c a .)We perform the fit to the Γ a data assuming that the data points for Γ a have Gaussian uncertainties σ = α Γa (c a ), where α is a hyperparameter and where Γa (c a ) is the model prediction for Γ a for a given choice of c a .To reduce statistical and systematic noise we perform a total of nine simulations with different initial states; see Appendix B for more details.The results of these simulations are similar and hence we combine them by averaging Γ a , Γ s , and ξ.We profile over α to compute the confidence interval for c a .The best-fit value for α at the best-fit value for c a is then used to construct the 1σ error bars on the data points shown in Fig. 6.
The quantity Γ a /(8H 3 ξπf 2 a ) is observed to rise linearly with log(m s /H), as expected.In contrast, we expect Γ a /(8H 3 ξπf 2 a ) = c, referring to (7), to be constant with log(m s /H).In the right panel of Fig. 6 we show the Γ s data from our fiducial choice of z, which is extracted directly from the simulation at each time-step.The best-fit value of c is shown in dashed black with the grey band indicating the 1σ confidence interval; we find c ≈ 0.33 ± 0.13.We use this value of c in our analyses below looking at the effect of radial-mode induced energy injection, though keep in mind that different UV completions may give slightly different values for c.
Our best-fit value of c has a mild dependence on the choice of z for the radial model emission.Choosing a constant z of 3.8 (3.5) (3.3) leads to a central value for c of ∼0.38 (0.25) (0.17).We do not consider this source of uncertainty further since it is subdominant compared to other sources of uncertainty in our analyses.

IV. RADIAL MODE DECAYS
While the semi-static axion strings are protected from decay by the topology of the configuration, the radiated radial modes, with production rate given in (7), will promptly decay to lighter states.In particular, the radial modes may decay to pairs of axions or to pairs of SM states.The amount of visible energy injected by the string network is proportional to the branching ratio with Γ s→SM (Γ s→aa ) the decay rate to SM final states (axion final states).This quantity depends on the UV completion of the theory.We consider KSVZ-type [73,74] axion-like-particle scenarios for illustration.(Note that the simplest implementations of the DFSZ [78,79] scenario give a coupling of the axion to gluons, which we want to avoid since we are considering axion-like particles and not the QCD axion.) In the KSVZ-type scenario, the Lagrangian is given by Y but is a singlet under SU (3) c , since we are focusing on axion-like-particles and not the QCD axion.After spontaneous symmetry breaking the vector-like fermion acquires a mass The vector-like quarks must be able to decay as otherwise the population produced in the early Universe, both thermally prior to the PQ phase transition and from the evolution of the string network, would become non-relativistic at later times and over-close the Universe.New interactions beyond those in (19) are needed for the vector-like quarks to decay (see [80] for a review).In this regard, for specific choices of SU (2) L ×U (1) Y charges for the vector-like quarks, additional interactions with the SM are possible.For example, dimension four operators of the form f QH, where f are SM fermions, Q are the vector-like quarks, and H is the SM Higgs, may lead to heavy quark decay.Dimension-five operators may also be responsible for such decays [80].For our purposes, the specific forms of such operators are not important; all that is required is that if a heavy quark is produced it will eventually decay to SM final states.
Without loss of generality let us assume that under a U (1) PQ transformation Φ → e iα Φ, Q L → e iα Q L , and Q R → Q R for constant α.This transformation leaves the Lagrangian in (19) invariant and allows us to construct operators that respect the PQ symmetry but that induce PQ fermion decay using Q R .
As an illustration, let us consider the possible dimension-four operators that give rise to heavy quark decay.We chose the convention for the weak-isospin such that U (1) EM is generated by Q = Y + T 3 , with T 3 = 1 2 σ 3 the third generator of SU (2) L and σ 3 the third Pauli matrix; the SM Higgs H then has weak hypercharge Y = 1 2 , right-handed leptons have Y = 1, and left-handed leptons have For the radial modes to decay to KSVZ fermions far away from the string cores we need 2m Q < m s , which implies y Q < √ λ Φ .If this inequality is not satisfied then radial modes that propagate sufficiently far from the string cores will not be able to decay to KSVZ fermions and must decay directly to axions or to SM final states, as we discuss below.On the other hand, near the string cores, the KSVZ fermion masses are greatly reduced and indeed the fermion masses vanish at the core centers.Thus, the radial modes are always able to decay to Q pairs near the string cores; these KSVZ fermions are then kinematically trapped near the string cores until they decay to SM final states.
We now consider several scenarios for the decays of the radial mode.Since the radial mode always decays to a pair of axions, we consider one SM channel at a time to compute the branching ratio into the SM.The expressions given below can be appropriately modified where N L = 1 if the fermions are SU (2) L singlets and N L = 2 if they are doublets.This decay channel contributes to Γ s→SM , since the KSVZ fermions decay completely to SM final states.On the other hand, because of the kinetic term, the radial mode may decay to axions at tree level, such that Thus, the branching ratio of the radial mode to SM final states is For N L = 1 (N L = 2) the branching ratio may be as large as B ≈ 0.16 (B ≈ 0.27).On the other hand, if the Yukawa coupling to KSVZ fermions is small relative to the PQ self-coupling, then the branching ratio is suppressed; for example, B ≈ 2 × 10 −4 for y = 10 −2 and N L = 2.
b. Decay into electroweak gauge bosons.In addition to the tree-level decays of the radial mode to KSVZ fermions and axions, there are one-loop decays to SM gauge bosons that may also be relevant.Taking the vector-like quarks to be in the fundamental representation of SU (2) L (N L = 2), for example, the decay rate to W -bosons is given by Assuming N L = 2 we thus find B ≈ 2 × 10 −5 , where we use the value of the weak fine structure constant α 2 (m s ≈ 10 16 GeV) ≈ 0.02 at energy scales of order the grand unified theory (GUT) scale.The field s can decay also to ZZ, γγ, and γZ, though these processes are roughly a factor of five smaller in total than the decay rate to W W . See App.A for more details.
c. Decay into SM Higgs.The UV Lagrangian may also contain the renormalizable terms connecting the PQ field with the Higgs field through the potential: Here, H is the SM Higgs doublet, λ H is the Higgs quartic, and −μ 2 H is related to the Higgs mass parameter.The stability of the electroweak-symmetry-breaking vacuum requires With this requirement, one can see by integrating out the radial mode that the Higgs acquires a VEV Expanding Φ about its VEV we see that the radial acquires an interaction with the Higgs field, which allows the radial mode to decay at tree-level to Higgs pairs.Since the radial mode mass is well above the Higgs mass, the decay rate s to Higgs pairs is given by HΦ )] and B < 1/[1 + λ Φ /(8λ H )]. As we discuss further below, the Higgs quartic famously runs to small and potentially even negative values at high energy scales.On the other hand, λ H does receive threshold corrections near the PQ scale that push it further positive in the UV (see, e.g., [81]).For definiteness, let us assume that λ H = 0.01 at the PQ scale, such that B ≲ 0.1, though smaller, positive values for λ H do not qualitatively change the branching ratio.(For example, B ≲ 0.01 if λ H = 10 −3 at the PQ scale.)On the other hand, if λ H is negative at the PQ scale this leads to runaway behavior for the axion strings, and so we do not consider that possibility further.
In summary, depending on the Lagrangian parameters the branching ratio to SM final states may be as large as B ∼ 0.3 or as small as B ∼ 10 −5 .On the other hand, there are many ways of achieving a large branching ratio: for example, if the radial mode is kinematically allowed to decay to PQ fermions then it will generically do so with a large branching ratio unless the Yukawa coupling to the fermions is small.The branching ratio B ∼ 10 −5 is irreducible if we insist on the axion coupling to Wbosons; it can be a factor of a few smaller still if all other s-decay channels to SM final states are removed and the axion only couples to hypercharge and not to SU (2) L .
Lastly, we note that even if y Q > λ Q , such that the radial mode is not kinematically allowed to decay to KSVZ fermions asymptotically far away from the string, the radial mode could still decay to fermions close to the string.In the presence of a background field Φ the KSVZ fermions have a mass ∼ y Q |Φ|; asymptotically far away from the strings Φ goes to its VEV f a / √ 2. However, at the string core Φ = 0, with the magnitude of Φ rising to f a over a distance of order m −1 s from the string core.However, we do not consider such decays further here and rather bracket the possible branching ratio to SM final states as being within the range B ∈ (10 −5 , 0.3).

V. AXION-HIGGS STRINGS
The quartic coupling between the PQ scalar and the SM Higgs field -the term parameterized by λ HΦ in ( 26) -is generically present in KSVZ models.This term could be present in the UV but otherwise, it is generated under the renormalization group through the KSVZ fermions at lower energy scales.We show in this section that the presence of this coupling leads to non-trivial, classical Higgs field profiles surrounding the strings.Below, we refer to these strings as axion-Higgs strings and to the Higgs profiles as "Higgs sheaths."These Higgs sheaths may have a number of important implications, but for our purposes, they provide efficient sources of SM radiation from the cosmologically-evolving axion-Higgs string network.In this section, we study the axion-Higgs strings and perform dedicated simulations to study their dynamics and the radiated Higgs fields that are shed during their cosmological evolution.
Axion-Higgs strings were studied previously in [82] in the context of the DFSZ model, where they were called electroweak axion strings.However, the DFSZ electroweak axion strings are fundamentally different from the axion-Higgs strings that we discuss in this section.In particular, the Higgs sheaths we discuss have no nontrivial winding around the string cores, while the electroweak axion strings exhibit non-trivial winding in the context of the two-Higgs doublet model (2HDM).In the DFSZ model, there are two Higgs doublets H 1 and H 2 that have non-trivial interactions with the complex PQ scalar Φ, including those of the form where κ is a coupling constant.If we let the PQ charge of Φ be unity and the PQ charges of H 1 and H 2 be X 1 and X 2 , respectively, then PQ invariance requires 2 − X 1 − X 2 = 0. Axion string solutions in this scenario may be of the form, for infinite, straight strings: Φ = f a e iθ g(r)/ √ 2, H 1 = v 1 e iθ (0, h 1 (r)) T , H 2 = v 2 e iθ (0, h 2 (r)) T [82].However, while the tension associated with the Φ profile is ∝ f 2 a , the tension associated with the H 1 and H 2 profiles is significantly smaller, of order the Higgs field VEVs v 2 1,2 .Given this reason, we do not consider the DFSZ scenario in further detail in this work. 4n contrast to the electroweak axion strings scenario discussed in [82], the Higgs sheath solutions that we present here contribute to the tension at order f 2 a .The Higgs sheath solutions appear in KSVZ-type scenarios, with single Higgs multiplets, with the Lagrangian as in (19) and (26).As we discuss below, the λ ΦH term generates non-trivial solutions for the Higgs field exterior to the strings with no winding.
Below, we first discuss the Higgs sheath profiles for infinite, straight strings, and then we confirm the semianalytic expectations for the Higgs profiles using numerical simulations of the axion-Higgs cosmology.
A. Semi-analytic solutions for infinitely straight axion-Higgs strings Consider an infinitely straight string in the ẑ direction in cylindrical coordinates (r, θ, z), in the theory consisting of the PQ scalar Φ and a single complex scalar field H, with potential given in (19) and (26).In reality, the Higgs is an SU (2) doublet, but by gauge symmetry, it is sufficient to work with a singlet Higgs field when computing the contribution to the string tension.
In cylindrical coordinates, the ansatz for an infinite, straight axion string along ẑ is for some function g(r).We hypothesize that in the presence of the λ HΦ interaction the Higgs field acquires a non-trivial profile that we may write as for a real function h(r).Then, in Minkowski space, a static solution (g, h) obeys the equations of motion and Here primes denote derivatives with respect to r. Writing the total energy density associated with the string solution as ρ tot = ρ H + ρ Φ + ρ int , the sub-components of the energy density associated with the different fields are The tension is defined as the total energy per unit length of the string Note that unlike for the PQ field Φ, which has non-trivial winding around the string, there is no topological protection for the Higgs profile.The stability of the Higgs profile can rather be determined by the following consideration.A trivial solution to the equation of motion (33) is given by h = 0 everywhere.We denote the tension corresponding to that solution as µ h=0 . 5If µ tot < µ h=0 , with µ tot including the non-trivial Higgs field profile, then the Higgs sheath is stable.As already discussed, stability of the PQ and electroweak vacua requires 0 < λ HΦ < √ 4λ Φ λ H , and λ H > 0. The PQ and Higgs VEVs are then We use this notation below in describing the behavior of the field profiles.In practice, we may approximate With the approximation v EW ≈ 0, and assuming all of the dimensionless coupling constants (λ H , λ HΦ , λ Φ ) are order unity, there is only one dimensionful scale, which is f a .At small r, the equations of motion enforce g(0) = 0 in addition to h ′ (0) = 0.The quantity h(0) cannot be computed analytically; it must be computed numerically.However, the fact that the equations of motion do not depend on any dimensionful parameters 5 More precisely, we define µ h=0 as the energy density per unit string length of the string solution with h = 0 fixed minus the energy density per unit string length of the Universe with h = 0 and g = 1 everywhere.ensures that h(0) ∼ 1, such that H(0) ∼ f a .Asymptotically far from the string core h → 0 and g → 1.The fields Φ and H thus have field excursions of order f a over distance scales of order f −1 a , since f a is the only dimensionful parameter in the problem.This implies that the contributions to the tension from both H and Φ are expected to be of order f 2 a .Note that both g and h approach their asymptotic values at r → ∞ through terms the fall of with r as 1/r 2 , which implies that both the radial mode and the Higgs field have IR-finite contributions to the tension, unlike the contribution from the axion field, which is logarithmically-divergent in the IR.
It is important to contrast the Higgs solution above with that found in the DFSZ electroweak axion string.In the latter case, the Higgs fields have non-trivial winding, which implies that regularity at r = 0 forces H 1 (0) = H 2 (0) = 0. On the other hand, H 1,2 (r = ∞) ∼ v 1,2 .Thus, the Higgs fields in the electroweak axion strings only have field excursions of order their VEVs v 1,2 .It is precisely because our Higgs field does not wind that it is able to have H(0) ∼ f a and thus contribute substantially to the string tension.
In Fig. 7 we illustrate the radial mode and Higgs field profiles found for the infinite, straight string by numerically integrating (33) and (34) using a 4 th order collocation method.We take the Higgs VEV to be zero since physically it is much less than f a , and we make the choices (λ Φ , λ H , λ HΦ ) = (1, 4, 2).Note that |Φ| = 0 at the string core, as required in order to remove the singularity associated with the axion winding, while at large distances from the string, the PQ mode asymptotes to The Higgs field has a non-zero value at the string core (|H|/f a ≈ 0.3), since it has trivial winding, and it asymptotes to zero infinitely far from the string.
Let us now verify that the Higgs sheaths are stable by computing the tension of these configurations and comparing to the solution with H = 0 everywhere.In Fig. 8 we illustrate µ tot −µ h=0 as a function of λ H for the choices (λ Φ , λ HΦ ) = (1, 2).Note that vacuum stability requires λ H > λ 2  HΦ /(4λ Φ ) = 1 in this case, which is indicated.For all λ H the differences in tension are negative, suggesting that the Higgs sheaths represent the energetically preferred solutions and are stable to decay.
A network of axion strings with Higgs sheaths may be expected to evolve by emitting classical axion, Higgs, and radial mode radiation.By energy conservation, the total rate of energy loss to radiation by the string network is equal to the time derivative of the energy difference between the evolving string network and the free string network (see, e.g., [20]).The axion emission rate should be mostly unaffected by the Higgs profiles since the Higgs profiles extend over a distance ∼m −1 s while the axion emission dominantly comes from longer wavelengths.The Higgs emission rate should then arise from the IR-finite part of the string tension µ ∼ πf 2 a .More precisely, assuming the Higgs emission arises from energy conservation associated with the IR-finite part of the string tension we expect where f (λ H , λ Φ , λ HΦ ) is a function of the dimensionless coupling constants.An analytic derivation of f (λ H , λ Φ , λ HΦ ) appears difficult, due to the non-linear nature of the equations of motion.In the following sub-section, we numerically calculate f for specific choices of the coupling constants by performing AMR simulations.

B. AMR simulations of axion-Higgs strings
We verify the development of axion-Higgs strings and the subsequent classical radiation of Higgs modes through AMR lattice simulations of the coupled equations of motion in the early Universe.The simulation setup is similar to that in Sec.III except that it uses an adaptive mesh; details are described in App. C. We use an AMR grid instead of a static lattice grid in order to access a larger dynamical range.Note that we simulate the coupled Higgs-PQ system on an adaptive lattice with comoving side length L = 46/(R 1 H 1 ) up to log(m s /H) = 7.2, with R 1 and H 1 as defined in Sec.III.We start the simulation with thermal initial conditions for both the Higgs field and the axion.Our fiducial choice for the couplings is λ Φ = 1, λ HΦ = 0.3, and λ H = 0.05, as discussed further below, though we also consider variations to this fiducial choice.A 2D projection of the simulation volume, for an example simulation, is shown in Fig. 9 at log(m s /H) ∼ 7 for both the Higgs radiation energy density and the axion radiation energy density.Higgs sheaths and Higgs emission regions are clearly visible surrounding axion strings.
The goal of the simulations is to measure the function f characterizing the Higgs emission rate that appears in (38).First, let us consider reasonable choices for the parameters λ H , λ Φ , λ HΦ .Without any prior for λ Φ , we simply take λ Φ = 1 for illustrative purposes.The value we take for λ H is more subtle since it is known that λ H runs to small and potentially even negative values at high energy scales, assuming no new heavy physics (see, e.g., [83]).We insist that λ H is positive for consistency.In fact, the PQ field itself should give a threshold correction to the Higgs quartic at energy scales of order f a , which is the relevant scale for considering the Lagrangian when solving the classical equations of motion at distances of order f −1 a from the string core.That threshold correction could push the Higgs quartic to values around λ H ∼ 5 × 10 −2 , depending on the PQ scale and on the top quark mass [84].For definiteness, let us then take λ H = 5 × 10 −2 . 6We may then vary λ HΦ from small values all the way to the vacuum stability limit.
For each simulation we compute Γ H through an analogous procedure to that used in Sec.III, restricting to log(m s /H) > 5.In particular, we extract Γ H analogously to Γ s and Γ a using (13) with X(x) = h(x) and (14) with z = 4.In Fig. 10 we illustrate the example data points for Γ H for our fiducial simulation, along with the resulting linear fit to extract f , as defined in (38).In Fig. 11 we show f (λ H , λ Φ , Λ HΦ ) = Γ H /(8H 3 ξπf 2 a ) as a function of λ HΦ and λ H for a sequence of five simulations.All other couplings are fixed to our fiducial choice.
Note that there are two distinct ways that the axion strings may create Higgs radiation: (i) the strings may radiate high-energy radial modes, which decay quantum mechanically to Higgs pairs with branching ratio B ≈ 2λ 2 HΦ /λ 2 Φ , in the limit λ HΦ ≪ λ Φ ; and (ii) the strings directly radiate classical Higgs radiation, with rate given in (38).Both contributions should be accounted for when computing the energy injection due to axion strings.In Fig. 11 we show, in addition to f (λ H , λ Φ , λ HΦ ), the radial-mode emission rate times the branching ratio B of the radial mode to Higgs particles, though for the purpose of illustration, we neglect the back-reaction of the Higgs field on the PQ field when calculating Γ s .That is, in Fig. 11 we use the PQ-only simulation results when computing Γ s .Still, this comparison suggests that while the direct emission of Higgs particles may dominate in certain regions of parameter space, the Higgs sheaths do not parametrically increase the energy-injection relative to what would naively be estimated based off of radial-mode emission alone.

VI. OBSERVATIONAL CONSTRAINTS ON AXION STRINGS FROM SM RADIATION
We now consider the observational constraints that arise from radial-mode-induced radiation into SM final states for string networks that survive until at least the epoch of BBN.Note that while we frame the discussion in terms of radial mode emission and decay, the following arguments also apply to high-energy Higgs emission from Higgs sheaths.A summary of all of the upper limits on f a derived in this section is provided in Tab.I.
As discussed in Sec.IV, radial modes decay to SM final states with model-dependent branching ratio B. The SM particles, which may be, e.g., Higgs boson pairs or heavy gauge bosons, then subsequently undergo a sequence of prompt decays to produce a spectrum of SM final states, with characteristic energy scale given by the radial mode mass m s .At the epochs of BBN and CMB decoupling those SM particles rapidly deposit their energy in the primordial plasma through scattering processes.Note that this also applies to neutrino final states.High-energy neutrinos (e.g., E ∼ 10 12 GeV) scatter off the primordial plasma at rates much faster than Hubble at the epochs of BBN and reionization (see [85], for example, for a discussion of the high-energy neutrino cross-sections).Thus, for the purpose of the following discussions we do not differentiate between neutrino and non-neutrino SM final states.
We derive constraints on the axion decay constant associated with axion strings at three different cosmological  epochs: (i) BBN, (ii) the dark ages between CMB decoupling and reionization, and (iii) today (redshift z = 0).All of these constraints arise from injecting additional energy into the Universe from the string network through radial mode decay into SM final states.At a given redshift z the energy injected per unit time per unit volume from radial mode decay is, referring to (7), The redshift dependence of ( 39) may be made more explicit by recalling that with Ω Λ (Ω m ) (Ω rad ) the present day relative abundances of the cosmological constant (matter) (radiation) relative to the critical density ρ c .
The expression in (39) is closely related to those for energy injection from DM decay and annihilation.Moreover, constraints exist already from the epochs of BBN, CMB decoupling, and today on annihilating and decaying DM models.We may thus reinterpret these constraints in the context of radial mode radiation from strings.The energy deposited per unit volume per unit Table I.A summary of the observational constraints derived in this work on axion strings from the high-energy SM radiation they emit in the scaling solution.The BBN and CMB constraints arise from primordial energy injection at these epochs, while the gamma-ray constraint is from present-day gamma-ray searches.We provide the approximate redshift of the constraint, the maximum axion mass m max a for the constraint to apply (though considering domain wall formation more non-trivial constraints could apply at higher axion masses), and the upper bound on the decay constant f max a for a given branching ratio B of the radial mode to SM final states.More formally all of the upper bounds apply to fa/ 0.33 c 1 B 25 ξ(z * ) , but for this table we fix c = 0.33, which parameterizes the energy injection into radial modes, and the number of strings-per-Hubble volume ξ(z * ) at the values given in the table.Note that these limits are quoted in terms of radial mode production and decay but also apply, with the appropriate modifications, to the scenario in which the strings directly produce high energy Higgs radiation, as discussed in Sec.V. time from DM decay and annihilation is, respectively, where in the top line Γ DM→SM is the DM decay rate to SM final states, Ω DM is the energy density fraction in DM today, and in the bottom line the DM with mass  m DM annihilates to the SM with velocity-averaged crosssection ⟨σv⟩.

A. Constraints from BBN
During the radiation-dominated epoch, the energy deposition for string-induced radial mode decay and for DM annihilation scale the same with redshift, up to the logarithmic dependence of ξ on z.Thus, we may determine the upper limit on f a by identifying where ξ(z * ) is the value at the epoch given by redshift z * where the DM annihilation constraint is evaluated.Note that above M pl ≈ 2.4 × 10 18 GeV is the reduced Planck mass.
We now consider the constraints on axion strings from BBN by reinterpreting the BBN constraints on DM annihilation.DM annihilation constraints during the epoch of BBN arise from two different mechanisms related to (i) hadronic energy injection, and (ii) photonic and leptonic energy injection.Hadronic energy injection may increase the neutron-to-proton ratio, which in turn increases the primordial 4 He mass fraction [86][87][88][89][90]. Electromagnetic energy injection, on the other hand, may photo-dissociate nuclei [87,88,[90][91][92].For example, photo-dissociation of 4 He may lead to the overproduction of 3 He.The hadronic energy injection constraints and the electromagnetic energy injection constraints scale differently with DM mass: for masses m DM well above a GeV, the electromagnetic constraints are for fixed ⟨σv⟩/m DM , since they are constraints on the total injected energy, while the hadronic constraints scale as ⟨σv⟩/m 3/2 DM since they are proportional to the number of injected nucleons [93].This implies that for very large DM masses the electromagnetic constraints are more powerful; thus, the electromagnetic constraints are the leading ones to use when constraining heavy radial mode decay.
We adopt the S-wave DM annihilation constraint from [88] on the 3 He/D ratio, which states where ϵ vis is the fraction of annihilation energy that goes into photons and e ± .For both W + W − and hh final states, ϵ vis ≈ 1/2, which is the value we adopt. 7ote that photo-dissociation only becomes efficient when T ≲ 0.3 keV [88], which implies that the BBN constraint requires m a ≲ 2 × 10 −23 eV.When the mass is less than this critical value, the decay constant is constrained to be less than: Note that above we use c ≈ 0.33 from Sec. III D. Recall from Sec. II that at this same epoch, constraints exist from not overproducing the observed value of N eff from the axion radiation emitted by the string network, constraining f a ≲ 9 × 10 14 GeV.The constraint in (44) is stronger than the N eff constraint at this epoch for B ≳ 0.01.

B. Constraints from the CMB
At redshifts z ≲ 3000, energy injected into the SM plasma by annihilating or decaying DM, or radial mode decay from strings, changes the ionization history of the ordinary matter.These changes, in turn, change the CMB angular power spectrum, which is accurately measured and modeled under standard cosmology.CMB constraints on annihilating and decaying DM have been extensively studied (see, e.g., [94][95][96][97][98][99][100][101]).Recently, Ref. [101] found that for W + W − final annihilation states8 the annihilation cross-section is constrained by the CMB power spectrum to be smaller than These constraints predominantly arise from energy injection at redshifts z ≈ 600.Note that the CMB constraints, as in the case of BBN, are a function of ⟨σv⟩/m DM since they constrain the total injected energy.
It is less straightforward to reinterpret the CMB angular power spectrum constraints in terms of strings.This is because the injected energy from string emission does not redshift the same way as either energy injection from DM annihilation or decay, as the Universe goes through matter radiation equality into the epoch of matter domination.Below, we carefully compute the upper limit on f a by performing a dedicated CMB power spectrum analysis for the specific form of the redshift-dependent energy injection appropriate for string emission.First, though, let us approximate the upper limit by assuming that the energy deposition happens instantaneously at z inj = 600.Then, we may translate the CMB constraints on DM annihilation9 to constraints on radial mode decay by equating (39) and dE dtdV DM ann from (41) at z inj .This leads to the relation where both f a and ⟨σv⟩ represent the upper limits.This then implies that Performing this calculation more carefully leads to a similar result, as we discuss below.Before describing our simulation framework for energy injection from axion strings, let us briefly comment on a crucial assumption in the above estimate.We assume that, apart from the redshift dependence, the energyinjection signal from axion strings has the same phenomenology as that from DM decay and annihilation.On the other hand, we know that in detail this assumption cannot be true, since the axion strings deposit their energy in narrow cylinders around the string cores, while DM annihilation and decay processes deposit energy in a relatively smooth fashion throughout the entire Hubble volume of interest.On small angular scales, we thus expect to see differences between the morphology of the DM-induced signals and the string-induced signals, for example as manifest through the angular power spectrum of the CMB.
Roughly speaking, the co-moving separation between strings at z inj is (R(z inj )H(z inj ) ξ(z inj )) −1 .From this estimate we may compute the typical angular separation between strings θ by dividing the length scale above by the co-moving distance from today to the redshift z inj ; then, identifying the angular multipole number ℓ through θ ∼ π/ℓ, we estimate that only for ℓ ≳ 800(z inj /600) 1/2 ξ(z inj )/30 should the morphology of the strings signal differ from that of DM annihilation and decay.The analyses we describe below make use of Planck data with ℓ ≲ 2500; thus, the high-ℓ modeling we perform likely underestimates the true anisotropy.We suspect that this implies the limits presented here are conservative, but this should be checked in the future with a dedicated analysis that accounts for the anisotropies on small angular scales arising from strings depositing their energy locally.On the other hand, we verify that limiting the analysis to ℓ ≤ 800, corresponding to z inj = 600, does not qualitatively change the sensitivity of the search, since the string-induced energyinjection signal predominantly appears at low ℓ; in particular, we check that limiting ℓ ≤ 800 actually makes the upper limit stronger by ∼20%, though at a level expected from statistical fluctuations alone.

Dedicated CMB bound on axion strings
To go beyond the approximation in (47) we perform a dedicated analysis that accounts for the unique redshift dependence of the axion-string-induced energy injection.We use CLASS to compute the CMB anisotropy with energy injection from radial mode decay, using the formalism described in [102].The steps of the computation are as follows: (i) the recombination histories of the matter temperature and ionization fraction of H are computed using a modified version of RECFAST [103], which allows for an exotic energy injection; (ii) the evolution of matter and metric perturbations is then determined following cosmological perturbation theory, by solving (in Fourier space) the coupled linearized Einstein plus fluid equations for the photons, baryons, DM, and neutrinos; (iii) finally, a line-of-sight integral is used to compute the angular power spectra for the temperature (TT), E-modes (EE), and their cross-spectrum (TE).
We now discuss the recombination modeling in more detail.In the following we assume the rate of energy injection from radial mode decays is given by (39), where the prefactor c is extracted from our simulations, and we take a constant value of ξ for simplicity.As the stable radial mode decay products (high-energy photons, electrons, and positrons) cool, they deposit their energy in various channels, including into the ionization of H and He, Ly-α excitations, free-streaming continuum photons, and heating of the intergalactic medium.These processes increase the ionization fraction of H. CLASS has built-in options to include exotic energy injection from, e.g., DM annihilation/decay.We perform straightforward modifications to the thermodynamics module to include energy injection from radial mode decay from axion strings.
To obtain the fraction of energy deposited into each channel for s → W W decay, 10 we use DarkHistory [104], which models the evolution of the ionization fraction of H and the gas temperature and the cooling of decay products.Importantly, it includes the back-reaction of changes in the ionization fraction and gas temperature on the various energy-loss mechanisms.As DarkHistory requires an initial electron and photon energy spectrum, but limits the mass of the annihilating particle to 10 5 GeV and the highest energy bin to ∼5.4 TeV for both the electron and photon, we assume a box function for these spectra spanning a single energy bin centered at 1 TeV with width ∼60 GeV for the electron and ∼80 GeV for the photon, with the integrated energy obtained from HDMSpectra [105], which computes decay spectra for GUT-scale masses.We verify that our results are insensitive to the width and mean energy of these box spectra, and thus, to a good approximation, depend only on the total energy injected; for example, injecting GeV energy particles instead of TeV energy particles leads to test statistic changes less than ∼10%.As a consistency check, we also verify that using the deposition fractions from DarkHistory in CLASS retrieves the electron ionization fraction computed by DarkHistory to the 10% level.The deposition fractions calculated by DarkHistory are not expected to be accurate below this level.
We set constraints on the injected energy using Planck 2018 CMB data [106].
We use cobaya [107] to interface between CLASS and the Planck 2018 likelihood code Plik (described in detail in [108]).In particular, we use the likelihood planck 2018 highl plik.[TT|TTTEEE]lite, a version of the Planck 2018 high-ℓ T T + EE + T E binned likelihood which is marginalized (in the Bayesian sense) over 47 nuisance parameters that model the foreground.
As noted in [109], modifications to the T T spectrum from energy injection are almost degenerate with the pri-mordial scalar spectral index n s and amplitude A s .This degeneracy is broken by the polarization information in the EE and T E spectra.To account for the degeneracy we profile (in the frequentist sense) over the ΛCDM parameters (h, Ω b , Ω cdm , A s , n s ), fixing all other cosmological parameters to their best-fit values from the Planck 2018 T T + EE + T E + lowℓ + lowE + lensing data analysis. 11When profiling, we also fix the deposition fractions to their values calculated with the Planck 2018 best-fit parameters for ease of computation.Note that for simplicity we do not profile over the optical depth τ and instead fix the reionization history to the default model of DarkHistory.As reionization only affects redshifts z < 20, our constraint would only differ negligibly under small perturbations to the reionization history.
More precisely, we construct the profile likelihood ratio where p is the Planck partially-marginalized likelihood, given the data d, the 5 ΛCDM nuisance parameters θ nuis , and the signal parameter f a .The quantities { θnuis , fa } represent those which maximize the marginalized likelihood, while θ denotes the nuisance parameters that maximize the likelihood at fixed f a .We then compute the test statistic which is illustrated in Fig. 12.We invoke Wilks' theorem and set the one-sided 95% upper limit as the value of f a > fa for which t(f a ) ≈ 2.71 (see [1] for details).Note that this analysis is formally a hybrid Bayesian-frequentist analysis, since the likelihood p has been marginalized over the 47 foreground nuisance parameters.Additionally, note we find no evidence in favor of the axion model, given that the test statistic difference between the best-fit point fa and the null hypothesis (f a = 0) is much less than unity.We find the 95% one-sided upper limit to be Comparing to (47) we see that this upper limit is similar to the naive estimate based off of translating the DM annihilation limit.
There are two caveats related to the bound in (50) that are important to consider related to ξ.First, we note that the energy injection does not take place instantaneously at z inj ∼ 600 but rather over a range of redshift values near this characteristic redshift.Thus, taking a constant number of strings-per-Hubble patch is not completely correct, but in practice since ξ varies logarithmically with time we verify that this approximation is valid to the precision quoted.Second, and more importantly, the time specified by z inj is within the epoch of matter domination, and the string network has a difference scaling solution during matter domination than during radiation domination, as we discuss in App.D. In particular, if we assume the string network follows the radiation-epoch scaling solution until z in then we expect a characteristic value ξ(z inj ) ∼ 30, while using the matter-dominated scaling solution in App.D we would infer ξ(z inj ) ∼ 7.In practice, we expect ξ(z inj ) to be between these two values, since z = 600 is only slightly below matter-radiation equality.Computing ξ(z inj ) directly through simulations is difficult because of the large log(m s /H) values where matter-radiation equality occurs.Thus, we simply note that the pre-factor in (50), accounting for the ξ(z inj ) dependence as well, may be as large as ∼ 4.5 × 10 12 GeV if we use the lower bound ξ(z inj ) ≳ 7.
To help unpack this analysis in Fig. 13 we illustrate the angular power spectra for an energy injection signal from radial mode decay corresponding to f a c 0.33 B 1 ξ 30 = 8.15 × 10 12 GeV.We observe a suppression of the power spectrum which is stronger at smaller scales.

C. Constraints from present-day gamma rays
Additional constraints appear if one assumes that the network persists until z = 0 (m a ≲ 10 −33 eV).In this case, we may reinterpret the results of searches for extragalactic DM decay.High-energy gamma-rays and e ± from DM decay are reprocessed to lower energy gammarays through a cascade of electron-positron pair production and inverse Compton scattering off of background radiation fields, and for high m s , well above the PeV scale, the spectrum of gamma-rays observed on Earth approaches a universal spectrum that peaks, in terms of the flux E 2 γ dΦγ dEγ (units of GeV/cm 2 /s/sr), between 10 and 100 GeV [110]. 12Ref.[111] constrained τ ≳ 1.3×10 27 s for DM decay to b b for m DM > 10 11 GeV; considering that the upper limit only depends on the energy injected into non-neutrino species, we may infer that the upper limit for hh or W + W − final states would be the same to within 10%.Ref. [112] used more aggressive modeling of the extragalactic gamma-ray background to constrain τ ≳ 10 28 s for m DM ≳ 10 9 GeV, for both W + W − and hh final states, with the limits insensitive at the less than 10% level to m DM for m DM ≳ 10 9 GeV [111].Since the energy injection is dominated by decays with z ≪ 1, we may translate these limits to limits on f a through the relation where ξ(z = 0) is the value at z = 0. Using the lifetime bound from [112] this then implies This upper limit is marginally stronger than the CMB upper limit in (50), though -as we discuss further in App.D, the number of strings-per-Hubble is likely more comparable to ξ(z = 0) ∼ 7 at this epoch.

VII. DISCUSSION
In this work, we set strong constraints on axion-like particle strings that survive to temperatures at or below that of BBN.We study the effects of primordial energy injection from the decays of massive radial modes, released during the evolution of the axion string network.The strength of the derived upper limits depends on how long the network persists and on the branching ratio of the radial modes to SM final states.For relatively generic branching ratios B ∼ 0.1 to SM final states, the upper limit from BBN (CMB) is around f a ≲ 3 × 10 14 GeV (f a ≲ 5 × 10 12 GeV).These upper limits rely crucially on understanding how the axion-string network sheds energy into radial modes and directly into SM final states, which we study through a combination of analytic arguments and dedicated numerical simulations.
QCD axion string simulations predict that the QCD axion decay constant should be within a factor of a few of 10 11 GeV in order to produce the correct DM abundance if the PQ symmetry is broken after inflation [23,25], assuming a standard cosmological history.It is thus also well motivated to consider axion-like particle strings with similar decay constants, as these axion-like particles may accompany the QCD axion in some realization of the axiverse paradigm.Interestingly, this region of parameter space should be probed by the next generation of CMB experiments [113], making CMB probes an exciting future possible discovery channel for axion-like particles.
Lastly, we note that while this work focuses on global strings it is possible that some of the results may be relevant to local strings, such as strings in the Abelian-Higgs model.The Abelian-Higgs model is obtained simply by taking the same Lagrangian used in this work to produce axion strings and gauging the U (1) PQ symmetry with an abelian gauge field (see, e.g., [114] for a review).After spontaneous symmetry breaking both the radial mode and the gauge field are heavy, with masses of order f a .The previously massless axion is now "eaten" by the massive gauge field, such that there are no light degrees of freedom.
It is typically assumed that local string networks evolve by emitting gravitational-wave radiation, with the direct production of massive modes exponentially suppressed and thus not relevant for dynamics [115].On the other hand, the question of whether heavy-state emission is truly exponentially suppressed is unresolved, with some works claiming that it is relevant for the network dynamics [39,40,47,48].Our work concludes that for global strings the production of heavy states is only logarithmically suppressed relative to the production of massless states.This may suggest that in the local string scenario the strings are able to produce heavy states with a similar efficiency as in the global case.Heavy-mode production from local strings would have a number of important implications, including making such strings susceptible to energy injection constraints along the lines of those discussed in this work.Dedicated local-string cosmological simulations are needed, however, to understand to what extent the results found here for global strings carry over to local strings.(4f 2 a s 2 + 4f a s 3 + s 4 ) Thus the mass of the radial mode is given by m s = √ 2λ Φ f a .From the above we can compute the relevant decay widths, where c w ≡ cos θ w , s w ≡ sin θ w , θ w is the Weinberg angle, and α w ≡ g 2 2 /(4π), with g 2 the coupling constant of SU (2) L .
Appendix B: Pre-evolution and adiabatic regime We follow the procedure described in [20] to avoid transient radial-mode excitations at small values of log(m s /H).Instead of starting the simulation from a thermal initial state and having strings form dynamically by explicitly simulating the PQ phase transition we start the simulation in this procedure after the PQ phase transition from a pre-evolved initial state that already contains strings.This pre-evolved initial state is generated by evolving a thermal initial state within a modified physics scenario where strings have a constant width.This width can be tuned to match the string width at the intended starting time of the actual simulation, η i ≈ 2.3 (log(m s /H) = 2).Additionally, this scenario contains a moderate amount of Hubble friction, which allows us to evolve strings for far longer such that they have sufficient time to de-excite.This is achieved by changing the relation R/R 1 = t/t 1 to R/R 1 = t/t 1 and forcing m s ∝ 1/R.These modified equations of motion read where the dimensionless fields ψ i = ψ 1 , ψ 2 are given by Φ = (ψ 1 + iψ 2 )f a / √ 2. The initial state for the pre-evolution stage is a thermal state with wavenumbers up to a certain threshold in each spatial direction, see [22] for more details.We perform a total of nine different simulations that differ in the wavenumber threshold: Two initial states are statistically independent realizations using the first 13 wavenumbers, the other initial states are based on the first 10, 15, 18, 20, 25, 30, and 35 wavenumbers, respectively.Due to the existence of an attractor solution [20] the impact of our choice of threshold is marginal at sufficiently large log(m s /H).The pre-evolution simulation is performed with 512 3 grid cells as the resolution is not an issue here due to the constant string width.The simulation starts at η = 1 and ends when the total string length is close to the attractor solution at η i , ξ(η i ) ≈ 0.18.
Furthermore, as the underlying physics changes instantaneously when starting the main simulation from a preevolved state, we introduce a short adiabatic period between the two regimes.We do this by computing ψ ′′ 1,2 in both scenarios and combining them, ψ

Appendix C: Axion-Higgs simulations
The equations of motion for the coupled PQ-Higgs system are derived analogously to the PQ-only case and read along with ) Here, given the ∼ 12 order of magnitude hierarchy between the SM Higgs VEV and f a , we set the Higgs VEV far away from strings to zero.
While we chose a simulation volume that is slightly larger than that in our PQ-only simulation it is nevertheless not large enough to avoid the necessity of preevolving the initial state to mitigate the effects of transient oscillations of the string cores.The corresponding pre-evolution equations of motion are derived analogously to the PQ-only simulation described in App.B and read ) Our pre-evolution procedure is identical to that described in App.B but with an increased static grid size of 1024 3 cells to accommodate the larger volume.The thermal initial state includes the first 18 wavenumbers in each spatial direction where the Higgs field is generated analogously to the PQ field but without an effective mass (see [22] for more details).
We use the final state of the pre-evolution as our initial state for the main simulation starting at log(m s /H) = 2.This state is taken as our coarse level but an extra refinement level with ∆x → ∆x/2 is introduced whenever the string core width would be resolved by less than four grid sites on the finest level.This means by the end of the simulation we will have three refinement levels on top of the coarse level.A static lattice simulation would have needed 8192 3 cells to match this dynamic range.In which part of our simulation volume the refinement level is placed is based on two criteria: (i) the location of string cores, and (ii) a data-driven convergence criterion.The refined region is re-adjusted frequently every ∆η = 0.12.
String cores are identified using the procedure outlined in [116].We ensure the refinement region around string cores is large enough that even a string segment moving with the speed of light will always be at least an entire string width away from any coarse-fine boundary until the grid is readjusted.This criterion ensures that the strings themselves are properly resolved at all times, however, the emission leaving the string may not be.To guarantee this emission is resolved appropriately as well we additionally employ a data-driven method that estimates the convergence at individual grid cells.
The basic idea behind the data-driven technique is to independently evolve a cell at two different resolutions for a short period of time.The difference between the results informs us about the size of the numerical truncation error due to the finite resolution.When this difference gets too large it means that numerical convergence is locally bad at the current resolution and refinement is needed.In practice, this is done easily within an AMR framework as we are already evolving the field at different resolutions.That is, at the end of every time step and before level synchronization we can compare the results of, let us say, the coarse level and the first refinement level.If this difference ∆X exceeds a threshold τ the area around this cell will be covered by the second refinement level during the next regrid.In order to identify problematic cells on the coarse level, however, an even coarser level with half the resolution of the coarse level is required.This setup is known as a self-shadow hierarchy.The parameter τ is chosen empirically and we find τ = 10 −3 to work well for X = ψ i , h i , and τ = 10 −3 /∆η ℓ for X = ψ ′ i , h ′ i .Our integration scheme is the same as that of our PQonly simulation.The simulations are performed on the NERSC Perlmutter GPU cluster and utilize 256 NVIDIA A100 GPUs and 64 AMD EPYC 7763 CPUs for about an hour per run.

Appendix D: String density in matter-domination
To study the evolution of the string density at times after matter-radiation equality, we simulate the string network in a matter-dominated cosmology.The simulations are such that the radial mode acquires its broken VEV and strings form when the Universe is already matterdominated.Of course, this does not correspond to the physical scenario of string network formation at the PQ phase transition in a radiation-dominated epoch, and the subsequent evolution of the Universe through the epoch of matter-radiation equality.However, the choice of initial conditions in the simulation is not important as here we are only interested in the scaling regime during the matter-dominated era.
Neglecting radiation energy density, we have at late times Ω m + Ω Λ ≈ 1, and the energy densities of matter and cosmological constant are equal when the scale factor is η eq = (Ω m /Ω Λ ) 1 3 .We perform two AMR simulations: (i) a matter-only simulation where we impose η eq → ∞ and which ends at log(m s /H) = 5.78, and (ii) a simulation entering the cosmological constant-dominated epoch with η eq = 20, which ends at log(m s /H) = 4.71.We use the same AMR simulation setup as for our axion-Higgs simulation in App. C.However, as we are not interested in measuring any emission spectrum we can safely skip the pre-evolution procedure.Instead, we simulate explicitly through the PQ phase transition by starting at η = 0.1 from a thermal initial state with the first 9 wavenumbers included.
Both simulations evolve the following equations of mo-tion [η 4 + η 3 eq.η]ψ ′′ i + [4η 3 + 2.5η 3 eq.]ψ ′ i − [1 + η 3 eq.] ∇2 ψ i + λψ i η 2 |ψ| 2 − 1 + where the thermal term ensures the PQ symmetry breaks early on in the simulation.These are the Euler-Lagrange equations of the Lagrangian (5) in an FRW metric with Hubble parameter given by (40) and Ω rad = 0. We also define the dimensionless fields ψ i as in App.B. Evolution of the string length per horizon ξ in the matter-dominated (black solid) and radiationdominated (black dotted) scenarios, and during a transition between matter-dominated and cosmological constantdominated epochs (dashed).For the matter-dominated case, the model ξ = c0 + c1 log is fit to the ξ data., (D3) using our convention for the Hubble scale for the simulations described in this Appendix, where by definition H = f a at η = 1.The string length per horizon is compared between the various cosmological scenarios in Fig. 14.For the matter-dominated simulation, we follow the method of Ref. [25] to fit ξ = c 0 +c 1 log(m s /H) in the interval log(m s /H) = [4, 5.48] to find c 1 = 0.0584 ± 0.0013.In the scenario where the Universe transitions from matter to cosmological constant domination, ξ begins to drop exponentially, as expected, when the energy densities in matter and cosmological constant are equal.The string length for the radiationdominated scenario in Fig. 14 is extracted from the AMR simulation discussed in [25], which was performed with 2048 3 grid-cells.
In Fig. 15 we vary the number K of wavenumbers included in the initial thermal state for the matterdominated simulation.The proximity of ξ at large values of log(m s /H) as we vary K = 9, 17, 30 indicates that, like in the radiation-dominated era, the string network approaches an attractor solution in the matter-dominated era as well.The attractor solution implies that at late times the number of strings per Hubble patch is the same regardless of the initial condition.In our case, we find c 1 = 0.0584±0.0013,0.0476±0.0003,and 0.0548±0.0003,for the simulations with initial conditions K = 9, 17, 30, respectively.Note that the error bars quoted above and shown in Fig. 15 are statistical.We may infer a systematic uncertainty from the variance in c 1 between the simulations with varying initial mode numbers, leading to the estimate c 1 = 0.0584 ± 0.0013 stat ± 0.0055 sys .Note, also, that in reality the network evolves smoothly from the radiation dominated scaling solution to the matter dominated scaling solution around matter-radiation equality; depending on the cosmological epoch of interest, the network may not be well approximated by the scaling solution in either epoch but could take on intermediate values.

Appendix E: Single loop spectrum
We simulate a single circular collapsing string with instantaneous radius R to study the spectral shape of the resulting instantaneous axion emission spectrum; we find evidence that the instantaneous emission spectrum scales with axion momentum k as 1/k between the characteristic frequency k ∼ 2π/R and k ∼ m s , supporting the consistent conclusions reached in [38].
Our simulation setup for this test is mostly the same as that of our axion-Higgs simulations but with an artificial initial state that creates a singular and perfectly circular string loop.To achieve this initial condition we generate the initial fields at rest with Φ 1 (i, j, k) = 1 − 2/(1 + e −0.05[|x−N/2|−N/4] ) Φ 2 (i, j, k) = sin(2πk/N ) , (E1) in index space {ijk} over the three Cartesian dimensions, with N cells in each direction (e.g., i = 0, 1, • • • , N − 1) and x = (i, j, k) T .The result is a circular string loop with a radius of N/4 and approximately periodic boundary conditions.This construction has the advantage that the initial radius R 0 can be controlled by choosing a box length L. The exact details of the initial state such as the overall field amplitude are not crucial as we start our simulation in the unbroken phase at η = 0.1.By explicitly simulating the PQ phase transition the string is generated dynamically.Since we study a single string we can use a relatively aggressive AMR setup.Our coarse level consists of merely 1024 3 grid sites, but this is compensated by five refinement levels at the end of our simulation at η ∼ 50 (log(m s /H) ∼ 8.2) to maintain at least four grid sites per string core width.On a static lattice such a simulation would require a grid of 32, 768 3 grid sites.We chose L = 100/(R 1 H 1 ) such that the initial string radius is R 0 = 25/(R 1 H 1 ).The evolution of R is shown in Fig. 16 with corresponding string length ξ in Fig. 17.We perform the simulation on the NERSC Perlmutter GPU cluster and utilize 256 NVIDIA A100 GPUs and 64 AMD EPYC 7763 CPUs.An illustration of the axion energy density is shown in Fig. 18.
We compute the instantaneous axion emission spectrum F ∝ (1/R 3 ) d dt (R 3 ∂ρ a /∂k) with ∂ρ a /∂k the timedependent differential axion energy density spectrum.Figure 19.Instantaneous axion emissions spectrum F for a collapsing circular string loop.We perform a power law fit to the regime between k ∼ 0.1ms and k ∼ ms (dashed line).Numerically, we obtain the time-derivative via finite differences uniform in log(m s /H) with ∆ log(m s /H) ∼ 0.25 (see, e.g., [20,25]).We fit a power-law model F ∼ 1/k q to the instantaneous spectrum between k ∼ 4π/R and k ∼ m s .We do not, however, extend the fit below k ∼ 0.1m s in case of large R.An example spectrum and the corresponding fit are presented in Fig. 19.The power-law index q for each fit is shown in Fig. 20.We perform a linear fit to these indices, analogously to the fit of Γ s /(8H 3 ξπf 2 a ), and find q = 1.06 ± 0.06.This result supports our claim of a 1/k scaling in this regime to within ∼6% accuracy.
In Sec.III B we argue that sub-horizon-size string loops are distributed as dn ℓ /dℓ ∼ 1 ℓ .On the other hand, above we find evidence that string loops of radius R emit axions with instantaneous spectra F ∝ 1/k between, roughly, 2π/R and m s .Let us now combine these two points to argue that the network as a whole should emit axions with instantaneous spectrum F ∝ 1/k.
The F ∝ 1/k scaling of the full network can be understood analytically in the following way.We can write the axion emission spectra by including the contributions from all loops of size ℓ as, The factor dn ℓ /dℓ governs the number density of loops as a function of their size ℓ, and as argued and seen in simulations in [25], we expect dn ℓ /dℓ ∝ 1/ℓ.The spectral function F ℓ (k) governs the instantaneous emission spectrum from a single loop with size ℓ and is normalized, without loss of generality, via dkkF ℓ (k) = 1.As an illustration, we first focus on the case where each loop of size ℓ emits axions at a single frequency k ∼ 1/ℓ, i.e., F ℓ (k) = δ(k − 1/ℓ).We then arrive at, ) log(m s /H), for some constant c, where Γ k≳ms a denotes the axion emission with k ≳ m s and Γ k≲ms a

Figure 2 .
Figure 2. 2D Fourier transform of the string energy density for an infinite, straight string, from (11).

Figure 3 .
Figure 3. (Top panels) 2D projection of the radial mode energy ṡ2 at the end of our 3D simulation investigating radial mode emission around log(ms/H) ∼ 6.5.The full simulation box, spanning ∼1.5 Hubble lengths, is shown on the right with a detailed view shown on the left.Axion strings stand out as bright closed loops with strong emissions in particular around kinks and recent string re-connections.(Bottom panels) The same state of the string network but illustrated for the axion energy density ȧ2 instead of that of the radial mode.The axion emission has more support at long wavelengths relative to that of the radial mode.

Figure 4 .
Figure 4. Energy densities ρ for the axion (black) and radial mode (grey) as extracted from the simulation.

Figure 5
Figure 5.The quantity z computed for the radial mode, defined in (14) and(15), describes how the instantaneous emission red-shifts at production.Completely non-relativistic (relativistic) radiation has z = 3 (z = 4), with the free particles scaling with the scale factor like ρs ∝ R −z .The emitted spectrum of radial modes is semi-relativistic, giving the intermediate z shown.
enforcing Y = −1 or Y = 0 for the hypercharge of Q R .If Y = 0, though, the axion does not have non-trivial interactions with gauge fields at late times.If instead the quantum numbers are (1, 2, Y ) then the Yukawa couplings must be to right-handed SM leptons via

Figure 8 .
Figure 8.The difference in tension between the infinitely straight Higgs sheath solution and the solution with no nontrivial Higgs field profile (h = 0).We illustrate this difference as a function of λH with the other dimensionless coupling constants fixed to the indicated values.Negative values indicate that the Higgs sheath profile is stable.

Figure 9 .
Figure 9.As in Fig. 3 but for the simulations including Higgs fields.(Top panels) 2D projection of the Higgs energy ḣ2 towards the end of our fiducial 3D simulation around log(ms/H) ∼ 7.0.The full simulation box, spanning ∼1.65 Hubble lengths, is shown on the right with a detailed view shown on the left.(Bottom panels) The same state of the string network but illustrated for the axion energy density ȧ2 instead of that of the Higgs.Animations available here.

Figure 13 . 33 B 1 ξ
Figure 13.We illustrate the binned CMB anisotropy data (T T , T E, and EE) from the Planck 2018 Data Release on top of the best-fit model for the CMB anisotropy without any energy injection computed from CLASS (fa = 0).We also show the best-fit model, profiled over nuisance parameters, with fa c 0.33 B 1 ξ 30 = 8.15 × 10 12 GeV.
Figure 14.Evolution of the string length per horizon ξ in the matter-dominated (black solid) and radiationdominated (black dotted) scenarios, and during a transition between matter-dominated and cosmological constantdominated epochs (dashed).For the matter-dominated case, the model ξ = c0 + c1 log is fit to the ξ data.

30 Figure 15 .
Figure 15.String length per horizon in the matter-dominated scenario, varying the number K of wavenumbers included in the initial thermal state.

Figure
Figure Evolution of the string loop radius R as a function of log(ms/H) as measured in our simulation.

Figure 17 .
Figure 17.String length ξ as a function of log(ms/H) for a circular decaying string as measured in our simulation.

Figure 18 . 6 F
Figure 18.2D projection of the axion energy density of a circular decaying string at log(ms/H) ∼ 7.6.The radius of the string loop at this snapshot is is R ∼ 12.8/(R1H1), i.e. about 0.34 Hubble lengths.

Figure 20 .
Figure 20.Power law index q from fits to the instantaneous axion emission spectra F from a collapsing circular axion at different log(ms/H).We perform a linear fit to those indices (dotted line) with 1σ uncertainty indicated by the grey band.The fit yields q = 1.06 ± 0.06.