# Probing Ultralight Bosons with Compact Eccentric Binaries

Boye Su<EMAIL_ADDRESS>Department of Physics, Tsinghua
University, Beijing 100084, China  Zhong-Zhi Xianyu<EMAIL_ADDRESS>Department of Physics, Tsinghua University, Beijing 100084, China  Xingyu
Zhang<EMAIL_ADDRESS>Department of Physics, Tsinghua
University, Beijing 100084, China

###### Abstract

Ultralight bosons can be abundantly produced through superradiance process by
a spinning black hole and form a bound state with hydrogen-like spectrum. We
show that such a “gravitational atom” typically possesses anomalously large
mass quadrupole and leads to significant orbital precession when it forms an
eccentric binary with a second compact object. Dynamically formed black hole
binaries or pulsar-black hole binaries are typically eccentric during their
early inspirals. We show that the large orbital precession can generate
distinct and observable signature in their gravitational wave or pulsar timing
signals.

A new era of gravitational wave (GW) astronomy has emerged following the LIGO
discovery of two coalescing black holes (BHs) Abbott  _et al._ (2016). Since
then it has been extensively discussed how we can use this new opportunity to
probe fundamental physics beyond the Standard Model (SM) Perkins  _et al._
(2021). Among many ideas put forward, using the compact binaries to search for
new light degrees of freedom has attracted lots of attentions.

Light bosons are ubiquitous in new physics models, and can be good dark matter
candidates Essig  _et al._ (2013). They are typically very weakly coupled to
SM and have escaped all direct experimental searches so far. However, they
most likely gravitate as known particles do. Strong gravity systems like BHs
are thus natural laboratories to study them. In particular, it is known that a
fast rotating BH of mass $M$ can copiously produce light bosons of mass $\nu$
due to superradiance instability, a pure gravitational effect Zeldovich
(1971); Starobinsky (1973); Brito  _et al._ (2015). When the “gravitational
fine structure constant” $\alpha\equiv GM\nu/\hbar c\ll 1$, the boson $\sigma$
can form nonrelativistic and metastable condensate around the BH with
hydrogen-like spectrum. We can thus call the formed condensate a
“superradiance cloud” (SC) and call the whole system of SC + central BH a
“gravitational atom” (GA). For stellar mass BHs, the corresponding mass of
$\sigma$ is around $\nu\sim 10^{-13}\text{eV}(\alpha/0.01)(10M_{\odot}/M)$.

The phenomenology of SC has long been under active study. See e.g. Arvanitaki
and Dubovsky (2011); Cardoso  _et al._ (2011); Yoshino and Kodama (2012);
Arvanitaki  _et al._ (2015). The LIGO/Virgo discovery of many coalescing BH
binaries added new opportunities Abbott  _et al._ (2019, 2021). Recent studies
have exploited the inspiral phase of a GA-borne binary, showing that a
resonant level transition can occur if the binary’s orbital period matches the
energy split of certain levels of a GA Baumann  _et al._ (2019a, b, 2020);
Ding  _et al._ (2021); Tong  _et al._ (2021). This is dubbed “gravitational
collider physics” Baumann  _et al._ (2020) as the induced resonant transition
is reminiscent of a resonant peak of on-shell particle production at a real
collider.

The level transition requires the match of frequencies of two different
origins, which is expected to be rare in generic inspiral systems.
Furthermore, the GW signals of level transition, while being a distinct
signature, could be challenging to search for in realistic GW experiments.
Therefore it will be useful to look for alternative signals of GAs that are
less restrictive on binary’s orbital parameters and easier to search for.

In this Letter, we point out such a signal for binaries with eccentric orbits.
The key point is that the GA typically possesses huge mass quadrupole Q. In
terms of a dimensionless quadrupole parameter $\kappa\equiv-QM/J^{2}$ (where
$J$ is the spin), a GA can have $\kappa\sim 10^{3}(0.1/\alpha)^{3}$. For
comparison, $\kappa=1$ for Kerr BHs and $\kappa\sim\mathcal{O}({1\text{-}10})$
for a neutron star Laarakkers and Poisson (1999); Pappas and Apostolatos
(2012). As a result, a GA can induce very significant apsidal precession of an
eccentric orbit for a wide range of orbital parameters, which is the main
signal of this Letter. We will show that this signal can be very significant
and potentially observable for BH-GA binaries and pulsar-GA binaries.

Our signal does not rely on any resonant transitions, and therefore can be
viewed as an “off-shell” signal of gravitational collider physics. Thus it
shares both the merit and the drawback of “off-shell” observables. It shows up
for wider range of parameters, but it does not probe the internal structure of
a GA. The signal can nevertheless be a useful probe of GAs, since we are
virtually unaware of any other astrophysical objects sharing similar masses
and mass quadrupoles. Thus the apsidal precession signal opens up a continuous
and large parameter space for detecting GA beyond the known narrow resonance
band, shown in Fig. 1.

LIGO-type binaries are expected to have observably large eccentricity at low
frequencies, as predicted by a wide class of dynamical formation scenarios
Samsing  _et al._ (2014); Samsing (2018); Samsing and D’Orazio (2018); Samsing
_et al._ (2020); Antonini and Perets (2012); Hoang  _et al._ (2018); Fragione
and Bromberg (2019); Randall and Xianyu (2018a, b, 2019a); Silsbee and
Tremaine (2017); Randall and Xianyu (2019b); Deme  _et al._ (2020). Such
eccentric binaries are important targets of future space GW telescopes Ruan
_et al._ (2020); Liu  _et al._ (2020); Toubiana  _et al._ (2020); Buscicchio
_et al._ (2021); Randall  _et al._ (2021). On the other hand, known pulsar
binaries typically possess nonzero eccentricity. No pulsar-BH binaries are
detected yet but they could be found in the future Faucher-Giguere and Loeb
(2011); Shao and Li (2018). Therefore we shall consider our signal in
eccentric pulsar-BH binaries as well.

The effects of the mass quadrupole on binary’s orbit and corresponding GWs
have long been studied Barker and O’Connell (1975); Poisson (1998); Barack and
Cutler (2007). But to our best knowledge, our work is the first to address GA-
induced apsidal precession in eccentric BH-BH and pulsar-BH binaries. See
Krishnendu  _et al._ (2017); Ferreira  _et al._ (2017); Hannuksela  _et al._
(2019); Zhang and Yang (2020); Amorim  _et al._ (2019) for related works.

Superradiance cloud. In this Letter we consider a real scalar field $\sigma$
of mass $\nu$. Generalization to other bosonic species will be presented in a
companion paper Su  _et al._ . The superradiance and the GA formation can be
studied by solving the field equation $(\square-\nu^{2})\sigma=0$ outside a
Kerr BH of mass $M$ and spin $J$. Here and below we take $\hbar=c=1$. When
$\alpha\ll 1$, one finds hydrogen-like bound states $|\mathsf{n\ell
m}\rangle$. At leading order in $\alpha$ the binding energy of $|\mathsf{n\ell
m}\rangle$ is given by $E_{\mathsf{n}}=-\alpha^{2}\nu/(2\mathsf{n}^{2})\ll\nu$
, justifying the nonrelativistic approximation.

One important difference from the hydrogen, though, is the ingoing boundary
condition for $\sigma$ at the BH horizon, which introduces imaginary parts in
the energy eigenvalues and thus can trigger production or depletion of states.
Detailed study shows that $|211\rangle$ and $|322\rangle$ are the two leading
states that a spinning BH can produce which then remain quasi-stable. The
leading depletion channel for them is GW radiation, which is typically very
slow Yoshino and Kodama (2014). Therefore, these states can exist with
astronomical lifetime. Other states might also be relevant but we will mainly
focus on $|211\rangle$ below.

The production and depletion time scales for $|211\rangle$ have been worked
out as (See e.g. Baumann  _et al._ (2020))

$\displaystyle T_{\text{prod}}(|211\rangle)\sim
10^{6}\text{yr}(M/10M_{\odot})(0.016/\alpha)^{9},$ (1) $\displaystyle
T_{\text{depl}}(|211\rangle)\sim
10^{8}\text{yr}(M/10M_{\odot})(0.062/\alpha)^{15}.$ (2)

This puts bounds on $\alpha$ for given masses, shown in Fig. 1. After a period
of time $T_{\text{prod}}$, a SC is formed in $|211\rangle$ state. Its energy
density at leading order in $\alpha$ is time independent,

$\rho_{211}(\mathbf{x})=\frac{1}{64\pi}\alpha^{5}\nu^{5}m_{\text{C}}\,r^{2}e^{-\alpha\nu
r}\sin^{2}\theta,$ (3)

where $m_{\text{C}}$ is the total mass of the SC. From this we can find the
spin $S_{\text{C}}$ and the scalar mass quadrupole
$Q_{\text{C}}\equiv\int{\mathrm{d}}^{3}\mathbf{x}\,\rho(\mathbf{x})x^{2}\text{P}_{2}(\cos\theta)$
as,

$\displaystyle S_{\text{C}}=\alpha m_{\text{C}}r_{\text{B}},$ $\displaystyle
Q_{\text{C}}=-6m_{\text{C}}r_{\text{B}}^{2},$ (4)

where $r_{\text{B}}\equiv 1/(\alpha\nu)$ is the Bohr radius. Assuming the BH
with almost maximal initial spin, the cloud mass can reach
$m_{\text{C}}=\alpha M$ when the production is saturated. Therefore we see
that the cloud carries $\mathcal{O}({\alpha})$ of total mass of the GA but
carries almost all the angular momentum and mass quadrupole of the GA when
$\alpha\ll 1$.

Figure 1: The allowed region (unshaded) for observing $\dot{\chi}_{\text{Q}}$
in $(\alpha,a)$ plane, for a binary of $(40+40)M_{\odot}$, assuming saturated
$|211\rangle$ state. The two shaded regions in gray are excluded by the
condition $a>4r_{\text{B}}$ and $\dot{\chi}_{\text{Q}}>\dot{\chi}_{\text{S}}$,
respectively. The blue belts mark several resonant transition bands Baumann
_et al._ (2020).

GA-induced orbital precessions. Now consider a GA-borne binary with eccentric
orbit. The GA’s oblate energy distribution produces a deviation from Keplerian
potential already at the Newtonian level. Consequently, the orbit will exhibit
both apsidal and nodal precessions which we shall review below.

Let the two binary members have masses $m_{1}$ and $m_{2}$. We shall also use
$m\equiv m_{1}+m_{2}$ and $\mu\equiv m_{1}m_{2}/m$. Let the first mass $m_{1}$
be a GA, consisting of a BH of mass $m_{B}$ and a saturated SC in
$|211\rangle$ of total mass $m_{\text{C}}\simeq\alpha m_{1}$, so that
$m_{1}=m_{\text{B}}+m_{\text{C}}$. We assume that the GA is almost maximally
spinning, with angular momentum $S_{1}=|\mathbf{S}_{1}|\simeq Gm_{1}^{2}$. The
mass quadrupole tensor of GA is given by $Q^{ij}=Q_{1}(\mkern
2.0mu\widehat{\mkern-2.0muS\mkern-2.0mu}\mkern 2.0mu_{1}^{i}\mkern
2.0mu\widehat{\mkern-2.0muS\mkern-2.0mu}\mkern 2.0mu_{1}^{j}-\delta^{ij}/3)$
where $Q_{1}=Q_{\text{C}}$ is the scalar quadrupole given in (4) and $\mkern
2.0mu\widehat{\mkern-2.0muS\mkern-2.0mu}\mkern 2.0mu_{1}^{i}\equiv
S_{1}^{i}/S_{1}$. The second object $m_{2}$ will be taken as a point mass. It
can possess finite spin but we will neglect it in most of discussions.

We focus on the early inspiral period of the binary evolution where the Post-
Newtonian (PN) approximation is easily valid. Our signal, the quadrupole-
induced precession, is a Newtonian effect. But it is also important to include
higher PN effects that contribute orbital precessions, including the GR-
induced apsidal precession at 1PN and the spin-orbit (SO) coupling at 1.5PN.
The spin-spin (SS) coupling at 2PN is always negligible in our system. The GR
precession can in principle be subtracted once we know the masses and orbital
parameters, while SO- and SS-induced precessions depend on more free
parameters, and thus should be regarded as contaminations of our signal.

The GA is also subject to tidal deformation and this can back react to the
orbit and affects GW signals Cardoso  _et al._ (2020); Takahashi and Tanaka
(2021); De Luca and Pani (2021). But this belongs to higher PN effects and we
will neglect them in this Letter. The validity of treating our system as a
well separated binary also puts a constraint on the orbital separation $a$.
Following Baumann  _et al._ (2019a) we take $a>4r_{\text{B}}$ for
$|211\rangle$ state, also shown in Fig. 1. We don’t include mass ratio
dependence in this constraint because we will never consider the case of
$m_{2}\gg m_{1}$ in this Letter.

With above considerations in mind, we can model the GA as a rigid symmetric
top, described by its mass $m_{1}$, spin $\mathbf{S}_{1}$, and quadrupole
$Q_{1}$. We will adopt the famous EIH Lagrangian at 1PN to describe the
barycenter motion of the two masses Einstein  _et al._ (1938), adding
appropriate terms to include rotational degrees of the GA. In particular, the
interaction potential $V=V_{\text{O}}+V_{\text{Q}}+V_{\text{SO}}$ contains the
point-mass term, the quadrupole term, and SO term, respectively,

$\displaystyle V_{\text{O}}=$
$\displaystyle-\displaystyle\frac{\,{Gm\mu}\,}{r}+\text{1PN terms},$ (5)
$\displaystyle V_{\text{Q}}=$
$\displaystyle~{}\displaystyle\frac{\,{Gm\mu}\,}{2r^{3}}q_{1}(1-3\cos^{2}n_{1}),$
(6) $\displaystyle V_{\text{SO}}=$
$\displaystyle~{}\displaystyle\frac{\,{2GL}\,}{r^{3}}(m_{1}+\mbox{$\frac{\,{3}\,}{4}$}m_{2})s_{1}\cos
n_{1},$ (7)

where $n_{1}$ is the angle between $\mathbf{S}_{1}$ and the orbital angular
momentum $\mathbf{L}$. We further define $q_{1}\equiv Q_{1}/m_{1}$,
$s_{1}\equiv S_{1}/m_{1}$, and $L\equiv|\mathbf{L}|$.

The computation of precession rates of various orbital elements from the above
potential is straightforward, of which we will review systematically in a
companion paper Su  _et al._ . As mentioned, there are both nodal precession
(precession of $L$) and the apsidal precession (precession of the periapsis
_within_ the orbital plane). We discuss them in turn.

First the nodal precession. Our system always has $L\gg S_{1}$, and so the
angular momentum conservation dictates that $\mathbf{L}$ precesses only within
a small cone of angular size $\sim S_{1}/L$. Therefore, the nodal precession
only leaves small modulation of GW amplitudes (and phases) in generic
situation. The only exception is when the line of sight overlaps with the
small precession cone. In this case the nodal precession mixes with the
orbital motion. But this is a rare situation. So the nodal precession is
unimportant for most GW sources. For pulsar-BH binaries, however, the nodal
precession can still be quite significant due to the high precision of pulsar
timing measurements.

The most important effect is the apsidal precession. Let $\chi$ be the
argument of periapsis in the orbital plane. As mentioned, the apsidal
precession $\dot{\chi}$ has three contributions, from the GR precession at
1PN, mass quadrupole coupling $V_{\text{Q}}$, and the SO coupling
$V_{\text{SO}}$, re-spectively:

$\displaystyle\dot{\chi}=$
$\displaystyle~{}\dot{\chi}_{\text{GR}}+\dot{\chi}_{\text{Q}}+\dot{\chi}_{\text{S}},$
(8) $\displaystyle\dot{\chi}_{\text{GR}}=$
$\displaystyle~{}\displaystyle\frac{\,{3Gm\omega}\,}{a(1-e^{2})},$ (9)
$\displaystyle\dot{\chi}_{\text{Q}}=$
$\displaystyle~{}\displaystyle\frac{\,{3\omega}\,}{4a^{2}(1-e^{2})^{2}}q_{1}(1-3\cos^{2}n_{1}),$
(10) $\displaystyle\dot{\chi}_{\text{S}}=$
$\displaystyle~{}\displaystyle\frac{\,{-3\omega^{2}}\,}{(1-e^{2})^{3/2}}s_{1}\Big{(}1+\displaystyle\frac{\,{m_{1}}\,}{3m}\Big{)}\cos
n_{1},$ (11)

where $a$ and $e$ are semi-major axis and the eccentricity of the orbit, and
$\omega=\sqrt{Gm/a^{3}}$ is the orbital frequency. Let us compare their sizes.
For maximally spinning and saturated GA in $|211\rangle$, we have $s_{1}\simeq
Gm_{1}$ and $q_{1}\simeq-6\alpha^{3}(Gm_{1})^{2}$. Then,

$\displaystyle\displaystyle\frac{\,{\dot{\chi}_{\text{Q}}}\,}{\dot{\chi}_{\text{GR}}}\simeq\displaystyle\frac{\,{3u_{1}}\,}{2\alpha^{3}}\displaystyle\frac{\,{Gm_{1}}\,}{a(1-e^{2})},$
(12)
$\displaystyle\displaystyle\frac{\,{\dot{\chi}_{\text{Q}}}\,}{\dot{\chi}_{\text{S}}}\simeq\displaystyle\frac{\,{3u_{1}}\,}{2\alpha^{3}(1+u_{1}/3)}\Big{[}\displaystyle\frac{\,{Gm_{1}}\,}{a(1-e^{2})}\Big{]}^{1/2},$
(13)

where $u_{1}\equiv m_{1}/m$. The factor $\alpha^{-3}$ from the large mass
quadrupole is an enhancement while the factor $Gm_{1}/a<\alpha^{2}/4$ is
always a suppression after imposing the perturbative condition
$a>4r_{\text{B}}$. (One can impose a stronger condition $a(1-e)>4r_{\text{B}}$
but we will neglect the $(1-e)$ factor since we do not consider very large
$e\simeq 1$.)

The GR precession can in principle be subtracted from the total $\dot{\chi}$
up to measurement errors once we know $a$, $e$, and $m$. On the other hand,
$\dot{\chi}_{\text{S}}$ contaminates our signal by introducing new free
parameter $s_{1}$. There is a similar contamination from the spin of $m_{2}$
which we assume to be no larger than that of $m_{1}$ and thus is suppressed
here. We impose a constraint on a by requiring
$|\dot{\chi}_{\text{Q}}|\gg|\dot{\chi}_{\text{S}}|$, to stay away from the SO-
induced contamination, shown in Fig. 1. The SS-induced precession is always
smaller than SO-precession and thus can be safely neglected. Then, in the
range of $|\dot{\chi}_{\text{Q}}|\gg|\dot{\chi}_{\text{S}}|$ and after
subtracting the GR precession, the quadrupole term dominates the apsidal
precession, leaving a distinct signal. Interestingly, the GR precession is
always prograde while the quadrupole precession can be either prograde or
retrograde.

The spin $\mathbf{S}_{1}$ also precesses due to the torque from quadrupole and
SO couplings, and one may hope to probe this from the time variation of $\cos
n_{1}=\mkern 2.0mu\widehat{\mkern-2.0mu\mathbf{L}\mkern-2.0mu}\mkern
2.0mu\cdot\mkern 2.0mu\widehat{\mkern-2.0mu\mathbf{S}\mkern-2.0mu}\mkern
2.0mu_{1}$. However, since $L$, $S$, and $\mathbf{L}+\mathbf{S}$ are all
conserved in our system, it follows that $n_{1}$ is also a constant in time.
Interestingly, $\dot{\chi}_{\text{Q}}$ remains constant even when $m_{2}$
carries nonzero quadrupole in which case $n_{1}$ is no longer a constant. We
will elaborate more on this in a companion paper Su  _et al._ .

In summary, by staying in the range of
$|\dot{\chi}_{\text{Q}}|\gg|\dot{\chi}_{\text{S}}|$, only
$\dot{\chi}_{\text{GR}}$ and $\dot{\chi}_{\text{Q}}$ contribute the apsidal
precession, and our signal, the quadrupole induced apsidal precession, depends
on a single constant parameter which we define as $q_{\text{eff}}$ for later
convenience:

$q_{\text{eff}}\equiv q_{1}(1-3\cos^{2}n_{1}).$ (14)

GWs from GA-BH Binaries. The GWs radiated from a binary essentially trace the
binary’s orbital motion. Thus one could hope to see the orbital precession
from the GW waveform. But this is challenging since we need to see the GW from
early inspiral phase rather than close to the merger. Early inspiral is
required both to have nonzero eccentricity and also to satisfy the
perturbative condition (i.e., the second body be well outside the cloud).
Requiring $a>4r_{\text{B}}$, the GW frequency (of the 2nd harmonic)
$f_{2}=\pi^{-1}\sqrt{Gm/a^{3}}<\pi^{-1}\sqrt{Gm/(4r_{\text{B}})^{3}}\sim
0.14\text{Hz}(\alpha/0.05)^{3}(10M_{\odot}/m_{1})$, assuming $m_{1}=m_{2}$.
Thus we see that for $10M_{\odot}$ GA the signal is already in the window of
LISA. The similar signal from massive BH binaries ($m>10^{4}M_{\odot}$) or
extreme-mass-ratio inspirals (typically a supermassive BH with a stellar-mass
BH) will be at even lower frequencies, challenging to resolve individually
even by space GW telescope. Therefore we will focus on stellar-mass binaries
in this section and leave more general possibilities to a future study.

The GW amplitudes from the binary in the transverse-traceless (TT) gauge is
computed by $h^{ij}=(2G/d)\Lambda^{ij,k\ell}\ddot{M}^{k\ell}$ where
$\Lambda^{ij,k\ell}$ is the projector to TT components, $M^{ij}=\mu
r^{i}r^{j}$ is the second mass moment of the binary, and $d$ is the source
distance. In our case $M^{ij}$ depends on $a$, $e$, the true anomaly
$\psi(t)$, and the three Euler angles describing the orbital orientation.
These angles can be chosen to be the argument of periapsis $\chi$, the
inclination $I$ and the longitude of ascending node $\Omega$, defined with
respect to a reference plane (e.g., the plane perpendicular to the line of
sight). In a GA-borne binary, all three angles could precess and this
precession could enter the GW signal. But as discussed above, in generic
situations, the nodal precession is tiny, and only the apsidal precession
$\dot{\chi}$ could potentially be observable.

To identify the signature of $\dot{\chi}$ in GW signals, we go to the orbital
plane, in which the position vector $r_{i}=r(\cos\mkern
2.0mu\overline{\mkern-2.0mu\psi\mkern-2.0mu}\mkern 2.0mu,\sin\mkern
2.0mu\overline{\mkern-2.0mu\psi\mkern-2.0mu}\mkern 2.0mu)$, with
$r=a(1-e^{2})/(1+e\cos\psi)$, and $\mkern
2.0mu\overline{\mkern-2.0mu\psi\mkern-2.0mu}\mkern 2.0mu=\psi+\chi$. (For low
inclination system $I\simeq 0$, one should also include the motion of
ascending node so that $\mkern
2.0mu\overline{\mkern-2.0mu\psi\mkern-2.0mu}\mkern 2.0mu=\psi+\chi+\Omega$.
However, once again, we consider such systems rare.) The important point is
that the radial motion $r(\psi)$ depends only on the true anomaly which is a
periodic function of time with frequency $\psi$, while the angular motion also
includes the apsidal precession, so that the angle $\mkern
2.0mu\overline{\mkern-2.0mu\psi\mkern-2.0mu}\mkern 2.0mu$ is a periodic
function with frequency $\omega+\dot{\chi}$.

Figure 2: The GW power spectrum of an eccentric binary with apsidal
precession. The split in each triplet is exaggerated.

To understand intuitively how to observe $\dot{\chi}$ through GWs, we note
that $\psi(t)$ executes non-uniform motion for eccentric orbit (Kepler’s 2nd
law), resulting overtones in GW spectrum with frequency $\omega_{n}=n\omega$
$(n=1,2,\cdots)$. We can assume $\dot{\chi}$ to be almost constant up to
slight orbital decay. Then, each of $M^{ij}$’s harmonic components splits into
a triplet with frequencies
$(\omega_{n},\omega_{n\pm})\equiv(n\omega,n\omega\pm 2\dot{\chi})$. It turns
out that the $\omega_{n+}$ component almost always dominates the GW power in a
triplet, as shown in Fig. 2. Therefore we can retain the $\omega_{n+}$
component only. Then the effect of apsidal precession enters through the
_anharmonic overtones_ , namely, $\omega_{n+}/\omega_{m+}$ is not integer for
any $n\neq m$. Of course the visibility of this non-integer ratio depends on
the frequency resolution, which requires $\dot{\chi}_{\text{Q}}>T_{O}^{-1}$
with $T_{O}$ the total observation time. Using (10) this translates very
roughly to $q_{1}/a^{2}>(\omega T_{O})^{-1}$. Further using $q_{1}\sim\alpha
r_{\text{B}}^{2}$ and the perturbativity condition $a\gtrsim r_{\text{B}}$, we
have $q_{1}/a^{2}<\alpha$. On the other hand, taking LISA as an example,
$\omega\sim$0.01Hz and $T_{O}\sim$ yr, so $\omega T_{O}\sim 10^{5}$. So there
is a finite range $10^{-5}<q_{1}/a^{2}<\alpha$ in which we can hope to see
$\dot{\chi}_{\text{Q}}$ through GWs.

We note that this effect is absent for circular orbit where $r(\psi)$ is a
constant, and $\mkern 2.0mu\overline{\mkern-2.0mu\psi\mkern-2.0mu}\mkern
2.0mu=(\omega+\dot{\chi})t$ is essentially a uniform function of time. So a
nonzero $\dot{\chi}$ amounts to a constant shift in $\omega$ and thus is
degenerate with other orbital parameters. Of course this degeneracy is broken
by PN effects including chirping. But we expect less significant effect from
the quadrupole in this case. Below we will show that increasing $e$ reduces
the error in $q_{\text{eff}}$.

Figure 3: The $1\sigma$ error contours of $q_{\text{eff}}$ with the total mass
$m$, the GW frequency of 2nd harmonic $f_{20}$, and the eccentricity $e_{0}$.
We choose N2A5 configuration of LISA and 4yr of total observation time. The
choice of parameters for this plot is a bit beyond the bound $a>4r_{\text{B}}$
so the result is only indicative.

To assess more quantitatively the observability of the quadrupole-induced
apsidal precession through GW, we perform a simplified Fisher analysis for a
binary of $(40+40)M_{\odot}$. A more complete analysis will be presented in Su
_et al._ . Here we consider a 4-parameter set that directly appears in
$\dot{\chi}$, including $q_{\text{eff}}$, the total mass $m$, the GW frequency
$f_{20}$ of 2nd harmonic and eccentricity $e_{0}$ at the starting time of
observation. This parameter set is useful to address potential degeneracy
between $\dot{\chi}_{\text{GR}}$ and $\dot{\chi}_{\text{Q}}$. We apply a time-
domain formula for Fisher matrix Barack and Cutler (2004) and assume N2A5
configuration of LISA noise curve Klein  _et al._ (2016) with 4yr of total
observation time. For our purpose it is enough to use the Newtonian orbit
augmented by apsidal precession (8). We apply Peters’s equation Peters (1964)
to account for orbital decay but treat Euler angles $(I,\Omega)$ as constant.
In Fig. 3 we show the $1\sigma$ contours for $e_{0}=(0.1,0.3,0.5)$, with all
other parameters indicated in the Figure. The signal-to-noise ratios (SNR) for
three choices of e are normalized to that of $e=0.5$ at $d=200$Mpc, which is
roughly SNR $\simeq 42$. One can see that the mass quadrupole $q_{\text{eff}}$
can be clearly identified for $e=0.3$ and well measured for $e=0.5$, while
large degeneracy and errors appear for $e=0.1$, just as expected.

Cloud-induced precession in pulsar-GA binaries. The effect of the
superradiance cloud can also be searched for in pulsar-BH binaries. Compared
to GWs, the very precise pulsar timing data usually allow more direct
measurement of orbital elements. The apsidal precession rate $\dot{\chi}$
becomes almost a direct observable readily to be extracted from the time delay
data.

In pulsar binaries the 5 “Keplerian” parameters describing the geometry and
kinetics of the orbit can often be very precisely measured, while in well
situated systems one can also measure several post-Keplerian (PK) parameters
with good precision Stairs (2003). In GR, all PK parameters can be calculated
from Keplerian parameters and the two masses. Therefore measuring two PK
parameters can determine the two masses, conventionally denoted by $m_{p}$ for
the pulsar and $m_{c}$ for the companion. (With our notation, $m_{1}=m_{c}$
and $m_{2}=m_{p}$.) Measuring more than 2 PK parameters then serves as
consistency check. It is often the case that the apsidal precession rate
$\dot{\chi}$, the Einstein parameter $\gamma$, and the orbital decay
$\dot{P}_{b}$ can be measured to high precision. For highly inclined system
the Shapiro delay parameters $r$ and $s$ can also be measured.

Figure 4: (a) The comparison between the cloud-induced apsidal precession rate
$\dot{\chi}_{\text{Q}}$ and the GR induced precession rate
$\dot{\chi}_{\text{GR}}$ for a range of $\alpha$ and orbital period $P_{b}$.
In this plot we take $m_{p}=1.4M_{\odot}$, $m_{c}=10M_{\odot}$, $e=0.6$, and
assume the cloud mass is saturated. (b) The PK parameters fit in
$(m_{p},m_{c})$ plane for a pulsar-GA binary marked by the yellow dot in (a).
The error bands, shown in parentheses, are taken with the same order of
magnitude from known pulsar binaries such as PSR B1913+16 and PSR B1534+12
Stairs (2003). The solid and dashed magenta curves show the $\dot{\chi}$
without and with subtracting the quadrupole contribution, respectively.

In our case, however, the orbital evolution depends on further parameters
including $q_{1}$ and $n_{1}$ which correct PK parameters. Most of them are
small corrections of order $q_{1}/a^{2}\ll 1$. But there are two post-
Keplerian parameters whose quadrupole correction is not suppressed by
$q_{1}/a^{2}$. The most important one is of course $\dot{\chi}$ given in (8)
(normally denoted by $\omega$ in pulsar binary literature). In Fig. 4(a) we
show the comparison between $\dot{\chi}_{\text{Q}}$ and
$\dot{\chi}_{\text{GR}}$ for a range of orbital period $P_{b}$ and $\alpha$.
We assume $m_{c}$ is a GA saturated in either $|211\rangle$ or $|322\rangle$.
We see that the correction from $\dot{\chi}_{\text{Q}}$ is actually huge
compared with measurement precision, which is easily better than $10^{-3}$. In
Fig. 4(a) we also put a yellow spot and show an imagined PK parameter fit in
the $(m_{p},m_{c})$ plane for this data point as is usually done for pulsar
binaries. Consistency requires that all curves meet at one point up to errors.
However, for the parameter we are choosing
($\dot{\chi}_{\text{Q}}/\dot{\chi}_{\text{GR}}\simeq 0.1$), the failure of the
consistency check is already dramatic, showing that there can be a good chance
to measure $\dot{\chi}_{\text{Q}}$ in this way.

The other PK parameter receiving significant correction is $\dot{x}$ where
$x\equiv a\sin I$. In ordinary system $\dot{x}$ comes mainly from the GW-
induced orbital decay $\dot{a}$. This is a tiny effect and $\dot{x}$ has not
yet been measured in known systems. In our system, however, the inclination
$\sin I$ becomes time dependent due to quadrupole coupling (6). One can show
that $\dot{x}$ has the same time scale as $\dot{\chi}_{\text{Q}}$ and so is
easily within the observation range. The total variation of $x$ remains small
as mentioned before, $\Delta x/x\sim\mathcal{O}({S/L})$, due to the angular
momentum conservation. But for the system plotted in Fig. 4(b), $S/L\simeq
0.06$ and so can easily be larger than observation precision of $x$, making
$\dot{x}$ potentially observable. We will present a more systematic study of
this effect in a future work.

Discussions. Fast rotating BHs are natural factories and habitat for
ultralight bosons. The resulting GAs are ideal labs for probing these
otherwise invisible particles. In this Letter we showed that the large mass
quadrupole of a GA can excite observable apsidal precession for eccentric
binaries. The effect does not require resonant transitions and thus extends
the reach of gravitational collider physics. It will be interesting to
consider more such non-resonant effects in similar systems, such as a binary
of two GAs, and also a “gravitational molecule” with one SC surrounding two
BHs Ikeda  _et al._ (2021); Liu and Lyu (2021). We leave these directions for
future works.

###### Acknowledgements.

We thank Jiandong Zhang for helpful discussions. We also thank Horng Sheng
Chia, Jing Ren, and Xi Tong for useful comments on the manuscript. ZZX is
supported by Tsinghua University Initiative Scientific Research Program.

## References

  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Perkins  _et al._ (2021) S. E. Perkins, N. Yunes, and E. Berti, Phys. Rev. D 103, 044024 (2021), arXiv:2010.09010 [gr-qc] .
  * Essig  _et al._ (2013) R. Essig _et al._ , in  _Community Summer Study 2013: Snowmass on the Mississippi_ (2013) arXiv:1311.0029 [hep-ph] .
  * Zeldovich (1971) Y. B. Zeldovich, JETP Lett. 14, 180 (1971).
  * Starobinsky (1973) A. A. Starobinsky, Sov. Phys. JETP 37, 28 (1973).
  * Brito  _et al._ (2015) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
  * Arvanitaki and Dubovsky (2011) A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011), arXiv:1004.3558 [hep-th] .
  * Cardoso  _et al._ (2011) V. Cardoso, S. Chakrabarti, P. Pani, E. Berti, and L. Gualtieri, Phys. Rev. Lett. 107, 241101 (2011), arXiv:1109.6021 [gr-qc] .
  * Yoshino and Kodama (2012) H. Yoshino and H. Kodama, Prog. Theor. Phys. 128, 153 (2012), arXiv:1203.5070 [gr-qc] .
  * Arvanitaki  _et al._ (2015) A. Arvanitaki, M. Baryakhtar, and X. Huang, Phys. Rev. D 91, 084011 (2015), arXiv:1411.2263 [hep-ph] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021), arXiv:2010.14527 [gr-qc] .
  * Baumann  _et al._ (2019a) D. Baumann, H. S. Chia, and R. A. Porto, Phys. Rev. D 99, 044001 (2019a), arXiv:1804.03208 [gr-qc] .
  * Baumann  _et al._ (2019b) D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, JCAP 12, 006 (2019b), arXiv:1908.10370 [gr-qc] .
  * Baumann  _et al._ (2020) D. Baumann, H. S. Chia, R. A. Porto, and J. Stout, Phys. Rev. D 101, 083019 (2020), arXiv:1912.04932 [gr-qc] .
  * Ding  _et al._ (2021) Q. Ding, X. Tong, and Y. Wang, Astrophys. J. 908, 78 (2021), arXiv:2009.11106 [astro-ph.HE] .
  * Tong  _et al._ (2021) X. Tong, Y. Wang, and H.-Y. Zhu, (2021), arXiv:2106.13484 [astro-ph.HE] .
  * Laarakkers and Poisson (1999) W. G. Laarakkers and E. Poisson, Astrophys. J. 512, 282 (1999), arXiv:gr-qc/9709033 .
  * Pappas and Apostolatos (2012) G. Pappas and T. A. Apostolatos, Phys. Rev. Lett. 108, 231104 (2012), arXiv:1201.6067 [gr-qc] .
  * Samsing  _et al._ (2014) J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, Astrophys. J. 784, 71 (2014), arXiv:1308.2964 [astro-ph.HE] .
  * Samsing (2018) J. Samsing, Phys. Rev. D 97, 103014 (2018), arXiv:1711.07452 [astro-ph.HE] .
  * Samsing and D’Orazio (2018) J. Samsing and D. J. D’Orazio, Mon. Not. Roy. Astron. Soc. 481, 5445 (2018), arXiv:1804.06519 [astro-ph.HE] .
  * Samsing  _et al._ (2020) J. Samsing, D. J. D’Orazio, K. Kremer, C. L. Rodriguez, and A. Askar, Phys. Rev. D 101, 123010 (2020), arXiv:1907.11231 [astro-ph.HE] .
  * Antonini and Perets (2012) F. Antonini and H. B. Perets, Astrophys. J. 757, 27 (2012), arXiv:1203.2938 [astro-ph.GA] .
  * Hoang  _et al._ (2018) B.-M. Hoang, S. Naoz, B. Kocsis, F. A. Rasio, and F. Dosopoulou, Astrophys. J. 856, 140 (2018), arXiv:1706.09896 [astro-ph.HE] .
  * Fragione and Bromberg (2019) G. Fragione and O. Bromberg, Mon. Not. Roy. Astron. Soc. 488, 4370 (2019), arXiv:1903.09659 [astro-ph.GA] .
  * Randall and Xianyu (2018a) L. Randall and Z.-Z. Xianyu, Astrophys. J. 853, 93 (2018a), arXiv:1708.08569 [gr-qc] .
  * Randall and Xianyu (2018b) L. Randall and Z.-Z. Xianyu, Astrophys. J. 864, 134 (2018b), arXiv:1802.05718 [gr-qc] .
  * Randall and Xianyu (2019a) L. Randall and Z.-Z. Xianyu, Astrophys. J. 878, 75 (2019a), arXiv:1805.05335 [gr-qc] .
  * Silsbee and Tremaine (2017) K. Silsbee and S. Tremaine, Astrophys. J. 836, 39 (2017), arXiv:1608.07642 [astro-ph.HE] .
  * Randall and Xianyu (2019b) L. Randall and Z.-Z. Xianyu, (2019b), arXiv:1902.08604 [astro-ph.HE] .
  * Deme  _et al._ (2020) B. Deme, B.-M. Hoang, S. Naoz, and B. Kocsis, Astrophys. J. 901, 125 (2020), arXiv:2005.03677 [astro-ph.HE] .
  * Ruan  _et al._ (2020) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Int. J. Mod. Phys. A 35, 2050075 (2020), arXiv:1807.09495 [gr-qc] .
  * Liu  _et al._ (2020) S. Liu, Y.-M. Hu, J.-d. Zhang, and J. Mei, Phys. Rev. D 101, 103027 (2020), arXiv:2004.14242 [astro-ph.HE] .
  * Toubiana  _et al._ (2020) A. Toubiana, S. Marsat, S. Babak, J. Baker, and T. Dal Canton, Phys. Rev. D 102, 124037 (2020), arXiv:2007.08544 [gr-qc] .
  * Buscicchio  _et al._ (2021) R. Buscicchio, A. Klein, E. Roebber, C. J. Moore, D. Gerosa, E. Finch, and A. Vecchio, (2021), arXiv:2106.05259 [astro-ph.HE] .
  * Randall  _et al._ (2021) L. Randall, A. Shelest, and Z.-Z. Xianyu, (2021), arXiv:2103.16030 [astro-ph.HE] .
  * Faucher-Giguere and Loeb (2011) C. A. Faucher-Giguere and A. Loeb, Mon. Not. Roy. Astron. Soc. 415, 3951 (2011), arXiv:1012.0573 [astro-ph.HE] .
  * Shao and Li (2018) Y. Shao and X.-D. Li, Mon. Not. Roy. Astron. Soc. 477, L128 (2018), arXiv:1804.06014 [astro-ph.HE] .
  * Barker and O’Connell (1975) B. M. Barker and R. F. O’Connell, Phys. Rev. D 12, 329 (1975).
  * Poisson (1998) E. Poisson, Phys. Rev. D 57, 5287 (1998), arXiv:gr-qc/9709032 .
  * Barack and Cutler (2007) L. Barack and C. Cutler, Phys. Rev. D 75, 042003 (2007), arXiv:gr-qc/0612029 .
  * Krishnendu  _et al._ (2017) N. V. Krishnendu, K. G. Arun, and C. K. Mishra, Phys. Rev. Lett. 119, 091101 (2017), arXiv:1701.06318 [gr-qc] .
  * Ferreira  _et al._ (2017) M. C. Ferreira, C. F. B. Macedo, and V. Cardoso, Phys. Rev. D 96, 083017 (2017), arXiv:1710.00830 [gr-qc] .
  * Hannuksela  _et al._ (2019) O. A. Hannuksela, K. W. K. Wong, R. Brito, E. Berti, and T. G. F. Li, Nature Astron. 3, 447 (2019), arXiv:1804.09659 [astro-ph.HE] .
  * Zhang and Yang (2020) J. Zhang and H. Yang, Phys. Rev. D 101, 043020 (2020), arXiv:1907.13582 [gr-qc] .
  * Amorim  _et al._ (2019) A. Amorim _et al._ (GRAVITY), Mon. Not. Roy. Astron. Soc. 489, 4606 (2019), arXiv:1908.06681 [astro-ph.GA] .
  * (48) B. Su, Z.-Z. Xianyu, and X. Zhang, “to appear,” .
  * Yoshino and Kodama (2014) H. Yoshino and H. Kodama, PTEP 2014, 043E02 (2014), arXiv:1312.2326 [gr-qc] .
  * Cardoso  _et al._ (2020) V. Cardoso, F. Duque, and T. Ikeda, Phys. Rev. D 101, 064054 (2020), arXiv:2001.01729 [gr-qc] .
  * Takahashi and Tanaka (2021) T. Takahashi and T. Tanaka, (2021), arXiv:2106.08836 [gr-qc] .
  * De Luca and Pani (2021) V. De Luca and P. Pani, (2021), arXiv:2106.14428 [gr-qc] .
  * Einstein  _et al._ (1938) A. Einstein, L. Infeld, and B. Hoffmann, Annals Math. 39, 65 (1938).
  * Barack and Cutler (2004) L. Barack and C. Cutler, Phys. Rev. D69, 082005 (2004), arXiv:gr-qc/0310125 [gr-qc] .
  * Klein  _et al._ (2016) A. Klein _et al._ , Phys. Rev. D93, 024003 (2016), arXiv:1511.05581 [gr-qc] .
  * Peters (1964) P. C. Peters, Phys. Rev. 136, B1224 (1964).
  * Stairs (2003) I. H. Stairs, Living Rev. Rel. 6, 5 (2003), arXiv:astro-ph/0307536 .
  * Ikeda  _et al._ (2021) T. Ikeda, L. Bernard, V. Cardoso, and M. Zilhão, Phys. Rev. D 103, 024020 (2021), arXiv:2010.00008 [gr-qc] .
  * Liu and Lyu (2021) T. Liu and K.-F. Lyu, (2021), arXiv:2107.09971 [astro-ph.HE] .


# How Can We Describe Density Evolution Under Delayed Dynamics?

Michael C. Mackey<EMAIL_ADDRESS>https://www.mcgill.ca/mathematical-physiology-lab/ Departments of Physiology,
Physics & Mathematics, McGill University, Montreal, Quebec, CANADA, H4X 2C1
Marta Tyran-Kamińska<EMAIL_ADDRESS>http://www.math.us.edu.pl/mtyran/
Institute of Mathematics, University of Silesia, Bankowa 14, 40-007 Katowice,
Poland

###### Abstract

Although the theory of density evolution in maps and ordinary differential
equations is well developed, the situation is far from satisfactory in
continuous time systems with delay. This paper reviews some of the work that
has been done numerically and the interesting dynamics that have emerged, and
the largely unsuccessful attempts that have been made to analytically treat
the evolution of densities in differential delay equations. We also present a
new approach to the problem and illustrate it with a simple example.

###### pacs:

02.30.Ks,02.30.Oz,02.50.Ey,05.40.-a,05.40.Jc

> In this paper we highlight an open problem in mathematics that has
> implications for any system whose dynamics are dependent on behavior in the
> past. Namely how can we describe the evolution of densities in such systems.
> We review what is known about the evolutions of densities in discrete time
> maps as well as in systems with dynamics described by ordinary differential
> equations or stochastic differential equations and then highlight the rather
> formidable mathematical problems that arise when one wishes to consider
> delayed, or hereditary, dynamics.

##  I Introduction

We are accustomed to thinking about the trajectories of dynamical systems and
their possible bifurcations as a parameter is varied. Typically one encounters
bifurcation sequences like stable steady state $\to$ simple limit cycle $\to$
complicated limit cycle $\to$ ‘chaotic’ solutions, but the definition of what
constitutes chaos is tricky Hunt and Ott (2015).

Here, we want to turn this around and think about the evolution of densities.
This is akin to the Gibbs’ notion of looking at an ensemble of dynamical
systems, and this ensemble is described by the corresponding density of
states. This just means that we are thinking about looking at a very large
number of copies of a dynamical system, under the assumption that each copy is
not interacting with any others.

##  II Density evolution in dynamical systems

Though the idea of density evolution may seem an unfamiliar one initially, in
point of fact many readers will find that they are really quite familiar with
it from other contexts.

From a purely formal standpoint we start with the definition of the Frobenius-
Perron (FP) operator $P^{t}:L^{1}\to L^{1}$

$\int_{A}P^{t}f(x)\,m(dx)=\int_{S_{t}^{-1}(A)}f(x)\,m(dx)$ (1)

which maps densities to densities. From a technical standpointLasota and
Mackey (1994), $(X,{\cal A},m)$ is a $\sigma$-finite measure space, and
$S_{t}:X\to X$ a measurable nonsingular transformation, _i.e._
$S_{t}^{-1}(A)\in{\cal A}$ for all $A\in{\cal A}$ and $m(S_{t}^{-1}(A))=0$
whenever $m(A)=0$.

This still may look rather unfamiliar, but some examples will smooth the way.
First of all note that if $A=[a,x]$ the Frobenius-Perron operator becomes

$\int_{a}^{x}P^{t}f(s)\,ds=\int_{S_{t}^{-1}([a,x])}f(s)\,ds,$

so

$P^{t}f(x)=\dfrac{d}{dx}\int_{S_{t}^{-1}([a,x])}f(s)\,ds.$ (2)

For example, with the tent (hat) map

$S(x)=\left\\{\begin{array}[]{ll}ax&\qquad\mbox{for}\quad
x\in\left[0,\frac{1}{2}\right)\\\ a(1-x)&\qquad\mbox{for}\quad
x\in\left[\frac{1}{2},1\right],\end{array}\right.$ (3)

$S_{n}$ is the $n$th iterate of $S$, $n\in\mathbb{N}$, and the corresponding
Frobenius-Perron operator is the $n$th iterate $P^{n}$ of the operator

$Pf(x)=\dfrac{1}{a}\left[f\left(\dfrac{x}{a}\right)+f\left(1-\dfrac{x}{a}\right)\right].$

In a more familiar vein if we have a system of ordinary differential equations

$\dfrac{dx_{i}}{dt}={\cal F}_{i}(x),\qquad i=1,\ldots,d,$

then from the definition of the Frobenius-Perron operator we can derive the
evolution equation for $f(x,t)=P^{t}f(x)$:

$\dfrac{\partial f}{\partial t}=-\sum_{i=1}^{d}\dfrac{\partial(f{\cal
F}_{i})}{\partial x_{i}},$ (4)

which is just the generalized Liouville equationLasota and Mackey (1994).

Finally if we have a stochastic differential equation

${dx}={\cal F}(x)dt+\sigma(x)dw(t)$

where $w$ is a Wiener process, then the evolution equation for the density
$f(x,t)\equiv P^{t}f_{0}(x)$ is the Fokker-Planck equationGardiner (1983)

$\frac{\partial f}{\partial t}=-\sum_{=1}^{d}\frac{\partial(f{\cal
F}_{i})}{\partial
x_{i}}+\frac{1}{2}\sum_{i,j=1}^{d}\frac{\partial^{2}(a_{ij}f)}{\partial
x_{i}\partial x_{j}}$ (5)

where $a_{ij}(x)=\sum_{k=1}^{d}\sigma_{ik}(x)\sigma_{jk}(x)$.

However, if we have a variable $x$ evolving under the action of some dynamics
described by a differential delay equation

$x^{\prime}(t)=\epsilon^{-1}{\cal F}(x(t),x(t-\tau)),\quad x(t)=\phi(t)\quad
t\in[-\tau,0],$ (6)

then things are not so clear. We would like to know how some initial density
of the variable $x$ will evolve in time.

Denote the ‘density’ of $x$ on the interval $[t-\tau,t]$ by
$\rho(x,t,t^{\prime}),\forall t^{\prime}\in[t-\tau,t]$. Then we would like to
be able to determine an evolution operator ${\cal U}$ such that the equation

${\cal U}^{t}\rho_{0}(x)=0,\quad\rho_{0}(x)=\rho(x,0,t^{\prime})$ (7)

describes the evolution of $\rho(x,t,t^{\prime})$ given a density of initial
functions $\rho_{0}(x,0,t^{\prime})$. Unfortunately, we don’t really know how
to do this, and that’s the whole point of this paper.

The reason that the problem is so difficult is embodied in (6) and the
infinite dimensional nature of the problem because of the necessity of
specifying the initial function $\phi(t)$ for $t\in[-\tau,0]$, and even
further the ‘density’ of initial functions $\rho_{0}(x,0,t^{\prime})$ for
$t^{\prime}\in[-\tau,0]$.

We know about what ${\cal U}^{t}$ should look like in various limiting cases.
For example, if ${\cal F}(x(t),x(t-\tau))=-x(t)+S(x(t-\tau))$ so (6) becomes

$\epsilon x^{\prime}(t)=-x(t)+S(x(t-\tau)),\quad x(t)=\phi(t)\quad
t\in[-\tau,0],$ (8)

then we expect that:

  1. 1.

If we let $\epsilon\to 0$ and restrict consideration to $t/\tau\in\mathbb{N}$,
then ${\cal U}^{t}$ should reduce to the Frobenius-Perron operator (2) for the
map $S$.

  2. 2.

If $\tau\to 0$ then we should recover the Liouville equation (4) from ${\cal
U}^{t}$.

  3. 3.

If $\epsilon\to 0$, then from ${\cal U}^{t}$ we should recover the operator
governing the evolution of densities in a function space under the action of
the functional map

$x(t)=S(x(t-\tau)),$ (9)

for $t\in\mathbb{R}^{+}$, though we don’t know what that should beSharkovskiĭ,
Maĭstrenko, and Romanenko (1986).

###  II.1 Classifying density evolution dynamics

Just as dynamists classify the different types of trajectory
behavioursKuznetsov (2013), ergodic theorists have also classified different
types of convergence of density evolution (Ref. Lasota and Mackey, 1994,
Chapter 4). In this classification, we always let $S\colon X\to X$ be a
nonsingular transformation on a $\sigma$-finite measure space $(X,{\cal A},m)$
that preserves a probability measure $\mu$ with density $f_{*}$.

The weakest type of convergence is contained in the property of ergodicity.
$S$ is _ergodic_ if every invariant set $A\in{\cal A}$111$A$ is invariant if
$S^{-1}(A)=A$ is such that either $\mu(A)=0$ or $\mu(X\setminus A)=0$. This is
equivalent to the existence of a unique stationary density $f_{*}$ so
$Pf_{*}\equiv f_{*}$.

Next in the hierarchy is the stronger property of mixing. $S$ is _mixing_ if

$\lim_{t\to\infty}\mu(A\cap S_{t}^{-1}(B))=\mu(A)\mu(B)\quad\mbox{for all
}\,A,B\in{\cal A}.$

Mixing is equivalent to

$\lim_{t\to\infty}\langle P^{t}f,g\rangle=\langle f_{*},g\rangle$

for every bounded measurable function $g$.

Then we have the property of asymptotic stability (or exactness). Assume $S$
is such that $S(A)\in{\cal A}$ for each $A\in{\cal A}$. $S$ is asymptotically
stable if

$\lim_{t\to\infty}\mu(S_{t}(A))=1\quad\mbox{for all }\,A,B\in{\cal A}.$

Asymptotic stability is equivalent to

$f(x,t)\equiv P^{t}f_{0}(x)\to f_{*}(x)$

for all initial densities.

Finally there is asymptotic periodicity (Ref. Lasota and Mackey, 1994, Chapter
5.3). In this case, for all initial densities $f_{0}(x)$, there exists a
sequence of basis densities $g_{1},\ldots,g_{r}$ and a sequence of bounded
linear functionals $\lambda_{1},\ldots,\lambda_{r}$ such that

$Pf_{0}(x)=\sum_{i=1}^{r}\lambda_{i}(f_{0})g_{i}(x)+Qf_{0}(x),$

where $Q$ is an operator such that $\parallel P^{t}Qf\parallel\rightarrow 0$
as $t\rightarrow\infty$ for all integrable $f$. The densities $g_{j}$ have
disjoint supports and $Pg_{j}=g_{\alpha(j)}$, where $\alpha$ is a permutation
of $(1,\ldots,r)$. An invariant density is given by

$f_{*}=\dfrac{1}{r}\sum_{j=1}^{r}g_{j}.$

### Example: General hat map

The hat map is perfect to illustrate these various types of dynamics, since it
is knownIto, Tanaka, and Nakada (1979a, b) that (3) is ergodic for $a>1$ and
we have an analytic expressionYoshida, Mori, and Shigematsu (1983) for the
stationary density $f_{*}$. Furthermore, (3) is asymptotically periodic
Provatas and Mackey (1991a) with period $r=2^{n}$, $n=0,1,\cdots$ for

$2^{1/2^{{n+1}}}<a\leq 2^{1/2^{{n}}}.$

Thus, for example, $\\{P^{t}f\\}$ has period $1$ for $2^{1/2}<a\leq 2$, period
$2$ for $2^{1/4}<a\leq 2^{1/2}$, period $4$ for $2^{1/8}<a\leq 2^{1/4}$, etc.
Finally it is known Lasota and Mackey (1994) that (3) is exact for $a=2$.

###  II.2 Can dynamical systems display a ’chaotic’ evolution of densities?

The short answer is that nobody knows–it’s an open problem!

As pointed out in the Introduction, the trajectory sequence of potential
solution behaviors through bifurcations in dynamical or semi-dynamical systems
is

$\boxed{\mbox{stable steady state}\to\mbox{simple limit cycle}\to\mbox{
complicated limit cycle}\to\mbox{`chaotic' solutions}}$ (10)

and a great deal is known about the possible transitions between different
qualitative behaviours Kuznetsov (2013).

Analogously, the bifurcation structure in the evolution of sequences of
densities under the action of a Frobenius-Perron operator isLasota and Mackey
(1994)

$\boxed{\mbox{asymptotically stable stationary density}\to\mbox{simple
asymptotic periodicity}\to\mbox{complicated asymptotic periodicity}}$ (11)

but the analysis of bifurcations of densities is only in a rudimentary state
of development (Ref. Arnold, 1998, Chapter 9).

Thinking about these two different sequences raises the immediate and obvious
question “How could (can) one construct an evolution operator for densities
that would display a ‘chaotic’ evolution of densities?”. Thus, is

$\boxed{\mbox{asymptotically stable stationary density}\to\mbox{asymptotic
periodicity}\to\mbox{'chaotic' density evolution}}$ (12)

possible?

The Frobenius-Perron operator (1) is a linear operator, so our suspicion is
that in order to have a chaotic density evolution it would be necessary to
have a non-linear evolution operator. We have speculated elsewhereMackey
(2009); Mackey  _et al._ (2012); Mackey (2016) that maybe the density
evolution operator might need to be density dependent, thus leading to
nonlinearity.

### Example: A density dependent extension of the hat map

Consider a density dependent hat map

$x_{n+1}=\left\\{\begin{array}[]{ll}a[f_{n}]x_{n}&x_{n}\in[0,\frac{1}{2}]\\\
a[f_{n}](1-x_{n})&x_{n}\in(\frac{1}{2},],\end{array}\right.$ (13)

where $f_{n}$ is the density of $x_{n}$ and the functional $a[f]$ is defined
by

$a[f]=1+\int_{A}^{A+\delta}f(x)dx.$ (14)

The corresponding nonlinear evolution (pseudo-Frobenius-Perron) operator is
(Mackey, 2009)

$P_{f}f(x)=\dfrac{1_{[0,a[f]/2]}(x)}{a[f]}\left\\{f\left(\dfrac{x}{a[f]}\right)+f\left(1-\dfrac{x}{a[f]}\right)\right\\}.$
(15)

##  III Density evolution in differential delay equations

###  III.1 Asymptotic periodicity in a deterministic differential delay
equation

Let $x_{\tau}\equiv x(t-\tau)$ with $\tau=1$ and consider the hat map (3)
turned into a delay equationLosson and Mackey (1995a) (see Eq. (8) in
particular)

$\frac{dx}{dt}=-\alpha x+\left\\{\begin{array}[]{ll}ax_{\tau}&\mbox{if
$x_{\tau}<1/2$}\\\ a(1-x_{\tau})&\mbox{if $x_{\tau}\geq
1/2$}\end{array}\;\;\;\;\;\;\;\frac{a}{\alpha}\in(1,2],\right.$ (16)

and examine the result of picking many different initial functions and
following the trajectories forward in time. At successive times $t$ we sample
across all of the trajectories and form a histogram of the values of $x(t)$
that is an approximation to a ‘density’.

See Figure 1 where there is clear numerical evidence for the existence of
periodicity in the evolution of the histograms along the trajectories, and
which the authors in Ref Losson and Mackey, 1995a argued was evidence for
asymptotic periodicity of densities in this system, supported by their
analytic calculations. Note in particular in Figure 1 that a change in the
distribution of the initial functions changes the temporal sequence of
densities, but not the period.

Figure 1: Illustration of asymptotic periodicity in the density evolution of a
hat map differential delay equation (16) with $\alpha=13$ and $a=10$. The
simulation extends from $t=400$ to $t=402.9$ and is based on the integration
of $22,500$ initial functions. (a) Each initial function was a random process
distributed uniformly on $[0.65,0.75]$. (b) Here the initial functions were
uniformly distributed on $[0.65,0.75]$ for $17,000$ of the cases and
$[0.35,0.45]$ for the remaining $5,500$. It was observed that the cycling was
not transient and persisted for as long as the simulation ran. Modified from
Ref. Losson and Mackey, 1995a.

###  III.2 Asymptotic periodicity in stochastically perturbed delay equations

Asymptotic periodicity can be induced by noiseLasota and Mackey (1987);
Provatas and Mackey (1991b) in a Keener map

$S(x)=(ax+b)\,\,\pmod{1},\,\,0<a,b<1,$ (17)

_i.e._ when the dynamics are given by

$x_{n+1}=(ax_{n}+b+\xi_{n})\,\,\pmod{1},\,\,0<a,b<1,$ (18)

and the noise source $\xi$ is distributed with a density $\tilde{f}$. Consider
the Keener map (18) with noise $\xi$ turned into a stochastic delay
equationLosson and Mackey (1995a) (again see Eq. (8))

$\frac{dx}{dt}=-\alpha x+[(ax_{\tau}+b+\xi)\quad\mbox{mod}\,1]\,\,0<a,b<1$
(19)

and examine the evolution of many initial functions as shown in Figure 2.

In this figure there are two noteworthy features. First, an alteration in the
distribution of the noise with the distribution of initial functions kept the
same leads to an apparent qualitative change in the ‘density’ dynamics, going
from asymptotically stable in (b) to asymptotically periodic in (c). Secondly,
a change in the distribution of the initial functions with the distribution of
the noise kept constant [(c) to (d)] leads to a change in the details of the
temporal density evolution without a change in the period.

Figure 2: Noise induced apparent asymptotic periodicity and stability in the
stochastic differential delay equation (19). As in Figure 1, each simulation
was performed with $22,500$ random initial functions. In all four panels, the
parameters of the equation were $a=0.5,\;b=0.567,\;{\alpha}=10$. For panels
(a)-(c) the initial density was as in Figure 1(a). (a) No noise in the system:
$f(x,t)$ is not a density, but rather a generalized function. (b) Noise
supported uniformly on $[0,0.1]$. Numerically, the system appears to be
asymptotically stable (exact). (c) Noise uniformly supported on $[0,0.2]$. (d)
Same noise as in (c), with an initial density as in Figure 1(b). Modified from
Ref. Losson and Mackey, 1995a.

###  III.3 Deterministic Brownian motion

A typical formulationGardiner (1983) of the Brownian motion of a particle of
mass $m$ with position $x$ and velocity $v$ subject to a frictional force
$\gamma v$ is

$\displaystyle\dfrac{dx}{dt}$ $\displaystyle=$ $\displaystyle v,$ (20)
$\displaystyle m\dfrac{dv}{dt}$ $\displaystyle=$ $\displaystyle-\gamma
v+\eta(t),$ (21)

where $\eta$ is a fluctuating “force” due to collisions of the particle with
others of much small mass, and is usually given by $\eta(t)=\sigma\xi(t)$, and
$\xi=\frac{dw}{dt}$ is a ‘white noise’ (and delta correlated) which is the
‘derivative’ of a Wiener process $w(t)$. $\xi(t)$ is normally distributed with
mean $\mu=0$ and variance $\sigma=1$.

The question of whether one can produce a similar Brownian motion using a
totally deterministic model is interesting, and has been answered in the
affirmativeLei and Mackey (2011); Mackey and Tyran-Kamińska (2006).

### Brownian motion from a differential delay equation

In Ref. Lei and Mackey, 2011, the authors studied numerically the system

$\displaystyle\dfrac{dx}{dt}$ $\displaystyle=$ $\displaystyle v$ (22)
$\displaystyle\dfrac{dv}{dt}$ $\displaystyle=$ $\displaystyle-\gamma
v+\sin(2\pi\beta v(t-1)),$ (23) $\displaystyle v(t)$ $\displaystyle=$
$\displaystyle\phi(t),\,\,-1\leq t\leq 0.$ (24)

In this system the ‘random’ force (the sinusoidal term) is oscillating ever
more rapidly as $\beta$ increases. It was shown Lei and Mackey (2011), for a
variety of numerical situations, that the mean square displacement of the
particle obeys $[\Delta x(t)]^{2}\sim t$, while the velocity is distributed as
a quasi-Gaussian $\sim e^{-Cv^{2}}$ for $v\in[-K,K]$ (see Figure 3). The
numerics indicated that the bound $K$ and the standard deviation $\sigma$ are
given by

$\displaystyle K(\beta,\gamma)$ $\displaystyle=$
$\displaystyle\dfrac{1}{\sqrt{\gamma}(0.68\sqrt{\beta}+0.60\sqrt{\gamma})}$
(25) $\displaystyle\sigma(\beta,\gamma)$ $\displaystyle=$
$\displaystyle\dfrac{0.32}{\sqrt{\beta\gamma}}.$ (26)

Figure 3: Velocity distribution for four different values of $\beta$ with
$\gamma=1$ when the dynamics are given by (23). See Ref Lei and Mackey, 2011
for more detail.

### Brownian motion induced by perturbation from a non-invertible and chaotic
map

In Ref. Mackey and Tyran-Kamińska, 2006, the authors examined the properties
of the system

$\displaystyle\dfrac{dx}{dt}$ $\displaystyle=$ $\displaystyle v$ (27)
$\displaystyle\dfrac{dv}{dt}$ $\displaystyle=$ $\displaystyle-\gamma
v+\eta(t)$ (28)

where the fluctuating “force” consists of a series of delta-function like
impulses given by

$\eta(t)=\kappa(\tau)\sum_{j=0}^{\infty}h(\xi(t))\delta(t-j\tau)$

with $\tau>0$ and a scaling parameter $\kappa(\tau)>0$.

The real valued function $h$ is defined on a probability space
$(X,\mathcal{A},\mu)$ and $\xi$ is a highly chaotic deterministic variable
defined by $\xi(t)=\xi_{j}$ for $j\tau\leq t<(j+1)\tau$, $j\geq 0$, with
$\xi_{j+1}=S(\xi_{j})$ and $\xi_{0}$ the identity on $X$, where $S\colon X\to
X$ is an ergodic map with invariant measure $\mu$, _e.g._ $S$ is the hat map
(3) on $[0,1]$ and $h$ is the identity. It was shown in Ref. Mackey and Tyran-
Kamińska, 2006, Section 4 that the limit $\tau\to 0$ reproduces the
characteristics of an Ornstein-Uhlenbeck process which is the solution of the
stochastic differential equation

$\displaystyle dX(t)$ $\displaystyle=V(t)dt$ $\displaystyle dV(t)$
$\displaystyle=-\gamma V(t)dt+\sigma dw(t),$

where $\\{w(t):t\geq 0\\}$ is a Wiener process, provided

  1. 1.

$\kappa(\tau)^{2}/\tau$ converges to $1$ as $\tau\to 0$ and

  2. 2.

there is $\beta>1/2$ such that

$\limsup_{t\to\infty}t^{2\beta}\|P^{t}h\|_{1}<\infty,$

where $\|\cdot\|_{1}$ is the norm in $L^{1}(X,\mathcal{A},\mu)$ and $P^{t}$ is
the Frobenius-Perron operator (1) corresponding to $S_{t}$.

For precise statements and further examples see Mackey and Tyran-Kamińska
(2006); Tyran-Kamińska (2014).

##  IV What do these examples show?

So what are we displaying in Figures 1 through 3? Reference to Figure 4 makes
it clear that what we are examining is not really the evolution of the density
$\rho(x,t,t^{\prime})$ but rather $f(x,t)=\rho(x,t,t^{\prime}\equiv t)$.

This is not necessarily a bad thing because $f(x,t)=\rho(x,t,t^{\prime}\equiv
t)$ is the quantity that would typically be measured in an experimental
situation: We make measurements of $x$ at either discrete times, or a
continuum of times, and from these measurements construct temporal histograms
of some state variable as in Figures 1 or 2, or maybe we look at the long time
limiting behaviour as in Figure 3.

In Figure 4 we are looking at a snapshot of the evolution of all of these
trajectories emanating from many initial functions just as we would do in an
experiment. However, the unanswered question is how this is related to the
density evolving under the delayed dynamics? That is, how do we get from
$\rho(x,t,t^{\prime})$ to $f(x,t)=\rho(x,t,t^{\prime}\equiv t)$? This is
precisely the problem raised in Section II.

Figure 4: A schematic illustration of the connection between the evolution of
an ensemble of initial functions and what would be measured in a laboratory.
In the case that the delay has been scaled to $\tau=1$, an ensemble of $N$
initial functions on $[-1,0]$ is allowed to evolve forward in time under the
action of the delayed dynamics. At time $t$ we sample the distribution of the
values of $x$ across all $N$ trajectories and form an approximation to a
density $f(x,t)$. Modified from Ref. Losson and Mackey, 1995a.

##  V How to formulate?

We now turn to a consideration of the fundamental problem posed in this paper.
Namely, given a system whose evolution is determined by a differential delay
equation of the form (6) for example, and whose initial phase point $\phi\in
C$ is distributed over many possible initial states, how does the probability
distribution for the phase point evolve in time? Or equivalently given an
ensemble of independent systems, each governed by (6), and whose initial
functions are distributed according to some density over $C$, how does this
ensemble density evolve in time?

The answer to this question seems deceptively simple and would appear to be
provided by the Frobenius-Perron operator formalism of Section II. Suppose the
initial distribution of phase points is described by a probability measure
$\mu_{0}$ on $C$ so the probability that the initial function $\phi$ is an
element of $A\subset C$ is given by $\mu_{0}(A)$. Then, after a time $t$, the
new distribution is described by the measure $\mu_{t}$ given by

$\mu_{t}=\mu_{0}\circ S_{t}^{-1}$ (29)

if $S_{t}$ is a measurable transformation on $C$. Thus after a time $t$ the
probability that the phase point is an element of $A\subset C$ is
$\mu_{t}(A)=\mu_{0}(S_{t}^{-1}(A))$.

If the initial distribution of states $\phi$ is described by a density
$f(\phi)$ with respect to some measure $m$ and $S_{t}$ is a nonsingular
transformation, then after time $t$ the density will have evolved to $P^{t}f$,
where the Frobenius-Perron operator $P^{t}$ corresponding to $S_{t}$ is
defined by

$\int_{A}P^{t}f(\phi)\,m(d\phi)=\int_{S_{t}^{-1}(A)}f(\phi)\,m(d\phi)\quad\forall\text{
measurable }A\subset C.$ (30)

Equations (29)–(30) apparently answer the question posed above about the
evolution of probability measures for delay differential equations. However,
this is illusoryLosson  _et al._ (2020) since they provide only a symbolic
restatement of the problem. Everything that is specific to a given
differential delay equation is contained in $S_{t}^{-1}$.

Although the differential delay equation can be expressed in terms of an
evolution semigroup, there is no apparent way to invert the resulting
transformation $S_{t}$. This inversion will most certainly be non-trivial,
since solutions of delay equations often cannot be uniquely extended into the
past (Driver, 1977), and thus $S_{t}$ will not have a unique inverse.
$S_{t}^{-1}$ may have numerous branches that need to be accounted for when
evaluating $S_{t}^{-1}(A)$ in the Frobenius-Perron equation (30). This is a
serious barrier to deriving a closed-form expression for the Frobenius-Perron
operator $P^{t}$.

There are other subtle issuesLosson  _et al._ (2020) raised by equations
(29)–(30).

  * •

The most apparent difficulty is that the integrals in (30) are over sets in a
function space, and it is unclear how such integrals can be carried out.

  * •

It is also unclear what family of measures we are considering, and in
particular what subsets $A\subset C$ are measurable (_i.e._ , what is the
relevant $\sigma$-algebra on $C$?).

  * •

Also, in equation (30) what should be considered a natural choice for the
measure $m$ with respect to which probability densities are to be defined?

  * •

And finally, does it make sense to talk about probability densities in the
function space $C$?

These, then, are the rather formidable mathematical problems that we see in
any attempt to formulate a framework to describe the evolution of a ‘density’
under the action of a dynamical systems containing delays. These have
bedeviled us, as well as many of our collaborators, for a number of years. We
now turn to a brief consideration of failed attempts that we have published
over the years in an attempt to solve this problem. A detailed presentation of
these can be found in Ref. Losson  _et al._ , 2020.

### Possible lines of attack

  * •

Considerations related to an examination of fluid flow, and turbulence in
particular, bear some superficial similarities to the problems involved in
thinking about density evolution in delay systems because both are infinite
dimensional. In Ref. Losson  _et al._ , 2020, Chap. 5 an extension of the work
of Ref. Hopf, 1952 involving Hopf functionals is examined as a possible way of
tackling delayed density evolution. This is a exposition of work first
published in Ref. Losson and Mackey, 1992, and though a Hopf-like functional
differential equation governing the evolution of these densities was derived
it was not possible to proceed further with this method. Interestingly this
approach leads to a suggestive potential connection with the functional
expansions of quantum field theory and Feynman diagrams.

  * •

The method of stepsDriver (1977) is a classic tool used in the proof of the
existence and uniqueness of solutions of differential delay equations and Ref.
Losson  _et al._ , 2020, Chap. 6 examined its potential utility for a
formulation of the delay differential equation density problem. This approach,
while promising and useful from a numerical perspective, seems to lead only to
a weak solution and to not be of great utility.

  * •

Another possible approach is to take a differential delay equation and,
through a process of discretization (e.g. with an Euler approximation) turn it
into a high dimensional mapLosson and Mackey (1995b). Using this approach it
was then possible to prove certain properties for the density evolution of the
high dimensional map using results from Ref. Ionescu Tulcea and Marinescu,
1950, but it was never possible to successfully make the transition from these
results to a consideration of the discretized differential delay equation
(Ref. Losson  _et al._ , 2020, Chaps. 7,8).

None of these approaches have yielded a satisfactory way of dealing with the
problems enumerated at the beginning of this section and it is clear from the
extensive considerations of Ref. Losson  _et al._ , 2020 that a radically new
approach is required. In the next section we outline a possible approach and
show its application to a very simple example.

##  VI A novel approach

In this section we outline a new approach and show its potential applicability
to a specific and relatively simple, tractable, example.

### Liouville-like formulation

Let $C=C([-\tau,0],\mathbb{R}^{d})$ denote the Banach space of all continuous
functions $\phi\colon[-\tau,0]\to\mathbb{R}^{d}$ equipped with the supremum
norm and the Borel $\sigma$-algebra. Consider the equation

$\begin{split}x^{\prime}(t)=&\mathcal{F}(x(t),x(t-\tau)),\quad t\geq 0,\\\
x(t)=&\phi(t),\quad t\in[-\tau,0],\end{split}$ (31)

where $\phi\in C$ and $\mathcal{F}$ is such that for each $\phi\in C$ there
exists a continuous function $x\colon[-\tau,\infty)\to\mathbb{R}^{d}$ such
that (31) has a unique global solution depending continuously on $\phi$. For
each $t\geq 0$ define the solution map $S_{t}\colon C\to C$ by

$S_{t}(\phi)(s)=x_{t}(s)=x(t+s)\quad\text{for }s\in[-\tau,0],$ (32)

where $x(t)$ is a solution of (31) with $x_{0}=\phi$. It is well knownHale and
Verduyn Lunel (1993) that $\\{S_{t}\\}_{t\in\mathbb{R}^{+}}$ is a semi-
dynamical system on $C$ and the transformation $(t,\phi)\mapsto S_{t}(\phi)$
is continuous.

Let $B(C)$ be the space of bounded Borel measurable functions $\psi\colon
C\to\mathbb{R}$ with the supremum norm

$\|\psi\|_{\infty}=\sup_{\phi\in C}|\psi(\phi)|,$

and let $\mathcal{M}(C)$ (resp. $\mathcal{M}_{1}(C)$) denote the space of
finite (resp. probability) Borel measures on $C$. For any $\psi\in B(C)$ and
$\mu\in\mathcal{M}(C)$ we use the scalar product notation

$\langle\psi,\mu\rangle=\int_{C}\psi(\phi)\mu(d\phi).$

Let $C_{b}\subset B(C)$ be the Banach space of all bounded uniformly
continuous functions $\psi\colon C\to\mathbb{R}$ with the supremum norm. We
say that a family $\psi_{t}\in C_{b}$, $t>0$ _converges weakly_ Dynkin (1965);
Mohammed (1984) to $\psi\in C_{b}$ as $t\to 0^{+}$ (denoted by
$\mathrm{w-}\lim_{t\to 0}\psi_{t}=\psi$) if

$\lim_{t\to 0}\langle\psi_{t},\mu\rangle=\langle\psi,\mu\rangle\quad\text{for
each }\mu\in\mathcal{M}(C).$

This is equivalent to the following two conditions: $\lim_{t\to
0}\psi_{t}(\phi)=\psi(\phi)$ for every $\phi\in C$ and
$\sup_{t}\left\|\psi_{t}\right\|_{\infty}<\infty$.

A semigroup $\\{T^{t}\\}_{t\geq 0}$ of linear operators on the space $C_{b}$
is defined by

$T^{t}\psi(\phi)=\psi(S_{t}(\phi)),\quad t\geq 0,\;\phi\in C,\psi\in C_{b},$
(33)

where $S_{t}$ is the solution map (32). Since the semidynamical system
$\\{S_{t}\\}_{t\geq 0}$ is continuous, we obtain

$\lim_{t\to 0}T^{t}\psi(\phi)=\psi(\phi)\quad\text{for all }\phi\in
C,\;\psi\in C_{b}.$

Note also that

$\sup_{t\geq 0}\|T^{t}\psi\|_{\infty}\leq\|\psi\|_{\infty}.$

Consequently, the semigroup $\\{T^{t}\\}_{t\geq 0}$ is weakly continuous at
$t=0$. Define the _weak generator_
$\mathcal{L}\colon\mathcal{D}(\mathcal{L})\subset C_{b}\to C_{b}$ of the
semigroup $\\{T^{t}\\}_{t\geq 0}$ byDynkin (1965); Mohammed (1984)

$\displaystyle\mathcal{D}(\mathcal{L})$ $\displaystyle=\\{\psi\in
C_{b}:\mathrm{w-}\lim_{t\to
0}\frac{1}{t}\bigl{(}T^{t}\psi-\psi\bigr{)}\;\text{exists}\\},$
$\displaystyle\mathcal{L}\psi$ $\displaystyle=\mathrm{w-}\lim_{t\to
0}\frac{1}{t}\bigl{(}T^{t}\psi-\psi\bigr{)}.$

In particular, for $\psi\in\mathcal{D}(\mathcal{L})$ and
$\mu\in\mathcal{M}(C)$ we have

$\langle
T^{t}\psi,\mu\rangle=\langle\psi,\mu\rangle+\int_{0}^{t}\langle\mathcal{L}(T^{r}\psi),\mu\rangle
dr,\quad t>0.$ (34)

Let $\mu_{0}\in\mathcal{M}_{1}(C)$. For each $t\geq 0$, define the probability
measure $\mu_{t}$ on the space $C$ by $\mu_{t}=\mu_{0}\circ S_{t}$ as in (29),
where $S_{t}$ is the solution map (32). Then we have

$\langle\psi,\mu_{t}\rangle=\langle T^{t}\psi,\mu_{0}\rangle,\quad t\geq
0,\psi\in C_{b}.$

It follows from (34) and the change of variables formula that

$\langle\psi,\mu_{t}\rangle=\langle\psi,\mu_{0}\rangle+\int_{0}^{t}\langle\mathcal{L}\psi,\mu_{r}\rangle
dr$ (35)

for all $t>0$ and $\psi\in\mathcal{D}(\mathcal{L})$. However, it is difficult
to identify the domain of the weak generator.

Thus, we introduce an _extended generator_ for the solution map $S_{t}$. It is
defined as a linear operator $\mathcal{L}$ from its domain $\mathcal{D}$ to
the set of all Borel measurable functions on $C$, where we say that
$\psi\in\mathcal{D}$ if for each $t>0$ we have

$\int_{0}^{t}\langle|\mathcal{L}\psi|,\mu_{r}\rangle dr<\infty$ (36)

and (35) holds. Then we say that $\\{\mu_{t}\\}_{t\geq 0}$ is the solution of
the equation

$\frac{\partial}{\partial t}\mu_{t}=\mathcal{L}^{*}\mu_{t}.$

Instead of defining the whole domain $\mathcal{D}$, one can start with a
sufficiently large subset of $\mathcal{D}$. A typical example is
$\mathcal{F}C_{c}^{\infty}(C)$ the set of smooth cylinder functions on $C$.
These are functions of the form
$\psi(\phi)=g(l_{1}(\phi),\ldots,l_{n}(\phi))$, $\phi\in C$, $n\geq 1$, for
some continuous linear functionals $l_{1},\ldots,l_{n}$ on $C$, and $g\in
C_{c}^{\infty}(\mathbb{R}^{n})$. Another one example of $\mathcal{D}$ is the
set of quasi-tame functions introduced by Mohammed Mohammed (1984).

The marginal distribution $\mu(t)$ of the measure $\mu_{t}=\mu_{0}\circ S_{t}$
is defined on $\mathbb{R}^{d}$ by

$\mu(t)(B)=\mu_{t}\\{\phi\in C:\phi(0)\in B\\},\quad
B\in\mathcal{B}(\mathbb{R}^{d}).$

This can be rewritten with the projection map $\pi_{0}\colon
C\to\mathbb{R}^{d}$ defined by $\pi_{0}(\phi)=\phi(0)$, $\phi\in C$, as
$\mu(t)=\mu_{t}\circ\pi_{0}^{-1}$. Note that $\mu(t)$ is the distribution of
$x(t)$ for all $t\geq 0$.

Let $C_{c}^{1}(\mathbb{R}^{d})$ denote the space of functions that are
continuously differentiable and have compact support. Consider the
differential operator from the space $C_{c}^{1}(\mathbb{R}^{d})$ to the set of
all Borel measurable functions on $C$ defined by

$Lg(\phi)=\sum_{i=1}^{d}\mathcal{F}_{i}(\phi(0),\phi(-\tau))\frac{\partial
g}{\partial x_{i}}(\phi(0)),\quad\phi\in C,g\in C_{c}^{1}(\mathbb{R}^{d}).$

Let $g\in C_{c}^{1}(\mathbb{R}^{d})$ be such that $Lg\in B(C)$. If
$\psi=g\circ\pi_{0}$, then $\psi\in\mathcal{D}$ and $\mathcal{L}\psi=Lg$.

We say that $\\{\mu_{t}\\}_{t\geq 0}$ is a solution of the equation

$\frac{\partial}{\partial t}\mu(t)=L^{*}\mu_{t}$ (37)

if for each $t>0$ and $g\in C_{c}^{1}(\mathbb{R}^{d})$ the following holds

$\int_{0}^{t}\langle|Lg|,\mu_{r}\rangle dr<\infty$ (38)

and

$\int_{\mathbb{R}^{d}}g(x)\mu(t)(dx)=\int_{\mathbb{R}^{d}}g(x)\mu(0)(dx)+\int_{0}^{t}\langle
Lg,\mu_{r}\rangle dr.$ (39)

Note that this is equivalent to requiring that $\psi\in\mathcal{D}$ and
$\mathcal{L}\psi=Lg$ for $\psi=g\circ\pi_{0}$ and all $g\in
C^{1}_{c}(\mathbb{R}^{d})$. Next, observe that if we introduce the measure

$\nu(t)=\mu_{t}\circ\pi_{0,-\tau}^{-1},\quad t\geq 0,$

where $\pi_{0,-\tau}\colon C\to\mathbb{R}^{d}\times\mathbb{R}^{d}$ is the
projection map $\pi_{0,-\tau}(\phi)=(\phi(0),\phi(-\tau))$, $\phi\in C$, then

$\langle
Lg,\mu_{t}\rangle=\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\sum_{i=1}^{d}\mathcal{F}_{i}(x,y)\frac{\partial
g}{\partial x_{i}}(x)\nu(t)(dx,dy),$

by the change of variables formula. The measure $\nu(t)$ is the distribution
of $(x(t),x(t-\tau))$.

Now suppose that the measure $\nu(t)$ has a density $f_{\nu}$ with respect to
the Lebesgue measure on $\mathbb{R}^{2d}$ _i.e._ ,
$\nu(t)(dx,dy)=f_{\nu}(x,y,t)dxdy$. Then the measure $\mu(t)$ has a density
$f$ with respect to the Lebesgue measure on $\mathbb{R}^{d}$ and

$f(x,t)=\int_{\mathbb{R}^{d}}f_{\nu}(x,y,t)dy.$

We also have

$f(y,t-\tau)=\int_{\mathbb{R}^{d}}f_{\nu}(x,y,t)dx,\quad t\geq\tau.$

We can rewrite (39) as

$\displaystyle\int_{\mathbb{R}^{d}}g(x)f(x,t)dx=\int_{\mathbb{R}^{d}}g(x)f_{0}(x)dx$
$\displaystyle+\int_{0}^{t}\int_{\mathbb{R}^{d}}\sum_{i=1}^{d}\mathcal{F}_{i}(x,y)\frac{\partial
g}{\partial x_{i}}(x)f_{\nu}(x,y,r)dxdydr,$

which is the weak form of the equation

$\frac{\partial}{\partial t}f(x,t)=-\sum_{i=1}^{d}\frac{\partial}{\partial
x_{i}}\int_{\mathbb{R}^{d}}\mathcal{F}_{i}(x,y)f_{\nu}(x,y,t)dy.$ (40)

Eq. (40) reduces to the Liouville equation if $\mathcal{F}$ does not depend on
$y$.

Finally, recall that a Borel probability measure $\mu$ on $C$ is called a
Gaussian measure (see Ref. Bogachev, 1998) if the measure $\mu\circ l^{-1}$ is
Gaussian on $\mathbb{R}$ for each continuous linear functional $l\in C^{*}$.
Suppose that the solution map $S_{t}$ is linear. Then if we take as $\mu_{0}$
a Gaussian measure on $C$, the measure $\mu_{t}$ will be again Gaussian. We
will look at such examples next.

### A tractable example: Gaussian processes

A Gaussian process is a family $\xi=\\{\xi(t)\\}_{t\in\mathbb{T}}$ of (real-
valued) random variables defined on some probability space
$(\Omega,\Sigma,\mathbb{P})$, indexed by a parameter set $\mathbb{T}$, such
that every finite linear combination $\sum c_{t_{i}}\xi(t_{i})$ is either
identically zero or has a Gaussian distribution on $\mathbb{R}$. Given a
Gaussian process $\xi$ its mean function is $a(t)=\mathbb{E}(\xi(t))$,
$t\in\mathbb{T}$, and its covariance function is the bivariate symmetric
function

$R(t,s)=\mathrm{cov}(\xi(t),\xi(s))=\mathbb{E}(\xi(t)-a(t))(\xi(s)-a(s)).$

The process is called centered if its mean function is zero. Note that a
covariance function is non-negative definite, i.e.

$\sum_{i,j=1}^{n}R(t_{i},t_{j})c_{i}c_{j}\geq 0$ (41)

for any $t_{1},\ldots,t_{n}\in\mathbb{T}$ and
$c_{1},\ldots,c_{n}\in\mathbb{R}$, $n\geq 1$. A Gaussian process is said to be
non-degenerate if its covariance function is positive definite, i.e.
inequality (41) is strict for all non-zero
$(c_{1},\ldots,c_{n})\in\mathbb{R}^{n}$.

For a centered Gaussian process the Gaussian distribution of a finite
combination $\sum c_{i}\xi(t_{i})$ is determined through the variance
$\mathbb{E}(\sum c_{i}\xi(t_{i}))^{2}$ that can be calculated with the help of
the covariance function $R$. Thus the covariance function of a centered
Gaussian process completely determines all of the finite-dimensional
distributions (that is, the joint distributions of
$(\xi(t_{1}),\ldots,\xi(t_{n}))$ for any $t_{1},\ldots,t_{n}\in\mathbb{T}$ and
$n\geq 1$). Consequently the distribution of the entire centered Gaussian
process $\xi$ is uniquely determined through its covariance function.

Given a symmetric non-negative definite function
$R\colon[-\tau,0]\times[-\tau,0]\to\mathbb{R}$ there exists a centered
Gaussian process with $R$ being its covariance function, by the Kolmogorov
extension theorem. If the process has continuous sample paths then it defines
a Gaussian measure on $C$, the distribution of the process $\xi$, by

$\mu(B)=\mathbb{P}(\xi\in B),\quad B\in\mathcal{B}(C).$

If we have two centered Gaussian processes with continuous sample paths and
the same covariance function then they have the same distribution.

Consider the Gaussian process

$\xi(t)=\zeta\cos(t-\theta),\quad t\in\mathbb{R},$

where $\zeta$ and $\theta$ are independent random variables, the amplitude
$\zeta$ has the Rayleigh distribution with density $xe^{-x^{2}/2}$, $x\geq 0$,
and $\theta$ is uniformly distributed on $[0,2\pi)$. Observe that the mean
$\mathbb{E}(\xi(t))$ is zero and the covariance function is

$R(t,s)=\mathrm{cov}(\xi(t),\xi(s))=\mathbb{E}(\xi(t)\xi(s))=\cos(t-s).$

Then $\xi(t)$ is normally distributed with mean $0$ and variance $1$. Since
$R(t,t-\tau)=0$ with $\tau=\pi/2$ and $(\xi(t),\xi(t-\tau))$ has a Gaussian
distribution on $\mathbb{R}^{2}$, we see that $\xi(t)$ and $\xi(t-\tau)$ are
uncorrelated, thus independent. Note that $\xi$ is a solution of

$\dfrac{dx}{dt}=-x\Big{(}t-\frac{\pi}{2}\Big{)}.$ (42)

We take as $\mu_{0}$ the Gaussian measure being the distribution of the
Gaussian process $\xi=\\{\xi(s)\\}_{s\in[-\tau,0]}$. Then $S_{t}\xi$ is a
centered Gaussian process with the same covariance function as $\xi$. Hence
$\mu_{t}=\mu_{0}\circ S_{t}^{-1}=\mu_{0}$ for all $t\geq 0$ implying that the
measure $\mu_{0}$ is invariant for the solution map $S_{t}$ of equation (42).
We also have $f(x,t)=f(x)$ and $f_{\nu}(x,y,t)=f(x)f(y)$, where $f(x)$ is the
density of the standard normal distribution. Hence (40) holds with $d=1$ and
$\mathcal{F}(x,y)=-y$. However, if we take $\xi(t)=\zeta\cos(t)$,
$t\in\mathbb{R}$, where now $\zeta$ has the standard normal distribution, then
$\xi$ is a Gaussian process with $R(t,s)=\cos(t)\cos(s)$ and $\xi(t)$ is also
a solution of (42), but $\xi$ is degenerate since the distribution of
$(\xi(t),\xi(t-\tau))=\zeta(\cos(t),\sin(t))$ is concentrated on a line.

We can extend (42) in the following way. Consider the linear differential
delay equation

$x^{\prime}(t)=ax(t)+bx(t-\tau),$ (43)

where $a,b$ are real constants and a Gaussian process
$\xi=\\{\xi(s)\\}_{s\in[-\tau,0]}$ with continuous sample paths on $[-\tau,0]$
as the initial condition $x_{0}=\xi$. Then the solution map $S_{t}$ is a
linear mapping of $C$, and thus $S_{t}\xi$ is a Gaussian process, so the
distribution of $S_{t}\xi$ is a Gaussian measure on $C$.

Thus we can start with a positive definite symmetric function $R_{0}$ such
that the corresponding Gaussian process $\xi$ has continuous sample paths on
$[-\tau,0]$. Then $S_{t}\xi$ will be a centered Gaussian process with a
covariance function $R_{t}$. The measure $\mu(t)$, being the distribution of
$x(t)=S_{t}\xi(0)$, is Gaussian with mean zero $0$ and variance
$\sigma^{2}(t)=R_{t}(0,0)$ and the measure $\nu(t)$, being the distribution of
$(x(t),x(t-\tau))=(S_{t}\xi(0),S_{t}\xi(-\tau))$, is Gaussian with covariance
matrix

$Q_{t}=\left(\begin{array}[]{cc}R_{t}(0,0)&R_{t}(-\tau,0)\\\
R_{t}(-\tau,0)&R_{t}(-\tau,-\tau)\\\ \end{array}\right).$

Note that $\sigma^{2}(t-\tau)=R_{0}(t-\tau,t-\tau)$ for $t\in[0,\tau)$ and
$\sigma^{2}(t-\tau)=R_{t}(-\tau,-\tau)$ for $t\geq\tau$. Now if
$\sigma^{2}(t)\neq 0$ then the density $f$ of $\mu(t)$ is

$f(x,t)=\frac{1}{\sqrt{2\pi\sigma^{2}(t)}}e^{-\frac{x^{2}}{2\sigma^{2}(t)}}$
(44)

and if $\det(Q_{t})\neq 0$ then the density $f_{\nu}$ of $\nu(t)$ is given by

$f_{v}(x,y,t)=\frac{1}{2\pi\sqrt{\det(Q_{t})}}e^{-\frac{\sigma^{2}(t-\tau)x^{2}+2R_{t}(-\tau,0)xy+\sigma^{2}(t)y^{2}}{2\det(Q_{t})}}.$
(45)

Observe that in this example Eq. (40) is of the form

$\frac{\partial}{\partial t}f(x,t)=-\frac{\partial}{\partial
x}(axf(x,t))-\frac{\partial}{\partial x}\int_{\mathbb{R}}byf_{\nu}(x,y,t)dy$

and that

$\int_{\mathbb{R}}yf_{\nu}(x,y,t)dy=\frac{R_{t}(-\tau,0)}{\sigma^{2}(t)}xf(x,t).$

Thus $f$ is a solution of

$\frac{\partial}{\partial
t}f(x,t)=-\Big{(}a+b\frac{R_{t}(-\tau,0)}{\sigma^{2}(t)}\Big{)}\frac{\partial}{\partial
x}(xf(x,t)).$ (46)

It is easily seen that $f$ satisfies Eq. (46) if and only if

$\frac{d}{dt}(\sigma^{2}(t))=2a\sigma^{2}(t)+2bR_{t}(-\tau,0),\quad t>0.$ (47)

It follows from (43) that (47) holds. We see that

$\sigma^{2}(t)=e^{2at}\sigma^{2}(0)+2b\int_{0}^{t}e^{2a(t-s)}R_{s}(-\tau,0)ds,\quad
t>0.$

Consequently, to obtain the densities $f$ and $f_{\nu}$ it is enough to
determine first $R_{t}(-\tau,0)$ and then $\sigma^{2}(t)$.

To find the covariance function $R_{t}$ observe that we can write the solution
of (43) as (see Section 1.6 in Ref. Hale and Verduyn Lunel, 1993)

$x(t)=X(t)\xi(0)+b\int_{-\tau}^{0}X(t-r-\tau)\xi(r)dr,$

where $X(t)$ is the fundamental solution of (43), i.e. $X(t)=0$ for $t<0$ and

$X(t)=\sum_{k=0}^{\lfloor
t/\tau\rfloor}\frac{b^{k}}{k!}(t-k\tau)^{k}e^{a(t-k\tau)},\quad t\geq 0,$ (48)

where $\lfloor s\rfloor=\max\\{k\in\mathbb{Z}:k\leq s\\}$. We rewrite $x(t)$
with the help of the Lebesgue-Stieltjes integral as

$x(t)=\int_{-\tau}^{0}\xi(r)dM_{t}(r),$

where the function $M_{t}\colon[-\tau,0]\to\mathbb{R}$ is defined by

$dM_{t}(r)=bX(t-r-\tau)dr+X(t)\delta_{0}(dr),\quad t\geq 0,$

and

$dM_{t}(r)=\delta_{t}(dr),\quad t<0,$

with $\delta_{t}$ denoting the point measure at $t$. Since

$S_{t}\xi(s)=\int_{-\tau}^{0}\xi(r)dM_{t+s}(r),$ (49)

we obtain

$S_{t}\xi(s_{1})S_{t}\xi(s_{2})=\int_{-\tau}^{0}\int_{-\tau}^{0}\xi(r_{1})\xi(r_{2})dM_{t+s_{1}}(r_{1})dM_{t+s_{2}}(r_{2}).$

Taking the expectation on both sides of the above leads to

$R_{t}(s_{1},s_{2})=\int_{-\tau}^{0}\int_{-\tau}^{0}R_{0}(r_{1},r_{2})dM_{t+s_{1}}(r_{1})dM_{t+s_{2}}(r_{2}).$
(50)

Since the covariance function
$R_{t}\colon[-\tau,0]\times[-\tau,0]\to\mathbb{R}$ is symmetric, we can assume
that $-\tau\leq s_{1}\leq s_{2}\leq 0$. If $t\in[0,-s_{2}]$ then

$R_{t}(s_{1},s_{2})=R_{0}(t+s_{1},t+s_{2}),$

for $t\in(-s_{2},-s_{1}]$ we have

$\displaystyle R_{t}(s_{1},s_{2})$
$\displaystyle=e^{a(t+s_{2})}R_{0}(0,t+s_{1})$
$\displaystyle\quad+b\int_{-\tau}^{t+s_{2}-\tau}e^{a(t+s_{2}-r)}R_{0}(r,t+s_{1})dr$

and if $t>-s_{1}$ then

$\displaystyle R_{t}(s_{1},s_{2})$
$\displaystyle=X(t+s_{1})X(t+s_{2})R_{0}(0,0)$
$\displaystyle\quad+bX(t+s_{1})\int_{-\tau}^{0}X(t+s_{2}-r-\tau)R_{0}(r,0)dr$
$\displaystyle\quad+bX(t+s_{2})\int_{-\tau}^{0}X(t+s_{1}-r-\tau)R_{0}(r,0)dr$
$\displaystyle\quad+b^{2}\int_{-\tau}^{0}\int_{-\tau}^{0}X(t+s_{1}-r_{1}-\tau)X(t+s_{2}-r_{2}-\tau)$
$\displaystyle\quad\quad\quad\times R_{0}(r_{1},r_{2})dr_{1}dr_{2}.$

Thus we can find the covariance function $R_{t}$ by specifying the covariance
function $R_{0}$ and using (48) in the above equation.

Define

$\alpha_{0}=\max\\{\mathrm{Re}\lambda:\lambda=a+be^{-\lambda\tau}\\}.$

By Theorem 5.2 in Chapter 1 of Ref. Hale and Verduyn Lunel, 1993 for each
$\alpha>\alpha_{0}$ there is a constant $c$ such that $|X(t)|\leq ce^{\alpha
t}$ for all $t>0$. In particular, if $\alpha_{0}<0$ then $X(t)$ converges to
zero exponentially fast. Since $R_{0}$ being a continuous function is bounded,
we see that $R_{t}$ approaches $0$ if $\alpha_{0}<0$. Based on the work of
Hayes (1950), see also Ref. Hale and Verduyn Lunel, 1993, Section 5.2 and Thm.
A.5, we have $\alpha_{0}<0$ if and only if

$\displaystyle a\tau<1,\quad b\tau+a\tau<0,\quad
b\tau+a\tau\cos\kappa+\kappa\sin\kappa>0,$

where $\kappa$ is the root of $\kappa=a\tau\tan\kappa$, $0<\kappa<\pi$ if
$a\neq 0$ and $\kappa=\pi/2$ if $a=0$. These are values of $(-a,-b)$ inside
the cusp like area of Ref. Hale and Verduyn Lunel, 1993, Figure 5.1. Then we
will have

$\lim_{t\to\infty}\sigma^{2}(t)=0$

leading to

$\lim_{t\to\infty}\int_{-\varepsilon}^{\varepsilon}f(x,t)dx=1\quad\text{for
all }\varepsilon>0,$

by Chebyshev’s inequality

$1-\int_{-\varepsilon}^{\varepsilon}f(x,t)dx=\mathbb{P}(|x(t)|>\varepsilon)\leq\frac{1}{\varepsilon^{2}}\sigma^{2}(t).$

If, for example $R_{0}$ is nonnegative with $\sigma^{2}(0)=R_{0}(0,0)>0$ and
$a\geq 0$, $b>0$ then $\sigma^{2}(t)\geq X(t)^{2}\sigma^{2}(0)$ and
$X(t)\to\infty$ as $t\to\infty$. Thus

$\lim_{t\to\infty}\sigma^{2}(t)=\infty$

implying that

$\lim_{t\to\infty}\int_{-\varepsilon}^{\varepsilon}f(x,t)dx=0,$

since

$f(x,t)\leq\frac{1}{\sqrt{2\pi\sigma^{2}(t)}},\quad x\in\mathbb{R},t>0.$

It might also happen that $\sigma^{2}(t)$ is a constant, as it was for Eq.
(42) and $R_{0}(s_{1},s_{2})=\cos(s_{2}-s_{1})$.

Suppose now that the covariance function $R_{0}$ can be written in the form

$R_{0}(s_{1},s_{2})=\int_{-\tau}^{0}\eta_{r}(s_{1})\eta_{r}(s_{2})dr,$ (51)

for some function $\eta_{r}\colon[-\tau,0]\to\mathbb{R}$ such that
$(s,r)\mapsto\eta_{r}(s)$ is Borel measurable. Then it follows from (50) and
(49) that the covariance $R_{t}$ is given by

$R_{t}(s_{1},s_{2})=\int_{-\tau}^{0}S_{t}\eta_{r}(s_{1})S_{t}\eta_{r}(s_{2})dr$
(52)

and in particular, we have

$\sigma^{2}(t)=\int_{-\tau}^{0}(S_{t}\eta_{r}(0))^{2}dr.$

One example of (51) is $R_{0}(s_{1},s_{2})=\min\\{s_{1},s_{2}\\}+\tau$ for
$s_{1},s_{2}\in[-\tau,0]$, since (51) holds with

$\eta_{r}(s)=1_{[r,0]}(s),\quad s,r\in[-\tau,0].$ (53)

Then we have $\xi(s)=W(s+\tau)$, where $W=\\{W(t)\\}_{t\geq 0}$ is the
standard Wiener process on $[0,\infty)$, Another one is

$R_{0}(s_{1},s_{2})=\left\\{\begin{array}[]{ll}u(s_{1})v(s_{2}),&s_{1}\leq
s_{2},\\\ u(s_{2})v(s_{1}),&s_{1}>s_{2},\end{array}\right.$

where $u,v$ are nonnegative differentiable functions with $u(-\tau)=0$ and
$v(s)>0$. Here we can define $\xi(s)=v(s)W(\frac{u(s)}{v(s)})$ and we get (52)
with

$\eta_{r}(s)=v(s)\sqrt{\frac{u(r)}{v(r)}}1_{[-\tau,s]}(r)\quad
s,r\in[-\tau,0].$

Finally, we calculate the covariance $R_{t}$ as in (52) for $t\in(0,\tau]$,
when $\eta_{r}$ is given by (53), so that $S_{t}\eta_{r}(s)=1_{[r,0]}(t+s)$
for $t+s<0$ and

$S_{t}\eta_{r}(s)=X(t+s)+b\int_{r}^{0}X(t+s-q-\tau)dq$

for $s,r\in[-\tau,0]$, $t>-s$. Observe that

$\displaystyle
S_{t}\eta_{r}(s)=\left\\{\begin{array}[]{ll}1+b(t+s-r-\tau)_{+},&a=0,\\\
e^{a(t+s)}+\frac{b}{a}(e^{a(t+s-r-\tau)_{+}}-1),&a\neq 0,\end{array}\right.$

for $t+s\in[0,\tau]$, where $(q)_{+}=\max\\{0,q\\}$. Let $-\tau\leq s_{1}\leq
s_{2}\leq 0$. We have

$R_{t}(s_{1},s_{2})=t+s_{1}+\tau,\quad t\in[0,-s_{2}],$

but if $t\in(-s_{2},-s_{1}]$ then

$R_{t}(s_{1},s_{2})=\left\\{\begin{array}[]{ll}t+s_{1}+\tau+\frac{b}{2}(t+s_{2})^{2},&a=0,\\\
e^{a(t+s_{2})}(t+s_{1}+\tau)&\\\
+\frac{b}{a^{2}}\big{(}e^{a(t+s_{2})}-1-a(t+s_{2})\big{)},&a\neq
0;\end{array}\right.$

and if $-s_{1}<t\leq\tau-s_{2}$, then for $a=0$ we obtain

$\displaystyle R_{t}(s_{1},s_{2})$
$\displaystyle=\tau+\frac{b}{2}(t+s_{2})^{2}+\frac{b}{2}(t+s_{1})^{2}$
$\displaystyle\quad+\frac{b^{2}}{3}(t+s_{1})^{3}+\frac{b^{2}}{2}(t+s_{1})^{2}(s_{2}-s_{1})$

while for $a\neq 0$ we have

$\displaystyle R_{t}(s_{1},s_{2})$
$\displaystyle=e^{a(2t+s_{1}+s_{2})}\tau+\frac{b}{a^{2}}e^{a(t+s_{1})}\big{(}e^{a(t+s_{2})}-1-a(t+s_{2})\big{)}$
$\displaystyle\quad+\frac{b}{a^{2}}\big{(}e^{a(t+s_{2})}-\frac{b}{a}\big{)}\big{(}e^{a(t+s_{1})}-1-a(t+s_{1})\big{)}$
$\displaystyle\quad+\frac{b^{2}}{2a^{3}}e^{a(s_{2}-s_{1})}\big{(}e^{a(t+s_{1})}-1\big{)}^{2}.$

In particular, we see that for $t\in[0,\tau]$

$R_{t}(-\tau,0)=\left\\{\begin{array}[]{ll}t+\frac{b}{2}t^{2},&a=0,\\\
e^{at}t+\frac{b}{a^{2}}\big{(}e^{at}-1-at\big{)},&a\neq 0,\end{array}\right.$

and

$\sigma^{2}(t)=\left\\{\begin{array}[]{ll}\tau+bt^{2}+\frac{b^{2}}{3}t^{3},&a=0,\\\
e^{2at}\tau+\frac{2b}{a^{2}}e^{at}\big{(}e^{at}-1-at\big{)}&\\\
+\frac{b^{2}}{2a^{3}}\big{(}(e^{at}-1)^{2}-2(e^{at}-1-at)\big{)},&a\neq
0.\end{array}\right.$

Since $\sigma^{2}(t)>0$ and $\det(Q_{t})>0$ for $t\in(0,\tau]$, the densities
$f$ in (44) and $f_{\nu}$ in (45) are well defined. However, it should be
noted that the Frobenius-Perron operator as in (30) will not be well defined
if we take as a reference measure $m$ on $C$ the Gaussian measure $\mu_{0}$.
Suppose, on the contrary, that $P^{t}$ is well defined. According to the
Feldman-Hájek theoremBogachev (1998), two Gaussian measures are either
equivalent (mutually absolutely continuous) or singular. The Gaussian measure
$\mu_{t}$, being the distribution of the process $S_{t}\xi$ with covariance
function $R_{t}$, is equivalent to the Gaussian measure $\mu_{0}$Shepp (1966);
Park (1972) if and only if

$R_{t}(s_{1},s_{2})=\min\\{s_{1},s_{2}\\}+\tau-\int_{-\tau}^{s_{1}}\int_{-\tau}^{s_{2}}K_{t}(r_{1},r_{2})dr_{1}dr_{2},$

where $K_{t}\in L^{2}([-\tau,0]^{2})$ is a symmetric function such that the
integral operator

$K_{t}\phi(s)=\int_{-\tau}^{0}K_{t}(s,r)\phi(r)dr$

does not have eigenvalue $1$. The kernel $K_{t}$ is unique and for a.e.
$(s_{1},s_{2})$ satisfies

$K_{t}(s_{1},s_{2})=-\frac{\partial}{\partial s_{1}}\frac{\partial}{\partial
s_{2}}R_{t}(s_{1},s_{2}).$

The formula that we have obtained for $R_{t}$ implies that the function
$K_{t}$ does not exist in this example, thus showing that $\mu_{t}$ is
singular with respect to $m=\mu_{0}$. Hence, there exists a measurable set
$A\subset C$ such that $m(A)=0$ and $m(S_{t}^{-1}(A))=1$, contradicting (30)
with $f(\phi)\equiv 1$. Consequently, the approach presented here in Section
VI may be more effective than the ones we briefly mentioned in Section V.

##  VII Summary

Here we have highlighted an open mathematical problem. Namely, how can one
formulate and study the evolution of densities in systems with dynamics that
contain delays. This is not simply an abstract mathematical problem devoid of
interest. Rather it is of prime scientific interest because of the increasing
prevalence of studies of systems whose dynamics contain significant time
delays, and the fact that we have no way of theoretically treating these
systems. Various approximations are available, and we have briefly discussed
some of these while mentioning that there is a much more extensive
consideration of theseLosson  _et al._ (2020). We have also presented a
tentative new approach to the problem in Section VI.

Tied into this problem is the related issue of the lack of a well developed
theory of bifurcation patterns of densities analogous to what exists for
trajectories of dynamical systems. In this vein we have also raised the
question of whether it is possible for semi-dynamical systems to display a
chaotic pattern of density evolution–currently an open question.

##  VIII Data Availability Statement.

Data sharing is not applicable to this article as no new data were created or
analyzed in this study.

###### Acknowledgements.

We would like to acknowledge research support from the NSERC (Natural Sciences
and Engineering Research Council of Canada), MITACS (Mathematics of
Information Technology and Complex Systems), the Alexander von Humboldt
Stiftung and the National Science Centre (Poland) with grant No.
2017/27/B/ST1/00100. Additionally we have greatly benefited from conversations
with our colleagues Jinzhi Lei, Jérôme Losson, Nicholas Provatas, and Richard
Taylor.

## References

  * Hunt and Ott (2015) B. Hunt and E. Ott, “Defining chaos,” Chaos 5, 097618 (2015).
  * Lasota and Mackey (1994) A. Lasota and M. C. Mackey,  _Chaos, fractals, and noise: Stochastic aspects of dynamics_ , Applied Mathematical Sciences, Vol. 97 (Springer-Verlag, New York, 1994).
  * Gardiner (1983) C. Gardiner,  _Handbook of Stochastic Methods_ (Springer Verlag, Berlin, Heidelberg, 1983).
  * Sharkovskiĭ, Maĭstrenko, and Romanenko (1986) A. N. Sharkovskiĭ, Y. L. Maĭstrenko, and E. Y. Romanenko,  _Raznostnye uravneniya i ikh prilozheniya_ (“Naukova Dumka”, Kiev, 1986).
  * Kuznetsov (2013) Y. A. Kuznetsov,  _Elements of applied bifurcation theory_ , Vol. 112 (Springer Science & Business Media, 2013).
  * Note (1) $A$ is invariant if $S^{-1}(A)=A$.
  * Ito, Tanaka, and Nakada (1979a) S. Ito, S. Tanaka, and H. Nakada, “On unimodal linear transformations and chaos. I,” Tokyo J. Math. 2, 221–239 (1979a).
  * Ito, Tanaka, and Nakada (1979b) S. Ito, S. Tanaka, and H. Nakada, “On unimodal linear transformations and chaos. II,” Tokyo J. Math. 2, 241–259 (1979b).
  * Yoshida, Mori, and Shigematsu (1983) T. Yoshida, H. Mori, and H. Shigematsu, “Analytic study of chaos of the tent map: band structures, power spectra, and critical behaviors,” J. Statist. Phys. 31, 279–308 (1983).
  * Provatas and Mackey (1991a) N. Provatas and M. C. Mackey, “Asymptotic periodicity and banded chaos,” Phys. D 53, 295–318 (1991a).
  * Arnold (1998) L. Arnold,  _Random dynamical systems_ , Springer Monographs in Mathematics (Springer-Verlag, Berlin, 1998).
  * Mackey (2009) M. C. Mackey, “Exploring the world with mathematics,” Ann. Math. Sil. 23, 11–42 (2009).
  * Mackey  _et al._ (2012) M. C. Mackey, M. Tyran-Kamińska, H.-O. Walther, _et al._ , “The mathematical legacy of Andrzej Lasota,” Wiad. Mat. 48, 143 (2012).
  * Mackey (2016) M. C. Mackey, “Adventures in Poland: Having fun and doing research with Andrzej Lasota,” Mat. Appl. (Warsaw) 35, 5–32 (2016).
  * Losson and Mackey (1995a) J. Losson and M. C. Mackey, “Coupled map lattices as models of deterministic and stochastic differential delay equations,” Phys. Rev. E 52, 115–128 (1995a).
  * Lasota and Mackey (1987) A. Lasota and M. C. Mackey, “Noise and statistical periodicity,” Phys. D 28, 143–154 (1987).
  * Provatas and Mackey (1991b) N. Provatas and M. C. Mackey, “Noise-induced asymptotic periodicity in a piecewise linear map,” J. Statist. Phys. 63, 585–612 (1991b).
  * Lei and Mackey (2011) J. Lei and M. C. Mackey, “Deterministic Brownian motion generated from differential delay equations,” Phys. Rev. E 84, 041105 (2011).
  * Mackey and Tyran-Kamińska (2006) M. C. Mackey and M. Tyran-Kamińska, “Deterministic Brownian motion: The effects of perturbing a dynamical system by a chaotic semi-dynamical system,” Phys. Rep. 422, 167–222 (2006).
  * Tyran-Kamińska (2014) M. Tyran-Kamińska, “Diffusion and deterministic systems,” Math. Model. Nat. Phenom. 9, 139–150 (2014).
  * Losson  _et al._ (2020) J. Losson, M. C. Mackey, R. Taylor, and M. Tyran-Kamińska,  _Density evolution under delayed dynamics: An open problem_ (Springer Verlag, 2020).
  * Driver (1977) R. D. Driver,  _Ordinary and delay differential equations_ , Vol. 20 (Springer Science & Business Media, 1977).
  * Hopf (1952) E. Hopf, “Statistical hydromechanics and functional calculus,” J. Rat. Mech. Anal. 1, 87–123 (1952).
  * Losson and Mackey (1992) J. Losson and M. C. Mackey, “A Hopf-like equation and perturbation theory for differential delay equations,” J. Statist. Phys. 69, 1025–1046 (1992).
  * Losson and Mackey (1995b) J. Losson and M. C. Mackey, “Evolution of probability densities in stochastic coupled map lattices,” Phys. Rev. E 52, 1403 (1995b).
  * Ionescu Tulcea and Marinescu (1950) C. T. Ionescu Tulcea and G. Marinescu, “Théorie ergodique pour des classes d’opérations non complètement continues,” Ann. of Math. (2) 52, 140–147 (1950).
  * Hale and Verduyn Lunel (1993) J. K. Hale and S. M. Verduyn Lunel,  _Introduction to functional-differential equations_ , Applied Mathematical Sciences, Vol. 99 (Springer-Verlag, New York, 1993).
  * Dynkin (1965) E. B. Dynkin,  _Markov processes. Vols. I, II_ , Die Grundlehren der Mathematischen Wissenschaften, Bände 121, Vol. 122 (Academic Press Inc., Publishers, New York; Springer-Verlag, Berlin-Göttingen-Heidelberg, 1965).
  * Mohammed (1984) S. E. A. Mohammed,  _Stochastic functional differential equations_ , Research Notes in Mathematics, Vol. 99 (Pitman Advanced Publishing Program, Boston, MA, 1984).
  * Bogachev (1998) V. I. Bogachev,  _Gaussian measures_, Mathematical Surveys and Monographs, Vol. 62 (American Mathematical Society, Providence, RI, 1998).
  * Hayes (1950) N. Hayes, “Roots of the transcendental equation associated with a certain difference-differential equation,” J. London Math. Soc. (2) 1, 226–232 (1950).
  * Shepp (1966) L. A. Shepp, “Radon-Nikodým derivatives of Gaussian measures,” Ann. Math. Statist. 37, 321–354 (1966).
  * Park (1972) W. J. Park, “On the equivalence of Gaussian processes with factorable covariance functions,” Proc. Amer. Math. Soc. 32, 275–279 (1972).