[15] Thomas Middelmann et al. "High Accuracy Correction of Blackbody Radiation Shift in an Optical Lattice Clock". In: *Phys. Rev. Lett.* 109 (26 Dec. 2012), p. 263004. DOI: 10.1103/PhysRevLett.109.263004. URL: https://link.aps.org/doi/10.1103/PhysRevLett.109.263004.

[16] T F Gallagher. "Rydberg atoms". In: *Reports on Progress in Physics* 51.2 (Feb. 1988), p. 143. DOI: 10.1088/0034-4885/51/2/001. URL: https://dx.doi.org/10.1088/0034-4885/51/2/001.

[17] J Millen et al. "Spectroscopy of a cold strontium Rydberg gas". In: *Journal of Physics B: Atomic, Molecular and Optical Physics* 44.18 (Sept. 2011), p. 184001. DOI: 10.1088/0953-4075/44/18/184001. URL: https://dx.doi.org/10.1088/0953-4075/44/18/184001.

[18] Elizabeth M. Bridge et al. "Tunable cw UV laser with &\#x0003C;35 kHz absolute frequency instability for precision spectroscopy of Sr Rydberg states". In: *Opt. Express* 24.3 (Feb. 2016), pp. 2281-2292. DOI: 10.1364/OE.24.002281. URL: https://opg.optica.org/oe/abstract.cfm?URI=oe-24-3-2281.
# Part IV

## Dish antennas detectors
# Chapter 11

## Search for dark photons using dish antennas

As a reminder, DP will induce an electric field everywhere in space due to its kinetic mixing with EM. In Chapter 10, we focused on an experiment to enhance this electric field with a cavity. In this chapter, we will focus on a different way to enhance this electric field, using a spherical dish antenna. Indeed, one can use a dish antenna to focus the small electric field in order to enhance it inside a detector located at the curvature center of the mirror [1–9]. One of the main assumptions used in existing studies is that the power generated by the dish is entirely focused on the curvature center, hence can be detectable by the detector [10].

In this chapter, we will show that this assumption is too optimistic and we will derive an improved modeling for such experiments. More precisely, we will derive analytically the EM field measured by the antenna detector which is highly dependent on the DP Compton frequency. This is done in two steps: (i) we compute the electric field at the location of the antenna and (ii) we determine how much of this electric power is measured by the antenna. This study is the object of an article [11] published in *Physical Review D*.

### 11.1 Electromagnetic emission from the dish

For our analysis, we consider a dish of radius $r$ and curvature radius $R$, as depicted in Fig. 11.1, which works as an EM reflector of the oscillating electric field Eq. (8.63c). Indeed, boundary conditions $\vec{E}_{tot,\parallel} = 0$ at the edge of the dish requires the emission of a regular EM wave, at every point of the dish, directed according to the normal vector of the dish, and with same frequency as the background DM oscillating field, i.e it has the form

$$ \vec{E}_D(\theta, \varphi) = i\chi\omega_U\vec{Y}_{\parallel,D}(\theta, \varphi)e^{-i\omega_U t}, \qquad (11.1) $$

with $\vec{Y}_{\parallel,D}(\theta, \varphi)$ the component of the DP polarization vector $\vec{Y}$ defined in Eq. (8.61a), parallel to the radial vector at cylindrical coordinates $(\theta, \varphi)$ such that the electric field parallel to the dish vanishes at the mirror’s surface. In addition, we choose a coordinate system whose origin is at the curvature center of the dish.

In the literature (see e.g [10]), one of the main assumptions used in this kind of experiment is that provided that the diameter of the dish $d_{dish} \gg \lambda$, the wavelength of the emitted field (i.e the Compton wavelength of the DP field), diffraction effects are negligible and the power
<figure><img src="image_9.png" /><figcaption>Figure 11.1: The oscillating standing DP electric field $\vec{E}_{DM}$ is shown with red arrows. The dish emitter in black, with curvature center shown in purple, acts as reflector and emits a classic propagating electric field $\vec{E}_D$ in green. The Kirchhoff integral boundary surfaces A considered for the computation of the electric field Eq. (11.6) are shown in blue (in addition to the dish, surface $A_1$). It consists of the junction of a semi-sphere centered on the curvature center of the dish with infinite radius (surface $A_3$), and one additional vertical infinite plane surfaces to close the boundary surface (surface $A_2$). This closed surface is chosen such that the contributions from surface $A_2$ and $A_3$ are zero at the curvature center of the dish.</figcaption></figure>

generated by the dish is entirely focused on the curvature center of the mirror [10]. If we further assume that the field is homogeneous over the surface of the spherical mirror $A_{\text{dish}}$ (i.e the DM field coherence length $l_{\text{DM}} = \hbar/m_{\text{DM}}v_{\text{DM}} \gg d_{\text{dish}}$), this total power emitted is [10]

$$P_{\text{emit}} = A_{\text{dish}} \langle |E_D(\theta, \phi)|^2 \rangle = \langle \beta^2 \rangle A_{\text{dish}} \chi^2 \rho_{\text{DM}} c \quad (11.2)$$

where $\beta$ is the projection of the DP field on the surface of the dish, similarly defined as in Eq. (10.12b) and $\langle \rangle$ represents the average value over the surface. Assuming random polarization scenario, $\beta = 2/3$, as we consider a planar emitter.

While all previous analysis were based on these assumptions, we will show in this chapter that they are over-optimistic. Indeed, we will provide a more detailed modeling of such experiments by deriving properly the propagation of this emitted field from the dish using Kirchhoff integral. We will show that the amount of power received at the curvature center does depend on a various number of parameters, mainly the Compton frequency of the DP field, and that in some cases, reaches only a few percent of the emitted power, even when the radius of the dish is larger than the wavelength by more than one order of magnitude, i.e when diffraction effects appear negligible.
## 11.2 Propagation of the field from the dish to an antenna

In this section, we will derive the expression of the total electric field induced by the reflection of the DP field by the dish at a given location. In Sec. 11.2.1, we will give a brief overview of the Kirchhoff method and show that it cannot be directly used to get an analytical expression of the electric field. Instead, we will decompose the problem into two sub-problems. First, we will propagate the field from the dish to the plane closing the dish (displayed in orange in Fig. 11.2) using an approximation of the Kirchhoff method valid for thin optical element as presented in Sec. 11.2.2. Subsequently, in Sec. 11.2.3, we will propagate the electric field from the plane to the position of the detector using exactly the Kirchhoff method.

### 11.2.1 Kirchhoff integral

The Kirchhoff theorem allows one to solve the wave equation at a given position $\vec{x}$ by computing an integral over a closed surface around $\vec{x}$ [12]. In practice, this means that one can derive the electric field at position $\vec{x}$ as function of the field at all points on this closed surface.

In this formalism, the temporal dependence of the field is separated from its spatial dependence. Therefore, we decompose the emitted electric field at the dish's surface as

$$ \vec{E}_D(\vec{x}', t) = \Re[\vec{U}_D(\vec{x}')e^{-i\omega t}], \qquad (11.3) $$

where the complex function $\vec{U}_D(\vec{x}')$ denotes the spatial part of the field and $\vec{x}'$ is a point on the dish's surface.

We now consider a closed surface, denoted A, which encloses the point $\vec{x}$ where the value of the field is calculated. From Kirchhoff integral theorem, the general expression of $\vec{U}_D$ is [12],

$$ \vec{U}_D(\vec{x}) = \int_A dS' \left( \vec{U}_D(\vec{x}')(\hat{n}' \cdot \nabla G(\vec{x}, \vec{x}')) - G(\vec{x}, \vec{x}')(\hat{n}' \cdot \nabla)\vec{U}_D(\vec{x}') \right), \qquad (11.4) $$

where $\hat{n}'$ is a unit vector normal to the surface A directed inwardly, the derivatives are with respect to the emission coordinates $x'$ and G is a Green function, appropriately defined for the situation.

In the situation depicted in the previous section, from Eq. (11.1), we know exactly the components of the electric field at the surface of the emitting dish. Therefore, an appropriate Green function is the Dirichlet Green function $G_D$ [12]

$$ G_D(\vec{x}, \vec{x}') = 0 \quad \forall \vec{x}' \in A. \qquad (11.5) $$

Then, Eq. (11.4) becomes

$$ \vec{U}_D(\vec{x}) = \int_A dS \left( \vec{U}_D(\vec{x}')(\nabla G_D(\vec{x}, \vec{x}') \cdot \hat{n}') \right). \qquad (11.6) $$

In the following of the chapter, unless otherwise specified, whenever we mention electric fields, we refer to $\vec{U}_D$, i.e its spatial part, following Eq. (11.6).
<figure><img src="image_10.png" /></figure>

Figure 11.2: To compute the electric field induced by the dish at a location *P* at cylindrical coordinates (*ρ*, *φ*, *z*), we first propagate the electric field emitted by the spherical dish (in light grey) into the fictional plane (in orange) following Eq. (11.7). Then, we propagate this electric field to the receiver using the Kirchhoff integral. This procedure is valid only for a spherical dish with low curvature, i.e whose radius *r* is much smaller than its curvature radius *R*.

In the present situation, the closed surface A is appropriately defined by the junction of the dish, a surface with radius $R \to \infty$ and one additional plane surface along $(x, y)$ plane to close the surface, such that only the field emitted by the dish will contribute. The Kirchhoff integral surface boundary is shown in blue in Fig. 11.1.

Our system is therefore composed of a spherical dish and Eq. (11.6) requires the derivation of a Dirichlet Green function for a portion of sphere. However, no analytical expression has been found for such geometry, which makes it impossible to find an exact analytical solution for our problem. Instead, we will decompose the problem into two sub-problems: using an approximation, we will propagate the electric field from the dish to the plane closing the dish (see Fig. 11.2) and then, using the Kirchhoff integral, from the plane up to the detector.

### 11.2.2 Propagation of the electric field from the dish to the plane

This situation is depicted in Fig. 11.2, where the distance between the fictional plane (orange on the figure) and the dish at coordinates $x = y = 0$ is noted *a*.

As mentioned above, it is not possible to find an analytical solution to Eq. (11.6) for this particular geometry. Instead, we use the thin optical element approximation presented in details in Sec. 2.2.7 from [13] which allows one to find an analytical expression for the field $\vec{U}(\vec{x})$ of Eq. (11.6) on this fictional plane. This approximation is valid if the two following conditions are fulfilled:

* The transverse propagation modes $p, q$ (which can be viewed as the coordinates of the wavevector in the $x-y$ plane [13]) need to be much smaller than the longitudinal one $k$. In our case, the transverse modes are determined by the spatial extent of the dish antenna in the $x-y$ plane $\Delta x, \Delta y$, with $p \sim q \sim 1/\Delta x \sim 1 \text{ m}^{-1} \ll \omega_U/c = k \sim 10^3 \text{ m}^{-1}$, in the GHz range. In addition, the galactic velocity distribution in the DM halo induces
a transverse contribution [14], of the order of $\omega_U v_{DM}/c^2 \sim 3 \times 10^{-3} \text{ m}^{-1} \ll \omega_U/c \sim k$. Therefore, this condition is fulfilled.

* The dish needs to have a low curvature, or in other words, the radius is much smaller than the curvature radius $r \ll R$. Therefore, in the following, we will restrict ourselves to low curvature dish emitter.

A straightforward but lengthy calculation presented in details in [13] shows that, under the conditions detailed above, the Kirchhoff integral Eq. (11.6) simply reduces to

$$ \vec{U}_P(x, y, z_{\text{plane}}) = -e^{ikf(x,y)} \vec{U}_D(x, y, z_{\text{dish}}(x, y)), \quad (11.7) $$

where $\vec{U}_D$ is the electric field emitted by the dish and $\vec{U}_P$ is the electric field at the location of the closing plane. Furthermore, the function $f(x, y) = z_{\text{dish}}(x, y) - z_{\text{plane}}$ is the surface equation of the dish. In other words, within the approximation of thin optical element, the curvature of the dish translates into simple phase factors determined by the distance between the closing plane and the dish [13].

We now introduce a cylindrical coordinate system $(\rho, \phi, z)$, see Fig. 11.2. Because of the cylindrical symmetry, the surface equation that writes

$$ f(x, y) = f(\rho) = \sqrt{R^2 - \rho^2} - R + a \approx \frac{r^2 - \rho^2}{2R}, \quad (11.8) $$

depends only on the cylindrical coordinate $\rho$, where we used $a \approx r^2/2R$, since the dish is a portion of sphere, with small curvature.

Inserting the spatial part of the field on the dish from Eq. (11.1) into Eq. (11.7), we can now estimate the expression of the field in this fictional plane as

$$ \vec{U}_P(\rho, \phi, z_{\text{plane}}) = -i\chi\omega_U e^{ikf(\rho)} \vec{Y}_{\parallel,D}(\rho, \phi, z_{\text{plane}} + f(\rho)), \quad (11.9) $$

where $\rho \le r$ and $z_{\text{plane}} = R - a$ is the z-coordinate of the plane. This expression is only valid in the thin optical element approximation, whose conditions are detailed at the beginning of this section.

### 11.2.3 Propagation of the electric field from the plane to the receiver

We can now focus on the propagation from an emitting plane, a subject that has been vastly treated in the literature, e.g in [12, 13]. The idea is to use the Kirchhoff integral from Eq. (11.4) using an appropriate Dirichlet Green function that vanishes on the plane. In the case of a plane emitter located at $z_{\text{plane}} = R - a$, such a Dirichlet Green function is given by

$$ G_D(\vec{x}, \vec{x}') = \frac{e^{ikL'}}{4\pi L'} - \frac{e^{ikL''}}{4\pi L''}, \quad (11.10) $$

where [12, 13]

$$ L' = |\vec{x} - \vec{x}'| = \sqrt{\rho^2 + \rho'^2 - 2\rho\rho' \cos(\phi - \phi') + (R - a + \Delta z - z')^2} \quad (11.11a) $$

$$ L'' = |\vec{x} - \vec{x}''| = \sqrt{\rho^2 + \rho'^2 - 2\rho\rho' \cos(\phi - \phi') + (R - a - \Delta z - z')^2}, \quad (11.11b) $$
with $\vec{x}''$ the point symmetrical to $\vec{x}'$ with respect to the plane. In Eq. (11.11), we defined

$$ \Delta z = z - (R - a), \qquad (11.12) $$

the difference of z coordinates between the receiving point $z < z'$ and the emitter on the plane. One can easily show that this property leads to $L' = L''$ if $x'$ is on the plane and therefore the required condition $G_D(\vec{x}, \vec{x}_{\text{plane}}) = 0$ is satisfied. To compute Eq. (11.6), we use the directional derivative since $\hat{n}' = \hat{z}'$

$$ \vec{\nabla} G_D(\vec{x}, \vec{x}') \cdot \hat{n}' = \frac{\partial G_D}{\partial z'} \qquad (11.13a) $$

$$ = - \frac{(ikL' - 1)(R - a + \Delta z - z')e^{ikL'}}{4\pi L'^3} + \frac{(ikL'' - 1)(R - a - \Delta z - z')e^{ikL''}}{4\pi L''^3}. \qquad (11.13b) $$

Now, plugging $z' = z_{\text{plane}}$ and considering reception at $\vec{x} = (\rho, \phi, z)$, we get

$$ \left. \frac{\partial G_D}{\partial z'} \right|_{z' \in \text{plane}} = - \frac{(ikL - 1)\Delta z}{2\pi L^3} e^{ikL}, \qquad (11.14) $$

$$ \text{with } L = \sqrt{\rho^2 + \rho'^2 - 2\rho\rho' \cos(\phi - \phi') + (\Delta z)^2}. $$

Then, the electric field on the fictional plane propagating towards $z < R - a$ reads

$$ \vec{U}_{\text{dish}}(\rho, \phi, z) \approx -\frac{i\chi\omega_U\Delta z}{2\pi} \int_0^r d\rho' \rho' e^{ikf(\rho')} \int_0^{2\pi} d\phi' \frac{ikL - 1}{L^3} e^{ikL} \vec{Y}_{\parallel,D}(\rho', \phi', f(\rho') + R - a), \qquad (11.15) $$

This result provides the expression of the electric field induced by the DP field reflected by the dish anywhere in space under the approximation that the dish is a thin optical element. In particular, this formula includes diffraction effects that were implicitly neglected in previous studies [10]. This integral is generally not solvable analytically and some approximations might be made to simplify it.

## 11.3 Detection of the electric field with a horn antenna

While the previous section was devoted to the emission and propagation of the electric field induced by the dark photon and enhanced by the dish, we will now focus on the detection of this electric field. In this section, we consider the detection system to be a horn antenna of long side A and small side B. At the output of the antenna is located a resistance $R_0$ such that when an oscillating electric field is applied to the horn antenna, it is translated into a measurable voltage.

We will compare two different approaches to estimate this measurement: the first one is based on an analytical calculation which computes the overlap integral between the electric field at the location of the antenna and the antenna mode; the second one consists in using the antenna gain factor provided by the antenna's manufacturer.
### 11.3.1 Computation using the modes overlap

Let us first consider an antenna of internal resistance $R_0$ as an emitter by applying an oscillating voltage $V(t) = V_0 \cos(\omega_U t)$ to its terminals. The reciprocity theorem ensures that this will be equivalent to considering the antenna as a receiver. As a result, the antenna will emit an electric field predominantly in the TE$_{10}$ mode, i.e. polarized parallel to the small side of the rectangular horn, which is the most widely used fundamental mode for pyramidal horn antennas [15].

<figure><img src="image_11.png" /><figcaption>Figure 11.3: Rectangular horn antenna and definition of its physical surface aperture.</figcaption></figure>

Therefore, in a coordinate system with the z-axis perpendicular to the surface aperture of the antenna, as shown in Fig. 11.3, the electric field generated by the antenna can be written as

$$ \vec{E}_{\text{ant}}(t, x, y, z) = \mathfrak{R} \left[ V_{\text{ant}} \vec{M}_{\text{ant}}(x, y) e^{i(kz - \omega_{\text{U}}t + \phi)} \right], \quad (11.16) $$

where $\phi$ is the phase of the electric field and $V_{\text{ant}}$ characterizes the amplitude of this field (Energy conservation implies that $V_0^2/2R_0 = V_{\text{ant}}^2/2Z_0$, where $Z_0 = 376.7 \, \Omega$ is the vacuum impedance.).

In general, horn antennas make an axial detection of the field, i.e. only on one single axis, the small axis of the rectangular horn, the $\hat{y}$ axis. The real antenna mode can be written as

$$ \vec{M}_{\text{ant}}(x, y) = m_{\text{TE}_{10}} \hat{y} \cos \left( \frac{\pi x}{A_{\text{eff}}} \right), \quad (11.17a) $$

where $A_{\text{eff}}$ is the effective width and $m_{\text{TE}_{10}}$ (dimension L$^{-1}$) is a normalisation factor that ensures the mode is normalized $\int dS_{\text{eff}} |\vec{M}_{\text{ant}}(x, y)|^2 = 1$ on the effective aperture of the antenna, i.e.

$$ m_{\text{TE}_{10}} = \sqrt{\frac{2}{S_{\text{eff}}}}. \quad (11.17b) $$

The effective aperture of the antenna $S_{\text{eff}}$ depends on the frequency $f_U = \omega_U/2\pi$ of the field and in general differs from the physical aperture of the antenna. Indeed, the effective area of the horn, which can be expressed as the product of an effective width $A_{\text{eff}}$ with an effective height $B_{\text{eff}}$, depends on the frequency as [15]

$$ S_{\text{eff}}(\omega_U) = A_{\text{eff}}(\omega_U) B_{\text{eff}}(\omega_U) = \frac{e_r \pi G(\omega_U) c^2}{\omega_U^2}, \quad (11.18) $$

where $e_r$ is the realistic efficiency of the antenna, which represents reflection losses inside the antenna, and where $G(\omega)$ is the frequency dependent gain of the antenna [15]. Note
that according to the IEEE Standards, the gain does not include reflection and polarization losses [15], this is why we separate $e_r$ from $G(\omega)$. The time-averaged electromagnetic power generated by the antenna is given by

$$P_{\text{out}} = \frac{1}{2Z_0} \int dS_{\text{eff}} |\vec{E}_{\text{ant}}|^2 = \frac{V_{\text{ant}}^2}{2Z_0}. \qquad (11.19)$$

We will now invert this reasoning and consider the antenna as a receiver that will output a voltage in response to an electric field

$$\vec{E}(t, x, y, z) = \Re[\vec{U}(x, y, z)e^{-i\omega_U t}] = |\vec{U}| \hat{e}_U \cos(\omega_U t + \varphi), \qquad (11.20)$$

where $\varphi$ is an irrelevant phase, $\hat{e}_U$ is the unit polarization vector of $\vec{U}$, $|\vec{U}| = (\vec{U} \cdot \vec{U}^\dagger)^{1/2}$ denotes the complex modulus. The antenna will measure electric fields that are propagating into the $\hat{z}$ direction and matching its mode $\vec{M}_{\text{ant}}(x, y)$. More precisely, by taking the dot product of Eq. (11.16) with $\vec{M}_{\text{ant}}$, one gets $V_{\text{ant}} \cos(\omega t - kz - \varphi) = \int dS_{\text{eff}} \vec{M}_{\text{ant}}(x, y) \cdot \vec{E}_{\text{ant}}(t, x, y, z)$. Using the reciprocity properties of electromagnetism, one can replace $\vec{E}_{\text{ant}}$ by the electric field we are trying to measure Eq. (11.20). This reasoning shows that the antenna will output a voltage proportional to

$$V_{\text{ant}}(z) \cos(\omega_U t + \varphi'(z)) = \cos(\omega_U t + \varphi) \int dS_{\text{eff}}(\omega_U) |\vec{U}(x, y, z)| \hat{e}_U \cdot \vec{M}_{\text{ant}}(x, y), \quad (11.21)$$

where $\varphi'(z) = -kz - \varphi$ and the 2D integral is carried out in the $x-y$ plane at the $z$-coordinate of the antenna.

Let us factorize the electric field generated by the dish at the location of the antenna provided by Eq. (11.15) in a voltage factor $V_{\text{dish}}$, a real mode factor $\vec{M}_{\text{dish}}$ and a plane wave factor, such that it takes the following form in the antenna aperture plane $(x, y)$

$$\vec{E}_{\text{dish}}(t, x, y, z) = V_{\text{dish}} \vec{M}_{\text{dish}}(x, y, z) \cos(\omega t + \varphi'(z)) \equiv |\vec{U}_{\text{dish}}(x, y, z)| \hat{e}_U \cos(\omega t + \varphi), \qquad (11.22a)$$

where $\vec{M}_{\text{dish}}(x, y, z)$ is normalized over the infinite surface of the antenna plane $S_\infty$ at $z = z_{\text{ant}}$ i.e

$$\int dS_\infty |\vec{M}_{\text{dish}}(x, y, z_{\text{ant}})|^2 = 1. \qquad (11.22b)$$

The constant amplitude $V_{\text{dish}}$ in Eq. (11.22a) can be computed through energy conservation: the total energy generated by the dish should equal the energy received on the infinite antenna plane $S_\infty$. Then, using Eqs. (11.2) and (11.22b), we have

$$P_{\text{dish}} = \int dS_{\infty} \frac{|\vec{U}_{\text{dish}}|^2}{2Z_0} = \frac{V_{\text{dish}}^2}{2Z_0} \Rightarrow V_{\text{dish}} = \sqrt{\frac{4Z_0 A_{\text{dish}} \rho_{\text{DM}} c}{3}} \chi. \qquad (11.23)$$
Then, in our experimental scheme, the time averaged power generated by the horn antenna is obtained by combining Eqs. (11.19), (11.21) and (11.22a) and is provided by

$$P_{\text{rec}}(z_{\text{ant}}, \omega_{\text{U}}) = \frac{V_{\text{dish}}^2}{2Z_0} \left( \int dS_{\text{eff}}(\omega_{\text{U}}) \vec{M}_{\text{dish}}(x, y, z_{\text{ant}}) \cdot \vec{M}_{\text{ant}}(x, y) \right)^2, \quad (11.24)$$

where $\vec{M}_{\text{dish}}$ depends explicitly on $z_{\text{ant}}$ and where the integral is performed at every position $(x, y)$ on the effective surface of the horn antenna.

Then, the ratio between the time averaged power measured by the horn antenna and the total power emitted by the dish is simply given by

$$\gamma(z_{\text{ant}}, \omega_{\text{U}})_{\text{Overlap}} = \left( \int dS_{\text{eff}}(\omega_{\text{U}}) \vec{M}_{\text{dish}}(x, y, z_{\text{ant}}) \cdot \vec{M}_{\text{ant}}(x, y) \right)^2, \quad (11.25)$$

where $\vec{M}_{\text{dish}}$ directly depends on the Kirchhoff integral Eq. (11.15).

### 11.3.2 Computation using the antenna factor

Another approach to estimate the output of the horn antenna is to consider the antenna factor (AF), a characteristic of the antenna provided by the manufacturer.

The AF is defined by considering an incident plane wave, or in other words, an incoming electromagnetic field whose mode is constant over the aperture of the antenna. It is defined by (see e.g. [16])

$$\text{AF}(\omega_{\text{U}}) = \frac{|\vec{U}_{\text{dish}}(\omega_{\text{U}})|}{V_0} = \frac{V_{\text{dish}} M_{\text{dish}}}{V_0}, \quad (11.26)$$

where $V_0$ is the generated tension, $\vec{U}_{\text{dish}}$ the amplitude of the incoming electric field, $\vec{M}_{\text{dish}} = M_{\text{dish}}\hat{e}_{\text{U}}$ from Eq. (11.22a), with $M_{\text{dish}}$ the value of the constant mode of the electric field on the effective antenna area. The AF depends on the frequency of the incident field. It is measured experimentally and therefore takes into account any loss inside the antenna, represented by the $e_r$ parameter.

Using energy conservation, we immediately get the expression of time-averaged generated power of the antenna

$$P_{\text{rec}}(z_{\text{ant}}, \omega_{\text{U}}) = \frac{V_0^2}{2R_0} = \frac{V_{\text{dish}}^2 M_{\text{dish}}^2(z_{\text{ant}}, \omega_{\text{U}})}{2R_0 A F(\omega_{\text{U}})^2}. \quad (11.27)$$

However, the definition of the antenna factor assumes perfectly aligned polarization modes [16], which is not necessarily true in our case. Indeed, we assume a antenna polarization mode along the $\hat{y}$ axis while the electric field emitted by the dish has its polarization lying in the $(x, y)$ plane.

If we assume the DP polarization to be isotropically distributed, the polarization of the electric field emitted by the dish has a cylindrical symmetry and a linearly polarized antenna will be sensitive only to half of the power from the electric field. Taking into account this additional factor, the ratio $\gamma$ becomes

$$\gamma(z_{\text{ant}}, \omega_{\text{U}})_{\text{AF}} = \frac{Z_0 M_{\text{dish}}^2(z_{\text{ant}}, \omega_{\text{U}})}{2 R_0 A F(\omega_{\text{U}})^2}. \quad (11.28)$$
From there, we have all the necessary ingredients to estimate the power received by the antenna as function of its position, as long as $z_{\text{ant}} < z_{\text{plane}}$. Indeed, Eqs. (11.25) and (11.28) give us what portion of the emitted power by the dish Eq. (11.2) is actually received in the antenna wires. Therefore, from a given constraint of an experiment which assumes 100% of power emitted is received by the antenna, one can infer how to correct this constraint using this ratio. In particular, in Chapter 18, we will use these results to derive the expected power received by the antenna in *SHUKET* [1], an experiment using a dish and a horn antennas to detect the small electric field Eq. (11.1). In [1], it is assumed that $\gamma = 1$, i.e all the power generated by the dish is transmitted to the antenna. Using Eqs. (11.25) and (11.28), we will show that this is an over optimistic result and that the power received by the antenna is multiple orders of magnitude lower.
[1] Pierre Brun, Laurent Chevalier, and Christophe Flouzat. "Direct Searches for Hidden-Photon Dark Matter with the SHUKET Experiment". In: Phys. Rev. Lett. 122 (20 May 2019), p. 201801. DOI: 10.1103/PhysRevLett.122.201801. URL: https://link.aps.org/doi/10.1103/PhysRevLett.122.201801.

[2] J. Suzuki et al. "Experimental search for hidden photon CDM in the eV mass range with a dish antenna". In: Journal of Cosmology and Astroparticle Physics 2015.09 (Sept. 2015), p. 042. DOI: 10.1088/1475-7516/2015/09/042. URL: https://dx.doi.org/10.1088/1475-7516/2015/09/042.

[3] Stefan Knirck et al. "First results from a hidden photon dark matter search in the meV sector using a plane-parabolic mirror system". In: Journal of Cosmology and Astroparticle Physics 2018.11 (Nov. 2018), p. 031. DOI: 10.1088/1475-7516/2018/11/031. URL: https://dx.doi.org/10.1088/1475-7516/2018/11/031.

[4] N. Tomita et al. "Search for hidden-photon cold dark matter using a K-band cryogenic receiver". In: Journal of Cosmology and Astroparticle Physics 2020.09 (Sept. 2020), p. 012. DOI: 10.1088/1475-7516/2020/09/012. URL: https://dx.doi.org/10.1088/1475-7516/2020/09/012.

[5] T. Mizumoto et al. "Experimental search for solar hidden photons in the eV energy range using kinetic mixing with photons". In: Journal of Cosmology and Astroparticle Physics 2013.07 (July 2013), p. 013. DOI: 10.1088/1475-7516/2013/07/013. URL: https://dx.doi.org/10.1088/1475-7516/2013/07/013.

[6] K. Ramanathan et al. "Wideband Direct Detection Constraints on Hidden Photon Dark Matter with the QUALIPHIDE Experiment". In: Phys. Rev. Lett. 130 (23 June 2023), p. 231001. DOI: 10.1103/PhysRevLett.130.231001. URL: https://link.aps.org/doi/10.1103/PhysRevLett.130.231001.

[7] S. Kotaka et al. "Search for Dark Photon Dark Matter in the Mass Range 74 - 110 μeV with a Cryogenic Millimeter-Wave Receiver". In: Phys. Rev. Lett. 130 (7 Feb. 2023), p. 071805. DOI: 10.1103/PhysRevLett.130.071805. URL: https://link.aps.org/doi/10.1103/PhysRevLett.130.071805.

[8] Stefan Knirck et al. First Results from a Broadband Search for Dark Photon Dark Matter in the 44 to 52 μeV range with a coaxial dish antenna. 2023. arXiv: 2310.13891 [hep-ex].

[9] Fayez Bajjali et al. "First results from BRASS-p broadband searches for hidden photon dark matter". In: Journal of Cosmology and Astroparticle Physics 2023.08 (Aug. 2023), p. 077. ISSN: 1475-7516. DOI: 10.1088/1475-7516/2023/08/077. URL: http://dx.doi.org/10.1088/1475-7516/2023/08/077.

[10] Dieter Horns et al. "Searching for WISPy cold dark matter with a dish antenna". In: Journal of Cosmology and Astroparticle Physics 2013.04 (Apr. 2013), pp. 016-016. ISSN: 1475-7516. DOI: 10.1088/1475-7516/2013/04/016. URL: http://dx.doi.org/10.1088/1475-7516/2013/04/016.

[11] Jordan Gué et al. "Improved modeling for dark photon detection with dish antennas". In: Phys. Rev. D 110.3 (2024), p. 035001. DOI: 10.1103/PhysRevD.110.035001. arXiv: 2403.13448 [hep-ex].

[12] John David Jackson. *Classical electrodynamics; 2nd ed.* New York, NY: Wiley, 1975. URL: https://cds.cern.ch/record/100964.
[13] J-Y Vinet. *The Virgo Physics Book: Optics and related Topics*. 2020.

[14] Haipeng An et al. “Direct Detection of Dark Photon Dark Matter Using Radio Telescopes”. In: *Phys. Rev. Lett*. 130 (18 May 2023), p. 181001. DOI: 10.1103/PhysRevLett.130.181001. URL: https://link.aps.org/doi/10.1103/PhysRevLett.130.181001.

[15] Constantine. A. Balanis. *Antenna theory: Analysis and Design, 3rd edition*. John Wiley and Sons, 2005.

[16] James S. Mclean, Robert Sutton, and Robert William Hoffman. “Interpreting Antenna Performance Parameters for EMC Applications: Part 3: Antenna Factor”. In: 2002. URL: https://api.semanticscholar.org/CorpusID:1097498.
# Part V

## Search through tests of the equivalence principle
# Chapter 12

# Classical tests of the equivalence principle

## 12.1 Experiments testing the equivalence principle

As described in Section 3.1, EP is not a fundamental symmetry of the universe, and many Beyond the Standard Model (BSM) theories predict its violation [1–6], in particular by the spacetime variation of fundamental constants such as the fine structure constant $\alpha$ or the electron mass $m_e$.

In order to quantify such violation, one can measure the acceleration of two test bodies A and B in an external gravitational field and compare them through the so-called Eötvös parameter defined as

$$ \eta = 2 \frac{(\vec{a}_A - \vec{a}_B)}{|\vec{a}_A + \vec{a}_B|} \cdot \hat{e}_{\text{meas}} \qquad (12.1) $$

where $\hat{e}_{\text{meas}}$ is the sensitive axis of the apparatus considered.

### 12.1.1 MICROSCOPE

The MICROSCOPE mission consists in a micro-satellite launched in circular orbit around Earth at a 710 km altitude [7]. It is the first mission aiming at measuring a violation of EP in space, taking advantage of a very quiet environment and long durations of free fall. It was launched on 25 April 2016 for a total duration of 120 orbits of the satellite around Earth [8], which as we shall see in Chapter 18 corresponds approximately to $8.1 \times 10^6$ s integration time. The satellite measured the difference of acceleration between a Platinum and a Titanium based test masses and it is sensitive to both Eqs. (8.12) and (8.40). Their final constraint on $\eta$ is [9]

$$ \eta = (-1.5 \pm 2.8) \times 10^{-15}, \qquad (12.2) $$

at a $1\sigma$ confidence-level.

### 12.1.2 Eötvös torsion balances

Torsion balance is an apparatus able to measure very weak forces, once used by Cavendish to measure the gravitational constant $G$, and now used for the measurement of the weak EP. In usual torsion balances experiments, a bar is suspended by a thin fiber, which acts as a torsion
spring. If a force is applied to one of the ends of the bar with direction perpendicular to it, the bar rotates and twists the fiber. This way, the magnitude of this unknown force can be derived by measuring the rotation angle of the bar. By setting two different test masses with two different compositions at each end of the bar (but with the same mass), one can measure the differential gravitational acceleration that act on them and therefore a violation of the equivalence principle. The current best constraint on the Eötvös parameter from Torsion balances experiments is [10]

$$ \eta = (0.3 \pm 1.8) \times 10^{-13}, \qquad (12.3) $$

at a $1\sigma$ confidence-level.

### 12.1.3 Lunar Laser Ranging

Thanks to retro-reflectors array installation on the surface of the Moon by several american and russian space missions, starting from 1969, it is now possible to measure the Earth-Moon distance by sending light signals to this array and measure the round trip time of photons. This is what is called Lunar Laser Ranging (LLR).

Both the Earth and the Moon have different compositions and are in the solar gravitational potential, a violation of the EP would imply the two systems fall at different rates towards the Sun. This would distort the lunar orbit, which could be measuring using LLR [11].

Since it would be unsatisfactory to model planetary and stellar bodies as point-like particles, we must include the self gravitation energy of each body, which is covered by the strong EP. In other words, LLR measurements provide us a constrain on a combination of strong and weak EP, but not on one of them individually.

## 12.2 Acceleration difference induced by oscillating rest mass

As derived in Sections 8.1 and 8.2.4, axion and dilaton fields can produce oscillations of the rest mass and transition frequencies of atoms from their different couplings with SM fields.

We make the calculations in the laboratory frame, where the gradient of the DM field is not neglected. The rest mass and transition frequency oscillations of an atom A are shown in Eq. (8.5), which leads to an UFF violating acceleration Eq. (8.12)

$$ \vec{a}_A(t, \vec{x}) = \left[ \omega_{\phi} \vec{v}_A - \vec{k}_{\phi} c^2 \right] \frac{\sqrt{16\pi G \rho_{\text{DM}} [Q_M^A]_d}}{\omega_{\phi} c} \sin(\omega_{\phi} t - \vec{k}_{\phi} \cdot \vec{x} + \Phi) \quad (12.4a) $$

$$ \approx \frac{\sqrt{16\pi G \rho_{\text{DM}} \vec{v}_{\text{DM}} [Q_M^A]_d}}{c} \sin(\omega_{\phi} t - \vec{k}_{\phi} \cdot \vec{x} + \Phi) \quad (12.4b) $$

where, as explained below Eq. (8.12), we immediately got rid off the transition frequency contribution to the acceleration, as it is suppressed by a factor $\hbar\omega_A^0/m_A^0c^2 \sim 10^{-10}$. At the second line, we used the fact that $-\vec{k}_{\phi}c^2 = \omega_{\phi}\vec{v}_{\text{DM}}$ and $v_{\text{DM}} \gg v_A$. It follows a differential acceleration between two macroscopic bodies of different composition A and B with the same initial velocities and position given by

$$ \Delta\vec{a}(t, \vec{x})_d \approx \frac{\sqrt{16\pi G \rho_{\text{DM}} \vec{v}_{\text{DM}}}}{c} ([Q_M^A]_d - [Q_M^B]_d) \sin(\omega_{\phi} t - \vec{k}_{\phi} \cdot \vec{x} + \Phi). \quad (12.5a) $$
Similarly, a differential acceleration between two bodies of different composition arises through the axion-gluon coupling, as derived in Section 8.2.4 of the form

$$ \Delta\vec{a}(t, \vec{x})_a \approx \frac{16\pi G \rho_{\text{DM}} \vec{v}_{\text{DM}} E_P^2}{f_a^2 \omega_a c^2} ([Q_M^A]_a - [Q_M^B]_a) \sin(2\omega_a t - 2\vec{k}_a \cdot \vec{x} + 2\Phi) . \quad (12.5b) $$

As we have derived it in Sections 8.1 and 8.2.4, the charges $[Q_M]_d$, $[Q_M]_a$ are atom dependent, as they depend on the mass and charge number of the atom. Therefore, the differential accelerations Eq. (12.5) are non zero, and lead to a violation of the UFF, i.e $\eta \neq 0$ in Eq. (12.1). The various experiments listed in Sec. 12.1 are sensitive to such non-zero differential acceleration, therefore they can put constraint on dilaton and axion couplings to SM fields. In Chapter 18, we will derive such sensitivities for MICROSCOPE. In this case, the two test-masses are not precisely at the same position but are approximately ~ 20 µm apart [7], which sets the typical length scale of the experiment. As we shall see in Chapter 17, MICROSCOPE is sensitive to oscillations up to 0.3 Hz. Following Eq. (7.28), this means that the propagation term of the plane wave in Eq. (12.5) can be neglected, and that we can approximate that both test-masses are at the same position.

In the next chapter, we will be interested in atom interferometry, another experimental tool which could detect such UFF violating accelerations.
# Chapter 13

## Atom interferometry as a quantum test of the equivalence principle

As derived in the previous chapter, the time-dependent mass and internal frequency pro-
duces a differential acceleration between two atoms. Standard UFF tests search for such
differential acceleration of macroscopic test masses. Atom interferometry (AI) constitutes
the quantum equivalent of such classical experiments. Therefore, as we shall see in this
chapter, AI experiments are also sensitive to oscillations from Eq. (12.5).

In this chapter, we will derive the observable signatures produced by oscillating atomic
rest mass and transition energy on various AI configurations. In the most generic case,
AI consists of generating interference between atomic wavepackets that are split and then
recombined using atom-laser interactions. We classify the various AI schemes into two
classes: (i) AI schemes using two photon transitions and (ii) AI schemes using single photon
transitions. In the former, the transition from a stable momentum state to another stable
momentum state is done using two photons, while in the latter, only one photon is needed
for such a transition. Those results are based on an article [12], published in *Physical Review*
*D.*

**13.1 Two-photon transitions : π/2 – π – π/2 setups**

**13.1.1 General principle**

In this section, we present the basics of atom interferometry and we will focus on the most
common AI setups : two-photon transition interferometers.

The first interferometric scheme considered is known as two-photon transition Raman
interferometry. Its sequence is depicted in Fig. 13.1. In this setup, two-level free falling atoms
A enter the interferometer with an initial momentum $\hbar k$ and in their internal energetic ground
state, noted $|g\rangle$ (i.e. their initial state is defined as $|g, \hbar k\rangle$). After entering the interferometer,
they interact with a pair of laser waves $L_1$ and $L_2$ with respective frequencies $\omega_{L1}$, $\omega_{L2}$, whose
energy difference is resonant with the transition of the $|g\rangle \rightarrow |e\rangle$ where $|e\rangle$ is the excited
state of the two-levels atom, i.e $\omega_A^0 = \omega_{L1} - \omega_{L2}$. This means that the atom first absorbs
a photon from $L_1$, which excites it to an intermediate state with energy $\hbar\omega_i \gg \hbar\omega_A^0$, and
afterwards, the atom stimulately emits a photon in $L_2$. This process splits the atoms into two
spatially distant wavepackets, which means the state of each atom becomes the superposition
of the two internal states: (i) the ground state that remains unchanged $|g, \hbar k\rangle$ and (ii) the excited state that has changed its momentum because of its interaction with the two lasers $|e, \hbar(k+k_{\text{eff}})\rangle$, where $\hbar \vec{k}_{\text{eff}} = \hbar(\vec{k}_{L1} - \vec{k}_{L2})$ denotes the effective momentum transfer, considering both lasers.

This whole experimental manoeuvre is called a $\pi/2$ laser pulse since the transition amplitude from ground state to excited state corresponds to a probability 1/2, i.e the state of the atoms is now a half-half superposition of the two different states.

The atom freely propagates inside the interferometer and at time $t = T$, both wavepackets undergo a $\pi$ pulse, which will invert the state of all atoms. In other words, depending on the internal state of the atom prior to the $\pi$ pulse, two state transitions happen: $|e, \hbar(k+k_{\text{eff}})\rangle \rightarrow |g, \hbar k\rangle$ and $|g, \hbar k\rangle \rightarrow |e, \hbar(k+k_{\text{eff}})\rangle$.

A final laser-atom interaction happens at time $t = 2T$, where a second $\pi/2$ pulse divides the two incoming wavepackets into four different ones: two of them are in the state $|g, \hbar k\rangle$ and the remaining two are in the state $|e, \hbar(k + k_{\text{eff}})\rangle$. The study of interference pattern between the wavepackets in the same state allows one to measure a phase shift difference. The full interferometric sequence is shown in Fig. 13.1.

The second two-photon transfer AI scheme considered is the two-photons transition Bragg-type interferometer. This scheme is similar to the Raman interferometry presented above except that the atoms remain in the same energy state during all the interferometric path, i.e, the laser pulses only change the external state of the atom (their momentum).

For both interferometers described above, the effective wavevector $k_{\text{eff}}$ depends on the setup of the experiment. If counterpropagating lasers are used, the atom absorbs a photon in one direction as a result of the interaction with the first laser, and emits another photon in the opposite direction during the interaction with the second laser, implying $k_{\text{eff}} = k_{L1} + k_{L2}$. In general, this effective wave vector is multiple orders of magnitude larger than the transition frequency of the atom, i.e $k_{\text{eff}} \gg \omega_A^0/c$, because both laser are usually operating at hundreds of THz frequencies (this is why the intermediate state is much higher in energy, see previous paragraphs and Fig. 13.1). On the opposite, if co-propagating laser waves are used in the experimental setup, the absorption and emission directions are the same, implying $k_{\text{eff}} = k_{L1} - k_{L2}$. In that case, $k_{\text{eff}} = \omega_A^0/c$.

### 13.1.2 Phase contributions

As in the previous chapter, we will make the detailed calculations considering a coupling between dilatons and SM fields, and we will generalize the main results to the axion case.

The calculation of the observable phase shift at the output of the interferometer follows closely the ones presented in [13]. This method relies on Feynman path integrals and can be used for Lagrangians which are at most quadratic in the position and velocity [13] and considering atomic plane waves at initial time $t = t_0$. This suggests that the calculations of the phase shift in our case can formally only be done in the galactic frame, where there is no $\cos(\omega t - \vec{k} \cdot \vec{x})$ terms<fnref n="1" />. Therefore, we will make all calculations in the galactic frame, where the field simply reads $\phi = \phi_0 \cos(\omega_\phi t + \Phi)$. In the next section, we will make a discussion on the possibility of making the calculations in the laboratory frame. Therefore, we consider an

<footnote n="1">¹As we shall see in the following, this framework can also work when the propagation is negligibly small, i.e when $kL \ll 1$, where L is the typical scale of the experiment.</footnote>
<figure><img src="image_12.png" /><figcaption>Figure 13.1: On the left, scheme of the two photon Raman transition considered in the present section. The atom starts in the ground state $|g\rangle$, absorbs a photon which makes it transition to a high energetic state $|i\rangle$, then emits a photon to transition to the excited state of the hyperfine transition $|e\rangle$. On the right, full Raman interferometric sequence, based on such two-photon transitions. At each atom-laser interaction, the laser pulses change both internal energy state and external momentum state of the wavepackets. The Bragg-type equivalent has the same spacetime diagram, but the atom wavepackets do not change their internal energy state. In black are shown the atomic paths without any perturbation, i.e straight lines (the gravitational field is neglected here). In green is shown the perturbed motion of the atoms induced by Eq. (13.6), with exaggerated amplitude of oscillation.</figcaption></figure>

atom A whose nominal rest mass and transition frequency $m_A^0, \omega_A^0$ are perturbed such that they oscillate in phase as in Eq. (8.5) but without the propagation term $\vec{k}_\phi \cdot \vec{x}$ in the phase.

$$m_A(t) = m_A^0 \left( 1 + \frac{\sqrt{16\pi G \rho_{\text{DM}} [Q_M^A]_d}}{\omega_\phi c} \cos(\omega_\phi t + \Phi) \right) \quad (13.1a)$$

$$\omega_A(t) = \omega_A^0 \left( 1 + \frac{\sqrt{16\pi G \rho_{\text{DM}} [Q_\omega^A]_d}}{\omega_\phi c} \cos(\omega_\phi t + \Phi) \right), \quad (13.1b)$$

By the same procedure as in Section 8.1.3, the macroscopic Lagrangian is

$$\mathcal{L}_A = - \left( m_A(t)c^2 + \hbar\omega_A(t) \right) \left( 1 - \frac{v_A^2}{2c^2} \right), \quad (13.2a)$$

which induces an acceleration on the atom A

$$\vec{a}_A(t) = \frac{\sqrt{16\pi G \rho_{\text{DM}} \vec{v}_A [Q_M^A]_d}}{c} \sin(\omega_\phi t + \Phi), \quad (13.2b)$$
<figure><img src="image_13.png" /></figure>

Figure 13.2: The various phase contributions along the AI path: in orange are the propagation phase shifts calculated along the wavepackets paths; in blue are the laser phase shifts, calculated at light-matter interaction vertices; and in green is the separation phase shift, calculated at the output of the interferometer.

where $\vec{v}_A$ is the velocity of the atom in the galactic frame, i.e $\vec{v}_A \approx \vec{v}_{DM}$ at leading order. One can already notice the similarity with Eq. (12.5a). Indeed, for situations where the gradient of the field is almost constant on the size of experiment $L$ considered, i.e $\lambda_{DM}^{\text{dB}} \gg L$, one can neglect the propagation term $\vec{k}_{\phi} \cdot \vec{x}$ in Eq. (12.5a) and the two accelerations are identical.

In practice, we will assume that the AI experiment is performed in a laboratory frame which has a velocity $\vec{v}_{DM}$ with respect to the DM rest frame where the Lagrangian takes the form given by Eqs. (13.1) and (13.2a), such that $\vec{v}_A = \vec{v}_{DM} + \vec{v}_0$ where the second term is the initial velocity of the atoms with respect to the laboratory reference frame. $\tilde{v}_0$ is the velocity impacting systematic effects such as gravity gradients or second-order Doppler shift of the atoms when interacting with light beams (see e.g [14] for a comprehensive list of such effects).

As derived in e.g. [13, 15, 16], there are three independent contributions to the total phase shift in an AI: (i) the separation phase noted $\Phi_u$ which corresponds to a spatial coincidence between the two output wavepackets, (ii) the laser phase $\Phi_\ell$ which gathers the additional phase factors of the laser, due to displacement of the light-matter interaction vertices and (iii) the propagation phase shift denoted $\Phi_s$ which is essentially the phase accumulated by the atom wavepacket over the full interferometric path. All of these contributions must be calculated accurately to predict the phase shift at the output of the interferometer. Those contributions are summarized in Fig. 13.2.

Finally, we will make all calculations assuming a Raman interferometer, such that the transition frequency of the atom is relevant. When relevant, we will show how the various Bragg interferometer phases differ from the Raman AI. We also make all calculations at first order in the perturbations.

We describe the quantum state of the different atom wavepackets by a wavefunction $\Psi$,
which we break down into two different wavefunctions $\Psi_I$ and $\Psi_{II}$, where I,II respectively the up/down paths in Fig. 13.1, depending on the classical path the atom followed. Considering free particles at the input of the interferometer, the atomic plane waves have the form [13]

$$ \Psi_{\text{init}}(t_0, \vec{x}_0) = \Psi_0 e^{i\Phi(t_0, \vec{x}_0)}, \qquad (13.3) $$

with the amplitude of the wavefunction $\Psi_0, \Phi(t_0, \vec{x}_0) = \vec{k} \cdot \vec{x}_0 - \omega t_0 - \varphi$ with $\omega, \vec{k}, \varphi$ respectively its angular frequency, wavevector and constant phase. At the output of the interferometer, at a time $T_f \ge 2T$, considering all the different phase contributions listed previously, the atomic plane wave along the trajectory j read

$$ \Psi_j(T_f, \vec{x}_f) = \Psi_{\text{init}}(T_f, \vec{x}_f) \times e^{i\Phi_j}, \qquad (13.4) $$

where $\Phi_j = \Phi_{sj} + \Phi_{lj} + \Phi_{uj}$ represents the trajectory dependent phase factor. Then, at some detection time $T_d \ge T_f$, a detector measures by fluorescence the number of atoms on each quantum state, which is essentially the measurement of the probability that the two wavepackets are in the same quantum state, i.e

$$ \int |\Psi_I(T_d, \vec{x}_d) + \Psi_{II}(T_d, \vec{x}_d)|^2 dS, \qquad (13.5) $$

where the integral is taken over the detector area S. Note that for Raman schemes, two output quantum states can differ by their internal energy and momentum state, while in Bragg, they only differ by their momentum. Plugging Eq. (13.4) into Eq. (13.5) and neglecting loss of contrast due to decoherence, we find that the measurement result is proportional to $(1 + \cos \Delta\Phi)$ where $\Delta\Phi = \Phi_I - \Phi_{II}$ is the phase difference between the two wave functions at $T_d$. For simplicity, we consider $T_d \equiv 2T$, because any additional phase shifts cumulated during $T_d - T_f$ is negligible (as in general $T_d - T_f \ll 2T$, and because the two wavefunctions on each output port have almost the same position, momentum and energy.).

### 13.1.3 Calculation of the phase shift observable

We first derive the perturbed equations of motion of the atom following the perturbation to the acceleration Eq. (13.2b), to get the motion followed by the different wavepackets along the trajectories presented in Fig. 13.1. To simplify the reading of the next set of calculations, we set $\sqrt{16\pi G\rho_{DM}/\omega_{\phi}c} \equiv X_{DM}$ (which is dimensionless). These equations of motion read

$$ \vec{v}_A(t, t_0) \approx \vec{v}_{DM} \left( 1 - X_{DM}[Q_M^A]_d \left( \cos(\omega_{\phi}t + \Phi) - \cos(\omega_{\phi}t_0 + \Phi) \right) \right), \qquad (13.6a) $$

$$ \vec{x}_A(t, t_0) \approx \vec{x}_0 + \vec{v}_{DM}(t - t_0) - \frac{X_{DM}\vec{v}_{DM}[Q_M^A]_d}{\omega_{\phi}} \left[ \sin(\omega_{\phi}t + \Phi) - \sin(\omega_{\phi}t_0 + \Phi) - \omega_{\phi}(t - t_0)\cos(\omega_{\phi}t_0 + \Phi) \right], \qquad (13.6b) $$

at first order in the perturbation $[Q_M^A]_d$, with $\vec{x}_0 = \vec{x}(t_0)$ the initial position of the atom when entering the interferometer. In the following, we set $x_0 = t_0 = 0$ and $\vec{v}_{DM} = v_{DM}\hat{v}$. Note that in Eq. (13.6), one of the terms arising in the position of the atom is $\propto \omega_{\phi}t$, which induces a temporal linear drift of the position of the atom, due to its initial inertia.
As already discussed in [17], the velocity kick experienced by the atom during its interac-
tion with the laser beams is perturbed by the effective mass of the atom at the time of the kick.
In addition to that effect, the beams being locked to a given frequency reference, their own
frequency oscillates as $\omega_L(t) = \omega_L^0(1 + (\sqrt{16\pi G\rho_{DM}}[Q_\omega^L]_d/\omega_\phi c) \cos(\omega_\phi t + \Phi))$, neglecting the
additional phase coming from the propagation of photons to the atom<fnref n="2" />. Therefore, if the
atom interacts with the laser beam at time $t$, still considering $[Q_M^A]_d, [Q_\omega^L]_d \ll 1$,

$$
v_{\text{kick}}^{A}(t) = \frac{\hbar k_{\text{eff}}(t)}{m_{A}(t)} \approx \frac{\hbar k_{\text{eff}}}{m_{A}^{0}} \left( 1 + X_{\text{DM}}([Q_{\omega}^{L}]_{d} - [Q_{M}^{A}]_{d}) \cos(\omega_{\phi} t + \Phi) \right) \equiv v_{\text{kick}, 0}^{A} + \delta v_{\text{kick}}^{A}(t), \quad (13.7)
$$

where $v_{\text{kick}, 0}^A = \hbar k_{\text{eff}}/m_A^0$ is the nominal kick velocity, i.e without perturbation, and $\delta v_{\text{kick}}^A(t)$ is the perturbed contribution to the total velocity kick imparted to the atom. In Appendix B, we compute explicitly the position and velocity of the wavepackets along the various interferometric paths, which will be useful for the computation of the various phase contributions.

**Propagation phase contribution** The first component of phase shift is the one coming from the phase accumulated by atoms throughout the whole interferometric paths taking into account modified equations of motion and perturbed kicks. In the special case of quadratic Lagrangian in the position and velocity of the atom at maximum, this phase is by the principle of least action the integral of the Lagrangian over the path from initial point $i$ to final point $f$ ($t_i, x_i \rightarrow t_f, x_f$)

$$
\Phi_s = \frac{1}{\hbar} \int_{t_i}^{t_f} L(x, \dot{x}) dt, \quad (13.8)
$$

where the Lagrangian is defined in Eq. (8.11b). The internal energy term of the Lagrangian
being associated with the oscillation of the transition energy, it only contributes when the
atom is on the excited state, i.e on paths 2 and 3 in Fig. 13.1. Then, the phase accumulated
by the atoms on the path I of the interferometer is

$$
\begin{equation}
\begin{split}
\Phi_{\text{sI}} = -\omega_A^0 & \int_0^T dt (1 + X_{\text{DM}}[Q_\omega^A]_d \cos(\omega_\phi t + \Phi)) \left(1 - \frac{|\vec{v}_A^{(2)}(t)|^2}{2c^2}\right) - \\
& \frac{m_A^0 c^2}{\hbar} \left[ \int_0^T dt (1 + X_{\text{DM}}[Q_M^A]_d \cos(\omega_\phi t + \Phi)) \left(1 - \frac{|\vec{v}_A^{(2)}(t)|^2}{2c^2}\right) + \right. \\
& \qquad \left. \int_T^{2T} dt (1 + [Q_M^A]_d \cos(\omega_\phi t + \Phi)) \left(1 - \frac{|\vec{v}_A^{(4)}(t)|^2}{2c^2}\right) \right],
\end{split}
\tag{13.9a}
\end{equation}
$$

<footnote n="2">²In the following sections, the travelling distance of the photon to reach the freely falling atoms will be of the order of 100 m maximum, which would induce a significant phase for oscillation frequencies $\omega_\phi \gtrsim 10^5$ rad/s, way above the DM frequencies of interest for this study, see Chapter 18.</footnote>
while the atom wavepacket on the path II accumulates a phase

$$ \Phi_{sII} = -\omega_A^0 \int_T^{2T} dt (1 + X_{DM}[Q_\omega^A]_d \cos(\omega_\phi t + \Phi)) \left(1 - \frac{|\vec{v}_A^{(3)}(t)|^2}{2c^2}\right) - \frac{m_A^0 c^2}{\hbar} \left[ \int_0^T dt \left(1 + X_{DM}[Q_M^A]_d \cos(\omega_\phi t + \Phi)\right) \left(1 - \frac{|\vec{v}_A^{(1)}(t)|^2}{2c^2}\right) + \int_T^{2T} dt \left(1 + X_{DM}[Q_M^A]_d \cos(\omega_\phi t + \Phi)\right) \left(1 - \frac{|\vec{v}_A^{(3)}(t)|^2}{2c^2}\right) \right], \quad (13.9b) $$

where $v^{(1)}, v^{(2)}, v^{(3)}, v^{(4)}$ are respectively the atom velocities along portions 1, 2, 3 and 4 in Fig. 13.1, which can all be recovered explicitly using the set of calculations described in Appendix B.

Then, the propagation phase shift between the two perturbed trajectories is

$$ \Phi_s = \Phi_{sI} - \Phi_{sII} = -\frac{4X_{DM}}{\omega_\phi} \left[ k_{\text{eff}} \left( v_{\text{DM}}\hat{e}_v \cdot \hat{e}_{\text{kick}} + \frac{\hbar k_{\text{eff}}}{2m_A^0} \right) [Q_M^A]_d + \omega_A^0 [Q_\omega^A]_d \right] \sin^2\left(\frac{\omega_\phi T}{2}\right) \sin(\omega_\phi T + \Phi) + 4k_{\text{eff}}TX_{DM}[Q_\omega^L]_d \left( v_{\text{DM}}\hat{e}_v \cdot \hat{e}_{\text{kick}} + \frac{\hbar k_{\text{eff}}}{2m_A^0} \right) \sin\left(\frac{\omega_\phi T}{2}\right) \sin\left(\frac{\omega_\phi T}{2} + \Phi\right) + \mathcal{O}\left(\left(\frac{v_{\text{DM}}}{c}\right)^2\right), \quad (13.10a) $$

where we used $m_A^0 v_{\text{kick}} = \hbar k_{\text{eff}}$ at zeroth order in the perturbation, following Eq. (13.7) and where we defined $\vec{v}_{\text{kick}} = v_{\text{kick}}\hat{e}_{\text{kick}}$. The terms $\propto (v_{\text{DM}}^2/c^2)$ arise from the small contribution of the atoms velocity to the internal state kinetic energy. The Bragg propagation phase can be simply recovered by remembering that in this configuration, one effectively manipulates one-level atom, i.e $\omega_A^0 = 0$.

**Laser phase contribution** We now consider the phase from light-matter interaction between the laser and the atoms. In both Bragg and Raman schemes, we assume that the atoms are freely falling inside a "vacuum tower" where the lasers are on ground, located on a mount at initial position $x_G(0) = 0$, while a retro-reflective mirror at initial position $x_M(0) = L$ is used to reflect the beams in order to create the counter-propagating scheme. In the Bragg case, only one laser beam is used and retro-reflected, while for the Raman case, two beams $L_1, L_2$ with respective frequencies $\omega_{L_1}, \omega_{L_2}$ (such that $\omega_{L_1} - \omega_{L_2} = \omega_A^0$) are retro-reflected and the atomic wavepackets interact only with $L_1$ going up and $L_2$ going down. Considering that the whole tower on Earth is freely falling during the entire interferometric process, its own mass composition is affected by the oscillating mass behavior Eq. (13.1), i.e $m_M = m_M^0 (1 + Q_M^M \cos(\omega_\phi t + \Phi))$, with $m_M^0, Q_M^M$ respectively its unperturbed mass and mass charge, implying it follows the same perturbed equations of motion as the atom Eq.(13.6). Therefore, both retro-reflective mirror and mount on which laser rest upon are oscillating together<fnref n="3" />. At each spacetime points of light-matter interaction is associated a phase, which contributes to the laser phase shift. These points are denoted A, B, C, D₂ in Fig. 13.1, since we

<footnote n="3">³At next-to-leading order, L is also oscillating as the interatomic bonds oscillate too.</footnote>
assume the measurement of the interference pattern between the two wavepackets ending up on the internal ground state $|g\rangle$.

Keeping in mind that $\Phi_L(t, \vec{x}) = \vec{k}_L \cdot \vec{x} - \omega_L t - \varphi$, with $\omega_L, \vec{k}_L$ the light angular frequency and wavevector respectively, and that $\vec{k}_L \cdot \vec{x} - \omega_L t = 0$ along a photon geodesic in a flat spacetime (which is the case at the surface of the Earth, see [18]), the total light phase felt by the atom along a path X is given by

$$ \Phi_{\ell X} = - \sum_{j=0}^{n} \varphi_j(t_i) = - \sum_{j=0}^{n} s_j \left( \varphi_j(t_i^{\text{down}}) - \varphi_j(t_i^{\text{up}}) \right), \qquad (13.11) $$

where we sum on the total number of interaction points of path X (occurring at times $\{0, T, 2T\}$), and where $t_i$ corresponds respectively to the time of laser at emission of the photon (the superscripts up and down refer to a photon coming from the up or down laser). At the second line, we took into account the fact that we have two-photon interactions, where the atom takes a photon from one of the laser and emits in the second one, and the $s_j = \pm 1$ parameter depends on the transition of the wavepacket at interaction j (it is +1 for a $|g\rangle \rightarrow |e\rangle$ transition and -1 for a $|g\rangle \rightarrow |e\rangle$ transition). Therefore, only the initial phase of the laser at the time of photon emission is needed. Since the laser frequency is locked on an atom ensemble whose frequency oscillates through Eq. (13.1), the phase, as the integral of the frequency, is also time dependent, i.e

$$ \varphi_j(t) = \int_0^t dt' \omega_L(t') \equiv \int_0^t dt' \omega_L^0 \left( 1 + X_{\text{DM}} Q_{\omega}^L \cos(\omega_{\phi} t' + \Phi) \right). \qquad (13.12) $$

<figure><img src="image_14.png" /><figcaption>Figure 13.3: Simple scheme on the computation of the time of emission of photons from both lasers.</figcaption></figure>

To compute the laser phase, we need to know the different times $t_i$ of emission of photons. For a given interaction time $t_{\text{int}}$, they are given by (see Fig. 13.3).

$$ t^u = t_{\text{int}} - \frac{2x_M(t_{\text{int}} - t'^u) - x_{\text{At}}(t_{\text{int}}) - x_G(t'^u - t^u)}{c}, \qquad (13.13a) $$

$$ t^d = t_{\text{int}} - \frac{x_{\text{At}}(t_{\text{int}}) - x_G(t_{\text{int}} - t^d)}{c}, \qquad (13.13b) $$

respectively for the beam retro-reflected (up) and not retro-reflected (down), where $x_M$ is the mirror coordinate, $x_G$ is the mount coordinate and $x_{\text{At}}$ is the perturbed vertex of the atomic path. In order to solve Eq. (13.13), the various $t^u, t'^u, t^d$ on the right side are treated as unperturbed times.
Then, for Raman AI, the total laser phase is

$$
\begin{align}
\Phi_{\ell} &= \Phi_{\ell I} - \Phi_{\ell II} \tag{13.14a} \\
&= \frac{4X_{\mathrm{DM}}}{\omega_{\phi}} \left( k_{\mathrm{eff}} v_{\mathrm{DM}} ([Q_{M}^{A}]_{d} - [Q_{M}^{M}]_{d}) \hat{e}_{v} \cdot \hat{e}_{\mathrm{kick}} + \frac{\hbar k_{\mathrm{eff}}^{2}}{m_{A}^{0}} [Q_{M}^{A}]_{d} + \omega_{A}^{0} [Q_{\omega}^{L}]_{d} \right) \sin^{2}\left(\frac{\omega_{\phi} T}{2}\right) \sin(\omega_{\phi} T + \Phi) + \right. \nonumber \\
&\qquad \left. 4 [Q_{\omega}^{L}]_{d} X_{\mathrm{DM}} \left( L \left( k_{\mathrm{eff}} - \frac{\omega_{A}^{0}}{c} \right) \sin^{2}\left(\frac{\omega_{\phi} T}{2}\right) \cos(\omega_{\phi} T + \Phi) - \frac{\hbar k_{\mathrm{eff}}^{2} T}{2m_{A}^{0}} \sin\left(\frac{\omega_{\phi} T}{2}\right) \sin\left(\frac{\omega_{\phi} T}{2} + \Phi\right) \right) \right). \tag{13.14b}
\end{align}
$$

The Bragg laser phase is obtained by setting $\omega_A^0 = 0$.

**Separation phase contribution** Following Eq. (13.5) and considering that, at the end of the classical paths at $t = T_d$, the wave packet of the path $i=\{I,II\}$ can be expressed as [16]

$$
\Psi_i = \Psi_0 e^{i(\Phi_{si} + \Phi_{ti} + \frac{\vec{p}_i}{\hbar}(\vec{x} - \vec{x}_i))}, \quad (13.15)
$$

where $\Phi_s, \Phi_\ell$ are respectively the propagation and laser phase contributions and where $\vec{p}_i, \vec{x}_i$ are respectively the momentum and position of the wavepacket just after the last $\pi/2$ pulse at $t = T_d$. Then, it can be shown easily that for small difference in momentum between the two wavepackets $\Delta\vec{p} = \vec{p}_I - \vec{p}_{II}$ (i.e $|\Delta\vec{p}|L_{det}/\hbar \ll 1$, with $L_{det}$, the size of the detector, which we show in Appendix B), the separation phase shift can be expressed as

$$
\Phi_u = \frac{\Delta \vec{p} \cdot \vec{x}_{\text{det,COM}}(2T)}{\hbar} - \frac{\Delta (\vec{p} \cdot \vec{x})(2T)}{\hbar}, \quad (13.16)
$$

where we compute the difference of momenta $\vec{p} = \hbar\vec{k}$ and position $\vec{x}$ between the two
wavepackets in the same energy state, i.e $|g\rangle$ in our calculation. $\vec{x}_{\text{det,COM}}$ represents the
detector center of mass position at time $t = T_d = 2T$, i.e in our case it is simply $\vec{x}_{\text{det,COM}}(2T) =
2\vec{v}_{\text{DM}}t$, since its equation of motion reads $\vec{x}_{\text{det,COM}}(t) = \vec{v}_{\text{DM}}t$ (its own oscillation would
induce a second order effect).

As we have shown in Appendix B, both positions and velocities of the wavepackets at
the end of paths 3 and 4 in Fig. 13.1 are different in our framework, therefore, the separation
phase is a function of both differences of positions and velocities at the end of the two paths

$$
\Delta \vec{x}(2T) = 4 \vec{v}_{\text{kick}} X_{\text{DM}} \left( [Q_{\omega}^{L}]_{d} T - \frac{[Q_{M}^{A}]_{d}}{\omega_{\phi}} \sin\left(\frac{\omega_{\phi} T}{2}\right) \sin\left(\frac{\omega_{\phi} T}{2}\right) \sin(\omega_{\phi} T + \Phi) \right) \quad (13.17\text{a})
$$

$$
\Delta \vec{v}(2T) = -4 \vec{v}_{\text{kick}} X_{\text{DM}} [Q_{\omega}^{L}]_{d} \sin^2\left(\frac{\omega_{\phi} T}{2}\right) \cos(\omega_{\phi} T + \Phi), \quad (13.17\text{b})
$$

implying that the difference in velocities is small (of order $[Q_\omega^L]_d$), justifying the form of the
separation phase Eq. (13.16). Then, the separation phase shift is

$$ \Phi_u = \frac{m_A^0}{\hbar} (\Delta\vec{v}(2T) \cdot \vec{x}_{\text{det.,COM}}(2T) - \Delta\vec{v}(2T) \cdot \vec{x}(2T) - \vec{v}(2T) \cdot \Delta\vec{x}(2T)) \quad (13.18a) $$

$$ = \frac{4k_{\text{eff}}X_{\text{DM}}}{\omega_{\phi}} \left( v_{\text{DM}}\hat{e}_v \cdot \hat{e}_{\text{kick}} \left[ [Q_M^A]_d \sin\left(\frac{\omega_{\phi}T}{2}\right) \sin(\omega_{\phi}T + \Phi) - \omega_{\phi}T[Q_{\omega}^L]_d \sin\left(\frac{\omega_{\phi}T}{2} + \Phi\right) \right] + \right. $$

$$ \left. \frac{\hbar k_{\text{eff}} \omega_{\phi} T}{m_A^0} [Q_{\omega}^L]_d \sin\left(\frac{\omega_{\phi}T}{2}\right) \cos(\omega_{\phi}T + \Phi) \right) \sin\left(\frac{\omega_{\phi}T}{2}\right), \quad (13.18b) $$

where $\vec{v}(2T), \vec{x}(2T)$ respectively correspond to the unperturbed velocity and position of the wavepackets at $t = 2T$.

### 13.1.4 Total phase shift

Adding all contributions of phase shift Eqs. (13.10b), (13.14a) and (13.18b), the total phase of the Raman AI is

$$ [\Phi_A^{\text{Raman}}]_d = \frac{4\sqrt{16\pi G \rho_{\text{DM}}}}{\omega_{\phi}^2 c} \left[ \left( k_{\text{eff}} v_{\text{DM}} \Delta[Q_M]_d \hat{e}_v \cdot \hat{e}_{\text{kick}} + \omega_A^0 \Delta[Q_{\omega}]_d \right) \sin^2\left(\frac{\omega_{\phi}T}{2}\right) \sin(\omega_{\phi}T + \Phi) + \right. \\ \left. \left( \frac{\hbar k_{\text{eff}}^2 [Q_M^A]_d}{2m_A^0} \sin(\omega_{\phi}T + \Phi) - k_{\text{eff}} \omega_{\phi} [Q_{\omega}^L]_d \left( L \left( 1 - \frac{\omega_A^0}{k_{\text{eff}}c} \right) + \frac{\hbar k_{\text{eff}} T}{m_A^0} \right) \cos(\omega_{\phi}T + \Phi) \right) \sin^2\left(\frac{\omega_{\phi}T}{2}\right) \right], \quad (13.19a) $$

where $\Delta[Q_M]_d = [Q_M^A]_d - [Q_M^L]_d$ and $\Delta[Q_{\omega}]_d = [Q_{\omega}^L]_d - [Q_{\omega}^A]_d$. For Bragg AI, one simply needs to set $\omega_A^0 = 0$ and the phase reads

$$ [\Phi_A^{\text{Bragg}}]_d = \frac{4\sqrt{16\pi G \rho_{\text{DM}}}}{\omega_{\phi}^2 c} \left[ k_{\text{eff}} v_{\text{DM}} \Delta[Q_M]_d \hat{e}_v \cdot \hat{e}_{\text{kick}} \sin^2\left(\frac{\omega_{\phi}T}{2}\right) \sin(\omega_{\phi}T + \Phi) + \right. \quad (13.19b) $$

$$ \left. \left( \frac{\hbar k_{\text{eff}}^2 [Q_M^A]_d}{2m_A^0} \sin(\omega_{\phi}T + \Phi) - k_{\text{eff}} \omega_{\phi} [Q_{\omega}^L]_d \left( L + \frac{\hbar k_{\text{eff}} T}{m_A^0} \right) \cos(\omega_{\phi}T + \Phi) \right) \sin^2\left(\frac{\omega_{\phi}T}{2}\right) \right], $$

at lowest order in $v_{\text{DM}}/c$. Similarly, using Eq. (8.40), we find the expected signal resulting from the axion-gluon coupling in a Raman AI, which reads

$$ [\Phi_A^{\text{Raman}}]_a = \frac{16\pi G \rho_{\text{DM}} E_P^2}{f_a^2 \omega_a^3 c^2} \left[ \left( k_{\text{eff}} v_{\text{DM}} \Delta[Q_M]_a \hat{e}_v \cdot \hat{e}_{\text{kick}} + \omega_A^0 \Delta[Q_{\omega}]_a \right) \sin^2(\omega_a T) \sin(2\omega_a T + \Phi) + \right. \\ \left. \left( \frac{\hbar k_{\text{eff}}^2 [Q_M^A]_a}{2m_A^0} \sin(2\omega_a T + \Phi) - k_{\text{eff}} \omega_a [Q_{\omega}^L]_a \left( L \left( 1 - \frac{\omega_A^0}{k_{\text{eff}}c} \right) + \frac{\hbar k_{\text{eff}} T}{m_A^0} \right) \cos(2\omega_a T + \Phi) \right) \sin^2(\omega_a T) \right], \quad (13.20a) $$
where $\Delta[Q_M]_a = [Q_M^A]_a - [Q_M^M]_a$ and $\Delta[Q_\omega]_a = [Q_\omega^L]_a - [Q_\omega^A]_a$. Again, for Bragg AI, one simply needs to set $\omega_A^0 = 0$ and the phase reads

$$[\Phi_A^{\text{Bragg}}]_a = \frac{16\pi G \rho_{\text{DM}} E_P^2}{f_a^2 \omega_a^3 c^2} \left[ k_{\text{eff}} v_{\text{DM}} \Delta[Q_M]_a \hat{e}_v \cdot \hat{e}_{\text{kick}} \sin^2(\omega_a T) \sin(2\omega_a T + \Phi) + \left( \frac{\hbar k_{\text{eff}}^2 [Q_M^A]_a}{2m_A^0} \sin(2\omega_a T + \Phi) - k_{\text{eff}} \omega_a [Q_\omega^L]_a \left( L + \frac{\hbar k_{\text{eff}} T}{m_A^0} \right) \cos(2\omega_a T + \Phi) \right) \sin^2(\omega_a T) \right]. \quad (13.20b)$$

As we derived it in Chapter 8, the mass and frequency charges in both dilaton and axion models are composition-dependent, therefore the terms appearing in the first line in Eqs. (13.19) and (13.20) are a signature of the violation of the Einstein EP. In the following, the terms proportional to $v_{\text{DM}}$, which are present in both Bragg and Raman configurations, will be denoted as the mass terms, while the terms proportional to $\omega_A^0$, only present in Raman configuration will be denoted as the frequency terms. Depending on the setup, the $Q_\omega^L$ and $Q_\omega^A$ charges can be the same, cancelling this term completely.

All the terms appearing at the second lines of the same equations, i.e the terms quadratic in the effective wavevector and respectively proportional to the mass charge of the atom A and the frequency charge of the laser and the terms proportional to the mirror coordinate L still remain, even in the universal charges case<fnref n="4" />. In fact, the terms quadratic in the effective wavevector can be understood as a non-local measurement, whose macroscopic counterpart would be to compare the free fall of two test masses with different velocities. In the following, we will neglect this term, since $v_{\text{DM}} \sim 10^5 \text{ m/s} \gg 10^{-2} \text{ m/s} \sim \hbar k_{\text{eff}}/m_A^0$. The argument is the same for the term $\propto L$. We will also drop this term in the following, as it is also much smaller than the leading order terms, in particular when computing a differential phase shift (see next paragraph). Therefore, in the Bragg case, the mass term is the dominant term.

Two-photon transition AI are commonly operating using hyperfine transitions of alkaline-Earth atoms (e.g the hyperfine transition of Rb, Cs, K atoms [19–21]). For both co-propagating and counter-propagating configurations, the laser beams are locked onto the optical transition of an atomic ensemble of the same species as the atoms in free fall inside the interferometer, and the frequency difference $\omega_{L1} - \omega_{L2} = \omega_A^0$ is provided by a radio-frequency source (in the GHz range). Then $Q_\omega^L$ is the charge of that source, which, depending on the experimental configuration, may or may not be the same as $Q_\omega^A$ the charge of the hyperfine transition of the atoms in the AI. As discussed previously, in the case of co-propagating lasers, the effective wavevector corresponds to the frequency transition of the atom, i.e $k_{\text{eff}} = \omega_A^0/c$. Therefore, the mass term in Eq. (13.19) and Eq. (13.20) is suppressed by a factor $v_{\text{DM}}/c$ compared to the frequency term. For counter-propagating laser beams, the effective wavevector $k_{\text{eff}} \gg \omega_A^0/c$, implying that the mass term in these equations is much bigger than in the co-propagating case. In that case, both mass and frequency terms are relevant and need to be taken into account. The mass term has already been derived in [17], in the case of oscillating mass coming from a coupling between matter and a classical oscillating dark matter field. However, the calculation in [17] was performed in the lab frame but only considering the velocity of

<footnote n="4">⁴At this point, the interpretation of such terms in terms of the Einstein Equivalence Principle starts to be subject to interpretation. One can argue that these would be a sign of violation of the UFF as they would still induce a non-zero phase shift when measuring a differential phase shift between two different atomic species. Another interpretation would be a violation of the Lorentz Position Invariance as those terms appear due to the variation of the laser phase, which now depends on the spacetime position of the interaction.</footnote>
the atoms in this frame, i.e $v_A \sim 10$ m/s, while we argue that another component, due to the galactic velocity $v_{DM} \sim 10^5$ m/s $\gg v_A$ should be taken into account, as we demonstrated it in the beginning of this section. Adding this contribution would improve the expected sensitivity to the charges $Q_M^A$ by several orders of magnitude. In addition, our calculations take into account the contribution from Earth oscillation (through $Q_M^M$) which was not the case in [17].

Note that in double diffraction interferometers, i.e when two pairs of laser beams transfer opposite momentum to the atom [22], such that the spatial separation between the two coherent wavepackets is twice as large, the total phase shift is the same with the change $k_{\text{eff}} \rightarrow 2k_{\text{eff}}$, as expected. In addition, in our calculation, the retro-reflective mirror is used to create the counter propagating waves, i.e we implicitly assumed a counter propagating scheme. For the co-propagating situation, one simply needs to set $L = 0$, as there is no mirror.

Dual atom interferometers using two atomic species A and B with different mass and frequency charges will measure the difference of the interferometric phases, whose amplitude is given by

$$ \left|\Delta\Phi_{AB}^{\text{Raman}}\right|_d = \frac{4\sqrt{16\pi G\rho_{\text{DM}}}}{\omega_{\phi}^2 c} \left( v_{\text{DM}}\hat{e}_v \cdot \hat{e}_{\text{kick}} \left( k_{\text{eff}}^A([Q_M^A]_d - [Q_M^M]_d) - k_{\text{eff}}^B([Q_M^B]_d - [Q_M^M]_d) \right) + \omega_A^0 ([Q_\omega^{L,A}]_d - [Q_\omega^A]_d) - \omega_B^0 ([Q_\omega^{L,B}]_d - [Q_\omega^B]_d) \right) \sin^2\left(\frac{\omega_\phi T}{2}\right) \quad (13.21a) $$

$$ \left|\Delta\Phi_{AB}^{\text{Raman}}\right|_a = \frac{16\pi G \rho_{\text{DM}} E_P^2}{f_a^2 \omega_a^3 c^2} \left( v_{\text{DM}}\hat{e}_v \cdot \hat{e}_{\text{kick}} \left( k_{\text{eff}}^A([Q_M^A]_a - [Q_M^M]_a) - k_{\text{eff}}^B([Q_M^B]_a - [Q_M^M]_a) \right) + \omega_A^0 ([Q_\omega^{L,A}]_a - [Q_\omega^A]_a) - \omega_B^0 ([Q_\omega^{L,B}]_a - [Q_\omega^B]_a) \right) \sin^2(\omega_a T), \quad (13.21b) $$

for Raman interferometers and

$$ \left|\Delta\Phi_{AB}^{\text{Bragg}}\right|_d = \frac{4\sqrt{16\pi G\rho_{\text{DM}}v_{\text{DM}}\hat{e}_v \cdot \hat{e}_{\text{kick}}}}{\omega_{\phi}^2 c} \left( k_{\text{eff}}^A([Q_M^A]_d - [Q_M^M]_d) - k_{\text{eff}}^B([Q_M^B]_d - [Q_M^M]_d) \right) \sin^2\left(\frac{\omega_{\phi}T}{2}\right) \quad (13.21c) $$

$$ \left|\Delta\Phi_{AB}^{\text{Bragg}}\right|_a = \frac{16\pi G \rho_{\text{DM}} E_P^2 v_{\text{DM}} \hat{e}_v \cdot \hat{e}_{\text{kick}}}{f_a^2 \omega_a^3 c^2} \left( k_{\text{eff}}^A([Q_M^A]_a - [Q_M^M]_a) - k_{\text{eff}}^B([Q_M^B]_a - [Q_M^M]_a) \right) \sin^2(\omega_a T), \quad (13.21d) $$

for Bragg interferometers, where we used Eqs. (13.19) and (13.20), neglected the subdominant terms $\propto k_{\text{eff}}^2$, $L$ and where $Q_\omega^{L,A}$, $Q_\omega^{L,B}$ are respectively the frequency charge of the beams used for the transition of the atomic species A and B in the Raman case. The mass terms of Eq. (13.21) are the quantum equivalent of the classical calculation Eq. (12.5). Indeed, for small $\omega_\phi T$, $\omega_a T$ and assuming $k_{\text{eff}}^A = k_{\text{eff}}^B = k_{\text{eff}}$, the Taylor expansion of the mass terms of Eq. (13.21) is related to Eq. (12.5) through $|\Delta\Phi_{AB}| = k_{\text{eff}}T^2|\Delta\vec{a}_{AB}|$, which is the usual relation between phase shift and acceleration in atom interferometers [13]. The second terms of Eq. (13.21) of Raman phase, proportional to the frequency charges $Q_\omega$ have no classical counterpart. Assuming mass and frequency charges of same order of magnitude, the frequency terms of Eq. (13.21) will dominate in co-propagating laser Raman interferometers, while for counter-propagating laser Raman interferometers, both terms contribute to the phase shift.
## 13.2 Are the calculations doable in the laboratory frame ?

In the previous section, the observable phase shift is computed in the galactic reference frame (the DM rest frame) where the perturbation to the rest mass/atomic frequency is proportional to $\cos(\omega_{\phi}t + \Phi)$ (in the dilatonic case, but the argument is the same for the axion). As it was shown in Chapter 7 (but this can be recovered by a simple Lorentz transformation), the field (and by extension the perturbation to the rest mass or atomic frequency) is proportional to $\cos(\omega_{\ell}t_{\ell} - \vec{k}_{\ell} \cdot \vec{x}_{\ell} + \Phi)$ in a lab-centric frame, moving at velocity $\vec{v}_{\text{DM}}$ with respect to the Galactocentric frame (see Eq. (13.22) for their explicit expression). The $\ell$ subscript denotes quantities expressed in the lab-centric frame. The reason why we performed the calculations in the galactocentric frame relies on the fact that the method presented in [13] is valid only for Lagrangians at most quadratic in the position and velocity, which is the case in the galactocentric case but not in the lab-centric one. For this reason, the quantities appearing in Eqs. (13.19), (13.20) and (13.21) (and as we shall see in the next sections, this works also for Eqs. (13.26) and (13.27)) are quantities evaluated in the galactocentric frame. The first goal of this section is to justify why one can safely replace the values of $\omega_{A}^{0}$, $k_{\text{eff}}$ and $L$ appearing in these equations by their lab-centric counterpart. The result from Eqs. (13.19) and (13.20) differs from the one obtained in [17] by the fact that in [17], the mass terms are proportional to the lab-centric initial velocity of the atoms while in the previous section, it corresponds to its galactic counterpart (i.e the galactic velocity, which is multiple orders of magnitude larger).

To strengthen our point, we will explicitly derive the classical equations of motion of the atoms in the lab-centric frame and show that the galactic velocity is indeed expected to appear in the phase shift. Finally, we will explain how a solution can be derived in the lab-centric frame (at first order in $v_{\text{DM}}/c$) using the formalism from [13] and show that this is also consistent with Eqs. (13.19).

First, the two reference frames are related to each other by a Lorentz transformation. This means that the angular frequency $\omega$ and wave vector $\vec{k}$ of the laser beam transform following

$$ \omega_{\ell} = \gamma \left( \omega - \vec{v}_{\text{DM}} \cdot \vec{k} \right) = \omega \left( 1 + O\left( \frac{v_{\text{DM}}}{c} \right) \right), \qquad (13.22a) $$

$$ \vec{k}_{\ell} = \vec{k} + \frac{1}{v_{\text{DM}}^2} (\gamma - 1) (\vec{v}_{\text{DM}} \cdot \vec{k}) \vec{v}_{\text{DM}} - \frac{1}{c^2} \gamma \omega \vec{v}_{\text{DM}} = \vec{k} \left( 1 + O\left( \frac{v_{\text{DM}}}{c} \right) \right), \qquad (13.22b) $$

where $\gamma = (1 - v_{\text{DM}}^2/c^2)^{-1/2}$ and $v_{\text{DM}}/c \sim 10^{-3}$. For this reason, although Eqs. (13.19) and (13.20) are expressed in the galactocentric reference frame, one can safely replace $\vec{k}_{\text{eff}}$, $\omega_{A}^{0}$ by their lab-centric counterpart. This would lead to a correction three orders of magnitude smaller than the leading order term. A similar argument applies for the other quantities such as $L, T$ the interrogation time, etc. This demonstrates that one can safely use lab-centric quantities in Eqs. (13.19) and (13.20) (and this applies to all the equations following them, including the ones derived in the next two sections).

Let us now convince ourselves that a derivation directly performed in the lab-centric frame would also lead to a phase shift whose main term is also proportional to the lab velocity with respect to the galactic reference frame, in agreement with Eqs. (13.19) and (13.20). In Chapter 12, we derived the acceleration felt by an atom A in the laboratory frame
Eq. (12.4b). Therefore, in the lab frame, the acceleration is proportional to the galactic velocity $v_{DM}$, at leading order. In addition, the AI schemes considered in Eqs. (13.19) and (13.20) act as accelerometers that provide, to first order, a phase shift $\Delta\Phi = k_{eff}aT^2$ where $a$ is the local acceleration. For this reason, as an output to the perturbative acceleration form Eq. (12.4b), one expects that the leading order term for $\Delta\Phi$ is $\propto k_{eff}v_{DM}[Q_M^A]_d T^2$, consistent for $\omega T \ll 1$ with the first term from Eqs. (13.19) and (13.20). This reasoning provides an argument showing that the velocity of the laboratory with respect to the galactocentric reference frame also appears when reasoning directly in the lab frame.

Formally, the method from [13] cannot be used to compute the AI phase shift from the lab frame Lagrangian provided by Eq. (8.11b) since it is not at most quadratic in the position. Nevertheless, it is possible to perform an approximate calculation in the lab frame in the case where the de Broglie wavelength of the DM field (which for demonstration we assume is the dilaton field) is much larger than the typical size of the experiment, i.e. if $\vec{k}_{\phi} \cdot \vec{x}_{\ell} \ll 1$ (which is the case for the AI experiments considered in this thesis (see Chapter 17) characterized by $\vec{k}_{\phi} \cdot \vec{x}_{\ell} \lesssim 5 \times 10^{-8} \ll 1^5$). In such a case, an expansion of the Lagrangian Eq. (8.11b) leads to

$$ \mathcal{L}_A \approx -m_A^0 c^2 \left( 1 - \frac{v_{A,\ell}^2}{2c^2} \right) \times \left( 1 - \frac{\sqrt{16\pi G\rho_{DM}}[Q_M^A]_d}{\omega_\phi c} \left( \cos(\omega_\phi t_\ell + \Phi) + \vec{k}_\phi \cdot \vec{x}_\ell \sin(\omega_\phi t_\ell + \Phi) \right) \right). \quad (13.23) $$

From this Lagrangian, one can derive the UFF violating acceleration in the laboratory frame, which reads

$$ \vec{a}_A \approx (\omega_\phi \vec{v}_{A,\ell} - \vec{k}_\phi c^2) \frac{\sqrt{16\pi G\rho_{DM}}[Q_M^A]_d}{\omega_\phi c} \sin(\omega_\phi t_\ell + \Phi) \quad (13.24a) $$

$$ = \frac{(\vec{v}_{A,\ell} + \vec{v}_{DM})}{c} \sqrt{16\pi G\rho_{DM}} [Q_M^A]_d \sin(\omega_\phi t_\ell + \Phi). \quad (13.24b) $$

This Lagrangian is now linear in the position and can be used to compute the phase shift using the method from [13]. Note that there is an important difference compared to the calculation performed in [17] where the second term has implicitly been neglected.

The equations of motion deriving from this Lagrangian consists in the ones from Eq. (12.4b) where one neglects the $\vec{k}_{\ell} \cdot \vec{x}_{\ell}$ within the sine function. The calculation of the phase shift directly in the lab frame follows exactly the equations presented in the previous section. There are mainly two differences in the derived equations. First, the perturbed trajectory and velocity of the atom in the lab frame now read (taking $t_0 = 0$ immediately and $t_{\ell} \to t$)

$$ \vec{v}_{A,\ell}(t) \approx \vec{v}_{\ell,0} - (\vec{v}_{\ell,0} + \vec{v}_{DM}) X_{DM}[Q_M^A]_d (\cos(\omega_\phi t + \Phi) - \cos(\Phi)), \quad (13.25a) $$

$$ \vec{x}_{A,\ell}(t) \approx \vec{x}_{\ell,0} + \vec{v}_{\ell,0}t - \frac{\vec{v}_{\ell,0} + \vec{v}_{DM}}{\omega_{\phi}} X_{DM}[Q_M^A]_d (\sin(\omega_{\phi} t + \Phi) - \sin(\Phi) - \omega_{\phi} t \cos(\Phi)), \quad (13.25b) $$

where $v_{\ell,0}$ is the initial velocity of the atom in the lab frame, which corresponds to the launch velocity. These equations replace Eqs. (13.6). Secondly, one needs to keep the $\vec{k}_{\phi} \cdot \vec{x}_{\ell}$ term in

<footnote n="5">⁵We consider DM frequencies characterized by $m_{\phi}c^2 \le 10^{-13}$ eV and size of experiments $L \le 100$ m).</footnote>
<figure><img src="image_15.png" /></figure>

Figure 13.4: Spacetime diagram of a single photon transition, as proposed in [23] with $n = 2$ (using the convention of [24], see text.) Two lasers are used, one located at $x = 0$, with emission (in red) towards the second laser located at $x = L$ (with blue emission). The frequencies of the lasers are selected carefully to interact only with the wavepacket of interest (see text). The black dots indicate the location of matter-wave interactions.

the Lagrangian when computing the propagation phase using Eq. (13.8). The full calculation following the method from the last section leads to a result that is consistent with Eqs. (13.19) and (13.20) to first order in $v_{DM}/c$, in Bragg and Raman AI respectively, which demonstrates the equivalence between the two frames.

## 13.3 Gradiometers

We will now focus on interferometric setup that involves single photon transitions. The first setup considered is known as a gradiometer, i.e a setup where two atom interferometers are stacked at different altitudes, for the study of e.g gravity gradients, and which has already been studied in [18, 23].

We are interested in the setup initially proposed by [23] and then studied in e.g. [18, 24]. Practically, we consider two ensembles of atoms $A$ (one for each interferometer) located respectively at $x_1$ and $x_2$, all initially in the state $|g, \hbar k\rangle$ and two lasers, one at coordinate $x=0$ with wavevector $\vec{k}_1$ and the other at coordinate $x=L$ with wavevector $\vec{k}_2$. In those single photon interactions configurations, $k_1 \approx k_2 = k = \omega_A^0/c$.

At an initial time $t_0$, the first laser sends a beam which interacts with both atom ensembles and which corresponds to a $\pi/2$ pulse. Then, the second laser beam sends a $\pi$ Doppler shifted laser pulse in order to interact only with the excited state wavepackets, whose motion induces a change in transition frequency, and to convert entirely this wavepacket to the ground state. After this sequence, all wavepackets are in the ground state, however one of them has gained momentum $2\hbar k$. Following the convention of [18, 24], at the end of the sequence, the fast wavepackets have received a Large Momentum Transfer (LMT) photon kick of order $2^6$. If

<footnote n="6">⁶Note that in the original proposal by [23], the LMT has a slightly different definition.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 13.5: Spacetime diagram of a gradiometer, where two interferometers as shown in Fig. 13.4 are stacked at different altitudes. In this figure, trajectories of the atom in its ground state are represented by dashed lines, while full lines represent paths where the atom is in its excited state.</figcaption></figure>

another pair of $\pi$ pulses (with the first one from the bottom laser and the second one from the top laser in the same way as before) is sent to the faster half of the atom, one makes a large momentum transfer (LMT) beam splitter of order 4. More generally, if $m$ pairs of $\pi$ pulses are sent, the order of the LMT beam splitter is $2(m + 1)$ and the faster wavepacket has gained total momentum $2(m + 1)\hbar k$. In other words, the order $n$ of the LMT is defined as $n = 2(m + 1)$.

Later at time $t = T$, a sequence of state inversion similar to the one used in two-photon transitions interferometers is performed, but this time, with three different $\pi$ pulses, the first one coming from the bottom laser, the second one from the top laser and the last one from the bottom laser again. However, in order to slow down the faster wavepacket, $m$ pairs of $\pi$ pulses are added before this sequence, such that it loses $2(m + 1)\hbar k$ momentum. Symmetrically, $m$ other pairs of LMT pulses are added after the state inversion to accelerate the other wavepacket, such that it gains $2(m + 1)\hbar k$ momentum.

Finally at time $t = 2T$, a sequence of pulses opposite to the one sent at the initial $t = 0$ is sent to the wavepackets, i.e $m$ LMT pairs of $\pi$ pulses are sent to the wavepackets before the final $\pi - \pi/2$ pulses used for recombination. This whole sequence is depicted in Fig. 13.4 with $n = 2$ (i.e $m = 0$) and only one interferometer.

In gradiometers, two such interferometers are stacked at different altitudes, separated by a distance $\Delta r$, and the same laser beam is used for the laser-atom interaction in both interferometers, which is depicted in Fig. 13.5. The main advantage of this setup is that laser phase noise is entirely cancelled when measuring the differential phase shift between the two interferometers.

Assuming $\omega_{\phi}L/c, \omega_aL/c \ll 1$ and $nL/c \ll T$ where the LMT kick is of order $n$, and $L$ is the baseline separation between the two lasers, which we assume for simplicity to be
the distance between the two interferometers, the differential phase shift can be computed following the exact same methodology as in Section 13.1.3 and reads

$$|\Delta\Phi_A^{\text{Grad}}|_d \approx \frac{4\sqrt{16\pi G\rho_{\text{DM}}n\omega_A^0\Delta r[Q_\omega^A]}_d}{\omega_\phi c^2} \sin^2\left(\frac{\omega_\phi T}{2}\right) \quad (13.26a)$$

$$|\Delta\Phi_A^{\text{Grad}}|_a \approx \frac{32\pi G \rho_{\text{DM}} E_p^2 n \omega_A^0 \Delta r [Q_\omega^A]_a}{f_a^2 \omega_a^2 c^3} \sin^2(\omega_a T) . \quad (13.26b)$$

We recover the leading order phase shift amplitude already derived in [18, 23].

## 13.4 Single Photon Isotope Differential interferometer

We study another setup employing single photon transitions, that we name SPID in the following, for Single Photon Isotope Differential AI, which is a variation of the interferometric sequence presented in the previous section. We will show that this experimental setup is more sensitive to oscillations in atom rest mass and transition frequency compared to regular gradiometers. This type of setup has already been proposed, without any detail, for ULDM detection in MAGIS-100 [25] (see end of this section for a discussion). The goal of this section is to show the expected signals of such a setup, in order to compare it directly with gradiometers, as the one expected to be used in AION-10 [18].

Contrary to usual gradiometers which only use one single species of atom, we consider two different atom isotopes, each of them undergoing individually the interferometric path described earlier in the previous section. The setup is presented in Fig. 13.6. The two interferometers overlap at the same elevation, so the experiment will test the universality of free fall between the two isotopes. This setup employs single photon transition (meaning we will consider only optical transitions) and measures the differential acceleration between two isotopes, so we will refer to this setup as Single Photon transition Isotope Differential (SPID).

Typically, an optical transition has a frequency of the order of ~10¹⁴ Hz while the typical frequency shift between two isotopes is of order ~10⁹ Hz (i.e. typically 5-6 orders of magnitude smaller than the nominal transition frequency, see [26] for Sr). Then, a unique laser source can be used in this setup and separated in two different beams: one which is directly used in the laser-atom interactions inside the interferometer of the first isotope; and the other one whose frequency is shifted, e.g. using an electro-optic modulator (EOM), in order to interact with the second isotope.

The calculation of the differential phase shift for such a setup can be computed following the exact same methodology as in Section 13.1.3. If we assume the optical transition frequencies to be close, i.e. $\omega_A^0 \approx \omega_B^0 \equiv \omega_0$ and the initial velocity ($\propto v_{\text{DM}}$) to be much larger than the velocity kick $\hbar k_{\text{eff}}/m^0$, the differential phase shift between the two interferometers reads

$$|\Delta\Phi_{AB}^{\text{SPID}}|_d \approx \frac{4\sqrt{16\pi G\rho_{\text{DM}}v_{\text{DM}}n\omega_0}}{\omega_\phi^2 c^2} |([Q_M^A]_d - [Q_M^B]_d)\hat{e}_v \cdot \hat{e}_{\text{kick}}| \sin^2\left(\frac{\omega_\phi T}{2}\right) \quad (13.27a)$$

$$|\Delta\Phi_{AB}^{\text{SPID}}|_a \approx \frac{16\pi G \rho_{\text{DM}} E_p^2 v_{\text{DM}} n \omega_0}{f_a^2 \omega_a^3 c^3} |([Q_M^A]_a - [Q_M^B]_a)\hat{e}_v \cdot \hat{e}_{\text{kick}}| \sin^2(\omega_a T), \quad (13.27b)$$
<figure><img src="image_17.png" /></figure>

Figure 13.6: Simplified setup for the *SPID* experiment. For the sake of simplicity, we assume two different isotopes with respective transition frequencies $f_1 = f_{tr} = 10^{14}$ Hz, $f_2 = f_{tr} + f_{EOM} = (1 + 10^{-5}) \times 10^{14}$ Hz. The laser is locked on the transition frequency of isotope 1 and split in three: the first two outputs are used for the isotope 1 AI (in red). The third output enters an electro-optic modulator (EOM) to shift its frequency by $f_{EOM}$ to account for the isotope shift in order to be used with the second isotope. As for the other experimental AI schemes, we do not show the various acousto-optic modulators (AOM) that shift the frequency of one of the laser input of each interferometer, to account for the Doppler shift of the freely falling atoms.

where we kept only lowest order terms in $\omega L/c$. In this result, we did not take into account the effect of an oscillating EOM frequency through $Q_{\omega}^{\text{EOM}} \neq 0$. Nonetheless, its effect is suppressed by a factor $\omega L/v_{\text{DM}} \times \Delta\omega^0/\omega^0 \le 10^{-7}$ compared to the leading term, where $\Delta\omega^0$ is the isotope shift. In such a case (and similarly for $L=0$), we recover the leading order Bragg phase shift derived in Eq. (13.21) with $n\omega_0/c = k_{\text{eff}}$, as expected.

The advantage of this setup compared to usual gradiometers can be immediately visualized by comparing Eq. (13.27) and Eq. (13.26): the *SPID* setup does not suffer from a small factor $\omega\Delta r/v_{\text{DM}}$. More precisely, the ratio of amplitude of signals between this variation and usual gradiometers is roughly $v_{\text{DM}}/(\omega\Delta r) \times (Q_M^A - Q_M^B)/Q_\omega^A \sim (1 \text{ rad}\cdot\text{s}^{-1}/\omega)$ for the following values: $\Delta r \sim 5 \text{ m}, v_{\text{DM}} \sim 10^{-3} \text{ c}, \hat{e}_v \cdot \hat{e}_{\text{kick}} \sim \mathcal{O}(1)$ and mass and frequency charges presented in Tables 8.1 and 8.3.). This implies that, at low angular frequency ($\omega < 1 \text{ rad/s}$), the signal of the single photon transition isotope differential AI will be larger than the one in a gradiometer. The simple reason for this difference in signal amplitudes is that the gradiometer leading order phase shift in Eq. (13.26) is proportional to the frequency charge $Q_\omega$, or in other words, it is proportional to the time in which the wavepackets are in their excited state. As can be noticed from Fig. 13.4, this happens for a limited time, of the order of $nL/c$. Conversely, the signal amplitude of the *SPID* variation is proportional to mass charges $Q_M$ whose effect is imprinted in the phase shift, whatever the internal state, i.e for a time $\sim T$.

As it was pointed out in the beginning of this section, *MAGIS-100* [25] will operate a similar mode as *SPID* for the search of ULDM and we want to derive the expected sensitivity of such a large scale experiment to axion and dilaton signals. More specifically, *MAGIS-100* will run a Bragg interferometer (i.e with two photons transitions instead of single-photon transition,
as in *SPID*). However, as mentioned in the text after Eq. (13.27), the leading order signal of both setup is the same, provided that the beams are locked on the same optical transition in both cases. The main difference between the two setup is that for *SPID*, the wavepackets spend some time in their excited state (which they do not in Bragg configuration), but as shown in Eq. (13.27), the additional phase shift is next-to-leading order, and thus negligible. In addition, as we shall see in Chapter 18, the phase noise levels will be equivalent in both setups. Therefore, the optimal choice between the two setup is a matter of practicality. In particular, we will be deriving the sensitivity of such setup considering optical transitions summarized in Table 8.2. In these choices of isotopes and transitions, we ignore the small lifetime of the various excited states (of hundreds of *ns* to hundreds of *µs*), which could limit the number of atoms detected at the end of the sequence due to spontaneous emission<fnref n="7" />

<footnote n="7">⁷This limitation could be overcome for optical transitions due to the large Rabi frequency, implying a very short π pulse duration [27].</footnote>
[1] Jean-Philippe Uzan. "Varying Constants, Gravitation and Cosmology". en. In: *Living Reviews in Relativity* 14.1 (Dec. 2011), p. 2. ISSN: 2367-3613, 1433-8351. DOI: 10.12942/lrr-2011-2.

[2] Thibault Damour. "Theoretical aspects of the equivalence principle". In: *Classical and Quantum Gravity* 29.18 (Aug. 2012), p. 184001. DOI: 10.1088/0264-9381/29/18/184001. URL: https://dx.doi.org/10.1088/0264-9381/29/18/184001.

[3] T. Damour and A.M. Polyakov. "The string dilation and a least coupling principle". In: *Nuclear Physics B* 423.2 (1994), pp. 532-558. ISSN: 0550-3213. DOI: https://doi.org/10.1016/0550-3213(94)90143-0. URL: https://www.sciencedirect.com/science/article/pii/0550321394901430.

[4] Aurélien Hees et al. "Violation of the equivalence principle from light scalar dark matter". In: *Phys. Rev. D* 98 (6 Sept. 2018), p. 064051. DOI: 10.1103/PhysRevD.98.064051. URL: https://link.aps.org/doi/10.1103/PhysRevD.98.064051.

[5] Gianfranco Bertone and Dan Hooper. "History of dark matter". In: *Rev. Mod. Phys.* 90 (4 Oct. 2018), p. 045002. DOI: 10.1103/RevModPhys.90.045002. URL: https://link.aps.org/doi/10.1103/RevModPhys.90.045002.

[6] Lam Hui et al. "Ultralight scalars as cosmological dark matter". In: *Phys. Rev. D* 95 (4 Feb. 2017), p. 043541. DOI: 10.1103/PhysRevD.95.043541. URL: https://link.aps.org/doi/10.1103/PhysRevD.95.043541.

[7] Pierre Touboul et al. "MICROSCOPE Mission: First Results of a Space Test of the Equivalence Principle". In: *Phys. Rev. Lett.* 119 (23 Dec. 2017), p. 231101. DOI: 10.1103/PhysRevLett.119.231101. URL: https://link.aps.org/doi/10.1103/PhysRevLett.119.231101.

[8] Joel Bergé. "MICROSCOPE's view at gravitation". In: *Reports on Progress in Physics* 86.6 (May 2023), p. 066901. DOI: 10.1088/1361-6633/acd203. URL: https://dx.doi.org/10.1088/1361-6633/acd203.

[9] Pierre Touboul et al. "MICROSCOPE Mission: Final Results of the Test of the Equivalence Principle". In: *Phys. Rev. Lett.* 129 (12 Sept. 2022), p. 121102. DOI: 10.1103/PhysRevLett.129.121102. URL: https://link.aps.org/doi/10.1103/PhysRevLett.129.121102.

[10] S. Schlamminger et al. "Test of the Equivalence Principle Using a Rotating Torsion Balance". In: *Phys. Rev. Lett.* 100 (4 Jan. 2008), p. 041101. DOI: 10.1103/PhysRevLett.100.041101. URL: https://link.aps.org/doi/10.1103/PhysRevLett.100.041101.

[11] JAMES G. WILLIAMS, SLAVA G. TURYSHEV, and DALE H. BOGGS. "LUNAR LASER RANGING TESTS OF THE EQUIVALENCE PRINCIPLE WITH THE EARTH AND MOON". In: *International Journal of Modern Physics D* 18.07 (July 2009), pp. 1129-1175. ISSN: 1793-6594. DOI: 10.1142/s021827180901500X. URL: http://dx.doi.org/10.1142/S021827180901500X.

[12] Jordan Gué, Aurélien Hees, and Peter Wolf. "Violation of the equivalence principle induced by oscillating rest mass and transition frequency, and its detection in atom interferometers". In: *Phys. Rev. D* 110.3 (2024), p. 035005. DOI: 10.1103/PhysRevD.110.035005. arXiv: 2401.14742 [hep-ph].
[13] Pippa Storey and Claude Cohen-Tannoudji. "The Feynman path integral approach to atomic interferometry. A tutorial". In: *Journal de Physique II* 4.11 (1994), pp. 1999-2027. DOI: 10.1051/jp2:1994103. URL: https://hal.archives-ouvertes.fr/jpa-00248106.

[14] Savas Dimopoulos et al. "Testing General Relativity with Atom Interferometry". In: *Phys. Rev. Lett.* 98 (11 Mar. 2007), p. 111102. DOI: 10.1103/PhysRevLett.98.111102. URL: https://link.aps.org/doi/10.1103/PhysRevLett.98.111102.

[15] Peter Wolf and Philippe Tourrenc. "Gravimetry using atom interferometers: Some systematic effects". In: *Physics Letters A* 251.4 (1999), pp. 241-246. ISSN: 0375-9601. DOI: https://doi.org/10.1016/S0375-9601(98)00881-0. URL: https://www.sciencedirect.com/science/article/pii/S0375960198008810.

[16] Peter Wolf and Christian J. Borde. *Recoil Effects in Microwave Ramsey Spectroscopy*. 2004. arXiv: quant-ph/0403194 [quant-ph].

[17] Andrew A. Geraci and Andrei Derevianko. "Sensitivity of Atom Interferometry to Ultralight Scalar Field Dark Matter". In: *Phys. Rev. Lett.* 117 (26 Dec. 2016), p. 261301. DOI: 10.1103/PhysRevLett.117.261301. URL: https://link.aps.org/doi/10.1103/PhysRevLett.117.261301.

[18] Leonardo Badurina, Diego Blas, and Christopher McCabe. "Refined ultralight scalar dark matter searches with compact atom gradiometers". In: *Phys. Rev. D* 105 (2 Jan. 2022), p. 023006. DOI: 10.1103/PhysRevD.105.023006. URL: https://link.aps.org/doi/10.1103/PhysRevD.105.023006.

[19] Holger Ahlers et al. *STE-QUEST: Space Time Explorer and QUantum Equivalence principle Space Test*. 2022. DOI: 10.48550/ARXIV.2211.15412. URL: https://arxiv.org/abs/2211.15412.

[20] A. Gauguet et al. "Off-resonant Raman transition impact in an atom interferometer". In: *Physical Review A* 78.4 (2008). DOI: 10.1103/physreva.78.043615. URL: https://doi.org/10.1103%2Fphysreva.78.043615.

[21] P. Gillot et al. "Limits to the symmetry of a Mach-Zehnder-type atom interferometer". In: *Phys. Rev. A* 93 (1 Jan. 2016), p. 013609. DOI: 10.1103/PhysRevA.93.013609. URL: https://link.aps.org/doi/10.1103/PhysRevA.93.013609.

[22] E. Giese et al. "Double Bragg diffraction: A tool for atom optics". In: *Physical Review A* 88.5 (Nov. 2013). DOI: 10.1103/physreva.88.053608. URL: https://doi.org/10.1103%2Fphysreva.88.053608.

[23] Peter W. Graham et al. "New Method for Gravitational Wave Detection with Atomic Sensors". In: *Phys. Rev. Lett.* 110 (17 Apr. 2013), p. 171102. DOI: 10.1103/PhysRevLett.110.171102. URL: https://link.aps.org/doi/10.1103/PhysRevLett.110.171102.

[24] Asimina Arvanitaki et al. "Search for light scalar dark matter with atomic gravitational wave detectors". In: *Phys. Rev. D* 97 (7 Apr. 2018), p. 075020. DOI: 10.1103/PhysRevD.97.075020. URL: https://link.aps.org/doi/10.1103/PhysRevD.97.075020.

[25] Mahiro Abe et al. "Matter-wave Atomic Gradiometer Interferometric Sensor (MAGIS-100)". In: *Quantum Science and Technology* 6.4 (July 2021), p. 044003. DOI: 10.1088/2058-9565/abf719. URL: https://dx.doi.org/10.1088/2058-9565/abf719.
[26] Tetsushi Takano, Ray Mizushima, and Hidetoshi Katori. "Precise determination of the isotope shift of 88Sr-87Sr optical lattice clock by sharing perturbations". In: *Applied Physics Express* 10.7 (June 2017), p. 072801. DOI: 10.7567/APEX.10.072801. URL: https://dx.doi.org/10.7567/APEX.10.072801.

[27] Jan Rudolph et al. "Large Momentum Transfer Clock Atom Interferometry on the 689 nm Intercombination Line of Strontium". In: *Phys. Rev. Lett.* 124 (8 2020), p. 083604. DOI: 10.1103/PhysRevLett.124.083604. URL: https://link.aps.org/doi/10.1103/PhysRevLett.124.083604.
# Part VI

## Effects on gravitational waves detectors
# Chapter 14

## Gravitational waves and LISA

### 14.1 Gravitational waves physics basics

As quickly mentioned in the introduction of this thesis, one of the main prediction of GR is the production of GW by massive objects. More precisely, only systems with time dependent quadrupolar moment (or higher order moment) produce such waves.

GW are solutions of the linearized Einstein equations. Assuming a static Minkowski background $\eta_{\mu\nu}$ perturbed by $h_{\mu\nu} \ll \eta_{\mu\nu}$ such that the spacetime metric is $g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu}$, around a system energy momentum tensor $T_{\mu\nu}$, one can show that [1]

$$ \Box \bar{h}_{\mu\nu} = \frac{-16\pi G}{c^4} T_{\mu\nu}, \qquad (14.1a) $$

where $\Box \equiv \eta^{\alpha\beta}\partial_{\alpha}\partial_{\beta}$ and

$$ \bar{h}_{\mu\nu} = h_{\mu\nu} - \frac{1}{2}\eta_{\mu\nu}h^{\alpha}_{\alpha}. \qquad (14.1b) $$

In presence of matter-energy, Eq. (14.1a) describes the emission of the GW sourced by $T_{\mu\nu}$, while in vacuum where $T_{\mu\nu} = 0$, it describes the propagation of the GW. The energy-momentum tensor being symmetric, it consists of 10 independent components, therefore Eq. (14.1a) is a set of 10 independent equations. Imposing Lorenz gauge $\partial_{\mu}\bar{h}^{\mu} = 0$ and the so-called TT (traceless-transverse) gauge $\bar{h}_{\mu}^{\mu} = 0, \partial^{i}\bar{h}_{ij} = 0$ leaves only two independent radiative degrees of freedom, which are the two GW polarizations. Therefore, in the TT gauge, we have $h_{\mu\nu}^{TT} = \bar{h}_{\mu\nu}^{TT}$. Assuming a GW wave with frequency $f = \omega/2\pi$ propagating in the z-direction, we have $h_{\mu\nu}^{TT} = \Re[A_{\mu\nu}^{TT}e^{-i(\omega t - kz)}]$ with $A_{11}^{TT} = A_+$ and $A_{12}^{TT} = A_-$, and the amplitude of the wave is

$$ \begin{aligned} A_{\mu\nu}^{TT} &= \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & A_{+} & A_{\times} & 0 \\ 0 & A_{\times} & -A_{+} & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} = A_{+} \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} + A_{\times} \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \\ &\equiv A_{+} \epsilon_{\mu\nu}^{+} + A_{\times} \epsilon_{\mu\nu}^{\times} \end{aligned} \qquad (14.2) $$

The amplitudes $A_+, A_\times$ can be found by solving Eq. (14.1a). Using Green's functions and the conservation of the energy-momentum tensor $\partial_\mu T^{\mu\nu} = 0$, one can show that the spatial
solutions of Eq. (14.1a) is the perturbed spacetime at coordinates $(t, \vec{x})$ induced by a GW emitted by a source located at spatial coordinates $\vec{y}$, and has the form [1]

$$h_{ij}^{\text{TT}}(t, \vec{x}) = \frac{2G}{c^4 R} \mathcal{P}_{ijab} \frac{\partial^2}{\partial t^2} \int d^3\vec{y} \left( y_a y_b - \frac{1}{3} \delta_{ab} |\vec{y}|^2 \right) T_{00} \left( t - \frac{|\vec{x} - \vec{y}|}{c}, \vec{y} \right), \quad (14.3a)$$

where the argument of the integral is the quadrupolar moment tensor of the source and where we introduced the projection operators

$$\mathcal{P}_{ijab} = \mathcal{P}_{ia}\mathcal{P}_{jb} - \frac{1}{2}\mathcal{P}_{ij}\mathcal{P}_{ab} \quad (14.3b)$$

$$\mathcal{P}_{ij} = \delta_i \delta_j - \hat{k}_i \hat{k}_j, \quad (14.3c)$$

that ensures that $h_{ij}^{\text{TT}}$ is traceless and transverse to the propagation vector $\hat{k}$. This result is derived assuming that the distance $|\vec{x} - \vec{y}| \gg |\vec{y}|$, the size of the source.

The effect of the GW is the modification of physical distances between test masses in the plane of polarization. Considering the same GW as above, assuming two freely falling bodies located at $z = 0$, separated on the x-axis by a coordinate distance $L_0$. When the GW passes by, this coordinate distance is unchanged, because the coordinates move with the GW [1]. The relevant and physical quantity is the proper distance, as the distance measured by a light signal and a clock, which is shown to change as [1]

$$L(t) = \int_{0}^{L_0} dx \sqrt{g_{xx}} \approx L_0 \left( 1 + \frac{1}{2} h_{xx}^{\text{TT}}(t, z=0) \right). \quad (14.4)$$

More generally, if we assume a ring of particles in the x-y plane, the + polarization of the GW makes the ring oscillating in the x-y directions while the × polarization acts on the ring along the x' - y' directions which are rotated by 45° with respect to x-y coordinates (i.e $\hat{e}_{x'} = \sqrt{2}(\hat{e}_x + \hat{e}_y)/2$ and $\hat{e}_{y'} = \sqrt{2}(-\hat{e}_x + \hat{e}_y)/2$). This is the reason for the +, × polarization names, as they respectively distort a circular ring of particles to make a + and × shape.

Since the GW strain (i.e the relative displacement $(L(t) - L_0)/L_0 \equiv \delta L/L_0$ in Eq. (14.4)) is expected to be extremely small, of the order of $10^{-21}$ or less, one must use the most precise method to measure distances, i.e optical interferometers. For example, the Earth-based GW detectors LIGO/VIRGO are operating a Michelson optical interferometer in a L-shape in order to maximize the effect of the GW, i.e stretch one arm and compress the other. In such a case, a photon of wavelength $\lambda$ that travels back and forth on the arm whose distance has changed acquires a phase shift $\Delta\phi$. Note that even if light is stretched by the passage of the GW (i.e its frequency changes), one can still use it as a "ruler" to detect any change in physical distance between two freely falling objects. This is because light is not used as a ruler but as a clock, i.e one does not measure the total number of crests that fits in a given distance (e.g. between reflecting mirror and beam splitter in an usual GW detector) but rather measures the number of crests of light that arrives at the beam splitter per unit of time. As the GW stretches (or compresses light), the crests will arrive at the beam splitter with a time delay compared to situations with no GW, and this creates an observable phase shift [2].
<figure><img src="image_18.png" /><figcaption>Figure 14.1: LISA constellation heliocentric orbit (from [4]).</figcaption></figure>

## 14.2 LISA, the first european space-based gravitational waves detector

After the success of the ground-based GW detectors such as LIGO and VIRGO, one of the next important mission aiming at detecting GW is LISA (Laser Interferometer Space Antenna), see [3] for a complete review of the mission and its scientific goals. In this manuscript, we will only describe the important features of the mission that will be needed for its understanding. As pictured in Fig. 14.1, LISA will consist of a constellation of three spacecrafts, in the following noted 1, 2, 3, which will follow Earth on an heliocentric orbit. In practice, LISA will be able to detect GW of frequencies in the $10^{-4} - 1$ Hz range, which is complementary to LIGO/VIRGO sensitivity bands (around $10 - 10^4$ Hz). However, the mean spatial separation between the spacecrafts being $L = 2.5 \times 10^9$ m, at frequencies $f > c/L \sim 0.1$ Hz, the sensitivity decreases due to the arm length penalty [5], i.e the GW wavelength is shorter than the arm length, which induces averaging of the signal. All spacecrafts will continuously produce two local laser signals through two local oscillators and exchange it with the other spacecrafts, and GW signals will be detected by interferometric measurement between those optical signals. Indeed, at the passage of a GW, the distance between spacecrafts changes and therefore, the observed phase of the laser signal that travelled in vacuum changes as well. More precisely, inside one given spacecraft (e.g. 1), interferometric measurements will be made between the local oscillator of 1 and the laser received from spacecraft 2 using a phasemeter. A very specific optical bench inside the spacecrafts allows to recombine the two optical signals and make the interference. A more precise description of such bench will be done in Chapter 16. Overall, in order to make interferometric measurements between all the optical signals, each spacecraft contains two optical benches. In addition, inside each of the spacecrafts is located a freely-falling test mass surrounded by metallic shields in order to protect from any exterior perturbations, e.g. radiation pressure from the Sun. This way, the test mass is insensitive to non-gravitational forces which allows it to stay on its geodesic. By measuring the distance between the test mass and the spacecraft around it, one knows when the spacecraft has left its geodesic, due to perturbations, and can replace the spacecraft on its geodesy using mini-rockets on board. Combining all these measurements in all spacecrafts, one can reconstruct the distance between the test masses.

Finally, let us discuss sources of noise in the apparatus, in particular laser noise. It is admitted that the laser noise amplitude spectral density (ASD) is $\sqrt{S_f} \sim 30$ Hz/$\sqrt{\text{Hz}}$ [6], which in terms of the nominal frequency of the laser $\nu_0 \sim 2.82 \times 10^{14}$ Hz [6], corresponds to
$\sqrt{S_f}/v_0 \sim 10^{-13}$ Hz$^{-1/2}$. In the same unit, the GW signal ASD is $\sim 10^{-21}$ Hz$^{-1/2}$, i.e at this level, the laser noise completely dominates the signal. Fortunately, there exists a method to effectively lower drastically the laser noise by combining signals from the different spacecrafts in such a way that the resulting observable is independent of laser noise, but still contains the GW signal. This is called *Time Delay Interferometry*, or TDI and were first introduced in [7–9]. The principle relies on the fact that each laser noise is measured several times.

It can be shown mathematically (for details, we refer to [6]) that some time retarded signals combinations between the three different arms completely reduces the laser noise to zero, thus the name of the method. Here, we will only present some of the TDI combinations of interest.

First, we will label the optical benches with two indices *ij*, where the spacecraft *i* receives light signal (host optical bench) from emitting spacecraft *j*. Defining the laser phase of the laser on spacecraft *i* by $\Phi_i$, the beatnote phase between $\Phi_i$ and $\Phi_j$ at time of reception $t_r$ is

$$ \Phi_{ij}(t_r) = \Phi_j \left( t_r - \frac{L_{ij}(t_r)}{c} \right) - \Phi_i(t_r), \qquad (14.5a) $$

where $L_{ij}$ is the arm length between the two spacecrafts and where we assume that the measurements inside different spacecrafts are synchronized. We define the relative phase fluctuation

$$ \delta\Phi_{ij}(t_r) = -\frac{\omega_j \delta L_{ij}}{c} \equiv \frac{\omega_j L_{ij} \gamma_{ij}}{c}, \qquad (14.5b) $$

where $\delta L_{ij}$ is the arm length fluctuation induced by the GW, $\omega_j$ is the local frequency of spacecraft *j* and where the interferometric measurement is expressed in terms of relative arm length fluctuation induced by the GW

$$ \gamma_{ij}(t_r) = -\frac{\delta L_{ij}}{L_{ij}}. \qquad (14.5c) $$

One can also define relative frequency fluctuation signals (i.e Doppler shifts) as

$$ y_{ij}(t_r) = \frac{1}{\omega_j} \frac{d\delta\Phi_{ij}(t_r)}{dt_r} = -\frac{1}{c} \frac{d\delta L_{ij}}{dt}. \qquad (14.6) $$

Using the sign convention of [10], we define the first (and second) generation Michelson $X_1$ ($X_2$) combination as [11]

$$ X_1(t) = y_{13} + D_{13}y_{31} + D_{131}y_{12} + D_{1312}y_{21} - (y_{12} + D_{12}y_{21} + D_{121}y_{13} + D_{1213}y_{31}) \qquad (14.7a) $$

$$ X_2(t) = (1 - D_{12131}) (y_{13} + D_{13}y_{31} + D_{131}y_{12} + D_{1312}y_{21}) - (1 - D_{13121}) (y_{12} + D_{12}y_{21} + D_{121}y_{13} + D_{1213}y_{31}), \qquad (14.7b) $$

where we used the compacted writing of TDI where $D_{ij}y_{jk}(t) = y_{jk}(t - L_{ij}/c)$ and $D_{ijk}y_{kl}(t) = y_{kl}(t - L_{ij}/c - L_{jk}(t - L_{ij}/c))$. It can be shown geometrically that the circulation of the laser beams in that case is very similar to an usual optical Michelson interferometers, thus the name of the combination. The first generation combinations are not sensitive to the laser noise in the constant armlength approximation while the second generation takes into account
effects of rotation of the detector (Sagnac effect) and of the time variation of armlengths. The residuals of laser noise in the second generation TDI is proportional to the difference in spacecraft velocities squared [12]. In the constant armlength approximation and considering that the Doppler shifts are harmonic functions of the frequency $\omega$, the second generation TDI X combination Eq. (14.7a) can be rewritten as

$$X_2 = \left(1 - e^{-\frac{4i\omega L}{c}}\right) \left(y_{13} + e^{-\frac{i\omega L}{c}} y_{31} + e^{-\frac{2i\omega L}{c}} y_{12} + e^{-\frac{3i\omega L}{c}} y_{21}\right) - \left(1 - e^{-\frac{4i\omega L}{c}}\right) \left(y_{12} + e^{-\frac{i\omega L}{c}} y_{21} + e^{-\frac{2i\omega L}{c}} y_{13} + e^{-\frac{3i\omega L}{c}} y_{31}\right) \quad (14.8a)$$

$$= -4 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \left(e^{-\frac{3i\omega L}{c}} (y_{13} - y_{12} + e^{-\frac{i\omega L}{c}} (y_{31} - y_{21}))\right), \quad (14.8b)$$

where we used the fact that one delay operator corresponds to a factor $\exp(-i\omega L/c)$. The other two Michelson Y, Z combinations can be obtained from Eq. (14.7a) by a cyclic permutation of the indices $1 \to 2 \to 3 \to 1$. Other important combinations are the A, E, T combinations which are linear combinations of X, Y, Z [13]

$$A(t) = \frac{1}{\sqrt{2}} (Z(t) - X(t)) \quad (14.9a)$$

$$E(t) = \frac{1}{\sqrt{6}} (X(t) - 2Y(t) + Z(t)) \quad (14.9b)$$

$$T(t) = \frac{1}{\sqrt{3}} (X(t) + Y(t) + Z(t)), \quad (14.9c)$$

which we will be using for the analysis. These combinations are useful because their cross noise power spectral density are vanishing in the equal armlength approximation [11]. For convenience, we also introduce the Sagnac combination $\alpha$ defined as [14]

$$\alpha_1(t) = -(y_{13} + D_{13}y_{32} + D_{132}y_{21} - (y_{12} + D_{12}y_{23} + D_{123}y_{31})) \quad (14.10a)$$

$$\alpha_2(t) = (1 - D_{1231}) \alpha_1(t), \quad (14.10b)$$

respectively for the first and second generations.

## 14.3 LISA's expected performance

As mentioned in the previous section, the Earth-based GW detectors LIGO/VIRGO are sensitive to GW in the ~ 10 - 10^4 Hz frequency range. Out of the ~ 90 various signals detected by the collaboration, all of them are transients signals [15], which are GW signals that last for a few seconds, and which correspond to the last periods of inspirals between the two bodies and their merge into a new body.

On the contrary, LISA is expected to be sensitive to various different signals [6]. In particular, it will be able to detect the inspiral phase of binary systems much before the merge. In such cases, the two bodies rotate around each other at an almost constant frequency, up to corrections (see below), for a time which exceeds largely the time of observation of the mission. Therefore, GW emitted by such sources are quasi-monochromatic. As we shall see
in the following, this is an important point as DM induces also a monochromatic oscillation on the test masses.

In terms of noise, we will be interested in the single-link optical metrology noise and single test mass acceleration noise which respectively read [5, 16]

$$S_{\text{oms}} = (1.5 \times 10^{-11} \text{ m})^2 \left(1 + \left(\frac{2 \text{ mHz}}{f}\right)^4\right) \text{Hz}^{-1} \qquad (14.11a)$$

$$S_{\text{acc}} = (3 \times 10^{-15} \text{ m.s}^{-2})^2 \left(1 + \left(\frac{0.4 \text{ mHz}}{f}\right)^2\right) \left(1 + \left(\frac{f}{8 \text{ mHz}}\right)^4\right) \text{Hz}^{-1}, \qquad (14.11b)$$

corresponding to the instrument performance requirements [4].
# Chapter 15

## Search for oscillations of rest mass in LISA

Many of the various couplings that we have seen so far violate the UFF, by inducing an oscillating acceleration on atoms that is element-dependent. Even though the different test masses used in LISA have the same composition, those oscillations could still be seen in gravitational waves detectors as a Doppler effect when photons are exchanged between one spacecraft and another.

### 15.1 Signatures of dark matter and gravitational waves on LISA arms

In this section, we derive the Doppler on photons exchanged on one arm of LISA induced by both ULDM fields and GW<fnref n="1" />. We make the calculations in the barycentric reference frame, where the spacecrafts positions are defined [6]. In this frame, the directions of DM wind and GW sources are defined with the ecliptic latitude $\beta$ and longitude $\lambda$.

#### 15.1.1 Oscillating rest mass

**One-arm Doppler shift through oscillating rest mass**

In this section, we derive the Doppler shift on the LISA interferometric measurements induced by the various ULDM candidates.

We start by the effects of scalar fields, i.e the UFF violating acceleration induced by the oscillation of the rest mass of a test mass A. In this part, we first consider only interactions between SM and a pure dilatonic field $\phi$, then the Doppler shift will be generalized to the axion field as well. From Eq. (8.12), the acceleration of a test mass A in a generic frame is

$$ \vec{a}_A(t, \vec{x}) = \left[ \omega_{\phi} \vec{v}_A - \vec{k}_{\phi} c^2 \right] \frac{\sqrt{16\pi G \rho_{\text{DM}}}}{\omega_{\phi} c} [Q_M^A]_d \sin(\omega_{\phi} t - \vec{k}_{\phi} \cdot \vec{x} + \Phi), \quad (15.1) $$

<footnote n="1">¹As we shall see in the following, we discuss the Doppler induced by the dynamical motion of the test masses induced by ULDM. As it is pointed out in [17], a change in the length of optical elements (like beam splitters) arises also through the variation of fundamental constants of Nature. While this effect is dominant in ground-based GW detectors (and was used for search for ULDM signals in LIGO, see e.g [18]), it is sub-dominant in space-based GW detectors.</footnote>
where $\vec{v}_A$ is the velocity of the test mass, $\omega_\phi$, $\vec{k}_\phi$ are respectively the Compton frequency and the wavenumber of the DM field, and $[Q_M^A]_d$ is the mass charge of the test mass, i.e the amplitude of oscillation of the rest mass of A. In the barycentric frame, the velocity of the test mass is roughly the same as the Earth velocity around the Sun, i.e $v_A \sim 3 \times 10^4$ m/s, which is one order of magnitude smaller than the galactic velocity, therefore the second term $\propto -\vec{k}_\phi c^2 = \omega_\phi v_{DM}$ dominates. One can treat the position $\vec{x}(t)$ in Eq. (15.1) as the unperturbed position $\vec{x}_0(t)$ at leading order and write it as $\vec{x}(t) = \vec{x}_0(t) + \delta\vec{x}_A(t)$, where the latter is $O([Q_M^A])$. Considering that $\vec{x}(t) = \vec{x}_{AU} \cos(\omega_E t)$ at leading order, where $|\vec{x}_{AU}| \sim 1.5 \times 10^{11}$m is one astronomical unit distance and $\omega_E$ is the Earth rotation frequency around the Sun, one can integrate twice this expression, which leads to the oscillation of the position of the test mass<fnref n="2" />

$$ \delta\vec{x}_A(t, t_0, \vec{x}) = \frac{\vec{k}_\phi c \sqrt{16\pi G \rho_{DM}}}{\omega_\phi^3} [Q_M^A]_d \left[ \sin(\omega_\phi t - \vec{k}_\phi \cdot \vec{x}_{AU} \cos(\omega_E t) + \Phi) - \sin(\omega_\phi t_0 - \vec{k}_\phi \cdot \vec{x}_{AU} \cos(\omega_E t_0) + \Phi) - \omega_\phi(t - t_0) \cos(\omega_\phi t_0 - \vec{k}_\phi \cdot \vec{x}_{AU} \cos(\omega_E t_0) + \Phi) \right]. \quad (15.2) $$

The three last terms are in general not taken into account in various analysis. Indeed, two of them have no time dependence. However, note that one of them $\propto \omega_\phi t$ induces a temporal drift of the test mass inside the spacecraft, which is due to its initial inertia. In the following, we will drop those terms and only work with the time oscillating signal, because they are outside the LISA band. In addition, we will not keep the explicit expression of $\vec{x}(t)$, but we will come back to it later. The slow oscillation at $\omega_E$ will produce side bands in the Fourier signal at frequencies $\omega_\phi \pm \omega_E$. However, note that in the following, we will assume a time of integration of one year $T_{obs} = 2\pi/\omega_E$, such that the size of one Fourier bin is exactly $f_E = \omega_E/2\pi$. Therefore, these side bands might not be visible but will most likely enlarge the signal at $\omega_\phi$.

Assuming that spacecraft e (at position $\vec{x}_e$) sends a signal to spacecraft r (at position $\vec{x}_r$, the one-way Doppler shift of light associated to this rest mass oscillation is obtained by the projection onto the corresponding LISA arm of the time derivative of the variation of the position $\delta\vec{x}$ i.e [19]

$$ y_{re} = -\frac{1}{c} \left( \hat{n}_{re} \cdot \frac{d\delta\vec{x}}{dt} \right), \quad (15.3a) $$

where $\hat{n}_{re}$ is a unit vector pointing from e to r. We can now express the Doppler shift induced by the dilaton-SM interaction Eq. (15.2), but also considering an axion-gluon interaction which induces an acceleration Eq. (8.40). We respectively note $y_{re}^d$ and $y_{re}^a$ these one-way Doppler shifts which read

$$ y_{re}^d = (\hat{n}_{re} \cdot \hat{e}_v) \frac{\sqrt{16\pi G \rho_{DM} v_{DM} [Q_M]_d}}{\omega_\phi c^2} \mathfrak{R} \left( e^{i(\omega_\phi t - \vec{k}_\phi \cdot \vec{x}_r + \Phi)} - e^{i(\omega_\phi t - \frac{\omega_\phi L}{c} - \vec{k}_\phi \cdot \vec{x}_e + \Phi)} \right), \quad (15.3b) $$

$$ y_{re}^a = (\hat{n}_{re} \cdot \hat{e}_v) \frac{8\pi G \rho_{DM} v_{DM} E_p^2 [Q_M]_a}{f_a^2 \omega_a^2 c^3} \mathfrak{R} \left( e^{2i(\omega_a t - \vec{k}_a \cdot \vec{x}_r + \Phi)} - e^{2i(\omega_a t - \frac{\omega_a L}{c} - \vec{k}_a \cdot \vec{x}_e + \Phi)} \right), \quad (15.3c) $$

<footnote n="2">²To make this integration, one needs to expand the $\sin(\omega_\phi t - \vec{k}_\phi \cdot \vec{x} + \Phi)$ in power of $\vec{k}_\phi \cdot \vec{x}_{AU} \cos(\omega_E t) \ll 1$, integrate each term of the series and then reconstruct the sine. We also considered $\omega_E \ll \omega_\phi$.</footnote>
respectively for dilaton and axion interactions, (using Eq. (8.40)), where all the test masses are made of same alloy, therefore they all share the same mass charge $[Q_M]_d$ (and similarly for the axion mass charge), and $\vec{k}_{DM} = -\omega_{DM}v_{DM}\hat{e}_v/c^2$ ($\hat{k} = -\hat{e}_v$, where $\hat{e}_v$ is the unit vector in the direction of the Sun velocity in the galactic halo).

As shown in Chapter 8.3, the light vector field induces an oscillating acceleration on test masses by a small Lorentz force that acts on them. From Eq. (8.68), the leading order contribution to the acceleration on the test mass A is

$$ \vec{a}_A(t, \vec{x}) = \epsilon e \omega_U \vec{Y} \frac{[Q_{B-L}^A]}{m_A} \sin(\omega_U t - \vec{k}_U \cdot \vec{x} + \Phi). \quad (15.4) $$

Using a similar approach as before with Eq. (7.26), we can derive the one-way Doppler shift associated to this oscillation as [19]

$$ y_{re}^{U} = -(\hat{n}_{re} \cdot \hat{e}_{Y}) \frac{\sqrt{2\mu_{0}\rho_{DM}}\epsilon e [Q_{B-L}]}{\omega_{U}m_{TM}} \mathfrak{R} \left( e^{i(\omega_{U}t - \vec{k}_{U} \cdot \vec{x}_{r} + \Phi)} - e^{i(\omega_{U}t - \frac{\omega_{U}L}{c} - \vec{k}_{U} \cdot \vec{x}_{e} + \Phi)} \right), \quad (15.5) $$

where we consider the test masses have all the same mass $m_{TM}$, $\vec{Y} = |\vec{Y}|\hat{e}_Y$ and, as in the scalar cases, the test masses have the same $B-L$ charge $[Q_{B-L}]$.

## Transfer function of the second generation TDI

Following [19], the transfer function of the one-arm Doppler shift induced by the pure DM scalar field is defined as the Fourier transform of $y_{re}(t)$ Eq. (15.3b) normalized by the constant amplitude

$$ \mathcal{T}_{re}^{\mathrm{DM}}(\omega) = (\hat{n}_{re} \cdot \hat{e}_v) \mathfrak{R} \left( e^{-i(\vec{k} \cdot \vec{x}_r - \Phi)} - e^{-i(\frac{\omega L}{c} + \vec{k} \cdot \vec{x}_e - \Phi)} \right). \quad (15.6) $$

Using Eq. (14.8b), we find the DM transfer function of the second generation X combination in the constant arm length approximation as

$$ \mathcal{T}_X^{\mathrm{DM}}(\omega) = -4 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathfrak{R} \left[ e^{\frac{-3i\omega L}{c}} \left( \mathcal{T}_{13}^{\mathrm{DM}} - \mathcal{T}_{12}^{\mathrm{DM}} + e^{\frac{-i\omega L}{c}} \left( \mathcal{T}_{31}^{\mathrm{DM}} - \mathcal{T}_{21}^{\mathrm{DM}} \right) \right) \right] \quad (15.7a) $$

$$ = -4 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathfrak{R} \left[ e^{\frac{-3i\omega L}{c}} e^{-i(\vec{k} \cdot \vec{x}_1 - \Phi)} \left(1 + e^{\frac{-2i\omega L}{c}}\right) (\hat{n}_{13} - \hat{n}_{12}) \cdot \hat{e}_v - 2e^{\frac{-4i\omega L}{c}} (\hat{n}_{13} \cdot \hat{e}_v e^{-i(\vec{k} \cdot \vec{x}_3 - \Phi)} - \hat{n}_{12} \cdot \hat{e}_v e^{-i(\vec{k} \cdot \vec{x}_2 - \Phi)}) \right] \quad (15.7b) $$

$$ = -8 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathfrak{R} \left[ e^{\frac{-4i\omega L}{c}} \left( (\hat{n}_{13} - \hat{n}_{12}) \cdot \hat{e}_v \cos\left(\frac{\omega L}{c}\right) e^{-i(\vec{k} \cdot \vec{x}_1 - \Phi)} - (\hat{n}_{13} e^{-i(\vec{k} \cdot \vec{x}_3 - \Phi)} - \hat{n}_{12} e^{-i(\vec{k} \cdot \vec{x}_2 - \Phi)}) \cdot \hat{e}_v \right) \right] \quad (15.7c) $$

$$ = -8 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathfrak{R} \left[ e^{-i\left(\frac{4\omega L}{c} + \vec{k} \cdot \vec{x}_1 - \Phi\right)} \left( \hat{n}_{13} \cdot \hat{e}_v \left( \cos\left(\frac{\omega L}{c}\right) - 1 + ikL\hat{n}_{13} \cdot \hat{e}_v \right) - \hat{n}_{12} \cdot \hat{e}_v \left( \cos\left(\frac{\omega L}{c}\right) - 1 + ikL\hat{n}_{12} \cdot \hat{e}_v \right) \right) \right] \quad (15.7d) $$
where we used $\vec{x}_2 = \vec{x}_1 - L\hat{n}_{12}$ and $\vec{x}_3 = \vec{x}_1 - L\hat{n}_{13}$ and we took the first order expansion of the exponent $\propto kL \ll 1$ for all frequencies of interest (for example, at $f = 1$ Hz, which is the maximum frequency of the LISA band, $kL \sim 0.05$, with the galactic velocity $v_{DM} = 10^{-3} c$). The amplitude of the transfer function is then

$$|\mathcal{T}_X^{\text{DM}}(\omega)| = 16 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \times \sqrt{(\hat{n}_{23} \cdot \hat{e}_v)^2 \sin^4\left(\frac{\omega L}{2c}\right) + ((\hat{n}_{13} \cdot \hat{e}_v)^2 - (\hat{n}_{12} \cdot \hat{e}_v)^2)^2 \left(\frac{\omega L |\vec{v}_{\text{DM}}|}{2c^2}\right)^2}, \quad (15.8)$$

where we used $\hat{n}_{13} - \hat{n}_{12} = \hat{n}_{23}$. As it can be noticed from Eq. (15.8), the first term dominates as long as $(\omega L/c)^2 > (|\vec{v}_{\text{DM}}|/c)^2$ (and $(\hat{n}_{23} \cdot \hat{e}_v)^2 \neq 0$), i.e for all frequencies $f \ge 5 \times 10^{-5}$ Hz, which fully contains the LISA band. Then, we can neglect the second term and the amplitude of the transfer function becomes

$$|\mathcal{T}_X^{\text{DM}}(\omega)| \approx 8 \left(\frac{\omega L}{c}\right)^4 |\hat{n}_{23} \cdot \hat{e}_v|, \quad (15.9)$$

i.e the transfer function scales as $f^4$. We do not recover the scaling of the transfer function of [19], because Eq. (15.8) corresponds to the amplitude of the transfer function of the second generation TDI, while [19] used the 1.5 generation for their calculation. However, using Eq. (D.12c), one can show easily that the transfer function of the first generation TDI is $|\mathcal{T}_X^{\text{DM}}(\omega)|/2 \sin(2\omega L/c)$, which allows us to recover the results of [19]$^3$.

## 15.1.2 Monochromatic gravitational waves

### One arm Doppler shift

We first model the GW with frequency $f_{GW} = \omega_{GW}/2\pi$ in the source frame propagating along the $\hat{k}$ direction as [20]

$$h_{\mu\nu}(\xi) = A_+ \cos(\varphi(\xi)) \epsilon_{\mu\nu}^{+} + A_- \sin(\varphi(\xi)) \epsilon_{\mu\nu}^{\times}, \quad (15.10a)$$

where $\xi = t - \hat{k} \cdot \vec{x}/c$ represents the surfaces of constant phase and

$$\varphi(\xi) = \omega_{\text{GW}}\xi + \frac{1}{2}\dot{\omega}_{\text{GW}}\xi^2 + \Phi_{\text{GW}}, \quad (15.10b)$$

where $\Phi_{GW}$ is the initial phase of the wave and where we assume a non zero time derivative of the frequency $\dot{\omega}_{GW}$. In the barycentric reference frame, we define the triad $(\hat{u}, \hat{v}, \hat{k})$ from the ecliptic latitude $\beta$ and ecliptic longitude $\lambda$ as [6, 20]

$$\hat{u} = \begin{pmatrix} \sin(\lambda) \\ -\cos(\lambda) \\ 0 \end{pmatrix}, \quad \hat{v} = \begin{pmatrix} -\sin(\beta) \cos(\lambda) \\ -\sin(\beta) \sin(\lambda) \\ \cos(\beta) \end{pmatrix}, \quad \hat{k} = \hat{u} \times \hat{v} = - \begin{pmatrix} \cos(\beta) \cos(\lambda) \\ \cos(\beta) \sin(\lambda) \\ \sin(\beta) \end{pmatrix}. \quad (15.11)$$

<footnote>$^3$To compute the transfer function, we did not average over all possible DM wind directions $\hat{k}$ as [19] did, since it is fixed to a single direction when the field is a coherent superposition of plane waves.</footnote>
We can further define two functionals $\xi_+(\hat{u}, \hat{v}, \hat{n})$ and $\xi_-(\hat{u}, \hat{v}, \hat{n})$ [6]

$$ \xi_+(\hat{u}, \hat{v}, \hat{n}) = (\hat{u} \cdot \hat{n})^2 - (\hat{v} \cdot \hat{n})^2, \qquad (15.12a) $$

$$ \xi_-(\hat{u}, \hat{v}, \hat{n}) = 2(\hat{u} \cdot \hat{n})(\hat{v} \cdot \hat{n}), \qquad (15.12b) $$

such that, in the barycentric frame (superscript SSB), the one-arm Doppler shift induced by the GW can be written as [6, 20, 21]

$$ y_{\text{re}}^{\text{GW}} = \frac{-\mathcal{A}}{2(1 - \hat{n}_{\text{re}} \cdot \hat{k})} \mathcal{R} \left[ \left( \hat{h}_{+}^{\text{SSB}} \xi_{+}(\hat{u}, \hat{v}, \hat{n}_{\text{re}}) + \hat{h}_{\times}^{\text{SSB}} \xi_{\times}(\hat{u}, \hat{v}, \hat{n}_{\text{re}}) \right) \left( e^{i\varphi(\xi_r)} - e^{i\varphi(\xi_e)} \right) \right] \quad (15.13a) $$

$$ = \frac{-\mathcal{A}\hat{n}_{re}^i \hat{n}_{re}^j}{2(1 - \hat{n}_{re} \cdot \hat{k})} \mathcal{R} \left[ \hat{h}_{ij}^{\text{SSB}} (e^{i\varphi(\xi_r)} - e^{i\varphi(\xi_e)}) \right], \qquad (15.13b) $$

where

$$ \hat{h}_{+}^{\text{SSB}} = \cos(2\Psi) (1 + \cos^2(t)) + 2i \sin(2\Psi) \cos(t), \qquad (15.14a) $$

$$ \hat{h}_{\times}^{\text{SSB}} = -\sin(2\Psi) (1 + \cos^2(t)) + 2i \cos(2\Psi) \cos(t), \qquad (15.14b) $$

$$ \hat{h}_{ij}^{\text{SSB}} = (u_i u_j - v_i v_j) \hat{h}_{+}^{\text{SSB}} + (u_i v_j + u_j v_i) \hat{h}_{\times}^{\text{SSB}}, \qquad (15.14c) $$

where $\Psi$ is the polarization angle, $t$ is the inclination of the source and $\mathcal{A}$ is the amplitude of the wave which depends only on the source parameters, and defined in the barycentric frame as

$$ A_+ = -\mathcal{A}(1 + \cos^2(t)), \quad A_\times = -2\mathcal{A}\cos(t). \qquad (15.15) $$

## Transfer function of the second generation TDI

We now compute the transfer function for a GW (assuming $\dot{f} = 0$ for simplicity), which reads (using Eqs. (15.13b) and (14.8b))

$$ \mathcal{T}_{\times}^{\text{GW}}(\omega) = 2 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathcal{R} \left[ \hat{h}_{ij}^{\text{SSB}} e^{-i\left(\frac{3\omega L}{c} + \vec{k} \cdot \vec{x}_1 - \Phi\right)} \times \right. \qquad (15.16a) $$

$$ \sum_{\ell=2,3} w_{\ell} \frac{\hat{n}_{1\ell}^{i} \hat{n}_{1\ell}^{j}}{1 - (\hat{n}_{1\ell} \cdot \hat{k})^2} \left( (1 + \hat{n}_{1\ell} \cdot \hat{k}) (1 - e^{-i \frac{\omega L}{c} (1 - \hat{n}_{1\ell} \cdot \hat{k})}) - e^{-i \frac{2\omega L}{c} (1 - \hat{n}_{1\ell} \cdot \hat{k})} (1 - e^{i \frac{\omega L}{c} (1 + \hat{n}_{1\ell} \cdot \hat{k})}) \right) \\ = 2 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathcal{R} \left[ \hat{h}_{ij}^{\text{SSB}} e^{-i\left(\frac{5\omega L}{c} + \vec{k} \cdot \vec{x}_1 - \Phi\right)} \times \right. \qquad (15.16b) $$

$$ \sum_{\ell=2,3} w_{\ell} \frac{\hat{n}_{1\ell}^{i} \hat{n}_{1\ell}^{j}}{1 - (\hat{n}_{1\ell} \cdot \hat{k})^2} \left( (\hat{n}_{1\ell} \cdot \hat{k}) (1 - 2e^{i\frac{\omega L}{c}(1+\hat{n}_{1\ell} \cdot \hat{k})}) + e^{\frac{2i\omega L}{c}} + 2ie^{\frac{i\omega L}{c}} \sin\left(\frac{\omega L}{c}\right) \right), $$
<figure><img src="image_19.png" /><figcaption>Figure 15.1: Amplitude of transfer functions of TDI X combinations for both scalar ULDM (in blue) and GW (in orange) from Eqs. (15.8) and (15.16a) respectively. For this plot, we neglected the impact of the geometric factors.</figcaption></figure>

where $w_ℓ = ±1$ respectively for $ℓ = 3, 2$. One can simplify this expression when $ωL/c \ll 2π$ to find the amplitude of the transfer function

$$|\mathcal{T}_X^{\text{GW}}(\omega)| = 8 \left(\frac{\omega L}{c}\right)^3 \left|\hat{n}_{ij}^{\text{SSB}} (\hat{n}_{13}^i \hat{n}_{13}^j - \hat{n}_{12}^i \hat{n}_{12}^j)\right|. \qquad (15.16c)$$

Therefore, at small $ωL/c$, $|\mathcal{T}_X^{\text{GW}}| \propto f^3$ which means that the DM signal will, on average (i.e neglecting the geometric factor), be suppressed by a factor $ωL/c$ compared to the GW signal, as it can be noticed in Fig. 15.1.

As it can be noticed from Eqs. (15.9) and (15.16c) and in Fig. 15.1, at low frequency, the DM and GW transfer function do not scale equally as function of the frequency, in particular it exists a third order term in $ωL/c$ in the GW transfer function and not in the DM one. There are two reasons for this difference. The first one is that the GW travels at speed $c$, while the DM field is non relativistic. Indeed, one can notice that assuming a "relativistic DM" which travels at speed $c$ instead of $10^{-3} c$, the transfer function Eq. (15.7d) scales as $(ωL/c)^3$ at low frequency, as the GW transfer function. This suggests that for a dipolar signal, the velocity of the wave impacts the scaling of the transfer function on frequency. The second reason is that DM is dipolar while GW is quadrupolar. Indeed, assuming now a non-relativistic GW (i.e neglecting terms $\propto kL \ll 1$), one can notice that the GW transfer function Eq. (15.16a) also scales as $(ωL/c)^3$ at low frequency, as the relativistic GW. This suggests that for a quadrupolar signal, the velocity of the wave does not influence the scaling of the transfer function on the frequency. Therefore, it is the fact that the DM produces a non-relativistic and dipolar signal which makes a difference in transfer function compared to the GW. Now, let us give a physical explanation. Assuming a low frequency and non-relativistic wave (i.e $ω \ll c/L \equiv 1/τ$, where $τ$ is the light travel time between spacecrafts; and with velocity $v \ll c$, such that the wavelength is much larger than the armlength $λ \gg L$), the phase of the wave at the location of the emitting and receiving spacecrafts is very close (i.e $ΔΦ \ll 2π$) and therefore, the one-way Doppler seen by the photon is small as well (from Eq. (15.6), it is $O(ωL/c)$). If we now assume a full round trip emitter-receiver-emitter, the wave has not propagated much and therefore, the outward effect is, at leading
order, the same as the inward effect (up to a sign). Therefore, for DM, since we are essentially looking at this inward-outward effect projected on a dipole, the signal decreases (for example, $\mathcal{T}_{13}^{\text{DM}} + \exp(-i\omega L/c)\mathcal{T}_{31}^{\text{DM}}$ is $O(\omega L/c)^2$). For the quadrupolar GW, the outward effect is still almost the same as the inward effect but they also share the same sign, therefore the signal is still $O(\omega L/c)$. This is why the polarization makes a difference (for non relativistic waves). Now, if we assume that both waves are relativistic, but the polarization still differs, the phase of oscillation between the emitting and receiving spacecrafts will be different for both dipolar and quadrupolar waves (even for low frequencies), and therefore, both effects are of the same order of magnitude. This is why DM, as a non relativistic wave and with a dipolar effect on LISA arm, has not the same transfer function as GW.

### 15.1.3 Summary of the expected one-arm signals

In summary, if we consider an ULDM or GW background that produces a Doppler shift at angular frequency $\omega$ (i.e $\omega = \omega_{\phi} = 2\omega_a = \omega_{\text{U}} = \omega_{\text{GW}}$), one can write the one-link Doppler shift in a very concise way, similarly as in [19]

$$y_{re} = \mathfrak{R} \left[ \mu_{re} \left( e^{i(\omega t - \vec{k} \cdot \vec{x}_r + \Phi)} - e^{i(\omega t - \frac{\omega L}{c} - \vec{k} \cdot \vec{x}_e + \Phi)} \right) \right], \quad (15.17a)$$

where the wavevector $\vec{k}$ depends on the velocity of the wave i.e $|\vec{k}| = \omega/c$ for GW and $|\vec{k}| = \omega v_{\text{DM}}/c^2$ for DM, and with the unitless $\mu_{re}$ are given by

$$\mu_{re} = \begin{cases} \frac{\sqrt{16\pi G \rho_{\text{DM}} v_{\text{DM}} [Q_M]_d}}{\omega c^2} & \text{for DM pure scalar field} \\ \frac{32\pi G \rho_{\text{DM}} v_{\text{DM}} E_P^2 [Q_M]_a}{f_a^2 \omega c^3} & \text{for DM pseudo scalar field} \\ -(\hat{n}_{re} \cdot \hat{e}_v) \frac{\sqrt{2\mu_0 \rho_{\text{DM}} \epsilon e [Q_{\text{B-L}}]}}{\omega m_{\text{TM}}} & \text{for DM vector field} \\ -\left(\frac{\hat{n}_{re}^i \hat{n}_{re}^j \hat{h}_{ij}^{\text{SSB}}(t, \Psi)}{2(1 - \hat{n}_{re} \cdot \hat{k})}\right) \mathcal{A} & \text{for gravitational wave} \end{cases} \quad (15.17b)$$

and where we neglected the frequency derivative in the case of the GW $\dot{f} = 0$, for simplicity and concise equations purposes. These $\mu_{re}$ factor are the product of an amplitude and a pure geometric factor which shows how the detector couples to DM or GW, and which gives us indication on the nature of the wave. Indeed, ULDM signal being dipolar, the LISA arm couples to a single direction (which corresponds to either the galactic velocity in the scalar cases, or to the vector polarization direction in the vector case) while GW signal being quadrupolar, the arm couples to the two orthogonal directions corresponding to the quadrupolar GW polarization.

On one hand, among other signals, LISA will be sensitive to (almost) monochromatic GW, i.e GW emitted from galactic binaries whose inspiral phase (with frequency which can be approximated as a constant) lasts longer than the total mission duration [6]. On the other hand, considering that the coherence time of the ULDM field exceeds the time of the mission, which we will consider to be 1 year in the following$^4$, the acceleration induced on the test

<footnote>$^4$Following Eq. (5.11), this implies that we need $f \le 10^{-2}$ Hz.</footnote>
masses would be monochromatic as well. This implies that both monochromatic signals from ULDM and GW Eq. (15.17a) would lead to similar shapes in Fourier space (essentially one Dirac at the wave frequency). This suggests that the two signals might resemble each other, which could lead to degeneracy between them.

In addition, one can find simple detection setup where the quadrupolar nature of the GW signal would be seen by the detector as a dipole, therefore where the ULDM and GW signals would look similar and therefore be degenerate. For example, let us assume a simple detection setup where the test mass A is static and the test mass B is located at a distance L and is rotating counterclockwise around A in the x-y plane with angular frequency $\omega_R$. Both test masses are sending light signals to each other. We assume that the DM field is a pure scalar field of frequency $\omega$ whose propagation direction is along the x-axis in this coordinate system. Let us further assume a monochromatic GW of the same angular frequency $\omega$ which propagates along the y-axis, and whose polarization is purely + in the x-z plane. One can easily deduce that the GW polarization in the z direction will not affect the light signals emitted between the two test masses, only the x polarization will be noticed. Therefore, the GW will be seen as a dipole by the detector.

However, in the case of a more complex detector like LISA, we might wonder if this difference in polarization would clearly lead to a different spectral shape of the signal. In addition to this polarization difference, DM and GW propagates at different speed, i.e $|\vec{k}_{GW}| = \omega_{GW}/c$ while $|\vec{k}_{DM}| = \omega_{DM}\nu_{DM}/c^2 \sim 10^{-3}\omega_{DM}/c$ in Eq. (15.17a), which might also lead to a difference in the spectral signature for both signals. Note that such a difference would not be visible for time scales much smaller than one year, where the geometric factor would simply be constant.

Therefore, in the following, we will try to answer the following question : using realistic orbits, would LISA be able to break the degeneracy between GW emitted from a galactic binary and (scalar) ULDM fields ?

To do so, the first step is to simulate one signal induced by a galactic binary and one signal induced by a (scalar) DM candidate. Then, we will model one galactic binary and one scalar DM candidate using Bayesian analysis on each of these data sets. Finally we will compare the efficiency of these models to see if, for each data set, one is preferred compared to the other.

## 15.2 Bayesian inference to discriminate between dark matter and gravitational waves

### 15.2.1 A short introduction on Bayesian inference

In this section, we will review how Bayesian inference works. It is based one Bayes' theorem, which for two random events A and B, allows to write the probability of A under the condition $B P(A|B)$ as

$$P(A|B) = \frac{P(A)}{P(B)} P(B|A), \qquad (15.18)$$

where $P(B|A)$ is the probability of B under the condition A and $P(A), P(B)$ respectively the probabilities of A and B. When considering a vector of data points $\vec{d}$ depending on a vector
of parameters $\vec{p}$, Eq. (15.18) becomes

$$ \mathcal{P}(\vec{p}|\vec{d}) = \frac{\mathcal{L}(\vec{d}|\vec{p}) \pi(\vec{p})}{\mathcal{N}(\vec{d})}. \qquad (15.19) $$

$\pi(\vec{p})$ is the prior distribution of the parameters $\vec{p}$. This represents what we know of the parameters, prior to looking at the data $\vec{d}$. It can e.g. come from prior measurements on another set of data. For example, the velocity distribution Eq. (5.7b) can be used as a prior distribution for the galactic velocity (this is what will be done in the following). On the contrary, if no information on the given parameter is available, one might use a flat uniform prior distribution. Therefore, in the Bayesian framework, the parameters $\vec{p}$ are treated as random variables and not as fixed constants.

$\mathcal{L}(\vec{d}|\vec{p})$ is the likelihood density function and represents the agreement between a given set of parameters $\vec{p}$ with the data $\vec{d}$. In the most simple case, the data (of length $j$) has white noise with variance $\sigma^2$, and we model the signal as $s(\vec{p})$ such that the likelihood reduces to

$$ \log \mathcal{L}(\vec{d}|\vec{p}) = -\frac{1}{2} \sum_{n=1}^{j} \frac{(d_n - s_n(\vec{p}))^2}{\sigma^2}. \qquad (15.20) $$

Maximizing the likelihood means minimizing the sum, which is very similar to $\chi^2$ method.

$\mathcal{P}(\vec{p}|\vec{d})$ is the posterior probability density function and represents the updated conditional probability of the parameters $\vec{p}$ after observing the data $\vec{d}$. For a given set of parameters $\vec{p}_w$, $\mathcal{P}(\vec{p}_w|\vec{d}) \to 0$ means that $\vec{p}_w$ is not compatible with the data. As the Bayes' theorem suggests it, the posterior distribution is the product of the prior distribution and the likelihood (up to $\mathcal{N}(\vec{d})$, a normalization constant, independent of the parameters), i.e it is the combination of our knowledge of the parameters before and after analyzing the data.

There exists two big classes of algorithms to reconstruct the full posterior probability distribution of the parameters $\vec{p}$ given a dataset $\vec{d}$. The first one gathers Markov Chain Monte Carlo (MCMC) algorithms. At each step $i$ of the algorithm, a random sample of parameters $\vec{p}_i$ is generated from the sample $\vec{p}_{i-1}$ and the acceptance of this new sample only depends on $\vec{p}_i$ and $\vec{p}_{i-1}$. The ratio of posteriors between parameters of $i-1$ and $i$ steps is computed in order to know if newly found parameters explain the data better than the previous one. In general, the distance in parameter space between $\vec{p}_i$ and $\vec{p}_{i-1}$ is small (compared to the parameter space volume), and therefore, the algorithm can be trapped in a local maximum, which we want to avoid. The second method is called Nested sampling [22] and is a way of avoiding being stuck in a local maximum likelihood. We first need to introduce the prior mass $X$ as the fraction of prior contained within an iso-likelihood contour. Since we expect large likelihood to be contained in a small region of the parameter space, the likelihood $\mathcal{L}(X)$ must be a decreasing function of $X$. The principle of nested sampling is to choose a set of $N$ random points covering uniformly the parameter space, and which have different likelihoods. At the iteration $i$, one deletes the point $N_i$ with lowest likelihood $\mathcal{L}(X_i)$ and finds a new point in the parameter space $N_{i+1}$ such that $\mathcal{L}(X_{i+1}) > \mathcal{L}(X_i)$. This way, the remaining prior mass decreases and converges to the region with largest likelihood.
Since one can store the value of $\mathcal{L}(X_i)$ at each iteration, one has access to the integral of $\mathcal{L}(X)$, which corresponds to the model evidence. When comparing two models, one uses the Bayes factor, which is the ratio of model evidence. More specifically, one can write Bayes' rule in terms of $\vec{p}$ which depends on $\vec{d}$ but also on a specific model $M$, i.e

$$ \mathcal{P}(\vec{p}|\vec{d}, M) = \frac{\mathcal{L}(\vec{d}|\vec{p}, M)\pi(\vec{p}|M)}{p(\vec{d}|M)}, \qquad (15.21) $$

where $p(\vec{d}|M)$ is known as the model evidence, which is the probability of the data $\vec{d}$ given $M$. This does not depend on the parameters $\vec{p}$, which have been integrated out (or marginalized out). Considering two given models $M_1, M_2$, the Bayes factor

$$ \mathcal{B} = \frac{p(\vec{d}|M_1)}{p(\vec{d}|M_2)}, \qquad (15.22) $$

indicates which model best explain the data. Therefore, nested sampling is very convenient to compare different models on the same data. In the following, we will be using a nested sampler algorithm, called *Nessai* [23] for the parameter estimation of both a galactic binary and a scalar ULDM field.

### 15.2.2 Fast likelihood modeling

We now focus on the analysis of TDI dataset induced by a galactic binary (GB) on one side and (scalar) ULDM on the other side. For a given model, we want to optimize the computation of the likelihood, in particular to make it faster. To do so, we work in the Fourier domain, because it is much more convenient for the signals that we are interested in. Indeed, for (quasi) monochromatic signals, the idea is to separate the fast oscillation (at the wave frequency) whose Fourier transform can be done analytically, from the slowly oscillating part which can be Fourier transformed numerically [24]. This speeds up the numerical integration necessary for the computation of the likelihood. As we shall see in the following, the TDI combination are very simple to implement in Fourier space.

In Appendix D, we show explicitly how to express the slow oscillating part of TDI $X_2$ combination for both scalar ULDM and GB signals. The two *Python* softwares that we use which follow this idea, are *FastGB* and *FastDM*, respectively for the modeling of GB and scalar DM.

### 15.2.3 Generation of signals and parameter space

We now explicitly describe the general process from the data simulation to the data modeling.

We first simulate a dataset which contains either GB or scalar DM signal. To do so, we simulate in time domain the second generation TDI Eq. (14.7a) induced by either a GB or a scalar ULDM signals through Doppler effects derived in the previous sections, see Eqs. (15.3b) and (15.13b). The positions of the spacecrafts are obtained using *LISAOrbits* [25], and we use Earth-trailing orbits for their orbits [26], which we simulate for one full year, i.e $T_{obs} = 365 \times 86400$ s. The positions of the spacecrafts are sampled with frequency $f_s = 0.1$ Hz,
i.e $N = T_{\text{obs}}f_s$ measurements are performed. For the generation of data, no noise is added i.e these are "perfect" dataset<fnref n="5" />.

For the data modeling, we fit the TDI A, E data combinations expressed in Fourier domain. More precisely, writing the TDI data in time domain $d^{A,E}(t)$ as a vector ($d_0^{A,E}$, $d_1^{A,E}$, ..., $d_{N-1}^{A,E}$) where $d_i^{A,E} = d^{A,E}(i/f_s)$ is the ith sampled element of $d^{A,E}(t)$, we perform the discrete Fourier transform (DFT) $\tilde{d}^{A,E}$ defined as

$$ \tilde{d}_k^{A,E} = \sum_{\ell=0}^{N-1} e^{-\frac{2\pi i \ell k}{N}} d_{\ell}^{A,E}, \qquad (15.23) $$

and we fit a model $\tilde{m}^{A,E}$ to such dataset in Fourier space. We use a Gaussian likelihood expressed in time domain as [27]

$$ \log \mathcal{L} = \sum_{x=A,E} \left( -\frac{1}{2} \log (\det(2\pi \mathbf{C}^x)) - \frac{1}{2} (\tilde{d}^x - \tilde{m}^x)^{\dagger} (\mathbf{C}^x)^{-1} (\tilde{d}^x - \tilde{m}^x) \right), \qquad (15.24a) $$

where $\tilde{d}$, $\tilde{m}$ represent respectively the time vector or data and model and where $\mathbf{C}^x$ is the noise covariance matrix of the x TDI combination. Using the DFT of the data and model, the likelihood can be written as [28]

$$ \log \mathcal{L} = \sum_{x=A,E} \left( -\frac{1}{2} \log (2\pi \Pi_k \tilde{C}_{kk}^x) - \sum_{j=0}^{N/2} \frac{|\tilde{d}_j^x - \tilde{m}_j^x|^2}{\tilde{C}_{jj}^x} \right), \qquad (15.24b) $$

where $\tilde{C}_{kk}^x$ is the diagonal element at Fourier bin k of the two-sided PSD matrix of the $x = A, E$ TDI combinations defined as [28]

$$ \tilde{C}_{kk}^x = N f_s \delta_{kl} N_x(f_l), \qquad (15.24c) $$

where $\delta_{kl}$ is the Kronecker delta and therefore $\tilde{C}_{kk}^x$ is diagonal because the noise is stationary [27], and where the $N_x(f_k)$ is the noise PSD of the second TDI A, E combinations at frequency bin k, which are defined as [11]

$$ N_{A,E}(f) = 32 \sin^2\left(\frac{2\pi f L}{c}\right) \sin^2\left(\frac{4\pi f L}{c}\right) \left(2\left(3 + 2\cos\left(\frac{2\pi f L}{c}\right) + \cos\left(\frac{4\pi f L}{c}\right)\right) S_{\text{acc}}(f) + \left(2 + \cos\left(\frac{2\pi f L}{c}\right)\right) S_{\text{oms}}(f)\right), \qquad (15.25) $$

assuming no correlation between noises, and which is obtained by extending the single-link noise PSD Eq. (14.11) on the full TDI combination. Since the noise PSD are equivalent for A, E combinations, we can rewrite Eq. (15.24b) as

$$ \log \mathcal{L} = -\frac{1}{2} \log \left( (2\pi)^2 f_s \Pi_k N_{A,E}(f_k) \right) - \sum_{x=A,E} \left( \frac{1}{N f_s} \sum_{j=0}^{N/2} \frac{|\tilde{d}_j^x - \tilde{m}_j^x|^2}{N_x(f_j)} \right). \qquad (15.26) $$

<footnote n="5">⁵Therefore, the results of the analysis, and its possible bias, will not depend on noise.</footnote>
The GB parameters that will be fitted (using Eqs. (D.6b), (D.6a) and (D.6c)) are the amplitude *A*, the frequency $f_{GW}$, the drift in frequency $\dot{f}_{GW}$, the sky localization $(\beta, \lambda)$, the inclination $i$, the polarization angle $\Psi$ and the phase $\Phi_{GW}$. There is no specific knowledge on the physics of one particular GB, therefore we will use isotropic prior distributions for each of those parameters.

The scalar DM parameters that will be fitted (using Eqs. (D.3), (D.5b) and (D.5c)) are the mass charge<fnref n="6" /> $[Q_M]_d \equiv \epsilon$, the frequency $f_\phi$, the direction of the DM 'wind' $(\beta_{DM}, \lambda_{DM})$, its amplitude, which is given by the galactic velocity amplitude $|\vec{v}_{DM}|$ and the phase $\Phi$. For most of these parameters, we will use isotropic priors, except for the galactic velocity amplitude prior that will be given by the amplitude of the galactic velocity distribution Eq. (5.7b). Even though we are working on time scales lower than the coherence time of the field, such that the DM wave is essentially monochromatic, the prior will help constraining the simulation to consistent regions of the parameter space<fnref n="7" />. Note that, by choosing $\vec{v}_{DM}$ as a random variable and not a fixed parameter, we will derive realistic sensitivity curves for ULDM of LISA. By doing so, we simplify the system by considering that the velocity amplitude $|\vec{v}_{DM}|$ is uncorrelated with the DM wind direction, while it is in reality not true as both parameters are constrained together by the same distribution Eq. (5.7a)<fnref n="8" />.

### 15.2.4 Model of galactic binary onto galactic binary signal

We first fit a dataset $d_{GW}$ from a galactic binary using FastGB whose parameters $D_{GW}$ are shown on the left of Table 15.1, and check for efficiency of the model in reproducing data. The 2D posterior distributions are shown in Fig. 15.2, with the joint posterior distribution for all parameters. The yellow lines represent the true values of the parameters, i.e the values injected in the simulation. One can notice strong correlations between the amplitude *A* and the inclination *i*, and between the ecliptic latitude *β* and the drift in frequency *f*. One can also notice that for the polarization angle *Ψ* and initial phase *Φ* parameters, two modes exist which are respectively shifted by *π*/2 and *π*. From Eq. (15.13b), it is straightforward to show that a change in polarization angle *Ψ* → *Ψ* + *π*/2 is equivalent to a change in phase *Φ* → *Φ* + *π*. Overall, the analysis recovers correctly the injected parameters.

<footnote n="6">⁶From Eq. (8.6a), we remind that $[Q_M]_d$ is a sum of partial dilatonic mass charges and dilatonic couplings.</footnote>

<footnote n="7">⁷As we shall see in the following, for most frequencies under consideration, the phase shift induced by DM propagation (through $|\vec{v}_{DM}|$) is small, and therefore hard to fit with accuracy.</footnote>

<footnote n="8">⁸As mentioned above, we will use isotropic priors for both ecliptic latitude and longitude of the DM wind, while the galactic velocity distribution Eq. (5.7a) provides prior information on such direction. First, we make this choice because there is no analytical expression for the angular distributions only. Second, by doing the numerical integration of Eq. (5.7a) over the velocity amplitude, one can show that the width of the angular distributions are much larger than the posterior distributions of these parameters (see below, Section 15.2.5), and therefore using isotropic prior will not impact the results.</footnote>
<figure><img src="image_20.png" /><figcaption>Figure 15.2: Results of the fit of a galactic binary signal by a galactic binary model. The yellow lines indicate the true values of the parameters. The analysis recovers correctly the injected parameters.</figcaption></figure>
The model best fit parameters $M_{\text{GW}}^{\text{GW}}$, that we define as the posterior with largest likelihood, in addition to the mean and standard deviation of the distribution of posteriors are shown in Table 15.1 (only the ~ 10 000 last iterations are taken into account by the algorithm). Let us now analyze the result of the fit, in particular, its efficiency. First, one can compare the signal to noise levels (i.e compute the signal-to-noise ratio (SNR)) for the GW using [16]$^9$

$$ \sqrt{\text{SNR}} = \sqrt{4\mathfrak{R}\left(\int df \frac{\tilde{d}_{\text{GW}}\tilde{d}_{\text{GW}}^{\dagger}}{N(f)}\right)}, \quad (15.27) $$

where $\tilde{d}_{\text{GW}}^{\dagger}$ is the conjugate transpose of $\tilde{d}_{\text{GW}}$, the GW data in Fourier domain, and $N(f)$ is the noise PSD. For the GW data from Table 15.1, one finds

$$ \sqrt{\text{SNR}}_{A}^{\text{GW}} \approx 632 \qquad (15.28a) $$

$$ \sqrt{\text{SNR}}_{E}^{\text{GW}} \approx 679, \qquad (15.28b) $$

for both A and E TDI combinations. The signal power is well above the noise, as it can be noticed in Fig. 15.3.

In a second step, in order to estimate the accuracy of the model, we simulate the TDI time series of a galactic binary with parameters given by $M_{\text{GW}}^{\text{GW}}$, which we call $m_{\text{GW}}^{\text{GW}}$ and we construct the residuals (in time domain)

$$ r_{\text{GW,GW}} = d_{\text{GW}} - m_{\text{GW}}^{\text{GW}}, \qquad (15.29) $$

which we compare with the LISA noise PSD. In Fig. 15.3, we compare the residuals to LISA noise in Fourier space. One can notice that the residuals are well below the noise, indicating that the model extracted all the necessary information out of the data$^{10}$. We now compute the residuals-to-noise ratio (RNR) between the residual in Fourier domain $\tilde{r}_{\text{GW,GW}}$ given by [16]

$$ \sqrt{\text{RNR}} = \sqrt{4\mathfrak{R}\left(\int df \frac{\tilde{r}\tilde{r}^{\dagger}}{N(f)}\right)}, \qquad (15.30) $$

<table><caption>Table 15.1: On the left column, we show the GW injected parameters in the simulation. On the middle, we show the GW model best fit, i.e the posterior with the largest likelihood. The model correctly recovers the injected parameters with less than 1% deviation. On the right column, we present the mean and standard deviation of the full distribution of posteriors.</caption><thead><tr><th></th><th><i>D</i><sub>GW</sub></th><th><i>M</i><sub>GW</sub><sup>GW</sup></th><th>Mean ± 1 σ</th></tr></thead><tbody><tr><td><i>A</i> (10<sup>-21</sup>)</td><td>1.397</td><td>1.398</td><td>1.397 ± 0.003</td></tr><tr><td><i>f</i> (mHz)</td><td>4.215</td><td>4.215</td><td>4.215 ± 0.000</td></tr><tr><td><i>f</i>' (aHz/s)</td><td>55.91</td><td>56.26</td><td>57.36 ± 18.20</td></tr><tr><td>β (rad)</td><td>0.817</td><td>0.8164</td><td>0.8165 ± 0.0005</td></tr><tr><td>λ (rad)</td><td>5.149</td><td>5.1485</td><td>5.1485 ± 0.0004</td></tr><tr><td>ι (rad)</td><td>1.047</td><td>1.048</td><td>1.047 ± 0.002</td></tr><tr><td rowspan="2">Ψ (rad)</td><td rowspan="2">0.785</td><td rowspan="2">0.785</td><td>0.785 ± 0.002</td></tr><tr><td>2.356 ± 0.002</td></tr><tr><td rowspan="2">Φ (rad)</td><td rowspan="2">4.890</td><td rowspan="2">4.888</td><td>4.889 ± 0.007</td></tr><tr><td>1.748 ± 0.007</td></tr></tbody></table>

<footnote>$^9$The definition of SNR is slightly changed compared to [16], in order to be consistent with its previous definition of ratio of signal to noise power.</footnote>

<footnote>$^{10}$Despite its efficiency, the residuals still present a peak at the GW frequency. This is due to the approximations done in FastGB (e.g. neglecting the spacecraft velocity during light travel time).</footnote>
where $\tilde{r}^\dagger$ is the conjugate transpose of $\tilde{r}$. For the GW model on GW data, one finds

$$ \sqrt{\mathrm{RNR}_A^{\mathrm{GW,GW}}} \approx 0.54 \qquad (15.31a) $$

$$ \sqrt{\mathrm{RNR}_E^{\mathrm{GW,GW}}} \approx 0.41. \qquad (15.31b) $$

The RNR are well below unity and the SNR values, indicating that the model extracts all the necessary information out of the data, as expected from what is shown in Fig. 15.3.

<figure><img src="image_21.png" /><figcaption>Figure 15.3: GW signal PSD (in lime) compared to LISA noise PSD of the TDI A combination (in red) in Fourier domain. The signal is well above the noise. We also show the $r_{GW,GW}$ residuals power in grey, which are below the noise.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 15.4: DM signal PSD (in lime) compared to LISA noise PSD of TDI A combination (in red) in Fourier domain. The signal is above the noise. We also show the $r_{DM,DM}$ residuals power in grey, which are below the noise.</figcaption></figure>

### 15.2.5 Model of scalar dark matter onto a scalar dark matter signal

Using *FastDM*, we now turn to the analysis of DM signal. For consistency with the galactic binary that was analyzed in the previous section, we consider a DM signal with the same frequency (in order to compare to the same noise level) and intrinsic amplitude of oscillation on the test masses. From Eq. (15.17b), this means that we choose the coupling $\epsilon$ such that

$$ \epsilon = \mathcal{A} \frac{2\pi f c^2}{\sqrt{16\pi G \rho_{\mathrm{DM}} v_{\mathrm{DM}}}}, \qquad (15.32) $$
where $f, \mathcal{A}$ are given in Table 15.1.

In addition, we will use the mean velocity direction and amplitude of the velocity distribution Eq. (5.7a). The mean direction is well known and it points towards ($\beta_{DM}, \lambda_{DM}$) = (59.91°, 24.67°) ~ (1.046 rad, 0.431 rad)<fnref n="11" />. This is the DM wind direction we will use. For the amplitude, we use $v_{DM} = 3 \times 10^5$ m/s for the velocity of the wave, i.e the mean of the galactic velocity distribution Eq. (5.7a). The result of the fit of a DM model onto a DM signal is presented in Fig. 15.5

One can notice that there exists a strong correlation between the coupling $\epsilon$ and the galactic velocity parameters $v_{DM}$, as expected from the amplitude of the Doppler Eq. (15.3b), which is proportional to both parameters. One can notice from Table 15.2, that the analysis correctly retrieve the correct value of the true parameter (at the 1% error level), except for the coupling $\epsilon$ and galactic velocity $v_{DM}$, whose best fit values differ from the true value of ~ 15%, and whose posterior distribution width is much larger than for the rest of the parameters ($\sigma_{v_{DM}}/v_{DM} \sim \sigma_{\epsilon}/\epsilon \sim 1$ while for the rest of the parameters $\sigma_{p}/p \ll 1$). The reason is that at this frequency, LISA cannot probe accurately the DM wind velocity, since it produces a phase in the Doppler Eq. (15.3b) which is too small to be seen. As a consequence, we have a large uncertainty on this parameter, which is translated to the coupling due to the large correlation between the two parameters. Indeed, one can notice from Table 15.2 that the typical width of the phase $\Phi$ is $\sigma_{\Phi} = 10^{-2}$ rad. As we have used a uniform prior on $\Phi$, the width of the posterior is essentially the width on the likelihood $\mathcal{L}(d_{DM}|\Phi)$, which is directly related to the noise PSD. At $f \sim 4 \times 10^{-3}$ Hz, the amplitude of the phase induced by the propagation of the field is $\sim \omega v_{DM}|\vec{x}|\cos(\beta)/c^2 \sim 6 \times 10^{-3} < \sigma_{\Phi}$ (where $\beta$ is the ecliptic latitude of the DM wind direction, considering $|\vec{x}| \sim 1.5 \times 10^{11}$ m $\equiv$ 1 AU, as mentioned previously and that the orbit of the spacecrafts are close to the ecliptic plane). In Section 15.3, we will analytically derive the typical width of the coupling $\epsilon$ parameter to be able to predict its value at any frequency in the LISA band.

For now, let us compute the SNR of such DM data

$$ \sqrt{\text{SNR}_{A}^{\text{DM}}} \approx 77 \qquad (15.33a) $$

$$ \sqrt{\text{SNR}_{E}^{\text{DM}}} \approx 65, \qquad (15.33b) $$

i.e the signal is still higher than the noise, but the SNR is smaller than in the case of the

<table><caption>Table 15.2: On the left column, we show the DM injected parameters in the simulation. On the middle, we show the DM model best fit, i.e the posterior with the largest likelihood. On the right column, we present the mean and standard deviation of the full distribution of posteriors.</caption><thead><tr><th></th><th>D<sub>DM</sub></th><th>M<sup>DM</sup><sub>DM</sub></th><th>Mean ± 1 σ</th></tr></thead><tbody><tr><td>ε (10<sup>-5</sup>)</td><td>2.396</td><td>2.785</td><td>2.430 ± 0.855</td></tr><tr><td>f (mHz)</td><td>4.215</td><td>4.215</td><td>4.215 ± 0.000</td></tr><tr><td>β (rad)</td><td>1.046</td><td>1.048</td><td>1.045 ± 0.010</td></tr><tr><td>λ (rad)</td><td>0.431</td><td>0.435</td><td>0.430 ± 0.021</td></tr><tr><td>v<sub>DM</sub> (10<sup>5</sup>) (m/s)</td><td>2.998</td><td>2.580</td><td>3.305 ± 1.079</td></tr><tr><td>Φ (rad)</td><td>3.444</td><td>3.444</td><td>3.444 ± 0.010</td></tr></tbody></table>

<footnote n="11">¹¹The Sun velocity in the galactic halo points towards $\alpha$ Cygni, the biggest star of the Cygnus constellation [29], which corresponds to a right ascension $\alpha_{DM} = 310.36°E$ and declination $\delta_{DM} = 45.28°N$ [30], in the equatorial frame. By transforming it to the ecliptic coordinate system [31], we find ($\beta_{DM}, \lambda_{DM}$) = (59.91°, 24.67°).</footnote>
GW data. This difference in $\sqrt{\text{SNR}}$ between DM and GW of a factor ~ 10 comes from the difference in transfer functions at low frequency Eqs. (15.9) and (15.16c)$^{12}$.

Eq. (15.33) is consistent with what is shown in Fig. 15.4, where we notice that the DM signal PSD is well above LISA noise$^{13}$. Calling $m_{\text{DM}}^{\text{DM}}$ the TDI time series of the model with parameters $M_{\text{DM}}^{\text{DM}}$, we can define the residuals as

$$r_{\text{DM,DM}} = d_{\text{DM}} - m_{\text{DM}}^{\text{DM}}, \qquad (15.34)$$

and compare them to the noise level. The RNR of the DM model on the DM data gives

$$\sqrt{\text{RNR}_{A}^{\text{DM,DM}}} \approx 0.35 \qquad (15.35a)$$

$$\sqrt{\text{RNR}_{E}^{\text{DM,DM}}} \approx 0.18, \qquad (15.35b)$$

smaller than 1 and the SNR values Eq. (15.33). As in the GW case above, this suggests, that the model extracts all the information out of the data. We can also check the residuals power in Fourier space Fig. 15.4, and we notice that it is below the noise, as expected.

<footnote>$^{12}$From those equations and Eq. (14.9), one can easily obtain the low frequency transfer functions of e.g. the TDI A combination. Then, using the true value of DM and GW parameters (i.e $f, \hat{c}_v$ for DM and $f, \Psi, \hat{k}$ for GW), one can compute numerically the various geometric factors $\hat{n} \cdot \hat{c}_v$ and $\hat{h}_{ij}^{\text{SSB}} \hat{n}^i \hat{n}^j$ from the spacecrafts orbits. We find $T_A^{\text{GW}}/T_A^{\text{DM}} \sim 8.3$ while the ratio of Eqs. (15.28) and (15.33) gives $\sim 8.2$ for TDI A combination. Note that from Fig. 15.1, the ratio of the two X transfer functions is $\sim 3$ at $f \sim 4.215$ mHz, but this is because we did not take into account the geometric factors for this plot.</footnote>

<footnote>$^{13}$The Fourier peak is narrower compared to the GW case, because the DM signal is completely monochromatic (up to small deviations from the orbit), while we implemented a small frequency drift for the GW.</footnote>
<figure><img src="image_23.png" /><figcaption>Figure 15.5: Results of the fit of the DM signal by a DM model. The yellow lines indicate the true values of the parameters. One can notice strong correlations between the coupling $\epsilon$ and the velocity $v_{DM}$.</figcaption></figure>
### 15.2.6 Can scalar dark matter be misinterpreted as a galactic binary ?

As it was mentioned in the beginning of this chapter, we now tackle one of the main goals of this study, i.e we try to answer the following question : can a real DM signal be misinterpreted as a GB ? To do so, we simulate the same DM signal as previously (the parameters are shown in Table 15.2) and we fit this signal with a GB model using *FastGB*. The fitted parameters are shown in Table 15.3. We use the same procedure as before, i.e we compute the residuals between the GW model time series $m_{DM}^{GW}$ (from the model parameters $M_{DM}^{GW}$, see Table 15.3) and the DM data $d_{DM}$

$$r_{GW,DM} = d_{DM} - m_{DM}^{GW}, \quad (15.36)$$

which allows us to compute the RNR

$$\sqrt{\mathrm{RNR}_{A}^{\mathrm{GW,DM}}} \approx 51 \qquad (15.37a)$$

$$\sqrt{\mathrm{RNR}_{E}^{\mathrm{GW,DM}}} \approx 52. \qquad (15.37b)$$

Those values must be compared with the SNR values of DM data Eq. (15.33).

One can notice that the SNR and RNR are of the same order of magnitude, and are larger than 1, meaning that the model does not fit the data (see Fig. 15.6). Now, we wonder if *LISA* will be able to non ambiguously detect DM. To do so, we use Bayesian model selection criteria which allow one to quantitatively compare the efficiency of the two different models, the GW $m_{DM}^{GW}$ and DM $m_{DM}^{DM}$ models, on the same dataset $d_{DM}$. More precisely, we use the Bayes factor $\mathcal{B}$ which gives the ratio of the probability of getting the dataset $d_{DM}$ under the condition of one of the model, e.g. $m_{DM}^{DM}$, to the probability of getting $d_{DM}$ under the condition of the second model, $m_{DM}^{GW}$. This is equivalent to the ratio of the evidences (see Eq. 15.22). For each model, the sampler provides us the evidence of the given model. For the dataset and models under consideration, the Bayes factor $\mathcal{B}$ is

$$\log \mathcal{B} = \log p(d_{\mathrm{DM}}|m_{\mathrm{DM}}^{\mathrm{DM}}) - \log p(d_{\mathrm{DM}}|m_{\mathrm{DM}}^{\mathrm{GW}}) = O(10^3) \gg 1, \quad (15.38)$$

i.e given the data that contains a simulated signal, the DM model $m_{DM}^{DM}$ is largely preferred compared to the GW model. Therefore, a priori, *LISA* is able to make the difference between one GB and scalar DM. This suggests that if a scalar DM field exists with coupling and oscillation frequency in *LISA* range, *LISA* will be able to uncorrelate it from single galactic binaries. In the following section, we try to answer the opposite question, namely if it is possible that a galactic binary is seen as DM by *LISA*, i.e if both GW and DM models efficiently fit a GB data.

<table><thead><tr><th></th><th><i>D</i><sub>DM</sub></th><th><i>M</i><sub>DM</sub><sup>GW</sup></th><th>Mean ± 1 σ</th></tr></thead><tbody><tr><td><i>A</i> (10<sup>-21</sup>)</td><td>1.397</td><td>0.129</td><td>0.128 ± 0.004</td></tr><tr><td><i>f</i> (mHz)</td><td>4.215</td><td>4.215</td><td>4.215 ± 0.000</td></tr><tr><td><i>f</i>' (aHz/s)</td><td>—</td><td>-8.040</td><td>38.14 ± 208.6</td></tr><tr><td>β (rad)</td><td>1.046</td><td>-1.520</td><td>-1.522 ± 0.004</td></tr><tr><td>λ (rad)</td><td>0.431</td><td>5.326</td><td>5.307 ± 0.021</td></tr><tr><td>ι (rad)</td><td>—</td><td>1.990</td><td>1.995 ± 0.024</td></tr><tr><td>Ψ (rad)</td><td>—</td><td>3.114</td><td>3.079 ± 0.046<br>1.571 ± 0.082<br>0.067 ± 0.053</td></tr><tr><td>Φ (rad)</td><td>3.444</td><td>4.740</td><td>4.740 ± 0.066<br>1.597 ± 0.069</td></tr></tbody></table>

Table 15.3: On the left, the GW parameters injected in the simulation equivalent to the DM parameters shown in Table 15.2 (see Eq. (15.32)). On the middle and on the right, the GW best fit and the mean and standard deviation of the full distribution of posteriors. One can notice different modes for $\Psi$ and $\Phi$ parameters.
<figure><img src="image_24.png" /><figcaption>Figure 15.6: DM signal PSD (in lime) compared to LISA noise PSD of TDI A combination (in red) in Fourier domain (as in Fig. 15.4). This time, the $r_{GW,DM}$ residuals power in grey are close to the signal power and above the noise. This indicates the model fits poorly the data.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 15.7: GW signal PSD (in lime) compared to LISA noise PSD of TDI A combination (in red) in Fourier domain (as in Fig. 15.3). This time, the $r_{DM,GW}$ residuals power in grey are close to the signal power and above the noise. This indicates the model fits poorly the data.</figcaption></figure>

### 15.2.7 Can a galactic binary be misinterpreted as scalar dark matter?

We now focus on the opposite situation compared to the previous section, namely we use a dataset from a GB which we fit with a DM model. For this situation, the results are shown in Table 15.4. Here, the $v$ parameter denotes the velocity of the wave, for the GW data, it corresponds to the speed of light, while the DM model is constrained by the DM velocity distribution, so it is driven towards $\sim 10^{-3}$ c.

Computing the RNR of the best fit DM model on GW data, one finds

$$ \sqrt{\mathrm{RNR}_{A}^{\mathrm{DM,GW}}} \approx 572 \qquad (15.39a) $$

$$ \sqrt{\mathrm{RNR}_{E}^{\mathrm{DM,GW}}} \approx 682. \qquad (15.39b) $$

The RNR values are still much larger than 1, and of the same order of magnitude as the GW data SNR Eq. (15.28), implying that the model is inefficient<fnref n="14" />. The Bayes factor in this case is

$$ \log \mathcal{B} = \log p(D_{\mathrm{GW}}|M_{\mathrm{GW}}^{\mathrm{GW}}) - \log p(D_{\mathrm{DM}}|M_{\mathrm{GW}}^{\mathrm{DM}}) = \mathcal{O}(10^5) \gg 1, \qquad (15.40) $$

which indicates that the GW model is by far preferred, given the GW dataset that we use. Therefore, we can conclude that LISA is able to distinguish between scalar DM and one single GB.

<footnote n="14">¹⁴Note that, here, we used the galactic velocity distribution as a prior, but for completeness, we can also relax this velocity constraint, i.e use a uniform prior in $[0, c]$. In this case, the best fit presents a much higher velocity (of order $10^7$ m/s), but the resulting RNR is of the same order of magnitude as before $\sim 500$.</footnote>
<table><caption>Table 15.4: On the left column, GW injected parameters. On the middle and right sides, we show respectively the DM model best fit and the mean and standard deviation of the full distribution of posteriors.</caption><thead><tr><th></th><th>D<sub>GW</sub></th><th>M<sup>DM</sup><sub>GW</sub></th><th>Mean ± 1 σ</th></tr></thead><tbody><tr><td>ε (10<sup>-5</sup>)</td><td>2.396</td><td>1.989</td><td>2.006 ± 0.014</td></tr><tr><td>f (mHz)</td><td>4.215</td><td>4.215</td><td>4.215 ± 0.000</td></tr><tr><td>β (rad)</td><td>0.817</td><td>-0.101</td><td>-0.098 ± 0.004</td></tr><tr><td>λ (rad)</td><td>5.149</td><td>0.806</td><td>0.804 ± 0.005</td></tr><tr><td>ν (10<sup>6</sup>) (m/s)</td><td>299.8</td><td>1.014</td><td>1.003 ± 0.010</td></tr><tr><td>Φ (rad)</td><td>4.890</td><td>5.521</td><td>5.521 ± 0.007</td></tr><tr><td>˙f (aHz/s)</td><td>55.91</td><td>—</td><td>—</td></tr><tr><td>ι (rad)</td><td>1.047</td><td>—</td><td>—</td></tr><tr><td>Ψ (rad)</td><td>0.785</td><td>—</td><td>—</td></tr></tbody></table>

### 15.2.8 What about two galactic binaries ?

In the beginning of this chapter, when discussing about the DM and GW transfer functions, we mentioned that the two main differences between DM and GW are the polarization and the velocity of the wave, i.e the former is effectively non propagating at the scale of TDI, while the latter travels at speed $c$. These are the reasons why LISA can a priori make the difference between both signals, as we have seen in the previous sections. We can wonder if, by considering the signals emitted by multiple GB simultaneously, this could change the story, i.e if LISA could be misled and interpret this effective signal as scalar DM.

It is not possible to construct a dipolar signal from one or multiple quadrupolar signals, therefore the polarizations will still differ. However, one can imagine situations where two GB with opposite sky localizations emit GW with same frequency and amplitude, which therefore would create a standing GW. In such (hypothetical) situations, both scalar DM wave and standing GW would not propagate at the time scale of one TDI combination (at leading order).

Let us now derive the transfer function of such system. Assuming same GW frequency, same amplitude $\mathcal{A}$, phase $\Phi$ and $\iota = \dot{f} = 0$ for both binaries, but opposite source localization $\hat{k}$, the transfer function of the two binaries combined is (following Eq. (15.16a))

$$ \mathcal{T}_X^{\text{Double GW}}(\omega) = 8 \sin\left(\frac{\omega L}{c}\right) \sin\left(\frac{2\omega L}{c}\right) \mathcal{R} \left[ i \hat{h}_{ij}^{\text{SSB}} e^{-i\left(\frac{6\omega L}{c} - \Phi\right)} \sum_{\ell=2,3} w_\ell \frac{\hat{n}_{1\ell}^i \hat{n}_{1\ell}^j}{1 - (\hat{n}_{1\ell} \cdot \hat{k})^2} \times \left( \cos(\vec{k} \cdot \vec{x}_1) \sin\left(\frac{\omega L}{c}\right) - (\hat{n}_{1\ell} \cdot \hat{k}) \left( \sin(\vec{k} \cdot \vec{x}_1) \cos\left(\frac{\omega L}{c}\right) - \sin\left(\vec{k} \cdot \vec{x}_1 - \frac{\omega L \hat{n}_{1\ell} \cdot \hat{k}}{c}\right) \right) \right) \right], \quad (15.41a) $$

i.e a standing wave, as expected. When $\omega L/c \ll 2\pi$, this reduces to

$$ |\mathcal{T}_X^{\text{Double GW}}(\omega)| \approx 16 \left(\frac{\omega L}{c}\right)^3 \left|\hat{h}_{ij}^{\text{SSB}} (\hat{n}_{13}^i \hat{n}_{13}^j - \hat{n}_{12}^i \hat{n}_{12}^j) \cos(\vec{k} \cdot \vec{x}_1)\right|. \quad (15.41b) $$

This transfer function is similar to the one of a single GW (up to a factor 2 as expected), and with extinction positions when $\hat{k}$ is orthogonal to $\hat{x}_1$. Note that the standing wave wavelength
does still correspond to $k_{GW} = \omega/c$, while in the case of DM, the de Broglie wavelength of the fields is $10^3$ larger (since $\lambda_{DM} = 2\pi/|\vec{k}_{DM}| = 2\pi c^2/\omega v_{DM} \sim 10^3 \lambda_{GW}$). If we assume that those GWs have the smallest frequency in LISA band, i.e $f = 10^{-4}$ Hz, $\lambda_{GW} \gg |\vec{x}_1|$, such that the double GW are seen as an homogeneous standing wave by LISA, as the scalar DM field.

We artificially generate a signal from the superposition of two GBs with opposite sky localizations. In Table 15.5, we show the parameters of only one of them, as well as the scalar DM best fit. The SNR is

$$ \sqrt{\mathrm{SNR}_{A}^{\mathrm{2 GW}}} \approx 251 \qquad (15.42a) $$

$$ \sqrt{\mathrm{SNR}_{E}^{\mathrm{2 GW}}} \approx 120, \qquad (15.42b) $$

while the RNR are

$$ \sqrt{\mathrm{RNR}_{A}^{\mathrm{DM},2 \mathrm{GW}}} \approx 243 \qquad (15.42c) $$

$$ \sqrt{\mathrm{RNR}_{E}^{\mathrm{DM},2 \mathrm{GW}}} \approx 119. \qquad (15.42d) $$

Therefore, it seems the model does not fit correctly the data. As we have not implemented multiple binaries fit, we cannot compare with a model of two GBs.

Therefore, even with such two GB, one should not expect degeneracy between DM and GW signals as seen by LISA. This conclusion is also strengthened by the footnote 13, where we stated that when trying to model GW data by scalar DM, even when relaxing the prior on the velocity, i.e with a best fit velocity close to the speed of light, the DM model still poorly reflects the GW data. This suggests that the main signature which allows one to make the difference between DM and GW is the polarization. This is why combining several GW of rank 2 polarization does not mimic the effect of a scalar DM with a rank 1 polarization.

<table><caption>Table 15.5: On the left column, we show the GW injected parameters in the simulation for the double GW fit. On the right, we show the DM model best fit.</caption><thead><tr><th></th><th>D<sub>2 GW</sub></th><th>M<sup>DM</sup><sub>2 GW</sub></th></tr></thead><tbody><tr><td>ε (10<sup>-4</sup>)</td><td>7.106</td><td>9.446</td></tr><tr><td>f (mHz)</td><td>0.121</td><td>0.121</td></tr><tr><td>β (rad)</td><td>1.046</td><td>0.482</td></tr><tr><td>λ (rad)</td><td>0.431</td><td>3.699</td></tr><tr><td>Φ (rad)</td><td>5.936</td><td>3.290</td></tr><tr><td>v (10<sup>6</sup>) (m/s)</td><td>299.8</td><td>1.149</td></tr><tr><td>ι (rad)</td><td>π/2</td><td>—</td></tr><tr><td>Ψ (rad)</td><td>0</td><td>—</td></tr><tr><td>˙f (aHz/s)</td><td>0</td><td>—</td></tr></tbody></table>

## 15.3 Realistic limit on the sensitivity of LISA to ultralight dark matter couplings

As mentioned in Section 15.2.5, the small velocity of the DM wave implies a large width of the $\varepsilon$ parameter posterior distribution. This broadening of the distribution due to this correlation implies that for the specific realization of the data considered in Section 15.2.5, the sensitivity of the detector to the coupling $\varepsilon$ decreases compared to a situation where $\varepsilon$ is decorrelated from the other parameters (i.e when $v_{DM}$ is fixed). In this section, we derive analytically the width of the $\varepsilon$ distribution, as function of the frequency, and we check for its consistency for the specific realization of the data used in the simulation. Then, in Chapter 18, we will use this analytical result to infer the sensitivity of LISA to ULDM couplings for any frequency in the LISA band, taking into account this correlation.
### 15.3.1 Analytical likelihood

As it was discussed in Section 15.2.3, we use the TDI A, E combinations to fit a model onto our dataset. As shown in Eq. (15.25), the two combinations have the same noise PSD, such that only the signal will differ between the two. In the following, we will work using continuous Fourier transforms (and not DFT), and we will express Fourier transform quantities with a hat. We start by expressing the (complex) Fourier signal of the A combination as a function of parameters ($\epsilon^D, f^D, \beta^D, \lambda^D, v_{DM}^D, \Phi^D$) at frequency $\omega^{D15}$

$$ \hat{d}_{A}^{\mathrm{DM}}(\omega^{D}) = \frac{\sqrt{16\pi G \rho_{\mathrm{DM}} v_{\mathrm{DM}}^{D} \epsilon^{D}}}{\omega^{D} c^{2}} \left( \frac{|\mathcal{T}_{Z}^{\mathrm{DM}}| e^{-i\left(\frac{4\omega_{L}^{D}}{c} + \vec{k}^{D} \cdot \vec{x}_{3} - \Phi^{D}\right)} - |\mathcal{T}_{X}^{\mathrm{DM}}| e^{-i\left(\frac{4\omega_{L}^{D}}{c} + \vec{k}^{D} \cdot \vec{x}_{1} - \Phi^{D}\right)}}{\sqrt{2}} \right) \quad (15.43\mathrm{a}) $$

$$ = \frac{16\sqrt{8\pi G\rho_{\mathrm{DM}}v_{\mathrm{DM}}^{D}\epsilon^{D}}}{\omega^{D}c^{2}} \sin\left(\frac{\omega^{D}L}{c}\right) \sin\left(\frac{2\omega^{D}L}{c}\right) \sin^{2}\left(\frac{\omega^{D}L}{2c}\right) e^{-i\left(\frac{4\omega_{L}^{D}}{c} - \Phi^{D}\right)} \quad (15.43\mathrm{b}) $$

$$ \left( \left( \hat{n}_{12} \cdot \hat{e}_{v}^{D} \right) e^{-i\vec{k}^{D} \cdot \vec{x}_{3}} - \left( \hat{n}_{23} \cdot \hat{e}_{v}^{D} \right) e^{-i\vec{k}^{D} \cdot \vec{x}_{1}} \right) \\ \approx \frac{16\sqrt{8\pi G\rho_{\mathrm{DM}}v_{\mathrm{DM}}^{D}\epsilon^{D}}}{\omega^{D}c^{2}} (\hat{n}_{12} - \hat{n}_{23}) \cdot \hat{e}_{v}^{D} \sin\left(\frac{\omega^{D}L}{c}\right) \sin\left(\frac{2\omega^{D}L}{c}\right) \sin^{2}\left(\frac{\omega^{D}L}{2c}\right) \quad (15.43\mathrm{c}) \\ e^{-i\left(\frac{4\omega_{L}^{D}}{c} - \frac{\omega_{v}^{D}v_{\mathrm{DM}}^{D}|\vec{x}_{\mathrm{AU}}|\cos(\beta^{D})}{c^{2}} - \Phi^{D}\right)} $$

where we have used the transfer function of X Eq. (15.8), which we also use for the Z combination by the simple $\hat{n}_{23} \rightarrow \hat{n}_{12}$ and $\vec{x}_{1} \rightarrow \vec{x}_{3}$ changes. At the last line, we approximated $\vec{k}^{D} \cdot \vec{x}_{1} \approx \vec{k}^{D} \cdot \vec{x}_{\mathrm{AU}} + O(kL)$ (as $kL \ll 1$ for all DM frequencies of interest, as discussed previously), where $|\vec{x}_{\mathrm{AU}}| \approx 1.5 \times 10^{11}$ m is the mean Earth-Sun distance, and similarly for $\vec{k}^{D} \cdot \vec{x}_{3}$. Therefore, both dot products are equal at leading order. In addition, we assume that spacecrafts are close to the ecliptic plane such that the dot product $\hat{k}^{D} \cdot \vec{x}_{\mathrm{AU}}$ evolves with time and can be approximated by $\sim -|\vec{x}_{\mathrm{AU}}|\cos(\beta^{D})\cos(\omega_{E}t - \lambda^{D})$, where $\omega_{E}$ is Earth rotation frequency around the Sun. As $\omega_{E} \ll \omega$ and in particular $\omega_{E}$ is outside of the LISA band, the oscillation can be neglected because it will only produce a small broadening of the Fourier peak at frequency $f = \omega/2\pi$, as mentioned previously. In a similar manner, the signal of the E combination reads

$$ \hat{d}_{E}^{\mathrm{DM}}(\omega^{D}) = \frac{16\sqrt{8\pi G\rho_{\mathrm{DM}}v_{\mathrm{DM}}^{D}\epsilon^{D}}}{\omega^{D}c^{2}} \frac{3\hat{n}_{13} \cdot \hat{e}_{v}^{D}}{\sqrt{3}} \sin\left(\frac{\omega^{D}L}{c}\right) \sin\left(\frac{2\omega^{D}L}{c}\right) \sin^{2}\left(\frac{\omega^{D}L}{2c}\right) \quad (15.43\mathrm{d}) \\ e^{-i\left(\frac{4\omega_{L}^{D}}{c} - \frac{\omega_{v}^{D}v_{\mathrm{DM}}^{D}|\vec{x}_{\mathrm{AU}}|\cos(\beta^{D})}{c^{2}} - \Phi^{D}\right)}, $$

i.e it differs from Eq. (15.43c) just by the geometric factor, at leading order $((\hat{n}_{12} - \hat{n}_{23}) \cdot \hat{e}_{v}^{D})/\sqrt{2} \rightarrow 3(\hat{n}_{13} \cdot \hat{e}_{v}^{D})/\sqrt{6}$. Similarly, in Fourier space, the model for both combinations reads

<footnote n="15">¹⁵The signal is purely monochromatic (up to a small correction, see next paragraph), and therefore we will be interested at the likelihood at a unique Fourier bin, which for simplicity, we assume is exactly the same as the data frequency $f^D = \omega^D/2\pi$.</footnote>
$$ \hat{m}_A^{\text{DM}}(\omega^M) = \frac{16\sqrt{8\pi G\rho_{\text{DM}}v_{\text{DM}}^M\epsilon^M}}{\omega^M c^2} (\hat{n}_{12} - \hat{n}_{23}) \cdot \hat{e}_v^M \sin\left(\frac{\omega^M L}{c}\right) \sin\left(\frac{2\omega^M L}{c}\right) \sin^2\left(\frac{\omega^M L}{2c}\right) e^{-i\left(\frac{4\omega^M L}{c} - \frac{\omega^M \epsilon_{\text{DM}}^M |\vec{x}_{\text{AU}}| \cos(\beta^M)}{c^2}\right) - \Phi^M} \quad (15.44\text{a}) $$

$$ \hat{m}_E^{\text{DM}}(\omega^M) = \frac{16\sqrt{8\pi G\rho_{\text{DM}}v_{\text{DM}}^M\epsilon^M}}{\omega^M c^2} \left(\sqrt{3}\hat{n}_{13} \cdot \hat{e}_v^M\right) \sin\left(\frac{\omega^M L}{c}\right) \sin\left(\frac{2\omega^M L}{c}\right) \sin^2\left(\frac{\omega^M L}{2c}\right) e^{-i\left(\frac{4\omega^M L}{c} - \frac{\omega^M \epsilon_{\text{DM}}^M |\vec{x}_{\text{AU}}| \cos(\beta^M)}{c^2}\right) - \Phi^M}, \quad (15.44\text{b}) $$

i.e they are the same as Eqs. (15.43c) and (15.43d), but where all quantities with M superscript are the unknown parameters.

Finally, we define the standard deviation of the TDI combinations A, E noise at frequency f as

$$ \hat{\sigma}_{A,E}(f) = \sqrt{\int_{f - \frac{1}{2T_{\text{obs}}}}^{f + \frac{1}{2T_{\text{obs}}}} df N_{A,E}(f)} \approx \sqrt{\frac{N_{A,E}(f)}{T_{\text{obs}}}} \equiv \hat{\sigma}(f), \quad (15.45) $$

where the $N_{A,E}$ are given by Eq. (15.25)). At the second to the last equality, we assumed that the noise level is constant over the frequency bandwidth of interest, which is $1/T_{\text{obs}}$ wide, for a monochromatic signal.

Let us now write the (monochromatic) joint Gaussian likelihood in Fourier space as (from Eq. (15.24a))

$$ \log \mathcal{L}(f) = \kappa - \frac{\left| \hat{d}_A^{\text{DM}}(f) - \hat{m}_A^{\text{DM}}(f) \right|^2 + \left| \hat{d}_E^{\text{DM}}(f) - \hat{m}_E^{\text{DM}}(f) \right|^2}{2\hat{\sigma}^2(f)}, \quad (15.46) $$

where f denotes the Fourier frequency and $\kappa$ is a normalization constant independent of the model, and therefore will not play a role in the following analysis. In the following, we will be interested in the joint likelihood of $\epsilon$ and $v_{\text{DM}}$ parameters. In practice, 1) we assume that all the other parameters are individually uncorrelated with the rest (see Fig. 15.5, where one can notice that the parameters that we want to marginalize over are not correlated with any other parameter); and 2) we use a Laplace approximation, which approximates the posterior distribution by a Gaussian around the best fit value, and we approximate the best fit value by the injected value<fnref n="16" />. Therefore, the marginalization consists essentially in taking $(f^M, \beta^M, \lambda^M, \Phi^M) = (f^D, \beta^D, \lambda^D, \Phi^D)$, such that the velocity-coupling joint likelihood

<footnote n="16">¹⁶This is a valid approximation if the data contains no noise, which is the case. This approximation can also be confirmed by Table 15.2.</footnote>
<figure><img src="image_26.png" /><figcaption>Figure 15.8: Joint likelihood $\mathcal{L}(v_{\text{DM}}^M, \epsilon^M)$ for the data parameters shown in Table 15.2.</figcaption></figure>

becomes

$$ \log \mathcal{L}(v_{\text{DM}}^M, \epsilon^M) = - \frac{\left| \hat{d}_A^{\text{DM}}(f) - \hat{m}_A^{\text{DM}}(f) \right|^2 + \left| \hat{d}_E^{\text{DM}}(f) - \hat{m}_E^{\text{DM}}(f) \right|^2}{2\hat{\sigma}^2(f)} \quad (15.47\text{a}) $$

$$ = -\frac{T_{\text{obs}}}{2N_{A,E}(f)} \mu_{\text{DM}} \left( \frac{16\sqrt{8\pi G\rho_{\text{DM}}}}{2\pi f c^2} \right)^2 \sin^2\left(\frac{2\pi f L}{c}\right) \sin^2\left(\frac{4\pi f L}{c}\right) \sin^4\left(\frac{2\pi f L}{2c}\right) \quad (15.47\text{b}) $$

$$ \left( \left(v_{\text{DM}}^D \epsilon^D\right)^2 + \left(v_{\text{DM}}^M \epsilon^M\right)^2 - 2v_{\text{DM}}^D \epsilon^D v_{\text{DM}}^M \epsilon^M \cos\left(\frac{2\pi f |\vec{x}_{\text{AU}}| \cos(\beta)}{c^2} \left(v_{\text{DM}}^D - v_{\text{DM}}^M\right)\right) \right), $$

with

$$ \mu_{\text{DM}} = \left( (\hat{n}_{12} - \hat{n}_{23}) \cdot \hat{e}_v \right)^2 + \left( \sqrt{3}\hat{n}_{13} \cdot \hat{e}_v \right)^2, \quad (15.47\text{c}) $$

where we dropped the first term of Eq. (15.46) since it is a normalization factor, where we consider a monochromatic signal at frequency $f$ and where we dropped the $D$ and $M$ superscripts for the marginalized parameters (i.e when $p^M = p^D$). By using the true value of parameters $\epsilon^D, f, \beta, v_{\text{DM}}^D$ from Table 15.2, one can obtain the joint likelihood as function of $\epsilon^M, v_{\text{DM}}^M$. This is shown in Fig. 15.8, where we set numerically $\mu_{\text{DM}} \sim 4.70^{17}$. We recover the banana-shaped likelihood as in Fig. 15.5.

<footnote n="17">¹⁷This factor is time dependent and evolves with the orbit. Numerically, we can recover the time evolution of the dot product between the DM wind direction and any LISA arm, and we find that $(\hat{n}_{12} \cdot \hat{e}_v) - (\hat{n}_{23} \cdot \hat{e}_v) \sim 1.29+0.43 \cos(\omega_E t + \varphi)$ and $\hat{n}_{13} \cdot \hat{e}_v \sim 0.79-0.03 \cos(\omega_E t + \varphi')$, where $\omega_E, \varphi, \varphi'$ are respectively the Earth rotation frequency around the Sun and two irrelevant phases. As we are interested in frequencies $\omega \gg \omega_E$, we assume that the second term is constant and therefore $(\hat{n}_{12} \cdot \hat{e}_v) - (\hat{n}_{23} \cdot \hat{e}_v) \sim 1.29+0.43 \sim 1.72$ and $\sqrt{3}(\hat{n}_{13} \cdot \hat{e}_v) \sim 1.32$, such that the sum of the squares is $\sim 4.70$.</footnote>
### 15.3.2 Coupling posterior width

From this joint likelihood, we can compute the joint posterior distribution which reads

$$ \mathcal{P}(v_{\text{DM}}^{M}, \varepsilon^{M}) = \frac{\mathcal{L}(v_{\text{DM}}^{M}, \varepsilon^{M})\pi(v_{\text{DM}}^{M})\pi(\varepsilon^{M})}{\int dv_{\text{DM}}^{M} d\varepsilon^{M} \mathcal{L}(v_{\text{DM}}^{M}, \varepsilon^{M})\pi(v_{\text{DM}}^{M})\pi(\varepsilon^{M})}, \quad (15.48) $$

where $\pi(v_{\text{DM}}^{M})$, $\pi(\varepsilon^{M})$ are respectively the prior distribution on the velocity, given by Eq. (5.7b) and the prior on the coupling, which is uniform. The denominator ensures that the posterior distribution integrates to 1 over the whole parameter space. Then, to obtain the posterior distribution of the $\varepsilon$ parameter, we marginalize over $v_{\text{DM}}^{M}$ i.e

$$ \mathcal{P}(\varepsilon^{M}) = \frac{\int dv_{\text{DM}}^{M} \mathcal{L}(v_{\text{DM}}^{M}, \varepsilon^{M}) \pi(v_{\text{DM}}^{M}) \pi(\varepsilon^{M})}{\int dv_{\text{DM}}^{M} d\varepsilon^{M} \mathcal{L}(v_{\text{DM}}^{M}, \varepsilon^{M}) \pi(v_{\text{DM}}^{M}) \pi(\varepsilon^{M})}. \quad (15.49) $$

From this, one can compute the mean $\mu_{\varepsilon}$ and standard deviation $\sigma_{\varepsilon}$ of $\mathcal{P}(\varepsilon^{M})$ by numerical integration. For the parameters that are correctly fitted (i.e the frequency $f$, the ecliptic latitude $\beta$), we use the true values for both data and model parameters (we neglect the small difference between the true values and the best fit values, see Table 15.2), and we do the same for the velocity $v_{\text{DM}}^{D}$ and the coupling $\varepsilon^{D}$. In practice, using *Mathematica*, we numerically integrate $\mathcal{P}(\varepsilon^{M})$ over $\varepsilon^{M}$ with a global adaptative method, where the region of integration is cut in various sub regions in which the numerical integration is performed. The boundaries of the integral correspond to the prior bounds used in the simulation ($[10^{-6}, 10^{-4}]$) and we find

$$ \mu_{\varepsilon} = \int d\varepsilon^{M} \varepsilon^{M} \mathcal{P}(\varepsilon^{M}) \sim 2.69 \times 10^{-5}, \quad (15.50a) $$

$$ \sigma_{\varepsilon} = \sqrt{\int d\varepsilon^{M} (\varepsilon^{M} - \mu_{\varepsilon})^{2} \mathcal{P}(\varepsilon^{M})} \sim 1.17 \times 10^{-5}. \quad (15.50b) $$

As it can be noticed from Table 15.2, the relative error on the mean and standard deviation is ~10% and ~36%, compared to what is found in the simulation. The main reason for this discrepancy is that low DM velocities $0 \le v_{\text{DM}}^{M} \le 1.5 \times 10^{5}$ m/s seem to be disfavoured by data (see Fig. 15.5), while this velocity range is permitted by our prior. As a consequence, "large" values of coupling $\varepsilon^{M}$ are not allowed in the analysis using *Nessai* since one still requires $v_{\text{DM}}^{M} \varepsilon^{M} \sim v_{\text{DM}}^{D} \varepsilon^{D}$ for both amplitudes to match. Therefore, the mean and standard deviation get shifted to smaller values than expected, as we have calculated it above. However, we argue that this velocity cut-off is an artefact of the sampler and is inconsistent with the data that is provided. First, the minimum value of velocities on the posterior distribution given by the simulation depends on the prior that we use for the coupling $\varepsilon^{18}$. Second, there is, to our knowledge, no physical reason why the low velocities would not be permitted by the data. Indeed, as mentioned previously, the velocity is not well constrained by the data because the phase induced by the propagation is too small to be detected, and therefore this induced a cut-off of large velocities. However, any arbitrary small velocity $v_{s}$ should be allowed (as

<footnote n="18">¹⁸We always use a uniform prior distribution, but the minimum and maximum allowed values are somewhat arbitrary.</footnote>
long as $\pi(v_s) \neq 0$ and that we can associate a "large" coupling $\epsilon_l$ such that $v_s \epsilon_l \sim v_{DM}^D \epsilon^D$ again).

In Chapter 18, where we will derive *LISA* sensitivities to scalar ULDM fields, we will use the previous simple semi-analytical method to derive $\sigma_\epsilon$ for all frequencies in *LISA* range.

## 15.4 Conclusion

The conclusions of this chapter are two-fold. The first one is that we managed to answer the question asked at the beginning of the chapter, namely, are GB and scalar ULDM signals degenerate for *LISA*, i.e is it possible to model accurately a GW signal with scalar DM and vice versa ? It seems that it is not the case, i.e through a heliocentric orbit of one year, *LISA* is able to decorrelate the two signals. This is encouraging with regard to the possible detection of DM by *LISA*. The second conclusion is that, as we have seen in Section 15.2.5, the fact that the ULDM field is non relativistic makes it non propagating as seen by *LISA*, especially at low frequency<fnref n="19" />. Therefore, the velocity of the DM wave is essentially invisible for *LISA*. The amplitude of the signal being proportional to the DM-SM coupling, but also to the DM velocity, this means that *LISA* is in reality sensitive to the product of the two parameters, but not to one or the other individually. In other words, if one assumes that the velocity is not fixed, but is a free parameter of the theory, the sensitivity of *LISA* to the coupling will decrease due to its correlation with the velocity. In Chapter 18, we will derive two sensitivity curves for *LISA*: a standard sensitivity curve where the coupling is assumed to be uncorrelated with the rest of the parameters (i.e we assume the velocity is fixed at $v_{DM} \sim 10^{-3}$ c), and a second one, where we take into account those correlations, i.e where we keep the velocity as a free parameter. The former has already been derived for one given dilaton-SM coupling in [19], but the latter, which is more realistic, has never been considered in the literature.

<footnote n="19">¹⁹As we shall see in Chapter 18, even at higher frequencies, when we consider a signal strength such that SNR = 1, the velocity is still badly constrained by data.</footnote>