© The Authors 2025

1. Introduction

The first direct detection of the gravitational wave (GW) source GW150914 opened an era of GW astronomy and added a new dimension to multi-messenger astrophysics (Abbott et al. 2016). Inspiraling binary neutron stars (BNS) are promising sources to be used as standard sirens (Abbott et al. 2017a), provided their host galaxies can be identified and their associated redshift are measured (Schutz 1986). More importantly, as GW detectors are operational and gathering data, it will be possible to test various aspects of General Relativity (GR) in a ways inaccessible to other techniques. For example, the speed of GWs has been measured using the time delay among GW detectors (Cornish et al. 2017) and the time delay between GW and electromagnetic observations (Abbott et al. 2017b). However, alternative theories of gravity predict v_g ≠ c due to the breaking of the weak equivalence principle or the existence of massive gravitons (Will 2005). In this Letter, we focus on a new approach to directly measuring v_g with the DECi-hertz Interferometer Gravitational-wave Observatory (DECIGO), a proposed Japanese space mission based on laser interferometer space satellites (Kawamura et al. 2011). We stress that the 0.1–10 Hz frequency band covered by DECIGO will fill the gap between LISA and ground-based detectors, thereby expanding the reach of nascent GW astronomy (Yagi et al. 2012; Cao et al. 2021, 2022). In particular, it will be possible to detect and follow the last stages of BNSs during the low-frequency adiabatic inspiral phase long before their final coalescence. This opens up the possibility for measuring the speed of GWs in the distant universe, especially in the 0.1–1 Hz band, which DECIGO is most sensitive to. In the context of the accelerating expansion of the Universe, effective field theories (EFT) provide a powerful theoretical framework to incorporate an additional scalar field to describe the dynamics of the universe and its possible couplings to gravity (Chen et al. 2015). In some EFT models (Joyce et al. 2015), transient deviations of v_g from the speed of light are allowed at frequencies well below ground-based detectors. Note that the massive graviton scenario can also lead to a frequency-dependent GW propagation speed (Ezquiaga et al. 2021), which is testable in any theory of gravity where the spectral dimension of spacetime changes with the probed scale (Yunes et al. 2016). The advantages of our method are the following. (I) It offers the first empirical assessment of the speed of GWs at times much earlier than the present. (II) Testing the consistency of v_g measurements across different redshifts and sky positions would be an important test of GR. (III) Measuring v_g(f_low) at different frequencies, which is nearly impossible at the present stage, would be of paramount importance for probing the landscape of dark energy and extended gravitational models. Throughout this analysis, we denote the speed of light by c and adopt a notation, c₀ = 299792.458 km s⁻¹, for its laboratory-measured value, which is also supported by recent astrophysical observations (Cao et al. 2017, 2020; Liu et al. 2021).

2. Methodology and simulation

Our proposed method, based on the theoretical framework described in Baker et al. (2022), relies on the ability of future GW interferometers to measure the luminosity distance D_L(z) and the expansion rate H(z). From the D_L(z), its derivative, and the expansion rate, it is possible to calculate the speed of GWs (v_g) (see Appendix A). DECIGO is the only proposed future detector that could achieve such a scientific goal (Seto et al. 2001). Compared to other ground-based and space-based GW detectors (Hou et al. 2025), the advantage of DECIGO lies in its ability to register a much larger number of GW cycles from BNSs, even at a redshift of z ∼ 5 (Kawamura et al. 2019). This would enable the discovery of a larger number of BNSs in their inspiral phase long before entering the frequency range of Laser Interferometer Gravitational-Wave Observatory (LIGO) (Kawamura et al. 2019). Therefore, the signal-to-noise ratio (S/N) of DECIGO is much higher than that of current ground-based GW detectors, so the measurement uncertainty of D_L(z) will be considerably reduced. In addition to the luminosity distance, we still need to assess the expansion rate, and the redshift-drift measurement meets our demands. It is given by the following expression: ${\mathrm{\Delta }}_{t}z=H_0\mathrm{\Delta }{t}_0(1+z-\frac{H(z)}{H_0})$, where Δt₀ and H₀ denote the observation period and the Hubble constant, respectively (Loeb 1998). To follow, the Hubble parameter at redshift z can be written as

$$\begin{array}{c}\hfill H(z)=(1+z)H_0-\frac{{\mathrm{\Delta }}_{t}z}{\mathrm{\Delta }{t}_0}.\end{array}$$(1)

The order of magnitude of the redshift drift is roughly given as the observation time divided by the cosmic age, which makes it difficult to measure Δ_tz with current technology (Yoo et al. 2011). On the other hand, DECIGO is the only proposed detector that can precisely measure a GW phase correction Ψ_acc(f) due to cosmic acceleration. The Fourier transform of the GW waveform from a coalescing binary system of masses m₁ and m₂ can be expressed as

$$\begin{array}{c}\hfill \stackrel{\sim }{h}(f){|}_{\text{no accel}}=\frac{\sqrt{3}}{2}\mathcal{A}{f}^{-7/6}{e}^{i\mathrm{\Psi }(f)}\left[\frac{5}{4}{A}_{\mathrm{pol},\alpha }(t(f))\right]e^{-i({\phi }_{\mathrm{pol},\alpha }+{\phi }_D)}.\end{array}$$(2)

The amplitude of the GW is

$$\begin{array}{c}\hfill \mathcal{A}=\frac{1}{\sqrt{30}{\pi }^{2/3}}\frac{{{\mathcal{M}}_z}^{5/6}}{D_L},\end{array}$$(3)

where the redshifted chirp mass, ℳ_z = (1 + z)η^3/5M_t, is related to the total mass, M_t = m₁ + m₂, and the symmetric mass ratio, η = m₁m₂/M_t². The function Ψ(f) represents the frequency-dependent phase arising from the orbital evolution, which can be given by the second post-Newtonian (2-PN) approximation (Maggiore 2008b). The polarisation amplitude A_pol, α(t), the polarisation phases φ_pol, α(t) (α represents the number of individual detectors), and the Doppler phase φ_D(t) are explicitly given in Yagi & Tanaka (2010).

One of the main objectives of DECIGO would be to directly measure the accelerated expansion of the Universe (Yagi et al. 2012). This is possible by measuring GWs coming from BNSs at far distances, since the phase of the waveform can shift if the sources recede due to acceleration. Let us first derive the correction to the GW phase caused by the accelerating expansion of the universe. We define h(Δt) as the observed GW waveform, where Δt = t_c − t denotes the time to coalescence measured in the observer’s frame with t_c representing the coalescence time. In our method, time to coalescence is equivalent to Δt_o used in Eq. (1). To follow, Δt can be related to ΔT ≡ (1 + z_c)Δt_e as (Takahashi & Nakamura 2005)

$$\begin{array}{c}\hfill \mathrm{\Delta }t=\mathrm{\Delta }T+X(z_c)\mathrm{\Delta }{T}^2,\end{array}$$(4)

where z_c is the source redshift at coalescence, Δt_e is the time to coalescence measured in the source frame, and X(z) is the acceleration parameter defined as $X(z)\equiv \frac{H_0}{2}(1-\frac{H(z)}{(1+z)H_0})$. In Friedmann-Lemaître-Robertson-Walker (FLRW) space-time, X(z) is related to the redshift drift expressed as Δ_tz = 2(1 + z)Δt_oX(z). The Fourier component of this waveform is

$$\begin{array}{c}\hfill \stackrel{\sim }{h}(f)=e^{i{\mathrm{\Psi }}_{\text{acc}}(f)}\stackrel{\sim }{h}{(f)|}_{\text{no accel}},\end{array}$$(5)

with

$$\begin{array}{cc}\hfill {\mathrm{\Psi }}_{\text{acc}}(f)& =-2\pi fX(z_c)\mathrm{\Delta }T{(f)}^2=-{\mathrm{\Psi }}_N(f)\frac{25}{768}X(z_c){\mathcal{M}}_z{x}^{-4},\hfill \end{array}$$(6)

where x = (πℳ_zf)^2/3, ${\mathrm{\Psi }}_N(f)=\frac{3}{128}{(\pi {\mathcal{M}}_{z}f)}^{-5/3}$ and $\stackrel{\sim }{h}{(f)|}_{\text{no}\text{accel}}$ corresponds to the gravitational waveform in the Fourier domain without cosmic acceleration. A term proportional to x⁻ⁿ represents the n-th PN order relative to the leading Newtonian phase Ψ_N(f); hence, this is a ‘4-PN’ correction. This means that lower-frequency GWs detectable with DECIGO are advantageous (Quartin & Amendola 2010). The measurement accuracy of X(z) has been estimated in Seto et al. (2001), using the so-called ultimate DECIGO, which is three orders of magnitude more sensitive than the realistic DECIGO detector (Takahashi & Nakamura 2005). The influence of the line-of-sight contamination could affect the phase shift measurement due to a weak lensing effect, especially for high-redshift GWs (Takahashi 2004). Therefore, we limit the following analysis to the redshift range z ≲ 2.0. We estimated the measurement accuracy of the parameter X_H ≡ X(z)/H₀, using deep learning techniques, based on the one-side noise power spectral density S_h(f) characterising the performance of the realistic GW detector (Yagi et al. 2011; Kawamura et al. 2019). The crucial part of our idea is to perform separate analyses in the low- and high-frequency bands for each BNS observed at redshift z. Let us now describe the division between these bands in more detail. Regarding the low-frequency band, the initial frequency registered, i.e. the frequency at Δt_o before the coalescence, is fixed at f_in = (256/5)^−3/8π⁻¹ℳ_z^−5/8Δt_o^−3/8, while the final frequency of this band is fixed at f_fin = 1/(1 + z)Hz, in order to guarantee a constant value of v_g(f) from the source to the observer (see Eq. (13) in Appendix A). In the high-frequency band, the initial frequency is fixed at f_in = 1 Hz, while the final cut-off frequency is taken as f_fin = 100 Hz, DECIGO’s higher frequency (Yagi & Tanaka 2010).

The waveform analysis allows us to fit the following parameter vector $\theta =(ln{\mathcal{M}}_z,ln\eta ,\beta ,t_c,{\varphi }_c,{\theta }_S,{\varphi }_S,{\theta }_L,{\varphi }_L,D_L,X_H,{\mathit{\theta }}_{\mathit{MG}})$, where β denotes the spin-orbit coupling parameter and ϕ_c is the coalescence phase. (θ_S, ϕ_S) and (θ_L, ϕ_L) denote, respectively, the direction of the source and the orientation of its orbital axis in the barycentric frame. In our simulations we set the fiducial values of m₁ = m₂ = 1.4 M_⊙ and t_c = ϕ_c = β = 0 without losing generality (Yagi & Tanaka 2010). Considering DECIGO’s angular resolution (∼arcsec²), which can uniquely identify the host galaxy of the binary, it is possible to determine the angular position (θ_S, ϕ_S) by pointing telescopes towards the location uncertainty box at the expected coalescence time from the chirp signal (Kawamura et al. 2019). We perform a Monte Carlo simulation by randomly distributing the orientation of sources (θ_L, ϕ_L). The effect of modified gravity θ_MG is quantified by two parameters: the transition frequency f_trans and the height of the v_g transition δv_g = v_g(f_high)−v_g(f_low). In our analysis the exact location of the transition frequency is not included as a free parameter to be fitted, provided that it lies between the low- and high-frequency bands of DECIGO (Baker et al. 2023). The uncertainties of D_L(z) and H(z) were estimated using a convolutional neural network (CNN), the network structure of which has been presented in detail in Sun et al. (2024, see our Appendix B). Since DECIGO is expected to detect BNS mergers at a rate of ∼10⁴ yr⁻¹ (Seto et al. 2001), two observational strategies were applied in the simulations: (Strategy I) a one-month observation of 10³ BNS, and (Strategy II) a one-year observation of 10⁴ BNS. The merger rate of double compact objects in a cosmological scenario is quantified by the conservative star formation rate (SFR) function (Dominik et al. 2013). We employed Gaussian processes (GP) to extract information from D_L(z) and its related errors σ_DL, which yield numerically reconstructed D_L(z) and D′_L(z) as smooth analytical functions (Seikel et al. 2012). Then, we evaluated the GP output function on the individual H(z) measurements, with the output of v_g(z) at different redshifts.

3. Numerical results and implication

We applied the techniques outlined above to the two observational strategies of DECIGO, focusing on the low-frequency and high-frequency GW signals from the BNS. A fiducial flat ΛCDM cosmology is applied, with Ω_m = 0.30 and H₀ = 70 km/s/Mpc. Our results show that tighter constraints on the speed of GWs are obtained in the low-frequency band, since the effect due to cosmic acceleration is greater on GW signals with lower frequencies (Seto et al. 2001; Kawamura et al. 2019). In order to analyse the performance of our method, we considered two scenarios, where the height of the transition δv_g(f_low)/c₀ = 10⁻⁴ and δv_g(f_low)/c₀ = 10⁻⁵, respectively, occurring at f = 1Hz (Baker et al. 2023). The results for both observational strategies are presented in Fig. 1. As shown, DECIGO will provide reasonably precise measurements of GW speeds from sources located at different redshifts. By benefiting from the accumulation of a large number of inspiral cycles in registered GW signals, which enables high-precision measurements of the phase shift associated with the redshift drift, such a space-based GW detector will more precisely measure D_L(z) and H(z), resulting in better measurements of v_g(f_low) at different redshifts. Disagreement between such individual v_g(z) measurements, if statistically significant, could be a probe of GR, considering modifications of the group velocity of GWs generated by Lorentz invariance violation in the gravity sector of the gravitational Standard-Model Extension (Kostelecký & Mewes 2016).

Fig. 1. Measurements of the speed of GWs based on a one-month (upper panel) and a one-year (lower panel) observation of BNS events by DECIGO. The scale of departure from the GR δvg(flow)/c0 is assumed to be 10−4 and 10−5, respectively.

When we summarised multiple v_g(f_low) measurements through the inverse variance weighting, the final constraint with 1σ uncertainty is v_g(f_low)/c₀ = 0.9997 ± 0.0005 for strategy I. This means that the one-month long observation enabled testing deviations of δv_g(f_low)/c₀ from zero (as predicted by GR) with a precision of 5 × 10⁻⁴. Hence, the fiducial transition height δv_g(f_low)/c₀ = 10⁻⁴ underlying the simulations of strategy I could be detected. It is also clear that a one-month observation by DECIGO would not be sufficient to detect a δv_g(f_low)/c₀ = 10⁻⁵ deviation, even if the GR were not the correct theory of gravity. We then proceeded by exploring whether a one-year observation of DECIGO could perform better. Summarising multiple v_g(f_low) measurements through the inverse variance weighting, the final constraint with 1σ uncertainty results in v_g(f_low)/c₀ = 0.99997 ± 0.00005. Therefore, if the speed GWs at low frequencies differs from c₀ by a factor of 10⁻⁵, DECIGO will succeed in detecting the signal of GR failure through a one-year observation of BNS events. Figure 2 gives us insight into the performance of our method and its relation to alternative methods of direct measurements of the speed of GWs. It should be stressed that we consider the propagation of GWs in an expanding universe under the far-field approximation. This assumption is equally valid for both f_low and f_high frequency limits relevant to our study. If the graviton is massive, it may cause backreaction on the source, which in turn could imprint timing uncertainties in measuring the v_g variation during later propagation (de Rham et al. 2011). However, the uncertainty or lacking information in the v_g variation near the source will not affect the conclusions derived from our simulations (see Appendix C for further discussion).

Fig. 2. Posterior distributions of vg determined by different techniques. DECIGO data under observing strategies I (one-month) and II (one-year) are compared with two other methods: a direct method based on transit time across a separated network of detectors, and an indirect method using the time delay of GW170817. Two strategies were considered: LIGO O2 – data from O2 of Advanced LIGO-Virgo – and LIGO O2 fixed – where the sky localization of the BNS source is fixed. Note that LIGO O2 and GW170817 only provide constraints on vg/c, and at low redshifts, we take c = c0 in this figure.

It would be appropriate to compare the above bounds with the recent observational studies of GWs. For example, Liu et al. (2020) proposed a method of determining the speed of GWs by measuring the transit time across a geographically separated network of detectors. By combing ten binary black hole (BBH) events and the BNS event from the second observing run of Advanced LIGO and Advanced Virgo, they constrained the speed of GWs as v_g(f_high)/c₀ to (0.97c, 1.05c). Such speed was narrowed down to (0.97c, 1.01c) by fixing the sky localization of the BNS source at the electromagnetic (EM) counterpart. This is actually another way to indirectly measure the speed of GWs. As clearly seen from the v_g assessment comparison between DECIGO and Advanced LIGO + Advanced Virgo, the second generation space-based GW detector will result in more stringent constraints on v_g. Moreover, very strong constraints on the difference between the speed of gravity and the speed of light (−3 × 10⁻¹⁵ < δtv_g(f_high)/c < +7 × 10⁻¹⁶) have also been obtained using the time delay between GW and EM observations of GW170817 (Abbott et al. 2017b). Since this result exceeds the yields of our method by orders of magnitude, some comments are in order. First, the bound relies on the particular interpretation of the time delay Δt = 1.74 ± 0.05s between the coalescence moment registered in GW band and the peak of gamma-ray fluence (Ciolfi & Siegel 2014). Relaxing model assumptions about the EM emission, or including certain exotic scenarios for the time difference could lead reduce the precision by two orders of magnitude (Abbott et al. 2017b). Moreover, this indirect method, which measures v_g relative to c at a specific redshift, depends on the prior assumption that c is constant. These effects would significantly reduce the effectiveness of time-delay based measurements, assuming much smaller numbers of associated BNS-GRB events probed at higher redshifts (Sathyaprakash et al. 2010). Moreover, our method offers a direct measurement of the speed of GWs - the only one proposed besides the time differences between detectors in the network. On the one hand, assuming additional BNS events observed at different redshifts and directions in the sky, the combination of both direct and indirect v_g measurements will place limits on the anisotropy of the speed of gravity. On the other hand, assessing v_g(f_low) from our analysis is of paramount importance for probing the different speed of GWs at high and low frequencies, which are nearly inaccessible at the present stage of GW astrophysics.

Acknowledgments

We thank Prof. Jin Li and Dr. Mengfei Sun for helpful discussion. This work is supported by Beijing Natural Science Foundation No. 1242021, National Natural Science Foundation of China (Nos. 12021003, 12433001); and the Fundamental Research Funds for Central Universities.

References

Appendix A: Theoretical framework of GW speed

In this appendix we introduce the general theoretical framework of frequency-dependent v_g(f) and its corresponding observables that were used in this Letter. This framework has been developed in the LISA science papers regarding tests of modified gravity theories with the future GW space-detectors (Baker et al. 2022). See Baker et al. (2023) for a comprehensive discussion of the speed of GWs in modified gravity theories. It has been shown that a non-trivial speed of GWs, affects both the phase and the amplitude of the signal from a coalescing binary, hence affecting the GW luminosity distance (Belgacem et al. 2019). Furthermore, the expression of the luminosity distance measured with GWs as a function of redshift is also affected by a non-trivial speed of GWs.

Specifically, a quadratic action for the linearised transverse-traceless GW modes was considered (Baker et al. 2022)

$$\begin{array}{c}\hfill S_T=\frac{M_{\mathrm{Pl}}^2}{8}{\displaystyle \int dt}{d}^{3}x{a}^3(t)\overline{\alpha }[{\dot{h}}_{\mathit{ij}}^2-\frac{v_g^2(f)}{a^2(t)}{(\mathbf{\nabla }{h}_{\mathit{ij}})}^2],\end{array}$$(A.1)

where M_Pl denotes the reduced Planck mass. This framework describes free GWs travelling freely in a flat FLRW metric (Qi et al. 2019; Cao et al. 2019)

$$\begin{array}{c}\hfill d{s}^2=v_g(f)\overline{\alpha }[-v_g^2(f)d{t}^2+a^2(t)d{\mathit{x}}^2].\end{array}$$(A.2)

with arbitrary speed of v_g(f). This is an effective metric which LISA Cosmology Working Group use for describing the propagation of the GW (Belgacem et al. 2019), with the Lagrangian density for a free spin-2 field (in the integrand of Eq. (A.1))

$$\begin{array}{cc}\hfill L_T& =\sqrt{-\stackrel{\sim }{g}}\left[{\stackrel{\sim }{g}}^{\mu \nu }{\partial }_{\mu }{h}_{\mathit{ij}}{\partial }_{\nu }{h}_{\mathit{ij}}\right]\hfill \\ \hfill & =a^3\overline{\alpha }[{\dot{h}}_{\mathit{ij}}^2-\frac{v_g^2}{a^2}{\left(\mathbf{\nabla }{h}_{\mathit{ij}}\right)}^2].\hfill \end{array}$$(A.3)

The comoving distance to the GW source whose signal was emitted at time t_e and registered by the observer at time t_o is (Maggiore 2008a)

$$\begin{array}{c}\hfill r_{\mathrm{com}}^{\mathrm{GW}}(t)={\displaystyle {\int }_0^r}d{r}^{\prime }={\displaystyle {\int }_{t_e}^{t_o}}\frac{v_g[f(t^{\prime })]}{a(t^{\prime })}d{t}^{\prime }.\end{array}$$(A.4)

In the discussion below, we use the following notations distinguishing several possible notions of time (and associated quantities, like frequency): t₀ - time measured by clocks of distant observer, t_s time measured by a clock near the source region (local wave zone, intrinsic source time) and t_e cosmic time when the signal was emitted (Baker et al. 2022). Now, the comoving distance defined above, is related to the physical distance to the source as

$$\begin{array}{c}\hfill r_{\mathrm{phys}}^{\mathrm{GW}}(t)=a(t){\left[\frac{v_g(f(t))}{v_g\left(f_s\right)}\right]}^{\frac{1}{2}}{r}_{\mathrm{com}}^{\mathrm{GW}}(t),\end{array}$$(A.5)

where f_s, according to the notation introduced, denotes the intrinsic frequency of the GW emitted. Considering the cosmological redshift effect, the relation between the source and observer frequencies (f_s and f_o) reads

$$\begin{array}{c}\hfill \frac{f_o}{v_g(f_o)}=\frac{f_s}{(1+)v_g(f_s)}.\end{array}$$(A.6)

Based on the luminosity at the source

$$\begin{array}{c}\hfill \mathcal{L}=\frac{d{E}_s}{d{t}_s}=\frac{{(1+z_e)}^2}{{(v_g(f_o)/v_g(f_s))}^2}\frac{d{E}_o}{d{t}_o},\end{array}$$(A.7)

and the energy flux at the observer

$$\begin{array}{c}\hfill \mathcal{F}=\frac{\mathcal{L}}{4\pi {(D_L^{\mathrm{GW}})}^2},\end{array}$$(A.8)

one could obtain the GW luminosity distance as (Maggiore 2008b)

$$\begin{array}{c}\hfill D_L^{\mathrm{GW}}=(1+z)\sqrt{\frac{v_g(f_s)}{v_g(f_o)}}{\displaystyle {\int }_0^z}\frac{v_g[f(z^{\prime })]}{H(z^{\prime })}d{z}^{\prime }.\end{array}$$(A.9)

As already emphasised in the main text, in order to be free from a particular modified gravity theory, the v_g(f) dependence is described by a step function. Therefore, we consider constant value of v_g(f) from the source to the observer v_g(f_s) = v_g(f_o). Then the information about the speed of GW could be derived from either low-frequency or high-frequency bands, with a sharp transition point (∼1Hz) between high-frequency and low-frequency DECIGO range (Baker et al. 2023). Now one has

$$\begin{array}{c}\hfill D_L^{\mathit{GW}}(z)=(1+z){\displaystyle {\int }_0^z}\frac{v_g(f)}{H(z^{\prime })}d{z}^{\prime }\end{array}$$(A.10)

and consequently

$$\begin{array}{c}\hfill \frac{\partial }{\partial z}\frac{D_L^{\mathit{GW}}(z)}{1+z}\equiv \frac{v_g(f)}{H^{\mathit{GW}}(z)}.\end{array}$$(A.11)

The methods and conclusions derived in this work apply beyond the scope of a illustrative example (Salzano 2017).

Appendix B: The derivation of observables with CNN

In this appendix we introduce how to use deep learning to provide a powerful tool for parameter estimation, in the analysis of DECIGO’s time-domain GW data. In the realm of time-domain GW data analysis, CNNs introduce many advantages for parameter estimation. Uniquely, CNNs enable the automated extraction of features, eliminating the need for manual intervention in feature design (Christensen & Meyer 2022). They also exhibit local perception abilities, through which they discern and extract features across different input data positions using convolutional filters, a crucial aspect for identifying local structures and temporal aspects in GW data. This capability is essential for the precise discrimination of various GW signals (Edwards 2021). Further, the multi-layered structure of CNNs, comprising convolutional and pooling layers, allows for a hierarchical representation of features, progressively abstracting the data and thereby augmenting the accuracy of parameter determination. The robustness and adaptability of CNNs, fostered through training on large datasets, play a pivotal role in ensuring resilience against noise and other real-world detection anomalies (George & Huerta 2018). Moreover, their proficiency in processing large-scale datasets is especially advantageous given the expected influx of extensive data from the forthcoming DECIGO detector, which will produce high-temporal-resolution data. This underscores CNNs’ suitability for handling significant data volumes effectively. Thus, the comprehensive benefits offered by CNNs solidify their role as a potent tool for parameter estimation in analysing GW data (Dreissigacker et al. 2019), particularly from DECIGO, leading us to choose a one-dimensional (1D) CNN model for our analysis.

In this analysis, we estimated the uncertainty of H(z) using CNN, with detailed descriptions of its network architecture and hyper-parameters available in Sun et al. (2024). For each BNS generated GW signal, the actual input data for the network was the time-domain data s(t) = h(t) + n(t), where h(t) is the inverse Fourier transform of the frequency-domain GW h(f), and n(t) represents the coloured Gaussian time-domain noise derived from the single-sided noise spectral density of DECIGO (Kawamura et al. 2006). The data length is 1 × 2,000,000, and prior to training the network, we employed a 1D convolutional layer to perform feature extraction and data reshaping on the raw data. This step was necessary because the original signal data length is excessively long, which could lead to extended training periods and difficulties in network convergence (Sun et al. 2024). The entire dataset was divided into a training set (70%) and a test set (30%). The input data were (x_train,y_train) and (x_test,y_test), where x_train and x_test denote the time-domain data for the training and test sets, respectively. y_train and y_test are the parameters that CNN needs as labels in the space θ. The results of the two observables (D_L, X_H) obtained from the neural network are displayed in Fig. A.1.

Fig. A.1. Comparison of DL(z) (left panel) and XH(z) (right panel) error estimates for BNSs observed during one-month and one-year DECIGO observations. The uncertainties of are evaluated with CNN.

Appendix C: Systematical analysis

In this appendix we introduce several sources of systematics based on the one-year observation of DECIGO, which were considered in the above analysis. Firstly, the peculiar acceleration of each binary source could possibly act as an additional noise when measuring the redshift drift. However, the magnitude of peculiar acceleration expected in typical clusters and galaxies is much smaller than the measurement errors of Δ_tz due to detector noise (Amendola et al. 2008; Uzan et al. 2008). Secondly, the performance of DECIGO is evaluated on the Fisher information matrix method for the uncertainty assessment. Such statistics estimates the uncertainty of binary parameters θ_i via $\mathrm{\Delta }{\theta }_i=\sqrt{{({\mathrm{\Gamma }}^{-1})}_{\mathit{ii}}}$ (Cutler & Flanagan 1994). Our findings demonstrated the effectiveness of the machine learning method in accurately inferring the necessary variables and uncertainties. This result is consistent with the one in Sun et al. (2024). Thirdly, comparing the v_g measurements on the ΛCDM cosmology from Planck 2018 results (Planck Collaboration VI 2020), we found that the systematics caused by different fiducial cosmologies are negligible. The fourth issue that needs clarification is the impact of degeneracies between mass and redshift. A recent review of the measurements from GWs of BNS (Landry & Read 2021), pulsars with radio timing, NS with high and low stellar mass companions using X-ray and optical observations (Alsing et al. 2018) showed that the mass distribution of neutron stars is very narrow (m = 1.27 ± 0.04 M_⊙) (You et al. 2025). The impact of such small prior variance on the inverse Fisher matrix for D_L(z) and X_H(z) (and also other parameters) is of the order 𝒪(σ_ℳBNS²). Therefore, introducing the prior σ_mNS ≈ 0.04 M_⊙ in our calculations would not significantly affect the results regarding the final measurement of v_g. Finally, modified gravity theories could possibly affect the estimation of BNS parameters and thus the determination accuracy of v_g. The easiest modification is to add some scalar degrees of freedom to gravity (Brans & Dicke 1961) and one of the first models of this kind is the Brans-Dicke scalar-tensor theory of gravity. Such gravity theory will reduce to GR in the limit of ω_BD → ∞ where ω_BD is the Brans-Dicke parameter (Yagi & Tanaka 2010). Our analysis indicated that the inclusion of such modified gravity theory will result in systematic uncertainties in the extraction of v_g. However, if we imposed the prior information of the source orientation, the constraint would become much stronger.

Finally, we emphasize that our analysis can be extended to the measurements of GW dispersion. Following Abbott et al. (2017c), we consider the following modified dispersion relation of E² = p²c² + Ap^αc^α, α ≥ 0, where E and p denote the energy and momentum of gravitational radiation. A change in the amplitude of dispersion relation A could modify the GW group velocity as v_g/c = 1 + (α − 1)AE^α − 2/2 (Yunes et al. 2016). By introducing the evolution of GW phase in the effective-precession waveform model (Mirshekari et al. 2012), three BBH events (GW170104, GW150914 and GW151226) provided the credible upper bounds of dispersion parameter |A| as 10⁻¹⁹ PeV^2 − α (Abbott et al. 2017c). However, this was achieved using only the high-frequency part of the GW signal. For several modified theories of gravity with various values of α, i.e. doubly special relativity (α = 3; Amelino-Camelia 2002), and extra-dimensional theory (α = 4; Sefiedgar et al. 2011), |A| is expected to be determined at the level of 10⁻⁵ TeV^2 − α, through one-year observation of BNS events. In this aspect, multi-frequency studies involving DECIGO can effectively differentiate between GR and modified gravity theories, which strengthens its probative power to inspire new observing programs in the moderate future.

All Figures

	Fig. 1. Measurements of the speed of GWs based on a one-month (upper panel) and a one-year (lower panel) observation of BNS events by DECIGO. The scale of departure from the GR δvg(flow)/c0 is assumed to be 10−4 and 10−5, respectively.
In the text

Fig. 2. Posterior distributions of vg determined by different techniques. DECIGO data under observing strategies I (one-month) and II (one-year) are compared with two other methods: a direct method based on transit time across a separated network of detectors, and an indirect method using the time delay of GW170817. Two strategies were considered: LIGO O2 – data from O2 of Advanced LIGO-Virgo – and LIGO O2 fixed – where the sky localization of the BNS source is fixed. Note that LIGO O2 and GW170817 only provide constraints on vg/c, and at low redshifts, we take c = c0 in this figure.

In the text

	Fig. A.1. Comparison of DL(z) (left panel) and XH(z) (right panel) error estimates for BNSs observed during one-month and one-year DECIGO observations. The uncertainties of are evaluated with CNN.
In the text