# Strong-coupling phases of trions and excitons in electron-hole bilayers at
commensurate densities

David D. Dai  Liang Fu  Department of Physics, Massachusetts Institute of
Technology, Cambridge MA 02139, USA

(August 28, 2024)

###### Abstract

We introduce density imbalanced electron-hole bilayers at a commensurate $2:1$
density ratio as a platform for realizing novel phases involving electrons,
excitons and trions. Three length scales are identified which characterize the
interplay between kinetic energy, intralayer repulsion, and interlayer
attraction. By a combination of theoretical analysis and numerical
calculation, we find a variety of strong-coupling phases in different
parameter regions, including quantum crystals of electrons, excitons, and
trions. We also propose an “excitonic supersolid” phase that features electron
crystallization and exciton superfluidity simultaneously. The material
realization and experimental signature of these phases are discussed in the
context of semiconductor transition metal dichalcogenide bilayers.

##  I Introduction

Recently, semiconductor transition metal dichalcogenide (TMD) heterostructures
have proven to be an ideal platform for exploring quantum phases of matter in
a controlled, tunable, and reproducible way. An extraordinarily rich variety
of quantum states have been predicted and observed, including Mott-Hubbard and
charge-transfer insulators Wu et al. (2018); Zhang et al. (2020a); Regan et
al. (2020); Tang et al. (2020); Li et al. (2021a); Park et al. (2023), Wigner
crystals Regan et al. (2020); Xu et al. (2020); Zhou et al. (2021); Jin et al.
(2021a); Li et al. (2021b); Huang et al. (2021); Xu et al. (2022); Padhi et
al. (2021); Jin et al. (2021b), itinerant ferromagnets Zhang et al. (2020b);
Davydova et al. (2023); Ciorciaro et al. (2023); Lee et al. (2023); Crépel and
Fu (2023a); Morera et al. (2023), interfacial ferroelectrics Wang et al.
(2022), heavy Fermi liquids Zhao et al. (2023a); Kumar et al. (2022); Guerci
et al. (2023), and spin-polaron liquids with pseudogap and magnetization
plateau Davydova et al. (2023); Zhang and Fu (2023); Tao et al. (2023), as
well as quantum spin Hall states Wu et al. (2019); Zhao et al. (2023b) and
quantum anomalous Hall states Devakul et al. (2021); Li et al. (2021c); Crépel
and Fu (2023b); Morales-Durán et al. (2023); Reddy et al. (2023); Wang et al.
(2023); Li et al. (2021d); Kou et al. (2014); Cai et al. (2023); Zeng et al.
(2023). It is remarkable that all these electronic phases realized in a single
material system are rooted in one common ground—the two-dimensional electron
(or hole) gas in monolayer TMDs. Here the large effective mass favors strong
interaction effects, and the presence of moiré bands in TMD heterostructures
further enriches the physics, leading to much of the observed phenomena.

In addition to moiré physics, TMD heterostructures provide a material
realization of electron-hole ($e$-$h$) bilayer systems featuring electrons and
holes on spatially separated layers, whose densities can be independently
tuned by top and bottom gate voltages Zeng and MacDonald (2020). Owing to the
coexistence of positively and negatively charged particles, both repulsion
between like charges and attraction between opposite charges are present,
enabling new phases of matter. When the densities of electrons and holes are
equal, interlayer excitons with intrinsic out-of-plane dipole moments form,
which may support high-temperature exciton superfluidity in an electrical
insulator. Recently, thermodynamic evidence of excitonic insulator ground
states has been observed in TMD bilayer WSe2/MoSe2 with WSe2 as the hole layer
and MoSe2 as the electron layer separated by insulating hBN layers Ma et al.
(2021a); Qi et al. (2023). Moreover, density imbalanced electron-hole bilayers
have attracted increasing interest owing to the interaction between charge
carriers and dipolar excitons Ma et al. (2021b).

In this work, we study the strong-coupling phases of an imbalanced electron-
hole bilayer at commensurate electron and hole densities $n_{e}/n_{h}=2$. Our
study is motivated by the following considerations. First, because there is a
net charge density $n=n_{e}-n_{h}\neq 0$, the Coulomb interaction between
charged particles is crucial in the low density regime, favoring strong-
coupling phases with crystalline order. This should be contrasted to balanced
electron-hole bilayers, where charge-neutral excitons condense into a
superfluid at low densities because their mutual dipole-dipole interaction is
parametrically weaker than the quantum kinetic energy. Second, the particular
choice of the electron-to-hole density ratio $n_{e}/n_{h}=2$ is motivated by
the prospect of three-body bound states known as trions, which are charge-$e$
composite particles made of two electrons on the same layer bound to a hole on
the other layer.

Through both theoretical analysis and numerical calculation, we find a number
of ordered phases driven by strong interactions. These include bilayer
electron-hole Wigner crystals, “composite crystals” of coexisting electrons
and excitons, as well as exotic quantum phases without classical counterparts.
In particular, we find a quantum supersolid made of electrons and excitons as
well as a quantum crystal of trions in the low density regime. Experimental
signatures of various phases are also discussed.

Before presenting our main results, it is instructive to first consider the
three-body problem of two electrons and one hole, which reside on spatially
separated layers ($z=d$ and $0$) and mutually interact through Coulomb forces.
Classically, the minimum energy configuration is an electron and a dipole that
are far apart due to their residual repulsion. This dipole (equivalent to a
classical exciton) consists of an electron and a hole sitting directly on top
of each other, separated only by the layer distance $d$. For comparison,
consider a “trion” charge cluster with the two electrons at
$(\pm{\mathbf{r}},d)$ and the hole at $(\mathbf{0},0)$. When
$r=1/\sqrt{2^{4/3}-1}\approx 0.811d$, the net force acting on each particle is
zero. However, this force-balanced configuration has energy $-0.937/d$, which
is higher than the energy $-1/d$ of a dipole plus an electron far away.
Additionally, this classical “trion” has an unstable normal mode, where the
hole moves towards one of the electrons along the $\mathbf{r}$ direction,
making the classical “trion” a saddle point of the energy instead of a local
minima.

Our above analysis of the classical three-body problem demonstrates that
quantum mechanics is crucial for the formation of a trion bound state. As
shown by numerical studies, the trion is the ground state of two electrons and
one hole in the bilayer when $a_{B}/d>0.065$ assuming equal electron and hole
masses Witham et al. (2018), where
$a_{B}=\frac{4\pi\epsilon\hbar^{2}}{e^{2}m}$ is the Bohr radius (which
vanishes in the classical limit $\hbar\rightarrow 0$ or $m\rightarrow\infty$).
We can understand the trion’s quantum origin heuristically. Start with the
electron and the dipolar exciton separated by a large distance $r\gg d$. The
electron exerts a repulsive force on the exciton $F\approx d^{2}/r^{4}$,
leading to a repulsive potential energy $U\propto d^{2}/r^{3}$. On the other
hand, the in-plane electric field $\propto 1/r^{2}$ polarizes the exciton and
lowers its energy through the second-order Stark effect by an amount $\delta
E_{s}\propto-a_{B}^{3}/r^{4}$ to zeroth order in $d$. Therefore, the
attraction due to the quantum Stark effect dominates the direct electron-
exciton repulsion over a range of distances below a crossover length
$r_{c}\sim a_{B}^{3}/d^{2}$. For sufficiently small $d$, this extended
attraction supports a bound state of the electron and the exciton, leading to
trion formation.

At finite charge density, an additional length scale appears: the average
inter-particle distance $a\equiv 1/\sqrt{\pi n}$, with $n=n_{h}$ for
$n_{e}/n_{h}=2$. With the presence of three length scales—the inter-particle
distance $a$, the layer distance $d$, and Bohr radius $a_{B}$, we expect that
the competition between intralayer repulsion, interlayer electron-hole
attraction, and quantum kinetic energy leads to a rich phase diagram for the
electron-hole bilayer, which we explore below by a combination of analytical
and numerical methods.

The Hamiltonian for the bilayer assuming equal electron and hole effective
masses $m$ and $1/r$ Coulomb interactions is:

$\begin{split}H&=\sum_{a=e,h}\sum_{s=\uparrow,\downarrow}\int\text{d}^{2}{\mathbf{r}}\bigg{[}\psi^{a\dagger}_{s}({\mathbf{r}})\left(\frac{-\nabla^{2}}{2}\right)\psi^{a}_{s}({\mathbf{r}})\bigg{]}\\\
&+\frac{1}{2}\int\text{d}^{2}{\mathbf{r}}\int\text{d}^{2}{\mathbf{r}^{\prime}}\bigg{[}\sum_{a=e,h}\frac{n^{a}({\mathbf{r}})n^{a}({\mathbf{r}^{\prime}})}{|\mathbf{r-r^{\prime}}|}\bigg{]}\\\
&-\int\text{d}^{2}{\mathbf{r}}\int\text{d}^{2}{\mathbf{r}^{\prime}}\bigg{[}\frac{n^{e}({\mathbf{r}})n^{h}({\mathbf{r}^{\prime}})}{\sqrt{|\mathbf{r-r^{\prime}}|^{2}+d^{2}}}\bigg{]}\end{split}$
(1)

where $\psi^{e,h}({\mathbf{r}})$ denotes the electron or hole operator,
$n^{a}({\mathbf{r}})=\sum_{s}\psi^{a\dagger}_{s}({\mathbf{r}})\psi^{a}_{s}({\mathbf{r}})$
is the density operator, and $s$ denotes the spin. We have already divided all
quantities by their appropriate atomic units, i.e. lengths by the Bohr radius
$a_{B}=\frac{4\pi\epsilon\hbar^{2}}{me^{2}}$ (we use this definition
throughout ) and energies by the Hartree energy
$E_{h}=\frac{\hbar^{2}}{ma_{B}^{2}}$. Unless stated otherwise, we assume for
simplicity that the electron and hole masses are equal, noting that our main
findings are qualitatively correct for a range of mass ratios.

##  II Classical Regime

We first consider the classical limit defined by taking $a_{B}\rightarrow 0$
while keeping $d$ and $a$ fixed, or equivalently, $a/a_{B}\rightarrow\infty$
and $d/a_{B}\rightarrow\infty$ for a fixed $d/a$. In the classical limit, the
phase diagram is characterized solely by the dimensionless ratio $d/a$. Both
limits $d/a\gg 1$ and $d/a\ll 1$ can be understood analytically. When $d/a\gg
1$, the interlayer coupling is negligible and a triangular lattice Wigner
crystal is formed independently in each layer. In the opposite regime $d/a\ll
1$, every hole pairs with an electron at the shortest possible distance $d$ to
form a small out-of-plane dipole, leaving an equal number of excess electrons.
These dipoles have a weak repulsion $1/r-1/\sqrt{r^{2}+d^{2}}\propto
d^{2}/r^{3}$ with the electrons, while electrons interact with each other
through a strong Coulomb repulsion $1/r$. Therefore, the electrons arrange
themselves essentially independently of the dipoles to form a triangular
lattice Wigner crystal. Once the electrons crystallize, the dipoles
crystallize in half of the voids between the electrons to minimize the
residual repulsion. The result is a “composite crystal” of electrons and
dipoles, taking the form of a honeycomb lattice with one sublattice occupied
by electrons and the other by dipoles.

Figure 1: The classical phase diagram of the electron-hole bilayer at
$n_{e}/n_{h}=2$ for $a/d>1$. The top panel depicts the checkerboard composite
crystal (left) and the honeycomb composite crystal (right). The bottom panel
shows the energy and energy difference of the two composite crystals, with a
dashed line indicating the phase transition.

To determine the classical ground state as a function of the charge density,
we minimized the classical electrostatic energy. To handle the Coulomb
interaction’s long-ranged tail, we used Ewald summation (details are provided
in the supplementary material). For the experimentally relevant parameter
range $a/d>1$, we find two distinct types of composite crystals. For
$a/d>5.42$, the numerical results confirm our expectation of a honeycomb
composite crystal in the low-density limit. Interestingly, for $a/d<5.42$, we
found a “checkerboard” composite crystal, which consists of two
interpenetrating square lattices of electrons and dipoles. This is likely due
to the increased electron-dipole repulsion at smaller $a/d$. Energy
minimization calculations with randomly initialized configurations generally
converged to a mixture of honeycomb and checkerboard structures alongside
various defects, with the former (latter) more predominant at large (small)
$a/d$. We did not observe the formation of classical trions in any of our
randomly initialized runs, which is consistent with the classical trion’s
instability. Once these two composite crystals were identified, we directly
calculated their electrostatic energy as a function of $a/d$ and determined
the phase transition point, which is shown in Fig. 1.

We now analyze the effect of quantum fluctuations around the classical
composite crystals. At sufficiently small $a_{B}$, the leading quantum effect
is the zero-point-motion of charges in the classical ground state. The root-
mean-square displacement of charges increases with $a_{B}$, and when it
becomes comparable to the lattice constant $a$, quantum melting of the crystal
takes place (Lindemann criterion). Since our composite crystals have three
charges per unit cell, there are a total of six phonon modes, including both
acoustic and optical branches. Notably, the optical phonons are associated
with the relative vibration of charges within the unit cell, which is absent
in the canonical electron Wigner crystal.

When the layer distance $d$ is small compared to the inter-particle distance
$a$, the repulsion between dipoles and electrons in the composite crystal is
much weaker than the repulsion between electrons. In this case, we expect that
the dipole’s center of mass has the largest zero-point displacement, denoted
as $\xi_{d}$. This is indeed confirmed by our direct calculation of the
optical phonon frequencies at zero wavevector. Two types of optical phonons
are present: the low-frequency one corresponds to the displacement of the
dipole’s center of mass relative to the electron, while the high-frequency one
is associated with the internal structure of the dipole. At small $d/a$, the
low-frequency optical phonon softens with
$\omega_{d}\propto\sqrt{d^{2}/ma^{5}}$, hence the zero-point displacement
increases and is given by:

$\displaystyle\xi_{d}\sim\sqrt{\frac{\hbar}{m\omega_{d}}}\propto
a(\frac{aa_{B}}{d^{2}})^{1/4}.$ (2)

In the classical limit $a_{B}/a\rightarrow 0,a_{B}/d\rightarrow 0$ with $a/d$
fixed, $\xi_{d}\ll a$ ensures the stability of the classical composite crystal
against quantum fluctuation.

For small $d/a$, Eq.(2) implies that quantum melting of the dipole
(=interlayer exciton) lattice in the composite crystal occurs first when
$a_{B}$ reaches the order of $d^{2}/a$. On the other hand, the stronger
electron-electron repulsion renders the electron lattice stable against
quantum fluctuations until $a_{B}$ reaches the order of the inter-particle
distance $a$, as in the case of the canonical Wigner crystal. This naturally
raises the question of what the ground state is at intermediate $a_{B}$
between $d^{2}/a$ and $a$, where both quantum effects and electrostatic
interactions play crucial roles.

Figure 2: Soft and hard optical phonon frequencies at the center of the
Brillouin zone as function of $a/d$ for the classical composite crystals. The
low-frequency branches correspond to the displacement of the dipole’s center
of mass relative to the electron, while the high-frequency branches are
associated with oscillations in the dipole’s internal structure.

##  III Dilute Quantum Regime

Of particular interest is the low-density limit defined by
$a\rightarrow\infty$ while keeping $d$ and $a_{B}$ fixed, especially if $d$
and $a_{B}$ are of the same order. With $a\rightarrow\infty$, $a_{B}$
necessarily lies between $d^{2}/a$ and $a$. Therefore, we expect partial
melting of the composite crystal due to the quantum motion of the exciton’s
center of mass.

The ground state in this low density limit depends crucially on the ratio of
the Bohr radius $a_{B}$ and layer distance $d$. For $d<d_{c}\approx
15.38a_{B}$ (assuming equal electron and hole masses), it is known that the
lowest-energy state of two electrons and one hole in a bilayer system is a
trion bound state with the two electrons in a spin-singlet Witham et al.
(2018). (Note that the spin-triplet trion is absent for $m_{e}=m_{h}$.) The
binding energy of a trion is on the order of $E_{h}$ at $d=0$ and decreases as
$d$ approaches the critical $d_{c}$ at which the trion unbinds. At a given
$d<d_{c}$ and in the dilute limit $a\rightarrow\infty$, the hierarchy of
energy scales is necessarily as follows: the trion binding energy is much
larger than the characteristic Coulomb energy $\propto 1/a$, which is much
larger than the kinetic energy $\propto 1/a^{2}$. Therefore, we conclude that
for $d<d_{c}$, the ground state of electron-hole bilayers in the low density
limit is a triangular lattice Wigner crystal of spin-singlet trions.

The physically realistic parameters for the TMD bilayer
$\text{WSe}_{2}$/$\text{MoSe}_{2}$ are electron effective mass
$m_{e}(\text{MoSe}_{2})=0.8m_{0}$ ($m_{0}$ is the bare mass), hole effective
mass $m_{h}(\text{WSe}_{2})=0.4m_{0}$, dielectric constant
$\epsilon=4.7\epsilon_{0}$, and interlayer spacing $d=3\text{ nm}$. With these
numbers, the trion binding energy predicted by numerical studies is on the
order of $15\text{ K}$ Witham et al. (2018).

Notably, trions can crystallize at significantly higher densities than
electrons. Because the trion mass is three times the electron mass, the
effective Bohr radius is three times smaller, hence the critical density for
crystal melting according to the Lindemann criterion is nine times larger.
Another distinctive feature of our trion Wigner crystal is that the electron
layer has zero total spin and a large spin gap on the order of $E_{h}$ due to
the trion binding energy, whereas the spins of localized holes interact with
each other through a weak exchange interaction which goes to zero in the low
density limit. Therefore, a small magnetic field can fully polarize the holes
while leaving the spin-singlet electrons intact. This sharp contrast between
the spin response of electron and hole layers is an indication of spin-singlet
trion formation.

Finally, we discuss the regime $d>d_{c}$. Here, the trion is not stable at the
three-particle level, i.e. the electron and the exciton unbind. On the other
hand, as we showed earlier, the exciton lattice in the composite crystal state
is necessarily unstable against quantum melting in the low density limit,
whereas the electron lattice is stable. Based on these considerations, we
propose that for $d>d_{c}$, the ground state at sufficiently low density is an
exciton superfluid that permeates through the electron crystal. This state is
remarkable as it simultaneously exhibits crystallization and superfluidity.
For this reason, we call this state an “excitonic supersolid”, a quantum
electron solid in which interstitial excitons Bose condense.

##  IV Hartree-Fock Calculations

Guided by our theoretical results in the low-density regimes, we perform
further numerical study to strengthen our prediction of singlet trion
crystals. Specifically, we apply the Hartree-Fock method to $N_{e}=72$
electrons and $N_{h}=36$ holes in a periodic rectangle. Our system consists of
four fermion species including spin $s=\uparrow,\downarrow$ and electron/hole
degrees of freedom, and our ansatz for the HF ground state is a product of
Slater determinants, one for each fermion species. No restrictions other than
orthonormality are applied to the Hartree-Fock orbitals
$\phi_{n,s}(\mathbf{r})$. As discussed previously, the coupling between hole
spins is weak and plays an unimportant role, so we assume that the holes are
fully spin-polarized, reducing the effective number of species to three. While
our ansatz is capable of describing the crystalline phases, it cannot capture
coherences between electrons and holes and therefore cannot describe the
proposed excitonic supersolid, which we defer to a future study.

Instead of the usual self-consistent field method, we use conjugate gradients
(any other unconstrained minimization algorithm would also work) to directly
minimize an analytic continuation of the Hartree-Fock energy functional onto
the space of non-orthonormal orbitals Arias et al. (1992). Although formally
equivalent to the self-consistent field method, this method is immune from
instabilities such as charge sloshing and is guaranteed to converge from any
starting point. Both of these features are convenient for exploratory
calculations with unknown ground states. We supply a detailed derivation of
our working equations in the supplementary materials.

Although our predictions are general, we choose experimentally relevant
parameters $m_{e}/m_{h}=1/4$, interlayer spacing $d=0.5a_{B,h}$, where
$a_{B,h}=\frac{4\pi\epsilon\hbar^{2}}{e^{2}m_{h}}$ is the hole Bohr radius,
and hole density parameter $r_{s}=7$ to provide a concrete example of our
results. We find that with equal electron spin populations
$n_{e,\uparrow}=36,n_{e,\downarrow}=36$, randomly initialized HF calculations
converge naturally to a spin-singlet state with equivalent electron $\uparrow$
and $\downarrow$ orbitals (apart from small imperfections). Additionally, the
electrons are concentrated on top of the holes, which are arranged on a
triangular lattice (with a few defects). The electron orbitals also spread out
more than the hole orbitals due to the electrons’ repulsion and smaller
effective mass.

Figure 3: The total (spin $\uparrow$ plus $\downarrow$) electron density is
shown in the left panel, and the hole density is shown in the right panel; the
coloring scheme is consistent for both densities and is shown in the colorbar
in units of $1/a_{B}^{2}$ for
$a_{B}=\frac{4\pi\epsilon\hbar^{2}}{e^{2}m_{h}}$. The calculation parameters
are $m_{e}=1/4$, $m_{h}=1$, $d=0.5a_{B}$, hole $r_{s}=7$, $n_{h}=36$, and
equal electron spin-populations.

Motivated by these randomly initialized calculations, we perform HF
calculations starting directly from singlet trions placed on a triangular
lattice. These calculations yield an energy per particle $-0.04300$, which is
slightly lower than the best outcomes from the randomly initialized runs. We
also find that the electron HF energy spectrum has a gap of $0.06082$, and the
hole HF energy spectrum has a gap of $0.21220$ (all energies are in terms of
the hole Hartree energy). The existence of a gap is consistent with the
stability of singlet trions.

For comparison, randomly initialized HF calculations in the fully spin-
polarized sector yield energies per particle in excess of $-0.02$,
significantly higher than that of the spin-unpolarized sector. While the hole
HF spectrum is still gapped, the electron HF spectrum appears gapless in the
fully spin-polarized sector. Therefore, our calculations clearly show that the
ground state of the electron-hole bilayer is an insulating quantum crystal of
spin-singlet trions for these parameters. Additionally, applying an external
magnetic field can melt the singlet trion crystal into a conducting state.

It should be noted that the Hartree-Fock method generally overestimates the
tendency for interactions to break symmetries. For example, Hartree-Fock
underestimates the critical $r_{s}$ at which the uniform electron gas forms a
Wigner crystal Trail et al. (2003); Drummond and Needs (2009). For our
electron-hole bilayer, it is important to combine an HF study with more
advanced numerical methods in order to accurately determine the parameter
regions for interaction-induced phases, including the trion crystal, the
composite crystals, and the putative excitonic supersolid. We leave a
comprehensive numerical study of the electron-hole bilayer at commensurate
$2:1$ density ratio to a future work.

##  V Conclusion

Figure 4: Schematic phase diagram for the electron-hole bilayer at
commensurate densities $n_{e}/n_{h}=2$, with analytically known boundaries
marked in red and limiting phases labelled.

Our work introduces the electron-hole bilayer at commensurate 2:1 density
ratio as a platform for realizing novel strong-coupling phases. We identify
two dimensionless ratios that govern the phase diagram: the ratio of the
average interparticle spacing to the Bohr radius $a/a_{B}$, and the ratio of
the interlayer spacing to the Bohr radius $d/a_{B}$. By a combination of
theoretical analysis and numerical calculation, we find a number of ordered
phases, shown schematically in Fig. 4. In the classical regime $a/a_{B}\gg
1,d/a_{B}\gg 1$, we identify three crystalline phases: decoupled crystals in
each layer for small $a/d$, a checkerboard “composite crystal” comprised of
electrons and bound electron-hole dipoles for intermediate $a/d$, and a
honeycomb composite crystal for large $a/d$. Focusing on the dilute regime
$a\gg d,a_{B}$, as $d$ is decreased, we propose that large zero-point-
fluctuations of the excitons partially melt the composite crystal into an
“excitonic supersolid”, where the unbound electrons remain crystalline but the
excitons condense into a superfluid. As $d$ is decreased further to the order
of the Bohr radius $a_{B}$, we show that the second-order Stark effect
mediates an effective attractive force between electrons and excitons. For
sufficiently small $d$, this attraction binds excitons to electrons to form
spin-singlet trions, which then subsequently crystallize into a trion crystal.

Notably, the trion crystal is stable throughout the low charge density regime
with a lattice constant varying continuously with the charge density
$n=n_{e}-n_{h}$, provided that the electron-to-hole density ratio is
maintained at the commensurate value $n_{e}/n_{h}=2$. Thus, the trion crystal
is charge-compressible with $\frac{\partial n}{\partial\mu}\neq 0$ with $\mu$
the charge chemical potential, but has an energy gap to adding excitons, i.e.,
exciton-incompressible. This should be contrasted with the excitonic insulator
at charge neutrality $n=0$ which is charge-incompressible and exciton-
compressible. Remarkably, recent capacitance and optical experiments on
WSe2/MoSe2 have shown that the charge and exciton compressibility can be
measured independently by varying the top and bottom gate voltages
concurrently with $\Delta V_{B}=\pm\Delta V_{T}$. Moreover, the formation of
trion is accompanied by a reduction of spin susceptibility in the electron
(majority carrier) layer, which can be detected by magnetic circular
dichroisum. Finally, it is noted that the binding energy of trions can be
further increased by applying an magnetic field, which we leave to future
study. We hope our theoretical work stimulates experimental study of TMD
electron-hole bilayers at commensurate electron-hole density ratio.

## Acknowledgements

We are grateful to Trithep Devakul and Aidan Reddy for helpful discussions. We
thank Kin Fai Mak, Jie Shan and Feng Wang for their interest and feedback.
This work was supported by the Simons Investigator Award from the Simons
Foundation. DDD was supported by the Undergraduate Research Opportunities
Program at MIT. The authors acknowledge the MIT SuperCloud and Lincoln
Laboratory Supercomputing Center for providing high performance computing
resources.

## References

  * Wu et al. (2018) F. Wu, T. Lovorn, E. Tutuc, and A. MacDonald, Physical Review Letters 121 (2018).
  * Zhang et al. (2020a) Y. Zhang, N. F. Yuan, and L. Fu, Physical Review B 102 (2020a).
  * Regan et al. (2020) E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, et al., Nature 579, 359–363 (2020).
  * Tang et al. (2020) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, et al., Nature 579, 353–358 (2020).
  * Li et al. (2021a) T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, et al., Nature 597, 350–354 (2021a).
  * Park et al. (2023) H. Park, J. Zhu, X. Wang, Y. Wang, W. Holtzmann, T. Taniguchi, K. Watanabe, J. Yan, L. Fu, T. Cao, et al., Nature Physics (2023).
  * Xu et al. (2020) Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, Nature 587, 214–218 (2020).
  * Zhou et al. (2021) Y. Zhou, J. Sung, E. Brutschea, I. Esterlis, Y. Wang, G. Scuri, R. J. Gelly, H. Heo, T. Taniguchi, K. Watanabe, et al., Nature 595, 48–52 (2021).
  * Jin et al. (2021a) C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, et al., Nature Materials 20, 940–944 (2021a).
  * Li et al. (2021b) H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, et al., Nature 597, 650–654 (2021b).
  * Huang et al. (2021) X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, et al., Nature Physics 17, 715–719 (2021).
  * Xu et al. (2022) Y. Xu, K. Kang, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, Nature Nanotechnology 17, 934–939 (2022).
  * Padhi et al. (2021) B. Padhi, R. Chitra, and P. W. Phillips, Physical Review B 103 (2021).
  * Jin et al. (2021b) C. Jin, Z. Tao, T. Li, Y. Xu, Y. Tang, J. Zhu, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, et al., Nature Materials 20, 940–944 (2021b).
  * Zhang et al. (2020b) Y. Zhang, H. Isobe, and L. Fu, arXiv preprint arXiv:2005.04238 (2020b).
  * Davydova et al. (2023) M. Davydova, Y. Zhang, and L. Fu, Physical Review B 107 (2023).
  * Ciorciaro et al. (2023) L. Ciorciaro, T. Smolenski, I. Morera, N. Kiper, S. Hiestand, M. Kroner, Y. Zhang, K. Watanabe, T. Taniguchi, E. Demler, et al., arXiv preprint arXiv:2305.02150 (2023).
  * Lee et al. (2023) K. Lee, P. Sharma, O. Vafek, and H. J. Changlani, Physical Review B 107 (2023).
  * Crépel and Fu (2023a) V. Crépel and L. Fu, Physical Review B 107 (2023a).
  * Morera et al. (2023) I. Morera, M. Kanász-Nagy, T. Smolenski, L. Ciorciaro, A. Imamoğlu, and E. Demler, Physical Review Research 5 (2023).
  * Wang et al. (2022) X. Wang, K. Yasuda, Y. Zhang, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, L. Fu, and P. Jarillo-Herrero, Nature Nanotechnology 17, 367–371 (2022).
  * Zhao et al. (2023a) W. Zhao, B. Shen, Z. Tao, Z. Han, K. Kang, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, Nature 616, 61–65 (2023a).
  * Kumar et al. (2022) A. Kumar, N. C. Hu, A. H. MacDonald, and A. C. Potter, Physical Review B 106 (2022).
  * Guerci et al. (2023) D. Guerci, J. Wang, J. Zang, J. Cano, J. H. Pixley, and A. Millis, Science Advances 9 (2023).
  * Zhang and Fu (2023) Y. Zhang and L. Fu, SciPost Physics Core 6 (2023).
  * Tao et al. (2023) Z. Tao, W. Zhao, B. Shen, P. Knüppel, K. Watanabe, T. Taniguchi, J. Shan, and K. F. Mak, arXiv preprint arXiv:2307.12205 (2023).
  * Wu et al. (2019) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. MacDonald, Physical Review Letters 122 (2019).
  * Zhao et al. (2023b) W. Zhao, K. Kang, L. Li, C. Tschirhart, E. Redekop, K. Watanabe, T. Taniguchi, A. Young, J. Shan, and K. F. Mak, arXiv preprint arXiv:2207.02312 (2023b).
  * Devakul et al. (2021) T. Devakul, V. Crépel, Y. Zhang, and L. Fu, Nature Communications 12 (2021).
  * Li et al. (2021c) H. Li, U. Kumar, K. Sun, and S.-Z. Lin, Physical Review Research 3 (2021c).
  * Crépel and Fu (2023b) V. Crépel and L. Fu, Physical Review B 107 (2023b).
  * Morales-Durán et al. (2023) N. Morales-Durán, J. Wang, G. R. Schleder, M. Angeli, Z. Zhu, E. Kaxiras, C. Repellin, and J. Cano, arXiv preprint arXiv:2304.06669 (2023).
  * Reddy et al. (2023) A. P. Reddy, F. F. Alsallom, Y. Zhang, T. Devakul, and L. Fu, arXiv preprint arXiv:2304.12261 (2023).
  * Wang et al. (2023) C. Wang, X.-W. Zhang, X. Liu, Y. He, X. Xu, Y. Ran, T. Cao, and D. Xiao, arXiv preprint arXiv:2304.11864 (2023).
  * Li et al. (2021d) T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, Z. Tao, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, et al., Nature 600, 641–646 (2021d).
  * Kou et al. (2014) A. Kou, B. E. Feldman, A. J. Levin, B. I. Halperin, K. Watanabe, T. Taniguchi, and A. Yacoby, Science 345, 55–57 (2014).
  * Cai et al. (2023) J. Cai, E. Anderson, C. Wang, X. Zhang, X. Liu, W. Holtzmann, Y. Zhang, F. Fan, T. Taniguchi, K. Watanabe, et al., Nature (2023).
  * Zeng et al. (2023) Y. Zeng, Z. Xia, K. Kang, J. Zhu, P. Knüppel, C. Vaswani, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, arXiv preprint arXiv:2305.00973 (2023).
  * Zeng and MacDonald (2020) Y. Zeng and A. H. MacDonald, Physical Review B 102 (2020).
  * Ma et al. (2021a) L. Ma, P. X. Nguyen, Z. Wang, Y. Zeng, K. Watanabe, T. Taniguchi, A. H. MacDonald, K. F. Mak, and J. Shan, Nature 598, 585–589 (2021a).
  * Qi et al. (2023) R. Qi, A. Y. Joe, Z. Zhang, Y. Zeng, T. Zheng, Q. Feng, E. Regan, J. Xie, Z. Lu, T. Taniguchi, et al., arXiv preprint arXiv:2306.13265 (2023).
  * Ma et al. (2021b) L. Ma, P. X. Nguyen, Z. Wang, Y. Zeng, K. Watanabe, T. Taniguchi, A. H. MacDonald, K. F. Mak, and J. Shan, Nature 598, 585–589 (2021b).
  * Witham et al. (2018) O. Witham, R. J. Hunt, and N. D. Drummond, Physical Review B 97 (2018).
  * Arias et al. (1992) T. A. Arias, M. C. Payne, and J. D. Joannopoulos, Physical Review Letters 69, 1077–1080 (1992).
  * Trail et al. (2003) J. R. Trail, M. D. Towler, and R. J. Needs, Physical Review B 68 (2003).
  * Drummond and Needs (2009) N. D. Drummond and R. J. Needs, Physical Review Letters 102 (2009).
  * Johnson and Ranganathan (2001) R. E. Johnson and S. Ranganathan, Physical Review E 63 (2001).
  * Sundararaman et al. (2017) R. Sundararaman, K. Letchworth-Weaver, K. A. Schwarz, D. Gunceler, Y. Ozhabes, and T. Arias, SoftwareX 6, 278–284 (2017).

Strong-coupling phases of trions and excitons in electron-hole bilayers at
commensurate densities: Supplementary Material
David D. Dai and Liang Fu

Department of Physics, Massachusetts Institute of Technology, Cambridge MA
02139, USA

##  I Ewald Summation

We use Ewald summation to efficiently calculate the classical electrostatic
energy and determine the classical phase diagram. Although the principle is
the same as in three dimensions, there are some additional subtleties in our
case, which we discuss here Johnson and Ranganathan (2001). Without Ewald
summation, the $v(r)=1/r$ potential converges very slowly or not at all in
real space, and likewise for its Fourier transform
$v(k)=\int\text{d}^{2}\textbf{r}[v(r)e^{-i\textbf{k}\cdot\textbf{r}}]=2\pi/k$
in reciprocal space. Consider the interlayer potential
$v(r)=1/\sqrt{d^{2}+r^{2}}$. Its Fourier transform $v(k)=2\pi e^{-kd}/k$
converges for finite $d$, but convergence becomes very slow for small $d$,
rending direct summation in reciprocal space inefficient. To solve this
problem, we split the potential into a short-range part and a long-range part
as shown below:

$\begin{split}v_{\text{short}}(r)&=\frac{1-\erf{(\frac{\sqrt{d^{2}+r^{2}}}{r_{0}})}}{\sqrt{d^{2}+r^{2}}},\\\
v_{\text{long}}(r)&=\frac{\erf{(\frac{\sqrt{d^{2}+r^{2}}}{r_{0}})}}{\sqrt{d^{2}+r^{2}}},\\\
\end{split}$ (S1)

where $r_{0}$ is a cutoff parameter that we may choose, and $r$ always denotes
the in-plane distance. It is very important that the argument of the error
function involves $\sqrt{d^{2}+r^{2}}$ instead of just $r$, because otherwise
the long-range potential would have “fine-detail” on the scale of $d$ as $r$
approaches $0$, rending its Fourier transform slowly convergent. The short-
range potential can be summed directly in real-space without any additional
trouble by summing over all interactions between a unit cell and its
neighbors, where we choose the neighborhood’s size based on the value of
$r_{0}$ relative to the unit cell dimensions and our desired accuracy.

Unlike the standard three-dimensional case, the Fourier transform of the long-
range potential does not have a nice closed-form expression (Mathematica was
unable to compute the Fourier transform analytically), so it is easiest to
approximate it numerically as the difference between the original potential’s
Fourier transform and the short-range potential’s Fourier transform:

$\begin{split}v_{\text{long}}(k)&=v(k)-v_{\text{short}}(k),\\\
&=\frac{2\pi}{k}e^{-kd}-v_{\text{short}}(k).\\\ \end{split}$ (S2)

The short-range potential’s Fourier transform can be expressed as:

$\begin{split}v_{\text{short}}(k)&=\int_{0}^{\infty}r\text{d}r\bigg{[}\frac{1-\erf{(\frac{\sqrt{d^{2}+r^{2}}}{r_{0}})}}{\sqrt{d^{2}+r^{2}}}\int_{0}^{2\pi}\text{d}\theta\big{[}e^{ikr\cos{\theta}}\big{]}\bigg{]},\\\
&=\int_{0}^{\infty}\text{d}r\bigg{[}\frac{1-\erf{(\frac{\sqrt{d^{2}+r^{2}}}{r_{0}})}}{\sqrt{d^{2}+r^{2}}}2\pi
rJ_{0}(kr)\bigg{]},\\\ \end{split}$ (S3)

where we have used the integral representation of the Bessel functions of the
first kind. Because the error function suppresses the integrand for large $r$,
the upper limit can be replaced with $6r_{0}$ ($6$ was chosen to make the
error from cutting off the integral negligible) and the resulting proper
integral evaluated numerically. To subtract the contribution of the uniform
background charge from the short-range sum, one more integral is needed:

$\begin{split}V_{\text{back}}&=\int\text{d}^{2}\textbf{r}\bigg{[}v_{\text{short}}(r)\bigg{]},\\\
&=\int_{0}^{\infty}\text{d}r\bigg{[}2\pi
r\frac{1-\erf{(\frac{\sqrt{d^{2}+r^{2}}}{r_{0}})}}{\sqrt{d^{2}+r^{2}}}\bigg{]},\\\
&=2\pi\bigg{[}\frac{r_{0}e^{-\frac{d^{2}}{r_{0}^{2}}}}{\sqrt{\pi}}-d\bigg{(}1-\erf{\frac{d}{r_{0}}}\bigg{)}\bigg{]},\end{split}$
(S4)

where the integral has been evaluated analytically using a computer algebra
system.

Figure S1: A local minima found by classical energy minimization with $a/d=2$
and $64$ holes is shown on the left. Because $a/d=2$ is denser than the
crossover $a/d=5.42$, the local minima is comprised predominantly of the
checkerboard composite crystal, although there are various defects present. On
the right, we show a classical “trion crystal” with balanced forces, where
there is clear grouping into linear $e-h-e$ trions. Although stable against
certain directions, such as a uniform rotation of all trions, this
configuration is unstable against arbitrarily small symmetry-breaking
perturbations.

The cutoff parameter $r_{0}$ was chosen to always be one-quarter of the
smallest unit cell dimension, and both the real-space and reciprocal-space
sums were converged to around eight decimal places. For our exploratory
calculations, we scanned through various $a/d$ and placed $128$ electrons and
$64$ holes in a large supercell to avoid artificially imposing any symmetry on
the result. For each value of $a/d$, we ran several calculations with
different random initializations to ensure that we discovered all possible
phases. An example from one of these calculations is shown in Fig. S1. Once
the honeycomb and checkerboard composite crystals were discovered, we used the
defect-free geometries to determine the structures’ precise energies.

One may wonder whether a lattice of classical trions may be stable. As with
the single trion, we were able to find force-balanced configurations that are
periodic arrays of trions, as shown in Fig. S1. However, even the slightest
symmetry-breaking perturbation destroys the classical trion crystal, and the
system relaxes into a composite crystal of electrons and dipoles, indicating
that the classical trion crystal is a saddle point and not a local minima of
the classical energy. Indeed, the trion crystal relaxes to the composite
crystal through the same instability that the single trion has, where one of
the electrons approaches the hole and binds with it to form a dipole while the
remaining electron leaves (shown in Fig. S2). These results are consistent
with the total absence of trion configurations in our randomly initialized
calculations.

Figure S2: Arrangement of electrons (red) and holes (blue) in the exciton
(top) and trion (bottom) configurations. The exciton is the unique ground
state of the system with energy $-1/d$ compared to the trion’s $-0.937/d$.
Additionally, the trion has an imaginary frequency normal mode shown by the
yellow arrows, making it a saddle point of the energy instead of a local
minima.

##  II Hartree-Fock Derivation

To study the quantum phases of the electron-hole bilayer, we introduce a mean-
field approximation that generalizes the Slater determinant of Hartree-Fock
theory. Consider a system of fermions of various species $1\dots f$, where
fermions of the same species are identical. For example, species could
distinguish between different spin polarizations, between electrons and holes,
or between valley degrees of freedom. A wavefunction describing this system
should be antisymmetric under the exchange of two same-species fermions, but
this condition does not apply to the exchange of different-species fermions,
as such an exchange does not really make sense. A simple ansatz for this
system is a product of Slater determinants, one for each species. Letting
$\mathbf{r}_{i}^{\alpha}$ denote the position of the $i$th fermion of species
$\alpha$ (throughout this paper, we use the Greek letters starting with
$\alpha$ and $\beta$ to index species) and $N_{\alpha}$ denote the number of
fermions of species $\alpha$, we can write this anstaz as:

$\displaystyle\psi(\mathbf{r}^{1}_{1},\ldots,\mathbf{r}^{1}_{N_{1}};\mathbf{r}^{2}_{1},\ldots,\mathbf{r}^{2}_{N_{2}};\ldots;\mathbf{r}^{f}_{1},\ldots,\mathbf{r}^{f}_{N_{f}})$
(S5)
$\displaystyle=\prod_{\alpha}\frac{1}{\sqrt{N_{\alpha}!}}\begin{vmatrix}\phi^{\alpha}_{1}(\mathbf{r}^{\alpha}_{1})&\phi^{\alpha}_{2}(\mathbf{r}^{\alpha}_{1})&\cdots&\phi^{\alpha}_{N_{\alpha}}(\mathbf{r}^{\alpha}_{1})\\\
\phi^{\alpha}_{1}(\mathbf{r}^{\alpha}_{2})&\phi^{\alpha}_{2}(\mathbf{r}^{\alpha}_{2})&\cdots&\phi^{\alpha}_{N_{\alpha}}(\mathbf{r}^{\alpha}_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\phi^{\alpha}_{1}(\mathbf{r}^{\alpha}_{N_{\alpha}})&\phi^{\alpha}_{2}(\mathbf{r}^{\alpha}_{N_{\alpha}})&\cdots&\phi^{\alpha}_{N_{\alpha}}(\mathbf{r}^{\alpha}_{N_{\alpha}})\end{vmatrix},$

where the $\phi^{\alpha}_{1\cdots N_{\alpha}}$’s are orthonormal single-
particle orbitals (orthonormality is only enforced between orbitals belonging
to the same species; there is no need for conditions on orbitals belonging to
different species). By comparing to the single-determinant ansatz, we can see
that our multi-determinant ansatz is indeed normalized. While it is true that
such a product of Slater determinants cannot capture correlations between
particles of the same or different species, this ansatz is a good place to
start for a first analysis, and we expect it to be adequate for describing
many strong-coupling phases that we are interested in.

To describe the most general system of $f$ interacting species, let species
$\alpha$ experience a single-particle Hamiltonian $h^{\alpha}$, and let
particles of species $\alpha$ and $\beta$ experience an interaction potential
$v^{\alpha\beta}$. The many-body Hamiltonian is:

$H=\sum_{\alpha}\sum_{i=1}^{N_{\alpha}}h^{\alpha}(\mathbf{r}_{i}^{\alpha})+\frac{1}{2}\sum_{\alpha\beta}\sum_{i=1}^{N_{\alpha}}\sum_{j=1}^{N_{\beta}}v^{\alpha\beta}(\mathbf{r}_{i}^{\alpha}-\mathbf{r}_{j}^{\beta}).$
(S6)

Technically, we should exclude $i=j$ in the sum when $\alpha=\beta$, but this
self-interaction will be exactly cancelled by the exchange energy, so we need
not worry about it. Following the same procedure used to calculate the
expected energy of a single Slater determinant, we can show that the expected
energy of our ansatz is:

$E=\sum_{\alpha}\sum_{i=1}^{N_{\alpha}}\langle\phi^{\alpha}_{i}|h^{\alpha}|\phi^{\alpha}_{i}\rangle+\frac{1}{2}\sum_{\alpha\beta}\sum_{i=1}^{N_{\alpha}}\sum_{j=1}^{N_{\beta}}\bigg{[}\langle\phi^{\alpha}_{i}\phi^{\beta}_{j}|v^{\alpha\beta}(\mathbf{r}_{1}^{\alpha}-\mathbf{r}_{2}^{\beta})|\phi^{\alpha}_{i}\phi^{\beta}_{j}\rangle-\delta^{\alpha\beta}\langle\phi^{\alpha}_{i}\phi^{\beta}_{j}|v^{\alpha\beta}(\mathbf{r}_{1}^{\alpha}-\mathbf{r}_{2}^{\beta})|\phi^{\beta}_{j}\phi^{\alpha}_{i}\rangle\bigg{]}.$
(S7)

Note that the second term is a generalization of the exchange interaction
between electrons of the same spin and is only between fermions of the same
species, while the first term in the two-particle sector generalizes the
Coulomb interaction and is between all species. To minimize this energy, we
express all orbitals as a linear combination of basis functions. Although we
could in principle use a different basis set for each species, it is simpler
to use a common plane-wave basis, which we will index with Greek letters
starting with $\mu$ and $\nu$. Defining the coefficient matrix as
$|\phi^{\alpha}_{i}\rangle=C^{\alpha}_{\mu i}|\mu\rangle$ and matrix elements
$h^{\alpha}_{\mu\nu}=\langle\mu|h^{\alpha}|\nu\rangle$ and
$v^{\alpha\beta}_{\mu\nu\kappa\lambda}=\langle\mu\nu|v^{\alpha\beta}|\kappa\lambda\rangle$,
we can write the Hartree-Fock energy expectation value as:

$E=\sum_{\alpha}h^{\alpha}_{\mu\nu}C^{\alpha*}_{\mu i}C^{\alpha}_{\nu
i}+\frac{1}{2}\sum_{\alpha\beta}(v^{\alpha\beta}_{\mu\kappa\nu\lambda}-\delta^{\alpha\beta}v^{\alpha\beta}_{\mu\kappa\lambda\nu})C^{\alpha*}_{\mu
i}C^{\alpha}_{\nu i}C^{\beta*}_{\kappa j}C^{\beta}_{\lambda j},$ (S8)

where Einstein summation is implied for all indices except the species indices
$\alpha$ and $\beta$. Defining the one-particle reduced density matrix as
$D^{\alpha}_{\mu\nu}=C^{\alpha*}_{\mu i}C^{\alpha}_{\nu i}$, we can further
simplify the energy expression into:

$\displaystyle
E=\bigg{[}\sum_{\alpha}h^{\alpha}_{\mu\nu}+\frac{1}{2}\sum_{\alpha\beta}\bigg{(}v^{\alpha\beta}_{\mu\kappa\nu\lambda}-\delta^{\alpha\beta}v^{\alpha\beta}_{\mu\kappa\lambda\nu}\bigg{)}D^{\beta}_{\kappa\lambda}\bigg{]}D^{\alpha}_{\mu\nu}.$
(S9)

In what follows, we demonstrate a way of directly minimizing the Hartree-Fock
energy functional using unconstrained minimization techniques instead of the
standard self-consistent field (SCF) method Arias et al. (1992). By
analytically continuing the energy functional onto the space of non-
orthonormal orbitals using Lowdin symmetric orthogonalization, no constraints
need to be enforced and no Lagrange multipliers are necessary. Although not
commonly used, direct minimization is implemented in some density functional
theory codes, and it is equivalent to the usual way of solving the Hartree-
Fock equations self-consistently Sundararaman et al. (2017). We chose this
approach because it is immune from instabilities such as charge sloshing and
is guaranteed to converge from any starting point, which is optimal for
exploratory calculations.

If the density matrix is varied by $\Delta D^{\alpha}_{\mu\nu}$, the energy
changes by an amount:

$\displaystyle\Delta
E=\sum_{\alpha}\Tr[\mathbf{F}^{\alpha}(\Delta\mathbf{D}^{\alpha})^{\mathsf{T}}],$
(S10) $\displaystyle
F^{\alpha}_{\mu\nu}[D]=h^{\alpha}_{\mu\nu}+\sum_{\beta}\bigg{(}v^{\alpha\beta}_{\mu\kappa\nu\lambda}-\delta^{\alpha\beta}v^{\alpha\beta}_{\mu\kappa\lambda\nu}\bigg{)}D^{\beta}_{\kappa\lambda}.$

In the top expression, bold font denotes that an object should be treated as a
matrix. Ordinarily, we would minimize the energy under the orthonormality
constraint using Lagrange multipliers. However, there is a way out that allows
us to convert the problem into unconstrained minimization. For any arbitrary
orbitals $|i\rangle$, we can perform Lowdin symmetric orthogonalization to
define the orthonormal orbitals $|\tilde{i}\rangle=|j\rangle(S^{-1/2})_{ji}$,
where $S^{-1/2}$ denotes the matrix inverse-square-root of $S_{ij}=\langle
i|j\rangle$. Using this procedure, the density matrix corresponding to an
arbitrary $C^{\alpha}_{\mu i}$ after orthogonalization is:

$D^{\alpha}_{\mu\nu}=C^{\alpha*}_{\mu j}C^{\alpha}_{\nu
k}(S^{\alpha,-1})_{kj},$ (S11)

where $\textbf{S}^{\alpha,-1}$ is the matrix inverse of
$S^{\alpha}_{ij}=C^{\alpha*}_{\mu i}C^{\alpha}_{\mu j}$. By first
orthogonalizing before plugging into the energy expression, the Hartree-Fock
energy can be analytically continued onto the space of non-orthonormal
orbitals, enabling unconstrained minimization to be used. Note that the matrix
inverse satisfies the relation
$\Delta\mathbf{S}^{-1}=-\mathbf{S}^{-1}(\Delta\mathbf{S})\mathbf{S}^{-1}$, and
the overlap matrix can expressed as
$\mathbf{S}=\textbf{C}^{\dagger}\textbf{C}$. Therefore, if we vary
$\mathbf{C}^{\dagger}$ by $\Delta\mathbf{C}^{\dagger}$ (we treat variables and
their complex conjugates as independent here, which we will justify later),
the change in the density matrix is (we have suppressed the species index for
brevity so everything here refers to a single species):

$\displaystyle\Delta\mathbf{D}^{\mathsf{T}}$
$\displaystyle=\Delta(\mathbf{C}\mathbf{S}^{-1}\mathbf{C}^{\dagger}),$ (S12)
$\displaystyle=\textbf{C}\Delta(\textbf{S}^{-1})\textbf{C}^{\dagger}+\textbf{C}\textbf{S}^{-1}\Delta\textbf{C}^{\dagger}$
$\displaystyle=-\textbf{C}\mathbf{S}^{-1}(\Delta\mathbf{S})\mathbf{S}^{-1}\textbf{C}^{\dagger}+\textbf{C}\textbf{S}^{-1}\Delta\textbf{C}^{\dagger},$
$\displaystyle=\textbf{C}\textbf{S}^{-1}\Delta\textbf{C}^{\dagger}(\textbf{I}-\textbf{C}\textbf{S}^{-1}\textbf{C}^{\dagger}).$

Substituting this into (S10) yields the variation in the energy:

$\displaystyle\Delta E$
$\displaystyle=\sum_{\alpha}\Tr[\mathbf{F}^{\alpha}\textbf{C}^{\alpha}\textbf{S}^{\alpha,-1}\Delta\textbf{C}^{\alpha,\dagger}(\textbf{I}-\textbf{C}^{\alpha}\textbf{S}^{\alpha,-1}\textbf{C}^{\alpha,\dagger})],$
(S13)
$\displaystyle=\sum_{\alpha}\Tr[\Delta\textbf{C}^{\alpha,\dagger}(\textbf{I}-\textbf{C}^{\alpha}\textbf{S}^{\alpha,-1}\textbf{C}^{\alpha,\dagger})\mathbf{F}^{\alpha}\textbf{C}^{\alpha}\textbf{S}^{\alpha,-1}].$

Therefore, we have the derivative of the analytically continued energy with
respect to $\textbf{C}^{\dagger}$:

$\frac{\partial E}{\partial
C^{\alpha*}_{\mu\nu}}=[(\textbf{I}-\textbf{C}^{\alpha}\textbf{S}^{\alpha,-1}\textbf{C}^{\alpha,\dagger})\mathbf{F}^{\alpha}\textbf{C}^{\alpha}\textbf{S}^{\alpha,-1}]_{\mu\nu}.$
(S14)

Using the expression for the Fock matrix $\textbf{F}^{\alpha}$, the energy
expression Eq. S9 can be rewritten to

$E=\frac{1}{2}\sum_{\alpha}(h^{\alpha}_{\mu\nu}+F^{\alpha}_{\mu\nu})D^{\alpha}_{\mu\nu},$
(S15)

allowing the computationally-expensive Fock matrix to be reused in both energy
and gradient computations. In general, if a real-valued function $g$ depends
on the complex-valued argument $z$, we have the following identities:

$\displaystyle\frac{\partial g}{\partial\real[z]}$
$\displaystyle=2\real[\frac{\partial{g}}{\partial z^{*}}],$ (S16)
$\displaystyle\frac{\partial g}{\partial\imaginary[z]}$
$\displaystyle=2\imaginary[\frac{\partial{g}}{\partial z^{*}}],$

where derivatives with respect to the real (imaginary) part of $z$ hold the
imaginary (real) part constant, and derivatives with respect to $z^{*}$ treat
$z$ and $z^{*}$ as independent. For the purposes of updating variables during
some optimization algorithm such as conjugate gradient, one can pretend that
$2\frac{\partial g}{\partial z^{*}}$ is the derivative of $g$ with respect to
$z$, as this produces the same variable update as treating the real and
imaginary parts of $z$ as separate independent variables. Therefore, because
we have computed both the value and derivative of the analytically continued
energy, off-the-shelf minimization algorithms such as conjugate gradient can
be used to solve for the Hartree-Fock orbitals.

For our calculations, we use a large $L_{x}$ by $L_{y}$ rectangular box with
periodic boundary conditions, which is equivalent to placing our system on a
torus. Implementing Hartree-Fock in a large periodic box is essentially
equivalent to enforcing a unit cell and sampling on a particular mesh in the
first Brillouin zone, where using a larger periodic box is equivalent to using
a finer mesh, all else being constant. As with using direct minimization,
using a large periodic box is very helpful when exploring a new system with
unknown symmetry. Our basis set is a rectangular grid of plane-waves
consistent with the periodic boundary conditions:
$\phi_{\mathbf{k}}(\mathbf{r})=e^{i\mathbf{k}\cdot\mathbf{r}}/\sqrt{L_{x}L_{y}}$,
$\mathbf{k}=(2\pi n_{x}/L_{x},2\pi n_{y}/L_{y})$ for integer $n_{x}$ and
$n_{y}$. The one-body contribution to the Hamiltonian is just the kinetic
energy and is diagonal in this basis:

$h^{\alpha}_{\mu\nu}=\delta_{\mu\nu}\frac{1}{2m_{\alpha}}k_{\mu}^{2},$ (S17)

where $\mathbf{k}_{\mu}$ is the wavevector corresponding to the plane-wave
with the flattened index $\mu$, and $m_{\alpha}$ is a species-dependent
effective mass. For a general interaction potential
$v(\mathbf{r}_{1}-\mathbf{r}_{2})$, the two-body matrix element is:

$\displaystyle\langle\mathbf{k}_{1}^{\prime}\mathbf{k}_{2}^{\prime}|v(\mathbf{r}_{1}-\mathbf{r}_{2})|\mathbf{k}_{1}\mathbf{k}_{2}\rangle=$
$\displaystyle\frac{1}{(L_{x}L_{y})^{2}}\int_{\text{box}}\text{d}^{2}\textbf{r}_{1}\int_{\text{box}}\text{d}^{2}\textbf{r}_{2}\bigg{[}v(\mathbf{r}_{1}-\mathbf{r}_{2})e^{-i\textbf{k}_{1}^{\prime}\cdot\textbf{r}_{1}-i\textbf{k}_{2}^{\prime}\cdot\textbf{r}_{2}+i\textbf{k}_{1}\cdot\textbf{r}_{1}+i\textbf{k}_{2}\cdot\textbf{r}_{2})}\bigg{]},$
(S18) $\displaystyle=$
$\displaystyle\frac{1}{(L_{x}L_{y})^{2}}\int_{\text{box}}\text{d}^{2}\textbf{r}_{1}\bigg{[}e^{i(-\textbf{k}_{1}^{\prime}-\textbf{k}_{2}^{\prime}+\textbf{k}_{1}+\textbf{k}_{2})\cdot\textbf{r}_{1}}\bigg{]}\int_{\text{box}}\text{d}^{2}(\textbf{r}_{2}-\textbf{r}_{1})\bigg{[}v(\mathbf{r}_{1}-\mathbf{r}_{2})e^{i(-\textbf{k}_{2}^{\prime}+\textbf{k}_{2})\cdot(\textbf{r}_{2}-\textbf{r}_{1})}\bigg{]},$
$\displaystyle=$
$\displaystyle\frac{\delta_{\textbf{k}_{1}^{\prime}+\textbf{k}_{2}^{\prime},\textbf{k}_{1}+\textbf{k}_{2}}}{L_{x}L_{y}}\int_{\text{box}}\text{d}^{2}\textbf{r}\bigg{[}v(\textbf{r})e^{i(-\textbf{k}_{2}^{\prime}+\textbf{k}_{2})\cdot\textbf{r}}\bigg{]},$
$\displaystyle=$
$\displaystyle\frac{\delta_{\textbf{k}_{1}^{\prime}+\textbf{k}_{2}^{\prime},\textbf{k}_{1}+\textbf{k}_{2}}}{L_{x}L_{y}}v(\textbf{k}_{2}^{\prime}-\textbf{k}_{2}),$

where all wavevectors are consistent with the periodic boundary conditions, we
have used the periodic boundary conditions to simplify the integral’s limits,
and
$v(\textbf{q})=\int\text{d}^{2}\textbf{r}[v(\textbf{r})e^{-i\textbf{k}\cdot\textbf{r}}]$
is the Fourier transform of the interaction potential. We can see that the
matrix element is only nonzero between initial and final two-particle states
of the same momentum and only depends on the magnitude of the momentum
transfer. As mentioned in the section discussing Ewald summation, the Fourier
transform of $v(\textbf{r})=1/\sqrt{r^{2}+d^{2}}$ is $v(\textbf{q})=2\pi
e^{-qd}/q$, which handles both the intralayer and interlayer interactions. The
Fourier transform is divergent at $\textbf{q}=\textbf{0}$ because our
interactions are long-ranged, but this is not an issue because the
contribution of $v(\textbf{q}=\textbf{0})$ to the Fock matrix is exactly
cancelled by the contribution of a neutralizing background charge. To see
this, we explicitly calculate the contribution of
$v^{\alpha\beta}(\textbf{q}=\textbf{0})$ to the Fock matrix:

$\Delta
F^{\alpha}_{\mu\nu}=\delta_{\mu\nu}\sum_{\beta}v^{\alpha\beta}(\textbf{q}=\textbf{0})\frac{D^{\beta}_{\kappa\kappa}}{L_{x}L_{y}}-v^{\alpha\alpha}(\mathbf{q}=\mathbf{0})\frac{D^{\alpha}_{\nu\mu}}{L_{x}L_{y}}.$
(S19)

The potential contributed by a uniform background charge $\rho^{\beta}$ of
species $\beta$ is simply
$\int\text{d}^{2}\textbf{r}[v^{\alpha\beta}(\textbf{r})\rho^{\beta}]=v^{\alpha\beta}(\textbf{q}=\textbf{0})\rho^{\beta}$.
The trace of the density matrix $D^{\beta}_{\kappa\kappa}$ is exactly the
number of particles of species $\beta$, and $L_{x}L_{y}$ is the system area,
so the first term in Eq. S19 is exactly canceled by incorporating a
neutralizing background charge for every species. The contribution to the
exchange term is slightly more subtle. The total exchange matrix is

$\displaystyle K^{\alpha}_{\mu\nu}$
$\displaystyle=v^{\alpha\alpha}_{\mu\kappa\lambda\nu}D^{\alpha}_{\kappa\lambda},$
(S20)
$\displaystyle=\frac{1}{L_{x}L_{y}}\sum_{\textbf{q}}v^{\alpha\alpha}(\textbf{q})D^{\alpha}_{\nu+q,\mu+q},$

where $\mu+q$ is shorthand for the plane-wave state corresponding to the
wavevector $\textbf{k}_{\mu}+\textbf{q}$. For a finite system size, the
divergent $v(\textbf{q}=\textbf{0})$ is indeed problematic at face value.
However, after replacing sum with integral in the infinite system size limit,
the measure $\text{d}^{2}\textbf{q}=2\pi q\text{d}q$ contributes a factor of
$q$ which cancels the singularity coming from the potential
$v(\textbf{q})\propto 1/q$. Cutting out the $\textbf{q}=\textbf{0}$
singularity simply results in a small finite-size error $\propto 1/L_{x}L_{y}$
in the exchange matrix. Therefore, it is completely safe simply to set
$v(\textbf{q}=\textbf{0})=0$, which both removes the singularity and
implements the neutralizing background charge.

Figure S3: A local minima found by a randomly initialized HF run at the same
parameters as in Fig. 3. This is the lowest energy state our randomly
initialized runs found, with an energy per particle higher than the trion-
initialized calculation by only $0.2\%$. Both the electrons and holes have
crystallized onto a nearly perfect triangular lattice, with some defects
visible on the right edge of the box. Additionally, the electrons are nearly
in a pure singlet, with a visually imperceptible difference between the spin
$\uparrow$ and $\downarrow$ densities.

For our calculations, we used a periodic rectangle with a height-to-width
ratio of $\sqrt{3}/2$ to allow a triangular lattice to fit cleanly. We also
performed randomly initialized runs with differently-shaped boxes and
different numbers of particles, but the results were qualitatively the same.
Our basis set was a $45$ by $39$ grid of plane waves consistent with the
periodic boundary conditions. We choose the specific numbers $45$ and $39$
because $39/45$ and $\sqrt{3}/2$ differ only by $0.07\%$, so the resolution in
each direction is essentially the same. The maximum occupancy of the plane-
waves on the borders of the basis grid is on the order of $10^{-8}$ times
smaller than the $k=0$ plane-wave, demonstrating that our basis is very well
converged. We wrote all of our code in Python and minimized the Hartree-Fock
energy using SciPy’s conjugate gradient with default settings (halt
optimization when maximum component of gradient is less than $10^{-5}$). To
provide further justification that the singlet trion crystal is the ground
state for the parameters discussed in the text, Fig. S3 presents the lowest
energy randomly initialized HF calculation we obtained at identical parameters
to that present in the main text.


Non-Existence of S-Integrable Three-Point Partial Difference Equations in the
Lattice Plane

Non-Existence of S-Integrable Three-Point Partial
Difference Equations in the Lattice Plane††This paper is a contribution to the
Special Issue on Symmetry, Invariants, and their Applications in honor of
Peter J. Olver. The full collection is available at
https://www.emis.de/journals/SIGMA/Olver.html

Decio LEVI a and Miguel A. RODRÍGUEZ b

D. Levi and M.A. Rodríguez

a) Mathematical and Physical Department, Roma Tre University,
a) Via della Vasca Navale, 84, I00146 Roma, Italy

b) Departamento de Física Teórica, Universidad Complutense de Madrid,
b) Pza. de las Ciencias, 1, 28040 Madrid, Spain<EMAIL_ADDRESS>
Received April 17, 2023, in final form October 23, 2023; Published online
November 01, 2023

Determining if an $(1+1)$-differential-difference equation is integrable or
not (in the sense of possessing an infinite number of symmetries) can be
reduced to the study of the dependence of the equation on the lattice points,
according to Yamilov’s theorem. We shall apply this result to a class of
differential-difference equations obtained as partial continuous limits of
3-points difference equations in the plane and conclude that they cannot be
integrable.

difference equations; integrability; Yamilov’s theorem

39A14; 39A36

This article is our contribution to the celebration
of Peter Olver’s 70th birthday. We want to express our admiration and
gratitude for his work and friendship.

##  1 Introduction

Partial difference equations have always played an important role in physics
and this has been noticeable the last decades. From one side, discrete systems
seem to be at the base of many important laws of physics (as in quantum
gravity [23]) and on the other side, with the increasing use of computers,
discretizations are playing a growing role in physical applications to
numerically solve differential equations [9] preserving some of the main
properties of these equations, in particular their symmetries [14].

The number of points involved in the equation characterizes these equations,
for instance, partial difference equations on four points in a plane have been
studied at length [3, 8, 13, 15], being the simplest class such that the
evolution can be invertible. The continuous limit corresponds to hyperbolic
partial differential equations. We can distinguish two classes of integrable
equations, those which are linearizable, C-integrable equations in Calogero’s
terminology [4], and equations whose integrability requires the solution of a
scattering problem (characterized by the compatibility of two linear problems
for an auxiliary wave function), S-integrable equations in Calogero’s
terminology.

Partial difference equations defined on three points, which appear as
particular situations in triangular lattices, have not received such detailed
attention (see [20] for a thoroughly study on the linearizable case) although
triangularization is a key tool in differential geometry and they allow to
define many discrete models [2, 7, 17, 18, 19]. Often, however, in the
references mentioned before, the partial difference equations defined on a
triangular lattice involve more than three lattice points.

Boundary value problems can be defined for equations involving only three
lattice points. In Figure 1, we present just an example of the initial data
and the corresponding set of points (a more complete set of possible
configurations is given in [20]).

$r_{1}$$u_{n+1,m}$$u_{n,m+1}$$u_{n,m}$

$r_{2}$$u_{n+1,m}$$u_{n,m+1}$$u_{n,m}$

$r_{3}$$u_{n+1,m}$$u_{n,m+1}$$u_{n,m}$

Figure 1: An example of points related by an equation defined on three points
and the choice of initial conditions (see the main text for a discussion of
these graphics).

We consider an equation relating the values of the field $u$ at the points,
$(m,n)$, $(m,n+1)$, and $(m+1,n)$ (given three points there is only one
relative position, disregarding orthogonality or step sizes). Generically, if
we write the equation as

$f(u_{n,m},u_{n+1,m},u_{n,m+1})=0,$

we can define the value in a point as a function of the values in the other
two points, for instance,

$u_{n,m+1}=g_{1}(u_{n,m},u_{n+1,m}),$

that is, if we give the initial values along $r_{1}$ (see Figure 1), we could
compute (solving the equation) the values of the field in the upper half
plane. We can (again generically) isolate $u_{n+1,m}$:

$u_{n+1,m}=g_{2}(u_{n,m},u_{n,m+1}).$

Then, giving the initial values along the line $r_{2}$ we obtain the values of
the field in the right half plane. Finally, we can also obtain $u_{n,m}$ as a
function of the values in the two other points:

$u_{n,m}=g_{3}(u_{n,m+1},u_{n+1,m}).$

If we provide the initial values along the line $r_{3}$, the solution of the
equation will be obtained in a left descending staircase. Other possible
schemes can be used for this kind of lattices and equations (see [22] for a
description of different initial-boundary value problems).

Multilinear equations can also be considered as particular cases. In fact,
nonlinear equations defined on three lattice points are simpler than those on
four points. The number of parameters in these equations is less than in the
four points case, since the nonlinearities are at most cubic. As shown in
[20], the classification problem for linearizable multilinear partial
difference equations can be carried out up to the end and provides many
examples of nonlinear three-point partial difference equations.

Let us write a generic equation on three lattice points as

$\mathcal{E}_{n,m}(u_{n,m},u_{n+1,m},u_{n,m+1})=0.$ (1.1)

If the independent variables are written as $x_{n}=nh_{n}$ and $t_{m}=mh_{m}$,
in terms of the discrete indices $n$ and $m$ and of the lattice parameters
$h_{n}$ and $h_{m}$, we can get from (1.1) in the continuous limit, when
$h_{n},h_{m}\rightarrow 0$ and $n,m\rightarrow\infty$ in the usual way,
equations of the form

$u_{t}=f(u)u_{x}+g(u)u_{xx},\qquad u=u(x,t),$

as, for instance, the continuous Burgers equation [11].

The analysis of the C-integrable nonlinear partial differential equations
(linearizable via a transformation) can be found in [5, 6]. In the continuous
case the class of possible transformations is richer than in the discrete one
as we can have transformations of the independent and dependent variables. In
the discrete case we have to restrict ourselves to just the transformations of
the dependent variable unless we consider transformable lattices, i.e.,
partial difference schemes (see [14] and references quoted therein).

Contrary to the case of C-integrable equations, the situation in the case of
S-integrable three-point partial difference equation is not known, up to our
knowledge.

In this sense, Yamilov’s theorem constitutes a useful tool in the discussion
of the integrability of differential-difference equations (see, for instance,
[12] for an example of its application in a particular case).

In fact, the classification of differential, difference or differential-
difference equations with higher order symmetries, the symmetry approach to
integrable equations, is a well known field with many important results (see
[1, 16, 21] as interesting examples of these results for differential
evolution equations). In this article, we will only focus on a negative result
on the existence of integrable systems of a particular type to show the role
of Yamilov’s theorem in this field.

The content of this article is as follows. In Section 2, we review the main
results on the linearization of partial difference equations defined on three
lattice points either by point transformations or by Hopf–Cole transformations
[20]. Then in Section 3, by using a theorem introduced by Yamilov (see [14]),
we show that a partial difference equation defined on three points cannot be
S-integrable. The general result is then confirmed by considering multilinear
equations where the calculations are explicit. Section 4 is devoted to some
concluding remarks and a proposal of future researches.

##  2 C-integrable three-point partial difference equations

Here we review the results on linearizable three-point partial difference
equations presented in [20]. To simplify the presentation, as we will limit
ourselves to autonomous equations, we will not need to index the complex
dependent variable by its lattice point but just by its relative position with
respect to the reference point $n,m$, i.e. $u_{n,m}=u_{00}$. Then, equation
(1.1) can be written as

$\mathcal{E}\left(u_{00},u_{10},u_{01}\right)=0.$ (2.1)

In some cases, we can use an autonomous point transformation (for a
nonconstant function $f$)

$\tilde{u}_{00}=f(u_{00})$ (2.2)

to linearize (2.1), providing the linear equation

$a\tilde{u}_{00}+b\tilde{u}_{10}+c\tilde{u}_{01}+d=0,$

with $a$, $b$, $c$ and $d$ being complex coefficients.

###  2.1 Linearizability

The following theorem, see [20], provides a solution to this problem:

######  Theorem 2.1 ([20, Theorem 2]).

Necessary and sufficient condition for the linearizability by a point
transformation (2.2) of an equation belonging to the class (2.1) is that the
equation can be written in the form

$u_{00}=\mathcal{H}^{-1}(-\beta\mathcal{H}(u_{10})-\alpha\mathcal{H}(u_{01})-\gamma),$

where $\mathcal{H}(x)$ is an arbitrary function of its argument,
$\mathcal{H}^{-1}$ is its inverse function and $\alpha$, $\beta$, $\gamma$ are
arbitrary integration constants. This equation is linearizable by the point
transformation $\tilde{u}_{00}=\mathcal{H}(u_{00})$ to the equation

$\tilde{u}_{00}+\beta\tilde{u}_{10}+\gamma\tilde{u}_{01}=0.$

We can also linearize (2.1) by the Cole–Hopf transformation:

$\tilde{u}_{01}=f\left(u_{00}\right)\tilde{u}_{00},$ (2.3)

which transforms the linear equation

$a\tilde{u}_{00}+b\tilde{u}_{10}+c\tilde{u}_{01}=0,$

into (2.1), where $a$, $b$ and $c$ are complex coefficients with
$(a,b,c)\not=(0,0,0)$. We refer to [20] for details on the conditions over the
equation to be linearizable under this transformation.

###  2.2 Classification of linearizable multilinear equations

We can use the above results to classify multilinear difference equations
depending on three points in a two-dimensional lattice:

$au_{00}+bu_{01}+cu_{10}+du_{00}u_{10}+eu_{00}u_{01}+fu_{10}u_{01}+gu_{00}u_{10}u_{01}+h=0,$
(2.4)

where $a,\ldots,h$ are arbitrary complex parameters. We can state the
following theorem [20]:

######  Theorem 2.2 ([20, Theorem 3]).

Apart from the equations which are Möbius equivalent to a linear equation, the
only equations belonging to the class (2.4) which are linearizable by a point
transformation are, up to a Möbius transformation of the dependent variable
$($eventually composed with an exchange of the independent variables
$n\leftrightarrow m)$, the following three

$\displaystyle u_{00}u_{10}u_{01}-1=0,$ (2.5a) $\displaystyle
u_{00}u_{01}-u_{10}=0,$ (2.5b) $\displaystyle u_{10}u_{01}-u_{00}=0.$ (2.5c)

Equations (2.5) linearize by the transformation $\tilde{u}_{00}=\log u_{00}$
with $\log$ always standing for the principal branch of the complex
logarithmic function, respectively, to the equations $(z\in\mathbb{Z})$:

$\displaystyle\tilde{u}_{00}+\tilde{u}_{10}+\tilde{u}_{01}=2\pi\mathrm{i}z,$
$\displaystyle\tilde{u}_{00}-\tilde{u}_{10}+\tilde{u}_{01}=2\pi\mathrm{i}z,$
$\displaystyle-\tilde{u}_{00}+\tilde{u}_{10}+\tilde{u}_{01}=2\pi\mathrm{i}z.$

Let us now consider the linearization via Hopf–Cole transformations (2.3). As
we classify up to a Möbius transformation, we can be always set
$g_{00}(u_{00})=u_{00}$. We can now state the following theorem:

######  Theorem 2.3 ([20, Theorem 5]).

The class of complex autonomous multilinear discrete equations defined on
three points which is linearizable to a homogeneous linear equation by a
Cole–Hopf transformation $\tilde{u}_{01}=u_{00}\tilde{u}_{00}$, is given, up
to a Möbius transformation of the dependent variable $($eventually composed
with an exchange of the independent variables $n\leftrightarrow m)$, by the
following equation:

$\displaystyle\frac{1+u_{00}}{u_{00}}-\frac{1+u_{01}}{u_{10}}=0.$

This equation linearizes to the equation

$\displaystyle\tilde{u}_{00}+\tilde{u}_{10}+\tilde{u}_{01}=0.$

##  3 Yamilov’s theorem and the S-integrability
of three-point partial difference equations

In [14] (see section “Why the shape of integrable equations on the lattice is
symmetric”), the following theorem is proposed and proved.

######  Theorem 3.1 (Yamilov).

If an equation of the form

$\displaystyle\dot{u}_{n}=f_{n}=f(u_{n+N},u_{n+N-1},\dots,u_{n+M}),\qquad
N\geq M,\qquad\frac{\partial f_{n}}{\partial u_{n+N}}\frac{\partial
f_{n}}{\partial u_{n+M}}\neq 0,$

where $f$ is a smooth enough function of its variables, possesses a
conservation law of order $m$, such that $m>\min(|N|,|M|)$, then

$N=-M,$ (3.1)

and $N\geq 0$.

This theorem states a necessary condition for a differential-difference
equation to be S-integrable based on the fact that while a differential-
difference equation with an infinity of symmetries is either S- or
C-integrable, if it has not a sufficiently high conservation law it cannot be
S-integrable. As the content of this theorem is a necessary condition, an
equation which satisfies it may still be not integrable. However, a
S-integrable equation has to satisfy necessarily this theorem.

Here in the following we will use this theorem to show that, while as we saw
in the previous Section partial difference equations defined on three lattice
points can be C-integrable, they may not be S-integrable.

To do so we carry out the partial continuous limit of (1.1). Equation (1.1) is
symmetric in the exchange of $n$ and $m$ so in the whole generality we can do
the continuous limit when $h_{m}\rightarrow 0$ and $m\rightarrow\infty$ so
that $t=mh_{m}$ remains finite. For convenience, we will call
$h_{m}=\varepsilon$ and $u_{n,m+1}=u_{n}(t+h_{m})=u_{n}(t+\varepsilon)$.
Assuming that $u_{n}(t)$ is an entire function of $t$, we can write

$u_{n}(t+\varepsilon)=u_{n}(t)+\varepsilon\dot{u}_{n}(t)+\frac{1}{2}\varepsilon^{2}\ddot{u}_{n}(t)+\mathcal{O}\big{(}\varepsilon^{3}\big{)}.$

So (1.1) becomes

$\mathcal{E}_{n,m}(u_{n,m},u_{n+1,m},u_{n,m+1})=\mathcal{E}_{n}(\varepsilon,u_{n}(t),u_{n+1}(t),u_{n}(t+\varepsilon)).$
(3.2)

Expanding the last result in (3.2) in $\varepsilon$, assuming that
$\mathcal{E}_{n}(\varepsilon,u_{n}(t),u_{n+1}(t),u_{n}(t+\varepsilon))$ is an
entire function, we have

$\displaystyle\mathcal{E}_{n}(\varepsilon,u_{n}(t),u_{n+1}(t),u_{n}(t+\varepsilon))=\mathcal{E}_{n}^{(0)}(u_{n}(t),u_{n+1}(t))\\!+\\!\varepsilon\mathcal{E}_{n}^{(1)}(u_{n}(t),u_{n+1}(t),\dot{u}_{n}(t))\\!+\\!\mathcal{O}\big{(}\varepsilon^{2}\big{)}.$

By a proper choice of the dependence of (1.1) from $u_{m,n}$, $u_{m,n+1}$ and
$u_{m+1,n}$, we can make $\mathcal{E}_{n}^{(0)}=0$ and then its semi-
continuous limit becomes

$\mathcal{E}_{n}^{(1)}(u_{n}(t),u_{n+1}(t),\dot{u}_{n}(t))=0.$ (3.3)

Equation (3.3) is not in the form of Yamilov’s theorem, since there is no
$u_{n-1}$ dependence in the equation, required in the theorem, and cannot be
reduced to it by any lattice re-parametrization as it depends just on two
points.

As an example let us consider the semi-continuous limit of (2.4) when the
complex parameters $a,\ldots,h$ are taken to be entire functions of
$\varepsilon$, the lattice spacing in the $m$ direction which we can assume to
be a constant along the lattice. We can expand the parameters of (2.4) in
powers of $\varepsilon$ and we have

$a=a^{(0)}+\varepsilon a^{(1)}+\varepsilon^{2}a^{(2)}+\cdots$

and similar expressions for the other parameters in the equation. The request
that the zero order in $\varepsilon$ of (2.4) be zero implies

$a^{(0)}=-b^{(0)},\qquad c^{(0)}=e^{(0)}=g^{(0)}=h^{(0)}=0,\qquad
d^{(0)}=-f^{(0)},$ (3.4)

and at first order in $\varepsilon$ we get

$\displaystyle\big{(}b^{(0)}+f^{(0)}u_{1}\big{)}\dot{u}_{0}+(a^{(1)}+b^{(1)})u_{0}+c^{(1)}u_{1}+e^{(1)}u_{0}^{2}+(d^{(1)}+f^{(1)})u_{0}u_{1}$
$\displaystyle\hphantom{\big{(}b^{(0)}+f^{(0)}u_{1}\big{)}\dot{u}_{0}}{}+g^{(1)}u_{0}^{2}u_{1}+h^{(1)}=0.$
(3.5)

Naturally (3.5) is of the form (3.3) and, by a proper choice of the parameters
(3.4), we can always make $\mathcal{E}_{n}^{(0)}=0$, as required in the
general case. In the notation of Yamilov’s theorem, (3.5) has $N=1$ and $M=0$
and consequently $N\geq 0$ but the condition (3.1) is not satisfied. Thus
(3.5) cannot be S-integrable as it does not satisfy the necessary
S-integrability conditions implied by Theorem 3.1. As (3.5) is not
S-integrable for any choice of the parameters so it will be also the
multilinear partial difference equation (2.4).

##  4 Conclusions

The application of Yamilov’s theorem to the three-points partial difference
equation yields the result that no integrable differential-difference equation
can be constructed by taking the limit of this class of equations.

We have just discussed the limits in the direct approach, taking the
continuous limit in one of the two indices of the equation. We have also
considered skew-limits (through a combination of both indices $n$, $m$),
since, as it is known [10, 12], in some cases one can obtain differential-
difference equations satisfying Yamilov’s theorem using this approach.
However, this is not possible in the three-points case.

We plan to extend the results presented here to the case of four-points
partial difference equations on the lattice plane. In this case we know there
are many S-integrable and C-integrable results. Since the classification of
C-integrable equations in the multilinear case is not complete, it would be
interesting to prove that, at least in this multilinear case, there is a
privileged shape of the equation which might contain S-integrable equations.

### Acknowledgements

We thank the anonymous referees for corrections, useful suggestions, and
constructive criticism that helped improve this article. MAR acknowledges the
support of Universidad Complutense de Madrid (Spain), under grant
G/6400100/3000.

Professor Decio Levi passed away at the time we were working on this article.
I will always miss him as my colleague and dearest friend.

## References

  * [1] Abellanas L., Galindo A., Conserved densities for nonlinear evolution equations. I. Even order case, J. Math. Phys. 20 (1979), 1239–1243.
  * [2] Adler V.E., Legendre transformations on a triangular lattice, Funct. Anal. Appl. 34 (2000), 1–9, arXiv:solv-int/9808016.
  * [3] Adler V.E., Bobenko A.I., Suris Yu.B., Classification of integrable equations on quad-graphs. The consistency approach, Comm. Math. Phys. 233 (2003), 513–543, arXiv:nlin.SI/0202024.
  * [4] Calogero F., Why are certain nonlinear PDEs both widely applicable and integrable?, in What is Integrability?, Springer Ser. Nonlinear Dynam., Springer, Berlin, 1991, 1–62.
  * [5] Calogero F., Eckhaus W., Necessary conditions for integrability of nonlinear PDEs, Inverse Problems 3 (1987), L27–L32.
  * [6] Calogero F., Eckhaus W., Nonlinear evolution equations, rescalings, model PDEs and their integrability. I, Inverse Problems 3 (1987), 229–262.
  * [7] Dynnikov I.A., Novikov S.P., Geometry of the triangle equation on two-manifolds, Mosc. Math. J. 3 (2003), 419–438, arXiv:math-ph/0208041.
  * [8] Garifullin R.N., Yamilov R.I., Generalized symmetry classification of discrete equations of a class depending on twelve parameters, J. Phys. A 45 (2012), 345205, 23 pages, arXiv:1203.4369.
  * [9] Hansen P.C., Discrete inverse problems: Insight and algorithms, Fundamentals of Algorithms, Vol. 7, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2010.
  * [10] Hietarinta J., Joshi N., Nijhoff F.W., Discrete systems and integrability, Cambridge Texts Appl. Math., Cambridge University Press, Cambridge, 2016.
  * [11] Levi D., Ragnisco O., Bruschi M., Continuous and discrete matrix Burgers’ hierarchies, Nuovo Cimento B 74 (1983), 33–51.
  * [12] Levi D., Rodríguez M.A., Yamilov’s theorem for differential and difference equations, Ufa Math. J. 13 (2021), 152–159.
  * [13] Levi D., Scimiterna C., Linearizability of nonlinear equations on a quad-graph by a point, two points and generalized Hopf–Cole transformations, SIGMA 7 (2011), 079, 24 pages, arXiv:1108.3648.
  * [14] Levi D., Winternitz P., Yamilov R.I., Continuous symmetries and integrability of discrete equations, CRM Monogr. Ser., Vol. 38, American Mathematical Society, Providence, RI, 2022.
  * [15] Levi D., Yamilov R.I., Generalized symmetry integrability test for discrete equations on the square lattice, J. Phys. A 44 (2011), 145207, 22 pages, arXiv:1011.0070.
  * [16] Mikhailov A.V., Shabat A.B., Yamilov R.I., The symmetry approach to the classification of nonlinear equations. Complete lists of integrable systems, Russian Math. Surveys 42 (1987), no. 4, 1–65.
  * [17] Moritz B., Schwalm W., Triangle lattice Green functions for vector fields, J. Phys. A 34 (2001), 589–602.
  * [18] Novikov S.P., Discrete connections on the triangulated manifolds and difference linear equations, arXiv:math-ph/0303035.
  * [19] Schwalm W., Moritz B., Giona M., Schwalm M., Vector difference calculus for physical lattice models, Phys. Rev. E 59 (1999), 1217–1233.
  * [20] Scimiterna C., Levi D., Three-point partial difference equations linearizable by local and nonlocal transformations, J. Phys. A 46 (2013), 025205, 13 pages.
  * [21] Svinolupov S.I., Sokolov V.V., Weak nonlocalities in evolution equations, Math. Notes 48 (1990), 1234–1239.
  * [22] van der Kamp P.H., Initial value problems for lattice equations, J. Phys. A 42 (2009), 404019, 16 pages.
  * [23] Williams R.M., Discrete quantum gravity, J. Phys. Conf. Ser. 33 (2006), 38–48.


11institutetext:  CNRS, CC-IN2P3, 21 avenue Pierre de Coubertin, CS70202,
69627 Villeurbanne CEDEX, France 22institutetext: UniversitÃ Clermont
Auvergne, CNRS/IN2P3, Laboratoire de Physique de Clermont, F-63000 Clermont-
Ferrand, France

# Experience deploying an analysis facility for the Rubin Observatoryâs Legacy
Survey of Space and Time (LSST) data

Gabriele Mainetti 11<EMAIL_ADDRESS>Fabio Hernandez 11 Fabrice
Jammes 22 Quentin Le Boulc’h 11

###### Abstract

The Vera C. Rubin Observatory is preparing for the execution of the most
ambitious astronomical survey ever attempted, the Legacy Survey of Space and
Time (LSST). Currently in its final phase of construction in the Andes
mountains in Chile and due to start operations in 2025 for 10 years, its
8.4-meter telescope will nightly scan the southern sky and collect images of
the entire visible sky every 4 nights using a 3.2 Gigapixel camera, the
largest imaging device ever built for astronomy. Automated detection and
classification of celestial objects will be performed by sophisticated
algorithms on high-resolution images to progressively produce an astronomical
catalog eventually composed of 20 billion galaxies and 17 billion stars and
their associated physical properties.

In this paper, we briefly present the infrastructure deployed at the French
Rubin data facility (operated by IN2P3 computing center, CC-IN2P3) to deploy
the Rubin Science Platform, a set of web-based services to provide effective
and convenient access to LSST data for scientific analysis. We describe the
main services of the platform, the components that provide those services and
our deployment model. We also present the Kubernetes-based infrastructure we
are experimenting with for hosting the LSST astronomical catalog, a petabyte-
scale relational database developed for the specific needs of the project.

##  1 Introduction

This article outlines the tools being deployed by the CNRS/IN2P3 Computing
Center111https://cc.in2p3.fr (CC-IN2P3) to provide researchers a convenient
way to explore and analyse the data that will be produced as part of the
Legacy Survey of Space and Time (LSST) by the Vera C. Rubin Observatory.

The Vera C. Rubin Observatory is in the final stages of its construction on
the Chilean Andes. Equipped with an 8.4-meter telescope and the largest
digital camera ever built (3.2 gigapixel), it will execute the most ambitious
survey ever attempted: the Legacy Survey of Space and Time (LSST) collecting
imagery covering the entire visible sky every four nights.

The four main science goals of LSST are to investigate dark matter and dark
energy, map the Milky Way, take an inventory of the solar system and study the
transient sky Ivezic:2019 . To achieve those goals, the Observatory will
produce, by the end of the survey after 10 years of operation, an astronomical
catalog containing 37 billion objects (20 billion galaxies, 17 billion stars)
as well as a set of 5.5 million images. Science-ready images and the
astronomical catalog will be regularly delivered by the Observatory to science
collaborations.

This paper is structured as follows. The section 2 introduces the role of CC-
IN2P3 as a data facility for the Rubin Observatory. Section 3 describes the
analysis facility, the technology behind it and the resources deployed by the
CC-IN2P3 to run it. The following sections cover the analysis facility main
components and, in particular, section 4 describes the database component and
section 5 the data analysis platform.

##  2 CC-IN2P3 as a Rubin Observatory Data Facility

Image data recorded by the Rubin Observatory will be processed by 3 centers,
one in the USA (US Data Facility - USDF222USDF is located at SLAC National
Accelerator Laboratory, https://www.slac.stanford.edu) and the other two in
Europe (UK Data Facility - UKDF333Supported by the IRIS Infrastructure
https://www.lsst.ac.uk and https://www.iris.ac.uk and France Data Facility -
FrDF), see Figure 1.

Figure 1: Images flow from the Summit Site where the telescope is located in
Chile to the three Rubin Data Facilities which collectively provide the
computational capacity for processing the images taken by the Observatory for
the duration of the survey.

FrDF is located in Lyon, France at the CC-IN2P3 and will provide the computing
and storage resources to annually process 40% of the raw images recorded since
the beginning of the survey, with 35% provided by the USDF and 25% by the
UKDF. In addition, CC-IN2P3 is preparing to provide long-term storage for the
entire raw images data set as well as for selected data products. By the end
of the survey, the size of the Rubin data sets could reach several hundred
petabytes, with 15 PB for the final version of the astronomical catalog alone.

##  3 The analysis facility

The Rubin Observatory’s data analysis facility deployed at CC-IN2P3 is
designed with three main objectives: to provide access to Rubin data, to
integrate the working environment with the rest of the site’s environment and
to reach high levels of scalability and resilience.

###  3.1 Data Access

Our main objective is to provide researchers with an intuitive platform
providing a convenient and effective way to access and to analyse the
extensive survey data (both images and catalog) generated by the Rubin
Observatory. This objective is achieved with two major components of the
platform:

  1. 1.

Qserv: Qserv qserv ; qservrepo  is an open source Massively Parallel
Processing database which hosts the Rubin astronomical catalog (see Section 4
below for details), and

  2. 2.

Rubin Science Platform (RSP): the RSP LSE-319  is an online interactive
analysis hub, allowing researchers to access and analyse the data via a user-
friendly web-based interface, see Section 5 for more details.

###  3.2 Integration with CC-IN2P3 Environment

We want the analysis facility to be seamlessly integrated into the existing
CC-IN2P3 working environment to make easy for the researchers the transition
between the platform and their familiar environment. This integration involves
many aspects, including authentication and transparent access to all the
storage areas individual users can usually access from a terminal session or
from the batch farm, including their $HOME directory and other project-
specific disk-based storage areas. Site-wide single sign on authentication at
CC-IN2P3 is based on Keycloak keycloak .

###  3.3 Scalability and Resilience

The platform is designed to be both scalable and resilient to ensure that it
can adapt to evolving demands linked to the large volumes of astronomical
data. The use of Kubernetes (K8S) k8s  as the orchestration tool for the
platform’s containers allows us to effectively manage and scale the components
as needed.

For this purpose, CC-IN2P3 provides two K8S clusters dedicated to the Qserv
and the Rubin Science Platform: one used as test-bench is based on OpenStack
openstack  virtual machines and another one used for production on top of
bare-metal444The bare metal cluster is composed of DELL PowerEdge R440 and
DELL PowerEdge R540 servers, see Table 1 for a summary.

Table 1: Details of the components of the Kubernetes clusters dedicated to the Rubin analysis platform at CC-IN2P3. Cluster | Production | Test-bench
---|---|---
Infrastructure | Bare-metal | OpenStack
Total Nodes | 25 | 8
Qserv Nodes | 17 | 3
RSP Nodes | 5 | 1
Local Storage | 50 TB | 1 TB

CC-IN2P3 provides also 4 data transfer nodes to share with the community
catalog data ready for ingestion into Qserv via Caddy caddy  HTTP servers;
these servers are also used internally to populate the Qserv databases
deployed at the site (see Section 4.2).

##  4 Qserv

The astronomical catalog produced by LSST will include the physical properties
of 20 billion galaxies and 17 billion stars, resulting in a total of 15
petabytes of data by the end of the survey. The Qserv database management
system qserv ; qservrepo  has been specifically developed by Rubin project
members at SLAC, with contributions from IN2P3, to handle this large volume of
data and to respond to the specific data extraction requests scientists
typically perform against astronomical catalogs.

Qserv is a shared-nothing Massively Parallel Processing Relational (SQL)
Database: this means the processing is split among servers and a leader node
handles communications with each of the individual nodes as a map-reduce
process. Nodes don’t share resources with other nodes. A representation of
this architecture is shown in Figure 2.

Figure 2: Shared-nothing database architecture, from DMTN-243

Developed as an open source software project, it includes data partitioning
over the celestial sphere’s region of equal area, data replication to ensure
resilience and high-availability, shared scanning to reduce the total I/O
load, and sciSQL scisql  User Defined Functions (UDFs) to simplify spherical
geometry, statistics and photometry queries.

###  4.1 Deployment: qserv-operator

Deployment of Qserv is based on Kubernetes operator SDK operator , a framework
to manage and automate the deployment of complex applications. This approach
makes the installation of Qserv really easy because only two commands are
necessary555A demo of the deployment is available at https://is.gd/FK62Wa.:

      $ kubectl apply -f manifest/operator.yaml
      $ kubectl apply -k manifest/<instance>

where <instance> refers to the customization needed for a specific
infrastructure. The qserv-operator qserv-operator  code is available in
GitHub.

###  4.2 Data ingestion: qserv-ingest

Qserv has a powerful distributed data ingestion algorithm integrated: qserv-
ingest.

The ingestion workflow’s schema is shown in Figure 3; qserv-ingest uses Argo
Workflow argowf  to orchestrate the ingestion tasks which concurrently execute
on all nodes in the cluster. The workflow’s input data to be loaded into the
catalog’s tables must be previously partitioned (using a partitioner included
in Qserv) and stored in CSV format.

Figure 3: Schematic view of the data ingestion process into Qserv database.

Each transaction can handle a significant number of files for ingestion.
Asynchronous REST requests are issued to drive the ingestion of each input
file, enabling recovery in case of errors. The workflow also performs
validation steps and executes benchmarks on the ingested database.

During a test campaign in 2022, qserv-ingest was able to ingest 22 million
files for an aggregated 40TB of data in 5 hours. At CC-IN2P3 we use qserv-
ingest as the workflow to ingest data in Qserv.

###  4.3 Status of Qserv at CC-IN2P3

The production Qserv instance at CC-IN2P3 is currently composed of 17 nodes,
15 qserv-worker and 2 qserv-master nodes. There are 5 astronomical catalogs
ingested for an aggregated size of 90 TB, as shown in Table 2.

Table 2: Status of Qserv-managed astronomical catalogs at CC-IN2P3. Catalog | Size (TB) | Number of rows (Billions)
---|---|---
idf-dp0.2-catalog | 36.6 | 139
frdf-dp0.2-catalog | 35.9 | 121
dp01_dc2_catalogs | 1.1 | 1.7
skysim5000_v1.1.1 | 13.6 | 20.5
cosmoDC2_v1.1.4 | 3.7 | 5.5

The DP0.1 and DP0.2 catalogs above are catalogs generated by the LSST data
processing campaigns RTN-001 ; pipecc  using simulated sky images and they are
a preview of what we expect once the Rubin Observatory will start to produce
real data.

##  5 The Rubin Science Platform

The Rubin Science Platform (RSP) is a coherent set of web-based applications
and services designed by the Rubin project’s Data Management team to offer
scientists a unified, near-to-the-data interactive analysis tool. It uses
Firefly firefly  for data visualisation and plotting and provides access to
both Qserv-managed catalogs as well as external catalogs. For advanced,
programmatic analysis using Python it integrates a Jupyter jupyter  notebook
platform configured to include the LSST Science Pipelines lsst-science-
pipelines . Finally, it acts also as gateway to Qserv catalogs for Virtual
Observatory666https://www.ivoa.net Table Access Protocol TAP  compatible
tools.

A view of the RSP is shown in Figure 4.

(a) Web portal access

(b) Tables and images viewer

(c) Jupyter notebook platform

Figure 4: A view of the RSP

An experimental instance of the Rubin Science Platform deployed at CC-IN2P3 is
online and accessible to registered project members at the address
http://data-dev.lsst.eu.

###  5.1 Deployment

Designed as a collection of applications and services, it has been developed
by the Rubin Observatory Science Quality and Reliability
Engineering777https://github.com/lsst-sqre (SQuaRE) team.

The RSP components run on top of Kubernetes clusters and each component is
configured via Helm helm  charts. Individual Helm charts are managed via the
Phalanx phalanx  git repository phalanxrepo .

The use of ArgoCD argocd  allows to synchronize these application deployment
manifests into the Kubernetes cluster of each environment.

Authentication is managed via the gafaelfawr (gafaelfawr , DMTN-224 , DMTN-234
) identity management application, which supports LDAP ldap  and OpenId
Connect opeind  providers (e.g. Github). CC-IN2P3’s RSP instance uses Keycloak
keycloak  as OpenID Connect (OIDC) provider and LDAP as provider of the user’s
group identifiers.

In addition to gafaelfawr, there are 4 core applications providing key
functionality for other applications:

  * •

argocd for deployment orchestration,

  * •

cert-manager for certificates management,

  * •

ingress-nginx for traffic routing,

  * •

vault-secrets-operator for secret management

Seventeen components are currently activated on the CC-IN2P3 RSP instance,
some of them requiring adaptations to the specifics of the CC-IN2P3
environment.

##  6 Summary

The Rubin Observatory will record a large number of images of the Southern sky
over several years. We outlined the role of CC-IN2P3 as one of the three main
data facilities that will process these images to regularly produce an updated
version of an astronomical catalog and we presented how we deployed an
analysis facility on top of Kubernetes, highlighting its two key components:
the scalable, shared and resilient database designed for serving the Rubin
astronomical catalog and the Rubin Science Platform to perform interactive
analysis of the data to be regularly released by the Observatory.

## References

  * (1) Z. Ivezic et al., Astrophys. J. 873, 111 (2019)
  * (2) D.L. Wang et al., _Qserv: A Distributed Shared-Nothing Database for the LSST Catalog_ , in _State of the Practice Reports_ (Association for Computing Machinery, 2011), SC ’11, ISBN 9781450311397, https://doi.org/fnk9nm
  * (3) _Qserv: petascale distributed database_ , https://github.com/lsst/qserv
  * (4) M. JuriÄ, D. Ciardi ,G.P. Dubois-Felsmann and L.P. Guy, _LSE-319 - LSST Science Platform Vision Document_ (2019), Vera C. Rubin Observatory Systems Engineering Technical Note, https://lse-319.lsst.io
  * (5) _KeyCloak_ , https://www.keycloak.org
  * (6) _Kubernetes_ , https://kubernetes.io
  * (7) _OpenStack_ , https://www.openstack.org
  * (8) _Caddy, an open source web server_ , https://caddyserver.com
  * (9) F. Mueller et al., _DMTN-243 - Qserv: A Distributed Petascale Database for the LSST Catalogs_ (2022), Vera C. Rubin Observatory Data Management Technical Note, https://dmtn-243.lsst.io/
  * (10) _sciSQL 0.3: Science tools for mysql_ , https://smonkewitz.github.io/scisql/
  * (11) _Operator Framework for kubernetes_ , https://operatorframework.io
  * (12) _qserv-operator: a qserv operator for Kubernetes_ , https://github.com/lsst/qserv-operator
  * (13) _Argo Workflow a workflow engine for orchestrating parallel jobs on kubernetes_ , https://argoproj.github.io/argo-workflows
  * (14) W. O’Mullane, _RTN-001 - Data Preview 0: Definition and planning._ (2021), Vera C. Rubin Observatory Technical Note, https://rtn-001.lsst.io/
  * (15) Q. Le Boulc’h, F. Hernandez, G. Mainetti, _The Rubin Observatoryâs Legacy Survey of Space and Time DP0.2 processing campaign at CC-IN2P3_ , in _Proc. of CHEP 2023_ (to appear)
  * (16) _Firefly a browser-based interactive particle visualization app_ , http://firefly-viz.com
  * (17) _Jupyter_ , https://jupyter.org
  * (18) J. Bosch et al., _An Overview of the LSST Image Processing Pipelines_ , in _Astronomical Data Analysis Software and Systems XXVII_ , edited by P.J. Teuben, M.W. Pound, B.A. Thomas, E.M. Warner (2019), Vol. 523 of _Astronomical Society of the Pacific Conference Series_ , p. 521, 1812.03248
  * (19) P. Dowler, G. Rixon, D. Tody, M. Demleitner, _Table Access Protocol Version 1.1_ , IVOA Recommendation 27 September 2019 (2019), https://10.5479/ADS/bib/2019ivoa.spec.0927D
  * (20) _Helm the package manager for kubernetes_ , https://helm.sh
  * (21) _Phalanx a gitops repository for rubin observatoryâs kubernetes environments_ , https://phalanx.lsst.io
  * (22) _Phalanx_ , https://github.com/lsst-sqre/phalanx
  * (23) _ArgoCD a declarative, gitops continuous delivery tool for kubernetes_ , https://argo-cd.readthedocs.io
  * (24) _Gafaelfawr rsp identity management_ , https://github.com/lsst-sqre/gafaelfawr
  * (25) R. Allbery, _DMTN-224 - RSP identity management implementation strategy_ (2023), Vera C. Rubin Observatory Data Management Technical Note, https://dmtn-224.lsst.io/
  * (26) R. Allbery, _DMTN-234 - RSP identity management design_ (2023), Vera C. Rubin Observatory Data Management Technical Note, https://dmtn-234.lsst.io/
  * (27) _LDAP_ , https://ldap.com
  * (28) _OpenId_ , https://openid.net