# Classical Sampling of Random Quantum Circuits with Bounded Fidelity

Gleb Kalachev<EMAIL_ADDRESS>Pavel Panteleev
<EMAIL_ADDRESS>Huawei 2012 Lab  Lomonosov Moscow State University.
PengFei Zhou<EMAIL_ADDRESS>Huawei 2012 Lab  Man-Hong Yung
<EMAIL_ADDRESS>Huawei 2012 Lab  Shenzhen Institute for Quantum Science
and Engineering, Southern University of Science and Technology, Shenzhen,
Guangdong 518055, China  International Quantum Academy, and Shenzhen Branch,
Hefei National Laboratory, Futian District, Shenzhen, Guangdong 518048, China
Guangdong Provincial Key Laboratory of Quantum Science and Engineering,
Southern University of Science and Technology, Shenzhen, Guangdong 518055,
China  Shenzhen Key Laboratory of Quantum Science and Engineering, Southern
University of Science and Technology, Shenzhen, Guangdong 518055, China

###### Abstract

Random circuit sampling has become a popular means for demonstrating the
superiority of quantum computers over classical supercomputers. While quantum
chips are evolving rapidly, classical sampling algorithms are also getting
better and better. The major challenge is to generate bitstrings exhibiting an
XEB fidelity above that of the quantum chips. Here we present a classical
sampling algorithm for producing the probability distribution of any given
random quantum circuit, where the fidelity can be _rigorously_ bounded.
Specifically, our algorithm performs rejection sampling after the introduced
very recently multi-tensor contraction algorithm. We show that the fidelity
can be controlled by partially contracting the dominant paths in the tensor
network and by adjusting the number of batches used in the rejection sampling.
As a demonstration, we classically produced 1 million samples with the
fidelity bounded by 0.2%, based on the 20-cycle circuit of the Sycamore
53-qubit quantum chip. Though this task was initially estimated to take 10,000
years on the Summit supercomputer, it took about 14.5 days using our algorithm
on a relatively small cluster with 32 GPUs (Tesla V100 16GB). Furthermore, we
estimate that for the Zuchongzhi 56-qubit 20-cycle circuit one can produce 1M
samples with fidelity 0.066% using the Selene supercomputer with 4480 GPUs
(Tesla A100 80GB) in about 4 days.

quantum simulation, quantum supremacy, tensor network

††preprint: APS/123-QED

One of the main motivations for quantum computing is a belief, shared by many
researches, that quantum systems are very hard to simulate classically. In the
language of complexity theory, it is usually expressed as a plausible but
currently unproven conjecture that $\mathbf{BPP}\subsetneq\mathbf{BQP}$, where
the complexity class $\mathbf{BPP}$ is, informally speaking, the class of
problems efficiently solvable by classical computers, and $\mathbf{BQP}$ is
the corresponding class for quantum computers. If this conjecture were true,
it would imply, at least theoretically, that quantum computers have a
_computational advantage_ (or _quantum supremacy_) over the classical ones in
the asymptotic regime, as the problem size grows to infinity. In fact, to
demonstrate such a quantum advantage on a real quantum hardware one should
also deal with the decoherence problem and other imperfections. While the
general consensus is that, under some reasonable assumptions on the noise
scaling, this problem can be solved by applying a fault-tolerant protocol [1],
there are some researchers who still doubt these assumptions [2].

From the practical point of view, a convincing argument in favor of the
quantum advantage would be an experiment demonstrating that some well-defined
problem can be efficiently solved on a real quantum computer, while the
solution of the same problem using the state-of-the-art algorithms cannot be
obtained within a reasonable amount of time even with the best classical
supercomputers in the world [3, 4, 5]. At the current stage, a leading
candidate for such a quantum advantage experiment, which was already performed
on a real quantum hardware [6, 7, 8], is the task of sampling from the output
probability distributions of _random quantum circuits_ (RQCs). Though
initially this task was estimated [6] to take thousands of years on the
fastest classical supercomputers, later the simulation time was significantly
reduced due to the recent progress on tensor network (TN) based quantum
simulation algorithms [9, 10, 11, 12, 13, 14, 15]. Another interesting idea,
which was used recently to significantly reduce the computational cost of RQCs
simulation, is to compress the wave function using matrix product states [16].

Usually TN based simulators for RQCs calculate single amplitudes by
contracting the tensor network representing the circuit. In fact, it was shown
in [17, 18] that TN contraction can also be used to produce not only single
amplitudes but also large batches of amplitudes (i.e., the amplitudes for a
collection of bitstrings that share some fixed bits), and the computational
cost of finding a batch of amplitudes is usually similar to the cost of
finding a single amplitude. This idea has been applied recently [13] to spoof
the  _linear cross-entropy benchmarking_ (linear XEB) test, which was used in
Google’s quantum advantage experiment [6] as a way to certify the fidelity.
Explicitly, the linear XEB fidelity $\mathcal{F}_{\mathrm{XEB}}$ for a
sequence of bitstrings $s_{1},\ldots,s_{k}$, produced in an experiment is
defined as

$\mathcal{F}_{\mathrm{XEB}}\equiv\frac{2^{n}}{k}\sum_{i=1}^{k}p_{C}(s_{i})-1\
,$ (1)

where $p_{C}(\cdot)$ is the theoretical output probability distribution for
the circuit $C$ used in the experiment. To spoof the linear XEB test and
produce $N$ samples with $\mathcal{F}_{\mathrm{XEB}}\geq f$ one can choose
$fN$ samples $s$ with the highest probabilities $p_{C}(s)$ from a large batch
of amplitudes found with 100% fidelity and then add $(1-f)N$ uniformly random
bitstrings. At first, it may seem that this very simple way of spoofing is
quite easy to detect since the produced samples are highly correlated (in $f$
fraction of samples the values for some bit positions are fixed). However, if
one uses a sufficiently large number of batches, the spoofing algorithm is
much harder to distinguish from the real simulation since the bit positions
are no longer fixed.

In fact, one can use a large number of batches not only for spoofing but also
for the _sampling_ task, where we want to produce $N$ independent random
samples according to the probability distribution $p_{C}(\cdot)$. A rather
standard approach here is to apply the _frugal rejection sampling_ algorithm
[9] to generate random bitstrings by calculating their probabilities and
accepting each bitstring with probability proportional to the calculated
value. This approach works for random quantum circuits where all bitstrings
have probabilities of the same order of magnitude. Unfortunately, for a large
number of samples the computational cost of this approach is quite high if the
batches are calculated independently one by one. At the same time, it was
shown recently [19] that by applying the _multi-tensor contraction_ algorithm,
which reuses the partial contraction results, one can reduce the computational
cost in this case by several orders of magnitude. In fact, it is shown in [19]
that not only the sampling task but also the much harder verification task for
RQCs, where one needs to find the _exact_ amplitudes for a large collection of
_uncorrelated_ bitstrings, can also be solved in several days on a modern
supercomputer though it was initially estimated in [6] to take millions of
years.

In the current work, we propose a slight modification of the frugal rejection
sampling from [9], which uses more batches than the actual number of the
random samples we need to produce. This allows us to give a rigorous analysis
of the proposed algorithm in terms of the fidelity and the statistical
variation distance to the ideal probability distribution. To produce $m$
samples we calculate $\alpha m$ random small batches (e.g., of size 64), and
then use our variant of the frugal rejection sampling algorithm to produce $m$
samples out of these $\alpha m$ batches. Our analysis indicates that if
$\alpha=2$, then the statistical variation distance between the probability
distributions of our algorithm and the ideal sampling is negligible, and
therefore in all our simulations we assume that $\alpha=2$. Note that the
computational cost grows less than linearly as $\alpha\to\infty$, and for
$\alpha=2$ the simulation time is less then two times larger than for
$\alpha=1$.

Moreover, to further reduce the computational cost of the simulations we also
apply the _partial slicing_ summations [9] in the tensor-network contraction.
This allows us to simulate RQCs with a target fidelity $f$. The partial
slicing is similar to the gate decomposition used in the Shrödinger-Feynman
algorithm from [6], where some 2-qubit gates are decomposed into sums of pairs
of $1$-qubit gates. If the target fidelity $f$ in a simulation is less than
$1$, one can speed up by skipping some of the terms in the summation [9, 6].
In our case, if we have $k$ sliced vertices in the tensor network, then we can
sum only over $f2^{k}$ out of $2^{k}$ slices and obtain the result with the
fidelity approximately equal to $f$. Compared with the method in Refs [9, 6],
the key feature of the current work is that our method can accurately predict
the obtained fidelity for a given set of slices and choose slices in order to
maximize the fidelity. Moreover, we will show that this fidelity can be found
by a contraction of some specifically designed tensor network, which gives a
new general way to control the fidelity when we apply the partial slicing
method.

In the current work, we demonstrate our experimental results using these new
algorithms, where we:

  1. 1.

produce samples for Google’s supremacy (ABCD) circuits from [6] up to $20$
cycles;

  2. 2.

show (using an approach similar to [13]) how to spoof the linear XEB test for
the hardest case in the Google’s experiment with fidelity $0.2\%$ in just $4$
hours on a personal computer with only one GPU.

All the experimental data produced in these experiments can be found in [20].

Update. Recently, after all the experiments in the current paper were already
finished we became aware of the work [15], where an approach, very similar to
the multi-tensor contraction algorithm from [19], was used in combination with
other techniques to significantly reduce the simulation time. Though the
computational cost of our algorithm is similar to the one from [15], the
former comes with a rigorous analysis of the fidelity, while the latter is
only justified by empirical estimates. Moreover, in the current paper, we also
independently confirm our analytical estimates of the fidelity using the
verification algorithm from [19].

##  I Simulation with target fidelity

In this section, we present a general method for simulating RQCs with given
target fidelity $f$. In general, for a quantum circuit $C$, by a standard
procedure proposed in [6], to produce multiple independent random samples, one
needs to calculate multiple independent batches of amplitudes, which is
usually a difficult computing task. However, the multi-tensor simulator from
[19] provides a much more efficient way to accomplish this task by utilizing a
global cache which can reuse some intermediate tensors to significantly save
the computing time. The algorithm from [19] also uses the simulated annealing
method to optimize the contraction tree and the list of sliced variables. Let
us remind that a  _contraction tree_ [21, 22, 10] encodes a particular way we
perform the contraction for a given tensor network. At the same time, the
_sliced variables_ (also called the _projected variables_) correspond to the
variables (i.e., the tensor legs) that we sum over at the very last step,
which allows us to control the memory budget of the contraction.

Note that if we do the full summation for all the sliced vertices, the target
fidelity of the obtained amplitudes is equal to $1$. However, in many cases we
cannot afford this, and our aim is to find the amplitudes with some target
fidelity $f\ll 1$. In this case, we divide the sliced variables into two
parts: the _partially sliced_ variables where we perform a partial summation
(i.e., sum only over $f$ fraction of slices), and the remaining _fully sliced_
variables, where perform the full summation. The former aim to find a balance
on the sampling time complexity and the fidelity of the produced samples,
while the latter only aim to reduce the intermediate memory of the tensor
contraction. Next, we will show how to choose the partially sliced vertices
and the particular $f$ fraction of slices for them.

###  I.1 Achieving target fidelity by choosing slices of maximal norms

Let us first show how to choose the $f$ fraction of slices we sum over on the
chosen partially sliced vertices. Let $C=C_{2}C_{1}$ be a quantum circuit
divided into the sub-circuits $C_{1}$ and $C_{2}$. Moreover, we assume that
the partially sliced vertices are already chosen in the cut of $C$ splitting
it into parts $C_{1}$ and $C_{2}$ (see Fig. 1). For the vector
$\left|\psi\right\rangle\equiv C\left|0\right\rangle$, we have

$\left|\psi\right\rangle=C_{2}C_{1}\left|0\right\rangle=\sum_{i\in\\{0,1\\}^{k}}\left|\psi_{i}\right\rangle=\sum_{i\in\\{0,1\\}^{k}}\sum_{j\in\\{0,1\\}^{n-k}}\left|\psi_{i}^{j}\right\rangle,$

where $n$ is number of qubits, $k$ is the number of the partially sliced
vertices in the cut of $C$. Here and below in this section we will assume that
$i,i^{\prime},i^{\prime\prime}\in\\{0,1\\}^{k}$, and
$j,j^{\prime}\in\\{0,1\\}^{n-k}$. The quantum state $\left|\psi\right\rangle$
can be expressed as:

$\left|\psi\right\rangle=C_{2}C_{1}\left|0\right\rangle=\sum_{i\in\\{0,1\\}^{k}}\underbrace{\sum_{j\in\\{0,1\\}^{n-k}}\overbrace{C_{2}\left|ji\right\rangle\left\langle
ji\right|C_{1}\left|0\right\rangle}^{\left|\psi_{i}^{j}\right\rangle}}_{\left|\psi_{i}\right\rangle}.$

Note that when $i\neq i^{\prime}$, the vectors $\left|\psi_{i}\right\rangle$
and $\left|\psi_{i^{\prime}}\right\rangle$ are orthogonal:

$\left\langle\psi_{i}\middle|\psi_{i^{\prime}}\right\rangle=\sum_{j,j^{\prime}}\left\langle
0\right|C_{1}^{*}\left|ji\right\rangle\underbrace{\left\langle
ji\right|\overbrace{C_{2}^{*}C_{2}}^{I}\left|j^{\prime}i^{\prime}\right\rangle}_{=0\mbox{
if }i\neq i^{\prime}}\left\langle
j^{\prime}i^{\prime}\right|C_{1}\left|0\right\rangle=0.$

For every set $X\subset\\{0,1\\}^{k}$ we define the vector
$\left|\psi_{X}\right\rangle=\sum_{i\in X}\left|\psi_{i}\right\rangle$. It is
the orthogonal projection of the vector $\left|\psi\right\rangle$ on the
vector
$\left|\overline{\psi_{X}}\right\rangle=\left|\psi_{X}\right\rangle/\|\left|\psi_{X}\right\rangle\|$.
The fidelity between 2 mixed states defined by density matrices $\rho$ and
$\sigma$ is defined by formula:

$F(\rho,\sigma)=\left(\mathop{\mathrm{tr}}\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right)^{2}.$

For the pure states
$\rho=\left|\psi_{\rho}\right\rangle\left\langle\psi_{\rho}\right|$ and
$\sigma=\left|\psi_{\sigma}\right\rangle\left\langle\psi_{\sigma}\right|$ this
formula can be simplified:

$F(\rho,\sigma)=\left|\left\langle\psi_{\rho}\middle|\psi_{\sigma}\right\rangle\right|^{2}.$

The vector $\left|\overline{\psi_{X}}\right\rangle$ has unit norm, hence it
can be interpreted as a pure state. Then we can calculate the fidelity between
the states $\left|\overline{\psi_{X}}\right\rangle$ and
$\left|\psi\right\rangle$ as

$\mathcal{F}=|\left\langle\overline{\psi_{X}}\middle|\psi\right\rangle|^{2}=\|\left|\psi_{X}\right\rangle\|^{2}=\sum_{i\in
X}\|\left|\psi_{i}\right\rangle\|^{2}.$

Assume we know all the norms $\|\left|\psi_{i}\right\rangle\|$, then we can
define $X$ to be the set of the indices $i\in\\{0,1\\}^{k}$ of the vectors
$\left|\psi_{i}\right\rangle$ with the maximal norms. In this case, we get

$\mathcal{F}=\sum_{i\in
X}\|\left|\psi_{i}\right\rangle\|^{2}\geq\frac{|X|}{2^{k}}\underbrace{\sum_{i}\|\left|\psi_{i}\right\rangle\|^{2}}_{=1}=\frac{|X|}{2^{k}}.$
(2)

Since the vectors $\left|\psi_{i}^{j}\right\rangle$ are also orthogonal, we
have

$\displaystyle\|\left|\psi_{i}\right\rangle\|^{2}=\sum_{j}\|\left|\psi_{i}^{j}\right\rangle\|^{2}$
$\displaystyle=\sum_{j}\underbrace{\|C_{2}\left|ji\right\rangle\|^{2}}_{=1}|\left\langle
ji\right|C_{1}\left|0\right\rangle|^{2}$ $\displaystyle=\sum_{j}|\left\langle
ji\right|C_{1}\left|0\right\rangle|^{2}.$

Let us mention that $\|\left|\psi_{i}\right\rangle\|^{2}$ can be interpreted
as the probability to obtain the state $\left|i\right\rangle$ after we measure
the qubits corresponding to the partially sliced vertices on the output of the
subcirciut $C_{1}$. Moreover, the norms of all vectors
$\left|\psi_{i}\right\rangle$ can be calculated simultaneously as the result
of a contraction for the tensor network shown on the right part of Fig. 1. It
is not hard to see that the result of the contraction for this tensor network
corresponds to

$\|\left|\psi_{i}\right\rangle\|^{2}=\sum_{i^{\prime},i^{\prime\prime}}\sum_{j}\delta_{i,i^{\prime},i^{\prime\prime}}\left\langle
0\right|C_{1}^{*}\left|ji^{\prime}\right\rangle\left\langle
ji^{\prime\prime}\right|C_{1}\left|0\right\rangle,$

where $\delta_{i,j,k}=\begin{cases}1,&\mbox{if }i=j=k,\\\
0,&\mbox{else}\end{cases}$. Recall that the contraction with the tensor
$\delta_{i,j,k}$ is equivalent to the identification of the variables $i$, $j$
and $k$.

circuit
$C$$C_{1}$$C_{2}$$=$$=$$=$$=$$=$$=$$=$$=$$=$$C_{1}$$C_{1}^{*}$$i$$j$$i^{\prime}$$i^{\prime\prime}$—
sliced vertices Figure 1: The tensor network for the calculation of
$\|\left|\psi_{i}\right\rangle\|$. The red dots represent the sliced vertices.
Combined with the blue vertices, they form a cut separating the full circuit.

Here, we point out that in recent paper [15], the authors propose cutting the
tensor network or the circuit into two parts (called the _big head_ and the
_small tail_), and choose 8 slicing vertices that are the input of 4 fSim
gates of the tail part on the cut interface. These sliced vertices on the cut
interface correspond to the partially sliced vertices $i$ in the current
paper. With the belief that all Feynman paths are almost orthogonal to each
other and have equal norms for random quantum circuits [9], in [15] the
authors propose to choose one slice $\left|00000000\right\rangle$ out of the
$2^{8}$ slices and estimate the fidelity as $1/2^{8}$. Note that our method
can be used to find the _exact_ fidelity for the distribution induced from the
slice $\left|00000000\right\rangle$.

###  I.2 Slice selection

In this part, we present an algorithm which shows how to choose the partially
sliced vertices from all sliced vertices and how to do the partial contraction
introduced above in (2).

The input data for the algorithm contains the quantum circuit $C$, the
corresponding contraction tree $T$, the list of sliced vertices $I$, the
target fidelity ${f<1}$. The algorithm produces the set of partially sliced
vertices $S$, the set of slices $X\subset\\{0,1\\}^{|S|}$, and the actual
fidelity $\mathcal{F}$. We proceed as follows.

  1. 1.

Select the set $S\subseteq I$ of $k$ partial summed sliced vertices such that
no vertex from $S$ inside the lightcone of another vertex from $S$ (i.e. all
$k$ vertices are outputs of some subcircuit of $C$). Moreover, vertices from
$S$ should be as close as possible to the inputs of the circuit $C$, and
$\lceil f2^{k}\rceil/(f2^{k})$ should be minimal. The value of $k$ can be
selected in each case depending on the preliminary estimate of the total
algorithm complexity.

  2. 2.

Let $C_{1}$ be the minimal subcircuit, with the same inputs as the inputs of
$C$ and the outputs containing the whole set $S$, i.e., the subcircuit $C_{1}$
consists of the union of the lightcones of all vertices from $S$.

  3. 3.

Calculate the norms $\|\left|\psi_{i}\right\rangle\|$ by contracting the
tensor network shown on Fig. 1.

  4. 4.

Select a set $X\subseteq\\{0,1\\}^{k}$ of minimal size such that $\sum_{i\in
X}\|\left|\psi_{i}\right\rangle\|^{2}\geq f$. It is enough to choose $X$
consisting of $i\in\\{0,1\\}^{k}$ with maximal norms
$\|\left|\psi_{i}\right\rangle\|$. It is easy to see that $|X|\leq\lceil
f2^{k}\rceil$. Together with the set $X$ we obtain the fidelity
$\mathcal{F}=\sum_{i\in X}\|\left|\psi_{i}\right\rangle\|^{2}$.

The computational cost of the simulation with fidelity $f$ can be estimated as

$C_{s}(f)=\frac{|X|}{2^{k}}C_{s}(1)\leq\frac{\lceil
2^{k}f\rceil}{2^{k}}C_{s}(1)<(f+2^{-k})C_{s}(1),$

where $C_{s}(1)$ is the corresponding computational cost with 100% fidelity.

In Appendix B, you can find a more detailed variant of the above slice
selection algorithm.

###  I.3 Partial contraction

Suppose we have one of the following simulation tasks: calculation of a single
amplitude, calculation of a batch of amplitudes or calculation of a set of
batches of amplitudes. In all these cases we need to calculate some components
of the full-state vector $\left|\psi\right\rangle=C\left|0\right\rangle$.
Suppose we also have a contraction tree $T$ for this task, a set of all sliced
vertices $I$, a set $S\subseteq I$ of partially sliced vertices, and a set
$X\subset\\{0,1\\}^{k}$ of slices such that $\|\psi_{X}\|^{2}=\mathcal{F}$. To
perform our simulation task with fidelity $\mathcal{F}$ it is enough to
calculate corresponding components of the vector
$\left|\overline{\psi_{X}}\right\rangle$ instead of $\left|\psi\right\rangle$.
This can be done as follows:

  1. 1.

Calculate the sum of all slices of the circuit $C$ with contraction tree $T$,
corresponding to the indices $X\times\\{0,1\\}^{|I|-|S|}$, i.e., when the
partially sliced vertices (from $S$) take values from $X$, the fully sliced
vertices (from $I\setminus S$) take all possible values.

  2. 2.

Normalize the result dividing by $\sqrt{\mathcal{F}}$.

##  II Sampling algorithm

###  II.1 Modified rejection sampling for RQCs

In the previous section, we described how to calculate some amplitudes or
batches of amplitudes for the state
$\left|\psi^{\prime}\right\rangle=\left|\overline{\psi_{X}}\right\rangle$ that
approximates the exact state $\left|\psi\right\rangle=C\left|0\right\rangle$
with the fidelity at least $f$. Here we describe the algorithm that we use to
simulate the RQCs from Google’s supremacy experiment. In this algorithm, we
assume that the output probability distribution does not have bitstrings with
very high probability. In this case, we can calculate some amplitudes and use
the modification of the frugal rejection sampling algorithm described below.

If target fidelity is $\ll 1$ as in Google’s experiment, then we cannot use
the approach proposed in [23] because it does not guarantee enough precision
to maintain the fidelity when it is already small. So, we need another
approach to reduce the number of batches needed to produce the given number of
samples. We suppose that the bitstring $b$ is sampled according to the
distribution corresponding to the measurement of the state
$\left|\psi^{\prime}\right\rangle$. We denote the corresponding measurement
operation by $M$. Let us separate the qubits into $2$ parts $A$ and $B$,
$N_{A}=2^{|A|}$, $N_{B}=2^{|B|}$. Then each bitstring $b$ can be represented
as a pair $(b_{A},b_{B})$, where $b_{A}\in[N_{A}]$, $b_{B}\in[N_{B}]$. Denote
the probability of the bitstring $(i,j)$ by
$p_{i,j}=P\\{M(\psi^{\prime})=(i,j)\\}$. Let $p_{j}$ be the probability that
$M_{B}(\psi^{\prime})=j$, then $p_{j}=\sum_{i\in[N_{A}]}p_{i,j}$.

At each step, we do the following.

  1. 1.

Select at random $j\in[N_{B}]$ and calculate the batch of amplitudes that
gives us the probabilities $p_{i,j}$ and hence $p_{j}$; $i\in[N_{A}]$,
$j\in[N_{B}]$.

  2. 2.

We accept this batch with probability
$t_{j}=\min\left(1,\frac{p_{j}N_{B}}{\alpha}\right)$, where $\alpha>1$ is a
parameter.

  3. 3.

If the batch is accepted, we sample one bitstring from this batch according to
the distribution where a bitstring $(i,j)$ has the conditional probability
$P\\{b=(i,j)\mid b_{B}=j\\}=P\\{b_{A}=i\mid b_{B}=j\\}=p_{i,j}/p_{j}$.

Since the average batch probability is $1/N_{B}$, the average number of
batches we need to get one sample is approximately equal to $\alpha$. If we
want to sample $k$ bitstrings, we should calculate approximately $\alpha k$
batches and then sample from them using the described algorithm. Note that we
should choose the parameter $\alpha$ using the knowledge about the bitstrings
distribution to be sure that $p_{j}<\alpha/N_{B}$ for almost all batches.

The probability $p^{\prime}_{i,j}$ that the sample is produced in one step and
it is the bitstring $b$ is the product of the probabilities $1/N_{B}$
($j=b_{B}$ at first step), $t_{j}$ (we accept batch) and $p_{i,j}/p_{j}$
($b=(i,j)$ conditioned on $b_{B}=j$), i.e.,

$p^{\prime}_{i,j}=\frac{1}{N_{B}}\min\left(1,\frac{p_{j}N_{B}}{\alpha}\right)\frac{p_{i,j}}{p_{j}}=\frac{1}{\alpha}\min\left(\frac{\alpha}{N_{B}},p_{j}\right)\frac{p_{i,j}}{p_{j}}.$

Let $p^{\prime}_{j}=\min(p_{j},\alpha/N_{B})$,
$\epsilon_{j}=p_{j}-p^{\prime}_{j}$,
$\epsilon=\sum_{j\in[N_{B}]}\epsilon_{j}$. The probability that a sample is
produced in one step is

$t=\sum_{i,j}p^{\prime}_{i,j}=\frac{1}{\alpha}\sum_{j\in[N_{B}]}p^{\prime}_{j}=\frac{1-\epsilon}{\alpha}.$

When some sample is produced, the probability that it is a bitstring $(i,j)$
is

$\tilde{p}_{i,j}=p^{\prime}_{i,j}/t=\frac{1}{1-\epsilon}p^{\prime}_{j}\frac{p_{i,j}}{p_{j}}.$

###  II.2 Statistical variational distance

Define
$\tilde{p}_{j}=\sum_{i=1}^{N_{A}}\tilde{p}_{i,j}=p^{\prime}_{j}/(1-\epsilon)$.
The _statistical variational distance_ between $p$ and $\tilde{p}$ is

$D(p,\tilde{p})=\frac{1}{2}\sum_{j=1}^{N_{B}}\sum_{i=1}^{N_{A}}|p_{i,j}-\tilde{p}_{i,j}|=\frac{1}{2}\sum_{j=1}^{N_{B}}|p_{j}-\tilde{p}_{j}|\\\
$ (3)

Since $\sum p_{j}=1=\sum\tilde{p}_{j}$ and $|x-y|=y-x+2\max(0,x-y)$, we have

$\sum_{j=1}^{N_{B}}|p_{j}-\tilde{p}_{j}|=2\sum_{j=1}^{N_{B}}\max(0,p_{j}-\tilde{p}_{j}).$
(4)

Taking into account $\tilde{p}_{j}\geq p^{\prime}_{j}$ and $p_{j}\geq
p^{\prime}_{j}$, from (3) and (4) we obtain

$D(p,\tilde{p})=\sum_{j=1}^{N_{B}}\max(0,p_{j}-\tilde{p}_{j})\\\
\leq\sum_{j=1}^{N_{B}}\max(0,p_{j}-p^{\prime}_{j})=\sum_{j=1}^{N_{B}}(p_{j}-p^{\prime}_{j})=\epsilon.$
(5)

We should choose the parameter $\alpha$ in order to make $\epsilon$ small
enough to obtain the given fidelity. Suppose we sample from state
$\left|\psi^{\prime}\right\rangle$ such that
$F(\left|\psi\right\rangle,\left|\psi^{\prime}\right\rangle)=|\left\langle\psi\middle|\psi^{\prime}\right\rangle|^{2}\geq
f$. Distance $D(p,\tilde{p})$ corresponds to trace distance between the
density matrix $M(\left|\psi^{\prime}\right\rangle)=\mathop{\mathrm{diag}}(p)$
and the density matrix $\mathop{\mathrm{diag}}(\tilde{p})$. We need to
estimate the fidelity
$f^{\prime}=F(M(\left|\psi\right\rangle),\mathop{\mathrm{diag}}(\tilde{p}))$.
Using the triangle inequality for Bures metric
$D_{B}(\rho,\sigma)=\sqrt{2\left(1-\sqrt{F(\rho,\sigma)}\right)}$, we have

$D_{B}(M(\left|\psi\right\rangle),\tilde{p})\leq
D_{B}(M(\left|\psi\right\rangle),p)+D_{B}(p,\tilde{p}).$ (6)

Since $\mathop{\mathrm{diag}}(p)=M(\left|\psi^{\prime}\right\rangle)$ and

$F(M(\left|\psi\right\rangle),M(\left|\psi^{\prime}\right\rangle))\geq
F(\left|\psi\right\rangle,\left|\psi^{\prime}\right\rangle)\geq f,$

we have

$D_{B}(M(\left|\psi\right\rangle),\tilde{p})\leq\sqrt{2(1-\sqrt{f})}.$ (7)

We will use following relation between fidelity and trace distance:
$1-\sqrt{F(p,\tilde{p})}\leq D(p,\tilde{p})$. For short denote
$d=D(p,\tilde{p})$, then

$D_{B}(p,\tilde{p})=\sqrt{2(1-\sqrt{F(p,\tilde{p})})}\leq\sqrt{2d}.$ (8)

Substituting (7), (8), and the definition of
$D_{B}(M(\left|\psi\right\rangle),\tilde{p})$ into (6) and dividing by
$\sqrt{2}$ we have

$\sqrt{1-\sqrt{f^{\prime}}}\leq\sqrt{1-\sqrt{f}}+\sqrt{d}.$

To obtain nontrivial estimation on $f^{\prime}$, we require that $d<f/16$. In
this case

$\sqrt{f^{\prime}}\geq\sqrt{f}-d-2\sqrt{d(1-\sqrt{f})}>\sqrt{f}-2\sqrt{d},$

Finally, we have an estimate

$f^{\prime}\geq f(1-4\sqrt{d/f}),$ (9)

where $d$ is the sampling algorithm error in terms of the trace distance.
Using (5) we can estimate $d$ from above by $\epsilon$ which we can estimate
in different ways.

For example, if we want to sample with fidelity $f^{\prime}=1\%$, then we can
set $f=1.1$ and choose $\alpha$ such that $\epsilon<5\cdot 10^{-6}$. Note that
when we estimate $\epsilon$, we always assume something about distribution $p$
(for example, for random quantum circuit with enough depth components of $p$
have Porter-Thomas distribution).

$0$$5$$10$$15$$20$$m=18$Single bitstrings

$0.5$$1$$1.5$$2$Batches

$0$$5$$10$$15$$20$$m=20$

$0.5$$1$$1.5$$2$

Figure 2: The distribution of the normalized probabilities for individual
bitstrings ($x$ axis is $p_{i,j}2^{n}$) and batches of size 64 ($x$ axis is
$p_{j}N_{B}$).

###  II.3 Estimation of $\epsilon$

One way to estimate $\epsilon$ is to use the assumption that the vector
$\left|\psi^{\prime}\right\rangle$ has Porter-Thomas distribution. In this
case, the probabilities $p_{i,j}$ are almost independent and have exponential
distribution $p_{i,j}\sim\mathrm{Exp}(2^{n})$. The probability $p_{j}$ of a
batch of size $N_{A}$ is sum of $N_{A}$ independent exponentially distributed
random variables and has gamma-distribution
$p_{j}\sim\mathrm{Gamma}(N_{A},2^{n})$. We can estimate error expectation

$\mathbb{E}\epsilon_{j}=\frac{\Gamma(N_{A},2^{n}\alpha/N_{B})}{\Gamma(N_{A})}=\frac{\Gamma(N_{A},\alpha
N_{A})}{\Gamma(N_{A})}$

where $\Gamma(s,x)=\int_{x}^{\infty}t^{s-1}\,e^{-t}\,{\rm{d}}t$ is the _upper
incomplete gamma function_. Hence,

$\mathbb{E}\epsilon=\sum_{j=1}^{N_{B}}\mathbb{E}\epsilon_{j}=N_{B}\frac{\Gamma(N_{A},\alpha
N_{A})}{\Gamma(N_{A})}.$ (10)

On figure 2 we compare actual distribution of bitstring probabilities with
exponential distribution and also compare distribution of batch probabilities
with gamma distribution.

Second way to estimate $\epsilon$ is to calculate $\epsilon_{J}=\sum_{j\in
J}\epsilon_{j}$ for set $J$ for which batches were calculated during
algorithm. Then calculate estimation $\tilde{\epsilon}=\epsilon_{J}/|J|$. This
method doesn’t use assumption about concrete distribution however works in
assumption that all probabilities have the same order of magnitude which is
true for RQCs of enough depth.

##  III Experiments

###  III.1 Sampling with fidelity

####  III.1.1 Experiment with $\left|\psi_{i}\right\rangle$ norm calculation

$0.5$$1$$1.5$$2$$k=10$$m=14$

$0.98$$1$$1.02$$m=16$

$0.95$$1$$1.05$$m=18$

$0.95$$1$$1.05$$m=20$

$0.5$$1$$1.5$$2$$k=15$

$0.98$$1$$1.02$

$0.95$$1$$1.05$

$0.95$$1$$1.05$

Figure 3: The distribution of $\|\left|\psi_{i}\right\rangle\|^{2}$. The
labels on $X$ axis shows the normalized values
$2^{k}\|\left|\psi_{i}\right\rangle\|^{2}$. Note that the mean values with
this normalization are equal to $1$.

The complexity of simulation with fidelity depends on how many slices should
be contracted to obtain given fidelity. The number of slices depends on how
does maximal $\|\left|\psi_{i}\right\rangle\|$ norm differ from mean
$\|\left|\psi_{i}\right\rangle\|$ norm. In Fig. 3 there are calculated norms
of $\|\left|\psi_{i}\right\rangle\|^{2}/2^{k}$ for schedules used for sampling
1M amplitudes for sycamore RQCs. There were 2 settings: $k=10$ and $k=15$
indices involved in partial slicing. In Table 1 there are complexity of
$\|\left|\psi_{i}\right\rangle\|$ calculation, normalized standard deviations
and range of $\|\left|\psi_{i}\right\rangle\|^{2}$. Normalization here is
multiplication by $2^{k}$ to make mean value equal to 1.

$m$ | $k$ | Complexity | $\sqrt{2^{k}D\|\left|\psi_{i}\right\rangle\|^{2}}$ | $\mathop{\mathrm{range}}(2^{k}\|\left|\psi_{i}\right\rangle\|^{2})$
---|---|---|---|---
20 | 10 | $1.2\times 10^{15}$ | 0.021 | [0.95,1.05]
20 | 15 | $1.6\times 10^{17}$ | 0.022 | [0.94,1.06]
18 | 10 | $2.1\times 10^{12}$ | 0.017 | [0.96,1.04]
18 | 15 | $2.3\times 10^{12}$ | 0.021 | [0.93,1.07]
16 | 10 | $2.7\times 10^{11}$ | 0.007 | [0.97,1.02]
16 | 15 | $5.6\times 10^{14}$ | 0.013 | [0.96,1.04]
14 | 10 | $8.7\times 10^{6}$ | 0.341 | [0.50,1.88]
14 | 15 | $5.4\times 10^{10}$ | 0.341 | [0.49,1.93]
12 | 8 | $3.6\times 10^{12}$ | 0.135 | [0.78,1.26]
Table 1: Complexity, deviation and range of
$\|\left|\psi_{i}\right\rangle\|^{2}$.

From Fig. 3 and Table 1 we see that the distribution of
$\|\left|\psi_{i}\right\rangle\|^{2}$ is far from normal is most cases, and
sometimes variance is very big, so it cannot be assumed that all
$\|\left|\psi_{i}\right\rangle\|^{2}$ are close to $2^{-k}$, and we should
find it directly. In some cases such as $m=14$ maximal value 2 times bigger
than average, this allows to calculate almost 2 times less slices to obtain
fidelity 0.02 (compare slicing ratio and fidelity $\mathcal{F}$ in Table 3).
For $m\geq 16$ variance is not so big, and slicing ratio is almost equal to
fidelity. Probably, this is because sliced vertices for $m\geq 16$ are closer
to the middle of circuit and subcircuit $C_{1}$ can be viewed as a random
circuit with enough depth.

####  III.1.2 Sampling algorithm validation on elided circuits

In this section, we verify our sampling algorithm on the elided circuits from
Google’s experiment [6]. Since our method adopts the frugal rejection
sampling, to work well it requires some special properties of the output
probability distribution. In particular, it works well when there is a very
small number of output bitstrings that have probabilities significantly bigger
than $1/2^{n}$, where $n$ is the number of qubits. Google’s team has already
provided some data in [6] that confirms this assumption for supremacy RQCs. In
this section, we demonstrate our algorithm on the elided circuits for which we
can calculate the exact amplitudes and $\mathcal{F}_{\mathrm{XEB}}$. We apply
our algorithm with target fidelity 1% and compare the obtained
$\mathcal{F}_{\mathrm{XEB}}$ with 0.01. The calculated fidelities are shown in
Fig. 2.

$m$ | 12 | 14 | 16 | 18 | 20
---|---|---|---|---|---
$\mathcal{F}$ | 0.0104 | 0.0103 | 0.0117 | 0.0104 | 0.0113
$\mathcal{F}_{\mathrm{XEB}}$ | 0.0108 | 0.0099 | 0.0121 | 0.0103 | 0.0114
Table 2: The fidelities obtained for the Google’s elided circuits when the
target fidelity is 1%. For each number of cycles $m$ and for all 10 elided
circuits provided by Google 1M samples was generated. The fidelities in the
table are averaged over these 10 circuits.

Experiment description. For each number of cycles $m=12,14,16,18,20$ and each
of 10 elided circuit instances provided by Google we calculate $2^{21}$ random
batches of size 64 with target fidelity 1%. Note that the actual fidelity
$\mathcal{F}=\sum_{i\in X}\|\left|\psi_{i}\right\rangle\|^{2}$, shown in table
2, is slightly higher. Then we apply our modified frugal rejection sampling
algorithm to calculated batches and obtain slightly more than $10^{6}$
samples. After we get all samples, we calculate the amplitudes for these
samples with 100% fidelity using the multi-tensor contraction algorithm from
[19]. After we get all amplitudes, we calculate $\mathcal{F}_{\mathrm{XEB}}$
for all 10 circuit instances. The average $\mathcal{F}_{\mathrm{XEB}}$ for
each number of cycles is shown in table 2.

####  III.1.3 Sampling supremacy circuits

In table 3 shown sampling complexity for Google’s RQC circuits with different
number of cycles where time is shown for 1 GPU Tesla V100.

$m$ |  $2^{2}$target fidelity$2_{q}$ | $\mathcal{F}$ |  slicing ratio | $C_{s}$ | efficiency | time
---|---|---|---|---|---|---
12$2^{2^{2}}$ | 0.02 | 0.0215 | 0.0136 | $1.1\cdot 10^{16}$ | 26% | 6.6 hours
14 | 0.02 | 0.0218 | 0.0117 | $5.4\cdot 10^{16}$ | 42% | 20.5 hours
16 | 0.02 | 0.0208 | 0.0205 | $4.0\cdot 10^{17}$ | 52% | 5 days
18 | 0.02 | 0.0201 | 0.0195 | $1.9\cdot 10^{18}$ | 32% | 40 days
20 | 0.002 | 0.0021 | 0.0021 | $2.2\cdot 10^{19}$ | 31% | 15 months
Table 3: Complexity of sampling 1M samples with fidelity. $C_{s}$ is
complexity measure, number of complex number multiplications during
contraction. Total number of flops is $8C_{s}$. Time and efficiency provided
for Tesla V100 GPU with 16 GB memory.

$0$$1$$2$$\cdot 10^{-2}$$0$$2$$4$$\cdot 10^{-2}$Ratio$m=12$

$0$$1$$2$$\cdot 10^{-2}$$0$$2$$4$$\cdot 10^{-2}$Ratio$m=14$

$0$$0.5$$1$$\cdot 10^{-2}$$0$$1$$2$$\cdot 10^{-2}$Ratio$m=16$

$0$$0.5$$1$$\cdot 10^{-2}$$0$$1$$2$$\cdot 10^{-2}$Ratio$m=18$

$0$$0.5$$1$$\cdot 10^{-2}$$0$$1$$2$$\cdot 10^{-2}$Ratiom=2010%20%36%50%

Figure 4: XEB spoofing with partial slicing

Experiment description. For each number of cycles $m=12,14,16,18,20$ we take
first full circuit instance provided by Google and calculate $2^{21}$ batches
of size 64 with target fidelity $2\%$ for $m=12,14,16,18$ and $0.2\%$ for
$m=20$, which is not less than as Google’s sycamore quantum computer has on
the same circuit (see target and actual fidelities in table 3). Then 1M
samples was generated for each $m$ using proposed modified frugal rejection
sampling algorithm. Experimental data [20] contains 5 text files, each
contains 1M bitstrings for corresponding $m$.

Simulation run on 4 servers, each has 8 GPU Tesla V100 16GB, total running
time for all cases is approximately 14.5 days. In table 3 there is detailed
information about each case complexity and running time normalized for one GPU
Tesla V100. Estimated time of generating 1M samples for $m=20$ on Summit
supercomputer is approximately 24 min.

For the recent experiment on Zuchongzhi quantum computer [7] the sampling task
is significantly more complex. We prepared a contraction tree for generating
2M batches of size 64 for a 56-qubit circuit with 20 cycles. The full
contraction complexity is $2.15\cdot 10^{25}$ Flops, the contraction requires
80 GB memory. We estimated the time for sampling with target fidelity 0.066%
on Selene supercomputer with 4480 Tesla A100 80GB. Note that the contraction
schedule has high arithmetic intensity. Hence in this estimation we assume
that the computational efficiency of our simulator will be at least 50%. The
time estimate in this case is:

$\frac{(2.15\cdot 10^{25})\mathrm{FLOPs}\times 0.066\%}{(79.2\cdot
10^{15})\mathrm{FLOP}/\mathrm{s}\times 50\%}\approx 3.6\cdot
10^{5}\,\mathrm{s}\approx 4\,\mathrm{days}.$

Taking into account a relatively high arithmetic intensity, some additional
optimization can be done using tensor cores with single precision. Therefore
potentially the simulation time can be reduced to 1 day or even less.

###  III.2 Spoofing Linear XEB

In this section, we show our experimental results on the spoofing Linear XEB
test. If we want to get a set of $N$ bitstrings with
$\mathcal{F}_{\mathrm{XEB}}\geq f$, then we need to do the following steps:

  1. 1.

Let $b=\lceil\log_{2}(10N)\rceil$.

  2. 2.

Choose the set $F$ consisting of $b$ free circuit outputs which gives the
minimal contraction complexity.

  3. 3.

Calculate the amplitudes for batch $B$ of $2^{b}$ bitstrings using partial
slicing with target fidelity $f$.

  4. 4.

Select the $N$ bitstrings from the batch $B$ with maximal absolute values of
amplitudes.

The computational cost of spoofing the XEB test in Google’s quantum supremacy
experiment for different number of cycles $m$ is shown in table 4. We also
estimated the time to spoof the Linear XEB test in the recent experiment with
56-qubit circuit on Zuchongzhi quantum computer [7]. Our estimates show that
in this case the spoofing can be done in $1$ month on one Tesla V100 16GB GPU.

Assume that we have random circuit $C$ on $n$ qubits, random variables
$p_{i}=|\left\langle i\right|C\left|0\right\rangle|^{2}$ have Porter-Thomas
distribution. Moreover, if circuit depth is enough, $C\left|0\right\rangle$ is
uniformly distributed on complex sphere $S^{2^{n}}$, for given set of
bitstrings $B$, $|B|\ll 2^{n}$ we can assume that random variables $p_{i},i\in
B$ are almost independent. If we calculated approximate values
$p^{\prime}_{i}$ with fidelity $f$ and take set $S$ of $N=r|B|$ bitstrings
with maximal values of $p^{\prime}_{i}$, then

$\mathop{\mathbb{E}}\mathcal{F}_{\mathrm{XEB}}(S)=-f\ln r+O(N^{-1}\ln N),$
(11)

see Appendix A for details. For real circuits assumption that $p_{i}$ are
independent is not always true especially when $B$ is a batch of corellated
bitstrings. In random subsets of a batch we can assume that all $p_{i}$ are
independent but distribution of $p_{i}$ is exponential with expectation that
depends on $\mathcal{F}_{\mathrm{XEB}}(B)$. Taking into account this fact we
have following heuristic equality

$\mathop{\mathbb{E}}(\mathcal{F}_{\mathrm{XEB}}(S)-\mathcal{F}_{\mathrm{XEB}}(B))\approx-f\ln
r.$ (12)

$m$ | XEB | slicing ratio |  $2^{25}$ amps batch complexity |  spoofing time
---|---|---|---|---
12 | 3.2% | 1.4% | $7.5\cdot 10^{13}$ | 1 sec
14 | 1.2% | 0.9% | $1.6\cdot 10^{14}$ | 2 sec
16 | 1.4% | 0.6% | $1.7\cdot 10^{17}$ | 17 min
18 | 0.83% | 0.4% | $1.1\cdot 10^{18}$ | 30 min
20 | 0.47% | 0.2% | $6.9\cdot 10^{18}$ | 4 hours
Table 4: Linear XEB spoofing (the time is for _one_ GPU Tesla V100 16GB)

Experiment description. For each number of cycles $m=12,14,16,18,20$ we
calculate batch $B_{m}(f)$ of $2^{25}$ amplitudes with different ratio of
slices $f$ from 0.1% to 2% and obtain partially calculated probability
distributions. From each batch of partially calculated amplitudes $B_{m}(f)$
we take the set $S_{m}(f,r)$ of $\lfloor r\cdot 2^{25}\rfloor$ samples with
maximal absolute values where $r\in\\{0.1,0.2,0.36,0.5\\}$ (0.36 is
approximate value of $e^{-1}$). Then we calculate precise amplitudes for the
same batch and which we use to calculate $\mathcal{F}_{\mathrm{XEB}}$ for all
sets of samples $S_{m}(f,r)$. Calculated $\mathcal{F}_{\mathrm{XEB}}$ for all
cases is shown on figure 4.

We can see that
$\mathcal{F}_{\mathrm{XEB}}(S_{m}(f,0.36))\approx\mathcal{F}_{\mathrm{XEB}}(S_{m}(0,1))+f$
where $\mathcal{F}_{\mathrm{XEB}}(S_{m}(0,1))$ corresponds to whole batch. In
the case $m=14$ $\mathcal{F}_{\mathrm{XEB}}$ of full batch $\approx-0.0094$
which means some weakness in entanglement of qubits in the output state of
circuit with 14 cycles. For bigger $m$ there are no such problems and
$\mathcal{F}_{\mathrm{XEB}}$ of full batch is close to 0.

## References

  * Shor [1996] P. Shor, in  _Proceedings of 37th Conference on Foundations of Computer Science_ (1996) pp. 56–65.
  * Kalai [2021] G. Kalai, The argument against quantum computers, the quantum laws of nature, and Google’s supremacy claims (2021), arXiv:2008.05188 [quant-ph] .
  * Preskill [2012] J. Preskill, arXiv:1203.5813 [cond-mat, physics:quant-ph] (2012).
  * Aaronson and Chen [2017] S. Aaronson and L. Chen, in  _Proceedings of the 32nd Computational Complexity Conference_ , CCC ’17 (Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, DEU, 2017) pp. 1–67.
  * Yung [2019] M.-H. Yung, National Science Review 6, 22 (2019).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, and others, Nature 574, 505 (2019).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, and et al., Physical Review Letters 127, 10.1103/physrevlett.127.180501 (2021).
  * Zhu  _et al._ [2021] Q. Zhu, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, M. Gong, C. Guo, C. Guo, S. Guo, L. Han, L. Hong, and others, Quantum computational advantage via 60-qubit 24-cycle random circuit sampling (2021), arXiv:2109.03494 [quant-ph] .
  * Markov  _et al._ [2018] I. L. Markov, A. Fatima, S. V. Isakov, and S. Boixo, arXiv:1807.10749 [quant-ph] (2018).
  * Gray and Kourtis [2021] J. Gray and S. Kourtis, Quantum 5, 410 (2021).
  * Huang  _et al._ [2020] C. Huang, F. Zhang, M. Newman, J. Cai, X. Gao, Z. Tian, J. Wu, H. Xu, H. Yu, B. Yuan, M. Szegedy, Y. Shi, and J. Chen, arXiv:2005.06787 [quant-ph] (2020).
  * Guo  _et al._ [2021] C. Guo, Y. Zhao, and H.-L. Huang, Phys. Rev. Lett. 126, 070502 (2021).
  * Pan and Zhang [2021] F. Pan and P. Zhang, arXiv:2103.03074 [physics, physics:quant-ph] (2021).
  * Liu  _et al._ [2021] Y. A. Liu, X. L. Liu, F. N. Li, H. Fu, Y. Yang, J. Song, P. Zhao, Z. Wang, D. Peng, H. Chen, C. Guo, H. Huang, W. Wu, and D. Chen, in  _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_, SC ’21 (Association for Computing Machinery, New York, NY, USA, 2021) pp. 1–12.
  * Pan  _et al._ [2021] F. Pan, K. Chen, and P. Zhang, arXiv preprint arXiv:2111.03011 (2021).
  * Zhou  _et al._ [2020] Y. Zhou, E. M. Stoudenmire, and X. Waintal, Physical Review X 10, 041038 (2020).
  * Pednault  _et al._ [2020] E. Pednault, J. A. Gunnels, G. Nannicini, L. Horesh, T. Magerlein, E. Solomonik, E. W. Draeger, E. T. Holland, and R. Wisnieff, arXiv:1710.05867 [quant-ph] (2020).
  * Schutski  _et al._ [2020] R. Schutski, D. Lykov, and I. Oseledets, Physical Review A 101, 042335 (2020), arXiv:1911.12242 .
  * Kalachev  _et al._ [2021] G. Kalachev, P. Panteleev, and M.-H. Yung, Recursive multi-tensor contraction for xeb verification of quantum circuits (2021), arXiv:2108.05665 [quant-ph] .
  * [20] https://gitee.com/Huawei-HiQ/supremacy.
  * Bienstock [1990] D. Bienstock, Journal of Combinatorial Theory, Series B 49, 103 (1990).
  * O’Gorman [2019] B. O’Gorman, in  _14th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2019)_, Leibniz International Proceedings in Informatics (LIPIcs), Vol. 135, edited by W. van Dam and L. Mancinska (Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany, 2019) pp. 10:1–10:19.
  * Villalonga  _et al._ [2019] B. Villalonga, S. Boixo, B. Nelson, C. Henze, E. Rieffel, R. Biswas, and S. Mandrà, npj Quantum Information 5, 1 (2019).

##  Appendix A Theoretical XEB estimation for spoofing

Assume that all $p_{i}$ are independent and $p_{i}\sim\mathrm{Exp}(\lambda)$
where $\lambda=2^{n}$. If we consider $N=|B|$ random variables
$x_{k}\sim\mathrm{Exp}(\lambda)$, $k=\overline{1,N}$, then for order
statistics $x_{(k)}$ we have

$\mathop{\mathbb{E}}x_{(N-k+1)}=\frac{1}{\lambda}(\ln N-\ln k+O(1/k)),$
$\mathop{\mathbb{E}}\sum_{j=1}^{k}x_{(N-j+1)}=\frac{1}{\lambda}\left(k\left(\ln\frac{N}{k}+1\right)+O(\ln
k)\right).$

Suppose we select the set $S$ of $rN$ samples with the maximal probabilities
from a batch of size $N$, then

$\displaystyle\mathop{\mathbb{E}}\mathcal{F}_{\mathrm{XEB}}(S)=\mathop{\mathbb{E}}\left(\frac{2^{n}}{rN}\sum_{j=1}^{rN}x_{(N-j+1)}-1\right)$
$\displaystyle\qquad=\frac{2^{n}}{\lambda}\left(\ln N-\ln rN+1+O((rN)^{-1}\ln
rN)\right)-1$ $\displaystyle\qquad=-\ln r+O(N^{-1}\ln N).$

For partial slicing we have orthogonal projection $v_{f}$ of state vector $v$
on some subspace with the norm $\|v_{f}\|^{2}=f$. We assume that $v_{f}$ has
the uniform distribution on the sphere of radius $\sqrt{f}$, so the components
$p_{f,i}=|v_{f,i}|^{2}$ have exponential distribution with the parameter
$2^{n}/f$. We have $v=v_{f}+v_{f}^{\bot}$, $\left\langle
v_{f}\middle|v_{f}^{\bot}\right\rangle=0$. Moreover, the random vectors
$v_{f}$ and $v_{f}^{\bot}$ are independent. It is not hard to see that

$\mathop{\mathbb{E}}(p_{i}\mid
p_{f,i})=\mathop{\mathbb{E}}p_{f,i}^{\bot}+p_{f,i}=p_{f,i}+\frac{1-f}{2^{n}}.$

Let $S=\\{i_{j}\mid j=1,...,rN\\}$ be the set of $rN$ indices $i$ with the
maximal $p_{f,i}$, then

$\displaystyle\mathop{\mathbb{E}}\mathcal{F}_{\mathrm{XEB}}(S)$
$\displaystyle=\mathop{\mathbb{E}}\left(\frac{2^{n}}{rN}\sum_{j=1}^{rN}p_{i_{j}}-1\Biggm{|}p_{f,i_{j}}=x_{(N-j+1)}\right)$
$\displaystyle=\frac{2^{n}}{rN}\sum_{j=1}^{rN}\mathop{\mathbb{E}}\left(p_{i_{j}}\bigm{|}p_{f,i_{j}}=x_{(N-j+1)}\right)-1$
$\displaystyle=\frac{2^{n}}{rN}\left(rN\frac{1-f}{2^{n}}+\sum_{j=1}^{rN}x_{(N-j+1)}\right)-1$
$\displaystyle=1-f+f\left(\ln\frac{N}{rN}+1+O\left(\frac{\ln
N}{N}\right)\right)-1$ $\displaystyle=-f\ln r+O(N^{-1}\ln N).$

##  Appendix B Formal slice selection algorithm

Let $C$ be quantum circuit, then by $G(C)$ we denote the set of all its gates.
Each gate inside the circuit is represented by a tuple $(\ell,A,q)$, where
$\ell$ is the index of a gate in the circuit, $A$ is a unitary matrix, $q$ is
a tuple of qubit indices. If $X$ is a set of gates, then by $\mathbf{C}(X)$ we
denote the circuit composed from these gates. By vertices of the circuit we
understand the tensor legs in the corresponding tensor network. The lightcone
$L_{C}(v)$ of a vertex $v\in V(C)$ contains all the gates the vertex $v$
depends on. Denote by $L^{\prime}_{C}(v)$ the set of all vertices that are the
inputs of the gates from $L_{C}(v)$. For a set $S\subset V(C)$ we can define
$L_{C}(S)=\cup_{v\in S}L_{C}(v)$, $L^{\prime}_{C}(S)=\cup_{v\in S}L_{C}(v)$.
Algorithm 1 selects the subset $S$ of partially sliced vertices from the set
$I$ of all sliced vertices.

$k_{0}:=\lceil 3-\log_{2}f\rceil$;

$S:=\mathrm{SlicedVertexSelect}(C,I,k_{0})$;

$k:=|S|$;

$C_{1}:=\mathbf{C}(L_{C}(S))$;

// Find tensor $R$, $R[i]=\|\left|\psi_{i}\right\rangle\|^{2}$ for
$i\in\\{0,1\\}^{k}$

$R:=\bigl{(}\left\langle
0\right|C_{1}^{*}\delta_{S=i}C_{1}\left|0\right\rangle\bigr{)}_{i\in\\{0,1\\}^{k}}$;

Sort $i\in\\{0,1\\}^{k}$ such that $R[i^{(1)}]\geq...\geq R[i^{(2^{k})}]$;

$j:=\min\\{t\in\\{1,...,2^{k}\\}:\sum_{j=1}^{t}R[i^{(j)}]\geq f\\}$;

return $(S,\\{i^{(1)},...,i^{(j)}\\})$;

Algorithm 1 Slice selection

The operator $\delta_{S=i}$ is a projector on the subspace where the qubit
corresponding to the vertex $s_{j}$ is $\left|i_{j}\right\rangle$. In the
tensor network it corresponds to fixing indices $s_{1}=i_{1}$,…,$s_{k}=i_{k}$.

The vertex subset $S$ is selected by the following algorithm that optimizes
the gate number in the lightcone $L(S)$ in a greedy fashion.

$S:=\varnothing$;

while  _$|S| <k$ and $I\setminus L^{\prime}_{C}(S)\neq\varnothing$_ do

$v:=\arg\min_{u\in I\setminus
L^{\prime}_{C}(S)}|L^{\prime}_{C}(S\cup\\{v\\})|$;

$S:=(S\setminus L^{\prime}_{C}(v))\cup\\{v\\}$;

return $S$;

Algorithm 2 SlicedVertexSelect($C,I,k$)


††thanks: Corresponding author

# The motion of twisted particles in a stellar gravitational field

Dan-Dan Lian<EMAIL_ADDRESS>Peng-Ming Zhang
<EMAIL_ADDRESS>School of Physics and Astronomy, Sun Yat-sen
University, 519082 Zhuhai, China

###### Abstract

A twisted particle possesses intrinsic orbital angular momentum (OAM), and the
dynamics of such a particle may challenge the Einstein Equivalence Principle.
In this study, we disregard the spin characteristic of the twisted particle,
modeling it as a massless complex twisted scalar wave packet to simplify its
interaction with gravitational fields. Building on this simplification, we
investigate the gravitational birefringence of the twisted particle by
analyzing the center of its energy density. We demonstrate that the
gravitational birefringence induced by OAM can potentially exceed that induced
by spin by several orders of magnitude, significantly enhancing its
detectability. Furthermore, we examine the influence of the nonminimal
coupling term $\lambda R|\phi|^{2}$ on the propagation of the twisted particle
through the internal gravitational field of a star. Contrary to the
predictions of the Mathisson-Papapetrou-Dixon equations, our findings show
that the trajectory of the twisted particle under nonminimal coupling differs
from that in the minimal coupling scenario. Specifically, we find that for a
positive nonminimal coupling constant, the trajectory of the twisted particle
is expected to deviate away from the stellar center, compared to the minimal
coupling scenario, and this deviation is independent of the particle’s OAM.
These findings could provide new avenues for testing the Einstein Equivalence
Principle.

##  I Introduction

Undoubtedly, in the field of modern physics, General Relativity (GR) holds its
place as the most widely accepted theory of gravity yet. The experimental
testing of its foundational postulate, the Einstein Equivalence Principle
(EEP), continues to be a vibrant area of research. The dynamics of spinning
particles are proposed to be accurately described by using the Mathisson-
Papapetrou-Dixon (MPD) equations Mathisson (2010); Papapetrou (1951); Dixon
(1970a, b); Puetzfeld  _et al._ (2015). However, these MPD equations seem to
challenge the EEP. And an intriguing phenomenon was first discovered by
Corinaldesi and Papapetrou, who found that particles with different angular
momentum would follow differing trajectories in a gravitational field
Corinaldesi and Papapetrou (1951). This phenomenon, due to its similarity to
the birefringence of light observed in certain materials, is commonly referred
to as gravitational birefringence. This suggests that the motion of particles
carrying angular momentum could violate the EEP, potentially providing a new
avenue for the experimental testing of GR.

This phenomenon, gravitational birefringence, has been extensively discussed
by using different methods. Several of these studies have approached this
phenomenon by solving the MPD equations Hartl (2003); Hackmann  _et al._
(2014); Duval and Schücker (2017); Duval  _et al._ (2019). In contrast, other
researches discuss this phenomenon by appling the Wentzel-Kramers-Brillouin
(WKB) approximation to the Maxwell equations in GR Frolov and Shoom (2012);
Oancea  _et al._ (2020); Deriglazov (2021); Harte and Oancea (2022); Andersson
and Oancea (2023). Moreover, this phenomenon has been further investigated
through the lens of the Energy-Momentum Tensor (EMT) Lian and Chen (2022). In
another study, this phenomenon for free-falling electrons in a static uniform
gravitational field has also been explored Wang . However, the primary focus
of these studies has been on spin angular momentum, with its associated
gravitational birefringence being too small to be detected.

It is noteworthy that particles such as electrons, photons, atoms, and
neutrons can exist in twisted states, which possess intrinsic orbital angular
momentum (OAM) Ivanov (2022). For photons in the lab environment, the OAM can
exceed $10^{4}\hbar$ per photon Fickler  _et al._ (2016), and even reach up to
$10^{6}\hbar$ per photon Zhang  _et al._ (2021). Furthermore, it has been
reported that light emitted from a Kerr black hole can also carry OAM
Tamburini  _et al._ (2011); Yang and Casals (2014). This OAM has been claimed
to be detected through the analysis of Event Horizon Telescope data Tamburini
_et al._ (2020, 2021). According to the MPD equations, gravitational
birefringence is proportional to the angular momentum $l$. Consequently, this
phenomenon induced by intrinsic OAM could be significantly larger than that
induced by spin. Therefore, it seems more practical to utilize the twisted
particles to test the EEP.

The primary focus of this work is to explore the influence of OAM on the
motion of a twisted particle. Accordingly, we introduce a simplification by
disregarding the spin characteristic of the twisted particle and representing
it as a massless twisted complex scalar wave packet, thus facilitating the
analysis of its interactions with gravitational fields. Building on this
simplification, we investigate the behavior of a free-falling twisted particle
in a stellar gravitational field. Within the framework of GR, the wave
function $\phi$ of the twisted particle is determined by the covariant Klein-
Gordon equation via the Minimal Coupling Procedure (MCP), in accordance with
the EEP. However, it has been pointed out that the EMT of a scalar field
cannot be finite at one-loop order unless a nonminimal coupling term $\lambda
R|\phi|^{2}$ is incorporated in its Lagrangian density Freedman  _et al._
(1974); Freedman and Weinberg (1974). Here $\lambda$ is a constant and $R$
denotes the Ricci scalar, and $|\phi|^{2}$ signifies the square of the modulus
of $\phi$. This nonminimal coupling term arises from quantum corrections Linde
(1982) and is anticipated when the spacetime curvature is large Faraoni
(1996). The associated EMT $T^{\mu\nu}$, defined as $\delta
I\equiv\frac{1}{2}\int d^{4}x\sqrt{g}T^{\mu\nu}\delta g_{\mu\nu}$, also
exhibits covariant conservation $\nabla_{\mu}T^{\mu\nu}=0$, where $I$ denotes
as the action of the scalar field. And this conservation satisfies the
requirement of the MPD equations. Therefore, the motion of the twisted
particle can be also described by the MPD equations, irrespective of the
constant $\lambda$. Importantly, the applicability of the MPD equations in
this nonminimal coupling case merely indicates the feasibility of defining a
consistent trajectory regardless of the value of the constant $\lambda$.
However, it does not negate the nonminimal coupling effects on the motion of
the twisted particle.

Notably, the nonminimal coupling term $\lambda R|\phi|^{2}$ vanishes when the
twisted particle moves outside gravitational sources, as the Ricci scalar
$R=0$. Therefore, this work concentrates on the motion of the twisted particle
inside stars, which are ubiquitous gravitational sources in the universe. We
exclude all non-gravitational couplings, considering them outside the scope of
this investigation. As previously mentioned, the MPD equations are not
suitable for discussing the effects of the nonminimal couplings. Consequently,
we employ the center of the EMT to describe the motion of the twisted
particle. However, the EMT has 16 components, each leading to different
trajectories as shown in Ref. Lian and Chen (2022). Additionally, the
gravitational birefringence given by the energy flux density is also affected
by another effect, named the Geometric Spin Hall Effect (GSHE) for its purely
geometric nature Aiello  _et al._ (2009); Wang and Chen (2022); Lian and Chen
(2022). However, the gravitational birefringence resulting from the energy
density remains unaffected, as noted in Ref. Lian and Chen (2022). Compared
with other components of the EMT, this unaffected property by GSHE suggests
that the energy density has an advantage in investigating the couplings
between the OAM and gravitational fields. Consequently, our focus in this
paper is primarily on the energy density. We reiterate that the topic of this
work is to investigate the motion of a twisted particle and discuss the
effects of the nonminimal coupling term $\lambda R|\phi|^{2}$ on its motion.

In this paper, we have set natural units where $\hbar=c=1$ and adopted the
metric signature $(-,+,+,+)$. The curvature tensor, defined in terms of the
Christoffel symbols $\Gamma^{\alpha}_{\mu\nu}$, is expressed as $R^{\alpha}_{\
\mu\nu\beta}\equiv\partial_{\beta}\Gamma^{\alpha}_{\mu\nu}-\partial_{\nu}\Gamma^{\alpha}_{\mu\beta}+\Gamma^{\gamma}_{\mu\nu}\Gamma^{\alpha}_{\beta\gamma}-\Gamma^{\gamma}_{\mu\beta}\Gamma^{\alpha}_{\nu\gamma}$.
Furthermore, the Ricci tensor and Ricci scalar are denoted by
$R_{\mu\beta}\equiv R^{\gamma}_{\ \mu\gamma\beta}$ and $R\equiv
g^{\mu\nu}R_{\mu\nu}$, respectively. Greek indices run over the four
coordinate labels in a general coordinate system, and Latin indices
$i,j,k,\dots$ are restricted to three spatial coordinate labels.

This paper is organized as follows: Section II provides an overview of the
methodology used to investigate the motion of a twisted particle in a
gravitational field. This section also introduces the configuration of the
physical system and the numerical method used to solve the equations of motion
(EoM) for a free-falling twisted particle. Section III calculates the center
of energy density of the twisted particle to describe its motion and details
the gravitational birefringence induced by its OAM. Section IV compares the
trajectory of the twisted particle under minimal and nonminimal coupling
scenarios and discusses the effects of nonminimal couplings. Section V offers
further discussion of our findings.

##  II Equation of motion of free-falling twisted particles

As mentioned in Section I, we represent twisted particles as massless complex
twisted scalar wave packets to simplify the analysis of their interactions
with gravitational fields. The simplest and most reasonable nonminimal
coupling term for complex scalar fields is given by $\lambda R|\phi|^{2}$. In
this case, the Lagrangian of a free-falling twisted particle becomes

$\mathcal{L}=-\frac{1}{2}\nabla_{\mu}\phi\nabla^{\mu}\phi^{*}-\frac{1}{2}\lambda
R|\phi|^{2},$ (1)

where $\phi^{*}$ denotes the complex conjugate of $\phi$. The EoM of the
twisted particle, determined by the principle of stationary action, is given
by

$\nabla_{\mu}\nabla^{\mu}\phi-\lambda R\phi=0,$ (2)

and its associated EMT takes the following form:

$T^{\mu\nu}=\frac{1}{2}\nabla^{\mu}\phi\nabla^{\nu}\phi^{*}+\frac{1}{2}\nabla^{\mu}\phi^{*}\nabla^{\nu}\phi-\frac{1}{2}g^{\mu\nu}\nabla^{\rho}\phi\nabla_{\rho}\phi^{*}+\lambda(G^{\mu\nu}-g^{\mu\nu}\nabla_{\rho}\nabla^{\rho}+\nabla^{\mu}\nabla^{\nu})(\phi\phi^{*}),$
(3)

where $G^{\mu\nu}\equiv R^{\mu\nu}-\frac{1}{2}g^{\mu\nu}R$ represents the
Einstein tensor. And the energy density center of the twisted particle in a
curved spacetime can be defined as

$\langle
x^{i}\rangle=\frac{\int^{+\infty}_{-\infty}\sqrt{g}x^{i}T^{00}dx^{3}}{\int^{+\infty}_{-\infty}\sqrt{g}T^{00}dx^{3}},$
(4)

where $g\equiv-\text{Det}(g_{\mu\nu})$ is the determinant of the metric
tensor.

Consider an arbitrary locally inertial coordinate system, denoted as
$\xi^{a}$, where the metric devolves to $g_{\mu\nu}\to\eta_{ab}$ and the
covariant derivative becomes to $\nabla_{\mu}\to\partial_{a}$. In the
coordinate system $\xi^{a}$, Eq. (2) transforms into
$\partial_{\mu}\partial^{\mu}\phi-\lambda R_{\xi}\phi=0$. Here, $R_{\xi}$
denotes the Ricci scalar in the $\xi^{a}$ coordinate system. Note, however,
that the Ricci scalar $R_{\xi}$ equals zero if and only if $R=0$ holds in the
global coordinates $x^{\mu}$. Consequently, if the nonminimal coupling
constant $\lambda\neq 0$, the EoM for a twisted particle propagating inside
gravitational sources, where $R\neq 0$, violates the EEP. This violation might
be observable when a twisted particle is emitted from the inside of the star.

Before delving into the specifics of this work, let us provide a brief
overview of the methodology employed to investigate the motion of the twisted
particle in a gravitational field. Initially, we configure the physical system
by establishing specific initial conditions and selecting an appropriate
gravitational background. Subsequently, we solve the EoM for the twisted
particle, as represented by Eq. (2), using the predetermined initial
conditions and gravitational background, a critical step that informs our
subsequent discussion on the dynamics of the twisted particle. Finally, we
determine the center of the energy density employing Eqs. (3) and (4), and
then investigate the motion of the twisted particle.

###  II.1 The Physical System and Initial Conditions

As mentioned in Sec. I, the nonminimal coupling term $\lambda R|\phi|^{2}$ is
activated if and only if the twisted particle is within a gravitational
source, prompting our study of particles emitted from within stars. As shown
in Fig. 1, a twisted particle with an OAM $l$ free-falls inside a star, shown
as a grey spherical region. Its initial position is located on the $x$-axis,
and its initial velocity is directed along the $z$-axis. Throughout this
paper, we set the radius of the star to $R_{\odot}=1$ and the initial position
of the twisted particle to $\vec{X}=(0.2,0,0)$.

$O$$\vec{v}$$R_{\odot}$$l$$X$$z$$x$$y$ Figure 1: Schematic illustration of the
physical system considered in this work. A twisted particle, shown as a blue
sphere and carrying an OAM $l$, free-falls inside a star. The star, shown as a
gray sphere, has a radius of $R_{\odot}$.

In a flat spacetime, the twisted particle possessing an OAM $l$, described by
the scalar wave packet, can be represented using the Fourier Transform as

$\bar{\phi}(\vec{x},t)=\frac{1}{(2\pi)^{3/2}}\int^{+\infty}_{-\infty}\exp(-\frac{(\vec{k}-\vec{k}_{0})^{2}}{2\Delta
k^{2}}+il\varphi)\exp(i\vec{k}\cdot\vec{x}-i\omega t)dk^{3},$ (5)

where $\omega=|\vec{k}|$ denotes the angular frequency and $\varphi$ is
azimuthal angle of $\vec{k}$. In Eq. (5), the initial average momentum of the
twisted particle ought to be $\vec{k}_{0}=(0,0,|\vec{k}_{0}|)$, considering
that its initial velocity $\vec{v}$ is directed along the $z$-axis.

For a specific wave packet, its wavelength $\lambda_{\phi}$ is typically far
smaller than its radius $d$. This characteristic implies that the magnitude of
its initial average momentum, denoted as $k_{0}\equiv|\vec{k}_{0}|$, is
significantly larger than the dispersion of its momentum $\Delta k$. The
presence of the Gaussian factor
$\exp\left(-\frac{(\vec{k}-\vec{k}_{0})^{2}}{2\Delta k^{2}}\right)$ ensures a
rapid decay of the function within the integral as
$|\vec{k}-\vec{k}_{0}|/\Delta k$ increases. This indicates that
$|\vec{k}-\vec{k}_{0}|\sim\Delta k$ is also substantially smaller than
$|\vec{k}_{0}|$. Retaining terms to the second order of
$|\vec{k}-\vec{k}_{0}|/k_{0}$, the angular frequency $\omega$ approximates to

$\omega\simeq k_{0}+\delta k_{\|}+\frac{\delta k^{2}_{\perp}}{k_{0}},$ (6)

where $\delta k_{\|}\equiv|\vec{k}\cdot\vec{k}_{0}-k_{0}^{2}|/k_{0}$ and
$\delta k_{\perp}\equiv|(\vec{k}-\vec{k}_{0})\times\vec{k}_{0}|/k_{0}$
represent the magnitudes of longitudinal and transverse components of
$\vec{k}-\vec{k}_{0}$ with respect to $\vec{k}_{0}$.

By substituting Eq. (6) into Eq. (5), the integration can be carried out,
yielding the scalar wave function $\bar{\phi}(\vec{x},t)$ as expressed by

$\begin{split}\bar{\phi}(\vec{x},t)&\simeq\frac{\pi^{2}r_{\perp}\Delta
k^{4}k_{0}^{\frac{3}{2}}}{(k_{0}+it\Delta
k^{2})^{\frac{3}{2}}}\exp\left(-\frac{r_{\perp}^{2}\Delta
k^{2}k_{0}^{2}}{4(k^{2}_{0}+t^{2}\Delta k^{4})}-\frac{(z-t)^{2}\Delta
k^{2}}{2}\right)\\\
&\times\left\\{I_{\frac{|l|-1}{2}}\left(\frac{r_{\perp}^{2}\Delta
k^{2}k_{0}}{4(k_{0}+it\Delta
k^{2})}\right)-I_{\frac{|l|+1}{2}}\left(\frac{r_{\perp}^{2}\Delta
k^{2}k_{0}}{4(k_{0}+it\Delta k^{2})}\right)\right\\}\\\
&\times\exp\left(i\vec{k}_{0}\cdot\vec{x}-ik_{0}t+il\varphi+it\frac{r_{\perp}^{2}\Delta
k^{4}k_{0}}{4(k^{2}_{0}+t\Delta k^{4})}\right),\\\ \end{split}$ (7)

where $r_{\perp}\equiv\sqrt{x^{2}+y^{2}}$ and $I_{a}(x)$ denotes the modified
Bessel function of the first kind. The wave function $\bar{\phi}(\vec{x},t)$
is an approximate solution to the massless Klein-Gordan equation in the flat
spacetime $\partial_{\mu}\partial^{\mu}\bar{\phi}(\vec{x},t)=0$ with the
following initial conditions:

$\begin{split}\bar{\phi}(\vec{x},t)|_{t=0}&=i^{|l|}\pi^{2}r_{\perp}\Delta
k^{4}\exp\left(-\frac{r^{2}_{\perp}\Delta k^{2}}{4}-\frac{z^{2}\Delta
k^{2}}{2}\right)\\\
&\times\left\\{I_{\frac{|l|-1}{2}}\left(\frac{r^{2}_{\perp}\Delta
k^{2}}{4}\right)-I_{\frac{|l|+1}{2}}\left(\frac{r^{2}_{\perp}\Delta
k^{2}}{4}\right)\right\\}\exp\left(i\vec{k}_{0}\cdot\vec{x}+il\varphi\right),\\\
\partial_{t}\bar{\phi}(\vec{x},t)|_{t=0}&=-\frac{i^{|l|+1}\pi^{2}\Delta
k^{4}r_{\perp}}{2\text{kz}}\left\\{c_{1}I_{\frac{|l|+3}{2}}\left(\frac{\Delta
k^{2}r^{2}_{\perp}}{4}\right)+\frac{c_{2}}{\Delta
k^{2}r^{2}_{\perp}}I_{\frac{|l|+1}{2}}\left(\frac{\Delta
k^{2}r^{2}_{\perp}}{4}\right)\right\\}\\\ &\times\exp\left(-\frac{\Delta
k^{2}\left(r^{2}_{\perp}+2z^{2}\right)}{4}+i\vec{k}_{0}\cdot\vec{x}+il\varphi\right),\\\
\end{split}$ (8)

where the two factors $c_{1}$ and $c_{2}$ are represented as

$\begin{split}c_{1}&=2k_{0}^{2}-\Delta
k^{4}r^{2}_{\perp}+(2ik_{0}z+|l|+2)\Delta k^{2},\\\
c_{2}&=8k_{0}^{2}(1+|l|)+\Delta k^{6}r^{4}_{\perp}-\Delta
k^{4}r^{2}_{\perp}(6+3|l|+2ik_{0}z)\\\ &+2\Delta
k^{2}\left\\{4+6|l|+2|l|^{2}-k_{0}^{2}r^{2}_{\perp}+4ik_{0}(1+|l|)z\right\\}.\\\
\end{split}$ (9)

The twisted particle, represented by the scalar wave function
$\bar{\phi}(\vec{x},t)$, is initially located at the origin of the coordinate
system $x^{\mu}$, which corresponds to the center of the star as shown in FIG.
1. However, we have set the twisted particle’s initial position to be located
on the $x$-axis, not at the center of the star. Consequently, we must adjust
the initial conditions via a translation along the $x$-axis. For convenience,
we define this position as $\vec{r}_{0}\equiv(|\vec{r}_{0}|,0,0)$. This
adjustment necessitates the replacement of $r_{\perp}$ and $\vec{x}$ with
$r^{\prime}_{\perp}\equiv\sqrt{(x-|\vec{r}_{0}|)^{2}+y^{2}}$ and
$\vec{x}^{\prime}\equiv\vec{x}-\vec{r}_{0}$, respectively. Upon this
adjustment, the aforementioned initial conditions should yield a wave function
$\phi$, as derived from solving the EoM given by Eq. (2). As suggested by Eq.
(7), this wave function $\phi$ represents a free-falling twisted particle with
an OAM of $l$ in a gravitational field.

###  II.2 The Gravitational Background

Despite excluding all non-gravitational couplings, it remains a considerable
challenge to accurately simulate the distribution of matter inside a star,
along with the corresponding gravitational field. However, our focus in this
study is not on the internal structure of the star, but on the motion of a
twisted particle. As such, we approximate the matter distribution within the
star using a polytrope. This simplification allows for a straightforward
calculation of the gravitational field, as demonstrated in the subsequent
subsection.

A polytrope is defined by its equation of state satisfying the given formula:

$p(r)=K\rho(r)^{\gamma},$ (10)

where $p(r)$ and $\rho(r)$ are the pressure and rest-mass density, $K$ is a
constant number for a specific star, and $r$ denotes the radial distance. The
process of solving the matter distribution of a polytrope and its
corresponding internal gravitational field is a standard topic extensively
covered in many textbooks on GR and stellar astrophysics. In the following
discussion, we merely provide a brief overview. For a more detailed
discussion, please refer to Chapter 11 in Ref. Weinberg (1972).

We introduce two variables, $\xi$ and $\theta(\xi)$, defined as follows:

$r\equiv\left(\frac{K\gamma}{4\pi
G(\gamma-1)}\right)^{\frac{1}{2}}\rho(0)^{(\gamma-2)/2}\xi,\quad\rho\equiv\rho(0)\theta^{1/(\gamma-1)}(\xi),$
(11)

where $\rho(0)$ is the rest-mass density at the center of this star.
Consequently, the rest-mass density $\rho(r)$ can be determined by solving the
Lane–Emden equation:

$\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta(\xi)}{d\xi}\right)+\theta^{1/(\gamma-1)}(\xi)=0$
(12)

with the following boundary conditions:

$\theta(0)=0\qquad\theta^{\prime}(0)=0,$ (13)

where $\theta^{\prime}\equiv\frac{d\theta(\xi)}{d\xi}$. The radius of this
star, denoted as $R_{\odot}$, can be obtained as:

$R_{\odot}=\left(\frac{K\gamma}{4\pi
G(\gamma-1)}\right)^{\frac{1}{2}}\rho(0)^{(\gamma-2)/2}\xi_{1},$ (14)

where $\xi_{1}$ is the corresponding value of $\xi$ when $\rho=0$.

In this study, we adopt an ad hoc value of $\gamma=4/3$. Upon solving Eq. (12)
and considering the definition of $\theta(\xi)$ from Eq. (11) , we can
determine the rest-mass density $\rho(r)$ and pressure $p(r)$ inside this
star. As shown in Fig. 2, the ratio of pressure $p(r)$ to rest-mass density
$\rho(r)$ is approximately $p(r)/\rho(r)\sim 10^{-6}\ll 1$. This indicates
that the pressure $p(r)$ is greatly overshadowed by the rest-mass density
$\rho(r)$ inside this star. Consequently, the pressure $p(r)$ can be neglected
in the metric $g_{\mu\nu}$.

For a specific value of $\gamma$, the internal metric $g_{\mu\nu}$ of a star
is determined by its surface gravitational potential, denoted as
$\varepsilon\equiv\frac{GM_{\odot}}{R_{\odot}}$. Here, $M_{\odot}$ represents
the total mass of the star. In this study, we concentrate on the case where
$\varepsilon\ll 1$, indicating a weak gravitational field. This approximation
is applicable to main-sequence stars, such as the Sun, for which
$\varepsilon\sim 10^{-6}\ll 1$. By retaining only the first order of
$\varepsilon$, the internal metric simplifies to

$g_{\mu\nu}\simeq\eta_{\mu\nu}+h_{\mu\nu},\quad
h_{00}=2\varepsilon+2\varepsilon f(r),\quad
h_{ij}=2\varepsilon\frac{g(r)x_{i}x_{j}}{r^{2}},\quad h_{0i}=0$ (15)

with the following definitions:

$M(r)\equiv\int^{r}_{0}4\pi r^{\prime
2}\rho(r^{\prime})dr^{\prime},\quad\varepsilon
g(r)\equiv\frac{GM(r)}{r},\quad\varepsilon
f(r)\equiv\int^{R_{\odot}}_{r}\frac{g(r^{\prime})}{r^{\prime}}dr^{\prime},$
(16)

where $\eta_{\mu\nu}=\text{diag}(-1,+1,+1,+1)$ is the Minkowski metric.

Figure 2: The rest-mass density and pressure of a polytrope when $\gamma=4/3$.
Here the radius of this star has been set to $R_{\odot}=1$ and $\bar{\rho}$ is
the average mass density.

###  II.3 Numerical Method

As demonstrated in Eqs. (15), the internal gravitational field of the star is
weak, with the gravitational component in the metric
$h_{\mu\nu}\propto\varepsilon\ll 1$. Retaining only the first order of
$\varepsilon$, the corresponding affine connection
$\Gamma^{\lambda}_{\rho\alpha}$ and Ricci scalar $R$ can be approximated as

$\Gamma^{\lambda}_{\rho\alpha}\simeq\frac{1}{2}\eta^{\lambda\gamma}(\partial_{\rho}h_{\gamma\alpha}+\partial_{\alpha}h_{\gamma\rho}-\partial_{\gamma}h_{\rho\alpha}),\qquad
R\simeq-\partial_{\lambda}\partial_{\alpha}h^{\lambda\alpha},$ (17)

where
$h^{\lambda\alpha}\equiv\eta^{\lambda\rho}\eta^{\alpha\beta}h_{\rho\beta}$.
Therefore, Eq. (2) simplifies to

$\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}\phi=h^{\mu\nu}\partial_{\mu}\partial_{\nu}\phi+\partial_{\alpha}h^{\alpha\beta}\partial_{\beta}\phi-\frac{1}{2}\eta^{\alpha\beta}\partial_{\alpha}h\partial_{\beta}\phi-\lambda\partial_{\alpha}\partial_{\beta}h^{\alpha\beta}\phi,$
(18)

where $h\equiv\eta^{\mu\nu}h_{\mu\nu}$. In accordance with Eqs. (15), the EoM
for the twisted particle inside the star can be expressed as

$\eta^{\mu\nu}\partial_{\mu}\partial_{\nu}\phi=2\varepsilon
f(r)\eta^{ij}\partial_{i}\partial_{j}\phi+2\varepsilon\frac{g(r)x^{i}x^{j}}{r^{2}}\partial_{i}\partial_{j}\phi+3\varepsilon\rho
x^{i}\partial_{i}\phi+2\varepsilon
g(r)\frac{x^{i}}{r^{2}}\partial_{i}\phi-6\lambda\varepsilon\rho^{\prime}(r)\phi,$
(19)

where $\varepsilon\rho^{\prime}(r)\equiv\rho(r)/\bar{\rho}$.

It is nearly impossible to find an analytical solution for the EoM given by
Eq. (19) with the initial conditions as represented by Eq. (8). Consequently,
we resort to numerical methods to solve the EoM in this work. It is worth
noting that the wave function $\phi$ is expected to exhibit rapid oscillations
in both space and time due to the term
$\exp(i\vec{k}_{0}\cdot\vec{x}-ik_{0}t)$, as suggested by Eq. (7). This
implies that the time and space step sizes in the numerical calculations must
be significantly smaller than the wavelength $\lambda_{\phi}=2\pi/k_{0}$.
However, reducing the time and space step sizes increases the computational
demand. Therefore, we must simplify Eq. (19) for efficient numerical
computations.

Let us introduce a new wave function $\phi^{\prime}$, defined as

$\phi\equiv\phi^{\prime}\exp(i\vec{k}_{0}\cdot\vec{x}-ik_{0}t).$ (20)

The EoM for $\phi^{\prime}$ can be derived from Eq. (19), which yields

$\begin{split}&\quad-\partial_{t}^{2}\phi^{\prime}+2ik_{0}\partial_{t}\phi^{\prime}+\eta^{ij}\partial_{i}\partial_{j}\phi^{\prime}+2ik_{0}\delta_{3}^{i}\partial_{i}\phi^{\prime}\\\
&=2\varepsilon\left(f(r)\eta^{ij}+\frac{g(r)x^{i}x^{j}}{r^{2}}\right)\left(\partial_{i}\partial_{j}\phi^{\prime}+ik_{0}\delta^{3}_{i}\partial_{j}\phi^{\prime}+ik_{0}\delta^{3}_{j}\partial_{i}\phi^{\prime}-k^{2}_{0}\delta^{3}_{i}\delta^{3}_{j}\phi^{\prime}\right)\\\
&+\varepsilon\left(3\rho+\frac{2g(r)}{r^{2}}\right)x^{i}(\partial_{i}\phi^{\prime}+ik_{0}\delta^{3}_{i}\phi^{\prime})-6\lambda\varepsilon\rho^{\prime}(r)\phi^{\prime}.\\\
\end{split}$ (21)

As suggested by Eq. (7), the rapid oscillations in the wave function
$\phi^{\prime}$ are expected to vanish. Consequently, the spatial step size
can significantly exceed the wavelength $\lambda_{\phi}$.

Significantly, Eq. (21) remains a second-order inhomogeneous partial
differential equation in four-dimensional spacetime, the complexity of which
poses substantial challenges to even numerical calculations. To simplify, we
apply a Fourier Transform to this equation in three-dimensional space, thereby
transforming Eq. (21) into an ordinary differential equation. Consequently,
the EoM of $\phi^{\prime}$ simplifies to

$\partial^{2}_{t}\phi^{\prime}_{f}-2ik_{0}\partial_{t}\phi^{\prime}_{f}-k_{0}^{2}\phi^{\prime}_{f}=-k^{\prime}_{i}k^{\prime
i}\phi^{\prime}_{f}+2\varepsilon\mathcal{F}^{ij}*(k^{\prime}_{i}k^{\prime}_{j}\phi^{\prime}_{f})-i\varepsilon\mathcal{G}^{i}*(k^{\prime}_{i}\phi_{f})-6\lambda\varepsilon\rho^{\prime}_{f}*\phi^{\prime}_{f},$
(22)

given the following definitions:

$\begin{split}&\vec{k}^{\prime}\equiv(k^{1},k^{2},k^{3}+k_{0}),\quad\mathcal{F}^{ij}\equiv\mathscr{F}\left(f(r)\eta^{ij}+\frac{g(r)x^{i}x^{j}}{r^{2}}\right),\quad\mathcal{G}^{i}\equiv\mathscr{F}\left(3\rho
x^{i}+\frac{2g(r)}{r^{2}}x^{i}\right)\\\
&\rho^{\prime}_{f}\equiv\mathscr{F}(\rho^{\prime}),\quad\phi^{\prime}_{f}\equiv\mathscr{F}(\phi^{\prime}),\\\
\end{split}$ (23)

where
$\mathscr{F}(\phi)\equiv\frac{1}{(2\pi)^{3/2}}\int^{+\infty}_{-\infty}\phi(x)\exp(-i\vec{k}\cdot\vec{x})dx^{3}$
is the Fourier Transform in three-dimensional space and the operator ’$*$’
signifies the convolution in three-dimensional space
$(f*g)(\vec{x})\equiv\int^{+\infty}_{-\infty}f(\vec{x}^{\prime})g(\vec{x}-\vec{x}^{\prime})dx^{\prime
3}$.

There exists a plethora of numerical methods for solving this ordinary
differential equation (22). In this study, we employ the midpoint method to
solve Eq. (22) with the initial conditions previously mentioned. Once
$\phi^{\prime}_{f}$ is obtained using this method, the wave function $\phi$
can be effortlessly calculated as

$\phi(\vec{x},t)=\exp(i\vec{k}_{0}\cdot\vec{x}-ik_{0}t)\mathscr{F}^{-1}(\phi^{\prime}_{f}(\vec{k},t)),$
(24)

where $\mathscr{F}^{-1}$ denotes the inverse Fourier Transform of
$\mathscr{F}$.

In this paper, unless otherwise specified, the numerical calculations are
performed with the following default parameters: the wavelength
$\lambda_{\phi}=\pi/2500$, the OAM $l=+1$, and the gravitational potential at
the surface of the star $\varepsilon=10^{-3}$.

##  III Gravitational Birefringence Induced By OAM

Numerous researchers have pointed out that the gravitational birefringence
occurs in the direction perpendicular to both the gravity and angular momentum
Duval  _et al._ (2019); Oancea  _et al._ (2020); Lian and Chen (2022); Wang .
Given the geometric configuration of the physical system shown in Fig. 1, this
phenomenon manifests along the direction parallel to the $y$-axis. Therefore,
this section primarily focuses on the trajectory along the $y$-axis under the
minimal coupling condition $\lambda=0$.

The trajectory along the $y$-axis, represented by the the temporal evolution
of $\langle y\rangle$ is shown in Fig. 3. The analysis reveals an approximate
linear dependence of $\langle y\rangle$ on the time $t$ for $t<0.1$. It is
also indicated that $\langle y\rangle$ is proportional to the wavelentgh
$\lambda_{\phi}$, surface gravitational potential of the star $\varepsilon$
and OAM $l$, respectively. Therefore, the trajectory along the $y$-axis can be
approximately given by the following relation:

$\langle y\rangle\simeq 2\varepsilon l\lambda_{\phi}t.$ (25)

Figure 3: The trajectories along the $y$-axis.

During the movement of the twisted particle from inside to outside of the
star, the wavelength $\lambda_{\phi}$, OAM $l$ and surface gravitational
potential of the star $\varepsilon$, are expected to remain constant.
Consequently, the proportionality $\langle
y\rangle\propto\varepsilon\lambda_{\phi}l$ is anticipated to be preserved.
However, due to significant spatial variations in the gravitational field
during the movement, the trajectory along the $y$-axis will inevitably deviate
from the simple linear function of time $t$. The numerical results are shown
in Fig. 4, which indicate that upon reaching the stellar surface, the velocity
along the $y$-axis of the twisted particle is approximately given by

$\frac{d\langle y\rangle}{dt}\arrowvert_{surf}\simeq 1.09\times 10^{-7}\simeq
0.08\varepsilon l\lambda_{\phi},$ (26)

where the parameters $\lambda_{\phi}=\pi/2500$, $\varepsilon=10^{-3}$, $l=\pm
1$ have been considered, along with the previously established proportionality
$\langle y\rangle\propto\varepsilon l\lambda_{\phi}$. Consequently, the angle
$\theta\simeq\frac{2d\langle y\rangle}{dt}$, representing the angular
separation along the $y$-axis between two twisted particles with opposite OAM
$l=\pm|l|$, is approximated as

$\theta\sim 0.16\times\varepsilon|l|\lambda_{\phi}.$ (27)

It is worth noting that the angle $\theta$ is proportional to both the
wavelength $\lambda_{\phi}$ and OAM $l$. Therefore, an increase in OAM leads
to a more pronounced gravitational birefringence. Under controlled laboratory
conditions, it has been demonstrated that light can carry an OAM as high as
$10^{6}\hbar$ per photon Zhang  _et al._ (2021). As such, the gravitational
birefringence induced by OAM could significantly exceed that induced by spin,
potentially facilitating its detection. Similarly, an increase in wavelength
also enhances the effect of gravitational birefringence. Therefore, detection
using radio waves might prove easier than using light or gamma rays.

Figure 4: The trajectories of two twisted particles with opposite OAM $l=\pm
1$ as they propagate from inside to outside of the star. The dot-dashed line
representes the surface of the star.

##  IV Nonminimal Coupling Effects

According to Eqs. (2) and (3), variations in the nonminimal coupling constant
$\lambda$ induce changes in both the wave function $\phi$ and EMT
$T^{\mu\nu}$, irrespective of identical initial conditions. As a result, the
energy density center of the twisted particle is expected to shift as
$\lambda$ varies. Consequently, The trajectory of the twisted particle,
defined by its energy density center, can alter with different values of
$\lambda$, as the twisted particle free falls inside a gravitational source.
As mentioned before, the EoM with the nonminimal coupling term $\lambda R\phi$
presents a challenge to the EEP, suggesting the possibility of a novel
phenomenon that could serve as a test for the EEP.

After numerical calculations, the trajectory of the twisted particle with
different nonminimal coupling constants $\lambda$ is shown in Fig. 5. It is
indicated that the nonminimal couplings can significantly affect the
trajectory along the $x$-axis, while leaving the trajectory along the $y$-axis
negligibly affected. Therefore, in this section, we mainly focus on the
trajectory along the $x$-axis which differs significantly with different
$\lambda$. In the remainder of this section, unless otherwise specified, the
numerical calculations are performed with the default nonminimal coupling
constant $\lambda=10^{4}$.

Figure 5: The trajectory of the free-falling twisted particle with different
nonminimal coupling constant $\lambda$. In the left panel, the initial
position of the twisted particle is aligned with the origin of the $x$-axis.

The numerical results, presented in Fig. 6, indicate that the deviation of the
trajectory along the $x$-axis, denoted as $\Delta x\equiv\langle
x\rangle_{\lambda}-\langle x\rangle_{\lambda=0}$, between the scenarios of
nonminimal and minimal couplings is approximately proportional to the
nonminimal coupling constant $\lambda$, the surface gravitational potential of
the star $\varepsilon$, and the square of the wavelength $\lambda_{\phi}^{2}$,
respectively. In contrast, the OAM $l$ appears to have a negligible influence
on the trajectory deviation $\Delta x$. Therefore, the relationship between
the deviation $\Delta x$ and the aforementioned parameters can be
approximately represented as

$\Delta x\propto\varepsilon\lambda\lambda_{\phi}^{2}.$ (28)

Figure 6: The deviation of the trajectory along the $x$-axis between the
nonminimal and minimal coupling cases.

The trajectory of the twisted particle along the $x$-axis during its movement
from inside to outside of the star is depicted in Fig. 7. Its velocity along
the $x$-axis, $v_{x}$, can be defined by $v_{x}\equiv\frac{d\langle
x\rangle}{dt}$. Correspondingly, the difference in $v_{x}$ between the
nonminimal and minimal coupling scenarios is represented by

$\Delta v_{x}\equiv-\left(\frac{d\langle
x\rangle_{\lambda}}{dt}-\frac{d\langle
x\rangle_{\lambda=0}}{dt}\right)=-\frac{d\Delta x}{dt},$ (29)

where the negative sign indicates that, with a positive nonminimal coupling
constant $\lambda$, the trajectory deviates along the direction away from the
stellar center when compared to the trajectory in the minimal coupling case.

Figure 7: The trajectory along the $x$-axis and its deviation as the twisted
particle propagates from inside to outside of the star within the scenarios of
nonminimal and minimal couplings. The dot-dashed line representes the surface
of the star.

Numerical results shown in Fig. 7 indicate that, upon the twisted particle
reaching the stellar surface, $\Delta v_{x}$ is approximately given by

$\Delta v_{x}|_{surf}\simeq-1.08\times 10^{-7}.$ (30)

As the twisted particle moves from inside to outside of the star, the
proportionality $\Delta x\propto\varepsilon\lambda\lambda_{\phi}^{2}$ is also
anticipated to be preserved, mirroring the reasons underlying the previously
established $\langle y\rangle\propto\varepsilon\lambda_{\phi}l$.
Correspondingly, the difference in $v_{x}$ at the stellar surface, denoted as
$\Delta v_{x}|_{surf}$, follows the established proportionality $\Delta
v_{x}|_{surf}\propto-\varepsilon\lambda\lambda_{\phi}^{2}$. Given the
parameter settings $\lambda=10^{3}$, $\varepsilon=10^{-3}$ and
$\lambda_{\phi}=\pi/2500$ in the numerical calculations to obtain Eq. (30),
$\Delta v_{x}|_{surf}$ can be approximately represented as

$\Delta v_{x}|_{surf}\simeq-0.07\times\varepsilon\lambda\lambda_{\phi}^{2}.$
(31)

Outside the star, the Ricci scalar $R$ becomes zero, leading to the vanishing
of the nonminimal coupling term $\lambda R|\phi|^{2}$. Consequently, the
velocity difference $\Delta v_{x}$ for the twisted particle becomes constant
upon exiting the stellar surface. Therefore, Eq. (31) suggests that the
trajectory in the case of nonminimal coupling deviates from the minimal
coupling case by an angle $\theta_{non}$ along the $x$-axis, which can be
approximated as

$\theta_{non}\sim-0.07\times\varepsilon\lambda\lambda_{\phi}^{2}.$ (32)

It has been established that the nonminimal coupling constant $\lambda$ should
be set to $1/6$ to ensure the finiteness of the EMT for scalar fields at one-
loop order Freedman  _et al._ (1974); Freedman and Weinberg (1974). In this
case, the order of magnitude for the angle $\theta_{non}$ is given by

$\theta_{non}\sim-10^{-2}\times\frac{\varepsilon\lambda_{\phi}^{2}}{R_{\odot}^{2}},$
(33)

where the radius of the star $R_{\odot}$ explicitly appears to maintain the
generality of the formula. As indicated by Eq. (33), the angle $\theta_{non}$
scales with $\lambda_{\phi}^{2}/R_{\odot}^{2}$. Given that the wavelength
$\lambda_{\phi}$ of a realistic wave packet is typically at least several
orders of magnitude smaller than the stellar radius $R_{\odot}$, the angle
$\theta_{non}$ induced by the nonminimal couplings is expected to be
drastically suppressed.

##  V Discussion

As previously established, the trajectories of two free-falling twisted
particles carrying opposite OAM exhibit a transverse separation perpendicular
to both the OAM and gravity, which is a manifestation of gravitational
birefringence. This separation adheres to the proportionality
$\theta\propto\varepsilon l\lambda_{\phi}/R_{\odot}$, which is in agreement
with the well-documented relationship between the separation and spin angular
momentum. Such agreement suggests that the gravitational birefringence induced
by OAM is expected to be analogous to that induced by spin angular momentum.

For a twisted particle, the OAM can significantly surpass its intrinsic spin,
potentially by orders of magnitude. This substantial difference suggests that
the gravitational birefringence induced by OAM might be much larger than that
induced by spin. To illustrate, consider the birefringence of sunlight at a
typical wavelength of $\lambda_{\phi}\sim 500\text{nm}$ as it passes through
the gravitational field of the Sun. The spin of light is limited to $\hbar$
per photon, while its OAM can exceed $10^{4}\hbar$ per photon, as reported in
Ref. Fickler  _et al._ (2016). The surface gravitational potential of the Sun
is approximated to be $\varepsilon\sim 2.43\times 10^{-6}$, with a radius of
$R_{\odot}\sim 6.96\times 10^{5}\text{km}$. Consequently, the separation
induced by the spin of light is estimated to be $\theta\sim
10^{-22}\text{rad}$, whereas the separation due to its OAM is expected to be
$\theta\sim 10^{-18}\text{rad}$, which is an increase by four orders of
magnitude over that induced by the spin. This substantial amplification
associated with OAM could significantly enhance the detectability of
gravitational birefringence, thereby improving the prospects for experimental
verification.

Unlike the MPD equations, the center of energy density approach suggests that
a free-falling twisted particle is expected to follow different trajectories
with variations in the nonminimal coupling constant $\lambda$. Specifically,
as the twisted particle propagates from inside to outside of a star along a
non-radial direction, its trajectory in the case of nonminimal coupling
deviates from the minimal coupling case by an angle
$\theta_{non}\sim-0.07\times\varepsilon\lambda\lambda_{\phi}^{2}/R_{\odot}^{2}$
along the direction of gravity. The negative sign here signifies that, for a
positive coupling constant $\lambda$, the deviation is directed away from the
stellar center. It is important to note that the angle $\theta_{non}$ does not
depend on the OAM, suggesting that the effect of nonminimal coupling is not
amplified by the OAM.

Due to the proportionality of $\lambda_{\phi}^{2}/R_{\odot}^{2}$, this angle
$\theta_{non}$ is expected to be many orders of magnitude smaller than the
deflection angle $\Delta\varphi$ predicted by GR, which can be approximated by
$\Delta\varphi\sim\varepsilon$. For a twisted particle emitted from the inside
of the Sun with a typical wavelength $\lambda_{\phi}\sim 500\text{nm}$, the
ratio of its wavelength to the radius of the Sun is approximated as
$\lambda_{\phi}/R_{\odot}\sim 7.18\times 10^{-16}$. Given that $1/6$ is
identified as the only reasonable value for the nonminimal coupling constant
$\lambda$, as indicated by Refs. Freedman  _et al._ (1974); Freedman and
Weinberg (1974), the estimated angle $\theta_{\text{non}}$ is approximated as
$\theta_{\text{non}}\sim-10^{-38}\,\text{rad}$. This value is $32$ orders of
magnitude smaller than the gravitational deflection angle $\Delta\varphi$,
which is on the order of $10^{-6}\text{rad}$. With current technology,
detecting an angle as small as $\theta_{non}$ is almost impossible.

It is important to note that both the angles, $\theta$ and
$\theta_{\text{non}}$, are proportional to the gravitational potential
$\varepsilon$. The angle $\theta$ is inversely proportional to the radius of
the gravitational source $R_{\odot}$, while the angle $\theta_{non}$ is
inversely proportional to its square, $R_{\odot}^{2}$. Therefore, with an
increase in the gravitational potential $\varepsilon$ and a decrease in the
radius of its source, the angles $\theta$ and $\theta_{\text{non}}$ are
expected to increase significantly. As an example, replacing the Sun with a
pulsar having a gravitational potential $\varepsilon\sim 0.1$ at a radius of
$R_{\odot}\sim 10$ km, the angle $\theta$ with the OAM $l=10^{4}$ is expected
to become $\theta\sim 10^{-8}\text{rad}$, which is nearly $10$ orders of
magnitude greater than that in the gravitational field of the Sun. Meanwhile,
the angle $\theta_{non}$ is expected to be
$\theta_{non}\sim-10^{-24}\text{rad}$, nearly $14$ orders of magnitude larger
than that in the gravitational field of the Sun. These results indicate that,
in the gravitational fields of compact objects such as black holes and neutron
stars, the angles $\theta$ and $\theta_{non}$ could potentially be large
enough to be detectable. Therefore, the gravitational birefringence and
nonminimal coupling effects of twisted particles could provide new
opportunities to test the EEP, despite the significant observational
challenges.

###### Acknowledgements.

We are grateful to Xiang-Song Chen for the helpful discussions, and we thank
Zheng-Cheng Liang for his suggestions during the numerical calculations. This
work is supported by the National Natural Science Foundation of China(Grant
No. 12375084).

## References

  * Mathisson (2010) M. Mathisson, Gen. Relativ. Gravit. 42, 1011 (2010).
  * Papapetrou (1951) A. Papapetrou, Proc. Roy. Soc. Lond. A 209, 248 (1951).
  * Dixon (1970a) W. G. Dixon, Proc. Roy. Soc. Lond. A 314, 499 (1970a).
  * Dixon (1970b) W. G. Dixon, Proc. Roy. Soc. Lond. A 319, 509 (1970b).
  * Puetzfeld  _et al._ (2015) D. Puetzfeld, C. Lämmerzahl, and B. Schutz,  _Equations of Motion in Relativistic Gravity_ (Springer International Publishing, Cham, 2015).
  * Corinaldesi and Papapetrou (1951) E. Corinaldesi and A. Papapetrou, Proc. Roy. Soc. Lond. A 209, 259 (1951).
  * Hartl (2003) M. D. Hartl, Phys. Rev. D 67, 024005 (2003).
  * Hackmann  _et al._ (2014) E. Hackmann, C. Lämmerzahl, Y. N. Obukhov, _et al._ , Phys. Rev. D 90, 064035 (2014).
  * Duval and Schücker (2017) C. Duval and T. Schücker, Phys. Rev. D 96, 043517 (2017).
  * Duval  _et al._ (2019) C. Duval, L. Marsot, and T. Schücker, Phys. Rev. D 99, 124037 (2019).
  * Frolov and Shoom (2012) V. P. Frolov and A. A. Shoom, Phys. Rev. D 86, 024010 (2012).
  * Oancea  _et al._ (2020) M. A. Oancea, J. Joudioux, I. Y. Dodin, _et al._ , Phys. Rev. D 102, 024075 (2020).
  * Deriglazov (2021) A. A. Deriglazov, Phys. Rev. D 104, 025006 (2021).
  * Harte and Oancea (2022) A. I. Harte and M. A. Oancea, Phys. Rev. D 105, 104061 (2022).
  * Andersson and Oancea (2023) L. Andersson and M. A. Oancea, Class. Quantum Grav. 40, 154002 (2023).
  * Lian and Chen (2022) D.-D. Lian and X.-S. Chen, Phys. Rev. D 105, 104008 (2022).
  * (17) Z.-L. Wang, arXiv:2310.17581 [gr-qc] .
  * Ivanov (2022) I. P. Ivanov, Prog. Part. Nucl. Phys. 127, 103987 (2022).
  * Fickler  _et al._ (2016) R. Fickler, G. Campbell, B. Buchler, _et al._ , Proc. Natl. Acad. Sci. U.S.A. 113, 13642 (2016).
  * Zhang  _et al._ (2021) H. Zhang, J. Zhao, Y. Hu, _et al._ , High Power Laser Sci. Eng. 9, e43 (2021).
  * Tamburini  _et al._ (2011) F. Tamburini, B. Thide, G. Molina-Terriza, _et al._ , Nature Phys. 7, 195 (2011).
  * Yang and Casals (2014) H. Yang and M. Casals, Phys. Rev. D 90, 023014 (2014).
  * Tamburini  _et al._ (2020) F. Tamburini, B. Thidé, and M. Della Valle, Mon. Not. Roy. Astron. Soc. 492, L22 (2020).
  * Tamburini  _et al._ (2021) F. Tamburini, F. Feleppa, I. Licata, _et al._ , Phys. Rev. A 104, 013718 (2021).
  * Freedman  _et al._ (1974) D. Z. Freedman, I. J. Muzinich, and E. J. Weinberg, Annals Phys. 87, 95 (1974).
  * Freedman and Weinberg (1974) D. Z. Freedman and E. J. Weinberg, Annals Phys. 87, 354 (1974).
  * Linde (1982) A. D. Linde, Phys. Lett. B 114, 431 (1982).
  * Faraoni (1996) V. Faraoni, Phys. Rev. D 53, 6813 (1996).
  * Aiello  _et al._ (2009) A. Aiello, N. Lindlein, C. Marquardt, _et al._ , Phys. Rev. Lett. 103, 100401 (2009).
  * Wang and Chen (2022) Z.-L. Wang and X.-S. Chen, Commun. Theor. Phys. 74, 105401 (2022).
  * Weinberg (1972) S. Weinberg,  _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_ (John Wiley and Sons, New York, 1972).