# A deep learning approach to resonant light transmission through single
subwavelength apertures

David Alonso-González  Departamento de Física Teórica de la Materia Condensada
and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
E-28049 Madrid, Spain  Michel Frising  Departamento de Física de la Materia
Condensada and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma
de Madrid, E-28049 Madrid, Spain  Ferry Prins  Departamento de Física de la
Materia Condensada and Condensed Matter Physics Center (IFIMAC), Universidad
Autónoma de Madrid, E-28049 Madrid, Spain  Jorge Bravo-Abad  Departamento de
Física Teórica de la Materia Condensada and Condensed Matter Physics Center
(IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain

###### Abstract

Resonant transmission of light is a surface-wave assisted phenomenon that
enables funneling light through subwavelength apertures milled in otherwise
opaque metallic screens. In this work, we introduce a deep learning approach
to efficiently compute and design the optical response of a single
subwavelength slit perforated in a metallic screen and surrounded by periodic
arrangements of indentations. First, we show that a semi-analytical framework
based on a coupled-mode theory formalism is a robust and efficient method to
generate the large training datasets required in the proposed approach.
Second, we discuss how simple, densely connected artificial neural networks
can accurately learn the mapping from the geometrical parameters defining the
topology of the system to its corresponding transmission spectrum. Finally, we
report on a deep learning tandem architecture able to perform inverse design
tasks for the considered class of systems. We expect this work to stimulate
further work on the application of deep learning to the analysis of light-
matter interaction in nanostructured metallic films.

The advance of the technological abilities for structuring metals in the
micro- and nanometer length scales led in 1998 to the experimental discovery
that the optical transmission through an array of subwavelength holes milled
in a metallic screen can, for certain wavelengths, be significantly boosted
with respect to what was predicted for independent holes eot1 . Later, it was
shown that a single aperture can also feature enhanced transmission, along
with remarkable beaming effects, when surrounded by periodic arrangements of
surface corrugations eot2 . The underlying physical origin of both phenomena
relies on the excitation of electromagnetic resonances (surface plasmons)
running along the corrugated metallic surfaces eot3 ; eot4 . Importantly, the
transmission can be precisely engineered by varying the geometry of the
surface corrugations, allowing for spectrally, directionally, and polarization
selective transmission properties laux2008 ; jun2011 ; davis2017 ; deleo2017 .

On a different front, during the last decades deep learning (DL), a sub-area
of machine learning goodfellow , has drawn a great deal of attention from the
scientific community due to its ability to deal with complex data-driven
problems cirac . DL techniques are based on _training_ artificial neural
networks (ANNs), a collection of nonlinear units or nodes called artificial
neurons murphy . Each artificial neuron receives a _signal_ , processes it,
and then transmits a new output to other neurons connected to it. The signal
goes from an input layer formed by a set of neurons to an output one, usually
after traversing several layers of neurons (so-called _hidden layers_).

In this work, we report on the application of ANNs to the simulation of the
optical response of a single subwavelength aperture surrounded by periodic
arrangements of indentations and perforated in an opaque metallic screen. In
our approach, we first use a semi-analytical theoretical framework, based on a
well-known coupled-mode theory, to generate large enough data sets than can be
used to train the studied ANNs. Then, we show that ANNs can indeed be used to
accurately obtain the non-linear mapping between the topology of the studied
class of nanophotonic structures and their corresponding transmission spectra.
In addition, we discuss a deep-learning based inverse design algorithm that
enables finding the geometrical parameters of the system that yield a desired
optical response.

We focus on the structure sketched in Fig. 1, which have been extensively
analyzed in the context of resonant transmission phenomena eot4 . It is based
on a single slit (of width $b$), perforated on a metallic film (of thickness
$W$) flanked to the left and right by finite arrays of one-dimensional
rectangular indentations. Each of these arrays is characterized by its
corresponding periodicity ($d_{1}$ and $d_{2}$ label the periodicity of the
left and right arrays, respectively), and the width and depth of each
indentation. In what follows, we assume that all indentations have the same
width ($a$) but different depths in the left and right arrays ($h_{1}$ and
$h_{2}$, respectively). We also assume the whole system is embedded in air
(refractive index $n=1$). The number and specific degrees of freedom chosen
here to describe the geometry allows having a moderate number of free
parameters, while at the same time capturing with their variation most of the
phenomenology present in this class of devices.

Figure 1: Schematics of the system under study, along with the geometrical
parameters defining it. Yellow arrows on top represent the plane wave
illumination incident normally onto the system.

One of the key ingredients of any deep learning approach is the dataset used
to train the algorithmic architecture. For studying our specific problem, we
have generated the dataset using a coupled mode theory (CMT) framework. This
theoretical framework has been widely applied to the studied systems and it is
well-known for its accuracy and computational efficiency eot4 . The latter
aspect is crucial for our purposes, as will be key for having access to large
training sets required in the proposed approach.

We consider the case of $p$-polarized incident light, for which the most
apparent resonant transmission effects have been reported. Within a CMT
framework, electromagnetic (EM) fields are expanded by a set of plane waves in
air and by a modal expansion inside the slit and indentations. As we are
mainly interested in subwavelength apertures, we consider that only the
fundamental propagating eigenmode is excited. We assume that the metal behaves
as a perfect metal, expelling the electric field from its bulk. This
assumption, ideally suited for applying a CMT as the one used here, captures
in a semi-quantitative fashion most enhanced transmission phenomena in the
optical regime eot3 . Within this formalism, the total transmittance of the
studied setup can be computed through the following compact expression:
$T=\textrm{Im}\left\\{E^{I}E^{O*}\right\\}/\sin(kW)$, where $E^{I}$ and
$E^{O}$ are complex amplitudes related the $x$-component of the electric field
at illuminated and non-illuminated the exit of the slit, respectively, and
$k=2\pi/\lambda$, with $\lambda$ being the wavelength in vacuum of the
incident light (see Ref. eot4 , and references therein, for details on these
magnitudes and a description of the re-illumination effects between the
indentations and the slit).

Figure 2: (a) Evolution of the cost function with epochs in the optimal ANN,
as obtained for different sizes of the dataset, ranging from 606,481 (blue
line) to 101,140 (red line). (b) Predictions of the ANN for a random set of
six transmission spectra that were contained in the test set (i.e., not used
to train the network). The values of the parameters corresponding to each case
(in nm) are indicated in the legend following the notation [$d_{1}$ $h_{1}$
$d_{2}$ $h_{2}$]. Solid lines and dashed correspond, respectively, to the
results computed with a coupled-mode theory formalism and the predictions from
the ANN.

First, we trained our neural network to accurately predict the transmission
spectra of the studied class of systems. In order to generate the
corresponding dataset, we varied the periodicity ($d_{1}$ and $d_{2}$) and the
depth ($h_{1}$ and $h_{2}$) of both left and right sides of the illuminated
part of the structure. In all cases, we fixed the number of indentations to
the left and right of the slit to $N=10$ (i.e., we have a total of 20
indentations surrounding the slit). Following Ref. soljacic , the initial
architecture of the network consisted of 3 hidden layers with 250 neurons
each. The input layer has 4 neurons corresponding to the four input parameters
($d_{1}$, $d_{2}$, $h_{1}$ and $h_{2}$), while the output layer has 81
neurons, corresponding to a sampling of the transmission spectra at 81
different wavelength values, equally spaced between $400$ nm and $800$ nm. Our
dataset is therefore composed by sets of the four parameters (playing the role
of _features_ in our deep learning approach) and their corresponding
transmission spectra (these spectra being the corresponding _labels_).

The detailed count of samples forming this dataset is obtained as follows. We
vary independently the periodicities of the left and right arrays of
indentations from $400$ nm to $800$ nm with a step of $10$ nm, yielding 41
possible values for the periodicity of each array. We also vary independently
the depths of the indentations of the left and right arrays, from $5$ nm to
$95$ nm with a step of $5$ nm, producing 19 possible configurations for each
side of the slit. Therefore, we have a total of $(41\times 19)^{2}=606,841$
spectra in the complete dataset. From this dataset, we take a 80% as the
training set, and a 10% each for the validation and the test sets.

Next, we performed a systematic set of calculations to find the optimal
hyperparameters of the described ANN. In all these studies, to quantify the
difference between the ANN’s predictions and the CMT results, during the
training process we used a cost function based on the mean average error (MAE)
of all sampling points assumed in the spectra. Regarding the architecture of
the network, we found that our initial guess (including 3 hidden layers and
250 neurons per layer) leads to an optimal performance of the network. On the
other hand, from all the activation functions, a _Sigmoid_ was the one
yielding the best accuracy (see SI  for details on this analysis). The best
performing optimizer among the series we tried resulted to be _Adam_. For that
optimizer we tried several learning rates, $\eta$. We checked that values of
$\eta$ lower than $10^{-4}$ show a very good performance, and found a
particularly good performance for $\eta=3\times 10^{-5}$. In this study of
hyperparameters, the variation in the efficiency of the ANN’s training with
the size of the training set is of particular interest. Figure 2(a) summarizes
the numerical results for the evolution of the cost function with the number
of epochs, for different sizes of the training set. As seen, as the size of
the training set decreases, the training process becomes slower and for the
maximum number of epochs considered in this work (5,000 epochs) the cost
function increases significantly when considering training sets smaller than
approximately 200,000 samples. This result shows that the ability to access
large enough datasets is a key ingredient to enable the proposed approach.

Figure 2(b) shows the predictions of the ANN featuring the optimal
hyperparameter set and for the largest number of epochs and training set sizes
considered in this work (5,000 and 606,841, respectively). A set of six
randomly selected transmission spectra from the test set is considered for
illustration. As observed, for the displayed set, the predictions of the ANN
(solid lines) show very good agreement with the results obtained with the CMT
theory (dashed lines), to the extent that in most cases the results from both
approaches are indistinguishable to the naked eye. A more quantitative way to
assess the performance of the network can be obtained through the evaluation
of the average value of the cost function over the main data subsets involved
in this problem (namely, the training, validation and test sets). Table I
summarizes the corresponding results, showing that an average error of 0.42%
is obtained for the three datasets (note that the fact that the validation and
test set features the same value of the cost function also shows that no
overfitting is occurring in this case).

|  | Average of the cost function (%)
---|---|---
Epochs | Samples | Train | Validation | Test
5000 | 606841 | 0.42 | 0.42 | 0.42
Table 1: Values of the cost function computed for the main data subsets of the
analyzed deep learning approach. Figure 3: (a) Schematics of the tandem
artificial neural network used for inverse design. Yellow circles to the left
and right represent the input and output layers, respectively (both holding a
sampling of 81 points of the spectra). Gray circles represent the hidden
layers of the proposed ANN architecture. Blue circles correspond to the inner
artificial neurons holding values that define the geometry of the structure.
The area in green marks the part of the deep learning system that corresponds
to a pretrained feed-forward neural network (FFNN). Trainable connections
among neurons are displayed in purple. (b) Comparison between a targeted
spectrum and the prediction of the ANN for than spectrum, obtained by using
the deep learning architecture displayed in (a).

Now, we turn to the _inverse design_ of nanophotonic devices so ; fan . The
main idea underlying deep learning for inverse design can be explained in
simple terms: it consists of ANNs able to design nanophotonic systems that can
match a predefined _targeted_ optical response. In this work, inspired by
Refs. so ; fan  (and references therein), we consider the tandem architecture
shown in Fig. 3(a), which we have specifically tailored to perform inverse
design tasks for the studied kind of metallic structures. The artificial
neurons of the input layer (shown in yellow in Fig. 3(a)) hold the sampling of
the targeted transmission spectrum (as in the previous ANN we assume a
81-point sampling of the spectrum). The output layer has a similar structure,
but in this case the neurons hold the sampling of the spectrum produced by the
whole deep learning system. The overall loss function of this deep learning
system is set up to minimize the difference between the targeted spectrum and
the one obtained from the output layer (see details in SI ). Right at the
middle of the system we have inserted a small hidden layer, consisting in only
4 neurons (shown in blue in Fig. 3(a)). Importantly, those neurons hold the
values of the geometrical parameters characterizing the structure (the values
of the set $[d_{1}\leavevmode\nobreak\ h_{1}\leavevmode\nobreak\
d_{2}\leavevmode\nobreak\ h_{2}]$) that, when operating in inference mode, the
proposed deep learning system predicts as optimal for achieving the targeted
response. The connection of this inner layer with the input and output layers
(to the left and right, respectively, in the sketch of Fig. 3(a)) is done
through the stacking of hidden layers (represented as columns of grey
circles). We assume the same topology for these hidden layers as the one used
in Fig. 2, i.e., 3 hidden layers of 250 artificial neurons are set to the left
and right of the inner layer.

To demonstrate the performance of the tandem architecture shown in Fig. 3(a),
we have selected as target spectrum the one shown in Fig. 3(b) (blue line). It
correspond to a physical setting, not included in the training set, defined by
the geometrical parameters: $[d_{1}\leavevmode\nobreak\
h_{1}\leavevmode\nobreak\ d_{2}\leavevmode\nobreak\
h_{2}]=[470\leavevmode\nobreak\ 75\leavevmode\nobreak\ 560\leavevmode\nobreak\
95]$ nm. As observed, the spectrum features two resonant transmission peaks at
very different wavelengths that are close in magnitude and width. We believe
that reaching a solution for a spectrum of these characteristics, which
combines simultaneously many different aspects of the physical mechanisms
underlying resonant transmission in the considered structures eot4 , is a
demanding proof-of-principle task that illustrates well the inverse design
capabilities of the system sketched in Fig. 3(a).

Remarkably, we have found that to deploy the above described tandem
architecture as efficient inverse design engine, it is not necessary to train
the whole deep learning system shown in Fig. 3(a). An alternative approach
consists in _transferring_ to the layers placed to the right of the central
inner layer of the tandem structure (green area in Fig. 3(a)) the weights of
the optimized predictive ANN discussed in the first part of this work (labeled
as _feed-forward neural network, FFNN_ in Fig. 3(a)). In addition, we transfer
to the layers located to the left of the central inner layer the same weights,
but with an important difference: we do not transfer the ones corresponding to
the connections between last hidden layer of that part of the architecture and
the inner central layer (these connections are shown in purple in Fig. 3(a)).
Following this route it is possible to avoid the large computational load
associated to a detailed training of the whole system (see SI  for the full
details of the implementation of this numerical process).

By running the well-known backpropagation algorithm goodfellow  in the above
discussed setup, we obtain that the system shown in Fig. 3(a) is able to
predict that the targeted spectrum corresponds to the set of geometrical
parameters $[469.96\leavevmode\nobreak\ \leavevmode\nobreak\
75.25\leavevmode\nobreak\ \leavevmode\nobreak\ 559.96\leavevmode\nobreak\
\leavevmode\nobreak\ 95.18]$ nm, which are very close to the ones associated
to the targeted spectrum ($[470\leavevmode\nobreak\ 75\leavevmode\nobreak\
560\leavevmode\nobreak\ 95]$ nm). Figure 3(b) shows the spectrum, also
produced by the same tandem structure, for the predicted optimal set of
parameters (orange dashed line). The very good agreement between predicted and
targeted spectra demonstrates the high degree of accuracy that can be obtained
with the proposed inverse design approach.

To sum up, in this work we have studied the application of deep learning
techniques to the simulation and design of the optical response of a single
subwavelength aperture milled in an opaque metallic screen and surrounded by
periodic arrangements of indentations. In particular, we have employed a
densely connected ANN to accurately obtain the non-linear mapping between the
topology of this class of structures and their corresponding transmission
spectra. To do that, we have used a semi-analytical coupled mode theory that
has allowed us to access the large dataset needed in the training process of
our ANN. We have found very good agreement between the transmission spectra
predicted by the proposed ANN and the spectra obtained with the coupled-mode
theory. We have also explored an important application of these kind of
methods, namely, the inverse design of the system. Specifically, we have
proposed a tailored ANN tandem architecture as solution for inverse design
problems in this context. We believe the present study can contribute to a
wider adoption of deep learning techniques as efficient tools for the analysis
of light-matter interaction in metallic films structured at the subwavelength
scale.

## References

  * (1) T.W. Ebbesen, H.J. Lezec, H.F. Ghaemi, T. Thio, P.A. Wolff, Extraordinary optical transmission through subwavelength hole arrays, Nature 391, 667 (1998).
  * (2) H.J. Lezec, A. Degiron, E. Devaux, R.A. Linke, L. Martín-Moreno, F.J. García-Vidal and T.W. Ebbesen, Beaming light from a Subwavelength Aperture, Science, 297, 820 (2002).
  * (3) L. Martín-Moreno, F.J. García-Vidal, H.J. Lezec, K.M. Pellerin, T. Thio, J.B. Pendry and T.W. Ebbesen, Theory of extraordinary optical transmission through subwavelength hole arrays, Phys. Rev. Lett. 86, 1114 (2001).
  * (4) F.J. García-Vidal, L. Martín-Moreno, T.W. Ebbesen and L. Kuipers, Light passing through subwavelength apertures, Rev. Mod. Phys. 82, 729 (2010).
  * (5) E. Laux, C. Genet, T. Skauli, T.W. Ebbesen, Plasmonic photon sorters for spectral and polarimetric imaging, Nat. Phot. 2, 161 (2008).
  * (6) Y.C. Jun, K.C.Y. Huang, M.L. Brongersma, Plasmonic beaming and active control over fluorescent emission, Nat. Comm. 2, 283 (2011).
  * (7) M.S. Davis, W. Zhu, T. Xu, J.K. Lee, H.J. Lezec, A. Agrawal, Aperiodic nanoplasmonic devices for directional colour filtering and sensing, Nat. Comm. 8, 1347 (2017) .
  * (8) E. De Leo, A. Cocina, P. Tiwari, L.V. Poulikakos, P. Marqués-Gallego, B. le Feber, D.J. Norris, F. Prins, Polarization Multiplexing of Fluorescent Emission Using Multiresonant Plasmonic Antennas, ACS Nano 11, 12167 (2017).
  * (9) I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, The MIT Press (Cambridge, 2016).
  * (10) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborova, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019).
  * (11) K.P. Murphy, Machine Learning: A Probabilistic Perspective, The MIT Press (Cambridge, 2012).
  * (12) J. Peurifoy, Y. Shen, L. Jing, Y. Yang, F. Cano-Renteria, B.G. DeLacy, J.D. Joannopoulos, M. Tegmark, M. Soljačić, Nanophotonic particle simulation and inverse design using artificial neural networks, Science Advances 4, eaar4206 (2018).
  * (13) S. So, T. Badloe, J. Noh, J. Bravo-Abad, J. Rho, Deep learning enabled inverse design in nanophotonics, Nanophotonics 9, 1041 (2020).
  * (14) J. Jiang, M. Chen, J.A. Fan, Deep neural networks for the evaluation and design of photonic devices, Nat. Rev. Mat. 1, 22 (2020).
  * (15) See detailed description in the Supplement accompanying this work.


# Gravity-induced entanglement between two massive microscopic particles in
curved spacetime: II. Friedmann- Lemaître-Robertson-Walker universe

Chi Zhang<EMAIL_ADDRESS>Department of Physics, Nanchang University,
Nanchang, 330031, China  Center for Relativistic Astrophysics and High Energy
Physics, Nanchang University, Nanchang, 330031, China  Fu-Wen Shu
<EMAIL_ADDRESS>Corresponding author Department of Physics, Nanchang
University, Nanchang, 330031, China  Center for Relativistic Astrophysics and
High Energy Physics, Nanchang University, Nanchang, 330031, China  GCAP-
CASPER, Physics Department, Baylor University, Waco, Texas 76798-7316, USA
Center for Gravitation and Cosmology, Yangzhou University, Yangzhou, China

###### Abstract

In our previous work Zhang:2023llk , we have explored quantum gravity induced
entanglement of mass (QGEM) in curved spacetime, observing entanglement
formation between particles moving along geodesics in a Schwarzschild
spacetime background. We find that long interaction time induces entanglement,
even for particles with microscopic mass, addressing decoherence concerns. In
this work, we build upon our previous work Zhang:2023llk  by extending our
investigation to a time-dependent spacetime. Specifically, we explore the
entanglement induced by the mutual gravitation of massive particles in the
Friedmann-Lemaître-Robertson-Walker (FLRW) universe. Through calculations of
the phase variation and the QGEM spectrum, our proposed scheme offers a
potential method for observing the formation of entanglement caused by the
quantum gravity of massive particles as they propagate in the FLRW universe.
Consequently, our research provides fresh insights into the field of
entanglement in cosmology.

##  I Introduction

There has been a long-standing debate about whether the gravitational field is
a quantum field, and to date we still have no effective way to give a clear
answer. In the past few years, we have gone through a boom of interest to use
desk-top experiments to explore quantum gravity effects. Many of these
experiment proposals date back to a thought experiment of Feynman feynman
that if one massive particle is put in a superposition of two locations by a
Stern-Gerlach device, then gravitational field would also form a superposition
state. Based on this idea, Bose and Marletto et al r2 ; r3  proposed a plan to
use quantum entanglement to detect quantum gravity effects, which is known as
quantum gravity induced entanglement of mass (QGEM). In the QGEM device, if
these two particles in the spatial superposition state are entangled under the
action of their gravitational fields, then it can be inferred that the
gravitational field has quantum properties based on the principle of quantum
information (see some recent progress Christodoulou:2018cmk ; Capolupo:2019peg
; Belenchia:2018szb ; Carney:2018ofe ; Marshman:2019sne ; Westphal:2020okx ;
Carlesso:2019cuh ; Danielson:2021egj ; Christodoulou:2022mkf ; Cho:2021gvg ;
Matsumura:2020law ; Miki:2020hvg ; Howl:2023xtf ; Yant:2023smr ; Feng:2023krm
; Schut:2023eux ; Fragolino:2023agd ; Li:2022yiy ; Bose:2022uxe ; He:2023hys
). However, due to the harsh experimental parameters and quantum states are
highly susceptible to decoherence, these experimental designs are difficult to
implement in practice.

Very recently, we developed an astronomical version of QGEM in curved
spacetime, focusing on the observation of entanglement formation between pairs
of particles moving along geodesics in a Schwarzschild spacetime background
Zhang:2023llk . We find that the long interaction time between particles can
induce entanglement even for particles with microscopic mass, providing a
possible solution of the issue of the decoherence. In this work Zhang:2023llk
, we also found that pairs of superposed particles exhibit varying degrees of
entanglement depending on different motion parameters. Specifically,
entanglement formation is more likely to occur when the proper motion time is
longer, the geodesic deviation rate is appropriate, and the mass of the galaxy
is smaller.

In the previous work, we considered the QGEM in the Schwarzschild spacetime,
which is static. In this work, we would like to investigate how QGEM evolves
in a time-dependent spacetime. To be more specific, we will investigate
entanglement induced by mutual gravitation between massive particles in the
FLRW universe.

There have been some previous studies on entanglement in cosmology. For
instance, the production of entanglement in quantum fields due to the
expansion of the underlying spacetime was shown in Martin-Martinez:2012chf .
The entanglement entropy of cosmological perturbations in the very early
universe has also been investigated Brahma:2020zpk . In the realm of photon
entanglement, various avenues to utilize entanglement in order to investigate
the nature of the universe have been explored. Such as schemes using entangled
photons to study the origins of cosmic microwave background (CMB) fluctuations
Kiefer:2008ku ; Martin:2007bw ; Martin:2015qta ; Martin:2016tbd  and
investigations involving entangled photons produced from the decay of cosmic
particles to test quantum mechanics Lello:2013bva ; Valentini:2006yj . In
addition, astronomers have even attempted to detect entangled photons that
have traveled different paths through the cosmos and reached Earth
Chen:2017cgw .

However, investigations into the entanglement induced by the gravitational
interaction of massive particles are still lacking. Our current knowledge does
not include an understanding of how QGEM evolves in the FLRW universe. Filling
this gap is one of our main motivations for conducting this research.
Intuitively, pairs of adjacent particles propagating in the FLRW universe are
expected to become entangled due to their gravitational field. However,
determining the specific factors influencing their entanglement requires
detailed calculations. The entanglement induced by quantum gravity, which we
will present in this article, offers a fresh perspective for studying
entanglement in cosmology.

This paper is organized as follows. In Sec. II, we present the overall method
for quantifying the entanglement formation induced by the QGEM mechanism in
the context of FLRW cosmology, and also study the influencing factors of the
entanglement phase. In Sec. III, incorporating more realistic observational
scenarios, we present the characteristic spectral lines of the entanglement
witness as varying functions of the initial motion parameters of the particle
pair and the observed kinetic energy. In Sec. IV, conclusions and discussions
are made.

Throughout this paper, we will adopt the natural units system, $c=G=1$, in
order to simplify the calculations. Except that the physical quantities with
units are given in the SI system of units.

##  II Entanglement generation in FLRW spacetime

For a flat FLRW universe, the metric is

$d{s^{2}}=d{t^{2}}-{a(t)^{2}}(d{r^{2}}+{r^{2}}d{\theta^{2}}+{r^{2}}{\sin^{2}}\theta{d^{2}}\varphi),$
(1)

where $a(t)$ is the scale factor of the FLRW universe and it obeys

$a^{\prime}{\left(t\right)^{2}}=\frac{8}{3}\pi
Ga{\left(t\right)^{2}}\rho\left(t\right).$ (2)

The general density in (2) is given by

$\rho\left(t\right)=\frac{{a_{0}^{3}{\rho_{{M_{0}}}}}}{{a{{\left(t\right)}^{3}}}}+\frac{{a_{0}^{4}{\rho_{{R_{0}}}}}}{{a{{\left(t\right)}^{4}}}}+{\rho_{{v_{0}}}},$
(3)

where

${\rho_{{v_{0}}}}=\frac{{3H_{0}^{2}{\Omega_{\Lambda}}}}{{8G\pi}},{\rho_{{M_{0}}}}=\frac{{3H_{0}^{2}{\Omega_{M}}}}{{8G\pi}},{\rho_{{R_{0}}}}=\frac{{3H_{0}^{2}{\Omega_{R}}}}{{8G\pi}}.$
(4)

Here $\Omega_{\Lambda}$, $\Omega_{M}$ and $\Omega_{R}$ are the density
parameters for dark energy, matter and radiation, respectively. $H_{0}$ is the
Hubble constant.

It was shown in r2 ; r3  that nearby neutral massive particle pairs exhibit
phase growth under their self-interaction of gravity. The phase growth in
general is given by Zhang:2023llk

$\delta\phi=-\frac{{{m_{0}}\delta\tau}}{\hbar}=\frac{{{m_{0}}^{2}}}{\hbar}\int{\frac{1}{{d\left(\tau\right)}}d\tau},$
(5)

where $m_{0}$ and $\tau$ are, respectively, the static mass and the proper
time of the particles. $d\left(\tau\right)$ represents the distance between
two particles and is mainly determined by the geodesic deviation. Eq. (5)
shows that the key to calculate phase change is to obtain the geodesic
deviation distance $d\left(\tau\right)$ of the particle pair.

To do this, we use the tetrad formalism and introduce the following orthogonal
normalized tetrad of flat FLRW universe

$\displaystyle{e_{(0)}}^{\mu}=(\partial t)^{\mu},\ \ \ \ \ \ \ \ \ \ \ \ \
{e_{(1)}}^{\mu}=\frac{1}{{a\left(t\right)}}(\partial r)^{\mu},$ (6)
$\displaystyle{e_{(2)}}^{\mu}=\frac{1}{{a\left(t\right)r}}(\partial\theta)^{\mu},\
\
{e_{(3)}}^{\mu}=\frac{1}{{a\left(t\right)r\sin\theta}}(\partial\varphi)^{\mu},$

so that they form an orthogonal base

$g_{\mu\nu}{e_{(a)}}^{\mu}{e_{(b)}}^{\nu}=\eta_{(a)(b)},$ (7)

where $(a=0,1,2,3)$ and $\eta_{(a)(b)}=\text{diag}.(1,-1,-1,-1)$.

Recalling the geodesic equations in terms of four-velocity $u^{\mu}(\tau)$ is

${u^{\mu}}\circ({e_{(a)}}_{\nu}{u^{\nu}})+({e^{(b)}}_{\nu}{u^{\nu}}){(d{e_{(b)}})_{\rho\sigma}}{u^{\rho}}{e_{(a)}}^{\sigma}=0.$
(8)

After substituting Eq. (6) into Eq. (8), it reduces to

$\begin{split}\ddot{t}\left(\tau\right)+a^{2}(t)H(t){{\dot{r}(\tau)}^{2}}=0,\\\
\ddot{r}\left(\tau\right)+2H(t)\dot{t}\left(\tau\right)\dot{r}\left(\tau\right)=0,\end{split}$
(9)

where the dot represents differentiation with respect to $\tau$. $H(t)$ in Eq.
(9) is defined as $H\equiv\frac{da(t)/dt}{a(t)}$, which is determined by the
Friedmann equation:

$H^{2}=H_{0}^{2}\left(\Omega_{\Lambda}+\Omega_{M}\left(\frac{a_{0}}{a}\right)^{3}+\Omega_{R}\left(\frac{a_{0}}{a}\right)^{4}\right),$
(10)

The full geodesics of the FLRW universe can be obtained by solving the
geodesic equation (9) and the Friedmann equation (10).

Now let us turn to calculate the geodesic deviation vectors, from which one
can read off values of $d(\tau)$, a key term in calculating $\delta\phi$ as
shown in Eq. (5). We obtain these vectors by solving the geodesic deviation
equations in tetrad formalism. Note that unlike the fixed tetrad used in Eq.
(6) for geodesics, the tetrad here must be parallelly transported along the
geodesics r7 , namely, ${{{\tilde{e}}_{(a)}}^{\mu}}{}_{;\nu}v^{\nu}=0$, where
$v^{\nu}$ is the tangent vector of the geodesics. This means that the
orientations of the axes are fixed and they are non-rotating as determined by
local dynamical experiments. It turns out that the following choice of tetrads
are appropriate

$\displaystyle{\tilde{e}_{(0)}}^{\mu}=\dot{t}(\partial
t)^{\mu}+\dot{r}(\partial r)^{\mu},$ (11)
$\displaystyle{\tilde{e}_{(1)}}^{\mu}=a\left(t\right)\dot{r}(\partial
t)^{\mu}+\frac{{\dot{t}}}{{a\left(t\right)}}(\partial r)^{\mu},$
$\displaystyle{\tilde{e}_{(2)}}^{\mu}=\frac{1}{{a\left(t\right)r}}(\partial\theta)^{\mu},$
$\displaystyle{\tilde{e}_{(3)}}^{\mu}=\frac{1}{{a\left(t\right)r\sin\left(\theta\right)}}(\partial\varphi)^{\mu}.$

Note that in order for the above tetrad satisfying the orthogonal normalized
condition (7), one extra condition should be imposed

${\dot{t}^{2}}-a{\left(t\right)^{2}}{\dot{r}^{2}}=1.$ (12)

The geodesic deviation equation in the above frames is given by

$\frac{{{d^{2}}{w^{(a)}}}}{{d{\tau^{2}}}}+{k^{(a)}}_{(b)}{w^{(b)}}=0,$ (13)

where $w^{(a)}$ is geodesic deviation vector in tetrad form

$w^{(a)}={\tilde{e}^{(a)}}_{\mu}w^{\mu},$ (14)

and

${k^{(a)}}_{(b)}=-{R^{\mu}}_{\nu\rho\sigma}{\tilde{e}^{(a)}}_{\mu}{v^{\nu}}{v^{\rho}}{\tilde{e}_{(b)}}^{\sigma}.$
(15)

Substituting the tetrad (11) and (12) into above definition, one finds that
the nonvanishing components of ${k^{(a)}}_{(b)}$ are

$\displaystyle{k^{(1)}}_{(1)}$ $\displaystyle=$
$\displaystyle-\frac{1}{a(t)}\frac{dH}{dt},$ (16)
$\displaystyle{k^{(2)}}_{(2)}$ $\displaystyle=$
$\displaystyle-\frac{1}{a(t)}\frac{dH}{dt}{\dot{t}}^{2}+H^{2}{\dot{r}}^{2}.$
(17)

To simplify our analysis, we assume that $w^{\mu}$ is orthogonal to
${e_{(0)}}^{\mu}$, giving us $w^{(0)}=0$. Since the background exhibits
spherical symmetry, the second and third components of $w^{(a)}$ are simply
symmetric angular components. Thus, we can choose a suitable frame where the
initial value of $w^{(3)}$ is zero. Notably, Eq. (13) is a homogeneous
equation. Consequently, if $\hat{w}^{(a)}$ is a solution, so is
$\kappa\hat{w}^{(a)}$, where $\kappa$ is an arbitrary nonzero constant.
Therefore, if we initially set $w^{(3)}$ and its first time derivative to
zero, it will remain zero throughout. Ultimately, we are left with two non-
zero components: $w^{(1)}$ and $w^{(2)}$.

By substituting these components into the geodesic deviation equation (13) and
performing numerical calculations with specified initial conditions, one can
determine the distance between the pair of particles

$d(\tau)=\sqrt{{w^{(1)}}{{\left(\tau\right)}^{2}}+{w^{(2)}}{{\left(\tau\right)}^{2}}}\Delta
l.$ (18)

The change in proper time is thus given by

$\delta\tau=-{m_{0}}\int_{0}^{\tau^{\prime}}{\frac{1}{{d\left(\tau\right)}}d\tau},$
(19)

which leads to the change in phase
$\delta\phi=-\frac{{{m_{0}}\delta\tau}}{\hbar}$.

In our numerical simulations, the particle’s static mass ${m_{0}}$ is also set
to be ${10^{-25}}\text{kg}$. And we adopt the following data:
$\Omega_{\Lambda}=0.6935$, $\Omega_{m}=0.3065$, $\Omega_{R}=0$ and
$H_{0}=67.6~{}\text{km}\cdot\text{s}^{-1}\cdot\text{Mpc}^{-1}$ as suggested in
Plank 2018  data Aghanim:2018eyx . The scale factor at present ${a_{0}}$ is
set to 1.

One point deserves mention is that solving the geodesic deviation equations of
present case needs assign both the initial temporal coordinate $t_{0}$ and the
initial spatial coordinate $r_{0}$. In the following simulations, we choose
four typical sets of geodesics and their initial values which are adopted in
the following numerical simulations can be found in Table 1. Note that
throughout the paper, $z$ represents the redshift of the source. It is not
difficult to generalize the simulations to the cases where $w_{0}^{(1)}\neq
w_{0}^{(2)}$ and $\dot{w}_{0}^{(1)}\neq\dot{w}_{0}^{(2)}$.

First of all, let’s focus on two assignments. One is that an initial time
slice ($t=t_{0}$) is chosen, while particles start at different radial
coordinates $r=r_{0}$ (that is, a snapshot of different place at time
$t_{0}$). The other is an initial radial coordinate $r_{0}$ is fixed, while
initial time coordinate is free (i.e., focusing on one fixed place for
different time). These two series of geodesics are illustrated in Fig.1.

Table 1: Sets of initial values for numerical simulations

Sets | redshift $z$ | $r_{0}$ | $w_{0}^{(1)}=w_{0}^{(2)}$ | $\dot{w}_{0}^{(1)}=\dot{w}_{0}^{(2)}$
---|---|---|---|---
S1 | 3.94 | $4.5\times{10^{17}}$ | $\frac{1}{3}\times 10^{-18}$ | 0
S2 | 3.94 | $4.5\times{10^{17}}$ | $\frac{2}{3}\times 10^{-18}$ | 0
S3 | 3.94 | $4.5\times{10^{17}}$ | $\frac{1}{3}\times 10^{-18}$ | $10^{-35}$
S4 | 3.94 | $4.5\times{10^{17}}$ | $\frac{1}{3}\times 10^{-18}$ | $-10^{-35}$

  * 1

Note that this is in natural unit system. In SI system, $r_{0}$, $w_{0}^{(i)}$
and $\dot{w}_{0}^{(i)}$ should be multiplied by $c=3\times 10^{8}\text{m/s}$.
For example, for S1, we have $r_{0}=1.35\times{10^{26}}\text{m}$,
$w_{0}^{(1)}=w_{0}^{(2)}=10^{-10}\text{m}$ ($\gg\lambda_{c}\sim 10^{-15}$m,
where $\lambda_{c}$ is the Compton wavelength of the particle), and
$\dot{w}_{0}^{(1)}=\dot{w}_{0}^{(2)}=0$.

Figure 1: Two different assignments of initial values for $r_{0}$ and $t_{0}$.
Solid one represents that particles start from different initial radial
coordinates $r_{0}$ at a fixed initial time slice $t=t_{0}$. Dashed line
represents that particles start from different time $t_{0}$ for a fixed radial
coordinate $r_{0}$. The black point marked ”Now” represents the spacetime
point where we receive particles.

(a)

(b)

Figure 2: Phase change $\delta\phi$ as a function of $t_{0}$ and $r_{0}$. (a)
$\delta\phi$ against $r_{0}$ with ${t_{0}}$ fixed at $z=3.94$. (b)
$\delta\phi$ against $z$ (or $t_{0}$) with ${r_{0}}$ fixed at
($r_{0}=4.5\times{10^{17}}$). In both cases, we take
${w^{(1)}}={w^{(2)}}=\frac{1}{3}\times{10^{-18}}$( $\approx 10^{-10}$m),
${{\dot{w}}_{0}^{(1)}}={{\dot{w}}_{0}^{(2)}}=0$ (in natural units system).

Figure 2 presents a plot illustrating the phase change of two series of
geodesics as a function of $t_{0}$ (with fixed $r_{0}$) or $r_{0}$ (with fixed
$t_{0}$). In Fig.2a, we observe that for geodesics S1 and S2, the phase change
initially increases to a maximum value and then gradually decreases. On the
other hand, for S3 and S4, the phase change diminishes as the radial
coordinate increases, accompanied by a decrease in descent speed. Fig.2b
reveals that the phase change initially grows with increasing initial time
coordinate, followed by a subsequent decrease. This behavior arises due to the
interplay between increased proper time and increased separation during the
movement as the initial time coordinate increases. Furthermore, we can
conclude that specific initial conditions, such as a smaller radial coordinate
and an appropriate time coordinate, yield the maximum phase change.

According to Christodoulou:2018cmk , strict ranges must be set for particle
mass, spacing distance, and action time. However, in our scenario, the
particle’s static mass can be significantly smaller, and the spacing distance
can be much larger. For instance, when the initial time coordinate is
$1.5\times 10^{17}$ and the radial coordinate is $4.5\times 10^{17}$, the
phase change of particle pairs with a static mass of $1\times 10^{-20}$ kg,
separated by one meter initially, and with no initial deviation rate, can
reach a substantial value of 0.361988. Such an experiment would be impractical
to conduct on Earth due to the extremely light particles and the large spacing
distance required. The relationship between phase change and initial geodesic
offset distance under these conditions is depicted in Fig.3.

Figure 3: Phase change $\delta\phi$ as a function of initial geodesic
deviation $w_{0}$. We also take
${t_{0}}=1.5\times{10^{17}},{r_{0}}=4.5\times{10^{17}},{w_{0}^{(1)}}={w_{0}^{(2)}}={w_{0}},{\dot{w}_{0}^{(1)}}={\dot{w}_{0}^{(2)}}=0$
(in natural units system).

We now investigate the influence of the initial velocity of the geodesic
deviation vector on the phase change. By keeping the starting point of the
particles fixed at a spacetime point and varying the initial velocity, we can
analyze this relationship. The specific initial conditions are provided in the
caption of the Fig.4.

In Fig.4a, we observe a concave phase change curve that indicates a decreasing
rate of phase change as the positive initial deviation velocity increases.
This behavior can be easily understood: the positive initial deviation
velocity causes an increase in the spacelike distance, as depicted in Eq.
(19), leading to a smaller phase change.

On the other hand, Fig.4b demonstrates that the phase change initially
increases with a negative initial deviation velocity and then gradually
decreases until it reaches nearly zero. This can be explained by considering
the following: with a relatively small negative velocity, the two particles
remain close neighbors for most of the time. However, with a relatively large
negative velocity, they move apart, resulting in a considerable separation
between them by the end.

(a)

(b)

Figure 4: Phase change $\delta\phi$ as a function of initial geodesic
deviation velocity $\dot{w}_{0}$. (a) Initial geodesic deviation velocity is
positive, ${\dot{w}_{0}^{(1)}}={\dot{w}_{0}^{(2)}}>0$. (b) Initial geodesic
deviation velocity is negative, ${\dot{w}_{0}^{(1)}}={\dot{w}_{0}^{(2)}}<0$.
In both cases, we take
${t_{0}}=5\times{10^{16}},{r_{0}}=4.5\times{10^{17}},{w_{0}^{(1)}}={w_{0}^{(2)}}=\frac{1}{3}\times{10^{-18}}$
(in natural units system).

##  III Characteristic spectrum

Entanglement between particle pairs might have been generated before geodesic
motion starts. Therefore, it becomes crucial to determine whether gravity
induces entanglement or if other physical processes are involved. To address
this challenge, we propose that QGEM during geodesic motion will exhibit a
characteristic spectrum as phase changes occur in a series of geodesic lines.
By analyzing the entangled patterns formed by different geodesics, we can
infer whether the entanglement originates from the gravitational field of
nearby particles or from alternative sources.

Phase itself is not directly observable. We often use entanglement witness
$\mathcal{W}$ as an experimental indicator to detect entanglement formation.
The definition of entanglement witness is

$\mathcal{W}=\left|{\left\langle{\sigma{1_{x}}\otimes\sigma{2_{z}}}\right\rangle+\left\langle{\sigma{1_{y}}\otimes\sigma{2_{y}}}\right\rangle}\right|.$
(20)

When it’s greater than 1, we can infer that there is entanglement between the
two particles.

Similar to the case of Zhang:2023llk , the entanglement witness $\mathcal{W}$
exhibits oscillatory variations with different initial motion parameters.
Again, let’s take two series of geodesics as mentioned in Fig. 1 for exampled.
In Fig. 5, we have converted different initial parameters $r_{0}$ into
corresponding initial velocity $u^{\mu}(0)$. While the initial values of four
velocity, $u^{\mu}(0)$, can be further transformed to the orthogonal frame
form as:

$u^{\mu}(\tau=0)=\gamma{e_{(0)}}^{\mu}-\gamma v_{0}{e_{(1)}}^{\mu},$ (21)

where $\gamma=\frac{1}{{\sqrt{1-{v_{0}^{2}}}}}$ and $v_{0}$ is the initial
velocity of the particle with which we replace $r_{0}$. It shows $v_{0}$
increases with $r_{0}$ increases and $\mathcal{W}$ oscillates faster in the
larger area of $v_{0}$.

Figure 5: Variation of entanglement witness $\mathcal{W}$ as a function of
$v_{0}$ in FLRW universe. Other initial parameters are set to the same as S1
in Table1

$\mathcal{W}$ oscillates very unevenly as a function of $r_{0}$. Only when
$r_{0}$ falls in some certain intervals, the entanglement phase that gives a
positive answer to whether to be entangled can be reached. This forces us to
limit the particle source to a reasonable spacetime region in order to ensure
that entangled particles can be observed when using the entanglement witness
(20).

Fig. 6 shows the entanglement witness as a function of the redshift $z$. From
where we find that due to the convexity of the phase transition curve S1 in
Fig. 2b, the entanglement witness curve has a valley in the middle of z that
is insufficient to confirm entanglement. At both ends of z, $\mathcal{W}$
oscillates fast enough to reach the threshold, 1, for confirming entanglement.

Figure 6: Variation of entanglement witness $\mathcal{W}$ as a function of $z$
in FLRW universe. Other initial parameters are set to the same as S1 in Table1

Compared to the initial conditions, particle’s energy is an important
observable. In observer’s frame, each particle emitted at different $r_{0}$,
$t_{0}$ will have a specific geodesic, and every instantaneous observer on the
geodesic will measure a corresponding kinetic energy.

We use the kinetic energy of the particles measured by observer on Earth to
label each particle. The observed kinetic energy for particles emitted from
different initial spacetime location is calculated as follows:

$\displaystyle{E_{kin}}={z_{a}}{p^{a}}-m,$ (22)

where ${z_{a}}$ is the four-velocity of an observer on Earth, ${p^{a}}$ is the
four-velocity of a particle upon reaching Earth. The variation of
$\mathcal{W}$ with respect to $E_{kin}$ is plotted in Fig. 7 and Fig. 8. From
these two figures, we observe that for both cases, $\mathcal{W}$ oscillates
faster in the low kinetic energy ${E_{kin}}$ region than in the high kinetic
energy region. Observing low-energy particle pairs will be more helpful for us
to detect entanglement.

Figure 7: Variation of entanglement witness $\mathcal{W}$ as a function of
$E_{kin}$ with $t_{0}$ fixed at $z=3.94$ in FLRW universe. The other kinematic
initial conditions are consistent with that in Fig. 5, and each ${v_{0}}$
corresponds to each ${E_{kin}}$ of this figure. Figure 8: Variation of
entanglement witness $\mathcal{W}$ as a function of $E_{kin}$ with fixed
${r_{0}}$ at $4.5\times{10^{17}}$ in FLRW universe. The kinematic initial
conditions are consistent with that in Fig. 6, and each $z$ corresponds to
each abscissa ${E_{kin}}$ of this figure.

##  IV Conclusions

In this paper, we generalize the QGEM experiment to be carried out in FLRW
cosmological spacetimes so as to verify quantum gravity after in Schwarzschild
spacetimes in Zhang:2023llk . Just like in Schwarzschild spacetime, particle
pairs propagating in FLRW cosmological spacetime can become entangled as a
result of the gravitational field’s influence, thus verifying quantum gravity.
The entanglement phase is also influenced by variations in its motion
parameters, such as the initial launch point position, initial geodesic
deviation and its rate.

Since the entangled phase cannot be directly observed, in order to connect it
more closely with actual cosmological observations, we drew the entanglement
witness-energy diagram. Similar to Zhang:2023llk , there is a characteristic
spectrum of QGEM. So we could distinguish whether the observed entanglement
arises from the quantum gravity effect of particle pairs or from other process
(e.g, the emission stage at the source as suggested by the Hawking radiation
of black holes).

In addition to verifying the quantum gravity effect, our solution can also
verify the extended equivalence principle in the spacetime background of FLRW
cosmology Giacomini:2020ahk . It also provides a new perspective for the
observation of cosmic rays, enriching the content of entanglement detection in
cosmology. We provide a new possible source of entanglement formation in
cosmology.

In the future, we would like to further study the impact of different
cosmological models on the quantum entanglement induced by the QGEM mechanism.
The entanglement induced in this way under different cosmological models is
very likely to have significantly different characteristics, which will
provide a new way for us to observe and study the universe. Moreover, in
addition to inducing entanglement between physical particles, quantum gravity
can also induce entanglement between quantum fields, such as scalar fields,
vector fields, and spinor fields, thus affecting the overall entangled
structure of the universe. This could be explored further in the future.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China
with the Grants Nos. 12375049, 11975116, and Key Program of the Natural
Science Foundation of Jiangxi Province under Grant No. 20232ACB201008.

## References

  * (1) R. Feynman, “The role of gravitation in physics”, in Chapel Hill Conference Proceedings, edited by C. M. DeWitt and D. Rickles (Edition Open Access, 1957) pp. 250–256.
  * (2) S. Bose, A. Mazumdar, G. W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A. Geraci, P. Barker, M. S. Kim and G. Milburn, “Spin Entanglement Witness for Quantum Gravity,” Phys. Rev. Lett. 119, no.24, 240401 (2017) [arXiv:1707.06050 [quant-ph]].
  * (3) C. Marletto and V. Vedral, “Gravitationally induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity,” Phys. Rev. Lett. 119, no.24, 240402 (2017) [arXiv:1707.06036 [quant-ph]].
  * (4) M. Christodoulou and C. Rovelli, “On the possibility of laboratory evidence for quantum superposition of geometries,” Phys. Lett. B 792, 64-68 (2019) [arXiv:1808.05842 [gr-qc]].
  * (5) A. Capolupo, G. Lambiase, A. Quaranta and S. M. Giampaolo, “Probing axion mediated fermion–fermion interaction by means of entanglement,” Phys. Lett. B 804, 135407 (2020) [arXiv:1910.01533 [hep-ph]].
  * (6) A. Belenchia, R. M. Wald, F. Giacomini, E. Castro-Ruiz, Č. Brukner and M. Aspelmeyer, “Quantum Superposition of Massive Objects and the Quantization of Gravity,” Phys. Rev. D 98, no.12, 126009 (2018) [arXiv:1807.07015 [quant-ph]].
  * (7) D. Carney, P. C. E. Stamp and J. M. Taylor, “Tabletop experiments for quantum gravity: a user’s manual,” Class. Quant. Grav. 36, no.3, 034001 (2019) [arXiv:1807.11494 [quant-ph]].
  * (8) R. J. Marshman, A. Mazumdar and S. Bose, “Locality and entanglement in table-top testing of the quantum nature of linearized gravity,” Phys. Rev. A 101, no.5, 052110 (2020) [arXiv:1907.01568 [quant-ph]].
  * (9) T. Westphal, H. Hepach, J. Pfaff and M. Aspelmeyer, “Measurement of gravitational coupling between millimetre-sized masses,” Nature 591, no.7849, 225-228 (2021) [arXiv:2009.09546 [gr-qc]].
  * (10) M. Carlesso, A. Bassi, M. Paternostro and H. Ulbricht, “Testing the gravitational field generated by a quantum superposition,” New J. Phys. 21, no.9, 093052 (2019) [arXiv:1906.04513 [quant-ph]].
  * (11) D. L. Danielson, G. Satishchandran and R. M. Wald, “Gravitationally mediated entanglement: Newtonian field versus gravitons,” Phys. Rev. D 105, no.8, 086001 (2022) [arXiv:2112.10798 [quant-ph]].
  * (12) M. Christodoulou, A. Di Biagio, M. Aspelmeyer, Č. Brukner, C. Rovelli and R. Howl, “Locally Mediated Entanglement in Linearized Quantum Gravity,” Phys. Rev. Lett. 130, no.10, 100202 (2023) [arXiv:2202.03368 [quant-ph]].
  * (13) H. T. Cho and B. L. Hu, “Quantum noise of gravitons and stochastic force on geodesic separation,” Phys. Rev. D 105, no.8, 086004 (2022) [arXiv:2112.08174 [gr-qc]].
  * (14) A. Matsumura and K. Yamamoto, “Gravity-induced entanglement in optomechanical systems,” Phys. Rev. D 102, no.10, 106021 (2020) [arXiv:2010.05161 [gr-qc]].
  * (15) D. Miki, A. Matsumura and K. Yamamoto, “Entanglement and decoherence of massive particles due to gravity,” Phys. Rev. D 103, no.2, 026017 (2021) [arXiv:2010.05159 [gr-qc]].
  * (16) R. Howl, N. Cooper and L. Hackermüller, “Gravitationally-induced entanglement in cold atoms,” [arXiv:2304.00734 [quant-ph]].
  * (17) J. Yant and M. Blencowe, “Gravitationally induced entanglement in a harmonic trap,” Phys. Rev. D 107, no.10, 106018 (2023) [arXiv:2302.05463 [quant-ph]].
  * (18) S. Feng, B. M. Gu and F. W. Shu, “Detecting Extra Dimension By the Experiment of the Quantum Gravity Induced Entanglement of Masses,” [arXiv:2307.11391 [gr-qc]].
  * (19) M. Schut, A. Grinin, A. Dana, S. Bose, A. Geraci and A. Mazumdar, “Relaxation of experimental parameters in a Quantum-Gravity Induced Entanglement of Masses Protocol using electromagnetic screening,” [arXiv:2307.07536 [quant-ph]].
  * (20) P. Fragolino, M. Schut, M. Toroš, S. Bose and A. Mazumdar, “Decoherence of a matter-wave interferometer due to dipole-dipole interactions,” [arXiv:2307.07001 [quant-ph]].
  * (21) P. Li, Y. Ling and Z. Yu, “Generation rate of quantum gravity induced entanglement with multiple massive particles,” Phys. Rev. D 107, no.6, 064054 (2023) [arXiv:2210.17259 [gr-qc]].
  * (22) S. Bose, A. Mazumdar, M. Schut and M. Toroš, “Mechanism for the quantum natured gravitons to entangle masses,” Phys. Rev. D 105, no.10, 106028 (2022) [arXiv:2201.03583 [gr-qc]].
  * (23) F. He and B. Zhang, “Generation of entanglement between two laser pulses through gravitational interaction,” Eur. Phys. J. Plus 138, no.2, 141 (2023) [arXiv:2302.06362 [gr-qc]].
  * (24) C. Zhang and F. W. Shu, “Gravity-induced entanglement between two massive microscopic particles in curved spacetime: I.The Schwarzschild background,” [arXiv:2308.16526 [gr-qc]].
  * (25) E. Martin-Martinez and N. C. Menicucci, “Cosmological quantum entanglement,” Class. Quant. Grav. 29, 224003 (2012) [arXiv:1204.4918 [gr-qc]].
  * (26) S. Brahma, O. Alaryani and R. Brandenberger, “Entanglement entropy of cosmological perturbations,” Phys. Rev. D 102, no.4, 043529 (2020) [arXiv:2005.09688 [hep-th]].
  * (27) C. Kiefer and D. Polarski, “Why do cosmological perturbations look classical to us?,” Adv. Sci. Lett. 2, 164-173 (2009) [arXiv:0810.0087 [astro-ph]].
  * (28) J. Martin, “Inflationary perturbations: The Cosmological Schwinger effect,” Lect. Notes Phys. 738, 193-241 (2008) [arXiv:0704.3540 [hep-th]].
  * (29) J. Martin and V. Vennin, “Quantum Discord of Cosmic Inflation: Can we Show that CMB Anisotropies are of Quantum-Mechanical Origin?,” Phys. Rev. D 93, no.2, 023505 (2016) [arXiv:1510.04038 [astro-ph.CO]].
  * (30) J. Martin and V. Vennin, “Bell inequalities for continuous-variable systems in generic squeezed states,” Phys. Rev. A 93, no.6, 062117 (2016) [arXiv:1605.02944 [quant-ph]].
  * (31) L. Lello, D. Boyanovsky and R. Holman, “Entanglement entropy in particle decay,” JHEP 11, 116 (2013) [arXiv:1304.6110 [hep-th]].
  * (32) A. Valentini, “Astrophysical and cosmological tests of quantum theory,” J. Phys. A 40, 3285-3303 (2007) [arXiv:hep-th/0610032 [hep-th]].
  * (33) J. W. Chen, S. H. Dai, D. Maity, S. Sun and Y. L. Zhang, “Towards Searching for Entangled Photons in the CMB Sky,” Phys. Rev. D 99, no.2, 023507 (2019) [arXiv:1701.03437 [quant-ph]].
  * (34) R. D. Inverno, Introducing Einstein’s Relativity (Clarendon Press, Oxford, 1992).
  * (35) N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (36) F. Giacomini and Č. Brukner, “Einstein’s Equivalence principle for superpositions of gravitational fields and quantum reference frames,” [arXiv:2012.13754 [quant-ph]].