# Supercell formation in epitaxial rare-earth ditelluride thin films

Adrian Llanos  Department of Applied Physics and Materials Science, California
Institute of Technology, Pasadena, California 91125, USA.  Institute for
Quantum Information and Matter, California Institute of Technology, Pasadena,
California 91125, USA.  Salva Salmani-Rezaie  School of Applied and
Engineering Physics, Cornell University, Ithaca, New York 14853, USA  Kavli
Institute at Cornell for Nanoscale Science, Cornell University, Ithaca, New
York 14853, USA.  Jinwoong Kim  Department of Physics and W. M. Keck
Computational Materials Theory Center, California State University,
Northridge, Northridge, California 91330, USA.  Nicholas Kioussis  Department
of Physics and W. M. Keck Computational Materials Theory Center, California
State University, Northridge, Northridge, California 91330, USA.  David A.
Muller  School of Applied and Engineering Physics, Cornell University, Ithaca,
New York 14853, USA  Kavli Institute at Cornell for Nanoscale Science, Cornell
University, Ithaca, New York 14853, USA.  Joseph Falson<EMAIL_ADDRESS>Department of Applied Physics and Materials Science, California Institute of
Technology, Pasadena, California 91125, USA.  Institute for Quantum
Information and Matter, California Institute of Technology, Pasadena,
California 91125, USA.

###### Abstract

Square net tellurides host an array of electronic ground states and commonly
exhibit charge-density-wave ordering. Here we report the epitaxy of DyTe2-δ on
atomically flat MgO (001) using molecular beam epitaxy. The films are single
phase and highly oriented as evidenced by transmission electron microscopy and
X-ray diffraction measurements. Epitaxial strain is evident in films and is
relieved as the thickness increases up to a value of approximately 20 unit
cells. Diffraction features associated with a supercell in the films are
resolved which is coupled with Te-deficiency. First principles calculations
attribute the formation of this defect lattice to nesting conditions in the
Fermi surface, which produce a periodic occupancy of the conducting Te square-
net, and opens a band gap at the chemical potential. This work establishes the
groundwork for exploring the role of strain in tuning electronic and
structural phases of epitaxial square-net tellurides and related compounds.

## Introduction

Layered compounds sharing the square-net structural motif support a myriad of
electronic phases including superconductivityWen-Hao et al. (2014),
topological protected modes,Klemenz et al. (2019) magnetism Sankar et al.
(2019); Lei et al. (2021) , and charge density waves (CDW)Ru et al. (2008);
Shin et al. (2005); Lei et al. (2019).Among this class of materials, rare-
earth tellurides (LnTex, where Ln is a lanthanide and $x$ = 2, 3 most
commonly), have attracted attention due to the interplay between itinerant
carriersLei et al. (2020) and magnetism, superconductivity and CDW-
order.Yumigeta et al. (2021) Sandwiched between corrugated quasi-ionic Ln-Te
layers, the electronic structure of the covalently bonded two-dimensional
square Te sheet can conceptualized by a partially filled p-orbital manifold
which obtains extra charge from an imbalance in valence states within the
corrugated layer Tremel and Hoffmann (1987); Hoffmann (1988). Modifying Ln
acts to exert chemical pressure and influences both the electronic properties
of the material. For example, in the case of the rare-earth tritellurides,
chemical pressure acts to modify the CDW wavevector $q$, the CDW-ordering
temperature, and the orthorhombicity of the crystal, which ultimately reveals
two orthogonal CDW with dissimilar ordering temperatures.Ru et al. (2008);
Yumigeta et al. (2021) Ultra-fast optical probes have demonstrated transient
expression of competing CDW orders, suggesting a complex energy landscape
governs their stability.Kogar et al. (2020) Studies on the effect of
hydrostatic pressure have revealed that CDW order can be further suppressed
and superconductivity being promoted, all in the presence of rare-earth-
derived (antiferromagnetic) ordering.Yumigeta et al. (2021).

Figure 1: Crystal structure of DyTe2 on MgO plotted along the (a) $a$ axis and
(b) $c$ axis directions. (c) Schematic of film growth recipe indicating
substrate temperature as a function of time. Color bars at the top of the
panel indicate when the Dy and Te effusion cells are open. (d) RHEED
oscillations from the specular reflection of films as a function of growth
time. (e) AFM image of a $T_{\mathrm{g}}=330^{\circ}$C film of 20 u.c. RHEED
of (f) thermally prepared substrate, (g) after buffer layer deposition, and
(h) at the conclusion of growth, all recorded along the [001] direction of the
MgO substrate.

The related LnTe2 display similar behaviors to the tritellurides, but differ
in key ways due to their electron filling of 7e-/Te in the square planar
layers Tremel and Hoffmann (1987); Kikuchi (1998). In contrast to the
tritellurides, which form stoichiometric crystals, the ditellurides have a
substantial width of formation and tend to form Te-deficient crystal
structures with a range of stoichimetries and resulting superstructure
modulations Niihara et al. (1972); Park et al. (1998); Shin et al. (2005);
Poddig et al. (2018). These modulations are often commensurate with the parent
structure and result from an ordering of Te vacancies within the Te sheet.
However, incommensurate modulations have also been observed in LaTe2-δ and
CeTe2-δ Shin et al. (2005). Electrical transport ranging from metallic Stöwe
(2000) to semiconducting Jung et al. (1997); Kwon and Min (2000) for similar
nominal compositions and photoemission studies reporting both pseudogap Shin
et al. (2005); Chung et al. (1999) and fully-gapped Garcia et al. (2007) Fermi
surface structures further demonstrate the complexity of the ditelluride
system. Superconductivity has also been reported in CeTe2 at 0.2 GPaJung et
al. (2003) and in Sb-doped LaTe2 at approximately 3 GPaChen et al. (2022). The
subtle dependence on defect stoichimetry and applied pressure motivate an
investigation into the role that epitaxial strain and dimensionality play in
the stabilization of lattice distortions and the resulting groundstate of the
electron-phonon system.

Here we report the epitaxial growth of rare-earth ditelluride DyTe2-δ on MgO
(001) substrates. While little is known about DyTe2-δ in the bulk, the closely
related structures of LaTe2-δ and CeTe2-δ Shin et al. (2005) are are reported
to either have tetragonal or orthorhomicDiMasi et al. (1996) symmetry
belonging to the $P4/nmm$ and $Pnma$ space groups, respectively. We have found
in all experiments that the ditelluride phase is stabilized preferentially
over the tritelluride phase despite synthesis being performed in Te-rich
conditions. We attribute this to the weak van-der Waals interactions between
stacked Te-sheets on the surface of samples during layer-by-layer growth.
While growth of the tritelluride phase has been reported in the ultra-thin
limit Xu et al. (2021) ($<$ 1 unit cell (u.c.)), this appears to not persist
into multilayer films, even as thin as 3 u.c. Our results identify that
epitaxial strain is evident in films and is dependent on the number of
quintuple layers (QL) deposited, being progressively relieved as the films are
made thicker. The films are tetragonal within the detection limit of the X-ray
diffraction apparatus. We additionally have resolved a commensurate
superlattice in a subset of films with wave vectors
$\vec{q_{1}}=\frac{2}{5}\vec{a}^{*}+\frac{1}{5}\vec{b}^{*}$ ,
$\vec{q_{2}}=-\frac{1}{5}\vec{a}^{*}+\frac{2}{5}\vec{b}^{*}$ textand
$\vec{q_{3}}=\frac{1}{2}\vec{c}^{*}$ where $\vec{a}^{*},\vec{b}^{*}\text{and}\
\vec{c}^{*}$ are the reciprocal lattice vectors of the undistorted structure.
In real space, this distortion can be thought of as
($\sqrt{5}\times\sqrt{5}$)R26.6∘ $\times$ 2 superlattice. In the case of
Ref.16, the distortion was associated with a reduction in symmetry of the
compound from the $P4/nmm$ to $P4/n$ space group. The superlattice emerges due
to a periodic Te-defect lattice, which also acts to open a gap in the
electronic spectrum and induce semiconducting transport behaviour. First
principles calculations point towards nesting conditions of the Fermi surface
at a $q$-vector that corresponds to the $\sqrt{5}\times\sqrt{5}$ condition,
suggesting that the formation of the defect lattice results from a similar
driving force to the conventional picture of CDW formation where sections of
the Fermi surface are gapped out by the formation of supercells with
periodicity corresponding to the nesting conditionLee and Foran (1996). It is
therefore a thermodynamically favored process rather than an outcome of
imperfect growth conditions. The calculations further reveal that a number of
competing atomic arrangements may be present and dependent on the degree of
strain in the crystal.

## Experimental Methods

Substrates are prepared by annealing in ultra-high vacuum (base pressure
$\approx 10^{-10}$ Torr) using a home-built CO2 laser heating apparatus,
inspired by Ref.27. We used MgO (001) (Crystec) substrates for growth without
back coating of an absorbing layer or additional cleaning steps. The
combination of MgO ($a$=0.421 nm) and DyTe2 ($a$=0.428 nm) was selected due to
the relatively good epitaxial relationship at approximately -1.6% mismatch, as
shown in Fig.1(a) and (b) where we highlight the DyTe2 unit cell and its
constituent atoms; sheet tellurium (Tes), the corrugated layer tellurium
(Tec), and the rare-earth ion. Being hygroscopic, the MgO surface without
annealing is of poor quality, but upon heating to in excess of 1000 ∘C we can
generate step-terraced structures within a few minutesBraun et al. (2020). The
substrate is then transferred to a dedicated telluride substrate holder in a
glove box, and then reloaded into the vacuum cluster and transferred to the
molecular beam epitaxy growth chamber. The reflected high energy electron
diffraction (RHEED) signal from the surface remains qualitatively the
identical upon arrival into the growth chamber. Substrate heating in the
growth chamber is provided by a conventional resistive SiC coil, and the
temperature reported is read from a thermocouple between the substrate and
heating coil. Te was supplied using a thermal cracker cell with tank
temperatures of approximately 450 ∘C and cracking zone temperature of 1000∘C.
Dysprosium was evaporated from a conventional effusion cell at a temperature
of 1075 ∘C. We have found that both thermal cracking of Te and thermal
preparation of substrates play a role in improving sample quality. Due to the
air sensitivity of the films, we cap each film before extracting them from the
vacuum chamber with either amorphous Ge or partially crystalline Te. This is
performed in-situ at $T\approx$90 K using a liquid nitrogen cooling stage.
X-ray diffraction data was obtained using a Rigaku Smartlab diffractometer
using a 2-bound Ge (220) monochromator. High resolution diffraction employed
copper $K\alpha 1$ radiation. Diffraction in the (hk0) plane (i.e. the
$\chi$-circle) was accessed using an in-plane diffraction arm in an incidence
angle $\omega~{}0.35^{\circ}$. A $0.5^{\circ}$ parallel slit collimator and
$0.5^{\circ}$ parallel slit analyzer were used to reduce in-plane beam
divergence.

The cross-sectional Scanning Transmission Electron microscopy (STEM) samples
are prepared using a standard lift-out process on Helios G4 UX DualBeam
focused ion beam system. We utilize the Spectra 300 X-CFEG operating at 200 kV
with a semi-convergence angle of 30mrad and a High-Angle Annular Dark-Field
(HAADF) detector with an angular range of 60-200 mrad to collect HAADF-STEM
images. The HAADF-STEM images are acquired as a series of 25 images (250 ns
per frame) and subsequently averaged to produce images with a high signal-to-
noise ratio. STEM- Energy dispersive X-ray spectroscopy (EDX) data is acquired
using a steradian Dual-X EDX detector, and the resulting spectrum is denoised
through the application of the principal component analysis (PCA). We use an
electron microscope pixel array detector (EMPAD) to collect Convergent beam
electron diffraction (CBED) patterns at a convergence angle of 1.2 mrad and
dwell time of 1 ms. For the CBED simulation, we use the Py4DSTEM package based
on the Bloch wave method.Savitzky et al. (2021)

## Calculation Methods

The density functional theory calculations are performed using the Vienna $ab$
$intio$ simulation packageKresse and Furthmüller (1996a, b) with the projector
augmented wave (PAW) method.Blöchl (1994); Kresse and Joubert (1999) The
Dy-4$f$ states are treated as core. All structures are optimized using the
PBEsol exchange correlation functionalPerdew et al. (2008) and the Bloch
states are wannierizedMostofi et al. (2014); Marzari et al. (2012) with the
modified Becke-Johnson (mBJ) potentialBecke and Johnson (2006); Tran and Blaha
(2009) which provides accurate band gaps, effective masses and frontier-band
ordering. The Fermi surfaces and nesting functions (supplementary information
(SI)) are obtained using the tetrahedron method with the Wannier Hamiltonian.

## Results and Discussion

The growth procedure is illustrated in Fig1(c). We begin with a low
temperature ($T$=200∘C) buffer layer on the annealed MgO surface. By observing
intensity oscillations of the RHEED pattern, precise thickness control can be
realized allowing for the buffer layer growth to be stopped after a single QL.
Attempts to grow films without a buffer layer resulted in poorly crystalline
films, as observed by weak, diffuse features in RHEED. Growth of the remainder
of the film takes place after annealing the buffer under Te flux to a final
growth temperature. Growth rate is set by the flux of Dy atoms and was set at
1.2 A/min for all growths under Te over pressure. Te/Dy flux ratio were
typically in the range of 10 to 20, corresponding to beam flux pressures of
7$\times 10^{-9}$ mbar for Dy and 1.4$\times 10^{-7}$ mbar for Te. We have
deposited films within a range of growth temperatures ($T_{\mathrm{g}}$), as
shown in 1(c). For films deposited at low temperature ($T_{g}<$300∘C), rough
island growth appears whereas those grown at elevated ($T_{g}>$ 350∘C) overall
suffer from poor crystallinity and rough surfaces, likely due to the low
sticking coefficient of Te. RHEED oscillations during growth for moderate
growth temperatures ($T_{g}=$ 315∘C, Fig.1(d)) indicate that the film grows
layer-by-layer with oscillations corresponding to the growth of a single unit
cell. These oscillations were observed to decay after approximately 10
oscillations and were completely damped by the completion of a 20uc thick
film. To estimate the growth rate, RHEED oscillations were observed for the
first several layers. The substrates are subsequently rotated at 5 RPM to
enhance sample uniformity. Fig.1(e) shows an atomic force microscopy image
(AFM) of the final surface of a 20 u.c. film, showing surface roughness on the
order of 2 3 u.c. RHEED images are shown of the pristine substrate (Fig.1(f)),
after buffer layer (g) and post-growth (h), the latter of which shows sharp,
streaky patterns with prominent Kikuchi bands.

Figure 2: HAADF imaging of the film from a wide angle (a) and at the interface
with MgO (b). (c-d) EDX mapping of the film showing occupancy of Dy and Te.
(f) Intensity of atoms in a vertical line cut of the HAADF data indicated as
the red bar in panel (b), showing the oscillatory intensity as a function of
distance. The Dy, Tes and Tec sites are noted.

Transmission electron microscopy analysis is presented in Fig.2 and shows
stacked layers of planar Te sheets and corrugated Dy-Te interlayers. The
stacking sequence confirms the presence of the DyTe2 phase and not DyTe3. We
attribute this to the low sticking coefficient of the Te-sheets, which
evidently cannot be stacked in pairs (as is the case in the tritellurides)
even with a Te/Dy flux ratio as high as 40. The film is phase pure with a
sharp interface formed between the MgO and Te-square net sheet (Fig.2(b)). The
square net tellurium sheet forms the first wetting layer on the substrate.
Rotated domains are evident in some areas of the film, but are suppressed
after a few monolayers of growth. The high-angle annular dark-field imaging
(HAADF) (Fig.2(a,b)) and energy-dispersive X-ray spectroscopy elemental map
(Fig.2(c-e)) shows nominal (P4/nmm) DyTe2 structure. The Dy/Te ratio with
absorption and ionization cross-section correction for the foil with a
thickness of 30 nm is 0.497. This ratio, at its best, has a 10-20% error
because of strong channeling.Lugg et al. (2015); MacArthur et al. (2017) When
the beam is along the zone axis, the intensity of x-ray spectra is not a
linear contribution of each element. The experimental CBED pattern (see SI)
matches the simulated CBED pattern of the ($P4/nmm$) point group when recorded
along the [010] direction. The high resolution images also capture dissimilar
intensity of Te sites along a line profile in the DyTe2 [001] direction, as
plotted in Fig.2(f). The intensity of the Te columns at the Te sheet (Tes) and
the Te columns at the corrugated layer (Tec) is different. We attribute this
to the phase space for vibrations of the atoms in the square sheet versus
those in the ionic Dy-Te layer; the higher thermal vibration in the Te sheets
is an indication that the atoms are loosely bound. This observation is in
agreement with our experience of being unable to stabilize the tritelluride
phase as a thin-film due to the poor sticking coefficient of Te atoms in the
square net layer.

Figure 3: (a) Out-of-plane $\theta$-2$\theta$ X-ray diffraction of a 20 u.c.
DyTe2-δ film on MgO (001), plotted in reciprocal lattice units of DyTe2-δ.
Peaks corresponding to the substrate are indicated by a star, with fundamental
peaks of the film labelled by a circle and superlattice peaks by a triangle.
(b) Rocking curve analysis of the DyTe2-δ (003) diffraction peak. (c) Out-of-
plane asymmetric scan of the film with diffraction from the MgO (204) and
DyTe2-δ (208) planes plotted as a function of momentum transfer in the
$Q_{x}$–$Q_{z}$ plane. (d) In-plane $a$ and (e) out-of-plane $c$ lattice
constants of films as a function of thickness. Dashed lines are fits intended
as guides for the eye.

The structural properties of grown films were further investigated using X-ray
diffraction. Figure 3(a) shows a $\theta$-$2\/\theta$ scan, plotted in
reciprocal lattice units, showing highly $c$-axis oriented films with intense
layer peaks accompanied by prominent Laue fringes indicating a high quality
film with a sharp interface with the underlying substrate pietsch. An $\omega$
rocking curve analysis of the DyTe2 (003) peak produces a full-width at half
maximum (FWHM) of 0.03∘ indicating low out-of-plane mosaicity. $\phi$ scan
data (not shown) around the (110) peaks taken in grazing incidence geometry
reveals the four-fold symmetry of the layer and further confirms the
$[100]_{\mathrm{DyTe_{2}}}||[100]_{\mathrm{MgO}}$ epitaxial relationship
observed with RHEED.

Figure 4: (a) In-plane (hk0) diffraction of a 20 u.c. film showing diffraction
peaks identified as fundamental (circle) and superlattice (triangle). (b)
Schematic of the Te square net layer showing the original square u.c. (upper
square) and the distorted $\sqrt{5}\times\sqrt{5}$ R26.6∘ superlattice (lower
rotated square) of partially occupied Te-sites, indicated by half-transparent
atoms.

One goal of this work is to investigate the presence of epitaxial strain and
its dependence on film thickness and growth parameters. To this end, we employ
high resolution reciprocal space mapping in the $Q_{x}$–$Q_{z}$ plane for
films with thickness between 3 and 35 u.c. A characteristic data set for a 5
unit cell film is plotted in Fig.3(c), showing the MgO (204) and DyTe2-δ (208)
diffraction spots. For the 3 unit cell film, intensity in the $Q_{x}$–$Q_{z}$
plane was too weak to fit a peak. In this case, measurements of the
$Q_{x}$–$Q_{y}$ RSM were performed about the $(hh0)$ for $odd\ h$ layer
diffraction peaks in grazing incidence geometry. These reflections were chosen
to ensure that the MgO substrate would not contribute to the measurement. The
$c$-axis lattice constants were determined from the $\theta-2\theta$
measurements. In Fig3(d) and (e) we plot the evolution of the measured $a$ and
$c$ lattice constants for films grown under a Te/Dy flux ratio of 20 and a
growth temperature $T_{g}$=310∘C. There is a continuous evolution of in-plane
compressive strain from -$1.23\%$ to $0.5\%$ relative to reported bulk values
for DyTe1.8. As will be described in the subsequent section, the presence of
the observed superlattice implies that the stoichiometry is closest to DyTe1.8
and the epitaxial strain estimates were referenced to these values Poddig et
al. (2018). The relaxed value of the in-plane lattice constant was found to be
4.26$\AA$ which. As expected for films showing in-plane compressive strain,
the $c$-axis lattice constant was observed to evolve from $+1.1\%$ at the
lowest thickness to $+0.1\%$ by 20nm. Films grown at higher temperatures
($T_{g}>$ 330∘) showed no evolution of strain with thickness.

It is known that in the case of La and Ce the ditelluride structure is stable
over a range of stoichiometries in the bulk, typically within the range of
${x}=1.8\sim 2.0$. The effects of such off-stoichiometries in the case of
LaTe1.8 leads to a gap opening at the chemical potential.Poddig et al. (2018)
The off-stoichiometries are known to produce a variety of structural
modulations depending on Te deficiency. The most commonly observed modulation,
for a Te deficiency of 20%, was a ($\sqrt{5}\times\sqrt{5}$$\times 2$)
superlattice of the $LnTe_{2}$ aristotype. The u.c. of this ordered defect
modulation is oriented parallel to the [210] direction the higher symmetry
structure. Such a defect modulation is reminiscent of that observed in the
insulating phase of Fe-deficient FeSe Liu and Zou (2020) and in certain oxygen
deficient surface reconstructions of SrTiO3 Newell et al. (2006).

Figure 5: Density functional theory calculations of DyTe2-δ: (a) Fermi surface
of pristine DyTe2 at a chemical potential corresponding to $\delta=0.2$ Te
deficiency, where q is the nesting vector. (b) In-plane nesting function,
$N_{\gamma}(q_{x},q_{y})$, (see SI) in the first Brillouin zone (solid lines),
where the eight yellow dots denote the nesting vector q, corresponding to the
$\sqrt{5}\times\sqrt{5}$ cell modulation. (c) Te bonding network for
$\delta=0.2$ where the Te dimers are colored as blue. (d) Band structure along
the symmetry directions of the $\sqrt{5}\times\sqrt{5}$ cell shown in (c)
revealing the emergence of the band gap. The orbital characters of the Dy-$d$,
Tec-$p$, and Tes-$p$ derived states in (a) and (d) are denoted by orange,
yellow, and blue, respectively.

Using grazing incidence X-ray diffraction we are able to analyze the
structural modulations present in DyTe2-δ thin films, as shown in Fig.4(a)
where we plot diffraction of a 20 unit cell film in the (hk0) plane. The
enhanced intensity in grazing incidence geometry for thin films enables use to
resolve the presence of the $\sqrt{5}\times\sqrt{5}$ modulation of the a-b
plane, indexed by the diffraction spots identified by triangles in the figure,
which are absent in the undistorted structure. Diffraction indices given by
triangles are with respect to the full supercell whereas those identified by
circles are referenced to the undistorted structure. The presence of the (310)
reflection in this region of reciprocal space indicates the presence of a
rotated domain with modulation wavevector
$\vec{q_{1}}=\frac{1}{5}\vec{a}^{*}+\frac{2}{5}\vec{b}^{*}$ ,
$\vec{q_{2}}=\frac{2}{5}\vec{a}^{*}-\frac{1}{5}\vec{b}^{*}$ textand
$\vec{q_{3}}=\frac{1}{2}\vec{c}^{*}$. As previously discussed,
$\theta-2\theta$ diffraction measurements revealed additional features at 0.5
reciprocal lattice units, which correspond to a doubling of the unit cell
along the $c$-axis. Additionally in RHEED, we are able to resolve weak streaks
associated with an in-plane modulation (see Fig.S1). Taken all together, these
diffraction features provide strong evidence for the
($\sqrt{5}\times\sqrt{5}$) R26.6∘ $\times$ 2 reconstruction being stabilized
in the films. The presence of the defect lattice would imply an off-
stoichiometry in the range of 10$\sim$20%, which is consistent with the EDX-
obtained chemical ratios within the uncertainty of the measurement. All films
with thickness greater than 10nm showed features consistent with the 2-fold
modulation of the $c$-axis. However, the comparably weaker diffracted
intensity from films in the thin-limit combined with the intrinsically weak
intensity from superlattice diffraction spots prevented the observation of the
superlattice below 20 u.c. with our lab-based diffractometer.

To further understand the underlying origin of the experimentally observed
supercell, we have performed first principles calculations with and without Te
deficiencies. Questions of interest include whether nesting within the Fermi
surface promotes the formation of the supercell, the atomic arrangement(s)
within the supercell, and the culminating band structure of the extended
compound. Figure 5 presents the main results of our calculations, and an
extended set of data is available in the SI. The Fermi surface of the
unmodulated pristine DyTe2-δ structure at the chemical potential corresponding
to $\delta=0.2$ Te deficiency is presented in Fig.5(a) showing the Tes $p$\-
and the Dy $d$-derived bands in blue and orange, respectively. The planar
shapes of the Tes derived Fermi surface confirm their 2D nature. The presence
of bands at the Fermi energy would be expected to produce metallic transport
in films. However, we observe semiconducting transport with an activation gap
on the order of 300 meV (see SI). The gapping of the Fermi surface can be
achieved through a number of means; in other square net compounds where CDW
ground states have been identified, the Fermi surface hosts a large number of
nesting vectors $q$ in momentum space, which produces enhancements in the
Lindhard response function $\chi(\mathbf{q})$ when $q\approx 2k_{\mathrm{F}}$.
The induced charge is a product of the response function and the time-
independent potential of the crystal. Thus, a highly nested Fermi surface, as
is the case in the quasi-2D compound studied here, may drive the formation of
CDW ground states.Zhu et al. (2015) The calculated nesting function, $N({\bf
q})$, (see SI) of the Fermi surface in Fig.5(a) is displayed in Fig.5(b),
where red colors indicate wavevectors with a strong nesting condition at the
Fermi energy. The yellow dots correspond to the $\mathbf{q}$-vector which
would produce a $\sqrt{5}\times\sqrt{5}$ modulation, and they are seen to
closely coincide with the nesting conditions of the Fermi surface. Therefore,
we hypothesize that the $\sqrt{5}\times\sqrt{5}$ modulation is promoted in the
structure to minimize energy. However, in contrast to purely electronic CDW-
order, it appears likely that the crystal preferentially forms a defect
lattice. Removal of a Te from the square net sheet effectively contributes
extra charge to the layer, which in turn may be localized on the periodic
defect site. The question remains what the defect lattice looks like locally
and how it influences the band structure. Consequently, we have carried out
calculations of the $\sqrt{5}\times\sqrt{5}$ supercell with Tes mono- and di-
vacancies corresponding to $\delta=0.1$ and $0.2$, respectively. We find that
the A-C divacancy configuration, consisting of Tes dimers and trimers, shown
in Fig. 5(c) has the lowest formation energy among the five considered
structures and is in agreement with previous X-ray diffraction measurements on
SmTe1.8.Ijjaali and Ibers (2006) This configuration is favored in the
calculations when the $a$-axis lattice parameter is compressed, which is a
trend that is resolved in diffraction in Fig.3(d). The band structure of the
A-C configuration along symmetry directions, shown in Fig. 5(d), demonstrates
the emergence of band gap (indirect gap of 143 meV and direct gap of 326 meV).
We also find a band gap for the mono-vacancy and A-B di-vacancy
configurations, while the A-D di-vacancy configuration with a higher formation
energy is metallic. Therefore, we speculate that these vacancy configurations,
especially A-C, are stabilized and are responsible for the semiconducting
transport behavior. Additional structure modulation along the $c$ direction is
also reproduced under the assumption of the dominant $\delta=0.2$ A-C
configuration. The calculated total energies of the
$\sqrt{5}\times\sqrt{5}\times 2$ structure (Table 1) show that the A-C
divacancy appears at a laterally displaced position $(0.5,0.5)$, relative to
the divacancy on the adjacent Te square net layer. This lateral displacement
lowers the total energy by approximately 10 meV/Dy compared to the
$\sqrt{5}\times\sqrt{5}\times 1$ structure, where the A-C dimers align as a
column.

## Conclusion

In summary, we have grown epitaxial DyTe2-δ on atomically flat MgO surfaces
which induces compressive strain in the films. The in-plane lattice constant
of the films increases as the thickness is increased, indicating a relieving
of strain. The high structural quality of films enables detection of a
$\sqrt{5}\times\sqrt{5}$$\times 2$ R26.6∘ superlattice in diffraction. First
principles calculations point towards nesting in the Fermi surface as a
possible mechanism for promoting the supercell modulation, which in turn
produces a band gap and semiconducting transport. This work sets the stage for
studying the role epitaxial strain has on the broader class of LnTe2-δ and
related square net compounds, and ultimately their heterointerfaces.

### Acknowledgements

We thank Kaveh Pezeshki for help in the early stages of this project and
appreciate discussions with Ian Fisher, Leslie Schoop and Anshul Kogar. J.F
acknowledges funding provided by the the Air Force Office of Scientific
Research (Contract number FA9550-22-1-0463), and partial funding provided by
the Gordon and Betty Moore Foundation’s EPiQS Initiative (Grant number
GBMF10638), and the Institute for Quantum Information and Matter, an NSF
Physics Frontiers Center (NSF Grant PHY-1733907). N.K and J.F. acknowledge
support from the NSF-PREP CSUN/Caltech-IQIM Partnership (Grant number 2216774)
and NSF-PREM (Grant number DMR-1828019). This work made use of a Helios FIB
supported by NSF (Grant No. DMR-1539918) and the Cornell Center for Materials
Research (CCMR) Shared Facilities, which are supported through the NSF MRSEC
Program (Grant No. DMR-1719875). The Thermo Fisher Spectra 300 X-CFEG was
acquired with support from PARADIM (NSF MIP DMR-2039380) and Cornell
University. S.S-R. acknowledges the support of the Kavli and PARADIM
fellowship.

##

## References

  * Wen-Hao et al. (2014) Z. Wen-Hao, S. Yi, Z. Jin-Song, L. Fang-Sen, G. Ming-Hua, Z. Yan-Fei, Z. Hui-Min, P. Jun-Ping, X. Ying, W. Hui-Chao, et al., Chinese Physics Letters 31, 017401 (2014), ISSN 0256-307X, eprint 1311.5370.
  * Klemenz et al. (2019) S. Klemenz, S. Lei, and L. M. Schoop, Annual Review of Materials Research 49, 185 (2019), eprint https://doi.org/10.1146/annurev-matsci-070218-010114, URL https://doi.org/10.1146/annurev-matsci-070218-010114.
  * Sankar et al. (2019) R. Sankar, I. P. Muthuselvam, K. R. Babu, G. S. Murugan, K. Rajagopal, R. Kumar, T.-C. Wu, C.-Y. Wen, W.-L. Lee, G.-Y. Guo, et al., Inorganic Chemistry 58, 11730 (2019), URL https://doi.org/10.1021/acs.inorgchem.9b01698.
  * Lei et al. (2021) S. Lei, S. M. L. Teicher, A. Topp, K. Cai, J. Lin, G. Cheng, T. H. Salters, F. Rodolakis, J. L. McChesney, S. Lapidus, et al., Advanced Materials 33, 2101591 (2021), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202101591, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/adma.202101591.
  * Ru et al. (2008) N. Ru, C. L. Condron, G. Y. Margulis, K. Y. Shin, J. Laverock, S. B. Dugdale, M. F. Toney, and I. R. Fisher, Phys. Rev. B 77, 035114 (2008), URL https://link.aps.org/doi/10.1103/PhysRevB.77.035114.
  * Shin et al. (2005) K. Y. Shin, V. Brouet, N. Ru, Z. X. Shen, and I. R. Fisher, Phys. Rev. B 72, 085132 (2005), URL https://link.aps.org/doi/10.1103/PhysRevB.72.085132.
  * Lei et al. (2019) S. Lei, V. Duppel, J. M. Lippmann, J. Nuss, B. V. Lotsch, and L. M. Schoop, Advanced Quantum Technologies 2, 1900045 (2019), ISSN 2511-9044.
  * Lei et al. (2020) S. Lei, J. Lin, Y. Jia, M. Gray, A. Topp, G. Farahi, S. Klemenz, T. Gao, F. Rodolakis, J. L. McChesney, et al., Science Advances 6, eaay6407 (2020), eprint https://www.science.org/doi/pdf/10.1126/sciadv.aay6407, URL https://www.science.org/doi/abs/10.1126/sciadv.aay6407.
  * Yumigeta et al. (2021) K. Yumigeta, Y. Qin, H. Li, M. Blei, Y. Attarde, C. Kopas, and S. Tongay, Advanced Science 8, 2004762 (2021), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/advs.202004762, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/advs.202004762.
  * Tremel and Hoffmann (1987) W. Tremel and R. Hoffmann, Journal of the American Chemical Society 109, 124 (1987), ISSN 0002-7863.
  * Hoffmann (1988) R. Hoffmann, _Solids and Surfaces: A Chemist’s View of Bonding in Extended Structures_ (VCH Publishers, 1988), ISBN 9783527269051, URL https://books.google.com/books?id=Y--LO5YSU24C.
  * Kogar et al. (2020) A. Kogar, A. Zong, P. E. Dolgirev, X. Shen, J. Straquadine, Y.-Q. Bie, X. Wang, T. Rohwer, I.-C. Tung, Y. Yang, et al., Nature Physics 16, 159 (2020).
  * Kikuchi (1998) A. Kikuchi, Journal of the Physical Society of Japan 67, 1308 (1998), ISSN 0031-9015.
  * Niihara et al. (1972) K. Niihara, T. Shishido, and S. Yajima, Bulletin of the Chemical Society of Japan 45, 1239 (1972), ISSN 0009-2673.
  * Park et al. (1998) S.-M. Park, S.-J. Park, and S.-J. Kim, Journal of Solid State Chemistry 140, 300 (1998), ISSN 0022-4596, URL https://www.sciencedirect.com/science/article/pii/S0022459698978914.
  * Poddig et al. (2018) H. Poddig, T. Donath, P. Gebauer, K. Finzel, M. Kohout, Y. Wu, P. Schmidt, and T. Doert, Zeitschrift für anorganische und allgemeine Chemie 644, 1886 (2018), ISSN 0044-2313.
  * Stöwe (2000) K. Stöwe, Journal of Solid State Chemistry 149, 155 (2000), ISSN 0022-4596.
  * Jung et al. (1997) M. Jung, Y. Kwon, and T. Suzuki, Physica B: Condensed Matter 240, 83 (1997), ISSN 0921-4526.
  * Kwon and Min (2000) Y. Kwon and B. Min, Physica B: Condensed Matter 281, 120 (2000), ISSN 0921-4526.
  * Chung et al. (1999) J. Chung, J. Park, J. G. Park, B.-H. Choi, S. J. Oh, E. J. Cho, H. D. Kim, and Y. S. Kwon, arXiv (1999), eprint cond-mat/9910443.
  * Garcia et al. (2007) D. R. Garcia, G.-H. Gweon, S. Y. Zhou, J. Graf, C. M. Jozwiak, M. H. Jung, Y. S. Kwon, and A. Lanzara, Physical Review Letters 98, 166403 (2007), ISSN 0031-9007, eprint cond-mat/0703535.
  * Jung et al. (2003) M. H. Jung, A. Alsmadi, H. C. Kim, Y. Bang, K. H. Ahn, K. Umeo, A. H. Lacerda, H. Nakotte, H. C. Ri, and T. Takabatake, Phys. Rev. B 67, 212504 (2003), URL https://link.aps.org/doi/10.1103/PhysRevB.67.212504.
  * Chen et al. (2022) X. Chen, P. han Sun, Z. Xie, F. Meng, C. Pei, Y. Qi, T. Ying, K. Liu, J. gang Guo, and X. Chen (2022), eprint 2209.07188.
  * DiMasi et al. (1996) E. DiMasi, B. Foran, M. C. Aronson, and S. Lee, Phys. Rev. B 54, 13587 (1996), URL https://link.aps.org/doi/10.1103/PhysRevB.54.13587.
  * Xu et al. (2021) Z. Xu, S.-H. Ji, L. Tang, J. Wu, N. Li, X. Cai, and X. Chen, Chinese Physics Letters 38, 077102 (2021), URL https://dx.doi.org/10.1088/0256-307X/38/7/077102.
  * Lee and Foran (1996) S. Lee and B. J. Foran, Journal of the American Chemical Society 118, 9139 (1996), ISSN 0002-7863.
  * Braun et al. (2020) W. Braun, M. Jäger, G. Laskin, P. Ngabonziza, W. Voesch, P. Wittlich, and J. Mannhart, APL Materials 8 (2020), ISSN 2166-532X, 071112, eprint https://pubs.aip.org/aip/apm/article-pdf/doi/10.1063/5.0008324/14564941/071112_1_online.pdf, URL https://doi.org/10.1063/5.0008324.
  * Savitzky et al. (2021) B. H. Savitzky, S. E. Zeltmann, L. A. Hughes, H. G. Brown, S. Zhao, P. M. Pelz, T. C. Pekin, E. S. Barnard, J. Donohue, L. Rangel DaCosta, et al., Microscopy and Microanalysis 27, 712 (2021), ISSN 1431-9276, eprint https://academic.oup.com/mam/article-pdf/27/4/712/48292107/mam0712.pdf, URL https://doi.org/10.1017/S1431927621000477.
  * Kresse and Furthmüller (1996a) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996a), URL https://link.aps.org/doi/10.1103/PhysRevB.54.11169.
  * Kresse and Furthmüller (1996b) G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996b), ISSN 0927-0256, URL http://www.sciencedirect.com/science/article/pii/0927025696000080.
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994), URL https://link.aps.org/doi/10.1103/PhysRevB.50.17953.
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999), URL https://link.aps.org/doi/10.1103/PhysRevB.59.1758.
  * Perdew et al. (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.100.136406.
  * Mostofi et al. (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014), ISSN 0010-4655, URL http://www.sciencedirect.com/science/article/pii/S001046551400157X.
  * Marzari et al. (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012), URL http://link.aps.org/doi/10.1103/RevModPhys.84.1419.
  * Becke and Johnson (2006) A. D. Becke and E. R. Johnson, The Journal of Chemical Physics 124 (2006), ISSN 0021-9606, 221101, eprint https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2213970/15385734/221101_1_online.pdf, URL https://doi.org/10.1063/1.2213970.
  * Tran and Blaha (2009) F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.102.226401.
  * Lugg et al. (2015) N. Lugg, G. Kothleitner, N. Shibata, and Y. Ikuhara, Ultramicroscopy 151, 150 (2015), ISSN 0304-3991, special Issue: 80th Birthday of Harald Rose; PICO 2015 – Third Conference on Frontiers of Aberration Corrected Electron Microscopy, URL https://www.sciencedirect.com/science/article/pii/S0304399114002459.
  * MacArthur et al. (2017) K. E. MacArthur, H. G. Brown, S. D. Findlay, and L. J. Allen, Ultramicroscopy 182, 264 (2017), ISSN 0304-3991, URL https://www.sciencedirect.com/science/article/pii/S0304399117302668.
  * Liu and Zou (2020) C. Liu and K. Zou, Physical Review B 101, 140502 (2020), ISSN 2469-9950, eprint 2002.06314.
  * Newell et al. (2006) D. T. Newell, A. Harrison, F. Silly, and M. R. Castell, Physical Review B 75, 205429 (2006), ISSN 1098-0121.
  * Zhu et al. (2015) X. Zhu, Y. Cao, J. Zhang, E. W. Plummer, and J. Guo, Proceedings of the National Academy of Sciences 112, 2367 (2015), eprint https://www.pnas.org/doi/pdf/10.1073/pnas.1424791112, URL https://www.pnas.org/doi/abs/10.1073/pnas.1424791112.
  * Ijjaali and Ibers (2006) I. Ijjaali and J. A. Ibers, Journal of Solid State Chemistry 179, 3456 (2006), ISSN 0022-4596, URL https://www.sciencedirect.com/science/article/pii/S0022459606003926.
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976), URL https://link.aps.org/doi/10.1103/PhysRevB.13.5188.
  * Shim et al. (2004) J. H. Shim, J.-S. Kang, and B. I. Min, Phys. Rev. Lett. 93, 156406 (2004), URL https://link.aps.org/doi/10.1103/PhysRevLett.93.156406.

### Supplementary Data

Figure S1: Off-axis RHEED image with saturated brightness showing weak
additional streaks which we tentatively associated with the
$\sqrt{5}\times\sqrt{5}$ modulation. Figure S2: RHEED images taken at the end
of growths performed in the indicated temperature windows. Figure S3:
Temperature dependent transport of a 20 u.c. film showing semiconducting
transport. Figure S4: (a) CBED pattern which matches the (b) simulated CBED
pattern of the (p4/nmm) point group. {0k0} (k=odd) reflections, marked by the
arrow, are kinematically forbidden.

####  .0.1 Calculation Methods

The density functional theory calculations are performed using the Vienna $ab$
$initio$ simulation package.Kresse and Furthmüller (1996a, b) The pseudo-
potentials are treated with the projector augmented wave (PAW) methodBlöchl
(1994); Kresse and Joubert (1999) with valence configurations of
$6s^{2}5p^{6}5d^{1}$ for Dy and $5s^{2}5p^{4}$ for Te. The plane-wave cut-off
energy is set to 300 eV. All structures are optimized with the PBEsol
functionalPerdew et al. (2008) with an energy convergence criteria of
($<10^{-6}$ eV) where the Brillouin zone is sampled with Monkhorst-Pack
meshesMonkhorst and Pack (1976) of $12\times 12\times 6$ for pristine DyTe2,
$6\times 6\times 6$ for the $\sqrt{5}\times\sqrt{5}\times 1$ supercell
structures, and $6\times 6\times 3$ for the $\sqrt{5}\times\sqrt{5}\times 2$
supercell structures. The last Monkhorst-Pack grid is shifted by half grid
point along $k_{z}$ direction to sample equivalent grid points. The Bloch
states are then wannierizedMostofi et al. (2014); Marzari et al. (2012) with
the modified Becke-Johnson (mBJ) potentialBecke and Johnson (2006); Tran and
Blaha (2009) in $\Gamma$-centered k-points meshes of $6\times 6\times 4$ for
pristine DyTe2 and $4\times 4\times 4$ for the $\sqrt{5}\times\sqrt{5}\times
1$ supercell structures. The effect of spin-orbit coupling is found to be
negligible and not included here. The first Brillouin zone is sampled in
$100\times 100\times 50$ grid points with the Wannier Hamiltonian and the
tetrahedron method is employed to obtain the Fermi surface and nesting
functions.

####  .0.2 Nesting function

Figure S5: Electronic structure of pristine DyTe2. (a) Band structure along
symmetry directions calculated with the PBEsol functional where the electron
pocket derived from Dy-$d$ orbital appears at $\Gamma$. (b) Band structure
calculated with Wannier Hamiltonian that is obtained with mBJ potential. The
Dy-$d$ electron pocket at $\Gamma$ is lifted above the Fermi level. Right
panel shows the density of states in units of (number of electrons per eV
$\cdot$ 2Dy). Horizontal dashed line indicates electron doping level of
$+0.2e$ per 2Dy. (c) Fermi surface and (d - f) three in-plane nesting
functions, $N^{z}_{0}$, $N^{z}_{P}$, and $N^{z}_{AP}$ at $E=0$ eV. (g) Fermi
surface and (h - j) three in-plane nesting functions, $N^{z}_{0}$,
$N^{z}_{P}$, and $N^{z}_{AP}$ at $+0.2e$ doping level. The Tes-$p$ and Dy-$d$
orbital characters in the Fermi surfaces in (c) and (g) are illustrated with
blue and orange colors, respectively.

The Fermi surface nesting is often determined by the susceptibility
$\chi(\mathbf{q})$.Shim et al. (2004) However, in this work we calculate the
nesting function, $N(\mathbf{q})$, simply from the Fermi surface using the
tetrahedron method,

$\displaystyle N_{\gamma}(\mathbf{q})$ $\displaystyle=$
$\displaystyle\sum_{i,j}{s_{i}s_{j}w_{\gamma}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{j})\delta(\mathbf{k}_{j}-\mathbf{k}_{i}-\mathbf{q})},$

where the $i$ and $j$ indices run over the triangular facets of Fermi surface,
$s_{i}$ is the area of $i$-th facet whose center is the $\mathbf{k}_{i}$
wavevector, and ${\gamma}$ denotes the type of weighting function
$w_{\gamma}$, of the normalized group velocities
$\hat{\mathbf{v}}_{i}=\mathbf{v}_{i}/|\mathbf{v}_{i}|$ of the $i$-th facets.
We consider three types of weighting functions,

$\displaystyle w_{0}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{j})$
$\displaystyle=$ $\displaystyle 1,$ $\displaystyle
w_{P}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{j})$ $\displaystyle=$
$\displaystyle\frac{1}{2}(1+\hat{\mathbf{v}}_{i}\cdot\hat{\mathbf{v}}_{j}),$
$\displaystyle w_{AP}(\hat{\mathbf{v}}_{i},\hat{\mathbf{v}}_{j})$
$\displaystyle=$
$\displaystyle\frac{1}{2}(1-\hat{\mathbf{v}}_{i}\cdot\hat{\mathbf{v}}_{j}),$

where mutual parallel and anti-parallel group velocities at facets $i$ and $j$
are distinguished with $w_{P}$ and $w_{AP}$ weighting functions, respectively.
It is noteworthy that the parallel contribution is robust under shifts of the
chemical potential whereas the anti-parallel contribution is relevant to the
band gap opening. The in-plane nesting function is obtained by integrating
along the $k_{z}$ direction as,

$\displaystyle N^{z}_{\gamma}(q_{x},q_{y})$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}{dq_{z}N_{\gamma}(\mathbf{q})}.$

The $\mathbf{q}$ vector corresponding to the $\sqrt{5}\times\sqrt{5}$
modulation of DyTe2 appears in the vicinity of mutual peaks of $N^{z}_{P}$ and
$N^{z}_{AP}$ that explains the gap opening after the distortion.

####  .0.3 Formation energy

Figure S6: Calculated structural stability. (a - e) Te-bonding networks in the
$\sqrt{5}\times\sqrt{5}$ supercell structures with various Te vacancy
configurations. (f) Total energies of the various configurations shown in (a -
e) as a function of epitaxial strain. All energies are relative to the lowest
energy for each stoichiometry. (g) Out-of-plane, $c$ lattice constant under
the epitaxial strain. (h) Calculated formation energy, referenced to Dy2Te3
and Dy2Te5 structures. (i) Formation energy as a function of epitaxial strain.
Closed symbols in all panels denote the values obtained from fully relaxed
structure without the epitaxial constraint ($a=b$). The two closed triangle
symbols (A-B and A-C) breaks tetragonal symmetry ($a\neq b$) and are marked at
$a^{\prime}=\sqrt{ab}$.

For the $\sqrt{5}\times\sqrt{5}\times 1$ supercell structures we consider Te
mono- and di-vacancies per supercell, corresponding to $\delta=0.1$ and $0.2$,
respectively. Without Te vacancies, the structural distortion is negligible
indicating no phonon instability at the corresponding $\mathbf{q}$ vectors in
the pristine DyTe2. Introducing a single Te vacancy ($\delta=0.1$) causes
rearrangement of the Tes ions and the formation of chain-like bonding
networks. Two Te vacancies ($\delta=0.2$) can be arranged in three different
configurations, referred to as A-B, A-C, and A-D shown in Fig. S6. The
calculated relative total energies versus in-plane lattice constant for
$\delta=0.0$ (pristine), $0.1$ (mono-vacancy), and $0.2$ (di-vacancy) are
shown in Fig. S6 (f). The variation of lattice parameter $c$ versus in-plane
lattice constant for the different Te-deficiency $\delta$ values are shown in
Fig. S6 (g). The closed symbols denote the ground state energy and the
equilibrium lattice constants for all cases obtained by structural relaxation
without epitaxial constraint. We find that the A-C di-vacancy configuration is
the most stable. We note that A-B and A-C configurations break tetragonal
symmetry and the total energies (closed triangle symbols) deviate from the
epitaxial strain curves. In order to determine the structural stability
between different stoichiometric phases, the formation energy is calculated
relative to the Dy2Te3 and Dy2Te5 phases,

$\displaystyle
E_{\textrm{Form}}^{(x,y)}=\frac{E_{\textrm{T}}^{(x,y)}-\frac{5x-2y}{16}E_{\textrm{T}}^{(8,12)}-\frac{-3x+2y}{8}E_{\textrm{T}}^{(4,10)}}{x+y}.$

Here, $E_{\textrm{T}}^{(x,y)}$ is the total energy of a structure consisting
of $x$ Dy and $y$ Te atoms. The calculated formation energy [Fig. S6 (h)]
predicts that the $\delta=0.2$ A-C configuration is the most favored structure
while DyTe2 and Dy4Te7 phases are found to be less favorable states.
Considering the distortion from pristine DyTe2 to Te-deficient states, it is
noteworthy that the formation energies of DyTe2, DyTe1.9, and
DyTe${}_{1.8,\textrm{A-C}}$ align on a line that indicates mixed occurrence
between them.

It is thus an interesting question whether epitaxial strain which varies with
film thickness can affect their formation energies due to the different
equilibrium in-plane lattice constants,
$a_{\delta=0.2}<a_{\delta=0.1}<a_{\delta=0.0}$. Figure S6 (i) shows the
formation energy as a function of epitaxial strain for the five structures
that predicts a structural transition from $\delta=0.2$ A-C to $\delta=0.1$
phase with increasing epitaxial strain. The calculated threshold lattice
constant is, however, larger than the maximum in-plane lattice constant in the
experiment.

Table 1: Calculated total energies for the $\sqrt{5}\times\sqrt{5}\times 2$ supercell structures with A-C divacancy configuration. One of the two Te square net layers is laterally shifted by $(d{{}_{1}},d{{}_{2}})=d{{}_{1}}\mathbf{a}^{{}^{\prime}}_{1}+d{{}_{2}}\mathbf{a}^{{}^{\prime}}_{2}$, where $\mathbf{a}^{{}^{\prime}}_{i}$ are the in-plane translation lattice vectors of the supercell. For example, the vacancy sites A and C in Fig. S6(a) correspond to $(0.0,0.0)$ and $(0.2,0.4)$, respectively. Lateral shift |  $\Delta E$ (meV/Dy)
---|---
(0.0, 0.0) | 0.0
(0.3, 0.1) | -6.0
(0.6, 0.2) | -8.4
(0.9, 0.3) | -5.8
(0.2, 0.4) | -7.3
(0.5, 0.5) | -10.8


# On the helicity conservation for Mie optical cavities

Jorge Olmos-Trigo<EMAIL_ADDRESS>Donostia International Physics Center
(DIPC), 20018 Donostia-San Sebastián, Spain  Xavier Zambrana-Puyalto
<EMAIL_ADDRESS>Istituto Italiano di Tecnologia, Via Morego 30, 16163
Genova, Italy.

###### Abstract

The study of helicity in the context of light-matter interactions is an
increasing area of research. However, some fundamental aspects of the helicity
content of light fields inside scatterers have been overlooked. In this work,
we demonstrate that the helicity of light fields inside lossless spherical
cavities cannot be either conserved or sign-flipped as a result of a
scattering event. The underlying reason is that the internal electric and
magnetic Mie coefficients cannot simultaneously oscillate with equal amplitude
and equal/opposite phase. Our analytical demonstration is fulfilled regardless
of the refractive index, sphere size, and multipolar order. In addition, we
show that the helicity of light fields inside lossy spheres can be conserved.
This fact is in striking contrast to the behavior of the scattered field,
whose helicity content cannot be conserved precisely when the sphere has
losses. Finally, we show that the helicity content of internal fields can be
flipped for materials with gain.

The study of light-matter interactions is ubiquitous across physical sciences.
In particular, the increasing reach of Photonic technologies Saleh and Teich
(1991) is due to the advances in controlling light-matter interactions.

Many different properties of light-matter interactions can be partially
controlled or tweaked. In recent years, a property of light-matter
interactions that had been overlooked for a long time has gained significant
attention: electromagnetic helicity Calkin (1965). The first fundamental
explanation of the role of electromagnetic helicity in Maxwell equations dates
back to 1965 when Calkin showed that the electromagnetic fields in vacuum
satisfy a continuous symmetry: the electromagnetic duality. All continuous
symmetries have a generator, and Calkin found helicity as the generator of
duality for the electromagnetic fields in vacuum. Now, helicity is defined as
the projection of the total angular momentum ($\bf{J}$) onto the linear
momentum of the wave ($\bf{p}$), namely,
$\Lambda=({\bf{J}\cdot\bf{p}})/{|\bf{p}|}$ Fernandez-Corbaton  _et al._
(2012). Its definition becomes simpler in the plane wave decomposition of an
electromagnetic field, as helicity is associated with the handedness of
circular polarization of each plane wave to its momentum vector. However, due
to the non-existence of magnetic monopoles in nature, which prevents Maxwell
equations from being dual-symmetric in the presence of matter, the interest of
the optical community in studying duality symmetry and helicity was residual
for a large number of decades with very few exceptions Zwanziger (1968);
Białynicki-Birula (1994); Bialynicki-Birula  _et al._ (1981).

The study of helicity upon light-matter interactions took a clear boost in
2013 when it was theoretically unveiled that helicity can be conserved upon
scattering for dual materials, regardless of their geometry Fernandez-Corbaton
_et al._ (2013). A dual material is such that the ratio between its relative
magnetic permeability and its electric permittivity happens to be constant
($\mu/\varepsilon=\text{constant}$). These materials restore duality symmetry
in the macroscopic approximation of Maxwell equations. Unfortunately, there
are no magnetic materials ($\mu\neq 1$) at optical frequencies. The
implication is that it should be impossible to experimentally observe helicity
conservation at optical frequencies. However, shortly after, it was
demonstrated Zambrana-Puyalto  _et al._ (2013a) that high refractive index
particles could approximately conserve helicity due to their electric and
magnetic resonances García-Etxarri  _et al._ (2011); Kuznetsov  _et al._
(2016). Since then, a lot of work has been done to characterize the role of
helicity in light-matter interactions Tischler  _et al._ (2014); Fernandez-
Corbaton  _et al._ (2016); Alpeggiani  _et al._ (2018); Forbes and Jones
(2021), as well as to study its relation with chirality Cameron  _et al._
(2017); Gutsche and Nieto-Vesperinas (2018); Poulikakos  _et al._ (2019);
Hanifeh  _et al._ (2020).

One of the first platforms that was used to study the role of helicity in
light-matter interactions is Mie Theory Mie (1908). Mie Theory studies the
interaction of a plane wave with a spherical scatterer. It is a widely used
platform to discover new analytical effects, as well as to predict
experimental measurements Gouesbet  _et al._ (2011). In particular, Mie Theory
has been used to understand the role of helicity in a variety of scattering
phenomena at the nanoscale, such as the so-called Kerker conditions Nieto-
Vesperinas  _et al._ (2011); Olmos-Trigo  _et al._ (2020a), enhanced optical
localization errors provided by optical mirages Olmos-Trigo  _et al._ (2019a),
or non-radiating sources Labate  _et al._ (2017) such as hybrid optical
anapoles Sanz-Fernández  _et al._ (2021). Moreover, helicity-dependent optical
forces Liu  _et al._ (2018); Nieto-Vesperinas and Xu (2021) have given rise to
a large plethora of interesting phenomena with applications in the biomedical
and pharmaceutical industries, ranging from chiral sensing Graf  _et al._
(2019); Lasa-Alonso  _et al._ (2020, 2022) chiral sorting and enantioselective
and enantiospecificdetection Hayat  _et al._ (2015); Toyoda  _et al._ (2013);
Nieto-Vesperinas (2015), to optical tweezers Tkachenko and Brasselet (2014);
Ali  _et al._ (2020), among others.

Most of the works dealing with helicity have studied scattering features, yet
in the very last few years, some groups have started to study the helicity of
fields trapped in a cavity Feis  _et al._ (2020); Scott  _et al._ (2020);
Voronin  _et al._ . Having the enhancement of chiral sensing in mind, the
authors of Feis  _et al._ (2020) came up with the design of a cavity that
tightly concentrates light in areas where the helicity of the field is
maintained. However, it is important to note that this intriguing and
significant effect relies on numerical methods under very specific
illumination conditions. Here, we show that a simple lossy sphere can also be
used to create a cavity that internally conserves the helicity content of the
incident beam. In addition, we analytically demonstrate that the internal
field of a lossless sphere cannot conserve the helicity content of an incoming
beam. Note that our findings of the conservation of helicity for the internal
field of a spherical cavity are in striking contrast with what is known about
the scattered field, i.e. the scattered field can only conserve the helicity
of the incoming beam when the spherical cavity is lossless Olmos-Trigo  _et
al._ (2020b). Last but not least, we also show that the helicity of the
internal field can be flipped if a sphere with optical gain is used.

Next, we lay out the framework that we will use to show how helicity works for
internal fields. First, let us consider an incoming field with well-defined
helicity $\sigma=\pm 1$ Olmos-Trigo  _et al._ (2020c, 2019b),

$\displaystyle{\bf E}^{\sigma}_{\text{inc}}(k{\bf r})$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{\infty}\sum_{m=-l}^{+l}C_{lm}^{\sigma}\bm{\Psi}_{lm}^{\sigma}(k{\bf
r}).$ (1)

Here $C_{lm}^{\sigma}$ denotes the incoming coefficients characterizing the
nature of the wave, $k$ is the radiation wavevector, and

$\displaystyle\bm{\Psi}_{lm}^{\sigma}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left[{\bm{N}}_{lm}+\sigma{\bm{M}}_{lm}\right],$
(2) $\displaystyle{\bm{M}}_{lm}$ $\displaystyle\equiv$ $\displaystyle
j_{l}(kr)\bm{X}_{lm},\quad{\bm{N}}_{lm}\equiv\frac{1}{k}\bm{\nabla}\times{\bm{M}}_{lm},$
(3) $\displaystyle\bm{X}_{lm}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{l(l+1)}}{\bf{L}}Y_{lm}(\theta,\varphi).$ (4)

Here $\bm{M}_{lm}$ and $\bm{N}_{lm}$ are the so-called Hansen’s multipoles
Jackson (1999), $j_{l}(kr)$ are the spherical Bessel functions,
$Y_{lm}(\theta,\varphi)$ are the spherical harmonics, $\theta$ and $\varphi$
being the polar and azimuthal angles, and ${\bf{L}}=\left\\{-{\rm i}{\bf
r}\times\bm{\nabla}\right\\}$ is the total angular momentum operator. Let us
recall that the multipoles $\bm{\Psi}_{lm}^{\sigma}$ are simultaneous
eigenvectors of the squared angular momentum $L^{2}$, the projection of the
angular momentum on one direction $L_{z}$, and the helicity operator
${\Lambda}$ Fernandez-Corbaton  _et al._ (2013), with eigenvalues $l(l+1)$,
$m$ and $\sigma$, respectively Jackson (1999).

Figure 1: Decomposition of the internal and scattered fields in the Riemann-
Silberstein representation, $\Lambda^{\sigma^{\prime}}=|{\bf E}+{\rm
i}\sigma^{\prime}Z{\bf H}|$, under illumination of a well-defined helicity
Gaussian beam ($\sigma=+1$). According to the colormap, black color represents
the absence of the $\Lambda^{\sigma^{\prime}}$ component. The white circle
represents a lossless sphere with radius $a$ and refractive index contrast
$m\in\mathbb{R}$. a) Emergence of the First Kerker condition provided by a
Germanium (Ge) sphere of $m=\sqrt{\epsilon}=4.2$ and $a=140$ nm under
illumination of a Gaussian beam focused with $\rm{NA}=0.15$ at $\lambda=1347$
nm. In this case, $|\Lambda_{\rm{sca}}^{-}|=0$ while
$|\Lambda_{\rm{int}}^{-}|\neq 0$. The plot has a dimension of $1.29\times
1.29\mu m^{2}$. b) Electromagnetic duality restoration provided by a sphere of
$\epsilon=\mu=4.2$ and $a=140$ nm under illumination of a Gaussian beam
focused with $\rm{NA}=0.15$ at $\lambda=1347$ nm. In this scenario,
$|\Lambda_{\rm{sca}}^{-}|=|\Lambda_{\rm{int}}^{-}|=0$. The plot has a
dimension of $1.29\times 1.29\mu m^{2}$. c) Emergence of the hybrid anapole
mediated by a sphere of $m=\sqrt{\epsilon}=3.33$ and $a=453$ nm under the
illumination of a Gaussian beam with $\rm{NA}=0.9$ at $\lambda=633$ nm. In
this scenario, $|\Lambda_{\rm{sca}}^{+}|=|\Lambda_{\rm{sca}}^{-}|=0$ while
$|\Lambda_{\rm{int}}^{-}|\neq 0$ and $|\Lambda_{\rm{int}}^{+}|\neq 0$. The
plot has a dimension of $2.01\times 2.01\mu m^{2}$.

At this point, let us expand the electric field inside a sphere (internal
electric field) in the same $\bm{\Psi}_{lm}^{\sigma}$ basis. The internal
field is obtained as the result of a scattering process involving the incoming
field ${\bf E}^{\sigma}_{\text{inc}}(k{\bf r})$ and a sphere of radius $a$:

$\displaystyle{\bf E}^{\sigma}_{\rm{int}}(k_{i}{\bf r})$ $\displaystyle=$
$\displaystyle\sum_{\sigma^{\prime}=\pm 1}{\bf
F}^{\sigma\sigma^{\prime}}_{\rm{int}}(k_{i}{\bf r}),$ (5) $\displaystyle{\bf
F}^{\sigma\sigma^{\prime}}_{\rm{int}}(k_{i}{\bf r})$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{\infty}\sum_{m=-l}^{+l}F^{\sigma\sigma^{\prime}}_{lm}\bm{\Psi}_{lm}^{\sigma^{\prime}}(k_{i}{\bf
r}),$ (6)

where
$F^{\sigma\sigma^{\prime}}_{lm}=C^{\sigma}_{lm}\left({d_{l}+\sigma\sigma^{\prime}c_{l}}\right)/{2}$.
Here $d_{l}$ and $c_{l}$ denote the internal electric and magnetic Mie
coefficients, respectively Bohren and Huffman (2008) and $k_{i}=mk$, with
$m=m_{p}/m_{h}$, $m_{p}$ and $m_{h}$ being the refractive index of the sphere
and medium, respectively.

From Eq. (6), we can notice that $d_{l}=c_{l}$ preserves the incoming helicity
($F^{\sigma\sigma^{\prime}}_{lm}=C^{\sigma}_{lm}d_{l}\delta_{\sigma\sigma^{\prime}}$).
Similarly, it is straightforward to infer that $d_{l}=-c_{l}$ flips the
incoming helicity
($F^{\sigma\sigma^{\prime}}_{lm}=C^{\sigma}_{lm}d_{l}\left[1-\delta_{\sigma\sigma^{\prime}}\right]$).
At this stage, let us draw our attention to the internal electric and magnetic
Mie coefficients expressed in phase-shifts notation Luk’yanchuk  _et al._
(2017),

$\displaystyle d_{l}=-\frac{{\rm i}m}{F^{(a)}_{l}+{\rm i}G^{(a)}_{l}},$
$\displaystyle c_{l}=\frac{{\rm i}m}{F^{(b)}_{l}+{\rm i}G^{(b)}_{l}},$ (7)

where

$\displaystyle F^{(a)}_{l}$ $\displaystyle=$ $\displaystyle
m\psi^{{}^{\prime}}_{l}(q)\psi_{l}(mq)-\psi_{l}(q)\psi^{{}^{\prime}}_{l}(mq),$
(8) $\displaystyle G^{(a)}_{l}$ $\displaystyle=$ $\displaystyle
mx^{{}^{\prime}}_{l}(q)\psi_{l}(mq)-\psi^{{}^{\prime}}_{l}(mq)x_{l}(q),$ (9)
$\displaystyle F^{(b)}_{l}$ $\displaystyle=$ $\displaystyle
m\psi^{{}^{\prime}}_{l}(mq)\psi_{l}(q)-\psi_{l}(mq)\psi^{{}^{\prime}}_{l}(q),$
(10) $\displaystyle G^{(b)}_{l}$ $\displaystyle=$ $\displaystyle
mx_{l}(q)\psi^{{}^{\prime}}_{l}(mq)-\psi_{l}(mq)x^{{}^{\prime}}_{l}(q).$ (11)

Here $\psi_{l}(q)=(\frac{\pi q}{2})^{1/2}J_{l+\frac{1}{2}}(q)$ and
$x_{l}(q)=(\frac{\pi q}{2})^{1/2}N_{l+\frac{1}{2}}(q)$ denote the Riccati-
Bessel functions, where $J_{l}(q)$ and $N_{l}(q)$ are the Bessel functions of
first and second kind, respectively Bohren and Huffman (2008); ′ denotes the
derivative with respect the argument, $q=2\pi a/\lambda$ is the size
parameter, and $\lambda$ is the radiation wavelength. Now, we impose
$d_{l}=\tau c_{l}$, where $\tau=\pm 1$. After doing some algebra we arrive
from Eq. (7) to

$\displaystyle F^{(a)}_{l}+{\rm i}G^{(a)}_{l}=-\tau\left(F^{(b)}_{l}+{\rm
i}G^{(b)}_{l}\right),$ (12)

with

$\displaystyle F^{(a)}_{l}=mA-B,$ $\displaystyle G^{(a)}_{l}=mC-D,$ (13)
$\displaystyle F^{(b)}_{l}=mB-A,$ $\displaystyle G^{(a)}_{l}=mD-C.$ (14)

Here, we have defined $A=\psi^{\prime}_{l}(q)\psi_{l}(mq)$,
$B=\psi_{l}(q)\psi^{\prime}_{l}(mq)$, $C=x^{\prime}_{l}(q)\psi_{l}(mq)$, and
$D=\psi^{\prime}_{l}(mq)x_{l}(q)$. Now, by taking into account Eqs.(13)-(14),
it can be shown that $d_{l}=\tau c_{l}$, given by Eq. (12), can be re-written
as

$\displaystyle A+\tau B=-{\rm i}\left(C+\tau D\right).$ (15)

Now, by inspecting Eq. (15), it is straightforward to re-write $d_{l}=\tau
c_{l}$ as

$\tau\psi^{\prime}_{l}(mq)[\psi_{l}(q)+{\rm
i}x_{l}(q)]+\psi_{l}(mq)[\psi^{\prime}_{l}(q)+{\rm
i}x^{{}^{\prime}}_{l}(q)]=0.$ (16)

At this stage and by identifying $h_{l}(q)=\psi_{l}(q)+{\rm i}x_{l}(q)$ as the
spherical Hankel function of second kind Bohren and Huffman (2008), we finally
get

$\psi_{l}(mq)h^{\prime}_{l}(q)+\tau\psi^{\prime}_{l}(mq)h_{l}(q)=0.$ (17)

Eq. (17) represents a notable simplification as it allows us to compute
$d_{l}=\tau c_{l}$ by making use of fundamental properties of just two
spherical Bessel functions. At this point, let us split the solutions of Eq.
(17) into two possible physical scenarios, namely, non-absorbing
($m\in\mathbb{R}$) and lossy materials ($m\in\mathbb{C}$ with
$|\Im\\{m\\}|\neq 0$). Let us first draw our attention to the lossless case
($m\in\mathbb{R}$), in which the complex Eq. (17) can be re-expressed as two
real-valued equations:

$\displaystyle F^{(a)}_{l}=-\tau
F^{(b)}_{l}\rightarrow\psi_{l}(mq)\psi^{{}^{\prime}}_{l}(q)+\tau\psi^{{}^{\prime}}_{l}(mq)\psi_{l}(q)=0,$
(18) $\displaystyle G^{(a)}_{l}=-\tau
G^{(b)}_{l}\rightarrow\psi_{l}(mq)x^{{}^{\prime}}_{l}(q)+\tau\psi^{{}^{\prime}}_{l}(mq)x_{l}(q)=0.$
(19)

Now, it is important to notice that Eqs. (18)-(19) need to be simultaneously
fulfilled to obtain $d_{l}=\tau c_{l}$. At this point, let us drive our
attention to Lemma 1 Watson (1995):

  1. 1.

The zeros of any cylinder function or its derivative are simple; If $\nu$ is
real, then $J_{\nu}(z)$, $J^{\prime}_{\nu}(z)$, $N_{\nu}(z)$,
$N^{\prime}_{\nu}(z)$ each have an infinite number of positive real zeros. The
$m$-th positive zeros of these spherical functions are denoted by $j_{\nu,m}$,
$j^{\prime}_{\nu,m}$, $n_{\nu,m}$, $n^{\prime}_{\nu,m}$, and interlace
according to
$\nu<j^{\prime}_{\nu,1}<n_{\nu,1}<n^{\prime}_{\nu,1}<j_{\nu,1}<j^{\prime}_{\nu,2}<n_{\nu,2}<...\;.$

Let us now use the Lemma 1 to prove that the zeros of the Riccati-Bessel
functions are also interlaced: Let be $\psi_{\nu}(z)=(\frac{\pi
z}{2})^{1/2}J_{\nu+\frac{1}{2}}(z)=0$. For $z\neq 0$, this solution is given
by the zeros of $J_{\nu+\frac{1}{2}}(z)$. Let us now take the first derivative
of the Riccati-Bessel to then set it to zero, namely,

$\psi^{\prime}_{\nu}(z)=\frac{1}{2}{J_{\nu+\frac{1}{2}}(z)}+zJ^{\prime}_{\nu+\frac{1}{2}}(z)=0.$
(20)

From Eq. (20) and by taking into account that we have first settled
$J_{\nu+\frac{1}{2}}(z)=0$, the only possible solution of
$\psi^{\prime}_{\nu}(z)=0$ might be given by
$J^{\prime}_{\nu+\frac{1}{2}}(z)=0$. However, the latter cannot be met due to
the interlaced property of the zeros of Bessel functions provided by Lemma 1.
Consequently, the possible solutions of Eqs. (18)-(19) are notably reduced
(see 111Notice that if $\psi_{l}^{\prime}(z)=0$ then $\psi_{l}(z)\neq 0$. On
other hand, if $\psi_{l}(z)=0$ then $\psi^{\prime}_{l}(z)\neq 0$. Notice that
if $x_{l}^{\prime}(z)=0$ then $x_{l}(z)\neq 0$. On other hand, if $x_{l}(z)=0$
then $x^{\prime}_{l}(z)\neq 0$. to get insight into all possible
combinations). That is, the only possibilities that make $d_{l}=\tau c_{l}$
when $m\in\mathbb{R}$ are:

  1. i.

First Kerker condition (helicity conservation in scattering) Nieto-Vesperinas
_et al._ (2011); Zambrana-Puyalto  _et al._ (2013b), mathematically expressed
as $\psi_{l}^{\prime}(mq)=0$ (node of the first kind) and $\psi_{l}(mq)=0$
(node of the second kind) Olmos-Trigo  _et al._ (2020c).

  2. ii.

Hybrid or Kerker anapoles (optical transparency) Luk’yanchuk  _et al._ (2017),
mathematically given by $\psi_{l}^{\prime}(mq)=0$ and $\psi_{l}^{\prime}(q)=0$
or $\psi_{l}(mq)=0$ and $\psi_{l}(q)=0$ Sanz-Fernández  _et al._ (2021).

First, the conservation of helicity in scattering for a lossless sphere
precludes the conservation of the internal helicity (see Fig 1a)).
Mathematically, we can demonstrate this by using the property that the zeros
of the Riccati-Bessel functions are interlaced. That is, it is impossible to
fulfill both Eqs (18)-(19) if $\psi_{l}^{\prime}(mq)=0$ or $\psi_{l}(mq)=0$.
On physics grounds, this is a result of the fact that we are dealing with
spheres with $\mu=1$ and $\epsilon^{2}=m\in\mathbb{R}$, with $\epsilon$ and
$\mu$ denoting the electric permittivity and magnetic permeability,
respectively. Now, it is known that electromagnetic duality can only be
restored when $\epsilon=\mu$ Fernandez-Corbaton  _et al._ (2013). When duality
is restored, helicity is conserved for both the internal and scattered field,
as it can be observed in Fig. 1b). However, for non-magnetic lossless
particles, this physical picture is precluded: if the scattered helicity is
preserved, then the internal helicity is not.

Second, hybrid anapoles, namely spectral points in which the scattering
associated with a given order $l$ vanishes, also prevent $d_{l}=\tau c_{l}$.
Mathematically, the hybrid anapole condition imposes that
$F^{(a)}_{l}=F^{(b)}_{l}=0$, which yields $c_{l}/d_{l}=m$ when
$\psi_{l}^{\prime}(mq)=0$ and $\psi_{l}^{\prime}(q)=0$, and $d_{l}/c_{l}=m$
when $\psi_{l}(mq)=0$ and $\psi_{l}(q)=0$ Sanz-Fernández  _et al._ (2021).
Physically, we can notice that hybrid anapoles are a particular solution of
the first Kerker condition with zero scattering, and as explained above: a
sphere with $\mu=1$ cannot be dual. As a result, we conclude that hybrid
anapoles of order $l$ not only inhibit the conservation of the internal
helicity, as can be inferred from Fig 1c), but also constrain the internal Mie
coefficients to have a very specific relation given by $c_{l}/d_{l}=m^{\pm
1}$. Remarkably, this phenomenon remains valid regardless of the refractive
index, optical size, making our proof general for anapoles associated with
higher multipolar orders.

Figure 2: Decomposition of the internal and scattered fields in the Riemann-
Silberstein representation, $\Lambda^{\sigma^{\prime}}=|{\bf E}+{\rm
i}\sigma^{\prime}Z{\bf H}|$, under illumination of a well-defined helicity
Gaussian beam ($\sigma=+1$). According to the colormap, black color represents
the absence of the $\Lambda^{\sigma^{\prime}}$ component. The white circle
represents the lossy sphere with radius $a$ and refractive index contrast
$m\in\mathbb{C}$. a) Conservation of the helicity content of light fields
inside a sphere of $m=4.2+0.54{\rm i}$ and $a=140$ nm under illumination of a
Gaussian beam with $\rm{NA}=0.15$ at $\lambda=407$ nm. In this scenario, the
scattered helicity is not preserved while $|\Lambda_{\rm{int}}^{-}|=0$. The
plot has a dimension of $1.04\times 1.04\mu m^{2}$. b) Flipping of the
helicity content of light fields provided by a sphere of
$m=\sqrt{\epsilon}=4.2-0.37{\rm i}$ and $a=140$ nm under illumination of a
Gaussian beam with $\rm{NA}=0.15$ at $\lambda=1038$ nm. In this case,
$|\Lambda_{\rm{int}}^{+}|=0$. The plot has a dimension of $1.32\times 1.32\mu
m^{2}$. c) Second Kerker condition provided by a sphere of
$m=\sqrt{\epsilon}=4.2-0.33{\rm i}$ and $a=140$ nm under illumination of a
Gaussian beam with $\rm{NA}=0.15$ at $\lambda=1038$ nm. In this case,
$|\Lambda_{\rm{sca}}^{+}|=0$. The plot has a dimension of $1.07\times 1.07\mu
m^{2}$.

Hitherto, we have found that the conservation of helicity in scattering and
the emergence of hybrid anapoles are related to the fact that the helicity of
light fields inside a lossless sphere cannot be either conserved or sign-
flipped as a result of a scattering event; namely, $d_{l}\neq\tau c_{l}$ for
$m\in\mathbb{R}$. In this vein, it has been analytically proven in Ref. Olmos-
Trigo  _et al._ (2020b) and Ref. Sanz-Fernández  _et al._ (2021) that the
conservation of helicity in scattering and the emergence of hybrid anapoles
can only occur for lossless particles. Next, we show that conservation of
helicity for the internal fields can only happen for the opposite case, i.e.
for spheres made of a lossy material where $m\in\mathbb{C}$ with
$|\Im\\{m\\}|\neq 0$. To show that, let us consider $d_{l}=\tau c_{l}$ without
assuming constraints on the real and imaginary parts in Eq. (17). Taking this
crucial fact into account, we can re-write $d_{l}=\tau c_{l}$ as

$-\tau=\frac{\psi_{l}(mq)h^{\prime}_{l}(q)}{\psi^{\prime}_{l}(mq)h_{l}(q)}.$
(21)

Eq. (21) is a transcendental equation that can only be satisfied for
$|\Im\\{m\\}|\neq 0$. As a matter of fact, it can be shown that for $\tau=+1$
the set of solutions can only be given for $\Im\\{m\\}>0$ while for $\tau=-1$
these can only be met for $\Im\\{m\\}<0$. This phenomenon is depicted in Fig.
2. In particular, in Fig. 2a), we show a Ge-like nanosphere in the visible
spectral range ($\lambda=407$ nm, $m=4.2+0.54{\rm i}$) whose internal field
conserves the helicity content of the incident beam. In Fig. 2 we also show
two spheres in the presence of optical gain. On the one hand, in Fig. 2b) we
show the case of a sphere where the helicity content of the internal field is
just the opposite with respect to the incoming beam. On the other hand, in
Fig. 2c) we show a sphere whose scattered field helicity content is just the
opposite of the incident beam.

In conclusion, in our work we have unraveled a fundamental property of the
internal Mie coefficients: we have demonstrated that the helicity of light
fields inside lossless spheres cannot be either conserved or sign-flipped as a
result of a scattering event. This demonstration is analytical and involves
special features of the scattered field, such as the conservation of helicity
and the emergence of hybrid anapoles. Our proof does not depend on the
incident polarization, optical size, and multipolar order. We have also shown
that, in striking contrast to the behavior of the scattered field, losses are
a compulsory requirement to conserve the helicity content of the internal
field. That is, $d_{l}=c_{l}$ can only happen for lossy materials. Finally, we
have shown that the helicity content of internal fields can be flipped for
materials with gain. Note that the main findings of helicity conservation
inside spherical cavities can be extrapolated to other geometries such as
disks or cylinders, finding potential applications in chiral sensing and
chiral spectroscopy techniques.

J. O. T. acknowledges support from the Spanish Ministerio de Ciencia e
Innovacion (PID2019- 109905GA-C2) and from Eusko Jaurlaritza (KK-2021/00082).
Moreover, J. O. T. acknowledges funding from the Basque Government’s IKUR
initiative on Quantum technologies (Department of Education).

## References

  * Saleh and Teich (1991) B. E. Saleh and M. C. Teich,  _Fundamentals of photonics_ , Vol. 22 (Wiley New York, 1991).
  * Calkin (1965) M. Calkin, Am. J. Phys 33, 958 (1965).
  * Fernandez-Corbaton  _et al._ (2012) I. Fernandez-Corbaton, X. Zambrana-Puyalto, and G. Molina-Terriza, Phys. Rev. A 86, 042103 (2012).
  * Zwanziger (1968) D. Zwanziger, Physical Review 176, 1489 (1968).
  * Białynicki-Birula (1994) I. Białynicki-Birula, Acta Physica Polonica A 1, 97 (1994).
  * Bialynicki-Birula  _et al._ (1981) I. Bialynicki-Birula, E. Newman, J. Porter, J. Winicour, B. Lukacs, Z. Perjes, and A. Sebestyen, Journal of Mathematical Physics 22, 2530 (1981).
  * Fernandez-Corbaton  _et al._ (2013) I. Fernandez-Corbaton, X. Zambrana-Puyalto, N. Tischler, X. Vidal, M. L. Juan, and G. Molina-Terriza, Phys. Rev. Lett. 111, 060401 (2013).
  * Zambrana-Puyalto  _et al._ (2013a) X. Zambrana-Puyalto, X. Vidal, M. L. Juan, and G. Molina-Terriza, Opt. Express 21, 17520 (2013a).
  * García-Etxarri  _et al._ (2011) A. García-Etxarri, R. Gómez-Medina, L. S. Froufe-Pérez, C. López, L. Chantada, F. Scheffold, J. Aizpurua, M. Nieto-Vesperinas, and J. J. Sáenz, Opt. Express 19, 4815 (2011).
  * Kuznetsov  _et al._ (2016) A. I. Kuznetsov, A. E. Miroshnichenko, M. L. Brongersma, Y. S. Kivshar, and B. Luk’yanchuk, Science 354, aag2472 (2016).
  * Tischler  _et al._ (2014) N. Tischler, I. Fernandez-Corbaton, X. Zambrana-Puyalto, A. Minovich, X. Vidal, M. L. Juan, and G. Molina-Terriza, Light: Science & Applications 3, e183 (2014).
  * Fernandez-Corbaton  _et al._ (2016) I. Fernandez-Corbaton, M. Fruhnert, and C. Rockstuhl, Physical Review X 6, 031013 (2016).
  * Alpeggiani  _et al._ (2018) F. Alpeggiani, K. Bliokh, F. Nori, and L. Kuipers, Physical review letters 120, 243605 (2018).
  * Forbes and Jones (2021) K. A. Forbes and G. A. Jones, Journal of Optics 23, 115401 (2021).
  * Cameron  _et al._ (2017) R. P. Cameron, J. B. Götte, S. M. Barnett, and A. M. Yao, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20150433 (2017).
  * Gutsche and Nieto-Vesperinas (2018) P. Gutsche and M. Nieto-Vesperinas, Scientific reports 8, 1 (2018).
  * Poulikakos  _et al._ (2019) L. V. Poulikakos, J. A. Dionne, and A. García-Etxarri, Symmetry 11, 1113 (2019).
  * Hanifeh  _et al._ (2020) M. Hanifeh, M. Albooyeh, and F. Capolino, ACS Photonics 7, 2682 (2020).
  * Mie (1908) G. Mie, Annalen der physik 330, 377 (1908).
  * Gouesbet  _et al._ (2011) G. Gouesbet, J. Lock, and G. Gréhan, Journal of Quantitative Spectroscopy and Radiative Transfer 112, 1 (2011).
  * Nieto-Vesperinas  _et al._ (2011) M. Nieto-Vesperinas, R. Gomez-Medina, and J. J. Saenz, J. Opt. Soc. Am. A 28, 54 (2011).
  * Olmos-Trigo  _et al._ (2020a) J. Olmos-Trigo, D. R. Abujetas, C. Sanz-Fernández, J. A. Sánchez-Gil, and J. J. Sáenz, Physical Review Research 2, 013225 (2020a).
  * Olmos-Trigo  _et al._ (2019a) J. Olmos-Trigo, C. Sanz-Fernández, A. García-Etxarri, G. Molina-Terriza, F. S. Bergeret, and J. J. Sáenz, Phys. Rev. A 99, 013852 (2019a).
  * Labate  _et al._ (2017) G. Labate, A. Alù, and L. Matekovits, Physical Review A 95, 063841 (2017).
  * Sanz-Fernández  _et al._ (2021) C. Sanz-Fernández, M. Molezuelas-Ferreras, J. Lasa-Alonso, N. de Sousa, X. Zambrana-Puyalto, and J. Olmos-Trigo, Laser & Photonics Reviews , 2100035 (2021).
  * Liu  _et al._ (2018) L. Liu, A. Di Donato, V. Ginis, S. Kheifets, A. Amirzhan, and F. Capasso, Physical review letters 120, 223901 (2018).
  * Nieto-Vesperinas and Xu (2021) M. Nieto-Vesperinas and X. Xu, Physical Review Research 3, 043080 (2021).
  * Graf  _et al._ (2019) F. Graf, J. Feis, X. Garcia-Santiago, M. Wegener, C. Rockstuhl, and I. Fernandez-Corbaton, ACS Photonics 6, 482 (2019).
  * Lasa-Alonso  _et al._ (2020) J. Lasa-Alonso, D. R. Abujetas, Á. Nodar, J. A. Dionne, J. J. Sáenz, G. Molina-Terriza, J. Aizpurua, and A. García-Etxarri, Acs Photonics 7, 2978 (2020).
  * Lasa-Alonso  _et al._ (2022) J. Lasa-Alonso, J. Olmos-Trigo, A. García-Etxarri, and G. Molina-Terriza, Mater. Adv. , (2022).
  * Hayat  _et al._ (2015) A. Hayat, J. B. Mueller, and F. Capasso, Proceedings of the National Academy of Sciences 112, 13190 (2015).
  * Toyoda  _et al._ (2013) K. Toyoda, F. Takahashi, S. Takizawa, Y. Tokizane, K. Miyamoto, R. Morita, and T. Omatsu, Physical review letters 110, 143603 (2013).
  * Nieto-Vesperinas (2015) M. Nieto-Vesperinas, Phys. Rev. A 92, 023813 (2015).
  * Tkachenko and Brasselet (2014) G. Tkachenko and E. Brasselet, Nature communications 5, 1 (2014).
  * Ali  _et al._ (2020) R. Ali, F. A. Pinheiro, R. S. Dutra, F. S. Rosa, and P. A. M. Neto, Nanoscale 12, 5031 (2020).
  * Feis  _et al._ (2020) J. Feis, D. Beutel, J. Köpfler, X. Garcia-Santiago, C. Rockstuhl, M. Wegener, and I. Fernandez-Corbaton, Physical Review Letters 124, 033201 (2020).
  * Scott  _et al._ (2020) P. Scott, X. Garcia-Santiago, D. Beutel, C. Rockstuhl, M. Wegener, and I. Fernandez-Corbaton, Applied Physics Reviews 7, 041413 (2020).
  * Voronin  _et al._ (0) K. Voronin, A. S. Taradin, M. V. Gorkunov, and D. G. Baranov, ACS Photonics 0, null (0), https://doi.org/10.1021/acsphotonics.2c00134 .
  * Olmos-Trigo  _et al._ (2020b) J. Olmos-Trigo, C. Sanz-Fernández, D. R. Abujetas, J. Lasa-Alonso, N. de Sousa, A. García-Etxarri, J. A. Sánchez-Gil, G. Molina-Terriza, and J. J. Sáenz, Physical Review Letters 125, 073205 (2020b).
  * Olmos-Trigo  _et al._ (2020c) J. Olmos-Trigo, D. R. Abujetas, C. Sanz-Fernández, X. Zambrana-Puyalto, N. de Sousa, J. A. Sánchez-Gil, and J. J. Sáenz, Physical Review Research 2, 043021 (2020c).
  * Olmos-Trigo  _et al._ (2019b) J. Olmos-Trigo, M. Meléndez, R. Delgado-Buscalioni, and J. J. Sáenz, Opt. Express 27, 16384 (2019b).
  * Jackson (1999) J. D. Jackson,  _Classical Electrodynamics_ (John Wiley & Sons, New York, 1999).
  * Bohren and Huffman (2008) C. F. Bohren and D. R. Huffman,  _Absorption and scattering of light by small particles_ (John Wiley & Sons, 2008).
  * Luk’yanchuk  _et al._ (2017) B. Luk’yanchuk, R. Paniagua-Domínguez, A. I. Kuznetsov, A. E. Miroshnichenko, and Y. S. Kivshar, Phys. Rev. A 95, 063820 (2017).
  * Watson (1995) G. N. Watson,  _A treatise on the theory of Bessel functions_ (Cambridge university press, 1995).
  * Note (1) Notice that if $\psi_{l}^{\prime}(z)=0$ then $\psi_{l}(z)\neq 0$. On other hand, if $\psi_{l}(z)=0$ then $\psi^{\prime}_{l}(z)\neq 0$. Notice that if $x_{l}^{\prime}(z)=0$ then $x_{l}(z)\neq 0$. On other hand, if $x_{l}(z)=0$ then $x^{\prime}_{l}(z)\neq 0$.
  * Zambrana-Puyalto  _et al._ (2013b) X. Zambrana-Puyalto, I. Fernandez-Corbaton, M. Juan, X. Vidal, and G. Molina-Terriza, Opt. Lett. 38, 1857 (2013b).


# Embedding a critical point in a hadron to quark–gluon crossover equation of
state

J. I. Kapusta  School of Physics & Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  C. Plumberg  Illinois Center for Advanced Studies
of the Universe, Department of Physics, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA  T. Welle  School of Physics & Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA

###### Abstract

Lattice QCD simulations have shown unequivocally that the transition from
hadrons to quarks and gluons is a crossover when the baryon chemical potential
is zero or small. Many model calculations predict the existence of a critical
point at a value of the chemical potential where current lattice simulations
are unreliable. We show how to embed a critical point in a smooth background
equation of state so as to yield the critical exponents and critical amplitude
ratios expected of a transition in the same universality class as the
liquid–gas phase transition and the 3D Ising model. The resulting equation of
state has parameters which may be inferred by hydrodynamic modeling of heavy
ion collisions in the Beam Energy Scan II at the Relativistic Heavy Ion
Collider or in experiments at other accelerators.

The QCD (Quantum Chromodynamics) equation of state has been a subject of
intense interest ever since the discovery of asymptotic freedom. At high
temperature $T$ and baryon chemical potential $\mu$ it is a weakly interacting
gas of quarks and gluons, while at low $T$ and $\mu$ it is a strongly
interacting gas of hadrons. Lattice QCD simulations have shown conclusively
that the transition from one phase to the other at $T\approx 155$ MeV and
$\mu=0$ is smooth on account of the fact that the up and down quark masses,
and consequently the pion mass, are not zero. However, diverse model
calculations predict the existence of a line of first order phase transition,
beginning at $T=0$ and $\mu_{0}$ and terminating in a critical point at
$T_{c}<155$ MeV and $\mu_{c}<\mu_{0}$ stephanov ; CPOD . Such a purported
critical point appears to be beyond the reach of reliable lattice
calculations. Experiments during the Beam Energy Scan II at RHIC (Relativistic
Heavy Ion Collider) may or may not support the existence of critical behavior
QMseries . The goal of this paper is to propose a general construction for the
equation of state which is consistent with (i) lattice QCD for all $T$ and
small $\mu$, (ii) perturbative QCD for large $T$ and/or large $\mu$, and (iii)
a critical point with critical exponents and amplitude ratios from the same
universality class as the liquid–gas phase transition and the 3D Ising model.
Parameters in this construction can be adjusted to best fit the experimental
data. The construction is different than the one proposed in Refs. attract ;
Taylor1 , which are based on the work of Refs. Zinnbook ; Guida . These two
constructions can perhaps be viewed as alternatives which provide some idea as
to the range of uncertainty in how to describe matter near the critical point
and the associated critical line of first order phase transition.

Without loss of generality the equation of state can be expressed as

$P(T,\mu)=P_{BG}(T,\mu)R(T,\mu)$ (1)

where $P_{BG}(T,\mu)$ is a judiciously chosen background equation of state
with no critical behavior; all critical behavior resides in the dimensionless
function $R$. Motivated by solutions to the cubic equation, which produce
S-shaped curves characteristic of the van der Waals equation of state or the
liquid-gas phase transition Kapusta1 ; Goodman , we consider the auxilliary
functions

$Q_{\pm}(T,\mu)=\left\\{\left[(\Delta^{2}(T))^{2}+r^{2}(T,\mu)\right]^{1/2}\pm
r(T,\mu)\right\\}^{k}$ (2)

where

$r(T,\mu)=\frac{\mu^{m}-\mu_{x}^{m}(T)}{\mu^{m}+\mu_{x}^{m}(T)}$ (3)

with $m$ a positive even integer. The function $\mu_{x}(T)$ represents the
chemical potential where the two phases are in coexistence when $T\leq T_{c}$,
but it must also be a smooth function for all $T\geq T_{c}$ to avoid undesired
discontinuities. Note that $-1\leq r<1$ and that it vanishes along the
coexistence curve. The function $\Delta^{2}(T)$ is expected to have the
functional form $|T/T_{c}-1|^{p}$ near $T_{c}$. The parameters $k$ and $p$
will determine the four critical exponents.

When $T>T_{c}$ we take

$\displaystyle R(T,\mu)=1-a(T)\left(\sqrt{\Delta^{4}+1}+1\right)^{k}$ (4)
$\displaystyle-a(T)\left(\sqrt{\Delta^{4}+1}-1\right)^{k}+a(T)(Q_{+}+Q_{-})$

where $a(T)$ is a smooth function. This has the property that $P\rightarrow
P_{BG}$ as $\mu\rightarrow 0$ and as $\mu\rightarrow\infty$ for any fixed
value of $T$. It is an even function of $\mu$. The density is

$n(T,\mu)=n_{BG}(T,\mu)R(T,\mu)+P_{BG}(T,\mu)\frac{\partial
R(T,\mu)}{\partial\mu}$ (5)

The critical exponent $\delta$ is determined by $P-P_{c}\sim{\rm
sgn}(n-n_{c})|n-n_{c}|^{\delta}$ as $n\rightarrow n_{c}$ along the critical
isotherm $\Delta^{2}=0$. Assuming $1<k<2$ the leading behavior is

$\displaystyle n-n_{c}$ $\displaystyle=$
$\displaystyle\frac{m^{k}ka(T_{c})}{\mu_{c}}P_{BG}(T_{c},\mu_{c}){\rm
sgn}(\mu-\mu_{c})\left|\frac{\mu-\mu_{c}}{\mu_{c}}\right|^{k-1}$
$\displaystyle P-P_{c}$ $\displaystyle=$
$\displaystyle(1-2^{k}a(T_{c}))n_{BG}(T_{c},\mu_{c})(\mu-\mu_{c})$ (6)

Thus the critical exponent $\delta=1/(k-1)$ or $k=1+1/\delta$.

The baryon number susceptibility is
$\chi_{B}=\partial^{2}P/\partial\mu^{2}\equiv\chi_{\mu\mu}$. It diverges like
$\chi_{B}=\chi_{+}(T/T_{c}-1)^{-\gamma}$ as $T\rightarrow T_{c}^{+}$ with
$|\mu/\mu_{c}-1|\ll T/T_{c}-1$. With $\Delta^{2}(T)=d_{+}(T/T_{c}-1)^{p}$ near
$T_{c}$ we find that the susceptibility diverges as

$\chi_{B}\rightarrow\frac{m^{2}k^{2}a(T_{c})}{2\mu_{c}^{2}}d_{+}^{\,k-2}P_{BG}(T_{c},\mu_{c})\left(\frac{T}{T_{c}}-1\right)^{-(2-k)p}$
(7)

Thus the critical exponent $\gamma=(2-k)p$ or $p=\gamma\delta/(\delta-1)$.

The heat capacity at fixed volume is

$c_{V}=T\frac{\partial s}{\partial
T}(T,n)=T\left(\chi_{TT}-\frac{\chi^{2}_{T\mu}}{\chi_{\mu\mu}}\right)$ (8)

The critical behavior is $c_{V}\rightarrow c_{+}(T/T_{c}-1)^{-\alpha}$. It can
be shown, albeit numerically, that approaching the critical point along $r=0$
gives the same result as approaching it at fixed $n_{c}$. The amplitude is

$c_{+}=2p(2kp-k-p)d_{+}^{k}\frac{a(T_{c})}{T_{c}}P_{BG}(T_{c},\mu_{c})$ (9)

with $\alpha=2-kp=2-\gamma(\delta+1)/(\delta-1)$.

When $T<T_{c}$ we take the pressure in the quark phase when $\mu\geq\mu_{x}$
to be

$P_{Q}(T,\mu)=P_{BG}(T,\mu)R_{Q}(T,\mu)$ (10)

with

$R_{Q}=1+a(T)Q_{+}(T,\mu)-a(T)\left(\sqrt{\Delta^{4}+1}+1\right)^{k}$ (11)

and in the hadron phase when $\mu\leq\mu_{x}$ to be

$P_{H}(T,\mu)=P_{BG}(T,\mu)R_{H}(T,\mu)$ (12)

with

$R_{H}=1+a(T)Q_{-}(T,\mu)-a(T)\left(\sqrt{\Delta^{4}+1}+1\right)^{k}$ (13)

We refer to these as quark and hadron phases because, even though the
background equation of state is a crossover, one is predominantly comprised of
quarks and gluons while the other is predominantly comprised of hadrons. Note
that the pressure along the critical isotherm is

$P(T_{c},\mu)=P_{BG}(T_{c},\mu)\left[1+2^{k}a(T_{c})\left(\left|\frac{\mu^{m}-\mu_{c}^{m}}{\mu^{m}+\mu_{c}^{m}}\right|^{k}-1\right)\right]$
(14)

no matter whether $T_{c}$ is approached from below or above. Hence the
critical exponent $\delta$ is well defined.

The density difference along the coexistence curve is

$\Delta
n(T)=\frac{mka(T)}{\mu_{x}(T)}P_{BG}(T,\mu_{x}(T))\left[\Delta^{2}(T)\right]^{k-1}$
(15)

The critical exponent $\beta$ is defined via $\Delta
n\sim(1-T/T_{c})^{\beta}$. With $\Delta^{2}(T)=d_{-}(1-T/T_{c})^{p}$ near
$T_{c}$ we find that $\beta=p(k-1)=\gamma/(\delta-1)$.

The susceptibility along the coexistence curve is

$\displaystyle\chi_{B}(T)$ $\displaystyle=$ $\displaystyle
P_{BG}(T,\mu_{x})\left[\frac{m^{2}k^{2}a}{4\mu_{x}^{2}}\left(\Delta^{2}\right)^{k-2}\mp\frac{mka}{2\mu_{x}^{2}}\left(\Delta^{2}\right)^{k-1}\right]$
(16) $\displaystyle+$
$\displaystyle\chi_{B,BG}(T,\mu_{x})\left[1+a(\Delta^{2})^{k}-a\left(\sqrt{\Delta^{4}+1}+1\right)^{k}\right]$
$\displaystyle\pm$
$\displaystyle\frac{mka}{\mu_{x}}n_{BG}(T,\mu_{x})\left(\Delta^{2}\right)^{k-1}$

Recalling that $\gamma=(2-k)p$ we write the critical part as
$\chi_{B}(T)=\chi_{-}(1-T/T_{c})^{-\gamma}$. When $T_{c}$ is approached from
above at $n_{c}$ the critical part is
$\chi_{B}(T)=\chi_{+}(T/T_{c}-1)^{-\gamma}$. From the above equations the
ratio of critical amplitudes is

$\frac{\chi_{+}}{\chi_{-}}=2\left(\frac{d_{-}}{d_{+}}\right)^{2-k}$ (17)

Mathematically the 2 arises because above $T_{c}$ the sum $Q_{+}+Q_{-}$ enters
whereas below $T_{c}$ only $Q_{+}$ or $Q_{-}$ does. For the universality class
which includes the liquid–gas phase transition and the 3D Ising model the
exact result is $\chi_{+}/\chi_{-}\approx 5$, whereas in mean field
approximation $\chi_{+}/\chi_{-}=2$. For the latter then $d_{+}=d_{-}$.

The last quantity to examine when $T<T_{c}$ is the heat capacity. The critical
behavior is $c_{V}(T)=c_{-}(1-T/T_{c})^{-\alpha}$ when the critical point is
approached along the coexistence curve. We find that

$c_{-}=p(2kp-k-p)d_{-}^{k}\frac{a(T_{c})}{T_{c}}P_{BG}(T_{c},\mu_{c})$ (18)

The critical exponent $\alpha$ is the same above and below $T_{c}$ as it
should be. The ratio of critical amplitudes is

$\frac{c_{+}}{c_{-}}=2\left(\frac{d_{+}}{d_{-}}\right)^{k}$ (19)

This model has two independent exponents, $k$ and $p$, in terms of which the
critical exponents $\alpha$, $\beta$, $\gamma$ and $\delta$ are expressed.
They obey the known relations $\alpha+2\beta+\gamma=2$ and
$\gamma=\beta(\delta-1)$. The true critical exponents for the universality
class which includes the liquid–gas phase transition and the 3D Ising model
are $\alpha\approx 0.1101$, $\beta\approx 0.3264$, $\gamma\approx 1.2371$, and
$\delta\approx 4.7898$ bootstrap1 ; bootstrap2  which results in $p\approx
1.564$ and $k\approx 1.209$. In the mean field approximation $\alpha=0$,
$\beta=1/2$, $\gamma=1$, and $\delta=3$ which results in $p=3/2$ and $k=4/3$.
In mean field approximation there is a discontinuity in $c_{V}$ but no
divergence. Using the true critical exponents and assuming $d_{+}/d_{-}=1/3$
yields the ratios of critical amplitudes $c_{+}/c_{-}=0.530$ and
$\chi_{+}/\chi_{-}=4.769$. Given that the amplitude ratios have an uncertainty
of a few percent these numbers are entirely consistent with published results
Zinnbook ; Guida ; Hasenbusch1 ; Hasenbusch2 .

For the background equation of state we choose the one described in Ref.
matchingpaper . It uses a switching function to transition smoothly from a
hadron resonance gas, with excluded volume interactions, to a perturbative
quark–gluon plasma. The switching function is

$S(T,\mu)=\exp\left[-\left(T^{2}/T_{s}^{2}+\mu^{2}/\mu_{s}^{2}\right)^{-2}\right]$
(20)

It ranges between 0 and 1 as $\mu$ and $T$ increase. It is an even function of
$\mu$ and infinitely differentiable so as not to introduce an artificial phase
transition of any order. The parameters $T_{s}$ and $\mu_{s}$ were adjusted to
give a good representation of lattice results latticeQCD , which results in
$T_{s}=195$ MeV and $\mu_{s}=1300$ MeV.

Only behavior in the immediate vicinity of a critical point is universal. For
heavy ion collisions it is not even clear whether one can probe it closely
enough to reveal the true critical exponents or whether mean field values are
more appropriate. Choosing the functions appearing in $R(T,\mu)$ is informed
guesswork. Here we choose the function $\mu_{x}(T)$ to follow a curve of
constant density. Other choices are possible, but this one produces an
inverted U-shape in the $T$ versus $n$ plane. It is determined implicitly via
the relation

$R(T,\mu_{x}(T))n_{BG}(T,\mu_{x}(T))=n_{c}$ (21)

This function is displayed in Figure 1. The critical point lies along this
curve. For illustration we choose $T_{c}=100$ MeV and $\mu_{c}=750$ MeV.

Figure 1: The choice of critical curve as described in the text. The critical
temperature is taken to be 100 MeV and the critical chemical potential to be
750 MeV.

As can be seen from Eq. (15), the strength of the transition is directly
proportional to the exponent $m$. We choose $m=4$ for illustration. The
remaining functions are parameterized as

$\displaystyle a(T)$ $\displaystyle=$ $\displaystyle a_{0}\exp(-T/T_{a})$
$\displaystyle\Delta^{2}(T)$ $\displaystyle=$ $\displaystyle
d_{+}(T/T_{c}-1)^{p}\exp(-T/T_{d})\;\;{\rm when}\;\;T\geq T_{c}$
$\displaystyle\Delta^{2}(T)$ $\displaystyle=$ $\displaystyle
d_{-}(1-T/T_{c})^{p}\exp(-T/T_{d})\;\;{\rm when}\;\;T\leq T_{c}$

The parameters $T_{a}$ and $T_{d}$ determine how fast $R\rightarrow 1$ as $T$
increases beyond $T_{c}$. We use the true critical exponents and amplitude
ratios with parameters $d_{-}=50$, $a_{0}=0.15$, $T_{a}=80$ MeV, and
$T_{d}=200$ MeV. Figure 2 is a contour plot of the function $R(T,\mu)$. Note
that $R\rightarrow 1$ as $\mu\rightarrow 0$ and for large $T$ and/or $\mu$.

Figure 2: The dimensionless function $R$ with parameters given in the text.
The critical point is indicated by the solid dot.

Isotherms of pressure versus density are shown in Fig. 3. Even at $T_{c}$
there is almost a plateau on account of the large critical exponent
$\delta\approx 4.79$ relative to the mean field value of 3.

Figure 3: Isotherms of pressure versus density. The numbers label the
temperature with $T_{c}=100$ MeV.

Figure 4 shows the phase transition in the temperature versus density plane.
It has the shape of an inverted U. Other shapes are possible using different
functions $\mu_{x}(T)$.

Figure 4: Temperature versus density separating the two phases.

In conclusion, we have proposed a novel way to embed critical behavior in
background equations of state which exhibit a smooth crossover from hadrons to
quarks and gluons. The goal is to use these equations of state in hydrodynamic
simulations of heavy ion collisions in order to infer whether there is
critical behavior. The approach has flexibility in selecting the location of
the critical point, the coexistence curve, and the reach of these into the
background equation of state. This flexibility is an advantage, as it allows
parameters to be adjusted to best fit experimental data. Examples were
provided; further details and exploration, such as using a background equation
of state which includes more realistic attractive and repulsive nuclear
interactions at low temperature, will be published elsewhere.

## Acknowledgments

The work of J. I. K. and T. W. was supported by the U.S. DOE Grant No. DE-
FG02-87ER40328. The work of C. P. was supported by the U.S. DOE Grant No. DE-
SC0020633.

## References

  * (1) M. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 (2004); Int. J. Mod. Phys. A 20, 4387 (2005); PoS(LAT2006)024.
  * (2) See the proceedings of the International Conference on Critical Point and Onset of Deconfinement series, the most recently being held in March 2021, to appear in PoS.
  * (3) See the proceedings of the Quark Matter Conference series, the most recently available being: Nucl. Phys. A 1005, (2021) ed. F. Liu, E. Wang, X.-N. Wang, N. Xu and B.-W. Zhang.
  * (4) C. Nonaka and M. Asakawa, Phys. Rev. C 71, 044904 (2005).
  * (5) P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang, J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Schäfer, and M. Stephanov, Phys. Rev. C 101, 034901 (2020).
  * (6) J. Zinn-Justin, Quantum Field Theory and Critical Phenomena, Clarendon Press, Oxford, 3rd edition, 1996.
  * (7) R. Guida and J. Zinn-Justin, Nucl. Phys. B489 [FS], 626 (1997).
  * (8) J. I. Kapusta, Phys. Rev. C 29, 1735 (1984).
  * (9) A. L. Goodman, J. I. Kapusta, and A. Z. Mekjian, Phys. Rev. C 30, 851 (1984).
  * (10) S. El-Showk, M. F. Paulos, D. Poland, S. Rychkov, D. Simmons-Duffin, A. Vichi, J. Stat. Phys. 157, 869 (2014).
  * (11) D. Simmons-Duffin, JHEP 06 (2015) 174.
  * (12) M. Caselle, M. Hasenbusch, J. Phys. A 30, 4963 (1997).
  * (13) M. Hasenbusch, Phys. Rev. B 82, 174434 (2010).
  * (14) M. Albright, J. Kapusta, and C. Young, Phys. Rev. C 90, 024915 (2014); 92, 044904 (2015).
  * (15) Sz. Borsányi, G. Endrödi, Z. Fodor, A. Jakovác, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabó, J. High Energy Phys. 11 (2010) 077.