# Multi-messenger astronomy in the new physics modality with GPS constellation

Arko P. Sen1 Kalia Pfeffer1 Paul Ries2 Geoffrey Blewitt1 Andrei Derevianko1
1University of Nevada, Reno, NV 89511, USA  2Jet Propulsion Laboratory,
California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA
<EMAIL_ADDRESS>
###### Abstract

We explore a novel, exotic physics, modality in multi-messenger astronomy. We
are interested in exotic fields emitted by the mergers and their direct
detection with a network of atomic clocks. We specifically focus on the
rubidium clocks onboard satellites of the Global Positioning System. Bursts of
exotic fields may be produced during the coalescence of black hole
singularities, releasing quantum gravity messengers. To be detectable such
fields must be ultralight and ultra-relativistic and we refer to them as
exotic low-mass fields (ELFs). Since such fields possess non-zero mass, the
ELF bursts lag behind the gravitational waves emitted by the very same merger.
Then the gravitational wave observatories provide a detection trigger for the
atomic clock networks searching for the feeble ELF signals. ELFs would imprint
an anti-chirp transient across the sensor network. ELFs can be detectable by
atomic clocks if they cause variations in fundamental constants. We report our
progress in the development of techniques to search for ELF bursts with clocks
onboard GPS satellites. We focus on the binary neutron star merger GW170817 of
August 17, 2017. We find an intriguing excess in the clock noise post LIGO
gravitational wave trigger. Potentially the excess noise could be explained
away by the increased solar electron flux post LIGO trigger.

##  1 Introduction

Since the initial discovery of gravitational waves (GW) by LIGO in 2015 [1],
there have been multiple detection of GWs. Most of these GWs result from
mergers of black hole binaries. However, on August 17, 2017, a new class of GW
events was discovered: the binary neutron star merger GW170817 [2]. That was
the first astrophysical source detected in both the GW and multiwavelength
electromagnetic radiation modalities and such is referred to as the multi-
messenger astronomy [3]. The messengers came from a host galaxy 40 megaparsecs
($\sim$100 million light-years) away. Ref. [4] extended the gravitational and
electromagnetic modalities of multi-messenger astronomy to exotic, beyond the
Standard Model of elementary particles, fields. We start by summarizing
salient points of the concept paper [4] and then report on our progress in the
search for such exotica in Global Positioning System (GPS) atomic clock
datastreams.

Fig. 1 summarizes the idea of multi-messenger astronomy with the addition of
exotica modality. A merger emits both the GW and an exotic low-mass fields
(ELFs) bursts. Since such fields possess non-zero mass, the burst propagates
at a group velocity $v_{g}$ smaller than the speed of light. Thereby, the ELF
pulse lags behind the GW burst. Because of the dispersion, the ELF burst tends
to spread out as it travels. More energetic components of the ELF wave-packet
travel faster, imprinting a universal anti-chirp signature in the atomic clock
data. We refer the reader to Ref. [4] for technical details and to Ref. [5]
for a more qualitative and informal discussion.

Figure 1:  An emitted ELF burst propagates with the group velocity
$v_{g}\lesssim c$ to the detector, lagging behind the GW burst. Since the more
energetic ELF components propagate faster, the arriving ELF wavepacket
exhibits a characteristic frequency anti-chirp. From Ref. [4].

As to the exotic physics modality, we focus on ultralight (yet non-zero mass)
bosonic fields as the messengers. We assume that bursts of such ELFs could be
generated by cataclysmic astrophysical events such as black hole or neutron
star mergers [6, 7]. Considering that much of the underlying physics remains
to be discovered, a pragmatic approach [4] is based on energy arguments: based
on the uncertainties for the LIGO events, $\Delta E\sim M_{\odot}c^{2}$ of
energy can be released in the form of ELFs from a black hole binary and
$\Delta E\sim 0.1\,M_{\odot}c^{2}$ in a neutron star merger. These are extreme
(compared to terrestrial experiments) amounts of energy potentially released
in the ELF modality. Atomic clocks are sensitive to gentle perturbations of
energy levels by coherent, classical waves. Due to copious numbers of emitted
ELFs this condition is indeed satisfied all the way to the sensors. Notice
that the atomic clocks sensitivities are such that we do not need require
$\Delta E\sim M_{\odot}c^{2}$ for the ELFs to be detectable. Much smaller
fractions of the solar mass would suffice for LIGO detectable events.

There is a wide variety of speculative scenarios [4] for coherent ELF
production in the mergers. One particularly intriguing scenario is due to
quantum gravity. Much of the underlying physics of coalescing singularities in
black hole mergers remains unexplored as it requires understanding of the as
yet unknown theory of quantum gravity [8]. Then the ELF burst would emerge
from the merger as a quantum gravity messenger.

Relativistic massive fields satisfy the energy-momentum dispersion relation

$\omega(k)=\sqrt{(ck)^{2}+\Omega_{c}^{2}},$ (1)

where the Compton frequency $\Omega_{c}=mc^{2}/\hbar$ depends on the ELF mass
$m$. The energy and momentum of ELF quanta are given by
$\varepsilon=\hbar\omega$ and $p=\hbar k$, where $k$ is the wave number. Since
we require a close time correlation (say delay less than a week) between the
arriving GW and ELF bursts and the vast burst travel distances, the ELFs must
be necessarily relativistic, i.e. their group velocity $v_{g}\approx c$.

A propagation of the initially emitted Gaussian ELF pulse (of central
frequency $\omega_{0}$ and of pulse duration $\tau_{0}$) leads to the solution
[4] for a field at the clock a distance $R$ away from the progenitor

$\displaystyle\phi(R,t)\approx$ $\displaystyle\frac{1}{R}\left(\frac{c\Delta
E}{2\pi^{3/2}\omega_{0}^{2}\tau}\right)^{1/2}\exp\left(-\frac{(t-t_{s})^{2}}{2\tau^{2}}\right)$
$\displaystyle\times\cos\left(\omega_{0}(t-t_{s})-\frac{\omega_{0}}{4\delta
t}(t-t_{s})^{2}\right)\,,$ (2)

where $t_{s}=t_{GW}+\delta t$ is the time of arrival of the center of the
pulse at the sensor. It is offset by $t_{GW}=R/c$, the time of arrival of the
GW messenger. $\Delta E$ is the total energy released in the ELF channel. The
signal duration $\tau$ and GW-ELF time delay $\delta t$ are related as

$\tau\approx 2\delta t\,/(\omega_{0}\tau_{0})\,.$

Note that the instantaneous ELF frequency is time-dependent,
$\omega(t)={\left(1-(t-t_{s})/(2\delta t)\right)}\omega_{0}$, exhibiting a
frequency “anti-chirp” at the sensor. The slope of the anti-chirp is given by

$d\omega/dt=-1/{\left(\tau\tau_{0}\right)}=-\omega_{0}/{\left(2\delta
t\right)}\,.$

This reflects a qualitative fact that the more energetic (or higher frequency)
components have a higher phase velocity $\omega/k$. The waveform (2) is shown
in Fig. 1.

A GPS clock sampling rate in our dataset is 1 Hz, see Sec. 4. This limits
$\omega_{0}$ from above. These frequencies fix energies
$\varepsilon_{0}=\hbar\omega_{0}$ of detectable ELFs to below
$10^{-14}\,\mathrm{eV}$. Practically, we are limited to a duration of the
pulse at the sensor to about a month; this limits $\omega_{0}$ from below. In
fact, our GPS data set, see Sec. 5.2, is effectively limited to $\sim$ a day,
requiring that $\omega_{0}\gtrsim 2\pi/\text{day}$, i.e. above $\sim
10^{-19}~{}\rm{eV}$. Moreover, it is desirable to observe at least several
oscillations, so the anti-chirp can be observed. To summarize,

$10^{-19}\,\mathrm{eV}\lesssim\varepsilon_{0}\lesssim
10^{-14}\,\mathrm{eV}\,.$ (3)

Now we know how an ELF pulse looks like at a clock. How can we detect them?
The most obvious way how the clock frequencies can be affected is through the
variation of fundamental constants, discussed in the following section.

##  2 Transient variation in fundamental constants

ELFs can couple to atomic clocks through interactions that lead to the
variation of fundamental constants. We are interested in the phenomenological
Lagrangian for the SM-ELFs interactions

$\displaystyle\mathcal{L}=\hbar
c{\left(-\sum_{f}\Gamma_{f}m_{f}c^{2}\bar{\psi}_{f}\psi_{f}+\frac{\Gamma_{\alpha}}{4}F_{\mu\nu}F^{\mu\nu}\right)}\phi^{2}\,.$
(4)

Here $\psi_{f}$ are the SM fermion fields (of masses $m_{f}$) and $F_{\mu\nu}$
is the electromagnetic Faraday tensor. $\Gamma_{X}$ are coupling constants
(with $X=\alpha,m_{f}$) quantifying the strength of the quadratic interaction
between the scalar $\phi$ and SM fields. We focus on the interactions that are
quadratic in the ELF field $\phi$ for a number of reasons, like the U(1)
invariance for complex-valued fields and the fact that the current limits on
the couplings being much weaker for the quadratic couplings than for the
linear ones.

The interaction Lagrangian (4) leads to the effective redefinition of
fundamental constants (FCs): fermion masses and the fine structure constant,

$\displaystyle
m^{\mathrm{eff}}_{f}(\textbf{r},t)=m_{f}{\left[1+\Gamma_{f}\hbar
c\phi^{2}(\textbf{r},t)\right]}\,,$ (5)
$\displaystyle\alpha^{\mathrm{eff}}(\textbf{r},t)=\alpha{\left[1+\Gamma_{\alpha}\hbar
c\phi^{2}(\textbf{r},t)\right]}\,.$ (6)

The variations in FCs modify clock frequencies and, thereby, the apparent time
as measured by the clock. In general, the sensitivity of an oscillator
resonance frequency $\nu$ to the variation in FC $X$ can be quantified by the
coefficient

$\displaystyle\kappa_{X}=\frac{\partial\,\log\nu}{\partial\log X}\,.$ (7)

In atomic clocks, the clock frequency is determined by frequencies of both the
atomic transition and the local oscillator (LO), see Sec. 3. The LOs in GPS
clocks are quartz oscillators with their fundamental mode frequencies
depending on FCs. We can introduce the LO sensitivity coefficient
$\kappa_{X}^{\mathrm{LO}}$, as in Eq. (7). Moreover, in GPS we compare the
satellite and terrestrial clocks, so, in general, we have to operate in terms
of differential sensitivity coefficients, $K_{X}$, that take sensitivities of
individual clocks into account.

As the ELFs are ultra-relativistic, their travel time from satellites to
terrestrial atomic clock is about 0.1 s, an order of magnitude smaller than
our sampling time interval, $\Delta_{t}=1$ s. Here we used the fact that the
GPS constellation diameter $D_{\rm{GPS}}\approx 5\times 10^{4}$ km. Moreover,
all the clocks would experience approximately the same value of the field if
the ELF wavelength $\lambda\approx{2\pi\hbar c}/{\varepsilon_{0}}\gg
D_{\rm{GPS}}$. For the range (3) of detectable ELF energies, the corresponding
wavelengths map into the range $10^{5}\lesssim\lambda\lesssim 10^{10}$ km.
This means that for $\varepsilon_{0}\lesssim 10^{-13}\,\mathrm{eV}$,
$\lambda\gg 10^{4}$ km and the entire GPS constellation measures nearly the
same ELF phase and can be considered as a composite yet point-like quantum
sensor. In the opposite limit of $\lambda\ll 10^{5}\,\mathrm{km}$, field
values at different satellite clocks differ and if we were to take an average
over all the clocks at a given time, the ELF signal would appear as an excess
noise.

The fractional frequency excursion due to quadratic interaction (4) reads

$\displaystyle
s(t)=\frac{\delta\nu(t)}{\nu_{\mathrm{clock}}}=\Gamma_{\mathrm{eff}}\hbar
c\phi(t)^{2}\,,$ (8)

where $\nu_{\mathrm{clock}}$ is the unperturbed or the nominal clock
frequency, $\phi(t)$ is the ELF field at the sensor (2), and
$\Gamma_{\mathrm{eff}}\equiv\sum_{X}K_{X}\Gamma_{X}$ is the effective coupling
constant expressed in terms of the differential sensitivity coefficients,
$K_{X}$.

The sensitivity coefficients $K_{X}$ depend on the internal composition and
design of the clock and data taking rate. Our GPS clock phase (or bias) data,
see Sec. 5.2, is sampled every second. These clock biases are the differences
between the satellite and the fixed reference terrestrial atomic clock times
as measured by the clocks. We will focus on ${}^{87}\mathrm{Rb}$ microwave GPS
clocks and fix an active hydrogen maser located in Kourou, French Guiana
(station ID: KOUR) as a reference clock. In the following section, we discuss
the operation of such atomic clocks.

If the characteristic time scale of the transient variation of FCs is longer
than the servo loop time constant, $K_{X}=\kappa_{X}^{\mathrm{atom}}$,
otherwise $K_{X}=\kappa_{X}^{\mathrm{LO}}$. From Ref. [9, 10, 11], we have
$\kappa_{\alpha}^{\mathrm{Rb}}\approx 4.34$, $\kappa_{m_{e}}^{\mathrm{Rb}}=2$,
$\kappa_{m_{q}}^{\mathrm{Rb}}\approx-0.069$,
$\kappa_{\alpha}^{\mathrm{H}}\approx 4$, $\kappa_{m_{e}}^{\mathrm{H}}=2$,
$\kappa_{m_{q}}^{\mathrm{H}}\approx-0.150$, $\kappa_{\alpha}^{\mathrm{LO}}=2$
and $\kappa_{m_{e}}^{\mathrm{LO}}=3/2$, where $m_{e}$ and $m_{q}$ are the
electron and average up and down quark masses respectively. As discussed in
Sec. 4, our data is acquired in the regime when
$K_{X}=\kappa_{X}^{\mathrm{atom}}$. Thereby, the effective coupling constant
for the Rb-H pairwise comparison becomes

$\displaystyle\Gamma_{\mathrm{eff}}=0.34\,\Gamma_{\alpha}+0.081\,\Gamma_{m_{q}}\,,$
(9)

where the dependence on $\Gamma_{m_{e}}$ cancels out. Notice that the ELF
signal would disappear in the data if we were to compare clocks of the same
kind such as Rb-Rb, since $\Gamma_{\mathrm{eff}}\equiv 0$; this is a powerful
test of systematic effects.

##  3 Atomic clocks of GPS network

GPS satellites host either 87Rb or 133Cs clocks. We focus on the Rb clocks
because of their abundance (24 out of total 26 clocks are Rb in our data set
for the GW170817 event) and because their frequency noise is primarily white
[12]. The 87Rb $5S_{1/2}$ electronic ground state splits into a manifold of
two, $F=1$ and $F=2$, hyperfine levels. The clock transition is between two
$M_{F}=0$ B-field insensitive Zeeman sub-levels of the hyperfine manifold with
the nominal frequency $\nu_{\mathrm{Rb}}\approx 6.835\,\mathrm{GHz}$. The
frequency of a quartz oscillator is locked to $\nu_{\mathrm{Rb}}$ through a
servo-loop operation, driving a frequency synthesizer to generate a 10.23 MHz
digital clock. This then drives the rest of the electronics generating the GPS
microwave signals which are transmitted on carrier waves precisely at
$154\times 10.23\rm{~{}MHz}=1.57542\rm{~{}GHz}$ (L1 band) and $120\times
10.23\rm{~{}MHz}=1.22760\rm{~{}GHz}$ (L2 band) [13]. (Actually, the frequency
of the digital clock is set slightly lower to allow for the mean combined
relativistic effects of time dilation and gravitational red shift, so that the
clock frequency is precisely 10.23 MHz according to an observer at the
gravitational potential at mean sea level on Earth.) The satellites broadcast
microwave signals with a carrier phase that is commensurate with the on-board
atomic clock system, through a phase locked loop frequency synthesizer. The
phases of the microwave signals are further measured by the terrestrial
receivers and through the data processing are interpreted as a measured GPS
clock time, see Sec. 4. To estimate the effect of the FC variations, we need
to gain a deeper understanding of how the GPS signal is generated, and how it
propagates to the terrestrial receivers. We start with describing the
operation of on-board Rb clocks in Sec. 3.1 and then, in Sec. 3.2, the
terrestrial H-masers.

###  3.1 Rb atomic clocks

The GPS Rb atomic clocks are vapour-cell type [14, 15]. Such clocks include
several packages which function together to ensure the frequency locking
process. The Rb atoms are a part of the physics package (PP) which is
primarily comprised of a 87Rb discharge lamp, a filter cell containing 85Rb
isotope, a resonance cell containing 87Rb, and a photo-detector [14]. The
resonance vapor cell is placed in a microwave cavity and a solenoid generating
quantizing magnetic field [15]. The PP is integrated with supporting
electronic circuits which control temperature of the discharge lamp, filter
cell and resonance cell and also control the solenoid coil’s current supply
[14].

Notice the use of two isotopes of Rb: 87Rb and 85Rb with nuclear spins $I=3/2$
and $I=5/2$ respectively. The 87Rb lamp emits on the transitions from the
$5\,^{2}\\!P_{1/2},\,5\,^{2}\\!P_{3/2}$ levels to the ground state manifold
(D1 and D2 lines). The emitted light passes through a filter cell containing
the 85Rb isotope. The fortuitous coincidence in the spectral line frequency
associated with the $F=2$ hyperfine level of 87Rb with those associated with
$F=3$ hyperfine level of 85Rb leads to the absorption of the $F=2$ spectral
lines of D1 and D2 emissions of 87Rb from the discharge lamp by 85Rb in the
filter cell [16].

The filtered light is used for quantum state preparation via optical pumping
to the $F=2$ clock state. When this $F=2$ state is subjected to a microwave
radiation, the microwaves drive the population between the $F=2$ and $F=1$
clock levels. As a result the intensity of the light coming through the filter
cell depends on the detuning of the microwave frequency from
$\nu_{\mathrm{Rb}}$ [15, 16]. The difference in light intensity is measured by
the photodetector. This is used by a servo-loop circuit to generate an error
signal which further controls the frequency of the LO that determines the
microwave frequency. The magnitude of the error signal is the difference
between the microwave frequency and $\nu_{\text{Rb}}$.

At the LO heart in GPS clocks are quartz crystals. Quartz is a piezo-electric
material and when an electrical field is applied across the crystal, it
undergoes mechanical deformations and vibrates at a fundamental frequency of
the crystal that remains stable over relatively short timescales. The
vibration modes depend on the geometry of the crystal. The relevant
sensitivity coefficients to the variation of FCs was given in Sec. 2.

The servo-loop time constant $\tau_{\mathrm{servo}}$ determines the duration
of the fastest perturbation that can be recorded by the clock. Military
operators manually adjust $\tau_{\mathrm{servo}}$ and its precise value is not
available to us at a given epoch, however it is known to be within 0.01 s to
0.1 s [17, 18]. This is sufficient to suppress the effect of occasional
frequency jumps in the quartz oscillator down to the level of $<10^{-11}$ s in
accumulated clock phase [14]. Our dataset, described in Sec. 5, is sampled
every second, giving enough time for LO frequency to be steered to the Rb atom
frequency. This explains our choice of sensitivity constants in Eq. (9). It is
worth pointing out that the electronics itself can be affected by the varying
FCs; we presume that all the derived frequencies follow the Rb atom frequency
at a sufficiently fast time scale, so that the sampled phases of the L1 and L2
GPS signals by the terrestrial GPS receivers are a faithful representation of
the behavior of the Rb atomic frequency.

###  3.2 H-maser atomic clock

The terrestrial reference clock used for satellite clocks bias data (see Sec.
5.2), is an active H-maser T4Science atomic clock housed in the station named
KOUR, located in Kourou, French Guiana [19, 20]. Below we focus on the general
principles of H-maser operations [21]. In these clocks, the quartz oscillator
is locked to the transition frequency between hyperfine levels attached to the
ground state of the hydrogen atom, $\nu_{\mathrm{H}}\approx 1.420$ GHz.

Such a clock uses a beam of atomic hydrogen, obtained by dissociating of
molecular hydrogen $\text{H}_{2}$ in a discharge chamber. The atomic beam
passes through an inhomogeneous state-selecting magnetic field, which focuses
hydrogens in the $F=1$ states into a Teflon-coated quartz storage bulb. This
“magic” coating preserves atomic coherence for $\sim 10^{6}$ consecutive
collisions with the bulb walls, critically enabling long interrogation times.

The bulb is placed within a radiofrequency (RF) cavity tuned to the transition
frequency, $\nu_{\mathrm{H}}$ between $F=1,M_{F}=0$ and $F=0,M_{F}=0$ states
of the H-atom [22, 23]. As the H-atoms undergo the transition from the $F=1$
state down to the $F=0$ state, the emitted microwave photons stimulate
transitions in other atoms [24], hence the maser qualifier in these clocks.
Cascade of stimulated transitions produce a highly coherent oscillation in the
RF cavity. This results in a microwave signal to which the quartz oscillator
is phase-locked.

The RF cavity is enveloped by magnetic shields to attenuate the ambient
fields. An inner solenoid coil creates a uniform magnetic field to separate
the different magnetic sublevels. There are two types of interactions that
occur in the RF cavity: the atomic relaxation process that affects the
population in $F=1$ state and the process that leads to the decay of atomic
phase coherence that affect the frequency stability [23]. The major purpose of
the bulb is to spatially confine the atoms, enabling interaction with the in-
phase components of the RF magnetic field in the cavity [23].

For longer sustenance of the oscillation either a sufficiently high H-atom
density is to be maintained or the RF cavity is to be interrogated with
external RF radiation at $\nu_{\mathrm{H}}$. Atomic clocks based on the former
process is termed as active H-maser clocks, whereas clocks based on the latter
are passive H-maser clocks. In active H-masers the cavity undergoes self-
sustained oscillation. Both types of H-masers eventually drive a quartz
oscillator which outputs a reference signal with a high stability and a low
drift rate by control voltage adjustment [25]. The principal circuit of
H-maser atomic clock has a crystal oscillator, an isolation amplifier circuit,
an up and down-conversion circuit, a digital servo, and a frequency
synthesizer circuit. There is an auxiliary circuit system that includes a
constant temperature and a high-voltage source circuit, and a constant current
source [25]. In passive H-maser clocks the interrogation microwave field is
produced by an auxiliary frequency generator. H-masers have a good short-term
stability, $\sigma_{y}{\left(\tau=1~{}\mathrm{s}\right)}<10^{-12}$, reaching a
noise floor of $\approx 10^{-15}$ after about an hour of averaging [24].

The cavity frequency is stabilized by a servo loop. The value of the time
constant of the servo-loop response is in the order of 1 s [26]. As with the
Rb clocks, the response to a sudden change in FCs is that of a quartz for
time-scales much shorter than the servo-loop time constant. For time scales
$\gg$ than the servo-loop time constant, it is the response of the hydrogen
hyperfine splitting that matters. Once again, we ignore the complicated
response of electronic circuits to the FC variations, as we believe, such a
response is only important for very short time scales. These considerations
govern our choice of sensitivity constants to FC variations in Sec. 2.

##  4 GPS data processing

As discussed in Sec. 3, GPS satellites broadcast microwave signals with a
carrier phase that is commensurate with the on-board atomic clock system,
through a phase locked loop frequency synthesizer. The carrier phase is
sampled at regular intervals by a global terrestrial network of GPS receiver
stations. Each GPS receiver measures the beat phase between the incoming
carrier phase signal mixed with an internal replica signal, which is locked in
phase to the receiver’s internal quartz oscillator. Thus the beat phase
measures the difference between time of transmission according to each GPS
satellite clock, and the time of reception according to each receiver clock.

The Jet Propulsion Laboratory (JPL) routinely analyzes carrier phase data
every 30 seconds from the global GPS network [27]. For the special purpose of
this ELF search, JPL processed data at 1Hz across multiple days from an
identical global 30-station network. Data processing includes modeling the
time of flight of the carrier signal traveling at the phase velocity of light,
accounting for relativistic Shapiro delay, rotational wind-up of circularly
polarized phase, atmospheric refraction, and deterministic station motion
caused by tidal and surface loading deformation. For each station, a linear
combination of carrier phases at two different carrier frequencies is computed
such that it eliminates the effect of ionospheric refraction $n_{\rm ion}$
assuming that $(n_{\rm ion}-1)\propto 1/\omega^{2}$ [28]. Note that this
elimination applies even if fundamental constants are perturbed. The linear
combination conserves the satellite clock variations, which for GPS is
designed to be identical on carrier signals of both frequencies. Automatic
data editing algorithms allow for arbitrary variation of clock phase, and so
unusual clock behavior is not a determining factor in data outlier deletion
[29].

Parameter biases in the model are simultaneously estimated using a square-root
estimation filter. Those parameters include all the satellite and receiver
clock biases every second, with the exception of one high quality receiver
clock which defines the reference coordinate time. The receiver reference
clock is selected such that it is synchronized to a stable H-maser frequency
standard. Other biases that are important to estimate but are not explicitly
used in our subsequent analysis include station coordinates, satellite orbit
parameters, atmospheric delays and their variation, and carrier phase biases
associated with an ambiguity in the integer number of cycles, which are
subsequently determined precisely and fixed [28, 30].

For the neutral atmosphere, the signal delay is non-dispersive at microwave
frequencies, thus it is independent of frequency variations in the carrier
signal that may be caused by FC variations. However, in principle there could
be a contribution of FC variations to delay arising from changes in refractive
index. Variations in refractive index due to weather are accounted for using
the GPT2 model [31], which was based on fits to numerical weather models of
the European Center for Medium-Range Weather Forecasts.

Any residual variations in refractive index, such those caused by lack of
spatio-temporal resolution in the GPT2 model, are estimated as part of the
multi-parameter bias estimation by the square root information filter.
Specifically, the variations are estimated as a random walk in the zenith
delay and in two gradient parameters for each station. These parameters are
separable from satellite clock variations owing to the so-called “mapping
function,” that is the modeled variation of atmospheric slant thickness with
elevation angle of the satellite [31]. Such estimation should also absorb
variation in refractive index caused by the FC variations we are seeking,
leaving the imprint of clock frequency variations intact. The assumption is
that any refractive index variations caused by ELFs do not greatly exceed the
level of variations expected from weather.

The resulting data set used in our analysis is a 1-second time series of
atomic clock biases from each of the GPS satellites that have Rb clocks. All
of these biases are relative to the clock of the selected reference station.
We also have access to receiver clock biases, however most GPS stations only
use quartz oscillators without an external atomic frequency standard, and so
are of limited use. Other useful metadata include satellite positions and
velocities.

One aspect of the GPS system that we aim to exploit is that, in the co-
rotating Earth-fixed reference frame of the global GPS station network, the
satellite positions are designed to repeat approximately every sidereal day
[32, 33]. The precise repeat time must account for the precession of each
satellite’s orbit, caused by Earth’s oblateness, and the orbit’s eccentricity,
inclination, and semi-major axis. This precession is slightly different for
each satellite, resulting in a repeat period of that varies in the range 1 day
minus $245\pm 5$s [33]. This repeat of the orbit geometry causes systematic
errors such as ground multipath reflections to also repeat, hence errors in
the carrier phase measurement tend also to repeat. Therefore, there is the
potential to enhance signal to noise by comparing measured carrier phase time
series from one repeat cycle to the next, as in the common mode error
subtraction techniques. Such methods have proven to greatly reduce scatter in
geodetic GPS positioning every second, with accuracy at the few mm level [33],
therefore signal to noise ratio in clock bias data should similarly be
improved.

##  5 Search for ELF signal and preliminary results

###  5.1 Effect of ELFs on the clock output

The Rb-locked LO frequency is measured by a narrow band circuit whose signal
can be related to a sinusoidal wave [34]. The output voltage of such a circuit
can be expressed as

$\displaystyle
V(t)={\left[V_{0}+\varepsilon(t)\right]}\sin{\left[2\pi\nu_{\mathrm{clock}}t+\varphi(t)\right]}\,,$
(10)

where $V_{0}$ and $\nu_{\mathrm{clock}}$ are nominal peak voltage amplitude
and the clock frequency respectively and $\varepsilon(t)$ and $\varphi(t)$ are
the deviation from the nominal peak voltage and nominal clock phase [34, 35].
Usually the oscillator frequency remains constant for a short time interval
but it can drift over a relatively larger period of time.

The instantaneous frequency is conventionally defined as the first derivative
of the phase, which is

$\displaystyle\nu(t)=\nu_{\mathrm{clock}}+\frac{1}{2\pi}\frac{d\varphi(t)}{dt}.$
(11)

The extra phase term can be due to either technical noise or due to the ELF
signal. Given a sampling time interval $\Delta_{t}$, the off-set phase is
explicitly

$\displaystyle\varphi(t_{j})$
$\displaystyle=2\pi{\left[\int_{t_{j-1}}^{t_{j}}\nu(t^{\prime})dt^{\prime}-\nu_{\mathrm{clock}}\int_{t_{j-1}}^{t_{j}}dt^{\prime}\right]}\,$
$\displaystyle\approx
2\pi{\left[\nu(t_{j})-\nu_{\mathrm{clock}}\right]}\Delta_{t}\,,$ (12)

where $t_{j}=j\Delta_{t}$ ; $j=\overline{1,N_{t}}$ and $N_{t}$ being to total
number of data points.

Note that the deviation from the nominal clock frequency is
$\delta\nu(t_{j})=\nu(t_{j})-\nu_{\mathrm{clock}}$. Then the fractional clock
frequency excursions are

$\displaystyle\frac{\delta\nu(t_{j})}{\nu_{\mathrm{clock}}}\approx\frac{1}{2\pi\Delta_{t}}\frac{\varphi(t_{j})}{\nu_{\mathrm{clock}}}\,,$
(13)

to be matched to the sought ELF waveform (8).

###  5.2 Data analysis and results

In our GPS.ELF group, we pursue the analysis of data from the GPS atomic
clocks. One can think of the GPS constellation as the largest human-built
$\sim 50,000\,\text{km}$-aperture sensor array. Combined with other satellite
positioning constellations and terrestrial clocks, the total number of sensors
$N_{c}\sim 100$. Employing networks is crucial for distinguishing ELF signals
from technical noise. Furthermore, by having baselines with the diameter of
the Earth or larger, one can resolve the sky position of the ELF source. This
is a critical feature for multi-messenger astronomy that enables correlation
with other observations of the progenitor. However, an ELF feature (say
leading edge) propagation time across the GPS constellation is
$0.2\,\text{s}$. The publicly available GPS data is sampled every
$\Delta_{t}=30\,\text{s}$. Our GPS.ELF group has developed techniques to
generate more frequently sampled, $\Delta_{t}=1\,\text{s}$, GPS datastreams.
Still the limited time resolution makes tracking the leading edge of the ELF
pulse across the GPS network challenging. Nevertheless, clock networks can
still act collectively, gaining $\sqrt{N_{c}}$ in sensitivity and vetoing
signals that do not affect all the sensors in the network.

We work with the GPS clock bias data for the day of GW170817 event detection
and the day before. The satellite clock biases are relative to a fixed
reference terrestrial clock. This is a non-stationary time series, dominated
by random walk noise. This discrete clock bias data stream $d^{(0)}_{a,j}$,
sampled every $\Delta_{t}=1$ s, can be expressed, for the $a^{\text{th}}$
satellite clock, in units of time as

$\displaystyle{d^{(0)}_{a,j}}=\frac{\varphi_{a}(t_{j})}{\omega_{a}}-\frac{\varphi_{r}(t_{j})}{\omega_{r_{0}}}\,,$
(14)

where the index $j=\overline{1,N_{t}}$, enumerates the epochs. $N_{t}$ is the
number of points sampled during the orbital period, as discussed below. The
frequencies $\omega_{a}$ and $\omega_{r_{0}}$ are the nominal clock angular
frequency of a satellite $a$ and the reference terrestrial clocks. The time
dependent discrete phases of the satellite and reference clocks are denoted by
$\varphi_{a}(t_{j})$ and $\varphi_{r}(t_{j})$ respectively. The reference
terrestrial clock for this work is a hydrogen maser at the KOUR station
located in Kourou, French Guiana, see Sec. 3.2.

Although the proposal paper [4] has used the formalism of excess noise
statistic and time-frequency decomposition, here we carry out a simplified
analysis of the excess noise to also include a possibility of high-
frequency/short wavelength ELFs, see Sec. 1.

As discussed in Sec. 4, the orbital geometry of the satellites repeats
relative to observer on Earth, leading to repeating multipath errors. To
account for the multipath errors, the GPS satellite orbital period based on
the Earth Centered Earth Fixed (ECEF) reference frame is set approximately to
24 hours minus 247 s = 86,153 s. The sampling interval of 1 s thus leads to
the total number of points per day $N_{t}=86,153$. To whiten the data, we
employ the first order differencing of the discrete clock bias data, the so-
called discretized pseudo-frequency [36],

$\displaystyle{d^{(1)}_{a,j}}={d^{(0)}_{a,j+1}-d^{(0)}_{a,j}}\,,$ (15)

where $a$ indexes individual datastreams and $j$ spans over sampling grid. The
differencing reduces the total number of discrete pseudo-frequency points to
$N\equiv N_{t}-1$. We remove the frequency drift and the offset for each clock
$a$

$\displaystyle u^{(1)}_{a,j}=d^{(1)}_{a,j}-f_{a,j}\,,$ (16)

by fitting the ${d^{(1)}_{a,j}}$ data to the straight line of the form,
$f_{a,j}=m_{a}t_{j}+b_{a}$, with $m_{a}$ and $b_{a}$ being fitting
coefficients. The weighted average of $u^{(1)}_{a,j}$ over the network of
$N_{c}$ satellite clocks is

$\displaystyle\bar{d}^{(1)}_{j}=\frac{\sum_{a}^{N_{c}}{u}^{(1)}_{a,j}/\sigma^{2}_{a}}{\sum_{a}^{N_{c}}1/\sigma^{2}_{a}}\,,$
(17)

where

$\displaystyle\sigma_{a}^{2}=\frac{1}{N-1}\sum_{j=1}^{N}{\left({u}^{(1)}_{a,j}-{\left\langle{u}^{(1)}_{a}\right\rangle}\right)}^{2}\,,$
(18)
$\displaystyle{\left\langle{u}^{(1)}_{a}\right\rangle}=\frac{1}{N}\sum_{j=1}^{N}{u}^{(1)}_{a,j}~{}.$
(19)

Here $\sigma_{a}^{2}$ and ${\left\langle{u}^{(1)}_{a}\right\rangle}$ are the
variance and mean of ${u}^{(1)}_{a,j}$ respectively for $N_{c}=24$ Rb atomic
clocks.

Further, the network weighted average pseudo-frequency as calculated using Eq.
(17), is partitioned into $M$ windows with each window consisting of
$N_{w}=N/M$ points. We define $\bar{d}^{(1)}_{w,k}\equiv\bar{d}^{(1)}_{j}$,
where $w=\overline{1,M}$ enumerates the windows and index
$k=\overline{1,N_{w}}$ here spans points in the window $w$, i.e.
$j={\left(w-1\right)}N_{w}+k$. The standard deviation in a window is
calculated as

$\displaystyle\sigma_{w}=\sqrt{\frac{1}{N_{w}-1}\sum_{k=1}^{N_{w}}{\left(\bar{d}^{(1)}_{w,k}-{\left\langle\bar{d}_{w}^{(1)}\right\rangle}\right)}^{2}}\,,$
(20)

with

$\displaystyle{\left\langle\bar{d}_{w}^{(1)}\right\rangle}=\frac{1}{N_{w}}\sum_{k=1}^{N_{w}}\bar{d}^{(1)}_{w,k}\,.$
(21)

Figure 2: Standard deviation in a window, $\sigma_{w}$, of pseudo-frequency
for the GPS network average of Rb clocks for August 16 and 17, 2017. GW
trigger detected at 45,682 s from the start of the day of August 17, 2017
(GPST), is shown with the black line. In computing standard deviation, we
partitioned each day dataset into $M=712$ windows, each containing $N_{w}=121$
data points.

In Fig. 2, the standard deviation $\sigma_{w}$ is plotted with respect to time
$t_{w}$ associated with the beginning of the $w^{\text{th}}$ window,
$t_{w}=(w-1)N_{w}\Delta_{t}$. The time is based on the GPS Time (GPST) scale
offset at the beginning of each GPST day. Note that the difference between
GPST and Coordinated Universal Time (UTC) time scale is 18 s for the year
2017. There are two traces: for the day before the GW170817 LIGO trigger
(August 16, 2017) and for the day of the event (August 17, 2017). Comparing
the two days, we see that the clock network exhibits excess noise level post
the GW trigger shown by a black vertical line. In particular, we observe
several spikes in the standard deviation of $\bar{d}^{(1)}_{j}$ occuring after
the GW LIGO trigger (on August 17, 2017). The spike with the highest magnitude
is delayed by 17,600 s ($\approx 4.9$ hours), from the GW LIGO trigger.
Although we see the spike with a high magnitude in both the days’ data, the
spike post GW trigger is larger by a factor of about 2.3. This factor was
computed by taking the average standard deviation for both days (at
$\approx{4.7}\times 10^{-12}\,\mathrm{s}$) as the floor. Generically, one
would expect both spikes of approximately the same magnitude due to repeating
multipath errors in geodetic receivers. It is plausible that this extra noise
is due to the putative ELF signal. To quantify this statement additional
statistical analysis is required, in particular by finding the optimal size of
the window $N_{w}$ and characterizing the clock statistics. Notice, however,
that we found that the solar electron flux was elevated during this period,
see Fig. 3, perhaps negating the ELF hypothesis. This is discussed in the
following section.

Figure 3: The average electron flux depicts a period of elevated electron flux
that occurs after the GW event, GW170817. The GW LIGO event occurred at 45,682
s on August 17, 2017 (GPST) and is marked with a black vertical line. The
elevated electron flux continues onto the following two days, maintaining the
same order of magnitude.

##  6 Systematic effects due to solar activity

The expulsion of solar energetic particles into the near-Earth space
environment during periods of enhanced solar activity can be responsible for
compromising GPS radio communication and positioning accuracy. Enhanced solar
activity is defined as periods of time in which the solar corona exhibits
increased magnetic field deformation. The acceleration of energetic particles
from the solar corona, with the total energy release of about $10^{24}$ to
$10^{32}$ ergs, typically characterizes a solar flare. Coronal mass ejections
(CME) often accompany solar flares, where charged plasma particles with an
average velocity ranging from $500$ to $800$ km/s are ejected into the
environment [37].

Energetic particles emitted during solar activity events are introduced to
Earth’s magnetosphere and can precipitate to the ionosphere through
acceleration along magnetic field lines. The energetic electrons and protons
can also be trapped within the Van Allen radiation belts, accessible to GPS
satellites. Exposure to high energetic particle fluxes within the radiation
belts can induce malfunctions to the satellites’ electric power system [38].
The precipitation of energetic particles to the ionosphere can introduce
ionization density variations, resulting in fluctuations to the amplitude and
phase of transmitted signals. This phenomenon is known as ionospheric
scintillation. GPS transmitted signals are vulnerable to strong ionospheric
scintillation associated with enhanced electron content [39].

The Los Alamos National Laboratory has released energetic particle data
recorded by the Combined X-Ray Dosimeter (CXD) aboard 21 operational GPS
satellites. The CXD is equipped with high-energy and low-energy particle
sensors that can record electron fluxes with an energy of 1 MeV to above 5 MeV
and proton fluxes with an energy of 16 MeV to above 60 MeV [40]. An increased
sunspot number associated with enhanced solar activity can correspond to
elevated energetic particles introduced to the GPS orbit environment.
Increased energetic particle exposure can be represented by spikes in electron
and proton flux measurements. The daily sunspot number for August of 2017
indicates elevated values occurring at the end of the month.

We examined the network average proton and electron fluxes from the CXD aboard
17 GPS satellites for 16 August and 17 August, 2017, c.f. Fig. 2. The plot of
the average electron flux is shown in Fig. 3. We observe a period of elevated
electron flux that occurs after the GW event, GW170817. The proton flux
exhibited no obvious spike or irregularity before this GW event. The spike in
the network weighted average of the pseudo-frequency depicted in Fig. 2 could
be related to elevated energetic particle exposure. Notice that the solar
electron flux leads to a charge build-up on a satellite that could randomly
discharge after reaching a certain level [41], which could potentially cause
frequency jumps in the quartz oscillator [14, 18], modifying GPS microwave
signals. This elevated electron flux starts around the time of the LIGO
trigger. One of the intriguing explanations might be that the GWs perturb the
Sun causing the solar activity.

##  7 Conclusions

The ELF proposal [4] has introduced a novel, exotic physics, modality in
multi-messenger astronomy. In this report we outlined some practical
considerations for the ELF search with the network of atomic clocks on board
GPS satellites. We have generated 1 Hz clock bias data specifically for this
analysis. Our search method focused on the excess noise in the clock network
weighted average pseudo-frequency. There is an intriguing excess in the clock
noise post LIGO gravitational wave trigger for GW170817 binary neutron star
merger of August 17, 2017. Potentially it could be explained away by the
increased solar electron flux during the observed spike in the clock noise. It
is also intriguing that solar activity is greatly enhanced around the time and
after the LIGO trigger.

## Acknowledgments

This work was supported in part by the U.S. National Science Foundation grant
PHY-2207546, by the Heising-Simons Foundation, and by the NASA Grant
80NSSC23M0104.

## References

## References

  * [1] B Abbott et al (2016) Physical Review Letters 116 061102 ISSN 0031-9007 URL https://link.aps.org/doi/10.1103/PhysRevLett.116.061102
  * [2] B Abbott et al (2017) Phys. Rev. Lett. 119 161101 ISSN 0031-9007 (Preprint 1710.05836) URL https://link.aps.org/doi/10.1103/PhysRevLett.119.161101
  * [3] B Abbott et al (2017) The Astrophysical Journal Letters 848 L13 URL https://dx.doi.org/10.3847/2041-8213/aa920c
  * [4] C Dailey et al (2021) Nature Astronomy 5 150–158 ISSN 2397-3366 (Preprint 2002.04352) URL http://www.nature.com/articles/s41550-020-01242-7
  * [5] Derevianko A (2023) (Preprint 2305.17138) URL http://arxiv.org/abs/2305.17138
  * [6] Bini D, Geralico A and Ortolan A (2017) Phys. Rev. D 95(10) 104044 URL https://link.aps.org/doi/10.1103/PhysRevD.95.104044
  * [7] D Baumann et al (2019) Phys. Rev. D 99 044001 ISSN 2470-0010 URL https://doi.org/10.1103/PhysRevD.99.044001 https://link.aps.org/doi/10.1103/PhysRevD.99.044001
  * [8] Loeb A (2018) (Preprint 1805.05865) URL https://arxiv.org/abs/1805.05865
  * [9] W Campbell et al 2021 Phys. Rev. Lett. 126(7) 071301 URL https://link.aps.org/doi/10.1103/PhysRevLett.126.071301
  * [10] V Flambaum et al 2004 Phys. Rev. D 69(11) 115006 URL https://link.aps.org/doi/10.1103/PhysRevD.69.115006
  * [11] B Roberts et al 2018 Phys. Rev. D 97 083009 ISSN 2470-0010 URL http://arxiv.org/abs/1803.10264 https://link.aps.org/doi/10.1103/PhysRevD.97.083009
  * [12] Danzy F and Riley W 1989 Stability test results for gps rubidium clocks Proceedings of the Annual Precise Time and Time Interval (PTTI) applications and Planning Meeting (21st), Held in Redondo Beach, California on November 28-30, 1989 URL https://apps.dtic.mil/sti/pdfs/ADA224769.pdf
  * [13] Blewitt G (2015) GPS and space-based geodetic methods Treatise on Geophysics (Oxford: Elsevier) pp 307–338 URL http://linkinghub.elsevier.com/retrieve/pii/B9780444538024000609
  * [14] Riley W (1981) A rubidium clock for GPS Proceedings of the 13th Annual Precise Time and Time Interval Systems and Applications Meeting URL https://apps.dtic.mil/sti/pdfs/ADA494259.pdf
  * [15] Camparo J 2007 Physics Today 60 33–39 ISSN 0031-9228 (Preprint https://pubs.aip.org/physicstoday/article-pdf/60/11/33/16696572/33_1_online.pdf) URL https://doi.org/10.1063/1.2812121
  * [16] Y Guo et al (2022) AIP Advances 12 045216 ISSN 2158-3226 URL https://doi.org/10.1063/5.0086523
  * [17] Griggs E, Kursinski E R and Akos D 2015 Radio Science 50 813–826 (Preprint https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1002/2015RS005667) URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015RS005667
  * [18] Dupuis et al R 2008 Frequency Control Symposium, 2008 IEEE International 655–660
  * [19] International GNSS Service (IGS) File Access URL https://files.igs.org/pub/station/log/kour_20230405.log
  * [20] International Terrestrial Reference Frame, (ITRF) URL https://itrf.ign.fr/en/station/KOUR-97301M210
  * [21] Major F G 2007 The Quantum Beat:Principles and Applications of Atomic Clocks; 2nd ed. (Springer-Verlag) URL "https://link.springer.com/book/10.1007/978-1-4757-2923-8"
  * [22] Kleppner D, Goldenberg H M and Ramsey N F 1962 Phys. Rev. 126(2) 603–615 URL https://link.aps.org/doi/10.1103/PhysRev.126.603
  * [23] Vessot R F C 2005 Metrologia 42 S80 URL https://dx.doi.org/10.1088/0026-1394/42/3/S09
  * [24] Lombardi M 2002 Time and Frequency pp 783–801 URL https://www.researchgate.net/publication/269420574_Time_and_Frequency
  * [25] Hu W, Shuai T, Xie Y, Chen P, Pei Y, Zhao Y and Wang R 2023 Sensors 23 ISSN 1424-8220 URL https://www.mdpi.com/1424-8220/23/22/9202
  * [26] Audoin, C 1981 Rev. Phys. Appl. (Paris) 16 125–130 URL https://doi.org/10.1051/rphysap:01981001603012500, https://hal.science/jpa-00244902/document
  * [27] W Bertiger et al 2020 Advances in Space Research 66 469–489 ISSN 0273-1177 URL https://www.sciencedirect.com/science/article/pii/S0273117720302532
  * [28] Blewitt G (1989) Journal of Geophysical Research: Solid Earth 94 10187–10203 (Preprint https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JB094iB08p10187) URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB094iB08p10187
  * [29] Blewitt G (1990) Geophysical Research Letters 17 199–202 (Preprint https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/GL017i003p00199) URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/GL017i003p00199
  * [30] Bertiger et al W 2010 Journal of Geodesy 84 327–337 URL https://doi.org/10.1007/s00190-010-0371-9
  * [31] J Böhm et al 2015 GPS Solutions 19 433–441 URL https://doi.org/10.1007/s10291-014-0403-7
  * [32] D Agnew et al 2007 GPS Solutions 11 71–76 URL https://doi.org/10.1007/s10291-006-0038-4
  * [33] K Larson et al 2007 Journal of Geophysical Research: Solid Earth 112 (Preprint https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2006JB004367) URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006JB004367
  * [34] D Howe et al 1981 Properties of Signal Sources and measurement methods Proceedings of the 35th Annual Symposium on Frequency Control URL https://tf.nist.gov/general/pdf/554.pdf
  * [35] Riley W and Howe D 2008 Handbook of frequency stability analysis URL https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=50505
  * [36] B Roberts et al (2017) Nature Comm. 8 1195 ISSN 2041-1723 URL http://www.nature.com/articles/s41467-017-01440-4
  * [37] Duffin R T, White S M, Ray P S and Kaiser M L 2015 Journal of Physics: Conference Series 642 012006 URL https://dx.doi.org/10.1088/1742-6596/642/1/012006
  * [38] Chen Y, Carver M R, Morley S K and Hoover A S 2021 Determining ionizing doses in medium earth orbits using long-term gps particle measurements 2021 IEEE Aerospace Conference (50100) pp 1–21 URL https://ieeexplore.ieee.org/document/9438516
  * [39] Kintner P M, Ledvina B M and de Paula E R 2007 Space Weather 5 (Preprint https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2006SW000260) URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2006SW000260
  * [40] J Distel et al 1999 The combined X-ray dosimeter (CXD) on GPS Block IIR satellites Tech. rep. Tech. Rep. LA-UR-99 URL https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-99-2280
  * [41] Frederickson A, Mullen E, Brautigam D, Kerns K, Robinson P and Holman E 1991 IEEE Transactions on Nuclear Science 38 1614–1621 URL https://ieeexplore.ieee.org/document/124153