Accelerated binary black holes in globular clusters: forecasts and detectability in the era of space-based gravitational-wave detectors

The motion of the center of mass of a coalescing binary black hole (BBH) in a gravitational potential imprints a line-of-sight acceleration (LOSA) onto the emitted gravitational wave (GW) signal. The acceleration could be sufficiently large in dense stellar environments, such as globular clusters (GCs), to be detectable with next-generation space-based detectors. In this work, we use outputs of the \textsc{cluster monte carlo (cmc)} simulations of dense star clusters to forecast the distribution of detectable LOSAs in DECIGO and LISA eras. We study the effect of cluster properties -- metallicity, virial and galactocentric radii -- on the distribution of detectable accelerations, account for cosmologically-motivated distributions of cluster formation times, masses, and metallicities, and also incorporate the delay time between the formation of BBHs and their merger in our analysis. We find that larger metallicities provide a larger fraction of detectable accelerations by virtue of a greater abundance of relatively lighter BBHs, which allow a higher number of GW cycles in the detectable frequency band. Conversely, smaller metallicities result in fewer detections, most of which come from relatively more massive BBHs with fewer cycles but larger LOSAs. We similarly find correlations between the virial radii of the clusters and the fractions of detectable accelerations. Our work, therefore, provides an important science case for space-based GW detectors in the context of probing GC properties via the detection of LOSAs of merging BBHs.

The ∼ 90 BBH detections reported by the LIGO-Virgo-KAGRA collaboration (Abbott et al. 2021) have started to shed some light on their origin (Abbott et al. 2023).However, making precise inferences on formation channels from data needs to take into consideration two factors.The first is that detected binaries could be coming from a combination of the aforementioned formation channels.Indeed, the data suggest that multiple formation sub-channels even within isolated evolution contribute to this spectrum, although the extent of these contributions from different channels is unknown and not straightforward to constrain, in part because of the systematics associated with the population synthesis simulations (Wong et al. 2021;Zevin et al. 2021).The second is that in general, the shape of the GW waveform itself contains no definite signatures that can conclusively ascertain the provenance of the binary on a single-event basis In the case of BBH mergers assembled dynamically, the binaries move on orbits determined by the star cluster gravitational poten-tial.As this motion could leave an imprint on the GW signal in the form of a Doppler shift, its detection would contribute to our ability to identify the binary formation channel2 .However, a binary orbit at constant velocity would produce a constant Doppler shift in the GW waveform, degenerate with the mass of the binary.On the other hand, accelerated motion (with a non-zero component of the acceleration along the observer's line-of-sight) could modulate the signal and, therefore, be detectable (Yunes et al. 2011;Bonvin et al. 2017;Tamanini et al. 2020;Vĳaykumar et al. 2023).Constraints on this line-of-sight acceleration inferred directly from the GW signal could hence carry information on the environment in which the binary merged (Vĳaykumar et al. 2023).
GCs are among the dense stellar environments expected to efficiently assemble BBH mergers.They are stable, spherically symmetric, gravitationally bound collections of ∼ 10 4 −10 6 stars with typical sizes of ∼ 1 − 10 pc (Harris 1996;Gratton et al. 2019;Baumgardt & Vasiliev 2021).BBHs merging in GCs are expected to present an acceleration reminiscent of the environment in which they formed.Thus, detecting signatures of (time-varying) Doppler shift could not only point towards identifying different formation environments, but could also provide crucial information about masses, density profiles, metallicities, and ages of GCs.
In this work, we calculate accelerations of BBHs in GCs, as a function of the cluster properties, using the catalogue pertaining to the large-scale cluster monte carlo (cmc) (Kremer et al. 2020c) simulation.We extract and determine the accelerations of all the BBH binaries that merge within a Hubble time from the cmc catalogue, and employ a GW Fisher analysis3 (Cutler & Flanagan 1994) to estimate whether such accelerations can be sufficiently well constrained with the proposed DECIGO (Sato et al. 2017) and LISA (Danzmann & Rüdiger 2003) space-based detectors.We construct distributions of accelerations as a function of GC properties, with appropriately chosen detectability, metallicity, and cluster-mass weights.We then study the imprint of GC properties on these distributions.
The rest of the paper is organised as follows.Section 2 describes the cmc catalog models and outlines the prescription we use to construct distributions of BBH accelerations in GCs.Section 3 presents the results and, in particular, the imprint of GC properties on the distribution of accelerated BBHs.Section 4 summarizes the paper, discusses this work in the context of other GW probes of GCs, and suggests the scope for future work.In the entirety of the paper, we assume the standard cosmological model with parameters fixed to the Planck 2018 values (Aghanim et al. 2020).

The cmc models
The cmc catalogue comprises 144 simulations of GCs.It uses a Hénon type Monte Carlo algorithm which enables a long-term evolution of the GC (Hénon 1971b,a;Joshi et al. 2000Joshi et al. , 2001;;Fregeau et al. 2003;Fregeau & Rasio 2007;Chatterjee et al. 2010Chatterjee et al. , 2013;;Pattabiraman et al. 2013;Rodriguez et al. 2015), assuming a set of initial conditions.Details of the cmc simulation can be found in Kremer et al. 2020c,a.Here, we briefly summarize some of the most important features of the models.
A number of fixed initial conditions are assumed for the whole set of simulations.The initial cluster potential is assumed to follow a King profile (King 1962), with concentration parameter  0 = 5.The stellar masses are drawn from a Kroupa initial mass function (IMF; Kroupa 2001), assuming a mass range of 0.08M ⊙ − 150M ⊙ and the stellar binary fraction is set to  b = 5%.For binaries, the primary component is drawn from a Kroupa IMF, while the secondary component is chosen by drawing from a uniform distribution of mass ratios  ∈ [0.1, 1].The initial orbital period of binaries is drawn from a log-uniform distribution, with a lower limit on the separation set such that this separation () does not fall below five times the sum of the stellar radii of the binary ( ≥ 5( 1 +  2 )), and an upper limit set by the hard/soft boundary.Each simulation is evolved across 14 Gyr or until the GC undergoes tidal disruption (see eg. Heggie & Hut 2003) or collisional runaway (see eg. Portegies Zwart & McMillan 2002).
A number of physical processes have been incorporated into the cmc simulations.These include stellar and binary evolution, neutron star formation, black hole (BH) formation, modeling of strong encounters, two-body relaxation, three-body binary formation, implementation of galactic tides, and stellar collisions.We refer the reader to Kremer et al. 2020c for details on all these processes; we briefly summarize the prescriptions used for BH formation below.
BHs are modeled to form via standard iron core-collapse supernovae (CCSNe) using the "rapid model" for stellar remnants (Fryer et al. 2012).The CCSNe impart natal kicks to the BH, with mass fallback decreasing the magnitude of the kick.The kick velocities of neutron stars,  NS , are assumed to be described by a Maxwellian distribution with a dispersion set to  = 265km/s.For BHs, the kick magnitudes are then modulated as a function of the fallback mass fraction  b such that  BH = (1 −  b ) NS , where this fraction pertains to the percentage of stellar envelope mass that falls back onto the collapsed core.Additionally, pulsational pair-instability (Belczynski et al. 2016b) is implemented, which results in the mapping of stars with helium core masses in the range 45 − 65M ⊙ to BHs of masses in the vicinity of 40M ⊙ (Woosley 2017), producing an excess in that region of the BH mass spectrum.Stars with helium cores in excess of 65M ⊙ are modeled to produce no remnants at all (Heger et al. 2003).

Extracting accelerations from the cmc catalog
We describe below the prescription used to evaluate the accelerations of merging BBHs in GCs.A flowchart summary of the prescription we use below is provided in Figure 1.
(i) For each merger in the cmc catalog, we determine the mass of the cluster  enc enclosed within a radius , where  is the distance of the BBH from the centre of the cluster when it merges.
(ii) The acceleration of the center of mass of the BBH divided A flowchart summary of the prescription to extract accelerations from the cmc simulations, and construct weighted distributions of identifiable (found) and non-identifiable (missed) BBH accelerations (/).A standard FRW metric is assumed, with cosmological parameters taken from "Planck18" (Aghanim et al. 2020).CoM is the "centre of mass",  lb,cl is the lookback time of the GC at formation,  delay is the time to coalescence from the formation of the binary, which is assumed to coincide with the formation time of the cluster, and  lb,merg is the lookback time at merger.
by the speed of light /, is then evaluated as (Binney & Tremaine 1987;Bovy 2023): (iii) For each BBH and corresponding acceleration,  = 5000 redshift samples are drawn following the cosmic star-formation rate density (SFRD) as given in the Madau-Fragos prescription (Madau & Fragos 2017;Madau & Dickinson 2014): (2) These samples correspond to cluster-formation redshifts.In essence, we assume that the history of cluster formation follows that of stars.
(iv) To evaluate the merger epochs, we first convert the clusterformation redshifts to lookback time  lb,cl .Then, using the time-delay values from the simulation,  bbh,delay , the lookback time of the BBH at merger is calculated as  lb,bbh =  lb,cl −  bbh,delay .If positive, this lookback time is now converted back to a redshift at the merger.Otherwise, the sample is rejected, since it implies that the BBH will not merge within the age of the universe.
(v) Converting the redshift at merger to a luminosity distance at merger, and using the intrinsic parameters of the BBH provided from the simulation, the signal-to-noise ratio (SNR)  in DECIGO and LISA are calculated as: where h is the Fourier transform of the GW waveform as seen by (projected onto) the detector and modulated by /,  n (  ) is the detector's sky-averaged noise power spectral density (PSD), and where A ∝ M 5/6 / L with M,  L as the chirp mass and luminosity distance of binary respectively, Ψ(  ) is given by Eq. (3.18) of Buonanno et al. 2009, and ΔΨ(  ) is given by Eq. ( 4) of Vĳaykumar et al. 2023.The BBHs are assumed to have face-on (inner) circular orbits 5 .(vi) If  ≥ 10 (8) for DECIGO (LISA), the BBH is considered to be detectable, and the sample is kept.Otherwise, it is rejected.
(vii) Each detected BBH sample is assigned a set of weights: a cluster-mass weight,  cl , and a metallicity weight,   , to account for the relative cosmological abundance of clusters with different properties.We also assign a detectability weight  det .These weights are computed following Fragione & Banerjee 2021.The clustermass weight is assigned following the cluster initial mass function as (Portegies Zwart et al. 2010): (5) 4 The detector antenna pattern will change appreciably over the inspiral timescale.However, while calculating the SNR, we do not account for the effects of this time-varying detector antenna pattern for computational ease.We do not expect a significant change in the obtained SNR due to this effect. 5It is worth mentioning that in general, the inner binary is not necessarily circular at sub-Hz frequencies (Breivik et al. 2016).This is especially true in GCs where a non-trivial fraction of binaries could be eccentric ( where  cl is the mass of the cluster at formation.The metallicity weight   is assigned using lognormal distribution with a 0.5 dex standard deviation, and redshift-dependent mean given by (Madau & Fragos 2017): The detectability weight is given by: where the detection probability  det is accounted for by setting an SNR threshold and rejecting samples that do not exceed that threshold.As mentioned before, the threshold is 8 for LISA and 10 for DECIGO.To the surviving samples, we then assign a weight determined by the product of the cosmological time dilation piece 1/1 +  and the differential comoving volume  c /.
(viii) To determine if the acceleration of a sample BBH is constrainable, we resort to a Fisher Matrix Analysis (FMA) which approximates the shape of the GW parameter estimation likelihood to be Gaussian in the source parameters (Cutler & Flanagan 1994).From the corresponding covariance matrix, a statistical r.m.s.error Δ(/) is calculated, assuming the Gaussian is centered on true /.If / < Δ(/), / = 0 is contained within the 68% errorbar, and the BBH's acceleration is said to be "missed" ie. the event cannot be confidently identified as accelerating.If / > Δ(/), / = 0 lies outside the 68% errorbar and the BBH's acceleration is said to be "found".In other words, the event can be identified as accelerating at 68% CL.We briefly describe the application of the FMA to the identification of found-missed accelerations in Section 2.3.
(ix) We construct various histograms, including histograms of found and missed accelerations, weighted by  t , the product of the mass, metallicity, and detectability weights (Fragione & Banerjee 2021): We point out here that the BBH is assumed to be optimally oriented in a way that maximizes the SNR and the magnitude of the LOSA.The latter is also assumed to be unchanging.The inner orbit is assumed to be face-on, while the outer orbit (i.e., the orbit of the BBH's center of mass in the potential of the globular cluster) is assumed to be edge-on.The fractions of found accelerations in this work should therefore be considered as upper limits.
To assess the drop in the fraction of measurable accelerations due to a randomized orientation of the outer orbit, as well as due to a more stringent metric for measurability (2 − , 3 −  confidence), see Apendix A3.

Identifying found-missed BBH accelerations
A constant line-of-sight velocity component of the center of mass of a BBH will produce a constant Doppler shift that is degenerate with the mass of the BBH.On the other hand, a BBH with a LOSA will result in a time-varying Doppler shift, which in turn will modulate the GW waveform with respect to one that is not accelerated.At leading order, a deviation ΔΨ(  ) in the GW phase Ψ(  ) is incurred at −4 Post Newtonian (PN) order, and is given by (Bonvin et al. 2017): where   = (   / 3 ) 1/3 ,  is the total (detector frame) mass of the binary, and  is the symmetric mass ratio.Vĳaykumar et al. 2023 calculated 3.5 PN corrections beyond the leading order to ΔΨ(  ), and also showed that including these higher-order corrections is necessary for unbiased source property inference.We hence use the full expression of ΔΨ(  ) from Vĳaykumar et al. 2023 to construct our waveform approximant ℎ(  ).
To calculate the r.m.s error Δ(/), the Fisher matrix  is first constructed as (Cutler & Flanagan 1994): where  ,  are the binary's intrinsic and extrinsic parameters that determine the shape of the GW, and (|) represents a noise-weighted inner product between two GW waveforms (  ), (  ): The covariance matrix is then evaluated as C =  −1 , and Δ(/) is read-off (and square-rooted) from the corresponding diagonal element in C.
Intuitively, one would expect that a signal whose phase difference can be tracked for a longer time in-band would provide a better measurement of the acceleration.Since we choose to model the GW phase in the frequency domain, the tracking of the phase in time is equivalently converted to tracking the phase in frequency and is automatically taken into account by the choice of the bandpass of the detector (i.e. min ,  max ).Since time in-band  ≈  −5/3  −8/3 min , this also means that less massive events would have the least measurement uncertainty owing to their long in-band times.

RESULTS
In this section, we provide weighted distributions of found and missed accelerations, as well as the corresponding fractions with respect to the total number of detected BBHs.We also evaluate the fractions and the distributions as a function of the cluster properties (metallicity, galactocentric radius, and virial radius) in the cmc models.We restrict our attention to DECIGO and LISA detectors, which, by virtue of their sensitivity in the low-frequency regime, are especially suited to detect LOSAs from GCs.

Aggregate distributions of found and missed accelerations in DECIGO
There are two competing effects that determine if a BBH is detectable and its acceleration is found.BBHs with heavier masses produce GW signals with larger amplitudes and are therefore relatively easier to detect, although increasing the mass eventually reduces the detectability due to a smaller number of GW cycles in the detector frequency band.On the other hand, BBHs with lighter masses produce GWs with a smaller amplitude and are therefore relatively more difficult to detect out to large distances.Given that LOSA modulations to the GW are a low-frequency effect, incurring corrections in the GW phase at −4PN, longer durations of the in-band inspirals enable stronger constraints on the acceleration, or, equivalently, allow probes of smaller accelerations.BBHs with lighter masses spend a longer time in-band relative to BBHs with heavier masses, and thus contribute more significantly to the distribution of found accelerations.Moreover, the metallicity and cluster-mass weights also contribute to the fraction of found accelerations, as well as the shape of their distributions (see Appendix A).
We show, in Figure 2, the distribution of detectable (total), found, , corresponding to a fiducial GC enclosing 10 4 M ⊙ within a sphere of radius 10 −2 pc.About 12% of the accelerations are found in DECIGO.Moreover, the found accelerations peak at about 2 orders of magnitude larger than missed accelerations, as well as the total detectable accelerations, which is consistent with the modest fraction of accelerations that are found.Lighter BBHs provide stronger constraints on the accelerations (specifically, LOSAs), due to the increased number of cycles in-band.This is reflected in the distribution of found total masses having relatively less support at higher masses.The median of the found distribution is also smaller than the median of the missed distribution since the former is ∼ 61M ⊙ while the latter is ∼ 108M ⊙ .
and missed accelerations.Of the total detectable BBHs in DECIGO, 12% are found.The detectable accelerations follow a distribution that peaks between 10 −17 s −1 and 10 −16 s −1 with the median value being 1.7 × 10 −16 s −1 and 90% CI being [4.7 × 10 −18 , 4.8 × 10 −15 ]s −1 .The missed acceleration distribution peaks roughly at a similar value, having a median value: 8.5 × 10 −17 s −1 and 90% CI: [2.5×10 −18 , 1.8×10 −15 ]s −1 with relatively smaller support between 10 −15 s −1 and 10 −14 s −1 .Conversely, the found acceleration distribution peaks at ∼ 10 −15 s −1 , having a median value: 6.3 × 10 −16 s −1 and 90% CI: We depict, in Figure 3, the distributions of found and missed detector-frame total masses  det =  source (1 + ), where  is the cosmological redshift.The distribution of found  det is shifted towards smaller values relative to the corresponding missed distribution.This can be explained as follows.Lower-redshift mergers are lower-mass because of the mass-dependence of delay times and since mass segregation in GCs favors higher-mass mergers at early times (e.g., Chatterjee et al. 2017b;Fragione & Rasio 2023).Smaller masses then enable better constraints on / by virtue of spending more cycles in the detector band.
Similarly, Figure 4 gives the distributions of found and missed BBH redshifts.Again, as with  det , the distribution of found  is shifted towards smaller values relative to the corresponding missed distribution.Smaller  correspond to larger SNRs, which reduces the r.m.s error approximately as Δ(/) ∝ 1/.This allows relatively smaller accelerations to also be identified within 68% confidence.

Effect of GC properties on the distributions of found and missed acceleration
Properties of the GC determine the population of BBHs in the GC and its spatial distribution, and thus, by extension, the distribution of BBH accelerations.Here, we break down the effect of metallicity, virial radius, and galactocentric radius on the distribution of found and missed accelerations, and corresponding distributions of  det and  (outer orbital radius/cluster-centric radius).The cmc catalog encompasses 3 distinct GC metallicities:  = 2 × 10 −4 , 2 × 10 −3 , 2 × 10 −2 .We extract BBH accelerations and construct distributions (weighted by the metallicity, cluster mass, and detectability weights) pertaining to found accelerations for each of these metallicities.We find that the majority of found accelerations, 93%, come from relatively higher metallicity GCs,  = 0.02(29%), 0.002(64%), with the fraction dropping to 7% for  = 0.0002.
The found distributions of accelerations, detector-frame masses, and orbital radii, are shown in Figure 5.The left panel shows a systematic preference for higher accelerations with decreasing metallicity.This can be understood from the fact that GCs with a larger metallicity prefer forming at low redshift and have a relatively larger frac- .Weighted distributions of accelerations, detector frame masses  det , and outer orbital radii () of the found BBHs in GCs, with varying GC initial metallicity  in DECIGO.Larger  clusters have a relatively larger fraction of lighter BBHs.This is reflected in the fact that the majority (∼ 93%) of found accelerations come from relatively higher  clusters (left panel).Moreover, the lighter masses enable probes of smaller accelerations due to the increased number of cycles in-band.Heavier masses require larger accelerations to be found.This is reflected in the shift of found accelerations towards the larger values with decreasing .It is further corroborated by the distributions of detector frame masses  det (centre panel), where decreasing  pushes the distributions of found  det to larger values.The distributions of outer orbital radii  (right panel) are also consistent with the left and centre panels, where the distributions are peaking at systematically larger values with decreasing .Smaller  implies BBHs in such GCs need larger accelerations to be found, and therefore need to be closer to the centre of the potential.tion of low-mass BBHs.This decreases the detector-frame mass and enables measurements of smaller accelerations.The larger detectorframe mass BBHs in low-metallicity (and high-redshift) clusters, need larger accelerations to be confidently identified (found) as accelerating.This explanation is further corroborated by the corresponding distributions in the center panel, which show a systematic shift of  det distributions to larger values with decreasing metallicitythe medians being ∼ 34M ⊙ , ∼ 99M ⊙ , and ∼ 125M ⊙ in descending order of the metallicity.The right panel is also consistent with this picture since the distribution of outer orbital radii shifts to decreasing values with decreasing metallicity 6 .Smaller radii yield larger accelerations, which are required by heavier masses to be identified confidently (found) as accelerating.
We study the effect of changing the virial radius on the distributions of found accelerations and corresponding  det and .The cmc catalog provides 4 discrete values:   /pc = 0.5, 1.0, 2.0, 4.0.The effect of changing   on the distributions is less pronounced than the effect of changing .This can be explained as follows.For a given mass distribution and location of BBH mergers in a GC, a smaller   leads to more compact GCs, which in turn leads to larger accelerations.However, while there is a direct correlation between   and acceleration, there is no such correlation between   and BBH mass.Thus, while smaller   yield larger accelerations in general, they do not necessarily yield smaller  det which enables stronger constraints on /.Nevertheless, we do find that the fraction of found accelerations, varies markedly with decreasing   :   = 0.5pc (70%),   = 1.0pc (23%),   = 2.0pc (6%),   = 4.0pc (1%).
We additionally study the effect of varying  and   on the  distribution pertaining to found accelerations.This is shown in Figure 6.The left panel shows the  distribution getting progressively larger support at larger  values with decreasing metallicity .This can be readily explained in terms of the age of clusters.Lower metallicity GCs are older and thus reside at larger  values.Conversely, higher metallicity GCs are younger and contain a larger fraction of lower-6 See Appendix A for more details.
mass BBHs.This results in fewer samples of found accelerations at larger  -both due to reduced SNR as well as poorer acceleration constraints from higher-mass BBHs.The effect of   on  distributions is less pronounced (right panel of Figure 6), although larger accelerations from smaller values of   imply increasing support at larger redshifts.
We do not find any significant effect of changing   on the distributions or the fraction of found accelerations.This is due to the fact that the accelerations extracted from the cmc simulations consider only the potential of the GC and not the potential of the galaxy in which the GC is hosted.The center of mass of the GC itself will have an acceleration, which depends on   but has not been considered in this work.Adding this effect will likely cause a systematic shift in the acceleration distributions; however, we do not expect the distributions to be impacted significantly if the GC is situated at typical locations (  ∼ kpc) in the galaxy 7 .
We refer the reader to Appendix A for a more detailed explanation of how the application of metallicity and cluster-mass weights to the intrinsic distribution of found accelerations impact the variation of the fraction of BBHs with these accelerations as a function of cluster properties.

Distributions of found and missed accelerations in LISA
LISA's sensitivity band covers a frequency range that is lower than DECIGO:  ∈ [10 −4 , 1]Hz.The BBHs, therefore, spend a significantly longer time within the LISA band than the DECIGO band, which should enable stronger constraints on acceleration.However, LISA's sensitivity to stellar mass BBHs is much lower as compared to BBHs.As a result, the majority of the lighter BBHs are not detectable ( < 8) in LISA, given that the Madau-Fragos SFRD peaks at  ∼ 2. Nevertheless, among those BBHs that are detectable, ∼ 14%  .Weighted distributions of detectable found and missed accelerations in LISA (left panel; vertical axis shows counts in arbitrary units), as well as the effect of metallicity  on the found distributions (right panel).The fraction of found accelerations is ∼ 14%, and the acceleration distributions span similar values as for DECIGO (see Figure 2).The imprint of  on the distribution of found accelerations is different from what was found for DECIGO (see Figure 5). = 0.002 dominates the fraction of measurable accelerations.This can be explained as the result of competing effects.Higher  has a larger fraction of lighter BBHs, many of which are not detectable in LISA.On the other hand, lighter BBHs enable better constraints on acceleration.Among the discrete metallicities considered, the one closest to the "sweet spot" that has detectable BBHs with measurable accelerations is  = 0.002.
are found, in part because the lower frequency reach of LISA enables binaries to spend longer times in-band.
In Figure 7, we provide the distribution of found and missed accelerations (left panel), and the variation of found acceleration distributions with metallicity (right panel), whose imprint was found to be the most pronounced in the DECIGO analysis.We once again see that the found acceleration's distribution peaks between 10 −15 s −1 and 10 −14 s −1 with the median value being 1.2 × 10 −15 s −1 and 90% CI being [3.6 × 10 −17 , 1.7 × 10 −14 ]s −1 .However, unlike DECIGO, we find that the majority of these accelerations (93%) come from  = 0.002 clusters.This can be explained as the consequence of competing effects.GCs with higher metallicities have a larger fraction of lighter BBHs, many of which are undetectable with LISA.On the other hand, BBHs with lighter masses enable more precise acceleration measurements.Among the discrete metallicities considered, the metallicity value closest to the "sweet spot" that has both detectable and measurable accelerations is  = 0.002.It should be noted that the metallicity weight (and to a lesser extent the cluster-mass weight) also contributes to enhancing the fraction of found accelerations for  = 0.002 (see Appendix A).Correlations of metallicity with  det and  are similar though less pronounced than what was found for DECIGO, and are therefore not plotted.
Given that the sensitivities of other proposed millihertz spacebased detectors such as TianQin (Luo et al. 2016) are similar to LISA, we do not expect our forecasts to differ significantly for those detectors.

SUMMARY AND OUTLOOK
GCs are one class of dense stellar environments expected to host BBH mergers.The ∼ 90 events detected by the LIGO-Virgo-KAGRA network cannot conclusively determine if a given BBH was hosted by a GC, although the merger rate in GCs can be estimated by comparing against GC simulations (Rodriguez et al. 2021), or by calculating the fraction of the BBH population that is consistent with having isotropic spin directions (Fishbach & Fragione 2023).Such rates, however, are limited by the sample size of the detected BBHs, as well as uncertainties in the models of GCs and their initial properties (size, metallicity, etc).
On the other hand, LOSAs of BBHs leave an imprint on their GW waveform at -4PN, and can therefore be potentially constrained by detectors sensitive at low frequencies (e.g: decihertz, millihertz bands) such as DECIGO and LISA.BBHs in GCs are expected to contain finite LOSAs, and their distribution could contain imprints of the properties of the GCs.LOSAs could therefore assist in identifying the provenance of BBHs.
In this work, we forecast the distribution of detectable BBHs in GCs in DECIGO and LISA eras, that also produce accelerations that are identifiable (found) at ≥ 68% confidence.To do so, we use the outputs of the cmc catalogue to extract distributions of BBH accelerations, following the scheme presented in Figure 1.We summarize our main results below.
(i) We find that ∼ 12% (∼ 14%) of detectable BBHs in the DE-CIGO (LISA) era have accelerations that are well-constrained away from zero.We also find that the distribution of measurable (found) accelerations peaks at 10 −15 s −1 in DECIGO and between 10 −15 s −1 and 10 −14 s −1 in LISA.
(ii) Among found accelerations, the majority (∼ 93% in DECIGO and LISA) come from relatively higher metallicity ( = 2×10 −2 , 2× 10 −3 ) clusters.This is clearly reflected in the mass spectrum of BBHs with found accelerations.Higher metallicity clusters form at low redshift and have a larger fraction of relatively low-mass BBHs, thus enabling better measurements of acceleration.Conversely, low metallicity ( = 2×10 −4 ) results in a larger fraction of high (detector frame) mass BBHs, and their accelerations need to be 1 − 2 orders of magnitude larger to be found.In LISA,  = 0.002 dominates the fraction of measurable accelerations due to competing effects of lighter masses being more difficult to detect while also enabling more precise acceleration measurements.
(iii) We observe correlations between the virial radius   of the cluster and the shape of the distributions, although these are less pronounced compared to the correlations with metallicity.Nevertheless, the majority of the found accelerations come from small   (e.g.70% of found accelerations come from   = 0.5pc).We find no appreciable dependence of the fraction of identifiable accelerations on the galactocentric radius   , likely because the accelerations extracted from the cmc simulations do not account for the galactic potential that hosts the GC.
(iv) Converting the percentage of found accelerations to a rate of found accelerations in the DECIGO/LISA eras requires estimates of BBH merger rates out to redshifts  > 1, which to date is poorly constrained.We instead plot the fraction of found accelerations in DECIGO 8 as a function of redshift in Figure 8.This fraction initially decreases, reaches its minimum value in the redshift bin [5.5, 6], and starts rising slightly again.This rise coincides with the redshift ( ∼ 6) at which  = 0.0002 clusters overtake  = 0.002 clusters in their contribution to the total number of detected events (in part due to high   ; see Figure A2).Since the source-frame masses in  = 0.0002 clusters are slightly higher than those in  = 0.002 clusters, and events in low-metallicity clusters have higher acceleration owing to their relative closeness to the center 9 , the number of found events increases slightly in comparison to the number of missed events above  ∼ 6 10 .
We note that the results mentioned above and in the rest of the work are contingent on our modeling assumptions for   ,  cl , and  det .For instance, our understanding of cosmic GC formation history is incomplete, and the assumption that GC formation follows star formation might not be a good one.Semi-analytic models of GC formation built using dark matter halo merger trees (El-Badry et al. 2019) show that the cluster formation rate density peaks at a higher redshift ( ∼ 4) and does not track the SFRD.However, these estimates are themselves model-dependent, and we prefer to use an observation-oriented (ie. the Madau-Fragos SFRD) prescription in our work.While we only focus on model-dependent forecasts of LOSAs, measurements of LOSAs can also be used to constrain host GC properties of BBHs independently or in tandem with methods in Fishbach & Fragione 2023.Binaries in GCs, especially those that merge in the cluster cores, are expected to have non-negligible orbital eccentricities.Although we do not consider the effect of orbital eccentricity in our SNR calculations or Fisher forecasts, the number of mergers with found accelerations could increase if we include these effects (Xuan et al. 2023).
Other dense stellar environments that could host BBHs include nuclear star clusters (Hoang et al. 2018) and AGNs (Ford & McKernan 2022).As follow-up work, we plan to study the distributions 8 Given that the intrinsic rate of detectable stellar mass BBHs is expected to be small in the LISA era, we do not plot the corresponding evolution of the fraction of found accelerations with .All events with found acceleration in LISA lie at  ≲ 0.2. 9 This is due to lower natal kicks in low-metallicity environments (Kremer et al. 2020c).See also Appendix A1 for a related discussion. 10 The slight rise around  ∼ 1.5 can be similarly attributed to the redshift beyond which binaries in  = 0.002 clusters dominate over those in  = 0.02 clusters.
of accelerations of BBHs in these environments, and the imprints of their properties on said distributions.We also plan to compare distributions of accelerations coming from these different dense stellar environments, which, in principle, could help in determining the provenance of the BBHs.
The accelerations of BBHs extracted from the cmc simulations consider only the effect of the GC gravitational potential.However, encounters of BBHs with a third body, when they lie within the band of the detectors, could impart an acceleration that is significantly larger than those provided by the GC potential.Accelerations of such in-encounter mergers could therefore be detectable even by future ground-based detectors, such as the XG network.We plan to investigate this as well in future work.

APPENDIX A: EFFECT OF WEIGHTS ON THE DISTRIBUTION OF MEASURABLE ACCELERATIONS A1 Initial Cluster Mass Weight
To construct the distribution and determine the fraction of found (measurable) accelerations, we apply a weight  cl that is inversely proportional to the square of the initial cluster mass (cf.Eq. 5).The effect of applying this weight is to enhance the fraction of found accelerations pertaining to low metallicity mergers.This can be explained as follows: (i) High metallicity environments form low mass pre-supernova cores due to higher line-driven winds (Vink et al. 2001).
(ii) Low mass cores get a larger supernova natal kick owing to lesser mass fallback (Fryer et al. 2012).This high kick displaces them from the center of the cluster, i.e. to higher , possibly also ejecting them from the cluster in the process.
(iii) The only way to then have appreciable acceleration for high metallicity mergers is by having a very dense environment, i.e. clusters with a higher mass.
(iv) Since massive clusters are down-weighted by  cl ∼ 1/ 2 cl , the total number of high-metallicity mergers is also down-weighted.This is illustrated in the scatter plots of found BBHs in DECIGO.Figure A1 shows scatter plots of found accelerations vs corresponding radii for different metallicities, with and without accounting for  cl .

A2 Metallicity Weight
Another weight that is applied to the distribution of found accelerations is the metallicity weight.The weight is evaluated using a log-normal distribution in the metallicity whose mean is redshift dependent (Madau & Fragos 2017).Since the BBH redshifts are drawn following the Madau-Fragos SFRD, the metallicity weight is (broadly) a result of convolving this distribution with the log-normal distribution.
We plot metallicity weights for found samples as a function of redshift.We see that  = 0.002 has the largest weights between  = 1−4, in comparison to the other metallicities  = 0.02 and 0.002.Since the Madau-Fragos SFRD has the largest support between  = 1 − 4, metallicity weights tend to enhance the fraction of found accelerations for low metallicities (say  = 0.002), relative to the other metallicities.This in part explains the fractions displayed in Figure 5. Furthermore,  = 0.0002 has the largest weights only at  ≳ 7.5, where the Madau-Fragos SFRD has negligible support.On the other hand, where the SFRD has the largest support ( = 1 − 4), this metallicity value has the smallest weight.This explains, in part, the small fraction of found accelerations assigned to  = 0.0002.

A3 Impact of randomized orbital inclination and stringent measurability criteria
The analysis presented in this paper assumes the acceleration vector to be aligned with the line of sight.The resulting fractions of found accelerations should therefore be thought of as upper limits.To assess the reduction in this fraction from a more realistic set-up, we allow the angle  with respect to the line of sight to vary uniformly in cos .We find that the fraction of measurable accelerations reduces to ∼ 6% in DECIGO and ∼ 7% in LISA (see Figure A3).We also present results for accelerations measurable at 2 −  and 3 −  confidence (see Figure A4).We find that the fraction of found accelerations drops to ∼ 7% and ∼ 5% for the 2 −  and 3 −  detections, respectively, in DECIGO while the same fractions in LISA drop to ∼ 8% and ∼ 5%.
Figure1.A flowchart summary of the prescription to extract accelerations from the cmc simulations, and construct weighted distributions of identifiable (found) and non-identifiable (missed) BBH accelerations (/).A standard FRW metric is assumed, with cosmological parameters taken from "Planck18"(Aghanim et al. 2020).CoM is the "centre of mass",  lb,cl is the lookback time of the GC at formation,  delay is the time to coalescence from the formation of the binary, which is assumed to coincide with the formation time of the cluster, and  lb,merg is the lookback time at merger.
min ,  max are frequency limits set by the detector bandwidth 4 .To choose the frequency limits for the mergers, assuming an observation time of 4 years, we follow Berti et al. (2005) with {  min ,  max } being {10 −2 , 10}Hz for DECIGO and {10 −4 , 1}Hz for LISA.To model h, we use the TaylorF2 prescription given by:

Figure 2 .
Figure2.Weighted distributions of detectable, found, and missed accelerations of BBHs in GCs.The dotted line in the histogram represents an acceleration, of magnitude 4.65 × 10 −14 s −1 , corresponding to a fiducial GC enclosing 10 4 M ⊙ within a sphere of radius 10 −2 pc.About 12% of the accelerations are found in DECIGO.Moreover, the found accelerations peak at about 2 orders of magnitude larger than missed accelerations, as well as the total detectable accelerations, which is consistent with the modest fraction of accelerations that are found.

Figure 3 .
Figure3.Weighted distributions of detector frame total mass ( det ) of detectable BBHs with found and missed accelerations in DECIGO.The vertical axis shows counts in arbitrary units.Lighter BBHs provide stronger constraints on the accelerations (specifically, LOSAs), due to the increased number of cycles in-band.This is reflected in the distribution of found total masses having relatively less support at higher masses.The median of the found distribution is also smaller than the median of the missed distribution since the former is ∼ 61M ⊙ while the latter is ∼ 108M ⊙ .

Figure 4 .
Figure 4. Weighted distributions ( vertical axis shows counts in arbitrary units) of redshifts of detectable BBHs with found and missed accelerations in DECIGO.The redshift distribution of found BBHs peaks at smaller values compared to the missed redshift distributions.This is because the r.m.s error on the acceleration estimate scales inversely with SNR, which in turn scales inversely with luminosity distance  L .
Figure5.Weighted distributions of accelerations, detector frame masses  det , and outer orbital radii () of the found BBHs in GCs, with varying GC initial metallicity  in DECIGO.Larger  clusters have a relatively larger fraction of lighter BBHs.This is reflected in the fact that the majority (∼ 93%) of found accelerations come from relatively higher  clusters (left panel).Moreover, the lighter masses enable probes of smaller accelerations due to the increased number of cycles in-band.Heavier masses require larger accelerations to be found.This is reflected in the shift of found accelerations towards the larger values with decreasing .It is further corroborated by the distributions of detector frame masses  det (centre panel), where decreasing  pushes the distributions of found  det to larger values.The distributions of outer orbital radii  (right panel) are also consistent with the left and centre panels, where the distributions are peaking at systematically larger values with decreasing .Smaller  implies BBHs in such GCs need larger accelerations to be found, and therefore need to be closer to the centre of the potential.

7
See Figure 5 of Vĳaykumar et al. 2023 for an estimate of acceleration due to the gravitational potential of a Milky Way-like galaxy.

Figure 6 .
Figure6.Weighted distributions of redshifts (s) of found BBHs in DECIGO as a function of metallicity  (left panel) and virial radius   (right panel).GCs with larger  are younger (i.e., they are formed at lower ).Moreover, they have a larger fraction of lighter BBHs, which are only detectable at lower s.On the other hand, smaller  clusters are older and have a relatively larger fraction of heavier BBHs, which are detectable at higher .Both of these contribute to redshift distributions in low  GCs having support at large .On the other hand, GCs with smaller   have BBHs with relatively larger accelerations, which enables them to be found at larger .This explains the systematic shift of the redshift distribution of found accelerations towards larger values with decreasing   .

Figure 8 .
Figure 8. Variation of the fraction of found accelerations with redshift in DECIGO.Each fraction corresponds is calculated in a redshift bin of width 0.5, with each point sitting at the centre of that bin.The found fraction initially decreases, reaches its minimum value in the redshift bin[5.5, 6], and starts slightly rising again due to increasing metallicity weight.

Figure A1 .Figure A2 .
Figure A1.Scatter plots of detectable accelerations vs corresponding distance of the merger from the center of the cluster, colored by metallicities.The left panel shows mergers in clusters of different metallicities when 2000 mergers (with detectable accelerations) are drawn randomly without any weights, while the right panel shows the same when mergers are drawn following the initial cluster mass weight  cl .Notably, as explained in section A1, binaries in high metallicity clusters are farther away from the center on average and the number of binaries detected in high metallicity environments decreases when the initial cluster mass weight is applied.In particular,  = 0.02 binaries make up 46% of the total number in the left panel as opposed to 27% in the right panel.

Figure A3 .Figure A4 .
Figure A3.The left panel shows a comparison between the distributions of found accelerations (absolute values) for the aligned orientation (AO) and random orientation (RO) cases in DECIGO, while the right panel shows the same in LISA.We observe that the found fraction drops to ∼ 6% and ∼ 7% in DECIGO and LISA, respectively.