Extracting electromagnetic signatures of spacetime fluctuations

The theoretical outcome of combining quantum mechanics and general relativity remains one of the open problems of physics. In 1957, Wheeler hypothesised a fluctuating spacetime metric [1] while Wigner and Salecker noted the existence of fundamental limits on the measurement of spacetime distances [2–4]. Conceptual aspects of such fluctuations and the measurability of the gravitational field were also discussed at the 1957 Chapel Hill Conference [5]. However, it was not until 1999 that an experimental route to observe such fluctuations using laser-interferometric gravitational wave detectors was suggested [6]. Complementary to terrestrial detectors, astronomical observations were also suggested for similar purposes [7].

The general strategy has since been to constrain hypothesised theoretical models using data from these detectors and observations [8, table 1]. The former motivated the derivation of a device-independent limit on the minimum length that can be measured [9]. The latter includes a variety of sources whose emissions have been exploited to constrain proposed models [10, 11]. The increasing sensitivity of laser-interferometric gravitational wave detectors has also placed constraints on some of these models. Indeed, the data from dedicated co-located and co-aligned Michelson interferometers have also been used to rule out certain hypothesised non-commutative and holographic models of spacetime [12, 13].

Models based on the hypothesised holographic principle of quantum gravity have, in particular, received some recent attention. While certain ways of invoking this holographic principle have seemingly been ruled out by astronomical observations [11], others have been studied for their signatures in Michelson interferometers [14–16]. These single-parameter models seek to provide a microscopic mechanism for spacetime fluctuations, and their validity is to be probed in forthcoming experiments with Michelson interferometers such as QUEST and GQuEST [17–19].

In this paper, we present a formalism that describes how fluctuations in the spacetime metric manifest as measurable signatures in electromagnetic radiation propagating through it. The fluctuations in the elements of the metric tensor are modelled using their strengths and scales. Signatures of such fluctuations are identified in general electromagnetic field correlation tensors after propagation. We illustrate the strength of our formalism by constraining models of spacetime fluctuations, in two-dimensional exclusion plots, using the experimental data from a single as well as two co-located Michelson interferometers such as the Holometer [20, 21] and aLIGO [22]. The experimental data from these interferometric setups typically capture the fluctuations in the output beam power from its average value. In our approach, we obtain these output power fluctuations using the above-mentioned electromagnetic field correlation tensors, and a subsequent comparison with the experimental data helps in constraining the spacetime fluctuations. We also apply our methods to the proposed QUEST [18] experiment at the Cardiff University [17], and relate our constraints to those obtained using some previous single-parameter models.

In principle, a generic random process could require infinitely many parameters to completely describe its infinitely many moments. However, our two-parameter model is motivated by the empirical observation that physical phenomena may typically be characterised solely by their strength and scale. These are the two parameters captured in our model. A single-parameter model [10, 11, 16] elides these two distinct physical aspects. In fact, the single-parameter models assume the spacetime fluctuations to be white noise, thereby implicitly assuming the scale to be zero [14] or Planck length [10, 11]. Evidently, the single-parameter models lose a degree of freedom by assuming the scale of the fluctuation to be a fixed value. Foregoing these implicit assumptions with our two-parameter model can thus become less restrictive and lead to a more faithful reflection of the physical Universe. Using both the strength and the scale to model the metric fluctuations, in contrast to prior single-parameter models, enables stricter constraints from the cross-correlation between the outputs of the two co-located Michelson interferometers.

We list the mathematical assumptions invoked for analytical tractability, and discuss their consequences. Indeed, invoking a subset of these assumptions even allows these fluctuations to be interpreted as spatial and temporal fluctuations in the refractive index of the medium.

Our work complements prior ones [11, 16, 23] by taking a phenomenological approach to relate spacetime fluctuations to their electromagnetic signatures. Our model is not derived from a fundamental quantum theory of gravity and is agnostic of the underlying mechanism of spacetime fluctuations. It rather captures the scenario where any residual fluctuations in spacetime measurements, after identifying known sources, may be ascribed to possible quantum origins. Operationally, it allows the strength and scale of spacetime fluctuations to be constrained by experimental data. This in turn helps our formalism provide an unambiguous way of judging the performance of experiments designed to constrain spacetime fluctuations and in formulating new experiments. Finally, it should also enable a unified treatment for constraining spacetime fluctuations from terrestrial experiments and astronomical observations.

More generally, our work should enable unification with explorations of other consequences on spacetime fluctuations, particularly on massive [24] and charged objects [25]. The former is hypothesised to lead to spatial decoherence, and has already been constrained using interferometric gravitational wave detectors [26] and parts thereof [27]. The latter suggests the emission of radiation due to fluctuating accelerations, and has also been constrained using terrestrial experiments and astronomical observations [28]. These efforts are part of a larger body of work studying the coupling of quantum matter with gravity [29].

The paper is structured as follows. Section 2 lists the attributes of our model of spacetime fluctuations. Section 3 uses these attributes to simplify the relativistic wave equation for the propagation of light in curved spacetime. Their solution, subject to certain assumptions, provides analytical expressions for the relevant correlation tensors of the electric field. Section 4 specialises the correlation functions to interferometric setups to obtain the spectral densities of the relevant observable in a single, and two co-located interferometers. These expressions are compared with the experimental data in section 5 to obtain two-dimensional exclusion plots constraining the scale and strength of the fluctuations. We use the experimental data from the Holometer experiment and aLIGO, as well as the proposed QUEST experiment. Section 6 relates our results to prior constraints obtained using single-parameter models of spacetime fluctuations. We conclude in section 7.

We propose a phenomenological model to capture fluctuations in the spacetime metric tensor $g^{\mu \nu}.$ We begin by listing the attributes of these fluctuations. Of the ten independent coefficients g^µν that can fluctuate subject to symmetry, we focus on one (g⁰⁰) to illustrate our formalism. Some of these attributes are invoked to capture empirical insights while others to enable a tractable analysis. Each of these can be modified or revoked to further the understanding of possible spacetime fluctuations.

We consider metric fluctuations that are

• (1)  
  of the form,

  $$\begin{align} g^{00} = - n^{2}\left(\boldsymbol{r}\right), \quad g^{ij} = \delta_{ij}, \quad g^{0i} = g^{i0} = 0, \end{align} \tag{ 1 }$$

  for $i,j = 1,2,3$ where ⁴ $n(\boldsymbol{r}),$ $\boldsymbol{r} \equiv (x,y,z,t)$ is a real function of its arguments ⁵ . The quadratic form of g⁰⁰ is not essential, but enables an appealing interpretation in section 3.
• (2)  
  such that $n(x,y,z,t) = n_{\small\textrm{S}}(x,y,z) n_{\small\textrm{T}}(t),$ where $n_{\small\textrm{S}}(x,y,z),n_{\small\textrm{T}}(t)$ are real functions of their arguments. This attribute is crucial for variable separation assumed in section 3.
• (3)  
  small, that is, $n_{\small\textrm{S}}(x,y,z) = 1 + w_{\small\textrm{S}}(x,y,z)$, $n_{\small\textrm{T}}(t) = 1 + w_{\small\textrm{T}}(t)$, and $w_{\small\textrm{S}}, w_{\small\textrm{T}} \ll 1$ are real functions of their arguments. This captures the empirical fact that any fluctuations of spacetime are yet to be detected in the most sensitive of experiments. It is invoked towards the end of section 3.
• (4)  
  zero mean, stationary random processes, that is,

  $$\begin{align} \overline{w_{\small\textrm{X}}} & = 0 \quad \left(\small\textrm{X} = \small\textrm{S,T}\right), \end{align} \tag{ 2 }$$

  $$\begin{align} \overline{w_{\small\textrm{S}}\left(x,y,z\right) w_{\small\textrm{S}}\left(x+\delta_{x},y+\delta_{y},z+\delta_{z}\right)} & = \Gamma_{\small\textrm{S}} \, \rho\left(\delta_{x},\delta_{y},\delta_{z}\right), \end{align} \tag{ 3 }$$

  and

  $$\begin{align} \overline{w_{\small\textrm{T}}\left(t\right) w_{\small\textrm{T}}\left(t+\delta_{t}\right)} & = \Gamma_{\small\textrm{T}} \, \varrho\left(\delta_{t}\right), \end{align} \tag{ 4 }$$

  where $\overline{\mathcal{O}}$ denotes the average of any quantity $\mathcal{O}$ over realisations of the random process.
• (5)  
  uncorrelated in space and time, that is, the average

  $$\begin{align} \overline{w_{\small\textrm{S}}\left(x,y,z\right)\,w_{\small\textrm{T}}\left(t\right)} & = \overline{w_{\small\textrm{S}}\left(x,y,z\right)} \: \: \overline{w_{\small\textrm{T}}\left(t\right)}. \end{align} \tag{ 5 }$$

  Attributes (4) and (5) enable analytic computation of the relevant observables in section 5.
• (6)  
  Gaussian in space, that is,

  $$\begin{align} \rho\left(\delta_{x},\delta_{y},\delta_{z}\right) & = \text{exp} \left( -\sum\limits_{i = x,y,z} \dfrac{\delta_{i}^{2}}{2 \ell_{i}^{2}} \right), \end{align} \tag{ 6 }$$

  with correlation lengths $\ell_{i}, i = x,y,z.$ This attribute allows analytical closed-form expressions of relevant observables such as the power spectral and cross-spectral densities in section 5. We leave the form of $\varrho(\delta_{t})$ unspecified, noting that it may have a characteristic time scale τ_t .

Our model is thus parameterised only by the strength $\Gamma_{\small\textrm{S}}$ and scale $N_i \equiv 2\ell_i/\ell_{\small\textrm{P}}$ of the fluctuations, where the Planck length $\ell_{\small\textrm{P}}$ is introduced to enable a fully dimensionless parameterisation. The factor of 2 in the parametrisation of N_i is chosen in anticipation of section 5 where light propagates along the two arms of a Michelson interferometer oriented along the x and z axes, and the relevant observables become proportional to $(\ell_x + \ell_z)$. The observables will then be proportional to the chosen scale, especially in an isotropic model where $N_{\parallel} \equiv N_{x} = N_{y} = N_{z}$.

In what follows, we first describe how the spacetime fluctuations affect the propagation of electromagnetic fields. Using data from optical experiments, we then set upper bounds on the strength and scale of the spacetime fluctuations.

The relativistic wave equation governing the propagation of light in curved spacetime is [30, 31]

$$\begin{align} \Box F_{\alpha \beta} + 2 R_{\alpha \gamma \beta \delta} F^{\gamma \delta} - R_{\alpha \gamma} F^{\gamma}_{\beta} + R_{\beta \gamma} F^{\gamma}_{\alpha} = 0, \end{align} \tag{ 7 }$$

where F_αβ is the electromagnetic field tensor, $\Box$ denotes the d'Alembertian operator with the covariant derivative, $R_{\alpha \beta \gamma \delta}$ is the Riemann tensor, and R_αβ is the Ricci curvature tensor defined in terms of the spacetime metric g^µν . Equation (7) is obtained from Maxwell's equations written in terms of the covariant derivative without further assumptions of gauge or form of the 4-potential [31].

In principle, all observable signatures of the propagation of electromagnetic fields through the fluctuating spacetime, can be extracted by solving equation (7) for a general tensor $g^{\mu\nu}.$ In practice, this seems intractable, even when the metric fluctuations are restricted to possess the six attributes listed in section 2. Further assumptions must thus be made to solve equation (7).

The Christoffel symbols contained in equation (7) involve derivatives of the metric tensor $g^{\mu\nu}.$ In particular, it contains derivatives of $g^{00} = - n^{2}(\boldsymbol{r})$, as per Attribute (1) of our model in section 2. The corresponding wave equation, provided in appendix A, eluded our attempts at a closed-form analytical solution. To advance we thus assume that for a field tensor F_αβ ,

$$\begin{align} \mathcal{O}\left( \partial_{\mu} n^{2} \right) &\ll \frac{\mathcal{O}\left( \partial_{\nu} F_{\alpha \beta} \right) }{\mathcal{O}\left( F_{\varphi \gamma} \right)}, \end{align} \tag{ 8a }$$

$$\begin{align} \quad \mathcal{O}\left( \partial_{\mu} n^{2} \right) &\ll \frac{\mathcal{O}\left( \partial_{\chi} \partial_{\nu} F_{\alpha \beta} \right)}{\mathcal{O}\left( \partial_{\sigma} F_{\varphi \gamma} \right)}, \quad \end{align} \tag{ 8b }$$

$$\begin{align} \mathcal{O}\left( \partial_{\mu} \partial_{\nu} n^{2} \right) &\ll \frac{\mathcal{O}\left( \partial_{\chi} \partial_{\sigma} F_{\alpha \beta} \right)}{\mathcal{O}\left( F_{\varphi \gamma} \right)}, \end{align} \tag{ 8c }$$

where $\mathcal{O}\left( \cdot \right)$ denotes an order-of-magnitude of the tensor quantity within the parentheses. These conditions suppose that any variation in n² across spacetime occurs at a rate much slower than the variation of the electromagnetic field. Therefore, in subsequent analysis, the wavelength and time period of the wave, that set the natural scales of the electromagnetic field variations, should be much smaller than the corresponding scales for variations ⁶ in n².

The inequalities in (8a) imply that the terms in equation (7) comprising either the Christoffel symbols or any derivative thereof are negligible compared to $g^{\gamma \delta} \partial_{\delta} \partial_{\gamma} F_{\alpha \beta}$. This simplifies equation (7) to

$$\begin{align} g^{\gamma \delta} \partial_{\delta} \partial_{\gamma} F_{\alpha \beta} = 0, \end{align} \tag{ 9 }$$

which reduces to

$$\begin{align} \nabla^{2} \mathbf{E} & = \frac{n^{2}\left(x,y,z,t\right)}{c^{2}} \, \frac{\partial^{2} \mathbf{E}}{\partial t^{2}} \quad \text{for } \alpha = 0, \, \beta = i \text{ or } \alpha = i, \, \beta = 0, \end{align} \tag{ 10a }$$

$$\begin{align} \nabla^{2} \mathbf{B} & = \frac{n^{2}\left(x,y,z,t\right)}{c^{2}} \, \frac{\partial^{2} \mathbf{B}}{\partial t^{2}}\quad \text{for } \alpha = i,\, \beta = j \neq i \,\left(i,j = 1,2,3\right), \end{align} \tag{ 10b }$$

where E and B denote 3-vector electric and magnetic fields. Using the Attribute (1) and the inequalities (8a), the fluctuating metric tensor thus appears as an effective refractive index for the electromagnetic wave equation in free, flat space. It is also evident from equation (10a) that the choice of the metric in Attribute (1) is crucial for interpreting $n(\boldsymbol{r})$ as a refractive index fluctuating randomly in space and time. Propagation of light through a medium with homogeneous isotropic fluctuations of the refractive index was studied in [32]. The stochastic nature of the random medium is shown to manifest itself in the average of the fourth-order correlation of the radiation field over an ensemble of possible fluctuations of the refractive index.

Motivated by [32], we examine general electric field correlation tensors for observable signatures of the underlying spacetime fluctuations. Electric field correlation tensors can be of the form [33, 34]

$$\begin{align} M^{m,m^{{\prime}}}_{\boldsymbol{c};\boldsymbol{c^{{\prime}}}} \left(\left\{\boldsymbol{r}\right\}; \left\{\boldsymbol{r^{{\prime}}}\right\} \right) = E_{c_1}\left(\boldsymbol{r}_1\right) E_{c_2}\left(\boldsymbol{r}_2\right) \cdots E_{c_m}\left(\boldsymbol{r}_m\right) E_{c^{{\prime}}_1}^*\left(\boldsymbol{r^{{\prime}}}_1\right) E_{c^{{\prime}}_2}^*\left(\boldsymbol{r^{{\prime}}}_2\right) \cdots E_{c^{{\prime}}_{m^{{\prime}}}}^*\left(\boldsymbol{r^{{\prime}}}_{m^{{\prime}}}\right) \end{align} \tag{ 11 }$$

where $c_i, c^{{\prime}}_j \in \{x,y,z\}$ for $i = 1,2,\dots,m$, $j = 1,2,\dots,m^{{\prime}}$, with m not necessarily equal to mʹ and $\{\boldsymbol{r}\} \equiv \{\boldsymbol{r}_1,\cdots \boldsymbol{r}_m\},$ $\boldsymbol{c} \equiv \{c_1,\cdots c_m\}.$ The primed variables are denoted similarly. Here $(m,m^{{\prime}})$ denote the order of the electric field correlation tensor. We choose to focus on the electric field for the rest of this paper. From equation (10a), it follows that the magnetic field will carry equivalent signatures.

The correlation tensors are crucial in studying the different orders of optical coherence of the light beam. The theory of partial coherence becomes relevant specifically when investigating the unknown stochastic nature of the intervening medium. In this paper, we study intensity–intensity correlations captured in an interferometer with classical sources of light. If and when captured, the higher-order correlation tensors could potentially reveal more about the nature of the spacetime fluctuations and may lead to stricter bounds. Another strength of these correlation tensors lies in their natural extendability to quantum sources of light [35, 36] such as entangled and squeezed states. These could also lead to newer and stricter bounds from future experiments using such quantum states [37–39].

Solving equation (10a) for any $n(\boldsymbol{r})$ to obtain the above field correlations in closed form analytically seems intractable, although numerical solutions may be possible. Seeking a closed form solution of equation (10a), we employ the following assumptions motivated by an electromagnetic wave propagating along a certain direction, say z.

• (i)  
  Separation of variables for the electric field. This relies on Attribute (2) of our model in section 2.
• (ii)  
  Paraxial and slowly varying envelope approximation (SVEA) [40].

Using both the ansatz, $E_{y}(x,y,z,t) = \mathcal{S}(x,y,z) \mathcal{T}(t)$ ($i = x,y,z$) and the Attribute (2), we invoke Assumption (i). This separation of variables splits equation (10a) into the Helmholtz-type equations,

$$\begin{align} \nabla^{2} \mathcal{S} = - k^{2} n_{\small\textrm{S}}^{2} \mathcal{S}, \end{align} \tag{ 12 }$$

$$\begin{align} \frac{\partial^{2} \mathcal{T}}{\partial t^{2}} = - \frac{k^{2} c^{2}}{n_{\small\textrm{T}}^{2}} \mathcal{T}, \end{align} \tag{ 13 }$$

with $k = 2 \pi/\lambda = \Omega/c,$. where λ is the wavelength of light and Ω the corresponding frequency. The spatial part of the field $\mathcal{S}(x,y,z)$ (respectively, the temporal part, $\mathcal{T}(t)$) is any function that satisfies equation (12) (respectively, equation (13)).

Subsequently using $\mathcal{S}(x,y,z) = E_{\small\textrm{S}}(x,y,z) e^{-i k z}$ and $\mathcal{T}(t) = E_{\small\textrm{T}}(t) e^{i \Omega t}$, and applying (ii), that is,

$$\begin{align} k \frac{\partial E_{\small\textrm{S}}}{\partial z} \gg \frac{\partial^{2} E_{\small\textrm{S}}}{\partial z^{2}}, \quad \mathrm{and} \quad \Omega \frac{\partial E_{\small\textrm{T}}}{\partial t} \gg \frac{\partial^{2} E_{\small\textrm{T}}}{\partial t^{2}}. \end{align} \tag{ 14 }$$

Equation (10a) simplifies to

$$\begin{align} \frac{\partial E_{\small\textrm{S}}}{\partial z} & = -\frac{i}{2 k} \left( \frac{\partial^{2} E_{\small\textrm{S}}}{\partial x^{2}} + \frac{\partial^{2} E_{\small\textrm{S}}}{\partial y^{2}} \right)- \frac{i k}{2} \left(n_{\small\textrm{S}}^{2}-1\right) E_{\small\textrm{S}}, \end{align} \tag{ 15a }$$

$$\begin{align} \frac{\partial E_{\small\textrm{T}}}{\partial t} & = - \frac{i \Omega}{2} E_{\small\textrm{T}} + \frac{i k^{2} c^{2}}{2 \Omega n_{\small\textrm{T}}^{2}} E_{\small\textrm{T}}. \end{align} \tag{ 15b }$$

Now invoking Attribute (3) and neglecting terms of order $w_{\small\textrm{X}}^{2}$ ($\small\textrm{X} = \small\textrm{S,T}$) and above, the corresponding solutions are

$$\begin{align} E_{\small\textrm{S}}\left(x,y,z\right) & = \iint_{-\infty}^{\infty} \text{d} k_{x} \text{d} k_{y} e^{-i \left(k_{x} x + k_{y} y\right)} \widetilde{E}_{\small\textrm{S}}\left(k_{x},k_{y},z\right), \end{align} \tag{ 16 }$$

$$\begin{align} E_{\small\textrm{T}}\left(t\right) & = E_{\small\textrm{T}}\left(0\right) \, \text{exp} \left( - i\, \Omega\, \int\limits_{0}^{t} \text{d} t^{{\prime}} \, w_{\small\textrm{T}}\left(t^{{\prime}}\right) \right), \end{align} \tag{ 17 }$$

where

$$\begin{align} \widetilde{E}_{\small\textrm{S}}\left(k_{x},k_{y},z\right)& = \widetilde{E}_{\small\textrm{S}}\left(k_{x},k_{y},0\right) \, e^{i z \left( k_{x}^{2} + k_{y}^{2} \right)/2 k} \text{exp} \left(- i k \int\limits_{0}^{z} \text{d} z^{{\prime}} \, w_{\small\textrm{S}}\left(x,y,z^{{\prime}}\right) \right), \end{align} \tag{ 18 }$$

and $k_x$, $k_y$ are the Fourier-conjugate variables corresponding to $x$ and $y$ respectively.

We note that the spatial component of the propagating field captures the cumulative fluctuations along the direction of traversal via $\smallint\nolimits_{0}^{z} \text{d} z^{{\prime}} \, w_{\small\textrm{S}}(x,y,z^{{\prime}})$ in the phase factor of equation (16). This remains unchanged even on implementing the inverse Fourier transform (equation (18)). This effective phase fluctuation with the $x,y$ dependence modifies the shape of the wavefront at each point as the light traverses along the z-axis. The temporal component also has effective phase fluctuation as in equation (17). As the fluctuations appear in the phase, it follows from equation (17) that some correlation tensors given by equation (11) do not contain their signatures, such as $M^{m,m}_{\boldsymbol{c}; \boldsymbol{c}} \left( \{\boldsymbol{r} \}; \{\boldsymbol{r} \} \right)$ which is an instantaneous intensity correlation function.

Others, such as the second-order correlation function

$$\begin{align} M^{1,1}_{y;y} \left(\boldsymbol{r}_1; \boldsymbol{r^{{\prime}}}_1 \right) = E_y\left(\boldsymbol{r}_1\right)E^*_y\left(\boldsymbol{r^{{\prime}}}_1\right) \end{align} \tag{ 19 }$$

do contain signatures of spacetime fluctuations. Indeed, any measurement that records some aspect of the wavefront or reconstructs it entirely should reveal observable signatures of the underlying spacetime fluctuations.

A typical measurement can involve measuring the second-order correlation function $M^{1,1}_{y;y} \left(\boldsymbol{r}_1; \boldsymbol{r^{{\prime}}}_1 \right)$ over a cross-section of a beam at $z = z_{\mathrm{det}}$ and time $t = t_{\mathrm{det}}$, such that $\boldsymbol{r}_1 = (x_{1},y_{1},z_{\mathrm{det}},t_{\mathrm{det}})$ and $\boldsymbol{r^{{\prime}}}_1 = (x_{2},y_{2},z_{\mathrm{det}},t_{\mathrm{det}})$. In this case, the temporal part of the correlation function, $E_{\small\textrm{T}}(t_{\mathrm{det}}) E^{\ast}_{\small\textrm{T}}(t_{\mathrm{det}})$, does not contain any signatures of spacetime fluctuations. Only the spatial part reveals signatures of interest. This is why it was not necessary to specify $\varrho(\delta_t)$ in Attribute (6). This argument is applicable even when considering the appropriate correlation tensor in the context of interferometric setups in subsequent sections.

There is a great variety of possible configurations of the field correlation tensors in equation (11) whose measurements can reveal signatures of spacetime fluctuations. We leave a fuller exploration for the future. In the rest of this paper, we focus on how interferometric measurements can identify signatures of spacetime fluctuations [6].

Given their extremely high sensitivities, laser interferometers have been proposed and used to probe spacetime fluctuations [6, 16, 17, 20]. As shown in figure 1(a), the power of the light at the output port B of such an interferometer, with the beam cross-section in the y-z plane, is measured by the photodetector. The electric field at the output is given by

$$\begin{align} E_{\mathrm{out}}\left(\boldsymbol{r}_0\right) = \frac{1}{\sqrt{2}} \left[ E^{\left(\small\textrm{C}\right)}_{y}\left(\boldsymbol{r}_0\right) - E^{\left(\small\textrm{D}\right)}_{y}\left(\boldsymbol{r}_0\right) e^{-2 i \varphi_{\text{off}}} \right], \end{align} \tag{ 20 }$$

where the superscripts $(\small\textrm{C})$ and $(\small\textrm{D})$ denote the y-component of the electric field along the respective arms of the interferometer, $\boldsymbol{r}_0 \equiv \{0,y,z,\tau_{0} \},$ and τ₀ is the time taken to complete a round trip in either of the arms. We note that $\varphi_{\mathrm{off}}$ is the phase difference due to the offset in the interferometer arm lengths between the two arms of each interferometer. As the arms are almost of equal length $\mathcal{L}$ with $\varphi_{\text{off}}/k \ll \mathcal{L}$, $\tau_{0} = 2 \mathcal{L} / c$. This deliberate small phase offset is present in interferometers that use the so-called DC-readout scheme (including all interferometers considered here), which makes the destructive interference at the output imperfect [41].

Zoom In Zoom Out Reset image size

For brevity, we define the correlation function $M_{\mathrm{out}}^{m,m^{{\prime}}} \left(\{\boldsymbol{r}\}; \{\boldsymbol{r^{{\prime}}}\} \right)$ for the output field of an interferometer in line with $M^{m,m^{{\prime}}}_{\boldsymbol{c};\boldsymbol{c^{{\prime}}}} \left(\{\boldsymbol{r}\}; \{\boldsymbol{r^{{\prime}}}\} \right)$, except with $E_{\mathrm{out}}$ replacing the field components E_c and $E_{c^{{\prime}}}$ in equation (11) ⁷ . In particular, the fourth-order correlation function $M_{\mathrm{out}}^{2,2} \left(\boldsymbol{R}; \boldsymbol{R} \right),$ $\boldsymbol{R} \equiv \{\boldsymbol{r}_{1},\boldsymbol{r}_{2}\}$ with $\boldsymbol{r}_{1} = (0,y_{1},z_{1},\tau_{1})$ and $\boldsymbol{r}_{2} = (\delta_{s},y_{2},z_{2},\tau_{2})$ will be of interest both for a single interferometer as in figure 1(a) as well as two co-located interferometers [17, 20] as in figure 1(b).

The time-averaged output power of a single interferometer measured at time τ is then

$$\begin{align} P_{\mathrm{out}}\left(\tau\right) = \frac{\epsilon_{0} c}{\Delta_{T}} \int\limits_{\tau}^{\tau+\Delta_T} \mathrm{d}\tau_{1} \int_{A} \mathrm{d}s \, M^{1,1}_{\mathrm{out}} \left(\boldsymbol{r}_0; \boldsymbol{r}_0 \right), \nonumber \end{align} \tag{ 21 }$$

where $M^{1,1}_{\mathrm{out}} \left(\boldsymbol{r}_0; \boldsymbol{r}_0 \right) = |E_{\mathrm{out}}(\boldsymbol{r}_0)|^{2},$ $\mathrm{d}s \equiv \mathrm{d}y \: \mathrm{d}z,$ $\Delta_{T}$ is integration time of the detector, and $\int_{A}$ denotes a surface integral over the beam cross-section. Similarly, the time-averaged power correlation is

$$\begin{align} P_{\mathrm{out}}^{i}\left(\tau\right) P_{\mathrm{out}}^{j}\left(\tau+\delta_{\tau}\right) = \left(\frac{\epsilon_{0} c}{\Delta_T} \right)^{2} \iint\limits_{\tau,\tau+ \delta_{\tau}}^{\tau+\Delta_T,\tau+ \delta_{\tau}+ \Delta_T} \!\!\! \! \mathrm{d}\tau_{1} \: \mathrm{d}\tau_{2} \iint\limits_A \!\! \mathrm{d}s_{1} \: \mathrm{d}s_{2} \: M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right), \end{align} \tag{ 22 }$$

where $i,j = \text{I,II}$ denote the interferometers which determine the different $E_{\mathrm{out}}$ fields appearing in $M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)\!\!.$ An explicit expression is given in appendix B, equation (B2).

In interferometric setups such as the Holometer, the control systems carrying out measurements monitor fluctuations in the output light and feed back the differential signals to actuated end mirror mounts aimed to hold the average output power constant [20]. The measurement of relevance is thus a time-series of the deviation from the average output power accumulated over periods of $T_{\text{acc}} = 1.8/700~\mathrm{s} = 2.6$ ms [20]. Within each $T_{\text{acc}},$ the average output power itself is obtained by averaging over $N_{\text{sp}} = T_{\text{acc}}/\Delta_{T}$ time measurements.

Equations (21) and (22) capture the experimental scenario where $n(x,y,z,t)$ (equivalently, $w_{\small\textrm{S}}$ and $w_{\small\textrm{T}}$) is sampled as the light propagates in space. However, our theoretical formalism requires assuming specific characteristics of the distribution of $n(x,y,z,t)$ to obtain analytical expressions. While the former can be used to obtain the necessary temporal averages, the latter can only be used to obtain ensemble averages. We thus use the stationarity of the time-series of the deviation from the average output power, to replace time averages with ensemble averages. This can be justified by the tests establishing the stationarity of the Holometer/aLIGO data which is necessary for the ergodicity of a process.

As the recorded data from experiments such as the Holometer are the power spectral density (PSD) from a single interferometer's output or the cross-spectral density (CSD) from the outputs of two co-located interferometers [17, 20, 22], we invoke ergodicity for the correlation functions of the fourth-order in the electric fields, leading to

$$\begin{align} \overline{P_{\mathrm{out}}^{i}\left(\tau\right) P_{\mathrm{out}}^{j}\left(\tau+\delta_{\tau}\right)} = \left(\epsilon_{0} c \right)^{2} \iint\limits_A \!\! \mathrm{d}s_{1} \: \mathrm{d}s_{2} \: \overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)}, \end{align} \tag{ 23 }$$

and the covariance

$$\begin{align} \mathrm{Cov}_{i,j}\left(P_{\mathrm{out}}\right) & = \overline{P^{i}_{\mathrm{out}}\left(\tau\right) P^{j}_{\mathrm{out}}\left(\tau+\delta_{\tau}\right)} -\overline{P^{i}_{\mathrm{out}}\left(\tau\right)} \:\:\: \overline{P^{j}_{\mathrm{out}}\left(\tau+\delta_{\tau}\right)}. \end{align} \tag{ 24 }$$

Note that $\mathrm{Cov}_{i,j}(P_{\mathrm{out}}) \neq \mathrm{Cov}_{j,i}(P_{\mathrm{out}})$, in general ⁸ , for i ≠ j. Extending ergodicity to higher-order moments will require commensurate tests of stationarity on corresponding data.

The PSD is thus

$$\begin{align} S\left(\,f\,\right) = \frac{1}{\pi} \int_{0}^{\infty} \text{d} \delta_{\tau} \, \mathrm{Cov}\left(\Delta x\right) \text{cos} 2 \pi f \delta_{\tau}. \end{align} \tag{ 25 }$$

and the CSD is

$$\begin{align} CS\left(\,f\,\right) = &\dfrac{1}{2 \pi} \int_{0}^{\infty} \text{d} \delta_{\tau} \, \left\{ e^{-2 \pi i f \delta_{\tau}} \, \mathrm{Cov}_{\small\textrm{I,II}}\left(\Delta x\right) + e^{2 \pi i f \delta_{\tau}} \, \mathrm{Cov}_{\small\textrm{II,I}}\left(\Delta x\right) \right\}, \end{align} \tag{ 26 }$$

where ⁹

$$\begin{align} \mathrm{Cov}_{i,j}\left(\Delta x\right) = \left(\frac{\lambda}{4 \pi \varphi_{\mathrm{off}} P_{0}}\right)^{2} \mathrm{Cov}_{i,j}\left(P_{\mathrm{out}}\right). \end{align} \tag{ 27 }$$

For brevity, $\mathrm{Cov}(\Delta x) \equiv \mathrm{Cov}_{i,i}(\Delta x)$.

In practice, the PSD and the CSD are obtained using discrete Fourier transforms of time series of the measured path difference. This corresponds to using discretized transforms in place of equations (25) and (26). In the next section, we compare the PSDs and CSDs obtained analytically by evaluating these discretized versions of equations (25) and (26) with the corresponding measured spectral densities. These analytical evaluations use Attributes (4), (5), and (6) to provide constraints on the strength ($\Gamma_{\small\textrm{S}}$) and scale ($N_{\parallel}$) of the spacetime fluctuations.

In the Holometer and QUEST experiments, the setup is a simple Michelson interferometer as represented in figure 1(b).

We show how constraints on the strength and scale of the spacetime fluctuations could be established from the PSD and then the CSD, both from completed and proposed experiments such as the Holometer [21], aLIGO [22], and QUEST [17]. These results rely on Attributes (4), (5), and (6), and assume the spacetime fluctuations to be isotropic in space with the scale $N_{\parallel} \equiv N_{x} = N_{y} = N_{z}$. We present the constraints as two-dimensional $\Gamma_{\small\textrm{S}}$-$N_{\parallel}$ exclusion plots. Importantly, these constraints only apply when spectral density produced by the spacetime fluctuations is present at a frequency at which the interferometer is sensitive. As explained below, this is not the case for our model from section 2 and presently operational interferometers.

5.1. Constraint from PSD

Constraints based on the Holometer and the QUEST experiments are presented in the first two subsections. The former is based on experimental data [21, figure 2], and the latter on proposed experimental parameters [17]. The third subsection reports similar constraints from aLIGO data, which is a Michelson interferometer with cavities in each arm which requires some additional analysis.

5.1.1. Holometer.

Using equation (B17) from appendix B in the expression for $\mathrm{Cov}_{i,j}(\Delta x)$ in equation (27), and setting $\delta_{s} = 0$ in the discretized version of equation (25), the PSD is

$$\begin{align} S\left(\,f\,\right) = \mathcal{A} \left[ 1 + \, \mathcal{B}^2 \text{cos} 4 \varphi_{\mathrm{off}} - 2 \mathcal{B} \text{cos}^{2} 2 \varphi_{\mathrm{off}} \right], \end{align} \tag{ 28 }$$

where

$$\begin{align} &\mathcal{A} = \frac{1}{4 \pi}\left(\frac{\lambda}{8 \pi \varphi_{\mathrm{off}}}\right)^{2} \, \Theta\left(\,f\,\right), \end{align} \tag{ 29 }$$

$$\begin{align} & \mathcal{B} = \text{exp}{\left(- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right)}, \end{align} \tag{ 30 }$$

$$\begin{align} &\mathcal{K} = k^{2} \ell_{\small\textrm{P}} \Gamma_{\small\textrm{S}} N_{\parallel}, \end{align} \tag{ 31 }$$

and the step function $\Theta(f)$ (in units of time) is

$$\begin{align} \Theta\left(\,f\,\right) = \begin{cases} 1 \,\text{s}, & \mathrm{if }\:\: 0 \unicode{x2A7D} f < f_{\mathrm{step}},\\ 0 \,\text{s}, & \mathrm{otherwise,} \end{cases} \end{align} \tag{ 32 }$$

where $f_{\mathrm{step}}$ denotes the frequency resolution of the discrete Fourier transform, i.e. the width of the frequency bins. Thus the signal has all its power in the zero frequency (DC) bin.

To first-order in Planck length, i.e. $\mathcal{B} = 1- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L}$, the PSD becomes

$$\begin{align} S\left(\,f\,\right) \approx \left(\frac{\lambda}{2 \pi} \right)^{2} \frac{\mathcal{K} \mathcal{L}}{2 \sqrt{2 \pi}} \Theta\left(\,f\,\right). \end{align} \tag{ 33 }$$

Here we have used the fact that $\varphi_{\mathrm{off}}\ll 1$, as in the setups considered $\varphi_{\mathrm{off}}\approx 10^{-5}$ rad (see table 1).

Table 1. Values of the system parameters as measured in the Holometer [21] and aLIGO [22]. Values of $\Delta_{T}$ and $\varphi_{\mathrm{off}}$ are up to an order-of-magnitude. System parameter values used in the simulation of the QUEST experiment [17] are also quoted.

Parameter	Holometer [21]	QUEST [17]	aLIGO [22]
P0 (kW)	2	10	750
$\mathcal{L}$ (m)	40	3	$4 \times 10^{3}$
$\Delta_{T}$ (s)	10−8	10−9	10−5
$\varphi_{\mathrm{off}}$ (rad)	$5\times 10^{-5}$	$5\times 10^{-5}$	$2.5 \times 10^{-5}$
λ ($\times 10^{-6}$ m)	1.064	1.064	1.064
PSD (m$^{2}\,$Hz−1)	$6.3\times 10^{-37}$	$3.2\times 10^{-38}$	$2.3\times 10^{-40}$

In our drastically simplified spacetime fluctuations model, the use of Attribute (2) and Assumption (i) results in an $M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)$ independent of time due to equation (17). This in turn leads to a covariance that is independent of δ_τ and thereby, a PSD with a Dirac delta function of the frequency centred on f = 0. For a discrete cosine transform, $\delta(f)$ is replaced by $\Theta(f).$ This predicted frequency dependence is very different from spectra theorised in other models for quantum spacetime fluctuations [23, 42].

Interferometers have limited or no sensitivity to signals at zero frequency. Therefore, a direct comparison between our predicted signal and interferometric observations is not possible, and our model in section 2 cannot be constrained by existing observations. Nevertheless, to illustrate the principle, we proceed to constrain our model by comparing the order of magnitude of the predicted signal spectral density and the noise of the interferometers, anticipating the possibility of the signal having a different frequency dependence but a similar magnitude.

The sensitive band of the Holometer data (i.e. at frequencies above $f_\mathrm{sens}\approx1$ MHz is limited by photon shot noise [21]—a white noise spectrum in principle [43], which is given by [17]

$$\begin{align} S_{\mathrm{meas}}\left(\,f>1\;\mathrm{MHz}\right) = \frac{\hbar c \lambda}{8 \pi P_{0}}. \end{align} \tag{ 34 }$$

We now consider spacetime fluctuations described by our model to produce a displacement spectral density at frequencies within this band but with a magnitude smaller than that of the shot noise. The model can thus be constrained at a signal-to-noise ratio of one ¹⁰ by setting $S_{\mathrm{meas}}(f_\mathrm{sens}\gt0) \unicode{x2A7E} S(f = 0)$, which gives

$$\begin{align} \Gamma_{\small\textrm{S}} N_{\parallel} \unicode{x2A7D} \frac{\hbar c }{2\sqrt{2 \pi} \, \ell_{\small\textrm{P}}} \: \frac{\lambda}{P_{0} \mathcal{L}}, \end{align} \tag{ 35 }$$

We use equation (35) and the parameter values listed in table 1 in the case of the Holometer, to obtain the allowed parameter space as shown in figure 3(a).

Furthermore, (35) and the subsequent inequalities are valid only between the limits $N_{\parallel} \ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$ and $N_{\parallel} \gg \lambda/\ell_{\small\textrm{P}}$. We note that the latter limit arises from equation (8a) as mentioned in section 3. We consider the case $N_{\parallel} = 2 \mathcal{L}/\ell_{\small\textrm{P}}$ in appendix C for the Holometer. While a constraint on $\Gamma_{\small\textrm{S}}$ is obtained in this case too, the product $\Gamma_{\small\textrm{S}} N_{\parallel} \unicode{x2A7D} 10^{-1}$. This is a weaker constraint than $\Gamma_{\small\textrm{S}} N_{\parallel} \unicode{x2A7D} 5 \times 10^{-3}$ obtained using equation (35) at $N_{\parallel}\ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$. On the other hand, the constraint from CSD becomes stronger than the corresponding constraint for $N_{\parallel}\ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$.

5.1.2. QUEST.

The simulation for QUEST considered a shot-noise limited sensitivity above $f_\mathrm{sens}\approx1$ MHz, and the use of squeezed vacuum states of light, aiming to attain 6 dB of squeezing in the output [17]. This leads to a measured PSD $S_{\small\textrm{QUEST}}(f)$, which is written in terms of $S_{\mathrm{meas}}(f)$ as,

$$\begin{align} S_{\small\textrm{QUEST}}\left(\,f>1\;\mathrm{MHz}\right) = 10^{-0.6} S_{\mathrm{meas}}\left(\,f>1\;\mathrm{MHz}\right). \end{align} \tag{ 36 }$$

As before, the model can thus be constrained at a signal-to-noise ratio of one by setting $S_{\mathrm{meas}}(f_\mathrm{sens}\gt0) \unicode{x2A7E} S(f = 0)$

$$\begin{align} \Gamma_{\small\textrm{S}} N_{\parallel} \unicode{x2A7D} 10^{-0.6} \, \frac{\hbar c }{2\sqrt{2 \pi} \, \ell_{\small\textrm{P}}} \: \frac{\lambda}{P_{0} \mathcal{L}}. \end{align} \tag{ 37 }$$

Using parameters listed in table 1, equation (37) yields the maximum limit on the product $\Gamma_{\small\textrm{S}} N_{\parallel}$ which is plotted in figure 3(a) as a line, indicating that an experimental run of the QUEST experiment could rule out the parameter space above this line. This constraint from the simulation is stricter than that from the Holometer when 6 dB of squeezing is used to achieve sub-shot-noise performance.

5.1.3. aLIGO.

aLIGO is a Michelson interferometer with cavities in each arm [22] as shown in figure 2. These arm cavities are formed by introducing a mirror in each arm. The electric field at the detector B is

$$\begin{align} E_{\mathrm{out}}\left(0,y,z,\tau\right) = \frac{1}{\sqrt{2}} \sum_{q = 1}^{N_{\mathrm{rt}}\left(\tau\right)} \sqrt{T_{\small\textrm{M}} R_{\small\textrm{M}}^{q-1}} \left[ E^{\left(\small\textrm{C}\right)}_{y}\left(d_{\small\textrm{C}},y,z,\tau\right) - E^{\left(\small\textrm{D}\right)}_{y}\left(d_{\small\textrm{D}},y,z,\tau\right) e^{- 2 i \varphi_{\text{off}}} \right]. \end{align} \tag{ 38 }$$

Here $N_{\mathrm{rt}}(\tau) = \lfloor\frac{c \tau}{2 \mathcal{L}}\rfloor$ and $T_{\small\textrm{M}} = 1-R_{\small\textrm{M}}$ is the transmission coefficient of the mirrors introduced to render arm cavities. The analogous expression without arm cavities corresponds to equation (20).

Zoom In Zoom Out Reset image size

The PSD to first-order in Planck length is

$$\begin{align} S\left(\,f\,\right) \approx \frac{\mathcal{K} \mathcal{L}}{2 \sqrt{2 \pi}} \left(\frac{\lambda T_{\small\textrm{M}}}{2 \pi}\right)^{2} \, \left(\frac{1- \left(\sqrt{R_{\small\textrm{M}}}\right)^{280}}{1- \sqrt{R_{\small\textrm{M}}}}\right)^{4} \Theta\left(\,f\,\right). \end{align} \tag{ 39 }$$

The derivation is in appendix D. We use $N_{\mathrm{rt}}(\tau) = 280$ based on $T_{\small\textrm{M}} = 0.014,$ whereby $R_{\small\textrm{M}}^{280} \approx 0.019 \lt 0.02$, i.e. less than 2$\%$ of the input light remains after 280 round-trips of the light beam in each arm.

Again, we compare $S(f)$ in equation (39) with the experimentally measured PSD at the frequency of maximum sensitivity ($f_\mathrm{sens} = 235$ Hz, attained at the LIGO Hanford site [44]), and we constrain the product $\Gamma_{\small\textrm{S}} N_{\parallel}$ as in figure 3(a).

Zoom In Zoom Out Reset image size

5.2. Constraint from CSD

The cross-correlation between the output power from the two co-located interferometers is captured using the cross-spectral density $CS(f)$ in equation (26).

We first compute the $CS(f)$ in the case of the Holometer, where the beamsplitters of the two co-located interferometers are separated by 0.91 m diagonally in 3 dimensions with $\delta_{s} = 0.91/\sqrt{3}$ m being the separation along each axis. We consider the two interferometers to have equal arm lengths $\mathcal{L}_{\small\textrm{I}} = \mathcal{L}_{\small\textrm{II}} = \mathcal{L}$. Defining

$$\begin{align} &\mathcal{D} = \text{exp} \left(-\frac{4 \delta_{s}^{2}}{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}}\right), \end{align} \tag{ 40 }$$

and using equations (26), (27) and (B17), we obtain the cross-spectral density, as

$$\begin{align} CS\left(\,f\,\right) = \mathcal{A} \left[\mathcal{B}^{1-\mathcal{D}} + \mathcal{B}^{1+\mathcal{D}} \, \text{cos} 4 \varphi_{\mathrm{off}} - 2 \mathcal{B} \text{cos}^{2} 2 \varphi_{\mathrm{off}} \right]. \end{align} \tag{ 41 }$$

As $\mathcal{D} \rightarrow 0, CS(f) \rightarrow 0$, and as $\mathcal{D} \rightarrow 1, CS(f) \rightarrow S(f)$ in equation (28). For intermediate values of $\mathcal{D}$, the CSD thus provides additional information beyond the PSD. It adds to the constraint on $\Gamma_{\small\textrm{S}} N_{\parallel}$ such as in equation (35) or (37) which appears as a straight line of the log-log plot in figure 3(a). This shows the additional strength of our two-parameter model as compared to a single-parameter one capturing only $\Gamma_{\small\textrm{S}} N_{\parallel}$.

As with the PSD, we use the measured cross-spectral density in the frequency range $f \gt f_{\mathrm{sens}}$, denoted by $|CS_\mathrm{meas}(f\gt f_\mathrm{sens})|$, to set the upper bound on $|CS(f)|$ at  $f$= 0. For the Holometer with $|CS_\mathrm{meas}(f\gt f_\mathrm{sens})| = 10^{-40}$ m² Hz⁻¹ [21], the constraint is shown in figure 3(b). Comparing figures 3 (a) and (b), it is evident that considering the cross-correlations has improved the bound on $\Gamma_{\small\textrm{S}}$ and $N_{\parallel}$. We see that the constraint is improved only where $\ell_{\small\textrm{P}} N_{\parallel}$ is comparable to $\delta_{s}.$ It is thus evident that stricter constraints can be obtained by reducing the separation between the co-located interferometers.

In the proposed QUEST experiment, the beamsplitters are to be separated by 0.3 m along the diagonal in the plane of the interferometer [17]. Using equation (41) with $\delta_{s} = 0.3/2$ m and an estimated $|CS_\mathrm{meas}(f\gt f_\mathrm{sens})| = 10^{-41}$ m² Hz⁻¹ for 6 dB of squeezing we obtain the improved bound on $\Gamma_{\small\textrm{S}}$ and $N_{\parallel}$ displayed in figure 3(b). The choice of $\delta_{s} = 0.3/2$ m in place of $\delta^{{\prime}}_{s} = 0.3/\sqrt{2}$ m is explained in appendix B using equation (B14), while $|CS_\mathrm{meas}(f\gt f_\mathrm{sens})| = 10^{-41}$ m² Hz⁻¹ is obtained by combining (A) shot-noise-limited PSD $S_{\small\textrm{QUEST}}(f)\approx 10^{-38}$ m² Hz⁻¹ for 6 dB of squeezing (evaluating equation (36)), (B) an assumed total measurement time of $T_{\mathrm{tot}} = 10^{6}$ s, and (C) $CS_{\mathrm{meas}}(f)$ is of the order-of-magnitude of $S_{\mathrm{meas}}(f) / \sqrt{T_{\mathrm{tot}}}$.

In the case of aLIGO, the two interferometers are located $\delta_{L} \approx 3030.13$ km apart. The corresponding decay factor $\mathcal{D}$ decreases exponentially and vanishes as $\frac{\delta_{L}}{\ell_{\small\textrm{P}} N_{\parallel}}\gg 1$ given that $N_{\parallel}\ll \mathcal{L}/\ell_{\small\textrm{P}}$. Therefore the limits are set only by the PSD measurements.

Our model of spacetime fluctuations and the constraints placed upon it by experimental data can be translated to some prior models of spacetime fluctuations. One such class is captured by [45]

$$\begin{align} \frac{\Delta x}{\mathcal{L}} \unicode{x2A7E} \left(\frac{\ell_{\small\textrm{P}}}{\mathcal{L}}\right)^{\alpha}, \end{align} \tag{ 42 }$$

where $\Delta x$ is the standard deviation of the path difference in an interferometric setup due to spacetime fluctuations. From equations (33) and (39), it is evident that the PSD is directly proportional to $\mathcal{K}$ and hence, to the product $\Gamma_{\small\textrm{S}} N_{\parallel}$ in all three interferometric setups considered. Therefore, we introduce the combined parameter $\Pi = \Gamma_{\small\textrm{S}} N_{\parallel}$ to relate our two-parameter approach to the single-parameter α model. Comparing $(\Delta x)^2$ with the corresponding measured $(\Delta x)_{\mathrm{meas}}^{2}$ due to noise, given by

$$\begin{align} \left(\Delta x\right)_{\mathrm{meas}}^{2} & = 4 \pi \int\limits_{0}^{\infty} \text{d} f \left[\left.S\left(\,f\,\right) \right\vert_{\Pi_{\mathrm{max}} \equiv \text{bound}\left(\Pi\right)} \right], \end{align} \tag{ 43 }$$

provides a lower bound on $\alpha.$ Here, $\text{bound}(\Pi)$ denotes the value of Π such that it gives the strictest bound $\Gamma_{\small\textrm{S}} N_{\parallel}\unicode{x2A7D} \Pi_{\text{max}}$, from Π computed at the boundary of any allowed parameter space. For instance, in the context of PSD measurements, $\Pi_{\mathrm{max}}$ is the value of Π corresponding to the straight-line boundary of the allowed parameter space in figure 3(a). As these allowed parameter spaces are themselves bounded using the measured data, equation (43) is valid and along with equation (42), it provides the link between $\Pi_{\text{max}}$ and α.

For the Holometer, equation (33) gives

$$\begin{align} \left(\Delta x\right)_{\mathrm{meas}}^{2} = \sqrt{\frac{\pi}{2}} \, \Pi_{\mathrm{max}} \ell_{\small\textrm{P}} \mathcal{L}. \end{align} \tag{ 44 }$$

Combining this with equation (42) gives

$$\begin{align} \alpha \unicode{x2A7E} \frac{1}{2} - \frac{\text{log} \left( \sqrt{\pi/2} \, \Pi_{\mathrm{max}} \right)}{2 \text{log} \left( \mathcal{L}/\ell_{\small\textrm{P}} \right)}. \end{align} \tag{ 45 }$$

Similarly, using equation (39) for aLIGO,

$$\begin{align} \alpha \unicode{x2A7E} \frac{1}{2} - \frac{\text{log} \left( \sqrt{\pi/2} \, \Pi_{\mathrm{max}} \, T_{\small\textrm{M}}^{2} \: \left(\frac{1- \left(\sqrt{R_{\small\textrm{M}}}\right)^{280}}{1- \sqrt{R_{\small\textrm{M}}}}\right)^{4}\right)}{2 \text{log} \left( \mathcal{L}/\ell_{\small\textrm{P}} \right)}. \end{align} \tag{ 46 }$$

Lower bounds on the value of α can thus be obtained from our constraints on $\Gamma_{\small\textrm{S}} N_{\parallel} \unicode{x2A7D} \Pi_{\mathrm{max}}$ obtained for different interferometric setups, such as in equation (35) or (37). We present the numerical estimates in table 2.

Table 2. Comparison of $\Pi_{\mathrm{max}}$ and $\alpha_{\mathrm{min}}$. For completeness, we also quote the $\alpha_{\mathrm{min}}$ values obtained from the astronomical observations [11, table 1].

	Setup	Constraint from	$\Pi_{\mathrm{max}}$	$\alpha_{\mathrm{min}}$
Interferometric	Holometer	PSD	$5 \times 10^{-3\phantom{1}}$	0.53
		CSD	$1 \times 10^{-6\phantom{1}}$	0.58
	QUEST	PSD	$3 \times 10^{-3\phantom{1}}$	0.53
		CSD	$1 \times 10^{-6\phantom{1}}$	0.58
	aLIGO	PSD	$4 \times 10^{-13}$	0.60
Astronomical	Optical	[11]	—	0.53
	x-ray	[11]	—	0.58
	GeV γ-ray	[11]	—	0.67
	TeV γ-ray	[11]	—	0.72

Stronger lower bounds can be placed on α from CSD measurements. This is done by taking the maximum value of $N_{\parallel}$ allowed for different allowed values of $\Gamma_{\small\textrm{S}}$ to compute a set of values $\Gamma_{\small\textrm{S}} N_{\parallel}$; the constraint is then set at $\Pi_{\mathrm{max}}$, which is the minimum value in that set. The corresponding $\alpha_{\mathrm{min}}$ is computed using equation (45) and reported in table 2. We have also presented the $\alpha_{\mathrm{min}}$ values obtained from the astronomical observations [10, 11]. For instance, in [10], spacetime fluctuations are shown to inhibit astronomical observations at sufficiently short wavelengths. This allows using x-ray and γ-ray observations to place stricter constraints as evident from the larger $\alpha_{\mathrm{min}}$. It is interesting to note that the $\alpha_{\mathrm{min}}$ obtained using the CSD measurements in our approach place stricter bounds than that obtained from optical astronomical observations.

Our model is not directly comparable to that by Zurek and coworkers [14–16]. This is because of the different assumptions in the two models of fluctuations. For instance, in [16], the metric itself is identified with a scalar field satisfying the massless free scalar wave equation. This leads to correlations, particularly angular ones, in the spectral densities that are absent in our model. Our Attribute (5) removes possible correlations in the spatial and temporal fluctuations, enables the separation of variables in Assumption (i), and thereby leads to spectral densities that are delta functions.

We have presented a theoretical formalism for evaluating electromagnetic signatures of spacetime fluctuations. It is based on measuring field correlation tensors at different spacetime locations, and enables to bound the strength and the scale of spacetime fluctuations by comparison with experimental data. To that end, we compare the predicted signal to the power and the cross-spectral densities of some of the most sensitive optical Michelson interferometers.

Crucially, our predicted spectral densities cannot be compared to the experimental data directly, as the predicted signal manifests at zero frequency and where interferometers are not sensitive. We thus illustrate the principle of constraining models of spacetime fluctuations by assuming that the model spectral density is present at a frequency at which the interferometer is sensitive. We also provide a correspondence between our bounds and prior single-parameter models of quantum spacetime fluctuations.

Our two-parameter model, capturing the strength and the scale of spacetime fluctuations is able to extract stricter constraints from the cross-correlation between the experimental data recorded in the two co-located interferometers, subject to the above caveat. This is evident not only from our two-dimensional exclusion plots in figure 3(b) but also from the stronger bounds on the α of single-parameter models in table 2. In fact, it is important to note that the bound on α obtained using the CSD measurements is more stringent than that obtained from optical astronomical observations. Subject to the aforementioned caveat, our bound matches that obtained from the shorter wavelength x-ray observations.

The two-dimensional exclusion plots that we introduce should facilitate an unambiguous comparison of various theoretical models of quantum spacetime fluctuations as well as the experiments being designed to detect them. These could include different laboratory-based interferometer configurations [13, 17, 23, 39] and detection schemes such as counting single photons [19] or detecting higher-order field correlations involving both electric and magnetic fields [34] as well as astronomical observations [11]. Finally, our approach is naturally suited to incorporate novel configurations of coupled interferometers using quantum light for quantum gravity tests [37–39]. These could provide new constraints on the structure and nature of spacetime.

Our formalism has already provided insights into experimental routes for detecting quantum spacetime fluctuations. These include the role of the separation between two interferometers as well as their arm lengths in constraining the fluctuations based on their CSD. The former determines the area of the two-dimensional parameter space that is additionally ruled out by the CSD. As to the latter, if the arm lengths of the co-located interferometers are not equal to each other, the spacetime fluctuations can manifest as a complex CSD. Complex CSDs at positive frequencies are typical in interferometric experiments which are limited by uncorrelated noise.

On the theoretical side, future work could identify the specific physical phenomenon that is ignored, due to each assumption or attribute considered when building our model of spacetime fluctuations. Further, our approach can be extended to identify the components of metric fluctuations that can accumulate over a distance vis-a-vis those that cannot. We could also examine the role played by the specific autocorrelation function of the spacetime fluctuations [46]. Further investigations with other motivated functional dependences of spacetime fluctuations can potentially yield better insight into the frequency dependence and angular correlations of quantum spacetime fluctuations. For instance, appendix E shows how changing Attribute (1) can alter the wave equations in our formalism, moving beyond the interpretation as a fluctuating refractive index.

We thank H Grote and K Dooley for extensive discussions and suggestions crucial to this work, as well as Y Shikano, O Kwon, R Roemer, K Zurek and D Steeghs for useful discussions. BS and AD acknowledge the UK STFC 'Quantum Technologies for Fundamental Physics' program (Grant No. ST/T006404/1) for support. SMV acknowledges the support of the Leverhulme Trust under research Grant RPG-2019-022.

No new data were created or analysed in this study.

We recall the relativistic wave equation in terms of the electromagnetic field tensor F_αβ [31]

$$\begin{align} \Box F_{\alpha \beta} + 2 R_{\alpha \gamma \beta \delta} F^{\gamma \delta} - R_{\alpha \gamma} F^{\gamma}_{\phantom{\gamma}\beta} + R_{\beta \gamma} F^{\gamma}_{\phantom{\gamma}\alpha} = 0, \end{align} \tag{ A1 }$$

where

$$\begin{align} \Box F_{\alpha \beta} = g^{\gamma \delta} \nabla_{\gamma} \nabla_{\delta} F_{\alpha \beta} = g^{\gamma \delta} \left( \partial_{\delta} X_{\alpha \beta \gamma} - \Gamma^{\eta}_{\alpha \delta} X_{\eta \beta \gamma} - \Gamma^{\eta}_{\beta \delta} X_{\alpha \eta \gamma} - \Gamma^{\eta}_{\gamma\delta} X_{\alpha \beta \eta} \right), \end{align} \tag{ A2 }$$

with

$$\begin{align} X_{\alpha \beta \gamma} & = \partial_{\gamma} F_{\alpha \beta} - \Gamma^{\eta}_{\beta \gamma} F_{\alpha \eta} - \Gamma^{\eta}_{\alpha \gamma} F_{\eta \beta}, \end{align} \tag{ A3 }$$

$$\begin{align} R_{\alpha \beta \gamma \delta} & = g_{\alpha \nu} \left[ \partial_{\gamma} \Gamma^{\nu}_{\delta \beta} - \partial_{\delta} \Gamma^{\nu}_{\gamma \beta} + \Gamma^{\nu}_{\gamma \eta} \Gamma^{\eta}_{\delta \beta} - \Gamma^{\nu}_{\delta \eta} \Gamma^{\eta}_{\gamma \beta}\right], \end{align} \tag{ A4 }$$

and

$$\begin{align} R_{\alpha \beta} = R^{\nu}_{\alpha \nu \beta}, \end{align} \tag{ A5 }$$

with the Christoffel symbol

$$\begin{align} \Gamma^{\alpha}_{\beta \gamma} = g^{\alpha \mu} \left(\partial_{\beta} g_{\mu \gamma}+\partial_{\gamma} g_{\beta \mu}-\partial_{\mu} g_{\beta \gamma}\right)/2. \end{align} \tag{ A6 }$$

Using the Attribute (1) of our model in section 2 in equation (A1) and defining $G_{\alpha} = \left( \frac{1}{2 \, n^{2} } \partial_{\alpha} n^{2} \right)$, we find

$$\begin{align} \begin{array}{*{20}{c}} {}&{{g^{\gamma \delta }}{\partial _\gamma }{\partial _\delta }{F_{0i}} - {n^2}{\partial _0}\left[ {{G_0}{F_{0i}}} \right] - 2{n^2}{G_0}{\partial _0}{F_{0i}} - 2{n^2}G_0^2{F_{0i}}}\\ {}&{ + \mathop \sum \limits_{j = 1}^3 \left[ { - {n^2}{\partial _0}\left[ {\frac{1}{{{n^2}}}{G_j}{F_{ji}}} \right] + {\partial _j}\left[ {{G_j}{F_{0i}}} \right] - {G_j}{\partial _0}{F_{ji}} - \frac{1}{{2{n^2}}}\left( {\partial _j^2{n^2}} \right){F_{0i}}} \right.}\\ {}&{\left. { - \,2{G_j}{G_0}{F_{ji}} + \left[ {{\partial _j}{G_i}} \right]{F_{0j}} + 4G_j^2{F_{0i}}} \right] = 0,\left( {i = 1,2,3} \right)} \end{array} \end{align} \tag{ A7a }$$

$$\begin{align} & g^{\gamma \delta} \partial_{\gamma} \partial_{\delta} F_{ij} - n^{2} \, \partial_{0} \left[ G_{j} F_{i0} \right] - n^{2} \, \partial_{0} \left[ G_{i} F_{0j} \right] - n^{2} G_{i}\, \partial_{0} F_{0j} \nonumber\\ & \quad - n^{2} G_{j} \, \partial_{0} F_{i0} - 2\, n^{2} \, G_{0} G_{i} \, F_{0j} - 2\, n^{2} \, G_{0} G_{j} \, F_{i0} - n^{2} G_{0} \, \partial_{0} F_{ij} \nonumber\\ & \quad - \sum_{q = 1}^{3} \left[ \left\{ \partial_{q} G_{i} \right\} F_{qj} + \left\{ \partial_{q} G_{j} \right\} F_{iq} + G_{q} \, \partial_{q} F_{ij} \right] = 0, \quad \left(i,j = 1,2,3\right) \quad \end{align} \tag{ A7b }$$

with $\alpha,\gamma,\delta = 0,1,2,3$.

Assuming the inequalities in (8a), all terms comprising G_α , any derivative of G_α or $\partial_{j}^{2} n^{2}$ in equation (A7a) are neglected. This reduces the wave equation to

$$\begin{align} g^{\gamma \delta} \partial_{\delta} \partial_{\gamma} F_{\alpha \beta} = 0. \end{align} \tag{ A8 }$$

As an aside, for better illustration of the different terms in equation (A7a), we consider the dimensions to reduce to just $(c t,x)$ in place of $(c t,x,y,z)$, with $\partial_{t}$ and $\partial_{x}$ denoting the partial derivatives with respect to x and t respectively. This simplifies equation (A7a) to

$$\begin{align} & -n^{2} \partial_{t}^{2} E_{x} + c^{2} \partial_{x}^{2} E_{x} - \frac{c^{2}}{2 \, n^{2}} \left(\partial_{x}^{2} n^{2} \right) E_{x} - n^{2} \, \left[ \partial_{t} \left( \frac{1}{2 \, n^{2} } \partial_{t} n^{2} \right)\right] E_{x} \nonumber\\ & \quad + 2 \, c^{2} \, \left[ \partial_{x} \left( \frac{1}{2 \, n^{2} } \partial_{x} n^{2} \right) \right] E_{x} + c^{2} \, \left( \frac{1}{2 \, n^{2} } \partial_{x} n^{2} \right) \, \partial_{x} E_{x} \nonumber\\ & \quad - 3 \, n^{2} \, \left( \frac{1}{2 \, n^{2} } \partial_{t} n^{2} \right) \, \partial_{t} E_{x} - 2 \, n^{2} \, \left( \frac{1}{2 \, n^{2} } \partial_{t} n^{2} \right)^{2} \, E_{x} +4 \, c^{2} \, \left( \frac{1}{2 \, n^{2} } \partial_{x} n^{2} \right)^{2} \, E_{x} = 0. \quad \end{align} \tag{ A9 }$$

Here $E_{x}(x,t)$ is the electric field along the only spatial dimension x. We note that equation (A7b) has no equivalent equation when considering only one spatial dimension as there can be no magnetic field in such a Universe. In order to better illustrate the terms involved, we simplify equation (A9) assuming $n(x,t) = 1+w(x,t)$ with $w(x,t)\ll 1$ and neglecting $w^q$ for $q\unicode{x2A7E} 2$. This reduces equation (A9) to

$$\begin{align} & \partial_{t}^{2} E_{x} - c^{2} \, \partial_{x}^{2} E_{x} + 2 \, c^{2} \, w \, \partial_{x}^{2} E_{x} + 3 \, \left(1-w\right) \, \left( \partial_{t} w \right) \, \left(\partial_{t} E_{x}\right) \nonumber\\ & \quad - c^{2} \left(1-3 w\right) \left(\partial_{x} w\right) \left(\partial_{x} E_{x}\right) - c^{2} \, \left(1- 3 w\right) \left(\partial_{x}^{2} w \right) E_{x}\nonumber\\ & \quad - c^{2} \, \left(1 - 12 w\right) \left(\partial_{x} w \right)^{2} \, E_{x} + \left(1 - w\right) \left( \partial_{t}^{2} w \right) \, E_{x} + \left(1-6 w\right) \, \left( \partial_{t} w \right)^{2} E_{x} = 0. \end{align} \tag{ A10 }$$

A detailed analysis of the effects of each of these terms involved could potentially reveal interesting signatures of spacetime fluctuation.

In this section, we find $\overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)}$ with $\boldsymbol{R} = \{(0,y,z,\tau), (\delta_{x},y+\delta_{y},z+\delta_{z},\tau+\delta_{\tau})\}$, where the electric field corresponding to the output of the p-interferometer ($p = \text{I,II}$) is

$$\begin{align} E^{p}_{\mathrm{out}}\left(\delta_{x},y+\delta_{y},z+\delta_{z},\tau+\delta_{\tau}\right) & = \frac{1}{\sqrt{2}} \left[ E^{\left(\small\textrm{C}_{p}\right)}_{y}\left(\delta_{x},y+\delta_{y},z+\delta_{z},\tau+\delta_{\tau}\right) \right.\nonumber\\ & \quad \left. -\ E^{\left(\small\textrm{D}_{p}\right)}_{y}\left(\delta_{x},y+\delta_{y},z+\delta_{z},\tau+\delta_{\tau}\right) e^{- 2 i \varphi_{\text{off}}} \right]. \quad \end{align} \tag{ B1 }$$

We note that we have invoked Attribute (4) and $\overline{\mathcal{O}}$ denotes the average of any quantity $\mathcal{O}$ over realisations of the underlying random process. The fourth-order correlation $M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)$ is written explicitly in terms of $M^{2,2}_{\{\small\textrm{X,Y}\};\{\small\textrm{X}^{{\prime}},\small\textrm{Y}^{{\prime}}\}} \left( \boldsymbol{R};\boldsymbol{R} \right)$ as

$$\begin{align} M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right) & = \frac{1}{4} \left[ \sum_{\text{X}\in \left\{\text{C}_{p},\text{D}_{p}\right\}}\sum_{\text{Y}\in \left\{\text{C}_{q},\text{D}_{q}\right\}}M^{2,2}_{\left\{\small\textrm{X,Y}\right\};\left\{\small\textrm{X,Y}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) - \sum_{\text{X}\in \left\{\text{C}_{p},\text{D}_{p}\right\}} M^{2,2}_{\left\{\small\textrm{X,C}_{q}\right\};\left\{\small\textrm{X,D}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) e^{2 i \varphi_{\text{off}}} \right.\nonumber\\ & - \sum_{\text{X}\in \left\{\text{C}_{p},\text{D}_{p}\right\}} M^{2,2}_{\left\{\small\textrm{X,D}_{q}\right\};\left\{\small\textrm{X,C}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) e^{- 2 i \varphi_{\text{off}}} - \sum_{\text{X}\in \left\{\text{C}_{q},\text{D}_{q}\right\}} M^{2,2}_{\left\{\small\textrm{C}_{p},\small\textrm{X}\right\};\left\{\small\textrm{D}_{p},\small\textrm{X}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) e^{2 i \varphi_{\text{off}}}\nonumber\\ & - \sum_{\text{X}\in \left\{\text{C}_{q},\text{D}_{q}\right\}} M^{2,2}_{\left\{\small\textrm{D}_{p},\small\textrm{X}\right\};\left\{\small\textrm{C}_{p},\small\textrm{X}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) e^{- 2 i \varphi_{\text{off}}} + M^{2,2}_{\left\{\small\textrm{C}_{p},\small\textrm{C}_{q}\right\};\left\{\small\textrm{D}_{p},\small\textrm{D}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) e^{4 i \varphi_{\text{off}}}\nonumber\\ & \left. +\ M^{2,2}_{\left\{\small\textrm{D}_{p},\small\textrm{D}_{q}\right\};\left\{\small\textrm{C}_{p},\small\textrm{C}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) e^{- 4 i \varphi_{\text{off}}} + M^{2,2}_{\left\{\small\textrm{D}_{p},\small\textrm{C}_{q}\right\};\left\{\small\textrm{C}_{p},\small\textrm{D}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) + M^{2,2}_{\left\{\small\textrm{C}_{p},\small\textrm{D}_{q}\right\};\left\{\small\textrm{D}_{p},\small\textrm{C}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right) {} \right]\!. \end{align} \tag{ B2 }$$

Here $p,q = \text{I, II}$. Further, in the case p = q, $\delta_{x} = \delta_{y} = \delta_{z} = 0$. We discuss the salient steps involved in computing $M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)$ by listing the steps in computing one of the terms in this moment, such as,

$$\begin{align} & \overline{M^{2,2}_{\left\{\small\textrm{D}_{p},\small\textrm{D}_{q}\right\};\left\{\small\textrm{C}_{p},\small\textrm{C}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right)}\nonumber\\ & \quad = \int_{-\infty}^{\infty} \text{d} k_{y_1} \int_{-\infty}^{\infty} \text{d} k_{z_1} \, \int_{-\infty}^{\infty} \text{d} k_{y_2} \int_{-\infty}^{\infty} \text{d} k_{z_2} \, \int_{-\infty}^{\infty} \text{d} k_{y_3} \int_{-\infty}^{\infty} \text{d} k_{z_3} \, \int_{-\infty}^{\infty} \text{d} k_{y_4} \int_{-\infty}^{\infty} \text{d} k_{z_4} \nonumber\\ & \qquad \times \overline{\widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{D}_{p}\right)}\left(k_{z_1},k_{y_1},0\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{C}_{p}\right) \ast}\left(0,k_{y_2},k_{z_2}\right)\widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{D}_{q}\right)}\left(k_{z_3},k_{y_3},\delta_{z}\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{C}_{q}\right) \ast}\left(\delta_{x},k_{y_4},k_{z_4}\right)} \nonumber\\ & \qquad \times \overline{E_{\small\textrm{T}}^{\left(\small\textrm{D}_{p}\right)}\left(\tau\right) E_{\small\textrm{T}}^{\left(\small\textrm{D}_{q}\right)}\left(\tau+\delta_{\tau}\right) E_{\small\textrm{T}}^{\left(\small\textrm{C}_{p}\right) \ast}\left(\tau\right) E_{\small\textrm{T}}^{\left(\small\textrm{C}_{q}\right) \ast}\left(\tau+\delta_{\tau}\right)} \nonumber\\ & \qquad \times e^{-i \left(k_{y_1} y + k_{z_1} z - k_{y_2} y - k_{z_2} z\right)} e^{-i \left(k_{y_3} \left(y+\delta_{y}\right) + k_{z_3} \left(z+\delta_{x}\right) - k_{y_4} \left(y+\delta_{y}\right) - k_{z_4} \left(z+\delta_{z}\right)\right)}. \end{align} \tag{ B3 }$$

We note here that we have used the Attribute (5) in obtaining the above expression. It is evident from equation (17) that $\overline{E_{\small\textrm{T}}^{(\small\textrm{D}_{p})}(\tau) E_{\small\textrm{T}}^{(\small\textrm{D}_{q})}(\tau+\delta_{\tau}) E_{\small\textrm{T}}^{(\small\textrm{C}_{p}) \ast}(\tau) E_{\small\textrm{T}}^{(\small\textrm{C}_{q}) \ast}(\tau+\delta_{\tau})} = |E_{\small\textrm{T}}(0)|^{4}$.

We also find using equations (16) and (18),

$$\begin{align} &\overline{\widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{D}_{p}\right)}\left(k_{z_1},k_{y_1},0\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{C}_{p}\right) \ast}\left(0,k_{y_2},k_{z_2}\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{D}_{q}\right)}\left(k_{z_3},k_{y_3},\delta_{z}\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{C}_{q}\right) \ast}\left(\delta_{x},k_{y_4},k_{z_4}\right)} \nonumber\\ & \quad = \frac{1}{4} \, \overline{\Phi_{z}\left(0,0,0\right) \Phi_{z}\left(\delta_{x},\delta_{y},\delta_{z}\right) \Phi_{x}^{\ast}\left(0,0,0\right) \Phi_{x}^{\ast}\left(\delta_{x},\delta_{y},\delta_{z}\right)} \: \widetilde{E}_{\mathrm{in}}\left(k_{z_1},k_{y_1},0\right) \, \widetilde{E}_{\mathrm{in}}^{\ast}\left(0,k_{y_2},k_{z_2}\right) \, \widetilde{E}_{\mathrm{in}}\left(k_{z_3},k_{y_3},0\right)\nonumber\\ & \qquad \times \widetilde{E}_{\mathrm{in}}^{\ast}\left(0,k_{y_4},k_{z_4}\right) \, \text{exp} \left[ \frac{i \mathcal{L}}{k} \left( k_{y_1}^{2} + k_{z_1}^{2} - k_{y_2}^{2} - k_{z_2}^{2} \right)\right] \text{exp} \left[ \frac{i \mathcal{L}}{k} \left( k_{y_3}^{2} + k_{z_3}^{2} - k_{y_4}^{2} - k_{z_4}^{2} \right)\right], \end{align} \tag{ B4 }$$

where

$$\begin{align} &\Phi_{z}\left(\delta_{x},\delta_{y},\delta_{z}\right) = \text{exp} \left[ -i k \int\limits_{0}^{2 \mathcal{L}} \text{d} z^{{\prime}} \, w_{\small\textrm{S}}\left(z+\delta_{x},y+\delta_{y},z^{{\prime}}+\delta_{z}\right) \right], \end{align} \tag{ B5a }$$

$$\begin{align} &\Phi_{x}\left(\delta_{x},\delta_{y},\delta_{z}\right) = \text{exp} \left[ -i k \int\limits_{0}^{2 \mathcal{L}} \text{d} x^{{\prime}} \, w_{\small\textrm{S}}\left(x^{{\prime}}+\delta_{x},y+\delta_{y},z+\delta_{z}\right) \right]. \end{align} \tag{ B5b }$$

Here we have considered the case $\mathcal{L}_{\small\textrm{I}} = \mathcal{L}_{\small\textrm{II}} = \mathcal{L}$, i.e. the co-located interferometers are, in principle, completely identical. This is the reason why the LHS of equation (B4) depends on $p,q$ while the RHS shows no such dependence.

Further, $\widetilde{E}_{\mathrm{in}}$ is the Fourier transformed field factor corresponding to the input beam at BS_p ($p = \text{I, II}$), while $\widetilde{E}_{\small\textrm{S}}^{(\small\textrm{C}_{p})}$ and $\widetilde{E}_{\small\textrm{S}}^{(\small\textrm{D}_{p})}$ correspond to the output measured at B_p . For instance, in the Holometer, we know that $\delta_{i}\approx 0.52$ m and $\mathcal{L} = 40$ m with $\delta_{i}\lt 0.1 \, \mathcal{L}$ ($i = x,y,z$). We also note that in the above expressions, arguments such as $z+\delta_{x}$ feature in the definition of $\Phi_{z}$, in the context of field contributions due to the light traversing along z-axis in the D_p arm. This unusual argument is because of the following reason. We know that the light in D_p suffers a reflection at BS_p before interference with the light from C_p arm. Specifically, the beam cross-section on the x − y plane in D_p arm on reflection becomes the y − z plane at the detector, especially with $x$  →  $z$ and $y$  →  $y$ in terms arising due to light propagation in the D_p arm.

Denoting the input electric field coming in from port A_p by $E_{\text{in}}(x,y,0)$, the electric field $E_{\text{in}}(x,y,0)$ on reflection at BS_p into arm C_p becomes $E_{\text{in}}(0,y,z)$. We emphasize here that in equation (B4), the factors $\widetilde{E}_{\mathrm{in}}(k_{z_1},k_{y_1},0)$ and $\widetilde{E}_{\mathrm{in}}(k_{z_3},k_{y_3},0)$ are Fourier transforms of $E_{\text{in}}(x,y,0)$ in $x$ and $y$ , while $\widetilde{E}_{\mathrm{in}}(0,k_{y_2},k_{z_2})$ and $\widetilde{E}_{\mathrm{in}}(0,k_{y_4},k_{z_4})$ are Fourier transforms of $E_{\text{in}}(0,y,z)$ in $y$ and $z$.

We then use the Attribute (4) to simplify $\overline{\Phi_{z}(0,0,0) \Phi_{z}(\delta_{x},\delta_{y},\delta_{z}) \Phi_{x}^{\ast}(0,0,0) \Phi_{x}^{\ast}(\delta_{x},\delta_{y},\delta_{z})}$ as follows. We first define

$$\begin{align} \sigma_{z}\left(\delta_{x},\delta_{y}\right) & = \frac{1}{2 \mathcal{L}} \int\limits_{0}^{2 \mathcal{L}} \text{d} \zeta_{+} \int\limits_{-\zeta_{+}}^{\zeta_{+}} \text{d} \zeta_{-} \, \rho\left(\delta_{x},\delta_{y},\zeta_{-}\right), \end{align} \tag{ B6a }$$

$$\begin{align} \sigma_{x}\left(\delta_{y},\delta_{z}\right) & = \frac{1}{2 \mathcal{L}}\int\limits_{0}^{2 \mathcal{L}} \text{d} \chi_{+} \,\int\limits_{-\chi_{+}}^{\chi_{+}} \text{d} \chi_{-} \, \rho\left(\chi_{-},\delta_{y},\delta_{z}\right), \end{align} \tag{ B6b }$$

$$\begin{align} \xi_{1}\left(z\right) & = \int\limits_{0}^{2 \mathcal{L}} \text{d} x^{{\prime}} \, \int\limits_{0}^{2 \mathcal{L}} \text{d} z^{{\prime}} \, \rho\left(x^{{\prime}}-z,0,z-z^{{\prime}}\right), \end{align} \tag{ B6c }$$

$$\begin{align} \xi_{2 x}\left(c_{0},c_{1},c_{2},c_{3}\right) & = \int\limits_{0}^{2 \mathcal{L}} \text{d} x \, \int\limits_{0}^{2 \mathcal{L}} \text{d} z \, \rho\left(x+\delta_{x}-c_{0},c_{2}-c_{1},c_{3}-z\right), \end{align} \tag{ B6d }$$

$$\begin{align} \xi_{2 z}\left(c_{0},c_{1},c_{2},c_{3}\right) & = \int\limits_{0}^{2 \mathcal{L}} \text{d} x \, \int\limits_{0}^{2 \mathcal{L}} \text{d} z \, \rho\left(x-c_{0},c_{2}-c_{1},c_{3}-z-\delta_{z}\right), \end{align} \tag{ B6e }$$

where the σ's arise from spatial correlations in the spacetime metric fluctuations within an arm of the interferometer, and the ξ's correspond to spatial correlations of the metric fluctuations between the two arms.

Using equations (B5a) and (B6a), we find

$$\begin{align} & \overline{\Phi_{z}\left(0,0,0\right) \Phi_{z}\left(\delta_{x},\delta_{y},\delta_{z}\right) \Phi_{x}^{\ast}\left(0,0,0\right) \Phi_{x}^{\ast}\left(\delta_{x},\delta_{y},\delta_{z}\right)}\nonumber\\ & \quad = \text{exp} \Bigg[ - k^{2} \Gamma_{\small\textrm{S}} \bigg\{\mathcal{L} \left[2 \sigma_{z}\left(0,0\right) + 2 \sigma_{z}\left(\delta_{x},\delta_{y}\right) \right] - 2 \xi_{1}\left(z\right) \nonumber\\ & \qquad + \mathcal{L} \left[ 2 \sigma_{x}\left(0,0\right) + 2 \sigma_{x}\left(\delta_{y},\delta_{z}\right) \right] - \xi_{2 z}\left(z+\delta_{x},y+\delta_{y},y,z\right) - \xi_{2 x}\left(z,y,y+\delta_{y},z+\delta_{z}\right) \bigg\} \Bigg]. \end{align} \tag{ B7 }$$

We recall here from Attribute (6) that $\rho(\delta_{x},\delta_{y},\delta_{z}) = \text{exp} ( -\sum\nolimits_{i = x,y,z} \frac{\delta_{i}^{2}}{2 \ell_{i}^{2}} )$. For instance, $\sigma_z$ can be simplified using Attribute (6), to

$$\begin{align} \sigma_{z}\left(0,0\right) & = \left[\sqrt{2 \pi} \,\ell_{z} \, \mathrm{erf} \left( \frac{\sqrt{2} \mathcal{L}}{\ell_{z}} \right) + \frac{\ell_{z}^{2}}{\mathcal{L}} \, \left\{\text{exp} \left( - \frac{2 \mathcal{L}^{2}}{\ell_{z}^{2}} \right)-1 \right\} \right], \end{align} \tag{ B8 }$$

where $\mathrm{erf}(z)$ is the error function. We recall that the fluctuation scale $N_{i} = 2\ell_i/\ell_{\small\textrm{P}}$ ($i = x,y,z$). We consider an isotropic fluctuation scale $N_{\parallel} \equiv N_{x} = N_{y} = N_{z}$. Using $N_{\parallel} \ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$,

$$\begin{align} \sigma_{z}\left(0,0\right)\approx \sqrt{\frac{\pi}{2}} \ell_{\small\textrm{P}} N_{\parallel}. \end{align} \tag{ B9 }$$

Similar simplifications can be carried out to find

$$\begin{align} \sigma_{x}\left(0,0\right)\approx \sqrt{\frac{\pi}{2}} \ell_{\small\textrm{P}} N_{\parallel}. \end{align} \tag{ B10 }$$

Further, all ξ-type correlations are negligible in comparison to the σ-type correlations in the limit $N_{\parallel} \ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$. This can be seen even from the qualitative argument that the correlation across two different arms should be significantly smaller than the correlations within an arm when the correlation scale $N_{\parallel} \ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$.

We also consider the input beam to be a perfect Gaussian beam at BS with the input electric field $E_{\text{in}}(x,y,0)$ coming in from port A_p ($p = \text{I, II}$), given by,

$$\begin{align} E_{\text{in}}\left(x,y,0\right) = \sqrt{\frac{2}{\pi} \,} \, \frac{z_{R}}{W_{0} \sqrt{z_{R}^{2} + z_{0}^{2}}} \text{exp} \left[ \left( -\frac{i k z_{0} W_{0}^{2} + 2 z_{R}^{2} }{2 W_{0}^{2} \left(z_{0}^{2}+z_{R}^{2}\right)} \right) \left(x^{2} + y^{2}\right) \right]. \end{align} \tag{ B11 }$$

Here W₀ is the beam waist, z₀ is the position of the beam waist, and $z_{R} = \pi W_{0}^{2}/\lambda$.

Using equations (B4), (B7), (B9), (B10) and (B11) in equation (B3), we find

$$\begin{align} \overline{M^{2,2}_{\left\{\small\textrm{D}_{p},\small\textrm{D}_{q}\right\};\left\{\small\textrm{C}_{p},\small\textrm{C}_{q}\right\}} \left( \boldsymbol{R};\boldsymbol{R} \right)} & = \, |\mu_{1}\left(\mathcal{L},y,z\right) \mu_{1}\left(\mathcal{L},y+\delta_{s},z+\delta_{s}\right)|^{2} \, \text{exp} \left[ - 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \left\{1 + \mathcal{D} \right\} \right]. \end{align} \tag{ B12 }$$

Here,

$$\begin{align} \mu_{1}\left(\mathcal{L},y,z\right) & = \sqrt{\frac{2}{\pi}} \frac{k z_{R} W_{0} \sqrt{z_{R}^{2} + z_{0}^{2}}}{\left(W_{0}^{2} k \left(z_{R}^{2} + z_{0}^{2}\right) - 2 i \mathcal{L} \left(2 z_{R}^{2} + i k W_{0}^{2} z_{0}\right) \right)} \nonumber\\ & \quad \times \text{exp} \left[ - \frac{k \left( 2 z_{R}^{2} + i k W_{0}^{2} z_{0} \right) \left( y^{2} + z^{2} \right) }{2 \left[ W_{0}^{2} k \left(z_{R}^{2} + z_{0}^{2}\right) - 2 i \mathcal{L} \left(2 z_{R}^{2} + i k W_{0}^{2} z_{0}\right) \right]} \right], \quad \end{align} \tag{ B13 }$$

$$\begin{align} \mathcal{D} & = \frac{1}{2} \left[ \text{exp} \left(-\frac{2 \left(\delta_{x}^{2}+\delta_{y}^{2}\right)}{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}}\right) + \text{exp} \left(-\frac{2 \left(\delta_{y}^{2}+\delta_{z}^{2}\right)}{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}}\right) \right]. \end{align} \tag{ B14 }$$

We recall that $\mathcal{K} = k^{2} \Gamma_{\small\textrm{S}} \ell_{\small\textrm{P}} N_{\parallel}$. When $\delta_{x} = \delta_{y} = \delta_{z} = \delta_{s}$, we find $\mathcal{D} = \text{exp} (-\frac{4 \delta_{s}^{2}}{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}})$. In the case $\delta_{x} = \delta_{z} = \delta^{{\prime}}_{s}$ and $\delta_{y} = 0$, we find $\mathcal{D} = \text{exp} (-\frac{2 \delta_{s}^{{\prime} 2}}{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}})$. However, we use the definition $\mathcal{D} = \text{exp} (-\frac{4 \delta_{s}^{2}}{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}})$ (recall equation (40)), and in the case when $\delta_{x} = \delta_{z} = \delta^{{\prime}}_{s}$ and $\delta_{y} = 0$, we consider $\delta_{s} = \delta^{{\prime}}_{s}/\sqrt{2}$.

Using arguments similar to those used in obtaining $\overline{M^{2,2}_{\{\small\textrm{D}_{p},\small\textrm{D}_{q}\};\{\small\textrm{C}_{p},\small\textrm{C}_{q}\}} ( \boldsymbol{R};\boldsymbol{R} )}$, we can compute the other terms in equation (B2) to obtain $\overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)}$. Therefore,

$$\begin{align} \overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)} & = \, |\mu_{1}(\mathcal{L}, y,z)|^{2} \, \left|\mu_{1}\left(\mathcal{L},y+\delta_{s},z+\delta_{s}\right)\right|^{2} \Bigg( 1 - 2 \text{exp} \left[ - \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right] \text{cos} 2 \varphi_{\mathrm{off}} \nonumber\\ & \quad + \frac{1}{2} \text{exp} \left[ - 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right] \left\{\text{exp} \left[ 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \mathcal{D} \right] + \text{exp} \left[ - 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \mathcal{D} \right] \text{cos} 4 \varphi_{\mathrm{off}} \right\} \Bigg). \end{align} \tag{ B15 }$$

The corresponding two-point correlation of power is obtained using

$$\begin{align} \overline{P_{\mathrm{out}}^{i}\left(\tau\right) P_{\mathrm{out}}^{j}\left(\tau+\delta_{\tau}\right)} = \left(\epsilon_{0} c\right)^{2} \, \int_{A} \mathrm{d}y_{1} \, \mathrm{d}z_{1} \, \int_{A} \mathrm{d}y_{2} \, \mathrm{d}z_{2} \: \overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)},\quad \end{align} \tag{ B16 }$$

where $\int_{A}$ denotes a surface integral over the beam cross-section. We also note that the co-located interferometers in the Holometer or the QUEST experiment are built such that $\delta_{s}\,\gg \sqrt{A}$. For instance, $\delta_{s} = 0.91/\sqrt{3}$ m, $\sqrt{A} = 5 \sqrt{\pi} \times 10^{-3}$ m in the case of the Holometer. Therefore, we see that $y_{2}-y_{1}\approx \delta_{s}$ and $z_{2}-z_{1}\approx\delta_{s}$ in the integrals considered above. Using this in addition to equations (B15), (B16) and (24), we find

$$\begin{align} \mathrm{Cov}_{i,j}\left(P_{\mathrm{out}}\right)& = \frac{P_{0}^{2}}{4} \, \left( \frac{e^{- \sqrt{2 \pi} \mathcal{K} \mathcal{L}}}{2} \left\{e^{- \sqrt{2 \pi} \mathcal{K} \mathcal{L} \left(1 - 2 \mathcal{D} \right) } + e^{- \sqrt{2 \pi} \mathcal{K} \mathcal{L} \left(1 + 2 \mathcal{D} \right) } \text{cos} 4 \varphi_{\mathrm{off}} \right\}\right.\nonumber\\ & \quad \left. -\ e^{- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} } \text{cos}^{2} 2 \varphi_{\mathrm{off}} {}\right). \end{align} \tag{ B17 }$$

Using equation (27) and carrying out the cosine transform in equation (25), we obtain the PSD taking $\mathcal{D} = 1$ (as $\delta_{s} = 0$) in the above expression. The CS (with $\delta_{s}\neq 0$) corresponding to equation (B17) is obtained easily using equations (27) and (26) when $\mathcal{L}_{\small\textrm{I}} = \mathcal{L}_{\small\textrm{II}} = \mathcal{L}$. We note here that when $\mathcal{L}_{\small\textrm{I}} = \mathcal{L}_{\small\textrm{II}} = \mathcal{L}$, $\mathrm{Cov}_{\small\textrm{I,II}}(\Delta x) = \mathrm{Cov}_{\small\textrm{II,I}}(\Delta x)$ in equation (26).

Considering $N_{\parallel} = 2 \mathcal{L}/\ell_{\small\textrm{P}}$ in the case of the Holometer and the QUEST experiment, we find that equations (B6a)–(B6e) give

$$\begin{align} &\sigma_{z}\left(\delta_{s},\delta_{s}\right) = \sigma_{x}\left(\delta_{s},\delta_{s}\right) = \mathcal{L} \mathcal{D} \mathcal{F}_{1}, \end{align} \tag{ C1 }$$

$$\begin{align} &\xi_{1}\left(z\right) = \mathcal{L}^{2} \mathcal{F}_{2}, \end{align} \tag{ C2 }$$

$$\begin{align} &\nonumber \xi_{2 x}\left(z,y,y+\delta_{s},z+\delta_{s}\right) = \nonumber\\ &\quad\xi_{2 z}\left(z+\delta_{s},y+\delta_{s},y,z\right) = \mathcal{L}^{2} \mathcal{F}_{2} \sqrt{\mathcal{D}}, \end{align} \tag{ C3 }$$

$$\begin{align} &\mathrm{with} \: \mathcal{F}_{1} = \sqrt{2 \pi} \, \mathrm{erf}\left(\sqrt{2}\right) + e^{-2} -1, \end{align} \tag{ C4 }$$

$$\begin{align} &\mathrm{and} \: \mathcal{F}_{2} = \frac{\pi}{2} \left( \mathrm{erf}\left(\sqrt{2}\right) \right)^{2}. \end{align} \tag{ C5 }$$

Using the steps illustrated in appendix B, the PSD in this limit $N_{\parallel} = 2 \mathcal{L}/\ell_{\small\textrm{P}}$ is

$$\begin{align} S\left(\,f\,\right) = \mathcal{A} \, &\left[\right. 1 + \, e^{- 2 \left(2 \mathcal{F}_{1} -\mathcal{F}_{2}\right) \mathcal{K} \mathcal{L}} \, \text{cos} 4 \varphi_{\mathrm{off}} - 2 e^{- \left(2 \mathcal{F}_{1} -\mathcal{F}_{2}\right) \mathcal{K} \mathcal{L}} \, \text{cos}^{2} 2 \varphi_{\mathrm{off}} \left.\right], \end{align} \tag{ C6 }$$

and the cross-spectral density is

$$\begin{align} CS\left(\,f\,\right) & = \mathcal{A} \, \left[ e^{- \left(2 \mathcal{F}_{1} -\mathcal{F}_{2}\right) \mathcal{K} \mathcal{L}} \, e^{\mathcal{K} \mathcal{L} \sqrt{\mathcal{D}} \left(2 \mathcal{F}_{1} \sqrt{\mathcal{D}} - \mathcal{F}_{2} \right)} + \, e^{- \left(2 \mathcal{F}_{1} -\mathcal{F}_{2}\right) \mathcal{K} \mathcal{L}} \, e^{-\mathcal{K} \mathcal{L} \sqrt{\mathcal{D}} \left(2 \mathcal{F}_{1} \sqrt{\mathcal{D}} - \mathcal{F}_{2} \right)} \text{cos} 4 \varphi_{\mathrm{off}} \right.\nonumber\\ & \left. \qquad -\ 2 e^{- \left(2 \mathcal{F}_{1} -\mathcal{F}_{2}\right) \mathcal{K} \mathcal{L}} \, \text{cos}^{2} 2 \varphi_{\mathrm{off}}\right]. \end{align} \tag{ C7 }$$

To first-order in Planck length,

$$\begin{align} S\left(\,f\,\right) & \approx \left(\frac{\lambda}{2 \pi} \right)^{2} \frac{\left(2 \mathcal{F}_{1}-\mathcal{F}_{2}\right) \mathcal{K} \mathcal{L}}{8 \pi} \delta\left(\,f\,\right), \end{align} \tag{ C8 }$$

$$\begin{align} CS\left(\,f\,\right) & \approx \frac{\sqrt{\mathcal{D}}}{2} \left(2 \sqrt{D} \mathcal{F}_{1} -\mathcal{F}_{2}\right) \Gamma_{\small\textrm{S}} \mathcal{L}^{2} \delta\left(\,f\,\right). \end{align} \tag{ C9 }$$

Comparing $S(f)$ with the measured PSD $S_{\mathrm{meas}}(f)$, we get

$$\begin{align} \Gamma_{\small\textrm{S}} \unicode{x2A7D} \frac{\hbar c \lambda}{2 \left(2 \mathcal{F}_{1}-\mathcal{F}_{2}\right) \mathcal{L}^{2} P_{0}}. \end{align} \tag{ C10 }$$

Therefore, $\Gamma_{\small\textrm{S}}\unicode{x2A7D} 2 \times 10^{-38}$ in the Holometer, while $\Gamma_{\small\textrm{S}}\unicode{x2A7D} 2 \times 10^{-37}$ in the QUEST experiment with 6 dB of squeezing. Comparing $|CS(f)|$ with $C_{\mathrm{max}}$, we find that $\Gamma_{\small\textrm{S}}\unicode{x2A7D} 10^{-43}$ in the Holometer, while $\Gamma_{\small\textrm{S}}\unicode{x2A7D} 10^{-42}$ in the QUEST experiment. It is clear that the cross-spectral density improves the bound on $\Gamma_{\small\textrm{S}}$ for the given $N_{\parallel}$.

In this section, we find $\overline{M^{2,2}_{\text{out}} \left(\boldsymbol{R}_{0};\boldsymbol{R}_{0}\right)}$ with $\boldsymbol{R}_{0} = \{(0,y,z,\tau), (0,y,z,\tau)\}$, where

$$\begin{align} &E_{\mathrm{out}}\left(0,y,z,\tau\right) = \frac{1}{\sqrt{2}} \sum\limits_{q = 1}^{N_{\mathrm{rt}}\left(\tau\right)} \sqrt{T_{\small\textrm{M}} R_{\small\textrm{M}}^{q-1}} \, \left[ E^{\left(\small\textrm{C}\right)}_{y}\left(d_{\small\textrm{C}},y,z,\tau\right) - E^{\left(\small\textrm{D}\right)}_{y}\left(d_{\small\textrm{D}},y,z,\tau\right) e^{- 2 i \varphi_{\text{off}}} \right]. \end{align} \tag{ D1 }$$

We note that we have invoked Attribute (4) and $\overline{\mathcal{O}}$ denotes the average of any quantity $\mathcal{O}$ over realisations of the underlying random process. As mentioned in section 5.1, we take $N_{\mathrm{rt}}(\tau) = 280$ using the fact that less than 2$\%$ of the input light remains after 280 round-trips of the light beam in each arm. The fourth-order correlation $M^{2,2}_{\text{out}} \left(\boldsymbol{R}_{0};\boldsymbol{R}_{0}\right)$ is written explicitly in terms of $M^{2,2}_{\{\small\textrm{X,Y}\};\{\small\textrm{X}^{{\prime}},\small\textrm{Y}^{{\prime}}\}} \left( \boldsymbol{R}_{\small\textrm{X,Y}};\boldsymbol{R}_{\small\textrm{X}^{{\prime}},\small\textrm{Y}^{{\prime}}} \right)$ as

$$\begin{align} \begin{array}{l} \begin{array}{*{20}{c}} {M_{{\textrm{out}}}^{2,2}\left( {{R_0};{R_0}} \right)}&{ = \frac{1}{4}\mathop \sum \limits_{\begin{array}{*{20}{c}} {{q_1},{q_2}}\\ {{q_3},{q_4}} \end{array} = 1}^{280} T_{\textrm{M}}^2{{\left( {\sqrt {{R_{\textrm{M}}}} } \right)}^{{q_1} + {q_2} + {q_3} + {q_4} - 4}}\left[ {\mathop \sum \limits_{X \in \left\{ {C,D} \right\}} \mathop \sum \limits_{Y \in \left\{ {C,D} \right\}} M_{\left\{ {{\textrm{X,Y}}} \right\};\left\{ {{\textrm{X,Y}}} \right\}}^{2,2}\left( {{R_{{\textrm{X,Y}}}};{R_{{\textrm{X,Y}}}}} \right)} \right.} \end{array}\\ \begin{array}{*{20}{c}} {}&{ - \mathop \sum \limits_{X \in \left\{ {C,D} \right\}} M_{\left\{ {{\textrm{X,C}}} \right\};\left\{ {{\textrm{X,D}}} \right\}}^{2,2}\left( {{R_{{\textrm{X,C}}}};{R_{{\textrm{X,D}}}}} \right){e^{2i{\varphi _{{\textrm{off}}}}}} - \mathop \sum \limits_{X \in \left\{ {C,D} \right\}} M_{\left\{ {{\textrm{X,D}}} \right\};\left\{ {{\textrm{X,C}}} \right\}}^{2,2}\left( {{R_{{\textrm{X,D}}}};{R_{{\textrm{X,C}}}}} \right){e^{ - 2i{\varphi _{{\textrm{off}}}}}}} \end{array}\\ \begin{array}{*{20}{c}} {}&{ - \mathop \sum \limits_{X \in \left\{ {C,D} \right\}} M_{\left\{ {{\textrm{C,X}}} \right\};\left\{ {{\textrm{D,X}}} \right\}}^{2,2}\left( {{R_{{\textrm{C,X}}}};{R_{{\textrm{D,X}}}}} \right){e^{2i{\varphi _{{\textrm{off}}}}}} - \mathop \sum \limits_{X \in \left\{ {C,D} \right\}} M_{\left\{ {{\textrm{D,X}}} \right\};\left\{ {{\textrm{C,X}}} \right\}}^{2,2}\left( {{R_{{\textrm{D,X}}}};{R_{{\textrm{C,X}}}}} \right){e^{ - 2i{\varphi _{{\textrm{off}}}}}}} \end{array}\\ \begin{array}{*{20}{c}} {}&{ + M_{\left\{ {{\textrm{C,C}}} \right\};\left\{ {{\textrm{D,D}}} \right\}}^{2,2}\left( {{R_{{\textrm{C,C}}}};{R_{{\textrm{D,D}}}}} \right){e^{4i{\varphi _{{\textrm{off}}}}}} + M_{\left\{ {{\textrm{D,D}}} \right\};\left\{ {{\textrm{C,C}}} \right\}}^{2,2}\left( {{R_{{\textrm{D,D}}}};{R_{{\textrm{C,C}}}}} \right){e^{ - 4i{\varphi _{{\textrm{off}}}}}}} \end{array}\\ \begin{array}{*{20}{c}} {}&{\left. { + \,M_{\left\{ {{\textrm{D,C}}} \right\};\left\{ {{\textrm{C,D}}} \right\}}^{2,2}\left( {{R_{{\textrm{D,C}}}};{R_{{\textrm{C,D}}}}} \right) + M_{\left\{ {{\textrm{C,D}}} \right\};\left\{ {{\textrm{D,C}}} \right\}}^{2,2}\left( {{R_{{\textrm{C,D}}}};{R_{{\textrm{D,C}}}}} \right)} \right].} \end{array} \end{array} \end{align} \tag{ D2 }$$

Here $\boldsymbol{R}_{\small\textrm{X,Y}} = \{(d_{\small\textrm{X}},y,z,\tau), (d_{\small\textrm{Y}},y,z,\tau)\}$ ($\small\textrm{X,Y}\in \{\small\textrm{C,D}\}$). We discuss the salient steps involved in computing $M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R};\boldsymbol{R}\right)$ by listing the steps in computing one of the terms in this moment, such as,

$$\begin{align} & \sum_{\substack{q_{1},q_{2}, q_{3},q_{4}} = 1}^{280} \left(\sqrt{R_{\small\textrm{M}}}\right)^{q_{1}+q_{2}+q_{3}+q_{4}-4} \: \overline{M^{2,2}_{\left\{\small\textrm{D,D}\right\};\left\{\small\textrm{C,C}\right\}} \left( \boldsymbol{R}_{\small\textrm{D,D}};\boldsymbol{R}_{\small\textrm{C,C}} \right)}\nonumber\\ & \quad = \int_{-\infty}^{\infty} \text{d} k_{y_1} \int_{-\infty}^{\infty} \text{d} k_{z_1} \, \int_{-\infty}^{\infty} \text{d} k_{y_2} \int_{-\infty}^{\infty} \text{d} k_{z_2} \, \int_{-\infty}^{\infty} \text{d} k_{y_3} \int_{-\infty}^{\infty} \text{d} k_{z_3} \, \int_{-\infty}^{\infty} \text{d} k_{y_4} \int_{-\infty}^{\infty} \text{d} k_{z_4} \nonumber\\ & \qquad \times \left(\sum_{\substack{q_{1},q_{2}, q_{3},q_{4}} = 1}^{280} \left(\sqrt{R_{\small\textrm{M}}}\right)^{q_{1}+q_{2}+q_{3}+q_{4}-4} \: \overline{\widetilde{E}_{\small\textrm{S}}^{(\small\textrm{D})}(k_{z_1},k_{y_1},0) \widetilde{E}_{\small\textrm{S}}^{(\small\textrm{C}) \ast}(0,k_{y_2},k_{z_2})\widetilde{E}_{\small\textrm{S}}^{(\small\textrm{D})}(k_{z_3},k_{y_3},0) \widetilde{E}_{\small\textrm{S}}^{(\small\textrm{C}) \ast}(0,k_{y_4},k_{z_4})} \right)\nonumber\\ & \qquad \times \overline{\left[E_{\small\textrm{T}}^{\left(\small\textrm{D}\right)}\left(\tau\right)\right]^{2} \, \left[E_{\small\textrm{T}}^{\left(\small\textrm{C}\right) \ast}\left(\tau\right)\right]^{2}} e^{-i \left(k_{y_1} y + k_{z_1} z - k_{y_2} y - k_{z_2} z\right)} e^{-i \left(k_{y_3} y + k_{z_3} z - k_{y_4} y - k_{z_4} z\right)}. \end{align} \tag{ D3 }$$

We note here that we have used the Attribute (5) in obtaining the above expression. It is evident from equation (17) that $\overline{[E_{\small\textrm{T}}^{(\small\textrm{D})}(\tau)]^{2} \, [E_{\small\textrm{T}}^{(\small\textrm{C}) \ast}(\tau)]^{2}} = |E_{\small\textrm{T}}(0)|^{4}$.

We also find using equations (16) and (18),

$$\begin{align} \nonumber &\sum_{\substack{q_{1},q_{2}, q_{3},q_{4}} = 1}^{280} \left(\sqrt{R_{\small\textrm{M}}}\right)^{q_{1}+q_{2}+q_{3}+q_{4}-4} \: \overline{\widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{D}\right)}\left(k_{z_1},k_{y_1},0\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{C}\right) \ast}\left(0,k_{y_2},k_{z_2}\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{D}\right)}\left(k_{z_3},k_{y_3},0\right) \widetilde{E}_{\small\textrm{S}}^{\left(\small\textrm{C}\right) \ast}\left(0,k_{y_4},k_{z_4}\right)} \\ \nonumber & \quad = \frac{1}{4} \, \sum_{\substack{q_{1},q_{2}, q_{3},q_{4}} = 1}^{280} \left(\sqrt{R_{\small\textrm{M}}}\right)^{q_{1}+q_{2}+q_{3}+q_{4}-4} \: \overline{\Phi_{z}^{\left(q_{1}\right)}\left(0,0,0\right) \Phi_{z}^{\left(q_{2}\right)}\left(0,0,0\right) \left[\Phi_{x}^{\left(q_{3}\right)}\left(0,0,0\right) \Phi_{x}^{\left(q_{4}\right)}\left(0,0,0\right)\right]^{\ast}} \\ \nonumber & \qquad \times \widetilde{E}_{\mathrm{in}}\left(k_{z_1},k_{y_1},0\right) \, \widetilde{E}_{\mathrm{in}}^{\ast}\left(0,k_{y_2},k_{z_2}\right) \, \widetilde{E}_{\mathrm{in}}\left(k_{z_3},k_{y_3},0\right) \widetilde{E}_{\mathrm{in}}^{\ast}\left(0,k_{y_4},k_{z_4}\right) \, \\ & \qquad \times \text{exp} \left[ \frac{i \mathcal{L}}{k} \left( k_{y_1}^{2} + k_{z_1}^{2} - k_{y_2}^{2} - k_{z_2}^{2} \right)\right] \, \text{exp} \left[ \frac{i \mathcal{L}}{k} \left( k_{y_3}^{2} + k_{z_3}^{2} - k_{y_4}^{2} - k_{z_4}^{2} \right)\right], \end{align} \tag{ D4 }$$

where

$$\begin{align} &\Phi_{z}^{\left(q\right)}\left(0,0,0\right) = \text{exp} \left[ -i k \int\limits_{0}^{2 q \mathcal{L}} \text{d} z^{{\prime}} \, w_{\small\textrm{S}}\left(z,y,z^{{\prime}}\right) \right], \end{align} \tag{ D5a }$$

$$\begin{align} &\Phi_{x}^{\left(q\right)}\left(0,0,0\right) = \text{exp} \left[ -i k \int\limits_{0}^{2 q \mathcal{L}} \text{d} x^{{\prime}} \, w_{\small\textrm{S}}\left(x^{{\prime}},y,z\right) \right]. \end{align} \tag{ D5b }$$

As pointed out in appendix B, $\widetilde{E}_{\mathrm{in}}$ is the Fourier transformed field factor corresponding to the input beam at BS, while $\widetilde{E}_{\small\textrm{S}}^{(\small\textrm{C})}$ and $\widetilde{E}_{\small\textrm{S}}^{(\small\textrm{D})}$ correspond to the output measured at B. As before, denoting the input electric field coming in from port A by $E_{\text{in}}(x,y,0)$, the electric field $E_{\text{in}}(x,y,0)$ on reflection at BS into arm C becomes $E_{\text{in}}(0,y,z)$. We emphasize here that in equation (D4), the factors $\widetilde{E}_{\mathrm{in}}(k_{z_1},k_{y_1},0)$ and $\widetilde{E}_{\mathrm{in}}(k_{z_3},k_{y_3},0)$ are Fourier transforms of $E_{\text{in}}(x,y,0)$ in $x$ and $y$ , while $\widetilde{E}_{\mathrm{in}}(0,k_{y_2},k_{z_2})$ and $\widetilde{E}_{\mathrm{in}}(0,k_{y_4},k_{z_4})$ are Fourier transforms of $E_{\text{in}}(0,y,z)$ in $y$ and $z$ .

We then use the Attribute (4) to simplify $\overline{\Phi_{z}^{(q_{1})}(0,0,0) \Phi_{z}^{(q_{2})}(0,0,0)[\Phi_{x}^{(q_{3})}(0,0,0)}$ $\overline{\Phi_{x}^{(q_{4})}(0,0,0) ]^{\ast}}$ as follows. We first define

$$\begin{align} &\sigma_{z}^{\left(p,q\right)} = \frac{1}{2 \mathcal{L}} \, \int\limits_{0}^{2 p \mathcal{L}} \text{d} z_{1} \int\limits_{0}^{2 q \mathcal{L}} \text{d} z_{2} \, \rho\left(0,0,z_{2}-z_{1}\right), \end{align} \tag{ D6a }$$

$$\begin{align} &\sigma_{x}^{\left(p,q\right)} = \frac{1}{2 \mathcal{L}} \, \int\limits_{0}^{2 p \mathcal{L}} \text{d} x_{1} \int\limits_{0}^{2 q \mathcal{L}} \text{d} x_{2} \, \rho\left(x_{2}-x_{1},0,0\right), \end{align} \tag{ D6b }$$

$$\begin{align} &\xi_{1}^{\left(p,q\right)}\left(z\right) = \int\limits_{0}^{2 p \mathcal{L}} \text{d} x^{{\prime}} \, \int\limits_{0}^{2 q \mathcal{L}} \text{d} z^{{\prime}} \, \rho\left(x^{{\prime}}-z,0,z-z^{{\prime}}\right), \end{align} \tag{ D6c }$$

where $p,q \in \mathbb{Z}_{+}$. We also note that the σ's arise from spatial correlations in the spacetime metric fluctuations within an arm of the interferometer, and the ξ's correspond to spatial correlations of the metric fluctuations between the two arms.

Using equations (D5a) and (D6a), we find

$$\begin{align} & \nonumber \overline{\Phi_{z}^{\left(q_{1}\right)}\left(0,0,0\right) \Phi_{z}^{\left(q_{2}\right)}\left(0,0,0\right) \left[\Phi_{x}^{\left(q_{3}\right)}\left(0,0,0\right) \Phi_{x}^{\left(q_{4}\right)}\left(0,0,0\right)\right]^{\ast}}\\ & \nonumber \quad = \text{exp} \Bigg[ - k^{2} \Gamma_{\small\textrm{S}} \bigg\{\mathcal{L} \left[ \sigma_{z}^{\left(q_{1},q_{1}\right)} + \sigma_{z}^{\left(q_{2},q_{2}\right)} + 2 \sigma_{z}^{\left(q_{1},q_{2}\right)} \right] \\ & \qquad -\ \left( \xi_{1}^{\left(q_{1},q_{3}\right)}\left(z\right) + \xi_{1}^{\left(q_{1},q_{4}\right)}\left(z\right) + \xi_{1}^{\left(q_{2},q_{3}\right)}\left(z\right) + \xi_{1}^{\left(q_{2},q_{4}\right)}\left(z\right)\right)\nonumber\\ & \qquad + \mathcal{L} \left[ \sigma_{x}^{\left(q_{3},q_{3}\right)} + \sigma_{x}^{\left(q_{4},q_{4}\right)}+ 2 \sigma_{x}^{\left(q_{3},q_{4}\right)} \right] \bigg\} \Bigg]. \end{align} \tag{ D7 }$$

We recall here from Attribute (6) that $\rho(\delta_{x},\delta_{y},\delta_{z}) = \text{exp} ( -\sum\nolimits_{i = x,y,z} \frac{\delta_{i}^{2}}{2 \ell_{i}^{2}} )$. We also use the fluctuation scale $N_{i} = 2\ell_i/\ell_{\small\textrm{P}}$ ($i = x,y,z$) with an isotropic fluctuation scale $N_{\parallel} \equiv N_{x} = N_{y} = N_{z}$. This helps in simplifying equation (D6a) as follows. For instance, $\sigma_{z}^{(p,q)}$ can be simplified using Attribute (6), to

$$\begin{align} \nonumber \sigma_{z}^{\left(p,q\right)} & = \sqrt{\frac{\pi}{2}} \,\frac{\ell_{\small\textrm{P}} N_{\parallel}}{2} \, \left[ p \, \mathrm{erf} \left( \frac{2\sqrt{2} p \mathcal{L}}{\ell_{\small\textrm{P}} N_{\parallel}} \right) + q \, \mathrm{erf} \left( \frac{2\sqrt{2} q \mathcal{L}}{\ell_{\small\textrm{P}} N_{\parallel}} \right) - \left(p-q\right) \, \mathrm{erf} \left( \frac{2\sqrt{2} \left(p-q\right) \mathcal{L}}{\ell_{\small\textrm{P}} N_{\parallel}} \right) \right]\\ & \quad - \frac{\ell_{\small\textrm{P}}^{2} N_{\parallel}^{2}}{8 \mathcal{L}} \, \left\{1 - \text{exp} \left( - \frac{8 p^{2} \mathcal{L}^{2}}{N_{\parallel}^{2} \ell_{\small\textrm{P}}^{2}} \right)- \text{exp} \left( - \frac{8 q^{2} \mathcal{L}^{2}}{N_{\parallel}^{2} \ell_{\small\textrm{P}}^{2}} \right)+ \text{exp} \left( - \frac{8 \left(p-q\right)^{2} \mathcal{L}^{2}}{N_{\parallel}^{2} \ell_{\small\textrm{P}}^{2}} \right) \right\}. \end{align} \tag{ D8 }$$

Using $\ell_{\small\textrm{P}} N_{\parallel} / 2 \mathcal{L} \ll 1$,

$$\begin{align} \sigma_{z}^{\left(p,q\right)}\approx \sqrt{\frac{\pi}{2}} \, \text{min}\left(p,q\right) \, \ell_{\small\textrm{P}} N_{\parallel}, \end{align} \tag{ D9 }$$

where $\text{min}(p,q)$ chooses the minimum of the two, namely, p and q. Similar simplifications can be carried out to find

$$\begin{align} \sigma_{x}^{\left(p,q\right)}\approx \sqrt{\frac{\pi}{2}} \, \text{min}\left(p,q\right) \, \ell_{\small\textrm{P}} N_{\parallel}. \end{align} \tag{ D10 }$$

As in appendix B, all ξ-type correlations are negligible in comparison to the σ-type correlations in the limit $N_{\parallel} \ll 2 \mathcal{L}/\ell_{\small\textrm{P}}$.

As before, we consider the input beam to be a perfect Gaussian beam at BS with the input electric field $E_{\text{in}}(x,y,0)$ coming in from port A, given by (recalling equation (B11)),

$$\begin{align} E_{\text{in}}\left(x,y,0\right) = \sqrt{\frac{2}{\pi} \,} \, \frac{z_{R}}{W_{0} \sqrt{z_{R}^{2} + z_{0}^{2}}} \text{exp} \left[ \left( -\frac{i k z_{0} W_{0}^{2} + 2 z_{R}^{2} }{2 W_{0}^{2} \left(z_{0}^{2}+z_{R}^{2}\right)} \right) \left(x^{2} + y^{2}\right) \right]. \end{align} \tag{ D11 }$$

Here we recall that W₀ is the beam waist, z₀ is the position of the beam waist, and $z_{R} = \pi W_{0}^{2}/\lambda$.

Using equations (D4), (D7), (D9), (D10), and (D11) in equation (D3), we find

$$\begin{align} & \nonumber \sum_{\substack{q_{1},q_{2}, q_{3},q_{4}} = 1}^{280} \left(\sqrt{R_{\small\textrm{M}}}\right)^{q_{1}+q_{2}+q_{3}+q_{4}-4} \: \overline{M^{2,2}_{\left\{\small\textrm{D,D}\right\};\left\{\small\textrm{C,C}\right\}} \left( \boldsymbol{R}_{\small\textrm{D,D}};\boldsymbol{R}_{\small\textrm{C,C}} \right)}\\ & \nonumber \quad = \sum_{\substack{q_{1},q_{2}, q_{3},q_{4}} = 1}^{280} \left(\sqrt{R_{\small\textrm{M}}}\right)^{q_{1}+q_{2}+q_{3}+q_{4}-4} \: |\mu_{1}\left(\mathcal{L},y,z\right)|^{4} \, e^{- 4 i \varphi_{\mathrm{off}}}\nonumber\\ & \qquad \times \text{exp} \left[ - \sqrt{\pi / 2} \, \mathcal{K} \mathcal{L} \left\{\sum\limits_{j = 1}^{4} q_{j} + 2 \, \text{min}\left(q_{1},q_{2}\right) + 2 \, \text{min}\left(q_{3},q_{4}\right) \right\} \right]. \end{align} \tag{ D12 }$$

$$\begin{align} & \quad = \, |\mu_{1}\left(\mathcal{L},y,z\right)|^{4} \, e^{- 4 i \varphi_{\mathrm{off}}} \, \text{exp}\left( -2 \sqrt{2 \pi} \, \mathcal{K} \mathcal{L} \right) \mathcal{R}_{1}^{2}, \end{align} \tag{ D13 }$$

where

$$\begin{align} \nonumber \mathcal{R}_{1} & = \left[ \frac{1 - \left(\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}} \right)^{280}}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)\left(1-\sqrt{R_{\small\textrm{M}}} e^{- 3 \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} \right. \\ & \left. \quad -\ \: \frac{\left(\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}} \right)^{280} \left[1 - \left(\sqrt{R_{\small\textrm{M}}} e^{-3 \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}} \right)^{280}\right]}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)\left(1-\sqrt{R_{\small\textrm{M}}} e^{- 3 \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)}\right. \nonumber\\ & \left. \quad +\ \sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi / 2} \, \mathcal{K} \mathcal{L}} \, \frac{\left[1 - \left(R_{\small\textrm{M}} e^{-2\sqrt{2 \pi} \, \mathcal{K} \mathcal{L}} \right)^{280}\right]}{\left(1-R_{\small\textrm{M}} e^{-2\sqrt{2 \pi} \, \mathcal{K} \mathcal{L}}\right)} \right. \nonumber\\ & \left. \quad \times \: \left\{\frac{1}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{- \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} - \frac{e^{-\sqrt{2 \pi} \, \mathcal{K} \mathcal{L}}}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{- 3 \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} \right\} \right]. \end{align} \tag{ D14 }$$

We also recall that the definition of $\mu_{1}(\mathcal{L},y,z)$ is given in equation (B13) and $\mathcal{K} = k^{2} \Gamma_{\small\textrm{S}} \ell_{\small\textrm{P}} N_{\parallel}$.

Using arguments similar to those used in obtaining equation (D13), we can compute the other terms in equation (D2) to obtain $\overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R}_{0};\boldsymbol{R}_{0}\right)}$. Defining

$$\begin{align} \nonumber \mathcal{R}_{2} & = \left[ \frac{1 - \left(\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}} \right)^{280}}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)\left(1-\sqrt{R_{\small\textrm{M}}} e^{\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} \right. \nonumber\\ & \left. - \: \frac{\left(\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}} \right)^{280} \left[1 - \left(\sqrt{R_{\small\textrm{M}}} e^{\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}} \right)^{280}\right]}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{-\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)\left(1-\sqrt{R_{\small\textrm{M}}} e^{\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} \right. \nonumber\\ & \left. +\ \sqrt{R_{\small\textrm{M}}} \, \frac{\left[1 - \left(R_{\small\textrm{M}}\right)^{280}\right]}{\left(1-R_{\small\textrm{M}}\right)} \: \left\{\frac{e^{-\sqrt{\pi / 2} \, \mathcal{K} \mathcal{L}}}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{- \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} - \frac{e^{\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}}{\left(1-\sqrt{R_{\small\textrm{M}}} e^{\sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)} \right\} \right], \end{align} \tag{ D15 }$$

$$\begin{align} \mathcal{R}_{3}& = \frac{1-\left(\sqrt{R_{\small\textrm{M}}}\right)^{280}}{1-\sqrt{R_{\small\textrm{M}}}}, \quad \text{and} \end{align} \tag{ D16 }$$

$$\begin{align} \mathcal{R}_{4}& = \frac{1-\left(\sqrt{R_{\small\textrm{M}}} \, e^{- \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}\right)^{280}}{1-\sqrt{R_{\small\textrm{M}}}\, e^{- \sqrt{\pi/2} \, \mathcal{K} \mathcal{L}}}, \end{align} \tag{ D17 }$$

we find

$$\begin{align} \nonumber \overline{M^{2,2}_{\mathrm{out}} \left(\boldsymbol{R}_{0};\boldsymbol{R}_{0}\right)} & = \, |\mu_{1}(\mathcal{L}, y,z)|^{4} \,T_{\small\textrm{M}}^{2} \, \Bigg( \mathcal{R}_{3}^{4} - 2 \text{exp} \left[ - \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right] \mathcal{R}_{3}^{2} \mathcal{R}_{4}^{2} \text{cos} 2 \varphi_{\mathrm{off}} \\ & \quad + \frac{1}{2} \text{exp} \left[ - 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right] \left\{\text{exp} \left[ 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right] \mathcal{R}_{2}^{2} \right.\nonumber\\ & \quad \left. +\ \text{exp} \left[ - 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right] \mathcal{R}_{1}^{2} \text{cos} 4 \varphi_{\mathrm{off}} \right\} \Bigg). \end{align} \tag{ D18 }$$

Using equations (D18), (23) and (24), we find

$$\begin{align} \mathrm{Cov}_{i,j}\left(P_{\mathrm{out}}\right)& = \frac{P_{0}^{2}}{4} \,T_{\small\textrm{M}}^{2} \, \left(\frac{e^{- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L}}}{2} \left\{e^{2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} } \, \mathcal{R}_{2}^{2} + e^{- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L}} \, \mathcal{R}_{1}^{2} \text{cos} 4 \varphi_{\mathrm{off}} \right\} \right.\nonumber\\ & \quad \left.-\ e^{- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} } \mathcal{R}_{4}^{4} \text{cos}^{2} 2 \varphi_{\mathrm{off}} {}\right). \end{align} \tag{ D19 }$$

Using equation (27) and carrying out the cosine transform, we obtain the PSD from the above expression as,

$$\begin{align} S\left(\,f\,\right) = \mathcal{A} T_{\small\textrm{M}}^{2} \left[ \mathcal{R}_{2}^{2} + \, \mathcal{B}^2 \mathcal{R}_{1}^{2} \text{cos} 4 \varphi_{\mathrm{off}} - 2 \mathcal{B} \mathcal{R}_{4}^{4} \text{cos}^{2} 2 \varphi_{\mathrm{off}} \right]. \end{align} \tag{ D20 }$$

Here,

$$\begin{align} &\mathcal{A} = \frac{1}{4 \pi}\left(\frac{\lambda}{8 \pi \varphi_{\mathrm{off}}}\right)^{2} \, \Theta\left(\,f\,\right), \end{align} \tag{ D21 }$$

$$\begin{align} & \mathcal{B} = \text{exp}{\left(- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L} \right)}. \end{align} \tag{ D22 }$$

To first-order in Planck length, i.e. $\mathcal{B} = 1- 2 \sqrt{2 \pi} \mathcal{K} \mathcal{L}$, the PSD becomes

$$\begin{align} S\left(\,f\,\right) \approx \frac{\mathcal{K} \mathcal{L}}{2 \sqrt{2 \pi}} \left(\frac{\lambda T_{\small\textrm{M}}}{2 \pi}\right)^{2} \, \left(\frac{1- \left(\sqrt{R_{\small\textrm{M}}}\right)^{280}}{1- \sqrt{R_{\small\textrm{M}}}}\right)^{4} \Theta\left(\,f\,\right). \end{align} \tag{ D23 }$$

Here we have used the fact that $\varphi_{\mathrm{off}}\ll 1$, as in the setups considered $\varphi_{\mathrm{off}}\approx 10^{-5}$ rad (see table 1). To first-order in Planck length, we have taken $\mathcal{R}_{1}$, $\mathcal{R}_{2}$, and $\mathcal{R}_{4}^{2}$ to be approximately equal to $\mathcal{R}_{3}^{2}$.

Attribute (1) is an important choice made on how the spacetime metric fluctuates. If we choose another metric, for instance,

$$\begin{align} \nonumber \widetilde{g}^{00} & = -1, \quad \widetilde{g}^{11} = n^{2}\left(x,y,z,t\right),\\ \widetilde{g}^{22} & = \widetilde{g}^{33} = 1, \end{align} \tag{ E1 }$$

with all off-diagonal elements set to zero. Here $n(x,y,z,t)$ is a random field in the 4-dimensional spacetime with $n(x,y,z,t) = n_{0}(y,z,t) n_{1}(x)$. The inequalities listed in (8a) are used in what follows. The Assumptions (i) and (ii) are suitably modified according to the chosen form of $n(x,y,z,t)$ and $E_{i}(x,y,z,t)$ in a straightforward manner. Now the SDEs are listed below for three different possible directions of light propagation.

Case X: If light propagates along x-axis and if the electric field is of the form $E_{i}(x,y,z,t) = E_{1}(x) E_{0}(y,z,t) e^{-i (k x - \Omega t)}$, the SDEs are

$$\begin{align} \frac{1}{c} \, \frac{\partial E_{0}}{\partial t} & = -\frac{i}{2 k} \left( \frac{\partial^{2} E_{0}}{\partial y^{2}} + \frac{\partial^{2} E_{0}}{\partial z^{2}} \right) + \frac{i k}{2} \left(n_{0}^{2}-1\right) E_{0}, \:\: \end{align} \tag{ E2 }$$

$$\begin{align} \frac{\partial E_{1}}{\partial x} & = \frac{i k}{2} E_{1} -\frac{i k}{2 n_{1}^{2}} E_{1}. \end{align} \tag{ E3 }$$

In this case, we clearly see that there is a phase accumulated as light propagates in time t. However, due to the similarity in the structure of the SDEs to those discussed in section 3, we can see that the phase difference becomes independent of metric fluctuations along x.

Case Y: If light propagates along y-axis and if the electric field is of the form $E_{i}(x,y,z,t) = E_{1}(x) E_{0}(y,z,t) e^{-i (k y - \Omega t)}$, the SDEs are

$$\begin{align} \frac{\partial E_{0}}{\partial y}+ \frac{1}{c} \, \frac{\partial E_{0}}{\partial t} & = -\frac{i}{2 k} \frac{\partial^{2} E_{0}}{\partial z^{2}} + \frac{i k}{2} n_{0}^{2} E_{0}, \:\: \end{align} \tag{ E4 }$$

$$\begin{align} \frac{\partial^{2} E_{1}}{\partial x^{2}} & = -\frac{k^{2}}{n_{1}^{2}} E_{1}. \end{align} \tag{ E5 }$$

Equation (E4) points to a phase fluctuation along y and in time t. In addition to this, $k^{2}/n_{1}^{2}$ in equation (E5) suggests a fluctuating wavelength. This agrees with the earlier work (see, for instance, [11]) that point to both fluctuations in phase difference between the arms of the interferometer and fluctuations in the wavelength of the light. However, the SDEs corresponding to this case are not readily solvable analytically, in contrast to those considered in section 3. So, in this preliminary investigation, we confine our computations to the metric assumed in Attribute (1).

Case Z: If light propagates along z-axis, the results are identical to that obtained in Case Y.

As expected, the choice of the metric is crucial on how the spacetime fluctuations affect the observables of interest. So, this investigation reveals the need to further examine the precise effects of the different metric fluctuations in detail.