Upper limit on scalar field dark matter from LIGO-Virgo third observation run

If dark matter is a light scalar field weakly interacting with elementary particles, such a field induces oscillations of the physical constants, which results in time-varying force acting on macroscopic objects. In this paper, we report on a search for such a signal in the data of the two LIGO detectors during their third observing run (O3). We focus on the mass of the scalar field in the range of $10^{-13}-10^{-11}~{\rm eV}$ for which the signal falls within the detectors' sensitivity band. We first formulate the cross-correlation statistics that can be readily compared with publically available data. It is found that inclusion of the anisotropies of the velocity distribution of dark matter caused by the motion of the solar system in the Milky Way Galaxy enhances the signal by a factor of $\sim 2$ except for the narrow mass range around $\simeq 3\times 10^{-13}~{\rm eV}$ for which the correlation between the interferometer at Livingston and the one at Hanford is suppressed. From the non-detection of the signal, we derive the upper limits on the coupling constants between the elementary particles and the scalar field for five representative cases. For all the cases where the weak equivalence principle is not satisfied, tests of the violation of the weak equivalence principle provide the tightest upper limit on the coupling constants. Upper limits from the fifth-force experiment are always stronger than the ones from LIGO, but the difference is less than a factor of $\sim 5$ at large-mass range. Our study demonstrates that gravitational-wave experiments are starting to bring us meaningful information about the nature of dark matter. The formulation provided in this paper may be applied to the data of upcoming experiments as well and is expected to probe much wider parameter range of the model.


Introduction
Although there is little doubt for the existence of dark matter, it's nature is still elusive.There is considerable uncertainty in the mass of constituting elementary particle or object: the allowed mass range is 10 −19 eV m 10 60 eV.Many models of dark matter such as Weakly Interacting Massive Particles (WIMPs), axions, and primordial black holes (PBHs) have been proposed in the literature.Different models predict different observational signals and various experiments suited for detecting specific type of dark matter have been conducted.Yet, no robust detection is reported (for a review of dark matter, see e.g.[1]).In light of this situation, it is important to consider utilizing even available experiments whose main target is not dark matter.
A dark matter model we consider in this paper is a light scalar field dark matter weakly interacting with elementary particles in the standard model of particle physics.Such a scalar field may be realized in string landscape (see e.g., [2,3]), although seeking the fundamental origin is not our focus and we remain at the phenomenological level in this paper.If the mass of the boson is m 10 eV, the dark matter behaves as a classical wave [4].Such a wave induces time variation of the physical constants and leads to detectable signals such as in the atomic transitions [3,[5][6][7][8][9] and the laser interferometers [3,[10][11][12][13][14]. Searches for vector field dark matter with data from advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO) [15] have already been carried out [16,17].The scalar wave can also manifest itself as the violation of the equivalence principle [18][19][20][21] and the fifth force [22].Since the amplitude of the scalar wave increases toward higher redshift, the scalar dark matter also affects the abundance of light elements produced by the big bang nucleosynthesis and the temperature anisotropies of the cosmic microwave background [23].
If the mass is in the range (10 −13 eV, 10 −11 eV), the frequency of the signal falls in the range relevant to the ground laser interferometers.In [3,12], it was suggested that the laser interferometers can probe the scalar field dark matter by measuring the motion of mirrors caused by the scalar field.When the scalar field interacts with elementary particles, the mass of the mirror varies in responding to the value of the scalar field at the mirror's location.As a result, the mirror under the influence of the scalar field accelerates toward a direction along which the mirror minimizes its energy (= mass).Since the mass of an atom dominantly originates from the sector of the strong interaction, the scalar force that accelerates the mirror is most sensitive to the coupling to the gluon fields.
In the non-relativisitc limit of the mirror's motion, the force appears only when the scalar field φ varies in space, which means that the scalar force is parallel to ∇φ.Such a force causes the modulation of the round-trip time of light between the front and the end mirrors, which eventually leads to the detectable signal in the interferometers.There are two physical effects that contribute to such modulation of the round-trip time: i) finiteness of speed of light which produces signal even if the front and the end mirrors are perfectly synchronized (denoted by δt 1 in [12]) and ii) difference of the displacement between the front and the end mirrors due to spatial variation of the scalar field (denoted by δt 2 in [12]).For aLIGO, the former/(latter) becomes dominant if the frequency of the scalar field oscillation is larger/(smaller) than 20 Hz.
The papers [3,12] give expected upper limit on the magnitude of interaction between the scalar field dark matter and elementary particles for the laser interferometers such as aLIGO.To the best of our knowledge, there is no study which derived the upper limit by confronting the predicted signal with the real data obtained by the currently operating laser interferometers.In this paper, we report the result of our analysis using the real data of two LIGO detectors taken during their third observation run.In the analysis, we take into account the stochastic nature of the scalar field, velocity dispersion of dark matter, and the anisotropic distribution of velocity of dark matter due to rotation of the solar system in the Milky Way Galaxy.Formulation to compute the signal that includes the above effects and can be readily compared with the data released by LVK Collaboration is developed in Sec. 2. The upper limits based on our analysis are given in Sec. 3.
Before closing this section, it should be mentioned that the scalar field not only exerts the force but also changes the size of the bodies such as beam splitters and mirrors.The latter phenomenon occurs due to the variation of the atomic size ( Bohr radius) which depends on the fine structure constant and the electron mass.It was recognized in [13] that the time variation of the beam splitter produces detectable signal which is more prominent in GEO600 interferometer than aLIGO.The reason why GEO600 is more powerful than aLIGO is because GEO600 has a direct sensitivity to the signal, while in aLIGO, the strength of the signal is attenuated by a factor of arm cavity finesse.Thus, GEO600 is suited for probing the couplings of the scalar field to the photons and electrons and plays a role complementary to aLIGO which is sensitive to the coupling to the gluons.In [24], the signal was searched for in the data of GEO600 and the upper limits on the couplings to photons and electrons were derived.

Formulation of the signal
In this section, we first give the Lagrangian of the model in which the dark matter is a classically oscillating massive scalar field weakly interacting with elementary particles such as electrons, quarks, photons, and gluons, and give a brief overview of how dark matter yields signal in the laser interferometers.We then develop the formulation which enables to directly compare the predicted signal with public data provided by the LVK collaborations.

Model considered in this paper
The mathematical framework of the model we consider is given in [18].In this model, the scalar field φ which comprises all dark matter has interactions with elementary particles.The Lagrangian density for the scalar field is given by Here κ ≡ √ 4π M pl , and GeV is the Planck mass.The first, second, and third term in L int represents interaction between φ and electromagnetic field, gluon fields, and fermions including electrons, up quarks, and down quarks, respectively.Magnitude of each interaction is parametrized by dimensionless constant d e , d g , d m i which we treat as free parameters.One interesting example in which these five parameters are parametrized by a single parameter as d g = d i , d e = 0 (i = e, u, d) is realized if the scalar field has the nonminimal coupling to gravity in the Jordan frame.Description of this model is given in the Appendix.β 3 is the QCD beta function and γ m i is the anomalous dimension.The interaction Lagrangian given above is defined at some high energy scale.We can translate the effect of L int into equivalent but more physically relevant representation by considering the renormalization-group flow: Here, ϕ ≡ κφ is the scalar field normalized by the Planck scale.α is the fine-structure constant, Λ 3 is the QCD mass scale, and m u,d are the quark masses defined at the QCD mass scale.Due to the presence of ϕ, these physical constants deviate from those in the absence of ϕ.We assume that the mass m is in the range 10 −13 eV m 10 −11 eV.In this case, the φ field as dark matter behaves as a classically oscillating wave since the corresponding occupation number is huge (e.g., [4]).In Fourier domain, frequencies are sharply localised at f As it is argued in [18], coupling to the fermion kinetic term can be eliminated by the field redefinition.To see this, let us suppose that the φ couples to the kinetic term of the fermion as where we leave f (φ) unspecified to keep generality of our argument.The second term is necessary to make the kinetic term be Hermitian.Then, by changing the variable as ψ → e − 1 2 f ψ, the kinetic term in terms of the new field becomes and the coupling f (φ) has been completely absorbed into the field redefinition.Since mass of an atom depends on physical constants, it depends on ϕ through Eqs.(2.3)-(2.6).Denoting the mass of an atom by m A (ϕ), m A in the small ϕ limit is approximately written as [18] and d * g is defined by where A, Z is the mass number and the atomic number, respectively, and As it is clear, d * g represents a part independent of the atomic species and the remaining terms are dependent on the atomic species.Thus, any model belonging to a class satisfying d m = d me = d g , d e = 0 does not violate the weak equivalence principle.A non-minimal coupling model described in the Appendix belongs to this class.
The time variation of the mass of atoms due to the oscillations of φ causes the time variation of the mass M (ϕ) of a macroscopic object.The equation of motion of the macroscopic object under the influence of ϕ is given by (e.g., [12]) where α A is α A averaged over different types of molecules constituting the macroscopic object.The mirrors installed in the interferometers are purely made of silica (=SiO 2 ).In this case, we have It is seen that the mirrors are mostly sensitive to the coupling to the gluon fields d g and the couplings to electrons and photons are suppressed.This is because the mass of nucleons is mainly determined by the strong interaction.

Scalar field in Galaxy's center-of-mass frame
Inside the Milky Way Galaxy, dark matter is virialized.The ocsillating scalar field at the location of the solar system in the galaxy's center-of-mass frame may be modeled as stochastic wave with stationary and isotropic velocity distribution with velocity dispersion v * 220 km/s.In this coordinate system, the scalar dark matter may be expressed as the superposition of plane waves as Here, η is defined by and Ω is the unit vector representing the propagation direction of the plane wave u(f, Ω).
From the reality of φ, we have u * (f, Ω) = u(−f, Ω).From the stationarity and isotropy condition, we assume that the ensemble average of the stochastic variable u(f, Ω) can be written as Here, P φ (f ) is the powerspectrum of φ whose shape may be determined by comparing the energy density of the scalar field evaluated in the wave picture with that in the particle picture in the following way.Let us first evaluate the energy density in the wave picture.
The energy density of the scalar field is given by (2.17) Using Eq. (2.14), the ensemble average of ρ φ becomes Next, we evaluate the same quantity in the particle picture.In this picture, we interpret the stochastic scalar field as a collection of randomly moving particles with mass m, for which we can introduce a notion of the distribution function g(v) in the velocity space.As it is usually done, we assume that the distribution of the particles obeys the Boltzmann distribution as

.19)
Here v = |v|, and v * is the velocity dispersion of dark matter.The normalization constant is fixed by requiring that the integral of g(v) coincides with the local number density of dark matter particles; Therefore, we have  In order to compare, we change the integration variable from v to f by using the relation In terms of f , Eq. (2.21) becomes By identifying this equation with Eq. (2.18), we can express P φ (f ) in terms of the Boltzmann distribution function as (2.24) Fig. 1 shows a plot of P φ (f ) for µ = 100Hz, v * = 220km/s.

Scalar field in the Earth's rest frame
Since observations by the terrestrial detectors are conducted on the Earth, we need to translate the scalar field represented in the galaxy's center-of-mass frame into the one in the Earth's rest frame.Although the Earth, which is orbiting in the Milky Way Galaxy, is accelerating with respect to the center-of-mass frame of the Galaxy, we ignore the effect of the acceleration and make an approximation that the Earth is moving with a constant velocity v 0 in the Galaxy.We also ignore the Earth's orbital motion around the Sun since the orbital velocity ( 30 km/s) is much smaller than |v 0 | 200 km/s.To be definite, we put prime on coordinates and relevant quantities when they are defined in the Galaxy's rest frame.For instance, spacetime coordinates are written as (t , x ), and frequency and the propagation direction are written as (f , Ω ).On the other hand, quantities without prime are defined in the Earth's rest frame.Based on our assumptions, both frames are inertial systems and quantities in different frames are related by the Lorentz transformation.Frequency f and the wavenumber vector k transform as Since motions are non-relativistic (i.e., |v 0 | v * = O(10 −3 ) 1), we retain only the relevant lowest order terms in v 0 or k/m in the above relations.After some manipulations, we obtain ) Since φ itself is Lorentz-invariant, Eq. (2.14), which is expressed in terms of the variables of the Galaxy's rest frame, may be directly used to represent φ in the Earth's rest frame.We note that a combination kdf d Ω, where k is the wavenumber, is Lorentzinvariant: (2.30) By exploiting this relation, Eq. (2.14) may be written as Note that (t, x) are coordinates in the Earth's rest frame.Plugging the transformation rules for f , k , and Ω , we have For notational convenience, we introduce a new function ū(f, Ω) by so that φ may be written as This is the expression of φ that is used in the rest of this paper.

Detector's signal of the scalar field
A mirror obeys the equation of motion given by Eq. (2.12) and undergoes oscillations caused by the scalar force ∇ϕ.Response of a mirror by the action of the scalar field was computed in [12].For the scalar field given by Eq. (2.34), mirror's motion is written as Then, using the result in [12], the signal s(t) of the Michelson-type interferometer located at x, which is the phase difference of the lasers propagating along each arm, can be computed.The result is given by where the detector pattern function is given by (2.37) The first/(second) term corresponds to δt 1 /(δt 2 ) defined in [12], respectively.This function satisfies F * ( Ω, f ) = F (− Ω, f ).
The Fourier transform of the signal is given by

Two-detectors correlation
In the observation, the detector's output S(t) consists of the signal s (if it exists) and the noise n as S(t) = s(t) + n(t). (2.39) Since ϕ is a stochastic variable, the signal s(t) also behaves in a stochastic manner.Practically, it is difficult to distinguish the signal from the noise.As it is usually done, we thus consider correlation of the outputs between two detectors in order to extract the signal from the noise.Using Eqs.(2.38), (2.16), and (2.33), we can compute the cross correlation of the outputs between the detector 1 and the detector 2 as (2.40) Here, it has been assumed that the correlation of the noise between the two detectors is zero.∆x = x 2 − x 1 is the separation between the two detectors.Defining the overlap function Γ φ (f ) in the absence of the Earth's motion as the cross correlation becomes Here ε(f ) defined by quantifies the change of the correlation signal caused by the motion of the solar system with respect to the rest frame of the Milky Way Galaxy.Notice that ε(f ) is a complex number.
From the relation F * ( Ω, f ) = F (− Ω, f ), we can verify that Γ φ (f ) is a real function.For LIGO interferometers, we have f L 1 and 2πf η| Ω • ∆x| 1.Thus, it is a good approximation to ignore the phase proportional to ∆x in Eq. (2.41).With this approximation, we can perform the angular integration analytically; (2.44) where a 1 , a 2 are defined by (2.46) Fig. 2 shows numerically computed Γ φ (f ) and the approximate one (2.44) for κα SiO 2 = 1 #1 .Clearly, the approximate formula nicely reproduces the exact one.
Our formulation expanded up to this point assumes infinitely long observation time T .However, in reality, the processed public data which we can use, is obtained for T = 192 s.Thus, we need to relate the real observable with what we have computed for T → ∞.First, suppose that the data S(t) was taken during the period (−T /2, T /2).We define the Fourier transformation of the data S(t) as (2.47) #1 Numerical values of the vectors m, n for the existing detectors can be found at https://lscsoft.docs.ligo.org/lalsuite/lal/_l_a_l_detectors_8h_source.htmlwhere W T (t) is the Hann Window function which is used in LVK Collaboration.On the side, we have put the subscript "T" to emphasize that it is directly computed from data.Meanwhile, the original Fourier transformation is defined by (2.48) From these equations, we find that ST (f ) and S(f ) are related as where WT (f ) is defined by and in the last equation we have given an explicit expression for the Hann Window function.Using this relation, we can derive the cross-correlation of ST (f ) in terms of P φ (f ) and Γ φ (f ) as To evaluate this integral approximately, we note P φ (f ) is sharply peaked at µ ≤ f µ + µv 2 * .Thus, the integration is contributed only from this short interval.Since µv 2 * 1/T , the Window function remains almost constant in this interval.Then, it is a good approximation to replace f with µ in WT (f − f ) and pull it out of the integration; Here Q(µ) defined by quantifies the effect caused by the motion of the solar system.Because ε(f ) is a complex number, Q(µ) is also a complex number.Fig. 3 shows a plot of Re Q(µ) and Im Q(µ).In making this plot, it is assumed that dark matter wind comes from the direction specified by ( , b) = (270 • , 0) in the galactic coordinates which correspond to (α, δ) = (3.70231,−0.81267) in the equilateral coordinates and the signal is averaged over α to take into account the Earth's daily rotation.We find that the signal is suppressed at µ 70 Hz (m 3 × 10 −13 eV) whereas it is enhanced by a factor of ∼ 2 at other frequencies.We expect that the imaginary part of Q may be used as additional observable to increase the sensitivity of the search of the scalar dark matter by using the laser interferometers.Since only the real part of the cross correlation is publically available, we do not make use of Im Q in our analysis.
It is possible to perform the integral ∞ µ Γ φ P φ df analytically to a good approximation.To this end, using Eqs.(2.24) and (2.44), we obtain (2.54) where we have replaced f which is not in η with µ.By changing variable from f to v by the relation (2.22), we have (2.55) Thus, our final expression of the cross correlation becomes The cross-correlation statistic Ĉ(f ) given by the LVK collaboration is defined as a summation over discretized frequency bins with its width 1/T ; where f i = f − ∆f /2 + i/T , γ T is the overlap reduction function for gravitational waves, and S 0 ≡ 3H 2 0 /(10π 2 f 3 ).In the presence of the scalar field dark matter, using Eq.(2.56), the expectation value of Ĉ(f ) becomes This is the actual signal that we seek in the GW data provided by the LVK Collaboration.

Results
We search for the signal given by Eq. (2.58) in the data of the two LIGO detectors during their third observing run (O3).The LVK Collaboration gives the cross correlation statistics which has been obtained by computing Ĉ for each segment containing data of the duration T and averaging them #2 .Thus, the cross correlation statistics obeys the Gaussian distribution, and the likelihood becomes [25] p( Ĉ|Θ) where Θ = {µ, α SiO 2 } is a set of parameters characterizing the scalar field model.
In our analysis, we employ the Bayesian inference.We first fix µ and compute the posterior on α SiO 2 by adopting two priors (uniform in α SiO 2 and uniform in ln α SiO 2 ).We then repeat this calculation by scanning the relevant range of µ.We did not find any conclusive evidence of a signal in data.As a result, we are able to place upper limit on α SiO 2 for various values of µ.
Fig. 4 is one of the main results.Each black dot in both panels shows upper limit for corresponding value of m.The left/(right) panel assumes uniform prior for α SiO 2 /(ln α SiO 2 ).#2 Data is available at https://dcc.ligo.org/cgi-bin/DocDB/ShowDocument?.submit=Identifier& docid=G2001287&version=5 Although the original data of the cross correlation statistics is given with frequency interval δf 0.03 Hz, we did not use data for all frequencies but picked up data only at frequency interval 10δf in order to reduce the file size of Fig. 4. Vertical white bands in both panels where there are no black dots are due to deficit of data.Spike at m 3 × 10 −13 eV reflects the effect of the motion of the Earth with respect to the rest frame of the Milky Way galaxy (see discussion regarding Fig. 3).Finite spread of the black dots (i.e., that difference of values of y-axis between neighboring dots fluctuate) is caused by the fluctuations of the cross correlation data that exist even among slightly different frequencies.We find that there is no significant difference between the left panel and the right panel.Thus, we use only the left panel for the following results.Fig. 5 shows translation of Fig. 4 into the upper limits on the four coupling constants (d g , d m, d me , d e ) for five representative cases by using Eq.(2.13).In addition to our constraint by LIGO O3 data, we also show constraints from another gravitational-wave (GW) experiment (GEO600) [24], which measures change of the size of the beam splitter and is sensitive only to d me and d e , as well as other non-GW experiments including test of violation of the weak equivalence principle (WEP) [21,26,27] and the fifth-force experiment [22].Experiments of WEP such as MICROSCOPE mission and the torsion balance are only sensitive to the difference of gravitational accelerations between different materials.The effect of the violation of WEP is normally expressed by the Eotvos parameter η defined by [27] Here a A /(a B ) is the gravitational acceleration of material A/(B) measured at distance r from the center of the Earth and α E is α A for materials constituting the Earth.R E is the mean Earth radius and Φ(x) ≡ 3(x cosh x − sinh x)/x 3 .In our analysis, we assume that the Earth is made of silica (50%) and magnesium oxide (50%).For the MICRO-SCOPE mission, A and B are beryllium (Be) and Platinum (Pt).For the torsion balance experiment, A and B are beryllium (Be) and Titanium (Ti).The fifth force experiments measure extra force operating between two materials in addition to the gravitational force.For a pair of material A and the Earth, the experiment can probe a direct product α A α E .
We make an approximation that the material used in [22] is purely aluminium (Al).The quantity α Al α E does not vanish for all the five cases in Fig. 5. Thus, the upper limit set by the fifth-force experiment is present in all the panels and is shown as a blue curve.
Let us now consider the upper limit for each case one by one.The fist case in which only d g is assumed to be non-vanishing is shown in the top left panel.Since both d me and d e are zero, GEO600 does not give any meaningful constraint in this case.Since the component-dependent part of Eq. (2.9) contains d g , WEP experiments can probe d g .The upper limit set by WEP is shown as a green curve.We find that WEP provides the most stringent constraint by several orders of magnitude stronger than the other constraints.Even the LIGO at design sensitivity will not be able to reach the current constraint by WEP.However, Cosmic Explorer, planned third GW detector, may be able to have better sensitivity than WEP [12].The second case in which only d m is assumed to be nonvanishing is shown in the top right panel.This case is similar to the first case except that all the constraints are weaker than the ones in the first case.This is because α E is less sensitive to d m than d g , namely α E 0.9d g /(0.08d m) for the first/(second) case.The third case in which only d me is assumed to be non-vanishing is shown in the middle left panel.In this case, GEO600 gives stronger constraint than LIGO for m 3 × 10 −13 eV.The WEP still provides the tightest upper limit on d me .The forth case in which only d e is assumed to be non-vanishing is shown in the middle right panel.This case is similar to the third case except that WEP is much stronger.This can be understood from Eq. (2.9) that component-dependent term proportional to d me is suppressed by a factor (A − 2Z)/A compared to d e .The fifth case in which d g = d m = d me and d e = 0 is assumed is shown in the bottom panel.In this case, WEP constraint is not present.The LIGO (O3) constraint is weaker than the fifth force experiment although the difference is less than a factor of ∼ 5 at large mass side.The upgraded LIGO observation is expected to exceed the fifth force constraint [12].

Conclusion
If dark matter is a light scalar field weakly interacting with elementary particles, such a field induces oscillations of the physical constants, which results in time-varying force acting on macroscopic objects.We searched for the signal of the scalar field in the latest data of LIGO observations (O3 run) and placed upper limits on the coupling strength between the mirrors in the laser interferometers and the scalar field.To this end, we first formulated the cross-correlation statistics that can be readily compared with publically available data.It has been found that inclusion of the anisotropies of the velocity distribution of dark matter enhances the signal by a factor of ∼ 2 except for the narrow mass range around 3 × 10 −13 eV for which the correlation between the interferometer at Livingston and the one at Hanford is blind to the signal.
From the non-detection of the signal, we also derived the upper limits on the coupling constants between the elementary particles and the scalar field for five representative cases.For all the cases where the weak equivalence principle is not satisfied, tests of the violation of the weak equivalence principle provides the tightest upper limit on the coupling constants.Upper limits from the fifth-force experiment are always stronger than the ones from LIGO, but the difference is less than a factor of ∼ 5 at large-mass range, which vividly demonstrates that the gravitational-wave experiments have been coming close to the sensitivity of the fifth force experiment.For the models where only the coupling to electrons or photos is non-vanishing, upper limits from GEO600 are stronger than those from LIGO at large-mass range while LIGO provides stronger limits for the case where the weak equivalence principle is not violated.
Our study presents that gravitational-wave experiments are starting to bring us meaningful information to constrain the nature of dark matter.The formulation provided in this paper may be applied to the data of upcoming experiments as well and is expected to probe much wider parameter range of the model.

m 2 10 − 6 , 3
2π within the width ∆f /f = v 2 * where v * = 220 km/s is the velocity dispersion of dark matter.Thus, oscillations of φ produces almost monochromatic time variation of the physical constants listed in Eqs.(2.3)-(2.6).The frequency range corresponding to the mass range 10 −13 eV m 10 −11 eV is 10 Hz m 2π 10 Hz which falls in the frequency band of the ground-based laser interferometers.

1 .Figure 5 :
Figure 5: Upper limits on the coupling constants for five representative cases which have been obtained by translating Fig. 4 using Eq.(2.13).Top left panel: upper limit on d g assuming d m = d me = d e = 0. Top right panel: upper limit on d m assuming d g = d me = d e = 0. Middle left panel: upper limit on d me assuming d g = d m = d e = 0. Middle right panel: upper limit on d e assuming d g = d m = d me = 0. Bottom panel: upper limit on d g assuming d g = d m = d me , d e = 0 which is the model of the gravitational non-minimal coupling discussed in the Appendix. .21)