# Two electrons in harmonic confinement coupled to light in a cavity

Chenhang Huang  Department of Physics and Astronomy, Vanderbilt University,
Nashville, Tennessee, 37235, USA  Alexander Ahrens  Department of Physics and
Astronomy, Vanderbilt University, Nashville, Tennessee, 37235, USA  Matthew
Beutel  Department of Physics and Astronomy, Vanderbilt University, Nashville,
Tennessee, 37235, USA  Kálmán Varga<EMAIL_ADDRESS>Department of
Physics and Astronomy, Vanderbilt University, Nashville, Tennessee, 37235, USA

###### Abstract

The energy and wave function of a harmonically confined two-electron system
coupled to light is calculated by separating the wave functions of the
relative and center of mass (CM) motions. The relative motion wave function
has a known quasi-analytical solution. The light only couples to the CM
variable and the coupled equation can be solved with diagonalization without
approximations. The approach works for any coupling strength. Examples of wave
functions of light-matter hybrid states are presented.

##  I Introduction

Analytically or numerically easily solvable systems (e.g. by ”exact
diagonalization”) have always been important test grounds for models and
approximations. Recently, there is an intense interest in strongly coupled
light-matter systems Buchholz et al. (2019); Schäfer et al. (2019);
Ruggenthaler et al. (2018); Flick et al. (2015, 2017); Rokaj et al. (2018);
Rivera et al. (2019); Flick and Narang (2018); Hoffmann et al. (2020); Tokatly
(2018); Galego et al. (2017); Mandal et al. (2020a); Cederbaum and Kuleff
(2021); Szidarovszky et al. (2018); Ashida et al. (2021). In these systems,
the light-matter coupling cannot be treated perturbatively. The electronic
excitations and the photons are superimposed, forming hybrid light-matter
excitations. In this regime, there are only a few analytical approaches
available to test and develop efficient numerical methods. Reviews of the
recent theoretical and experimental development can be found in Refs.Rivera
and Kaminer (2020); Le Boité (2020); Garcia-Vidal et al. (2021).

In this paper, we consider a two-electron system interacting via the Coulomb
interaction, confined by a harmonic oscillator interaction coupled to light in
a cavity. The system is described on the level of the Pauli-Fierz (PF)
nonrelativistic QED Hamiltonian. The two-electron system in harmonic
oscillator confinement is a quasi-exactly solvable (QES) problem. The wave
function can be written as a product of the wave functions of the relative and
CM motion. The relative motion wave function can be expanded into infinite
series. For certain oscillator parameters, this infinite series can be reduced
to a recursion Taut (1993). The wave function of the CM motion is a simple
harmonic oscillator eigenfunction. We will show that the photons only couple
to the CM coordinate and the coupled CM photon system can be solved by exact
diagonalization.

The two-particle systems have long been investigated due to their analytic and
quasi-exact solvability, which provides straightforward intuition for the
physical system under scrutiny as well as an excellent benchmark test for
numerical computations. Examples of QES quantum systems are the two-
dimensional (2D) harmoniumsTaut (1993, 1994) and the hydrogen-like atoms in
homogeneous magnetic fieldsTaut (1995). These QES problems have been
generalized to relativistic cases as well Karwowski and Pestka (2007);
Villalba and Pino (1998). For harmonium systems, the separability condition
guarantees the quasi-exact solvability for the Schrödinger equationKarwowski
(2008), and linearly coupled oscillators have been studied under this
conditionTurbiner (1988). For hydrogen-like models, solutions have been found
for particular forms of the inhomogeneous magnetic fields Liu and Hao (2015);
Downing and Portnoi (2016). Examples of other known QES models include the
planar Dirac electron in hydrogen-like atomsHo and Khalilov (2000); Chiang and
Ho (2002), one-body problems in power-law central potentialsAgboola (2012);
Karwowski and Witek (2016), relativistic 2D pion in constant magnetic
fieldsAkhmedov and Guseinova (2009), and 1D and 3D regularized Calogero
modelsPont et al. (2018); Downing (2017). QES models with different forms of
confinements, e.g. two electrons in oneLoos and Gill (2012) or twoGuo et al.
(2012) 1D rings, two electrons on the surface of the n-sphere (spherium)Loos
and Gill (2009, 2010), have also been studied.

The exact or even the numerical solution for light-matter coupled systems is
very difficult even on the level of a minimal coupling Hamiltonian in the
long-wavelength limit Jestädt et al. (2019), because the photons substantially
increase the number degrees of freedom of the system. Theoretical approaches
have been developed to tame the light-matter coupled systems using
approximations and transformations Frisk Kockum et al. (2019); Schäfer et al.
(2018); Flick et al. (2018); Ashida et al. (2021); Sidler et al. (2020, 2020);
Jestädt et al. (2019); Lacombe et al. (2019); Ashida et al. (2021); Andolina
et al. (2019); Schuler et al. (2020); Settineri et al. (2021); Sidler et al.
(2020); Rokaj et al. (2021). In Refs. Schäfer et al. (2018); Flick et al.
(2018), an electron in a 2D potential coupled to a single photon mode is used
as a numerical benchmark test. The spatial part of the wave function is
represented on a real space grid and coupled to the Fock space of the photons.
The Hamiltonian of the system can be diagonalized in this representation and
the light coupled wave function can be studied. In Ref. Sidler et al. (2020),
the spatial wave function of the He, HD+, and H${}_{2}^{+}$ three-particle
system is represented using a 3D product of pseudospectral basis functions,
and a few Fock spaces states of a single photon mode are coupled to the
spatial part. The energy and wave function is calculated by exact
diagonalization of the PF Hamiltonian and the Jaynes–Cummings limit for
electronic and ro-vibrational transitions are studied. One-dimensional model
systems of atoms and molecules Buchholz et al. (2019); Lacombe et al. (2019)
often using the Shin-Metiu potential Shin and Metiu (1995) are also useful to
describe potential energy surfaces in cavities and test numerical approaches.

The free electron gas also allows analytical treatment Rokaj et al. (2021). In
Ref. Rokaj et al. (2021), the free electron gas in cavity is analytically
solved in the long-wavelength limit for an arbitrary number of non-interacting
electrons. It is found that the electron-photon ground state is a Fermi liquid
containing virtual photons.

Approaches to reformulating the problem have also been proposed. In Ref.
Ashida et al. (2021), the light and matter degrees of freedom are decoupled
using a unitary transformation. In the transformed frame, both the light and
the matter Hilbert spaces can be truncated systematically to facilitate an
efficient solution. In Ref. Rivera et al. (2019) a variational formulation is
developed and the semianalytical formula is derived for the ground and excited
state energies.

##  II Formalism

We consider two particles with positions $\mathbf{r}_{1}$, $\mathbf{r}_{2}$
and charges $q_{1}$, $q_{2}$. Later we show that an analytical approach only
works for $q_{1}=q_{2}$, but it is useful to consider the general case to show
the origin of the coupling to the center of mass. The Hamiltonian of the
system is

$H=H_{e}+H_{ph}=H_{e}+H_{p}+H_{ep}+H_{d}.$ (1)

$H_{e}$ is the electronic Hamiltonian, $H_{ph}$ describes the electron-photon
interaction, which is a sum of three terms, the photon Hamiltonian $H_{p}$,
the electron-photon coupling $H_{ep}$, and the dipole self-interaction
$H_{d}$. The electron-photon interaction can be described by using the PF
nonrelativistic QED Hamiltonian. The PF Hamiltonian can be rigorously derived
Ruggenthaler et al. (2018); Rokaj et al. (2018); Mandal et al. (2020a, b);
Tokatly (2018) by applying the Power-Zienau-Woolley gauge transformation Power
et al. (1959), with a unitary phase transformation on the minimal coupling
($p\cdot A$) Hamiltonian in the Coulomb gauge,

$H_{ph}={1\over
2}\sum_{\alpha=1}^{N_{p}}\left[p_{\alpha}^{2}+\omega_{\alpha}^{2}\left(q_{\alpha}-{\boldsymbol{\lambda_{\alpha}}\over\omega_{\alpha}}\cdot\mathbf{D}\right)^{2}\right],$
(2)

where $\mathbf{D}=\sum_{i=1}^{N}q_{i}\mathbf{r}_{i}$ is the dipole operator.
The photon fields are described by quantized oscillators.
$q_{\alpha}={1\over\sqrt{2\omega_{\alpha}}}(\hat{a}^{+}_{\alpha}+\hat{a}_{\alpha})$
is the displacement field and $p_{\alpha}=-i\sqrt{{\omega_{\alpha}}\over
2}(\hat{a}_{\alpha}-\hat{a}^{+}_{\alpha})$ is the conjugate momentum. This
Hamiltonian describes $N_{p}$ photon modes with frequency $\omega_{\alpha}$
and coupling $\boldsymbol{\lambda}_{\alpha}$. The coupling term is usually
written as Ruggenthaler et al. (2014)

$\boldsymbol{\lambda}_{\alpha}=\sqrt{4\pi}\,S_{\alpha}(\mathbf{r})\mathbf{e}_{\alpha},$
(3)

where $S_{\alpha}(\mathbf{r})$ is the mode function at position $\mathbf{r}$
and $\mathbf{e}_{\alpha}$ is the transversal polarization vector of the photon
modes.

The three components of the electron-photon interaction are as follows: The
photonic part is

$H_{p}=\sum_{\alpha=1}^{N_{p}}\left(\frac{1}{2}p_{\alpha}^{2}+\frac{\omega_{\alpha}^{2}}{2}q_{\alpha}^{2}\right)=\sum_{\alpha=1}^{N_{p}}\omega_{\alpha}\left(\hat{a}_{\alpha}^{+}\hat{a}_{\alpha}+\frac{1}{2}\right).$
(4)

By using the creation and annihilation operators, the photon states
$|n_{\alpha}\rangle$ can be generated by multiple applications of the creation
operators on the vacuum state
$|n_{\alpha}\rangle=(\hat{a}_{\alpha}^{+})^{n}|0\rangle$. All other photon
operations can be done by using $\hat{a}_{\alpha}$ and $\hat{a}_{\alpha}^{+}$.
The interaction term is

$H_{ep}=-\sum_{\alpha=1}^{N_{p}}\omega_{\alpha}q_{\alpha}\boldsymbol{\lambda_{\alpha}}\cdot\mathbf{D}=-\sum_{\alpha=1}^{N_{p}}\sqrt{\omega_{\alpha}\over
2}(\hat{a}_{\alpha}+\hat{a}_{\alpha}^{+})\boldsymbol{\lambda_{\alpha}}\cdot\mathbf{D}.$
(5)

Only photon states $|n_{\alpha}\rangle$, $|n_{\alpha}\pm 1\rangle$ are
connected by $\hat{a}_{\alpha}$ and $\hat{a}_{\alpha}^{+}$. The matrix
elements of the dipole operator $\mathbf{D}$ are only nonzero between spatial
basis functions with angular momentum $l$ and $l\pm 1$ in 3D or $m$ and $m\pm
1$ in2D. The strength of the electron-photon interaction can be characterized
by the effective coupling parameter

$g_{\alpha}=\left|\boldsymbol{\lambda_{\alpha}}\right|\sqrt{\frac{\omega_{\alpha}}{2}}.$
(6)

The dipole self-interaction is

$H_{d}={1\over
2}\sum_{\alpha=1}^{N_{p}}\left(\boldsymbol{\lambda_{\alpha}}\cdot\mathbf{D}\right)^{2},$
(7)

which describes the effects of the polarization of the electrons back on the
photon field. The importance of this term for the existence of a ground state
is discussed in Ref Rokaj et al. (2018).

###  II.1 Separation of the relative and center of mass equations

For simplicity we only consider a single photon mode. The formalism can be
easily extended to many photon modes as described in Appendix A and Appendix
B. We will define the coupling strength as
$\boldsymbol{\lambda}=(\lambda,\lambda,0)$. A more general case is described
in Appendix B. In this section, we consider the Hamiltonian that acts only in
the electron space

$\displaystyle H_{e}+$ $\displaystyle H_{d}$
$\displaystyle=-\frac{1}{2}\nabla_{1}^{2}+\frac{1}{2}\omega_{0}^{2}\mathbf{r}_{1}^{2}-\frac{1}{2}\nabla_{2}^{2}+\frac{1}{2}\omega_{0}^{2}\mathbf{r}_{2}^{2}$
(8) $\displaystyle+$
$\displaystyle\frac{q_{1}q_{2}}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}+{1\over
2}(q_{1}\boldsymbol{\lambda}\cdot\mathbf{r}_{1}+q_{2}\boldsymbol{\lambda}\cdot\mathbf{r}_{2})^{2}.\
$

Atomic units $\hbar=m=e=1$ are used throughout and unit charges are assumed.

Defining relative and CM coordinates as

$\begin{array}[]{l}\mathbf{r}=\mathbf{r}_{2}-\mathbf{r}_{1},\\\
\mathbf{R}=\frac{1}{2}\left(\mathbf{r}_{1}+\mathbf{r}_{2}\right),\end{array}$
(9)

the Hamiltonian decouples into a relative and CM Hamiltonian

$\displaystyle H_{e}+H_{d}$ $\displaystyle=$
$\displaystyle-\nabla_{\mathbf{r}}^{2}+\frac{1}{4}\omega_{0}^{2}\mathbf{r}^{2}+\frac{q_{1}q_{2}}{r}-\frac{1}{4}\nabla_{\mathbf{R}}^{2}+\omega_{0}^{2}\mathbf{R}^{2}$
(10) $\displaystyle+$ $\displaystyle{1\over
2}\left(\boldsymbol{\lambda}\cdot\left((q_{1}+q_{2})\mathbf{R}+{1\over
2}(q_{1}-q_{2})\mathbf{r}\right)\right)^{2}$ $\displaystyle\equiv$
$\displaystyle H_{\mathbf{r}}+H_{\mathbf{R}},$

and the corresponding eigenvalue problem is

$\left(H_{\mathbf{r}}+H_{\mathbf{R}}\right)\Phi(\mathbf{r},\mathbf{R})=(\epsilon+\eta)\Phi(\mathbf{r},\mathbf{R}),$
(11)

and $E=\epsilon+\eta$ is the eigenenergy. Note that for like charges, the last
term only contributes to $H_{\mathbf{R}}$, otherwise, it only contributes to
$H_{\mathbf{r}}$ and there is no cross term between $\mathbf{R}$ and
$\mathbf{r}$.

####  II.1.1 $q_{1}=-q_{2}$

In this case the photon only couples to $\mathbf{r}=(x,y,z)$. The CM wave
function is a harmonic oscillator eigenfunction with frequency $2\omega_{0}$.
By introducing $u={x+y\over\sqrt{2}}$, and $v={-x+y\over\sqrt{2}}$, the
relative motion Hamiltonian takes the form

$\displaystyle H_{\mathbf{r}}=-$
$\displaystyle\nabla_{u}^{2}-\nabla_{v}^{2}-\nabla_{z}^{2}+{1\over
2}\omega_{u}^{2}u^{2}+{1\over 2}\omega_{v}^{2}v^{2}+{1\over
2}\omega_{z}^{2}z^{2}$ (12)
$\displaystyle-{1\over{(u^{2}+v^{2}+z^{2})}^{1/2}},$

where ${\omega_{u}}^{2}=2\lambda^{2}+\frac{1}{2}{\omega_{0}}^{2},\
{\omega_{v}}^{2}={\omega_{z}}^{2}=\frac{1}{2}{\omega_{0}}^{2}$. This is a
single particle Coulomb problem in an anisotropic harmonic potential. The
derivation is detailed in the next section. We are not aware of any existing
analytical solutions to this system. One can, in principle, solve this problem
using a product basis of the $u-v-z$ harmonic oscillators, but we do not
pursue this case any further in this paper.

####  II.1.2 $q_{1}=q_{2}$

In the following, we will consider $q_{1}=q_{2}$ because, in this case, the
equation for the relative motion can be analytically found for certain
frequencies as mentioned before. After multiplying the relative part by 1/2
and the CM part by 2 to bring the equations in a more convenient form, we have

$\left[-\frac{1}{2}\nabla_{\mathbf{r}}^{2}+\frac{1}{2}\omega_{\mathbf{r}}^{2}\mathbf{r}^{2}+\frac{1}{2}\frac{1}{r}\right]\varphi(\mathbf{r})=\varepsilon^{\prime}\varphi(\mathbf{r}),$
(13)

where $\omega_{\mathrm{r}}=\frac{1}{2}\omega_{0}$ and
$\varepsilon^{\prime}=\frac{1}{2}\varepsilon$, and

$\left[-\frac{1}{2}\nabla_{\mathbf{R}}^{2}+\frac{1}{2}\omega_{\mathbf{R}}^{2}\mathbf{R}^{2}+4(\boldsymbol{\lambda}\cdot\mathbf{R})^{2}\right]\xi(\mathbf{R})=\eta^{\prime}\xi(\mathbf{R}),$
(14)

where $\omega_{\mathbf{R}}=2\omega_{0}$ and $\eta^{\prime}=2\eta$. The total
wave function can be written as

$\Phi(\mathbf{r},\mathbf{R})=\varphi(\mathbf{r})\xi(\mathbf{R}).$ (15)

In this case, the CM motion in the $z$-direction is described by a harmonic
oscillator eigenfunction, and we drop this part from now.

In 2D, using $\mathbf{R}=(X,Y)$ one can rewrite $H_{\mathbf{R}}$ as (in 3D one
simply has to multiply the CM wave function with a harmonic oscillator
function of frequency $2\omega_{0}$ in the $Z$ direction)

$H_{\mathbf{R}}=-{1\over 2}{\partial^{2}\over\partial X^{2}}-{1\over
2}{\partial^{2}\over\partial Y^{2}}+{1\over 2}\omega_{X}^{2}X^{2}+{1\over
2}\omega_{Y}^{2}Y^{2}+{1\over 2}\omega_{XY}^{2}XY,$ (16)

where

$\omega_{X}^{2}=\omega_{Y}^{2}=\omega_{\mathbf{R}}^{2}+8\lambda^{2},\ \ \ \
\omega_{XY}^{2}=16\lambda^{2}.$ (17)

Using a unitary transformation (a generalized version is presented in Appendix
A)

$U={X+Y\over\sqrt{2}},\ \ \ \ V={-X+Y\over\sqrt{2}},$ (18)

we have

$\displaystyle H_{\mathbf{R}}$ $\displaystyle=$ $\displaystyle-{1\over
2}{\partial^{2}\over\partial U^{2}}-{1\over 2}{\partial^{2}\over\partial
V^{2}}+{1\over 2}\omega_{U}^{2}U^{2}+{1\over 2}\omega_{V}^{2}V^{2}$ (19)
$\displaystyle\equiv$ $\displaystyle H_{U}+H_{V},$

where

$\omega_{U}^{2}=\frac{1}{2}(\omega_{X}^{2}+\omega_{XY}^{2}+\omega_{Y}^{2})=\omega_{R}^{2}+16\lambda^{2},$
(20)
$\omega_{V}^{2}=\frac{1}{2}(\omega_{X}^{2}-\omega_{XY}^{2}+\omega_{Y}^{2})=\omega_{R}^{2}.$
(21)

This Hamiltonian is analytically solvable: the lowest energy is

$\eta={1\over 2}\left(\omega_{0}+\sqrt{\omega_{0}^{2}+4\lambda^{2}}\right).$
(22)

$H_{\mathbf{r}}$ is also analytically solvable, in this case, for certain
frequencies Taut (1993); Sup (2021). For example, for $\omega_{0}=1$ one gets
$\epsilon=2$ (see the Table in Ref. Sup (2021)) and the total energy is
$E=2+{1\over 2}+{1\over 2}\sqrt{1+4\lambda^{2}}$.

The wave function of the CM motion now can be written as

$\xi(\mathbf{R})=\phi_{k}({U})\phi_{l}(V),$ (23)

where $\phi_{k}$ is the $k$th eigenfunction of the one-dimensional harmonic
oscillator,

$\phi_{k}(U)=\left(\sqrt{\omega_{U}}\over\sqrt{\pi}\,2^{k}k!\right)^{1\over
2}{\rm e}^{-{\omega_{U}\over 2}U^{2}}H_{k}(\sqrt{\omega_{U}}\,U),$ (24)

and the eigenfunctions are similarly defined for $V$.

###  II.2 Photon-electron coupling

The coupling term Eq. (5) takes the form

$H_{ep}=-\sqrt{{\omega\over 2}}(\hat{a}+\hat{a}^{+})\lambda D,\ \ \ \
D=2\sqrt{2}U,$ (25)

so only the $U$ harmonic oscillators are coupled with photons. The Hamiltonian
that we have to solve is reduced to a single one-dimensional electronic
Hamiltonian coupled to light:

$H_{c}=H_{U}+\omega\left(\hat{a}^{+}\hat{a}+{1\over
2}\right)-2\omega\sqrt{2}\lambda Uq.$ (26)

This Hamiltonian can be solved by exact diagonalization using the basis states

$\phi_{k}(U)|n\rangle.$ (27)

For the diagonalization, one needs the matrix elements of the Hamiltonian
which are readily available. The operators $H_{U}$ and $U$ act on the real
space, and $\hat{a}+\hat{a}^{+}$ acts on the photon space. For the coupling
term in the photon space:

$\displaystyle q|n\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega}}\left(\hat{a}+\hat{a}^{+}\right)|n\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\omega}}\left(|\sqrt{n}|n-1\rangle+\sqrt{n+1}|n+1\rangle\right),$

and the matrix elements of $q$ are

$\langle m|q|n\rangle=\frac{1}{\sqrt{2\omega}}D_{mn},$ (29)

where

$D_{mn}=\left(\begin{array}[]{cccccc}0&\sqrt{1}&0&0&0&\ldots\\\
\sqrt{1}&0&\sqrt{2}&0&0&\ldots\\\ 0&\sqrt{2}&0&\sqrt{3}&0&\ldots\\\
0&0&\sqrt{3}&0&\sqrt{4}&\ldots\\\ 0&0&0&\sqrt{4}&0&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right).$ (30)

The Hamiltonian $H_{U}$ is diagonal in the harmonic oscillator bases

$\langle\phi_{i}|H_{U}|\phi_{j}\rangle=(j+{1\over 2})\omega_{U}\delta_{ij}.$
(31)

The matrix elements of the photon Hamiltonian are

$\langle n|\omega\left(\hat{a}^{+}\hat{a}+{1\over
2}\right)|m\rangle=(n+{1\over 2})\omega\delta_{nm}.$ (32)

The last piece is the matrix elements of the position operator in harmonic
oscillator bases:

$\langle\phi_{i}|U|\phi_{j}\rangle=\frac{1}{\sqrt{2\omega_{U}}}D_{ij}.$ (33)

Thus, the matrix elements of $H_{c}$ are

$\langle
m,\phi_{i}|H_{c}|n,\phi_{j}\rangle=\delta_{mn}\delta_{ij}(j+\frac{1}{2})\omega_{U}+\delta_{mn}\delta_{ij}(n+\frac{1}{2})\omega+\sqrt{\frac{2\omega}{\omega_{U}}}\,\lambda\,D_{mn}D_{ij}.$
(34)

This is a very sparse matrix and can be diagonalized with sparse matrix
approaches even for very large dimensions. In practice, a few dozen photon
bases $|n\rangle$ and harmonic oscillator bases $\phi_{i}$ give converged
energies. This matrix is generalized for $N_{p}$ photon modes in Appendix A.

After the diagonalization, we have the eigenenergies $\eta_{j}^{\prime}$ and
the eigenfunctions by defining the spatial wave function in the photon
subspace $n$ as

$\psi_{j}(\mathbf{R})=\phi_{0}(V)\phi_{j}(U),$ (35)

the eigenfunction for the CM motion is

$\displaystyle\xi_{k}(\mathbf{R})$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{K_{n}}\left(\sum_{j=0}^{K_{U}}c^{k}_{j,n}\psi_{n}(\mathbf{R})\right)|n\rangle$
(36) $\displaystyle=$
$\displaystyle\sum_{j=0}^{K_{U}}\left(\sum_{n=0}^{K_{n}}c^{k}_{j,n}|n\rangle\right)\psi_{j}(\mathbf{R}),$

where $K_{U}$ and $K_{n}$ are some suitably chosen upper limits that control
the convergence of the eigenvalues. For the $V$ part of the CM motion, we have
chosen the lowest state. The first line in Eq. (36) emphasizes the coupling of
the spatial part to photon spaces; the second line emphasizes the coupling of
the linear combination of photon states to a given CM eigenfunction.

##  III Results and discussion

A few examples will be presented in this section. For these calculations, we
have picked an oscillator frequency $\omega_{0}$ from the Table of Ref. Sup
(2021), calculated the radial part of the relative wave function as described
in Refs. Sup (2021); Taut (1993, 1994), and multiplied with the corresponding
spherical function. This function is then multiplied by $\xi_{k}(\mathbf{R})$
calculated using Eq. (36).

First, we show the wave function for the different CM excitation. In this
case, two variables determine the behavior: the confining strength
$\omega_{0}$ and the coupling parameter $\lambda$. We show 2D examples because
they are easier to visualize. The 3D cases are very similar, with the only
difference being that the wave function is multiplied by the lowest harmonic
oscillator function with frequency $2\omega_{0}$ in the $Z$ direction.

First, we show the spin-singlet case using $\omega_{0}=1$. The energy of the
relative motion is $\epsilon=1$ a.u. in this case (see the Table in Ref. Sup
(2021)). The state with $j=0$ CM wave function is spherically symmetric for
small $\lambda$ (Fig. 1a), and as the anharmonicity of the harmonic oscillator
dominates ($\omega_{V}<<\omega_{U}$), a slightly ellipsoidal structure appears
(Fig. 1b). For $j=1$, the CM state is multiplied by $U$
($H_{1}(\sqrt{\omega_{U}}\,U)=2\sqrt{\omega_{U}}\,U$) and becomes elongated in
the diagonal direction (Fig. 1c). This direction is set by the choice of
$\boldsymbol{\lambda}=(\lambda,\lambda$), and other values would change the
direction (see Appendix A). For larger $\lambda$, the confinement by
$\omega_{U}$ is much stronger and the elongation disappears (Fig. 1d). For
higher $j$ values the elongation increases due to the higher
$H_{j}(\sqrt{\omega_{U}}U)$ polynomials (Figs. 1e and 1g). Higher $\lambda$
values decrease the elongation (Figs. 1f and 1h). This trend continues for
even higher $j$ values as well. Solutions with other $\omega_{0}$ values show
very similar behaviors.

Figure 1: Two-dimensional densities
$\left(\phi(\mathbf{r})\psi_{j}(\mathbf{R})\right)^{2}$ of two electrons
confined by a harmonic potential with $\omega_{0}=1$ a.u. and spin $S=0$.
First row: j=0, (a) $\lambda=0.5$, (b) $\lambda=2$. Second row: j=1, (c)
$\lambda=0.5$, (d) $\lambda=2$. Third row: j=2, (e) $\lambda=0.5$, (f)
$\lambda=2$. Fourth row: j=5, (g) $\lambda=0.5$, (h) $\lambda=2$.

In the spin-triplet case in 2D we choose $\omega_{0}=1/3$ a.u., and the energy
of the relative motion is $\epsilon=1$ a.u. Fig. 2 shows the densities for
this case. This system is more sensitive to the choice of $\lambda$ and we use
three different $\lambda$ values (0.01, 0.5, 2) to illustrate that. In this
case, the spin function is symmetric, the spatial part is antisymmetric
($m=\pm 1$ in Eq. (1.2) in Ref. Sup (2021)) and the two peaks appear in the
density plot for $j=0$ as shown in Figs. 2a, 2b and 2c. By increasing
$\lambda$ the $U$ oscillator squeezes the electrons closer and the separation
between the two peaks is more visible (the density between the peaks being
lower). There are three peaks for $j=1$ for $\lambda=0.01$ and $\lambda=0.5$,
but as the $U$ confinement gets stronger the two peak structure returns (Figs.
2d, 2e, and 2f). The three-peak structure can be a nontrivial consequence
because, unlike the simple spherical structure in the singlet state, the
relative motion function, in this case, is in an $m=1$ state and multiplied by
$U$. If we neglect the Coulomb interaction, then the relative motion wave
function is a ground state harmonic oscillator for the first electron (one
density peak) and the first excited harmonic oscillator state for the second
electron (two peaks). These three peaks are magnified when the relative wave
function is multiplied by the center of mass wave function, which is now
proportional to $H_{1}(\sqrt{\omega_{U}}\,U)=2\sqrt{\omega_{U}}\,U$.

For higher $j$ states, the elongation caused by $H_{j}$ continues (see Figs.
2g and 2k), and the nodal structure of $H_{j}$ also contributes to the density
structure. Overall it seems that the $\lambda=0.01$ case captures the general
trend very well. For larger $\lambda$ values the same structures appear later
as $j$ increases. For a given $j$, increasing $\lambda$ squeezes the
elongation due to the $\omega_{U}$ confinement, as in the singlet case.

Figure 2: Two-dimensional densities
$\left(\phi(\mathbf{r})\psi_{j}(\mathbf{R})\right)^{2}$ of two electrons
confined by a harmonic potential with $\omega_{0}=1/3$ a.u. and spin $S=1$.
First row: $j=0$, (a) $\lambda=0.01$, (b) $\lambda=0.5$, (c) $\lambda=2$.
Second row: $j=1$, (d) $\lambda=0.01$, (e) $\lambda=0.5$, (f) $\lambda=2$.
Third row: $j=2$, (g) $\lambda=0.01$, (h) $\lambda=0.5$, (j) $\lambda=2$.
Third row: $j=5$, (k) $\lambda=0.01$, (l) $\lambda=0.5$, (m) $\lambda=2$.

The total wave function will be a linear combination of the
$\phi(\mathbf{r})\psi_{j}(\mathbf{R})|n\rangle$ components. The probability of
a given component is given by

$P_{j}(n)=|c_{jn}^{0}|^{2}$ (37)

and depends on $\omega_{0},\omega$ and $\lambda$. An example for $P_{j}(n)$ is
given in Fig. 3a. First, we note that due to the structure in Eq. (34), the
probabilities follow a checkerboard-like structure: odd photon numbers couple
to odd $j$ and even photon numbers couple to even $j$. The probabilities
decrease for large photon numbers. The low CM excitations, $j=0,1,3$, are the
most dominant terms for low photon numbers. The probability of the higher CM
excitations ($j=2,3,4,5$) first increases with the photon number, then reaches
a maximum and starts to decrease.

Fig. 3b shows the sum $P_{n}=\sum_{j}P_{n}(j)$. By increasing $\lambda$, the
higher photon spaces are coupled and the occupation of lower photon numbers
increases. However, the effect of $\omega$ is more complicated. The coupling
increases as $\sqrt{\omega}$ but with larger photon frequency the photon
harmonic oscillator states move higher in energy ($n\hbar\omega$) and their
occupation decreases. This latter effect seems to be dominant for smaller
$\lambda$. In Fig. 3b, in case of $\lambda=0.5$ a.u., $P_{n}$ is the same for
$n=0,1$ for $\omega=0.5$ a.u. and $\omega=5$ a.u., but for higher $n$, $P_{n}$
is much smaller for $\omega=5$ a.u. For higher $\lambda$ values this effect
becomes less important. The oscillations (the even states have higher
occupation than the odd states) in the case of $\omega$=0.5, $\lambda$=5 a.u.
always appear when $\omega$ is much smaller than $\lambda$ and probably due to
the checkerboard-like coupling.

Fig. 3c shows the sum $P_{j}=\sum_{n}P_{n}(j)$. By increasing $\omega$, the
occupations of the low photon number states decrease and the occupations of
the higher states increase. Increasing $\lambda$ increases $\omega_{U}$ and
pushes the CM states higher, and those states do not couple with the low $j$
sector, so increasing $\lambda$ decreases $P_{j}$. The effect of $\lambda$ is
similar to $\omega$ in the previous case: $\lambda$ increases the coupling,
but larger $\lambda$ means larger $\omega_{u}$ and the CM states are pushed
higher. For low $\omega$, increasing $\lambda$ relaxes the occupation, but for
large $\omega$, the coupling dominates and the $\lambda$ increases the
occupation of the higher $n$ states.

Figure 3: Top: Probability of the occupation of a $jn$ subspace ($\lambda=0.5$
a.u. and $\omega=0.5$ a.u.), $P_{j}(n)$. Middle: $P_{n}=\sum_{j}P_{n}(j)$.
Bottom :$P_{j}=\sum_{n}P_{n}(j)$. $\omega_{0}=1/3$ a.u. is used in the
calculations.

Fig. 4 shows the energy of the singlet state as a function of photon
frequency. The ground state energy is 3 a.u. in this case. Infinitely many
photon states and infinitely many CM states can couple to this state. Without
coupling of the photons to the center of mass, the energy of the photon states
increases linearly with $\omega$ and the energy of the CM states increases
linearly with $\omega_{U}$. Fig. 4 shows the lowest 20 states with coupling,
and we use $\lambda=\alpha\sqrt{\omega}$, which means that the coupling $g$ is
proportional to $\omega$. For $\alpha=1$ (Fig. 4a), some states (primarily
photon states) move linearly up with $\omega$ for small frequencies, while
other states (primarily CM states) only slowly increase with $\omega$ and
converge to a horizontal line. To magnify the behavior we redo the calculation
with $\alpha=1/20$ (Fig. 4b). In this case, $\lambda$ increases much less
while we increase $\omega$. $\omega_{U}$ barely changes while $\omega$ ascends
from 0 to 5. The lowest state does not change (it barely couples to photons
$n>0$) and remains a horizontal line, which is just the lowest CM state. The
second state increases with $\omega$, but then it reaches the first excited CM
state of energy $\hbar\omega_{U}\approx\hbar\omega_{0}$ and becomes a
horizontal line. The third state also increases until it reaches the energy of
the first excited state and continues on that line until meeting the second
state. To avoid crossing, it moves to the second excited state of the center
of mass and so on. Of course, these states do not lie exactly on horizontal
lines but rise gradually with lambda. For much higher $\omega$, one can
recover a similar picture to Fig. 4a.

To complete the energy spectrum of the system, one has to include the excited
states of the relative motion. As those states are orthogonal, the complete
spectrum can be obtained by shifting the energy levels in Fig. 4 by the
energies of the excited states.

The wavefunction of the system will be a linear combination of wave functions
shown in Figs. 1 and 2, with coefficients defined in the second line of Eq.
(36). These coefficients depend on the values of $\omega$, $\omega_{0}$, and
$\lambda$. For a single photon mode, the lowest states often dominate and it
is hard to pick parameters that favor a single $j$ CM mode or higher $j$
values. In Fig. 5a we present an example for the triplet case where the square
of the linear coefficients are 0.55, 0.18, 0.10, 0.06, 0.04, 0.03 (j=0,…5), so
a few $j\neq 0$ contribute to the density. Figs. 5b, 5c and 5d show the square
of the wave function in the $n=1,3$ and 5 spaces. The $n$=0 density is very
similar to Fig. 2b. The squares of the linear coefficients in $n$ space are
0.49, 0.13, 0.11, 0.07, 0.05 ($n=0,...,5$). This example shows that the
spatial wave functions in different photon subspaces can be very different and
quantum mechanical methods have to look for accurate wave functions in
different photon spaces.

The calculation can be extended to many photon modes as it is shown in
Appendix A. Examples of two-photon mode calculations are shown in Figs. 6 and
7. In particular, Fig. 6 shows the photon occupation numbers for the two-
photon modes, $\omega$ and $2\omega$. The occupation probability tilts toward
the $\omega$ axis, showing that the $\omega$ modes have higher probabilities
than the $2\omega$ ones. Fig. 7 is the same calculation as is shown in Fig.
4a, but with two-photon modes. Overall, the two figures are very similar. The
two-photon case reaches higher energies and there are more level crossings.
This is because some of the states shown in Fig. 7 are $2\omega$ states and
move higher faster. Increasing the number of photon modes helps to reach
higher $j$ states and multiphoton modes might be a way to select higher $j$
states or single out a desired $j$ value.

Figure 4: Energy levels as a function of $\omega$ for (a)
$\lambda=\sqrt{\omega}$ and (b) $\lambda=\sqrt{\omega}/20$, $\omega_{0}$=1
a.u. in both cases.

Figure 5: Two-dimensional density $\Phi(\mathbf{r},\mathbf{R})^{2}$ for the
$S=1$ case, (a) total density, (b) density in the $n=1$ space, (c) density in
the $n=3$ space, (d) density in the $n=5$ space. ($\omega_{0}=0.18055$ a.u.,
$\omega=1$ a.u. and $\lambda=1$ a.u.). Figure 6: Occupation numbers for the
photon modes. The photon frequencies are $\omega$ and $2\omega$. The coupling
vectors are $\boldsymbol{\lambda}$ and $-\boldsymbol{\lambda}$. The vertical
axis is the photon number for $\omega$, the horizontal axis is the photon
number for $2\omega$ ($\omega_{0}$=1 a.u., $\omega=1$ a.u., $\lambda=1$ a.u.).
Figure 7: Energy levels as a function of $\omega$ for the two photon modes.
The photon frequencies are $\omega$ and $2\omega$, the coupling vectors are
$\boldsymbol{\lambda}$ and $-\boldsymbol{\lambda}$, with
$\lambda=\sqrt{\omega}$ ($\omega_{0}$=1 a.u.)

##  IV Summary

In a harmonically confined two-electron system, the light couples to the
dipole moment which is proportional to the CM coordinate. By separating the
relative and center of mass motion, we have shown that the coupled photon
center of mass system can be solved by diagonalization and the relative motion
part has analytical solutions for certain frequencies.

The coupling of the light to the center of mass coordinate leads to elongated
wave functions. The symmetry axis of the electron density is determined by the
polarization direction. The density has several peaks depending on the center
of mass excitation and the symmetry axis of the density is determined by the
polarization direction. The competition between the confinement due to the
coupling to light and the node structure of the center of mass excitation
influences the location of the density peaks.

We have shown that the spatial wave functions belonging to different photon
spaces are very different, and this means that quantum mechanical approaches
solving coupled light-matter problems have to determine the wave functions in
each photon subspace, which might be a difficult task.

The approach can be extended to many photon modes and the only limitation is
the dimension of the Hamiltonian matrix. As this matrix is very sparse, one
can easily diagonalize it even for very large matrices.

As there are only very few light-matter coupled systems with analytical
solutions, the present work might be useful to test and develop efficient
approximations.

A similar approach can be used for a larger electron number, but then the
relative motion part has to be solved numerically.

##  Appendix A $N_{p}$ photon modes

Consider the same system as in Section II.1.2, except that here $n_{p}$
photons are coupled. Hence, there are $N_{p}=2n_{p}$ photon modes involved.

$|\vec{n}\rangle=|n_{1},n_{2},{\ldots},n_{N_{p}}\rangle.$ (38)

Define the vector Kronecker delta as

$\delta_{\vec{n}\vec{m}}=\prod_{k=1}^{N_{p}}\delta_{n_{k}m_{k}},$ (39)

$\delta_{\vec{n}\vec{m}}^{l}=\prod_{k=1,k\neq l}^{N_{p}}\delta_{n_{k}m_{k}}.$
(40)

It is straightforward to generalize Eq. (34)

$\langle\vec{m},\phi_{i}|H|\vec{n},\phi_{j}\rangle=\delta_{\vec{m}\vec{n}}\delta_{ij}(j+\frac{1}{2})\omega_{U}+\delta_{\vec{m}\vec{n}}\delta_{ij}\sum_{k=1}^{N_{p}}(n_{k}+\frac{1}{2})\omega_{k}+\sum_{k=1}^{N_{p}}\sqrt{\frac{2\omega_{k}}{\omega_{U}}}\,\lambda\,D_{n_{k}m_{k}}D_{ij}\delta_{\vec{n}\vec{m}}^{k}.$
(41)

##  Appendix B Center-of-mass motion for many photons

We assume $N_{p}=2n_{p}$ photon modes, and $\boldsymbol{\lambda_{\alpha}}$’s
are not necessarily isotropic in the $x$, $y$-directions. Thus, the
Hamiltonian becomes

$\displaystyle H$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\nabla_{1}^{2}+\frac{1}{2}\omega_{0}^{2}\mathbf{r}_{1}^{2}-\frac{1}{2}\nabla_{2}^{2}+\frac{1}{2}\omega_{0}^{2}\mathbf{r}_{2}^{2}$
(42) $\displaystyle+$
$\displaystyle\frac{q_{1}q_{2}}{\left|\mathbf{r}_{1}-\mathbf{r}_{2}\right|}+{1\over
2}\sum_{\alpha=0}^{N_{p}}(q_{1}\boldsymbol{\lambda}_{\alpha}\cdot\mathbf{r}_{1}+q_{2}\boldsymbol{\lambda}_{\alpha}\cdot\mathbf{r}_{2})^{2}.$

Still imposing $q_{1}=q_{2}$, the radial part remains unchanged and can be
solved by RefTaut (1993, 1994); Sup (2021). Now we solve the CM part. Eq. (14)
becomes

$\left[-\frac{1}{2}\nabla_{\mathbf{R}}^{2}+\frac{1}{2}\omega_{\mathbf{R}}^{2}\mathbf{R}^{2}+4\sum_{\alpha=0}^{N_{p}}(\boldsymbol{\lambda}_{\alpha}\cdot\mathbf{R})^{2}\right]\xi(\mathbf{R})=\eta^{\prime}\xi(\mathbf{R}).$
(43)

Suppose $\boldsymbol{\lambda}_{\alpha}=(\lambda_{\alpha 1},\lambda_{\alpha
2},0)$. Further define

$\displaystyle\tilde{\lambda}_{1}$
$\displaystyle=\sum_{\alpha=0}^{N_{p}}\lambda_{\alpha 1}^{2},$ (44)
$\displaystyle\tilde{\lambda}_{2}$
$\displaystyle=\sum_{\alpha=0}^{N_{p}}\lambda_{\alpha 2}^{2},$
$\displaystyle\tilde{\lambda}_{12}$
$\displaystyle=\sum_{\alpha=0}^{N_{p}}\lambda_{\alpha 1}\lambda_{\alpha 2},$

Eq. (16) and (17) now read

$H_{\mathbf{R}}=-{1\over 2}{\partial^{2}\over\partial X^{2}}-{1\over
2}{\partial^{2}\over\partial Y^{2}}+{1\over 2}\omega_{X}^{2}X^{2}+{1\over
2}\omega_{Y}^{2}Y^{2}+{1\over 2}\omega_{XY}XY,$ (45)

where

$\displaystyle\omega_{X}^{2}=\omega_{\mathbf{R}}^{2}+8{\tilde{\lambda}_{1}},$
(46)
$\displaystyle\omega_{Y}^{2}=\omega_{\mathbf{R}}^{2}+8{\tilde{\lambda}_{2}},$
$\displaystyle\omega_{XY}=16\tilde{\lambda}_{12}.$

This linearly coupled Hamiltonian can be easily decoupled with the following
unitary transformation

$\displaystyle U=\frac{1}{(1-ab)^{1/2}}(X+aY),$ (47) $\displaystyle
V=\frac{1}{(1-ab)^{1/2}}(bX+Y),$

where

$\displaystyle
a=\frac{(\tilde{\lambda}_{1}-\tilde{\lambda}_{2})+\sqrt{\,(\tilde{\lambda}_{1}-\tilde{\lambda}_{2})\,^{2}+4\tilde{\lambda}_{12}^{2}}}{2\tilde{\lambda}_{12}},$
(48) $\displaystyle
b=-\frac{(\tilde{\lambda}_{1}-\tilde{\lambda}_{2})+\sqrt{\,(\tilde{\lambda}_{1}-\tilde{\lambda}_{2})\,^{2}+4\tilde{\lambda}_{12}^{2}}}{2\tilde{\lambda}_{12}}.$

In this case, the decoupled Hamiltonian reads

$H_{R}(U,V)=-\frac{1}{2}\frac{\partial^{2}}{\partial
U^{2}}-\frac{1}{2}\frac{\partial^{2}}{\partial
V^{2}}+\frac{1}{2}\omega_{U}^{2}U^{2}+\frac{1}{2}\omega_{V}^{2}V^{2},$ (49)

where

$\displaystyle\omega_{U}=\sqrt{\omega_{R}^{2}+4(\tilde{\lambda}_{1}+\tilde{\lambda}_{2})+4\sqrt{(\tilde{\lambda}_{1}-\tilde{\lambda}_{2})^{2}+4\tilde{\lambda}_{12}^{2}}},$
(50)
$\displaystyle\omega_{V}=\sqrt{\omega_{R}^{2}+4(\tilde{\lambda}_{1}+\tilde{\lambda}_{2})-4\sqrt{(\tilde{\lambda}_{1}-\tilde{\lambda}_{2})^{2}+4\tilde{\lambda}_{12}^{2}}}.$

Same as in Eq. (19), this is just the Hamiltonian for two non-interacting
harmonic oscillators. Hence, the energies for the CM part are

$\eta=\frac{1}{2}\eta^{\prime}=\frac{1}{2}(n_{U}+\frac{1}{2})\omega_{U}+\frac{1}{2}(n_{V}+\frac{1}{2})\omega_{V}\,\text{,
}n_{U},\ n_{V}=0,1,2...$ (51)

and the ground state energy is

$\eta_{0}=\frac{1}{4}(\omega_{U}+\omega_{V}).$ (52)

Finally, the corresponding wave function is just the product of that of the
two independent harmonic oscillators.

###### Acknowledgements.

This work has been supported by the National Science Foundation (NSF) under
Grant No. IRES 1826917.

## References

  * Buchholz et al. (2019) F. Buchholz, I. Theophilou, S. E. B. Nielsen, M. Ruggenthaler, and A. Rubio, ACS Photonics 6, 2694 (2019).
  * Schäfer et al. (2019) C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 116, 4883 (2019), ISSN 0027-8424, eprint https://www.pnas.org/content/116/11/4883.full.pdf, URL https://www.pnas.org/content/116/11/4883.
  * Ruggenthaler et al. (2018) M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, Nature Reviews Chemistry 2, 0118 (2018), ISSN 2397-3358, URL https://doi.org/10.1038/s41570-018-0118.
  * Flick et al. (2015) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 112, 15285 (2015), ISSN 0027-8424, eprint https://www.pnas.org/content/112/50/15285.full.pdf, URL https://www.pnas.org/content/112/50/15285.
  * Flick et al. (2017) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proceedings of the National Academy of Sciences 114, 3026 (2017), ISSN 0027-8424, eprint https://www.pnas.org/content/114/12/3026.full.pdf, URL https://www.pnas.org/content/114/12/3026.
  * Rokaj et al. (2018) V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, Journal of Physics B: Atomic, Molecular and Optical Physics 51, 034005 (2018), URL https://doi.org/10.1088/1361-6455/aa9c99.
  * Rivera et al. (2019) N. Rivera, J. Flick, and P. Narang, Phys. Rev. Lett. 122, 193603 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.193603.
  * Flick and Narang (2018) J. Flick and P. Narang, Phys. Rev. Lett. 121, 113002 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.121.113002.
  * Hoffmann et al. (2020) N. M. Hoffmann, L. Lacombe, A. Rubio, and N. T. Maitra, The Journal of Chemical Physics 153, 104103 (2020).
  * Tokatly (2018) I. V. Tokatly, Phys. Rev. B 98, 235123 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.98.235123.
  * Galego et al. (2017) J. Galego, F. J. Garcia-Vidal, and J. Feist, Phys. Rev. Lett. 119, 136001 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.136001.
  * Mandal et al. (2020a) A. Mandal, S. Montillo Vega, and P. Huo, The Journal of Physical Chemistry Letters 11, 9215 (2020a), pMID: 32991814.
  * Cederbaum and Kuleff (2021) L. S. Cederbaum and A. I. Kuleff, Nature Communications 12, 4083 (2021), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-021-24221-6.
  * Szidarovszky et al. (2018) T. Szidarovszky, G. J. Halász, A. G. Császár, L. S. Cederbaum, and A. Vibók, The Journal of Physical Chemistry Letters 9, 6215 (2018).
  * Ashida et al. (2021) Y. Ashida, A. m. c. İmamoğlu, and E. Demler, Phys. Rev. Lett. 126, 153603 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.153603.
  * Rivera and Kaminer (2020) N. Rivera and I. Kaminer, Nature Reviews Physics 2, 538 (2020), ISSN 2522-5820, URL https://doi.org/10.1038/s42254-020-0224-2.
  * Le Boité (2020) A. Le Boité, Advanced Quantum Technologies 3, 1900140 (2020), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/qute.201900140, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qute.201900140.
  * Garcia-Vidal et al. (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373 (2021), ISSN 0036-8075, eprint https://science.sciencemag.org/content/373/6551/eabd0336.full.pdf, URL https://science.sciencemag.org/content/373/6551/eabd0336.
  * Taut (1993) M. Taut, Phys. Rev. A 48, 3561 (1993), URL https://link.aps.org/doi/10.1103/PhysRevA.48.3561.
  * Taut (1994) M. Taut, Journal of Physics A: Mathematical and General 27, 1045 (1994), URL https://iopscience.iop.org/article/10.1088/0305-4470/27/3/040/pdf.
  * Taut (1995) M. Taut, Journal of Physics A: Mathematical and General 28, 2081 (1995), URL https://iopscience.iop.org/article/10.1088/0305-4470/28/7/026/pdf.
  * Karwowski and Pestka (2007) J. Karwowski and G. Pestka, Theoretical Chemistry Accounts 118 (2007), URL https://link.springer.com/content/pdf/10.1007/s00214-007-0362-y.pdf.
  * Villalba and Pino (1998) V. M. Villalba and R. Pino, Physics Letters A 238, 49 (1998), URL https://www.sciencedirect.com/science/article/pii/S0375960197008918.
  * Karwowski (2008) J. Karwowski, Journal of Physics: Conference Series 104 (2008), URL https://iopscience.iop.org/article/10.1088/1742-6596/104/1/012033/pdf.
  * Turbiner (1988) A. Turbiner, Communications in Mathematical Physics 118, 467 (1988), URL https://projecteuclid.org/journals/communications-in-mathematical-physics/volume-118/issue-3/Quasi-exactly-solvable-problems-and-rm-sl2-algebra/cmp/1104162094.full.
  * Liu and Hao (2015) L. Liu and Q. Hao, Theoretical and Mathematical Physics 183 (2015), URL https://link.springer.com/content/pdf/10.1007/s11232-015-0291-1.pdf.
  * Downing and Portnoi (2016) C. A. Downing and M. E. Portnoi, Physical Review B 94 (2016), URL https://journals.aps.org/prb/pdf/10.1103/PhysRevB.94.045430.
  * Ho and Khalilov (2000) C.-L. Ho and V. R. Khalilov, Phys. Rev. A 61 (2000), URL https://journals.aps.org/pra/pdf/10.1103/PhysRevA.61.032104.
  * Chiang and Ho (2002) C.-M. Chiang and C.-L. Ho, Journal of Mathematical Physics 43, 43 (2002), URL https://aip.scitation.org/doi/abs/10.1063/1.1418426.
  * Agboola (2012) Y.-Z. Agboola, Davis; Zhang, Modern physics letters A 27 (2012), URL https://www.worldscientific.com/doi/pdf/10.1142%2FS021773231250112X.
  * Karwowski and Witek (2016) J. Karwowski and H. Witek, Molecular Physics 114 (2016), URL https://www.tandfonline.com/doi/pdf/10.1080/00268976.2015.1115565.
  * Akhmedov and Guseinova (2009) K. Akhmedov and N. Guseinova, Russian Physics Journal 52, 321 (2009), URL https://link.springer.com/content/pdf/10.1007/s11182-009-9230-7.pdf.
  * Pont et al. (2018) F. M. Pont, O. Osenda, and P. Serra, Journal of Physics A: Mathematical and General 51 (2018), URL https://iopscience.iop.org/article/10.1088/1751-8121/aab85e/pdf.
  * Downing (2017) C. A. Downing, Physical Review A 95 (2017), URL https://journals.aps.org/pra/pdf/10.1103/PhysRevA.95.022105.
  * Loos and Gill (2012) P.-F. Loos and P. M. W. Gill, Physical Review Letters 108 (2012), URL https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.108.083002.
  * Guo et al. (2012) G.-J. Guo, Z.-Z. REN, B. ZHOU, and X.-Y. GUO, International Journal of Modern Physics B 26 (2012), URL https://www.worldscientific.com/doi/pdf/10.1142%2FS0217979212502013.
  * Loos and Gill (2009) P.-F. Loos and P. M. W. Gill, Physical Review Letters 103 (2009), URL https://doi.org/10.1103/PhysRevLett.103.123008.
  * Loos and Gill (2010) P.-F. Loos and P. M. W. Gill, Molecular Physics 108 (2010), URL https://www.tandfonline.com/doi/pdf/10.1080/00268976.2010.508472.
  * Jestädt et al. (2019) R. Jestädt, M. Ruggenthaler, M. J. T. Oliveira, A. Rubio, and H. Appel, Advances in Physics 68, 225 (2019).
  * Frisk Kockum et al. (2019) A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nature Reviews Physics 1, 19 (2019), ISSN 2522-5820, URL https://doi.org/10.1038/s42254-018-0006-2.
  * Schäfer et al. (2018) C. Schäfer, M. Ruggenthaler, and A. Rubio, Phys. Rev. A 98, 043801 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.043801.
  * Flick et al. (2018) J. Flick, C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, ACS Photonics 5, 992 (2018).
  * Sidler et al. (2020) D. Sidler, M. Ruggenthaler, H. Appel, and A. Rubio, The Journal of Physical Chemistry Letters 11, 7525 (2020), pMID: 32805122, URL https://doi.org/10.1021/acs.jpclett.0c01556.
  * Lacombe et al. (2019) L. Lacombe, N. M. Hoffmann, and N. T. Maitra, Phys. Rev. Lett. 123, 083201 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.123.083201.
  * Andolina et al. (2019) G. M. Andolina, F. M. D. Pellegrino, V. Giovannetti, A. H. MacDonald, and M. Polini, Phys. Rev. B 100, 121109 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.121109.
  * Schuler et al. (2020) M. Schuler, D. D. Bernardis, A. M. Läuchli, and P. Rabl, SciPost Phys. 9, 66 (2020), URL https://scipost.org/10.21468/SciPostPhys.9.5.066.
  * Settineri et al. (2021) A. Settineri, O. Di Stefano, D. Zueco, S. Hughes, S. Savasta, and F. Nori, Phys. Rev. Research 3, 023079 (2021), URL https://link.aps.org/doi/10.1103/PhysRevResearch.3.023079.
  * Rokaj et al. (2021) V. Rokaj, M. Ruggenthaler, F. G. Eich, and A. Rubio, _The free electron gas in cavity quantum electrodynamics_ (2021), eprint 2006.09236.
  * Shin and Metiu (1995) S. Shin and H. Metiu, The Journal of Chemical Physics 102, 9285 (1995), URL https://doi.org/10.1063/1.468795.
  * Mandal et al. (2020b) A. Mandal, T. D. Krauss, and P. Huo, The Journal of Physical Chemistry B 124, 6321 (2020b), pMID: 32589846.
  * Power et al. (1959) E. A. Power, S. Zienau, and H. S. W. Massey, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 251, 427 (1959), eprint https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.1959.0008, URL https://royalsocietypublishing.org/doi/abs/10.1098/rsta.1959.0008.
  * Ruggenthaler et al. (2014) M. Ruggenthaler, J. Flick, C. Pellegrini, H. Appel, I. V. Tokatly, and A. Rubio, Phys. Rev. A 90, 012508 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.90.012508.
  * Sup (2021) (2021), see Supplemental Material at.