# Magnetic-Field-Induced Wigner Crystallization of Charged Interlayer Excitons
in van der Waals Heterostructures

Igor V. Bondarev<EMAIL_ADDRESS>Department of Mathematics & Physics,
North Carolina Central University, Durham, NC 27707, USA

Yurii E. Lozovik  Institute of Spectroscopy, Russian Academy of Sciences,
142190 Troitsk, Moscow Region, Russia  National Research University ”Higher
School of Economics”,Tikhonov Moscow Institute of Electronics & Mathematics,
123458 Moscow, Russia

###### Abstract

We develop the theory of the magnetic-field-induced Wigner crystallization
effect for charged interlayer excitons (CIE) discovered recently in
transition-metal-dichalcogenide (TMD) heterobilayers Lius-PKim19 . We derive
the ratio of the average potential interaction energy to the average kinetic
energy for the many-particle CIE system subjected to the perpendicular
magnetic field of an arbitrary strength, analyze the weak and strong field
regimes, and discuss the ’cold’ crystallization phase transition for the CIE
system in the strong field regime. We also generalize the effective g-factor
concept previously formulated for interlayer excitons Nagler17 , to include
the formation of CIEs in electrostatically doped TMD heterobilayers. We show
that magnetic-field-induced Wigner crystallization and melting of CIEs can be
observed in strong-field magneto-photoluminescence experiments with TMD
heterobilayes of systematically varied electron-hole doping concentrations.
Our results advance the capabilities of this new family of transdimensional
quantum materials.

Coherent states of semiconductors with high excitation density have been at
the forefront of semiconductor materials research since its inception back in
the early sixties KeldyshKozlov68 ; LozovikYudson ; Ogawa90 ; LozovikPRL07 ;
Kotthaus13 ; Kezer14 ; Fogler14 ; Suris16 ; LozovikUFN18 ; Lozovik19 ;
JonFinley ; FWang20 ; FWang21 ; MakShan21 ; CommPhys-bond . Nowadays, a new
family of ultrathin quasi-2D van der Waals (vdW) materials of controlled
thickness (also called transdimensional materials BoltSha19 ; BondMoSha20 ),
particularly bilayer transition metal dichalcogenide (TMD) heterostructures
MakShan16 , provide unique opportunities to study such states experimentally.
Recent experiments have revealed that not only neutral interlayer excitons
(IE) and biexcitons can form of electrons ($e$) in one layer and holes ($h$)
in another layer of these bilayers due to the dimensionality reduction
Rivera2015 ; Ross2017 ; Baranowski2017 ; Miller17 ; BondVlad18 . Predicted
theoretically BondVlad18 ; Drummond18 ; Thygesen18 ; SnokeBond21  and then
discovered experimentally SnokeBond21 ; Lius-PKim19 ; Geim20 ; Xu21  are also
bound three-particle fermion and four-particle boson $e$-$h$ complexes of
nonzero charge such as charged interlayer excitons (CIE, or interlayer trion,
a charge carrier bound to an IE Lius-PKim19 ; Geim20 ; Xu21 ) and doubly
charged quaternions (two like-charge carriers bound to an IE SnokeBond21 ).
Both of them feature the permanent electric dipole moment directed
perpendicular to the layers and so can be effectively manipulated by an
external electromagnetic field. In this way, the CIEs enable the optical
control of spin-1/2 states they carry, which is advantageous for applications
such as spintronics and quantum information processing, while the
multiparticle collective (Bose-condensed) quaternion states offer no less than
a new path to superconductivity without Cooper pairing due to their nonzero
charge. For TMD bilayers under intense irradiation, therefore, collective
cooperative phenomena such as crystallization and condensation phase
transitions must be studied for these bound few-particle $e$-$h$ complexes
theoretically and experimentally to advance potential capabilities of this new
family of quantum materials.

Figure 1: Schematic of the bilayer TMD structure under study, with negative
(a) and positive (b) charged interlayer exciton complexes (CIE, or interlayer
$X^{\pm}$-trions) subjected to a perpendicular homogeneous magnetostatic
field. The distance $\Delta\rho_{X^{\pm}}$ is between the centers of mass of
the $e$-$h$ pairs forming the CIE.

Recently, there were two scenarios predicted for the long-range ordered
crystal phases of the collective multiparticle CIE system in TMD bilayer
nanostructures in the absence of an external magnetic field CommPhys-bond .
They are the unlike-charge CIE crystallization and the Wigner crystallization
of the like-charge CIEs, which can be selectively realized by fabricating TMD
bilayers with appropriately chosen $e$-$h$ effective mass ratio and interlayer
distance in addition to using the standard electrostatic doping technique. The
unlike-charge CIEs form the periodic 1D chains of the two interpenetrating
sublattices of positive and negative CIEs with collinear permanent dipole
moments in each of the two, whereas the dipole moments of the sublattices are
anticollinear. The like-charge CIEs form the 2D Wigner crystal (WC) lattice
instead provided that temperature is below and the CIE density is above their
respective critical values, with the latter controlled by the interlayer
distance. The WC phase has been one of the longest anticipated exotic
correlated phases originally thought of as a periodic array of electrons held
in place when their Coulomb repulsion energy exceeds the Fermi and thermal
fluctuation energies Wigner34 ; Platzman74 ; Morf79 . Lozovik and Yudson were
the first to predict the magnetic-field-induced crystallization of two-
dimensional (2D) electron gas LozYud75 , having noticed that the strong
perpendicular magnetic field should turn the in-plane translational electron
motion into small-radius circular orbits, which the Coulomb repulsion will
then be ordering on a periodic lattice to minimize the Helmholtz free energy
of the system. The Lindemann’s melting criterion Khrapak20  and the simplest
single-mode (Einstein) model of lattice vibrations were used to derive the
upper bound of $\sim\\!100$ kG (10 T) for the critical magnetic field for the
’cold’ melting ($T\\!=\\!0$ K) of the Wigner solid of 2D electrons LozYud75 .
The more realistic Debye vibrational model, still based on the Lindemann’s
melting criterion, was used to study the $T$-dependence of the stability
regions of the magnetic-field-induced WC MusLoz81 . The model was then
extended to include the lattice next-nearest-neighbor anharmonicity effects
based on the properly modified Lindemann’s melting criterion BedGaLoz85pla ,
whereby the crystal phase shear modulus was shown to suddenly drop at melting
LozFarAbd85  and the melting phase diagram was obtained LozFar85 . The Wigner
lattice instability point was shown to be very close to the point of the
Kosterlitz-Thouless topological transition in 2D systems Kosterlitz73 , in
agreement with both earlier numerical studies Morf79  and most current
thermodynamical interpretation of the Lindemann’s law Khrapak20 . The
microscopic anharmonic theory of melting is confirmed by molecular dynamics
simulations BedGaLoz85jetp  and is consistent with experimental studies
GrimesAdams79 ; GrimesAdams80 ; Mendez ; Kajita ; Monarkha ; Shayegan19 ;
Shayegan20 . Recent experimental measurements of the magnetic-field-induced
Wigner crystallization in semiconductor 2D electron systems can be found in
Refs. Shayegan19 ; Shayegan20 .

Here, we present our theory for the Wigner crystallization of like-charge
CIEs, which is induced by the strong external magnetostatic field directed
perpendicular to the TMD bilayer structure. The sketch of the system is shown
in Fig. 1. With magnetic field treated nonperturbatively, we derive the ratio
of the potential interaction energy to the kinetic energy that controls the
crystallization process for the multiparticle CIE system in the presence of
the external perpendicular magnetostatic field. We analyze the weak and strong
field strength limits and obtain the critical temperature for the ’cold’
crystallization phase transition in the strong field regime where only the
lowest Landau level is occupied. We also explain how the effects of
crystallization and melting can be detected in photoluminescence (PL)
experiments. We show that they can be clearly seen in magneto-PL spectra of
bilayer TMD heterostructures under electrostatic doping. Finally, we remark on
the universality of the Wigner crystallization phenomenon for neutral and
charged quasiparticle excitations of bosonic and fermionic nature.

###  .1 Zero magnetic field like-charge trion Wigner crystallization

In the absence of an external field, an ensemble of repulsively interacting
particles (or quasiparticles, structureless or compound) forms a WC lattice
when its average potential interaction energy $\langle V\rangle$ exceeds
average kinetic energy $\langle K\rangle_{0}$, that is Platzman74

$\frac{\langle V\rangle}{\langle K\rangle_{0}}=\Gamma_{0}>1.$ (1)

For like-charge CIEs (or interlayer trions denoted as $X^{{}^{\pm}}$ in the
following) shown in Fig. 1, the Coulomb repulsion at large $R$
($\gg\\!r_{ee,hh}$) is strengthened at shorter $R$ by the dipole-dipole
repulsion of their collinear permanent dipole moments directed perpendicular
to the heterostructure plane and controlled by the interalayer distance $d$ of
the TMD bilayer.

For trion-trion separation distances $R$ greater than the trion size
($R\\!\gg\\!r_{ee,hh}$ in Fig. 1), the first order power series expansion of
the average repulsive trion-trion interaction potential takes the form

$\langle V\rangle=\frac{1}{R}\Big{(}1+\frac{d^{2}}{R^{2}}\Big{)}=\sqrt{\pi
n}\,\big{(}1+d^{2}\pi n\big{)},$ (2)

where $n\\!=\\!1/\pi R^{2}$ is the trion surface density such that $d^{2}\pi
n\\!\ll 1$ for our series expansion to hold true. With $d\\!\sim\\!3-6\,$Å
typical of vdW systems, one then has $n\\!\ll\\!10^{14}\,$cm-2 which is always
the case in experiments with TMD materials where the densities are at least
two orders of magnitude lower to reduce exciton nonradiative recombination
Wang17 . Also, it was shown recently CommPhys-bond  that
$r_{ee,hh}\\!\lesssim\\!2\;$nm and the actual inter-trion repulsion potential
rapidly decreases with distance due to screening, whereby
$R/r_{ee,hh}\\!\gtrsim\\!3$ and our CIEs are separated in space well enough to
make their description in terms of free particle occupation numbers

$n_{\mathbf{k}}=\frac{1}{e^{\,\beta(E_{\mathbf{k}}-\bar{\mu})}+1}$ (3)

legitimate. Here, $\beta\\!=\\!1/k_{B}T$,
$E_{\mathbf{k}}\\!=\\!\hbar^{2}k^{2}/2M$, $M\\!=\\!M_{X^{\pm}}$ and
$\bar{\mu}$ being the trion total effective mass and chemical potential,
respectively. At zero $T$, this turns into a unit-step function to give $n$ in
Eq. (2) in the form

$n=\frac{\langle
N\rangle}{S}=\frac{2}{S}\sum_{\mathbf{k}}n_{\mathbf{k}}=\frac{2}{S}\frac{S}{(2\pi)^{2}}\,2\pi\\!\\!\int_{0}^{k_{F}}\\!\\!\\!dk\,k=\frac{k_{F}^{2}}{2\pi}$
(4)

with $S$ and $k_{F}$ being the surface area of the sample and the trion Fermi-
momentum, respectively, whereby the translational kinetic energy per particle
takes the form

$\langle K\rangle_{0}=\frac{2}{\langle
N\rangle}\sum_{\mathbf{k}}E_{\mathbf{k}}n_{\mathbf{k}}=\frac{\hbar^{2}\pi
n}{2M}=\frac{\hbar^{2}\big{(}\langle k_{x}^{2}\rangle+\langle
k_{y}^{2}\rangle\big{)}}{2M}=\frac{\langle p_{x}^{2}\rangle+\langle
p_{y}^{2}\rangle}{2M}$ (5)

with $\langle p_{x}^{2}\rangle\\!=\\!\langle
p_{y}^{2}\rangle\\!=\\!\hbar^{2}\pi n/2$. This should be added by the
rotational kinetic energy contribution,
$K^{\texttt{(r)}}_{\\!X^{{}^{\pm\\!}}}\\!=\hbar^{2}l(l+1)/2I_{\\!X^{{}^{\pm}}}$,
in which $l\\!=\\!0,1,2,...$ is the orbital quantum number and
$I_{\\!X^{{}^{\pm}}}\\!=\\!m_{h,e}r_{ee,hh}^{2}/2$ is the moment of inertia
for the CIE rotation about its permanent dipole moment direction (see Fig. 1).
The low-$T$ statistical averaging over $l$ leads to the characteristic
rotational motion ”freezing” temperature
$T^{\texttt{(r)}}_{X^{{}^{\pm\\!}}}\\!=\\!\hbar^{2}\\!/k_{B}I_{\\!X^{{}^{\pm}}}$
(see, e.g., Ref. Pathria ). By analogy with hydrogen molecular ion (see, e.g,
Ref. AbersQM ) this can be written in terms of the trion binding energy as
$T^{\texttt{(r)}}_{\\!X^{{}^{+\\!}}}\\!\approx\\!\sigma|E_{\\!X^{{}^{+\\!}}}|/k_{B}$
and
$T^{\texttt{(r)}}_{\\!X^{{}^{\\!-\\!}}}\\!\approx\\!|E_{\\!X^{{}^{\\!-\\!}}}|/k_{B}\sigma$
with $\sigma\\!=\\!m_{e}/m_{h}$, whereby the rotational degrees of freedom are
to be frozen out (at least for $\sigma$ close to unity typical of TMDs
TMDmass18 ) as long as the CIEs are stable against thermal fluctuations.

With no rotational term contribution and no external fields, it is
straightforward to obtain a qualitative picture of the like-charge trion
Wigner crystallization by performing an analysis analogous to that reported in
Ref. Platzman74  for the 2D electron gas. In view of Eqs. (2) and (5), the
ratio (1) now takes the form

$\frac{\langle V\rangle}{\langle
K\rangle_{0}}=\frac{2M}{\hbar^{2}}\left(\frac{1}{\sqrt{\pi n}}+d^{2}\sqrt{\pi
n}\right)=\Gamma_{0}>1,$ (6)

whereby it can be seen that the Wigner crystallization process of CIEs can be
controlled by the TMD interlayer distance $d$, an extra parameter absent from
the electron WC problem. The exact expressions for the zero-$T$ critical
density $n_{c}$ and the critical temperature $T^{\texttt{(W)}}_{c}$ of the CIE
WC phase transition come out as CommPhys-bond

$\displaystyle n_{cX^{\pm}}\\!=\\!\frac{2}{\pi a_{B}^{\ast
2}d^{2}}\\!\left(\\!\frac{g_{\pm}\Gamma_{0}}{4d}\\!\right)^{\\!\\!2}\\!\left[1-\frac{1}{2}\\!\left(\frac{4d}{g_{\pm}\Gamma_{0}}\right)^{\\!\\!2}\\!-\sqrt{1\\!-\\!\left(\frac{4d}{g_{\pm}\Gamma_{0}}\right)^{\\!\\!2}}\,\right]\\!,$
$\displaystyle
k_{B}T^{\texttt{(W)}}_{cX^{\pm}}=\frac{4Ry^{\ast}}{g_{\pm}\Gamma_{0}^{2}}\,,\hskip
85.35826pt$ (7) $\displaystyle
g_{\pm}(\sigma)=\left(3+\Big{\\{}\\!\begin{array}[]{c}1\\\\[-7.11317pt]
2\end{array}\\!\Big{\\}}\,\sigma+\Big{\\{}\\!\begin{array}[]{c}2\\\\[-7.11317pt]
1\end{array}\\!\Big{\\}}\,\frac{1}{\sigma}\right)^{\\!\\!-1}\\!\\!.\hskip
45.52458pt$ (12)

Here the 3D ”atomic units” are used LandauQM , with distance and energy
measured in the units of exciton Bohr radius
$a^{\ast}_{B}\\!\\!=\\!0.529\,\mbox{\AA}\,\varepsilon/\mu$ and exciton Rydberg
energy $Ry^{\ast}\\!\\!=\\!\hbar^{2}\\!/(2\mu\,m_{0}a_{B}^{\ast
2})\\!=\\!e^{2}\\!/(2\varepsilon
a_{B}^{\ast})\\!=\\!13.6\,\mbox{eV}\,\mu/\varepsilon^{2}$, respectively,
$\varepsilon$ represents the _effective_ average dielectric constant of the
bilayer heterostructure and $\mu\\!=\\!m_{e}/(\lambda\,m_{0})$ stands for the
exciton reduced effective mass (in the units of free electron mass $m_{0}$)
with $\lambda\\!=\\!1+m_{e}/m_{h}\\!=\\!1+\sigma$. The zero magnetic field
quantities $n_{cX^{\pm}}$ and $T^{\texttt{(W)}}_{cX^{\pm}}$ were analyzed in
great detail in Ref. CommPhys-bond .

###  .2 Non-zero field like-charge trion Wigner crystallization

In the presence of the perpendicular magnetostatic field (see Fig. 1), the
non-interacting-particle energy contribution to the total energy of the system
of identical charged particles with non-zero spin (the many-particle CIE
system in our case) can be written as Mahan

$\frac{\hbar^{2}}{2M}\sum_{j}\Big{[}\Big{(}k_{\\!jx}+\frac{qB}{2\hbar
c}\,y_{j}\Big{)}^{2}+\Big{(}k_{\\!jy}-\frac{qB}{2\hbar
c}\,x_{j}\Big{)}^{2}\Big{]}-\sum_{j}{\bm{\mu}}_{j}\cdot\textbf{B}$ (13)

(summation over particles). Here,
$\textbf{B}=\\!{\bm{\nabla}}\\!\times\textbf{A}$ is the perpendicular
homogeneous magnetic field with $\textbf{A}\\!=\\!(-y,x,0)B/2$ representing
the vector potential in the symmetric gauge. The trion effective mass
$M\,(=\\!M_{X^{\pm}})$ in the first term (Landau energy) is the sum of the
$e$-$h$ band-structure masses $m_{e,h}$. The $j$-th trion magnetic moment
${\bm{\mu}}_{j}$ in the second term (Pauli interaction) is the vector sum of
those of electrons and holes it is composed of, which can be written as
${\bm{\mu}}_{j}\\!=\\!\mbox{g}_{X^{\pm}}\mu_{B\,}\textbf{s}_{j}$ with
$\mbox{g}_{X^{\pm}}$ and $\textbf{s}_{j}$ being the trion g-factor and
spin-1/2 operator, respectively, and $\mu_{B}$ being the Bohr magneton
constant. We focus on the Landau energy term first. The Pauli interaction term
will be added later.

Figure 2: The 2D Fermi surface center shift to the point
$(k_{bx},k_{by})\\!=\\!(-y,x)/(2l^{2})$ for the positive CIE subjected to the
perpendicular homogeneous magnetostatic field B as shown in Fig. 1 (a).

Since there is no interparticle interactions in Eq. (13), the many-particle
Landau energy states can be composed of the single-particle states. The
single-particle Hamiltonian is

$\frac{\hbar^{2}}{2M}\Big{[}\Big{(}k_{x}+\frac{qB}{2\hbar
c}\,y\Big{)}^{2}+\Big{(}k_{y}-\frac{qB}{2\hbar c}\,x\Big{)}^{2}\Big{]}.$ (14)

This can be seen to include the domain of the reciprocal space $(k_{x},k_{y})$
inside of the 2D Fermi surface whose center is shifted out of the point
$k_{x}\\!=\\!k_{y}\\!=\\!0$ due to the perpendicular magnetic field B applied.
For the many-particle Landau state at zero $T$, the average kinetic energy per
particle can then be obtained by summing up over all $k$-states bounded by the
displaced 2D Fermi sphere as shown for the positive CIE in Fig. 2 and given by

$\langle K\rangle=\frac{2}{\langle
N\rangle}\sum_{\mathbf{k}}E_{\mathbf{k}}n_{\mathbf{k}}=\frac{2}{\langle
N\rangle}\frac{S}{(2\pi)^{2}}\,\frac{\hbar^{2}}{2M}\;4\\!\int_{k_{bx}}^{k_{bx}+k_{F}}dk_{x}\int_{k_{by}}^{k_{by}+\sqrt{k_{F}^{2}-(k_{x}-k_{bx})^{2}}}dk_{y}\;\big{(}k_{x}^{2}+k_{y}^{2}\big{)}.$
(15)

Here, $\langle N\rangle$ is still given by Eq. (4), being unaffected by the
Fermi surface shift, $k_{bx}\\!=\\!-\pi By/\phi_{0}$ and $k_{by}\\!=\\!\pi
Bx/\phi_{0}$ are the coordinates of the center of the 2D Fermi surface
shifted, with $\phi_{0}\\!=\\!2\pi\hbar c/|q|$ representing the magnetic flux
quantum. To include both positive and negative CIE case, this can generally be
written as

$\displaystyle k_{bx}=-\mbox{sign}(q)\frac{\pi
By}{\phi_{0}}=-\mbox{sign}(q)\frac{y}{2l^{2}}=-\mbox{sign}(q)\frac{M\omega_{c}y}{2\hbar}\,,$
(16) $\displaystyle k_{by}=\mbox{sign}(q)\frac{\pi
Bx}{\phi_{0}}=\mbox{sign}(q)\frac{x}{2l^{2}}=\mbox{sign}(q)\frac{M\omega_{c}x}{2\hbar}\,,\hskip
21.33955pt$

where

$l^{2}=\frac{\hbar c}{|q|B}=\frac{\hbar}{M\omega_{c}}=\frac{\phi_{0}}{2\pi
B}~{}~{}~{}\mbox{and}~{}~{}~{}\omega_{c}=\frac{|q|B}{Mc}=\frac{\hbar}{Ml^{2}}=\frac{2\pi\hbar
B}{M\phi_{0}}$ (17)

are the magnetic localization length squared and the cyclotron frequency,
respectively.

The integration in Eq. (15) can be done analytically, whereupon using Eq. (4)
to eliminate $k_{F}$ one obtains the exact expression as follows

$\langle K\rangle=\langle K\rangle_{0}f(B,x,y).$ (18)

Here, $\langle K\rangle_{0}\\!=\\!\hbar^{2}\pi
n/2M\\!=\\!\hbar^{2}k_{F}^{2}/4M\\!=\\!E_{F}/2$ with $E_{F}$ being the Fermi
energy. This is the average kinetic energy per particle of Eq. (5) with no
magnetic field, the one that together with Eqs. (1) and (2) yields the set of
Eqs. (7) reported previously CommPhys-bond . With Eq. (17) this can now also
be written as $\langle K\rangle_{0}\\!=\\!\hbar\omega_{c\,}l^{2}\pi n/2$. The
dimensionless function

$f(B,x,y)=1-\frac{64}{9\pi^{2}}+\frac{\hbar\omega_{c}}{2E_{F}}\Big{[}\Big{(}\frac{x-x_{0}}{\sqrt{2}l}\Big{)}^{2}+\Big{(}\frac{y-y_{0}}{\sqrt{2}l}\Big{)}^{2}\Big{]}$
(19)

with

$x_{0}(n,l)=-y_{0}(n,l)=-\mbox{sign}(q)\frac{8}{3}\sqrt{\frac{2n}{\pi}}\,l^{2}$
(20)

is responsible for the magnetic field effect of an arbitrary strength. From
Eq. (19) the motion of particles in our system can be seen to be divided in
translational and rotational parts. The latter is given by the last term
describing the in-plane rotation at frequency $\omega_{c}$ about the point
$(x_{0},y_{0})$ in circles with the localization radius $\sqrt{2}\,l\,$. The
magnitude of this term is controlled by its prefactor, the rotational-to-
translational energy ratio. If this ratio is much less than unity, then the
translational motion dominates, otherwise the rotational motion does, which
corresponds to weak and strong magnetic fields, respectively. More insight
into this can be gained by plugging Eq. (20) in Eq. (19) and regrouping the
terms therein to obtain an equivalent form as follows

$f(B,x,y)=1+\frac{\hbar\omega_{c}}{2E_{F}}\Big{[}\frac{\langle
L_{z}\rangle}{\hbar}+\Big{(}\\!\frac{x}{\sqrt{2}l}\Big{)}^{\\!2}+\Big{(}\\!\frac{y}{\sqrt{2}l}\Big{)}^{\\!2}\Big{]},$
(21)

whereby the rotational motion in our system can be seen more explicitly. Here,

$\langle
L_{z}\rangle=\mbox{sign}(q)\frac{16}{3\pi}\big{(}x\bar{p}_{y}-y\bar{p}_{x}\big{)}$
(22)

is the averaged out-of-plain angular momentum with
$\bar{p}_{x,y}\\!=\\!\hbar\bar{k}_{x,y}\\!=\\!\sqrt{\langle
p_{x,y}^{2}\rangle}\\!=\\!\hbar\sqrt{\pi n/2}$ being the in-plain
translational momentum components statistically averaged. For weak magnetic
fields $\hbar\omega_{c}/2E_{F}\\!\ll\\!1$ with $l\\!\rightarrow\\!+\infty$,
which gives $f(B,x,y)\\!=\\!1$ in Eq. (21) regardless of what $x$ and $y$ are.
In Eq. (19), additionally, the absolute values of $x_{0}$ and $y_{0}$ tend to
infinity per Eq. (20), yielding $(x-x_{0})^{2}\\!\sim\\!x_{0}^{2}$ and
$(y-y_{0})^{2}\\!\sim\\!y_{0}^{2}$ with the same end result. For very strong
magnetic fields, on the other hand, one has
$x-x_{0}\\!=\\!l\sqrt{2}\cos\omega_{c}t$ and
$y-y_{0}\\!=\\!l\sqrt{2}\sin\omega_{c}t$ with $l\\!\rightarrow\\!+0$, whereby
both equations take the universal coordinate independent form

$f(B)=1-\frac{64}{9\pi^{2}}+\frac{\hbar\omega_{c}}{2E_{F}}\approx\frac{\hbar\omega_{c}}{2E_{F}}=\frac{M\omega_{c}^{2}l^{2}}{2E_{F}}=\frac{1}{2l^{2}\pi
n}=\frac{1}{\nu(n,B)}\gg 1\,.$ (23)

Here, $\nu$ is the filling factor function that determines the number of
occupied states within the circular area of radius $\sqrt{2}l$ by charged
particles of density $n$ bound inside Mahan ,

$\nu=\pi\big{(}\sqrt{2}l\big{)}^{2}n=\frac{\phi_{0}}{B}\,n=\\!\sum_{N,m_{s}}\\!n_{F}(E_{N,m_{s}}),$
(24)

where $n_{F}$ is the standard Fermi-Dirac distribution of Eq. (3) with

$E_{N,m_{s}}=\hbar\omega_{c}\Big{(}N+\frac{1}{2}\Big{)}-\mbox{g}_{X^{\pm}}m_{s}\mu_{B}B$
(25)

representing the quantized energies (Landau levels) that are characterized by
nonnegative integers $N\,(=\\!0,1,2,...)$ corresponding to the classical
rotation in quantized circular orbits. The Pauli interaction term, where
$m_{s}\,(=\\!\pm 1/2)$ is the spin-1/2 projection quantum number, is hereby
added for completeness. The lowest Landau level has $N\\!=\\!0$. All Landau
states are highly degenerate as the centers of orbits can be located at any
point inside the finite-size area $S\\!=\\!L_{x}L_{y}$ occupied by the sample.
The last equality in Eq. (24) can be obtained by writing the particle density
$n$ in terms of the density of states $\rho(E)$ per unit area, whereby

$n=\\!\int\\!\\!dE\,n_{F}(E)\rho(E)=\\!\int\\!\\!dE\,n_{F}(E)\frac{1}{S}\sum_{\bar{k}_{x}=0}^{L_{y}/l^{2}}\sum_{N,m_{s}}\\!\delta(E\\!-\\!E_{N,m_{s}})=\frac{1}{2\pi
l^{2}}\sum_{N,m_{s}}\\!n_{F}(E_{N,m_{s}}).$

This takes into account that in a finite-size sample, for the quantized
$\bar{k}_{x}\\!=\\!2\pi\\{0,1,2,\dots\\}/L_{x}$ the inequality
$0\\!\leq\\!\bar{k}_{x}l^{2}\\!=\\!3\pi|y_{0}|/16\\!\leq\\!L_{y}$ is to be
fulfilled as per Eq. (20) (assuming $l\\!\ll\\!L_{y}$). Plugging this in the
first part of Eq. (24) yields the last part.

From Eqs. (21), (23) and (25) it follows that the lowest Landau level
$\hbar\omega_{c}/2$ is the only one to contribute to the function $f(B)$ in
the strong field regime. In this regime, a charge carrier is localized within
the area $\pi(\sqrt{2}l)^{2}$ due to its circular motion with out-of-plane
angular momentum of absolute value $M(l\sqrt{2})^{2}\omega_{c}/2\\!=\\!\hbar$,
as per Eq. (17). If the magnetic localization length $l$ is much less than the
interparticle distance $a$, then the repulsive interparticle interaction can
promote the long-range (Wigner) crystal phase development. For the triangular
lattice typical of highly correlated many-particle systems Kleiner64 , one has
$n\\!=\\!1/S_{0}$ with $S_{0}\\!=\\!a^{2}\sqrt{3}/2$ being the area per
particle, whereby the condition $2\sqrt{2}\,l\\!<\\!a$ leads to an approximate
upper bound $\pi\sqrt{3}/6\\!\approx\\!0.907$ for the filling factor $\nu$ not
to exceed in order for the system to stay in the crystal phase. Rigorous
microscopic quantum theories of melting yield the 2D magnetic-field-induced
crystallization requirement in the form LozFarAbd85

$\nu=\frac{\phi_{0}}{B}\,n=\frac{n}{n_{0}}=\frac{B_{0}}{B}\leq\nu_{c}\approx
0.13,~{}~{}~{}0\leq\nu\leq 1,$ (26)

mostly because of the lattice anharmonicity, for both positive and negative
charge carriers regardless of their repulsive interparticle interaction
strength or the crystal lattice type endnote1 . Here,
$n_{0}(B)\\!=\\!B/\phi_{0}\\!=\\!1/\pi(\sqrt{2}l)^{2}$ is the number of
magnetic flux quanta per unit area (or the lowest Landau level degeneracy per
unit area to set up its maximum occupancy) and $B_{0}(n)\\!=\\!\phi_{0}n$ is
the magnetic field lower bound to be exceeded for the system to enter the
strong field regime. Thus, no large density is required, the strong field is
instead, for the magnetic-field-induced WC phase to occur. This conclusion can
also be reached from the zero-field expressions in Eq. (7). Plugging Eq. (18)
in Eq. (1) redefines the energy ratio of interest to the form
$\Gamma_{0}/f(B)\\!=\\!\nu\Gamma_{0}\\!=\\!\Gamma$, as per Eq. (23). Solving
this for $\Gamma_{0}$ and plugging the result in the critical density
expression of Eq. (7) yields

$n_{cX^{\pm}}(\nu)\approx\frac{\nu^{2}}{4\pi(g_{\pm}\Gamma)^{2}a_{B}^{\ast
2}}=\frac{1}{4\pi(g_{\pm}\Gamma_{0})^{2}a_{B}^{\ast
2}}\Big{(}\frac{\nu}{\nu_{c}}\Big{)}^{2}$

to the lowest nonvanishing order in $\nu d/\Gamma\,(=\\!\nu
d/\nu_{c}\Gamma_{0}\\!\ll\\!1$; the $d$ dependence comes out of the next order
term), whereby per Eq. (26) no nonzero lower limit exists for $n_{cX^{\pm}}$
in the strong field regime. Similarly, for the critical temperature of Eq. (7)
one obtains

$k_{B}T^{\texttt{(W)}}_{cX^{\pm}}(\nu)=\frac{\nu_{c}^{2}Ry^{\ast\\!}}{4g_{\pm}\Gamma^{2}}=\frac{Ry^{\ast}}{4g_{\pm}\Gamma_{0}^{2}}\Big{(}\frac{\nu_{c}}{\nu}\Big{)}^{2}=\frac{0.13^{2}Ry^{\ast}}{4g_{\pm}\Gamma_{0}^{2}}\left(\frac{B}{B_{0}}\right)^{2},$
(27)

presenting $T^{\texttt{(W)}}_{cX^{\pm}}$ as an inverse quadratic function of
$\nu$ and as a quadratic function of $B$ in the strong field regime.

###  .3 Optical detection of the magnetic-field-induced trion WC phase

Recently, a significant enhancement of the electronic transition coupling
strength to an external magnetostatic field was demonstrated experimentally
for TMD heterobilayers Nagler17 . Here, the type II (staggered gap) band
alignment at stacking angles $0^{\circ}$ (AA stacking) or $60^{\circ}$ (AB
stacking) enables the formation of IEs with electrons and holes in different
layers. As opposed to monolayers, in TMD bilayers the valley-conserving
requirement is lifted for optical $e$-$h$ transitions, allowing the IE-forming
charge carriers at the $K$-points of valleys with different index (of
different layers) to recombine radiatively. In an external magnetic field,
their recombination yields two-component circularly polarized PL emission,
which in the case of AB-stacked heterobilayers comes from the valley splitting
significantly enhanced as compared to monolayers, due to the addition of
conduction-band and valence-band valley magnetic moments. The effect allows
for efficient external manipulation by the spin-valley degrees of freedom in
the TMD heterobilayers, opening the door to new fundamental physics only
achievable with these transdimensional quantum materials.

In magneto-PL experiments Nagler17 , the TMD heterobilayer sample is first
excited by linearly polarized light to equally populate both valleys of
individual monolayers. Then, the optical detection and analysis follow of the
two circularly polarized PL emission components ($\sigma^{+}$ and
$\sigma^{-}$) resulted from the valley-selective splitting of the IE
recombination transitions in the magnetic field. The total transition energy
splitting is

$\Delta
E=E^{\sigma+}\\!-E^{\sigma-}\\!=\mbox{g}_{\texttt{eff}}\mu_{B}B,~{}~{}~{}\mbox{g}_{\texttt{eff}}=\pm\Big{(}4+\frac{2}{m_{h}}\mp\frac{2}{m_{e}}\Big{)},$
(28)

where masses are in units of the free electron mass, and the upper (lower)
sign in the effective g-factor is for the AA-stacked (AB-stacked)
heterobilayer Nagler17 .

Figure 3: Sketch (top view) of the crystal phase (a) and the liquid phase (b)
for the negative CIEs in the strong perpendicular magnetic field (directed
upward). In (a), the red periodically aligned circles with arrows show the
trajectories for the negative CIE circular motion on the magnetically-induced
WC lattice. In (b), the circular trajectories are increased in size to
overlap, which breaks down Eq. (26) either due to the field strength $B$
decreased or due to the doping level $n$ increased. The overlap areas are
shown in white to indicate that the CIEs can now hop from circle to circle
randomly, which is the liquid phase with overall rotational quasi-free CIE
motion.

The TMD bilayer heterostructures can also be doped electrostatically by using
two optically transparent electrical gates at the top and bottom monolayers
Lius-PKim19 . In this case, positive or negative CIEs form in the sample as
per the sign of the bias voltage applied. In the sufficiently strong
perpendicular magnetic field, they will be trapped by the magnetic-field-
induced parabolic potential and (Wigner) crystallized due to their Coulomb
repulsion. That is, the CIEs formed will be constrained to move on
periodically ordered circular orbits of area $\pi(\sqrt{2}l)^{2}$ with angular
momentum $M(l\sqrt{2})^{2}\omega_{c}/2\\!=\\!\hbar$, as shown in Fig. 3 (a),
to contribute the respective magnetic moment to $\mbox{g}_{\texttt{eff}}$ in
Eq. (28). Increasing the doping level with magnetic field strength unchanged
will initiate the melting of this WC of CIEs as prescribed by Eq. (26). Some
of the CIEs will then be forced to leave their respective Wigner lattice sites
for a relatively free (quasi-free) random motion on overall circular orbits of
much greater radii as shown in Fig. 3 (b), being still driven by the magnetic
field. Clearly, such a quasi-free motion will reduce (but will not eliminate)
the $\mbox{g}_{\texttt{eff}}$ angular momentum contribution that comes from
the CIE crystallization.

In twisted TMD bilayers, a combination of structural rippling and electronic
coupling gives rise to periodic moiré potentials that can confine neutral and
charged excitations and thus impede their motion MoireIXdiff20 . At the same
time, experiments have revealed that the largest moiré potentials are realized
for quite a narrow range of angles that need to be controlled precisely, so
that engineering of high-quality moiré lattice structures is a challenge
MoireColumbia21 . Even bilayer samples grown by chemical vapor deposition
(CVD) or created by the mechanical exfoliation and transfer (MET) with a well-
controlled twist angle present different results MoireIXdiff20 . The IE
diffusion length in the rotationally aligned CVD-grown heterobilayers exceeds
the size of the sample (a few microns), whereas no IE diffusion is observed in
the MET stacked heterobilayers with $\sim\\!1^{\circ}$ twist angle, suggesting
the IE localization by the moiré potential. Mobile excitons and trions are
preferred in our case. They are more likely found in either CVD-grown samples
or in MET heterobilayer samples with small moiré periods where they can tunnel
between the moiré potential minima, or in commensurate homo- and heterobilayer
structures with zero and near-zero moiré periods.

The role of a periodic moiré potential can be qualitatively understood in
terms of the effective mass model. A moiré potential with period $a_{m}$
greater than the intrinsic period $a$ of individual monolayers (always the
case), adds to the CIE effective mass to shrink the CIE band and also
contracts the first Brillouin zone of the reciprocal space as sketched in Fig.
4. Being independent of the effective mass, the ratio $\hbar\omega_{c}/2E_{F}$
remains unchanged, so that
$\hbar\omega_{c}/2E_{F}\\!=\\!\hbar\omega_{cm}/2E_{Fm}$, see Fig. 4, but the
respective magnetic length $l_{m}$ increases and the critical field
$B_{0m}(n_{m})\\!=\\!\phi_{0}n_{m}$ drops to keep Eq. (23) constant. This
leads to the $T_{c}$ increase as per Eq. (27). If $B\\!<\\!B_{0m}$ (weak field
regime), then $\hbar\omega_{cm}/2E_{Fm}\\!<\\!1$ and the narrow-band
translational motion is dominant for the CIEs trapped by the periodic moiré
potential to form a ’generalized’ Wigner crystal similar to that observed for
electrons in twisted TMD bilayers recently FWang20 ; FWang21 . If
$B\\!>\\!B_{0m}$ (strong field regime), then
$\hbar\omega_{cm}/2E_{Fm}\\!>\\!1$ and the rotational motion dominates for the
CIEs about the centers of their respective moiré spots, provided though that
the number of magnetic flux quanta
$\alpha\\!=\\!S_{0m}B/\phi_{0}\\!=\\!S_{0m}/\pi(\sqrt{2}l_{m})^{2}$ traversing
the unit cell of the moiré lattice is not integer (the Hofstadter parameter
Hofstadter ), where $S_{0m}\\!=\\!a_{m}^{2}$ and $a_{m}^{2}\sqrt{3}/2$ for
square and triangular moiré lattices, respectively (see also Refs. Azbel64 ;
ClaroWannier79 ; ChangNiu96  for advanced theory and Ref. FrancoNori16  for
review of recent experimental work). The magnetic-field-induced CIE Wigner
crystallization is to proceed as prescribed by Eq. (26) with $B_{0}$ replaced
by $B_{0m}$ in this case. If $\alpha$ is integer, then the magnetic field
effect vanishes identically Hofstadter , whereby the CIEs form a narrow-band
’generalized’ Wigner crystal on the moiré lattice, with no rotation present in
the system and so no $\mbox{g}_{\texttt{eff}}$ change in Eq. (28). An example
we consider below is the commensurate AA-stacked MoSe2/WSe2 heterobilayer
system studied in Ref. Lius-PKim19 .

Figure 4: Effective mass model sketch (not to scale) for the first Brillouin
zones of the CIEs without and with periodic moiré potential present. The
parabolic bands (black and blue lines, respectively) are chosen to be half-
filled (shaded in grey and light blue) with $E_{F}$ and $E_{Fm}$ indicating
their Fermi levels. The moiré potential adds to the CIE effective mass to
narrow the band. The moiré pattern period $a_{m}$ is greater than the
monolayer period $a$, to contract the zone as compared to no moiré case. The
lowest Landau levels (red lines) are such that
$\hbar\omega_{c}/2E_{F}\\!=\hbar\omega_{cm}/2E_{Fm}$ and the magnetic length
$l_{m}$ in the presence of moiré is greater than $l$ of no moiré case.

It has been recently shown experimentally Xu21  that the g-factors obtained
using Eq. (28) are indistinguishable for IEs and CIEs. Indeed, the Zeeman
effect of magnetic field on the trion energy is given by the Pauli interaction
term in Eq. (25), and can be viewed as the sum of the IE Zeeman shift and that
of the extra carrier. Clearly, the extra electron (hole) in a negative
(positive) trion only shifts the transition energy Xu21 , but does not
contribute to the spectroscopic Zeeman splitting given by Eq. (28). This does
not include the strong field regime, however, where in addition to the Pauli
interaction term the Landau energy term of Eq. (25) contributes the angular
momentum of Eq. (22) to the effective g-factor of the entire charge system due
to the circular motion of the CIE as a compound charged particle, which leads
to the WC formation of CIEs provided that Eq. (26) is fulfilled.

A simple phenomenological model allows one to understand the generic behavior
of the effective g-factor in Eq. (28) in the case of the Wigner solid-liquid
phase transition of the CIEs in doped TMD heterobilayers in the strong
magnetic field regime. Let the CIE trapping rate $\gamma_{t}$ to and
detrapping rate $\gamma_{f}$ from the magnetic-field-induced parabolic
localization potential be of the standard exponential form as follows

$\gamma_{t}=A_{t}e^{a_{t}U},~{}~{}~{}\gamma_{f}=A_{f}e^{a_{f}(U-{\cal{E}})}.$
(29)

Here, ${\cal{E}}$ stands for the CIE crystal-liquid phase internal energy
difference, $U$ represents the charge carrier potential energy change with
doping (proportional to the gate voltage that controls the dopant
concentration), both in units of $k_{B}T$, and $A_{t,f}$ are the CIE trapping-
detrapping amplitudes. The constants $a_{t,f}$ in the exponents are necessary
to convert the charge carrier energy added into the sample with doping into
that of the CIEs formed due to doping since not all of the doped charge
carriers participate in the CIE formation. The statistical detailed balance
further requires the condition

$\gamma_{t}B_{f}=\gamma_{f}B_{t}e^{a_{f}{\cal{E}}+(a_{t}-a_{f})U}$

to be fulfilled, where $B_{t}$ and $B_{f}$ are the effective numbers of
trapped and quasi-free states available for quasi-free and trapped CIEs to get
trapped and get quasi-free, respectively. Plugging Eq. (29) in here gives
$A_{f}\\!=\\!B_{f}A_{t}/B_{t}=CA_{t}$, whereby the detrapping amplitude comes
out being proportional to the trapping amplitude. Next, we introduce the
relative weights of the $e$-$h$ recombination components to take part in the
magneto-PL emission under doping. They are, respectively, $w_{t}^{+}$,
$w_{f}^{+}$, $w_{h}$ for trapped $X^{+}$, quasi-free $X^{+}$ and holes when
doped by holes, and $w_{t}^{-}$, $w_{f}^{-}$, $w_{e}$ for trapped $X^{-}$,
quasi-free $X^{-}$ and electrons when doped by electrons, all defined
explicitly as follows

$w_{t,f}^{\pm}=\frac{\gamma_{t,f}(X^{\pm})}{1+\gamma_{t}(X^{\pm})+\gamma_{f}(X^{\pm})},~{}~{}w_{h,e}=\frac{1}{1+\gamma_{t}(X^{\pm})+\gamma_{f}(X^{\pm})},~{}~{}w_{h,e}+w_{t}^{\pm}+w_{f}^{\pm}=1.$
(30)

Here, $\gamma_{t,f}(X^{+})$ and $\gamma_{t,f}(X^{-})$ are given by Eq. (29)
with $A_{t,f}$ and $a_{t,f}$ specified for the positive and negative CIEs
individually. These weight functions of $\cal{E}$ and $U$ now require that
$\mbox{g}_{\texttt{eff}}$ of Eq. (28) be redefined consistently to take into
account the positive and negative CIEs formed under hole and electron doping,
respectively. We also note that in the vicinity of the CIE crystal-liquid
phase transition, taking into account Eq. (27), the following is true:

${\cal{E}}=\frac{\delta E}{k_{B}T}\approx\frac{\delta
E}{k_{B}T^{\texttt{(W)}}_{cX^{\pm}}}=\kappa\Big{(}\frac{\nu}{\nu_{c}}\Big{)}^{2},~{}~{}~{}\kappa=\delta
E\,\frac{4g_{\pm}\Gamma_{0}^{2}}{Ry^{\ast}}\,.$ (31)

Here, $\delta E$ is the crystal-liquid phase internal energy difference in the
physical units and $\kappa$ is the dimensionless constant. With this in view,
one finally obtains the effective g-factor as a function of the two
independent variables $U$ and $\nu/\nu_{c}$. The former is controlled by the
gate voltage and varies the doped charge carrier concentration to result in
the trion formation. The latter, treated as a function of $B$ only as per Eqs.
(26) and (27), can be viewed as that setting up the _initial_ pre-doped value
of the filling factor $\nu=\bar{\nu}=\\!B_{0}/B$ for $e$-$h$ pairs created by
absorption of light (and controlled by light intensity) in an _undoped_ TMD
bilayer (with no trions formed) subjected to the strong external magnetic
field.

The effective g-factor of Eq. (28) is now to be redefined for the hole
($\nu_{h}$) and electron ($\nu_{e}$) doping individually as follows:

$\mbox{g}_{\texttt{eff}}^{+}=\pm\Big{[}4+2\Big{(}\frac{w_{h}}{m_{h}}+\frac{w_{t}^{+}+w_{f}^{+}}{M_{X^{+}}}\Big{)}\mp\frac{2}{m_{e}}\Big{]}$
(32)

with $M_{X^{+}}\\!=\\!2m_{h}+m_{e}$ for doping with holes and

$\mbox{g}_{\texttt{eff}}^{-}=\pm\Big{[}4+\frac{2}{m_{h}}\mp
2\Big{(}\frac{w_{e}}{m_{e}}+\frac{w_{t}^{-}+w_{f}^{-}}{M_{X^{-}}}\Big{)}\Big{]}$
(33)

with $M_{X^{-}}\\!=\\!2m_{e}+m_{h}$ for doping with electrons.

Figure 5: The effective g-factor as given by Eqs. (29)–(33) of the CIE Wigner
crystallization model for the positive (top left) and negative (top right)
CIEs in a generic AA-stacked TMD heterobilayer with doped carrier energy $U$
and initial filling factors $\bar{\nu}_{e,h}/\nu_{c}$ varied in the strong
perpendicular magnetic field regime. Bottom panel shows the crosscuts of the
top graphs by the vertical planes $\bar{\nu}_{e,h}/\nu_{c}=0.9$, $1.0$, $1.1$
(line thickness increase indicates stronger field applied). See text for
details.

Figure 5 shows our calculations of the effective g-factor as given by the CIE
Wigner solid-liquid model in Eqs. (29)–(33) for the positive (top left panel)
and negative (top right panel) CIEs in a generic AA-stacked TMD heterobilayer
taken as an example. For definiteness, we use the $e$-$h$ effective masses of
the MoSe2/WSe2 heterostructure studied experimentally in Ref. Lius-PKim19 .
They are $m_{h}\\!=\\!0.36$ for the holes of the WSe2 monolayer Kormanyos  and
$m_{e}\\!=\\!0.8$ for the electrons of the MoSe2 monolayer Larentis18 . This
gives $\mbox{g}_{\texttt{eff}}\\!\approx\\!7$ for the undoped heterobilayer
described by Eq. (28), in which case no CIEs are present in the system. Other
dimensionless model parameters we use are chosen to be consistent with the
detailed balance condition above (to be obtained from fitting our model to the
experimental data when available). They are as follows: $A_{t}\\!=\\!4\times
10^{-5}$, $a_{t}\\!=\\!1.7$, $A_{f}\\!=\\!10^{3}A_{t}$, $a_{f}\\!=\\!3a_{t}$,
$\kappa\\!=\\!5.5$ for the positive CIE and $A_{t}\\!=\\!10^{-6}$,
$a_{t}\\!=\\!1.7$, $A_{f}\\!=\\!10^{3}A_{t}$, $a_{f}\\!=\\!3a_{t}$,
$\kappa\\!=\\!6.5$ for the negative CIE. The graphs in both panels can be seen
to develop the valleys as $U$ departs from zero with doping under increasing
gate voltage, whereby the CIEs start forming at small concentrations and
getting crystallized as long as the condition $\nu_{h,e}\\!<\\!\nu_{c}$ (or
$\nu_{h,e}/\nu_{c}\\!<\\!1$) is fulfilled for increasing $\nu_{h,e}$ as they
depart from their initial values $\bar{\nu}_{e,h}$ per Eq. (26). The above-
mentioned $\kappa$ parameters have been adjusted to position the planes
$\bar{\nu}_{e,h}/\nu_{c}\\!=\\!1$ (light grey color in the figure) right where
these ’crystallization’ valleys start developing with $\bar{\nu}_{h,e}$
decreasing, which allows one to evaluate $\delta E$ of Eq. (31). With $e$-$h$
effective mass ratio
$\sigma\\!=\\!m_{e}/m_{h}\\!=\\!0.8/0.36\\!\approx\\!2.2$, Eq. (7) gives
$g^{+}\\!\approx\\!1/6.3$ and $g^{-}\\!\approx\\!1/7.9$, whereby solving Eq.
(31) for $\delta E$ one obtains $5.5\times 6.3\\!\approx\\!35$ and
$6.5\\!\times 7.9\\!\approx\\!51$ in units of $Ry^{\ast}/4\Gamma_{0}^{2}$ for
$X^{+}$ and $X^{-}$, respectively. Next, using $Ry^{\ast}\\!\approx\\!0.1$ eV
recently reported for TMDs CommPhys-bond  and $\Gamma_{0}\\!=\\!55$ reported
previously for 2D parallel dipole systems from molecular dynamics simulations
BedGaLoz85jetp , we arrive at $\delta E\\!\approx\\!0.3$ meV for $X^{+}$ and
$\delta E\\!\approx\\!0.4$ meV for $X^{-}$, which is consistent with earlier
numerical studies of the Wigner crystallization phenomena in 2D $e$-$h$
systems BedGaLoz85jetp . To give an idea of temperatures, densities and fields
involved in our case, we start with $n\\!\approx\\!5\\!\times\\!10^{11}$ cm-2
reported for IEs in the experiments with MoSe2/WSe2 heterostructures Lius-
PKim19 . The CIE density is expected to be significantly lower than that since
only a fraction of doped charge carriers couples to IEs to form CIEs. Taking
$n\\!\sim\\!10^{10}$ cm-2 for CIEs, Eq. (26) yields $B_{0}\\!\approx\\!0.4$ T
requiring that $B\\!>\\!B_{0}/\nu_{c}\\!\approx\\!0.4/0.13\\!\approx\\!3$ T
for the WC phase of CIEs to occur. If $B\\!=\\!45$ T, for example, then
$B/B_{0}\\!=\\!\nu_{c}/\nu\\!=\\!15$ to give
$T^{\texttt{(W)}}_{cX^{+}}\\!\approx\\!2$ K and
$T^{\texttt{(W)}}_{cX^{-}}\\!\approx\\!3$ K as per Eq. (27) with other
parameters as mentioned above. These values are all practically achievable in
real experiments.

If the initial values $\bar{\nu}_{h,e}$ are less than $\nu_{c\,}$, then with
$U$ further departing from zero the CIE densities increase together with their
respective doped charge carrier densities, pushing the _current_ $\nu_{h,e}$
values closer and closer to the $\nu_{c}$ value to finally exceed it, which
breaks the WC phase condition of Eq. (26). This is why in the top left and top
right panels of Fig. 5 after passing the crystallization drop-off the
effective g-factors raise up and plateau at values quite different from
$\mbox{g}_{\texttt{eff}}\\!\approx\\!7$, the no-doping value. The CIE WC melts
here as the crystallization condition (26) is no longer met when the _then-
current_ $\nu_{h,e}$ values start exceeding the $\nu_{c}$ value. To give a
better view of how this happens, at the bottom of Fig. 5 we show the crosscuts
of the top graphs by the vertical planes $\bar{\nu}_{e,h}/\nu_{c}\\!=\\!0.9$,
$1.0$, and $1.1$, where thicker lines indicate smaller
$\bar{\nu}_{e,h}/\nu_{c}$ (or stronger magnetic fields applied). For
$\bar{\nu}_{e,h}/\nu_{c}\\!=\\!0.9$ the trion liquid phase follows the WC
phase with $U$ increasing under doping, whereas for
$\bar{\nu}_{e,h}/\nu_{c}\\!=\\!1.0$ and $1.1$ the WC phase does not occur and
trions formed under charge doping enter the liquid phase directly. Either of
these scenarios will always be the case for the initially undoped bilayer
system in the strong field regime with charge carrier doping level smoothly
increasing from zero. They should in principle be distinguishable in the CIE
transport measurements similar to those reported in Ref. Lius-PKim19 , not
only in the magneto-PL measurements we are focusing on here. As opposed to the
CIE liquid phase, the WC phase is not expected to show in-plane current due to
being pinned on circular orbits inside the sample [cf. Fig. 3 (a) and Fig. 3
(b)]. If, on the other hand, the bilayer system is first sufficiently doped to
form the CIEs and then is subjected to a smoothly increasing perpendicular
magnetic field, neither of the two CIE phases will occur until the field gets
strong enough to approach the Wigner crystallization condition. Thus, the
strong-field magneto-PL experiments with TMD heterobilayes of systematically
varied $e$-$h$ doping concentrations make it possible to study this unique
strongly-correlated many-particle phenomenon of CIE Wigner crystallization.

###  .4 Conclusion and final remarks

We develop the theory for the magnetic-field-induced Wigner crystallization of
charged interlayer excitons that were discovered recently in TMD
heterobilayers Lius-PKim19 . We derive the ratio of the average potential
interaction energy to the average kinetic energy for the many-particle CIE
system subjected to the perpendicular magnetic field of an arbitrary strength
at zero temperature. We analyze the weak and strong field regimes and obtain
the critical temperature for the ’cold’ crystallization phase transition in
the strong field regime where the lowest Landau level alone is occupied. We
also generalize the effective g-factor concept previously formulated for
interlayer excitons Nagler17 , to include the formation of CIEs under
electrostatic doping of heterobilayers in the strong perpendicular magnetic
field. We show that the magnetic-field-induced Wigner crystallization and
melting of the many-particle system of CIEs can be detected and monitored in
magneto-PL experiments that measure the effective g-factor of bilayer TMD
heterostructures under controlled electrostatic doping. The CIE crystal-liquid
internal energy difference is estimated to be $\sim\\!0.3$ meV for a typical
case of the AA-stacked MoSe2/WSe2 heterobilayer.

On a final note, we wish to emphasize that the magnetic-field-induced 2D
Wigner crystallization effect in the strong perpendicular field is a universal
phenomenon. In this regime, the only internal parameter that remains there to
control the properties of the system is the cyclotron frequency, which exceeds
all other internal parameters including crystal lattice vibrational
frequencies and interparticle coupling strength. Hence, any many-particle (or
quasiparticle) system with repulsive interparticle interactions, be it
fermionic or bosonic, be it charged such as electrons, holes, CIEs,
quaternions, or even neutral such as IEs, will experience the Wigner
crystallization phase transition (which we will give the full account of
elsewhere) when the external magnetostatic field is strong enough to isolate
the lowest Landau level and make the cyclotron frequency the leading parameter
of the system. That is why the Wigner crystallization condition (26) we use
here is universal, being independent of the crystal lattice type and the sort
of (quasi-)particles crystallized. Despite this universality, however, only
for CIEs the Wigner crystallization effect can be detected and monitored in
the effective g-factor measurement experiments we discuss here, as the
features we predict for the g-factor are all coming from the additional
magnetic moment of the _charged_ IEs that is associated with their rotational
angular momentum in the crystallized phase.

##  I Data Availability

The data that supports the findings of this study are available within the
article.

##  II Code Availability

Wolfram Mathematica code used to generate the graphs in Fig. 5 is available
from the corresponding author upon email request.

##  III Acknowledgments

I.V.B. acknowledges discussions with Luis Jauregui (UC Irwine), Andrew Joe
(Harvard), and Philip Kim (Harvard) at the beginning of this work. This
research is supported by the U.S. Department of Energy, Office of Science,
Office of BES under award No. DE-SC0007117 (I.V.B.) and by the RFBR grants No.
20-02-00410 and No. 20-52-00035 (Yu.E.L.). Yu.E.L. is supported by the Basic
Research Program at the National Research University HSE. I.V.B. acknowledges
hospitality of the Kavli Institute for Theoretical Physics (KITP), UC Santa
Barbara, where this work was completed during his invited visit as a KITP
Fellow 2022–23 (the research program supported by the U.S. National Science
Foundation under Grant No. PHY-1748958).

##  IV Author contributions

Both authors conceived the project. I.V.B. developed the theory, carried out
numerical calculations, and wrote the final version of the manuscript. Yu.E.L.
provided expertise in many-particle electron-hole correlations and
crystallization phenomena. The authors discussed the results, their
significance, and the ways to best represent them in the manuscript.

##  V Competing interests

The authors declare no competing financial interests.

## References

  * (1) L.V.Keldysh and A.N.Kozlov, Collective properties of excitons in semiconductors, Phys. JETP 27, 521 (1968).
  * (2) Yu.E.Lozovik and V.I.Yudson, A new mechanism for superconductivity: pairing between spatially separated electrons and holes, Sov. Phys. JETP 44, 389 (1976).
  * (3) T.Fukuzawa, S.Kano, T.Gustafson, and T.Ogawa, Possibility of coherent light emission from Bose condensed states of SEHPs, Surf. Sci. 228, 482 (1990).
  * (4) G.E.Astrakharchik, J.Boronat, I.L.Kurbakov, and Yu.E. Lozovik, Quantum phase transition in a two-dimensional system of dipoles, Phys. Rev. Lett. 98, 060405 (2007).
  * (5) G.J.Schinner, J.Repp, E.Schubert, A.K.Rai, D.Reuter, A.D.Wieck, A.O.Govorov, A.W.Holleitner, and J.P. Kotthaus, Confinement and interaction of single indirect excitons in a voltage-controlled trap formed inside double InGaAs quantum wells, Phys. Rev. Lett. 110, 127403 (2013).
  * (6) O.L.Berman, R.Ya.Kezerashvili, and S.M.Tsiklauri, Trions in coupled quantum wells and Wigner crystallization, J. Mod. Phys. B 28, 1450064 (2014).
  * (7) M.M.Fogler, L.V.Butov, and K.S.Novoselov, High-$T$ superfluidity with indirect excitons in van der Waals heterostructures, Nature Commun. 5, 4555 (2014).
  * (8) R.A.Suris, Gas-crystal phase transition in a 2D dipolar exciton system, JETP 122, 602 (2016).
  * (9) Y.E.Lozovik, New effects in and the control of the exciton system in quasi-two-dimensional structures, Physics–Uspekhi 61, 1094 (2018).
  * (10) N.A.Asriyan, I.L.Kurbakov, A.K.Fedorov, and Yu.E.Lozovik, Optical probing in a bilayer dark-bright condensate system, Phys. Rev. B 99, 085108 (2019).
  * (11) M.Kremser, M.Brotons-Gisbert, J.Knörzer, J.Gückelhorn, M.Meyer, M.Barbone, A.V.Stier, B.D.Gerardot, K.Müller, and J.J.Finley, Discrete interactions between a few interlayer excitons trapped at a MoSe2-WSe2 heterointerface, npj 2D Mater Appl 4, 8 (2020).
  * (12) E.C.Regan, D.Wang, C.Jin, M.I.Bakti Utama, B.Gao, X.Wei, S.Zhao, W.Zhao, Z.Zhang, K.Yumigeta, M.Blei, J.D.Carlström, K.Watanabe, T.Taniguchi, S.Tongay, M.Crommie, A.Zettl, and F.Wang, Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices, Nature 579, 359 (2020).
  * (13) H.Li, S.Li, E.C.Regan, D.Wang, W.Zhao, S.Kahn, K.Yumigeta, M.Blei, T.Taniguchi, K.Watanabe, S.Tongay, A.Zettl, M.F.Crommie, and F.Wang, Imaging two-dimensional generalized Wigner crystals, Nature 597, 650 (2021).
  * (14) L.Ma, P.X.Nguyen, Z.Wang, Y.Zeng, K.Watanabe, T.Taniguchi, A.H.MacDonald, K.F.Mak, and J.Shan, Strongly correlated excitonic insulator in atomic double layers, Nature 598, 585 (2021).
  * (15) I.V.Bondarev, O.L.Berman, R.Ya.Kezerashvili, and Y.E.Lozovik, Crystal phases of charged interlayer excitons in van der Waals heterostructures, Communications Physics 4, 134 (2021).
  * (16) A.Boltasseva and V.M.Shalaev, Transdimensional photonics, ACS Photonics 6, 1 (2019).
  * (17) I.V.Bondarev, H.Mousavi, and V.M.Shalaev, Transdimensional epsilon-near-zero modes in planar plasmonic nanostructures, Phys. Rev. Research 2, 013070 (2020).
  * (18) K.F.Mak and J.Shan, Photonics and optoelectronics of 2D semiconductor transition metal dichalcogenides, Nature Photon. 10, 216 (2016).
  * (19) P.Rivera, J.Schaibley, A.M.Jones, J.S.Ross, S.Wu, G.Aivazian, P.Klement, K.Seyler, G.Clark, N.J.Ghimire, J. Yan, D.G.Mandrus, W.Yao, and X.Xu, Observation of long-lived interlayer excitons in monolayer MoSe2-WSe2 heterostructures, Nature Commun. 6, 6242 (2015).
  * (20) J.S.Ross, P.Rivera, J.Schaibley, E.Lee-Wong, H.Yu, T.Taniguchi, K.Watanabe, J.Yan, D.Mandrus, D.Cobden, W.Yao, and X.Xu, Interlayer exciton optoelectronics in a 2D heterostructure $p$-$n$ junction, Nano Lett. 17, 638 (2017).
  * (21) M.Baranowski, A.Surrente, L.Klopotowski, J.M.Urban, N.Zhang, D.K.Maude, K.Wiwatowski, S.Mackowski, Y.C.Kung, D.Dumcenco, A.Kis, and P.Plochocka, Probing the interlayer exciton physics in a MoS2/MoSe2/ MoS2 van der Waals heterostructure, Nano Lett. 17, 6360 (2017).
  * (22) B.Miller, A.Steinhoff, B.Pano, J.Klein, F.Jahnke, A.Holleitner, and U.Wurstbauer, Long-lived direct and indirect interlayer excitons in van der Waals heterostructures, Nano Lett. 17, 5229 (2017).
  * (23) I.V.Bondarev and M.R.Vladimirova, Complexes of dipolar excitons in layered quasi-two-dimensional nanostructures, Phys. Rev. B 97, 165419 (2018).
  * (24) O.Witham, R.J.Hunt, and N.D.Drummond, Stability of trions in coupled quantum wells modeled by two-dimensional bilayers, Phys. Rev. B 97, 075424 (2018).
  * (25) T.Deilmann and K.S.Thygesen, Interlayer trions in the MoS2/WS2 van der Waals heterostructure, Nano Lett. 18, 1460 (2018).
  * (26) Z.Sun, J.Beaumariage, Q.Wan, H.Alnatah, N.Hougland, J.Chisholm, Q.Cao, K.Watanabe, T.Taniguchi, B.M.Hunt, I.V.Bondarev, and D.Snoke, Charged bosons made of fermions in bilayer structures with strong metallic screening, Nano Lett. 21, 7669 (2021).
  * (27) L.A.Jauregui, A.Y.Joe, K.Pistunova, D.S.Wild, A.A. High, Y.Zhou, G.Scuri, K.De Greve, A.Sushko, C.-H.Yu, T.Taniguchi, K.Watanabe, D.J.Needleman, M.D.Lukin, H.Park, and P.Kim, Electrical control of interlayer exciton dynamics in atomically thin heterostructures, Science 366, 870 (2019).
  * (28) E.V.Calman, L.H.Fowler-Gerace, D.J.Choksy, L.V.Butov, D.E.Nikonov, I.A.Young, S.Hu, A.Mishchenko, and A.K.Geim, Indirect excitons and trions in MoSe2/WSe2 van der Waals heterostructures, Nano Lett. 20, 1869 (2020).
  * (29) X.Wang, J.Zhu, K.L.Seyler, P.Rivera, H.Zheng, Y.Wang, M.He, T.Taniguchi, K.Watanabe, J.Yan, D.G.Mandrus, D.R.Gamelin, W.Yao, and X.Xu, Moiré trions in MoSe2/WSe2 heterobilayers, Nature Nanotechn. 16, 1208 (2021).
  * (30) E.Wigner, On the interaction of electrons in metals, Phys. Rev. 46, 1002 (1934).
  * (31) P.M.Platzman and H.Fukuyama, Phase diagram of the two-dimensional electron liquid, Phys. Rev. B 10, 3150 (1974).
  * (32) R.H.Morf, Temperature dependence of the shear modulus and melting of the two-dimensional electron solid, Phys. Rev. Lett. 43, 931 (1979).
  * (33) Y.E.Lozovik and V.I.Yudson, Crystallization of a two-dimensional electron gas in a magnetic field, JETP Lett. 22, 11 (1975).
  * (34) S.A.Khrapak, Lindemann melting criterion in two dimensions, Phys. Rev. Research 2, 012040(R) (2020).
  * (35) D.R.Musin and Yu.E.Lozovik, Stability region of the two-dimensional electron crystal in the strong magnetic fields, Phys. Stat. Sol. B 104, 712 (1981).
  * (36) V.M.Bedanov, G.V.Gadiyak, and Yu.E.Lozovik, On a modified Lindemann-like criterion for 2D melting, Phys. Lett. A 109, 289 (1985).
  * (37) Yu.E.Lozovik, V.M.Farztdinov, and B.Abdullaev, 2D electron crystal in quantized magnetic fields: melting induced by zero-point oscillations, J. Phys. C 18, L807 (1985).
  * (38) Yu.E.Lozovik and V.M.Farztdinov, Oscillation spectra and phase diagram of 2D electron crystal: ”New” (3+4)-self-consistent approximation, Solid State Commun. 54, 725 (1985).
  * (39) J.M.Kosterlitz and D.J.Thouless, Ordering, metastability and phase transitions in two-dimensional systems, J. Phys. C: Solid State Phys. 6, 1181 (1973).
  * (40) V.M.Bedanov, G.V.Gadiyak, and Yu.E.Lozovik, Melting of two-dimensional crystals, Sov. Phys. JETP 61, 967 (1985).
  * (41) C.C.Grimes and G.Adams, Evidence for a liquid-to-crystal phase transition in a classical two-dimensional sheet of electrons, Phys. Rev. Lett. 42, 795 (1979).
  * (42) C.C.Grimes and G.Adams, Crystallization of electrons on the surface of liquid helium, Surf. Sci. 98, 1 (1980).
  * (43) E.E.Mendez, M.Heiblum, L.L.Chang, and L.Esaki, High-magnetic-field transport in a dilute two-dimensional electron gas, Phys. Rev. B 28, 4886(R) (1983).
  * (44) K.Kajita, Wigner crystallization of two-dimensional electrons formed on the surface of solid neon, J. Phys. Soc. Jpn. 54, 4092 (1985).
  * (45) Yu.P.Monarkha and V.E.Syvokon, A two-dimensional Wigner crystal (Review Article), Low Temp. Phys. 38, 1067 (2012).
  * (46) H.Deng, L.N.Pfeiffer,1 K.W.West, K.W.Baldwin, L.W.Engel, and M.Shayegan, Probing the melting of a two-dimensional quantum Wigner crystal via its screening efficiency, Phys. Rev. Lett. 122, 116601 (2019).
  * (47) M.K.Ma, K.A.Villegas Rosales, H.Deng, Y.J.Chung, L.N.Pfeiffer, K.W.West, K.W.Baldwin, R.Winkler, and M.Shayegan, Thermal and quantum melting phase diagrams for a magnetic-field-induced Wigner solid, Phys. Rev. Lett. 125, 036601 (2020).
  * (48) G.Wang, A.Chernikov, M.M.Glazov, T.F.Heinz, X.Marie, T.Amand, and B.Urbaszek, Colloquium: Excitons in atomically thin transition metal dichalcogenides, Rev. Mod. Phys. 90, 021001 (2018).
  * (49) R.K.Pathria and P.D.Beale, _Statistical Mechanics_ (Elsevier, Amsterdam, 2011).
  * (50) E.S.Abers, _Quantum Mechanics_ (Addison-Wesley, New York, 2004).
  * (51) F.A.Rasmussen and K.S.Thygesen, Computational 2D materials database: electronic structure of transition-metal dichalcogenides and oxides, J. Phys. Chem. C 119, 13169 (2015).
  * (52) L.D.Landau and E.M.Lifshitz, _Quantum Mechanics. Non-Relativistic Theory_ (Pergamon, Oxford, 1991)
  * (53) G.D.Mahan, _Many-particle physics_ (Kluwer Academic, New York, 2000).
  * (54) N.H.Kleiner, L.M.Roth, and S.H.Autler, Bulk solution of Ginzburg-Landau equations for Type-II superconductors: Upper critical field region, Phys. Rev. 133, A1226 (1964).
  * (55) We use the convention by which $\nu_{e}$ of electrons and $\nu_{h}$ of holes are counted from $0$ upward and from $1$ downward, respectively, whereby $0\leq\nu_{e,h}\leq 1$ and $\nu_{e}+\nu_{h}=1$ as required by the electron-hole interchange symmetry.
  * (56) P.Nagler, M.V.Ballottin, A.A.Mitioglu, F.Mooshammer, N.Paradiso, C.Strunk, R.Huber, A.Chernikov, P.C.M.Christianen, C.Schüller, and T.Korn, Giant magnetic splitting inducing near-unity valley polarization in van der Waals heterostructures, Nature Commun. 8, 1551 (2017).
  * (57) J.Choi, W.-T.Hsu, L.-S.Lu, L.Sun, H.-Y.Cheng, M.-H.Lee, J.Quan, K.Tran, C.-Y.Wang, M.Staab, K.Jones, T.Taniguchi, K.Watanabe, M.-W.Chu, S.Gwo, S.Kim, C.-K.Shih, X.Li, and W.-H.Chang, Moiré potential impedes interlayer exciton diffusion in van der Waals heterostructures, Science Adv. 6, eaba8866 (2020).
  * (58) S.Shabani, D.Halbertal, W.Wu, M.Chen, S.Liu, J.Hone, W.Yao, D.N.Basov, X.Zhu, and A.N.Pasupathy, Deep moiré potentials in twisted transition metal dichalcogenide bilayers, Nature Physics 17, 720 (2021).
  * (59) D.R.Hofstadter, Energy levels and wave functions in rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).
  * (60) M.Ya.Azbel, Energy spectrum of a conduction electron in a magnetic field, Sov. Phys. JETP 19, 634 (1964).
  * (61) F.H.Claro and G.H.Wannier, Magnetic subband structure of electrons and hexagonal lattices, Phys. Rev. B 19, 6068 (1979).
  * (62) M.-C.Chang and Q.Niu, Berry phase, hyperorbits, and the Hofstadter spectrum: Semiclassical dynamics in magnetic Bloch bands, Phys. Rev. B 53, 7010 (1996).
  * (63) A.V.Rozhkov, A.O.Sboychakov, A.L.Rakhmanov, and F.Nori, Electronic properties of graphene-based bilayer systems, Phys. Reports 648, 1 (2016).
  * (64) A.Kormányos, G.Burkard, M.Gmitra, J.Fabian, V.Zólyomi, N.D.Drummond, and V.Fal′ko, Corrigendum: k.p theory for two-dimensional transition metal dichalcogenide semiconductors (2015 2D Mater. 2, 022001), 2D Materials 2, 049501 (2015).
  * (65) S.Larentis, H.C.P.Movva, B.Fallahazad, K.Kim, A.Behroozi, T.Taniguchi, K.Watanabe, S.K.Banerjee, and E.Tutuc, Large effective mass and interaction-enhanced Zeeman splitting of K-valley electrons in MoSe2, Phys. Rev. B 97, 201407R (2018).