# Scale of Dirac leptogenesis and left-right symmetry in the light of recent
PTA results

Basabendu Barman<EMAIL_ADDRESS>Institute of Theoretical
Physics, Faculty of Physics, University of Warsaw,
ul. Pasteura 5, 02-093 Warsaw, Poland  Debasish Borah<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology Guwahati, Assam 781039,
India  Suruj Jyoti Das<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Guwahati, Assam 781039, India  Indrajit Saha
<EMAIL_ADDRESS>Department of Physics, Indian Institute of Technology
Guwahati, Assam 781039, India

###### Abstract

Motivated by the recent release of new results from five different pulsar
timing array (PTA) experiments claiming to have found compelling evidence for
primordial gravitational waves (GW) at nano-Hz frequencies, we study the
consequences for two popular beyond the Standard Model (SM) framework, where
such nano-Hz GW can arise due to annihilating domain walls (DW). Minimal
framework of Dirac leptogenesis, as well as left-right symmetric model (LRSM)
can lead to formation of DW due to spontaneous breaking of $Z_{2}$ symmetry.
Considering the NANOGrav 15 yr data, we show that the scale of Dirac
leptogenesis should be above $10^{7}$ GeV for conservative choices of Dirac
Yukawa couplings with fine-tuning at the level of the SM. The scale of minimal
LRSM is found to be more constrained $M_{\rm LR}\sim 10^{6}$ GeV in order to
fit the NANOGrav 15 yr data.

Introduction: Recently, four different pulsar timing array (PTA) experiments
namely NANOGrav Agazie et al. (2023), European Pulsar Timing Array (EPTA)
together with the first data release from Indian Pulsar Timing Array (InPTA)
Antoniadis et al. (2023), PPTA Reardon et al. (2023), all part of the
consortium called International Pulsar Timing Array (IPTA) have released their
latest findings hinting at a significant evidence for stochastic gravitational
waves (GW) background at nano-Hz frequencies. Similar evidence with larger
statistical significance has also been reported by the Chinese Pulsar Timing
Array (CPTA) collaboration Xu et al. (2023). While such a signal can be
generated by supermassive black hole binary (SMBHB) mergers though with a mild
tension, presence of exotic new physics alone or together with SMBHB can make
the fit better Afzal et al. (2023). Several follow-up papers have also studied
the possible origin or implications of this observation from the point of view
of dark matter Ghoshal and Strumia (2023); Shen et al. (2023), axions or
axion-like particles Yang et al. (2023); Guo et al. (2023), SMBHB Ellis et al.
(2023a), first order phase transition Megias et al. (2023); Fujikura et al.
(2023); Han et al. (2023); Zu et al. (2023), primordial black holes
Franciolini et al. (2023a), primordial magnetic field Li et al. (2023), domain
walls Kitajima et al. (2023); Bai et al. (2023), inflation Vagnozzi (2023),
cosmic strings Ellis et al. (2023b); Wang et al. (2023), astrophysical
neutrino oscillation Lambiase et al. (2023) and QCD crossover Franciolini et
al. (2023b).

While GW from domain walls (DW) has already been studied as a possible new
physics explanation for PTA results Afzal et al. (2023); Kitajima et al.
(2023); Bai et al. (2023), we consider the consequence for two popular beyond
standard model (BSM) scenarios namely, the minimal Dirac leptogenesis and the
left-right symmetric model (LRSM). The first model is a type I seesaw
realisation for light Dirac neutrino mass with the heavy vector-like neutral
fermions being responsible for generating baryogenesis via leptogenesis
Fukugita and Yanagida (1986) with light Dirac neutrinos, known as the Dirac
leptogenesis scenario Dick et al. (2000); Murayama and Pierce (2002). GW probe
of high scale leptogenesis models have received considerable attention in
recent times. In most of these works Dror et al. (2020); Blasi et al. (2020);
Fornal and Shams Es Haghi (2020); Samanta and Datta (2021), cosmic string (CS)
origin of GW has been studied by considering a $U(1)_{B-L}$ framework with in-
built heavy Majorana fermions responsible for generating Majorana mass of
light neutrinos as well as leptogenesis. The scale of leptogenesis or
$U(1)_{B-L}$ breaking scale then decides the amplitude of the CS generated GW
spectrum. However, in view of the latest PTA results preferring a positive
slope of the GW spectrum, stable CS in such models no longer provide a good
fit Afzal et al. (2023). This raises the prospects for a Dirac leptogenesis
model whose minimal version must have a softly broken $Z_{2}$ symmetry leading
to formation of DW followed by generation of GW due to annihilation or
collapse. While a general study related to GW probe of minimal Dirac
leptogenesis was carried out in Barman et al. (2022), here we consider the
implications of recent PTA findings on the scale of Dirac leptogenesis. On the
other hand, GW probe of LRSM considering DW as the source have been studied in
earlier works Craig et al. (2021); Borah and Dasgupta (2022). DW arise due to
spontaneous breaking of parity in such models. While earlier works considered
the detection aspects of this model, we now constrain the scale of left-right
symmetry considering the latest PTA data. While both the models can explain
the latest PTA data, the allowed parameter space remains squeezed to a tiny
window, which should face more scrutiny with future data.

Domain walls as source of GW: Domain wall is a two-dimensional topological
defect arising from spontaneous breaking of discrete symmetries Zeldovich et
al. (1974); Kibble (1976); Vilenkin (1981). With the expansion of the
universe, the energy density of DW falls slower compared to that of radiation
or ordinary matter, having the potential to start dominating the energy
density of the universe and ruin the successful predictions of standard
cosmology. Such a disastrous situation can be prevented if DW are made
unstable or diluted or if they have asymmetric initial field fluctuations
Coulson et al. (1996); Krajewski et al. (2021).

In minimal model of Dirac leptogenesis Barman et al. (2022) as well as left-
right symmetric models Borah and Dasgupta (2022), such DW arises due to the
spontaneous breaking of a $Z_{2}$ symmetry. If we consider a $Z_{2}$-symmetric
potential of a scalar field $\varphi$, it is straightforward to show the
existence of two different vacua $\langle\varphi\rangle=\pm u$. It is also
possible to find a static solution of the equation of motion given the two
vacua to be realized at $x\to\pm\infty$,

$\displaystyle\varphi({\bf
x})=u\,\tanh\left(\sqrt{\frac{\lambda_{\varphi}}{2}}\,u\,x\right)\,,$ (1)

representing a domain wall extended along the $x=0$ plane. Here
$\lambda_{\varphi}$ is the quartic self-coupling of the scalar field. The DW
width is $\delta\sim m_{\varphi}^{-1}=(\sqrt{2\lambda_{\varphi}}\,u)^{-1}$.
Another key parameter, known as the DW tension is given by

$\displaystyle\sigma_{w}=\int_{-\infty}^{\infty}dx\,\rho_{\varphi}=\frac{2\sqrt{2}}{3}\,\sqrt{\lambda_{\varphi}}\,u^{3}=\frac{2}{3}\,m_{\varphi}\,u^{2}\sim
u^{3}\,,$ (2)

where $\rho_{\varphi}$ denotes (static) energy density of $\varphi$ and in the
last step, $m_{\varphi}\sim u$ is used.

Assuming the walls to be formed after inflation, the simplest way to make them
disappear is to introduce a small pressure difference Zeldovich et al. (1974);
Vilenkin (1981); Sikivie (1982); Gelmini et al. (1989); Larsson et al. (1997),
a manifestation of a soft $Z_{2}$-breaking term. Such a pressure difference or
equivalently, a bias term in the potential $\Delta V$ needs to be sufficiently
large to ensure DW disappearance prior to the epoch of big bang
nucleosynthesis (BBN) that is, $t_{\rm BBN}>t_{\rm
dec}\approx\sigma_{w}/\Delta V$. It is also important to take care of the fact
that the DW disappear before dominating the universe, requiring
$t_{\text{dec}}<t_{\text{dom}}$, where $t_{\text{dom}}\sim
M_{P}^{2}/\sigma_{w}$. Both of these criteria lead to a lower bound on the
bias term $\Delta V$. However, $\Delta V$ can not be arbitrarily large as it
would otherwise prevent the percolation of both the vacua separated by DW.
Such decaying DW therefore can emit GW Kadota et al. (2015); Hiramatsu et al.
(2014); Krajewski et al. (2016); Nakayama et al. (2017); Dunsky et al. (2021);
Babichev et al. (2022); Ferreira et al. (2022); Deng et al. (2021); Gelmini et
al. (2021); Saikawa (2017). The amplitude the spectrum at peak frequency
$f_{\rm peak}$ can be estimated as Kadota et al. (2015); Hiramatsu et al.
(2014)111Here we are ignoring the friction effects between the walls and the
thermal plasma Nakayama et al. (2017); Gal’tsov et al. (2018).

$\displaystyle\Omega_{\rm GW}h^{2}\,(t_{0})\rvert_{\rm peak}$
$\displaystyle\simeq 5.2\times 10^{-20}\,\tilde{\epsilon}_{\rm
gw}\,A_{w}^{4}\left(\frac{10.75}{g_{*}}\right)^{1/3}$
$\displaystyle\times\left(\frac{\sigma_{w}}{1\,{\rm
TeV}^{3}}\right)^{4}\,\left(\frac{1\,{\rm MeV}^{4}}{\Delta V}\right)^{2}\,,$
(3)

with $t_{0}$ being the present time. Away from the peak, the amplitude varies
as

$\displaystyle\Omega_{\rm GW}\simeq\Omega_{\rm GW}\rvert_{\rm
peak}\times\begin{cases}\displaystyle{\left(\frac{f_{\rm
peak}}{f}\right)}&{\rm for}~{}~{}f>f_{\rm peak}\\\ &\\\
\displaystyle{\left(\frac{f}{f_{\rm peak}}\right)^{3}}&{\rm for}~{}~{}f<f_{\rm
peak}\end{cases}\,,$ (4)

where the peak frequency is given by

$\displaystyle f_{\rm peak}(t_{0})$ $\displaystyle\simeq 3.99\times
10^{-9}\,{\rm Hz}\,A^{-1/2}$ $\displaystyle\times\left(\frac{1\,{\rm
TeV}^{3}}{\sigma_{w}}\right)^{1/2}\,\left(\frac{\Delta V}{1\,{\rm
MeV}^{4}}\right)^{1/2}\,.$ (5)

In the above expressions, $A_{w}$ is the area parameter Caprini et al. (2018);
Paul et al. (2021) $\simeq 0.8$ for DW arising from $Z_{2}$ breaking, and
$\tilde{\epsilon}_{\rm gw}$ is the efficiency parameter $\simeq$ 0.7 Hiramatsu
et al. (2014). Note that the above spectrum can be obtained from a general
parametrisation $S(f/f_{\rm peak})$

$S(x)=\frac{(a+b)^{c}}{(bx^{-a/c}+ax^{b/c})^{c}}$ (6)

for $a=3$ (required by causality) and $b\approx c\approx 1$ (suggested by
simulation Hiramatsu et al. (2014)). However, as noted in Afzal et al. (2023);
Ferreira et al. (2022), the values of $b,c$ may depend upon the specific DW
annihilation mechanism or regime all of which have not been explored in
numerical simulations yet. This allows one to vary $b,c$ to get a better fit
with the PTA data Afzal et al. (2023).

When the GW production is ceased after the annihilation of the domain walls,
the energy density of GW redshifts mimicking that of the SM radiation. As a
result, GW itself acts as an additional source of radiation with the potential
to alter the prediction of BBN. Thus, an excess of the GW energy density
around $T\lesssim\mathcal{O}(\text{MeV})$, can be restricted by considering
the limits on the number of relativistic degrees of freedom from CMB and BBN,
encoded in $\Delta N_{\text{eff}}$. This, in turns, puts a bound on the
amplitude of GW spectrum, demanding $\Omega_{\text{GW}}\,h^{2}\lesssim
5.6\times 10^{-6}\,\Delta N_{\text{eff}}$ Maggiore (2000); Caprini and
Figueroa (2018). Here we consider several projected limits on $\Delta
N_{\text{eff}}$ on top of the existing limit from Planck:
$N_{\text{eff}}=2.99\pm 0.34$ at 95% CL Aghanim et al. (2020). This bound is
shown by the solid gray horizontal line in Fig. 2. Once the baryon acoustic
oscillation (BAO) data are included, the measurement becomes more stringent:
$N_{\text{eff}}=2.99\pm 0.17$. A combined BBN+CMB analysis shows
$N_{\text{eff}}=2.880\pm 0.144$, as computed in Ref. Yeh et al. (2022). This
constraint is denoted by the dashed horizontal line. On the other hand,
upcoming CMB experiments like CMB-S4 Abazajian et al. (2019) and CMB-HD Aiola
et al. (2022) will be able to probe $\Delta N_{\text{eff}}$ as small as $\sim
0.06$ and $\sim 0.027$, respectively. These are indicated by dot-dashed and
dotted lines respectively. The next generation of satellite missions, such as
COrE Bouchet et al. (2011) and Euclid Laureijs et al. (2011), leads to $\Delta
N_{\text{eff}}\lesssim 0.013$, as shown by the large dashed line.

In Fig. 1 we summarize theoretical bounds on the VEV $u$ and the bias term
$\Delta V$, where all the shaded regions are disallowed from (i) decay of the
DWs post BBN (in gray), where $t_{\text{dec}}>1$ sec, (ii) DW domination (in
red) or $t_{\text{dom}}<t_{\text{dec}}$ and (iii) percolation of the vacua
separated by DW (in orange), i.e., an arbitrarily large $\Delta V>u^{4}$. This
leaves us with the white region in-between that is allowed, from where we
choose our benchmark points (BP), as indicated in Tab. 1.

|  $u$ (GeV) | $\Delta V(\text{MeV})^{4}$
---|---|---
BP1 | $2\times 10^{5}$ | $10^{8}$
BP2 | $3\times 10^{5}$ | $10^{8}$
BP3 | $4\times 10^{5}$ | $10^{8}$
Table 1: Details of the benchmark points (BPs) used in Fig. 2 and Fig. 3.

The GW spectrum corresponding to the BPs in Tab. 1 is illustrated in Fig. 2.
As explained before, we distinctively see a blue-tilted pattern for
$f<f_{\text{peak}}$, while the spectrum is red-tilted in the opposite limit.
Here we project limits from BBO Yagi and Seto (2011), LISA Amaro-Seoane et al.
(2012), DECIGO Kawamura et al. (2006), ET Punturo et al. (2010), CE Abbott et
al. (2017), THEIA Garcia-Bellido et al. (2021), HL (aLIGO) Aasi et al. (2015),
$\mu$ARES Sesana et al. (2021) and SKA Weltman et al. (2020). In this plot,
the range of GW spectrum from NANOGrav results Agazie et al. (2023) is shown
by the red points. The gray-shaded region is completely disallowed from
$\Delta N_{\text{eff}}$ bound on overproduction of $\Omega_{\text{GW}}$ as
discussed before, depending on the sensitivity of a particular experiment. As
one can already notice, BP3 is already ruled out from PLANCK bound, while BP1
is beyond the reach of any future experiments proposed so far.

Figure 1: Bound on the size of VEV and the bias term. All shaded regions are
disallowed (see text for details). Figure 2: Spectrum of gravitational wave
from DW decay, where we show sensitivities of several GW experiments. The
black curves correspond to the chosen benchmark points. The gray region marked
as “$\Delta N_{\text{eff}}$ constraint” is disallowed from overproduction of
GW energy density (see text for details). The red points correspond to the
NanoGrav 15 yr observation Agazie et al. (2023). Figure 3: Viable parameter
space in the bi-dimensional plane of $T_{\text{ann}}$-$r_{w}$, based on the
NanoGrav 15 yr dataset Agazie et al. (2023). The black solid contour
corresponds to the $\Delta V$-value relevant for the BPs in Tab. 1. The red
and blue contours correspond to 1- and 2$\sigma$ CL respectively.

Corresponding to the two epochs $t_{\text{dom}}$ and $t_{\text{ann}}$ stated
before, we define two temperatures $T_{\text{dom}}$ and $T_{\text{ann}}$. We
are typically interested in the regime $T_{\text{ann}}>T_{\text{dom}}$, i.e.,
the DWs disappear before they dominant the Universe. Following Bai et al.
(2023), $T_{\text{ann}}$ reads

$\displaystyle T_{\text{ann}}\simeq 120\,\text{MeV}\,\sqrt{\frac{\Delta
V/\text{MeV}^{4}}{10^{8}}}\,\left(\frac{A_{w}}{0.8}\right)^{-1/2}$
$\displaystyle\left(\frac{\sigma_{w}/\text{GeV}^{3}}{10^{16}}\right)^{-1/2}\,\left(\frac{g_{*}(T_{\text{ann}})}{10}\right)^{-1/4}\,,$
(7)

implying, for larger surface tension, it takes longer the walls to collapse,
while for larger bias the opposite happens. We also define another quantity

$\displaystyle r_{w}=\frac{\rho_{r}}{1+\rho_{r}}\,,$ (8)

where

$\displaystyle\rho_{r}=\frac{\rho_{w}(T_{\text{ann}})}{\rho_{R}(T_{\text{ann}})}\simeq
0.14\,\left(\frac{A_{w}}{0.8}\right)^{2}\,\left(\frac{\sigma_{w}/\text{GeV}^{3}}{10^{16}}\right)^{2}\,\left(\frac{10^{8}}{\Delta
V\,\text{MeV}^{4}}\right)\,,$ (9)

which quantifies the energy density contained within the DW compared to that
of radiation. We show the compatibility of our relevant model parameters,
namely, the bias term $\Delta V$ and the VEV $u$ (equivalently, the strain
$\sigma$) in Fig. 3 with the NANOGrav data, utilizing Eq. (Scale of Dirac
leptogenesis and left-right symmetry in the light of recent PTA results) and
Eq. (8). We superimpose the 1 and 2-sigma contours (shown by red and blue
solid curves) provided by the NANOGrav result Agazie et al. (2023). We see,
BP1 lies well within the $1\sigma$ contour, while the other two BPs are well
off. It is possible to derive a lower and upper bound on the VEV for a fixed
$\Delta V$, as denoted by the gray dashed horizontal lines for $\Delta
V=10^{8}\,\text{MeV}^{4}$. Thus, for $\Delta V=10^{8}\,\text{MeV}^{4}$, we
find $189\lesssim u/\text{TeV}\lesssim 225$, in compliance with NANOGrav
$2\sigma$ contour. On the other hand, the viable range of bias term, that lies
within $2\sigma$ CL of NANOGrav result, turns out to be $4\times
10^{6}\lesssim\Delta V/\text{MeV}^{4}\lesssim 5\times 10^{11}$, as shown by
the green dotted and orange dashed curves. Depending on the choice of $\Delta
V$, the upper limit on $u$ can be pushed to larger values. Note that, the
constraint on VEV mentioned above (together with the limit on $\Delta V$),
satisfies the bounds shown in Fig. 1.

Consequence for Dirac leptogenesis: In the minimal model of Dirac leptogenesis
or Dirac neutrino seesaw Barman et al. (2022), the standard model (SM) is
extended by three copies of vector-like neutral singlet fermions $N_{L,R}$ and
three copies of right chiral part $(\nu_{R})$ of light Dirac neutrinos. A real
singlet scalar field $\varphi$ is introduced to couple $\nu_{R}$ with $N$. A
$Z_{2}$ symmetry under which $\varphi,\nu_{R}$ are odd, prevents direct
coupling of the SM lepton doublet $L$ with $\nu_{R}$ via SM Higgs $H$. After
the neutral components of $H$ and $\varphi$ acquire VEV $v,\,u$ respectively,
light Dirac neutrino mass arises from the Type-I seesaw equivalent for Dirac
neutrino as $m_{\nu}\propto y^{2}\,u\,v/M_{N}$ with $M_{N}$ being the scale of
Dirac seesaw and $y$ being Dirac Yukawa couplings of $N$ with
$\nu_{L},\nu_{R}$ considered to be equal. The same heavy fermions $N_{L,R}$
can have out-of-equilibrium decay to achieve successful Dirac leptogenesis.

Now, for $u\gtrsim 190$ TeV, preferred from NANOGrav 2023 data as discussed
before, and considering light neutrino mass $m_{\nu}\leq 0.1$ eV, we get
$M_{N}>y^{2}\times 10^{17}$ GeV. This implies that, for order one Yukawa
couplings, the scale of Dirac leptogenesis is even above the upper limit on
reheating temperature, disfavouring the possibility of thermal Dirac
leptogenesis. If the Yukawa couplings are made as low as electron Yukawa
coupling, we have $M_{N}>10^{7}$ GeV, keeping it at intermediate scale. To
summarise, the possibility of low scale Dirac leptogenesis is disfavoured
unless we tune Yukawa couplings involved in Dirac seesaw more than what we
have in the SM.

Consequence for left-right symmetry: Left-right symmetric models Pati and
Salam (1974); Mohapatra and Pati (1975a, b); Senjanovic and Mohapatra (1975);
Senjanović (1979); Mohapatra et al. (1978); Mohapatra and Marshak (1980);
Mohapatra and Senjanovic (1981); Lim and Inami (1982); Gunion et al. (1989);
Deshpande et al. (1991); Duka et al. (2000); Fileviez Perez (2009) have been
one of the most well-studied BSM frameworks where the SM gauge symmetry is
extended to $\rm SU(3)_{c}\times SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$.
In addition to the gauge symmetry, the model also has an in-built $Z_{2}$
symmetry or parity symmetry under which the left and right sector fields are
interchanged, keeping the framework parity or left-right symmetric. While a
detailed discussion related to GW probe of LRSM can be found in Craig et al.
(2021); Borah and Dasgupta (2022), here we briefly summarise the implications
of 2023 PTA data on the scale of left-right symmetry.

Figure 4: Same as Fig. 3, but for the minimal LRSM. Along the black solid line
the left-right symmetry breaking scale $u$ is varying.

Unlike in the $Z_{2}$-odd scalar singlet model discussed above, the LRSM gauge
symmetry does not allow arbitrary bias terms. It is possible to generate such
bias term via higher dimensional operators invariant under gauge symmetry but
explicitly breaking the parity symmetry. While Planck scale effects are
expected to break any global symmetries like parity Abbott and Wise (1989);
Kallosh et al. (1995); Hawking (1975), the corresponding bias term can lead to
DW disappearance Rai and Senjanovic (1994); Lew and Riotto (1993). As shown in
Borah and Dasgupta (2022), the bias term in the minimal model is $\Delta V\sim
u^{6}/M_{P}^{2}$, where $u$ is the $SU(2)_{R}\times U(1)_{B-L}$ as well as
parity breaking scale. Due to the dependence of the bias term on the scale of
left-right breaking, the constraint on left-right symmetry breaking scale is
stronger than what we had on the scale of Dirac leptogenesis. As shown in Fig.
4, the scale of left-right symmetry should be approximately around $\sim
10^{6}$ GeV (shown by the black solid horizontal line), in order to be in
agreement with NANOGrav 15 yr data at $2\sigma$.

Conclusion: We have investigated the consequence of the recent PTA results on
the scale of Dirac leptogenesis and left-right symmetric model. In minimal
version of both these scenarios, domain walls arise due to spontaneous
breaking of a discrete $Z_{2}$ symmetry. While the bound on the scale of Dirac
leptogenesis from PTA data depend upon the size of Dirac Yukawa couplings, for
conservative choice of such couplings with fine-tuning at the level of the SM,
we find a lower bound $M_{N}>10^{7}$ GeV, keeping leptogenesis at intermediate
scales. However, for order one Yukawa coupling, this bound is much stronger
$M_{N}>10^{17}$ GeV keeping only non-thermal Dirac leptogenesis option viable.
Due to the constrained structure of the LRSM, we get much tighter constraint
on the scale of left-right breaking namely $M_{\rm LR}\sim 10^{6}$ GeV, in
order to satisfy NANOGrav 15 yr data, keeping the model out of reach from
direct search experiments like colliders. Future data from PTA or other GW
experiments are expected to shed more light on the parameter space of this
model, by constraining the spectrum at higher frequencies.

###### Acknowledgements.

The work of D.B. is supported by the science and engineering research board
(SERB), Government of India grant MTR/2022/000575.

## References

  * Agazie et al. (2023) G. Agazie et al. (NANOGrav) (2023), eprint 2306.16213.
  * Antoniadis et al. (2023) J. Antoniadis et al. (2023), eprint 2306.16214.
  * Reardon et al. (2023) D. J. Reardon et al. (2023), eprint 2306.16215.
  * Xu et al. (2023) H. Xu et al. (2023), eprint 2306.16216.
  * Afzal et al. (2023) A. Afzal et al. (NANOGrav) (2023), eprint 2306.16219.
  * Ghoshal and Strumia (2023) A. Ghoshal and A. Strumia (2023), eprint 2306.17158.
  * Shen et al. (2023) Z.-Q. Shen, G.-W. Yuan, Y.-Y. Wang, and Y.-Z. Wang (2023), eprint 2306.17143.
  * Yang et al. (2023) J. Yang, N. Xie, and F. P. Huang (2023), eprint 2306.17113.
  * Guo et al. (2023) S.-Y. Guo, M. Khlopov, X. Liu, L. Wu, Y. Wu, and B. Zhu (2023), eprint 2306.17022.
  * Ellis et al. (2023a) J. Ellis, M. Fairbairn, G. Hütsi, J. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe (2023a), eprint 2306.17021.
  * Megias et al. (2023) E. Megias, G. Nardini, and M. Quiros (2023), eprint 2306.17071.
  * Fujikura et al. (2023) K. Fujikura, S. Girmohanta, Y. Nakai, and M. Suzuki (2023), eprint 2306.17086.
  * Han et al. (2023) C. Han, K.-P. Xie, J. M. Yang, and M. Zhang (2023), eprint 2306.16966.
  * Zu et al. (2023) L. Zu, C. Zhang, Y.-Y. Li, Y.-C. Gu, Y.-L. S. Tsai, and Y.-Z. Fan (2023), eprint 2306.16769.
  * Franciolini et al. (2023a) G. Franciolini, A. Iovino, Junior., V. Vaskonen, and H. Veermae (2023a), eprint 2306.17149.
  * Li et al. (2023) Y. Li, C. Zhang, Z. Wang, M. Cui, Y.-L. S. Tsai, Q. Yuan, and Y.-Z. Fan (2023), eprint 2306.17124.
  * Kitajima et al. (2023) N. Kitajima, J. Lee, K. Murai, F. Takahashi, and W. Yin (2023), eprint 2306.17146.
  * Bai et al. (2023) Y. Bai, T.-K. Chen, and M. Korwar (2023), eprint 2306.17160.
  * Vagnozzi (2023) S. Vagnozzi (2023), eprint 2306.16912.
  * Ellis et al. (2023b) J. Ellis, M. Lewicki, C. Lin, and V. Vaskonen (2023b), eprint 2306.17147.
  * Wang et al. (2023) Z. Wang, L. Lei, H. Jiao, L. Feng, and Y.-Z. Fan (2023), eprint 2306.17150.
  * Lambiase et al. (2023) G. Lambiase, L. Mastrototaro, and L. Visinelli (2023), eprint 2306.16977.
  * Franciolini et al. (2023b) G. Franciolini, D. Racco, and F. Rompineve (2023b), eprint 2306.17136.
  * Fukugita and Yanagida (1986) M. Fukugita and T. Yanagida, Phys. Lett. B174, 45 (1986).
  * Dick et al. (2000) K. Dick, M. Lindner, M. Ratz, and D. Wright, Phys. Rev. Lett. 84, 4039 (2000), eprint hep-ph/9907562.
  * Murayama and Pierce (2002) H. Murayama and A. Pierce, Phys. Rev. Lett. 89, 271601 (2002), eprint hep-ph/0206177.
  * Dror et al. (2020) J. A. Dror, T. Hiramatsu, K. Kohri, H. Murayama, and G. White, Phys. Rev. Lett. 124, 041804 (2020), eprint 1908.03227.
  * Blasi et al. (2020) S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Res. 2, 043321 (2020), eprint 2004.02889.
  * Fornal and Shams Es Haghi (2020) B. Fornal and B. Shams Es Haghi, Phys. Rev. D 102, 115037 (2020), eprint 2008.05111.
  * Samanta and Datta (2021) R. Samanta and S. Datta, JHEP 05, 211 (2021), eprint 2009.13452.
  * Barman et al. (2022) B. Barman, D. Borah, A. Dasgupta, and A. Ghoshal (2022), eprint 2205.03422.
  * Craig et al. (2021) N. Craig, I. Garcia Garcia, G. Koszegi, and A. McCune, JHEP 09, 130 (2021), eprint 2012.13416.
  * Borah and Dasgupta (2022) D. Borah and A. Dasgupta, Phys. Rev. D 106, 035016 (2022), eprint 2205.12220.
  * Zeldovich et al. (1974) Y. B. Zeldovich, I. Y. Kobzarev, and L. B. Okun, Zh. Eksp. Teor. Fiz. 67, 3 (1974).
  * Kibble (1976) T. W. B. Kibble, J. Phys. A 9, 1387 (1976).
  * Vilenkin (1981) A. Vilenkin, Phys. Rev. D 23, 852 (1981).
  * Coulson et al. (1996) D. Coulson, Z. Lalak, and B. A. Ovrut, Phys. Rev. D 53, 4237 (1996).
  * Krajewski et al. (2021) T. Krajewski, J. H. Kwapisz, Z. Lalak, and M. Lewicki, Phys. Rev. D 104, 123522 (2021), eprint 2103.03225.
  * Sikivie (1982) P. Sikivie, Phys. Rev. Lett. 48, 1156 (1982).
  * Gelmini et al. (1989) G. B. Gelmini, M. Gleiser, and E. W. Kolb, Phys. Rev. D 39, 1558 (1989).
  * Larsson et al. (1997) S. E. Larsson, S. Sarkar, and P. L. White, Phys. Rev. D 55, 5129 (1997), eprint hep-ph/9608319.
  * Kadota et al. (2015) K. Kadota, M. Kawasaki, and K. Saikawa, JCAP 10, 041 (2015), eprint 1503.06998.
  * Hiramatsu et al. (2014) T. Hiramatsu, M. Kawasaki, and K. Saikawa, JCAP 02, 031 (2014), eprint 1309.5001.
  * Krajewski et al. (2016) T. Krajewski, Z. Lalak, M. Lewicki, and P. Olszewski, JCAP 12, 036 (2016), eprint 1608.05719.
  * Nakayama et al. (2017) K. Nakayama, F. Takahashi, and N. Yokozaki, Phys. Lett. B 770, 500 (2017), eprint 1612.08327.
  * Dunsky et al. (2021) D. Dunsky, L. J. Hall, and K. Harigaya, JHEP 01, 125 (2021), eprint 2007.12711.
  * Babichev et al. (2022) E. Babichev, D. Gorbunov, S. Ramazanov, and A. Vikman, JCAP 04, 028 (2022), eprint 2112.12608.
  * Ferreira et al. (2022) R. Z. Ferreira, A. Notari, O. Pujolas, and F. Rompineve (2022), eprint 2204.04228.
  * Deng et al. (2021) X. Deng, X. Liu, J. Yang, R. Zhou, and L. Bian, Phys. Rev. D 103, 055013 (2021), eprint 2012.15174.
  * Gelmini et al. (2021) G. B. Gelmini, S. Pascoli, E. Vitagliano, and Y.-L. Zhou, JCAP 02, 032 (2021), eprint 2009.01903.
  * Saikawa (2017) K. Saikawa, Universe 3, 40 (2017), eprint 1703.02576.
  * Gal’tsov et al. (2018) D. Gal’tsov, E. Melkumova, and P. Spirin, JHEP 01, 120 (2018), eprint 1711.01114.
  * Caprini et al. (2018) C. Caprini, M. C. Guzzetti, and L. Sorbo, Class. Quant. Grav. 35, 124003 (2018), eprint 1707.09750.
  * Paul et al. (2021) A. Paul, U. Mukhopadhyay, and D. Majumdar, JHEP 05, 223 (2021), eprint 2010.03439.
  * Maggiore (2000) M. Maggiore, Phys. Rept. 331, 283 (2000), eprint gr-qc/9909001.
  * Caprini and Figueroa (2018) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, 163001 (2018), eprint 1801.04268.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Yeh et al. (2022) T.-H. Yeh, J. Shelton, K. A. Olive, and B. D. Fields, JCAP 10, 046 (2022), eprint 2207.13133.
  * Abazajian et al. (2019) K. Abazajian et al. (2019), eprint 1907.04473.
  * Aiola et al. (2022) S. Aiola et al. (CMB-HD) (2022), eprint 2203.05728.
  * Bouchet et al. (2011) F. R. Bouchet et al. (COrE) (2011), eprint 1102.2181.
  * Laureijs et al. (2011) R. Laureijs et al. (EUCLID) (2011), eprint 1110.3193.
  * Yagi and Seto (2011) K. Yagi and N. Seto, Phys. Rev. D 83, 044011 (2011), [Erratum: Phys.Rev.D 95, 109901 (2017)], eprint 1101.3940.
  * Amaro-Seoane et al. (2012) P. Amaro-Seoane, H. Audley, S. Babak, J. M. Baker, E. Barausse, P. L. Bender, E. Berti, P. Binétruy, M. Born, D. Bortoluzzi, et al. (2012).
  * Kawamura et al. (2006) S. Kawamura et al., Class. Quant. Grav. 23, S125 (2006).
  * Punturo et al. (2010) M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia, M. Beker, et al., Classical and Quantum Gravity 27, 194002 (2010), URL https://doi.org/10.1088/0264-9381/27/19/194002.
  * Abbott et al. (2017) B. P. Abbott et al. (LIGO Scientific), Class. Quant. Grav. 34, 044001 (2017), eprint 1607.08697.
  * Garcia-Bellido et al. (2021) J. Garcia-Bellido, H. Murayama, and G. White (2021), eprint 2104.04778.
  * Aasi et al. (2015) J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), eprint 1411.4547.
  * Sesana et al. (2021) A. Sesana et al., Exper. Astron. 51, 1333 (2021), eprint 1908.11391.
  * Weltman et al. (2020) A. Weltman et al., Publ. Astron. Soc. Austral. 37, e002 (2020), eprint 1810.02680.
  * Pati and Salam (1974) J. C. Pati and A. Salam, Phys. Rev. D10, 275 (1974), [Erratum: Phys. Rev.D11,703(1975)].
  * Mohapatra and Pati (1975a) R. N. Mohapatra and J. C. Pati, Phys. Rev. D11, 566 (1975a).
  * Mohapatra and Pati (1975b) R. N. Mohapatra and J. C. Pati, Phys. Rev. D11, 2558 (1975b).
  * Senjanovic and Mohapatra (1975) G. Senjanovic and R. N. Mohapatra, Phys. Rev. D12, 1502 (1975).
  * Senjanović (1979) G. Senjanović, Nucl.Phys. B153, 334 (1979).
  * Mohapatra et al. (1978) R. N. Mohapatra, F. E. Paige, and D. P. Sidhu, Phys. Rev. D17, 2462 (1978).
  * Mohapatra and Marshak (1980) R. N. Mohapatra and R. E. Marshak, Phys. Rev. Lett. 44, 1316 (1980), [Erratum: Phys. Rev. Lett.44,1643(1980)].
  * Mohapatra and Senjanovic (1981) R. N. Mohapatra and G. Senjanovic, Phys. Rev. D23, 165 (1981).
  * Lim and Inami (1982) C. S. Lim and T. Inami, Prog. Theor. Phys. 67, 1569 (1982).
  * Gunion et al. (1989) J. F. Gunion, J. Grifols, A. Mendez, B. Kayser, and F. I. Olness, Phys. Rev. D40, 1546 (1989).
  * Deshpande et al. (1991) N. G. Deshpande, J. F. Gunion, B. Kayser, and F. I. Olness, Phys. Rev. D44, 837 (1991).
  * Duka et al. (2000) P. Duka, J. Gluza, and M. Zralek, Annals Phys. 280, 336 (2000), eprint hep-ph/9910279.
  * Fileviez Perez (2009) P. Fileviez Perez, JHEP 03, 142 (2009), eprint 0809.1202.
  * Abbott and Wise (1989) L. F. Abbott and M. B. Wise, Nucl. Phys. B325, 687 (1989).
  * Kallosh et al. (1995) R. Kallosh, A. D. Linde, D. A. Linde, and L. Susskind, Phys. Rev. D52, 912 (1995), eprint hep-th/9502069.
  * Hawking (1975) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Rai and Senjanovic (1994) B. Rai and G. Senjanovic, Phys. Rev. D 49, 2729 (1994), eprint hep-ph/9301240.
  * Lew and Riotto (1993) H. Lew and A. Riotto, Phys. Lett. B 309, 258 (1993), eprint hep-ph/9304203.


aainstitutetext: Department of Physics & State Key Laboratory of Surface
Physics, Fudan University, Shanghai 200433, China

# Digital Pseudoquantum Simulation of $\mathbb{Z}_{2}$ Gauge Higgs Model

Yiming Ding a Xiaopeng Cui a,1 Yu Shi 111Corresponding author.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We perform a digital pseudoquantum simulation of $\mathbb{Z}_{2}$ gauge Higgs
model on a $3\times 3$ lattice. First we propose the quantum algorithm for the
digital quantum simulation, based on Trotter decomposition, quantum adiabatic
algorithm and its circuit realization. Then we classically demonstrate it in a
GPU simulator, obtaining useful results, which indicate the topological
properties of deconfined phase and clarify the phase diagram. Especially, our
work suggests that the tricitical point, where the two critical lines of
second-order transitions meet, lies on the critical line of the first-order
transition rather than its end.

##  1 Introduction

Lattice gauge theory is a nonperturbative approach to gauge theory, and is
important in particle physics, condensed matter physics and even topological
quantum computing PhysRevD.10.2445 . It is usually implemented in terms of
Monte Carlo simulation, which, however, lacks real-time dynamics and may
suffer the well known sign problem of fermions.

Recently, quantum simulation has appeared to be able to dissolve these
problems Martinez2016 ; Gonz_lez_Cuadra_2017 ; PhysRevX.10.021041 ;
JHEP08(2020)160 . As done for pure $\mathbb{Z}_{2}$ gauge theory
JHEP08(2020)160 ; S0217979220502926 , for a quantum simulation involving only
dozens of qubits, it is very useful to make classical demonstration in a high-
performance platform, which we call pseudoquantum simulation. It not only
serves as benchmark for real quantum simulation and facilitates the
development of quantum algorithms, but is also a new numerical method for
computational problems.

An extension of the $\mathbb{Z}_{2}$ pure gauge theory is the $\mathbb{Z}_{2}$
gauge Higgs model, which contains coupling between matter and gauge fields.
There exists duality between the gauge field and the matter field in this
model, which is equivalent to the transverse-field toric code model
PhysRevB.82.085114 . This model has been widely studied analytically
PhysRevD.19.3682 ; Banks:1980gs ; Nussinov:2004ns ; Vidal:2008uy  and
numerically Jongeward:1980wx ; PhysRevD.21.1006 ; PhysRevB.82.085114 ;
2012.15845 . It has been known that in this model (cf. Fig. 1), there is a
deconfined phase, separated from a confined phase on one hand, and from a so-
called Higgs phase on the other. The phase transition between the deconfined
and the confined or Higgs phase is second-order, leading to a topological
region surrounded by two second-order lines on the phase diagram. These two
lines meet at a self-dual point. The confined and Higgs phases are separated
by a finite dual line of a first-order transition, beyond which the two phases
are continuously connected.

However, with strong competition between matter and gauge fields, how these
critical lines are connected and where the two second-order lines meet were
not yet clear previously. It was pointed out that there are three
possibilities PhysRevB.82.085114 . A quantum Monte Carlo (QMC) study provides
the evidence that the tricritical point, where the three critical lines meet,
is scale-invariant and second-order 2012.15845 .

In this paper, we report a digital pseudoquantum simulation of
$\mathbb{Z}_{2}$ gauge Higgs model. It is digital as it is based on Trotter
decomposition of the unitary evolution. It realizes the quantum adiabatic
algorithm in terms of quantum circuits. Our pseudoquantum simulation is its
classical demonstration by using a GPU simulator called Quantum Exact
Simulation Toolkit (QuEST) in a NVIDIA GeForce RTX 3090 GPU server.

We have investigated the quantum phase transition, the dynamic process along
the dual line on the phase diagram, and the behavior near the tricritical
point. Especially, our work suggests that the tricritical point, where the two
second-order lines end, lies on the line of the first-order transition, rather
than at the lower end of it (see Fig. 1).

Figure 1: The phase diagram from our calculation. The red points represent the
second-order transitions and the blue line represents first-order transitions.
The first-order transition line has not only a part outside the deconfined
phase, but also a small part within the deconfined phase.

The rest of the paper is arranged as the following. In Sec. 2, we briefly
introduce the $\mathbb{Z}_{2}$ gauge Higgs model. In Sec. 3, we elaborate the
digital quantum adiabatic algorithm and its realization in terms of quantum
circuits. In Sec. 4, we describe the preparation of initial state, as well as
the topological phase. In Sec. 5, we analyze the critical points according to
the analysis of the density of states (DOS). The Trotter error is analyzed in
Sec. 6, with some details in the Appendix. A summary is made in Sec. 7.

##  2 Model description

The Hamiltonian of the $\mathbb{Z}_{2}$ guage-Higgs model is

$\begin{split}H=-J\sum_{v}\tau^{x}_{v}-h\sum_{p}B_{p}^{z}-g\sum_{l}\sigma^{x}_{l}-\lambda\sum_{\langle
j,k\rangle}\tau^{z}_{j}\sigma^{z}_{\langle
j,k\rangle}\tau_{k}^{z},\end{split}$ (1)

where $\sigma$ denotes gauge fields on links and $\tau$ denotes Ising matter
fields at vertices,

$B^{z}_{p}=\prod_{j\in\partial p}\sigma^{z}_{j},$ (2)

is the tensor product of four $\sigma^{z}$’s on the sides of a plaquette $p$.
The Gauss’ law is given by

$Q_{v}^{x}=\tau^{x}_{v}\prod_{j\in\partial
v}\sigma^{x}_{j}\equiv\tau^{x}_{v}A_{v}^{x}.$ (3)

Under a mathematical mapping, Eq. (1) is equivalent to a toric code model with
two transverse fields PhysRevB.82.085114 , which is used in our work. We focus
on the parameter subspace with $J=h=1/2$, that is,

$H=-\frac{1}{2}\sum_{v}A_{v}^{x}-\frac{1}{2}\sum_{p}B_{p}^{z}-g\sum_{l}\sigma^{x}_{l}-\lambda\sum_{l}\sigma^{z}_{l}.$
(4)

The total number of qubits is 19, one for ancilla and the others for a
$3\times 3$ lattice model on torus.

##  3 Digital quantum simulation based on quantum adiabatic algorithm

Our purpose is to obtain the energy of the system as a function of parameters
$\lambda$ and $g$. Since this model possesses an intrinsic duality, we only
need to investigate the behavior below the dual line $\lambda=g$ on the
parameter plane $\lambda-g$. The $g$ axis in this model represents pure
$\mathbb{Z}_{2}$ gauge theory, for which we prepare the initial ground state
on a point $P(p,0)$ on the $g$ axis.

Then two kinds of paths of parameter variation are used for the adiabatic
algorithm (See FIG. 2). On a red path, the parameters vary first to $R(p,r)$,
then to $R^{\prime}(0,r)$. On a blue path, the parameters vary from $P(p,0)$
to $P^{\prime}(0,p)$ on the straight line.

Figure 2: Two paths in our simulations.

In this section, we introduce the algorithm for the digital quantum
simulation, which is implemented by using quantum circuit, and use quantum
adiabatic algorithm to evolve the ground state along paths in the parameter
space. For convenience, we write

$H=H_{1}+H_{2}+H_{3},$ (5)

where

$\begin{split}&H_{1}=-g\sum_{l}\sigma^{x}_{l},\\\
&H_{2}=-\frac{1}{2}\sum_{v}A_{v}^{x}-\frac{1}{2}\sum_{p}B_{p}^{z},\\\
&H_{3}=-\lambda\sum_{l}\sigma^{z}_{l}\end{split}.$ (6)

We decompose the evolution operator $e^{-iHt}$ with second-order Trotter-
Suzuki formula Hatano2005 ,

$\begin{split}e^{-iHt}\approx
e^{-iH_{1}\frac{t}{2}}e^{-iH_{2}\frac{t}{2}}e^{-iH_{3}t}e^{-iH_{2}\frac{t}{2}}e^{-iH_{1}\frac{t}{2}},\end{split}$
(7)

where

$\begin{split}&e^{-iH_{1}\frac{t}{2}}=\prod_{l}e^{ig\sigma^{x}_{l}\frac{t}{2}},\\\
&e^{-iH_{2}\frac{t}{2}}=\prod_{v}e^{i\frac{1}{2}A^{x}_{v}\frac{t}{2}}\prod_{p}e^{i\frac{1}{2}B^{z}_{p}\frac{t}{2}},\\\
&e^{-iH_{3}t}=\prod_{l}e^{i\lambda\sigma^{z}_{l}t}.\\\ \end{split}$ (8)

This decomposition does not generate errors, since the summands in $H_{j}$
commute with each other.

$e^{-iH_{1}\frac{t}{2}}$ and $e^{-iH_{3}{t}}$ can be realized by using two
sets of rotation gates $\mathrm{R}_{x}$’s and $\mathrm{R}_{z}$’s

$\begin{split}&\prod_{l}e^{ig\sigma^{x}_{l}\frac{t}{2}}=\prod_{l}\mathrm{R}_{x}^{l}(-gt),\\\
&\prod_{l}e^{i\lambda\sigma^{z}_{l}{t}}=\prod_{l}\mathrm{R}_{z}^{l}(-2\lambda
t).\\\ \end{split}$ (9)

By introducing an ancilla $a$, each $e^{i\frac{1}{2}B_{p}^{z}\frac{t}{2}}$ can
be realized as

$\begin{split}&e^{i\frac{1}{2}B_{p}^{z}\frac{t}{2}}=\mathrm{U}_{1}[\mathrm{R}_{z}(-\frac{t}{2})]_{a}\mathrm{U}_{1}^{\dagger},\end{split}$
(10)

where

$\begin{split}\mathrm{U}_{1}=\prod_{j\in p}\mathrm{CNOT}_{j\to a}.\end{split}$
(11)

$e^{i\frac{1}{2}A_{v}^{x}\frac{t}{2}}$ can be realized similarly, and we only
need 4 additional Hadamard gates to switch into $z$ basis the 4 spins on the
sides connected at $v$, that is

$\begin{split}&e^{i\frac{1}{2}A_{v}^{x}\frac{t}{2}}=\mathrm{U_{2}}^{\dagger}[\mathrm{R}_{z}(-\frac{t}{2})]\mathrm{U_{2}},\\\
\end{split}$ (12)

where

$\begin{split}\mathrm{U_{2}}=[\prod_{k\in v}\mathrm{CNOT}_{k\to
a}][\prod_{k\in v}{\mathrm{H}}_{k}].\\\ \end{split}$ (13)

Note that we use the convention for the time order of operators, that is, from
right to left. We omit drawing the circuits, which is straightforward.

To ensure the adiabaticity in the variation of the parameters, evolution on
each path is divided into numerous tiny steps, each with the same time
duration.

##  4 Preparation of initial state

The initial ground state at $P(p,0)$ is the same as an intermediate state in
our previous work on the pure $\mathbb{Z}_{2}$ gauge theory JHEP08(2020)160 ,
as the mere addition of $A_{v}$ does not affect the state.

To prepare the ground state from $|00\cdots 00\rangle$, we use the circuit of

$\prod_{p}\\{[\prod_{j\in p}\mathrm{CNOT}_{j\to
a}][\mathrm{P}_{a=|0\rangle}][\prod_{j\in p}\mathrm{CNOT}_{j\to
a}]\\}\prod_{l}H_{l}.$ (14)

The idea is to first generate the equal superposition of all the computational
basis states, by using a set of hadamard gates on all the qubits, then for
each plaquette, to apply four CNOT gates to transfer the relevant information
on each basis state onto an ancilla. P represents the projection of the
ancilla to be $|0\rangle$, which is feasible in the pseudoquantum computing.
It leads to a superposition of all possible configurations with $B_{p}=1$ for
every $p$, which is also the ground state of a toric code at $O(0,0)$ on
$\lambda-g$ plane. The state prepared by using (14) is denoted as $\psi_{1}$.

Then we adiabatically evolve $\psi_{1}$ along $g$ axis, by using the circuit-
based digital quantum adiabatic algorithm in Sec. 3. Finally we reach the
ground state at $P(p,0)$.

Note that if we initially prepare a ground state at
$P^{\prime}(0,p^{\prime})$, $\psi_{1}$ should be replaced correspondingly
according to the duality.

As the ground state of $\mathbb{Z}_{2}$ gauge theory can be described by ’t
Hooft loop operator $V_{\mu}^{x}$ and non-contractible Wilson loop operator
$W_{\nu}^{z}$, the eigenstates of $V_{\mu}^{x}$ on $g$ axis are dual to the
eigenstates of $W_{\nu}^{z}$ on the $\lambda$ axis. To obtain the ground state
at $P^{\prime}$, the state prepared at $O(0,0)$ is
$\psi_{1}^{\prime}=(\psi_{1}+\psi_{2}+\psi_{3}+\psi_{4})/2$, where $\psi_{2}$,
$\psi_{3}$ and $\psi_{4}$ are degenerate with $\psi_{1}$ at $O$.

Similarly, the quantum circuit for $\psi_{1}^{\prime}$ is

$\prod_{v}\\{[\mathrm{H}_{k}\prod_{k\in v}\mathrm{CNOT}_{k\to
a}][\mathrm{P}_{a=|0\rangle}][\prod_{k\in v}\mathrm{CNOT}_{k\to
a}\mathrm{H}_{k}]\\}.$ (15)

$\psi_{1}$ and $\psi_{1}^{\prime}$ are also degenerate at $O$. If the
evolution is restricted in a topological phase, they cannot evolve to each
other because of the topological protection. In other words, the evolution
from $P(0,p)$ to $P^{\prime}(0,p^{\prime})$ should be different from that on
the reversed path.

We prepare each of these two ground states on the corresponding parameter
point and evolve it to the other parameter point along the blue path in Fig.
2. The result shows that these two energy curves cross if the path is chosen
within the deconfined phase, and do not cross if the path is outside of the
deconfined phase, as indicated in Fig. 3.

Figure 3: Energy as a function of $\lambda$ along the blue path. (a) For
$p=p^{\prime}=0.12$, $E^{\prime}=-9.16003$ and $E=-9.15039$, where
$E^{\prime}$ is the energy for the state evolved from that prepared at $P$,
while $E$ is the energy for the state prepared at $P^{\prime}$. Two energy
curves cross, indicating a topological phase. (b) For $p=p^{\prime}=0.50$,
$E^{\prime}=-14.06550$ and $E=-14.06553$. Two energy curves nearly coincide
and the energy is symmetric with respect to the center of the blue path,
indicating a nontopological phase.

##  5 Phase diagram

We now determine the critical points of quantum phase transitions as the
extreme points in the second derivatives of the energy.

After numerous trials, we choose $p=0.5$ for the red path in Fig. 2, of which
the $\lambda$ value denoted as $r$. The results show that for $r<0.19$, the
critical value of $g$ is $g_{c}\approx 0.19$, while for $r>0.19$,
$g_{c}\approx r$, as depected in Fig. 4. In this way, we find that the two
critical lines surrounding the deconfined phase meet at $g\approx 0.24$, as
shown in Fig. 1.

Figure 4:  For the $\lambda$ value $r=0.05,0.15,0.25$ on red paths (see Fig.
2) , (a) energy, (b) first derivative of the energy, (c) second derivative of
the energy, as functions of $g$. The minima in (c) are at $0.1920$, $0.1916$
and $0.2489$ respectively.

As shown in Fig. 5, we determine the order of the quantum phase transition
according to the DOS (density of states) of the eigenstates of
$\tilde{Z}=-\sum_{l}\sigma^{z}_{l}$, as in our previous work. When $r<0.16$,
the DOS of $\tilde{Z}$ shows only one maximum before and after the critical
point, respectively. This is same as in the pure $\mathbb{Z}_{2}$ gauge
theory. Hence the quantum phase transition is of second-order. This verifies
that the two critical lines surrounding the deconfined phase are second-order.

Such kind of DOS becomes more and more disordered when $r>0.16$, which
represent the region close to the meeting point $T$ of the two second-order
critical lines on the phase diagram. Similar phenomena can also be observed on
blue paths.

Figure 5:  (a) (b) DOS of $\tilde{Z}$ for $g=0.19$ and $r\in[0.01,0.19]$ on
red paths. (c) DOS of $\tilde{Z}$ for $r=0.14$. (d) DOS of $\tilde{Z}$ for
$r=0.18$.

A characteristic feature of the first-order phase transition is the
coexistence of different phases. The conversion from single-maximum to multi-
maxima of the DOS of $\tilde{Z}$ is then regarded as the conversion from
second-order phase transition to first-order phase transition.

As our system is small, although the multi-maxima feature appears in DOS when
$r>0.16$ in red paths, the first-order phase transition may not actually
appear when $0.16<g<0.24$. It is likely that it appears only after the two
second-order lines cross at $g\approx 0.24$.

Besides, no matter what value $p$ we choose for the blue path, we can always
observe an extreme point of $\mathrm{d}^{2}E/dg^{2}$ right on the dual line,
as shown in Fig. 7. Yet, an extreme point does not have to be a critical
point. To determine the end of the first-order phase transition line outside
the deconfined phase, we investigate the DOS of $\tilde{Z}$.

Figure 6:  DOS of $\tilde{Z}$ on the dual line. The disorder recedes as $g$
increase. (a) At $g\approx 0.30$, the multi-maxima is reduced to two maxima.
(b) At $g\approx 0.65$, we detect a single maximum.

Since the Higgs phase and confined phase are continuously connected when $g$,
$\lambda\to\infty$ PhysRevD.19.3682 , the first-order line should vanish
somewhere. Fig. 6 shows the variation of the DOS on the dual line and gives
two special points, $g_{1}=0.30$ and $g_{2}=0.65$. We consider $g_{1}$ as the
end of the line because the maxima at $\tilde{Z}=-10$ indicates that the
system starts to mainly permit the excitation of the smallest loop, which is
composed of four spins, and other possible configurations are generally
disfavored. In other word, the system turns into confinement.

Figure 7:  Energy as a function of $\lambda$, for blue paths with different
values of $p$. (a) $p=0.36$; (b) $p=0.32$; (c) $p=0.28$; (d) $p=0.24$. Figure
8: The second derivative of energy on the dual line, with an extreme point at
$g=0.1798<0.24$.

Besides, we should observe three extreme points when $0.19<p<0.48$ with blue
paths in principle according to Fig. 1. However, such kind of diagram can only
be clearly observed when $p<0.28$. For $p>0.28$, we observe a single extreme
point on the dual line. We regard this phenomenon as a limitation of the size
of the lattice, which makes the extreme point on the dual line
indistinguishable with the two on both sides.

The extreme points on the dual line make it natural for us to see the energy
curve right on the dual line, so we choose a new path from $(0,0.5)$, through
$(0.5,0.5)$, and to $O(0,0)$. Surprisingly we find an extreme point at
$g_{3}=0.1798$, which is in the deconfined phase (See Fig. 8). As the DOS is
also disordered in the corner of deconfined phase, we suspect that the first-
order line will extend for a while in deconfined phase and ends at $g_{3}$
(See Fig. 1).

The phase diagram we obtain is similar to the diagram in Vidal:2008uy , which
is derived by perturbation theory. The difference lies in the extension of the
first-order line in deconfined phase. The region near the tricritical point
has not been clear. Our result corresponds to one of the possibilities
proposed in PhysRevB.82.085114 .

##  6 Error analysis

As the Trotter-Suzuki decomposition is used in each step, the upper bound of
the total error should be a summation of the errors in each step.

We take the sub-path from $R(p,r)=(0.5,0.05)$ to $R^{\prime}(0,r)=(0,0.05)$ as
an example. As derived in Appendix, the accumulated Trotter error $\epsilon$
is

$\begin{split}\epsilon\lesssim&{(\Delta
t)^{3}}N(13.5p+54pr^{2}+54pr+54r^{2}\\\ &+9p^{2}+18p^{2}r+13.5r),\\\
\end{split}$ (16)

where $\Delta t$ is the time for each step of varying the parameter values,
$N$ is the number of steps.

For variation of $g$ from $0$, $\Delta g$ is the variation in each step.
$\Delta g\propto\Delta t$. We choose $\Delta g=10^{-5}$ and $\Delta t=2\times
10^{-2}$, giving $\epsilon\approx 4.581$, which seems to be too large for us.
However, this error is achieved by adding the positive errors in all steps,
the real errors are not always positive and may cancel with each other
PRXQuantum.2.010323 . We demonstrate it in the following.

Rewrite $\Delta g=0.0001/n$ and $\Delta t=0.2/n$, then $n=n_{0}=10$ under the
parameter values above. Define

$\begin{split}&\delta_{n}=\sum_{g=0}^{p}|E_{n_{0}}(g)-E_{n}(g)|\\\
&\alpha_{n}=\max_{g\in[0,p]}|E_{n_{0}}(g)-E_{n}(g)|\\\ \end{split}$ (17)

where $\delta_{n}$ is a 1-norm of the difference between two energy functions
$E_{n_{0}}(g)$ and $E_{n}(g)$, $\alpha_{n}$ represents the largest value of
the difference. The results in TALBLE 1 and Fig. 9 show that even
$\varepsilon$ at $n_{0}=10$ is around 6363 times more than that at $n=800$,
$\delta_{n}$ is merely $0.04075$ and the maximal deviation $\alpha_{n}$ is
only $0.005\%$. This is a direct evidence that the error at each step cannot
always be positive, thus our choice of $n=10$ is feasible.

$n$ | $\Delta g$ | $\Delta t$ | $\varepsilon$ | $\delta_{n}$ | $\alpha_{n}$
---|---|---|---|---|---
10 | $1.00\times 10^{-5}$ | $2.0\times 10^{-2}$ | 4.58100 | 0 | 0
200 | $5.00\times 10^{-7}$ | $1.0\times 10^{-3}$ | 0.01145 | 0.03982 | $4.79\times 10^{-5}$
500 | $2.00\times 10^{-7}$ | $4.0\times 10^{-4}$ | 0.00183 | 0.04009 | $4.99\times 10^{-5}$
800 | $1.25\times 10^{-7}$ | $2.5\times 10^{-4}$ | 0.00072 | 0.04075 | $4.98\times 10^{-5}$
Table 1: Different values of $\delta_{n}$ and $\alpha_{n}$ for
$n=10,200,500,800$. Figure 9: The diagram of difference
$E_{n_{0}}(g)-E_{n}(g)$ for different values of $n$.

##  7 Summary

We have first designed the quantum simulation of $\mathbb{Z}_{2}$ gauge Higgs
model. The quantum simulation is digital, as it is based on Trotter
decomposition of the unitary evolution. It is also based on quantum adiabatic
algorithm. For each parameter value, which is varied slowly, the Trotter
decomposition is used in execute the unitary transformation. Within each step
in the Trotter decomposition, the unitary transformations are realized in
terms of simple quantum circuits.

Then we make a thorough classical demonstration of this quantum simulation by
using QuEST simulator in a NVIDIA GeForce RTX 3090 GPU server. We call such an
approach pseudoquantum simulation. We haved perform a psudoquantum simulation
of $\mathbb{Z}_{2}$ gauge Higgs model on a $3\times 3$ lattice. The lattice is
small, as limited by computational time.

We have obtained some clear results on the topological properties of
deconfined phase. Especially, our work suggests that the two lines of second-
order transitions meet on the line of the first-order transition rather than
its end.

##  Appendix A Calculation of the Trotter errors

It is useful to calculate $\lVert H_{j}\rVert$ and $\lVert[H_{k},H_{l}]\rVert$
in advance where $\lVert\cdot\rVert$ denotes the spectral norm of an operator,
or the largest singular value of it. We use the estimation that $\lVert
O\rVert\sim 1$, where $O$ is a Pauli operator or a tensor product of Pauli
operators.

Figure 10: Lattice used in the simulations.

For our specific $3\times 3$ lattice model (See Fig. 10), we have

$\begin{split}&\lVert H_{1}\rVert\sim 18g,\quad\lVert H_{3}\rVert\sim
18\lambda\\\ &\lVert
H_{2}\rVert\leq\frac{1}{2}\lVert\sum_{p}B_{p}^{z}\rVert+\frac{1}{2}\lVert\sum_{v}A^{x}_{v}\rVert\sim
9\end{split}$ (18)

and

$\begin{split}&\lVert[H_{1},H_{2}]\rVert\\\
=&\lVert[-g\sum_{l}\sigma^{x}_{l},-\frac{1}{2}\sum_{p}B_{p}^{z}-\frac{1}{2}\sum_{v}A_{v}^{x}]\rVert\\\
=&\frac{g}{2}\lVert[\sum_{l}\sigma^{x}_{l},\sum_{p}B_{p}^{z}]\rVert\\\
\leq&\frac{g}{2}\sum_{l}\lVert[\sigma^{x}_{l},\sum_{p}B_{p}^{z}]\rVert\\\
=&18\frac{g}{2}\lVert[\sigma^{x}_{4},\sum_{p}B_{p}^{z}]\rVert\\\
=&18\frac{g}{2}\lVert[\sigma^{x}_{4},\sigma^{z}_{1}\sigma^{z}_{10}\sigma^{z}_{4}\sigma^{z}_{13}]+[\sigma^{x}_{4},\sigma^{z}_{4}\sigma^{z}_{11}\sigma^{z}_{7}\sigma^{z}_{14}]\rVert\\\
\sim&18{g}\end{split}$ (19)

where the coefficient $18$ comes up due to the symmetry of the lattice.
Similarly, we have $\lVert[H_{2},H_{3}]\rVert\sim 18\lambda$ and

$\begin{split}&\lVert[H_{1},H_{3}]\rVert\\\
=&\lVert[-g\sum_{l}\sigma^{x}_{l},-\lambda\sum_{s}\sigma_{s}^{z}]\rVert\\\
=&g\lambda\lVert[\sum_{l}\sigma^{x}_{l},\sum_{s}\sigma_{s}^{z}]\rVert\\\
=&g\lambda\lVert\sum_{l}[\sigma_{l}^{x},\sigma^{z}_{l}]\rVert\\\
\sim&36g\lambda\end{split}$ (20)

Multiplying a minus sign does not change the eigenvalues of a matrix, thus
$\lVert[H_{2},H_{1}]\rVert\sim 18g$, $\lVert[H_{3},H_{1}]\rVert\sim
36g\lambda$ and $\lVert[H_{3},H_{2}]\rVert\sim 18\lambda$.

Suppose $H=\sum_{\gamma=1}^{\Gamma}H_{\gamma}$, then the tight error bound for
the second-order decomposition is PhysRevX.11.011020

$\begin{split}\lVert S_{2}(\Delta t)-e^{-i(\Delta t)H}\rVert\leq&\frac{(\Delta
t)^{3}}{12}\sum_{\gamma_{1}=1}^{\Gamma}\lVert[\sum_{\gamma_{3}=\gamma_{1}+1}^{\Gamma}H_{\gamma_{3}},[\sum_{\gamma_{2}=\gamma_{1}+1}^{\Gamma}H_{\gamma_{2}},H_{\gamma_{1}}]]\rVert\\\
&+\frac{(\Delta
t)^{3}}{24}\sum_{\gamma_{1}=1}^{\Gamma}\lVert[H_{\gamma_{1}},[H_{\gamma_{1}},\sum_{\gamma_{2}=\gamma_{1}+1}^{\Gamma}H_{\gamma_{2}}]]\rVert\\\
\end{split}$ (21)

Our Hamiltonian is $H=H_{1}+H_{2}+H_{3}$, thus the error bound is

$\begin{split}\lVert S_{2}(t)-e^{-itH}\rVert\leq&\frac{(\Delta
t)^{3}}{12}\sum_{l=1}^{3}\lVert[\sum_{n=l+1}^{3}H_{n},[\sum_{m=l+1}^{3}H_{m},H_{l}]]\rVert\\\
&+\frac{(\Delta
t)^{3}}{24}\sum_{l=1}^{3}\lVert[H_{l},[H_{l},\sum_{m=l+1}^{3}H_{m}]]\rVert\\\
\leq&\frac{(\Delta
t)^{3}}{12}\sum_{l=1}^{3}\sum_{n=l+1}^{3}\sum_{m=l+1}^{3}\lVert
H_{n}[H_{m},H_{l}]-[H_{m},H_{l}]H_{n}\rVert\\\ &+\frac{(\Delta
t)^{3}}{24}\sum_{l=1}^{3}\sum_{m=l+1}^{3}\lVert
H_{l}[H_{l},H_{m}]-[H_{l},H_{m}]H_{l}\rVert\\\ \leq&\frac{(\Delta
t)^{3}}{6}\sum_{l=1}^{3}\sum_{n=l+1}^{3}\sum_{m=l+1}^{3}\lVert
H_{n}\rVert\lVert[H_{m},H_{l}]\rVert\\\ &+\frac{(\Delta
t)^{3}}{12}\sum_{l=1}^{3}\sum_{m=l+1}^{3}\lVert
H_{l}\rVert\lVert[H_{l},H_{m}\rVert\\\ =&\frac{(\Delta
t)^{3}}{6}W_{1}+\frac{(\Delta t)^{3}}{12}W_{2}\end{split}$ (22)

where

$\begin{split}&W_{1}=\sum_{l=1}^{3}\sum_{n=l+1}^{3}\sum_{m=l+1}^{3}\lVert
H_{n}\rVert\lVert[H_{m},H_{l}]\rVert\\\
&W_{2}=\sum_{l=1}^{3}\sum_{m=l+1}^{3}\lVert
H_{l}\rVert\lVert[H_{l},H_{m}]\rVert\\\ \end{split}$ (23)

According to (18), we have

$\begin{split}W_{1}=&\sum_{n=2}^{3}\sum_{m=2}^{3}\lVert
H_{n}\rVert\lVert[H_{m},H_{1}]\rVert+\lVert
H_{3}\rVert\lVert[H_{3},H_{2}]\rVert\\\ =&\lVert
H_{2}\rVert(\lVert[H_{2},H_{1}]\rVert+\rVert\lVert[H_{3},H_{1}]\rVert)\\\
&+\lVert
H_{3}\rVert(\lVert[H_{2},H_{1}]\rVert+\lVert[H_{3},H_{1}]\rVert+\lVert[H_{3},H_{2}]\rVert)\\\
\sim&9(18{g}+36g\lambda)+18\lambda(18{g}+36g\lambda+18{\lambda})\\\
=&162({g}+2g\lambda)+324\lambda({g}+2g\lambda+{\lambda})\\\
=&162g+648g\lambda^{2}+648g\lambda+324\lambda^{2}\end{split}$ (24)

and

$\begin{split}W_{2}=&\sum_{m=2}^{3}\lVert
H_{1}\rVert\lVert[H_{1},H_{m}\rVert+\lVert
H_{2}\rVert\lVert[H_{2},H_{3}\rVert\\\ =&\lVert
H_{1}\rVert(\lVert[H_{1},H_{2}\rVert+\lVert[H_{1},H_{3}\rVert)+\lVert
H_{2}\rVert\lVert[H_{2},H_{3}\rVert\\\
\sim&18g(18g+36g\lambda)+9(18{\lambda})\\\
=&324g^{2}+648g^{2}\lambda+162\lambda\end{split}$ (25)

The evolution in Sec. 6 is from $(0,r)$ to $(p,r)$ on $l:\lambda=r$ with $g$
increasing from $0$ to $p$ in $N=p/\Delta g$ steps. Denote $g=n\Delta g$, then
the total error is

$\begin{split}\epsilon\leq&\frac{(\Delta
t)^{3}}{6}\sum^{N}_{n=1}W_{1,n}+\frac{(\Delta
t)^{3}}{12}\sum_{m=1}^{N}\Delta_{m}W_{2,m}\\\ \sim&\frac{(\Delta
t)^{3}}{6}\sum^{N}_{n=1}(162g+648g\lambda^{2}+648g\lambda+324\lambda^{2})\\\
&+\frac{(\Delta
t)^{3}}{12}\sum_{m=1}^{N}(324g^{2}+648g^{2}\lambda+162\lambda)\\\ =&{(\Delta
t)^{3}}\sum^{N}_{n=1}(27g+108g\lambda^{2}+108g\lambda+54\lambda^{2})\\\
&+{(\Delta t)^{3}}\sum_{m=1}^{N}(27g^{2}+54g^{2}\lambda+13.5\lambda)\\\
=&{(\Delta
t)^{3}}\sum^{N}_{n=1}(27g+108g\lambda^{2}+108g\lambda+54\lambda^{2}\\\
&+27g^{2}+54g^{2}\lambda+13.5\lambda)\\\ =&{(\Delta
t)^{3}}\sum^{N}_{n=1}(27g+108gr^{2}+108gr+54r^{2}\\\
&+27g^{2}+54g^{2}r+13.5r)\\\ \end{split}$ (26)

Besides, we have

$\begin{split}&\sum_{n=1}^{N}g=\sum_{n=1}^{N}(n\Delta g)\approx\Delta
g\frac{N^{2}}{2}=\frac{1}{2}pN\\\ &\sum_{n=1}^{N}g^{2}=\sum_{n=1}^{N}(n\Delta
g)^{2}\approx\Delta g^{2}\frac{N^{3}}{3}=\frac{1}{3}p^{2}N\\\ \end{split}$
(27)

Rewrite (26) with (LABEL:g^1), we have

$\begin{split}\epsilon\lesssim&{(\Delta
t)^{3}}(27\frac{1}{2}pN+108\frac{1}{2}pNr^{2}+108\frac{1}{2}pNr+54r^{2}N\\\
&+54\frac{1}{3}p^{2}N+54\frac{1}{3}p^{2}Nr+13.5rN)\\\ =&{(\Delta
t)^{3}}N(13.5p+54pr^{2}+54pr+54r^{2}\\\ &+9p^{2}+18p^{2}r+13.5r)\\\
\end{split}$ (28)

## Acknowledgments

This work was supported by National Science Foundation of China (Grant No.
12075059).

## References

  * [1] K. G. Wilson.  Confinement of quarks.  Phys. Rev. D 10 (1974) 2445.
  * [2] E. A. Martinez et al.  Real-time dynamics of lattice gauge theories with a few-qubit quantum computer.  Nature 534 (2016) 516.
  * [3] E. Zohar, D. González-Cuadra, and J. I. Cirac.  Quantum simulation of the abelian-higgs lattice gauge theory with ultracold atoms.  New Journal of Physics 19 (2017) 063038.
  * [4] F. M. Surace et al.  Lattice gauge theories and string dynamics in rydberg atom quantum simulators.  Phys. Rev. X. 10 (2020) 021041.
  * [5] X. Cui, Y. Shi, and J. Yang.  Circuit-based digital adiabatic quantum simulation and pseudoquantum simulation as new approaches to lattice gauge theory.  J. High Energ. Phys. 160 (2020) .
  * [6] X. Cui and Y. Shi.  Trotter errors in digital adiabatic quantum simulation of quantum $\mathbb{Z}_{2}$ lattice gauge theory.  International Journal of Modern Physics B 34 (2020) 2050292.
  * [7] I. S. Tupitsyn, A. Kitaev, N. V. Prokof’ev, and P. C. E. Stamp.  Topological multicritical point in the phase diagram of the toric code model and three-dimensional lattice gauge higgs model.  Phys. Rev. B 82 (2010) 085114.
  * [8] E. Fradkin and S. H. Shenker.  Phase diagrams of lattice gauge theories with higgs fields.  Phys. Rev. D 19 (1979) 3682.
  * [9] J. L. Banks and D. K. Sinclair.  Hamiltonian Calculations for the Z2 Lattice Gauge Theory With Matter Fields.  Phys. Rev. D 23 (1981) 2962.
  * [10] Z. Nussinov.  A Derivation of the Fradkin-Shenker result from duality: Links to spin systems in external magnetic fields and percolation crossovers.  Phys. Rev. D 72 (2005) 054509.
  * [11] J. Vidal, S. Dusuel, and K. P. Schmidt.  Low-energy effective theory of the toric code model in a parallel field.  Phys. Rev. B 79 (2009) 033109.
  * [12] G. A. Jongeward, J. D. Stack, and C. Jayaprakash.  Monte Carlo calculations on $\mathbb{Z_{2}}$ gauge-Higgs theories.  Phys. Rev. D 21 (1980) 3360.
  * [13] M. Creutz.  Phase diagrams for coupled spin-gauge systems.  Phys. Rev. D 21 (1980) 1006.
  * [14] A. M. Somoza, P. Serna, and A. Nahum.  Self-dual criticality in three-dimensional $\mathbb{Z}_{2}$ gauge theory with matter.  arXiv:2012.15845.
  * [15] N. Hatano and M. Suzuki.  Finding Exponential Product Formulas of Higher Orders.  Springer Berlin Heidelberg, Berlin, Heidelberg (2005).
  * [16] M. C. Tran, Y. Su, D. Carney, and J. M. Taylor.  Faster digital quantum simulation by symmetry protection.  PRX Quantum 2 (2021) 010323.
  * [17] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. C. Zhu.  Theory of trotter error with commutator scaling.  Phys. Rev. X 11 (2021).


2022

These authors contributed equally to this work.

These authors contributed equally to this work.

[2,3]Haitan Xu

[1,4,5]Zheng Li

1]State Key Laboratory for Mesoscopic Physics and Collaborative Innovation
Center of Quantum Matter, School of Physics, Peking University, Beijing,
100871, China

2]School of Materials Science and Intelligent Engineering, Nanjing University,
Suzhou, 215163, Jiangsu, China 3]Shenzhen Institute for Quantum Science and
Engineering, Southern University of Science and Technology, Shenzhen, 518055,
Guangdong, China

4]Peking University Yangtze Delta Institute of Optoelectronics, Nantong, China

5]Collaborative Innovation Center of Extreme Optics, Shanxi University, ,
Taiyuan, 030006, Shanxi, China

6]Physics Department, Stanford University, CA, 94305, USA

# Entangled X-ray Photon Pair Generation by Free Electron Lasers

Linfeng Zhang  Zunqi Li  Dongyu Liu  Chengyin Wu<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [  [

Einstein, Podolsky and Rosen’s prediction einstein1935can  on incompleteness
of quantum mechanics was overturned by experimental tests on Bell’s inequality
bell1966problem  that confirmed the existence of quantum entanglement. In
X-ray optics, entangled photon pairs can be generated by X-ray parametric down
conversion (XPDC), which has certain wavelength window shwartz2012x .
Meanwhile, free electron laser (FEL) has successfully lased at X-ray
frequencies recently emma2010first ; ishikawa2012compact ; allaria2013two ;
wang2021free ; pompili2022free ; duris2020tunable ; duris2021controllable .
However, FEL is usually seen as a classical light source, and its quantum
effects are considered minor corrections to the classical theory. Here we
investigate entangled X-ray photon pair emissions in FEL. We establish a
theory for coherently amplified entangled photon pair emission from
microbunched electron pulses in the undulator. We numerically demonstrate the
properties of entangled emission, and provide a scheme to generate highly
entangled X-ray photon pairs, which is of great importance in X-ray quantum
optics. Our work shows a unique advantage of FELs over synchrotrons in
entangled X-ray photon pair generation.

Being an ultrabright and ultrashort X-ray laser source, FEL has shed light on
various research areas including single particle imaging seibert2011single ;
loh2012fractal ; gomez2014shapes , ultrafast X-ray spectroscopy
kern2013simultaneous ; erk2014imaging ; rudek2012ultra  and high energy
density physics vinko2012creation ; fletcher2015ultrabright ;
ciricosta2016measurements . The FEL emission was mostly seen as a classical
phenomenon with minor quantum corrections pellegrini2016physics . However, a
critical quantum effect, i.e., the quantum entanglement, has never been
investigated for the X-ray emission from the FEL.

In this work, we study the entangled X-ray photon pair emission from FEL,
which can be enhanced by electron microbunching in the undulator. We start by
presenting the quantum electrodynamics (QED) treatment of single-electron
photon pair emission process in the undulator and calculate the cross section
and the entanglement degree. Then we apply the Feynman rules to microbunched
electrons, which takes into account the many-body effect, to analyze the
coherent amplification and establish the condition for coherence. We further
numerically investigate the enhancement of entangled photon pair emission by
microbunched electrons.

## Entangled photons from FEL

Figure 1: Mechanism of coherent entangled X-ray photon pair emission in SASE
FEL. a, Schematic in the lab frame, where a relativistic electron pulse (red)
passes through an undulator, causing photon emission (yellow) while being
microbunched. The fundamental frequency of the emission can reach X-ray
region, with single-photon energy of hundreds of $\mathrm{eV}$s. The schematic
of the helical undulator is drawn according to the configuration of the Delta
undulator at LCLS lutman2016polarization . b, Schematic in the electron frame,
where the electron pulse is initially at rest. The strong quasi-EM wave (blue)
of the relativistic undulator is scattered from the electrons. The resulting
emission contains multiple components, including single-photon emission,
double-photon emission, etc. c, The Feynman diagrams corresponding to double-
photon emission in the electron frame. The black lines refer to the Volkov
states of electron. Inset: a typical constraint of four-momenta of the
particles. $n$ is the net number of Floquet photons absorbed. d, Coherent
summation of Feynman diagrams with $n=1$, including permutations of $k_{1}$
and $k_{2}$.

We present in Fig. 1 the quantum mechanism of entangled X-ray photon pair
emission. Fig. 1a shows a typical SASE FEL process in the lab frame. A
relativistic electron pulse travels through an undulator, interacts with the
electromagnetic (EM) field of the undulator, resulting in electron
microbunching and amplified emission.

To facilitate the study of entangled photon emission from the FEL, we use the
electron frame (EF), where the incident electron pulse has a zero initial
average velocity. We adopt the frame transformation formalism from the
Weizsäcker-Williams method madey1971stimulated , and use the natural units
$\hbar=c=\epsilon_{0}=1$. For a helical undulator,
$A_{\mathrm{Lab}}^{\mu}=a(0,\mathrm{cos}(k_{u}x_{\mathrm{Lab}}^{3}),\mathrm{sin}(k_{u}x_{\mathrm{Lab}}^{3}),0)$,
where $a=B_{0}/k_{u}$, $B_{0}$ is the amplitude of the undulator field, and
$k_{u}=2\pi/\lambda_{u}$ is the spatial frequency of the undulator with a
period of $\lambda_{u}$ in the lab frame temnykh2008delta . The undulator
field is transformed to a quasi-EM wave in the EF by the Lorentz boost as

$\begin{split}A_{\mathrm{EF}}^{\mu}&=a(0,\mathrm{cos}[k_{u}(\gamma\beta
x_{\mathrm{EF}}^{0}+\gamma x_{\mathrm{EF}}^{3})],\\\
&\phantom{=}\mathrm{sin}[k_{u}(\gamma\beta x_{\mathrm{EF}}^{0}+\gamma
x_{\mathrm{EF}}^{3})],0)\\\ &=a(0,\mathrm{cos}(k\cdot
x_{\mathrm{EF}}),\mathrm{sin}(k\cdot x_{\mathrm{EF}}),0)\;,\end{split}$ (1)

where $\beta$ and $\gamma$ are the average velocity and Lorentz factor of the
electron pulse. $k=\gamma k_{u}(\beta,0,0,-1)$ is the wave vector of the
undulator quasi-EM wave in the EF, which is close to the photon mass shell for
$\beta\to 1$. The helical undulator is practically available, such as the
Delta undulator at LCLS lutman2016polarization .

The SASE emission can be approximated as the scattering of the quasi-EM wave
from the electron pulse in the EF. The emission takes part in the
microbunching of the electron pulse, and is significantly enhanced at the
fundamental frequency and its harmonics. In the lab frame, the fundamental
frequency $\omega_{\mathrm{fd}}=2\gamma^{2}(2\pi/\lambda_{u})/(1+K^{2})$ for
the helical undulator, where $K=eB_{0}/mk_{u}=ea/m$ is the undulator parameter
emma2010first ; pellegrini2016physics ; huang2007review .

For strong undulator quasi-EM wave, nonlinear scattering can be observed in
the emission. We focus on the double-scattering processes, as shown in Fig.
1b. In this work we use the FEL parameters of Linac Coherent Light Source
(LCLS) emma2010first , as shown in Table 1.

For a representative electron initially at rest in the EF, it satisfies the
Dirac equation $(\not{p}-e\not{A}-m)\Psi=0$, where
$\not{f}=\gamma^{\mu}f_{\mu}$, $\gamma^{\mu}$ are the Dirac matrices, $A$ is
the EM four-potential, $e$ is the electron charge, and $m$ is the rest
electron mass.

Electron energy | $5~{}\mathrm{GeV}$ | Undulator period number | $1392$
---|---|---|---
Charge per pulse | $180~{}\mathrm{pC}$ | Pierce parameter | $2\times 10^{-3}$
Peak current | $3.4~{}\mathrm{kA}$ | Fundamental wavelength | $2.30~{}\mathrm{nm}$
Repetition rate | $120~{}\mathrm{Hz}$ | Saturation peak power | $10~{}\mathrm{GW}$
Undulator peak field | $1.32~{}\mathrm{T}$ | Emission pulse FWHM | $230~{}\mathrm{fs}$
Undulator period length | $30~{}\mathrm{mm}$ |  |
Table 1: FEL parameters used in this work, based on LCLS emma2010first .

In the undulator quasi-EM wave, the electron is in a Volkov state
$\Psi^{\mathrm{V}}_{p}(x)=(1-\frac{e\not{A}\not{k}}{2k\cdot
p})u_{p}\times\mathrm{exp}\\{-i[p\cdot x+\int_{0}^{\phi}{(\frac{ep\cdot
A(\phi^{\prime})}{k\cdot p}-\frac{e^{2}A^{2}(\phi^{\prime})}{2k\cdot
p})d\phi^{\prime}}]\\}$, where $\phi=k\cdot
x=k^{\mu}x_{\mu}=k^{0}x^{0}-\boldsymbol{k}\cdot\boldsymbol{x}$, $u_{p}$ is the
free electron spinor, and $A\simeq A_{\mathrm{EF}}$ [see Section I of
Supplementary Information (SI) for details]. For an electron of initial four-
momentum $p$, the corresponding Volkov state has a time-varying four-momentum
whose average is the quasi-momentum $q^{\mu}=p^{\mu}+\frac{e^{2}a^{2}}{2k\cdot
p}k^{\mu}$, satisfying $|q|^{2}=m^{2}+e^{2}a^{2}=m_{*}^{2}$.

In the Floquet picture, the Volkov state corresponds to the superposition
state of the electron dressed by $\tilde{n}$ Floquet photons, i.e.,
$\Psi^{\mathrm{V}}_{p}(x)=\sum_{\tilde{n}=-\infty}^{+\infty}{F_{\tilde{n}}\mathrm{exp}\\{-i(q\cdot
x+\tilde{n}k\cdot x)}\\}u(p)$ (the expression of $F_{\tilde{n}}$ can be found
in SI). The electron in the Volkov state can spontaneously decay into another
Volkov state with a different quasi-momentum and thus radiate photons, the
amplitude of which can be evaluated by the Feynman diagrams (see details in
Section II of SI). Especially, the entangled photon pair emission is
equivalent to non-perturbative double Compton scattering of the quasi-EM wave
from the electron in the EF ritus1985quantum ; lotstedt2009nonperturbative ,
which satisfies the four-momentum conservation
$q_{\mathrm{f}}+k_{1}+k_{2}-q_{\mathrm{i}}-nk=0$ (see Fig. 1c). $k_{1}$ and
$k_{2}$ are the momentum of entangled photons, and $n$ corresponds to the
number of Floquet photons absorbed. As the energies of photons
($k^{0},k_{1}^{0},k_{2}^{0}$) are much less than the electron mass ($m$) and
the quasi-momenta ($q_{\mathrm{i}},q_{\mathrm{f}}$) in the EF when $n$ is
reasonably small, the four-momentum conservation leads to
$k_{1}^{0}+k_{2}^{0}\simeq nk^{0}$. We focus on the representative case of
$n=1$, which corresponds to photon pairs with total frequency of
$\omega_{\mathrm{fd}}$ in the lab frame, and can be distinguished from other
cases with different $n$ by detecting the emitted photon pair energies (see
details in Section III of SI).

With the above treatment, we obtain the scattering matrix element
$S_{\mathrm{fi}}$ for the process from an initial Volkov state $\mathrm{i}$ to
any final Volkov state $\mathrm{f}$ emitting a pair of photons in plane wave
states with arbitrary polarization or helicity. The differential double-photon
emission rate is given by

$\begin{split}&\frac{d\dot{W}}{dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}}=\frac{e^{4}}{(2\pi)^{5}}\frac{m^{2}}{2q_{\mathrm{i}}^{0}}\\\
&\phantom{=}\times\sum_{n=1}^{+\infty}{\frac{k_{1}^{0}(k_{2}^{0})^{2}\Theta(q_{\mathrm{i}}^{0}+nk^{0}-k_{1}^{0}-k_{2}^{0})M}{2|(q_{\mathrm{i}}+nk-
k_{1})\cdot k_{2}|}}\bigg{|}_{\mathrm{DE}(n)}\;,\end{split}$ (2)

where $\Omega_{k_{1}}$ and $\Omega_{k_{2}}$ are the solid angles of the
emitted photon pair, the notation $|_{\mathrm{DE}(n)}$ represents the
constraint of $(q_{\mathrm{i}}+nk-k_{1})\cdot k_{2}=q_{\mathrm{i}}\cdot nk-
q_{\mathrm{i}}\cdot k_{1}-nk\cdot k_{1}$, and the helicity-relevant part $M$
is explained in SI (see details in Methods and Section II of SI).

In order to quantify the entanglement degree of the emitted X-ray photon
pairs, we calculate the density matrix $\rho_{\mathrm{f}}$ in the helicity
basis from the scattering matrix elements $S_{\mathrm{fi}}$ as
$\rho_{\mathrm{f}}=\frac{1}{2}\sum_{r_{\mathrm{i}},r_{\mathrm{f}}=1}^{2}N\big{(}S_{\mathrm{fi},j_{1}}S_{\mathrm{fi},j_{2}}^{*}\big{)}_{4\times
4}$, where $S_{\mathrm{fi},j_{1}}$ and $S_{\mathrm{fi},j_{2}}$ represent the
scattering matrix elements with the emitted photon pair at helicity
eigenstates $j_{1}$ and $j_{2}$, and $j_{1},j_{2}=1,2,3,4$ correspond to the
helicity eigenstates $|++\rangle,|+-\rangle,|-+\rangle$, and $|--\rangle$,
where $\pm$ corresponds to the right/left-handed helicity
lotstedt2009correlated . $N$ is chosen to keep
$\mathrm{Tr}[\rho_{\mathrm{f}}]=1$. $r_{\mathrm{i}}$ and $r_{\mathrm{f}}$
represent the spin of the electron at initial and final states, respectively,
which have been traced out. The scattering matrix elements can be calculated
for any choice of $k_{1},k_{2}$ satisfying the constraint
$(k_{1}+k_{2}-q_{\mathrm{i}}-nk)^{2}-m_{*}^{2}=0$.

In order to characterize the quantum entanglement of the emitted photon pair,
we use the so-called concurrence as a measure of the entanglement
wootters1998entanglement . The concurrence can be calculated from the density
matrix as

$\mathcal{C}(\rho_{\mathrm{f}})=\mathrm{max}(0,\sqrt{\zeta_{1}}-\sqrt{\zeta_{2}}-\sqrt{\zeta_{3}}-\sqrt{\zeta_{4}})\;,$
(3)

where $\zeta_{1,2,3,4}$ are the 4 eigenvalues of the matrix
$Q=\rho_{\mathrm{f}}(\sigma^{2}\otimes\sigma^{2})\rho_{\mathrm{f}}^{*}(\sigma^{2}\otimes\sigma^{2})$
in descending order, and $\sigma^{2}$ is one of the Pauli matrices.
$\mathcal{C}\in[0,1]$, with $\mathcal{C}=0$ corresponding to a unentangled
state and $\mathcal{C}=1$ corresponding to a fully entangled state. As the
helicity of massless particles like photons are Lorentz invariant, the density
matrix and the concurrence are also Lorentz invariant.

We have represented the density matrix and the concurrence of the 2-qubit
states of the emitted photon pairs with the basis of the helicity eigenstates
$|++\rangle,|+-\rangle,|-+\rangle,|--\rangle$. The choice of basis is not
unique, and we can alternatively choose other basis states, such as linear
polarization eigenstates lotstedt2009correlated , which result in a unitary
transformation of the density matrix, and the concurrence remains invariant
(see Section VII of SI for details). In addition, we can choose other
entanglement monotones instead of the concurrence to quantify the entanglement
degree, such as entanglement of formation and negativity. A comparison showing
qualitative agreement between these entanglement measures and the concurrence
is given in Section IX of SI.

## Microbunch enhancement

In the high-gain regime of FEL, a coherent enhancement due to the electron
microbunching plays a vital role. For a sufficiently bright electron beam and
a sufficiently long undulator, the combination of the undulator field and the
radiation field will induce a ponderomotive potential modulating the energy of
the electron beam huang2007review . This energy modulation eventually forces
the electrons to form periodic microbunches along the undulator axis with a
modulation wavelength $\lambda_{1}=2\pi/\omega_{\mathrm{fd}}$ equal to the
fundamental wavelength.

The conventional analysis of the coherence property due to the microbunching
relies on paraxial wave equation huang2007review , which cannot explain the
quantum phenomenon of entangled photon pair emission. To overcome this
problem, we develop a quantum collective double emission (QCDE) theory (see
details in Section IV of SI). This theory provides a QED description of the
microbunched state of electrons, and deals with the interference of the
Feynman diagrams corresponding to double emissions from different electrons.
The collective double emission rate in the QCDE theory is

$\displaystyle\bigg{(}\frac{d\dot{W}}{dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}}\bigg{)}_{\mathrm{C}}=\frac{d\dot{W}}{dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}}F_{\mathrm{MB}}(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})\;,$
(4)

where $F_{\mathrm{MB}}(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})$ is the
enhancement factor on the scattering rate introduced by microbunching.
$(d\dot{W}/dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}})_{\mathrm{C}}$ represents
the collective differential scattering rate for a microbunched electron pulse.
In addition to the rigorous QCDE theory presented in SI, here we elaborate the
physical picture leading to the coherent enhancement factor
$F_{\mathrm{MB}}(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})$ via an analogy with
the quasi-phase-matching in nonlinear optics (see details in Section IV of
SI).

In the electron frame, we can view the quasi-EM wave generated by the
relativistic undulator as the pump and the emitted photon pair as signal and
idler. The electron microbunches can be regarded as a nonlinear medium with
periodic structure. Therefore, the phase difference between photons emitted by
adjacent microbunches should be equal to
$\Delta\phi=\Delta\boldsymbol{l}\cdot\Delta\boldsymbol{k}$, where
$\Delta\boldsymbol{l}=|\Delta\boldsymbol{l}|\hat{\boldsymbol{e}}^{3}=(\lambda_{u}/\gamma)(1+K^{2})/(2+K^{2})|\hat{\boldsymbol{e}}^{3}$
is the average relative displacement between adjacent microbunches, with
$\hat{\boldsymbol{e}}^{3}$ being the unit vector along the undulator, and
$\Delta\boldsymbol{k}=\boldsymbol{k_{1}}+\boldsymbol{k_{2}}-n\boldsymbol{k}$
is the momentum transfer between the pump photon and the emitted photon pair.

If the momenta of the signal and idler photons are both along the $x^{3}$
axis, the total momentum of the photon pair is
$\boldsymbol{k_{1}}+\boldsymbol{k_{2}}=-n\boldsymbol{k}/(1+K^{2})$. Therefore
the phase difference is exactly $\Delta\phi=2\pi n$, resulting in a fully
constructive interference between emissions from different electron
microbunches and hence an $N_{\mathrm{e}}^{2}$ enhancement of the collective
emission rate. In a practical FEL system, the situation deviates from the
ideal one, and only adjacent microbunches within the coherence length
contribute to the constructive interference, the number of which is estimated
to be about $N_{\mathrm{c}}=22$ for the parameters shown in Table 1.

Meanwhile, if the directions of the emitted photon pair are away from the
$x^{3}$ axis, the phase difference will deviate from $2\pi n$ with a
nontrivial angular dependence

$\begin{split}\Delta\phi(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})&=2\pi\Big{[}(1+K^{2})\\\
&\phantom{=}\times\frac{1+\mathrm{cos}(Z_{2}^{\prime})-(\mathrm{cos}(Z_{2}^{\prime})-\mathrm{cos}(Z_{1}^{\prime}))k_{1}^{0}/k^{0}}{2+K^{2}(1+\mathrm{cos}(Z_{2}^{\prime}))}\Big{]}\end{split}$
(5)

(see details in Section IV of SI), where $Z_{1}^{\prime},Z_{2}^{\prime}$ are
the zenith angles of the two emitted photons relative to the $x^{3}$ axis in
the electron frame. The enhancement parameter on the scattering amplitude by
electron microbunches within one coherent length is given by
$H(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})=\sum_{j=1}^{N_{\mathrm{c}}}N_{j}e^{i[j\Delta\phi(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})]}$,
where $N_{j}$ is the number of electrons within the ${j}^{\mathrm{th}}$
electron microbunch. The collective emission rate of such coherent
microbunches is the product of the emission rate from a single electron and
the enhancement parameter $|H(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})|^{2}$.
For the entire electron pulse, the coherent enhancement factor in Eq. 4 is
given by

$\displaystyle
F_{\mathrm{MB}}=|H(Z_{1}^{\prime},Z_{2}^{\prime},k_{1}^{0})|^{2}N_{\mathrm{i}}\;,$
(6)

where $N_{\mathrm{i}}\simeq N_{\mathrm{e}}/(\sum_{j=1}^{N_{\mathrm{c}}}N_{j})$
represents the approximate number of coherent sections.

In addition, using the Lorentz transformation
$\gamma\mathrm{tan}(Z_{l})=\mathrm{sin}(Z_{l}^{\prime})/(\beta+\mathrm{cos}(Z_{l}^{\prime}))$,
we can determine the zenith angles $Z_{1}$ and $Z_{2}$ in the lab frame for
observation of the entangled photon pair emission. As the density matrix of
the emitted photon pairs from microbunched electrons has an identical form as
that from the single electron emission, the entanglement degree is retained
under microbunching enhancement.

## Numerical results

Figure 2: Numerical simulation of entangled X-ray photon pair emission with
microbunch enhancement. a, Illustration of the angles in b and c. $Z_{1,2}$
refer to the zenith angle between the three-momentum
$\boldsymbol{k}_{1},\boldsymbol{k}_{2}$ of the emitted photons and the $x^{3}$
axis. $A_{2}-A_{1}$ is the difference between the azimuth angles of the
emitted photons at the transverse far-field plane. Photons are treated as
being emitted from the point $S$ at the center of undulator for far-field
measurement. b, Normalized double-angular distribution of the differential
scattering rate $d\dot{W}/dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}$ in the lab
frame. c, Double-angular distribution of the concurrence $\mathcal{C}$. In
both b and c, the polar angles refers to $A_{2}\in[0,2\pi]$, and the radii
refers to $\gamma\mathrm{tan}(Z_{2})\in[0,1]$ for photon 2. The red points
represent the directions of photon 1. Here we choose
$\omega_{1}\simeq\omega_{\mathrm{fd}}/3$.

We firstly calculate the total scattering rate $\dot{W}_{\mathrm{fi}}$. We set
practical ranges of photon energies and angular directions for both photon
detectors, and focus on the case of $n=1$. The relevant angles are illustrated
in Fig. 2a. For the helical undulator, the dependence on the azimuthal angles
$A_{1},A_{2}$ of the emitted photon pair is simplified to $A_{2}-A_{1}$. Thus,
$d\dot{W}/dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}$ can be displayed in a
three-dimensional space of the zenith angles $Z_{1},Z_{2}$ and $A_{2}-A_{1}$.
Choosing $\omega_{1}\simeq\omega_{\mathrm{fd}}/3$ and $\omega_{2}\simeq
2\omega_{\mathrm{fd}}/3$ in the lab frame, the double angular distribution of
the differential scattering rate
$d\dot{W}/dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}$ is shown in Fig. 2b.

To demonstrate the quantum entanglement, we show the concurrence $\mathcal{C}$
of the emitted photon pairs in Fig. 2c. Especially, $\mathcal{C}$ can reach
its upper limit of $1$ for certain angles, meaning that fully-entangled two-
photon states may be generated. Concurrence close to $1$ can be approached at
multiple angles. However, for most regions in the double-angular parameter
space, the concurrence is relatively low. Thus, without selection or
purification, the overall entanglement of the emitted photon pairs would
vanish.

Figure 3: Density matrices of the emitted photon pairs. a, Geometry of the
inserted circular aperture, where $Z$ describes the angular position of the
circular hole relative to the equivalent source point $\mathrm{S}$. b,c,d,
Real parts of the density matrices of the resulting entangled states for the
aperture shown in a with $\gamma\mathrm{tan}(Z)=0.51,0.60,0.68$, which leads
to $\mathcal{C}=0.90,0.99,0.89$, respectively. The size of the circular
aperture is chosen to be $\sim 1/300$ of the far-field speckle. The imaginary
parts of the density matrices are zero due to the specific choice of the
helicity basis and the identical emitting directions of the entangled photons
(see Section VIII of SI for details).

To obtain higher concurrence, a circular aperture can be inserted into the
system at far field as shown in Fig. 3a, and only photon pairs that are
emitted in the direction $Z_{2}\simeq Z_{1}\simeq Z$ can pass through the
aperture. In light of the numerical calculation for the angular distribution
of concurrence shown in Fig. 2b, the circular aperture can be chosen to
maximize the concurrence of the photon pairs that pass it. As demonstrated in
Fig. 3b-d, concurrence close to $1$ can be achieved at multiple angles such as
$\gamma\mathrm{tan}(Z)=0.51,0.60,0.68$ with a sufficiently small aperture,
whose size is chosen to be $\sim 1/300$ of the far-field speckle.

We can now estimate the emission rate of the entangled photon pairs in certain
energy ranges that pass through the aperture in the lab frame. The energy
ranges of the emitted photons are chosen to be
$[\omega_{\mathrm{fd}}/3-10~{}\mathrm{eV},~{}\omega_{\mathrm{fd}}/3]$ and
$[2\omega_{\mathrm{fd}}/3-10~{}\mathrm{eV},~{}2\omega_{\mathrm{fd}}/3]$. The
angular position of the circular hole is chosen to be
$\gamma\mathrm{tan}(Z)=0.60$.

For a single electron, the emission rate of the entangled photon pairs is
$R_{\mathrm{SE}}=6\times 10^{-21}$. For an electron pulse in a synchrotron
system, with the same amount of electrons $N_{\mathrm{e}}=1.12\times 10^{9}$
as in Table 1, the total emission rate is
$R_{\mathrm{Sync}}=N_{\mathrm{e}}R_{\mathrm{SE}}=7\times 10^{-12}$. For a
microbunched electron pulse in the FEL system where the number of electrons in
each microbunch is chosen to be a typical value of $10^{6}$, with Eq. 4, we
obtain a rate of $R_{\mathrm{FEL}}=1.5\times 10^{-4}$. The FEL system can thus
provide an enhancement of order of $10^{7}$ for the entangled photon pair
emission rate compared to the synchrotron system.

## Conclusion

In conclusion, we have theoretically demonstrated the entangled photon pair
generation with FEL. In the electron frame, the double photon emission can be
understood as the scattering between the electrons and the quasi-EM wave of
the undulator field, which can be enhanced by the electron microbunching in
FEL. We demonstrate that the FEL can generate highly entangled photon pairs
with enhanced emission rate, which is impossible from a classical perspective.
The FEL can thus be utilized as an entangled photon pair source for X-ray
quantum optics, especially for applications in the soft X-ray regime where it
is challenging to use nonlinear effects like spontaneous parametric down
conversion due to the strong absorption in nonlinear crystals.

## Methods

We use the Floquet expansion of Volkov state in the the calculation of double
emission scattering matrix $S_{\mathrm{fi}}$. As the dependence of the Volkov
wave function on the phase $\phi$ (thus four-position $x$) shows a periodicity
of $2\pi$, a Fourier expansion can be performed such that

$\Psi^{\mathrm{V}}_{p}(x)=\sum_{\tilde{n}=-\infty}^{+\infty}{F_{\tilde{n}}e^{-i(q\cdot
x+\tilde{n}k\cdot x)}}u(p)\;.$ (7)

Thus, the contributions to $S_{\mathrm{fi}}$ can be divided according to the
number of Floquet photons absorbed, $n$. E.g., for the Feynman diagram shown
in Fig. 1c,

$\begin{split}S_{\mathrm{fi}}&=\int d^{4}x_{1}\int
d^{4}x_{2}\bar{\Psi}^{\mathrm{V}}_{\mathrm{f}}(x_{2})(-ie\not{A}_{2}(x_{2}))\\\
&\phantom{=}\times(iG_{\mathrm{V}}(x_{2},x_{1}))(-ie\not{A}_{1}(x_{1}))\Psi^{\mathrm{V}}_{\mathrm{i}}(x_{1})\\\
&\phantom{=}+\int d^{4}x_{1}\int
d^{4}x_{2}\bar{\Psi}^{\mathrm{V}}_{\mathrm{f}}(x_{2})(-ie\not{A}_{1}(x_{2}))\\\
&\phantom{=}\times(iG_{\mathrm{V}}(x_{2},x_{1}))(-ie\not{A}_{2}(x_{1}))\Psi^{\mathrm{V}}_{\mathrm{i}}(x_{1})\\\
&=-i\frac{e^{2}}{2V}\sqrt{\frac{1}{k_{2}^{0}k_{1}^{0}}}\sum_{n=-\infty}^{+\infty}{\bigg{\\{}\sum_{n_{1}=-\infty}^{+\infty}{M_{\mathrm{S}}(n,n_{1})\bigg{|}_{\mathrm{MC}(n_{1})}}}\\\
&\phantom{=}\times(2\pi)^{4}\delta^{(4)}(q_{\mathrm{f}}+k_{1}+k_{2}-q_{\mathrm{i}}-nk)\bigg{\\}}\;,\end{split}$
(8)

where the helicity-relevant part
$M_{\mathrm{S}}(n,n_{1})|_{\mathrm{MC}(n_{1})}$ is explained in Section II of
SI, which only survives when momentum conservation is fulfilled at each vertex
of the Feynman diagrams.

The differential scattering rate is then calculated by considering the total
scattering probability $\dot{W}$ of double emission that satisfies

$\displaystyle\dot{W}=$
$\displaystyle\int{V\frac{d^{3}\boldsymbol{q}_{\mathrm{f}}}{(2\pi)^{3}}V\frac{d^{3}\boldsymbol{k}_{1}}{(2\pi)^{3}}V\frac{d^{3}\boldsymbol{k}_{2}}{(2\pi)^{3}}\frac{|S_{\mathrm{fi}}|^{2}}{VT}V}$
(9) $\displaystyle=$
$\displaystyle\int\frac{4V^{4}}{(2\pi)^{9}}d^{3}\boldsymbol{q}_{\mathrm{f}}\int{d\Omega_{k_{1}}}\int{d\Omega_{k_{2}}}\int_{0}^{+\infty}{dk_{1}^{0}}\int_{0}^{+\infty}{dk_{2}^{0}}\,(k_{1}^{0})^{2}(k_{2}^{0})^{2}\frac{|S_{\mathrm{fi}}|^{2}}{VT}\;,$

where $VT$ is the space-time volume of the interaction zone. Combining Eq. 8
and Eq. 9, and completing the integral over variables except for $k_{1}^{0}$,
$\Omega_{k_{1}}$ and $\Omega_{k_{2}}$, we obtain the Eq. 2 that

$\begin{split}&\frac{d\dot{W}}{dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}}=\frac{e^{4}}{(2\pi)^{5}}\frac{m^{2}}{2q_{\mathrm{i}}^{0}}\\\
&\phantom{=}\times\sum_{n=1}^{+\infty}{\frac{k_{1}^{0}(k_{2}^{0})^{2}\Theta(q_{\mathrm{i}}^{0}+nk^{0}-k_{1}^{0}-k_{2}^{0})M}{2|(q_{\mathrm{i}}+nk-
k_{1})\cdot k_{2}|}}\bigg{|}_{\mathrm{DE}(n)}\;,\end{split}$

where $|_{\mathrm{DE}(n)}$ means the constraint of $(q_{\mathrm{i}}+nk-
k_{1})\cdot k_{2}=q_{\mathrm{i}}\cdot nk-q_{\mathrm{i}}\cdot k_{1}-nk\cdot
k_{1}$, and the helicity-relevant part $M$ is explained in Section II of SI.

In order to explicitly show the dependence of the entangled emission
probability on the key FEL parameters, we can reduce the formula of the
differential entangled double emission probability for a single electron in
the lab frame based on Eq. 2. Focusing on $n=1$, the entangled emission
probability $U_{\mathrm{pair}}$ is given by

$U_{\mathrm{pair}}\simeq
F_{\mathrm{L}}F_{\mathrm{S}}\frac{d\dot{W}}{dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}}\bigg{|}_{\mathrm{EF},k_{1},k_{2}}\;,$
(10)

where
$\frac{d\dot{W}}{dk_{1}^{0}d\Omega_{k_{1}}d\Omega_{k_{2}}}\bigg{|}_{\mathrm{EF},k_{1},k_{2}}$
is the value of scattering probability in Eq. 2 with the two photon four-
momenta in EF chosen as $k_{1},k_{2}$, and $F_{\mathrm{L}}$, $F_{\mathrm{S}}$
are the prefactors related to the Lorentz transformation and solid angle,
respectively.

$U_{\mathrm{pair}}$ can be expressed with key FEL parameters $\gamma$,
$\lambda_{u}$ and $B_{0}$, which represent the Lorentz factor of electrons,
undulator period length and magnetic field amplitude of the undulator, in the
form of

$\displaystyle U_{\mathrm{pair}}$
$\displaystyle\propto\frac{Q(B_{0}\lambda_{u})}{\lambda_{u}^{2}(1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2})^{2}(1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2}/2)[1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2}(1+\mathrm{cos}Z_{2}^{\prime})/2]}\;,$
(11)

where

$\displaystyle Q(B_{0}\lambda_{u})$
$\displaystyle=\sum_{t_{1},t_{2},t_{3},t_{4}=0}^{2}P_{t_{1},t_{2},t_{3},t_{4}}(B_{0}\lambda_{u})J_{t_{1}}(\mathrm{sin}Z_{1}^{\prime}\frac{(c_{\mathrm{K}}B_{0}\lambda_{u})}{1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2}})$
$\displaystyle\times
J_{t_{2}}(\mathrm{sin}Z_{1}^{\prime}\frac{(c_{\mathrm{K}}B_{0}\lambda_{u})}{1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2}})J_{t_{3}}(\mathrm{sin}Z_{2}^{\prime}\frac{(c_{\mathrm{K}}B_{0}\lambda_{u})}{1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2}})$
$\displaystyle\times
J_{t_{4}}(\mathrm{sin}Z_{2}^{\prime}\frac{(c_{\mathrm{K}}B_{0}\lambda_{u})}{1+(c_{\mathrm{K}}B_{0}\lambda_{u})^{2}})\;,$

and $c_{\mathrm{K}}=e/2\pi m$. $P_{t_{1},t_{2},t_{3},t_{4}}(B_{0}\lambda_{u})$
are polynomials of $B_{0}\lambda_{u}$ with the highest order between 4 and 8,
and $J_{t}(x)$ refers to the Bessel function of the first kind with order $t$.
The detailed derivation of Eq. 11 is presented in Section V of SI.

At the vicinity of a specific choice of $k_{1},k_{2}$, we can obtain
approximated power function for $Q(B_{0}\lambda_{u})$. For example, consider
$k_{1},k_{2}$ that satisfy $\omega_{1}=\omega_{\mathrm{fd}}/3$ and
$\gamma\mathrm{tan}Z_{1,2}\simeq 0.51,0.60,0.68$ in the lab frame, with the
corresponding $\omega_{2}\simeq 2\omega_{\mathrm{fd}}/3$ given by the
conservation of four-quasi-momentum. These three cases correspond to the three
circumstances shown in Fig. 3 of the main text. we obtain by fitting that

$\displaystyle Q(B_{0}\lambda_{u})\propto$
$\displaystyle(B_{0}\lambda_{u})^{1.90}\quad\mathrm{for~{}~{}\gamma\mathrm{tan}Z_{1,2}\simeq
0.51}\;,$ (12) $\displaystyle Q(B_{0}\lambda_{u})\propto$
$\displaystyle(B_{0}\lambda_{u})^{2.00}\quad\mathrm{for~{}~{}\gamma\mathrm{tan}Z_{1,2}\simeq
0.60}\;,$ (13) $\displaystyle Q(B_{0}\lambda_{u})\propto$
$\displaystyle(B_{0}\lambda_{u})^{2.10}\quad\mathrm{for~{}~{}\gamma\mathrm{tan}Z_{1,2}\simeq
0.68}\;.$ (14)

The details are shown in Section V of SI.

The helical undulator simplifies the dependence of the double-emission process
on the emitting azimuth angles $A_{1},A_{2}$ to $A_{2}-A_{1}$ and reduces the
parameter space, making it more convenient to optimize the detection scheme
for the entangled photon pairs in experiment. A detailed calculation for the
case of linearly polarized undulator is given in Section VI in SI.

## Data availability

The data that support the findings of this study are available from the
corresponding authors upon reasonable request.

## Code availability

The code for the simulation is available from the corresponding authors upon
reasonable request.

## Acknowledgement

The authors are thankful to Zhirong Huang and Zhaoheng Guo for helpful
discussions. This work was supported by the National Natural Science
Foundation of China (No. 12174009, 11974031, 12174011).

## References

  *   * (1) Einstein, A., Podolsky, B., Rosen, N.: Can quantum-mechanical description of physical reality be considered complete? Phys. Rev. 47(10), 777 (1935)
  * (2) Bell, J.S.: On the problem of hidden variables in quantum mechanics. Rev. Mod. Phys. 38(3), 447 (1966)
  * (3) Shwartz, S., Coffee, R.N., Feldkamp, J.M., Feng, Y., Hastings, J.B., Yin, G.Y., Harris, S.E.: X-ray parametric down-conversion in the Langevin regime. Phys. Rev. Lett. 109(1), 013602 (2012)
  * (4) Emma, P., Akre, R., Arthur, J., Bionta, R., Bostedt, C., Bozek, J., Brachmann, A., Bucksbaum, P., Coffee, R., Decker, F.J., et al.: First lasing and operation of an ångstrom-wavelength free-electron laser. Nature Photon. 4(9), 641–647 (2010)
  * (5) Ishikawa, T., Aoyagi, H., Asaka, T., Asano, Y., Azumi, N., Bizen, T., Ego, H., Fukami, K., Fukui, T., Furukawa, Y., et al.: A compact X-ray free-electron laser emitting in the sub-ångström region. Nature Photon. 6(8), 540–544 (2012)
  * (6) Allaria, E., Castronovo, D., Cinquegrana, P., Craievich, P., Dal Forno, M., Danailov, M.B., D’Auria, G., Demidovich, A., De Ninno, G., Di Mitri, S., et al.: Two-stage seeded soft-X-ray free-electron laser. Nature Photon. 7(11), 913–918 (2013)
  * (7) Wang, W., Feng, K., Ke, L., Yu, C., Xu, Y., Qi, R., Chen, Y., Qin, Z., Zhang, Z., Fang, M., et al.: Free-electron lasing at 27 nanometres based on a laser wakefield accelerator. Nature 595(7868), 516–520 (2021)
  * (8) Pompili, R., Alesini, D., Anania, M.P., Arjmand, S., Behtouei, M., Bellaveglia, M., Biagioni, A., Buonomo, B., Cardelli, F., Carpanese, M., et al.: Free-electron lasing with compact beam-driven plasma wakefield accelerator. Nature 605(7911), 659–662 (2022)
  * (9) Duris, J., Li, S., Driver, T., Champenois, E.G., MacArthur, J.P., Lutman, A.A., Zhang, Z., Rosenberger, P., Aldrich, J.W., Coffee, R., et al.: Tunable isolated attosecond X-ray pulses with gigawatt peak power from a free-electron laser. Nature Photon. 14(1), 30–36 (2020)
  * (10) Duris, J.P., MacArthur, J.P., Glownia, J.M., Li, S., Vetter, S., Miahnahri, A., Coffee, R., Hering, P., Fry, A., Welch, M.E., et al.: Controllable x-ray pulse trains from enhanced self-amplified spontaneous emission. Phys. Rev. Lett. 126(10), 104802 (2021)
  * (11) Seibert, M.M., Ekeberg, T., Maia, F.R.N.C., Svenda, M., Andreasson, J., Jönsson, O., Odić, D., Iwan, B., Rocker, A., Westphal, D., et al.: Single mimivirus particles intercepted and imaged with an X-ray laser. Nature 470(7332), 78–81 (2011)
  * (12) Loh, N.D., Hampton, C.Y., Martin, A.V., Starodub, D., Sierra, R.G., Barty, A., Aquila, A., Schulz, J., Lomb, L., Steinbrener, J., et al.: Fractal morphology, imaging and mass spectrometry of single aerosol particles in flight. Nature 486(7404), 513–517 (2012)
  * (13) Gomez, L.F., Ferguson, K.R., Cryan, J.P., Bacellar, C., Tanyag, R.M.P., Jones, C., Schorb, S., Anielski, D., Belkacem, A., Bernando, C., et al.: Shapes and vorticities of superfluid helium nanodroplets. Science 345(6199), 906–909 (2014)
  * (14) Kern, J., Alonso-Mori, R., Tran, R., Hattne, J., Gildea, R.J., Echols, N., Glöckner, C., Hellmich, J., Laksmono, H., Sierra, R.G., et al.: Simultaneous femtosecond X-ray spectroscopy and diffraction of photosystem II at room temperature. Science 340(6131), 491–495 (2013)
  * (15) Erk, B., Boll, R., Trippel, S., Anielski, D., Foucar, L., Rudek, B., Epp, S.W., Coffee, R., Carron, S., Schorb, S., et al.: Imaging charge transfer in iodomethane upon x-ray photoabsorption. Science 345(6194), 288–291 (2014)
  * (16) Rudek, B., Son, S., Foucar, L., Epp, S.W., Erk, B., Hartmann, R., Adolph, M., Andritschke, R., Aquila, A., Berrah, N., et al.: Ultra-efficient ionization of heavy atoms by intense X-ray free-electron laser pulses. Nature Photon. 6(12), 858–865 (2012)
  * (17) Vinko, S.M., Ciricosta, O., Cho, B.I., Engelhorn, K., Chung, H.-K., Brown, C.R.D., Burian, T., Chalupskỳ, J., Falcone, R.W., Graves, C., et al.: Creation and diagnosis of a solid-density plasma with an X-ray free-electron laser. Nature 482(7383), 59–62 (2012)
  * (18) Fletcher, L.B., Lee, H.J., Döppner, T., Galtier, E., Nagler, B., Heimann, P., Fortmann, C., LePape, S., Ma, T., Millot, M., et al.: Ultrabright X-ray laser scattering for dynamic warm dense matter physics. Nature Photon. 9(4), 274–279 (2015)
  * (19) Ciricosta, O., Vinko, S.M., Barbrel, B., Rackstraw, D.S., Preston, T.R., Burian, T., Chalupskỳ, J., Cho, B.I., Chung, H.-K., Dakovski, G.L., et al.: Measurements of continuum lowering in solid-density plasmas created from elements and compounds. Nat. Commun. 7(1), 1–6 (2016)
  * (20) Pellegrini, C., Marinelli, A., Reiche, S.: The physics of x-ray free-electron lasers. Rev. Mod. Phys. 88(1), 015006 (2016)
  * (21) Lutman, A.A., MacArthur, J.P., Ilchen, M., Lindahl, A.O., Buck, J., Coffee, R.N., Dakovski, G.L., Dammann, L., Ding, Y., Dürr, H.A., et al.: Polarization control in an X-ray free-electron laser. Nature Photon. 10(7), 468–472 (2016)
  * (22) Madey, J.M.J.: Stimulated emission of bremsstrahlung in a periodic magnetic field. J. Appl. Phys. 42(5), 1906–1913 (1971)
  * (23) Temnykh, A.B.: Delta undulator for Cornell energy recovery linac. Phys. Rev. ST Accel. Beams 11(12), 120702 (2008)
  * (24) Huang, Z., Kim, K.: Review of x-ray free-electron laser theory. Phys. Rev. ST Accel. Beams 10(3), 034801 (2007)
  * (25) Ritus, V.I.: Quantum effects of the interaction of elementary particles with an intense electromagnetic field. J. Sov. Laser Res. 6, 497–617 (1985)
  * (26) Lötstedt, E., Jentschura, U.D.: Nonperturbative treatment of double Compton backscattering in intense laser fields. Phys. Rev. Lett. 103(11), 110404 (2009)
  * (27) Lötstedt, E., Jentschura, U.D.: Correlated two-photon emission by transitions of Dirac-Volkov states in intense laser fields: QED predictions. Phys. Rev. A 80(5), 053419 (2009)
  * (28) Wootters, W.K.: Entanglement of formation of an arbitrary state of two qubits. Phys. Rev. Lett. 80(10), 2245 (1998)