# Symmetry breaking for semiconductor excitons induced by Coulomb coupling
between heavy and light holes

Shiue-Yuan Shiau1, Benoit Eble2 and Monique Combescot2 (1) Physics Division,
National Center for Theoretical Sciences, Taipei, 10617, Taiwan  (2) Institut
des NanoSciences de Paris, Sorbonne Université, CNRS, 75005 Paris, France

###### Abstract

Semiconductor excitons are commonly seen as hydrogen atom. This analogy
requires a unique hole mass. In reality, this is not so due to the complexity
of the semiconductor band structure. The precise consequences on the exciton
physics of the Coulomb coupling between heavy and light holes remain a tricky
open problem. Through an “optimized perturbative” approach that uses excitons
with a flexible hole mass as a basis, we show that for zero exciton wave
vector, the heavy-light hole mass difference does not split the $(2\times 4)$
exciton degeneracy in zinc-blende-like semiconductors, the hole mass for
binding energy being close to the average mass inverse. By contrast, for
nonzero exciton wave vector, that physically breaks the crystal symmetry, the
exciton degeneracy splits into two branches quantized along the exciton wave
vector, with nontrivial center-of-mass dependence not only on the heavy and
light hole masses, but also on the electron mass.

The impact of excitons on linear and nonlinear optics of semiconductors is of
paramount importanceKirabook  for modern technology. Powered by mature
conventional heterostructures and novel low-dimensional materials, a new
momentum has recently gathered to exploit excitons for sensingFeierabend2017 ,
quantum simulationLagoin2022 , quantum memory and quantum circuitsKrenner2008
; Cong2020 .

The exciton is a composite quantum particle usually seen as one negatively-
charged conduction electron and one positively-charged valence hole correlated
by Coulomb attraction, in this way sharing similarity with the hydrogen atom.
Driven by this attraction, bound-state excitons show up as pronounced narrow
absorption lines below the band gap.

Such a simplified exciton picture however has to be questioned when facing
cutting-edge challenges in the emerging field of “Excitonics”Kis_2019  because
the hole definitely is a very tricky quantum object: it fundamentally
corresponds to an electron absence in the full valence band, with spin and
spatial degeneracies mixed by the spin-orbit interaction. The resulting
valence band complexity defies the gross exciton reduction to hydrogen-like
eigenstates. Even within the spherical approximation valid for zinc-blende-
like cubic semiconductors, there are still two holes with different masses,
coupled by Coulomb scatterings in a far from trivial wayShiauprb21 . These two
masses prevent using the hydrogen procedure — that fundamentally transforms
the two-body problem into two separate one-body problems — to solve the
resulting Schrödinger equation exactly, as required to derive bound states. A
procedure, completely different from the very first line, must be found.

The fact that heavy holes can turn light under Coulomb interaction, leads us
to anticipate that for zero exciton center-of-mass wave vector, the heavy and
light hole masses should appear in the exciton binding energy, through a
unique averaged value because in the absence of cubic-symmetry breaking, there
is no physical reason for the exciton energy to split. By contrast, a nonzero
exciton wave vector that produces such a symmetry breaking, should lead to a
splitting of the exciton energy into two branches with different center-of-
mass masses.

We here show that the exciton binding energy is driven by a single hole mass
value, so that the exciton degeneracy is preserved despite the hole mass
difference. This single mass is close to $m^{\ast}_{h}$, obtained by averaging
the hole mass inverses (see Eq. (1)), not only for small hole mass difference,
but also for light conduction electron. Yet, this single hole mass does not
describe the exciton motion: when the exciton wave vector ${\bf K}$ differs
from zero, the exciton splits into two branches in which the heavy and light
holes appear through the two ($\pm 3/2$) and ($\pm 1/2$) linear combinations
quantized along ${\bf K}$ in a way similar to the heavy-light hole splitting
induced by the spin-orbit interaction in the valence band.

The major difficulty of this problem is to handle the Coulomb interaction
exactly — as required for bound states — while the extremely complicated
heavy-light Coulomb scatteringsShiauprb21  deprive any hope to analytically
solve the corresponding Schrödinger equation. The trick we have found to
overcome this difficulty, is to introduce an exciton basis constructed on
holes having a single, flexible hole mass $m_{h}$. Its states are obtained
from a hydrogen-like procedure, with the Coulomb interaction handled exactly.
We look for the consequences of having two hole masses different from $m_{h}$,
by treating the residual one-body hole kinetic term through a perturbative
approach. The $m_{h}$ hole mass is ultimately adjusted for the resulting
exciton ground-state energy to be minimum. The procedure we propose, that can
be qualified as an “optimized perturbative approach”, is similar to the
optimized $\delta$ expansionexpansion1 ; expansion2  that aims at obtaining
better than conventional perturbative results. This procedure allows us to
catch the physics of the heavy-light hole problem. The result we find reads in
terms of the hole mass difference
$(m_{{}_{H}}-m_{{}_{L}})/(m_{{}_{H}}+m_{{}_{L}})$, within a
$m_{e}/(m_{e}+m^{\ast}_{h})$ prefactor, for the electron mass has to show up
in some way when dealing with exciton.

The problem – We consider a zinc-blende-like cubic semiconductor with valence
band having a threefold spatial degeneracy labeled as $\mu=(x,y,z)$ along the
crystal axes. The spin-orbit interaction mixes the $(3\times 2)$ hole states
$|\mu\rangle\otimes|s\rangle$ and splits their subspace into a fourfold level
$\mathcal{J}_{\textbf{z}}=(\pm 3/2,\pm 1/2)$ quantized along the crystal axis
z, and a twofold level that we here neglect as it is far below in
energyCardona . The ${\bf k}\cdot{\bf p}$ coupling of the fourfold valence
electrons to the conduction levelsFishman  produces two energy dispersions
associated with heavy and light masses $(m_{{}_{H}},m_{{}_{L}})$ in the
spherical approximation that neglects the warpingDresselhausPR1955 ;
Dresselhaus1956 ; Luttinger ; Baldereschi_prl ; Baldereschi_prb . The
resulting heavy and light holes are labeled by $\mathcal{J}_{{\bf k}_{h}}=(\pm
3/2)$ and $\mathcal{J}_{{\bf k}_{h}}=(\pm 1/2)$ indices that are quantized
along the hole wave vector ${\bf k}_{h}$.

Figure 1: A “in” heavy hole, $\mathcal{J}_{{\bf k}_{h}}=\pm 3/2$, can stay
heavy, $\mathcal{J}_{{\bf k}^{\prime}_{h}}=\pm 3/2$, or turn light,
$\mathcal{J}_{{\bf k}^{\prime}_{h}}=\pm 1/2$, under a Coulomb interaction. The
associated scattering depends on the angle between the “in” and “out” hole
wave vectors $({\bf k}_{h},{\bf k}^{\prime}_{h})$ in a very complicated
wayShiauprb21 .

The Coulomb interaction, diagonal for holes in the spatial Bloch states
labeled by $\mu$, stays diagonalMonicbook  for spin-orbit holes labeled by
$\mathcal{J}_{\textbf{z}}$, but does not stay diagonalShiauprb21  for
$\mathcal{J}_{{\bf k}_{h}}$ holes because the hole wave vector ${\bf k}_{h}$
changes in a Coulomb scattering (Fig. 1). This is why working with different
heavy and light masses makes the exciton problem in bulk samples tremendously
difficult. Note that this difficulty does not exist in quantum wells because
the confinement energy brings the light holes far below in energy, so that we
can neglect their roledense ou puits .

It is then not a surprise that when dealing with bulk exciton, a unique mass
is assigned to the hole, its common valueDresselhausPR1955 ; Dresselhaus1956
reading in terms of heavy and light masses as

$\frac{1}{m^{\ast}_{h}}=\frac{1}{2}\left(\frac{1}{m_{{}_{H}}}+\frac{1}{m_{{}_{L}}}\right)$
(1)

While averaging mass inverses is reasonable for exciton binding because this
binding depends on the hole mass as $1/m_{h}$, such an averaging is hard to
accept for exciton center of mass because its hole part would then tend to
$2m_{{}_{L}}$ for $m_{{}_{H}}\gg m_{{}_{L}}$.

As the Coulomb interaction couples heavy and light holes, a heavy hole can
turn light in a Coulomb scattering, with an amplitude that heavily
dependsShiauprb21  on the angle between the incoming and outgoing wave
vectors, ${\bf k}_{h}$ and ${\bf k}^{\prime}_{h}$ (Fig. 1). For holes having
different masses, the exciton that results from the repeated Coulomb
scatterings between electron and hole is not analytically solvable for a very
simple reason: it is no more possible to reduce the two-body electron-hole
problem to two one-body problems. Indeed, the kinetic energy of an electron-
hole pair with masses $(m_{e},m_{h})$ can be rewritten as

$\frac{\hbar^{2}{\bf k}_{e}^{2}}{2m_{e}}+\frac{\hbar^{2}{\bf
k}_{h}^{2}}{2m_{h}}=\frac{\hbar^{2}{\bf
K}^{2}}{2(m_{e}+m_{h})}+\frac{\hbar^{2}{\bf
k}^{2}}{2}\left(\frac{1}{m_{e}}+\frac{1}{m_{h}}\right)$ (2)

for ${\bf k}_{e}={\bf k}+\gamma_{e}{\bf K}$ and ${\bf k}_{h}=-{\bf
k}+\gamma_{h}{\bf K}$ with $\gamma_{e}=1-\gamma_{h}=m_{e}/(m_{e}+m_{h})$.
Depending on which hole is scattered, $m_{h}$ can be $m_{{}_{H}}$ or
$m_{{}_{L}}$. So, although the pair center-of-mass wave vector, ${\bf
k}_{e}+{\bf k}_{h}={\bf K}$, stays constant in a Coulomb scattering, it is
impossible to define a pair relative-motion wave vector ${\bf k}$ and follow
its change at each Coulomb interaction.

To approach this problem, we introduce the $\hat{H}_{0}$ Hamiltonian in which
all $\mathcal{J}_{\bf k}=(\pm 3/2,\pm 1/2)$ holes are taken with the same mass
$m_{h}$, not necessarily equal to $m^{*}_{h}$. The resulting ground-state
exciton, obtained from the standard hydrogen-like procedure, is $(2\times
4)$-fold degenerate. The consequences of the hole mass difference follow from
$\hat{\Delta}=\hat{H}-\hat{H}_{0}$ acting in this eightfold exciton subspaceSM
, for $\hat{H}$ being the Hamiltonian when $m_{{}_{H}}\neq m_{{}_{L}}$.

Table 1: Values for $m_{h}^{\ast}$ in Eq. (1), $m_{h}^{\ast\ast}$ in Eq. (3), $E^{(\pm)}_{\bf K}=E_{\bf K}(1\pm r)$ in Eq. (4) and $T_{K}=E_{\bf K}/k_{{}_{B}}$, for various semiconductors. The mass unit is the free electron mass $m_{0}$. The exciton wave vector ${\bf K}$ is obtained from $E_{gap}=\hbar(c/n_{r})|{\bf K}|$ where $n_{r}$ is the semiconductor refractive index.  |  GaAs |  GaSb |  InP |  InAs
---|---|---|---|---
$m_{e}$ |  0.063 |  0.041 |  0.08 |  0.023
$m_{{}_{H}}$ |  0.51 |  0.4 |  0.6 |  0.42
$m_{{}_{L}}$ |  0.082 |  0.05 |  0.089 |  0.026
$E_{gap}$ (eV) |  1.52 |  0.81 |  1.42 |  0.42
$m_{h}^{\ast}$ |  0.14 |  0.089 |  0.16 |  0.049
$m_{h}^{\ast\ast}/m_{h}^{\ast}$ |  1.15 |  1.17 |  1.17 |  1.22
$r$ |  0.5 |  0.53 |  0.49 |  0.6
$n_{r}$ |  3.26 |  3.71 |  3.08 |  3.51
$E_{\bf K}$ ($\mu$eV) |  117 |  68.5 |  79.4 |  28.9
$T_{K}$ (K) |  1.36 |  0.8 |  0.92 |  0.34

Main results –

$\bullet$ For zero exciton center-of-mass wave vector, the $\hat{\Delta}$
difference is diagonal in the eightfold ground-exciton subspace, all states
suffering the same energy shift. So, the hole mass difference does not break
the exciton degeneracy. The effective average hole mass associated with the
ground exciton binding, as obtained from Eq. (25), reads

$m_{h}^{\ast\ast}\simeq
m_{h}^{\ast}\left[1+\beta\frac{m_{e}}{m_{e}+m_{h}^{\ast}}\left(\frac{m_{{}_{H}}-m_{{}_{L}}}{m_{{}_{H}}+m_{{}_{L}}}\right)^{2}\right]$
(3)

with $\beta\simeq 0.90$. It deviates from $m^{\ast}_{h}$ at second order only
in hole mass difference, with a larger deviation for large electron mass.
Table 1 for zinc-blende-like semiconductors, shows that the hole mass
$m_{h}^{\ast\ast}$ is heavier than $m_{h}^{\ast}$ by $15\%$ to $20\%$.

$\bullet$ For nonzero exciton wave vector ${\bf K}$, the $\hat{\Delta}$
difference splits the eightfold ground exciton into two branches in a tricky
way: the resulting heavy and light excitons are not made of heavy or light
holes as naïvely thought, but of their spin-orbit-like combinations,
$\mathcal{J}_{\bf K}=\pm 3/2$ and $\mathcal{J}_{\bf K}=\pm 1/2$, with
$\mathcal{J}_{\bf K}$ quantized along ${\bf K}$. The associated center-of-mass
energies, given in Eq. (29), read (Table 1)

$E_{{\bf K}}^{(\pm)}=\frac{\hbar^{2}{\bf
K}^{2}}{2M_{{}_{X}}^{(\pm)}}\simeq\frac{\hbar^{2}{\bf
K}^{2}}{2(m_{e}+m^{\ast}_{h})}\left[1\pm\frac{m_{{}_{H}}-m_{{}_{L}}}{m_{{}_{H}}+m_{{}_{L}}}\,\,\frac{m^{\ast}_{h}}{m_{e}{+}m^{\ast}_{h}}\right]$
(4)

When the hole mass ratio $m_{{}_{L}}/m_{{}_{H}}$ decreases from 1 to 0, the
center-of-mass ratio $(m_{e}+m^{\ast}_{h})/M_{{}_{X}}^{(\pm)}$ varies from 1
to $1\pm m^{\ast}_{h}/(m_{e}+m^{\ast}_{h})=1\pm
2m_{{}_{L}}/(m_{e}+2m_{{}_{L}})$: the smaller the electron mass compared to
light hole, the larger the change in center of mass, a result hard to
anticipate.

Relevant Hamiltonians – The Hamiltonian $\hat{H}$ for one electron and one
hole reduces to the kinetic energies of the spatially nondegenerate electron
and the threefold hole, plus their Coulomb interaction,
$\hat{H}=\hat{K}^{(e)}+\hat{K}^{(h)}+\hat{V}^{(eh)}$.

$\bullet$ The kinetic part for $m_{e}$ conduction electron reads

$\hat{K}^{(e)}=\sum_{{\bf k}_{e}}\sum_{s=\pm 1/2}\frac{\hbar^{2}{\bf
k}_{e}^{2}}{2m_{e}}\hat{a}^{\dagger}_{{\bf k}_{e},s}\hat{a}_{{\bf k}_{e},s}$
(5)

where $\hat{a}^{\dagger}_{{\bf k}_{e},s}$ creates a ${\bf k}_{e}$ electron
with spin $s$ quantized along an arbitrary axis.

In the spherical approximation that neglects the warping, the upper fourfold
valence level consists of holes having different masses
$(m_{{}_{H}},m_{{}_{L}})$, labeled by $\mathcal{J}_{{\bf k}_{h}}=(\pm 3/2,\pm
1/2)$ index quantized along the hole wave vector ${\bf k}_{h}$. The hole
kinetic part reads

$\hat{K}^{(h)}=\sum_{{\bf k}_{h}}\sum_{\mathcal{J}=(\pm 3/2,\pm
1/2)}\frac{\hbar^{2}{\bf k}_{h}^{2}}{2m_{{}_{\mathcal{J}_{{\bf
k}_{h}}}}}\hat{b}^{\dagger}_{{\bf k}_{h},\mathcal{J}_{{\bf
k}_{h}}}\hat{b}_{{\bf k}_{h},\mathcal{J}_{{\bf k}_{h}}}$ (6)

where $\hat{b}^{\dagger}_{{\bf k}_{h},\mathcal{J}_{{\bf k}_{h}}}$ creates a
heavy or light hole with mass $(m_{\pm 3/2},m_{\pm
1/2})=(m_{{}_{H}},m_{{}_{L}})$.

$\bullet$ The electron-hole Coulomb interaction is diagonal Monicbook  when
written in terms of holes in Bloch states labeled by $\mu=(x,y,z)$, or in
spin-orbit states labeled by $\mathcal{J}_{\textbf{z}}$ quantized along the
crystal axis z

$\displaystyle\hat{V}^{(eh)}$ $\displaystyle=$ $\displaystyle-\sum_{s=\pm
1/2}\,\,\,\sum_{\mathcal{J}=(\pm 3/2,\pm 1/2)}$ $\displaystyle\sum_{{\bf
q}\not=0}\frac{4\pi e^{2}}{\epsilon_{sc}L^{3}q^{2}}\sum_{{\bf k}_{e}{\bf
k}_{h}}\hat{a}^{\dagger}_{{\bf k}_{e}+{\bf
q},s_{\textbf{z}}}\hat{b}^{\dagger}_{{\bf k}_{h}-{\bf
q},\mathcal{J}_{\textbf{z}}}\hat{b}_{{\bf
k}_{h},\mathcal{J}_{\textbf{z}}}\hat{a}_{{\bf k}_{e},s_{\textbf{z}}}$

for a sample volume $L^{3}$ and a dielectric constant $\epsilon_{sc}$.

By contrast, this interaction does not stay diagonal, with complicated Coulomb
scatteringsShiauprb21 , when written in terms of heavy/light holes with
$\mathcal{J}_{{\bf k}_{h}}$ index quantized along the hole wave vector ${\bf
k}_{h}$. Indeed, the creation operators for the two sets of holes are related
by the following basis change

$\hat{b}^{\dagger}_{{\bf k}_{h},\mathcal{J}_{{\bf
k}_{h}}}=\sum_{\mathcal{J}^{\prime}_{\textbf{z}}=(\pm 3/2,\pm
1/2)}\hat{b}^{\dagger}_{{\bf
k}_{h},\mathcal{J}^{\prime}_{\textbf{z}}}\,\,\,\,{}_{\textbf{z}}\langle\mathcal{J}^{\prime}|\mathcal{J}\rangle_{{\bf
k}_{h}}$ (8)

where the ${}_{\textbf{z}}\langle\mathcal{J}^{\prime}|\mathcal{J}\rangle_{{\bf
k}_{h}}$ overlap dependsShiauprb21  on the $(\theta_{{\bf
k}_{h}},\varphi_{{\bf k}_{h}})$ angles of the ${\bf k}_{h}$ vector in the
$(x,y,z)$ crystal axes, $\mathcal{J}_{{\bf k}_{h}}$ reducing to
$\mathcal{J}_{\textbf{z}}$ for ${\bf k}_{h}$ along z.

Exciton basis – To derive the effects of the hole mass difference on exciton,
we introduce as exciton basis, the eigenstates of

$\displaystyle\hat{H}_{0}=\hat{K}^{(e)}+\hat{K}^{(h)}_{0}+\hat{V}^{(eh)}$ (9)

where $\hat{K}^{(h)}_{0}$ reads as $\hat{K}^{(h)}$ with
$(m_{{}_{H}},m_{{}_{L}})$ replaced by $m_{h}$, which is determined by
minimizing the resulting exciton energy in the presence of different hole
masses.

$\bullet$ To calculate the $\hat{H}_{0}$ eigenstates, we first note that

$\sum_{\mathcal{J}=(\pm 3/2,\pm 1/2)}|\mathcal{J}\rangle_{{\bf k}_{h}}{}_{{\bf
k}_{h}}\langle\mathcal{J}|=\sum_{\mathcal{J}=(\pm 3/2,\pm
1/2)}|\mathcal{J}\rangle_{\textbf{z}}{}_{\textbf{z}}\langle\mathcal{J}|$ (10)

since they both form a complete set. This allows us to replace
$\mathcal{J}_{{\bf k}_{h}}$ by $\mathcal{J}_{\textbf{z}}$ in
$\hat{K}^{(h)}_{0}$, since the hole masses in $\hat{K}^{(h)}_{0}$ are the
same. The $\hat{H}_{0}$ Hamiltonian then splits as

$\hat{H}_{0}=\sum_{s_{\textbf{z}}=\pm
1/2}\,\,\,\sum_{\mathcal{J}_{\textbf{z}}=(\pm 3/2,\pm
1/2)}\hat{h}_{s_{\textbf{z}},\mathcal{J}_{\textbf{z}}}$ (11)

where $\hat{h}_{s_{\textbf{z}},\mathcal{J}_{\textbf{z}}}$ is a hydrogen-like
Hamiltonian.

$\bullet$ The resulting eigenenergies split into a center-of-mass part and a
relative-motion part

$E_{{\bf K},\nu;s_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=\frac{\hbar^{2}{\bf
K}^{2}}{2(m_{e}+m_{h})}+\varepsilon_{\nu}$ (12)

The creation operators for the corresponding $(2\times 4)$-fold excitons
$|{\bf K},\nu;s_{\textbf{z}},\mathcal{J}_{\textbf{z}}\rangle$ read as

$\hat{B}^{\dagger}_{{\bf
K},\nu;s_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=\sum_{\bf
k}\hat{a}^{\dagger}_{-{\bf k}+\gamma_{e}{\bf
K},s_{\textbf{z}}}\hat{b}^{\dagger}_{{\bf k}+\gamma_{h}{\bf
K},\mathcal{J}_{\textbf{z}}}\langle{\bf k}|\nu\rangle$ (13)

where $\langle{\bf k}|\nu\rangle$ is the hydrogen-like wave function, with
energy $\varepsilon_{\nu}$ that scales in Rydberg unitLandau

$R_{{}_{X}}=\frac{\mu_{{}_{X}}e^{4}}{2\hbar^{2}\epsilon_{sc}^{2}}=\frac{\hbar^{2}}{2\mu_{{}_{X}}a_{{}_{X}}^{2}}$
(14)

for $\mu_{{}_{X}}^{-1}=m_{e}^{-1}+m_{h}^{-1}$.

These $\hat{H}_{0}$ eigenstates are used as a basis to approach excitons with
different hole masses.

Effects of heavy and light holes – The hole mass difference is concentrated
into the one-body operator

$\displaystyle\hat{\Delta}=\hat{H}-\hat{H}_{0}=\hat{K}^{(h)}-\hat{K}^{(h)}_{0}$
(15)

$\bullet$ To better catch its effect in the $\hat{H}_{0}$ exciton subspace, we
split $\hat{\Delta}$ as $\hat{\Delta}^{(h)}+\hat{\Delta}^{(HL)}$. The part
$\hat{\Delta}^{(h)}=\hat{K}_{0}^{(h^{*})}-\hat{K}_{0}^{(h)}$ for
$\hat{K}_{0}^{(h^{*})}$ reading as $\hat{K}_{0}^{(h)}$ with $m_{h}$ replaced
by $m^{*}_{h}$, concentrates on using $m_{h}$ instead of the commonly accepted
$m^{*}_{h}$ value. The part
$\hat{\Delta}^{(HL)}=\hat{K}^{(h)}-\hat{K}_{0}^{(h^{*})}$, which cancels for
$m_{{}_{H}}=m_{{}_{L}}$, concentrates on the mass difference.

$\bullet$ The $\hat{\Delta}$ matrix elements in the $\nu=\nu_{0}$ ground
exciton subspace

$\langle{\bf
K}^{\prime},\nu_{0};s^{\prime}_{\textbf{z}},\mathcal{J}^{\prime}_{\textbf{z}}|\hat{\Delta}|{\bf
K},\nu_{0};s_{\textbf{z}},\mathcal{J}_{\textbf{z}}\rangle=\delta_{{\bf
K}^{\prime},{\bf
K}}\delta_{s^{\prime}_{\textbf{z}},s_{\textbf{z}}}\Delta_{{\bf
K};\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}$ (16)

differ from zero for $s^{\prime}_{\textbf{z}}=s_{\textbf{z}}$ and ${\bf
K}^{\prime}={\bf K}$ because $\hat{\Delta}$, that comes from difference in
hole kinetic energies, does not act on spin nor on wave vector. This feature
holds true for higher-order terms in $\hat{\Delta}$.

Change in binding energy – We first focus on binding energy, i.e.,
$\textbf{K}=\bf 0$ exciton.

$\bullet$ We find that the first-order contribution in $\hat{\Delta}^{(HL)}$
cancels (Eq. (S39) of SM ),

$\Delta^{(HL)}_{\textbf{0};\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=0$
(17)

while the $\hat{\Delta}^{(h)}$ part leads to (Eq. (S41) of SM )

$\Delta^{(h)}_{\textbf{0};\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=\delta_{\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}\frac{\hbar^{2}}{2a_{{}_{X}}^{2}}\left(\frac{1}{m_{h}^{*}}-\frac{1}{m_{h}}\right)$
(18)

Being diagonal in $\mathcal{J}_{\textbf{z}}$, the $\Delta^{(h)}$ difference
thus produces the same energy shift

$R_{{}_{X}}\mu_{{}_{X}}\left(\frac{1}{\mu^{*}_{{}_{X}}}-\frac{1}{\mu_{{}_{X}}}\right)$
(19)

to all eightfold ground-exciton states, that is, no splitting. The resulting
exciton energy then reads

$R_{{}_{X}}\left[-1+\mu_{{}_{X}}\left(\frac{1}{\mu^{*}_{{}_{X}}}-\frac{1}{\mu_{{}_{X}}}\right)\right]=R_{{}_{X}}^{*}\frac{\mu_{{}_{X}}}{\mu^{*}_{{}_{X}}}\left[-2+\frac{\mu_{{}_{X}}}{\mu^{*}_{{}_{X}}}\right]$
(20)

where $R_{{}_{X}}^{*}$ is given by Eq. (14) for $m_{h}=m_{h}^{*}$. Its minimum
value, obtained for $\mu_{{}_{X}}=\mu^{*}_{{}_{X}}$, gives a maximum ground-
exciton binding equal to $R_{{}_{X}}^{*}$. So, up to first order in
$\hat{\Delta}$, the hole mass difference does not break the ground-exciton
degeneracy, the appropriate hole mass value for binding being the one given in
Eq. (1).

$\bullet$ Similar but somewhat heavier calculationsSM  performed for the
product $\hat{\Delta}^{(h)}\,\hat{\Delta}^{(h)}$, show that the
$(m_{h}-m^{*}_{h})$ difference leads to (Eq. (S65) of SM )

$\Delta^{(h,h)}_{{\bf
0},\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=\delta_{\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}\left[\frac{\hbar^{2}}{2a_{{}_{X}}^{2}}\left(\frac{1}{m_{h}^{*}}-\frac{1}{m_{h}}\right)\right]^{2}\frac{\beta_{{{}_{1}}}}{R_{{}_{X}}}$
(21)

with $\beta_{{}_{1}}\simeq-0.50$. In contrast to first order, the
$(m_{{}_{H}}-m_{{}_{L}})$ difference now brings a nonzero contribution through
$\hat{\Delta}^{(HL)}\,\hat{\Delta}^{(HL)}$, equal to (Eq. (S56) of SM )

$\Delta^{(HL,HL)}_{{\bf
0},\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=\delta_{\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}\left[\frac{\hbar^{2}}{2a_{{}_{X}}^{2}}\left(\frac{1}{m_{{}_{H}}}-\frac{1}{m_{{}_{L}}}\right)\right]^{2}\frac{\beta_{{{}_{2}}}}{R_{{}_{X}}}$
(22)

with $\beta_{{}_{2}}\simeq-0.22$, while the mixed term in
$\hat{\Delta}^{(HL)\,}\hat{\Delta}^{(h)}$ gives no contribution.

$\bullet$ Since the above contributions also are diagonal in
$\mathcal{J}_{\textbf{z}}$, the exciton degeneracy remains unbroken, but the
energy shift they produce now depends on the hole mass difference. When
combined with Eq. (20), the exciton energy now reads as

$\displaystyle
R_{{}_{X}}\left[-1+\mu_{{}_{X}}\left(\frac{1}{\mu_{{}_{X}}^{*}}{-}\frac{1}{\mu_{{}_{X}}}\right)+\mu_{{}_{X}}^{2}\left(\frac{1}{\mu_{{}_{X}}^{*}}{-}\frac{1}{\mu_{{}_{X}}}\right)^{2}\beta_{1}\right.$
$\displaystyle\left.+\mu_{{}_{X}}^{2}\left(\frac{1}{m_{{}_{H}}}{-}\frac{1}{m_{{}_{L}}}\right)^{2}\beta_{2}\right]\hskip
28.45274pt$ (23)

Its minimum value, obtained for

$\mu_{{}_{X}}=\mu_{{}_{X}}^{*}\left(1-\frac{3\beta_{2}}{2+2\beta_{1}}\mu_{{}_{X}}^{\ast\,2}\left(\frac{1}{m_{{}_{H}}}-\frac{1}{m_{{}_{L}}}\right)^{2}\right)$
(24)

gives the ground exciton energy as (Eq. (S71) of SM )

$R_{{}_{X}}^{*}\left[-1+\beta_{2}\mu_{{}_{X}}^{\ast\,2}\left(\frac{1}{m_{{}_{H}}}-\frac{1}{m_{{}_{L}}}\right)^{2}\right]\equiv-\frac{\mu_{{}_{X}}^{\ast\ast}e^{4}}{2\hbar^{2}\epsilon_{sc}^{2}}$
(25)

that corresponds, for
$\mu_{{}_{X}}^{\ast\ast\,-1}=m_{e}^{-1}+m_{h}^{\ast\ast\,-1}$, to the hole
effective mass $m_{h}^{\ast\ast}$ given in Eq. (3).

Change in center-of-mass energy – We now look for the effect of the
$\hat{\Delta}$ difference on the eightfold ground exciton subspace when its
center-of-mass wave vector ${\bf K}$ differs from zero. This will tell how the
hole mass difference affects the exciton motion. By taking $m_{h}=m_{h}^{*}$,
which produces the correct binding at first order in $\hat{\Delta}$, this
difference reduces to $\hat{\Delta}^{(HL)}$. Its contribution to the exciton
energy appears as (Eq. (S73) of SM )

$\Delta^{(HL)}_{{\bf
K};\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}=\sum_{{\bf
k}}|\langle{\bf k}-\gamma^{\ast}_{h}{\bf
K}|\nu_{0}\rangle|^{2}\frac{\hbar^{2}{\bf
k}^{2}}{2m_{{}_{HL}}}\,{}_{\textbf{z}}\langle\mathcal{J}^{\prime}|\hat{D}_{{\bf
k}}|\mathcal{J}\rangle_{\textbf{z}}$ (26)

with $\gamma^{\ast}_{h}=m_{h}^{\ast}/(m_{e}+m_{h}^{\ast})$ and
$1/m_{{}_{HL}}=(1/m_{{}_{H}}-1/m_{{}_{L}})/2$, for $\hat{D}_{\bf k}$ given by

$\displaystyle\hat{D}_{\bf k}=\Big{(}\sum_{\mathcal{J}=\pm
3/2}-\sum_{\mathcal{J}=\pm 1/2}\Big{)}|\mathcal{J}\rangle_{\bf k}\,{}_{\bf
k}\langle\mathcal{J}|$ (27)

As $\Delta^{(HL)}_{{\bf
K};\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}$ is equal to
zero for ${\bf K}={\bf 0}$, its lowest-order term in ${\bf K}$, obtained from
the $|\langle{\bf k}-\gamma^{\ast}_{h}{\bf K}|\nu_{0}\rangle|^{2}$ expansion,
leads after summing over ${\bf k}$, to

$\Delta^{(HL)}_{{\bf
K};\mathcal{J}^{\prime}_{\textbf{z}},\mathcal{J}_{\textbf{z}}}\simeq\frac{\hbar^{2}\gamma^{\ast\,2}_{h}{\bf
K}^{2}}{2m_{{}_{HL}}}\,{}_{\textbf{z}}\langle\mathcal{J}^{\prime}|\hat{D}_{{\bf
K}}|\mathcal{J}\rangle_{\textbf{z}}$ (28)

To go further, we note that changing the hole quantization axis from z to
${\bf K}$, renders diagonal the above matrix elements. This leads to (Eq.
(S80) of SM )

$\Delta^{(HL)}_{{\bf K};\mathcal{J}^{\prime}_{\bf K},\mathcal{J}_{\bf
K}}=\delta_{\mathcal{J}^{\prime}_{\bf K},\mathcal{J}_{\bf
K}}\frac{\hbar^{2}\gamma_{h}^{\ast\,2}{\bf
K}^{2}}{2m_{{}_{HL}}}\Big{(}\delta_{\mathcal{J}_{\bf
K},\pm\frac{3}{2}}{-}\delta_{\mathcal{J}_{\bf K},\pm\frac{1}{2}}\Big{)}$ (29)

For different hole masses, the exciton center-of-mass energy thus splits into
a heavy branch, $\mathcal{J}_{\bf K}=\pm 3/2$, and a light branch,
$\mathcal{J}_{\bf K}=\pm 1/2$, as said in Eq. (4), the energy difference being
more significant for light electron, $m_{e}\ll m^{\ast}_{h}$. Note that these
heavy and light exciton masses differs from the naïve $(m_{e}+m_{{}_{H}})$ and
$(m_{e}+m_{{}_{L}})$ values the excitons would have if they were simply made
of heavy or light hole.

Experimental consequences –

$\bullet$ The most direct way to evidence this exciton splitting would be to
see two distinct absorption lines associated with the two exciton branches.
This appears difficult in the case of GaAs because the exciton linewidth in
usual samples ($\simeq 320\,\mu$eV in Gopal2000 ) is large compared with the
$E_{{\bf K}}^{(\pm)}$ difference $(\simeq 117\,\mu$eV in Table 1), the
photocreated exciton having a wave vector equal to the photon wave vector; so,
$E_{gap}\simeq\hbar(c/n_{r})|{\bf K}|$ where $n_{r}$ is the semiconductor
refractive index. As a result, the heavy and light exciton lines would appear
as a single broad line.

$\bullet$ Another idea is to note that due to their energy differences, the
thermal populations of heavy and light excitons evolve differently

$N_{{\bf K},T}^{(\pm)}=N_{{\bf K},\infty}\exp(-E_{{\bf K}}(1\pm
r)/k_{{}_{B}}T)\propto\exp(\pm rT_{K}/T)$ (30)

for $T_{K}=E_{{\bf K}}/k_{{}_{B}}\simeq 1.4\,$K. Since heavy excitons have a
lower energy, they contribute to the low-energy side of the absorption line,
with a number that is not only larger than the light exciton number but that
also increases faster when $T$ decreases; so, the shape of the absorption line
should show an enhancement on the low-energy side.

Another way to quantify the difference in the heavy and light exciton
populations, is through the total weight of the photoluminescence line as a
function of temperature: being an overlap of two lines with amplitudes that
vary differently, the $\exp(2rT_{K}/T)$ difference must show up in a Log plot
of the temperature evolution in the subKelvin temperature range.

$\bullet$ Such temperature dependences can be directly compared to the
predicted exciton splitting driven by heavy-light hole Coulomb couplings,
provided that there is no extrinsic effect associated with residual crystal
strain1 ; 3 ; 5 . Regarding this point, it has been shown7  that epitaxial
GaAs samples on a Si substrate are essentially strain-free beyond a critical
thickness $\simeq 2.3\,\mu$m.

Another major problem is that, for such low temperatures, bright excitons are
going to mostly condense into dark statesBEC ; dubin2017 , making their
observation quite stringent.

Conclusion – We here study the consequences of valence holes having two
different masses, ($m_{{}_{H}},m_{{}_{L}}$), on the exciton energy, due to
their nondiagonal intraband Coulomb couplingearlywork .

We show that for zero center-of-mass wave vector, the ground-state exciton
stays $(2\times 4)$-fold degenerate, its binding energy reading in terms of a
single hole mass that is close to the commonly-used average value
$m_{h}^{\ast}$ given in Eq. (1). We expect this conclusion, mathematically
proved up to second order in hole mass difference, to stay valid up to any
order for a fundamental reason: in the absence of exciton wave vector, there
is no symmetry-breaking axis to possibly split the exciton degeneracy. This
strong argument leads us to anticipate that, when included, the valence band
warping will not split the exciton energywarping .

By contrast, a nonzero center-of-mass wave vector K splits the exciton
degeneracy into a heavy branch and a light branch, with $(\pm 3/2)$ and $(\pm
1/2)$ indices quantized along the symmetry-breaking axis K. The associated
center-of-mass masses depend on the ($m_{{}_{H}}-m_{{}_{L}}$) difference and
the electron mass $m_{e}$ as given in Eq. (4).

This splitting raises fundamental questions that deserve further investigation
on using a unique hole mass for many-body effects in which excitons with
finite wave vector play a key role, like the exciton Bose-Einstein
condensationBEC ; BEC2 ; BEC3 , the two-component condensatesKuklov ; Cabrera
, and the BEC-BCS crossoverBEC-BCS ; Li .

## References

  * (1) M. Kira, and S. W. Koch, Semiconductor Quantum Optics, Cambridge University Press (Cambridge, 2012).
  * (2) M. Feierabend, G. Berghäuser, A. Knorr, and E. Malic, Nat. Commun. 8, 14776 (2017).
  * (3) C. Lagoin, U. Bhattacharya, T. Grass, R. W. Chhajlany, T. Salamon, K. Baldwin, L. Pfeiffer, M. Lewenstein, M. Holzmann, and F. Dubin, Nature 609 485 (2022).
  * (4) H. J. Krenner, C. E. Pryor, J. He, and P. M. Petroff, Nano Letters 8, 1750 (2008).
  * (5) K. Cong, W.-W. Jiang, B. E. Anthonio, G. T. Noe, H.-P. Liu, H. Kataura, M. Kira, and J. Kono, ACS Photonics 7, 1382 (2020).
  * (6) D. Unuchek, A. Ciarrocchi, A. Avsar, Z. Sun, K. Watanabe, T. Taniguchi, and A. Kis, Nat. Nanotechnol. 14, 1104 (2019).
  * (7) S.-Y. Shiau, and M. Combescot, Phys. Rev. B 104, 045203 (2021).
  * (8) C. M. Bender, A. Duncan, and H. F. Jones, Phys. Rev. D 49, 4219 (1994).
  * (9) R. Guida, K. Konishi, and H. Suzuki, Ann. Phys. 241, 152 (1995).
  * (10) P. Y. Yu, and M. Cardona, Fundamentals of Semiconductors, 3rd ed. (Springer, Berlin, 2005).
  * (11) G. Fishman, Semiconducteurs: Les Bases de la Théorie ${\bf k}\cdot{\bf p}$ (Ecole Polytechnique, Paris, 2010).
  * (12) G. Dresselhaus, A. F. Kip, and C. Kittel, Phys. Rev. 98, 368 (1955).
  * (13) G. Dresselhaus, J. Phys. Chem. Solids, 1, 14 (1956).
  * (14) J. M. Luttinger, Phys. Rev. 102, 1030 (1956).
  * (15) A. Baldereschi, and N. O. Lipari, Phys. Rev. Lett. 25, 373 (1970).
  * (16) A. Baldereschi, and N. O. Lipari, Phys. Rev. B 3, 439 (1971).
  * (17) M. Combescot and S.-Y. Shiau, Excitons and Cooper Pairs, Oxford University Press (Oxford, 2015).
  * (18) Light holes also play a minor role for dense electron-hole plasmas because their number is very small (see M. Combescot, Phys. Rev. Lett. 32, 15 (1974)).
  * (19) See the Supplemental Information for detailed calculations of the first- and second-order contributions to the exciton binding energy and center-of-mass energy, and the $(\beta_{{}_{1}},\beta_{{}_{2}})$ values, which includes Refs. Shiauprb21  and Landau .
  * (20) L. D. Landau, and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory, Vol. 3, 3rd ed. (Pergamon Press, 1977).
  * (21) A. V. Gopal, R. Kumar, A. S. Vengurlekar, A. Bosacchi, S. Franchi, and L. N. Pfeiffer, J. Appl. Phys. 87, 1858 (2000).
  * (22) M. Joschko, M. Woerner, T. Elsaesser, E. Binder, T. Kuhn, R. Hey, H. Kostial, and K. Ploog, Phys. Rev. Lett. 78, 737 (1997).
  * (23) N. H. Bonadeo, D. G. Steel, and R. Merlin, Phys. Rev. B 60, 8970 (1999).
  * (24) K. Wundke, U. Neukirch, J. Gutowski, and D. Hommel, Phys. Rev. B 53, 10973 (1996).
  * (25) Y. Huang, P. Y. Yu, M.-N. Charasse, Y. Lo, and S. Wang, Appl. Phys. Lett. 51, 192 (1987).
  * (26) R. Anankine, M. Beian, S. Dang, M. Alloing, E. Cambril, K. Merghem, C. G. Carbonell, A. Lemaître, and F. Dubin, Phys. Rev. Lett. 118, 127402 (2017).
  * (27) Early worksBaldereschi_prl ; Baldereschi_prb  that mainly focus on exciton binding corrections do not address these consequences because they stay with the effective-mass Hamiltonian, instead of the spherical Hamiltonian.
  * (28) The effect of the hole warping on the exciton binding and center-of-mass motion will be considered elsewhere.
  * (29) M. Alloing, M. Beian, M. Lewenstein, D. Fuster, Y. Gonzalez, L. Gonzalez, R. Combescot, M. Combescot, and F. Dubin, Europhys. Lett. 107, 10012 (2014).
  * (30) M. Combescot, R. Combescot, and F. Dubin, Rep. Prog. Phys. 80, 066501 (2017).
  * (31) Y. Morita, K. Yoshioka, and M. Kuwata-Gonokami, Nat. Commun. 13, 5388 (2022)
  * (32) A. Kuklov, N. Prokof’ev, and B. Svistunov, Phys. Rev. Lett. 92, 030403 (2004).
  * (33) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, Science 359, 301 (2017).
  * (34) Y. Nakagawa, Y. Kasahara, T. Nomoto, R. Arita, T. Nojima, and Y. Iwasa, Science 372, 190 (2021).
  * (35) J. I. A. Li, T. Taniguchi, K. Watanabe, J. Hone, and C. R. Dean, Nature Phys. 13, 751 (2017).


# Variance Minimisation of the Lipkin-Meshkov-Glick Model
on a Quantum Computer

I. Hobday  School of Mathematics and Physics, University of Surrey,
Guildford, Surrey, GU2 7XH, United Kingdom  P. D. Stevenson  School of
Mathematics and Physics, University of Surrey,
Guildford, Surrey, GU2 7XH, United Kingdom  AWE, Aldermaston, Berkshire, RG7
4PR, United Kingdom  J. Benstead  AWE, Aldermaston, Berkshire, RG7 4PR, United
Kingdom  School of Mathematics and Physics, University of Surrey,
Guildford, Surrey, GU2 7XH, United Kingdom

###### Abstract

Quantum computing can potentially provide advantages for specific
computational tasks. The simulation of fermionic systems is one such task that
lends itself well to quantum computation, with applications in nuclear physics
and electronic systems. Here we present work in which we use a variance
minimisation method to find the full spectrum of energy eigenvalues of the
Lipkin-Meshkov-Glick model; an exactly-solvable nuclear shell model-type
system.

We perform these calculations using both quantum simulators and real quantum
hardware accessed via IBM cloud-based quantum computers. Using these IBM
quantum computers we are able to obtain all eigenvalues for the cases of three
and seven fermions (nucleons) in the Lipkin-Meshkov-Glick model.

##  I Introduction

The use of quantum algorithms to perform calculations in nuclear physics is a
rapidly developing field [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21]. The exponential scaling of Hilbert space with the
number of quantum bits (qubits) and the ability of multiple qubits to exhibit
highly entangled wave functions give quantum computers the potential to have a
great impact in simulating many-body quantum systems. In particular, one can
expect a quantum advantage where quantum computation outperforms classical
computation as the size of the system under study becomes sufficiently large.

Richard Feynman first proposed using quantum computers to study quantum
systems in the 1980s [22]. Since then, various algorithms have been developed
that are able to perform calculations on many-body quantum systems, such as
Quantum Phase Estimation (QPE) [23, 24, 25], and Quantum Imaginary Time
Evolution [26, 27, 28]. However, these algorithms are often too complex to be
fully implemented on current quantum computers.

Current quantum computers are said to be in the Noisy, Intermediate-Scale
Quantum (NISQ) era. This is due to the low numbers of qubits, and large
amounts of noise and error present within the devices. Current devices also
have low coherence times, restricting the amount of time a qubit can maintain
its state. This limits the length of the quantum circuit, or circuit depth,
that quantum computers can use and still produce meaningful results [29].

###  I.1 Variational Algorithms

Variational quantum algorithms have established themselves due to the relative
ease of running them on NISQ-era quantum computers [30, 31, 32]. These are
hybrid algorithms, meaning they use both quantum and classical processes in
tandem, reducing the amount of computational work needed to be performed on
the quantum computer. Working in conjunction with a classical computer allows
for reduced circuit depths and improved error rates by reducing the number of
gates in the quantum circuit. This reduction in error allows for meaningful
results to be obtained on current quantum hardware.

The Variational Quantum Eigensolver (VQE) algorithm uses the variational
principle of quantum mechanics to approximate the ground state of some
Hamiltonian. In the work presented here, we apply a version of the VQE which
targets any eigenvalue of a quantum system by minimizing the variance of the
Hamiltonian rather than the expectation value. Application of the method is
made to the simplified nuclear model due to Lipkin, Meshkov, and Glick, which
we summarize in the next section. The presentation proceeds with a description
of our implementation of the algorithm, followed by results on simulated and
real quantum computers.

###  I.2 Lipkin-Meshkov-Glick Model

The Lipkin-Meshkov-Glick (LMG) model is an exactly solvable nuclear model
introduced in the 1960s [33]. It is a simple shell model consisting of two
levels, separated by an energy $\epsilon$ with a permutation-symmetric
potential. The two levels are each $N$-fold degenerate and a standard
treatment (which we adopt here) considers $N$ fermions in the system. The
Hamiltonian for the model is given as

$\displaystyle H=$ $\displaystyle\frac{1}{2}\epsilon\sum_{p\sigma}\sigma
a_{p\sigma}^{\dagger}a_{p\sigma}+\frac{1}{2}V\sum_{pp^{\prime}\sigma}a_{p,\sigma}^{\dagger}a_{p^{\prime},\sigma}^{\dagger}a_{p^{\prime},-\sigma}a_{p,-\sigma}$
(1)
$\displaystyle+\frac{1}{2}W\sum_{pp^{\prime}\sigma}a_{p,\sigma}^{\dagger}a_{p^{\prime},-\sigma}^{\dagger}a_{p^{\prime},\sigma}a_{p,-\sigma}.$

Here the labels $p$ and $p^{\prime}$ run over the $N$ degenerate states in
each level, while $\sigma=\pm 1$ labels the level. The strength $V$ controls
pair (de-)excitations between the two levels, while the $W$ term scatters one
particle up and another down. The model can be simplified by writing the
Hamiltonian in the quasi-spin basis, as

$H=\epsilon
J_{z}+\frac{1}{2}V(J_{+}^{2}+J_{-}^{2})+\frac{1}{2}W(J_{+}J_{-}+J_{-}J_{+}),$
(2)

where $J_{z}$ and $J_{\pm}$ are given by

$J_{z}=\frac{1}{2}\sum_{p\sigma}\sigma a_{p\sigma}^{\dagger}a_{p\sigma}$ (3)

and

$J_{\pm}=\sum_{p}a_{p\pm 1}^{\dagger}a_{p\mp 1}$ (4)

respectively.

These quasi-spin operators satisfy angular momentum commutation relations. The
total quasi-spin operator $J^{2}=\frac{1}{2}(J_{+}J_{-}+J_{-}J_{+})+J_{z}^{2}$
commutes with the Hamiltonian, so that each $J$ value may be considered
separately. This symmetry allows one to reduce the maximum Hamiltonian matrix
dimension from $2^{N}$ to $N+1$. Further simplification is possible since the
Hamiltonian also exhibits good “parity” [34], which splits each matrix of a
given $J$ into two block-diagonal parts.

The LMG model is applicable as a test-platform for quantum computing
algorithms applied to nuclear physics due to its exactly solvable nature, and
the ability to scale the model indefinitely by increasing the number of
particles in the system. The model has already been studied using real and
simulated quantum computers, particularly for calculation of its ground state
[13, 16], but also for excited states using a quantum-assisted algorithm to
prepare a generalised eigenvalue problem for solution on a classical computer
for the excited state spectrum [17]. For our work, we use an LMG model
Hamiltonian with the parameters $\frac{V}{\epsilon}=0.5$,
$\frac{W}{\epsilon}=0$, making the common choice to ignore the $W$ term which
does not affect ground state correlations.

##  II Encoding Methodology

The LMG-model Hamiltonian must be encoded in such a way as to allow it to be
processed by a quantum computer. For circuit-based quantum computers, this
form is conveniently expressed as linear combinations of Pauli spin matrices.
Often a representation of a Hamiltonian in second-quantized form, such as Eq.
(1), leads to a mapping from creation and annihilation operators to Pauli
strings, keeping the second quantized representation where solutions of any
particle number can be encoded. Such encodings as Jordan-Wigner [35] or
Bravyi-Kitaev [36, 37] fall in this category. Having made use of the quasispin
representation, and specialising to the standard Lipkin Model in which the
number of particles is fixed at $N$ \- the same as the degeneracy of each of
the two unperturbed levels - we make a more specialised and efficient encoding
starting from a representation of the Hamiltonian in the quasispin basis [13]
in which only fixed-particle-number states are considered.

For a given specific LMG Hamiltonian, $H_{N}$ with fixed $N$ we make a matrix
representation in the quasispin basis and consider (as usual) the maximum
quasispin value, that being the one containing the unperturbed ground state.
The matrix dimension is $N+1$, but owing to the parity symmetry, it can be
written as two block diagonal parts which can be treated independently. When
$N$ is odd, the two submatrices are of equal dimension, $(N+1)/2$. In the
present work we consider $N=3,7$ so that $(N+1)/2$ is of the form $2^{m}$ and
each matrix whose eigenvalues are sought has a power-of-two dimension. For a
Hamiltonian $H_{N}$ for which $2^{m}=(N+1)/2$ we can decompose in terms of
Pauli matrices as

$H_{N}=\sum_{i=0}^{2^{2m}-1}\beta_{i}\bigotimes_{q=1}^{m}\sigma_{{i_{4}}[q]}$
(5)

where $\beta_{i}$ is a coefficient and the tensor product $\bigotimes$
produces all Pauli strings from the set
$\\{\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3}\\}=\\{I,X,Y,Z\\}$ and the
notation $\sigma_{{i_{4}}[q]}$ means the $q^{th}$ digit from the right of $i$
when represented in base 4.

For the given numerical Hamiltonians with particular values for the
interaction strengths $V/\epsilon$ and $W/\epsilon$, the coefficients
$\beta_{i}$ were found using a custom code [38], whose functionality is also
found in common quantum computing libraries (e.g. the pauli_decompose method
of PennyLane [39]). Specific examples of the encodings are found in subsequent
sections.

##  III Finding Excited States

In a typical Variational Quantum Algorithm (VQA), an ansatz is made for a wave
function representing the ground state of the problem at hand. The ansatz
contains parameters which can be adjusted to optimise the wave function and
bring it as close as possible to the true ground state. Usually the
expectation value of the Hamiltonian, expressed in terms of Pauli operators,
is evaluated on the quantum computer, while a classical algorithm manages the
parameter optimisation.

Here, we adapt the typical VQA to seek the minimum of the variance of the
Hamiltonian,

$\sigma^{2}=\langle H^{2}\rangle-\langle H\rangle^{2}.$ (6)

The variance is a positive semi-definite function which is zero when the wave
function is an eigenfunction of the Hamiltonian. Hence, its set of equally
deep global minima will correspond to the set of eigenstates of the
Hamiltonian, except possibly for accidental zeros of the variance which can be
checked for.

This method allows the use of the same quantum circuit ansatzes as when using
the VQE to find the ground state of $H$ directly, but requires additional
circuit measurements to be performed for the terms in $\langle H^{2}\rangle$.
This increases the time taken to perform calculations but does not introduce
any additional circuit depth or variational parameters to the calculations
compared with an energy-minimizing ground state VQE. Since no additional
circuit complexity is introduced, using the variance to find excited states of
the nuclear system is possible on NISQ-era quantum hardware for problems where
the standard VQE is applicable.

In our implementation, we make use of the IBM qiskit environment to perform
our VQE on real quantum hardware, as noted later where individual results are
presented, using qiskit’s classical COBYLA solver for the ansatz parameter
optimisation.

##  IV  $N=3$ Simulation and Measurement

###  IV.1 $N=3$ Hamiltonian

The Hamiltonian for the $N=3$ LMG model is split into its two parity
submatrices labelled as $A$ and $B$. For our chosen values of $V/\epsilon=0.5$
and $W/\epsilon=0$ the submatrices can be expressed as

$H_{N=3A}=\begin{bmatrix}-1.5&-0.866\\\ -0.866&0.5\end{bmatrix}$ (7)

and

$H_{N=3B}=\begin{bmatrix}-0.5&-0.866\\\ -0.866&1.5\end{bmatrix}.$ (8)

When decomposed into the form of Pauli matrices, the Hamiltonians can be
represented as

$H_{N=3A}=-0.5-1.0Z_{0}-0.8660254X_{0}$ (9)

and

$H_{N=3B}=0.5-1.0Z_{0}-0.8660254X_{0}$ (10)

respectively, where we have kept more significant figures than given in the
matrix representation.

We note that the two submatrices differ only by the signs of all the diagonal
elements, while the off-diagonal elements are indetical. This means that the
eigenspectra of the two submatrices differ only by their sign. This is a
general result (for odd $N$) and one needs only compute the spectrum of one of
the submatrices in order to know the spectrum for both, but we test the method
in the N=3 case for both submatrices.

For the evaluation of the variance, Eq. (6), the square of the Hamiltonian
operator is needed. From the square of the matrix representations

$H^{2}_{N=3A}=\begin{bmatrix}1.0&0.866\\\ 0.866&3.0\end{bmatrix}$ (11)

and

$H_{N=3B}^{2}=\begin{bmatrix}1.0&-0.866\\\ -0.866&3.0\end{bmatrix}$ (12)

their Pauli representations are found as

$H^{2}_{N=3A}=2.0+0.8660254X_{0}+1.0Z_{0}$ (13)

and

$H^{2}_{N=3B}=2.0-0.8660254X_{0}-1.0Z_{0}$ (14)

###  IV.2 Circuit Ansatz

From the submatrices presented in Eqs. (7) and (8), or their Pauli form in
Eqs. (9) and (10), it is clear that the wavefunctions of these submatrices can
be expressed using a single qubit ansatz. Thus we use a one-parameter circuit,
shown in Fig. 2, to cover the Hilbert space to find both eigenvalues
associated with the respective submatrix being measured.

###  IV.3 Simulation

We perform variance minimisations, starting from random initial parameters,
using noiseless simulations. The results of these simulations are plotted in
Fig. 1 for the submatrix $H_{N=3A}$, which show the energy values converging
to their respective exact eigenvalues.

Figure 1: Iterative variance minimisation to converge on energy eigenstates
(dashed lines) of the $H_{N=3A}$ sector of the $N=3$ LMG model with parameters
$\frac{V}{\epsilon}=0.5$, $\frac{W}{\epsilon}=0$, using a simulator.

###  IV.4 Quantum Hardware Results

We perform the calculations of the $N=3$ LMG excited state spectrum using a
parameter sweep, since this is feasible for a single parameter and will help
visualise the results. For this, we use IBM’s cloud-based quantum computer,
IMBQ_manila. This is a 5-qubit quantum computer using superconducting circuits
and a “nearest neighbour” entanglement scheme (though for this $N=3$ case we
use only a single qubit). We perform the calculations using 20,000
measurements per circuit (shots), and one circuit for each term in the
Hamiltonians. We also introduce additional circuits to mitigate some of the
measurement readout bias that is innate to the quantum computer. These
additional circuits prepare a state that represents each possible output of
the quantum computer, using Pauli X gates. For a single qubit, this method of
mitigation adds two additional circuits, one measuring the state $|0\rangle$
and the other measuring the state $|1\rangle$, per iteration step of the VQE.
This allows for up-to-date corrections to be applied to our counts [40]. We
also perform VQE variance minimisations, starting from random initial
parameters using IBMQ_manila as described in section III.

We present the results of the calculations in Fig. 3, which shows the
parameter sweep, in Fig. 4 which shows the variance minimisations, and in
Table 1, which shows the resulting states for which the Hamiltonian variance
is zero. These results obtained on a quantum computer are close to the exact
values for the respective eigenstates, within error.

0 & R_y(θ)

Figure 2: Single-qubit circuit ansatz. Figure 3: Parameter sweep of the block-diagonal submatrices, $A$ and $B$, for the $N=3$ LMG model, using 50 iteration steps, using quantum hardware. Figure 4: Iterative variance minimisation to converge on energy eigenstates (dashed lines) of the $H_{N=3A}$ sector of the $N=3$ LMG model with parameters $\frac{V}{\epsilon}=0.5$, $\frac{W}{\epsilon}=0$, using quantum hardware. Random starting parameters are used, and two illustrative results are shown, resulting in the full spectrum. Table 1: $N=3$ LMG model block-diagonal submatrix results calculated using the IBMQ_manila quantum computer using the circuit in Fig. 2 using 20000 shots. The Ground and second excited state results are achieved using submatrix A, with submatrix B used for the first and third excited states.  Eigenstate | Exact Value ($\epsilon$) | Variance | QC Result ($\epsilon$)
---|---|---|---
Ground | $-1.823$ | 0.073 | -1.788 $\pm$ 0.062
2nd | $0.823$ | 0.001 | 0.826 $\pm$ 0.064
1st | $-0.823$ | -0.004 | -0.816 $\pm$ 0.063
3rd | $1.823$ | 0.001 | 1.810 $\pm$ 0.063

##  V  $N=7$ Simulation and Measurement

###  V.1 Ansatz

The $N=7$ submatrix featuring the lowest unperturbed eigenvalue in the
diagonal for our choice of $V/\epsilon=0.5$, $W=0$ is

$H_{N=7\text{ Submatrix}}=\begin{bmatrix}-3.5&-2.291&0&0\\\
-2.291&-1.5&-3.873&0\\\ 0&-3.873&0.5&-3.354\\\ 0&0&-3.354&2.5\\\
\end{bmatrix}.$ (15)

When encoded into the form of Pauli spin matrices, this becomes

$\displaystyle H_{N=7}=$
$\displaystyle-0.5-2.8225X_{1}-1.0Z_{1}-1.9365X_{0}X_{1}$
$\displaystyle-1.9365Y_{0}Y_{1}-2.0Z_{0}+0.5315Z_{0}X_{1}.$ (16)

As with the $N=3$ case the eigenvalues of the opposite parity submatrix can be
obtained from those of the above matrix through multiplication by $-1$.

From Eq. (V.1), a two-qubit circuit is required to represent the wave function
of the system. A sufficiently general circuit for this case, with three
parameters, is shown in Fig. 5.

0 & R_y(θ_1)
0 R_y(θ_0) -1 -1 R_y(θ_2)

Figure 5: Quantum circuit to calculate the eigenvalues of LMG model $N=7$
Hamiltonian.

For the calculation of the variance, the square of $H_{N=7}$ is needed. Matrix
multiplication gives

$H^{2}_{N=7}=\begin{bmatrix}17.498681&11.455&8.873043&0.\\\
11.455&22.49881&3.873&12.990042\\\ 8.873043&3.873&26.499445&-10.062\\\
0.&12.990042&-10.062&17.499316\\\ \end{bmatrix},$ (17)

from which Pauli decomposition produces the form needed for quantum
computation as

$\displaystyle H^{2}_{N=7}=$ $\displaystyle
20.999063+0.6965X_{1}+1.0Z_{1}+10.9315425X_{0}$
$\displaystyle+1.9365X_{0}X_{1}-2.0584995X_{0}Z_{1}+1.93654Y_{0}Y_{1}$
$\displaystyle-1.0003175Z_{0}+10.7585Z_{0}X_{1}-3.5Z_{0}Z_{1}.$ (18)

###  V.2 Simulation

Similarly to the $N=3$ simulation in Section IV.3, we perform noiseless
simulations of the variance minimisation to find each eigenvalue. The results
of these simulations for the $N=7$ LMG model are plotted in Fig. 6, which
shows convergence across each of the four simulations to their respective
exact eigenstate. These simulations are also started from randomised initial
parameters.

Figure 6: Energy values during iteration of the simulated variance
minimisation method for the $N=7$ LMG model with parameters
$\frac{V}{\epsilon}=0.5$, $\frac{W}{\epsilon}=0$. Known exact eigenvalues are
shown with dashed lines. Four random starting points are shown corresponding
to runs which find the four eigenstates.

###  V.3 Quantum Hardware Results

Measurements for the $N=7$ LMG excitation spectra were performed using the
IBM_nairobi quantum computer, again using 20,000 shots per circuit. We apply
the same method of readout error mitigation as with the $N=3$ case, described
in Section IV.4, as well as applying a method of CNOT mitigation. This
mitigation technique adds pairs of CNOT gates to the circuit where there is a
CNOT present, allowing for a linear extrapolation to the theoretical case
where there are no CNOT gates in the circuit. Adding additional pairs of CNOT
gates does not alter the value of the measurement, other than introducing
additional error from the added CNOT pairs [41]. The results of the iterative
variance minimizations on IBM_nairobi are presented in Fig. 7 and in Table 2.
The results are close to the known, exact solutions. However, these results
are further from their respective exact solutions than the $N=3$ results. This
is likely due to the increase in circuit depth required to perform the
calculations for the $N=7$ case, as well as the additional qubit, resulting in
greater gate error and noise within the circuit.

Figure 7: Energy values during iteration of the variance minimisation method on quantum hardware for the $N=7$ LMG model with parameters $\frac{V}{\epsilon}=0.5$, $\frac{W}{\epsilon}=0$. Known exact eigenvalues are shown with dashed lines. Four random starting points are shown corresponding to runs which find the four eigenstates. Table 2: $N=7$ LMG model results calculated using the IBMQ_manila quantum computer using 20000 shots.  Eigenstate | Exact Value ($\epsilon$) | Variance | QC Result ($\epsilon$)
---|---|---|---
Ground | $-6.208$ | 0.139 | -6.067 $\pm$ 0.901
1st | $-2.944$ | 0.016 | -3.151 $\pm$ 0.503
2nd | $1.208$ | 0.010 | 1.184 $\pm$ 0.484
3rd | $5.944$ | 0.114 | 5.902 $\pm$ 0.660

In order to explore further the quality of the solutions obtained we calculate
the fidelity, or overlap of the wavefunctions obtained by the quantum computer
with the known analytic wavefunctions for each eigenvalue of the $N=7$ LMG
model. These overlaps are given by

$\langle\psi_{A}|\psi_{QC}\rangle^{2}$ (19)

where $\psi_{A}$ is the wavefunction of the analytic solution, and $\psi_{QC}$
is the wavefunction determined by the quantum computer. The results of the
overlap calculations are shown in Table 3.

Table 3: Overlaps of the $N=7$ LMG model result wavefunctions with analytic solution wavefunctions. Measured Eigenstate | $-6.2$ | $-2.94$ | $1.2$ | $5.94$
---|---|---|---|---
Ground | $0.924$ | $0.073$ | $0.002$ | $0.001$
1st | $0.020$ | $0.947$ | $0.033$ | $<0.001$
2nd | $0.001$ | $0.009$ | $0.991$ | $<0.001$
3rd | $0.001$ | $<0.001$ | $0.026$ | $0.972$

Here we see that all results found using variance minimization on quantum
hardware find states that closely overlap with their respective known analytic
eigenvector. The ground state introduces slightly more overlap with the first
excited state, which could be due to only achieving a partial convergence in
the time available for the run on real quantum hardware, as can be seen in
Fig. 7.

##  VI Conclusion

We adapted the Variational Quantum Eigensolver (VQE) to target any eigenstate
of a Hamiltonian through minimization of the variance of the Hamiltonian. The
cost is in the form of increased circuit evaluations compared with a ground
state VQE, but not circuit depth or number of qubits. States close to zero
variance could be found on presently-available quantum computing hardware,
corresponding to the known eigenstates of a model Hamiltonian. Further
exploration of larger problems would be needed to check how reliably one could
extract a full spectrum with the necessarily increasingly complex circuits
needed to parameterize larger problems.

###### Acknowledgements.

This work was funded by AWE. We acknowledge the use of IBM Quantum services
for this work. The views expressed are those of the authors, and do not
reflect the official policy or position of IBM or the IBM Quantum team. In
this paper we used IBMQ_manilla and IBMQ_nairobi, which are IBM Quantum
Processors. PDS acknowledges support from the UK STFC under grant numbers
ST/V001108/1 and ST/W006472/1. UK Ministry of Defence ©Crown owned copyright
2024/AWE.

## References

  * Višňák [2015] J. Višňák, Quantum algorithms for computational nuclear physics, EPJ Web of Conferences 100, 01008 (2015).
  * Višňák and Veselý [2017] J. Višňák and P. Veselý, Quantum algorithms for computational nuclear physics revisited, particular case of second quantized formulation, EPJ Web of Conferences 154, 01030 (2017).
  * Dumitrescu  _et al._ [2018] E. F. Dumitrescu, A. J. McCaskey, G. Hagen, G. R. Jansen, T. D. Morris, T. Papenbrock, R. C. Pooser, D. J. Dean, and P. Lougovski, Cloud Quantum Computing of an Atomic Nucleus, Physical Review Letters 120, 210501 (2018).
  * Roggero  _et al._ [2020a] A. Roggero, C. Gu, A. Baroni, and T. Papenbrock, Preparation of excited states for nuclear dynamics on a quantum computer, Physical Review C 102, 1 (2020a).
  * Roggero  _et al._ [2020b] A. Roggero, A. C. Li, J. Carlson, R. Gupta, and G. N. Perdue, Quantum computing for neutrino-nucleus scattering, Physical Review D 101, 74038 (2020b).
  * Lacroix [2020] D. Lacroix, Symmetry-Assisted Preparation of Entangled Many-Body States on a Quantum Computer, Physical Review Letters 125, 230502 (2020).
  * Siwach and Lacroix [2021] P. Siwach and D. Lacroix, Filtering states with total spin on a quantum computer, Physical Review A 104, 062435 (2021).
  * Baroni  _et al._ [2022] A. Baroni, J. Carlson, R. Gupta, A. C. Y. Li, G. N. Perdue, and A. Roggero, Nuclear two point correlation functions on a quantum computer, Phys. Rev. D 105, 074503 (2022).
  * Robin  _et al._ [2021] C. Robin, M. J. Savage, and N. Pillet, Entanglement rearrangement in self-consistent nuclear structure calculations, Physical Review C 103, 034325 (2021).
  * Zhang  _et al._ [2021] D.-B. Zhang, H. Xing, H. Yan, E. Wang, and S.-L. Zhu, Selected topics of quantum computing for nuclear physics, Chinese Physics B 30, 020306 (2021).
  * Robbins and Love [2021] K. Robbins and P. J. Love, Benchmarking near-term quantum devices with the variational quantum eigensolver and the lipkin-meshkov-glick model, Phys. Rev. A 104, 022412 (2021).
  * Pérez-Fernández  _et al._ [2022] P. Pérez-Fernández, J.-M. Arias, J.-E. García-Ramos, and L. Lamata, A digital quantum simulation of the Agassi model, Physics Letters B , 137133 (2022).
  * Cervia  _et al._ [2021] M. J. Cervia, A. B. Balantekin, S. N. Coppersmith, C. W. Johnson, P. J. Love, C. Poole, K. Robbins, and M. Saffman, Lipkin model on a quantum computer, Physical Review C 104, 024305 (2021).
  * Kiss  _et al._ [2022] O. Kiss, M. Grossi, P. Lougovski, F. Sanchez, S. Vallecorsa, and T. Papenbrock, Quantum computing of the 6Li nucleus via ordered unitary coupled clusters, Physical Review C 106, 034325 (2022).
  * Romero  _et al._ [2022] A. M. Romero, J. Engel, H. L. Tang, and S. E. Economou, Solving nuclear structure problems with the adaptive variational quantum algorithm, Physical Review C 105, 064317 (2022).
  * Chikaoka and Liang [2022] A. Chikaoka and H. Liang, Quantum computing for the Lipkin model with unitary coupled cluster and structure learning ansatz, Chinese Physics C 46, 024106 (2022).
  * Hlatshwayo  _et al._ [2022] M. Q. Hlatshwayo, Y. Zhang, H. Wibowo, R. LaRose, D. Lacroix, and E. Litvinova, Simulating excited states of the Lipkin model on a quantum computer, Physical Review C 106, 024319 (2022).
  * Ruiz Guzman and Lacroix [2022] E. A. Ruiz Guzman and D. Lacroix, Accessing ground-state and excited-state energies in a many-body system after symmetry restoration using quantum computers, Physical Review C 105, 024324 (2022).
  * Hobday  _et al._ [2022] I. Hobday, P. D. Stevenson, and J. Benstead, Quantum computing calculations for nuclear structure and nuclear data, in  _Quantum Technologies 2022_, Vol. 12133, International Society for Optics and Photonics (SPIE, 2022) p. 109.
  * Stetcu  _et al._ [2022] I. Stetcu, A. Baroni, and J. Carlson, Variational approaches to constructing the many-body nuclear ground state for quantum computing, Phys. Rev. C 105, 064308 (2022).
  * Li  _et al._ [2023] Y. H. Li, J. Al-Khalili, and P. Stevenson, A quantum simulation approach to implementing nuclear density functional theory via imaginary time evolution (2023), arXiv:2308.15425 .
  * Feynman [1982] R. P. Feynman, Simulating physics with computers, International Journal of Theoretical Physics 21, 467 (1982).
  * O’Brien  _et al._ [2019] T. E. O’Brien, B. Tarasinski, and B. M. Terhal, Quantum phase estimation of multiple eigenvalues for small-scale (noisy) experiments, New Journal of Physics 21, 023022 (2019).
  * Zhou  _et al._ [2013] X.-Q. Zhou, P. Kalasuwan, T. C. Ralph, and J. L. O’brien, Calculating unknown eigenvalues with a quantum algorithm, Nature photonics 7, 223 (2013).
  * Du  _et al._ [2021] W. Du, J. P. Vary, X. Zhao, and W. Zuo, Ab initio nuclear structure via quantum adiabatic algorithm (2021), arXiv:2105.08910 .
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nature Physics 16, 205 (2020).
  * Turro  _et al._ [2022] F. Turro, A. Roggero, V. Amitrano, P. Luchi, K. A. Wendt, J. L. Dubois, S. Quaglioni, and F. Pederiva, Imaginary-time propagation on a quantum chip, Phys. Rev. A 105, 022440 (2022).
  * Jouzdani  _et al._ [2022] P. Jouzdani, C. W. Johnson, E. R. Mucciolo, and I. Stetcu, Alternative approach to quantum imaginary time evolution, Phys. Rev. A 106, 062435 (2022).
  * Preskill [2018] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 4213 (2014).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, The variational quantum eigensolver: a review of methods and best practices (2022).
  * Lipkin  _et al._ [1965] H. Lipkin, N. Meshkov, and A. Glick, Validity of many-body approximation methods for a solvable model: (i). exact solutions and perturbation theory, Nuclear Physics 62, 188 (1965).
  * Agassi  _et al._ [1966] D. Agassi, H. J. Lipkin, and N. Meshkov, Validity of many-body approximation methods for a solvable model: (IV). The deformed Hartree-Fock solution, Nuclear Physics 86, 321 (1966).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Über das Paulische Äquivalenzverbot, Zeitschrift für Physik 47, 631 (1928).
  * Bravyi and Kitaev [2002] S. B. Bravyi and A. Y. Kitaev, Fermionic quantum computation, Annals of Physics 298, 210 (2002).
  * Seeley  _et al._ [2012] J. T. Seeley, M. J. Richard, and P. J. Love, The Bravyi-Kitaev transformation for quantum computation of electronic structure, The Journal of Chemical Physics 137, 224109 (2012).
  * Pesce and Stevenson [2021] R. M. N. Pesce and P. D. Stevenson, H2ZIXY: Pauli spin matrix decomposition of real symmetric matrices (2021), arXiv:2111.00627 .
  * Bergholm  _et al._ [2022] V. Bergholm, J. Izaac, M. Schuld, C. Gogolin, S. Ahmed, V. Ajith, M. S. Alam, G. Alonso-Linaje, B. AkashNarayanan, A. Asadi, J. M. Arrazola, U. Azad, S. Banning, C. Blank, T. R. Bromley, B. A. Cordier, J. Ceroni, A. Delgado, O. D. Matteo, A. Dusko, T. Garg, D. Guala, A. Hayes, R. Hill, A. Ijaz, T. Isacsson, D. Ittah, S. Jahangiri, P. Jain, E. Jiang, A. Khandelwal, K. Kottmann, R. A. Lang, C. Lee, T. Loke, A. Lowe, K. McKiernan, J. J. Meyer, J. A. Montañez-Barrera, R. Moyard, Z. Niu, L. J. O’Riordan, S. Oud, A. Panigrahi, C.-Y. Park, D. Polatajko, N. Quesada, C. Roberts, N. Sá, I. Schoch, B. Shi, S. Shu, S. Sim, A. Singh, I. Strandberg, J. Soni, A. Száva, S. Thabet, R. A. Vargas-Hernández, T. Vincent, N. Vitucci, M. Weber, D. Wierichs, R. Wiersema, M. Willmann, V. Wong, S. Zhang, and N. Killoran, Pennylane: Automatic differentiation of hybrid quantum-classical computations (2022), arXiv:1811.04968 .
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).