# Enhanced spin coherence while displacing electron in a 2D array of quantum
dots

Pierre-André Mortemousque  Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut
Néel, 38000 Grenoble, France  Univ. Grenoble Alpes, CEA, Leti, F-38000
Grenoble, France  Baptiste Jadot  Univ. Grenoble Alpes, CNRS, Grenoble INP,
Institut Néel, 38000 Grenoble, France  Univ. Grenoble Alpes, CEA, Leti,
F-38000 Grenoble, France  Emmanuel Chanrion  Univ. Grenoble Alpes, CNRS,
Grenoble INP, Institut Néel, 38000 Grenoble, France  Vivien Thiney  Univ.
Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble, France
Christopher Bäuerle  Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel,
38000 Grenoble, France  Arne Ludwig  Lehrstuhl für Angewandte
Festkörperphysik, Ruhr-Universität Bochum, Universitätsstraße 150, D-44780
Bochum, Germany  Andreas D. Wieck  Lehrstuhl für Angewandte Festkörperphysik,
Ruhr-Universität Bochum, Universitätsstraße 150, D-44780 Bochum, Germany
Matias Urdampilleta  Univ. Grenoble Alpes, CNRS, Grenoble INP, Institut Néel,
38000 Grenoble, France  Tristan Meunier  Univ. Grenoble Alpes, CNRS, Grenoble
INP, Institut Néel, 38000 Grenoble, France

###### Abstract

The ability to shuttle coherently individual electron spins in arrays of
quantum dots is a key procedure for the development of scalable quantum
information platforms. It allows the use of sparsely populated electron spin
arrays, envisioned to efficiently tackle the one- and two-qubit gate
challenges. When the electrons are displaced in an array, they are submitted
to site-dependent environment interactions such as hyperfine coupling with
substrate nuclear spins. Here, we demonstrate that the electron multi-
directional displacement in a $3\times 3$ array of tunnel coupled quantum dots
enhances the spin coherence time via the motional narrowing phenomenon. More
specifically, up to 10 configurations are explored by the electrons to study
the impact of the displacement on spin dynamics. An increase of the coherence
time by a factor up to 10 is observed in case of fast and repetitive
displacement. The physical mechanism responsible for the loss of coherence
induced by displacement is quantitatively captured by a simple model and its
implications on spin coherence properties during the electron displacement are
discussed in the context of large-scale quantum circuits.

The control over the flow of electrons in semiconductor circuits is core to
the success of micro and nanoelectronics. In their quantum counterparts,
similar processes are investigated at the single-particle level, to preserve
the fragile quantum information stored in individual electron spins [1].
Indeed, displacing these electron spins coherently opens possibilities to
convey on-chip quantum information and increase the connectivity of spin-based
semiconductor quantum circuits [2, 3, 4, 5, 6, 7]. Several strategies have
been recently demonstrated to preserve spin coherence while shuttling charged
particles. A first example consists of the confinement of electrons in moving
quantum dots [8, 9]. Another protocol is based on repetitive coherent spin
tunnelling between adjacent dots. Such a strategy is particularly appealing
for architecture based on large 2D arrays of quantum dots [10, 11] and has
been demonstrated for linear [12, 13, 14] and circular [15] arrays. Therefore,
understanding the mechanisms at play affecting the fidelity of the
displacement in 2D arrays is an important task to optimize this quantum
information conveyer procedure and increase as much as possible the distance
over which the spin transfer is coherent.

For an electron spin in a quantum dot, the loss of coherence is driven by spin
dephasing mechanisms that arise from intrinsic properties of the
semiconducting nanostructures like 1/f noise coupled to spin-orbit interaction
[16] or hyperfine interaction [17]. In both cases, the electron spins are
experiencing slow effective magnetic field fluctuations inducing uncertainty
on the Larmor precession of individual electron spins. Therefore, refocusing
techniques are extremely efficient on spin systems [18, 19, 20, 21]. With this
property in mind, displacing electron spins on fast timescales is expected to
enhance the spin coherence times in virtue of the averaging of the magnetic
field fluctuations, as demonstrated in NMR via a phenomenon called motional
narrowing effect [22, 23]. The unprecedented level of control recently
demonstrated in a 2D array of quantum dots [24] permits to displace the
electron through an important set of possible dot configurations on a
timescale faster than the coherence time. This array turns out to be a
relevant platform to investigate and reveal the consequences of electron
displacement on spin coherence properties.

Here we study the coherence of individual electron spins when they are
displaced via tunnelling through a 2D array of five quantum dots defined in an
AlGaAs heterostructure. In this displacement regime with an electron speed
below $100$ $\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, the hyperfine interaction
with the nuclei of the semiconductor nanostructure is responsible for the loss
of coherence [22]. This hyperfine interaction results in an effective
fluctuating magnetic field, called the Overhauser field, with a distribution
imposed by the number of nuclei the electron spins are interacting with. The
larger the number of nuclei, the longer is the coherence time, with a typical
square-root dependence [25]. During the electron displacement, the system
visits a number of dot configurations imposed by the gate voltage sequence. We
relate to each dot configuration an effective local magnetic field. Two
principal parameters are influencing the coherence time of the displaced
electron spins: (i) the number of effective local magnetic field
configurations explored which can be varied by changing the number of
accessible quantum dots in the array, and (ii) the typical interaction time
fixed by the time spent in each charge configuration. The demonstrated control
of the dot system enables us to vary precisely both parameters and to study
their impacts on the coherent spin transfer.

## Electron shuttling

Figure 1:  Sample and spin loading and readout sequence. a, Electron
micrograph of a sample similar to the one used in this work. The nine circles
indicate the $3\times 3$ array of quantum dots (QDs) arising from the gate-
induced potential landscape (see text). The five dots L, B, R, T, C in a
crossbar configuration (solid circles) are used for the coherent electron spin
shuttling in 2D. The square lattice of QDs is electrostatically defined thanks
to the red gate voltages, and its outer edges by the green, light blue and
blue gate voltages. The couplings of the array to electron reservoirs are
achieved through the corner QDs TL, BL, BR, and TR and by controlling the
light blue and blue gates. The purple gates are used to define four quantum
point contacts operated as local electrometers, whose conductances set the
measured currents $i_{\mathrm{TL,BL,BR,TR}}$. Scale bar (white) is $200$
$\mathrm{nm}$. b,c, Schematics of the sequence used to load a singlet spin in
C (b), and to read out the final spin state after manipulation (c).

The sample [24], shown in Fig. 1a, is an array of $3\times 3$ laterally
defined quantum dots (Methods). The dots are labelled according to their
position in the array (e.g. TL stands for Top Left dot). The system is
operated in a mode where the three corner dots TR, BR, and BL are
intentionally inaccessible by the electrons (tuned at a higher potential than
the crossbar dots L, B, R, T, C). In this configuration, it is possible to
load and isolate two electrons from the reservoir in the TL dot and to
transfer them into the dot C by increasing the $V_{\mathrm{L,B,R,T}}$ voltages
before rising the potential of the TL dot (Fig. 1b). First, we prove that all
the 15 expected charge states for two electrons in five dots ($\binom{5}{1}=5$
where both electrons are in the same dot, $\binom{5}{2}=10$ where electrons
are separated in different dots) are accessible by tuning the 4 linear
combinations $\delta V^{\pm}_{\mathrm{X,Y}}$ of the $V_{\mathrm{L,B,R,T}}$
voltages (Methods). Figure 2a-c shows charge stability diagrams of two
electrons recorded as functions of $\delta V^{-}_{\mathrm{X,Y}}$, which act
like in-plane electric dipole gates, for different values of $\delta
V^{+}_{\mathrm{X,Y}}$. By adjusting simultaneously $\delta
V^{+}_{\mathrm{X,Y}}$, it is possible to change the relative potential of the
dot C so that the central region of the diagrams either contains two electrons
(Fig. 2a), or no electron revealing the regions where the electrons can be
separated in (T, B) (Fig. 2b) and (L, R) (Fig. 2c). Figure 2d-f shows
simulations (Methods) of the stability diagrams shown in Fig. 2a-c,
respectively, in which the identification of the charge distributions is
confirmed. In our precedent report on the coherent control of electron spin in
this structure [24], we have shown that all charge states exhibited in the
stability diagrams of Fig. 2 can be accessed while setting sufficiently high
tunnel couplings for the coherent spin transfers between dots. Up to 10 charge
configurations where the electrons are separated in different dots are
explored, and multi-directional and complex one- and two-electron
displacements are performed.

Figure 2:  Two-electron stability diagrams. a-c, Charge stability diagrams of
two electrons in the five dots L, B, R, T, C (a), or in the four dots L, B, R,
T (b, c). The charge state of the array of dots is controlled by sweeping
$\delta V_{\mathrm{X,Y}}^{-}$. The signal is recorded as linear combinations
of the quantum point contact current derivatives: (a, c)
$\partial_{V_{\mathrm{Y}}^{-}}i_{\mathrm{TL}}+i_{\mathrm{BL}}$, and (b)
$\partial_{V_{\mathrm{Y}}^{-}}i_{\mathrm{TL}}+i_{\mathrm{BL}}-i_{\mathrm{BR}}$.
The dashed white lines in (a) are guide for the eye for weak intensity
degeneracy lines. The stability diagrams are obtained for different gate
voltage tunings. d-f, Simulations of the stability diagrams (Methods). The
labels indicate the position of the isolated electron in the QD array.

## Coherent electron spin shuttling in a 2D array

We investigate the two-electron spin coherence while the electrons are
displaced within the (L, B, R, T, C) subset of quantum dots. Before all
coherent spin shuttling, two electrons in a singlet state (Methods) are loaded
in TL before being transferred to L then C (Fig. 1b). Next, a sequence of
voltage pulses $\delta V_{\mathrm{L,B,R,T}}$ is used to separate and displace
the electrons coherently in the array. By shuttling one or two electrons
within the dot array, the separated electrons will experience decoherence,
resulting in the spin mixing with the triplet states. Finally, the two
electrons are recombined in C and transferred to L and TL for spin readout
(Fig. 1c). To consider the three vectorial components of the hyperfine field
equally, the experiments are performed under zero external magnetic field [25,
15].

Figure 3:  Quasi-static electron spin coherence. a-d, Schematic pictures of
the different electron spin displacement sequences inside the 2D array of five
QDs employed to probe the electron spin coherence shown in (e). e, Singlet
probability plotted as a function of the total electron separation time
$\tau_{\mathrm{S}}$ spent in one (T-C, blue, a), two (T-C and T-L, orange, b),
three (T-L, T-C, and T-B, green, c), and four (T-C, T-L, T-B, and T-R, red, d)
different charge configurations. The two electrons are initialised in the
singlet state in the C-dot. Then, one electron is first transferred to T and
the second electron visits each other dot only once over an equal time of
$\tau_{R}$. As a consequence, the second electron spends for the different
sequences a total separation time $\tau_{S}=1\times\tau_{R}$ in C,
$\tau_{S}=2\times\tau_{R}$ in C and L, $\tau_{S}=3\times\tau_{R}$ in L, C, and
B, or $\tau_{S}=4\times\tau_{R}$ in C, L, B, and R, where $\tau_{R}$ is in
integer value of the AWG clock period. The data are fitted with Gaussian
decays (solid lines) with characteristic times $T_{2}^{\star}$ equal to
$7.9\pm 0.6$ (blue), $12.3\pm 0.7$ (orange), $15.8\pm 1.5$ (green), and
$21.9\pm 1.6$ $\mathrm{ns}$ (red). f, $\overline{T_{2}^{\star}}$ averaged over
different possible charge configurations plotted as a function
$N_{\textrm{config}}$ (black circles). The data are fitted with a square-root
function (red solid line)
$\overline{T_{2}^{\star}(N_{\mathrm{config}}=1)}\sqrt{N_{\mathrm{config}}}$,
with $\overline{T_{2}^{\star}(N_{\mathrm{config}}=1)}=9.6(0.7)$. The
experimental data of Fig. 4c,d are plotted as green squares
($T_{\textrm{coh}}$ of $32\pm 7$ and $44\pm 4$ $\mathrm{ns}$, respectively),
and their projection on the square root curve as open blue squares (equivalent
$N_{\textrm{config}}$ of $11$ and $21$, respectively).

We first investigate the quasi-static case where the electrons are separated
in different dots, one electron remaining in one dot, the other visiting a
single time an increasing number of dots. For a first electron fixed in dot T,
it results in $N_{\textrm{config}}=1$ to 4 distinct possible configurations
(Fig. 3a-d, respectively). The corresponding voltage pulse sequences are
depicted in Supp. Fig. S1. For each experiment, the time spent in each
configuration, referred to as the resting time $\tau_{R}$, is varied. The
resulting singlet probabilities are plotted (filled circles) in Fig. 3e as a
function of the total separation time
$\tau_{S}=N_{\textrm{config}}\times\tau_{R}$. The coherence time is directly
extracted by fitting a Gaussian decay
$\textrm{e}^{-(\tau_{S}/T_{\textrm{coh}})^{2}}$ (solid lines) and increases
with $\sqrt{N_{\mathrm{config}}}$. Indeed, the separated electrons explore a
large ensemble of nuclear spins during the displacement, averaging the
hyperfine interaction [26, 25, 27]. We have repeated this set of experiments
for different permutations of QDs to extract a statistical dependence of
$T_{\textrm{coh}}$ on $N_{\textrm{config}}$ and average out the influence of
the non-regular coupling between the different dots [24]. For each size of
visited subsets, the average and the standard deviation of $T_{\textrm{coh}}$
is plotted in Fig. 3f as a function of $N_{\textrm{config}}$. They
quantitatively fit with the expected square root law (solid line). It confirms
the importance of the hyperfine interaction and the number of nuclei the
electrons interact with. It also demonstrates the successful coherent
displacement within the 2D array of dots.

## Enhancement of the spin coherence time via electron shuttling

Figure 4:  Enhanced spin coherence time for periodic electron displacements.
a-b, Schematic pictures of the different electron spin periodic displacement
sequences inside the 2D array of five QDs employed to probe the electron spin
coherence shown in (c-f). c, d, Singlet probabilities plotted (black dots) as
a functions of the total electron separation time $\tau_{\mathrm{S}}$, c for
the single electron periodic displacement pattern (a), and d for the two-
electron periodic displacement pattern (b). The data were acquired for
configuration durations $\tau_{\mathrm{R}}$ of $1.7$ $\mathrm{ns}$ (c) and
$2.1$ $\mathrm{ns}$ (d). The experimental data are fitted using simulations of
the spin dynamics (see text) as solid blue lines. The fit error functions (see
text) are shown in insets. e, Singlet probability plotted (dots) as a function
of the total electron separation time $\tau_{\mathrm{S}}$, for increasing
(bottom to top) values of $\tau_{R}$: 1, 2, 3, 5, 7, 9, 11, and 13
$\mathrm{ns}$, for the two-electron periodic displacement case (b). The solid
lines are the simulated singlet probabilities (see text) computed for the
corresponding $\tau_{R}$ times, for a common remanence value $\mathcal{R}=0$.
f, Characteristic coherence times $T_{\textrm{coh}}$ plotted as functions of
$\tau_{R}$ and fitted with an inverse function
$T_{\textrm{coh}}^{0}/\tau_{R}+T_{2}^{\star}$, with
$T_{\textrm{coh}}^{0}=130\pm 24$ $\mathrm{ns}$.

We now analyse the impact of the electron tunnelling on their spin coherence.
In comparison with the previous case where the electrons explore a single time
distinct charge configurations, we study the situation where the charge
configurations are periodically explored. For each experiment, the resting
time spent in each charge configuration $\tau_{R}$ is set to a constant value.
We use two periodic pulse sequences corresponding to two-electron displacement
patterns: (i) single- (ii) two-electron periodic displacement. In the single-
electron periodic displacement case, the first electron is maintained in L,
and the second is displaced among B, C, T, and R, so that four isochronous dot
configurations are periodically visited (Fig. 4a and Supp. Fig. S2a). In the
two-electron periodic displacement case, an electron is first moved in B and
the second in R. From this point, each electron is sequentially displaced
along a two-step trajectory, either vertically, or horizontally. Eight
isochronous QD configurations are periodically visited with C the only dot
visited by both electrons (Fig. 4b and Supp. Fig. S2b). The experimental
singlet probabilities (black points) recorded for the single-electron
($\tau_{R}=1.7$ $\mathrm{ns}$), and for the two-electron ($\tau_{R}=2.1$
$\mathrm{ns}$), are plotted in Fig. 4c,d, respectively. The curves are fitted
by exponential decays (Supp. Figure S3a,b), with characteristic times of
$32\pm 7$ and $44\pm 4$ $\mathrm{ns}$, corresponding to coherence length of
$1.9\pm 0.4$ and $2.2\pm 0.2$ $\mathrm{\SIUnitSymbolMicro m}$ (assuming a
conservative interdot distance of 100 $\mathrm{nm}$), respectively. These
decay times are plotted as green squares in Fig. 3f, and their projection on
the square root curve as open blue squares ($N_{\textrm{config}}$ equivalent
to $11$ and $21$, respectively). A significant improvement is observed
compared to the quasi-static case: they are 1.6 times higher than the expected
coherence time for four and eight static dot configurations, respectively.
Moreover, the decoherence law changes from Gaussian-like (Fig. 3e) to
exponential (Fig. 4c,d) decay when visiting the charge configuration multiple
times, which confirms the significant impact of the electron dynamics on their
spin coherence.

By varying the resting time $\tau_{R}$, we can explore different shuttling
regimes for the two-electron periodic displacement case. The resulting singlet
probability is plotted in Fig. 4e as a function of
$\tau_{S}=(8N_{\textrm{cycle}}+1)\tau_{R}$ for different values of $\tau_{R}$
going from 1 $\mathrm{ns}$ to 13 $\mathrm{ns}$ (bottom to top). No singlet
probability decay is observed for $\tau_{R}$ longer than the static spin
coherence time due to complete mixing between the singlet and triplet states.
By progressively reducing $\tau_{R}$, we observe the emergence of a single
exponential decay of the singlet probability, demonstrating preservation of
the spin coherence during displacement. We extract the characteristic decay
time $T_{\textrm{coh}}$ by fitting with exponential decay (Supp. Fig. S3c)
which we plot against $\tau_{R}$ (Fig. 4f). The points are fitted with an
inverse function $\alpha/\tau_{R}$, typical of a motional narrowing
phenomenon, with a scale factor $\alpha=130\pm 24$ $\mathrm{n}\mathrm{s}^{2}$.
This value is in agreement with the expected square of the dephasing time
[23], estimated in Fig. 3f to be
$\overline{T_{2}^{\star}(N_{\mathrm{config}}=1)}^{2}=9.6^{2}=93(13)$
$\mathrm{n}\mathrm{s}^{2}$.

Figure 5:  Electron path randomisation consequences. a,b, Simulated
decoherence decay curves for one- (a) and two- (b) displaced electrons for
different values of the remanence parameter $\mathcal{R}$: 0 (blue), 0.5
(orange), 0.7 (green), 0.9 (red), and 1 (purple). c, Estimated remanence
$\mathcal{R}$ plotted as a function of the distance between the dot location
of two consecutive periods, for a charging energy of 1 (blue), 3 (orange), and
10 $\mathrm{meV}$ (green). d, Cross-sectional view of two-electron densities
computed using the fundamental Fock-Darwin state [28], and distant of 33
$\mathrm{nm}$ ($\mathcal{R}=0.5$), assuming a charging energy of 3
$\mathrm{meV}$.

## Path randomisation

If only a definite set of dot configurations are cyclically explored, no
difference in the spin dynamics should be observed between the quasi-static
and the periodic displacement protocols. Then, the only relevant parameter is
the time spent in each configuration, and a gaussian decay of the singlet
probabilities would be expected. However, the difference observed in the
experiment forces us to reconsider the initial assumptions made about the
electron displacement: the electrons are not periodically exploring the same
set of random magnetic field configurations. To better understand the impact
of the electron shuttling on spin coherence time, we model the system and
study the spin dynamics (Supp. Fig. S4 and Methods). We introduce the notion
of randomization of the effective hyperfine field during the coherent electron
displacements. It takes the form of a remanent parameter $0\leq\mathcal{R}\leq
1$, and it describes the change in effective hyperfine field, between two
consecutive displacement periods $i$ and $i+1$, at equivalent dot $j$. It is
expressed as
$\vec{B}_{i+1}^{j}=\mathcal{R}\vec{B}_{i}^{j}+(1-\mathcal{R})\vec{\delta
B}_{i+1}^{j}$, where $\vec{\delta B}_{i+1}^{j}$ is a random occurence of the
magnetic field distribution. For example, a full remanence $\mathcal{R}=1$
means that the effective hyperfine interaction for an electron visiting
multiple times the same dot is constant over the entire displacement sequence
of a single shot. Assuming (i) the evolution of the effective hyperfine field
much slower than the coherent spin displacement [29, 30], (ii) the nuclear
spin immune to the electrostatic manipulation, and (iii) not affected by the
electron motion, the remanence value is then interpreted as a small
electrostatic potential reorganization, leading to fluctuations in the dot
positions. For $\mathcal{R}<0.5$, the decoherence decay curves are
exponential-like (see Fig. 5a,b) and only weakly dependent on $\mathcal{R}$,
which is the fingerprint of a continuous motional narrowing effect. On the
other side, for $\mathcal{R}>0.5$, the motional narrowing process is limited
by the number of visited sites, which leads to Gaussian or Gaussian-like [25]
coherence decays. Such a transition from Gaussian to exponential dynamics is a
clear manifestation of the motional narrowing effect that has been observed in
liquid NMR.

The simulation well reproduces the experimental data in the single- and two-
electron displacement with a remanence parameter $\mathcal{R}\approx 0$ (solid
blue curve in Fig. 4c,d). The evaluation of the best fitting simulated
coherence decay curves (calculated to minimize the sum of the square
differences with the experimental data, insets) show that the data of the two
experiments can be reproduced with a value of remanence
$0\lesssim\mathcal{R}<0.5$. One can note that the difference in decay times
between the single- and the two-electron displacement cases is not significant
regarding the fit uncertainties (Supp. Figure S3a,b). Moreover, this parameter
also provides simulated coherence decays very similar to the experimental data
of Fig. 4e, for the various values $\tau_{R}$ of the displacement dynamics
(solid lines), as well as for the data obtained in the quasi-static
experiments (Supp. Fig. S5).

Therefore, the cyclic nature of the evolution has to be randomized along the
path of the electrons to reproduce the experimental data. Since the hyperfine
interaction is a contact interaction, a displacement of the size of the dot
will be sufficient to change the effective magnetic field completely. The
relationship between $\mathcal{R}$ and the electron centre of mass separation
between two consecutive periods $d$ is given for various charging energies in
Fig. 5c. It is calculated as the overlap between the electron density at two
consecutive periods $\mathcal{R}=\int\psi^{2}_{i}\psi^{2}_{i+1}dr$. Therefore,
between two consecutive displacement periods, the electron centre of mass is
expected to be distant of about $d\approx 33$ $\mathrm{nm}$ for a charging
energy of 3 $\mathrm{meV}$ (Fig. 5d and Supp. Fig. S6). The schematic electron
centre of mass trajectories in the single- and the two-electron displacement
experiments are superimposed on the gate pattern and the simulated
electrostatic potential of Supp. Fig. S6a and b, respectively.

Here we speculate that these fluctuations in dot positions may occur because
of the large voltage ranges employed to displace the electron through the
quantum dot array (as large as few hundreds of millivolts, Supp. Fig. S2).
Such energy is expected to induce perturbations of the semiconductor
nanostructure resulting in a change in the disorder potential. Besides, the
electrostatic potential simulations of Supp. Fig. S6a,b confirm that the
confinement potentials remain shallow during the electron shuttling, making
the dot positions very sensitive to fluctuations of charges in the substrate
(e.g. donors in the doping layer of the heterostructure, Supp. Fig. S6c).
Moreover, the typical timescale of this alteration is comparable to the one of
the excitation potential and therefore can be as fast as 1 nanosecond. It
implies that the electrons are displaced along a path that is fluctuating on a
timescale comparable to the time needed to explore few dot configurations. As
a consequence, the effective magnetic field configuration is not precisely
cyclic anymore and become randomized.

A similar assumption on the electron path has been reported recently by our
group in a different displacement regime [9]. The electrons were then
displaced in moving quantum dots at a speed of 3000
$\mathrm{m}\mathrm{/}\mathrm{s}$. Indeed, at this speed, the motional
narrowing process is extremely efficient for hyperfine interaction, and the
main mechanism for decoherence is mediated by spin-orbit interaction. Evidence
of change in the disorder potential was imprinted in the coherence of the
transported electrons initially prepared in a superposition of antiparallel
spin states.

## Conclusions

We have performed the one- and two-electron coherent shuttling in a two-
dimensional array and explored various displacement trajectories through
different sets of quantum dots. We report an increase of the coherence time,
which corresponds to a coherence length well above one micron for both one-
and two-electron displacements. Furthermore, we specifically study the zero
magnetic field regime where all the three components of the hyperfine
interaction contribute equally to the decoherence of the two electron-spin
singlet states. This hyperfine interaction is expected to be averaged during
the electron displacement. A signature of a motional narrowing process is
observed with a decoherence rate inversely proportional to the time spent in
the static phase. Indeed, reducing the time where the electrons are
effectively static and averaging faster than the spin dynamics over many
nuclear spin configurations increase the observed coherence time [22]. An
important consequence of the observation is the signature that the path of the
electron is slightly modified from one period of displacement to the next. It
results in the randomization of the effective field distribution experienced
by the electrons. Experimental results are reproduced with a simple model.
Therefore, only the single electron displacement during the gate movement from
one dot to another is relevant to explain the increase in spin coherence time.

## Methods

The methods about the spin initialisation in the singlet state, the
electrostatic potential, the stability diagram simulations and the virtual
gates employed in this study are essentially similar to [24].

Acquisition protocol. The voltage control apparatus is programmed to execute a
list of single-shot experiments (about a few thousand different sequences).
The total duration of each single-shot (comprising the electron spin
initialisation, manipulation, readout and the instrumental over-head) is about
50 $\mathrm{ms}$. The singlet spin probabilities are computed by averaging
1000 single-shot (repetition of the programmed list of sequences). Therefore,
the typical time between the repetition of two equivalent sequences is about a
few minutes.

Materials and set-up. Our device was fabricated using a Si-doped AlGaAs/GaAs
heterostructure grown by molecular beam epitaxy, with a two-dimensional
electron gas 100 $\mathrm{nm}$ below the crystal surface which has a carrier
mobility of $10^{6}$ $\mathrm{c}\mathrm{m}^{2}\mathrm{V}^{-1}\mathrm{s}^{-1}$
and an electron density of $2.7\times 10^{11}$ $\mathrm{c}\mathrm{m}^{2}$. The
quantum dots are defined using electrostatic confinement generated using Ti/Au
Schottky gates. It is anchored to the cold finger, which is in turn
mechanically attached to the mixing chamber of a homemade dilution
refrigerator with a base temperature of $60$ $\mathrm{m}\mathrm{K}$. It is
placed at the centre of a superconducting solenoid generating the static out-
of-plane magnetic field. Quantum dots are defined and controlled by applying
negative voltages on Ti/Au Schottky gates deposited on the surface of the
crystal. Homemade electronics ensure fast changes of both chemical potentials
and tunnel couplings with voltage pulse rise times approaching 100
$\mathrm{ns}$ and refreshed every 16 $\mathrm{\SIUnitSymbolMicro s}$.

A Tektronix 5014C with a typical channel voltage rise time (20% - 80%) of
$0.9$ $\mathrm{ns}$ are used to rapidly change the $V_{\mathrm{L}}$,
$V_{\mathrm{B}}$, $V_{\mathrm{R}}$, and $V_{\mathrm{T}}$ gate voltages. The
charge configurations can be read out by four quantum point contacts, tuned to
be sensitive local electrometers, and independently biased with 300
$\mathrm{\SIUnitSymbolMicro V}$. The resulting currents $i_{\text{TL}}$,
$i_{\text{BL}}$, $i_{\text{BR}}$ and $i_{\text{TR}}$ are measured using
current-to-voltage converters with a typical bandwidth of 10 $\mathrm{kHz}$.
The single-shot repetition rate is about 20 Hz.

Nuclear spin dynamics and experimental timescales The experimental singlet
probabilities are calculated as the average spin readout results of either 150
(Fig. 5a), or single 1000 shots (Figs. 3b and 4c,d). Sequences of multiple
displacement patterns are repeated to get a large enough number of shots. A
sequence has a typical duration of 1 day, with a sequence length of $\lesssim
4000$ different shots, and a shot duration of about 50 $\mathrm{ms}$.
Therefore, the time between two equivalent shots is about 200 $\mathrm{s}$,
which is much shorter than the nuclear spin relaxation time in the dipole-
dipole field of other nuclear spins is in the range of 100
$\mathrm{\SIUnitSymbolMicro s}$ [25, 31]. It allows the nuclear spin
distribution to be randomised between each equivalent single-shot experiment.
Consequently, there is no preferential quantification axis at zero magnetic
fields, and the spin-orbit interactions are averaged out.

Simulation of the electron spin dynamics. The simulations are performed with
an initial singlet spin state. A four-level system is considered to model the
2 electron spins, each submitted to a specific magnetic field (random
hyperfine interaction with the substrate nuclear spins). The numerical time-
integration of the Schrödinger equation is computed for 1000 different couples
of magnetic field vectors, following a centered Gaussian distribution with a
standard deviation of 3.9 $\mathrm{mT}$ (see Fig. S4). Magnetic field
remanence parameters are added to the model to consider the spatial evolution
of the effective dot locations during the different displacement patterns.
Finally, the average singlet probability is calculated [32, 33].

## References

  * Burkard and Loss [2002] G. Burkard and D. Loss, Electron spins in quantum dots as qubits for quantum information processing, in  _Semiconductor Spintronics and Quantum Computation_ (Springer-Verlag Berlin Heidelberg, 2002).
  * Taylor  _et al._ [2005] J. M. Taylor, H.-A. Engel, W. Dur, A. Yacoby, C. M. Marcus, P. Zoller, and M. D. Lukin, Fault-tolerant architecture for quantum computation using electrically controlled semiconductor spins, Nature Phys. 1, 177 (2005).
  * Greentree  _et al._ [2004] A. D. Greentree, J. H. Cole, A. R. Hamilton, and L. C. L. Hollenberg, Coherent electronic transfer in quantum dot systems using adiabatic passage, Phys. Rev. B 70, 235317 (2004).
  * Skinner  _et al._ [2003] A. J. Skinner, M. E. Davenport, and B. E. Kane, Hydrogenic spin quantum computing in silicon: A digital approach, Phys. Rev. Lett. 90, 087901 (2003).
  * Barnes  _et al._ [2000] C. H. W. Barnes, J. M. Shilton, and A. M. Robinson, Quantum computation using electrons trapped by surface acoustic waves, Phys. Rev. B 62, 8410 (2000).
  * Vandersypen  _et al._ [2017] L. M. K. Vandersypen, H. Bluhm, J. S. Clarke, A. S. Dzurak, J. Ishihara, A. Morello, D. J. Reilly, L. R. Schreider, and M. Veldhorst, Interfacing spin qubits in quantum dots and donors-hot, dense, and coherent, npj Quant. Inf. 3, 34 (2017).
  * Mills  _et al._ [2019] A. R. Mills, D. M. Zajac, M. J. Gullans, F. J. Schupp, T. M. Hazard, and J. R. Petta, Shuttling a single charge across a one-dimensional array of silicon quantum dots, Nature Communications 10, 1063 (2019).
  * Bertrand  _et al._ [2016] B. Bertrand, S. Hermelin, S. Takada, M. Yamamoto, S. Tarucha, A. Ludwig, A. D. Wieck, C. Bäuerle, and T. Meunier, Fast spin information transfer between distant quantum dots using individual electrons, Nature Nanotech. 11, 672 (2016).
  * Jadot  _et al._ [2020] B. Jadot, P. A. Mortemousque, E. Chanrion, V. Thiney, A. Ludwig, A. D. Wieck, M. Urdampilleta, C. Bauerle, and T. Meunier, arXiv:2004.02727 (2020).
  * Vinet  _et al._ [2018] M. Vinet, L. Hutin, B. Bertrand, S. Barraud, J. Hartmann, Y. Kim, V. Mazzocchi, A. Amisse, H. Bohuslavskyi, L. Bourdet, A. Crippa, X. Jehl, R. Maurand, Y. M. Niquet, M. Sanquer, B. Venitucci, B. Jadot, E. Chanrion, P.-A. Mortemousque, C. Spence, M. Urdampilleta, S. D. Franceschi, and T. Meunier, Towards scalable silicon quantum computing, in  _2018 IEEE International Electron Devices Meeting (IEDM)_ (2018) pp. 6.5.1–6.5.4.
  * Li  _et al._ [2018] R. Li, L. Petit, D. P. Franke, J. P. Dehollain, J. Helsen, M. Steudtner, N. K. Thomas, Z. R. Yoscovits, K. J. Singh, S. Wehner, _et al._ , A crossbar network for silicon quantum dot qubits, Sci. Adv. 4, eaar3960 (2018).
  * Fujita  _et al._ [2017] T. Fujita, T. A. Baart, C. Reichl, W. Wegscheider, and L. M. K. Vandersypen, Coherent shuttle of electron-spin states, npj Quant. Inf. 3, 22 (2017).
  * [13] J. Yoneda, W. Huang, M. Feng, C. H. Yang, K. W. Chan, T. Tanttu, W. Gilbert, R. C. C. Leon, F. E. Hudson, K. M. Itoh, A. Morello, S. D. Barlett, A. Laucht, A. Saraiva, and A. S. Dzurak, Coherent spin qubit transport in silicon, arXiv:2008.04020 .
  * Ginzel  _et al._ [2020] F. Ginzel, A. R. Mills, J. R. Petta, and G. Burkard, Spin shuttling in a silicon double quantum dot, Phys. Rev. B 102, 195418 (2020).
  * Flentje  _et al._ [2017a] H. Flentje, P. A. Mortemousque, R. Thalineau, A. Ludwig, A. D. Wieck, C. Bäuerle, and T. Meunier, Coherent long-distance displacement of individual electron spins, Nat. Commun. 8, 501 (2017a).
  * Yoneda  _et al._ [2018] J. Yoneda, K. Takeda, T. Otsuka, T. Nakajima, M. R. Delbecq, G. Allison, T. Honda, T. Kodera, S. Oda, Y. Hoshi, N. Usami, K. M. Itoh, and S. Tarucha, A quantum-dot spin qubit with coherence limited by charge noise and fidelity higher than $99.9\%$, Nature Nanotechnology 13, 102 (2018).
  * Veldhorst  _et al._ [2014] M. Veldhorst, J. C. C. Hwang, C. H. Yang, A. W. Leenstra, B. de Ronde, J. P. Dehollain, J. T. Muhonen, F. E. Hudson, K. M. Itoh, A. Morello, and A. S. Dzurak, An addressable quantum dot qubit with fault-tolerant control-fidelity, Nature Nanotechnology 9, 981 (2014).
  * Tyryshkin  _et al._ [2011] A. M. Tyryshkin, S. Tojo, J. J. L. Morton, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, T. Schenkel, M. L. W. Thewalt, K. M. Itoh, and S. A. Lyon, Electron spin coherence exceeding seconds in high-purity silicon., Nature materials 11, 143 (2011).
  * Shulman  _et al._ [2014] M. D. Shulman, S. P. Harvey, J. M. Nichol, S. D. Bartlett, A. C. Doherty, V. Umansky, and A. Yacoby, Suppressing qubit dephasing using real-time hamiltonian estimation, Nature Communications 5, 5156 (2014).
  * Wolfowicz  _et al._ [2016] G. Wolfowicz, P.-A. Mortemousque, R. Guichard, S. Simmons, M. L. W. Thewalt, K. M. Itoh, and J. J. L. Morton, 29si nuclear spins as a resource for donor spin qubits in silicon, New Journal of Physics 18, 023021 (2016).
  * Bluhm  _et al._ [2011] H. Bluhm, S. Foletti, I. Neder, M. Rudner, D. Mahalu, V. Umansky, and A. Yacoby, Dephasing time of gaas electron-spin qubitscoupled to a nuclear bath exceeding 200 microsecond, Nature Physics 7, 109 (2011).
  * Huang and Hu [2013] P. Huang and X. Hu, Spin qubit relaxation in a moving quantum dot, Phys. Rev. B 88, 075301 (2013).
  * Slichter [1978] C. P. Slichter,  _Principles of Magnetic Resonance_ , 3rd ed., edited by Springer-Verlag (Springer-Verlag Berlin Heidelberg, 1978).
  * Mortemousque  _et al._ [2020] P.-A. Mortemousque, E. Chanrion, B. Jadot, H. Flentje, A. Ludwig, A. D. Wieck, M. Urdampilleta, C. Bäuerle, and T. Meunier, Coherent control of individual electron spins in a two-dimensional quantum dot array, Nature Nanotechnology 10.1038/s41565-020-00816-w (2020).
  * Merkulov  _et al._ [2002] I. A. Merkulov, A. L. Efros, and M. Rosen, Electron spin relaxation by nuclei in semiconductor quantum dots, Physical Review B 65, 205309 (2002).
  * Flentje  _et al._ [2017b] H. Flentje, B. Bertrand, P. A. Mortemousque, V. Thiney, A. Ludwig, A. D. Wieck, C. Bäuerle, and T. Meunier, A linear triple quantum dot system in isolated configuration, Appl. Phys. Lett. 110, 233101 (2017b).
  * Bloembergen  _et al._ [1948] N. Bloembergen, E. M. Purcell, and R. V. Pound, Relaxation effects in nuclear magnetic resonance absorption, Phys. Rev. 73, 679 (1948).
  * Kouwenhoven  _et al._ [2001] L. P. Kouwenhoven, D. G. Austing, and S. Tarucha, Rep. Prog. Phys. 64, 701 (2001).
  * Reilly  _et al._ [2008] D. J. Reilly, J. M. Taylor, E. A. Laird, J. R. Petta, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Measurement of temporal correlations of the overhauser field in a double quantum dot, Phys. Rev. Lett. 101, 236803 (2008).
  * Nakajima  _et al._ [2020] T. Nakajima, A. Noiri, K. Kawasaki, J. Yoneda, P. Stano, S. Amaha, T. Otsuka, K. Takeda, M. R. Delbecq, G. Allison, A. Ludwig, A. D. Wieck, D. Loss, and S. Tarucha, Coherence of a driven electron spin qubit actively decoupled from quasistatic noise, Phys. Rev. X 10, 011060 (2020).
  * Malinowski  _et al._ [2017] F. K. Malinowski, F. Martins, L. Cywiński, M. S. Rudner, P. D. Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, Spectrum of the nuclear environment for gaas spin qubits, Phys. Rev. Lett. 118, 177702 (2017).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Comp. Phys. Comm. 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, Comp. Phys. Comm. 184, 1234 (2013).
  * Bautze  _et al._ [2014] T. Bautze, C. Süssmeier, S. Takada, C. Groth, T. Meunier, M. Yamamoto, S. Tarucha, X. Waintal, and C. Bäuerle, Theoretical, numerical, and experimental study of a flying qubit electronic interferometer, Phys. Rev. B 89, 125432 (2014).

## Acknowledgements

We would like to thank Maud Vinet, Xuedong Hu, and Lieven M. K. Vandersypen
for enlightening discussions. We acknowledge technical support from the Poles
of the Institut Néel, and in particular, the NANOFAB team who helped with the
sample realisation, as well as E. Eyraud, T. Crozes, P. Perrier, G. Pont, H.
Rodenas, D. Lepoittevin, C. Hoarau and C. Guttin. M.U. acknowledges the
support of project CODAQ (ANR-16-ACHN-0029). A.L. and A.D.W. acknowledge
gratefully the support of DFG-TRR160, BMBF-Q.com-H 16KIS0109, and the DFH/UFA
CDFA-05-06. T.M. acknowledges ﬁnancial support from ERC QSPINMOTION, ERC
QUCUBE, ANR CMOSQSPIN (Grant No. ANR-17-CE24-0009), ANR SiQuBus and UGA IDEX
(Grant No. ANR-15-IDEX-02).

## Author contributions statement

P.-A. M. fabricated the sample and performed the experiments with the help of
T.M. and C.B.. P.-A. M. and T.M. interpreted the data. P.-A. M. and T.M. wrote
the manuscript with the input of all the other authors. H.F. contributed to
the experimental setup. A.L. and A.D.W. performed the design and molecular-
beam-epitaxy growth of the high mobility heterostructure. All authors
discussed the results extensively, as well as the manuscript.

Additional information Correspondence and requests for materials should be
addressed to<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
Supplemental Materials:

Figure S1:  Single-period electron displacement sequences. a-d, Sequences of
nanosecond voltage pulses used to perform the single period electron
displacement of Fig. 3. a, The pulses are applied during $\tau_{S}$ to
displace the electrons in the TC configuration. b, The pulses are applied
during $\tau_{S}/2$ to displace the electrons in the TC and LT configurations.
c, The pulses are applied during $\tau_{S}/3$ to displace the electrons in the
LT, TC and BT configurations. d, The pulses are applied during $\tau_{S}/4$ to
displace the electrons in the TC, LT, BT and RC configurations. The x-axes
labels denote the location of the dot occupied by the electron. For example,
TC stands for one electron in dot T, and one in dot C. The initial and final
configurations correspond to the two electrons in the C dot. The schematic
representations of the charge configurations are superimposed. Only one period
of displacement is shown for each sequence. The voltage rising times are given
as illustrations.

Figure S2:  Periodic electron displacement sequences. a,b, Sequences of
nanosecond voltage pulses used to perform the periodic electron displacements.
a, One electron remains in the left quantum dot while the second electron is
displaced in the B, C, T, R dots. b, One electron periodically visits L, C, R,
and the second T, C, B. The x-axes labels denote the location of the dot
occupied by the electron. For example, LB stands for one electron in dot L,
and one in dot B. The schematic representations of the charge configurations
are superimposed. Only one period of displacement is shown for each sequence.
The voltage rising times are given as illustrations.

Figure S3:  Exponential spin coherence fits of the periodically displaced
electron. In order to extract the coherence time $\textrm{T}_{\mathrm{coh}}$
that are discussed in the manuscript, we fit the experimental data using a
single exponential decay law. a, b Same dataset as in Fig. 4c and d,
respectively. The data are fitted with an exponential decay (solid lines) with
characteristic times of (a) $32\pm 7$, and (b) $44\pm 4$ $\mathrm{ns}$ (red).
c Same dataset as in Fig. 4e. The data are fitted with an exponential decay
(solid lines) with the characteristic times: $91.3\pm 10$ (blue), $44.8\pm 7$
(orange), $51.7\pm 8$ (green), $31.3\pm 9$ (red), $31.2\pm 15$ (purple), and
$8.9\pm 9$ (brown).

Figure S4:  Simulated spin coherence time in a static double quantum dot. The
electron spin decoherence curves can be simulated for different standard
deviation $\sigma(B_{\textrm{Z}})$ values of a Gaussian distribution in
hyperfine fields (Methods). The coherence times extracted from Gaussian fits
are plotted as functions of $\sigma(B_{\textrm{Z}})$. The field corresponding
to the spin coherence time $T_{2}^{\star}(N_{\textrm{config}}=1)=9.6$
$\mathrm{ns}$ (extracted from the $T_{2}^{\star}(\sqrt{N_{\textrm{config}}})$
fit of Fig. 3f) is pinpointed in dashed red line ($B_{\textrm{Z}}=3.9$
$\mathrm{mT}$).

Figure S5:  Simulated decays of the quasi-static electron spin coherence. Same
dataset as in Fig. 3e. The experimental data (dots) are compared with the
simulations (solid lines) of the coherence spin in the quasi-static case
(remanence parameter is irrelevant here as each configuration is explored only
a single time). The overshoot in singlet probability described in [25] is
observed in the simulations only for the double quantum dot case (blue), and
is not observed in the experimental data. For the two electrons in three
(orange), four (green), or five (red) dots, the simulated decays are Gaussian-
like.

Figure S6:  Electrostatic potential simulations. a, b, Electrostatic potential
simulations performed with a gate voltage configurations corresponding to a
confinement of two electrons (a) in L-B, and (b) in T-R. Hypothetical electron
trajectories (dots) are overlapped for the two displacement patterns
considered. The solid lines are guides for the eye. a, A single electron is
periodically displaced between dots (B, C, T, and R), the other electron
remains in the L dot. Four dot configurations are explored, but the dot
positions evolve at each period (two periods are drawn, e.g. $1\rightarrow
1^{\prime}\rightarrow 1"$). b, Two electrons are displaced, one between (L, C,
R) and the other between (T, C, B). Eight dot configurations are explored and
here again, the dot positions evolve at each displacement period (a single
period is drawn, e.g. $1\rightarrow 1^{\prime}$). The electrostatic potential
simulations were performed by solving the Laplace equation [34]. The orange
circles have a radius of 20 $\mathrm{nm}$ corresponding to the Bohr radius of
an electron with a charging energy $E_{\mathrm{C}}$ of 3 $\mathrm{meV}$. The
dashed circles have a radius of 33 $\mathrm{nm}$, value computed to reach a
remanence $\mathcal{R}=0.5$. c, Electrostatic potential simulation of a random
donor distribution in the Si-doped layer of the heterostructure (Methods). The
induced potential can localise the dots in different spatial positions than
the one intended by the gate geometry. d-f, Schematic representations of the
Bohr radii and corresponding dot displacements to reach a remanence
$\mathcal{R}=0.5$ for different electron charging energies.


# Gravitational waves from rapid structure formation on microscopic scales
before matter-radiation equality

Marcos M. Flores  Department of Physics and Astronomy, University of
California, Los Angeles
Los Angeles, California, 90095-1547, USA  Alexander Kusenko  Department of
Physics and Astronomy, University of California, Los Angeles
Los Angeles, California, 90095-1547, USA  Kavli Institute for the Physics and
Mathematics of the Universe (WPI), UTIAS
The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Theoretical Physics
Department, CERN, 1211 Geneva 23, Switzerland  Misao Sasaki  Kavli Institute
for the Physics and Mathematics of the Universe (WPI), UTIAS
The University of Tokyo, Kashiwa, Chiba 277-8583, Japan  Center for
Gravitational Physics and Quantum Information, Yukawa Institute for
Theoretical Physics,
Kyoto University, Kyoto 606-8502, Japan  Leung Center for Cosmology and
Particle Astrophysics, National Taiwan University,
Taipei 10617, Taiwan

###### Abstract

The existence of scalar fields can be probed by observations of stochastic
gravitational waves. Scalar fields mediate attractive forces, usually stronger
than gravity, on the length scales shorter than their Compton wavelengths,
which can be non-negligible in the early universe, when the horizon size is
small. These attractive forces exhibit an instability similar to the
gravitational instability, only stronger. They can, therefore, lead to the
growth of structures in some species. We identify a gravitational waves
signature of such processes and show that it can be detected by the future
gravitational waves experiments.

††preprint: IPMU22-0019††preprint: YITP-22-95

Most of the elementary particles discovered in the past century have Compton
wavelengths that are much smaller than the size of the atom, and, therefore,
they cannot mediate any long-range forces between atoms. However, in the early
Universe with the total energy density $\rho$, any force mediated by particles
with masses smaller than $\rho^{1/2}/M_{\rm Pl}$, where $M_{\rm Pl}$ is the
reduced Planck mass, acts as a long-range force on the scales of the Hubble
horizon at that time. Of particular interest are the forces that mediate
attractive interactions between particles, such as the Yukawa forces. The
Yukawa interactions are always attractive, and, unless the coupling is
unusually small, such forces are much stronger than gravity. This has profound
effect on the dynamics of any particles participating in such interactions.
The attractive nature of the force implies an instability similar to
gravitational instability, and the growth of structures in the distribution of
density of some particle species is possible even in the radiation dominated
universe Amendola  _et al._ (2018); Savastano  _et al._ (2019); Casas  _et
al._ (2016); Flores and Kusenko (2021a); Domènech and Sasaki (2021); Flores
and Kusenko (2021b). This growth of isocurvature overdensities can lead to
bound states Wise and Zhang (2014); Gresham  _et al._ (2018); Gresham and
Zurek (2019) or, accompanied by the radiative cooling caused by the same
Yukawa interaction, can lead to formation of primordial black holes Flores and
Kusenko (2021a, b).

We will explore a generic observable manifestation of the structure formation
due to the presence of the attractive forces in the early Universe, namely,
stochastic gravitational waves (GWs). We will show that such an event in the
early Universe during either a radiation dominated (RD) era or an intermediate
matter dominated (IMD) era can produce a signal detectable by the future GW
detectors. While previous work has demonstrated that gravitational forces
alone can generate GWs Dalianis and Kouvaris (2021); Jedamzik  _et al._
(2010); Domènech  _et al._ (2022), here we will explore the effects of a force
significantly stronger than gravity.

The formation of dark matter overdensities due to a long range scalar force
has been examined in many different contexts Amendola  _et al._ (2018); Flores
and Kusenko (2021a); Domènech and Sasaki (2021); Savastano  _et al._ (2019).
As in Ref. Flores and Kusenko (2021a), we consider a heavy fermion $\psi$
interacting with a scalar field $\chi$:

$\mathcal{L}\supset\frac{1}{2}m_{\chi}^{2}\chi^{2}-y\chi\bar{\psi}\psi+\cdots.$
(1)

We will assume that the scalar is either massless or very light, $m_{\chi}\ll
m_{\psi}^{2}/M_{\rm Pl}$, where $M_{\rm Pl}\approx 2\times 10^{18}$ GeV.
Further, we assume that the fermions $\psi$ are either stable or have a total
decay width $\Gamma\ll m_{\psi}^{2}/M_{\rm Pl}$ which ensures there is a
cosmological epoch where the $\psi$ particles can become nonrelativistic,
decoupled from equilibrium and interact via long-range force mediated by the
$\chi$ field.

During a RD era, gravity alone only allows for logarithmic growth of density
perturbations whereas in an IMD era, gravitational forces cause
$\delta(x,t)=\delta\rho/\bar{\rho}$ to grow linearly with scale factor
Mukhanov (2005) i.e. $\delta\propto a$. In either case, the presence of a
long-range “fifth force” stronger than gravity causes the fluctuations to grow
at a much faster rate. We note that the scalar force does not couple to the
mass density, but to the $\psi$ number density. The strength of this force can
be characterized by its relative strength with gravity, $\beta\equiv y(M_{\rm
Pl}/m_{\psi})\gg 1$. This causes fluctuations $\Delta(x,t)=\Delta
n_{\psi}/\bar{n}_{\psi}$ to grow more rapidly so long as the $\psi$ is
decoupled from radiation so that pressure can be neglected.

In Fourier space, the growth of these perturbations is described by the system
of coupled equations Gradwohl and Frieman (1992); Gubser and Peebles (2004);
Nusser  _et al._ (2005); Amendola  _et al._ (2018); Domènech and Sasaki
(2021),

$\displaystyle\ddot{\delta}_{k}+2H\dot{\delta}_{k}-\frac{3}{2}H^{2}(\Omega_{r}\delta_{k}+\Omega_{m}\Delta_{k})$
$\displaystyle=0$ (2)
$\displaystyle\ddot{\Delta}_{k}+2H\dot{\Delta}_{k}-\frac{3}{2}H^{2}\left[\Omega_{r}\delta_{k}+\Omega_{m}(1+\beta^{2})\Delta_{k}\right]$
$\displaystyle=0.$ (3)

where $\Omega_{r}=\rho_{r}/(\rho_{r}+\rho_{\psi})$ and
$\Omega_{m}=\rho_{\psi}/(\rho_{r}+\rho_{\psi})$ are the radiation and matter
components respectively such that $\Omega_{r}+\Omega_{m}=1$ and $H$ is the
Hubble parameter $H\equiv\dot{a}/a$. Assuming that radiation and the $\psi$
particles are the only relevant energy components, then the time dependence of
the energy fractions are

$\Omega_{r}=\frac{1}{1+(t/t_{\rm eq})^{s}},\qquad\Omega_{m}=\frac{1}{1+(t_{\rm
eq}/t)^{s}}$ (4)

where $t_{\rm eq}$ is the matter-radiation equality time when the $\psi$
component comes to dominate and $s$ is defined as the power $a\propto t^{s}$
such that $s=1/2$ for radiation domination and $s=2/3$ for matter domination.
We can examine the solution of (3) in the limit that $\beta\gg 1$ in either
era.

$\Delta_{k}\propto a^{p},\quad p=\begin{cases}\sqrt{6(1+\beta^{2})}&{\rm
RD}\\\\[7.11317pt] \frac{1}{4}\left(\sqrt{1+24(1+\beta^{2})}-1\right)&{\rm
MD}\end{cases}$ (5)

where $p$ can be much larger than one. This growth continues until the density
contrast $\Delta_{k}$ enters the nonlinear regime when $\Delta_{k}\simeq 1$.

Our discussion of the evolution of the overdensities and their subsequent
collapse will closely follow the work of Dalianis and Kouvaris Dalianis and
Kouvaris (2021). Rather than utilizing the formalism of induced GWs from a
decaying gravitational potential with time, Dalianis and Kouvaris (2021)
applies an alternative analytic approach based on the Zel’dovich approximation
for the nonlinear evolution of density perturbations Zel’dovich (1970). As
described by Zel’dovich, deviations in spherical symmetry will cause
instability as time goes on, leading to flattening in one particular
direction. These objects, often referred to as Zel’dovich pancakes, may lead
to the formation of black holes if the final configuration satisfies hoop
conjectures Thorne (1972); Misner  _et al._ (1973). If not, the system will
continue to collapse and go through a phase crossings and oscillations which
lead to a virialized halo. In the absence of an asymmetry in the $\psi$
component, this leads to the destruction of the halo itself through
$\bar{\psi}\psi$ annihilation. This process halts the further emission of GWs
and, in the case of an IMD era, reestablish the radiation era. Our approach,
similar to Ref. Dalianis and Kouvaris (2021), will only focus on the nonlinear
collapse process within the Zel’dovich approximation and the gravitational
signal which may exist due to the relaxation process will be left for future
work.

Within the Zel’dovich formalism the coordinate of a particle is written as

$r_{i}=a(t)q_{i}+b(t)p_{i}(q_{i})$ (6)

where $a(t)$ is the usual scale factor encoding the expansion of the Universe,
$q_{i}$ is the comoving coordinate, $b(t)$ is the growing mode describing the
instability of the overdensity and $p_{i}$ are deviation vectors that depend
on the initial perturbation.

To study the motion of a group of particles around $q_{i}$, a deformation or
strain tensor $D_{ik}$ is defined as

$\begin{split}D_{ik}&=\frac{\partial r_{i}}{\partial
q_{k}}=a(t)\delta_{ik}+b(t)\frac{\partial p_{i}}{\partial
q_{k}}\\\\[7.11317pt] &={\rm diag}(a-\alpha b,a-\beta b,a-\gamma
b)\end{split}$ (7)

where we have chosen a basis so that $\partial p_{i}/\partial q_{k}$ is
diagonal,

$\frac{\partial p_{i}}{\partial q_{k}}=-{\rm diag}(\alpha,\beta,\gamma).$ (8)

The mass contained within the Lagrangian volume is,

$M_{\psi}=\int\rho_{\psi}\ d^{3}r=\bar{\rho}_{\psi}a^{3}\int d^{3}q$ (9)

where the difference between the two right-most expressions is the Jacobian
determinant, i.e.

$\rho_{\psi}(a-\alpha b)(a-\beta b)(a-\gamma b)=\bar{\rho}_{\psi}a^{3}.$ (10)

We will also assume that the initial deviations are small, so that the
perturbations are nearly spherical. Specifically, when a perturbation of size
$q$ and wavenumber $k=q^{-1}$ enters the horizon at $t_{q}$,

$a(t_{q})q=H^{-1}(t_{q})$ (11)

we assume that $\alpha b(t_{q})/a(t_{q})\ll 1$, $\beta b(t_{q})/a(t_{q})\ll 1$
and $\gamma b(t_{q})/a(t_{q})\ll 1$ so that we can ensure that the entire
ellipsoid is within the Hubble volume. To leading order, the density contrast
is given by

$\delta_{L}\equiv\left(\frac{\rho_{\psi}-\bar{\rho}_{\psi}}{\bar{\rho}_{\psi}}\right)_{L}=(\alpha+\beta+\gamma)\frac{b}{a}.$
(12)

Note that for $b>0$, we have $\delta_{L}>0$ if and only if
$\alpha+\beta+\gamma>0$. As discussed previously, perturbations at most grow
linearly with scale factor. The presence of a long-range scalar force enhances
this growth and so, the key assumption of this work is that $\delta_{L}\propto
a^{p}$, $p\geq 1$. Crucially, the $p=1$, $s=2/3$ should reproduce precisely
the work of Ref. Dalianis and Kouvaris (2021), which we shall demonstrate
later.

Under this assumption, (12) implies that $b\propto a^{p+1}$ and therefore
grows faster than the scale factor. The perturbation will collapse along at
least one of the three axes. Without loss of generality we assume that

$\alpha>0,\quad-\infty<\gamma\leq\beta\leq\alpha<\infty,\quad\alpha+\beta+\gamma>0.$
(13)

Additionally, we can determine the timescale $t_{q}$ utilizing (9) and (11) to
give

$t_{q}=2sGM_{\psi}\left(1+\frac{\Omega_{\rm
rad}(t_{q})}{\Omega_{\psi}(t_{q})}\right)$ (14)

where $\Omega_{X}=\rho_{X}/\rho_{\rm crit}$. The horizon entry timescale is
only one of three important events in our scenario. After entering the
horizon, the perturbation grows until reaching some maximum size at
$t_{\max}$. The overdensity subsequently collapses, at $t_{\rm col}$, due to
the dominant scalar force. In this work, we will only examine the effects of
scalar-induced collapse, though scalar cooling might dramatically decrease the
time needed for the halo to collapse Flores and Kusenko (2021a, c).

The derivation of $t_{\max}$ and $t_{\rm col}$ follows a similar methodology
as in Ref. Dalianis and Kouvaris (2021) the details of which have been
included in the supplementary material. We find,

$t_{\max}=\left(\frac{\alpha+\beta+\gamma}{(p+1)\alpha\delta_{L}(t_{q})}\right)^{1/ps}t_{q}.$
(15)

and

$t_{\rm col}=(p+1)^{1/ps}t_{\max}.$ (16)

The deviations $\alpha$, $\beta$, $\gamma$ are determined by the Doroshkevich
probability function. Specified as Doroshkevich (1970),

$\begin{split}&\mathcal{F}_{D}(\alpha,\beta,\gamma,\sigma_{3})=-\frac{27}{8\sqrt{5}\pi\sigma_{3}^{6}}(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)\\\
&\times\exp\left[-\frac{3}{5\sigma_{3}^{2}}\left((\alpha^{2}+\beta^{2}+\gamma^{2})-\frac{1}{2}(\alpha\beta+\beta\gamma+\gamma\alpha)\right)\right]\end{split}$
(17)

for $\alpha\geq\beta\geq\gamma$. The probability function is normalized to
one,

$\int_{-\infty}^{\infty}d\alpha\int_{-\infty}^{\alpha}d\beta\int_{-\infty}^{\beta}d\gamma\
\mathcal{F}_{D}(\alpha,\beta,\gamma,\sigma_{3})=1.$ (18)

We will denote the standard deviation of $\delta_{L}(t_{q})$ as $\sigma$,
given explicitly as

$\sigma^{2}=\langle\delta_{L}^{2}(t_{q})\rangle=\langle(\alpha+\beta+\gamma)^{2}\rangle\left.\left(\frac{b}{a}\right)\right|_{t=t_{q}}=5\sigma_{3}^{2}.$
(19)

For this expression, we fixed the normalization such that
$b(t)/a(t)=(a(t)/a(t_{q}))^{p}$. In addition to our restrictions (13) we will
also require that the maximum expansion time $t_{\rm max}$ is larger than the
time of entry $t_{q}$. This defines a subspace $\mathcal{S}$,

$\begin{split}0<\alpha&<\infty\\\
-\frac{\alpha}{2}[1-(p+1)\sigma]<\beta&<\alpha\\\
-\beta-\alpha[1-(p+1)\sigma]<\gamma&<\beta\end{split}$ (20)

over which we will limit our calculations.

The stochastic GW background is traditionally characterized by the normalized
energy density,

$\Omega_{\rm GW}(t,f)\equiv\frac{1}{\rho_{\rm crit}}\frac{d\rho_{\rm GW}}{d\ln
f}.$ (21)

Following a framework similar to Ref. Dalianis and Kouvaris (2021), the
present day GW energy density fraction reads

$\begin{split}\Omega_{\rm GW}(t_{0},f_{0})&=\frac{1}{\rho_{\rm
crit}(t_{0})}\frac{2\pi
G}{5c^{5}}\int_{\mathcal{S}}d^{3}{\bm{\alpha}}\sum_{n=1}^{N}\frac{1}{1+z_{n}}\\\\[7.11317pt]
&\times\sum_{ij}|\tilde{\dddot{Q}}_{ij}^{n}(2\pi
f_{0}(1+z_{n}))|^{2}\\\\[7.11317pt] &\times(2\pi
f_{0}(1+z_{n}))\left(\frac{4\pi}{3}q^{3}\right)^{-1}\\\\[7.11317pt]
&\times\mathcal{F}_{D}({\bm{\alpha}},\sigma)\Theta(t_{\rm CO}-t_{\rm
col}({\bm{\alpha}})).\end{split}$ (22)

where the sum is over partitioned time intervals from the initial time to the
final time of emission and $\tilde{\dddot{Q}}_{ij}$ is the Fourier transform
of the third time derivative of the quadrupole moment. The Heaviside function
enforces the condition that collapse occurs before some “cut-off time.” In the
context of a RD era, this time could correspond to lifetime of the constituent
$\psi$ particles where in the IMD era $t_{\rm CO}$ is simply the time when
radiation domination is reestablished.

As in Ref. Dalianis and Kouvaris (2021), we will assume that all of the GWs
are emitted in a single interval i.e. $N=1$. Partitioning the interval was
motivated by the fact that GWs emitted in different subintervals will be
redshifted different amounts as they propagate. For simplicity, we assume that
the GWs are emitted instantaneously at $t_{\rm col}$. In view of the fact that
$t_{\rm col}=(p+1)^{1/ps}t_{\max}$, with $p\geq 1$, we see that this
approximation not only is reasonable, but becomes better as $p$ increases. In
this approximation the final expression for $\Omega_{\rm GW}(t_{0},f_{0})$ is,

$\begin{split}&\Omega_{\rm GW}(t_{0},f_{0})=\frac{2\pi f_{0}}{\rho_{\rm
crit}(t_{0})}\left(\frac{4\pi}{3}q^{3}\right)^{-1}\frac{2\pi
G}{5c^{5}}\\\\[7.11317pt]
&\times\int_{\mathcal{S}}d^{3}{\bm{\alpha}}\sum_{ij}|\tilde{\dddot{Q}}_{ij}(2\pi
f_{0}(1+z_{\rm col}({\bm{\alpha}})))|^{2}\\\\[7.11317pt] &\times(2\pi
f_{0}(1+z_{\rm
col}({\bm{\alpha}})))\mathcal{F}_{D}({\bm{\alpha}},\sigma)\Theta(t_{\rm
CO}-t_{\rm col}({\bm{\alpha}})).\end{split}$ (23)

Here, the redshift $z_{\rm col}$ is related to the collision timescale through

Figure 1: The GW signal for varying $p$ in both matter domination and
radiation domination. Here, we take $T_{\rm CO}=10^{4}$ GeV, $\sigma=10^{-1}$
and, for radiation domination, $\Omega_{\rm rad}/\Omega_{\psi}=100$. We assume
a monochromatic mass function, such that all the overdensities have mass
$M_{\psi}=10^{-12}\ M_{\odot}$. Sensitivity curves for LIGO Amaro-Seoane  _et
al._ (2017), DECIGO Seto  _et al._ (2001); Sato  _et al._ (2017) and BBO
Corbin and Cornish (2006).

$1+z_{\rm col}=\frac{1+z_{\rm CO}}{t_{\rm col}^{s}}\left(\frac{3s^{2}}{8\pi
G\rho_{\rm CO}}\right)^{s/2}$ (24)

where $z_{\rm CO}$ is the redshift at $T_{\rm CO}$ and $\rho_{\rm
CO}=\pi^{2}g_{*}T_{\rm CO}^{4}/30$. We may also use these quantities to
express the $q$, the comoving radius of the perturbation at the time of
horizon entry, $t_{q}$, Dalianis and Kouvaris (2021)

$\begin{split}&q^{-1}(M,T_{\rm CO})=\left(\frac{3M(1+z_{\rm
CO})^{3}}{4\pi\rho_{\rm CO}}\right)^{-1/3}\\\ &\simeq 1.2\times 10^{10}\
\text{Mpc}^{-1}\left(\frac{T_{\rm CO}}{10^{10}\ {\rm
GeV}}\right)^{\frac{1}{3}}\left(\frac{M}{M_{\odot}}\right)^{-\frac{1}{3}}\end{split}$
(25)

where $z_{\rm CO}$ was determined using conservation of entropy. Our selection
of $T_{\rm CO}$ will be informed by the Heaviside function in (23). In
particular,

$\begin{split}T_{\rm CO}&\lesssim\frac{0.2\ {\rm
GeV}}{\left(1+\frac{\Omega_{\rm
rad}(t_{q})}{\Omega_{\psi}(t_{q})}\right)^{1/2}}\left(\frac{M_{\psi}}{M_{\odot}}\right)^{-1/2}\left(\frac{g_{*}}{106.75}\right)^{-1/4}\times\\\
&\quad\times\left(\frac{\alpha\sigma}{\alpha+\beta+\gamma}\right)^{1/(2ps)}.\end{split}$
(26)

In the RD case, $\Omega_{\rm rad}\gg\Omega_{\psi}$ greatly restrict possible
values of $T_{\rm CO}$.

Fundamental to (23) is the mass quadruple moment. We have included a full
derivation of $Q_{ij}$ in the supplementary material. For simplicity, we will
simply present the Fourier transform:

$\tilde{\dddot{Q}}_{ii}(\omega)=\frac{1}{2\pi}\int_{t_{1}}^{t_{2}}\dddot{Q}_{ii}(t)e^{-i\omega
t}\ dt\\\\[7.11317pt]
=\frac{2M}{30\pi(p+1)s}\frac{t_{q}^{2(1-s)}}{t_{\max}^{2ps}}\left[\mathcal{A}_{ps}v_{i}({\bm{\alpha}})\
t_{\max}^{ps}\ \mathcal{I}_{1}(t_{1},t_{2},\omega)\right.\\\
\left.-\mathcal{B}_{ps}w_{i}({\bm{\alpha}})\
\mathcal{I}_{2}(t_{1},t_{2},\omega)\right]$ (27)

where

$\begin{split}\mathcal{I}_{1}(t_{1},t_{2},\omega)&=\int_{t_{1}}^{t_{2}}t^{(p+2)s-3}e^{-i\omega
t}\ dt,\\\\[7.11317pt]
\mathcal{I}_{2}(t_{1},t_{2},\omega)&=\int_{t_{1}}^{t_{2}}t^{2(p+1)s-3}e^{-i\omega
t}\ dt.\end{split}$ (28)

and the expressions for the constants $\mathcal{A}_{ps}$, $\mathcal{B}_{ps}$
and the vectors ${\bf v}({\bm{\alpha}})$ and ${\bf w}({\bm{\alpha}})$ are
defined alongside their derivation in the supplementary material. The square
of (27) with $(t_{1},t_{2})\to(t_{\max},t_{\rm col})$, and $\omega\to 2\pi
f_{0}(1+z_{\rm col})$, completes our derivation of $\Omega_{\rm
GW}(t_{0},f_{0})$.

Estimating $\Omega_{\rm GW}(t_{0},f_{0})$ requires numerical integration over
$(\alpha,\beta,\gamma)$ subject to the restrictions specified in (20). To
begin, the left side of Fig. 1 illustrates the enhanced amplitude resulting
from varying values of $p$ in an IMD era. The enhancement of the peak signal
$\Omega_{\rm GW}^{\rm peak}\equiv\Omega_{\rm GW}(t_{0},f_{\rm peak})$ is quite
dramatic, even for modest values of $p>1$, with $p=1$ being consistent with
the result obtained in Ref. Dalianis and Kouvaris (2021). We note that the
sizable enhancement may conflict with indirect bounds derived for the
integrated quantity, $\Omega_{\rm GW}\equiv\int\Omega_{\rm GW}(f)\ d\ln f$
Lasky  _et al._ (2016); Abbott  _et al._ (2017); Caldwell  _et al._ (2022).

Given the complexity of $\Omega_{\rm GW}(t_{0},f_{0})$, we can determined an
approximate parameterization of $\Omega_{\rm GW}^{\rm peak}$ and $f_{\rm
peak}$ after computing the full signal for various combinations of parameters.
We select a parameterization of the form,

$\begin{split}\Omega_{\rm GW}^{\rm peak}&\approx\Omega_{\rm GW,0}^{\rm peak}\
\sigma^{\delta_{A}/p}\left(\frac{p}{3}\right)^{\alpha_{A}}\left(\frac{T_{\rm
CO}}{10^{4}\ {\rm GeV}}\right)^{\beta_{A}}\times\\\\[7.11317pt]
&\times\left(\frac{M}{10^{-12}\
M_{\odot}}\right)^{\gamma_{A}}\left(\frac{\Omega_{\rm
rad}(t_{q})/\Omega_{\psi}(t_{q})}{100}\right)^{\varepsilon_{A}}\end{split}$
(29)

and

$\begin{split}f_{\rm peak}&\approx f_{\rm peak,0}\
\left(\frac{p}{3}\right)^{\alpha_{f}}\left(\frac{T_{\rm CO}}{10^{4}\ {\rm
GeV}}\right)^{\beta_{f}}\\\\[7.11317pt] &\times\left(\frac{M}{10^{-12}\
M_{\odot}}\right)^{\gamma_{f}}\left(\frac{\Omega_{\rm
rad}(t_{q})/\Omega_{\psi}(t_{q})}{100}\right)^{\varepsilon_{f}}\end{split}$
(30)

where the coefficients are for each formula are presented in Table 1 and 2
respectively.

| $\Omega_{\rm GW,0}^{\rm peak}$ | $\alpha_{A}$ | $\beta_{A}$ | $\gamma_{A}$ | $\delta_{A}$ | $\varepsilon_{A}$
---|---|---|---|---|---|---
MD | $6.8\times 10^{-8}$ | 5 | $4/3$ | $2/3$ | 1 | 0
RD | $2.2\times 10^{-9}$ | 7 | $1$ | $1/2$ | $p/5$ | $-1/2$
Table 1: Coefficients for parmeterization of $\Omega_{\rm GW}^{\rm peak}$. | $f_{\rm peak,0}\ [{\rm Hz}]$ | $\alpha_{f}$ | $\beta_{f}$ | $\gamma_{f}$ | $\varepsilon_{f}$
---|---|---|---|---|---
MD | $1.1\times 10^{-2}$ | 5/6 | $1/3$ | $-1/3$ | 0
RD | $2.0\times 10^{-3}$ | 2/5 | $0$ | $-1/2$ | $-1/2$
Table 2: Coefficients for parmeterization of peak frequency $f_{\rm peak}$.

For a RD era, illustrated on the right of Fig. 1 there is a notable
suppression in amplitude compared to the IMD case. We also note in this new
context, $\delta\propto a^{p}$, $p\geq 1$ always exceeds the logarithmic
growth enabled by gravity only during a radiation era. In either era, the peak
frequency at emission corresponds roughly to $f_{\rm peak}^{\rm emit}\sim
1/t_{\max}$.

The GW spectrum calculated scales as

$\Omega_{\rm GW}(t_{0},f_{0})\propto\begin{cases}f_{0}&\text{for}\ f_{0}\ll
1\\\ f_{0}^{-1}&\text{for}\ f_{0}\gg 1\\\ \end{cases}$ (31)

as also demonstrated in Dalianis and Kouvaris (2021), however, there is reason
to potentially doubt this scaling far from the peak amplitude. By examining a
finite time interval, $[t_{\max},t_{\rm col}]$, we essentially are windowing
the quadrupole moment with a step function. The net result leads to a Fourier
transform proportional to a constant in the IR, and scaling as $1/\omega^{2}$
in the UV regardless of the functional form of $Q_{ij}(t)$. Furthermore, this
result is inconsistent with the expectation that the GW amplitude should scale
as $\omega^{3}$ for small frequencies during emission during a RD era Cai  _et
al._ (2020). While the solution near the peak is reliable, in future, it would
be advantageous to investigate scaling behavior through alternative methods
including numerical simulations and cosmological perturbation theory.

In this work, we investigated the gravitational waves produced from the
evolution of non-spherical density perturbations in the presence of a long-
range scalar force. We consider the scenario where heavy, non-relativistic
fermions interact via an attractive Yukawa interaction. This force, which can
be significantly stronger than gravity, causes structure to form rapidly. The
subsequent collapse of these aspherical overdensities, which could lead to
PBHs, has a non-vanishing quadrupole moment allowing for the emission of GWs.
This signal, which fits inside the sensitivity curve of LISA, can be
significantly stronger than the gravity-only scenario and would be
distinguishable from other signals which occur during a RD era Mollerach  _et
al._ (2004); Ananda  _et al._ (2007); Baumann  _et al._ (2007); De Luca  _et
al._ (2019).

###### Acknowledgements.

We thank I. Dalianis, G. Domènech and D. Inman, for helpful discussions. This
work was supported by the U.S. Department of Energy (DOE) Grant No. DE-
SC0009937. The work of A.K. was also supported by the World Premier
International Research Center Initiative (WPI), MEXT, Japan and by Japan
Society for the Promotion of Science (JSPS) KAKENHI grant No. JP20H05853. MS
is supported by the JSPS KAKENHI grant Nos. 19H01895, 20H04727, and 20H05853.

## References

  * Amendola  _et al._ (2018) L. Amendola, J. Rubio, and C. Wetterich, Phys. Rev. D 97, 081302 (2018), arXiv:1711.09915 [astro-ph.CO] .
  * Savastano  _et al._ (2019) S. Savastano, L. Amendola, J. Rubio, and C. Wetterich, Phys. Rev. D 100, 083518 (2019), arXiv:1906.05300 [astro-ph.CO] .
  * Casas  _et al._ (2016) S. Casas, V. Pettorino, and C. Wetterich, Phys. Rev. D 94, 103518 (2016), arXiv:1608.02358 [astro-ph.CO] .
  * Flores and Kusenko (2021a) M. M. Flores and A. Kusenko, Phys. Rev. Lett. 126, 041101 (2021a), arXiv:2008.12456 [astro-ph.CO] .
  * Domènech and Sasaki (2021) G. Domènech and M. Sasaki, JCAP 06, 030 (2021), arXiv:2104.05271 [hep-th] .
  * Flores and Kusenko (2021b) M. M. Flores and A. Kusenko, (2021b), arXiv:2108.08416 [hep-ph] .
  * Wise and Zhang (2014) M. B. Wise and Y. Zhang, Phys. Rev. D 90, 055030 (2014), [Erratum: Phys.Rev.D 91, 039907 (2015)], arXiv:1407.4121 [hep-ph] .
  * Gresham  _et al._ (2018) M. I. Gresham, H. K. Lou, and K. M. Zurek, Phys. Rev. D 97, 036003 (2018), arXiv:1707.02316 [hep-ph] .
  * Gresham and Zurek (2019) M. I. Gresham and K. M. Zurek, Phys. Rev. D 99, 083008 (2019), arXiv:1809.08254 [astro-ph.CO] .
  * Dalianis and Kouvaris (2021) I. Dalianis and C. Kouvaris, JCAP 07, 046 (2021), arXiv:2012.09255 [astro-ph.CO] .
  * Jedamzik  _et al._ (2010) K. Jedamzik, M. Lemoine, and J. Martin, JCAP 04, 021 (2010), arXiv:1002.3278 [astro-ph.CO] .
  * Domènech  _et al._ (2022) G. Domènech, S. Passaglia, and S. Renaux-Petel, JCAP 03, 023 (2022), arXiv:2112.10163 [astro-ph.CO] .
  * Mukhanov (2005) V. Mukhanov,  _Physical Foundations of Cosmology_ (Cambridge University Press, Oxford, 2005).
  * Gradwohl and Frieman (1992) B.-A. Gradwohl and J. A. Frieman, Astrophys. J. 398, 407 (1992).
  * Gubser and Peebles (2004) S. S. Gubser and P. Peebles, Phys. Rev. D 70, 123510 (2004), arXiv:hep-th/0402225 .
  * Nusser  _et al._ (2005) A. Nusser, S. S. Gubser, and P. Peebles, Phys. Rev. D 71, 083505 (2005), arXiv:astro-ph/0412586 .
  * Zel’dovich (1970) Y. B. Zel’dovich, Astron. Astrophys. 5, 84 (1970).
  * Thorne (1972) K. S. Thorne, “Non-spherical Gravitational Collapse: A Short Review,” (1972).
  * Misner  _et al._ (1973) C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman, San Francisco, 1973).
  * Flores and Kusenko (2021c) M. M. Flores and A. Kusenko, Phys. Rev. D 104, 063008 (2021c), arXiv:2106.03237 [astro-ph.CO] .
  * Doroshkevich (1970) A. Doroshkevich, Astrophysics 6, 320 (1970).
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Seto  _et al._ (2001) N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001), arXiv:astro-ph/0108011 .
  * Sato  _et al._ (2017) S. Sato _et al._ , J. Phys. Conf. Ser. 840, 012010 (2017).
  * Corbin and Cornish (2006) V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435 (2006), arXiv:gr-qc/0512039 .
  * Lasky  _et al._ (2016) P. D. Lasky _et al._ , Phys. Rev. X 6, 011035 (2016), arXiv:1511.05994 [astro-ph.CO] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 118, 121101 (2017), [Erratum: Phys.Rev.Lett. 119, 029901 (2017)], arXiv:1612.02029 [gr-qc] .
  * Caldwell  _et al._ (2022) R. Caldwell _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07972 [gr-qc] .
  * Cai  _et al._ (2020) R.-G. Cai, S. Pi, and M. Sasaki, Phys. Rev. D 102, 083528 (2020), arXiv:1909.13728 [astro-ph.CO] .
  * Mollerach  _et al._ (2004) S. Mollerach, D. Harari, and S. Matarrese, Phys. Rev. D 69, 063002 (2004), arXiv:astro-ph/0310711 .
  * Ananda  _et al._ (2007) K. N. Ananda, C. Clarkson, and D. Wands, Phys. Rev. D 75, 123518 (2007), arXiv:gr-qc/0612013 .
  * Baumann  _et al._ (2007) D. Baumann, P. J. Steinhardt, K. Takahashi, and K. Ichiki, Phys. Rev. D 76, 084019 (2007), arXiv:hep-th/0703290 .
  * De Luca  _et al._ (2019) V. De Luca, V. Desjacques, G. Franciolini, and A. Riotto, JCAP 09, 059 (2019), arXiv:1905.13459 [astro-ph.CO] .
  * Harada  _et al._ (2016) T. Harada, C.-M. Yoo, K. Kohri, K.-i. Nakao, and S. Jhingan, Astrophys. J. 833, 61 (2016), arXiv:1609.01588 [astro-ph.CO] .

##  Appendix A Supplementary Material

###  A.1 Derivation of time scales

We begin by first deriving $t_{\max}$. To do so, will utilize the continuity
of $\alpha$, $\beta$ and $\gamma$ so that we can locally take the coordinates
$q_{i}$ to be Dalianis and Kouvaris (2021); Harada  _et al._ (2016)

$\begin{cases}r_{1}=(a-\alpha b)q_{1}\\\ r_{2}=(a-\beta b)q_{2}\\\
r_{3}=(a-\gamma b)q_{3}.\end{cases}$ (32)

By definition, $t_{\max}$ is the time that satisfies

$\dot{r}_{1}(t_{\rm max})=0$ (33)

which implies $\dot{b}(t_{\rm max})/\dot{a}(t_{\rm max})=1/\alpha$ and
therefore,

$\frac{b(t_{\max})}{a(t_{\max})}=\frac{1}{(p+1)\alpha}.$ (34)

Combining (12) evaluated at $t_{q}$ and (34), we find that

$t_{\max}=\left(\frac{\alpha+\beta+\gamma}{(p+1)\alpha\delta_{L}(t_{q})}\right)^{1/ps}t_{q}.$
(35)

Next, we focus our attention on $t_{\rm col}$. Collapse takes place when
$r_{1}(t_{\rm col})=0$. This naturally implies $b(t_{\rm col})/a(t_{\rm
col})=1/\alpha$. Using this result in combination with (34) implies

$t_{\rm col}=(p+1)^{1/ps}t_{\max}.$ (36)

###  A.2 Derivation of the quadrupole moment

We can express $Q_{ij}$ in terms of the moment of inertia,

$Q_{ij}(t)=-I_{ij}(t)+\frac{1}{3}\delta_{ij}{\rm Tr}I(t).$ (37)

In the reference frame of the principle axes, the moment of inertia is

$I_{ij}=\frac{M}{5}\begin{pmatrix}r_{2}^{2}+r_{3}^{2}&0&0\\\
0&r_{1}^{2}+r_{3}^{2}&0\\\ 0&0&r_{1}^{2}+r_{2}^{2}\\\ \end{pmatrix}.$ (38)

Using (32) and (34) we find,

$\begin{split}r_{1}(t)&=\frac{1}{s}t_{q}^{1-s}t^{s}\left(1-\frac{1}{(p+1)}\left(\frac{t}{t_{\max}}\right)^{ps}\right),\\\\[7.11317pt]
r_{2}(t)&=\frac{1}{s}t_{q}^{1-s}t^{s}\left(1-\frac{\beta}{(p+1)\alpha}\left(\frac{t}{t_{\max}}\right)^{ps}\right),\\\\[7.11317pt]
r_{3}(t)&=\frac{1}{s}t_{q}^{1-s}t^{s}\left(1-\frac{\gamma}{(p+1)\alpha}\left(\frac{t}{t_{\max}}\right)^{ps}\right).\end{split}$
(39)

Combining (37) - (39) we obtain

$\begin{split}\dddot{Q}_{ii}(t)=&\frac{2M}{15(p+1)s}\frac{t_{q}^{2(1-s)}}{t_{\max}^{2ps}}t^{(p+2)s-3}\\\\[7.11317pt]
&\hskip
28.45274pt\times\left[\mathcal{A}_{ps}v_{i}({\bm{\alpha}})t_{\max}^{ps}-\mathcal{B}_{ps}w_{i}({\bm{\alpha}})t^{ps}\right].\end{split}$
(40)

The constants $\mathcal{A}_{ps}$ and $\mathcal{B}_{ps}$ utilized in (27) are

$\begin{split}\mathcal{A}_{ps}&=(p+2)(s(p+2)-2)(s(p+2)-1)\\\\[7.11317pt]
\mathcal{B}_{ps}&=2(ps+s-1)(2s(p+1)-1).\end{split}$ (41)

we note that for $p=1$ and $s=2/3$, $\mathcal{A}_{ps}=0$, thus simplifying
much of the expressions above. The vector expressions ${\bf v}({\bm{\alpha}})$
and ${\bf w}({\bm{\alpha}})$ are defined as,

$\begin{split}{\bf
v}({\bm{\alpha}})&=\left\langle-2+\frac{\beta}{\alpha}+\frac{\gamma}{\alpha},\
1-\frac{2\beta}{\alpha}+\frac{\gamma}{\alpha},\
1+\frac{\beta}{\alpha}-\frac{2\gamma}{\alpha}\right\rangle,\\\\[7.11317pt]
{\bf
w}({\bm{\alpha}})&=\left\langle-2+\left(\frac{\beta}{\alpha}\right)^{2}+\left(\frac{\gamma}{\alpha}\right)^{2},\
1-2\left(\frac{\beta}{\alpha}\right)^{2}+\left(\frac{\gamma}{\alpha}\right)^{2},1+\left(\frac{\beta}{\alpha}\right)^{2}-2\left(\frac{\gamma}{\alpha}\right)^{2}.\right\rangle\end{split}$
(42)

Lastly, we may also rewrite the integral expressions (28) in terms of the the
generalized exponential integral

$E_{n}(x)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\ dt.$ (43)

Using this definition integrals (28) read,

$\begin{split}\mathcal{I}_{1}(t_{1},t_{2},\omega)&=\frac{1}{t_{1}^{2-(p+2)s}}\left[E_{3-(p+2)s}(i\omega
t_{1})-\left(\frac{t_{1}}{t_{2}}\right)^{2-(p+2)s}E_{3-(p+2)s}(i\omega
t_{2})\right]\\\\[7.11317pt]
\mathcal{I}_{2}(t_{1},t_{2},\omega)&=\frac{1}{t_{1}^{2-2(p+1)s}}\left[E_{3-2(p+1)s}(i\omega
t_{1})-\left(\frac{t_{1}}{t_{2}}\right)^{2-2(p+1)s}E_{3-2(p+1)s}(i\omega
t_{2})\right].\end{split}$ (44)