Recompute the radiation-era plateau of a cosmic-string background

Read From a loop population to the background spectrum first; it derives every equation used here. The paper is Auclair et al., 2020, arXiv version 1, published in JCAP 04 (2020) 034.

Get the materials and run them

Download the scripts, inputs and saved outputs and extract the archive. It contains six Python files, the three saved JSON files and a short README. The scripts use only the Python standard library. The saved run used Python 3.9.6 and took 37.4 seconds.

cd cosmic-string-background/scripts
python3 reproduce_benchmarks.py ../my-run

The run writes inputs.json, outputs.json and failure-cases.json into the folder you name. Compare them with the saved copies in the archive's benchmarks folder; the same files are also served as inputs, outputs and failure cases. Run the command from the scripts folder: the main script imports its five siblings from there.

Fix the inputs and units

Every time and length is in units of $1/H_0$ with $c=1$, frequencies are in hertz, and densities are reported as $h^2\Omega_{\mathrm{gw}}$. The fixed inputs are:

The model blends the radiation and matter loop amplitudes and birth sizes linearly in $\rho_{\mathrm r}/(\rho_{\mathrm r}+\rho_{\mathrm m})$. That blend is a teaching choice. The paper solves the one-scale equations through the transition instead, so any result that depends on the transition is a property of this model only. NANOGrav Collaboration, 2023, Section 5.4 describes the same Model I family with $\alpha_{\mathrm L}\simeq0.37$ and $F=0.1$.

Two independent routes to the same spectrum

The first route is the redshift integral of eqs. 4.4 and 4.5. It evaluates $C_1$ by Simpson's rule in $u=\ln(1+z)$ with 800 intervals, from today back to the birth boundary of the loops that can emit at the requested frequency. It uses the loop density $\mathsf n(\ell,t)$ and its Jacobian.

The second route never builds a loop density. It follows cohorts of loops from birth, lets each loop shrink at $\Gamma G\mu$, deposits the energy it radiates in its fundamental mode into logarithmic frequency bins after redshifting, and divides by today's critical density. An error in the Jacobian or in the $(1+z)$ powers of the first route would not appear in the second.

from loop_population import LoopPopulation
from background_by_redshift_integral import fundamental_only_density
import string_cosmology as cosmology

table = cosmology.CosmicTimeTable()
loops = LoopPopulation(table, 1e-10, 0.1)
print(fundamental_only_density(loops, 1e4))

This prints about $4.537\times10^{-10}$, the fundamental-mode plateau for $G\mu=10^{-10}$ and $\alpha=0.1$.

What the recomputation gives

The pure-radiation integral of the simulation-based loop density, eqs. 3.18 and 3.19, gives a plateau coefficient of 8.0425. The paper prints 8.04 in eq. 5.2. An analytic evaluation of the same integral gives 8.0425 as well.

Eq. 6.1 matches the numerical fundamental-mode plateau at $10^4\,\mathrm{Hz}$ within $0.06\,\%$ for 24 of 25 cases spanning $G\mu=10^{-8}$ to $10^{-17}$ and loop sizes from $0.1$ down to $10^{-3}\,\Gamma G\mu$. The one exception, $G\mu=10^{-17}$ with $\alpha=5\times10^{-19}$, is $4.6\,\%$ low. There the loops heard at $10^4\,\mathrm{Hz}$ can only have been emitting at $1+z$ below about $6\times10^4$, close to matter–radiation equality at $z\approx3366$, so that frequency is not yet deep in the radiation plateau.

The cohort histogram and the redshift integral agree within $0.9\,\%$ at 27 frequencies from $1.1\times10^{-12}$ to $4\times10^{3}\,\mathrm{Hz}$, for $G\mu=10^{-10}$ and $10^{-17}$. Summing all the energy radiated by every cohort, without any spectrum, agrees with the frequency integral of the spectrum to $0.3\,\%$ for both tensions.

As the tension falls toward zero the plateau divided by $\sqrt{G\mu}$ stays at $4.537\times10^{-5}$ from $G\mu=10^{-12}$ to $10^{-16}$. At $10^{-18}$ and $10^{-20}$ the ratio rises to $4.549\times10^{-5}$ and $5.16\times10^{-5}$, because the spectrum reaches its plateau at higher frequency when the tension is lower and $10^4\,\mathrm{Hz}$ is no longer on it. At one cycle per year the background falls as $(G\mu)^2$ at small tension, reaching $8.0\times10^{-27}$ at $G\mu=10^{-20}$. The background goes to zero with the tension, as it must.

Two printed values the model does not reproduce

Eq. 6.1 and eq. 6.2 print the prefactor $(128\pi/9)A_{\mathrm r}\Omega_{\mathrm r}h^2$ as about $1.02\times10^{-2}$. With the paper's own inputs it is $1.0146\times10^{-3}$, a factor of $10.05$ smaller. Eq. 6.5 shows the same factor of ten. The symbolic form agrees with eq. 5.2 and with both numerical routes, so the scripts use it. They do not edit the paper; the discrepancy is recorded as unresolved in outputs.json under printedPlateauCoefficient.

Section 6.1, in the text after figure 8, quotes a plateau of $h^2\Omega\approx6.04\times10^{-14}$ for $G\mu=10^{-17}$. Eq. 5.2 with constant particle species gives $1.51\times10^{-13}$, so the quoted value is 0.40 of it. The paper's figure 4 shows that the Standard Model's changing species count alters the high-frequency part of the spectrum, which may account for some of the difference. The statement is recorded as unresolved, and left uncorrected, until that effect is computed.

Four ways to get it wrong

The saved failure-cases.json records four mistakes and their size.

1. Dropping the Jacobian, using $1/\alpha$ for $\mathrm dt_{\mathrm i}/\mathrm d\ell$ when $\alpha=\Gamma G\mu$, doubles the plateau: the ratio is 2.000.
2. Using the plateau formula at one cycle per year for $G\mu=10^{-10}$ underestimates the fundamental-mode background by a factor of 2.43.
3. Assuming $\sum P_n=\Gamma$ while summing only $10^5$ cusp harmonics misses $1.79\,\%$ of the power.
4. Starting the network at $z=10^{10}$ instead of $10^{26}$ removes the high-frequency plateau: at $1\,\mathrm{Hz}$ the result falls from $4.54\times10^{-10}$ to $2.6\times10^{-14}$.

What this does and does not establish

The recomputation shows that one model agrees with itself: its loop density, its spectrum and its energy budget agree with each other and with two published closed forms. It does not test whether cosmic strings exist, whether Model I describes them, or whether constant particle species is a good approximation. It fits no data. Comparing it with pulsar-timing results needs the full likelihood analysis described in the guide to reading pulsar-timing and LISA results.