# On the observability of field-assisted birefringent Delbrück scattering

N. Ahmadiniaz  Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400,
01328 Dresden, Germany  M. Bussmann  Helmholtz-Zentrum Dresden-Rossendorf,
Bautzner Landstraße 400, 01328 Dresden, Germany  Center for Advanced Systems
Understanding (CASUS), Görlitz, Germany  T. E. Cowan  Helmholtz-Zentrum
Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany  Institut
für Kern-und Teilchenphysik, Technische Universität Dresden, 01062 Dresden,
Germany  A. Debus  Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße
400, 01328 Dresden, Germany  T. Kluge  Helmholtz-Zentrum Dresden-Rossendorf,
Bautzner Landstraße 400, 01328 Dresden, Germany  R. Schützhold  Helmholtz-
Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany
Institut für Theoretische Physik, Technische Universität Dresden, 01062
Dresden, Germany

###### Abstract

We consider the scattering of an x-ray free-electron laser (XFEL) beam on the
superposition of a strong magnetic field $\bf{B}_{\rm ext}$ with the Coulomb
field $\bf{E}_{\rm ext}$ of a nucleus with charge number $Z$. In contrast to
pure Delbrück scattering (Coulomb field only), the magnetic field $\bf{B}_{\rm
ext}$ introduces an asymmetry (i.e., polarization dependence) and renders the
effective interaction volume quite large, while the nuclear Coulomb field
facilitates a significant momentum transfer $\Delta\bf k$. For a field
strength of $B_{\rm ext}=10^{6}~{}\rm T$ (corresponding to an intensity of
order $10^{22}~{}\rm W/cm^{2}$) and an XFEL frequency of 24 keV, we find a
differential cross section $d\sigma/d\Omega\sim 10^{-25}~{}Z^{2}/(\Delta{\bf
k})^{2}$ in forward direction for one nucleus. Thus, this effect might be
observable in the near future at facilities such as the Helmholtz
International Beamline for Extreme Fields (HIBEF) at the European XFEL.

#### Introduction

According to classical electrodynamics, electromagnetic waves in vacuum obey
the superposition principle and thus do not influence each other. Quantum
electrodynamics (QED), on the other hand, predicts that they do interact via
their coupling to the fermionic degrees of freedom euler-35 ; euler+heisenberg
; karplus+neuman . The difficulties of observing this interaction can be
understood by recalling the characteristic scales of QED. First, the electron
mass $m\approx 0.51~{}{\rm MeV}/c^{2}$ sets an energy scale where the
associated length scale $\lambdabar=\hbar/(mc)\approx 386~{}\rm fm$ is the
Compton length. Second, the elementary charge $q$ can be used to construct the
Schwinger critical field sauter-31 ; schwinger-51

$\displaystyle E_{\rm crit}=\frac{m^{2}c^{3}}{\hbar q}\approx 1.3\times
10^{18}\,\frac{\rm V}{\rm m}\,,$ (1)

where the corresponding magnetic field $B_{\rm crit}=E_{\rm crit}/c$ is given
by $B_{\rm crit}\approx 4.4\times 10^{9}~{}\rm T$. It is extremely hard to
reach such field strengths in the laboratory, but even stronger fields exist
in extra-terrestrial environments, cf. neutron-star ; neutron-star-comment ;
mesz ; ruffini .

As an exception, the Coulomb field of a nucleus with charge $Zq$ exceeds the
field strength (1) very close to the nucleus, i.e., on a distance of order
$\,{\cal O}(\sqrt{Z}\lambdabar)$. Such a high field strength helps to observe
the interaction of electromagnetic fields and the scattering of photons at the
nuclear Coulomb field (usually referred to as Delbrück scattering meitner-33 ;
bethe+rohrlich ; costantini-71 ) has been observed in several experiments,
see, e.g., cheng1 ; cheng2 ; milshtein-1 ; milshtein-2 ; schumacher-75 ;
falkenberg-92 ; DESY-exp ; jackson-69 ; akhm-98 ; moreh ; kahane ; rev1 ; rev2
; rev3 ; rev4 ; rull-81 . In order to probe the high field in the close
vicinity of the nucleus, these photons had an energy of the order of the
electron mass or above. As another example, the interaction of the Coulomb
fields of two nuclei almost colliding with each other at ultra-high energies
and the resulting emission of a pair of photons has also been observed atlas .

Figure 1: Exemplary Feynman diagram of the considered process (wavy/straight
lines represent photons/electrons).

In contrast to the cases mentioned above, the interaction of electromagnetic
fields (as predicted by QED) in a regime well below the QED scales $mc^{2}$,
$E_{\rm crit}$ and $B_{\rm crit}$ has not been observed yet. Prominent
proposals for ongoing and planned experiments include the interaction of an
optical (or near-optical) laser with a (quasi) static magnetic field of a few
Tesla, see, e.g., della-16 ; pvlas ; zava-12 ; zava-06 ; zava-08 ; OVAL ; BMV
; hartman-17 ; battesti-18 , the interaction of x-ray free electron laser
(XFEL) beams among each other or with optical lasers, see, e.g., heinzl-06 ;
inada-14 ; yamaji-16 ; Schlenvoigt-2016 ; Karbstein-2016 ; inada-2017 , or the
interaction of several optical lasers, see, e.g., dipiazza-07 ; tommasini-09 ;
tommasini-10 ; koga-17 ; king-12 ; gies1 ; gies2 ; gies3 ; gies-09 ; grote-15
; ahm-20 .

Here, we consider a mixed set-up where an XFEL beam is scattered at the
combination of a strong magnetic field $\bf{B}_{\rm ext}$ superimposed by the
Coulomb field $\bf{E}_{\rm ext}$ of a nucleus (as schematically depicted in
Fig. 1), see also dipiazza-08 . This scenario offers several advantages: The
nuclear Coulomb field facilitates a significant momentum transfer $\Delta\bf
k$ while we still obtain a birefringent signal (where the polarization of the
XFEL photons flips). The crux of our proposal and the prime difference to the
aforementioned scenarios involving nuclear fields studied previously is that
we obtain a very large interaction volume whose length and time scales are set
by the momentum transfer $\Delta\bf k$ and thus well above the Compton length
$\lambdabar$. Hence, all involved field strengths are sub-critical, i.e., well
below (1).

#### Euler-Heisenberg Lagrangian

Since we are considering slowly varying electric $\mathfrak{E}$ and magnetic
$\mathfrak{B}$ fields well below $E_{\rm crit}$ and $B_{\rm crit}$, we may
start with the lowest-order Euler-Heisenberg Lagrangian
($\hbar=c=\epsilon_{0}=\mu_{0}=1$)

$\displaystyle\mathfrak{L}=\frac{1}{2}\left(\mathfrak{E}^{2}-\mathfrak{B}^{2}\right)+\xi\left[\left(\mathfrak{E}^{2}-\mathfrak{B}^{2}\right)^{2}+7\left(\mathfrak{E}\cdot\mathfrak{B}\right)^{2}\right],$
(2)

with a non-linearity set by the parameter

$\displaystyle\xi=\frac{q^{4}}{360\pi^{2}m^{4}}=\frac{2\alpha_{\rm
QED}^{2}}{45m^{4}}=\frac{\alpha_{\rm QED}}{90\pi E_{\rm crit}^{2}}\,,$ (3)

where $\alpha_{\rm QED}\approx 1/137$ is the fine-structure constant
euler+heisenberg , see also review-karbstein-2020 ; review-king+heinzl-2016 ;
review-karbstein-2016 ; review-dipiazza-12 ; review-dunne-12 ; akhm-02 ;
Bialynicki-Birula-1970 ; Toll-PhD ; adler70 ; adler71 ; adler96 . This is a
great advantage because we do not have to construct the electron propagator
whose explicit form is known in special cases only foot1 , e.g., in an
electromagnetic plane-wave background which facilitates Volkov solutions
volkov ; red-65 , see also milstein-05 ; diphatkei-08 ; dipiazza-12 ;
dipiazza-14 .

In the weak-field limit considered here, we neglect quadratic terms $\,{\cal
O}(\xi^{2})$. Furthermore, we assume that the magnetic field $\bf{B}_{\rm
ext}$ is (approximately) constant. In addition to this magnetic field
$\bf{B}_{\rm ext}$ and the static Coulomb field $\bf{E}_{\rm ext}$ of the
nucleus, we have the space-time dependent XFEL fields ${\bf E}$ and ${\bf B}$.
Inserting this split $\mathfrak{E}={\bf E}_{\rm ext}+{\bf E}$ and
$\mathfrak{B}={\bf B}_{\rm ext}+{\bf B}$ into (2), we obtain the effective
Lagrangian for the XFEL fields (cf. ahm-20 )

$\displaystyle\mathfrak{L}_{\rm XFEL}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[{\bf
E}\cdot(\mathbb{1}+\delta\epsilon)\cdot{\bf E}-{\bf
B}\cdot(\mathbb{1}-\delta\mu)\cdot{\bf B}\right]$ (4) $\displaystyle+{\bf
E}\cdot\delta\Psi\cdot{\bf B}\,,$

with the symmetric permittivity/permeability tensors $\delta\epsilon^{ij}=8\xi
E_{\rm ext}^{i}E_{\rm ext}^{j}+14\xi B_{\rm ext}^{i}B_{\rm
ext}^{j}+4\xi\delta^{ij}({\bf E}_{\rm ext}^{2}-{\bf B}_{\rm ext}^{2})$ and
$\delta\mu^{ij}=8\xi B_{\rm ext}^{i}B_{\rm ext}^{j}+14\xi E_{\rm
ext}^{i}E_{\rm ext}^{j}-4\xi\delta^{ij}({\bf E}_{\rm ext}^{2}-{\bf B}_{\rm
ext}^{2})$ plus the symmetry-breaking contribution

$\displaystyle\delta\Psi^{ij}$ $\displaystyle=$ $\displaystyle-8\xi E_{\rm
ext}^{i}B_{\rm ext}^{j}+14\xi B_{\rm ext}^{i}E_{\rm ext}^{j}$ (5)
$\displaystyle+14\xi\delta^{ij}({\bf E}_{\rm ext}\cdot{\bf B}_{\rm ext})\,,$

which describe the polarizability of the QED vacuum. Note that the latter
tensor is not symmetric $\delta\Psi^{ij}\neq\delta\Psi^{ji}$.

The equations of motion stemming from (4) can be cast into the same form as
the macroscopic Maxwell equations in a medium $\nabla\cdot{\bf D}=0$,
$\nabla\cdot{\bf B}=0$, $\nabla\times{\bf E}=-\partial_{t}{\bf B}$, and
$\nabla\times{\bf H}=\partial_{t}{\bf D}$, provided that we introduce the
electric ${\bf D}=(\mathbb{1}+\delta\epsilon)\cdot{\bf E}+\delta\Psi\cdot{\bf
B}$ and magnetic displacement fields ${\bf H}=(\mathbb{1}-\delta\mu)\cdot{\bf
B}-\delta\Psi^{\rm T}\cdot{\bf E}$.

#### Scattering Theory

Now we may calculate the scattering of the XFEL beam with standard approaches,
see, e.g., jackson . Combining the above Maxwell equations to

$\displaystyle\Box{\bf D}=\nabla\times\left[\nabla\times({\bf D}-{\bf
E})\right]+\partial_{t}\left[\nabla\times({\bf H}-{\bf B})\right]={\bf J}^{\rm
eff}\,,$ (6)

where the effective source term ${\bf J}^{\rm eff}$ on the right-hand side is
small, allows us to employ the Born approximation. To this end, we split the
XFEL field ${\bf D}$ into an ingoing plane wave ${\bf D}^{\rm in}$ plus a
small scattering contribution ${\bf D}^{\rm out}$ induced by vacuum
polarizability $\delta\epsilon$, $\delta\mu$, and $\delta\Psi$. Assuming a
stationary time-dependence $e^{-i\omega t}$ for the XFEL field
($\delta\epsilon$, $\delta\mu$, and $\delta\Psi$ are static), we find

$\displaystyle\Box{\bf D}_{\omega}^{\rm out}$ $\displaystyle=$
$\displaystyle-\left(\nabla^{2}+\omega^{2}\right){\bf D}_{\omega}^{\rm
out}={\bf J}_{\omega}^{\rm eff}$ (7) $\displaystyle=$
$\displaystyle\nabla\times\left[\nabla\times(\delta\epsilon\cdot{\bf
E}_{\omega}^{\rm in}+\delta\Psi\cdot{\bf B}_{\omega}^{\rm in})\right]$
$\displaystyle+i\omega\nabla\times(\delta\mu\cdot{\bf B}_{\omega}^{\rm
in}+\delta\Psi^{\rm T}\cdot{\bf E}_{\omega}^{\rm in})\,,$

where we may approximate ${\bf E}_{\omega}^{\rm in}\approx{\bf
D}_{\omega}^{\rm in}$ and ${\bf B}_{\omega}^{\rm in}\approx{\bf
H}_{\omega}^{\rm in}$ on the right-hand side within the Born approximation. As
usual, solving this Helmholtz differential equation for ${\bf D}_{\omega}^{\rm
out}$ with the standard Greens function and considering the behavior at large
spatial distances, we find the differential cross section
$d\sigma/d\Omega=|{\mathfrak{A}}|^{2}$ with the scattering amplitude jackson

$\displaystyle{\mathfrak{A}}=\frac{1}{4\pi|{\bf D}_{\omega}^{\rm in}|}\,{\bf
e}_{\rm out}\cdot\int d^{3}r\,\exp\\{-i{\bf k}_{\rm out}\cdot{\bf r}\\}\,{\bf
J}_{\omega}^{\rm eff}\,,$ (8)

where ${\bf k}_{\rm out}$ is the wave-number and ${\bf e}_{\rm out}$ the
polarization unit vector of the scattered (outgoing) XFEL radiation.

#### Field-Assisted Scattering

The effective source term ${\bf J}_{\omega}^{\rm eff}$ contains contributions
from $\delta\epsilon$, $\delta\mu$, and $\delta\Psi$ which add up to give the
full amplitude

$\displaystyle{\mathfrak{A}}$ $\displaystyle=$
$\displaystyle\frac{\omega^{2}}{4\pi}\int d^{3}r\,\exp\\{i\Delta{\bf
k}\cdot{\bf r}\\}\big{[}{\bf e}_{\rm out}\cdot\delta\epsilon\cdot{\bf e}_{\rm
in}$ (9) $\displaystyle+{\bf e}_{\rm out}\cdot\delta\Psi\cdot({\bf n}_{\rm
in}\times{\bf e}_{\rm in})+{\bf e}_{\rm in}\cdot\delta\Psi\cdot({\bf n}_{\rm
out}\times{\bf e}_{\rm out})$ $\displaystyle+({\bf n}_{\rm out}\times{\bf
e}_{\rm out})\cdot\delta\mu\cdot({\bf n}_{\rm in}\times{\bf e}_{\rm
in})\big{]}\,.$

Here ${\bf n}_{\rm in}$ and ${\bf n}_{\rm out}$ are the initial and final
propagation directions (${\bf k}_{\rm in}=\omega{\bf n}_{\rm in}$ and ${\bf
k}_{\rm out}=\omega{\bf n}_{\rm out}$), respectively, while ${\bf e}_{\rm in}$
and ${\bf e}_{\rm out}$ denote their polarizations. As usual in scattering
theory, the oscillating phase is governed by the momentum transfer $\Delta{\bf
k}={\bf k}_{\rm in}-{\bf k}_{\rm out}$.

In the following, we only consider the terms stemming from the combined impact
of the magnetic field ${\bf B}_{\rm ext}$ and the nuclear Coulomb field ${\bf
E}_{\rm ext}$. Keeping only those cross terms, just the symmetry-breaking
contribution $\delta\Psi$ in (5) survives, i.e., we focus on the second line
of Eq. (9).

Since ${\bf B}_{\rm ext}$ is (nearly) constant, the spatial integral yields
the Fourier transform of the nuclear Coulomb field ${\bf E}_{\rm ext}({\bf
r})={\bf e}_{r}Q/(4\pi{\bf r}^{2})$ with the charge $Q=Zq$

$\displaystyle\int d^{3}r\,e^{i\Delta{\bf k}\cdot{\bf r}}\,\frac{Q}{4\pi{\bf
r}^{2}}\,{\bf e}_{r}=iQ\,\frac{\Delta{\bf k}}{(\Delta{\bf k})^{2}}\,.$ (10)

An important point here is that the $1/r^{2}$ scaling from the Coulomb field
cancels the $r^{2}$ volume factor in the $d^{3}r$ integration. As a
consequence, the spatial integration is cut off by the momentum transfer
$\Delta{\bf k}$ resulting in a very large interaction volume – which may even
span many XFEL wavelengths for small scattering angles $|{\bf n}_{\rm in}-{\bf
n}_{\rm out}|\ll 1$, i.e., in forward direction. Of course, at some point the
approximation of a constant ${\bf B}_{\rm ext}$ breaks down.

#### Forward Scattering

The large interaction volume mentioned above goes along with a peak in forward
direction, i.e., for small $\Delta{\bf k}$, where the signal is enhanced as
$\sim 1/|\Delta{\bf k}|$. Since the pre-factor is quite small, let us focus on
this leading-order contribution $\sim 1/|\Delta{\bf k}|$. Thus, we approximate
${\bf n}_{\rm in}\approx{\bf n}_{\rm out}\to{\bf n}$ in order to simplify the
expressions. In this limit, the isotropic contribution $14\xi\delta^{ij}({\bf
E}_{\rm ext}\cdot{\bf B}_{\rm ext})$ in (5) cancels and we are left with the
anisotropic terms. For the birefringent signal, where ${\bf e}_{\rm in}$ and
${\bf e}_{\rm out}$ are (nearly) orthogonal, we may approximate ${\bf e}_{\rm
out}\approx\pm{\bf n}\times{\bf e}_{\rm in}$ and ${\bf e}_{\rm
in}\approx\mp{\bf n}\times{\bf e}_{\rm out}$ which simplifies the integrand in
(9) to $\pm[{\bf e}_{\rm out}\cdot\delta\Psi\cdot{\bf e}_{\rm out}-{\bf
e}_{\rm in}\cdot\delta\Psi\cdot{\bf e}_{\rm in}]$. Altogether, we get the
birefringent amplitude

$\displaystyle{\mathfrak{A}}_{\delta\Psi}^{\perp}$ $\displaystyle=$
$\displaystyle\pm 6i\xi\,\frac{Q}{(\Delta{\bf
k})^{2}}\,\frac{\omega^{2}}{4\pi}\left[({\bf e}_{\rm out}\cdot{\bf B}_{\rm
ext})({\bf e}_{\rm out}\cdot\Delta{\bf k})\right.$ (11)
$\displaystyle\left.-({\bf e}_{\rm in}\cdot{\bf B}_{\rm ext})({\bf e}_{\rm
in}\cdot\Delta{\bf k})\right]\,.$

Thus, one way to obtain a maximum birefringent signal would be to align the
momentum transfer $\Delta{\bf k}$ with ${\bf B}_{\rm ext}$ as well as either
${\bf e}_{\rm in}$ or ${\bf e}_{\rm out}$, for example.

In this case, we find the amplitude (up to a sign)

$\displaystyle{\mathfrak{A}}_{\delta\Psi}^{\perp}=6i\xi B_{\rm
ext}\,\frac{Q}{|\Delta{\bf k}|}\,\frac{\omega^{2}}{4\pi}=i\,\frac{\alpha_{\rm
QED}^{2}}{15\pi}\,\frac{qB_{\rm
ext}}{m^{2}}\,\frac{\omega^{2}}{m^{2}}\,\frac{Z}{|\Delta{\bf k}|}\,.$ (12)

Apart from a numerical pre-factor, this amplitude scales with the ratios of
the optical laser field strength $B_{\rm ext}$ over the critical field
strength $B_{\rm crit}\approx 4.4\times 10^{9}~{}\rm T$ and the square of the
XFEL frequency $\omega$ in comparison to the electron mass $m$.

#### Detectability

Since the amplitude (12) for field-assisted scattering is proportional to
$B_{\rm ext}$ and $\omega^{2}$, it is desirable to have a strong magnetic
field and a large XFEL frequency. Ultra-high field strengths $B_{\rm
ext}=10^{6}~{}\rm T$ can be reached in the focus of an ultra-strong optical or
near optical laser with an intensity of order $10^{22}~{}\rm W/cm^{2}$ in a
colliding-beam set-up. Then, inserting an XFEL frequency of $\omega=24~{}\rm
keV$, which is still within the range of the European XFEL xfel-data ;
Schlenvoigt-2016 ; xfel-URL , we obtain an amplitude of around
${\mathfrak{A}}_{\delta\Psi}^{\perp}\approx 5\times 10^{-13}Z/|\Delta{\bf k}|$
which yields the birefringent differential cross section in forward direction

$\displaystyle\frac{d\sigma_{\delta\Psi}^{\perp}}{d\Omega}=\left|{\mathfrak{A}}_{\delta\Psi}^{\perp}\right|^{2}\sim
10^{-25}\,\frac{Z^{2}}{(\Delta{\bf k})^{2}}\,.$ (13)

For a lower frequency of $\omega=1~{}\rm keV$, the amplitude would be reduced
to ${\mathfrak{A}}_{\delta\Psi}^{\perp}\approx 10^{-15}Z/|\Delta{\bf k}|$
corresponding to a cross section of $10^{-30}Z^{2}/(\Delta{\bf k})^{2}$.

The suppression in (13) by more than twenty orders of magnitude is roughly
comparable to other proposals for vacuum birefringence experiments, see, e.g.,
heinzl-06 ; Schlenvoigt-2016 . Yet, this suppression does not imply that the
effect is beyond reach. In order to demonstrate that, let us discuss two main
enhancement factors. The first enhancement factors is a large number $\,{\cal
O}(10^{11})$ of polarized XFEL photons, analogous to other vacuum
birefringence proposals, see, e.g., heinzl-06 ; Schlenvoigt-2016 . In
addition, a large number $N$ of nuclei represents another enhancement factor
of our set-up. Note that the peak $\propto 1/(\Delta{\bf k})^{2}$ in forward
direction, after integrating over the solid angle $d\Omega$, does only yield a
weak (i.e., logarithmic) enhancement. In addition, as explained after (10),
this $1/(\Delta{\bf k})^{2}$ behavior is only valid up to minimum values of
$|\Delta{\bf k}|$ set by the size of the laser focus. For even smaller
$|\Delta{\bf k}|$, one would have to include the Fourier transform of the
dependence of ${\bf B}_{\rm ext}({\bf r})$ or ${\bf B}_{\rm ext}(t,{\bf r})$
as an effective form factor.

Figure 2:  In an exemplary 3D particle-in-cell (PIC) simulation (PIConGPU), a
cubic carbon cluster target with 100 nm side length and an (initial) electron
density $n_{e}$ of $5\times 10^{23}~{}\rm cm^{-3}$ (i.e., $n_{e}=290~{}n_{c}$)
was irradiated by a short pulse laser with 800 nm wavelength, 30 fs width,
$1.946\times 10^{21}~{}\rm W/cm^{2}$ intensity (i.e., $a_{0}=30$), incident
from the left. In order to assess the number of electrons and ions a
transversely oriented XFEL probe beam would observe, we integrated their
respective densities perpendicular to the plane shown in the lower panels in a
200 nm $\times$ 200 nm square area. The resulting normalized temporal
evolution is plotted in the top: electrons blue, ions orange and the laser
electric field is shown in green for reference. The vertical gray dashed lines
indicate the times for which the lower panels show the electron density and
laser electric field. From top left to bottom right these are $-117$ fs, $-80$
fs, $-43$ fs, $+32$ fs, and $+101$ fs. We see that the electron number in the
interaction volume drops by five orders of magnitude while the ion number only
decreases moderately (due to the Coulomb explosion). The simulation box has
periodic boundaries in the dimensions transverse to the laser; the black scale
bar in the bottom right is 500 nm long.

As one possible scenario (see Fig. 2), one may consider a cubic cluster with
an edge length of 100 nm with typical solid-state density made of carbon, for
example. Then, applying a pre-pulse with a high intensity a bit below $\,{\cal
O}(10^{22}~{}\rm W/cm^{2})$, one may blow out almost all electrons, leaving
behind $N=\,{\cal O}(10^{8})$ ionized nuclei. Shortly afterwards, before the
Coulomb explosion of the remaining nuclei, one would have them interact with
the XFEL superimposed by the main pulse (in the form of a colliding beam set-
up, for example). To exploit the peak in forward direction, let us assume a
small momentum transfer $\Delta\bf{k}$ in the eV regime (corresponding to
scattering angles $\vartheta$ of order millirad). In this case, the amplitudes
from the $N=\,{\cal O}(10^{8})$ nuclei would have basically the same phase and
thus add up coherently, effectively acting as one giant nucleus with charge
$Z_{\rm eff}=NZ$. For much larger $\Delta\bf{k}$, one would have to include
the spatial distribution of the nuclei in analogy to (10).

For smaller clusters (e.g., 50 nm), the electrons would already be blown out
before the laser pulse reaches its peak (see the Appendix). Thus, in this
second scenario, one could avoid the “pre-pulse – main pulse” sequence and use
the same laser pulse for ionization and assisting Delbrück scattering.
However, in this case one should also take into account the electric field
component of the optical laser. This electric field of the laser would then
generate additional contributions in $\delta\epsilon$ and $\delta\mu$ after
combining it with the nuclear Coulomb field. Inserting those additional
contributions in $\delta\epsilon$ and $\delta\mu$ into the full amplitude (9),
we find that they exactly cancel the terms from $\delta\Psi$ if the optical
laser and the XFEL propagate in the same direction. This cancellation does
also occur in the case without the nuclear Coulomb field and demonstrates an
important difference between a propagating plane wave (crossed fields) and a
pure magnetic field. To avoid this cancellation and obtain a birefringent
signal, the XFEL should propagate at a finite angle (e.g., perpendicular) to
the propagation direction as well as the magnetic field component of the
optical laser.

As a third scenario, one could envisage an optical laser focus (expelling the
electrons) co-propagating with the XFEL pulse through a less dense medium, in
analogy to laser wake-field acceleration Esarey2009 . In the usual set-up, the
optical laser would co-propagate with the created fully blown-out plasma
cavity and thus also with the XFEL. However, this would again lead to the
cancellation problem discussed above. To overcome this problem, one could use
an optical laser with a propagation direction different from that of the XFEL,
whose focus is co-propagating with the XFEL pulse. This can be achieved by
laser pulse-front shaping techniques Debus2019 ; Steiniger2019 ; Steiniger2014
; Debus2010 . Inserting typical numbers such as a density of order
$10^{19}~{}\rm cm^{-3}$, the number $N$ of ions within the interaction region,
i.e., the laser focus with a 2.5 $\mu$m spot size, travelling over a distance
of order millimeter, would be about two to three orders of magnitude larger
$N=\,{\cal O}(6\times 10^{10})$ compared to the first scenario. However,
presumably not all of these ions would contribute coherently to the scattering
amplitude in view of the larger spatial extent (again compared to the first
and second scenario). As an advantage, some of the background processes (see
the Appendix) due to the residual electrons and the residual radiation are
minimized in the third scenario. In summary, the three scenarios offer
different advantages and drawbacks, which should be compared carefully for
designing an experimental realization.

Finally, the measurement of the birefringent signal, i.e., the detection of
the XFEL photons with flipped polarization could be achieved in complete
analogy to other vacuum birefringence proposals, see, e.g., heinzl-06 ;
Schlenvoigt-2016 .

#### Conclusions

As an example for the QED vacuum nonlinearity, we calculate the scattering of
XFEL photons at the combined field of a nucleus plus an external magnetic
field (e.g., generated by an optical laser focus), see Fig. 1. In contrast to
previous work involving nuclear fields, such as pure Delbrück scattering, this
scenario yields a large interaction volume, which goes along with a peak of
the differential cross section in forward direction $d\sigma/d\Omega\sim
1/(\Delta{\bf k})^{2}$.

As another distinction, the scales relevant to our scenario are well below the
characteristic scales of QED mentioned in the Introduction, i.e., the critical
field strength (1) and the electron mass $m$. Thus, our scenario is more
within the regime where the approximation of a classical electromagnetic field
applies – especially for the coherent superposition of the signal from many
nuclei – than within the particle (photon) picture often associated with pure
Delbrück scattering.

In addition to the normal polarization conserving scattering (see the
Appendix), we obtain a birefringent signal (12) whose amplitude
${\mathfrak{A}}_{\delta\Psi}^{\perp}$ is just a little bit smaller. Besides
the anisotropy induced by the magnetic field and the aforementioned peak in
forward direction, this birefringence may be used to distinguish the process
of field-assisted Delbrück scattering considered here from other background
processes discussed in the Appendix. Using a large number of XFEL photons
xfel-URL ; xfel-data  and a large number of nuclei, we show that it might be
possible to overcome the suppression of the signal (12) by more than twenty
orders of magnitude such that this effect might be observable in the near
future at facilities such as the Helmholtz International Beamline for Extreme
Fields (HIBEF) at the European XFEL hibef-URL .

###### Acknowledgements.

We would like to thank C. Schubert as well as R. Sauerbrey, R. Shaisultanov,
G. Torgrimsson, and other colleagues from the HZDR for helpful discussions.
This work was partially funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – Project-ID 278162697 – SFB 1242; as well as by
the Center of Advanced Systems Understanding (CASUS) which is financed by
Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon
Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis
of the budget approved by the Saxon State Parliament.

## References

  * (1) H. Euler and B. Kockel, Über die Streuung von Licht an Licht nach der Diracschen Theorie, Naturwissensechaften 23, 246 (1935).
  * (2) W. Heisenberg, H. Euler, Folgerungen aus der Diracschen Theorie des Positrons, Z. Physik 98, 714 (1936).
  * (3) R. Karplus and M. Neuman, The Scattering of Light by Light, Phys. Rev. 83, 776 (1951).
  * (4) F. Sauter, Über das Verhalten eines Elektrons im homogenen elektrischen Feld nach der relativistischen Theorie Diracs, Z. Physik 69, 742 (1931).
  * (5) J. Schwinger, On Gauge Invariance and Vacuum Polarization, Phys. Rev. 82, 664 (1951).
  * (6) R. P. Mignani et al, Evidence for vacuum birefringence from the first optical-polarimetry measurement of the isolated neutron star RX J1856.5-3754, Monthly Notices of the Royal Astronomical Society 465, 492 (2017).
  * (7) L. M. Capparelli et al, A note on polarized light from magnetars, Eur. Phys. J. C 77, 754 (2017).
  * (8) P. Mészáros, High Energy Radiation from Magnetized Neutron Stars, (Chicago Univeristy Press, Chicago, 1992).
  * (9) R. Ruffini, G. Vereshchagin and S.S. Xue, Electron-positron pairs in physics and astrophysics: from heavy nuclei to black holes, Phys. Rept. 487, 1 (2010).
  * (10) L. Meitner and H. Kösters (with addition from M. Delbrück), Über die Streuung kurzwelliger $\gamma$-Strahlen, Zeitschrift für Physik 84, 137 (1933).
  * (11) H. A. Bethe and F. Rohrlich, Small Angle Scattering of Light by a Coulomb Field, Phys. Rev. 86, 10 (1952).
  * (12) V. Costantini, B. De Tollis and G. Pistoni, Nonlinear effects in quantum electrodynamics, Nuovo Cimento Soc. Ital. Fis. A 2, 733 (1971).
  * (13) G. Jarlskog et al, Measurement of Delbrück scattering and observation of photon splitting at high energies, Phys. Rev. D 8, 3813 (1973).
  * (14) H. Cheng and T.T. Wu, High-Energy Elastic Scattering in Quantum Electrodynamics, Phys. Rev. Lett. 22, 666 (1969).
  * (15) H. Cheng and T.T. Wu, High-Energy Collision Processes in Quantum Electrodynamics. I-IV, Phys. Rev. 182, 1852 (1969); 182, 1868 (1969); 182, 1873 (1969); 182, 1899 (1969).
  * (16) A. I. Milstein and V. M. Strakhovneko, Quasiclassical approach to the high-energy Delbrück scattering, Phys. Rev. A 95, 135 (1983).
  * (17) A. I. Milstein and V. M. Strakhovneko, Coherent scattering of high-energy photons in a Coulomb field, Soviet Physics JETP 58, 8 (1983).
  * (18) M. Schumacher, I. Borchert, F. Smend, P. Rullhusen, Delbrück scattering of 2.75 MeV photons by lead, Phys. Lett. B 59 134 (1975).
  * (19) H. Falkenberg et al, Amplitudes for Delbrück scattering, Atomic Data and Nuclear Data Tables 50, 1 (1992).
  * (20) H. E. Jackson and K. J. Wetzel, Delbrück Scattering of 10.8-MeV $\gamma$ Rays, Phy. Rev. Lett. 22, 1008 (1969).
  * (21) R. Moreh and S. Kahana, Delbrück scattering of 7.9 MeV photons, Phys. Lett. B 47, 351 (1973).
  * (22) S. Kahane, O. Shahal, R. Moreh, Rayleigh and Delbrück scattering of 6.8-11.4 MeV photons at $\theta=1.5^{\circ}$, Phy. Lett. B 66, 229 (1977).
  * (23) Sh. Zh. Akhmadaliev et al, Delbrück scattering at energies of 140-450 MeV, Phys. Rev. C 58, 2844 (1998).
  * (24) A. I. Milstein, M. Schumacher, Present status of Delbrück scattering, Phys. Rept. 243 183 (1994).
  * (25) M. Schumacher, Delbrück scattering, Radiation Physics and Chemistry 56, 101 (1999).
  * (26) P. Papatzacos and K. Mork, Delbrück scattering calculations, Phys. Rev. D 12, 206 (1975).
  * (27) P. Papatzacos and K. Mork, Delbrück scattering, Phys. Rept. 21, 81 (1975).
  * (28) P. Rullhusen et al, Test of vacuum polarization by precise investigation of Delbrück scattering, Phys. Rev. C 23, 1375 (1981).
  * (29) ATLAS collaboration, Evidence for light-by-light scattering in heavy-ion collisions with the ATLAS detector at the LHC, Nature Phys. 13, 852 (2017).
  * (30) D. Valle et al, The PVLAS experiment: measuring vacuum birefringence and dichroism with a birefringent Fabry Perot cavity, Eur. Phys. J. C 76, 24 (2016).
  * (31) E. Zavattini et al, Measuring the Magnetic Birefringence of Vacuum: The PVLAS Experiment, Int. J. Mod. Phys. A 27, 1260017 (2012).
  * (32) D. Valle et al, Towards a direct measurement of vacuum magnetic birefringence: PVLAS achievements, Optics Communications 283, 4194 (2010).
  * (33) E. Zavattini et al, New PVLAS results and limits on magnetically induced optical rotation and ellipticity in vacuum, Phys. Rev. D 77, 032006 (2008).
  * (34) E. Zavattini et al, Experimental Observation of Optical Rotation Generated in Vacuum by a Magnetic Field, Phys. Rev. Lett. 96, 110406 (2006).
  * (35) X. Fan et al, The OVAL experiment: a new experiment to measure vacuum magnetic birefringence using high repetition pulsed magnets, Eur. Phys. J. D 71, 308 (2017).
  * (36) M. T. Hartman, R. Battesti and C. Rizzo, Characterization of the Vacuum Birefringence Polarimeter at BMV: Dynamical Cavity Mirror Birefringence, IEEE Transactions on Instrumentation and Measurement 68, 2268 (2019).
  * (37) M. T. Hartman,A. Rivère, R. Battesti and C. Rizzo, Noise characterization for resonantly-enhanced polarimetric vacuummagnetic-birefringence experiments, Review of Scientific Instruments, 88, 123114 (2017).
  * (38) R. Battesti et al, High magnetic fields for fundamental physics, Phys. Rep. 765, 1 (2018).
  * (39) T. Heinzl et al, On the observation of vacuum birefringence, Opt. Commun. 267, 318 (2006).
  * (40) T. Inada et al, Search for photon-photon elastic scattering in the X-ray region, Phys. Lett. B 732, 356 (2014).
  * (41) T. Yamaji et al, An experiment of X-ray photon-photon elastic scattering with a Laue-case beam collider, Phys. Lett. B 763, 454 (2016).
  * (42) H.-P. Schlenvoigt, T. Heinzl, U. Schramm, T. Cowan, R. Sauerbrey, Detecting vacuum birefringence with X-ray free electron lasers and high-power optical lasers: A feasibility study, Phys. Scr. 91, 023010 (2016).
  * (43) F. Karbstein, C. Sundqvist, Probing vacuum birefringence using X-ray free electron and optical high-intensity lasers, Phys. Rev. D 94, 013004 (2016).
  * (44) T. Inada et al, Probing Physics in Vacuum Using an X-ray Free-Electron Laser, a High-Power Laser, and a High-Field Magnet, Appl. Sci. 7, 671 (2017).
  * (45) A. Di Piazza, A.I. Milstein, and C.H. Keitel, Photon splitting in a laser field, Phys. Rev. A 76, 032103 (2007).
  * (46) D. Tommasini et al, Precision tests of QED and non-standard models by searching photon-photon scattering in vacuum with high power lasers, JHEP 11, 043 (2009).
  * (47) D. Tommasini and H. Michinel, Light by light diffraction in vacuum, Phys. Rev. A 82, 011803 (2010).
  * (48) J. K. Koga and T. Hayakawa, Possible Precise Measurement of Delbrück Scattering Using Polarized Photon Beams, Phys. Rev. Lett. 118, 204801 (2017).
  * (49) B. King and Ch. Keitel, Photon-photon scattering in collisions of intense laser pulses, New J. Phys. 14, 103002 (2012).
  * (50) H. Gies, F. Karbstein and N. Seegert, Quantum reflection as a new signature of quantum vacuum nonlinearity, New J. Phys. 15, 083002 (2013).
  * (51) H. Gies, F. Karbstein and C. Kohlfürst, Photon-photon scattering at the high-intensity frontier, Phys. Rev. D 97, 076002 (2018).
  * (52) H. Gies, F. Karbstein and C. Kohlfürst, All-optical signatures of strong-field QED in the vacuum emission picture, Phys. Rev. D 97, 036022 (2018).
  * (53) B. Döbrich and H. Gies, Interferometry of light propagation in pulsed fields, Europhys. Lett. 87, 21002 (2009).
  * (54) H. Grote, On the possibility of vacuum QED measurements with gravitational wave detectors, Phys. Rev. D 91, 022002 (2015).
  * (55) N. Ahmadiniaz, T.E. Cowan, R. Sauerbrey, U. Schramm, H.-P. Schlenvoigt, R. Schützhold, On the Heisenberg limit for detecting vacuum birefringence, Phys. Rev. D 101, 116019 (2020).
  * (56) A. Di Piazza and A. I. Milstein, Delbrück scattering in combined Coulomb and laser fields, Phys. Rev. A 77, 042102 (2008).
  * (57) F. Karbstein, Probing Vacuum Polarization Effects with High-Intensity Lasers, Particles 3, 39 (2020).
  * (58) B. King and T. Heinzl, Measuring Vacuum Polarisation with High Power Lasers, High Power Laser Sci. Eng. 4 (2016).
  * (59) F. Karbstein, The quantum vacuum in electromagnetic fields: From the Heisenberg-Euler effective action to vacuum birefringence, in Proceedings of the Helmholtz International Summer School 2016 (HQ 2016), Dubna, Russia, 18–30 July 2016; pp. 44–57, doi:10.3204/DESY-PROC-2016-04.
  * (60) A. Di Piazza, C. Müller, K. Z. Hatsagortsyan and C.H. Keitel, Extremely high-intensity laser interactions with fundamental quantum systems, Rev. Mod. Phys. 84, 1177 (2012).
  * (61) G. V. Dunne, The Heisenberg-Euler Effective Action: 75 years on, Int. J. Mod. Phys. A 27, 1260004 (2012).
  * (62) Sh. Zh. Akhmadaliev et al, Experimental Investigation of High-Energy Photon Splitting in Atomic Fields, Phys. Rev. Lett. 89, 061802 (2002).
  * (63) Z. Bialynicka-Birula, I. Bialynicki-Birula, Nonlinear effects in Quantum Electrodynamics. Photon propagation and photon splitting in an external field, Phys. Rev. D 2, 2341 (1970).
  * (64) S. L. Adler, J. N. Bahcall, C. G. Callan and M.N. Rosenbluth, Photon splitting in a strong magnetic field, Phys. Rev. Lett. 25, 1061 (1970).
  * (65) S. L. Adler, Photon splitting and photon dispersion in a strong magnetic field, Ann. of Phys. 67, 599 (1971).
  * (66) S. L. Adler and C. Schubert, Photon splitting in a strong magnetic field: Recalculation and comparison with previous calculations, Phys. Rev. Lett. 77, 1695 (1996).
  * (67) J. S. Toll, The Dispersion Relation for Light and Its Application to Problems Involving Electron Pairs, Ph.D. Thesis, Princeton University, Princeton, NJ, USA, 1952, unpublished.
  * (68) One could derive the same result from QED Feynman diagrams such as in Fig. 1. Since three vertices (the two XFEL photons and the external magnetic field) correspond to low energies and momenta (i.e., well below the electron mass $m$), energy-momentum conservation implies that the fourth vertex (representing the Coulomb field) does also involve low energies and momenta. In fact, for the case of forward scattering considered here, these energies and momenta are extremely low. Thus, we may use the low-energy limit of the QED Feynman diagrams, which are equivalent to the Euler-Heisenberg Lagrangian (2). Note that the situation is different for pure Delbrück scattering where only two vertices (the two XFEL photons) correspond to low energies and momenta such that the other two (representing the Coulomb field) can involve high momenta.
  * (69) D. M. Volkov, Über eine Klasse von Lösungen der Diracschen Gleichung, Z. Physik 94, 250 (1935).
  * (70) P. J. Redmond, Solution of the Klein-Gordon and Dirac Equations for a Particle with a Plane Electromagnetic Wave and a Parallel Magnetic Field, J. Math. Phys. 6, 1163 (1965).
  * (71) A. I. Milstein, I. S. Terekhov, U. D. Jentschura, and C. H. Keitel, Laser-dressed vacuum polarization in a Coulomb field, Phys. Rev. A 72, 052104 (2005).
  * (72) A. Di Piazza, K. Z. Hatsagortsyan and C. H. Keitel, Nonperturbative Vacuum-Polarization Effects in Proton-Laser Collisison, Phys. Rev. Lett. 100, 010403 (2008).
  * (73) A. Di Piazza and A. I. Milstein, Quasiclassical approach to high-energy QED processes in strong laser and atomic fields, Phys. Lett. B 717, 224 (2012).
  * (74) A. Di Piazza and A. I. Milstein, Ultrarelativistic quasiclassical wave functions in strong laser and atomic fields, Phys. Rev. A 89, 062114 (2014).
  * (75) J. D. Jackson, Classical electrodynamics, John Wiley & Sons, 2007.
  * (76) T. Tschentscher et al, Photon Beam Transport and Scientific Instruments at the European XFEL, Appl. Sci. 7, 592 (2017).
  * (77) https://www.xfel.eu
  * (78) E. Esarey, C.B. Schroeder, W.P. Leemans, Physics of laser-driven plasma-based electron accelerators, Rev. Mod. Phys. 81, 1229 (2009).
  * (79) A. Debus et al, Traveling-wave Thomson scattering and optical undulators for high-yield EUV and X-ray sources, App. Phys. B 100, 61 (2010).
  * (80) K. Steiniger et al, Optical free-electron lasers with Traveling-Wave Thomson-Scattering, J. Phys. B: Atomic, Molecular and Optical Physics 47, 234011 (2014).
  * (81) A. Debus et al, Circumventing the Dephasing and Depletion Limits of Laser-Wakefield Acceleration, Phys. Rev. X 9, 031044 (2019).
  * (82) K. Steiniger et al, Building an Optical Free-Electron Laser in the Traveling-Wave Thomson-Scattering Geometry, Front. Phys 6, 155 (2019).
  * (83) http://www.hibef.eu
  * (84) G. Baldwin and G. Klaiber, Photo-Fission in Heavy Elements, Phys. Rev. 71, 3 (1947).
  * (85) This estimate can also be obtained from the lowest-order Feynman diagram. The four vertices of the fermion box yield $q^{4}$ and the two external Coulomb field lines contribute $(Zq)^{2}$. Furthermore, scaling everything with respect to the electron mass (as the natural QED scale), the two external XFEL photon lines scale with $\omega^{2}/m^{2}$.
  * (86) O. Klein and Y. Nishina, Über die Streuung von Strahlung durch freie Elektronen nach der neuen relativistischen Quantendynamik von Dirac, Zeitschrift für Phys. 52, 853 (1929).
  * (87) O. Klein and Y. Nishina, The Scattering of Light by Free Electrons according to Dirac’s New Relativistic Dynamics, Nature 122, 398 (1928).
  * (88) A. Wightman, Note on polarization effects in Compton scattering, Phys. Rev. 74, 1813 (1948).
  * (89) W. Heitler, The quantum theory of radiation, Courier Corporation, 1984.
  * (90) R. Schützhold, H. Gies and G. Dunne, Dynamically assisted Schwinger mechanism, Phys. Rev. Lett. 101, 130404 (2008).
  * (91) G. Torgrimsson, C. Schneider and R. Schützhold, Sauter-Schwinger pair creation dynamically assisted by a plane wave, Phys. Rev. D 97, 096004 (2018).
  * (92) M. Bussmann et al, Radiative signature of the relativistic Kelvin-Helmholtz Instability, Proc. Int. Conf. High Perform. Comput. Networking, Storage Anal. - SC 13, 1 (2013).
  * (93) H. Burau et al, PIConGPU: a fully relativistic particle-in-cell code for a GPU cluster, IEEE Transactions on Plasma Science 38, 2831 (2010).
  * (94) A. Huebl et al, PIConGPU 0.4.3: System Updates and Bug Fixes, 10.5281/zenodo.2565503 (2020).
  * (95) A. Huebl et al, ComputationalRadiationPhysics/picongpu: Perfectly Matched Layer (PML) and Bug Fixes, 10.5281/zenodo.3875374 (2020).
  * (96) K. Yee, Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media, IEEE Transactions on Antennas and Propagation, 14, 302 (1966).
  * (97) J. P. Boris, Relativistic plasma simulation-optimization of a hybrid code, Proc. Fourth Conf. Num. Sim. Plasmas, 3-67 (1970).
  * (98) T. Zh. Esirkepov, Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor, Computer Physics Communications 135,144 (2001).
  * (99) R. W. Hockney and J. W. Eastwood. Computer simulation using particles, crc Press, 1988.
  * (100) N. B. Delone and V. P. Krainov,Tunneling and barrier-suppression ionization of atoms and ions in a laser radiation field, Phys. Usp. 41, 469 (1998).
  * (101) E. Clementi and D. Raimondi, Atomic Screening Constant from SCF Functions, The Journal of Chemical Physics 38, 2686 (1963).
  * (102) E. Clementi and D. Raimondi, Atomic Screening Constant from SCF Functions. II. Atoms with 37 to 86 Electrons, The Journal of Chemical Physics 47, 1300 (1967).
  * (103) R. M. More, Pressure Ionization, Resonances, and the Continuity of Bound and Free States, Advances in Atomic, Molecular and Optical Physics 21 C, 305 (1985).

##  I Appendix: Backward Scattering

In contrast to the case of forward scattering discussed above, let us study
the opposite case of backward scattering ${\bf n}_{\rm in}=-{\bf n}_{\rm out}$
for comparison. Again focusing on the birefringent signal ${\bf e}_{\rm
in}\perp{\bf e}_{\rm out}$, we find

$\displaystyle{\mathfrak{A}}_{\delta\Psi}^{\perp}=i\,\frac{7\alpha_{\rm
QED}^{2}}{45\pi}\,\frac{q{\bf B}_{\rm ext}\cdot{\bf n}_{\rm
out}}{m^{2}}\,\frac{Z\omega}{m^{2}}\,\,.$ (14)

In contrast to forward scattering, here the isotropic part
$14\xi\delta^{ij}({\bf E}_{\rm ext}\cdot{\bf B}_{\rm ext})$ in (5) contributes
while the anisotropic terms cancel. As a consequence, the maximum signal is
obtained if the magnetic field ${\bf B}_{\rm ext}$ is aligned with the XFEL
propagation direction ${\bf n}_{\rm in}$ and thus perpendicular to both ${\bf
e}_{\rm in}$ and ${\bf e}_{\rm out}$.

##  II Appendix: Parallel Polarizations

To complete our analysis, we also consider the case of parallel polarizations,
even though this signal is probably much harder to distinguish from the
background processes (see below). Interestingly, for backward scattering ${\bf
n}_{\rm in}=-{\bf n}_{\rm out}$, the field-assisted contribution from
$\delta\Psi$ vanishes. In forward direction, we get a $1/|\Delta{\bf k}|$ peak
in analogy to the birefringent signal (11)

$\displaystyle{\mathfrak{A}}_{\delta\Psi}^{\|}$ $\displaystyle=$
$\displaystyle i\xi\,\frac{Q}{(\Delta{\bf
k})^{2}}\,\frac{\omega^{2}}{4\pi}\left[28({\bf e}\cdot{\bf B}_{\rm ext})({\bf
n}\times{\bf e})\cdot\Delta{\bf k}\right.$ (15) $\displaystyle\left.-16({\bf
e}\cdot\Delta{\bf k})({\bf n}\times{\bf e})\cdot{\bf B}_{\rm ext}\right]\,,$

where we have set ${\bf e}_{\rm in}$ and ${\bf e}_{\rm out}$ to the overall
XFEL polarization unit vector ${\bf e}$. For both polarizations (perpendicular
and parallel), the leading-order amplitudes (11) and (15) in forward direction
vanish identically if the XFEL propagation direction $\bf n$ is parallel to
the magnetic field ${\bf B}_{\rm ext}$. In contrast to the birefringent signal
(11), the above amplitude (15) assumes its maximum contribution if the
polarization $\bf e$ (i.e., the electric field) of the XFEL is parallel to the
magnetic field ${\bf B}_{\rm ext}$ (and thus the XFEL propagation direction
$\bf n$ perpendicular to ${\bf B}_{\rm ext}$), while the momentum transfer
$\Delta{\bf k}$ is perpendicular to both, i.e., in the “equatorial plane”.

##  III Appendix: Background Processes

Of course, finding such small differential cross sections (13) necessitates
the careful estimate of potential background processes.

###  III.1 Nuclear Thomson Scattering

Perhaps the most obvious candidate is ordinary Thomson scattering of the XFEL
photons by the nuclei themselves. This process is suppressed for heavy nuclei
as its amplitude scales inversely proportional to the mass $M$ of the nucleus

$\displaystyle{\mathfrak{A}}_{\rm T}=\alpha_{\rm QED}\,\frac{Z^{2}}{M}\,.$
(16)

Inserting a typical mass of 0.9 GeV per nucleon, we find that
${\mathfrak{A}}_{\rm T}$ actually becomes comparable to the parallel
polarization signal ${\mathfrak{A}}_{\delta\Psi}^{\|}$ in (15) for a momentum
transfer in the eV regime. However, Thomson scattering ${\mathfrak{A}}_{\rm
T}$ does not display the characteristic (anisotropic) dependence on
$\Delta{\bf k}$ and ${\bf B}_{\rm ext}$, for example, and does also not
generate a birefringent signal.

###  III.2 Nuclear Resonances

Similar arguments apply to the scattering of photons by the nuclei via
coupling to their internal degrees of freedom, usually referred to as nuclear
resonances. Of course, this process does depend on the concrete nuclear
structure, such as the nuclear resonance frequencies, which are typically way
above the XFEL frequency. As detailed in rev4 , for example, the giant dipole
resonance (see, e.g., baldwin-47 ) dominates in such a case – but it is still
negligible in comparison to ${\mathfrak{A}}_{\rm T}$ for energies well below
the MeV regime. Furthermore, these nuclear resonances do not generate a
birefringent signal to lowest order rev4 .

###  III.3 Delbrück Scattering

Of course, another background process is pure Delbrück scattering, i.e., the
interaction of the XFEL beam with the nuclear Coulomb field only. If we
estimate the amplitude with the approach described above (in the absence of
the optical laser), we find that $\delta\Psi$ vanishes while $\delta\epsilon$
and $\delta\mu$ scale as $1/r^{4}$. Thus, in contrast to
${\mathfrak{A}}_{\delta\Psi}$ considered here, the effective interaction
volume would not be large and the spatial integral would be dominated at small
(instead if large) radii. Formally, this integral would even diverge for $r\to
0$, but this is an artifact of the used approach based on the leading-order
Euler-Heisenberg Lagriangian (2), which is no longer applicable for very small
radii since the Coulomb field ${\bf E}_{\rm ext}({\bf r})$ and its gradients
become large. Even though it would be a very interesting and non-trivial
problem to study Delbrück scattering in this regime, let us obtain a very
rough estimate foot2  by cutting off the $r$-integral at the Compton length
$\lambdabar$

$\displaystyle{\mathfrak{A}}_{\rm D}\sim Z^{2}\frac{\alpha_{\rm
QED}^{3}}{m}\,\frac{\omega^{2}}{m^{2}}\,.$ (17)

Thus, pure Delbrück scattering is negligible in comparison to the signal
${\mathfrak{A}}_{\delta\Psi}$ considered here.

###  III.4 Magnetic Field

Going to the other limit, the interaction of the XFEL beam with an optical
laser focus only (i.e., without the nucleus) has already been considered in
several works, see, e.g., heinzl-06 ; tommasini-09 ; tommasini-10 ; king-12 .
Since a constant ${\bf B}_{\rm ext}$ would generate constant $\delta\epsilon$
and $\delta\mu$ such that the spatial integral would yield
$\delta^{3}(\Delta{\bf k})$, the possible momentum transfer $\Delta{\bf k}$ is
strongly limited by the spatial and temporal inhomogeneities of the optical
laser focus, see also dipiazza-08 ; milstein-05 ; dipiazza-07 ; dipiazza-12 ;
dipiazza-14 . As another option for discriminating, this process would be
independent of the charge $Z$ and the number $N$ of nuclei. Note that the
optical laser focus can also generate a birefringent signal.

###  III.5 Electronic Compton Scattering

Another important background process is Compton scattering from the residual
electrons. According to the particle-in-cell (PIC) simulations displayed in
Fig. 2, for example, the electron number in the interaction volume drops by a
factor of $10^{-5}$ or more, but there may still be a few $\,{\cal O}(10^{4})$
electrons left. Thus, even assuming that the amplitudes from the electrons add
up coherently, their lowest-order Thomson scattering signal would already be
below the nuclear Thomson scattering ${\mathfrak{A}}_{\rm T}$ as well as the
signal ${\mathfrak{A}}_{\delta\Psi}$ considered here. However, this Thomson
scattering ${\mathfrak{A}}_{\rm T}=\alpha_{\rm QED}/m$ from the residual
electrons does not generate a birefringent signal to lowest order. A
birefringent signal can be generated by quantum effects due to the electron
spin or by coupling the electrons to the external magnetic field. However,
these higher-order processes are suppressed in comparison the the lowest-order
Thomson scattering ${\mathfrak{A}}_{\rm T}=\alpha_{\rm QED}/m$.

For example, the birefringent signal from free electrons can be estimated from
the Klein-Nishina formula KN1 ; KN2 ; Wig ; Hei

$\displaystyle\frac{d\sigma_{\rm KN}^{\perp}}{d\Omega}=\frac{\alpha_{\rm
QED}^{2}}{m^{2}}\,\frac{\omega^{2}}{m^{2}}\,\frac{\vartheta^{4}}{16}\,,$ (18)

where $\vartheta\ll 1$ denotes the small scattering angle. The birefringent
signal induced by coupling to the external magnetic field can be estimated by
the insertion of an additional vertex into the lowest-order Feynman diagram
for Compton scattering. The resulting amplitude is suppressed by the inverse
combined Keldysh parameter ass1 ; ass2

$\displaystyle\frac{1}{\gamma}=\frac{qB_{\rm ext}}{m\omega}\,,$ (19)

which is on the percent level for our parameters. Thus, these effects are
sufficiently suppressed in comparison with the signal
${\mathfrak{A}}_{\delta\Psi}^{\perp}$ considered here – unless we have too
many residual electrons.

Even for a larger number of residual electrons, one could still observe the
desired signal ${\mathfrak{A}}_{\delta\Psi}^{\perp}$ if the unwanted
contributions from the electrons would not add up coherently. This would be
quite natural since the light electrons are much more prone to de-phasing and
decoherence than the heavy nuclei. They experience a stronger recoil and
typically have a broader distribution in velocity and position. Further, if
the birefringent amplitude depends on the initial electron spin state of the
(unpolarized) electrons or if their spin state is changed by the birefringent
scattering event, one would get an incoherent superposition.

###  III.6 Residual Radiation

Further sources for x-ray photons which might lead to a false signal could be
radiation processes, such as from the electrons accelerated by the optical
laser during ionization or bremsstrahlung effects etc. However, these photons
will typically have frequencies different from the XFEL frequency and thus can
be filtered out efficiently by the very sharp energy resolution (in the sub-eV
regime) of the x-ray polarization filters based on multiple Bragg scattering.

###  III.7 Interferences

Note that the amplitude ${\mathfrak{A}}_{\delta\Psi}$ derived here is
imaginary. This is related to the parity $\mathcal{P}$ symmetry of QED as the
Coulomb field ${\bf E}_{\rm ext}({\bf r})$, the symmetry-breaking tensor
$\delta\Psi({\bf r})$, and thus also the effective current ${\bf
J}_{\omega}^{\rm eff}$ generating ${\mathfrak{A}}_{\delta\Psi}$ are all odd
functions of space. In contrast, $\delta\epsilon$ and $\delta\mu$ are even
functions, such that the amplitude for pure Delbrück scattering
${\mathfrak{A}}_{\rm D}$ is real – as is the amplitude for Thomson scattering
${\mathfrak{A}}_{\rm T}$. As a result, there are no interference terms (to
leading order) between the field-assisted Delbrück scattering
${\mathfrak{A}}_{\delta\Psi}$ considered here and Thomson scattering
${\mathfrak{A}}_{\rm T}$ or Delbrück scattering ${\mathfrak{A}}_{\rm D}$,
while there can be interference terms between the latter two
${\mathfrak{A}}_{\rm T}$ and ${\mathfrak{A}}_{\rm D}$.

##  IV Appendix: PIC simulations

The 3D particle-in-cell simulations shown in Fig. 2 and supplementary Figs. 3
and 4 were performed using the code PIConGPU buss1 ; buss2 , version $0.4.3$
buss3  (Figs. 2 and 3) $0.4.2$ buss4  (Fig. 4). For Figs. 2 and 3, the grid
size is $640\times 1280\times 640$, the spatial resolution is $256$ cells per
$\lambda_{0}$ ($\lambda_{0}=800$ nm) and the temporal resolution is $6\times
10^{-3}$ fs. For Fig. 4 we employed a moving-window frame with a total size of
$256\times 2560\times 1024$ cells which propagates for $50,000$ iterations.
The spatial resolution here is $4.5/\lambda_{0}\times 18/\lambda_{0}\times
18/\lambda_{0}$ with a temporal resolution of $103$ as. The electromagnetic
field evolution is simulated with the Yee solver yee , while the particle
motion is computed using the Boris pusher boris . Particles influence the
fields via the Esirkepov current deposition scheme esi  with a TSC macro-
particle shape hock . Ionization was taken into account for Figs. 2 and 3 by
field ionization via the tunneling and barrier suppression del  with orbital
structure in shielding of the ion charge according to cle1 ; cle2 , as well as
by collisions via a modified Thomas-Fermi model more .

###  IV.1 $50~{}\rm nm$ Carbon Cluster

In analogy to Fig. 2, the PIC simulation for a cubic carbon cluster target at
50 nm side length is shown in Fig. 3. For this smaller cluster, the number of
electrons within the interaction volume drops by more than four orders of
magnitude before the peak of the laser pulse is reached. In addition, the
effect of the Coulomb explosion of the remaining ions is weaker (as expected).

Figure 3: Same as in Fig. 2, but for a cubic carbon cluster target at 50 nm
side length.

###  IV.2 Co-moving Laser Focus

Laser-wakefield acceleration Esarey2009  can sustain a cavity with full
electron blowout until reaching laser pump depletion. If an XFEL laser
propagates in this blowout in the region for the entire duration where the
laser is most intense, electron background is minimal, so that background
processes are minimized. For such an approach, lower densities of several
$10^{19}~{}\rm cm^{-3}$ are required, because the drive laser pulse dimensions
have to be smaller than the plasma wavelength $2\pi c/\omega_{p}$. Compared to
high-density targets, see Figs. 2 and 3, this reduction in density is
compensated by a much longer interaction length, which can be in the mm-range.
However for laser-wakefield acceleration, the drive laser propagation needs to
be parallel to XFEL propagation, greatly diminishing the cross-section of
laser-assisted Delbrück scattering, so that a different scattering geometry is
chosen.

A simplified version of Traveling-wave electron acceleration (TWEAC) Debus2019
, i.e., using one drive laser only, provides a focal region co-moving at the
vacuum speed of light with the XFEL laser. Locally in the interaction region
of the XFEL laser, the ions see an incoming plane wave coming at some incident
angle $\phi$. Globally, the laser is cylindrically focused to a line collinear
with the XFEL direction of propagation. Also, in order to maintain a
stationary overlap with the XFEL beam, the laser pulse shape has to be
prepared to feature a pulse-front tilt of $\phi/2$ Steiniger2019 .

The example shown in Fig. 4 is based on a 3D-PIC simulation, displaying the
electron density time-averaged over the full 1 mm interaction length, showing
a sizable region devoid of residual electrons (at least 3 orders of magnitudes
electron density reduction) in the vicinity of the XFEL pulse and close to
maximum laser intensity.

Neither the depletion limit, nor laser pulse guiding constraints limit
aforementioned geometry and interaction length, but the laser pulse energy
available. For intensities at several $10^{21}~{}\rm W/cm^{2}$ ($a_{0}=30$),
laser pulse energies only remain in a technically feasible range if the
incidence angle does not become too large and features a tight focus.
Furthermore, for minimal differences between initial onset of interaction and
the later stationary interaction conditions, it is useful to choose an
interaction region close to a steep edge of the gas target. This facilitates
full plasma blowout around an XFEL laser over the entire interaction with a
gas target.

Figure 4: In a 3D-PIC simulation (_PIConGPU_), an electron blow-out cavity is
driven in a 1 mm long $H_{2}$ plasma target at $10^{19}~{}\rm cm^{-3}$ (black-
dotted on density scale) using a cylindrically focused, pulse-front tilted
laser (800 nm, 30 fs, 335 J, $1.9\times 10^{21}~{}\rm W/cm^{2}$, $a_{0}=30$,
2.5 $\mu$m focus height and 1 mm interaction length) incident at $10^{\circ}$
with respect to the XFEL direction of propagation and $5^{\circ}$ pulse-front
tilt. For assessing the amount of residual electrons over the entire
interaction length the electron density is time-averaged over the full 1 mm
interaction length, showing a sizable region cleared of electrons in the
vicinity of the XFEL pulse and close to maximum laser intensity.


# Permutational symmetry for identical multi-level systems: a second quantized
approach

Rui E. F. Silva<EMAIL_ADDRESS>Instituto de Ciencia de Materiales de
Madrid, Consejo Superior de Investigaciones Científicas (ICMM-CSIC), Madrid,
Spain  Departamento de Física Teórica de la Materia Condensada and Condensed
Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, E-28049
Madrid, Spain  Johannes Feist<EMAIL_ADDRESS>Departamento de Física
Teórica de la Materia Condensada and Condensed Matter Physics Center (IFIMAC),
Universidad Autónoma de Madrid, E-28049 Madrid, Spain

###### Abstract

We develop a framework that provides a straightforward approach to fully
exploit the permutational symmetry of identical multi-level systems. By taking
into account the permutational symmetry, we outline a simple scheme that
allows to map the dynamics of $N$ identical $d$-level systems to the dynamics
of $d$ bosonic modes with $N$ particles, achieving an exponential reduction on
the dimensionality of the problem in a simple and straightforward way. In
particular, we consider the Lindblad dynamics of several identical multi-level
systems interacting with a common subsystem under the action of collective
dissipation terms.

##  I Introduction

When dealing with a collection of $N$ $d$-level systems, a well-known problem
is the so-called curse of dimensionality, i.e., the fact that the dimension of
the Hilbert space scales exponentially as $d^{N}$. However, in many different
physical phenomena, such as lasing (Scully and Zubairy, 1997; Richter  _et
al._ , 2015), phase transitions (Wang and Hioe, 1973; Walls  _et al._ , 1978;
Gegg  _et al._ , 2018), superradiance (Dicke, 1954; Garraway, 2011), strong
coupling with organic molecules (Herrera and Spano, 2016; Feist  _et al._ ,
2018) and microwave photonics (Putz  _et al._ , 2017), the theoretical
modeling usually assumes that the emitters are identical. In these situations,
permutational symmetry of the $N$ $d$-level systems can be used to greatly
reduce the complexity of the problem. This was addressed in the works of Gegg
_et al._ (Gegg and Richter, 2016, 2017; Gegg, 2017), Shammah _et al._ (Shammah
_et al._ , 2018) and Kirton _et al._ (Kirton and Keeling, 2017, 2018). In
these works, the dynamics of an open quantum system composed of several
identical emitters interacting with a common subsystem under the action of
individual, but identical, collapse operators is considered. By exploiting the
permutational symmetry of the density matrix in the symmetrized Liouville
space, a huge reduction in the complexity of the problem is achieved, allowing
calculations for larger numbers of emitters than possible otherwise. These
efforts were conducted for an ensemble of multi-level systems (Gegg and
Richter, 2016, 2017; Gegg, 2017) and specialized for the case of two-level
systems (Shammah  _et al._ , 2018; Kirton and Keeling, 2017, 2018). In the
absence of individual dephasing operators and for appropriate initial states,
one can further restrict the Hilbert space to the totally symmetric subspace
(Gegg and Richter, 2016; Gegg, 2017). In the case of $2$-level systems, the
construction of the totally symmetric subspace can be achieved by using the
Dicke basis, restricting the Hilbert space to the highest super-spin subspace
(Dicke, 1954; Garraway, 2011).

In this work, we notice that working in the totally symmetric subspace is
completely equivalent to restricting the possible states to bosonic many-body
states. Therefore, by applying the rules of second quantization for bosons, we
achieve the reduction to the totally symmetric subspace in a simple and
straightforward way, mapping the dynamics of $N$ identical $d$-level systems
to the dynamics of $d$ bosonic modes with $N$ particles.

##  II Theory

Figure 1: (a,b) Schematic representation of two product states of four 3-level
systems. (a)
$\left|1\right\rangle_{1}\left|2\right\rangle_{2}\left|3\right\rangle_{3}\left|2\right\rangle_{4}$
and (b)
$\left|2\right\rangle_{1}\left|3\right\rangle_{2}\left|1\right\rangle_{3}\left|2\right\rangle_{4}$.
If we restrict dynamics to the totally symmetric subspace, these two states
will be always in a symmetric superposition that can be correctly taken into
account by using many-body bosonic states, (c)
$\frac{1}{\sqrt{2!}}b_{3}^{\dagger}b_{2}^{\dagger}b_{2}^{\dagger}b_{1}^{\dagger}\left|\mathrm{vac}\right\rangle$.

We start by considering the dynamics of $N$ identical $d$-level systems
interacting with a common subsystem and under the action of collective
dissipation terms, described by the Lindblad master equation,

$\dot{\rho}=-i\left[H,\rho\right]+\sum_{i}\mathcal{L}_{C_{i}}\left[\rho\right],$
(1)

where $H$ is the Hamiltonian, possibly time-dependent, and
$\mathcal{L}_{C_{i}}\left[\rho\right]=C_{i}\rho
C_{i}^{\dagger}-\frac{1}{2}\left(C_{i}^{\dagger}C_{i}\rho+\rho
C_{i}^{\dagger}C_{i}\right)$ is the Lindblad dissipator for the collapse
operator $C_{i}$. In the following, we assume the Hamiltonian, $H$, to be
invariant under any permutation of the $d$-level systems. We also restrict the
collapse operators, $C_{i}$, to collective operators that are also invariant
under any permutation of the $d$-level systems. For the case where the
collapse operators may act locally on each $d$-level system, we must construct
and work on the symmetrized Liouville space and this was taken into account in
(Gegg and Richter, 2016, 2017; Gegg, 2017; Shammah  _et al._ , 2018; Kirton
and Keeling, 2017, 2018). In this work, we restrict to the case where both the
Hamiltonian and the collapse operators are invariant under any permutation of
the emitters.

We may define the symmetrization operator,

$S=\frac{1}{N!}\sum_{\pi}P_{\pi}$ (2)

where $P_{\pi}$ is a permutation operator and $\pi$ runs over all possible
permutations of the $d$-level systems. A permutationally invariant operator,
$O$, then satisfies $[S,O]=0$. We also assume that the initial state,
$\rho_{0}=\sum_{i}p_{i}\left|\psi_{i}\right\rangle\left\langle\psi_{i}\right|$,
is a totally symmetric state, i.e.,
$S\left|\psi_{i}\right\rangle=\left|\psi_{i}\right\rangle$ for all
$\left|\psi_{i}\right\rangle$. Since both the Hamiltonian, $H$, and the
collapse operators, $C_{i}$, are permutationally invariant operators, when
solving the Lindblad dynamics, the density matrix will always remain in the
totally symmetric subspace. This may be used to substantially reduce the
dimensionality of the problem. As already noticed in (Gegg and Richter, 2016;
Gegg, 2017), this reduction of dimensionality is even larger than the one
obtained by working in the symmetrized Liouville space. We must stress that
the difference in the restrictions imposed in this work and the symmetrized
Liouville space approach is that here, each collapse operator $C_{i}$ must be
permutationally invariant, whereas in Refs. (Gegg and Richter, 2016, 2017;
Gegg, 2017; Shammah  _et al._ , 2018; Kirton and Keeling, 2017, 2018) only the
sum of all collapse operators, $\sum_{i}C_{i}$, must be permutationally
invariant.

We can write any permutationally invariant $M$-body operator acting solely on
the emitters as

$O_{\mathrm{em}}^{M}=\frac{1}{M!}\sideset{}{{}^{\prime}}{\sum}_{i_{1}\ldots
i_{M}}^{N}\sum_{\begin{subarray}{c}\alpha_{1}\ldots\alpha_{M}\\\
\beta_{1}\ldots\beta_{M}\end{subarray}}^{d}V_{\alpha_{1}\ldots\alpha_{M}}^{\beta_{1}\ldots\beta_{M}}\sigma_{\beta_{1},\alpha_{1}}^{i_{1}}\ldots\sigma_{\beta_{M},\alpha_{M}}^{i_{M}}$
(3)

where
$\sigma_{\beta,\alpha}^{j}=\left|\beta\right\rangle_{j}\left\langle\alpha\right|_{j}$
is an operator acting on emitter $j$, and the primed sum indicates that all
indices $i_{1},\ldots,i_{M}$ have to be distinct. Using this definition, any
permutationally invariant operator acting on the emitters and containing up to
$M$-body terms can be written as

$O_{\mathrm{em}}=\sum_{J=1}^{M}O_{\mathrm{em}}^{J}.$ (4)

In general, a permutationally invariant operator that may act on the emitters
and on a common subsystem can be written as the sum of three terms

$O=O_{\mathrm{em}}+O_{\mathrm{em-sub}}+O_{\mathrm{sub}},$ (5)

where $O_{\mathrm{em}}$ ($O_{\mathrm{sub}}$) is an operator acting only on the
emitters (subsystem). The interaction term, $O_{\mathrm{em-sub}}$, can be
written as $O_{\mathrm{em-sub}}=\sum_{q}A_{q}B_{q}$, where $A_{q}$ acts solely
on the emitters and must have the form of Eq. (missing) 4 and $B_{q}$ is an
operator acting on the common subsystem.

At this point, one may realize that if the state of the system is restricted
the totally symmetric Hilbert space, one can take advantage of all the
formalism of second quantization for bosons, for which the many-body states
are automatically restricted to the totally symmetric Hilbert space. This can
be done by applying the rules of second quantization (Bruus and Flensberg,
2004) and mapping all the operators to a second quantized form. For instance,
Eq. (missing) 3 becomes

$O_{\mathrm{em}}^{M}=\frac{1}{M!}\sum_{\begin{subarray}{c}\alpha_{1}\ldots\alpha_{M}\\\
\beta_{1}\ldots\beta_{M}\end{subarray}}^{d}V_{\alpha_{1}\ldots\alpha_{M}}^{\beta_{1}\ldots\beta_{M}}b_{\beta_{1}}^{\dagger}\ldots
b_{\beta_{M}}^{\dagger}b_{\alpha_{1}}\ldots b_{\alpha_{M}}$ (6)

where $b_{\alpha}^{\dagger}$ and $b_{\alpha}$ are the bosonic creation and
annihilation operators of an emitter in state $\alpha$. The recipe then simply
consists in applying these rules to all the relevant operators, i.e., the
Hamiltonian, the collapse operators and all the desired observables, giving a
Hamiltonian that can be easily implemented with standard quantum optics
packages such as QuTiP Johansson  _et al._ (2012). For the typical cases where
the number of emitters $N$ is fixed, one can restrict the Hilbert space to the
$N$-particle subspace, for which
$\left\langle\sum_{\alpha=1}^{d}b_{\alpha}^{\dagger}b_{\alpha}\right\rangle=N$.
Similarly, for the initial state,
$\rho_{0}=\sum_{i}p_{i}\left|\psi_{i}\right\rangle\left\langle\psi_{i}\right|$,
each $\left|\psi_{i}\right\rangle$ has to be mapped to its second quantized
version, i.e., expressed in the Fock space.

This approach thus solves the dynamics of $N$ permutationally invariant
$d$-level systems by treating it as the dynamics of $N$ bosons in a system
with $d$ modes. It has the usual advantages of a second quantized formulation.
In particular, it is not necessary to explicitly construct a totally symmetric
subspace, and the correct symmetry enhancement factors are automatically
encoded within and obtained from the bosonic operator algebra.

##  III Results

The study of few-level emitters interacting with light is at the core of our
understanding of light-matter interaction. In principle, the problem of light-
matter interaction can be fully understood within the laws of quantum
electrodynamics (QED) (Cohen-Tannoudji  _et al._ , 1987). However, for
practical applications in fields such as cavity QED, quantum optics, quantum
nanophotonics, and quantum plasmonics, a very common assumption is that matter
degrees of freedom can be described using only a few levels and that the
interaction with light is dominated by a single mode of the electromagnetic
field. In the case of two-level systems, this leads to the well known Rabi
(Rabi, 1937), Dicke (Dicke, 1954), Jaynes-Cummings (Jaynes and Cummings, 1963)
and Tavis-Cummings (Tavis and Cummings, 1968) models, which differ in the
number of emitters and the use of the rotating wave approximation.
Nevertheless, even when dealing with an ensemble of few-level systems, the
exponential scaling of the Hilbert space dimension quickly makes the problem
intractable and symmetry considerations must be taken into account to reduce
the dimension of the problem. In the context of the interaction of an ensemble
of identical two-level systems with a cavity mode, superradiance can be
observed. To study superradiance, it is useful to rewrite the Hamiltonian
using spin operators and use the so-called Dicke basis (Garraway, 2011). In
the following, we demonstrate our approach for several examples within this
context.

###  III.1 Tavis-Cummings model

As a first example to illustrate our approach, we will apply it to the Tavis-
Cummings model (Tavis and Cummings, 1968). For this relatively simple example,
we show explicitly that the state space and matrix elements within the second
quantized picture are the same as in conventional approaches. The Tavis-
Cummings Hamiltonian is given by

$\displaystyle H_{\mathrm{TC}}$
$\displaystyle=\sum_{j=1}^{N}\frac{\omega_{0}}{2}\left(\sigma_{\mathrm{e},\mathrm{e}}^{j}-\sigma_{\mathrm{g},\mathrm{g}}^{j}\right)+\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}$
$\displaystyle+g\sum_{j=1}^{N}\left(a_{\mathrm{c}}\sigma_{\mathrm{e},\mathrm{g}}^{j}+a_{\mathrm{c}}^{\dagger}\sigma_{\mathrm{g},\mathrm{e}}^{j}\right)$
(7)

where $\mathrm{g}$ ($\mathrm{e}$) stands for the ground (excited) state, $N$
is the number of two-level systems, $\omega_{0}$ ($\omega_{\mathrm{c}}$) is
the two-level system (cavity) energy and $a_{c}$ is the bosonic annihilation
operator for the cavity. A standard approach is to rewrite the above
Hamiltonian using spin operators and work in the Dicke basis (Garraway, 2011),
where $S_{+}=\sum_{j=1}^{N}\sigma_{\mathrm{e},\mathrm{g}}^{j}$,
$S_{-}=\left(S_{+}\right)^{\dagger}$ and
$S_{z}=\sum_{j=1}^{N}\frac{1}{2}\left(\sigma_{\mathrm{e},\mathrm{e}}^{j}-\sigma_{\mathrm{g},\mathrm{g}}^{j}\right)$.
The Hamiltonian can then be written as

$H_{\mathrm{TC}}=\omega_{0}S_{z}+\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}+g\left(a_{\mathrm{c}}S_{+}+a_{\mathrm{c}}^{\dagger}S_{-}\right),$
(8)

and the emitter states are $\left|s,m\right\rangle$, where $s$ and $m$ are the
quantum numbers associated to $S^{2}$ and $S_{z}$. The totally symmetric
subspace is then the highest spin subspace, where $s=N/2$.

If we instead apply our approach and second quantize Eq. (missing) 7, we
obtain

$H_{\mathrm{TC}}=\frac{\omega_{0}}{2}\left(b_{\mathrm{e}}^{\dagger}b_{\mathrm{e}}-b_{\mathrm{g}}^{\dagger}b_{\mathrm{g}}\right)+\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}+g\left(a_{\mathrm{c}}b_{\mathrm{e}}^{\dagger}b_{\mathrm{g}}+a_{\mathrm{c}}^{\dagger}b_{\mathrm{g}}^{\dagger}b_{\mathrm{e}}\right).$
(9)

In order to see that both approaches are completely equivalent in the totally
symmetric subspace ($s=N/2$), we examine the matrix elements of $S_{-}$. The
action of $S_{-}$ on a Dicke state is

$S_{-}\left|s,m\right\rangle=\sqrt{s(s+1)-m(m-1)}\left|s,m-1\right\rangle.$
(10)

Here, $m=N_{\mathrm{exc}}-N/2$ is directly related to the number of excited
emitters, $N_{\mathrm{exc}}$.

For the second quantized version, $S_{-}$ maps to $b_{g}^{\dagger}b_{e}$,
which acts on the Fock states $\left|n_{g},n_{e}\right\rangle$, where
$n_{g}=N-N_{\mathrm{exc}}$ and $n_{e}=N_{\mathrm{exc}}$, as

$b_{g}^{\dagger}b_{e}\left|n_{g},n_{e}\right\rangle=\sqrt{(n_{g}+1)n_{e}}\left|n_{g}+1,n_{e}-1\right\rangle.$
(11)

The Dicke state $\left|s=N/2,m=N_{\mathrm{exc}}-N/2\right\rangle$ is equal to
the Fock state $\left|N-N_{\mathrm{exc}},N_{\mathrm{exc}}\right\rangle$, and
comparing Eq. (missing) 10 and Eq. (missing) 11 shows that the matrix elements
are indeed equal. The equality can also easily be checked for
$S_{z}\equiv\frac{1}{2}\left(b_{\mathrm{e}}^{\dagger}b_{\mathrm{e}}-b_{\mathrm{g}}^{\dagger}b_{\mathrm{g}}\right)$.
Therefore, within the totally symmetric subspace, it is completely equivalent
to work with either of the two Hamiltonians, Eq. (missing) 7 or Eq. (missing)
9.

###  III.2 Holstein-Tavis-Cummings model

The field of molecular polaritonics and polaritonic chemistry (Hertzog  _et
al._ , 2019; Herrera and Owrutsky, 2020; Climent  _et al._ , 2021; Garcia-
Vidal  _et al._ , 2021; Fregoni  _et al._ , ; Sánchez-Barquilla  _et al._ , )
studies how to manipulate and use the changes in electronic and vibrational
structure and dynamics of molecules under strong coupling with confined modes
of light. Since molecules are complex systems with significant internal
structure due to rovibrational (nuclear) motion, describing them as two-level
systems is often not a good approximation. At the same time, the influence of
individual collapse operators acting on each molecule can often be neglected.
On the one hand, their individual radiative decay (on scales of nanoseconds)
is often much slower than the dynamics of interest. On the other hand, the
influence of the vibrational modes that is sometimes included through a pure-
dephasing Lindblad term (which has to be replaced by a more careful treatment
under strong light-matter coupling to prevent unphysical effects del Pino  _et
al._ (2015)) can be much better described by treating some vibrational modes
(or superpositions of them corresponding to so-called reaction coordinates)
explicitly, which allows neglecting the other ones at reasonably short
timescales Silva  _et al._ (2020); Zhao  _et al._ (2020). These considerations
apply especially for organic molecules interacting with a plasmonic nanocavity
(Chikkaraddy  _et al._ , 2016; Ojambati  _et al._ , 2019), since their
ultrafast loss is typically the dominant decay channel in the system. Explicit
inclusion of nuclear degrees of freedom also allows to represent many effects
that cannot be understood within a two-level system description (Herrera and
Spano, 2016; Galego  _et al._ , 2015). A workhorse in this field is the so-
called Holstein-Tavis-Cummings model (Herrera and Spano, 2016), in which the
molecule is approximated using the Holstein model, i.e., two displaced
harmonic oscillators for the electronic ground and excited states. Therefore,
when dealing with molecular polaritonics, it is common to face situations
where one needs to solve the dynamics of identical multi-level systems without
any individual collapse operator.

The Holstein-Tavis-Cummings Hamiltonian can be written as

$\displaystyle
H_{\mathrm{HTC}}=\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}+\sum_{i=1}^{N_{\mathrm{mol}}}H_{\mathrm{mol}}^{(i)}+\sum_{i=1}^{N_{\mathrm{mol}}}H_{\mathrm{cav-
mol}}^{(i)},$ (12) $\displaystyle
H_{\mathrm{mol}}^{(i)}=\omega_{\mathrm{e}}\sigma_{i}^{+}\sigma_{i}^{-}+\omega_{\mathrm{v}}c_{i}^{\dagger}c_{i}-\lambda_{\mathrm{v}}\sigma_{i}^{+}\sigma_{i}^{-}\left(c_{i}^{\dagger}+c_{i}\right)$
(13) $\displaystyle H_{\mathrm{cav-
mol}}^{(i)}=g\left(\sigma_{i}^{+}a_{\mathrm{c}}+a_{\mathrm{c}}^{\dagger}\sigma_{i}^{-}\right),$
(14)

where $\sigma_{i}^{+}$ ($\sigma_{i}^{-}$) is the raising (lowering) operator
for the electronic state in molecule $i$ with excitation energy
$\omega_{\mathrm{e}}$, whereas $c_{i}$ is the annihilation operator for the
vibrational mode in molecule $i$, with frequency $\omega_{\mathrm{v}}$ and
exciton–phonon coupling strength $\lambda_{\mathrm{v}}$. The cavity is
described through the photon annihilation (creation) operators
$a_{\mathrm{c}}$ ($a_{\mathrm{c}}^{\dagger}$), with photon energy
$\omega_{\mathrm{c}}$. In addition to the coherent dynamics described by the
Hamiltonian, the cavity mode decays with rate $\gamma_{\mathrm{c}}$, described
by a standard Lindblad decay operator
$C=\sqrt{\gamma_{\mathrm{c}}}a_{\mathrm{c}}$.

The Holstein-Tavis-Cummings Hamiltonian can be rewritten in terms of the
eigenstates of the single-molecule Hamiltonian,

$H_{\mathrm{mol}}^{(i)}=\sum_{s=\mathrm{g},\mathrm{e}}\sum_{\nu}\omega_{\mathrm{s},\nu}|\mathrm{s},\nu\rangle_{i}\langle\mathrm{s},\nu|_{i},$
(15)

which are labeled as $|\mathrm{g},\nu\rangle_{i}$ and
$|\mathrm{e},\nu\rangle_{i}$ for vibrational sublevel $\nu$ in the electronic
ground and excited state, respectively. Their corresponding energies are
$\omega_{\mathrm{g},\nu}=\omega_{\mathrm{v}}\nu$ and
$\omega_{\mathrm{e},\nu}=\omega_{\mathrm{e}}+\omega_{\mathrm{v}}\nu-\lambda_{\mathrm{v}}^{2}/\omega_{\mathrm{v}}$.
In this basis, the light-matter interaction operator is given by

$H_{\mathrm{cav-
mol}}^{(i)}=g\sum_{\nu\nu^{\prime}}\left(\text{$a_{\mathrm{c}}F_{\nu\nu^{\prime}}\left|\mathrm{e},\nu\right\rangle_{i}\left\langle\mathrm{g},\nu^{\prime}\right|_{i}$}+\mathrm{H.c.}\right),$
(16)

where
$F_{\nu\nu^{\prime}}=\left\langle\mathrm{e},\nu\left|\sigma_{i}^{+}\right|\mathrm{g},\nu^{\prime}\right\rangle$
is a vibrational overlap integral or Franck-Condon factor and can be
analytically obtained.

Figure 2: (a) Normalized cavity population $\left\langle
a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}\right\rangle N_{\mathrm{mol}}^{-1}$ in
the Holstein-Tavis-Cummings model, for $N_{\mathrm{mol}}=1,2,3,4,5$ calculated
both within the second-quantization approach (full lines) and without
exploiting permutational symmetry (dashed lines). See main text for
parameters. (b) The same for the electronic excited state population,
$\sum_{i}^{N_{\mathrm{mol}}}\left\langle\sigma_{i}^{+}\sigma_{i}^{-}\right\rangle
N_{\mathrm{mol}}^{-1}$.

Since both the Hamiltonian and collapse operators are permutationally
invariant, we can map Eq. (missing) 12 to its second quantized form as long as
the initial state is fully symmetric. This gives

$H_{\mathrm{HTC}}=\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}+\sum_{s=\mathrm{g},\mathrm{e}}\sum_{\nu}\omega_{\mathrm{s},\nu}b_{\mathrm{s},\nu}^{\dagger}b_{\mathrm{s},\nu}\\\
+g\sum_{\nu\nu^{\prime}}\left(\text{$a_{\mathrm{c}}F_{\nu\nu^{\prime}}b_{\mathrm{e},\nu}^{\dagger}b_{\mathrm{g},\nu^{\prime}}$}+\mathrm{H.c.}\right).$
(17)

In the following, we choose parameter values typical for organic molecules
such as anthracene coupled to nanoplasmonic cavities Silva  _et al._ (2020),
with $\omega_{\mathrm{e}}=3.5\,$eV, $\omega_{\mathrm{v}}=0.182\,$eV,
$\lambda_{\mathrm{v}}=0.096\,$eV, $\gamma_{\mathrm{c}}=0.2\,$eV,
$g=0.035\,$eV. We set the cavity photon energy to be on resonance with the
emission peak of the molecule,
$\omega_{\mathrm{c}}=\omega_{\mathrm{e}}-2\lambda_{\mathrm{v}}^{2}/\omega_{\mathrm{v}}$.
The initial state is chosen to be the fully inverted state, i.e., the state
where all molecules are instantaneously excited to the electronic excited
state by a vertical Franck-Condon transition.

In Figure 2(a,b), we show the numerical results for the dynamics for the
cavity and excited state population, respectively. In the basis truncation for
the single-molecule Hilbert space, we include the 6 lowest vibrational states
for the ground state and the 4 lowest vibrational states for the electronic
excited state, which gives converged results. Within the second quantization
approach, we show results up to $N_{\mathrm{mol}}=5$, while for the results
without resorting to any permutational symmetry, we show results up to
$N_{\mathrm{mol}}=3$.

The exciton population in Figure 2(b) displays a clear enhancement of the
spontaneous emission due to Dicke superradiance Dicke (1954) as the number of
emitters is increased. Furthermore, a modulation of the decay rate with a
period of about $22\,$fs can be observed. This modulation is more clearly
visible in the cavity population, see Figure 2(a), and is a signature of the
vibrational motion Silva  _et al._ (2020), which has a period of
$T_{\mathrm{v}}=2\pi/\omega_{\mathrm{v}}=22.7\,$fs.

As expected, both approaches are completely equivalent. However, while the
Hilbert space for the brute-force approach reaches size
$N_{\mathrm{Hilb}}=4000$ for three molecules, it only has size
$N_{\mathrm{Hilb}}=220$ within the second quantization approach exploiting the
permutational symmetry. For five molecules, this advantage improves to
$N_{\mathrm{Hilb}}=12012$ versus $N_{\mathrm{Hilb}}=600000$. Here, it should
be noted that the size of the density matrix that is propagated in the
Lindblad master equation is $N_{\mathrm{Hilb}}\times N_{\mathrm{Hilb}}$, while
the Liouvillian superoperator describing this evolution can be formally
treated as a $N_{\mathrm{Hilb}}^{2}\times N_{\mathrm{Hilb}}^{2}$ matrix.

###  III.3 Three-level systems

Figure 3:  (a) Dynamics for $N=5$ three-level emitters calculated within the
second quantized approach (full lines) and without using the permutational
symmetry (dashed lines). See main text for parameters. (b) The same for $N=17$
emitters, which is only possible with reasonable effort when using the second
quantized approach.

Figure 4:  (a) Expectation value of the dipole-dipole interaction term,
$\langle\mu_{i}\mu_{j}\rangle$ for $i\neq j$ for the simulation shown in
Figure 3(a) (with $N=5$), within the second quantized approach (full lines)
and without using the permutational symmetry (dashed lines). (b) The same for
the simulation in Figure 3(b) with $N=17$ emitters, only using the second
quantized approach.

To give another numerical example to illustrate this mapping, we formulate a
simple model Hamiltonian of $N$ $d$-level systems, where the levels of each
emitter are equally separated in energy by $\omega_{\mathrm{e}}$, coupled to a
cavity mode with frequency
$\omega_{\mathrm{c}}=\omega_{\mathrm{e}}=1\,\mathrm{eV}$. The transition
operator of each emitter is defined as $\mu_{i}=\mu_{i}^{+}+\mu_{i}^{-}$,
where $\mu_{i}^{+}=\sum_{\nu=1}^{d-1}\sigma_{\nu,\nu+1}^{i}$ and
$\mu_{i}^{-}=\left(\mu_{i}^{+}\right)^{\dagger}$ and each emitter is coupled
to the cavity mode by a coupling strength of
$g=\frac{0.15}{\sqrt{N}}\,\mathrm{eV}$ in the rotating wave approximation. We
also include an all-to-all dipole-dipole interaction term,
$H_{\mathrm{d-d}}=D\sum_{k,j\neq k}^{N}\mu_{i}\mu_{j}$, where
$D=0.1\,\mathrm{eV}$. The Hamiltonian is then given by

$H=\sum_{i=1}^{N}\sum_{\nu=1}^{d}\omega_{\mathrm{\nu}}\sigma_{\nu,\nu}^{i}+\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}+H_{\mathrm{d-d}}\\\
+g\sum_{i=1}^{N}\left(\mu_{i}^{-}a_{\mathrm{c}}^{\dagger}+\mathrm{H.c.}\right),$
(18)

where $\omega_{\mathrm{\nu}}=\nu\omega_{\mathrm{e}}$. The system is under the
action of the incoherent decay of the cavity,
$C_{\mathrm{cav}}=\sqrt{\gamma_{\mathrm{c}}}a_{c}$ where
$\gamma_{\mathrm{c}}=0.15\,$eV, as well as collective spontaneous emission,
$C_{\nu}=\sqrt{\Gamma_{\downarrow}}\sum_{i=1}^{N}\sigma_{\nu,\nu+1}^{i}$,
where $\nu$ runs from $1$ to $d-1$ and
$\Gamma_{\downarrow}=0.05\,\mathrm{eV}$. The initial state is chosen to be the
fully inverted state,
$|\psi_{0}\rangle=\prod_{j=1}^{N}|d\rangle_{j}\otimes|\chi\rangle_{\mathrm{cav}}$,
where all emitters are in the most excited state and the cavity is in the
vacuum state.

The above Hamiltonian and collapse operators are clearly invariant under any
permutation of the emitters. Also, the initial state belongs to the totally
symmetric subspace. Therefore, we can again proceed with the mapping by second
quantizing all relevant operators. In particular, the Hamiltonian can be
written as

$H=\sum_{\nu=1}^{d}\omega_{\mathrm{\nu}}b_{\nu}^{\dagger}b_{\nu}+\omega_{\mathrm{c}}a_{\mathrm{c}}^{\dagger}a_{\mathrm{c}}+H_{\mathrm{d-d}}\\\
+g\sum_{i=1}^{N}\sum_{\nu=1}^{d}\left(b_{\nu}^{\dagger}b_{\nu+1}a_{\mathrm{c}}^{\dagger}+\mathrm{h.c.}\right),$
(19)

where

$H_{\mathrm{d-d}}=D\sum_{\nu=1}^{d-1}\sum_{\mu=1}^{d-1}\left(b_{\nu}^{\dagger}b_{\mu}^{\dagger}b_{\nu+1}b_{\mu+1}+b_{\nu+1}^{\dagger}b_{\mu+1}^{\dagger}b_{\nu}b_{\mu}\right.\\\
\left.+b_{\nu+1}^{\dagger}b_{\mu}^{\dagger}b_{\nu}b_{\mu+1}+b_{\nu}^{\dagger}b_{\mu+1}^{\dagger}b_{\nu+1}b_{\mu}\right).$
(20)

Note that this operator is expressed using normal ordering. The $d-1$
collective spontaneous emission collapse operators can be rewritten as
$C_{\nu}=\sqrt{\Gamma_{\downarrow}}b_{\nu}^{\dagger}b_{\nu+1}$. Finally, the
initial state is just
$|\psi_{0}\rangle=\frac{1}{\sqrt{N!}}(b_{d}^{\dagger})^{N}|\mathrm{vac}\rangle_{\mathrm{em}}|\mathrm{vac}\rangle_{\mathrm{cav}}$.

In Figure 3 and Figure 4, we show the results of the dynamics for the case of
3-level systems, i.e. $d=3$. The total population of the different levels,
$\left\langle\sum_{i=1}^{N}\sigma_{\nu,\nu}^{i}\right\rangle$, can be mapped
in the second quantized approach to $\left\langle
b_{\nu}^{\dagger}b_{\nu}\right\rangle$. Figure 3(a) shows the results for
$N=5$. In Figure 3(b), we show the results for $N=17$ emitters. In Figure 4,
we show the time-dependent expectation value of the dipole-dipole interaction
term, $\langle\mu_{i}\mu_{j}\rangle$ for $i\neq j$.

Again, as expected, the second quantized approach is completely equivalent to
the direct solution. In this case, the brute-force approach is numerically
intractable, as the number of entries in the density matrix is
$d^{2N}N_{\mathrm{c}}^{2}$, where $N_{\mathrm{c}}$ is the dimension of the
cavity Hilbert space. When working only with the totally symmetric subspace,
the number of entries in the density matrix is reduced to

$\left(\frac{\left(N+d-1\right)!}{N!\left(d-1\right)!}\right)^{2}N_{\mathrm{c}}^{2},$
(21)

greatly reducing the size of the dynamical object. When comparing with the
approach that uses the symmetrized Liouville space, where the number of
entries in the density matrix is
$\frac{\left(N+d^{2}-1\right)!}{N!\left(d^{2}-1\right)!}N_{\mathrm{c}}^{2}$
(Gegg and Richter, 2016; Gegg, 2017), we also get a substantial reduction. As
an example, for $d=3$ and $N=17$, we have a reduction of the number of entries
in the density matrix by a factor of $37$.

Since we start in the fully inverted state, the dipole-dipole interaction,
$H_{\mathrm{d-d}}$, starts to transfer population from the highest excited
emitter state to the intermediate excited state and to a smaller extent to the
emitter ground state. After this first moment, the cavity starts to become
populated and due to its decay, drives the system to its overall ground state,
see Figure 3. It is important to notice that due to the dipole-dipole
interaction, $H_{\mathrm{d-d}}$, the ground state of the system is not the
state where all emitters are in their bare ground state. Consequently, there
is a residual population of the intermediate excited state for long times, see
Figure 3. This is also the reason why the dipole-dipole interaction goes to
negative values for long times, see Figure 4.

###  III.4 $N$-excitation subspace

When working with the dynamics of emitters coupled to cavity modes, there are
situations in which we are not interested in working with the full excitation
subspace. Indeed, in many common cases, restricting to the first or second
excitation subspace is enough (Galego  _et al._ , 2015; Feist  _et al._ ,
2018). Implementing such a restriction within the current approach is rather
simple, as one only needs to define an operator that determines the number of
excitations in terms of creation and annihilation operators of the emitter
levels. As an example, if one is working with the Holstein-Tavis-Cummings
model, where each emitter is described as having two electronic states, ground
and excited, with one vibrational mode, one could define a subspace where
restrictions are imposed on either the electronic or nuclear excitations, or
both.

As a concrete example, we discuss vibrational strong coupling for the case
where a single (approximately harmonic) vibrational mode per molecule is in
resonance with a cavity mode. The simplest Hamiltonian to model collective
vibrational strong coupling is (del Pino  _et al._ , 2015)

$H=\omega_{\mathrm{c}}a_{c}^{\dagger}a_{c}+\sum_{i=1}^{N_{\mathrm{mol}}}\omega_{\mathrm{v}}c_{i}^{\dagger}c_{i}+\sum_{i=1}^{N_{\mathrm{mol}}}g\left(a_{c}^{\dagger}c_{i}+\mathrm{H.c.}\right),$
(22)

where $a_{\mathrm{c}}$ is the annihilation operator for the cavity mode with
frequency $\omega_{\mathrm{c}}$, and $c_{i}$ is the annihilation operator of
the optically active vibrational mode of molecule $i$, characterized by its
frequency $\omega_{\mathrm{v}}$. $N_{\mathrm{mol}}$ is the number of
molecules, and the cavity-phonon interaction is given by $g$. Rewriting the
vibrational operators using the eigenstates of the harmonic oscillator,
$c_{i}^{\dagger}=\sum_{n=0}^{\infty}\sqrt{n+1}\left|n+1\right\rangle_{i}\left\langle
n\right|_{i}$, Eq. (missing) 22 can be written as

$H=\omega_{\mathrm{c}}a_{c}^{\dagger}a_{c}+\sum_{i=1}^{N_{\mathrm{mol}}}\sum_{n=0}^{\infty}n\omega_{\mathrm{v}}\left|n\right\rangle_{i}\left\langle
n\right|_{i}\\\
+\sum_{i=1}^{N_{\mathrm{mol}}}\sum_{n=0}^{\infty}g\left(a_{c}\sqrt{n+1}\left|n+1\right\rangle_{i}\left\langle
n\right|_{i}+\mathrm{h.c.}\right).$ (23)

This Hamiltonian is permutationally invariant under the exchange of any two
molecules. If the initial state is in the totally symmetric subspace, we can
map the Hamiltonian to

$H=\omega_{\mathrm{c}}a_{c}^{\dagger}a_{c}+\sum_{n=0}^{\infty}n\omega_{\mathrm{v}}b_{n}^{\dagger}b_{n}\\\
+\sum_{n=0}^{\infty}g\left(a_{c}\sqrt{n+1}b_{n+1}^{\dagger}b_{n}+\mathrm{h.c.}\right),$
(24)

where $b_{n}$ is the bosonic annihilation operator for the state
$\left|n\right\rangle$ and the states of interest are restricted to the
subspace where
$\left\langle\sum_{n=0}^{\infty}b_{n}^{\dagger}b_{n}\right\rangle=N_{\mathrm{mol}}$.
For regimes in which $\omega_{\mathrm{v}}\approx\omega_{\mathrm{c}}$, it is
reasonable to work in the $N_{\mathrm{exc}}$-excitation subspace (Campos-
Gonzalez-Angulo  _et al._ , 2021). In this formalism, this additional
restriction can be simply formulated as
$\left\langle\sum_{n=0}^{\infty}nb_{n}^{\dagger}b_{n}+a_{c}^{\dagger}a_{c}\right\rangle=N_{\mathrm{exc}}$.

##  IV Conclusion

To conclude, we have proposed a scheme to fully exploit the permutational
symmetry of identical, but arbitrary emitters when only collective dissipation
operators are considered. This scheme relies on the fact that the totally
symmetric subspace is equivalent to a bosonic many-body state. After mapping
all relevant operators to a second quantized picture using a simple procedure,
the explicit construction of the totally symmetric subspace from direct state
products is not required anymore. This approach thus provides a
straightforward and easily implemented way to treat such systems while fully
exploiting their permutational symmetry to significantly reduce the size of
the Hilbert space. We discuss several examples, such as the Tavis-Cummings
model, the Holstein-Tavis-Cummings model and a model Hamiltonian where two-
body operators are taken into account, and explicitly demonstrate the
equivalence of the second quantized approach to direct solution.

We expect that this work will be helpful for simulations that can fully
exploit the permutational symmetry of emitters in totally symmetric cases in a
very simple way. This can be especially useful for situations where each
emitter must be considered as having an internal structure that goes beyond
the two-level approximation, such as necessary in the field of molecular
polaritonics. For such systems, the current approach can provide a significant
reduction of the numerical complexity for very little effort.

## Acknowledgments

This work has been funded by the European Research Council through grant
ERC-2016-STG-714870 and by the Spanish Ministry for Science, Innovation, and
Universities-AEI through grants RTI2018-099737-B-I00, PCI2018-093145 (through
the QuantERA program of the European Commission), and CEX2018-000805-M
(through the María de Maeztu program for Units of Excellence in R&D). R. E. F.
S. also acknowledges support from the fellowship LCF/BQ/PR21/11840008 from “La
Caixa” Foundation (ID 100010434) and from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No. 847648.

## References

  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Richter  _et al._ (2015) Marten Richter, Michael Gegg, T Sverre Theuerholz, and Andreas Knorr, “Numerically exact solution of the many emitter–cavity laser problem: Application to the fully quantized spaser emission,” Physical Review B 91, 035306 (2015).
  * Wang and Hioe (1973) Yo K Wang and FT Hioe, “Phase transition in the dicke model of superradiance,” Physical Review A 7, 831 (1973).
  * Walls  _et al._ (1978) DF Walls, PD Drummond, SS Hassan, and HJ Carmichael, “Non-equilibrium phase transitions in cooperative atomic systems,” Progress of Theoretical Physics Supplement 64, 307–320 (1978).
  * Gegg  _et al._ (2018) Michael Gegg, Alexander Carmele, Andreas Knorr, and Marten Richter, “Superradiant to subradiant phase transition in the open system dicke model: Dark state cascades,” New Journal of Physics 20, 013006 (2018).
  * Dicke (1954) Robert H Dicke, “Coherence in spontaneous radiation processes,” Physical Review 93, 99 (1954).
  * Garraway (2011) Barry M Garraway, “The dicke model in quantum optics: Dicke model revisited,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369, 1137–1155 (2011).
  * Herrera and Spano (2016) Felipe Herrera and Frank C Spano, “Cavity-controlled chemistry in molecular ensembles,” Physical Review Letters 116, 238301 (2016).
  * Feist  _et al._ (2018) Johannes Feist, Javier Galego, and Francisco J Garcia-Vidal, “Polaritonic chemistry with organic molecules,” ACS Photonics 5, 205–216 (2018).
  * Putz  _et al._ (2017) Stefan Putz, Andreas Angerer, Dmitry O Krimer, Ralph Glattauer, William J Munro, Stefan Rotter, Jörg Schmiedmayer, and Johannes Majer, “Spectral hole burning and its application in microwave photonics,” Nature Photonics 11, 36–39 (2017).
  * Gegg and Richter (2016) Michael Gegg and Marten Richter, “Efficient and exact numerical approach for many multi-level systems in open system cqed,” New Journal of Physics 18, 043037 (2016).
  * Gegg and Richter (2017) Michael Gegg and Marten Richter, “Psiquasp–a library for efficient computation of symmetric open quantum systems,” Scientific Reports 7, 1–14 (2017).
  * Gegg (2017) Michael Gegg,  _Identical emitters, collective effects and dissipation in quantum optics_ , Ph.D. thesis, Technische Universität Berlin (2017).
  * Shammah  _et al._ (2018) Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Simone De Liberato, and Franco Nori, “Open quantum systems with local and collective incoherent processes: Efficient numerical simulations using permutational invariance,” Physical Review A 98, 063815 (2018).
  * Kirton and Keeling (2017) Peter Kirton and Jonathan Keeling, “Suppressing and restoring the dicke superradiance transition by dephasing and decay,” Physical Review Letters 118, 123602 (2017).
  * Kirton and Keeling (2018) Peter Kirton and Jonathan Keeling, “Superradiant and lasing states in driven-dissipative dicke models,” New Journal of Physics 20, 015009 (2018).
  * Bruus and Flensberg (2004) Henrik Bruus and Karsten Flensberg,  _Many-Body Quantum Theory in Condensed Matter Physics: An Introduction_ (Oxford University Press, 2004).
  * Johansson  _et al._ (2012) J Robert Johansson, D Nation Paul, and Franco Nori, “Qutip 2: An open-source python framework for the dynamics of open quantum systems,” Computer Physics Communications 184, 1234 (2012).
  * Cohen-Tannoudji  _et al._ (1987) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms: Introduction to Quantum Electrodynamics_ (Wiley-Interscience, 1987).
  * Rabi (1937) Isidor Isaac Rabi, “Space quantization in a gyrating magnetic field,” Physical Review 51, 652 (1937).
  * Jaynes and Cummings (1963) Edwin T Jaynes and Frederick W Cummings, “Comparison of quantum and semiclassical radiation theories with application to the beam maser,” Proceedings of the IEEE 51, 89–109 (1963).
  * Tavis and Cummings (1968) Michael Tavis and Frederick W Cummings, “Exact solution for an n-molecule—radiation-field hamiltonian,” Physical Review 170, 379 (1968).
  * Hertzog  _et al._ (2019) Manuel Hertzog, Mao Wang, Jürgen Mony, and Karl Börjesson, “Strong Light–Matter Interactions: A New Direction within Chemistry,” Chem. Soc. Rev. 48, 937 (2019).
  * Herrera and Owrutsky (2020) Felipe Herrera and Jeffrey Owrutsky, “Molecular Polaritons for Controlling Chemistry with Quantum Optics,” J. Chem. Phys. 152, 100902 (2020).
  * Climent  _et al._ (2021) C Climent, FJ Garcia-Vidal, and J Feist, “Cavity-modified chemistry: Towards vacuum-field catalysis,” in  _Effects of Electric Fields on Structure and Reactivity_ (2021) pp. 343–393.
  * Garcia-Vidal  _et al._ (2021) Francisco J. Garcia-Vidal, Cristiano Ciuti, and Thomas W. Ebbesen, “Manipulating Matter by Strong Coupling to Vacuum Fields,” Science 373, eabd0336 (2021).
  * (27) Jacopo Fregoni, Francisco J. García-Vidal, and Johannes Feist, “Theoretical Challenges in Polaritonic Chemistry,” arXiv:2111.08394 .
  * (28) Mónica Sánchez-Barquilla, Antonio I. Fernández-Domínguez, Johannes Feist, and Francisco J. García-Vidal, “A Theoretical Perspective on Molecular Polaritonics,” arXiv:2201.02827 .
  * del Pino  _et al._ (2015) Javier del Pino, Johannes Feist, and Francisco J Garcia-Vidal, “Quantum theory of collective strong coupling of molecular vibrations with a microcavity mode,” New Journal of Physics 17, 053040 (2015).
  * Silva  _et al._ (2020) R E F Silva, Javier Del Pino, Francisco J García-Vidal, and Johannes Feist, “Polaritonic molecular clock for all-optical ultrafast imaging of wavepacket dynamics without probe pulses,” Nature Communications 11, 1–8 (2020).
  * Zhao  _et al._ (2020) Dongxing Zhao, R E F Silva, Clàudia Climent, Johannes Feist, Antonio I Fernández-Domínguez, and Francisco J García-Vidal, “Impact of vibrational modes in the plasmonic purcell effect of organic molecules,” ACS Photonics 7, 3369–3375 (2020).
  * Chikkaraddy  _et al._ (2016) Rohit Chikkaraddy, Bart De Nijs, Felix Benz, Steven J Barrow, Oren A Scherman, Edina Rosta, Angela Demetriadou, Peter Fox, Ortwin Hess, and Jeremy J Baumberg, “Single-molecule strong coupling at room temperature in plasmonic nanocavities,” Nature 535, 127–130 (2016).
  * Ojambati  _et al._ (2019) Oluwafemi S Ojambati, Rohit Chikkaraddy, William D Deacon, Matthew Horton, Dean Kos, Vladimir A Turek, Ulrich F Keyser, and Jeremy J Baumberg, “Quantum electrodynamics at room temperature coupling a single vibrating molecule with a plasmonic nanocavity,” Nature Communications 10, 1–7 (2019).
  * Galego  _et al._ (2015) Javier Galego, Francisco J Garcia-Vidal, and Johannes Feist, “Cavity-induced modifications of molecular structure in the strong-coupling regime,” Physical Review X 5, 041022 (2015).
  * Campos-Gonzalez-Angulo  _et al._ (2021) Jorge Arturo Campos-Gonzalez-Angulo, Raphael Florentino Ribeiro, and Joel Yuen Zhou, “Generalization of the tavis-cummings model for multi-level anharmonic systems.” New Journal of Physics 23, 063081 (2021).