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

Kouki Fukusumi$\,{}^{1}$, Soichiro Morisaki$\,{}^{2}$, Teruaki
Suyama$\,{}^{1}$

1Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama,
Meguro-ku, Tokyo 152-8551, Japan
2Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan

###### Abstract

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.

##  1 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}~{}{\rm eV}\lesssim m\lesssim 10^{60}~{}{\rm 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\ll 10~{}{\rm 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 [5, 3, 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}~{}{\rm eV},10^{-11}~{}{\rm 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 $\phi$ varies in space, which means that the scalar
force is parallel to ${\bm{\nabla}}\phi$. 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 $\delta 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 $\delta
t_{2}$ in [12]). For aLIGO, the former/(latter) becomes dominant if the
frequency of the scalar field oscillation is larger/(smaller) than $20~{}{\rm
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 ($\simeq$ 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.

##  2 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.

###  2.1 Model considered in this paper

The mathematical framework of the model we consider is given in [18]. In this
model, the scalar field $\phi$ which comprises all dark matter has
interactions with elementary particles. The Lagrangian density for the scalar
field is given by

$\displaystyle{\cal L}$
$\displaystyle=-\frac{1}{2}{(\partial\phi)}^{2}-\frac{m^{2}}{2}\phi^{2}+{\cal
L}_{\rm int},$ (2.1) $\displaystyle{\cal L}_{\rm
int}=-\kappa\phi\bigg{[}\frac{d_{e}}{4e^{2}}F_{\mu\nu}F^{\mu\nu}-\frac{d_{g}\beta_{3}}{2g_{3}}G_{\mu\nu}^{A}G^{A\mu\nu}-\sum_{i=e,u,d}(d_{m_{i}}+\gamma_{m_{i}}d_{g})m_{i}{\bar{\psi}_{i}}\psi_{i}\bigg{]}.$
(2.2)

Here $\kappa\equiv\frac{\sqrt{4\pi}}{M_{\rm pl}}$, and $M_{\rm
pl}=1/\sqrt{G}\simeq 1.22\times 10^{19}~{}{\rm GeV}$ is the Planck mass. The
first, second, and third term in ${\cal L}_{\rm int}$ represents interaction
between $\phi$ 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. $\beta_{3}$ is the QCD beta function and $\gamma_{m_{i}}$ is the
anomalous dimension. The interaction Lagrangian given above is defined at some
high energy scale. We can translate the effect of ${\cal L}_{\rm int}$ into
equivalent but more physically relevant representation by considering the
renormalization-group flow:

$\displaystyle\alpha(\varphi)=(1+d_{e}\varphi)\alpha,$ (2.3)
$\displaystyle\Lambda_{3}(\varphi)=(1+d_{g}\varphi)\Lambda_{3},$ (2.4)
$\displaystyle m_{e}(\varphi)=(1+d_{m_{e}}\varphi)m_{e},$ (2.5) $\displaystyle
m_{u,d}(\varphi)=(1+d_{m_{u},m_{d}}\varphi)m_{u,d}(\Lambda_{3}).$ (2.6)

Here, $\varphi\equiv\kappa\phi$ is the scalar field normalized by the Planck
scale. $\alpha$ is the fine-structure constant, $\Lambda_{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 $\varphi$, these physical constants deviate from those in
the absence of $\varphi$. We assume that the mass $m$ is in the range
$10^{-13}~{}{\rm eV}\lesssim m\lesssim 10^{-11}~{}{\rm eV}$. In this case, the
$\phi$ 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\simeq\frac{m}{2\pi}$ within the width
$\Delta f/f=\frac{v_{*}^{2}}{2}\simeq 10^{-6}$, where $v_{*}=220~{}{\rm km/s}$
is the velocity dispersion of dark matter. Thus, oscillations of $\phi$
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}~{}{\rm eV}\lesssim m\lesssim 10^{-11}~{}{\rm eV}$ is $10~{}{\rm
Hz}\lesssim\frac{m}{2\pi}\lesssim 10^{3}~{}{\rm Hz}$ which falls in the
frequency band of the ground-based laser interferometers.

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
$\phi$ couples to the kinetic term of the fermion as

${\cal
L}=-e^{f(\phi)}{\bar{\psi}}\gamma^{\mu}\partial_{\mu}\psi-\frac{1}{2}e^{f(\phi)}\partial_{\mu}f{\bar{\psi}}\gamma^{\mu}\psi,$
(2.7)

where we leave $f(\phi)$ 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 $\psi\to e^{-\frac{1}{2}f}\psi$, the kinetic term in
terms of the new field becomes

${\cal L}=-{\bar{\psi}}\gamma^{\mu}\partial_{\mu}\psi,$ (2.8)

and the coupling $f(\phi)$ has been completely absorbed into the field
redefinition.

Since mass of an atom depends on physical constants, it depends on $\varphi$
through Eqs. (2.3)-(2.6). Denoting the mass of an atom by $m_{A}(\varphi)$,
$m_{A}$ in the small $\varphi$ limit is approximately written as [18]

$\displaystyle\alpha_{A}$ $\displaystyle\equiv\frac{d\ln
m_{A}(\varphi)}{d\varphi}$ $\displaystyle\simeq
d_{g}^{*}+\bigg{[}-\frac{0.036}{A^{\frac{1}{3}}}-0.02\frac{{(A-2Z)}^{2}}{A^{2}}-1.4\times
10^{-4}\frac{Z(Z-1)}{A^{\frac{4}{3}}}\bigg{]}(d_{\hat{m}}-d_{g})$
$\displaystyle~{}~{}~{}-2.7\times
10^{-4}\frac{A-2Z}{A}(d_{m_{e}}-d_{g})+\left(-4.1\frac{A-2Z}{A}+7.7\frac{Z(Z-1)}{A^{\frac{4}{3}}}\right)\times
10^{-4}d_{e},$ (2.9)

and $d_{g}^{*}$ is defined by

$d_{g}^{*}=d_{g}+0.093(d_{\hat{m}}-d_{g})+2.7\times
10^{-4}(d_{m_{e}}-d_{g})+2.7\times 10^{-4}d_{e},$ (2.10)

where $A,Z$ is the mass number and the atomic number, respectively, and

$d_{\hat{m}}\equiv\frac{d_{m_{u}}m_{u}+d_{m_{d}}m_{d}}{m_{u}+m_{d}}.$ (2.11)

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_{\hat{m}}=d_{m_{e}}=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 $\phi$
causes the time variation of the mass $M(\varphi)$ of a macroscopic object.
The equation of motion of the macroscopic object under the influence of
$\varphi$ is given by (e.g., [12])

$\frac{d^{2}{\bm{x}}}{dt^{2}}\simeq-\langle\alpha_{A}\rangle{\bm{\nabla}}\varphi,$
(2.12)

where $\langle\alpha_{A}\rangle$ is $\alpha_{A}$ averaged over different types
of molecules constituting the macroscopic object. The mirrors installed in the
interferometers are purely made of silica (=${\rm SiO_{2}}$). In this case, we
have

$\alpha_{SiO_{2}}\simeq d_{g}+0.083(d_{\hat{m}}-d_{g})+2.7\times
10^{-4}(d_{m_{e}}-d_{g})+3.1\times 10^{-3}d_{e}.$ (2.13)

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.

###  2.2 Scalar field in the Milky Way Galaxy

####  2.2.1 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_{*}\simeq 220~{}{\rm km/s}$. In this
coordinate system, the scalar dark matter may be expressed as the
superposition of plane waves as

$\phi(t,{\bm{x}})=\int
d{\hat{\bm{\Omega}}}\int_{-\infty}^{\infty}df~{}e^{-2\pi
if(t-\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}})}u(f,{\hat{\bm{\Omega}}}).$ (2.14)

Here, $\eta$ is defined by

$\eta=\sqrt{1-\frac{\mu^{2}}{f^{2}}},~{}~{}~{}~{}~{}~{}~{}~{}~{}\left(\mu=\frac{m}{2\pi}\right),$
(2.15)

and ${\hat{\bm{\Omega}}}$ is the unit vector representing the propagation
direction of the plane wave $u(f,{\hat{\bm{\Omega}}})$. From the reality of
$\phi$, we have $u^{*}(f,{\hat{\bm{\Omega}}})=u(-f,{\hat{\bm{\Omega}}})$. From
the stationarity and isotropy condition, we assume that the ensemble average
of the stochastic variable $u(f,{\hat{\bm{\Omega}}})$ can be written as

$\langle
u^{*}(f,{\hat{\bm{\Omega}}})u(f^{\prime},{\hat{\bm{\Omega}^{\prime}}})\rangle=\delta(f-f^{\prime})\frac{1}{4\pi}\delta({\hat{\bm{\Omega}}},{\hat{\bm{\Omega}^{\prime}}})\times\frac{1}{2}P_{\phi}(f).$
(2.16)

Here, $P_{\phi}(f)$ is the powerspectrum of $\phi$ 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

$\rho_{\phi}=\frac{1}{2}{\dot{\phi}}^{2}+\frac{1}{2}{({\bm{\nabla}}\phi)}^{2}+\frac{m^{2}}{2}\phi^{2}.$
(2.17)

Using Eq. (2.14), the ensemble average of $\rho_{\phi}$ becomes

$\langle\rho_{\phi}\rangle=\int_{\mu}^{\infty}df~{}{(2\pi f)}^{2}P_{\phi}(f).$
(2.18)

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({\bm{v}})$ in the velocity space. As it is usually
done, we assume that the distribution of the particles obeys the Boltzmann
distribution as

$g({\bm{v}})=\frac{3^{\frac{3}{2}}\rho_{\rm
loc}}{{(2\pi)}^{3/2}mv_{*}^{3}}\exp\left(-\frac{3v^{2}}{2v_{*}^{2}}\right).$
(2.19)

Here $v=|{\bm{v}}|$, and $v_{*}$ is the velocity dispersion of dark matter.
The normalization constant is fixed by requiring that the integral of
$g({\bm{v}})$ coincides with the local number density of dark matter
particles;

$\int d{\bm{v}}~{}g({\bm{v}})=\frac{\rho_{\rm loc}}{m}.$ (2.20)

Therefore, we have

$\int d{\bm{v}}\frac{3^{\frac{3}{2}}\rho_{\rm
loc}}{{(2\pi)}^{3/2}v_{*}^{3}}\exp\left(-\frac{3v^{2}}{2v_{*}^{2}}\right)=\rho_{\rm
loc}.$ (2.21)

In order to compare, we change the integration variable from $v$ to $f$ by
using the relation

$f^{2}=\mu^{2}(1+v^{2}).$ (2.22)

In terms of $f$, Eq. (2.21) becomes

$\int_{\mu}^{\infty}df~{}{\left(\frac{3}{2}\right)}^{\frac{3}{2}}\frac{4\rho_{\rm
loc}}{\sqrt{\pi}\mu^{3}v_{*}^{3}}f^{2}\eta\exp\left(-\frac{3f^{2}\eta^{2}}{2v_{*}^{2}\mu^{2}}\right)=\rho_{\rm
loc}$ (2.23)

By identifying this equation with Eq. (2.18), we can express $P_{\phi}(f)$ in
terms of the Boltzmann distribution function as

$P_{\phi}(f)={\left(\frac{3}{2}\right)}^{\frac{3}{2}}\frac{4\rho_{\rm
loc}}{\pi^{5/2}\mu^{3}v_{*}^{3}}\eta\exp\left(-\frac{3f^{2}\eta^{2}}{2v_{*}^{2}\mu^{2}}\right).$
(2.24)

Fig. 1 shows a plot of $P_{\phi}(f)$ for $\mu=100{\rm Hz},v_{*}=220{\rm
km/s}$.

Figure 1: Plot of $P_{\phi}(f)$ for $\mu=100~{}{\rm Hz}$, $v_{*}=220{\rm
km/s}$.

####  2.2.2 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
${\bm{v}_{0}}$ in the Galaxy. We also ignore the Earth’s orbital motion around
the Sun since the orbital velocity ($\simeq 30~{}{\rm km/s}$) is much smaller
than $|{\bm{v}_{0}}|\simeq 200~{}{\rm 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^{\prime},{\bm{x}^{\prime}})$, and frequency and the propagation direction
are written as $(f^{\prime},{\hat{\Omega}^{\prime}})$. 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 ${\bm{k}}$ transform as

$\displaystyle
f^{\prime}=\gamma\left(f+\frac{{\bm{v}_{0}}\cdot{\bm{k}}}{2\pi}\right),~{}~{}~{}~{}~{}\gamma\equiv\frac{1}{\sqrt{1-{\bm{v}_{0}}^{2}}}$
(2.25) $\displaystyle{\bm{k}^{\prime}}=\gamma\left({\bm{k}}+2\pi
f{\bm{v}_{0}}\right).$ (2.26)

Since motions are non-relativistic (i.e., $|{\bm{v}_{0}}|\simeq v_{*}={\cal
O}(10^{-3})\ll 1$), we retain only the relevant lowest order terms in
${\bm{v}_{0}}$ or ${\bm{k}}/m$ in the above relations. After some
manipulations, we obtain

$\displaystyle f^{\prime}$
$\displaystyle=f+\frac{\mu}{2}{\bm{v}_{0}}^{2}+\mu\eta{\bm{v}_{0}}\cdot{\hat{\bm{\Omega}}},$
(2.27) $\displaystyle k^{\prime}$
$\displaystyle=\Big{|}{\hat{\bm{\Omega}}}+\frac{\bm{v}_{0}}{\eta}\Big{|}k,$
(2.28) $\displaystyle{\hat{\bm{\Omega}^{\prime}}}$
$\displaystyle={\Big{|}{\hat{\bm{\Omega}}}+\frac{\bm{v}_{0}}{\eta}\Big{|}}^{-1}\left({\hat{\bm{\Omega}}}+\frac{\bm{v}_{0}}{\eta}\right).$
(2.29)

Since $\phi$ 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 $\phi$ in the Earth’s rest frame. We note that a combination
$kdfd{\hat{\Omega}}$, where $k$ is the wavenumber, is Lorentz-invariant:

$k^{\prime}df^{\prime}d{\hat{\bm{\Omega}^{\prime}}}=kdfd{\hat{\bm{\Omega}}}.$
(2.30)

By exploiting this relation, Eq. (2.14) may be written as

$\phi(t,{\bm{x}})=\int d{\hat{\bm{\Omega}}}df~{}\frac{k}{k^{\prime}}e^{-2\pi
if(t-\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}})}u(f^{\prime},{\hat{\bm{\Omega}^{\prime}}}).$
(2.31)

Note that $(t,{\bm{x}})$ are coordinates in the Earth’s rest frame. Plugging
the transformation rules for $f^{\prime},k^{\prime},$ and
${\bm{\Omega}^{\prime}}$, we have

$\phi(t,{\bm{x}})=\int d{\hat{\bm{\Omega}}}df~{}e^{-2\pi
if(t-\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}})}{\Big{|}{\hat{\bm{\Omega}}}+\frac{\bm{v}_{0}}{\eta}\Big{|}}^{-1}u\left(f+\frac{\mu}{2}{\bm{v}_{0}}^{2}+\mu\eta{\bm{v}_{0}}\cdot{\hat{\bm{\Omega}}},\frac{\eta{\hat{\bm{\Omega}}}+{\bm{v}_{0}}}{|\eta{\hat{\bm{\Omega}}}+{\bm{v}_{0}}|}\right).$
(2.32)

For notational convenience, we introduce a new function
${\bar{u}}(f,{\hat{\bm{\Omega}}})$ by

${\bar{u}}(f,{\hat{\bm{\Omega}}})={\Big{|}{\hat{\bm{\Omega}}}+\frac{\bm{v}_{0}}{\eta}\Big{|}}^{-1}u\left(f+\frac{\mu}{2}{\bm{v}_{0}}^{2}+\mu\eta{\bm{v}_{0}}\cdot{\hat{\bm{\Omega}}},\frac{\eta{\hat{\bm{\Omega}}}+{\bm{v}_{0}}}{|\eta{\hat{\bm{\Omega}}}+{\bm{v}_{0}}|}\right),$
(2.33)

so that $\phi$ may be written as

$\phi(t,{\bm{x}})=\int d{\hat{\bm{\Omega}}}df~{}e^{-2\pi
if(t-\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}})}{\bar{u}}(f,{\hat{\bm{\Omega}}}).$
(2.34)

This is the expression of $\phi$ that is used in the rest of this paper.

###  2.3 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 ${\bm{\nabla}}\varphi$. 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

$\delta{\bm{x}}(t)=-\kappa\alpha_{SiO_{2}}\int
d{\hat{\bm{\Omega}}}\int_{-\infty}^{\infty}\frac{df}{2\pi if}\eta e^{-2\pi
if(t-\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}})}{\bar{u}}(f,{\hat{\bm{\Omega}}}){\hat{\bm{\Omega}}}.$
(2.35)

Then, using the result in [12], the signal $s(t)$ of the Michelson-type
interferometer located at ${\bm{x}}$, which is the phase difference of the
lasers propagating along each arm, can be computed. The result is given by

$s(t)=\int
d{\hat{\bm{\Omega}}}\int_{-\infty}^{\infty}df~{}F({\hat{\bm{\Omega}}},f)e^{-2\pi
if(t-\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}})}{\bar{u}}(f,{\hat{\bm{\Omega}}}),$
(2.36)

where the detector pattern function is given by

$F({\hat{\bm{\Omega}}},f)=-\frac{\kappa\alpha_{SiO_{2}}}{\pi
iLf}\left({\hat{\bm{m}}}\cdot{\hat{\bm{\Omega}}}-{\hat{\bm{n}}}\cdot{\hat{\bm{\Omega}}}\right)\eta\sin^{2}(\pi
fL)-\kappa\alpha_{SiO_{2}}\left({({\hat{\bm{m}}}\cdot{\hat{\bm{\Omega}}})}^{2}-{({\hat{\bm{n}}}\cdot{\hat{\bm{\Omega}}})}^{2}\right)\eta^{2}.$
(2.37)

The first/(second) term corresponds to $\delta t_{1}$/$(\delta t_{2})$ defined
in [12], respectively. This function satisfies
$F^{*}({\hat{\bm{\Omega}}},f)=F(-{\hat{\bm{\Omega}}},f)$.

The Fourier transform of the signal is given by

${\tilde{s}}(f)=\int d{\hat{\bm{\Omega}}}~{}F({\hat{\bm{\Omega}}},f)e^{2\pi
if\eta{\hat{\bm{\Omega}}}\cdot{\bm{x}}}{\bar{u}}(f,{\hat{\bm{\Omega}}})$
(2.38)

####  2.3.1 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 $\varphi$ 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

$\langle{\tilde{S}_{1}^{*}}(f){\tilde{S}_{2}}(f^{\prime})\rangle=\int\frac{d{\hat{\bm{\Omega}}}}{4\pi}F_{1}^{*}({\hat{\bm{\Omega}}},f)F_{2}({\hat{\bm{\Omega}}},f)e^{2\pi
if\eta{\hat{\bm{\Omega}}}\cdot{\bm{\Delta}x}}\times\frac{1}{2}\delta(f-f^{\prime})P_{\phi}(f+\frac{\mu}{2}{\bm{v}_{0}}^{2}+\mu\eta{\bm{v}_{0}}\cdot{\hat{\bm{\Omega}}}).$
(2.40)

Here, it has been assumed that the correlation of the noise between the two
detectors is zero. ${\bm{\Delta}x}={\bm{x}_{2}}-{\bm{x}_{1}}$ is the
separation between the two detectors. Defining the overlap function
$\Gamma_{\phi}(f)$ in the absence of the Earth’s motion as

$\Gamma_{\phi}(f)\equiv\int\frac{d{\hat{\bm{\Omega}}}}{4\pi}F_{1}^{*}({\hat{\bm{\Omega}}},f)F_{2}({\hat{\bm{\Omega}}},f)e^{2\pi
if\eta{\hat{\bm{\Omega}}}\cdot{\bm{\Delta}x}},$ (2.41)

the cross correlation becomes

$\langle{\tilde{S}_{1}^{*}}(f){\tilde{S}_{2}}(f^{\prime})\rangle=\frac{\delta(f-f^{\prime})}{2}\Gamma_{\phi}(f)P_{\phi}(f)\varepsilon(f).$
(2.42)

Here $\varepsilon(f)$ defined by

$\varepsilon(f)\equiv\frac{1}{\Gamma_{\phi}(f)P_{\phi}(f)}\int\frac{d{\hat{\bm{\Omega}}}}{4\pi}F_{1}^{*}({\hat{\bm{\Omega}}},f)F_{2}({\hat{\bm{\Omega}}},f)e^{2\pi
if\eta{\hat{\bm{\Omega}}}\cdot{\bm{\Delta}x}}P_{\phi}(f+\frac{\mu}{2}{\bm{v}_{0}}^{2}+\mu\eta{\bm{v}_{0}}\cdot{\hat{\bm{\Omega}}})$
(2.43)

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 $\varepsilon(f)$ is a complex number.

From the relation $F^{*}({\hat{\bm{\Omega}}},f)=F(-{\hat{\bm{\Omega}}},f)$, we
can verify that $\Gamma_{\phi}(f)$ is a real function. For LIGO
interferometers, we have $fL\ll 1$ and $2\pi
f\eta|{\hat{\bm{\Omega}}}\cdot{\bm{\Delta}x}|\ll 1$. Thus, it is a good
approximation to ignore the phase proportional to ${\bm{\Delta}x}$ in Eq.
(2.41). With this approximation, we can perform the angular integration
analytically;

$\Gamma_{\phi}(f)\approx{(\kappa\alpha_{SiO_{2}}\eta)}^{2}\left(a_{1}{(fL)}^{2}+a_{2}\eta^{2}\right).$
(2.44)

where $a_{1},a_{2}$ are defined by

$\displaystyle a_{1}$
$\displaystyle=\frac{\pi^{2}}{3}({\hat{m}_{1}}-{\hat{n}_{1}})\cdot({\hat{m}_{2}}-{\hat{n}_{2}}),$
(2.45) $\displaystyle a_{2}$
$\displaystyle=\frac{2}{15}\big{[}\left(({\hat{m}_{1}}+{\hat{n}_{1}})\cdot{\hat{m}_{2}}\right)\left(({\hat{m}_{1}}-{\hat{n}_{1}})\cdot{\hat{m}_{2}}\right)-({\hat{m}_{2}}\leftrightarrow{\hat{n}_{2}})\big{]}.$
(2.46)

Fig. 2 shows numerically computed $\Gamma_{\phi}(f)$ and the approximate one
(2.44) for $\kappa\alpha_{SiO_{2}}=1$ #1#1#1 Numerical values of the vectors
${\hat{m}},{\hat{n}}$ for the existing detectors can be found at
https://lscsoft.docs.ligo.org/lalsuite/lal/_l_a_l_detectors_8h_source.html.
Clearly, the approximate formula nicely reproduces the exact one.

Figure 2: Plot of $\Gamma_{\phi}(f)$ for $\mu=100~{}{\rm Hz}$,
$\kappa\alpha_{SiO_{2}}=1$.

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~{}{\rm s}$. Thus, we need to relate the real observable
with what we have computed for $T\to\infty$. 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

${\tilde{S}}_{\rm T}(f)=\int_{-\frac{T}{2}}^{\frac{T}{2}}~{}dt~{}e^{2\pi
ift}S(t)W_{T}(t),$ (2.47)

where $W_{T}(t)$ is the Hann Window function which is used in LVK
Collaboration. On the left-hand side, we have put the subscript “${\rm T}$” to
emphasize that it is directly computed from data. Meanwhile, the original
Fourier transformation is defined by

${\tilde{S}}(f)=\int_{-\infty}^{\infty}~{}dt~{}e^{2\pi ift}S(t).$ (2.48)

From these equations, we find that ${\tilde{S}}_{\rm T}(f)$ and
${\tilde{S}}(f)$ are related as

${\tilde{S}}_{\rm
T}(f)=\int_{-\infty}^{\infty}df^{\prime}~{}{\tilde{S}}(f^{\prime}){\tilde{W}}_{T}(f-f^{\prime}),$
(2.49)

where ${\tilde{W}}_{T}(f)$ is defined by

${\tilde{W}}_{T}(f)=\int_{-\frac{T}{2}}^{\frac{T}{2}}~{}dt~{}e^{2\pi
ift}W_{T}(t)=\sqrt{\frac{8}{3}}\frac{\sin(\pi fT)}{2\pi f(1-f^{2}T^{2})},$
(2.50)

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
${\tilde{S}}_{\rm T}(f)$ in terms of $P_{\phi}(f)$ and $\Gamma_{\phi}(f)$ as

$\langle{\tilde{S}}_{1,\rm T}^{*}(f){\tilde{S}}_{2,\rm
T}(f)\rangle=\frac{1}{2}\int_{-\infty}^{\infty}{|{\tilde{W}}_{T}(f-f^{\prime})|}^{2}\Gamma_{\phi}(f^{\prime})P_{\phi}(f^{\prime})\varepsilon(f^{\prime})df^{\prime}.$
(2.51)

To evaluate this integral approximately, we note $P_{\phi}(f^{\prime})$ is
sharply peaked at $\mu\leq f^{\prime}\lesssim\mu+\mu v_{*}^{2}$. Thus, the
integration is contributed only from this short interval. Since $\mu
v_{*}^{2}\ll 1/T$, the Window function remains almost constant in this
interval. Then, it is a good approximation to replace $f^{\prime}$ with $\mu$
in ${\tilde{W}}_{T}(f-f^{\prime})$ and pull it out of the integration;

$\displaystyle\langle{\tilde{S}}_{1,\rm T}^{*}(f){\tilde{S}}_{2,\rm
T}(f)\rangle$
$\displaystyle\approx\frac{1}{2}{|{\tilde{W}}_{T}(\mu-f)|}^{2}\left(\int_{\mu}^{\infty}\Gamma_{\phi}(f^{\prime})P_{\phi}(f^{\prime})df^{\prime}\right)Q(\mu).$
(2.52)

Here $Q(\mu)$ defined by

$Q(\mu)\equiv\frac{\int_{\mu}^{\infty}\Gamma_{\phi}(f)P_{\phi}(f)\varepsilon(f)df}{\int_{\mu}^{\infty}\Gamma_{\phi}(f)P_{\phi}(f)df}$
(2.53)

quantifies the effect caused by the motion of the solar system. Because
$\varepsilon(f)$ is a complex number, $Q(\mu)$ is also a complex number. Fig.
3 shows a plot of ${\rm Re}~{}Q(\mu)$ and ${\rm Im}~{}Q(\mu)$. In making this
plot, it is assumed that dark matter wind comes from the direction specified
by $(\ell,b)=(270^{\circ},0)$ in the galactic coordinates which correspond to
$(\alpha,\delta)=(3.70231,-0.81267)$ in the equilateral coordinates and the
signal is averaged over $\alpha$ to take into account the Earth’s daily
rotation. We find that the signal is suppressed at $\mu\simeq 70~{}{\rm Hz}$
($m\simeq 3\times 10^{-13}~{}{\rm eV}$) whereas it is enhanced by a factor of
$\sim 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
${\rm Im}~{}Q$ in our analysis.

Figure 3: Plot of ${\rm Re}~{}Q(\mu)$ and ${\rm Im}~{}Q(\mu)$.

It is possible to perform the integral
$\int_{\mu}^{\infty}\Gamma_{\phi}P_{\phi}df$ analytically to a good
approximation. To this end, using Eqs. (2.24) and (2.44), we obtain

$\int_{\mu}^{\infty}\Gamma_{\phi}(f)P_{\phi}(f)df\approx{(\kappa\alpha_{SiO_{2}})}^{2}\int_{\mu}^{\infty}df~{}\eta^{2}\left(a_{1}\mu^{2}L^{2}+a_{2}\eta^{2}\right)\times{\left(\frac{3}{2}\right)}^{\frac{3}{2}}\frac{4\rho_{\rm
loc}\eta}{\pi^{5/2}\mu^{3}v_{*}^{3}}\exp\left(-\frac{3\eta^{2}}{2v_{*}^{2}}\right).$
(2.54)

where we have replaced $f$ which is not in $\eta$ with $\mu$. By changing
variable from $f$ to $v$ by the relation (2.22), we have

$\int_{\mu}^{\infty}\Gamma_{\phi}(f^{\prime})P_{\phi}(f^{\prime})df^{\prime}\approx{(\kappa\alpha_{SiO_{2}})}^{2}\frac{\rho_{\rm
loc}v_{*}^{2}}{2\pi^{2}\mu^{2}}\left(2a_{1}\mu^{2}L^{2}+\frac{10}{3}a_{2}v_{*}^{2}\right).$
(2.55)

Thus, our final expression of the cross correlation becomes

$\langle{\tilde{S}}_{1,\rm T}^{*}(f){\tilde{S}}_{2,\rm
T}(f)\rangle\approx{(\kappa\alpha_{SiO_{2}})}^{2}\frac{\rho_{\rm
loc}v_{*}^{2}}{2\pi^{2}\mu^{2}}\left(a_{1}\mu^{2}L^{2}+\frac{5}{3}a_{2}v_{*}^{2}\right){|{\tilde{W}}_{T}(\mu-f)|}^{2}Q(\mu).$
(2.56)

The cross-correlation statistic ${\hat{C}}(f)$ given by the LVK collaboration
is defined as a summation over discretized frequency bins with its width
$1/T$;

${\hat{C}}(f)=\frac{1}{T^{2}\Delta f}\sum_{i=0}^{5}\frac{{\rm
Re}[{\tilde{S}_{1}^{*}}(f_{i}){\tilde{S}_{2}}(f_{i})]}{\gamma_{T}(f_{i})S_{0}(f_{i})},~{}~{}~{}~{}~{}(\Delta
f=0.03125~{}{\rm Hz},~{}T=192~{}{\rm s})$ (2.57)

where $f_{i}=f-\Delta f/2+i/T$, $\gamma_{T}$ is the overlap reduction function
for gravitational waves, and $S_{0}\equiv 3H_{0}^{2}/(10\pi^{2}f^{3})$. In the
presence of the scalar field dark matter, using Eq. (2.56), the expectation
value of ${\hat{C}}(f)$ becomes

$\langle{\hat{C}}(f)\rangle=\frac{{(\kappa\alpha_{SiO_{2}})}^{2}\rho_{\rm
loc}v_{*}^{2}}{2\pi^{2}\mu^{2}T^{2}\Delta
f}\left(a_{1}\mu^{2}L^{2}+\frac{5}{3}a_{2}v_{*}^{2}\right)\real
Q(\mu)\sum_{i=0}^{5}\frac{{|{\tilde{W}}_{T}(\mu-
f_{i})|}^{2}}{\gamma_{T}(f_{i})S_{0}(f_{i})}.$ (2.58)

This is the actual signal that we seek in the GW data provided by the LVK
Collaboration.

##  3 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
${\hat{C}}$ for each segment containing data of the duration $T$ and averaging
them #2#2#2Data is available at https://dcc.ligo.org/cgi-
bin/DocDB/ShowDocument?.submit=Identifier&docid=G2001287&version=5. Thus, the
cross correlation statistics obeys the Gaussian distribution, and the
likelihood becomes [25]

$p({\hat{C}}|{\bm{\Theta}})\propto\exp\left(-\sum_{f}\frac{{({\hat{C}}(f)-\langle{\hat{C}}(f)\rangle)}^{2}}{2\sigma^{2}(f)}\right),$
(3.1)

where ${\bm{\Theta}}=\\{\mu,\alpha_{SiO_{2}}\\}$ is a set of parameters
characterizing the scalar field model.

In our analysis, we employ the Bayesian inference. We first fix $\mu$ and
compute the posterior on $\alpha_{SiO_{2}}$ by adopting two priors (uniform in
$\alpha_{SiO_{2}}$ and uniform in $\ln\alpha_{SiO_{2}}$). We then repeat this
calculation by scanning the relevant range of $\mu$. We did not find any
conclusive evidence of a signal in data. As a result, we are able to place
upper limit on $\alpha_{SiO_{2}}$ for various values of $\mu$.

---
Figure 4: Upper limit on $\alpha_{SiO_{2}}$ obtained by LIGO O3 data. The
left/right panel assumes uniform prior for
$\alpha_{SiO_{2}}/\ln\alpha_{SiO_{2}}$.

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 $\alpha_{SiO_{2}}/(\ln\alpha_{SiO_{2}})$. Although the original data
of the cross correlation statistics is given with frequency interval $\delta
f\simeq 0.03~{}{\rm Hz}$, we did not use data for all frequencies but picked
up data only at frequency interval $10\delta 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\simeq 3\times 10^{-13}~{}{\rm
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.

|
---|---
|
|
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_{\hat{m}}=d_{m_{e}}=d_{e}=0$. Top
right panel: upper limit on $d_{\hat{m}}$ assuming $d_{g}=d_{m_{e}}=d_{e}=0$.
Middle left panel: upper limit on $d_{m_{e}}$ assuming
$d_{g}=d_{\hat{m}}=d_{e}=0$. Middle right panel: upper limit on $d_{e}$
assuming $d_{g}=d_{\hat{m}}=d_{m_{e}}=0$. Bottom panel: upper limit on $d_{g}$
assuming $d_{g}=d_{\hat{m}}=d_{m_{e}},d_{e}=0$ which is the model of the
gravitational non-minimal coupling discussed in the Appendix.

Fig. 5 shows translation of Fig. 4 into the upper limits on the four coupling
constants $(d_{g},d_{\hat{m}},d_{m_{e}},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_{m_{e}}$ and $d_{e}$, as well as other non-GW experiments including
test of violation of the weak equivalence principle (WEP) [26, 27, 21] 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 $\eta$ defined by [27]

$\eta=\frac{2(a_{A}-a_{B})}{a_{A}+a_{B}}=\alpha_{E}(\alpha_{A}-\alpha_{B})\Phi(mR_{E})\left(1+mr\right)e^{-mr}\propto(\alpha_{A}-\alpha_{B}).$
(3.2)

Here $a_{A}/(a_{B})$ is the gravitational acceleration of material A/(B)
measured at distance $r$ from the center of the Earth and $\alpha_{E}$ is
$\langle\alpha_{A}\rangle$ for materials constituting the Earth. $R_{E}$ is
the mean Earth radius and $\Phi(x)\equiv 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 MICROSCOPE 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
$\alpha_{A}\alpha_{E}$. We make an approximation that the material used in
[22] is purely aluminium (Al). The quantity $\alpha_{\rm Al}\alpha_{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_{m_{e}}$ and $d_{e}$ are zero, GEO600 does not give any
meaningful constraint in this case. Since the component-dependent part of Eq.
(2.1) 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_{\hat{m}}$ is assumed to be non-vanishing 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 $\alpha_{E}$ is less
sensitive to $d_{\hat{m}}$ than $d_{g}$, namely $\alpha_{E}\simeq
0.9d_{g}/(0.08d_{\hat{m}})$ for the first/(second) case. The third case in
which only $d_{m_{e}}$ is assumed to be non-vanishing is shown in the middle
left panel. In this case, GEO600 gives stronger constraint than LIGO for
$m\gtrsim 3\times 10^{-13}~{}{\rm eV}$. The WEP still provides the tightest
upper limit on $d_{m_{e}}$. 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.1) that component-dependent term proportional to $d_{m_{e}}$ is
suppressed by a factor $(A-2Z)/A$ compared to $d_{e}$. The fifth case in which
$d_{g}=d_{\hat{m}}=d_{m_{e}}$ 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 $\sim 5$ at large mass side. The upgraded LIGO observation is
expected to exceed the fifth force constraint [12].

##  4 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 $\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 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 $\sim 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.

## Acknowledgements

We would like to thank Yuta Michimura for helpful comments on the manuscript.
This work is supported by the MEXT KAKENHI Grant Number 17H06359 (TS),
JP21H05453 (TS), and the JSPS KAKENHI Grant Number JP19K03864 (TS).

##  Appendix A Non-minimal coupling model

We assume that the scalar field dark matter does not directly interact with
standard model fields and is only non-minimally coupled to gravity through the
kinetic term mixing between the gravitons and the scalar field. In the Jordan
frame, the action is defined as

$S[{\tilde{g}}_{\mu\nu},\Phi,\Psi]=\int
d^{4}x~{}\sqrt{-{\tilde{g}}}\left(\frac{1}{16\pi
G}{\tilde{R}}-\frac{1}{2}{\tilde{g}}^{\mu\nu}\partial_{\mu}\Phi\partial_{\nu}\Phi-{\kappa}^{-1}\zeta{\tilde{R}}\Phi-\frac{M^{2}}{2}\Phi^{2}\right)+S_{m}[{\tilde{g}}_{\mu\nu},\Psi].$
(A.1)

Here ${\tilde{g}}_{\mu\nu}$ is the Jordan frame metric, scalar field,
respectively. The standard model fields are symbolically denoted by $\Psi$.
The mass term is parametrized by $M$ and dimensionless parameter $\zeta$
represents the strength of the kinetic term mixing.

In the standard manner [28], we can move to the Einstein frame where the
gravitational action takes the Einstein-Hilbert form by the following
conformal transformation:

${\tilde{g}}_{\mu\nu}=\frac{1}{1-4\zeta\kappa\Phi}g_{\mu\nu}.$ (A.2)

In the Einstein frame, $\Phi$ is not a canonical field (i.e. coefficient of
the kinetic term is not $-\frac{1}{2}$) and the canonical field $\phi$ is
obtained from $\Phi$ by

$\phi=\int_{0}^{\Phi}\frac{\sqrt{1+12\zeta^{2}-4\zeta\kappa x}}{1-4\zeta\kappa
x}dx=\sqrt{1+12\zeta^{2}}\Phi+{\cal O}(\Phi^{2}).$ (A.3)

Then, at the leading order in $\phi$, the action in the Einstein frame becomes

$S[g_{\mu\nu},\phi,\Psi]=\int d^{4}x~{}\sqrt{-g}\left(\frac{1}{16\pi
G}R-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-\frac{m^{2}}{2}\phi^{2}\right)+S_{m}[\Omega^{2}(\phi)g_{\mu\nu},\Psi],$
(A.4)

where the mass $m$ of $\phi$ and the conformal function $\Omega(\phi)$ are
given by

$m=\frac{M}{\sqrt{1+12\zeta^{2}}},~{}~{}~{}~{}~{}~{}\Omega(\phi)=1+\frac{2\zeta}{\sqrt{1+12\zeta^{2}}}\kappa\phi+{\cal
O}(\phi^{2}).$ (A.5)

The presence of the conformal function leads to the dependence of the energy
scale $E$ inherent in $S_{m}$ on the value of $\phi$ as
$E(\phi)=\Omega(\phi)E$. Consequently, the low-energy physical constants
listed in Eqs. (2.3)-(2.6) depends on $\varphi$ as

$\displaystyle\alpha(\varphi)=\alpha,$ (A.6)
$\displaystyle\Lambda_{3}(\varphi)=\left(1+\frac{2\zeta}{\sqrt{1+12\zeta^{2}}}\varphi\right)\Lambda_{3},$
(A.7) $\displaystyle
m_{e}(\varphi)=\left(1+\frac{2\zeta}{\sqrt{1+12\zeta^{2}}}\varphi\right)m_{e},$
(A.8) $\displaystyle
m_{u,d}(\varphi)=\left(1+\frac{2\zeta}{\sqrt{1+12\zeta^{2}}}\varphi\right)m_{u,d}(\Lambda_{3}).$
(A.9)

Thus, the model defined by Eq. (A.1) corresponds to a set of parameter values

$d_{g}=d_{m_{e}}=d_{m_{u,d}}=\frac{2\zeta}{\sqrt{1+12\zeta^{2}}},~{}~{}~{}~{}~{}~{}~{}d_{e}=0.$
(A.10)

## References

  * [1] A. Arbey and F. Mahmoudi, Dark matter and the early Universe: a review, Prog. Part. Nucl. Phys. 119 (2021) 103865, [arXiv:2104.11488].
  * [2] D. B. Kaplan and M. B. Wise, Couplings of a light dilaton and violations of the equivalence principle, JHEP 08 (2000) 037, [hep-ph/0008116].
  * [3] A. Arvanitaki, J. Huang, and K. Van Tilburg, Searching for dilaton dark matter with atomic clocks, Phys. Rev. D 91 (2015), no. 1 015015, [arXiv:1405.2925].
  * [4] L. Hui, Wave Dark Matter, Ann. Rev. Astron. Astrophys. 59 (2021) 247–289, [arXiv:2101.11735].
  * [5] A. Derevianko and M. Pospelov, Hunting for topological dark matter with atomic clocks, Nature Phys. 10 (2014) 933, [arXiv:1311.1244].
  * [6] K. Van Tilburg, N. Leefer, L. Bougas, and D. Budker, Search for ultralight scalar dark matter with atomic spectroscopy, Phys. Rev. Lett. 115 (2015), no. 1 011802, [arXiv:1503.06886].
  * [7] A. Hees, J. Guéna, M. Abgrall, S. Bize, and P. Wolf, Searching for an oscillating massive scalar field as a dark matter candidate using atomic hyperfine frequency comparisons, Phys. Rev. Lett. 117 (2016), no. 6 061301, [arXiv:1604.08514].
  * [8] N. Leefer, A. Gerhardus, D. Budker, V. V. Flambaum, and Y. V. Stadnik, Search for the effect of massive bodies on atomic spectra and constraints on Yukawa-type interactions of scalar particles, Phys. Rev. Lett. 117 (2016), no. 27 271601, [arXiv:1607.04956].
  * [9] E. Savalle, A. Hees, F. Frank, E. Cantin, P.-E. Pottie, B. M. Roberts, L. Cros, B. T. Mcallister, and P. Wolf, Searching for Dark Matter with an Optical Cavity and an Unequal-Delay Interferometer, Phys. Rev. Lett. 126 (2021), no. 5 051301, [arXiv:2006.07055].
  * [10] Y. V. Stadnik and V. V. Flambaum, Enhanced effects of variation of the fundamental constants in laser interferometers and application to dark matter detection, Phys. Rev. A 93 (2016), no. 6 063630, [arXiv:1511.00447].
  * [11] A. Pierce, K. Riles, and Y. Zhao, Searching for Dark Photon Dark Matter with Gravitational Wave Detectors, Phys. Rev. Lett. 121 (2018), no. 6 061102, [arXiv:1801.10161].
  * [12] S. Morisaki and T. Suyama, Detectability of ultralight scalar field dark matter with gravitational-wave detectors, Phys. Rev. D 100 (2019), no. 12 123512, [arXiv:1811.05003].
  * [13] H. Grote and Y. V. Stadnik, Novel signatures of dark matter in laser-interferometric gravitational-wave detectors, Phys. Rev. Res. 1 (2019), no. 3 033187, [arXiv:1906.06193].
  * [14] S. Morisaki, T. Fujita, Y. Michimura, H. Nakatsuka, and I. Obata, Improved sensitivity of interferometric gravitational wave detectors to ultralight vector dark matter from the finite light-traveling time, Phys. Rev. D 103 (2021), no. 5 L051702, [arXiv:2011.03589].
  * [15] LIGO Scientific Collaboration, J. Aasi et al., Advanced LIGO, Class. Quant. Grav. 32 (2015) 074001, [arXiv:1411.4547].
  * [16] H.-K. Guo, K. Riles, F.-W. Yang, and Y. Zhao, Searching for Dark Photon Dark Matter in LIGO O1 Data, Commun. Phys. 2 (2019) 155, [arXiv:1905.04316].
  * [17] LIGO Scientific, KAGRA, Virgo Collaboration, R. Abbott et al., Constraints on dark photon dark matter using data from LIGO’s and Virgo’s third observing run, Phys. Rev. D 105 (2022), no. 6 063030, [arXiv:2105.13085].
  * [18] T. Damour and J. F. Donoghue, Equivalence Principle Violations and Couplings of a Light Dilaton, Phys. Rev. D 82 (2010) 084033, [arXiv:1007.2792].
  * [19] J. G. Williams, S. G. Turyshev, and D. Boggs, Lunar Laser Ranging Tests of the Equivalence Principle, Class. Quant. Grav. 29 (2012) 184004, [arXiv:1203.2150].
  * [20] T. A. Wagner, S. Schlamminger, J. H. Gundlach, and E. G. Adelberger, Torsion-balance tests of the weak equivalence principle, Class. Quant. Grav. 29 (2012) 184002, [arXiv:1207.2442].
  * [21] MICROSCOPE Collaboration, P. Touboul et al., MICROSCOPE Mission: Final Results of the Test of the Equivalence Principle, Phys. Rev. Lett. 129 (2022), no. 12 121102, [arXiv:2209.15487].
  * [22] E. G. Adelberger, B. R. Heckel, and A. E. Nelson, Tests of the gravitational inverse square law, Ann. Rev. Nucl. Part. Sci. 53 (2003) 77–121, [hep-ph/0307284].
  * [23] Y. V. Stadnik and V. V. Flambaum, Can dark matter induce cosmological evolution of the fundamental constants of Nature?, Phys. Rev. Lett. 115 (2015), no. 20 201301, [arXiv:1503.08540].
  * [24] S. M. Vermeulen et al., Direct limits for scalar field dark matter from a gravitational-wave detector, Nature 600 (2021), no. 7889 424–428, [arXiv:2103.03783].
  * [25] V. Mandic, E. Thrane, S. Giampanis, and T. Regimbau, Parameter Estimation in Searches for the Stochastic Gravitational-Wave Background, Phys. Rev. Lett. 109 (2012) 171102, [arXiv:1209.3847].
  * [26] S. Schlamminger, K. Y. Choi, T. A. Wagner, J. H. Gundlach, and E. G. Adelberger, Test of the equivalence principle using a rotating torsion balance, Phys. Rev. Lett. 100 (2008) 041101, [arXiv:0712.0607].
  * [27] J. Bergé, P. Brax, G. Métris, M. Pernot-Borràs, P. Touboul, and J.-P. Uzan, MICROSCOPE Mission: First Constraints on the Violation of the Weak Equivalence Principle by a Light Scalar Dilaton, Phys. Rev. Lett. 120 (2018), no. 14 141101, [arXiv:1712.00483].
  * [28] K.-i. Maeda, Towards the Einstein-Hilbert Action via Conformal Transformation, Phys. Rev. D 39 (1989) 3159.


# On solutions of an ill-posed Stefan problem

Evgeny Yu. Panov
Yaroslav-the-Wise Novgorod State University,
Veliky Novgorod, Russian Federation

###### Abstract

We study multi-phase Stefan problem with increasing Riemann initial data and
with generally negative latent specific heats for the phase transitions. We
propose the variational formulation of self-similar solutions, which allows to
find precise conditions for existence and uniqueness of the solution.
Bibliography$:$ $2$ titles.

##  1 Introduction

In a half-plane $\Pi=\\{\ (t,x)\ |\ t>0,x\in{\mathbb{R}}\ \\}$ we consider the
multi-phase Stefan problem for the heat equation

$u_{t}=a_{i}^{2}u_{xx},\quad u_{i}<u<u_{i+1},$ (1)

where $u_{-}=u_{0}<u_{1}<\cdots<u_{n}<u_{n+1}=u_{+}$, $u_{i}$, $i=1,\dots,n$
being the temperatures of phase transitions, $a_{i}>0$, $i=0,\dots,n$, are the
diffusivity constants. On the unknown lines $x=x_{i}(t)$ of phase transitions
where $u=u_{i}$ the following Stefan condition

$d_{i}x_{i}^{\prime}(t)+k_{i}u_{x}(t,x_{i}(t)+)-k_{i-1}u_{x}(t,x_{i}(t)-)=0$
(2)

is postulated, where $k_{i}>0$ are the thermal conductivity of the $i$-th
phase, while $d_{i}$ is the Stefan number (the latent specific heat) for the
$i$-th phase transition. In the case when the temperature $u(t,x)$ non-
decreases with respect to the spatial variable $x$ the Stefan numbers $d_{i}$
should be nonnegative by physical reasons. In this case the problem (1), (2)
is well-posed and reduces to a degenerate nonlinear diffusion equation, see
[3, Chapter 5]. In the present paper we consider the case of arbitrary $d_{i}$
when the problem is generally ill-posed. This was firstly demonstrated in [4]
for a one-phase Stefan problem with a negative Stefan number. Our aim is to
find general conditions which guarantees the correctness of our problem. We
will study the Cauchy problem with the Riemann initial data

$u(0,x)=\left\\{\begin{array}[]{lr}u_{-},&x<0,\\\
u_{+},&x>0.\end{array}\right.$ (3)

By the invariance of our problem under the transformation group
$(t,x)\to(\lambda^{2}t,\lambda x)$, $\lambda\in{\mathbb{R}}$, $\lambda\not=0$,
it is natural to seek a self-similar solution of problem (1), (2), (3), which
has the form $u(t,x)=v(\xi)$, $\xi=x/\sqrt{t}$. For the heat equation
$u_{t}=a^{2}u_{xx}$ a self-similar solution must satisfy the linear ODE
$a^{2}v^{\prime\prime}=-\xi v^{\prime}/2$, the general solution of which is

$v=C_{1}F(\xi/a)+C_{2},\ C_{1},C_{2}={\rm const},\mbox{ where
}F(\xi)=\frac{1}{2\sqrt{\pi}}\int_{-\infty}^{\xi}e^{-s^{2}/4}ds.$

This allows to write our solution in the form

$\displaystyle
v(\xi)=u_{i}+\frac{u_{i+1}-u_{i}}{F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})}(F(\xi/a_{i})-F(\xi_{i}/a_{i})),$
(4) $\displaystyle\xi_{i}<\xi<\xi_{i+1},\ i=0,\ldots,n,$

where $-\infty=\xi_{0}<\xi_{1}<\cdots<\xi_{n}<\xi_{n+1}=+\infty$ and we agree
that $F(-\infty)=0$, $F(+\infty)=1$. The parabolas $\xi=\xi_{i}$,
$i=1,\ldots,n$, where $u=u_{i}$, are free boundaries. They must be determined
by conditions (2). In the variable $\xi$ these conditions have the form

$d_{i}\xi_{i}/2+\frac{k_{i}(u_{i+1}-u_{i})F^{\prime}(\xi_{i}/a_{i})}{a_{i}(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))}-\frac{k_{i-1}(u_{i}-u_{i-1})F^{\prime}(\xi_{i}/a_{i-1})}{a_{i-1}(F(\xi_{i}/a_{i-1})-F(\xi_{i-1}/a_{i-1}))}=0,\
i=1,\ldots,n.$ (5)

To investigate this nonlinear system, we notice that (5) coincides with the
condition $\nabla E(\bar{\xi})=0$, where the function

$\displaystyle
E(\bar{\xi})=-\sum_{i=0}^{n}k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))+\sum_{i=1}^{n}d_{i}\xi_{i}^{2}/4,$
(6) $\displaystyle\bar{\xi}=(\xi_{1},\ldots,\xi_{n})\in\Omega,$

the open convex domain $\Omega\subset{\mathbb{R}}^{n}$ is given by the
inequalities $\xi_{1}<\cdots<\xi_{n}$. Observe that $E(\bar{\xi})\in
C^{\infty}(\Omega)$.

##  2 Coercivity of $E$

Let us introduce the sub-level sets

$\Omega_{c}=\\{\ \bar{\xi}\in\Omega\ |\ E(\bar{\xi})\leq c\ \\},\quad
c\in{\mathbb{R}}.$

We are going to prove that under some exact conditions on the parameters these
sets are compact. This means the coercivity of function $E(\bar{\xi})$. First,
we investigate when the sets $\Omega_{c}$ are bounded.

######  Lemma 1.

The sets $\Omega_{c}$ are bounded if and only if the following two conditions
are satisfied

$\displaystyle
S^{j}\doteq\sum_{i=1}^{j}[\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}]\geq 0,\
j=1,\ldots,n;$ (7) $\displaystyle
S_{j}\doteq\sum_{i=j}^{n}[\kappa_{i}(u_{i+1}-u_{i})+d_{i}]\geq 0,\
j=1,\ldots,n,$ (8)

where $\kappa_{i}=k_{i}/a_{i}^{2}$.

###### Proof.

Assume that conditions (7), (8) hold. We have to prove that the set
$\Omega_{c}$ is bounded. Assuming the contrary, we can find the sequence
$\bar{\xi}_{m}\in\Omega_{c}$, $m\in{\mathbb{N}}$, such that
$|\bar{\xi}_{m}|\to+\infty$ as $m\to\infty$ (here and in the sequel we denote
by $|v|$ the Euclidean norm of a finite-dimensional vector $v$). Passing to a
subsequence if necessary, we can assume that the coordinates $\xi_{mi}$,
$i=1,\ldots,n$, of the points $\bar{\xi}_{m}$ have finite or infinite limits
in $[-\infty,+\infty]$. Let $l$ be the minimal $i$ such that
$\displaystyle\lim_{m\to\infty}\xi_{mi}=+\infty$. If such $i$ does not exist
we set $l=n+1$. Similarly, we define $r$ as the maximal $i$ such that
$\displaystyle\lim_{m\to\infty}\xi_{mi}=-\infty$. In the case when there is no
such $i$ we set $r=0$. By the condition $|\bar{\xi}_{m}|\to+\infty$ we can
claim that either $l\leq n$ or $r\geq 1$. Since the coordinates $\xi_{mi}$
increase with respect to $i$ then

$\lim_{m\to\infty}\xi_{mi}=+\infty\ \forall i,l\leq i\leq
n;\quad\lim_{m\to\infty}\xi_{mi}=-\infty\ \forall i,1\leq i\leq r$

and, in particular, $r<l$. Now we observe that by L’Hôpital’s rule and the
identities $F^{\prime}(\xi)=\frac{1}{2\sqrt{\pi}}e^{-\xi^{2}/4}$,
$F^{\prime\prime}(x)=-\xi F^{\prime}(\xi)/2$,

$\displaystyle\lim_{\xi\to+\infty}\frac{1-F(\xi)}{2F^{\prime}(\xi)/\xi}=\lim_{\xi\to+\infty}\frac{F^{\prime}(\xi)\xi^{2}}{2(F^{\prime}(\xi)-\xi
F^{\prime\prime}(\xi))}=\lim_{\xi\to\infty}\frac{\xi^{2}}{\xi^{2}+2}=1,$
$\displaystyle\lim_{\xi\to-\infty}\frac{-F(\xi)}{2F^{\prime}(\xi)/\xi}=\lim_{\xi\to-\infty}\frac{F^{\prime}(\xi)\xi^{2}}{2(F^{\prime}(\xi)-\xi
F^{\prime\prime}(\xi))}=\lim_{\xi\to\infty}\frac{\xi^{2}}{\xi^{2}+2}=1.$

By these relations

$F(-\xi)=1-F(\xi)=\frac{1}{\sqrt{\pi}|\xi|}e^{-\xi^{2}/4}(1+o(1))\ \mbox{ as
}\xi\to+\infty.\\\ $ (9)

Since

$-\ln(F(\xi_{m(i+1)}/a_{i})-F(\xi_{mi}/a_{i}))\geq\max(-\ln(1-F(\xi_{mi}/a_{i})),-\ln
F(\xi_{m(i+1)}/a_{i}))\geq 0,$

we have

$\displaystyle
E(\bar{\xi}_{m})=\sum_{i=0}^{r-1}[-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{m(i+1)}/a_{i})-F(\xi_{mi}/a_{i}))+d_{i+1}\xi_{m(i+1)}^{2}/4]+$
$\displaystyle\sum_{i=l}^{n}[-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{m(i+1)}/a_{i})-F(\xi_{mi}/a_{i}))+d_{i}\xi_{mi}^{2}/4]-$
$\displaystyle\sum_{i=r}^{l-1}k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{m(i+1)}/a_{i})-F(\xi_{mi}/a_{i}))+\sum_{i=r+1}^{l-1}d_{i}\xi_{mi}^{2}/4\geq$
$\displaystyle\sum_{i=0}^{r-1}[-k_{i}(u_{i+1}-u_{i})\ln
F(\xi_{m(i+1)}/a_{i})+d_{i+1}\xi_{m(i+1)}^{2}/4]+$
$\displaystyle\sum_{i=l}^{n}[-k_{i}(u_{i+1}-u_{i})\ln(1-F(\xi_{mi}/a_{i}))+d_{i}\xi_{mi}^{2}/4]+p_{m}=$
$\displaystyle\sum_{i=1}^{r}[-k_{i-1}(u_{i}-u_{i-1})\ln
F(\xi_{mi}/a_{i-1})+d_{i}\xi_{mi}^{2}/4]+$
$\displaystyle\sum_{i=l}^{n}[-k_{i}(u_{i+1}-u_{i})\ln(1-F(\xi_{mi}/a_{i}))+d_{i}\xi_{mi}^{2}/4]+p_{m},$
(10)

where

$p_{m}=\sum_{i=r+1}^{l-1}d_{i}\xi_{mi}^{2}/4$

is a bounded sequence, because the sequences $\xi_{mi}$ have finite limits as
$m\to\infty$ for all $i=r+1,\ldots l-1$. Remark that the set of such $i$ may
be empty, in this case $p_{m}=0$. In view of (9), it follows from (2) that

$\displaystyle c\geq
E(\bar{\xi}_{m})\geq\sum_{i=1}^{r}[\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}]\xi_{mi}^{2}/4+\sum_{i=l}^{n}[\kappa_{i}(u_{i+1}-u_{i})+d_{i}]\xi_{mi}^{2}/4+q_{m},$
(11)

where, up to a vanishing sequence,

$q_{m}=\sum_{i=1}^{r}k_{i-1}(u_{i}-u_{i-1})\ln(\sqrt{\pi}|\xi_{mi}|/a_{i-1})+\sum_{i=l}^{n}k_{i}(u_{i+1}-u_{i})\ln(\sqrt{\pi}|\xi_{mi}|/a_{i})+p_{m}.$

In particular, the sequence $q_{m}\to+\infty$ as $m\to\infty$ and it follows
from (11) that

$\sum_{i=1}^{r}[\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}]\xi_{mi}^{2}/4+\sum_{i=l}^{n}[\kappa_{i}(u_{i+1}-u_{i})+d_{i}]\xi_{mi}^{2}/4\leq
c-q_{m}<0$ (12)

for all sufficiently large $m\in{\mathbb{N}}$. On the other hand,
$\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}=S^{i}-S^{i-1}$,
$\kappa_{i}(u_{i+1}-u_{i})+d_{i}=S_{i}-S_{i+1}$ and applying the summation by
parts formula, we find that

$\displaystyle\sum_{i=1}^{r}[\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}]\xi_{mi}^{2}/4=S^{r}\xi_{mr}^{2}/4-\sum_{i=1}^{r-1}S^{i}(\xi_{m(i+1)}^{2}-\xi_{mi}^{2})/4\geq$
$\displaystyle-\sum_{i=1}^{r-1}S^{i}(\xi_{m(i+1)}-\xi_{mi})(\xi_{m(i+1)}+\xi_{mi})/4\geq
0,$
$\displaystyle\sum_{i=l}^{n}[\kappa_{i}(u_{i+1}-u_{i})+d_{i}]\xi_{mi}^{2}/4=S_{l}\xi_{ml}^{2}/4+\sum_{i=l+1}^{n}S_{i}(\xi_{mi}^{2}-\xi_{m(i-1)}^{2})/4\geq$
$\displaystyle\sum_{i=l+1}^{n}S_{i}(\xi_{mi}-\xi_{m(i-1)})(\xi_{mi}+\xi_{m(i-1)})/4\geq
0$

for sufficiently large $m$ because $S_{i},S^{i}\geq 0$ by assumptions (7), (8)
while $\xi_{mi}<0$ for $i=1,\ldots,r$, $\xi_{mi}>0$ for $i=l,\ldots,n$, if $m$
is sufficiently large. But this contradicts to (12). Hence, the set
$\Omega_{c}$ is bounded for every real $c$.

Conversely, assume that at least one of conditions (7), (8) fails. For
definiteness we suppose that (7) fails. This means that for some $1\leq r\leq
n$

$S^{r}=\sum_{i=1}^{r}[\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}]<0.$ (13)

For $\sigma>0$ we define the point $\bar{\xi}(\sigma)\in\Omega$ with
coordinates $\xi_{i}=-\sigma+i-r$, $i=1,\ldots,r$; $\xi_{i}=i-r$,
$i=r+1,\ldots,n$. Since for $i=1,\ldots,r-1$

$F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})=\frac{1}{2\sqrt{\pi}}\int_{\xi_{i}/a_{i}}^{\xi_{i+1}/a_{i}}e^{-s^{2}/4}ds\geq\frac{\xi_{i+1}-\xi_{i}}{2\sqrt{\pi}a_{i}}e^{-\xi_{i}^{2}/4a_{i}^{2}}=\frac{1}{2\sqrt{\pi}a_{i}}e^{-(\sigma+r-i)^{2}/4a_{i}^{2}},$

and, in view of (9),

$F(\xi_{1}/a_{0})=\frac{a_{0}}{\sqrt{\pi}|\sigma+r-1|}e^{-(\sigma+r-1)^{2}/4a_{0}^{2}}(1+o(1))\
\mbox{ as }\sigma\to+\infty$

then for some constant $C$

$\displaystyle
E(\bar{\xi}(\sigma))\leq-\sum_{i=0}^{r-1}k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))+$
$\displaystyle\sum_{i=1}^{r}d_{i}\xi_{i}^{2}/4+C\leq\sum_{i=1}^{r}[\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}]\sigma^{2}/4+O(\sigma).$
(14)

By (13) it follows from (2) that

$E(\bar{\xi}(\sigma))\to-\infty\ \mbox{ as }\sigma\to+\infty.$

Therefore, $\bar{\xi}(\sigma)\in\Omega_{c}$ for all sufficiently large
$\sigma$ and the set $\Omega_{c}$ is unbounded. Thus, if the set $\Omega_{c}$
is bounded (for some $c$) then both conditions (7), (8) are satisfied. This
completes the proof. ∎

Now, we are going to demonstrate that these conditions are necessary and
sufficient for coercivity of the function $E$.

######  Theorem 1.

The sets $\Omega_{c}$ are compact for each $c\in{\mathbb{R}}$ if and only if
conditions (7), (8) are satisfied. In particular, under these conditions the
function $E(\bar{\xi})$ reaches its minimal value.

###### Proof.

If the set $\Omega_{c}$ is compact for some $c\in{\mathbb{R}}$ then it is
bounded and by Lemma 1 conditions (7), (8) are satisfied.

Conversely, if conditions (7), (8) hold then the set $\Omega_{c}$ is bounded
(for every fixed $c$). Therefore, there exists $R>0$ such that
$|\bar{\xi}|\leq R$ for all $\bar{\xi}\in\Omega_{c}$. This implies the
estimate

$\displaystyle-\sum_{i=0}^{n}k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))\leq
E(\bar{\xi})-\sum_{i=1}^{n}d_{i}\xi_{i}^{2}/4\leq$ $\displaystyle c_{1}\doteq
c+\frac{R^{2}}{4}\sum_{i=1}^{n}\max(-d_{i},0).$

Since $-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))>0$, we
claim that for all $i=0,\ldots,n$

$-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i}))\leq c_{1}.$

It follows that

$F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})\geq\delta\doteq\exp(-c_{1}/\alpha)>0,$
(15)

where $\displaystyle\alpha=\min_{i=0,\ldots,n}k_{i}(u_{i+1}-u_{i})>0$. Since
$F^{\prime}(\xi)=\frac{1}{2\sqrt{\pi}}e^{-\xi^{2}/4}<1$, the function $F(\xi)$
is Lipschitz with constant $1$, and it follows from (15) that

$(\xi_{i+1}-\xi_{i})/a_{i}\geq
F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})\geq\delta,\quad i=1,\ldots,n-1,$

and we obtain the estimates $\xi_{i+1}-\xi_{i}\geq\delta_{1}=\delta\min
a_{i}$. Thus, the set $\Omega_{c}$ is contained in a compact

$K=\\{\ \bar{\xi}=(\xi_{1},\ldots,\xi_{n})\in{\mathbb{R}}^{n}\ |\
|\bar{\xi}|\leq R,\ \xi_{i+1}-\xi_{i}\geq\delta_{1}\ \forall i=1,\ldots,n-1\
\\}.$

Since $E(\bar{\xi})$ is continuous on $K$, the set $\Omega_{c}$ where
$E(\bar{\xi})\leq c$ is a closed subset of $K$ and therefore is compact. For
$c>N\doteq\inf E(\bar{\xi})$, this set is not empty and the entropy
$E(\bar{\xi})$ reaches on it a minimal value, which is evidently equal $N$. ∎

We have established the existence of minimal value $E(\bar{\xi}_{0})=\min
E(\bar{\xi})$. At the point $\bar{\xi}_{0}$ the required condition $\nabla
E(\bar{\xi}_{0})=0$ is satisfied, and $\bar{\xi}_{0}$ is a solution of system
(5). The coordinates of $\bar{\xi}_{0}$ determine the solution (4) of our
Stefan-Riemann problem. Thus, we establish the following existence result.

######  Theorem 2.

Under conditions (7), (8) there exists a self-similar increasing solution of
problem (1), (2), (3).

##  3 Convexity of the Function $E$ and Uniqueness of the Solution

In this section we find a condition of strict convexity of the function
$E(\bar{\xi})$. Since a strictly convex function can have at most one critical
point (and it is necessarily a global minimum) the system (5) has at most one
solution, that is, a self-similar solution (4) of problem (1), (2), (3) is
unique. We introduce the functions

$E_{i}(\bar{\xi})=-k_{i}(u_{i+1}-u_{i})\ln(F(\xi_{i+1}/a_{i})-F(\xi_{i}/a_{i})),\quad
i=0,\ldots,n,$

which are exactly summands of the first sum in expression (6). They depend
only on two coordinates $\xi_{i},\xi_{i+1}$ if $i=1,\ldots,n-1$ or on one
coordinate $\xi_{1}$, $\xi_{n}$ if $i=0,n$, respectively. We compute the
second derivatives

$\frac{\partial^{2}E_{i}(\bar{\xi})}{\partial\xi_{i}^{2}}=\kappa_{i}(u_{i+1}-u_{i})\left(\frac{F^{\prime\prime}(y_{i})}{F(x_{i})-F(y_{i})}+\frac{(F^{\prime}(y_{i}))^{2}}{(F(x_{i})-F(y_{i}))^{2}}\right),$

where we denote $x_{i}=\xi_{i+1}/a_{i}$, $y_{i}=\xi_{i}/a_{i}$. Since
$F^{\prime\prime}(x)=-\frac{x}{2}F^{\prime}(x)$ we can transform the above
relation

$\displaystyle\frac{\partial^{2}E_{i}(\bar{\xi})}{\partial\xi_{i}^{2}}=\frac{\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(y_{i})}{F(x_{i})-F(y_{i})}\left(-\frac{y_{i}}{2}-\frac{F^{\prime}(x_{i})-F^{\prime}(y_{i})}{F(x_{i})-F(y_{i})}\right)+$
$\displaystyle\frac{\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(y_{i})F^{\prime}(x_{i})}{(F(x_{i})-F(y_{i}))^{2}}=\beta_{i}^{-}+\gamma_{i},$
(16)

where

$\displaystyle\beta_{i}^{-}=\frac{\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(y_{i})}{F(x_{i})-F(y_{i})}\left(-\frac{y_{i}}{2}-\frac{F^{\prime}(x_{i})-F^{\prime}(y_{i})}{F(x_{i})-F(y_{i})}\right),\
i=1,\ldots,n;$
$\displaystyle\gamma_{i}=\frac{\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(y_{i})F^{\prime}(x_{i})}{(F(x_{i})-F(y_{i}))^{2}},\
i=0,\ldots,n.$

Notice that $x_{n}=+\infty$, $y_{0}=-\infty$ and therefore $F(x_{n})=1$,
$F(y_{0})=0$, $F^{\prime}(y_{0})=F^{\prime}(x_{n})=0$. In particular,
$\gamma_{0}=\gamma_{n}=0$. By the similar computations we find

$\frac{\partial^{2}E_{i}(\bar{\xi})}{\partial\xi_{i+1}^{2}}=\beta_{i}^{+}+\gamma_{i},\quad\frac{\partial^{2}E_{i}(\bar{\xi})}{\partial\xi_{i}\partial\xi_{i+1}}=-\gamma_{i},$
(17) $\mbox{ where }\
\beta_{i}^{+}=\frac{\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(x_{i})}{F(x_{i})-F(y_{i})}\left(\frac{x_{i}}{2}+\frac{F^{\prime}(x_{i})-F^{\prime}(y_{i})}{F(x_{i})-F(y_{i})}\right),\quad
i=0,\ldots,n-1.$

To estimate the values $\beta_{i}^{\pm}$, we use the following simple lemma.

######  Lemma 2.

For each $x>y$ there exist such $\theta_{1},\theta_{2}\in(y,x)$ that

$\displaystyle\frac{1}{F(x)-F(y)}\left(-\frac{y}{2}-\frac{F^{\prime}(x)-F^{\prime}(y)}{F(x)-F(y)}\right)=\frac{1}{4F^{\prime}(\theta_{1})},$
(18)
$\displaystyle\frac{1}{F(x)-F(y)}\left(\frac{x}{2}+\frac{F^{\prime}(x)-F^{\prime}(y)}{F(x)-F(y)}\right)=\frac{1}{4F^{\prime}(\theta_{2})}.$
(19)

###### Proof.

Observe that the left hand side of (18) can be written in the form

$\frac{-\frac{y}{2}(F(x)-F(y))-(F^{\prime}(x)-F^{\prime}(y))}{(F(x)-F(y))^{2}}=\frac{g(x)-g(y)}{h(x)-h(y)},$

where

$g(y)=-\frac{y}{2}(F(x)-F(y))-(F^{\prime}(x)-F^{\prime}(y)),\quad
h(y)=(F(x)-F(y))^{2}.$

By the Cauchy mean value theorem there exists such a value
$\theta_{1}\in(y,x)$ that

$\frac{g(x)-g(y)}{h(x)-h(y)}=\frac{g^{\prime}(\theta_{1})}{h^{\prime}(\theta_{1})}.$

Using again the identity $F^{\prime\prime}(y)=-y/2F^{\prime}(y)$, we find that

$g^{\prime}(y)=-(F(x)-F(y))/2+yF^{\prime}(y)/2+F^{\prime\prime}(y)=-(F(x)-F(y))/2,\quad
h^{\prime}(y)=-2(F(x)-F(y))F^{\prime}(y).$

Therefore,
$g^{\prime}(\theta_{1})/h^{\prime}(\theta_{1})=1/(4F^{\prime}(\theta_{1}))$
and (18) follows. Similarly, the left hand side of (19) can be represented as
the ratio $(g_{1}(x)-g_{1}(y))/(h(x)-h(y))$, where

$g_{1}(x)=\frac{x}{2}(F(x)-F(y))+(F^{\prime}(x)-F^{\prime}(y)),\quad
h(x)=(F(x)-F(y))^{2}.$

Applying again the Cauchy mean value theorem, we arrive at (19). Notice also
that (19) reduce to (18) after the change $x^{\prime}=-y$, $y^{\prime}=-x$. ∎

######  Remark 1.

As is easy to realize from the proof, relations (18), (19) remain valid,
respectively, for $x=+\infty$, $y=-\infty$.

######  Corollary 1.

There exist values $\theta_{1i},\theta_{2i}\in(y_{i},x_{i})$ such that

$\beta_{i}^{-}=\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(y_{i})/(4F^{\prime}(\theta_{1i})),\quad\beta_{i}^{+}=\kappa_{i}(u_{i+1}-u_{i})F^{\prime}(x_{i})/(4F^{\prime}(\theta_{2i})).$
(20)

In particular, $\beta_{i}^{\pm}>0$ and for $i=1,\ldots,n$

$\beta_{i}^{-}+\beta_{i-1}^{+}>\frac{1}{4}\min(\kappa_{i}(u_{i+1}-u_{i}),\kappa_{i-1}(u_{i}-u_{i-1})).$
(21)

###### Proof.

Inequalities (20) directly follows from Lemma 1 and the definitions of
$\beta_{i}^{\pm}$. Further, if $\xi_{i}\geq 0$ then
$\theta_{1i}>y_{i}=\xi_{i}/a_{i}\geq 0$, which implies that
$F^{\prime}(y_{i})/F^{\prime}(\theta_{1i})>1$. Therefore,
$\beta_{i}^{-}>\kappa_{i}(u_{i+1}-u_{i})/4$ and (21) follows. Similarly, in
the case $\xi_{i}\leq 0$ we have $\theta_{2(i-1)}<x_{i-1}=\xi_{i}/a_{i-1}\leq
0$ and therefore $F^{\prime}(x_{i-1})/F^{\prime}(\theta_{2(i-1)})>1$. Hence
$\beta_{i-1}^{+}>\kappa_{i-1}(u_{i}-u_{i-1})/4$ and (21) is again fulfilled. ∎

Now we are ready to formulate a sufficient condition for convexity of the
function $E$.

######  Theorem 3.

Suppose that for all $i=1,\ldots,n$

$\min(\kappa_{i}(u_{i+1}-u_{i}),\kappa_{i-1}(u_{i}-u_{i-1}))+2d_{i}\geq 0.$
(22)

Then the function $E(\bar{\xi})$ is strictly convex on $\Omega$.

###### Proof.

We fix a point $\bar{\xi}_{0}=(\xi_{01},\ldots,\xi_{0n})\in\Omega$ and define
values $d_{i}^{-}$, $d_{i}^{+}$ such that $d_{i}^{-}+d_{i}^{+}=d_{i}$ and that

$\beta_{i}^{-}+d_{i}^{-}/2>0,\quad\beta_{i-1}^{+}+d_{i}^{+}/2>0,$ (23)

where $\beta_{i}^{\pm}$ correspond to $\bar{\xi}_{0}$. Since, in view of (21),
(22)

$\beta_{i}^{-}+\beta_{i-1}^{+}+d_{i}/2>0,$

such values really exist, for instance we can take

$d_{i}^{-}=\frac{d_{i}\beta_{i}^{-}}{\beta_{i}^{-}+\beta_{i-1}^{+}},\quad
d_{i}^{+}=\frac{d_{i}\beta_{i-1}^{+}}{\beta_{i}^{-}+\beta_{i-1}^{+}}.$

We introduce the functions

$\displaystyle
P_{i}(\bar{\xi})=E_{i}(\xi_{i},\xi_{i+1})+d_{i}^{-}\xi_{i}^{2}/4+d_{i+1}^{+}\xi_{i+1}^{2}/4,\
i=1,\ldots,n-1;$ $\displaystyle
P_{0}(\bar{\xi})=E_{0}(\xi_{1})+d_{1}^{+}\xi_{1}^{2}/4,\quad
P_{n}(\bar{\xi})=E_{n}(\xi_{n})+d_{n}^{-}\xi_{n}^{2}/4.$

Observe that

$\sum_{i=0}^{n}P_{i}(\bar{\xi})=\sum_{i=0}^{n}E_{i}(\bar{\xi})+\sum_{i=1}^{n}(d_{i}^{-}+d_{i}^{+})\xi_{i}^{2}/4=\sum_{i=0}^{n}E_{i}(\bar{\xi})+\sum_{i=1}^{n}d_{i}\xi_{i}^{2}/4=E(\bar{\xi}).$

First, we observe that by relations (3), (17), (23)

$\frac{\partial^{2}P_{0}(\bar{\xi}_{0})}{\partial\xi_{1}^{2}}=P_{0}^{\prime\prime}(\xi_{01})=\beta_{0}^{+}+d_{1}^{+}/2>0,\quad\frac{\partial^{2}P_{n}(\bar{\xi}_{0})}{\partial\xi_{n}^{2}}=P_{n}^{\prime\prime}(\xi_{0n})=\beta_{n}^{-}+d_{n}/2>0.$
(24)

For $i=1,\ldots,n-1$ the function $P_{i}$ depend only on two variables
$\xi_{i},\xi_{i+1}$, $P_{i}(\bar{\xi})=P_{i}(\xi_{i},\xi_{i+1})$. Denote by
$D_{2}^{2}P_{i}$ the Hessian of $P_{i}$ as a function of two variables. In
view of (3), (17) we find

$\frac{\partial^{2}P_{i}(\bar{\xi}_{0})}{\partial\xi_{i}^{2}}=\beta_{i}^{-}+d_{i}^{-}/2+\gamma_{i},\
\frac{\partial^{2}P_{i}(\bar{\xi}_{0})}{\partial\xi_{i+1}^{2}}=\beta_{i}^{+}+d_{i+1}^{+}/2+\gamma_{i},\
\frac{\partial^{2}P_{i}(\bar{\xi}_{0})}{\partial\xi_{i}\partial\xi_{i+1}}=-\gamma_{i},$

and the matrix $D_{2}^{2}P_{i}$ can be represented as the sum
$R_{1}+\gamma_{i}R_{2}$, where $R_{1}$ is a diagonal matrix with the positive
(in view of (23)) diagonal elements $\beta_{i}^{-}+d_{i}^{-}/2$,
$\beta_{i}^{+}+d_{i+1}^{+}/2$ while $R_{2}=\left(\begin{smallmatrix}1&-1\\\
-1&1\end{smallmatrix}\right)$. Since $R_{1}>0$, $R_{2}\geq 0$, then the matrix
$D_{2}^{2}P_{i}>0$ (strictly positive definite). This, together with (24),
implies that the “full“ Hessian $D^{2}P_{i}(\bar{\xi}_{0})\geq 0$ for all
$i=0,\ldots,n$ and therefore

$D^{2}E(\bar{\xi}_{0})=\sum_{i=0}^{n}D^{2}P_{i}(\bar{\xi}_{0})\geq 0.$

Let us prove that $D^{2}E(\bar{\xi}_{0})>0$. Assume that for some
$\zeta=(\zeta_{1},\ldots,\zeta_{n})\in{\mathbb{R}}^{n}$

$D^{2}E(\bar{\xi}_{0})\zeta\cdot\zeta=\sum_{i,j=1}^{n}\frac{\partial^{2}E(\bar{\xi}_{0})}{\partial\xi_{i}\partial\xi_{j}}\zeta_{i}\zeta_{j}=0.$
(25)

Since $D^{2}E(\bar{\xi}_{0})$ is a sum of nonnegative matrices
$D^{2}P_{i}(\bar{\xi}_{0})$, we find that

$D^{2}P_{i}(\bar{\xi}_{0})\zeta\cdot\zeta=0\quad\forall i=0,\ldots,n.$ (26)

Taking $i=0,n$ we derive with the help of (24) that $\zeta_{1}=\zeta_{n}=0$.
In the case $i=1,\ldots,n-1$ (26) reduces to the equality

$D_{2}^{2}P_{i}(\bar{\xi}_{0})v_{i}\cdot v_{i}=0,\quad
v_{i}\doteq(\zeta_{i},\zeta_{i+1})\in{\mathbb{R}}^{2}.$

Since $D_{2}^{2}P_{i}(\bar{\xi}_{0})>0$, this implies that $v_{i}=0$. In
particular, $\zeta_{i}=0$ and we claim that $\zeta=0$. Thus, (26) may hold
only for $\zeta=0$, that is, the Hessian $D^{2}E(\bar{\xi}_{0})>0$. Taking
into account that a point $\bar{\xi}_{0}\in\Omega$ is arbitrary, we conclude
that the function $E(\bar{\xi})$ is strictly convex. ∎

######  Theorem 4.

Under requirement (22) there exists a unique increasing self-similar solution
(4) of problem (1), (2), (3).

###### Proof.

Obviously, condition (22) implies that

$\displaystyle\kappa_{i}(u_{i+1}-u_{i})+d_{i}=\kappa_{i}(u_{i+1}-u_{i})/2+(\kappa_{i}(u_{i+1}-u_{i})+2d_{i})/2>0,$
$\displaystyle\kappa_{i-1}(u_{i}-u_{i-1})+d_{i}=\kappa_{i-1}(u_{i}-u_{i-1})/2+(\kappa_{i-1}(u_{i}-u_{i-1})+2d_{i})/2>0$

and, in particular, conditions (7), (8) are satisfied. By Theorem 2 solution
(4) of problem (1), (2), (3) exists. Uniqueness of this solution follows from
strict convexity of the function $E(\bar{\xi})$ claimed in Theorem 3. ∎

## Acknowledgments

The research was supported by the Russian Science Foundation, grant
22-21-00344.

## References

  * [1]
  * [2]
  * [3] O. A. Ladyzhenskaya, V. A. Solonnikov and N. N. Ural’tseva, Linear and Quasi-Linear Equations of Parabolic Type, AMS, Providence (1968).
  * [4] V. V. Pukhnachov, Generation of a singularity in the solution of a Stefan-type problem [in Russian], Differ. Uravn., 16, No. 3, 492–500 (1980).