Searching for stochastic background of ultra-light fields with atomic sensors

We propose a cross-correlation method for the searches of ultra-light fields, in particular, with a space network of atomic sensors. The main motivation of the approach is cancellation of uncorrelated noises in the observation data and unique pattern the fields leave on the cross-spectrum, depending on their nature (i.e., scalar, vector or tensor). In particular, we analytically derive a dependence of the cross-spectrum on the angle between two pairs of detectors. We then confirm obtained angular curves with a numerical simulation. We imply application of the method to the detection of dark matter and gravitational waves.


I. INTRODUCTION
Search for new fields has been the central part of modern high-energy physics. While most of the effort is concentrated in the mass range of GeVs and above [1], where the new physics can be tested with particle accelerators and cosmic rays, relatively little is done in the sub-eV region. In that region, most active are the axion studies [2], usually spanning the interval of µeV-meV. Fields with even lower masses have gained more attention only very recently, mostly, as possible candidates for the dark matter (DM), see Refs. [3,6] for an overview. Large Compton wavelengths for such fields (λ > 1 m for m < 1µeV) require a methodology beyond the traditional particle physics but can be tested with atomic physics experiments and/or large scale experiments in space. One example is a comparison of frequencies of two atomic clocks being affected by slowly oscillating scalar field background. Another would be a test of the weak equivalence principle with two isotopes in free fall, as an additional acceleration might be induced by a new ultra-light vector field. Example of a large-scale experiment would be a fifth-force type of searches, seeking for anomalies in the trajectories of celestial bodies and spacecrafts [7] or searches for timing anomalies in the GPS data [8]. For an overview of possible spatial configurations of ultra-light fields, see Ref. [6] and Refs. therein. We recall that the current literature provides the following most common configurations of the fields: waves with the frequency at the field mass [9][10][11], clumps (e.g., topological defects [12][13][14]), caustics [15] and simple static distributions [16].
In this paper, we consider a new case: a stochastic background of waves of light fields, which we propose to measure by means of a network of precise atomic sensors. Some of the past methods assumed that the entire energy density of DM is carried by a monochromatic (or quasi-monochromatic [17][18][19][20]) wave with a frequency fixed at m φ /2π. However, if the total energy density is distributed over a range of frequencies, then the limits summarized in Ref. [6] will be significantly reduced. Therefore, it makes sense to put limits not only on the DM couplings, but also on the spectrum of DM excitations. Stochastic background of waves may appear in the context of Bose-Einstein Condensate (BEC) and superfluid dark matter models [21][22][23]. Technically, it may also be easier to search for a stochastic background of waves rather that single waves, because one does not have to form a bank of templates and the signal-to-noise ratio can be substantially improved by increasing the observation time (we discuss it later in the paper). We first analyze the cases of scalar and vector fields in the context of DM, and then also show how similar methodology can be used for tensor fields in the context of gravitational waves.

II. STOCHASTIC SCALAR FIELDS
In this section, we consider a new hypothetical scalar field, which we, for convenience, associate with DM to be able to refer to the existing literature. The scalar field φ with mass m φ will be coupled linearly to the Standard Model (SM) operators (see, e.g., dilaton dark matter studies [10,11,24]), where F µν and G µν standard electromagnetic and gluon field strength tensors, β 3 is the beta-function of the gauge coupling g 3 , γ m f is the anomalous dimension of the fermion (electron, u-, d-quark) mass operator and we use natural units, = c = 1. Parameters Λ a have dimension of mass and play a role of (unknown) inverse coupling constants. For other possible couplings, see Refs. [3][4][5][6]. The chosen interaction Lagrangian, Eq. (1), introduces local changes in values of fundamental constants, such as [24] where α is the fine-structure constant and Λ QCD is the QCD scale. These variations can be studied by the changes in the atomic clock frequency [25] ν where m q = (m u + m d )/2 and exponents K a are tabulated for the most common types of atomic clocks, see Refs. [10,25,26] and references therein. As an example, for 133 Cs, K α = 2.83, K q = 0.07, K e = 1. For optical clocks, only K α is nonzero, i.e., the clocks are only sensitive to the changes in the fine structure constant. Using (2), we obtain For our study, it will be crucial that the fractional frequency variation depends linearly on the scalar field strength φ, allowing us to express the Fourier transform of the frequency variation through the Fourier transform of φ. This, together with the general requirements, such as the gauge-and Lorentz-invariance, is the main motivation behind the chosen form of the Lagrangian (1). There exist strict experimental limits on the parameters Λ a based on this (and similar) methods, see Refs. [6,9,11]. However, these limits are based on the assumption of DM wave being monochromatic and can be relaxed in the case of a more general spectrum of DM waves. We refer the reader to the Appendix A for a relevant example. It is important to emphasize that the limits on Λ a should be always considered in the context of chosen DM configurations (e.g., waves, lumps, constant field, etc.). Since stochastic DM backgrounds were not considered before, we do not contrast possible Λ a sensitivities of our method to the projected sensitivities or ruled out regions of parameter space attributed to other methods.
We consider a triangle-shaped configuration of atomic clocks with frequencies ν i , i = 01, 02, 1, 2, placed at positions x i , see Fig. 1(a). Each clock responds to the scalar field by a shift in the frequency, δν i such as δν i /ν i = κ i φ, cf. Eq. (4). We will measure the difference in the relative frequency shifts, X 1 (t) = δν 01 /ν 01 − δν 1 /ν 1 and X 2 (t) = δν 02 /ν 02 − δν 2 /ν 2 . In a real experiment, one compares absolute clock frequencies in pairs, so we will assume that the clocks are identical in each pair, and X i (t) represents the relative frequency difference between clocks in a pair with the gravitational redshift difference taken into account. We have chosen two pairs of atomic clocks instead of three identical clocks because the common reference clock would introduce a possible unwanted correlated noise that can be mistaken for the signal. We will be interested in the cross-spectrum between X 1 (t) and X 2 (t) and derive a universal dependency of this cross-spectrum on the angle between the pairs of atomic clocks. The approach for the derivation is similar to the calculation of the Hellings-Downs curve [27] for pulsar-timing arrays used in the gravitational wave searches. The scalar field in consideration does not have to be DM and can be any neutral light scalar field obeying the energy density limit on the hidden matter content in the solar system (or, more generally, in the confining volume of the performed experiment), see, e.g., Ref [32]. In what follows, we assume that the data collection time T is finite but large, giving the frequency resolution of the signal ∆f ∼ 1/T . The finite-time Fourier transform for all time-dependent quantities in this article is given byX We will be replacing the integration limits by infinity, whenever it does not lead to a confusion. The power spectral density S X (f ) of X for one of the clock pairs is defined via S X (f ) = |X(f )| 2 and has the property (Parseval's theorem) We have chosen a two-sided version of the power spectral density to make it easier to change the order of integration, when needed. Notice that the given form of Fourier transform (together with the property of stationarity) leads to where we defined the finite-time delta function as If the two frequencies coincide, then δ(0) = T cancels the factor T in the denominator and we recover the definition of the power spectrum. Let us consider a stationary isotropic background of scalar field waves. We represent the scalar field as where the dispersion relation, in natural units, is 2πf = (m 2 φ + k 2 ) 1/2 for massive matter waves or f = const · k for gapless modes (such as phonons). Next, we introduce the cross-spectrum for two pairs of clocks, We will be using two identities (the second one comes from the isotropic property of the background), to calculate the cross-spectrum from the power spectrum of the scalar field, where the response function R c (f, ζ) is given by By chosing the appropriate coordinate system (see, e.g., Ref. [29]), one can integrate each term in the sum and obtain where The response function R c (f, ζ) relates the power spectrum of the scalar field fluctuations to the cross-spectrum of the measured signals, If the scalar field φ is identified with DM, then its power spectrum density, S φ (f ), has the following integral characteristic, where ρ DM is the average DM energy density in the vicinity of the experiment. A common model for the DM distribution in the Milky Way is a non-rotating isothermal spherical halo with velocities of the DM objects following Maxwell distribution. For a solar system experiment, it is believed that ρ DM ≈ 0.4 GeV/cm 3 [30] and that the DM objects are moving with (virial) velocities of v b ≈ 270 km/s and the velocity dispersion is δv b ≈ v b [31]. The upper bound on ρ DM is currently set to ρ DM < 10 5 GeV/cm 3 , based on positional observations of planets and spacecraft [32]. Since the exact shape of S φ (f ) is not known a priori, from the practical point of view, it makes sense to eliminate it by normalizing the crossspectrum with the power-spectrum of the signal from one pair of atomic clocks, Fig. 2. It is clear from the plot that the sensitivity for a pair of identical atomic clocks reaches its maximum for the separation D i being larger than the scalar field wavelength. After the normalization of the cross-spectrum, the final result becomes There are several natural limiting cases that can greatly simplify this expression due to the properties of the sinc function, where ζ is the angle between ( x 1 − x 0 ) and ( x 2 − x 0 ). The second limit is geometrically restricted to a situation when D 1 ≈ D 2 and ζ ≈ D 12 /D i 1 or, in other words, when 2πD 12 is much smaller that the wavelength of the scalar field excitation. It is important to notice that expressions in Eq. (18) do not depend on the DM wave-vectors, frequencies and masses explicitly. Here we used the property sinc(a) ≈ 0 when a 1 and sinc(a) ≈ 1 − a 2 /6, when a 1. Another nontrivial observation is that the first two limits give us nonzero constants [29], while an uncorrelated noise would have given a vanishing cross-correlator.
Eq. (18) is the main result of this section. One can use it to identify the presence of a scalar signal in the cross-correlation data, either directly or by constructing an optimal filter to extract the signal from noise. In the described procedure, the frequency of the signal is bounded from below by several factors. First, the frequency should be larger than m φ /(2π), due to the dispersion relation for the scalar field. Second, it should be larger than the inverse data collection time, f 1/T , due to the property of the Fourier transform. Finally, it is limited by the time T ζ ∼ ∆ζ/ζ by which the angular distance between sensors changes by the value equal to the uncertainty in the angle ζ, so f ζ /∆ζ. Value of ∆ζ should be chosen such that there is enough recorded data for the extraction of the power spectrum density at the given frequency.
The maximal frequency is limited by the Nyquist-Shannon-Kotelnikov theorem, f max = 1/(2∆t), where ∆t is the time interval between consecutive measurements. Value of ∆t depends on the averaging interval giving desirable performance of the atomic sensors (typically, ∆t ≥ 1 s for high-performance atomic clocks). Frequencies larger than f max can be studied with the help of the temporal aliasing effect [6], however, such analysis goes beyond the scope of this paper.
The method described here is general in nature and can be also applied beyond the use of atomic sensors. For instance, one could search for the long-wavelength curve, Eq. (18), using instruments for the gravitational wave studies. First, one could analyze pulsar timing array data. When passing through a pulsar, the DM field can alter the neutron mass, size of the pulsar and hence its moment of inertia [13]. This will lead to the variation of the pulsar rotation period. For the sensitivity estimates in the case of a single monochromatic DM wave, see Ref. [3]. Another method of DM detection with pulsar timing is presented in Ref. [33]. Second, the DM wave could lead to the periodic displacement of test masses in LIGO and LISA experiments (for instance, LIGO mirrors), see the sensitivity estimates for Λ g and a single monochromatic scalar wave in Ref. [10]. Third, one could use binary pulsars for studies of stochastic backgrounds, see, e.g., Refs. [34,35].
Dark matter "wind" DM direct detection experiments should take into account the possibility of DM moving with a constant velocity v b with respect to the observer (DM "wind"), see Ref. [31] and refs. therein. Since the solar system moves through the galactic halo of the DM, such velocity may be given by the velocity of the Sun with respect to the galactic rest frame. If the solar system has its own DM halo, then the velocity of DM with respect to the near-Earth experiment would be given by the Earth orbital speed. In this section, we consider the effect of the constant motion of the configuration of detectors through the DM rest frame without making assumptions on the particular direction and absolute value of such constant velocity. Our goal is to see if there is a change in the angular part of cross-spectra when the velocities of DM waves are shifted by a constant vector v b . Another approach to the DM "wind" detection with atomic clocks is proposed in Ref. [17].
For scalar excitations, due to the properties of the Fourier transform and Eq. (6), where is the Dopplershifted frequency of the stochastic excitations. The cross-spectrum will be given then by where we consider the configuration given by Fig. 1(a), and x 0 = 0, for the sake of simplicity. After the integration, the expression becomes giving us familiar results (with the only difference of k being replaced by k ), To reiterate the strategy, if the normalized cross-spectra give a constant (e.g., 1 or 1/2, for identical clocks) or cosine dependence on the angle between two pairs of atomic clocks, then the measured signal may be dominated by the presence of a new scalar field. In our derivation, we ignored the contribution of noise, assuming it to be weak. If the noise is strong, however, then one has to apply a different methodology, where the obtained analytic results, Eq. (18) or Eq. (22), are used in the combination with a statistical analysis for the extraction of a weak signal from noise. We describe this approach in the next section.
Signal-to-noise ratio and statistical inference In this section, we provide the signal-to-noise ratio (SNR) for stochastic background measurements, an optimal filter, as well as the probability of the presence of the signal in the noise. This will allow us to estimate measurement parameters (such as the measurement duration) to achieve a given measurement sensitivity. Our derivations follow the gravitational wave detection methods [36][37][38]; for the basics of the signal extraction from noise see, e.g., Ref. [39]. We focus on the case of scalar fields, the generalization on other cases is straightforward. We consider the measured observable as a sum of a (weak) signal φ(t) and a noise n(t), both having zero expectation values. For each pair of atomic clocks the noise is characterized by the (one-sided) power spectrum S ni , and the noises are assumed to be uncorrelated and stationary. In order to solve the problem of optimal filtering, we consider the signal to be If the filter Q(t − t ) = δ(t − t ), then the signal is simply S = φ 1 (t)φ 2 (t) . This signal can be, however, buried under noise and we need to construct a filter Q(t) that will maximize SNR, eventually making it large enough at the expense of long observation time T . It will be shown that SNR ∼ √ T for our measurement. The Fourier space representation of the signal is where we assumed that Q(τ ) decays quickly with τ → ±∞. Next, statistical properties of the field amplitudes and the noise can be expressed as and, in addition to the property δ T (0) = T , lead to The noise is defined in the standard way, N ≡ S − S . Assuming the noise contributions dominate the signal, we can write and further obtain where we used a standard way of expressing the 4th moment through the covariances and took into account that Comparing this expression to Eq. (27), one can see that and gives Notice that SNR ∝ √ T and the noise decreases with more data points collected, so, in ideal conditions, an arbitrarily small signal can be extracted from noise for long enough duration of the observation, which is, indeed, the power of the cross-correlation method.
In order to claim detection of a weak signal, one has to assume a certain shape of the spectrum S φ (f ), false alarm rate α and false dismissal rate β. Consider a set of n statistically independent measurements {S i }, each of duration T . In order to test the null hypothesis (no signal of scalar waves in the data), one can form a random variableX which is normally distributed with unit variance, in assumption of a large enough n. The Neyman-Pearson decision criterion (maximizing probability of the detection with fixed alarm rate α) allows us to use it as a test statistic and choose the null hypothesis ifX < X * , where X * is fixed by the choice of the false alarm rate, and reject the null hypothesis otherwise. Here erfc −1 (x) is the inverse of the complementary error function. In other words, one can claim the presence of the DM signal of unknown amplitude if √ nŜNR ≤ √ 2 erfc −1 (2α). At the next step, after the null hypothesis is rejected, i.e., assuming the signal is present in the data, the theoretical SNR required for the detection of the DM background in at least (1 − β) × 100% of measurements is given by [36] √ If the spectrum is flat, S φ (f ) =S φ = const, in order to detect the signal, one would require the following minimal value ofS φ , The same technique can be applied to a network of more that two pairs of clocks, see Ref. [36] discussing it in the context of gravitational wave detectors. For the anisotropic backgrounds and real data complications see, e.g., Ref. [38]. For illustration purposes, let us consider a system of identical clocks in the short-wavelength regime (F (ζ) = 1/2), with the scalar field background being localized around frequency f 0 in a narrow band ∆f and, for the particular statistics α = β = 0.05. Then the mentioned above expressions are simplified to where the power spectrum density is for the 5% false alarm and false dismissal rates. Again, one see the advantage of the large number of measurement sessions and their long duration, with fixed levels of noise.

Numerical simulations
In this section, we check validity of some of our results numerically. Necessity of the numerical test is due to the finite number of sources of DM waves, finite measurement time and number of measurement sessions (i.e., statistics). The numerical analysis provides some guidance with respect to the choice of measurement parameters and the impact of imperfections in the stochastic background. We choose N = 1000 sources of DM waves randomly positioned in the sky (Fig. 3) with random amplitudes, phases and random frequencies normally distributed around a given mean value. For the angular positions,  Fig. 4(a) and its spectrum is presented in Fig. 4(b). The cross-correlator for the two sets of data and the cross-spectrum are shown in Figs. 4(c) and 4(d), respectively. The angular curves are obtained within the frequency band given by the spread of the spectrum of the sources and shown in Figs. 5(a) and 5(b). The long-wavelength limit gives the identical result to the analytic formula, Eq. (18), here D 2 = 0.1 D 1 by choice. The short wavelength limit agrees with Eq. (18) within at least one standard deviation.
We also made preliminary measurements of the SNR with the optimal filter derived in the previous section. In the presence of a weak and strong white Gaussian noise, these measurements are consistent with the SNR ∝ √ T dependence. More precise quantitative statement requires additional computing power and will be made elsewhere. The effect of noise suppression with increase of observation time is demonstrated in Fig. 6.

III. STOCHASTIC VECTOR FIELDS
Another popular DM candidate is the massive B − L vector field that exerts an additional force on neutral atomic species with different number of neutrons (e.g., different isotopes). This addtional force can be probed with atom interferometer (AI) measurements in acccelerometer configurations. Such B −L field detection has been mostly discussed in the context of experimental tests of the weak equivalence principle (WEP) and manifests itself in an additional relative acceleration between two particle species a and b [3,40], where g B−L is the coupling constant to the B − L field W , Z a,b and A a,b are the atomic numbers and weights, respectively, m N is the neutron mass, and we also introduced the constant κ for the sake of compactness. It has been proposed to measure such acceleration difference with the dual-species atom interferometers (AI) [3,40]. For the review on principles of atom interferometry and use of AI as accelerometers, see, e.g, Ref. [41]. Let us consider four AI in a triangle-shaped configuration, such that two of them are at the position x 0 = 0 and the other two are at the positions x 1 and x 2 , respectively, see Fig. 1(b). All AI are oriented in the radial direction with respect to the center at x 0 . The quantities we  are interested in are defined by the radial accelerations (along the lines connecting AIs in a pair), whereê i are unit vectors directed along the AIs. We represent the vector field W as whereˆ i (k) are polarization vectors, such asˆ i (k) ·ˆ j (k) = δ ij and the dispersion relation is 2πf = (m 2 W + k 2 ) 1/2 with m W being the mass of the B − L field [57]. We further consider a stochastic, isotropic and unpolarized background of vector waves characterized by the power spectrum density S W (f ), Notice the factor 3 in the denominator, which shows the number of physical polarizations. In analogy with the previous section, we calculate the cross-spectrum, To compute this integral we first fix the orientations of AI in the Cartesian (x, y, z)-coordinates, e 1 =ẑ,ê 2 = sin ζx + cos ζẑ , and then fix the orthonormal basis of polarization vectors in polar coordinates (k, θ, ϕ), 1 (k) = cos θ cos ϕx + cos θ sin ϕŷ − sin θẑ =θ, 3 (k) = sin θ cos ϕx + sin θ sin ϕŷ + cos θẑ =k.
In the short-wavelength approximation, D i , D 12 1/k, we can neglect the exponents in the right-hand side of (41) and obtain R c (f, ζ) = κ1κ2 3 cos ζ. If D i 1/k, D 12 1/k, then R c (f, ζ) = 2κ1κ2 similar to the previous section. Finally, in the long-wavelength limit, D i 1/k, taking into account that x i =ê i and expanding the exponent in series, we get R c (f, ζ) = κ1κ2 9 k 2 D 1 D 2 cos 2 ζ. Considering the power spectrum for each individual pair of AI, we summarize the results of this section: so the strategy of detecting the DM footprint in the data is to search for a cosine or cosine squared modulation of the normalized cross-spectra.
There are several remarks in order. First, if one expects a contribution of the DM "wind", then the result repeats Eq. (46) with k being replaced by k in the limits. Second, if the vector field excitation is gapless, then the number of physical polarizations is reduced to two and the problem resembles the case of electromagnetic stochastic background [29]. Third, if all conventional sources of acceleration were known (if, e.g., the setup was in deep space), can be controlled or are spectrally uncorrelated, then one could consider, for simplicity, four spatially separated single-species AI in configuration of longbaseline gradiometer [42]. In this case, the provided derivation of the angular curve can be repeated after subtraction of the all conventional accelerations. Finally, the triangular measurement configuration for B − L vector field discussed above resembles a two-arm gravitational wave detector such as LIGO and LISA. We speculate that either of them can be used for the B − L vector DM detection because of the differential acceleration between test masses and highly correlated laser noises in each arm.

IV. STOCHASTIC TENSOR FIELDS
In this section, we discuss an application of the method to the searches of isotropic unpolarized stationary gravitational wave (GW) background. Even though the method of cross-correlation of data from several detectors is widely used in the pulsar astronomy community (see, e.g., Refs [36][37][38]), we are not aware of the use of the method in direct detection experiments with atomic sensors. An example of the configuration is depicted in Fig. 1(b), where the atomic clocks are placed in spacecrafts or attached to celestial bodies in free fall. Frequency of an electromagnetic signal (e.g., laser) is locked to the frequency of one of the clocks in the pair and being compared with the frequency of the other clock, see, e.g., proposals in Refs. [43][44][45]. Presence of GW will introduce a relative Doppler shift between the two frequencies. One can also use atom interferometers as gravitational wave detectors [46][47][48][49], with test masses being atoms in excited and ground states. Even though in AI experiments the acceleration is measured rather than velocity, fundamentally, the observables are equivalent to the ones of the atomic clock detectors [50]. We focus, for simplicity, on the case of clock comparison. The frequency difference due to the presence of the GW is denoted X a = δν a /ν a . Our goal, as before, is to extract the angular part of the cross-spectrum, R c (f, ζ). The small perturbation of the metric around the Minkowski metric η µν is given by g µν = η µν + h µν . The metric perturbation h µν in a spatial transverse traceless gauge has h 0µ = 0 and can be represented by the plane wave expansion [36] where f is the frequency of the wave, unit vectorΩ points in the direction of the propagation of the planewave component. Dispersion relation for the gravitational waves is simply 2πf = k. The polarization tensors can be defined through unit vectorsn andm orthogonal toΩ, where in the standard spherical coordinateŝ Ω = (sin θ cos φ, sin θ sin φ, cos θ), n = (cos θ cos φ, cos θ sin φ, − sin θ).
We begin from considering an effect of a single plane wave propagating in directionΩ on a signal sent between atomic clocks separated by distance D a . The unit vectorp (a) is pointing from the observation point to the singnal source. In order to find the Doppler shift X a , one can consider the null vector [51] σ µ (a) = s µ (a) − at the moments when the signal is emmited and received, with the unperturbed value given by s µ (a) = ν(1, −p (a) ). The null geodesics can be found by solving the null condition σ µ(a) σ µ (a) = 0 together with the standard condition σ µ V µ (a) = const (a) , where V µ (a) are the Killing vectors for the perturbed geometry. The final result is given by where we take into account the isotropic nature of the radiation and consider the perturbation amplitudes as functions of light-cone coordinates. The power spectrum density for the induced Doppler shifts, S Xa (f ) = X a (f )X * a (f ) is given by averaging over time, directionsΩ and polarizations of the incoming radiation. It can be shown that S Xa (f ) = R a (f )S h (f ), where the GW power spectral density is defined by h P (f )h * P (f ) = 1 2 S h (f )δ P P and the response function is given by see Fig. 7(a). In a long wavelength limit, 2πf D a 1, the response function reduces to R a (f ) = 4π 2 15 (D a f ) 2 , while in the short wavelength limit it is simply R a (f ) = 1/3. By considering Eq. (53) in the Fourier space, one can derive the cross-spectrum, In order to perform the direction averaging for the radiation background, we chose the z-axis to be alonĝ p (1) andp (2) to have an angle ζ with respect top (1) , similar to the previous section, p (1) = (0, 0, 1),p (2) = (sin ζ, 0, cos ζ).
The angular dependency on ζ in the power spectrum density S c (f, ζ) can be factorized in the shortwavelength limit, D a , D 12 1/k and corresponds to the Hellings-Downs curve [27,52] used in the pulsar timing studies, To calculate this expression, one can neglect exponents in Eq. (55), since their arguments are quickly oscillating functions of spherical angles. When D a 1/k, but D 12 1/k, one will get the same result multiplied by factor 2, as in the previous sections. In the opposite limit, D a 1/k, on can expand the exponents in series and obtain which reproduces the response function of one pair of atomic clocks in this limit, when ζ = 0. One can notice that the angular function above is proportional to the second Legendre polynomial, P 2 (cos ζ), which one can expect from the expansion of the background by spherical harmonics [53,54]. Finally, introducing, again, the angular curve F (f, ζ) ≡ S c (f, ζ)/S Xai (f ) we conclude this section, , D a 1/k, D 12 1/k, We plot the first expression and the third expression (with D 1 = D 2 ) from that equation in Fig. 7(b). Sensitivity of a subsystem of a pair of atomic sensors is discussed in, e.g., Ref. [44].

V. CONCLUSIONS
In this article, we proposed a general method of experimental measurements of ultra-light (or even massless) fields with atomic sensors. Such fields include various dark matter candidates (e.g., a dilaton or B-L field) and the gravitational field. We assumed that the fields form an isotropic stationary background and are being measured by two pairs of detectors in a triangular-shaped configuration (with or without variable angle). Depending on the nature of the fields (scalar, vector or tensor) they will leave a characteristic angular dependence on the cross-correlation of the data obtained from each pair. If the instrumental noise is weak, then the angular dependence of the cross-spectrum can be measured directly. If the noise is stronger than the expected signal, then the analysis should be done in a different way, that involves a statistical analysis with optimal filtering (we used the frequentist approach in our paper). Such analysis allows to recover the signal from a strong noise by performing long time observations. If the signal is not seen but is expected to be present, then one could put limits on the parameters of the model (such as DM couplings or characteristics of S φ ). This, however, would require further assumptions on the shape of power spectrum density of the stochastic background (see, e.g., Eq. (36) and Appendix A), which has to be motivated by the underlying theory of DM or GW generation. Comprehensive review of such theories is beyond the scope of this paper.