JAXGB Model

The sections below introduce the waveform and response model used in jaxgb.jaxgb.JaxGB.

Strain

Let \(h_{ij}\) denote the (transverse-traceless) strain of a plane gravitational wave (GW), also referred to as the metric perturbation. This strain can be decomposed into its “plus” and “cross” polarisations by means of the polarisation tensors \(e_{ij}^+\) and \(e_{ij}^\times\) associated with the source frame of the GW source:

\[h_{ij}(\xi) = h_+(\xi) e_{ij}^+ + h_\times(\xi) e_{ij}^\times,\]

where

\[\xi = t - \frac{\hat{\mathbf k} \cdot \mathbf{x}(t)}{c}\]

expresses the surfaces of constant phase. In the above expression, the (unit) vector \(\hat{\mathbf k}\) is the GW direction of propagation, while \(\mathbf{x}\) is the position and \(t\) the time.

More specifically, for a mildly-chirping, circular Galactic binary, the polarisations can be expressed as

\[ \begin{align}\begin{aligned}h_+(\xi) &= -\mathcal A (1 + \cos^2\iota) \cos\Psi(\xi),\\h_\times(\xi) &= -2\mathcal A\,\cos\iota \sin\Psi(\xi),\end{aligned}\end{align} \]

with

\[\Psi(\xi) = -\phi_0 + 2\pi f_0 \xi + \pi \dot{f}_0 \xi^2.\]

The GW amplitude \(\mathcal A\), inclination \(\iota\), initial phase \(\phi_0\), GW frequency \(f_0\), and GW frequency first derivative \(\dot f_0\) follow the usual parameterisation.

Aside from the propagation vector \(\hat{\mathbf k}\), the unit vectors \(\hat{\mathbf u}\) and \(\hat{\mathbf v}\) can be defined to form the orthonormal triad

\[ \begin{align}\begin{aligned}\hat{\mathbf u} &= (\sin\alpha, -\cos\alpha, 0),\\\hat{\mathbf v} &= (-\sin\delta\cos\alpha, -\sin\delta\sin\alpha, \cos\delta),\\\\hat{\mathbf k} &= (-\cos\delta\cos\alpha, -\cos\delta\sin\alpha, -\sin\delta).\end{aligned}\end{align} \]

Here, the right ascension and declination (\(\alpha\), \(\delta\)) express the sky position of the GW source. Defining then the GW polarisation angle \(\psi\), the polarisation tensors \(e^+_{ij}\) and \(e^\times_{ij}\) can be expanded as

\[ \begin{align}\begin{aligned}e_{ij}^+ &= \cos 2 \psi \epsilon_{ij}^+ + \sin 2 \psi \epsilon_{ij}^\times,\\e_{ij}^\times &= -\sin 2 \psi \epsilon_{ij}^+ + \cos 2 \psi \epsilon_{ij}^\times,\end{aligned}\end{align} \]

where

\[ \begin{align}\begin{aligned}\epsilon_{ij}^+ &= (\hat{\mathbf u}\otimes\hat{\mathbf u} - \hat{\mathbf v}\otimes\hat{\mathbf v})_{ij},\\\epsilon_{ij}^\times &= (\hat{\mathbf u}\otimes\hat{\mathbf v} + \hat{\mathbf v}\otimes\hat{\mathbf u})_{ij}.\end{aligned}\end{align} \]

These polarisation tensors, associated with the Solar System barycenter (SSB), enable another expansion of the GW strain:

\[h_{ij}(\xi) = h_+^\mathrm{SSB}(\xi)\epsilon_{ij}^+ + h_\times^\mathrm{SSB}(\xi)\epsilon_{ij}^\times,\]

where

\[ \begin{align}\begin{aligned}h_+^\mathrm{SSB} &= h_+ \cos 2 \psi - h_\times \sin 2 \psi,\\h_\times^\mathrm{SSB} &= h_+ \sin 2 \psi + h_\times \cos 2 \psi.\end{aligned}\end{align} \]

To end this section on the strain definition, it can be mentioned that the 8 above GW parameters (\(\mathcal A\), \(\iota\), \(\phi_0\), \(f_0\), \(\dot f_0\), \(\psi\), \(\alpha\), and \(\delta\)) uniquely define any given circular Galactic binary.

Response

To calculate the response to the above strain, the single-link response function can be introduced. Let the indices \(i\) and \(j\) denote the receiving and the emitting spacecraft associated with a given laser link. The relative frequency fluctuations caused by a GW can be obtained as follows:

\[y_{ij}(t_i) = -\frac{1}{2} \frac{H_{ij}(\xi_i) - H_{ij}(\xi_j)}{1 - \hat{\mathbf{k}} \cdot \hat{\mathbf{x}}_{ij}},\]

where

\[H_{ij}(\xi) = \hat{\mathbf x}_{ij}^k \hat{\mathbf x}_{ij}^l h_{kl}(\xi) = \left((\hat{\mathbf u}\cdot\hat{\mathbf x}_{ij})^2 - (\hat{\mathbf v}\cdot\hat{\mathbf x}_{ij})^2\right) h_+^\mathrm{SSB}(\xi) + 2(\hat{\mathbf u}\cdot\hat{\mathbf x}_{ij})(\hat{\mathbf v}\cdot\hat{\mathbf x}_{ij}) h_\times^\mathrm{SSB}(\xi).\]

In the above expression, the vector

\[\hat{\mathbf x}_{ij} = \frac{\mathbf x_i(t_i) - \mathbf x_j(t_j)}{\lVert \mathbf x_i(t_i) - \mathbf x_j(t_j) \rVert}\]

is the unit vector pointing from the emitting spacecraft to the receiving spacecraft.

Next, a more specific expression to calculate the response for a circular Galactic binary strain can be obtained. It is convenient to express the strain in the complex domain:

\[h_{ij}(\xi) = \mathfrak{R}\left\{\left(A_+^\mathrm{SSB}\epsilon_{ij}^+ + A_\times^\mathrm{SSB}\epsilon_{ij}^\times\right)e^{i\Psi(\xi)}\right\},\]

where

\[A_+^\mathrm{SSB} = -\mathcal A\left((1+\cos^2\iota)\cos 2\psi + 2i\cos\iota\sin 2\psi\right),\]

\[A_\times^\mathrm{SSB} = -\mathcal A\left((1+\cos^2\iota)\sin 2\psi - 2i\cos\iota\cos 2\psi\right).\]

In a next step, define

\[\mathcal A_{ij} = A_+^\mathrm{SSB}\left((\hat{\mathbf u}\cdot\hat{\mathbf x}_{ij})^2 - (\hat{\mathbf v}\cdot\hat{\mathbf x}_{ij})^2\right) + 2A_\times^\mathrm{SSB}(\hat{\mathbf u}\cdot\hat{\mathbf x}_{ij})(\hat{\mathbf v}\cdot\hat{\mathbf x}_{ij}),\]

so that

\[y_{ij}(t_i) = -\frac{1}{2}\,\mathfrak R\left\{\frac{\mathcal A_{ij}\left(e^{i\Psi(\xi_i)} - e^{i\Psi(\xi_j)}\right)}{1 - \hat{\mathbf{k}}\cdot \hat{\mathbf{x}}_{ij}}\right\}.\]

The following step is to rewrite the difference between the two exponentials in the above equation:

\[e^{i\Psi(\xi_i)} - e^{i\Psi(\xi_j)} = e^{i\Psi(\xi_i)}\, e^{-\frac{i}{2}(\Psi(\xi_i) - \Psi(\xi_j))} \left( e^{\frac{i}{2}(\Psi(\xi_i) - \Psi(\xi_j))} - e^{-\frac{i}{2}(\Psi(\xi_i) - \Psi(\xi_j))} \right)\]

\[= 2i\, e^{i\Psi(\xi_i)}\, e^{-\frac{i}{2}(\Psi(\xi_i) - \Psi(\xi_j))} \sin\left(\frac{1}{2}(\Psi(\xi_i) - \Psi(\xi_j))\right).\]

The difference between the phases \(\Psi\) can be evaluated further:

\[ \begin{align}\begin{aligned}\frac{1}{2}\big(\Psi(\xi_i) - \Psi(\xi_j)\big) &= \frac{1}{2}\big( -\phi_0 + 2\pi f_0 \xi_i + \pi \dot{f}_0 \xi_i^2 - (-\phi_0 + 2\pi f_0 \xi_j + \pi \dot{f}_0 \xi_j^2) \big)\\&= \pi (\xi_i - \xi_j)\left(f_0 + \dot f_0 \frac{1}{2}(\xi_i + \xi_j)\right).\end{aligned}\end{align} \]

where the definition of \(\xi\) can be used to write

\[\xi_i - \xi_j = t_i - \frac{\hat{\mathbf k}\cdot \mathbf{x}_i(t_i)}{c} - \left( t_j - \frac{\hat{\mathbf k}\cdot \mathbf{x}_j(t_j)}{c} \right) = L_{ij}(t_i) - \frac{\hat{\mathbf k}\cdot (\mathbf{x}_i(t_i) - \mathbf{x}_j(t_j))}{c}.\]

In the above expression, \(L_{ij}(t_i) = t_i - t_j\) is the light-travel time from spacecraft \(j\) to \(i\), which is the difference between the time of reception \(t_i\) and the time of emission \(t_j\).

To simplify the notation, the following definition is introduced:

\[\alpha_{ij}(t_i) = \pi L_{ij}(t_i) \left(1 - \hat{\mathbf k}\cdot \hat{\mathbf{x}}_{ij}\right) (f_0 + \dot f_0 \xi_i)\]

(note the dependency of the “hat” variable \(\hat{\mathbf x}_{ij}\), introduced earlier on).

Making the approximation that

\[L_{ij}(t_i) \simeq \frac{\lVert \mathbf x_i(t_i) - \mathbf x_j(t_j) \rVert}{c}\]

and that

\[\dot f_0(\xi_i - \xi_j)/2 \ll f_0\]

(this is the “adiabitic approximation”, assuming no significant change of the gravitational-wave strain frequency during the light travel time), it follows that

\[\frac{1}{2}\big(\Psi(\xi_i) - \Psi(\xi_j)\big) \simeq \alpha_{ij}(t_i).\]

With the above definition and assumption, the response can be rewritten (or rather, approximated) as

\[y_{ij}(t_i) = -\mathfrak R\left\{ \frac{i \mathcal A_{ij} \sin \alpha_{ij}} {1 - \hat{\mathbf{k}}\cdot \hat{\mathbf{x}}_{ij}} e^{-i\big(\alpha_{ij} - (-\phi_0 + 2\pi f_0 \xi_i + \pi \dot f_0 \xi_i^2)\big)} \right\}.\]

Heterodyned response

For typical data generation scenarios, the response is evaluated on a discrete grid for a finite observation time \(T_\mathrm{obs}\). In the corresponding frequency-domain representation of this signal, frequencies are spaced by integer numbers of \(T_\mathrm{obs}^{-1}\). The (integer) index of the closest frequency bin \(q\) to a given GW frequency \(f_0\) can be calculated by the following rounding-to-closest-integer operation:

\[q = \mathrm{round}[f_0 T_\mathrm{obs}], \qquad q \in \mathbb{N}.\]

Using \(q\), the response can be rewritten as

\[y_{ij}(t_i) = \mathfrak R\left\{\mathcal G_{ij} e^{2\pi i q t_i / T_\mathrm{obs}}\right\},\]

where

\[\mathcal G_{ij} = -\mathfrak R\left\{ \frac{i}{2} \frac{f_0 + \dot f_0 \xi_i}{f_{\star,ij}(t_i)} \mathcal A_{ij}\,\mathrm{sinc}\,\alpha_{ij} e^{-i\big( \alpha_{ij} + 2\pi f_0 \frac{\hat{\mathbf k}\cdot \mathbf{x}_i(t_i)}{c} - (-\phi_0 + 2\pi (f_0 - q/T_\mathrm{obs}) t_i + \pi \dot f_0 \xi_i^2) \big)} \right\}.\]

The above expression introduces \(f_{\star,ij} = (2 \pi L_{ij}(t_i))^{-1}\).

Under some conditions, \(\mathcal G_{ij}\) is a slowly varying function over time, allowing the single-link response to be decomposed into a fast component (\(e^{2\pi i q t_i / T_\mathrm{obs}}\)) and a slowly varying component (\(\mathcal G_{ij}\)). In general, these conditions are: (1) observation times \(T_\mathrm{obs}\) need to be sufficiently long, such that \((f_0 - q/T_\mathrm{obs}) \ll q/T_\mathrm{obs}\), (2) \(\dot f_0 \xi_i \ll f_0\), i.e., the first-order change in strain frequency needs to stay well below the zeroth order gravitational-wave strain frequency, during the relevant time interval.

The Fourier transforms of \(y\) and \(\mathcal G\) are now related as follows:

\[\widetilde{y}_{ij}(f) = \widetilde{\mathcal G}_{ij}(f - q / T_\mathrm{obs}).\]

In words, the frequency-domain response can be obtained by computing \(\mathcal G_{ij}\), transforming it into the frequency domain, and shifting it to higher frequencies. Given that \(\mathcal G_{ij}\) is slowly varying over time, its frequency content only covers a narrow part of the low-frequency domain.

In principle, to obtain \(y_{ij}\), it is sufficient to have a sparse time-domain evaluation of \(\mathcal G_{ij}\), transform to the Fourier domain, shift the signal as in the above equation, and, optionally, transform back to time domain if needed.

Windowing

For specific purposes, it is useful to window (i.e., taper) the function \(\mathcal G_{ij}\) over the (finite) time domain over which it is calculated. For example, analysing a finite stretch of data usually requires a similar window function to be applied on the model response.

Defining a generic windowing function \(W(t)\), the heterodyned approach is updated as follows:

\[\mathcal G_{ij}(t_i) \rightarrow W(t_i)\,\mathcal G_{ij}(t_i).\]

Examples of window functions supported in JAXGB include:

Boxcar/rectangular window: this window is 1 over the full time domain over which the slow response is calculated, i.e., \(W(t_i) = 1\).
`Planck-taper window <https://iopscience.iop.org/article/10.1088/0264-9381/27/8/084020>`_: with an additional parameter \(\alpha\) setting the amount of tapering (smaller values giving sharper transitions).
“Tukey4”: a classical Tukey (cosine-tapered) window raised to the 4th power, as defined in Riegger et al. (2024). This window function also has an additional parameter \(\alpha\) defining the amount of tapering (\(\alpha = 0\) is equivalent to boxcar, values closer to 1 have smoother edges).

Trade-off between performance and accuracy

JAXGB introduces two implementation variants:

jaxgb.jaxgb.JaxGB: baseline implementation focused on performance, using the above model and approximations. This implementation is intended for fast data analysis.
jaxgb.jaxgb.JaxGBaccurate: this mostly relies on the jaxgb.jaxgb.JaxGB implementation but (1) avoids most of the above simplifications, (2) supports the extension of the simulated time domain to reduce edge artefacts and obtain a response that is accurrate over the full intended time domain. This implementation is intended for accurate time-domain injection.

JaxGB

The baseline implementation makes the following simplifications:

The emitting spacecraft position \(\mathbf{x}_j\) at time of emission \(t_j\) can be approximated as the position at time of reception \(t_i\): \(\mathbf{x}_j(t_j) \simeq \mathbf{x}_j(t_i)\).
Light-travel times are equal and constant; they are approximated as the average LTT value across the 6 links and over the full time domain of interest.

A consequence of the first point is that reversing the direction of transfer on a single link only brings a sign flip (and no other direction or magnitude change):

\[\mathbf x_{21} = -\mathbf x_{12}, \quad \mathbf x_{32} = -\mathbf x_{23}, \quad \mathbf x_{13} = -\mathbf x_{31},\]

which simplifies the calculation of any dependent variable (e.g., \(\mathcal A_{ij} = \mathcal A_{ji}\)). The second point reduces \(f_{\star, ij}\) to a constant (scalar) quantity. This lowers the computational effort associated with the response calculation while still preserving (typically percentage-level) accuracy.

JaxGBaccurate

The more accurate implementation drops the above assumptions. It evaluates emitting spacecraft positions at times of emission \(\mathbf{x}_j(t_j)\) and full link vectors for which, in general, \(\mathbf{x}_{ij}(t_i) \neq \mathbf{x}_{ji}(t_j)\). Also, \(f_{\star, ij}\) values are calculated using actual, time-dependent light-travel times \(L_{ij}(t_i)\).

To avoid edge artefacts, the windowing approach introduced above is exploited. Following Riegger et al. (2024), the time domain of the slow response \(\mathcal G_{ij}\) calculation is extended from \([0, T_\mathrm{obs}]\) to \([-\varepsilon T_\mathrm{obs}, (1+\varepsilon)T_\mathrm{obs}]\), where \(\varepsilon \ge 0\) is an extension factor.

The “Tukey4” window function \(W_a(t)\) used for tapering the time-domain slow response is defined as:

\[\begin{split}W_a(t) = \begin{cases} 0, & t < -\varepsilon T_\mathrm{obs}, \\ \frac{1}{16}\left(1 - \cos\left(\frac{\pi(t + \varepsilon T_\mathrm{obs}(1-a))}{2a\varepsilon T_\mathrm{obs}}\right)\right)^4, & -\varepsilon T_\mathrm{obs} \le t < -\varepsilon T_\mathrm{obs}(1-a), \\ 1, & -\varepsilon T_\mathrm{obs}(1-a) \le t < T_\mathrm{obs} + \varepsilon T_\mathrm{obs}(1-a), \\ \frac{1}{16}\left(1 - \cos\left(\frac{\pi(T_\mathrm{obs} + \varepsilon T_\mathrm{obs}(1-a) - t)}{2a\varepsilon T_\mathrm{obs}}\right)\right)^4, & T_\mathrm{obs} + \varepsilon T_\mathrm{obs}(1-a) \le t < (1+\varepsilon)T_\mathrm{obs}, \\ 0, & (1+\varepsilon)T_\mathrm{obs} \le t. \end{cases}\end{split}\]

with \(0 < a < 1\).

Ultimately, this leads to the following expression for the discretised, time-domain single-link response \(y_{ij}\):

\[y_{ij}(t_k) = e^{-2\pi i \varepsilon q T_\mathrm{obs}/T_\mathrm{ext}} \,\mathrm{IFFT}\left[ \mathrm{FFT}[\mathcal G_{ij} W_a] \left(f - \frac{q}{T_\mathrm{ext}}\right) \right][k],\]

where

\[t_k = -\varepsilon T_\mathrm{obs} + k\,\Delta t, \qquad k = 0, \ldots, N_\mathrm{ext} - 1.\]

In the above:

\(\Delta t\) is the desired time step of the response
\(N_\mathrm{ext} = T_\mathrm{ext} / \Delta t\) is the size of the extended time-domain response array
\(T_\mathrm{ext} = (1 + 2\varepsilon) T_\mathrm{obs}\) is the duration of the extended time domain. The right-hand-side exponential factor depending on the ratio \(T_\mathrm{obs}/T_\mathrm{ext}\) corrects the response phase related to the time-domain extension prior to \(t=0\).

In words:

(A sparse evaluation of) the slow response is calculated and windowed using the “Tukey4” window; the number of time samples is referred to as \(n\).
This response is Fourier transformed and modulated to a higher frequency (cf. the \(q/T_\mathrm{ext}\) shift in the frequency domain), further upsampled by zero padding up to the frequency \(1/\Delta t\) (i.e., \(N_\mathrm{ext}\) samples with frequency step \(1/T_\mathrm{ext}\)).
The result is transformed back to the time domain, with an additional multiplicative phase factor \(e^{-2\pi i \varepsilon q T_\mathrm{obs}/T_\mathrm{ext}}\) related to the time-domain extension, to yield the single-link response at the desired \(\Delta t\).
Any samples in the above time series beyond the time domain \([0, T_\mathrm{obs}]\) can be dropped.

Practical implementation for Mojito

The above formulation requires the specification of three parameter values:

For \(a\), associated with the window function, the value \(a = 0.45\) is picked (in line with Riegger et al.).
For \(\varepsilon\), which is the extension factor of the time domain, the value \(1/2\) is picked. Some care is needed to choose this value to ensure the original time grid points are included in the extended time grid. Picking \(\varepsilon = 0.5\) makes \(T_\mathrm{ext} = (1 + 2\varepsilon) T_\mathrm{obs}\) an integer multiple of \(T_\mathrm{obs}\), while only extending the time domain by a factor 2.
For \(n\), the number of slow-response points to be evaluated in the time domain, the value \(n = 4096\) is chosen. This value is sufficient to capture the evolution of the slow response for all \((f, \dot{f})\) pairs in the Mojito population.