Will LISA Detect Harmonic Gravitational Waves from Galactic Cosmic String Loops?

For the toy models, we adopt the $\mathrm{\Lambda CDM}$ cosmology with Planck 2015 cosmological parametersAde et al. (2016):

Neutrinos are regarded as massive, and $\Omega_{\mathrm{\nu,0}}$ is estimated by refs. Robitaille, Thomas P. et al. (2013); Price-Whelan et al. (2018) using the method described in ref. Komatsu et al. (2011). Due to considerations of computational performance, they are treated as relativistic particles contributing to the radiation density, resulting in a further overall enhancement of the loop number density by $\sim 10$ compared to treating them as massless. This has no negative impact on the toy models. Physical results computed with massless neutrinos are presented in section IV.2.

The distribution of CDM in the Galaxy is well modeled by the Navarro-Frenk-White (NFW) profileMcMillan (2011); Klypin et al. (2002). We therefore adopt the NFW profileNavarro et al. (1996) for the Galactic dark matter halo, which also models the number density of loops clumped in the Galaxy. Consistent with ref. DePies and Hogan (2009), we adopt $R_{s}=21.5\ \mathrm{kpc}$ and $C=10$, with the common definitionDehnen et al. (2006); D’Onghia et al. (2010); Shull (2014) of $R_{\mathrm{vir}}\sim R_{200}$. The turnaround radius $R_{\mathrm{ta}}=1.1\ \mathrm{Mpc}$ is adopted from ref. Chernoff (2009).

With plane wave solutions of linearized gravity in the short wavelength approximation, the total signal as a result of superposition of GW from all loops in the Galaxy in the frequency bin $f$ is simply

by the harmonic addition theoremOo and Gan (2012). Define $A\equiv\sum_{i=1}^{N}h_{i}\sin\delta_{i}$ and $B\equiv\sum_{i=1}^{N}h_{i}\cos\delta_{i}$, where $h_{i}$ given by eq. 23 and $\delta_{i}\in\left(-\pi,\pi\right]$ are individual strain amplitudes and uncorrelated phases, respectively, and $N$ is the total number of loops in the bin. The total strain amplitude is

with the total phase

The simulation computes the superposition for each frequency bin, as the bins are independent of each other.

When loop decay is not significant, $N(f)\propto f^{\frac{3}{2}}$ from cosmology. From fig. 3, this is generally the case at lower frequency, where loops are bigger and are created later, and therefore have not yet experienced significant decay. Loops with lower $G\mu$ experience slower decay and therefore enter the decay regime at higher frequencies. For the LISA frequency range, the decay regime is entirely absent for $G\mu\lesssim 10^{-20}$. The flat region represents the decay regime. For $G\mu\gtrsim 10^{-14}$, loops emitting GW in the entire frequency range have decayed. Loops in the decay regime have similar $t_{\mathrm{c}}$, and therefore similar number densities and sizes initially. However, fractional size differences are amplified by loop decay over time, resulting in these loops having very different sizes today, while their number densities remain similar. In the decay regime, loops with higher $G\mu$ are created much later in time, and hence the lower number densities.

The important point from fig. 3 is that loops can be extremely abundant in the Galaxy, especially at lower $G\mu$. This is clearly not only true in the toy models, but is also the case in the physical model, as rebinning decreases $N(f)$ by a factor of $\frac{1}{T_{\mathrm{obs}}}/0.1f\lesssim 3\times 10^{-7}$. Tracking all loops in the entire Galaxy is certainly far beyond current computational capabilities.

The solution is to only track loops closest to the solar system. For our purpose, it is then possible to recover the statistics for the entire Galaxy based on the outcome of the limited-scale simulation. We call this the “closest loops method”. We verify this method both theoretically (below) and through simulations (section III.2.1).

Our simulations have a signal side focusing on accounting for signal from individual nearby loops, and a noise side focusing on generating GW from unresolved background loops in the Galaxy, i.e. the Galactic confusion. The latter at a given frequency is the result of superposition of GW emitted by all loops in the bin. This can vary significantly across frequency bins due to the uncorrelated phases of GW from background loops, even though the theoretical Galactic background is a slowly varying function of frequency as shown in fig. 1. This variation in the Galactic confusion can mimic the signal to a certain extent. Thus, there are two components to the Galactic confusion for which our simulations need to account for, the Galactic confusion background and the Galactic confusion noise. The former is the mean of the Galactic confusion, $\bar{h}$, which is essentially the theoretically calculated Galactic background. The latter, $h_{\mathrm{n,c}}$, is the amount of variation in the Galactic background across frequency bins.

On the signal side, we will see in section III.3.1 that the likelihood of detection at a given frequency, even in the toy models, is so low that most detections are attributed to just the single closest loop. Thus, simulating just a few closest loops at each frequency is sufficient to account for the signal.

On the noise side, we have both $\bar{h}$ and $h_{\mathrm{n,c}}$. Since overall the power of GW emission from background loops adds. Then from eq. 23,

As $h_{\mathrm{n,c}}$ is a result of the superposition of GW from background loops with uncorrelated phases,

because the interference is effectively a 1D random walk on the strain amplitude.

Both $\bar{h}$ and $h_{\mathrm{n,c}}$ are modulated by the geometry of the spatial distribution of background loops in the Galaxy. Geometry then cancels when comparing eqs. 35 and 36, and we simply have

As $\bar{h}$ is known through theoretical calculations, eq. 37 enables us to recover $h_{\mathrm{n,c}}$ for the entire Galaxy by tracking only the closest subset of background loops. The closest loops method proves to be especially powerful for the interesting region of the parameter space where loop abundance is very high.

We vary both $G\mu$ and $\alpha$ with and without the rocket effect to cover the parameter space. Simulations are performed with $\alpha=0.1$ and $10^{-5}$, with $G\mu$ varying from $10^{-12}$ down to $10^{-20}$, producing 24 parameter combinations.

Simulations cover the peak LISA sensitivity frequency range $f\in\left[10^{-5},1\right]\ \mathrm{Hz}$, with the frequency resolution $\Delta f=\frac{1}{T_{\mathrm{obs}}}$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{month}$ to conserve computational resources. For the signal side, we run 100 realizations for each simulation set including the closest 10 loops at each frequency to generate detection candidates. For the noise side, we run 10 realizations for each set including the closest 1000 loops at each frequency to estimate the Galactic confusion utilizing the closest loops method.

Additionally, we run 100 signal realizations with the rocket effect at $G\mu=10^{-18}$ and $\alpha=0.1$, with $\Delta f$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{year}$. This simulation set effectively constitutes a different toy model from the 24 sets above. It helps clarify physical implications of simulating the toy models, and also serves as a consistency test when we renormalize results from the toy models in section III.3.8.

We consider three main sources of noise:

1. 1.

   The Galactic confusion;

2. 2.

   The stochastic background;

3. 3.

   The LISA instrumental noise with the Galactic WD binaries background.

We set $T_{\mathrm{obs}}=1\ \mathrm{year}$ for the noise level regardless of the $\Delta f$ that a given set employs, as our goal is to make predictions for the possibility of detection over a 1-year mission period.

The full effects of signal modulations due to the LISA orbital motion averaged over a complete orbit are incorporated, by generating the location of each simulated loop in the Galaxy $(r,\theta,\phi)$ according to the distribution of loop number density in the Galaxy, with random $\iota\in[0,\pi]$ and $\psi\in[0,\pi]$.

The cliff-like feature at $f\sim 0.1\ \mathrm{Hz}$ in fig. 4 is located at $f_{\odot}$, and is the same break barely visible in fig. 1 due to the rocket effect. With the absence of the majority of loops in the Galaxy, the feature becomes very pronounced here because the entire population of loops is effectively shifted further away for $f>f_{\odot}$, where distances to these loops become bound by $r_{\mathrm{tr}}$.

The tall vertical lines in $0.01\ \mathrm{Hz}\lesssim f\lesssim 0.1\ \mathrm{Hz}$ in fig. 4 are in fact signal from resolved loops. As they can be several orders of magnitude stronger than $h_{\mathrm{n,c}}$, they need to be carefully excluded on the noise side to avoid skewing estimates of the Galactic confusion. As $\bar{h}$ is a slowly varying function of frequency, we compute $\bar{h}$ and $h_{\mathrm{n,c}}$ at a particular frequency bin by considering the neighborhood of $1\%$ of the total number of frequency bins. The quantities are then estimated from Gaussian fits to the histograms of the neighborhoods. This procedure eliminates the effects of signals located in the far-off tails of the histograms. The estimated $h_{\mathrm{n,c}}$ may appear low visually because of the great abundance of frequency bins at higher frequency, to the point that the true distribution of the frequency bins cannot be adequately resolved. The Galactic confusion for the simulation set is computed by averaging over all noise side realizations in the set.

We directly test the closest loops method by computing $h_{\mathrm{n,c}}$ and $\bar{h}$ for $G\mu=10^{-18}$ and $\alpha=0.1$, for the simulation sets including the closest $\left[10^{1},10^{2},10^{3},10^{4},10^{5}\right]$ loops, averaged over $[96,96,96,96,10]$ realizations in each set, respectively. We present an example at $f=0.05\ \mathrm{Hz}$ in fig. 5, because this is in the frequency range where the presence of signals is common. The 5 data points from lower to higher strain amplitudes corresponding to the growing number of loops included display a good linear relationship as required by eq. 37. By utilizing the closest loops method, we essentially extrapolate the line at each frequency to compute $h_{\mathrm{n,c}}$ for the entire Galaxy. This linear relationship holds down to a subset of just 10 loops and does start breaking down for even smaller subsets. Thus, our noise side simulation sets including 1000 loops are indeed sufficient for recovering the full Galactic $h_{\mathrm{n,c}}$.

Note that the cliff in fig. 4 is in a sense an artificial feature caused by including only a subset of loops. It therefore must go away if we correctly recover the full Galactic $h_{\mathrm{n,c}}$ through the closest loops method. This is indeed the case in fig. 6 showing the full $h_{\mathrm{n,c}}$ of the Galaxy with the rocket effect and $\alpha=0.1$ computed using the closest loops method. Clearly, the small breaks at $f_{\odot}$ are fully recovered, and the shapes of the curves resemble those in fig. 1 as expected.

The stochastic background of GW from loops can be computed by redshifting and integrating the power of GW emitted by all loops over all times, from the formation time of the earliest loops whose GW becomes relevant to the LISA frequency range, $t_{i}$DePies and Hogan (2007),

As a background energy density, it is customarily expressedDePies and Hogan (2007); Moore et al. (2015) as an energy density spectrum in terms of the critical density $\rho_{\mathrm{c}}=\frac{3H_{0}^{2}}{8\pi G}$,

The fact that eq. 38 concentrates the total power of emission in the fundamental mode is unimportant for our purpose, as the stochastic background is not significantly changed by GW emission spectraDePies and Hogan (2007). Higher frequency modes are incorporated when necessary.

Among the various quantities commonly used for expressing GW, the characteristic strain $h_{\mathrm{c}}$ is usually the most convenient to work with for data analysis, as it automatically incorporates the integration time over which the signal is observed. Eq. 39 can be converted usingMoore et al. (2015)

Since the stochastic background is really an unresolved GW background energy density, its actual strain amplitude cannot be simply defined, nor is that necessary as we always work with $h_{\mathrm{c}}$ in our data analysis. However, as we will see in section III.2.4, it is convenient practically to define an effective strain amplitude for the stochastic background

We plot $h_{\mathrm{n,u}}$ in fig. 7 with $\alpha=0.1$ and $T_{\mathrm{obs}}=1\ \mathrm{year}$. The change in slope indicates the transition into the high frequency flat region when plotted as $\Omega_{\mathrm{g}}(f)$, where loop decay is significant.

As a space-based GW detector, the orbiters of LISA experience considerable variations in their relative positions capable of mimicking GW signalAmaro-Seoane et al. (2017), setting the floor of LISA sensitivity. The level of signal necessary so that a signal-to-noise ratio (SNR) of unity can be obtained from the LISA observation is defined to be its sensitivity curve, which is then interpreted as a source of noise, i.e. the instrumental noise, for data analysis. We adopt the all-sky LISA instrumental noiseLarson et al. (2000) averaged over polarizations with a standard triangular detector configuration. This LISA instrumental noise accounts for most of the features in fig. 8 which is plotted as the effective strain amplitude $h_{\mathrm{n,d}}$ with $T_{\mathrm{obs}}=1\ \mathrm{year}$, including the overall trend and the ringing at high frequency. The same as in section III.2.2, $h_{\mathrm{n,d}}$ is defined in terms of $h_{\mathrm{c,d}}$ as

A major astrophysical source of noise in the LISA frequency range consists of the WD binaries population which is extremely abundant in the Galaxy, with an abundance in the halo estimated to beHiscock et al. (2000) $N_{\mathrm{WD}}\sim 10^{12}$ in the LISA frequency range. They are therefore an important source of GW in the Galaxy while remaining mostly unobservable otherwise. Far from the end of their inspirals, WD binaries in this frequency range are essentially monochromatic GW emitters, making them particularly relevant for our purpose. Their GW forms an unresolved confusion background primarily contributing to $f\lesssim 10^{-3}\ \mathrm{Hz}$ as can be seen in fig. 8. The fact that Galactic WD binaries mostly affect the low frequency region diminishes their relevance for our purpose, because of the relatively low loop abundance in that region.

The general optimal SNR for data analysis of GW detection isMaggiore (2007)

where $\tilde{h}(f)$ is the Fourier transform of the time-domain signal strain, and $S_{\mathrm{n}}(f)$ is the noise root spectral density. For monochromatic sources, the Dirac delta function in $\tilde{h}(f)$ kills the integral, and the SNR greatly simplifies toMaggiore (2007)

To proceed further, we turn our attention to $h_{\mathrm{c}}$ of monochromatic sources. Binary inspirals are quasi-monochromatic sources of GW withMoore et al. (2015); Finn and Thorne (2000)

where $N_{\mathrm{cyc}}$ is the number of wave cycles the binary inspiral emits GW in the frequency bin $f$, reflecting the limit on signal integration. For monochromatic sources, signal integration is instead limited by $T_{\mathrm{obs}}$, and therefore the correct identification of $h_{\mathrm{c}}$ is thenMoore et al. (2015)

Now the SNR for monochromatic sources can be derived by substituting eq. 46 into eq. 44 and converting $S_{\mathrm{n}}$ into $h_{\mathrm{c,n}}$,

The SNR for monochromatic sources is simply given by the ratio of the signal and the noise $h_{\mathrm{c}}$. Dividing both the numerator and the denominator by $N_{\mathrm{obs}}$, schematically, the SNR in our simulations is computed as

where the sum is over all sources of noise considered. As loops emitting harmonic GW are persistent sources in the Galaxy, we require a relatively low

for an event to be counted.

We present here results illustrating the validity of the closest loops method on the signal side. Snapshots from a realization in the signal side simulation set with the rocket effect at $G\mu=10^{-18}$ and $\alpha=0.1$ are plotted in fig. 9. For visual clarity, we select the realization with the lowest number of events satisfying the SNR requirement for this purpose. The total power of GW emission is concentrated in the fundamental mode to avoid unnecessary complications irrelevant to the central issue of this section. The upper panel reflects the normal method with which our simulation is performed and its results analyzed, where the closest 10 loops for each frequency bin are included in the signal side. There are 13 events in total. The data analysis procedure is repeated in the lower panel for the same realization, with the signal for each frequency bin coming from only the closest single loop. There are again 13 events in total, as expected.

A more conclusive validation can be obtained by comparing statistically results for all realizations in the simulation set. The histograms in fig. 10 cover the same simulation set with the rocket effect at $G\mu=10^{-18}$ and $\alpha=0.1$, with the difference that one (magenta) is generated by including contribution from the closest 10 loops for each frequency bin for the signal side, while the other (cyan) only includes signal from the closest single loop. Visually, it is clear that they are essentially identical, with the exception of a single event in a single realization which appears to be absent in the latter case. Statistically, the numbers of events from both histograms are also identical,

This comparison serves as a clear illustration that most events are in fact due to just a single loop being located extremely closely, and not from the combined contribution of GW emitted by multiple nearby loops.

The simulation set chosen is in the relevant region of the parameter space, where loop abundance in the Galaxy is very high, and the number of events is plentiful, ensuring the relevance and generality of the result obtained. Physically, the result means that even with such high loop abundance in the toy model, the probability that a single loop happens to be located sufficiently nearby to be resolved is still quite low, to the point that one has to be very lucky to have just one such occurrence. As the probability of multiple such occurrences in the same frequency bin decreases exponentially, tracking the closest 10 loops for each frequency bin for the signal side safely ensures that our simulation accounts for all potential event candidates in the Galaxy. Coupled with the noise side validation from section III.2.1, this result makes sure that despite constraints on computational capabilities, our simulation is a reasonable representation of loops in the Galaxy for the purpose of this study.

The computationally viable means of incorporating higher harmonic modes of GW emission from loops into our simulation is through postprocessing. The procedure is implemented from lower to higher frequencies, and consists of first renormalizing the fundamental strain amplitude according to eq. 11, then distributing the power into higher harmonics according to eqs. 15 and 17, and lastly computing the superposition of these higher harmonics with existing GW at those frequencies. This procedure avoids the burden of tracking higher harmonics for each loop, enabling us to incorporate different GW emission spectra easily and quickly.

A disadvantage of the above method is that the very low frequency region is missing contribution from higher harmonics of GW emitted by loops whose fundamental modes are below the frequency range, with the primary effect being a slight reduction in the Galactic confusion in this region. This is unimportant for our purpose for several reasons. First, the reduction is not very significant because the fundamental mode dominates in any GW emission spectrum, and after distributing power to higher harmonics, the fundamental only $h$ itself is only reduced by a factor of $\sim\sqrt{2}$. Second, $h_{n}$ decreases rapidly as $n$ becomes large, making the underestimation to the Galactic confusion very slight, as shown in ref. DePies and Hogan (2009). Third, this rapid decrease of $h_{n}$ means that any effect is only relevant for a narrow band in the low frequency range $1\times 10^{-5}\ \mathrm{Hz}\leq f\lesssim 5\times 10^{-5}\ \mathrm{Hz}$. As will become clear in section III.3.5, this low frequency band is irrelevant for the simulation sets in the important part of the parameter space, due to a combination of low loop abundance and the fact that in these sets, the instrumental noise with the Galactic WD binaries background greatly dominates over the Galactic confusion at those frequencies. One final argument for the unimportance of the missing contribution at very low frequency is the fact that the effects of higher harmonics themselves are actually minimal for the purpose of our simulation, as will become clear below.

An example illustrating any effects GW emission spectra may have on our simulation is presented in fig. 11, showing snapshots of the same realization selected in fig. 9. The top panel displays the signal when all the power of GW emission is concentrated in the fundamental mode, and is identical to the upper panel of fig. 9, with 13 events in total. The lower two panels display snapshots of the same realization when some power of GW emission is distributed in higher harmonics, with loops in the middle and bottom panels having the cusps- and the kinks- dominated spectra, respectively. We do not require coincident detection of higher harmonics of a fundamental mode signal here for a fair comparison. Higher harmonics of events whose fundamental modes have already been detected are carefully avoided, which are clearly visible at high frequency in the lower two panels. There are 11 and 12 events in total for signals with the cusps- and the kinks-dominated spectra, respectively. The slight differences in the total numbers of events for the three panels can be understood quite simply. The SNR is of course maximized when all the power is concentrated in the fundamental mode, and hence the highest number of events in that case. The SNR is reduced when some power is distributed in higher harmonics, and the effect is more pronounced for the shallower cusps-dominated spectrum, giving the lowest number of events.

A full comparison of the three cases for all realizations in the simulation set can be made statistically, with histograms shown in fig. 12 for the same set as that shown in fig. 10. Minor differences in the detailed features notwithstanding, the histograms for the three GW emission spectra have similar overall shapes and statistical properties, with the total numbers of events

Unsurprisingly, $E_{\mathrm{F}}$ is the highest as the SNR is maximized there. But the important point here is that GW emission spectra do not have significant effects on our simulation, which statistically is not at all sensitive to them when no coincident detection is required.

Since the harmonic signature is unique to GW from loops not replicated by any other sources, it is important to take advantage of it in the search of harmonic signal from resolved loops. As discussed in section III.2.3, WD binaries in the Galaxy are the most capable of causing confusion with loops. While the confusion background formed by these sources have been incorporated, with an abundance in the Galaxy as high as $N_{\mathrm{WD}}\sim 10^{12}$ in the LISA frequency rangeHiscock et al. (2000), they possibly have the potential to mimic the signal as well, should one happen to be located extremely nearby. While a detailed study of the possibility of detecting such signal from WD binaries would be beneficial for its own merit, it is not necessary for our purpose, because WD binaries are not harmonic sources. While they are also persistent sources of essentially monochromatic GW, they are unable to consistently mimic the harmonic signature of loops, and therefore can easily be excluded by requiring coincident detection of higher harmonics for candidate events. Even though $N_{\mathrm{WD}}$ appears high, our results (fig. 9) have already shown that even with a total number of loops in the Galaxy that is orders of magnitude higher in the toy model, the probability that a loop can be resolved out of even just the confusion of other loops in the Galaxy is still very low. This means that we would be lucky to individually resolve a few Galactic WD binaries, and the likelihood of multiple such sources consistently mimicking the harmonic signature is negligible.

While requiring coincident detection of higher harmonics can effectively eliminate false events caused by potentially confusing signal, such as that from Galactic WD binaries, genuine events should not be significantly affected as long as we do not require coincident detection of an excessive number of harmonic modes. We verify this by comparing statistical results for the same simulation set as that shown in fig. 10, with different coincident detection requirements. The histograms with the cusps- and the kinks-dominated GW emission spectra are presented in figs. 13 and 14, respectively, with the requirements for coincident detection set as no requirement (grey), at most the 1st harmonic mode (cyan), at most the 2nd harmonic mode (magenta). Histograms of the first case correspond to cyan and magenta histograms in fig. 12. The SNR requirement is modified according to the relative strengths of higher harmonics.

Minor differences in the detailed features notwithstanding, histograms in each figure all have similar shapes and statistical properties. It is somewhat noticeable that magenta histograms appear to have peaks shifted slightly to the left relative to the others. They give us the following event numbers

The results show that requiring coincident detection of the first harmonic appears to be the sweet spot, where we can reap the benefit of being able to eliminate false events without significantly affecting legitimate ones. We therefore adopt this requirement for the toy model unless otherwise specified, with the cusps-dominated GW emission spectrum, as cusps are generally dominant.

To analyze the effects varying $G\mu$ has on the results, we compare statistical results of event numbers obtained from histograms similar in nature to those presented earlier, for the simulation sets at $\alpha=0.1$ without the rocket effect and with $G\mu$ varying from $10^{-12}$ to $10^{-20}$. They are summarized in the second column in table 1. As will become clear in sections III.3.5 and III.3.6, this selection of simulation sets maximizes the total numbers of events, thereby giving the most informative results, while precluding complications not pertinent to the main issue discussed presently.

The most prominent feature of these statistical results is the relatively large number of events for $G\mu=10^{-12}$. Then the number of events just falls off a cliff as $G\mu$ decreases, before it starts increasing again and peaks at $G\mu=10^{-18}$, after which it again drops. To explain this feature, we note that a large number of events can be generated in two ways. The first way is to have strong GW emission from each loop. This effectively makes sources of noise other than the Galactic confusion less relevant, and the harmonic signal from a nearby loop can be detected without requiring it to be located in extreme proximity, as the signal only has to stand out over the Galactic confusion. Since the probability that a loop can overwhelm the Galactic confusion is not dependent on the strength of GW emission from each loop, and that $h\propto G\mu$ from eq. 23, this way favors higher $G\mu$ to ensure that other sources of noise are subdominant. The second way is to have high loop abundance in the Galaxy, which increases the likelihood that a loop happens to be present in extreme proximity that it overwhelms all sources of noise. From fig. 3, this way favors lower $G\mu$ up to a point. These two ways mean that the event number is a result of the tension between the strength of GW emission of each loop and loop abundance in the Galaxy, and that the two sides favor higher and lower $G\mu$, respectively. The ideal scenario for detection would of course be to have both. But fig. 3 tells us that we do not have that case, not even in the toy model. This tension explains features observed when $G\mu$ is varied.

At $G\mu=10^{-12}$, the Galactic confusion dominates the noise, and GW emission from each loop is the strongest. Thus, despite low loop abundance in the Galaxy, detecting harmonic signal from nearby loops is not very difficult as they just have to be nearby relative to the rest of the loop population at their frequency, and hence the large event number. As $G\mu$ decreases, the strength of GW emitted by each loop decreases, and other sources of noise take over. It now becomes impossible to detect harmonic signal from nearby loops that are just relatively nearby, while loop abundance is still far too low to enable an appreciable number of loops to be located in extreme proximity, and hence the very low event numbers at $G\mu=10^{-14}$ and $10^{-16}$. Loop abundance increases as $G\mu$ decreases further, contributing to a higher number of events and reaching a peak at $G\mu=10^{-18}$. Such loops with very low $G\mu$ decay slowly, with those in the LISA frequency range not significantly affected by loop decay. The resultant great loop abundance enables many to be located extremely closely, to the point that their harmonic signal becomes detectable despite the relatively weak GW emission. The number of events drops again as $G\mu$ decreases even further, due to less significant enhancement to loop abundance. At a given frequency, loop abundance does not depend on $G\mu$ intrinsically as can be seen from eq. 4, and stops increasing as $G\mu$ decreases once loop decay ceases being significant. This means that while the requirement of proximity for detection gets higher at $G\mu=10^{-19}$ and $10^{-20}$ due to weaker GW emission, these loops are not much more likely to be located more closely, lowering the event numbers. More importantly, the event number will only become even lower should $G\mu$ decrease more.

These results mean that when the rocket effect is not considered, $G\mu\lesssim 10^{-17}$ appears to be the best region in the parameter space, while $G\mu\sim 10^{-12}$ also seems viable.

As the rocket effect ejects all loops in the Galaxy beyond $r_{\mathrm{tr}}$, our expectation for its effects on simulation results is simply that no detection should be possible when $f>f_{\odot}$, because no loop can be located particularly closely as the entire population is bound by $r_{\mathrm{tr}}$. This strict elimination of signal should severely constrain detection for simulation sets with all but the lowest $G\mu$.

We verify this by analyzing statistical results of the simulation sets with the rocket effect, again at $\alpha=0.1$ so that they can be compared with results presented earlier, with $G\mu$ also varying from $10^{-12}$ to $10^{-20}$. They are summarized in the first column in table 1.

As expected, all events are eliminated by the rocket effect at $G\mu=10^{-12}$. This is especially significant given that the event number is high without the rocket effect, and is a direct result of the fact that here $f>f_{\odot}\forall f$ in the frequency range. For sets with lower $G\mu$, $f_{\odot}$ increases, and the impact of the rocket effect diminishes, as larger portions of the frequency range become viable for detection. We therefore start recovering significant numbers of events with $G\mu\lesssim 10^{-18}$. Recall from section II.3, that $f_{\odot}$ tracks the beginning of the decay regime as frequency increases. Thus, just as is the case for loop decay, the rocket effect also starts becoming irrelevant for $G\mu\lesssim 10^{-19}$. For this reason, event numbers at $G\mu\lesssim 10^{-19}$ are similar with and without the rocket effect. The small enhancement at $G\mu=10^{-20}$ with the rocket effect may be attributed to the fact that when $r_{\mathrm{tr}}\gtrsim r_{\odot}$, the rocket effect is actually beneficial to detection due to the reduction in the Galactic confusion.

The rocket effect eliminates $G\mu=10^{-12}$ as a viable region in the parameter space, and the sweet spot for detection is therefore $G\mu\lesssim 10^{-18}$, which should be the focus for physical results.

Eq. 4 appears to suggest that loop abundance should increase with a smaller $\alpha$, which is true for loops with a given $t_{\mathrm{c}}$. But when the frequency range is fixed, decreasing $\alpha$ increases $t_{\mathrm{c}}$, and actually reduces loop abundance because of the larger $H\left(t_{\mathrm{c}}\right)^{-3}$. This fact is clearly reflected by comparing loop abundance with $\alpha=10^{-5}$ shown in fig. 15 to that with $\alpha=0.1$ in fig. 3. The former is consistently lower than the latter by about two orders of magnitude. We again exclude the rocket effect here to eliminate complications from effects not pertinent to the issue discussed presently. Another feature to note from the comparison is the fact that the decay regimes remain the same, because they are really determined by the relative values of $\alpha t_{\mathrm{c}}$ and $\Gamma G\mu\left(t_{0}-t_{\mathrm{c}}\right)$. Recall from section II.3, the first quantity is independent of $\alpha$. For the second quantity, even though overall loops with $\alpha=10^{-5}$ are younger, the difference in $t_{\mathrm{c}}$ is negligible when compared to $t_{0}$. For the same reason, $f_{\odot}$ also remains similar and tracks the decay regime. The primary difference between the simulation sets with different $\alpha$ for our purpose manifests in differences in loop abundance in the Galaxy.

We now ask the question, for a given frequency bin, how does changing loop abundance in the Galaxy affect the probability of detection? A detailed answer is likely quite complicated, and is not essential to our study. Heuristically, there are three scenarios depending on the relative dominance of the Galactic confusion among all sources of noise. The first scenario applies when other sources of noise strongly dominate, and detections are very low probabilistic events, with any contribution being essentially due to the shell of closest loops, which has $N^{\frac{2}{3}}$. The probability of detection should be proportional to the mean scattering of the distribution of distances to these loops, which is a random walk with the mean proportional to $\sqrt{N^{\frac{2}{3}}}=N^{\frac{1}{3}}$. Since the overall noise level remains fixed as $N$ varies, we need to take into account the contribution from geometry, and the probability should also be inversely proportional to the average distance to these loops, i.e.

where $k$ is the frequency bin. The second scenario applies when other sources of noise only weakly dominate over or are comparable to the Galactic confusion. It is now easier for nearby loops to be detected, and instead of a shell, contribution comes from the volume of nearby loops. Thus, the scattering in the distribution of distances now scales as $\sqrt{N}$, and the probability of detection becomes

The third scenario applies when the Galactic confusion dominates. Detections are relatively easy, with contribution still coming from the nearby volume with scaling $\sqrt{N}$. Unlike the other two scenarios however, the overall noise now changes with $N$. Combining eqs. 35 and 36, the noise also scales as $\sqrt{N}$ modulated by geometry. We therefore disregard geometry, and see that we do not expect $p_{k}$ to depend strongly on loop abundance in this case.

We analyze effects of changing $\alpha$ on the simulation by computing statistical results for the simulation sets without the rocket effect and with $\alpha=10^{-5}$, summarized in the fourth column in table 1.

Though the simulation sets with $\alpha=10^{-5}$ have more events than those with $\alpha=0.1$ when $G\mu$ is high, by a factor of $\sim 3$, they are still consistent with the third scenario of the heuristic argument above, given that loop abundance differs by about two orders of magnitude. Moreover, as these sets contain relatively few loops, with only $N\sim\mathcal{O}(60)$ at each frequency when $G\mu=10^{-12}$, the Galactic confusion itself becomes less well-defined. In this case, conceivably it is relatively easy for a loop to become detectable by virtue of being located just somewhat more closely relative to the others, without actually being very close, thereby generating some minor enhancement in the event number.

Just as with results with $\alpha=0.1$, as $G\mu$ decreases, the event number first drops and then stages a rebound. However, here the rebound is far less significant in comparison, which comes as a direct result of the lower loop abundance. Thus, the only potentially viable region in the parameter space with $\alpha=10^{-5}$ is $G\mu\sim 10^{-12}$. We can already sense that when coupled with the rocket effect, this result implies that having smaller $\alpha$ is disadvantageous for detection.

We summarize in table 1 numbers of events for all 24 simulation sets of the toy model with $\Delta f$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{month}$. The third column shows that as expected, the rocket effect eliminates all events with $\alpha=10^{-5}$, including at $G\mu=10^{-12}$ which has a very high event number without the rocket effect. Recall from section II.3 that assuming perfect clumping of loops in the Galaxy is an overestimation for $\alpha=10^{-5}$, implying that detection is actually even less likely than what our results suggest. This also means that detection is not possible for loops created from the small-loop paradigm. Our attention should be focused on big loops, $\alpha=0.1$, where the hope of detection rests upon. The viable region of the parameter space for LISA detection of harmonic signal from individual loops is therefore $G\mu\lesssim 10^{-18}$ and $\alpha=0.1$.

One more thing to note is that despite the fact that the 24 simulation sets above have $\Delta f$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{month}$, the guidance they provide on the interesting region of the parameter space for further analysis should be applicable to our goal of making predictions for a 1-year LISA mission. To see this, we first note that in the limit when $p_{k}$ is low, directly decreasing $\Delta f$ without rebinning loop abundance just increases the number of events correspondingly. This fact is verified by the statistical result from the simulation set with the rocket effect at $G\mu=10^{-18}$ and $\alpha=0.1$, in the alternative toy model performed with $\Delta f$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{year}$, where the number of events is

computed without incorporating higher harmonics, and should therefore be compared with eq. 50. Clearly, the difference is consistent with a factor of 12. Recall from the description of the toy models at the beginning of this section, rebinning loop abundance in the Galaxy for each frequency bin from the toy models into the physical model results in reductions by orders of magnitude as shown in eq. 64, especially at higher frequency where most events are found. From the discussion in section III.3.6, for the simulation sets in the interesting region of the parameter space, the corresponding reduction in $p_{k}$ should be by a significant power-law of that factor which can be up to $\sim 10^{-7}$. This reduction more than offsets the enhancement shown above, and therefore results from the 24 sets in the toy model with $\Delta f$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{month}$ are already overestimations even for the physical simulation sets with $T_{\mathrm{obs}}=1\ \mathrm{year}$.

Arguments above already outline the basics for renormalizing results obtained from the toy models into those for the physical Galaxy. We now offer a concrete description of the renormalization procedure. The binning in $t_{\mathrm{c}}$ in eq. 4 translates into a binning in $f$,

When loop decay is not significant,

and eq. 61 simplifies to the expected relation

In the decay regime where the condition eq. 62 is not satisfied, natural frequency bins are no longer regular log bins, and can become enormous according to eq. 61, meaning that the already reduced population of loops in the Galaxy of the decay regime actually emits GW spanning a wide range of frequencies. The relation eq. 63 is largely applicable to all simulation sets in the interesting region of the parameter space from previous discussions. Then, rebinning from the toy models into the physical model reduces the number of loops in the frequency bin by a bin reduction factor

Recall from section III.1.3, a given frequency bin $i$ in the toy models can register either no event or one event with some small probability $p_{i}$. This means that the total number of events for a simulation set can be interpreted as the summation of outcomes from a collection of Bernoulli trials (random trials with two possible outcomes, such as flipping a coin), with probabilities $p_{i}$ and $1-p_{i}$ for outcomes 1 and 0, respectively. Though the actual $p_{i}$ is unknown, from eqs. 58 and 59, we know that $p_{i}$ is roughly reduced by a power-law of loop abundance which can be represented by the factor $\eta(f)^{\kappa}$ after rebinning, with the scaling parameter $\kappa=\frac{2}{3}$ when other sources of noise strongly dominate, and $\kappa=\frac{5}{6}$ when they only dominate weakly. From results of sets without the rocket effect and with $G\mu=10^{-18}$ for both values of $\alpha$, we see that the set with smaller $\alpha$ has a lower event number by a factor of $\sim 33$. Given that the difference in loop abundance is by about two orders of magnitude, we see that the actual simulation sets in fact contain a mixture of scenarios corresponding to different values of $\kappa$. We therefore compute predictions for physical detection probabilities for both cases, which should be interpreted as estimates of upper and lower bounds.

A set of $n$ Bernoulli trials follows the binomial distribution with the expectation value $np$. Adopting a frequentist interpretation, we can therefore effectively account for the rebinning without knowing $p_{i}$ by changing the success outcome for a frequency bin from 1 to $\eta(f)^{\kappa}$, effectively altering $n$. Then the predicted detection number for a simulation set in the physical model is

where the sum is over bins with events in the toy model. The more straightforward but less rigorous way of interpreting this relation is to think of $\eta\left(f_{i}\right)^{\kappa}$ as effective physical detection numbers for frequency bins.

Since events mostly occur at higher frequency, $\eta\left(f_{i}\right)$ is very small, and $N_{\mathrm{P}}$ is less than unity. It is therefore more suitable to interpret $N_{\mathrm{P}}$ as the total physical detection probability, $p_{\mathrm{det}}$, defined as the probability that at least one detection is made in the entire observed frequency range over the course of the mission. The more proper way of estimating $p_{\mathrm{det}}$ is

This interpretation with $\eta\left(f_{i}\right)^{\kappa}$ in the equation is well-defined in the limit $\eta\left(f_{i}\right)^{\kappa}\ll 1$ to the point that $N_{\mathrm{P}}$ is less than unity, because

after expansion and only keeping terms up to linear order in $\eta\left(f_{i}\right)^{\kappa}$. We should clarify that in writing down eq. 66, we are not taking $\eta\left(f_{i}\right)^{\kappa}$ to be the actual physical detection probability for the frequency bin, $p_{\mathrm{p},k}$. If the real $p_{\mathrm{p},k}$ is known, then the expression for $p_{\mathrm{det}}$ will be

where the index $k$ is over all physical frequency bins. Clearly, $\eta\left(f_{i}\right)^{\kappa}\gg p_{\mathrm{p},k}$, as otherwise

from eq. 64, and the condition for eq. 66 is not satisfied. In writing eq. 66 with $\eta\left(f_{i}\right)^{\kappa}$, we are in fact implicitly estimating $p_{\mathrm{p},k}$ through the overall statistics of events in the toy model modulated by $\eta\left(f_{i}\right)^{\kappa}$ as a result of rebinning. The more intuitive and less rigorous way of interpreting eq. 67 is simply that if the expected $N_{\mathrm{P}}$ is less than unity, say 0.1, then we reasonably expect to repeat the simulation set 10 times before a detection is made, i.e. the probability that a detection is made by performing the simulation set once can be said to be 0.1.

Eq. 66 assumes that the $T_{\mathrm{obs}}$ used to bin the simulation set in the toy model is the actual LISA mission duration, and is therefore only applicable to the single simulation set in the alternative toy model. Rebinning results from the toy model with $\Delta f$ corresponding to $T_{\mathrm{obs}}=1\ \mathrm{month}$ requires some further manipulation. To obtain a correct estimate, we first compute $\eta(f)$ with $T_{\mathrm{obs}}=1\ \mathrm{year}$, producing the correct $\eta(f)$ for the physical model. The problem now is just a lack of frequency bins, i.e. a lack of Bernoulli trials. The solution is to effectively duplicate each frequency bin 12 times in eq. 66. This procedure produces reasonable estimates of $p_{\mathrm{det}}$, because loop abundance in the Galaxy is a smooth and slowly varying function of frequency, and therefore $p_{\mathrm{p},k}$ and $\eta(f_{k})$ in nearby bins do not differ significantly. Then $p_{\mathrm{det}}$ can be estimated as