# Possible Topological Superconductivity in a Topological Crystalline
Insulator (Pb1-xSnx)1-yInyTe

I. Pletikosić  Condensed Matter Physics and Materials Science Department,
Brookhaven National Lab, Upton, New York 11973, USA  Department of Physics,
Princeton University, Princeton, NJ 08544, USA  T. Yilmaz  National
Synchrotron Light Source, Brookhaven National Lab, Upton, New York 11973, USA
B. Sinkovic  Department of Physics, University of Connecticut, Storrs,
Connecticut 06269, USA
A. P. Weber  National Synchrotron Light Source, Brookhaven National Lab,
Upton, New York 11973, USA  ICFO-Institut de Ciencies Fotoniques: Barcelona,
Spain  G. D. Gu  Condensed Matter Physics and Materials Science Department,
Brookhaven National Lab, Upton, New York 11973, USA  T. Valla  Condensed
Matter Physics and Materials Science Department, Brookhaven National Lab,
Upton, New York 11973, USA  Donostia International Physics Center, 20018
Donostia-San Sebastian, Spain<EMAIL_ADDRESS>
###### Abstract

Superconductivity in topological insulators is expected to show very
unconventional features such as $p+ip$ order parameter, Majorana fermions etc.
However, the intrinsic superconductivity has been observed in a very limited
number of materials in which the pairing symmetry is still a matter of debate.
Here, we study the topological crystalline insulator (Pb1-xSnx)1-yInyTe, where
a peculiar insulator to superconductor transition was previously reported near
the gap inversion transition, when the system is nearly a 3-dimensional Dirac
semimetal. Both the existence of superconductivity near the 3-dimensional
Dirac semimetal and the occurrence of insulator to superconductor transition
in an isotropic material is highly unusual. We suggest that the observed
phenomena are related to an intrinsic instability of a 3-dimensional Dirac
semimetal state in (Pb1-xSnx)1-yInyTe and ”flattening” of the bulk valence and
conduction bands as they acquire a Mexican hat-like dispersion on the inverted
side of the phase diagram. This favors the pairing instability if the chemical
potential is pinned to these flat regions.

###### pacs:

74.25.Kc, 71.18.+y, 74.10.+v, 74.72.Hs

A 3-dimensional (3D) Dirac point with linear dispersion in all three
directions can be created at a time-reversal invariant momentum (TRIM) by
tuning an external parameter $m$, if the valence and conduction bands have the
opposite parities [1, 2, 3]. This occurs, for example, at a quantum critical
point ($m=m_{c}$) between a normal insulator and a Z2 topological insulator,
or a topological crystalline insulator. In the gap-inverted (topological)
phase, the topological states that unavoidably exist at boundaries [4, 5],
have been in the focus of both theoretical and experimental studies over the
last $\sim$ 2 decades [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].
The approach to the 3D Dirac quantum critical point requires a fine-tuning of
an external parameter ($m$) such as pressure, or material composition that
drives the transition [20, 21]. In general, the 3D Dirac point with four-fold
degeneracy does not have a topological number and it is therefore unstable
against small perturbations. However, it was later realized that particular
lattice symmetries may result in the band structure with stable 3D Dirac cones
and shortly after the theoretical prediction [22, 23] several materials
including Na3Bi and Cd3As2 were experimentally identified as 3D Dirac semi-
metals (DSM) [24, 25, 26]. While the 3D Dirac points in these 3D DSMs are
protected by crystal symmetries, the generally unprotected 3D Dirac points at
the band inversion transitions could potentially display a very rich physics.
A particularly interesting example is In-doped Pb1-xSnxTe a topological
crystalline insulator (TCI), that shows an intriguing transition from a
relatively good insulator to a robust superconductor as a function of In
concentration when in the vicinity of a transition between a trivial and
topological crystalline insulator at $x=x_{c}\approx 0.36$ [27, 28, 29, 30,
31, 32]. The insulator to superconductor transitions have been usually
observed in thin films and in some quasi-2D materials as a consequence of
strong quantum phase fluctuations due to reduced dimensionality [33, 34]. In a
relatively isotropic 3D material, such transition is highly unexpected. The
appearance of superconductivity itself is also highly unusual in a system
where the low-energy quasiparticles resemble 3D Dirac fermions with linearly
vanishing density of states.

Another interesting aspect of In-doped Pb1-xSnxTe is the posibility that its
superconducting phase represents an intrinsic topological superconductor [35].
While for a long time the focus has been on Sr2RuO4, considered to be an
intrinsic chiral $p$-wave triplet superconductor [36, 37, 38], more-recent
results excluded that possibility [39]. Currently, the number of bulk
topological superconductors is very limited, with CuxBi2Se3 and Sn(1-x)InxTe
being the rare candidates, with $T_{c}$ around $2-4$ K [40, 41, 42]. CuxBi2Se3
is found to be very inhomogeneous, with a small superconducting volume and
therefore is not ideal for spectroscopic studies that could determine the
order parameter or the existence of Majorana fermions. In more-homogeneous
Sn(1-x)InxTe, spectroscopic studies have provided evidence for the odd-parity
pairing and topological superconductivity [42]. More recently, several
candidates for topological superconductors are discovered in materials that
make natural superlattices involving rock-salt Pb1-xSnxTe layers, in which the
specific heat suggests the existence of gap nodes [43, 44].

Figure 1: Electronic structure of a thick (111) Pb0.85Sn0.15Te film. (a)
Valence band dispersion along the $\bar{\Gamma}-\bar{K}$ momentum line of the
surface Brillouin zone (SBZ) for the Te- terminated surface. (b) Constant
energy contour of the ARPES intensity at $E=-60$ meV (marked in (c) with the
line) as a function of the in-plane momentum. (c) Close-up of the spectral
region marked by the rectangle in (a). (d) and (e) show the corresponding
spectra for the Pb-Sn terminated surface. (f) Energy distribution curve (EDC)
at $k_{x}=0$ from (e). The raw intensity and that divided by the Fermi
distribution is plotted as dot-markers and solid line, respectively. The
spectra in (a-f) were recorded at $\hbar\nu=17$ eV, corresponding to
$k_{z}=2.6$ Å($L$ point in between the 2nd and 3rd bulk BZ). The arrow in (e)
represents the gap between the VB and CB, $\Delta_{L}\approx$ 130 meV. (g)
ARPES intensity for the Te-terminated surface at $E=-60$ meV as a function of
the $k_{z}$ along the $\Gamma-L-\Gamma$ of the BBZ and $k_{x}$ along the
$\bar{\Gamma}\bar{K}$ line of the SBZ. (h) In-plane dispersion at $k_{z}=2.6$
Å(same as (a)). (i) $k_{z}$ dispersion of the valence band along the
$\Gamma-L-\Gamma$ line of the bulk BZ. Intensity maps in (g,i) were produced
from spectra taken at photon energies ranging from $\hbar\nu=12.5$ eV to
$\hbar\nu=24$ eV. All the spectra were taken at $\sim 50$ K. (j) Projection of
the bulk BZ to the (111) surface and the SBZ (dashed hexagon) (left). The
small circle at the $\bar{\Gamma}$ is the same contour from (b). The 1st and
2nd bulk BZ with the probed region of $k_{z}$ marked by a thick line.  Figure
2: 3D Electronic structure of a thick (111) Pb0.6Sn0.4Te film. (a) ARPES
intensity at the Fermi level (Fermi surface) as a function of the $k_{z}$
along the $\Gamma-L-\Gamma$ of the BBZ and $k_{x}$ along t he
$\bar{\Gamma}\bar{K}$ line of the SBZ. (b) $k_{z}$ dispersion of the valence
band along the $\Gamma-L-\Gamma$ line, near the $L$ point between the $2^{nd}$
and $3^{rd}$ BBZ. (c) Dispersion of the valence band along the
$\bar{\Gamma}\bar{K}$ line of the SBZ, taken at 17 eV photon energy
(corresponding to the $W-L-W$ line in the bulk BZ). (d) The same as in (c) but
recorded at $\sim 300$ K and divided by the corresponding Fermi distribution.
(e) Constant energy contour at $E=-100$ meV as a function of $k_{x}$ and
$k_{z}$ ($k_{y}=0$). (f) Constant energy contours of in-plane ARPES intensity
at $E=-100$ meV, taken at $\hbar\nu=15$ eV (bottom), $\hbar\nu=17$ eV,
(middle) and $\hbar\nu=19$ eV (top), corresponding to three different $k_{z}$,
as indicated by dashed lines in (e). All the spectra were taken at 80 K,
except for the spectrum in (d).

An alternative avenue towards topological superconductivity in
heterostructures involving superconductors and topological insulators has been
also explored [35, 45, 46], including the TI films on cuprate high-temperature
superconductors (HTSC). However, the studies claiming a large superconducting
gap on the topological surface state on Bi2Se3 grown on Bi2212 [47, 48] were
later shown to be false [49]. The complete absence of proximity effect seems
to be intrinsic to HTSC - the films grown on cuprate superconductors donate
electrons, thus making the interface so underdoped that superconductivity
completely vanishes [50, 51], making HTSCs inadequate for inducing topological
superconductivity.

Here, we show that the topological crystalline insulator Pb1-xSnxTe might be a
good candidate for an intrinsic topological superconductor when doped with
Indium. When near the gap inversion point, In-doped Pb1-xSnxTe shows a
transition between a real bulk TCI and a bulk superconductor [31]. By studying
the details of electronic structure in angle-resolved photoemission
spectroscopy (ARPES), we suggest that this peculiar transition is caused by
the transformation of bulk bands from the linear, 3D-Dirac-like to ”Mexican
hat”-like ones on the topological side, where an inevitable flattening at the
onset of that transformation must occur. This transformation could support
superconductivity when the flat regions are brought near the chemical
potential via In-doping.

The ARPES experiments were carried out on a Scienta SES-R4000 electron
spectrometer at beamline U5UA at the National Synchrotron Light Source (NSLS),
in the photon energy range from 12.5 to 62 eV and at the OASIS-laboratory at
Brookhaven National Laboratory, using monochromatized He Iα radiation [52].
The core level spectra were taken at NSLS using 110-150 eV photons. The total
instrumental energy resolution was $\sim$ 6 meV at 12-24 eV, $\sim$ 12 meV at
50-62 eV and $\sim$ 30 meV at 110-150 eV photon energy. Angular resolution was
better than $\sim 0.15^{\circ}$ and $0.4^{\circ}$ along and perpendicular to
the slit of the analyzer, respectively.

Figure 3: 3D Electronic structure of Pb0.7Sn0.3Te and
(Pb0.7Sn0.3)0.84In0.16Te. (a) ARPES intensity at the Fermi level (Fermi
surface) as a function of the $k_{z}$ along the $\Gamma-L-\Gamma$ of the BBZ
and $k_{x}$ along the $\bar{\Gamma}\bar{K}$ line of the SBZ for Pb0.7Sn0.3Te.
(b) $k_{z}$ dispersion of the valence band along the $\Gamma-L-\Gamma$ line,
near the $L$ point between the $2^{nd}$ and $3^{rd}$ BBZ. (c) Dispersion of
the valence band along the $\bar{\Gamma}\bar{K}$ line of the SBZ, taken at 17
eV photon energy (corresponding to the $W-L-W$ line of the bulk BZ). (d)
Close-up of the spectral region marked with the rectangle in (c). (e-h) The
corresponding panels for In-doped sample, (Pb0.7Sn0.3)0.84In0.16Te. The arrow
in (h) represents the gap between the VB and CB, $\Delta_{L}\approx$ 65 meV.
(i) EDC at $k_{x}=0$ from (h). The raw intensity and that divided by the Fermi
distribution is plotted as dot-markers and solid line, respectively. The
spectra in (a-d) and (e-i) were taken at 80 and 50 K, respectively.

Thick, bulk-like (2-5 $\mu$m) films of (111) (Pb1-xSnx)1-yInyTe were grown
either on BaF2 or Bi2Te3 substrates, using the open hot wall epitaxy method
[53] and a single-source evaporators loaded with crushed single crystals of
(Pb1-xSnx)1-yInyTe [29]. The BaF2 substrates were cleaved in air and annealed
in ultra high vacuum (UHV) at 600∘C for one hour, while the Bi2Te3 substrates
were cleaved in UHV before the film deposition. During the growth, the
substrate was kept at $\sim 300^{\circ}$C. The thickness was determined by the
quartz crystal thickness monitor. Composition of the films was checked by
measuring photoemission from the shallow core levels (In $4d$, Sn $4d$, Te
$4d$ and Pb $5d$) and comparing them with spectra of the source material that
was broken in situ. After ARPES experiments, films were checked for
superconductivity by SQUID magnetometry. In several cases, the $T_{c}$ in the
film was slightly higher than in the source material, with the maximum $T_{c}$
of $\sim$6 K.

##  .1 Pure Pb1-xSnxTe

Fig. 1 shows the ARPES data from the film on the trivial side of the phase
diagram: Pb0.85Sn0.15Te. The (111) surfaces can be either Te- or Pb-Sn-
terminated and in the Pb-rich films ($x<0.3$), both terminations can be
distinguished in ARPES, as shown in Fig. 1. The panels (a-c) represent the in-
plane electronic structure of the mostly Te-terminated surface, while the
panels (d-e) represent mostly Pb-Sn termination. The two terminations could be
distinguished as they display slightly different ratio between the Te $4d$ and
Pb $5d$ (Sn $4d$) core levels. Both terminations show almost identical (aside
from the $\sim 0.2$ eV energy shift) and fairly isotropic in-plane dispersion.
Both also display bands split in a Rashba-like manner, near their extrema. The
splitting is more obvious on the Te-terminated surface, indicating a stronger
inversion symmetry breaking surface field that would result in larger spin-
orbit interaction. However, the bands disperse rather strongly with $k_{z}$
(panels (g,i) and it is not clear if the Rashba-like split features are the
surface states. Indeed, a suggestion that the bulk bands could also display
the Rashba splitting when the breaking of the inversion symmetry due to the
surface goes deep enough into the crystal was recently discussed in relation
to a new candidate for spin-triplet superconductor, UTe2 [54, 55]. Here, we
note that the momentum width of the valence band along the $k_{z}$ ($\Delta
k_{\perp}\sim 0.06$Å-1) is comparable to the in-plane width, implying the
surprisingly large ARPES probing depth at these photon energies. The more-
conventional possibility is that the observed Rashba-like split features are
trivial precursors of the topological surface states. Their apparent
dispersion with $k_{z}$, or alternatively, their existence only inside an
extremely limited $k_{z}$ interval corresponding to the very vicinity of the
$L$ point of the BZ, would suggest their strong mixing with the bulk states.
Similar effects were observed on the topological side of Pb1-xSnxSe [56].

On the surface terminated mostly by Pb-Sn, the bottom of the bulk conduction
band is visible at $\hbar\nu=17$ eV, enabling the determination of the direct
band-gap $\Delta_{L}\approx 130$ meV at 50 K, in agreement with the early
optical studies at similar compositions [57, 58]. The photon-energy dependence
that gives the access to the out-of-plane electronic structure is shown in
panels (g,i) for the mostly Te-terminated surface. The set of $k_{y}=0$
spectra measured at photon energies ranging from 12.5 to 23 eV is converted to
$k_{z}$ using the free electron approximation for the final photoemission
state, with $V_{0}=13$ eV for the “inner” potential. The emission at
$\hbar\nu=17$ eV corresponds to $k_{z}=2.6$ Å, the $L$ point between the
second and the third bulk Brillouin zone (BBZ). The constant energy contours,
shown for the Te-terminated surface (in-plane: panel (b), out-of-plane: panel
(g)) correspond to $E=-60$ meV as there is no intensity at the Fermi level.
The contour is elongated in the $k_{z}$ direction, in agreement with the band
structure calculations [20]. It also shows a limited span of the Rashba-like
features in this direction. The in-plane dispersion (panel (h)) is
approximately 2-3 times faster than the out-of plane one (panel (i)), and
aside from a $\sim 130$ meV gap, the electronic structure seems already very
close to the 3D DSM.

Fig.2 shows the 3D electronic structure of a film on the TCI side, but very
near the inversion transition: Pb0.6Sn0.4Te. At this composition, we no longer
observe either Rashba-like split states or two distinct terminations. The
reason is probably in the fact that with increasing $x$, the samples become
significantly more hole-doped and for $x\geq 0.4$, we only observe hole
pockets with a very consistent sizes. The higher carrier concentration
provides more effective screening of the surface potential, reducing the
surface band bending. Consequently, the differences for different terminations
should diminish. The Fermi surface (panel (a)) is a small ellipsoidal hole
pocket, with $\sim 3:1$ ratio in $k_{z}$ over the in-plane size. If the pocket
is the bulk feature, as the strong $k_{z}$ dispersion (panel (b)) indicates,
and by assuming the ellipsoidal shape, we can estimate the bulk carrier
concentration to be $n_{h}\sim 1.2\times 10^{18}$ cm-3 (four $L$ pockets in
the BBZ containing $3\times 10^{17}$ cm-3 holes each). The Fermi velocities,
deduced from momentum distribution curves (MDC) fitting are 4.4 eVÅ and 1.7
eVÅ for in- and out-of-plane directions, respectively. The TSS probably
overlaps with these (bulk) states as was shown in $n$-type surfaces of
Pb1-xSnxSe [59, 18, 19].

At higher energies, the contours become distorted and strongly trigonally
warped, both below and above the $L$ point (panel (f)), with the direction of
warping being opposite across the zone boundary. This was predicted in the
band structure calculations, but never observed experimentally [20]. It seems
from the spectra in Fig. 2 that this composition is indeed very close to the
critical one at which the band inversion occurs. The bands are nearly linear
and the Fermi surface is very small. Unfortunately, as the films at this
composition are always $p$-type, we do not have an access to the conduction
band at low temperature. However, at room temperature, we could access it, and
the spectrum (panel (d)) indicates that the gap is indeed negligible.

##  .2 Effects of In Doping

Figure 4: 3D Electronic structure of Pb0.5Sn0.5Te and (Pb0.5Sn0.5)0.7In0.3Te.
(a) ARPES intensity at the Fermi level (Fermi surface) as a function of the
$k_{z}$ along the $\Gamma-L-\Gamma$ of the BBZ and $k_{x}$ along the
$\bar{\Gamma}\bar{K}$ line of the SBZ for Pb0.5Sn0.5Te. (b) $k_{z}$ dispersion
of the valence band along the $\Gamma-L-\Gamma$ line, near the $L$ point
between the $2^{nd}$ and $3^{rd}$ BBZ. (c) Dispersion of the valence band
along the $\bar{\Gamma}-\bar{K}$ line of the SBZ, taken at 17 eV photon energy
(corresponding to the $W-L-W$ line of the BBZ). (d) Constant energy contour of
in-plane ARPES intensity at $E=0$ (Fermi surface) taken at $\hbar\nu=17$ eV,
corresponding to $k_{z}=2.6$ Å. (e-g) The same as (a-c), but for In-doped
sample, (Pb0.5Sn0.5)0.7In0.3Te. All the spectra were taken at 80 K. (h) Close-
up of the spectral region marked with the rectangle in (g). The white points
indicate dispersion (peak position) of the bulk valence band from fitting the
EDCs in (i). (i) EDCs of the spectral intensities from (g, h) (black) with
selected fitting curves (red). The thick spectrum corresponds to the
$k_{x}=0$, or the $L$ point in the bulk BZ.

Fig. 3 shows the ARPES results for the trivial side of the
(Pb0.7Sn0.3)1-yInyTe phase diagram, $x=0.3$, with two different In
concentrations, $y=0$ and $y=0.16$. The valence band in Fig. 3(c, d) again
displays two states split in $k$ in the Rashba-like manner. If, as previously
discussed, these are the surface states peaking in intensity at the $L$ point
of the bulk BZ, they would indicate that the material is still on the trivial
side of the phase diagram, and although they might have a significant
contribution in transport of thin samples [30], they are not protected by any
topological invariance. In that sense, they are very similar to the
conventional Rashba surface states on Au(111). Increasing the In concentration
to $y=0.16$ in Fig. 3(b) clearly shifts the Fermi level towards the bulk
conduction band and now the Rashba-like surface states can be seen on both the
conduction side and the valence side, forming a double-cone like shape, with a
small gap still present between the two sides. This shows that at this
composition, the material is very close to the critical point, at which the
two cones should collapse into a single, un-gapped one—and the topological
surface states should be fully formed. A previous ARPES study has provided
evidence for topological surface states at $x=0.4$ (without In) [15]. The
crucial point here is that the (Pb0.7Sn0.3)1-yInyTe does not become
superconducting for any amount of In-doping [27, 30, 31] up to the In
solubility limit. On the contrary, initially metallic, Pb0.7Sn0.3Te becomes a
bulk insulator with In-doping. However, at low temperature, the resistivity
becomes dominated by the contribution from the surface states and saturates in
thin samples [30, 31], a clear sign that even the trivial surface states can
significantly participate in transport in certain cases.

Fig. 4 shows the effects of In doping on the electronic properties of
Pb0.5Sn0.5Te, on the gap inverted, TCI side of the phase diagram. Without In,
ARPES indicates a heavily hole-doped system with the carrier concentration
$n_{h}\sim 5\times 10^{19}$ cm-3 (4 times $1.25\times 10^{19}$ cm-3), almost
two orders of magnitude larger than in the Pb0.6Sn0.4Te sample. In-plane Fermi
velocities show significant anisotropy, ranging from $\sim 2.2$ eVÅ to $\sim
6$ eVÅ. The slower states are likely bulk, while the fast dispersing ones are
the TSS, as their velocity is similar to those spanning the gap in the In-
doped case from Fig. 4(g.h). The out of plane Fermi velocity is $\sim 2$ eVÅ.
When doped with In, the $n_{h}$ dramatically decreases and, for
(Pb0.5Sn0.5)0.7In0.3Te, the chemical potential is shifted by $\sim 200$ meV,
from inside of the valence band to near the bottom of conduction band (note
that at this composition without indium, the inverted gap would be around 100
meV [60]). This shift in chemical potential seems to be in a nice agreement
with transport in (Pb0.5Sn0.5)1-yInyTe which displays rather interesting
behavior: initially metallic ($y=0$), the material becomes less conductive and
shows an insulating bulk conductance at low In concentration ($y\approx
0.06$). However, at even higher In concentrations ($y\geq 0.1$), the material
becomes metallic again, this time with the superconducting ground state [29,
31]. This non-monotonic behavior, together with our spectroscopic results,
would suggests that the carriers, initially holes, turn into electrons with
sufficient In doping. This is similar to what has been observed in Bi-doped
material where the measured Hall coefficients support the transition from the
hole to the electron conduction, albeit at much lower nominal dopant
concentrations [61].

Furthermore, the bulk valence band, now fully occupied, displays a local
minimum at $L$ point and two nearby maxima as it disperses along the $k_{x}$
(Fig. 4(h, i)). Its dispersion is obtained from the peak position of the
higher binding energy peaks in EDCs, as represented by white points in panel
(h).

So, the question is why is the In-doped Pb1-xSnxTe superconducting only on the
gap-inverted, TCI side of the phase diagram [27]? After all, pure Pb has a
higher $T_{c}$ than Sn and one could naively expect that the superconductivity
will be stronger near the $x=0$. We think that the answer lies in the
different shape of the bands before and after the gap inversion, as presented
in Fig. 5. On the trivial side, the valence and conduction bands have a single
extremum along all three $k$-axis (Fig1). On approaching the critical point at
$x\approx 0.36$, these extrema approach each other and get lighter until they
touch in a single (3D Dirac) point. On the topological side, and far from the
inversion point, the bands have so called a Mexican hat band (MHB) shape [62,
63, 17, 64]. Specifically, in the case of SnTe, both the valence and
conducting bands have a single maximum (minimum) at L point as they disperse
along the $\Gamma-L$ line, whereas in the perpendicular direction, or within
the (111) plane, the bands are MHB, as realized early on by Tung et al and
Rabii [65, 62]. This is consistent with our experimental results (Fig. 4
(e-j)). The density of states (DOS) near the tops (bottoms) of such MHB can be
dramatically enhanced relative to the band extrema on the trivial side, with
the details depending on the models and dimensionality of the system [63, 17,
64]. Moreover, the transition from a 3D Dirac (linear) bands to the MHB
involves not only mass generation (gap opening), but also a flattening related
with a transition from a single extremum to a MHB, that can further enhance
the DOS and lead to singularities even in 3D systems [63]. We argue that this
will necessarily have an augmenting effect on the properties of a material
dependent on DOS, as it will enhance electronic correlations, especially if
the enhanced DOS can be brought near the chemical potential. In particular, we
can expect that the electron-phonon coupling would be enhanced, eventually
leading to superconductivity. The flat bands near the Fermi level are
generally unstable against other symmetry-breaking phases, such as charge
density waves and magnetism, but in In-doped Pb1-xSnxTe, the superconductivity
is the ground state tied to the MHB. Moreover, since the MHBs are inverted,
the materials that possess them are excellent candidates for intrinsic
topological superconductors.

Figure 5: Phase diagram of (Pb1-xSnx)1-yInyTe. (a) Experimentally obtained in-
plane dispersions of occupied (tinted) bulk states in Pb0.85Sn0.15Te (left),
Pb0.6Sn0.4Te (middle) and (Pb0.5Sn0.5)0.7In0.3Te (right). The fitting of
momentum distribution curves (MDCs) was performed for the first two cases,
while EDC fitting was used in the third case. The unoccupied dispersions are
inverted copies of the occupied ones. (b) The phase diagram of
(Pb1-xSnx)1-yInyTe. The solid (dashed) line represents the fundamental gap,
$E_{g}$, at $T=0$ ($T=50$ K) in Pb1-xSnxTe from Heremans et al [60]. The black
circles are gap estimates from our measurements, as discussed in the main
text. Superconducting $T_{c}$ data points for In-doped Pb1-xSnxTe (red
symbols) are taken from ref. [27].

We note that both the gap magnitudes and the details in dispersions of
valence/conduction bands shown in Fig. 5 depend on the analysis. They also
depend on temperature, both in the theoretical and experimental studies. Our
clearest cases are Pb0.85Sn0.15Te and (Pb0.7Sn0.3)0.84In0.16Te, where both
valence and conduction bands are visible, suggest that the most appropriate
method for gap determination is by the steep onsets of emission at band
extrema (Fig. 1(e,f) and Fig. 3(h,i)). For the remaining samples, the gap
determination is less reliable. For example, in Pb0.6Sn0.4Te, the room
temperature measurements displayed no detectable gap and we set the value of
$\sim 30$ meV as the upper limit for low-temperature gap. In the pure
Pb0.3Sn0.7Te sample, the conduction band is clearly out of the ARPES reach
(Fig. 3(c,d)). However, we were able to detect its threshold at room
temperature after depositing potassium on the surface. The extracted gap was
$\sim$150 meV, much larger than the one seen in In doped sample at low
temperature, but when adjusted for temperature, by using the known thermal
coefficients, the gaps were almost identical ($\sim 65$ meV at 50 K). Also,
the room temperature spectrum of the Pb0.4Sn0.6Te sample shows essentially no
gap, but we have set the 30 meV as an upper limit at the low temperature. In
the (Pb0.5Sn0.5)0.7In0.3Te case, the gap is determined from the dispersion of
the peak maximum of the valence band and the assumption that the chemical
potential is near the bottom of the conduction band (Fig.4(h,i)). From the
Pb0.85Sn0.15Te case, we know that the distance between the sharp onset of
emission and the peak maximum can be quite significant ($\sim 50$ meV) which
would lead to a large overestimate for the gaps determined from peak
positions. Therefore, in Fig. 5, we have adjusted the gap magnitudes obtained
from peak dispersions to those that would correspond to the onsets of
emission, both for $x=0.15$ and $x=0.5$. Also, there is some ambiguity in how
the gap is actually defined in the inverted case, where the bands display the
MHB dispersion. However, the most important finding of our study, that the In-
doped Pb1-xSnxTe displays the MHB dispersion on the TCI side of the phase
diagram, where superconductivity occurs, does not depend on those details.

On the other hand, the question whether the superconductivity in
(Pb1-xSnx)1-yInyTe is conventional or topological is still open. Topological
superconductors are accompanied by gapless states at their boundaries, and
they may be composed of Majorana fermions. The point-contact spectroscopy on
In-doped SnTe suggests the odd-parity pairing and topological
superconductivity [42], but the specific heat reveals a fully gapped
superconductivity that would argue against it [66]. Similarly, the STM
measurements of the superconducting gap in (Pb0.5Sn0.5)0.7In0.3Te on the (001)
surface suggest that the sample seems to be fully gapped, without in-gap
states—contrary to what would be expected for a topological superconductor
[67]. However, we note that those measurements were performed on (100)
surfaces of bulk crystals that poorly cleave. In addition, topological surface
states on the (100) surfaces are different from those on the (111) surfaces
studied here and their ability to support the zero modes might also be
different [16]. Therefore, it would be highly desirable to perform STM studies
on the (111) surfaces. Our ARPES results unambiguously show that the
superconductivity in In-doped Pb1-xSnxTe occurs not only on the topological
side of the phase diagram, but also because of it, with the inverted MHBs
playing the decisive role. Therefore, the In-doped Pb1-xSnxTe is a TCl (for
$x\geq 0.36$) and is also a superconductor, but it remains to be seen whether
it is a topological superconductor with an odd order parameter and in-gap
states. We hope that our results will stimulate further studies capable of
answering these questions directly.

###### Acknowledgements.

This work was supported by the US Department of Energy, Office of Basic Energy
Sciences, contract no. DE-SC0012704 and ARO MURI program, grant
W911NF-12-1-0461.

## References

  * Shuichi Murakami [2007] Shuichi Murakami, New Journal of Physics 9, 356 (2007), ISSN 1367-2630, URL http://stacks.iop.org/1367-2630/9/i=9/a=356.
  * Murakami et al. [2007] S. Murakami, S. Iso, Y. Avishai, M. Onoda, and N. Nagaosa, Physical Review B 76, 205304 (2007), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.76.205304.
  * Yang and Nagaosa [2014] B.-J. Yang and N. Nagaosa, Nature Communications 5, 4898 (2014), ISSN 2041-1723, URL http://www.nature.com/articles/ncomms5898.
  * Volkov and Pankratov [1985] B. A. Volkov and O. A. Pankratov, JETP Letters 42, 178 (1985), URL http://www.jetpletters.ac.ru/ps/1420/article_21570.shtml.
  * Fradkin et al. [1986] E. Fradkin, E. Dagotto, and D. Boyanovsky, Physical Review Letters 57, 2967 (1986), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.57.2967.
  * Murakami et al. [2004] S. Murakami, N. Nagaosa, and S.-C. Zhang, Physical Review Letters 93, 156804 (2004), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.93.156804.
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Physical Review Letters 95, 146802 (2005), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.95.146802.
  * Bernevig et al. [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science (New York, N.Y.) 314, 1757 (2006), ISSN 1095-9203, URL http://www.ncbi.nlm.nih.gov/pubmed/17170299.
  * Konig et al. [2007] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X. L. Qi, and S. C. Zhang, Science 318, 766 (2007).
  * Fu et al. [2007] L. Fu, C. Kane, and E. Mele, Physical Review Letters 98, 106803 (2007), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.98.106803.
  * Noh et al. [2008] H.-J. Noh, H. Koh, S.-J. Oh, J.-H. Park, H.-D. Kim, J. D. Rameau, T. Valla, T. E. Kidd, P. D. Johnson, Y. Hu, et al., EPL (Europhysics Letters) 81, 57006 (2008), ISSN 0295-5075, URL http://stacks.iop.org/0295-5075/81/i=5/a=57006?key=crossref.3e71f59820d47157a6627f67beccd3c3.
  * König et al. [2008] M. König, H. Buhmann, L. W. Molenkamp, T. Hughes, C.-X. Liu, X.-L. Qi, and S.-C. Zhang, Journal of the Physical Society of Japan 77, 031007 (2008), ISSN 0031-9015, URL http://jpsj.ipap.jp/link?JPSJ/77/031007/.
  * Hsieh et al. [2008] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Nature 452, 970 (2008).
  * Fu [2011] L. Fu, Physical Review Letters 106, 106802 (2011), ISSN 0031-9007, URL https://link.aps.org/doi/10.1103/PhysRevLett.106.106802.
  * Xu et al. [2012] S.-Y. Xu, C. Liu, N. Alidoust, M. Neupane, D. Qian, I. Belopolski, J. D. Denlinger, Y. J. Wang, H. Lin, L. A. Wray, et al., Nature communications 3, 1192 (2012), ISSN 2041-1723, URL http://dx.doi.org/10.1038/ncomms2191.
  * Liu et al. [2013] J. Liu, W. Duan, and L. Fu, Physical Review B 88, 241303 (2013), ISSN 1098-0121, URL http://link.aps.org/doi/10.1103/PhysRevB.88.241303.
  * Yang et al. [2014] G. Yang, J. Liu, L. Fu, W. Duan, and C. Liu, Physical Review B 89, 085312 (2014), ISSN 1098-0121, URL http://link.aps.org/doi/10.1103/PhysRevB.89.085312.
  * Pletikosić et al. [2014] I. Pletikosić, G. D. Gu, and T. Valla, Physical review letters 112, 146403 (2014), ISSN 1079-7114, URL http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.112.146403.
  * Gyenis et al. [2013] A. Gyenis, I. K. Drozdov, S. Nadj-Perge, O. B. Jeong, J. Seo, I. Pletikosić, T. Valla, G. D. Gu, and A. Yazdani, Physical Review B 88, 125414 (2013), ISSN 1098-0121, URL http://journals.aps.org/prb/abstract/10.1103/PhysRevB.88.125414.
  * Svane et al. [2010] A. Svane, N. E. Christensen, M. Cardona, A. N. Chantis, M. van Schilfgaarde, and T. Kotani, Physical Review B 81, 245120 (2010), ISSN 1098-0121, URL http://link.aps.org/doi/10.1103/PhysRevB.81.245120.
  * Suchalkin et al. [2020] S. Suchalkin, M. Ermolaev, T. Valla, G. Kipshidze, D. Smirnov, S. Moon, M. Ozerov, Z. Jiang, Y. Jiang, S. P. Svensson, et al., Applied Physics Letters 116, 032101 (2020), ISSN 0003-6951, URL http://aip.scitation.org/doi/10.1063/1.5128634.
  * Wang et al. [2012a] Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Physical Review B 85, 195320 (2012a), ISSN 1098-0121, URL http://link.aps.org/doi/10.1103/PhysRevB.85.195320.
  * Wang et al. [2013a] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Physical Review B 88, 125427 (2013a), ISSN 1098-0121, URL http://link.aps.org/doi/10.1103/PhysRevB.88.125427.
  * Liu et al. [2014] Z. K. Liu, B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D. Prabhakaran, S.-K. Mo, Z. X. Shen, Z. Fang, X. Dai, et al., Science (New York, N.Y.) 343, 864 (2014), ISSN 1095-9203, URL http://www.sciencemag.org/content/343/6173/864.abstract.
  * Borisenko et al. [2014] S. Borisenko, Q. Gibson, D. Evtushinsky, V. Zabolotnyy, B. Büchner, and R. J. Cava, Physical Review Letters 113, 027603 (2014), ISSN 0031-9007, URL https://link.aps.org/doi/10.1103/PhysRevLett.113.027603.
  * Neupane et al. [2014] M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin, et al., Nature Communications 5, 3786 (2014), ISSN 2041-1723, URL http://www.nature.com/articles/ncomms4786.
  * Bushmarina et al. [1991] G. Bushmarina, I. Drabkin, D. Mashovets, R. Parfeniev, D. Shamshur, and M. Shachov, Physica B: Condensed Matter 169, 687 (1991), ISSN 09214526, URL http://www.sciencedirect.com/science/article/pii/092145269190389V.
  * Kozub et al. [2006] V. I. Kozub, R. V. Parfen’ev, D. V. Shamshur, D. V. Shakura, A. V. Chernyaev, and S. A. Nemov, JETP Letters 84, 35 (2006), ISSN 0021-3640, URL http://link.springer.com/10.1134/S002136400613008X.
  * Zhong et al. [2014] R. D. Zhong, J. A. Schneeloch, T. S. Liu, F. E. Camino, J. M. Tranquada, and G. D. Gu, Physical Review B 90, 020505 (2014), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.90.020505.
  * Zhong et al. [2015] R. Zhong, X. He, J. A. Schneeloch, C. Zhang, T. Liu, I. Pletikosić, T. Yilmaz, B. Sinkovic, Q. Li, W. Ku, et al., Physical Review B 91, 195321 (2015), ISSN 1098-0121, URL http://journals.aps.org/prb/abstract/10.1103/PhysRevB.91.195321.
  * Zhong et al. [2017] R. Zhong, J. Schneeloch, Q. Li, W. Ku, J. Tranquada, and G. Gu, Crystals 7, 55 (2017), ISSN 2073-4352, URL http://www.mdpi.com/2073-4352/7/2/55.
  * Sapkota et al. [2020] A. Sapkota, Y. Li, B. L. Winn, A. Podlesnyak, G. Xu, Z. Xu, K. Ran, T. Chen, J. Sun, J. Wen, et al., Physical Review B 102, 104511 (2020), ISSN 2469-9950, URL https://link.aps.org/doi/10.1103/PhysRevB.102.104511.
  * Lin et al. [2015] Y. H. Lin, J. Nelson, and A. M. Goldman, Physica C: Superconductivity and its Applications 514, 130 (2015), ISSN 0921-4534.
  * Chudnovskiy and Kamenev [2022] A. Chudnovskiy and A. Kamenev, Physical Review Letters 129, 266601 (2022), ISSN 0031-9007, URL https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.129.266601.
  * Fu and Kane [2008] L. Fu and C. Kane, Physical Review Letters 100, 096407 (2008), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.100.096407.
  * Rice and Sigrist [1995] T. M. Rice and M. Sigrist, Journal of Physics: Condensed Matter 7, L643 (1995), ISSN 0953-8984, URL https://iopscience.iop.org/article/10.1088/0953-8984/7/47/002.
  * Ishida et al. [1998] K. Ishida, H. Mukuda, Y. Kitaoka, K. Asayama, Z. Q. Mao, Y. Mori, and Y. Maeno, Nature 396, 658 (1998), ISSN 0028-0836, URL http://www.nature.com/articles/25315.
  * Das Sarma et al. [2006] S. Das Sarma, C. Nayak, and S. Tewari, Physical Review B 73, 220502 (2006), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.73.220502.
  * Pustogow et al. [2019] A. Pustogow, Y. Luo, A. Chronister, Y.-S. Su, D. A. Sokolov, F. Jerzembeck, A. P. Mackenzie, C. W. Hicks, N. Kikugawa, S. Raghu, et al., Nature 574, 72 (2019), ISSN 0028-0836, URL http://www.nature.com/articles/s41586-019-1596-2.
  * Hor et al. [2010] Y. S. Hor, A. J. Williams, J. G. Checkelsky, P. Roushan, J. Seo, Q. Xu, H. W. Zandbergen, A. Yazdani, N. P. Ong, and R. J. Cava, Physical Review Letters 104, 057001 (2010), ISSN 0031-9007, URL http://prl.aps.org/abstract/PRL/v104/i5/e057001.
  * Erickson et al. [2009] A. S. Erickson, J.-H. Chu, M. F. Toney, T. H. Geballe, and I. R. Fisher, Physical Review B 79, 024520 (2009), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.79.024520.
  * Sasaki et al. [2012] S. Sasaki, Z. Ren, A. A. Taskin, K. Segawa, L. Fu, and Y. Ando, Physical Review Letters 109, 217004 (2012), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.109.217004.
  * Sasaki et al. [2014] S. Sasaki, K. Segawa, and Y. Ando, Physical Review B 90, 220504 (2014), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.90.220504.
  * Luo et al. [2016] H. Luo, K. Yan, I. Pletikosic, W. Xie, B. F. Phelan, T. Valla, and R. J. Cava, Journal of the Physical Society of Japan 85, 064705 (2016), ISSN 0031-9015, URL https://journals.jps.jp/doi/10.7566/JPSJ.85.064705.
  * Wang et al. [2012b] M.-X. Wang, C. Liu, J.-P. Xu, F. Yang, L. Miao, M.-Y. Yao, C. L. Gao, C. Shen, X. Ma, X. Chen, et al., Science 336, 52 (2012b), ISSN 0036-8075, URL https://www.science.org/doi/10.1126/science.1216466.
  * Yi et al. [2022] H. Yi, L.-H. Hu, Y. Wang, R. Xiao, J. Cai, D. R. Hickey, C. Dong, Y.-F. Zhao, L.-J. Zhou, R. Zhang, et al., Nature Materials pp. 1–7 (2022), ISSN 1476-1122, URL https://www.nature.com/articles/s41563-022-01386-z.
  * Zareapour et al. [2012] P. Zareapour, A. Hayat, S. Y. F. Zhao, M. Kreshchuk, A. Jain, D. C. Kwok, N. Lee, S.-W. Cheong, Z. Xu, A. Yang, et al., Nature communications 3, 1056 (2012), ISSN 2041-1723, URL http://dx.doi.org/10.1038/ncomms2042.
  * Wang et al. [2013b] E. Wang, H. Ding, A. V. Fedorov, W. Yao, Z. Li, Y.-F. Lv, K. Zhao, L.-G. Zhang, Z. Xu, J. Schneeloch, et al., Nature Physics 9, 621 (2013b), ISSN 1745-2473, URL http://dx.doi.org/10.1038/nphys2744.
  * Yilmaz et al. [2014] T. Yilmaz, I. Pletikosić, A. P. Weber, J. T. Sadowski, G. D. Gu, A. N. Caruso, B. Sinkovic, and T. Valla, Physical review letters 113, 067003 (2014), ISSN 1079-7114, URL http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.113.067003.
  * Kundu et al. [2020] A. K. Kundu, Z. Wu, I. Drozdov, G. Gu, and T. Valla, Advanced Quantum Technologies 3, 2000038 (2020), ISSN 2511-9044, URL https://onlinelibrary.wiley.com/doi/10.1002/qute.202000038.
  * Kundu et al. [2021] A. K. Kundu, G. Gu, and T. Valla, ACS Applied Materials and Interfaces 13, 33627 (2021), ISSN 1944-8244, URL https://pubs.acs.org/doi/10.1021/acsami.1c06821.
  * Kim et al. [2022] C. K. Kim, I. K. Drozdov, K. Fujita, J. C. S. Davis, I. Božović, and T. Valla, Journal of Electron Spectroscopy and Related Phenomena 257, 146775 (2022), ISSN 0368-2048, URL https://www.sciencedirect.com/science/article/pii/S0368204818300884?via%3Dihub.
  * Strauss [1967] A. Strauss, Physical Review 157, 608 (1967), ISSN 0031-899X, URL http://link.aps.org/doi/10.1103/PhysRev.157.608.
  * Ran et al. [2019] S. Ran, C. Eckberg, Q. P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I. L. Liu, M. Zic, H. Kim, J. Paglione, et al., Science 365, 684 (2019), ISSN 10959203, URL https://www.science.org/doi/10.1126/science.aav8645.
  * Yu et al. [2022] Y. Yu, V. Madhavan, and S. Raghu, Physical Review B 105, 174520 (2022), ISSN 2469-9950, URL https://link.aps.org/doi/10.1103/PhysRevB.105.174520.
  * Polley et al. [2014] C. M. Polley, P. Dziawa, A. Reszka, A. Szczerbakow, R. Minikayev, J. Z. Domagala, S. Safaei, P. Kacman, R. Buczko, J. Adell, et al., Physical Review B 89, 075317 (2014), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.89.075317.
  * Dimmock et al. [1966] J. Dimmock, I. Melngailis, and A. Strauss, Physical Review Letters 16, 1193 (1966), ISSN 0031-9007, URL http://link.aps.org/doi/10.1103/PhysRevLett.16.1193.
  * Takaoka et al. [1989] S. Takaoka, S. Shimomura, H. Takahashi, and K. Murase, Physical Review B 40, 5642 (1989), ISSN 0163-1829, URL http://link.aps.org/doi/10.1103/PhysRevB.40.5642.
  * Turowski et al. [2023] B. Turowski, A. Kazakov, R. Rudniewski, T. Sobol, E. Partyka-Jankowska, T. Wojciechowski, M. Aleszkiewicz, W. Zaleszczyk, M. Szczepanik, T. Wojtowicz, et al., Applied Surface Science 610, 155434 (2023), ISSN 0169-4332.
  * Heremans et al. [2012] J. P. Heremans, B. Wiendlocha, and A. M. Chamoire, Energy and Environmental Science 5, 5510 (2012), ISSN 1754-5706, URL https://pubs.rsc.org/en/content/articlehtml/2012/ee/c1ee02612ghttps://pubs.rsc.org/en/content/articlelanding/2012/ee/c1ee02612g.
  * Volobuev et al. [2017] V. V. Volobuev, P. S. Mandal, M. Galicka, O. Caha, J. Sánchez-Barriga, D. Di Sante, A. Varykhalov, A. Khiar, S. Picozzi, G. Bauer, et al., Advanced Materials 29, 1604185 (2017), ISSN 09359648, URL http://doi.wiley.com/10.1002/adma.201604185.
  * Tung and Cohen [1969] Y. W. Tung and M. L. Cohen, Physics Letters A 29, 236 (1969), ISSN 0375-9601.
  * He et al. [2014] X.-G. He, X. Xi, and W. Ku (2014), eprint 1410.2885, URL https://arxiv.org/abs/1410.2885v1.
  * Xu et al. [2020] H. Xu, J. Zhou, H. Wang, and J. Li, The Journal of Physical Chemistry Letters 11, 6119 (2020), ISSN 1948-7185, URL https://pubs.acs.org/doi/10.1021/acs.jpclett.0c01552.
  * Rabii [1969] S. Rabii, Physical Review 182, 821 (1969), ISSN 0031899X, URL https://journals.aps.org/pr/abstract/10.1103/PhysRev.182.821.
  * Novak et al. [2013] M. Novak, S. Sasaki, M. Kriener, K. Segawa, and Y. Ando, Physical Review B 88, 140502 (2013), ISSN 1098-0121, URL http://link.aps.org/doi/10.1103/PhysRevB.88.140502.
  * Du et al. [2015] G. Du, Z. Du, D. Fang, H. Yang, R. D. Zhong, J. Schneeloch, G. D. Gu, and H.-H. Wen, Physical Review B 92, 020512 (2015), ISSN 1098-0121, URL https://link.aps.org/doi/10.1103/PhysRevB.92.020512.


}k=1,\ldots,3\\\ 0&\text{ if }k=4,\ldots,9,\end{array}\right.$

Moreover,

$\|H(\mathcal{Z}_{i-1})\|=\sum_{k=1}^{9}\beta_{k}^{*}|\mathcal{Z}_{i-1,k}|+o(\|\mathcal{Z}_{i-1}\|),\qquad\text{
as }\;\left\|\mathcal{Z}_{i-1}\right\|\rightarrow\infty,$

where $\beta_{k}^{*}=|\beta_{2}|^{1/2}$ for $k=1,\ldots,3$ and zero otherwise.
Thus, from condition (B4) in Lu and Jiang [28], a sufficient condition to
guarantee the geometric ergodicity is

$\max_{1\leq k\leq
9}\left\\{\gamma_{k}^{*}+\gamma_{k}^{**}+\beta_{k}^{*}\mathbb{E}[|\upsilon_{i,k}|]\right\\}<1,$
(44)

where the $\upsilon_{i,k}$ stand for the components of $\upsilon_{i}$.

With our simulation design, for $k=4,\ldots,9$,

$\gamma_{k}^{*}+\gamma_{k}^{**}+\beta_{k}^{*}\mathbb{E}[|\upsilon_{i,k}|]=|\rho|+0+0=\frac{1}{4}<1.$

Hence, condition (44) rewrites as

$\max_{1\leq k\leq
3}\left\\{|\gamma_{11}|+|\gamma_{12}|+|\beta_{2}|^{1/2}\mathbb{E}[|\upsilon_{i,k}|]\right\\}<1$

Besides, noting that for each $1\leq k\leq 3$,
$\mathbb{E}^{2}[|\upsilon_{i,k}|]\leq\mathbb{E}[\upsilon_{i,k}^{2}]=\mathbb{E}[\zeta_{i}^{2}]=1$,
we obtain that for any $k=1,\ldots,3$,

$|\gamma_{11}|+|\gamma_{12}|+|\beta_{2}|^{1/2}\mathbb{E}[|\upsilon_{i,k}|]\leq\frac{1}{10}+0+\frac{1}{\sqrt{10}}.$

Therefore, condition (44) is fulfilled.

###  B.3 A conditional variance semiparametric model

In this section we extend the scope of our models. As mentioned in section 6,
$Y_{i}$ could be observed with some error. For illustration, consider a time
series $(R_{i})$, solution of the AR(1) equation

$R_{i}=\rho_{0}R_{i-1}+u_{i},\qquad i\in\mathbb{Z},$ (45)

Consider the PLSIM

$Y_{i}=u_{i}^{2}=\mu(V_{i};\gamma,m)+\varepsilon_{i}\;\textrm{ with
}\;\mu(V_{i};\gamma,m)=l(X_{i};\gamma_{1})+m(W_{i}^{\top}\gamma_{2})>0,$ (46)

with $\mathbb{E}[\varepsilon_{i}\mid V_{i},\mathcal{F}_{i-1}]=0$ a.s. and
$(V_{i}^{\top},\varepsilon_{i})^{\top}\in\mathbb{R}^{d_{X}+d_{W}}\times\mathbb{R}$
a strictly stationary and strongly mixing sequence. Moreover,
$\mathbb{E}[u_{i}\mid V_{i},\mathcal{F}_{i-1}]=0$ a.s. A more common way to
write model (46) is

$u_{i}=\sqrt{\mu_{i}}\;\nu_{i},$ (47)

with $(\nu_{i})$ a strong white noise process with unit variance, and
$\mu_{i}$ a positive function of the past values of $u_{i}$. The ARCH model is
a typical example. When covariates are also allowed to enter the expression of
$\mu_{i}$, one obtains a particular example of the so called GARCH-X models.
See Han and Kristensen (2014). Here we allow for a flexible semiparametric
form $\mu_{i}=\mu(V_{i};\gamma,m)$ and our additive error term is
$\varepsilon_{i}=\mu(V_{i};\gamma,m)(\nu_{i}^{2}-1).$

Although the conditional variance of $Y_{i}$ is not constant and the $Y_{i}$
are not directly observed, the PLSIM is still applicable, as we will briefly
justify in the following. Instead of $Y_{i}$, one has

$\widetilde{Y}_{i}=\left(R_{i}-\widetilde{\rho}R_{i-1}\right)^{2}=Y_{i}+R_{i-1}^{2}(\widetilde{\rho}-\rho_{0})^{2}-2u_{i}R_{i-1}(\widetilde{\rho}-\rho_{0}).$

Here, $\widetilde{\rho}$ is the least-squares estimator of $\rho_{0}$. Let
$\widehat{\widetilde{\eta}}_{\gamma}$ be the vector of nonparametric
estimators defined in (18) obtained with $\widetilde{Y}_{i}$ instead of
$Y_{i}$. Only the components $\widehat{\widetilde{\eta}}_{\gamma,m}$ and
$\widehat{\widetilde{\eta}}_{\gamma,m^{\prime}}$ are affected by the fact that
the $Y_{i}$’s are not available. Given the expression of
$\widetilde{Y}_{i}-Y_{i}$ we deduce

$\widehat{\widetilde{\eta}}_{\gamma,m}(t)=\widehat{\eta}_{\gamma,m}(t)-\frac{2(\widetilde{\rho}-\rho_{0})}{nh}\sum_{i=1}^{n}u_{i}R_{i-1}K\left(\frac{W_{i}^{\top}\gamma_{2}-t}{h}\right)+(\widetilde{\rho}-\rho_{0})^{2}O_{\mathbb{P}}(1),$

uniformly with respect to $t$, and a similar representation holds true for
$\widehat{\widetilde{\eta}}_{\gamma,m^{\prime}}$. Using the fact that
$\widetilde{\rho}-\rho_{0}=O_{\mathbb{P}}(n^{-1/2})$ and $\mathbb{E}[u_{i}\mid
V_{i},\mathcal{F}_{i-1}]=0$ a.s., the arguments in the proof of Theorem 1
remain valid and the limit of the ELR is still a chi-square distribution. The
technical details are quite straightforward and thus are omitted. Instead, we
propose an illustration using simulation data.

We generated data from model (45)-(47). First, we generate the covariates
$W_{i}=(W_{i1},W_{i2},W_{i3})^{\top}$ from a multivariate Gaussian
distribution with mean $W_{i-1}/4$ and covariance given by
$\textrm{cov}(W_{ik},W_{i\ell})=0.5^{|k-\ell|}$. Then, we set
$X_{i}=(U_{i-1}^{2},U_{i-2}^{2})^{\top}$ and we generate the $U_{i}^{2}$ from
(46)- (47) with

$\ell(X_{i};\gamma_{1})=\gamma_{11}U_{i-1}^{2}+\gamma_{12}U_{i-2}^{2}\quad\text{
and }\quad m(u)=\frac{1}{4}+\frac{3}{4}\sin^{2}(u\pi),$ (48)

with $\gamma_{1}=(0.1,0)^{\top}$, $\gamma_{2}=(1,1,1)^{\top}$. Finally, the
variables $R_{i}$ can be computed from (45) with $\rho_{0}=0.1$. Based on the
observed $R_{i}$’s, we can compute $\tilde{Y}_{i}$ and use the proposed EL
procedure for testing the order of the lagged values to consider in the
parametric function. Hypothesis testing is based on Wilks’ Theorem in Section
4.3 (results related to this method are named _estim_), along with the
unfeasible EL approach that previously learns the nonparametric estimators on
a sample of size $10^{4}$ (this case mimic the situation where $m$,
$m^{\prime}$ and the density of the index are known; results related to this
method are named _ref_). The nonparametric elements are estimated by the
Nadaraya-Watson method with Gaussian kernel and bandwidth $h=C^{-1}n^{-1/5}$
where $C$ is the standard deviation of the index. Thus, we consider the tests
introduced in Section 5. The empirical probabilities of rejection are
presented in Table 5 for a nominal level of $0.05$.

Table 5: Empirical probabilities of rejection obtained from 5000 replications using the PLSIM for testing the order for the lagged values of $Y_{i}$ in the parametric part $\ell(\cdot;\gamma_{1})$. Test | $n=1000$ | $n=2000$ | $n=4000$ | $n=8000$
---|---|---|---|---
| ref. | estim. | ref. | estim. | ref. | estim. | ref. | estim.
Lag(1) | 0.117 | 0.120 | 0.089 | 0.097 | 0.081 | 0.078 | 0.064 | 0.063
Lag(0) | 0.422 | 0.402 | 0.678 | 0.654 | 0.942 | 0.932 | 1.000 | 0.999
Lag(2) | 0.538 | 0.531 | 0.695 | 0.687 | 0.887 | 0.874 | 1.000 | 1.000

References

[1] Han, H. and D. Kristensen (2014). Asymptotic for the QMLE in GARCH-X
Models With Stationary and Nonstationary Covariates. _Journal of Business &
Economic Statistics 32_, 416–429.


# Using DFTB to Model Photocatalytic Anatase-Rutile TiO2 Nanocrystalline
Interfaces and their Band Alignment

Verena Kristin Gupta  Bremen Center for Computational Materials Science,
University of Bremen, P.O. Box 330440, D-28334 Bremen, Germany  Bálint Aradi
Bremen Center for Computational Materials Science, University of Bremen, P.O.
Box 330440, D-28334 Bremen, Germany  Kyoung Kweon  [  Nathan Keilbart  [  Nir
Goldman  [  Thomas Frauenheim  Computational Science Research Center, No.10
East Xibeiwang Road, Beijing 100193 and Computational Science and Applied
Research Institute Shenzhen, China<EMAIL_ADDRESS>Jolla Kullgren  Department of Chemistry, Structural Chemistry, Angström
Laboratory, Uppsala University, Box 538, 752 21, Uppsala, Sweden
<EMAIL_ADDRESS>
###### Abstract

Band alignment effects of anatase and rutile nanocrystals in TiO2 powders lead
to an electron hole separation, increasing the photo catalytic efficiency of
these powders. While size effects and types of possible alignments have been
extensively studied, the effect of interface geometries of bonded nanocrystal
structures on the alignment is poorly understood. In order to allow conclusive
studies of a vast variety of bonded systems in different orientations, we have
developed a new density functional tight binding parameter set to properly
describe quantum confinement in nanocrystals. Applying this set we found a
quantitative influence of the interface structure on the band alignment.

###### keywords:

DFTB, TiO2, nanocrystals, Interface, Band alignment

llnl]Physical and Life Sciences Directorate, Lawrence Livermore National
Laboratory, Livermore, California 94550, United States llnl]Physical and Life
Sciences Directorate, Lawrence Livermore National Laboratory, Livermore,
California 94550, United States llnl]Physical and Life Sciences Directorate,
Lawrence Livermore National Laboratory, Livermore, California 94550, United
States Department of Chemical Engineering, University of California, Davis,
California 95616, United States Bremen Center for Computational Materials
Science, University of Bremen, P.O. Box 330440, D-28334 Bremen, Germany
DFT,DFTB

##  1 Introduction

TiO2 powders commonly used for photo catalytic applications such as the
Degussa P25 (DP25) consist of anatase and rutile nanocrystals in different
compositions (commonly $\sim 80$% anatase and $\sim 20$% rutile). This
composition of the two different TiO2 phases is shown to have a great impact
on the photo catalytic efficiency due to band alignment effects that lead to
an effective electron-hole pair separation. There have been several
theoretical attempts to predict the type of band alignment and thus which
crystal acts as a hole or electron trap1, 2, 3, 4. Experiments indicate that
for the case of DP25 anatase acts as an electron trap, while rutile as a hole
trap5. Theoretical works have so far been focused on isolated bulk1, slab 2 or
nanocrystal3, 4 systems. While bulk studies predict a trapping of electrons on
anatase and a trapping of the holes on rutile in accordance with the
experiment, nanocrystal studies predict a significant size dependence for the
alignment until bulk-like behavior is reached. However, neither of those
approaches consider the potential impact of different types of additional
degrees of freedom in materials such as geometrical interface arrangements.
The role of such interface arrangements were investigated in Ref. 2 using slab
models. However, such models can only give limited insights in the interface
effects as the slabs of the two phases must be commensurate, allowing only for
a few possible orientations and surface geometries. This makes it impossible
to explore the degrees of freedom in the interface formation, which are
present under experimental conditions at the microscale.

The vast phase space of possible nanocrystal orientations and sizes thus makes
computational studies of anatase-rutile nanocrystal interfaces very
challenging, even for very small nanocrystal models where a conclusive study
on the role of the interfaces would be unfeasible using DFT based calculations
alone. Consequently, we have developed a density functional tight binding
(DFTB) parameter set to efficiently simulate TiO2 nanocrystals and anatase-
rutile interfaces. DFTB, being 2-3 orders of magnitude faster than comparable
DFT calculations, provides the possibility of studying a vast number of
systems and system sizes. Additionally, with the right choice of parameters,
DFTB retains most of the accuracy of DFT and thus yields the possibility of
making a closer to one-to-one comparison to experiments. Another benefit of
the DFTB approach is its localized basis set which simplifies the computation
of non-periodic structures such as nanocrystals.

In the present work we identify several issues in describing nanocrystals with
the previously published DFTB parameter sets for titanium dioxide 6, 7, e.g.
problems in describing the quantum confinement effects correctly and yielding
wrong relaxations when under-coordinated species are involved. The latter
problem has been also reported in Reference 8. To remedy those issues, we have
developed a new parameter set, named tio2nano 9, that focuses on the correct
description of TiO2 nanocrystals and interfaces. It is based on the 3ob-
parameter set 10. The electronic parameters of the new Ti species were
tailored to ensure the proper description of the quantum confinement effects
in TiO2 nanocrystals. Instead of the traditional two-center potentials, we
have used a many-body force field 11, 12 to represent the repulsive
contributions to the total energy, following the same approach as in Ref. 13.
The three-center contributions of the force field allowed us to obtain
improved agreement with DFT optimized geometries, compared to results from
two-center repulsive potentials only.

For the band alignment studies, we considered explicit anatase-rutile
interface structures. We are not aware of any previous attempts to determine
the electronic structure using such models, but similar CeO2 nanoparticle
models were addressed in Ref. 14. By applying the newly developed DFTB
parameter set on the generated interface structures, we were able to
demonstrate the impact of the geometric alignment of the crystals on the band
alignment.

##  2 Methods

DFTB calculations were performed with the DFTB+ code 15. As described below in
detail, the 3ob parameter set 10 has been extended with the element Ti by
using a density compression radius of 8.5 Bohr and wave compression radii of
5.609, 3.958 and 7.0 Bohr for the s-, p- and d- orbitals of the Ti-atom,
respectively. The energy and the chemical hardness (a.k.a. the Hubbard-U
value) of the virtual 4p orbital of the Ti-atom have been set to 0.206 Ha and
-0.08 Ha, respectively.

The energies were created using the Chebyshev Interaction Model for Efficient
Simulation (ChIMES)11, 12, a reactive many-body molecular dynamics (MD) force
field. ChIMES creates many-body interactions by projecting DFT computed data
(e.g., forces, stress tensors, energies) onto linear combinations of many-body
Chebyshev polynomials16. Briefly, this begins with an N-body expansion of the
total energy for a system:

$E_{\mathrm{rep}}=E_{\mathrm{ChIMES}}=\sum_{i}^{N}E_{i}+{\sum_{i}^{N}}\sum_{j>i}^{N}E_{ij}+\sum_{i}^{N}\sum_{j>i}^{N}\sum_{k>j}^{N}E_{ijk}+\mathcal{O}(n).$
(1)

Here, $E_{i}$ is the single atom energy for a given element, $E_{ij}$ and
$E_{ijk}$ represent the 2–body and 3–body interaction energy, respectively,
$N$ is the total number of atoms in the system, and $\mathcal{O}(n)$
corresponds to higher order terms.

Specifically, the 2–body interactions are expressed as a linear combination of
Chebyshev polynomials of the first kind:

$E_{ij}=f_{\mathrm{P}}^{ij}(r_{ij})+f_{\mathrm{C}}^{ij}(r_{ij})\sum_{n=1}^{\Theta_{2}}\mathbbm{c}^{e_{i}e_{j}}_{n}T_{n}(s^{e_{i}e_{j}}_{ij})$
(2)

In this case, $T_{n}(s_{ij})$ is the Chebyshev polynomial of the first kind of
$n^{\mathrm{th}}$ order, $e_{i}$ and $e_{j}$ are the element types of atoms
$i$ and $j$, and $s_{ij}$ is a transformation of the inter-atomic distance
$r_{ij}$ over the Chebyshev interval of [-1,1] using a Morse-like function. In
addition, $\Theta_{2}$ corresponds to the two-body polynomial order,
$f_{\mathrm{C}}^{ij}(r_{ij})$ is the cut-off function that ensures the
potential and its derivative vary smoothly to zero beyond a specified
distance, and $f_{\mathrm{P}}^{ij}(r_{ij})$ is a penalty function 13, 17 that
helps to prevent sampling of inter-atomic distances below those seen in the
training set (see Ref. 18 for further details). The
$\mathbbm{c}^{e_{i}e_{j}}_{n}$ are the set of permutationally invariant
coefficients of linear combination for a given atom pair-type that are
determined via a linear least-squares method.

Similarly, 3–body interaction energies are expressed as a product of Chebyshev
polynomials for the constituent atom pairs of a given triplet, yielding an
orthogonal 3–body polynomial with $\binom{3}{2}=3$ total interactions for a
given triplet:

$E_{ijk}=f_{\mathrm{C}}^{ij}(r_{ij})f_{\mathrm{C}}^{ij}(r_{ik})f_{\mathrm{C}}^{ij}(r_{jk})\sum_{m=0}^{\Theta_{3}}\sum_{p=0}^{\Theta_{3}}\sum_{q=0}^{\Theta_{3}}{{}^{\prime}}\mathbbm{c}^{e_{i}e_{j}e_{k}}_{mpq}T_{m}(s^{e_{i}e_{j}}_{ij})T_{p}(s^{e_{i}e_{k}}_{ik})T^{e_{j}e_{k}}_{q}(s^{e_{j}e_{k}}_{jk}).$
(3)

In this case, $\Theta_{3}$ is used to label the three-body polynomial order
with a single permutationally invariant coefficient
$\mathbbm{c}^{e_{i}e_{j}e_{k}}_{mpq}$ for each set of triplet atom types. To
guarantee that only 3-body interactions between $i,j,k$ are counted towards
the sum, only terms for which at least two of the three $m,p,q$ indices are
greater than zero are included in the sum (indicated by the prime in Eq. 3).
In this way, many-body interactions can be included in the repulsive energy
($E_{\mathrm{rep}}$) by solving a linear least-squares optimization problem
where optimal coefficients of the linear combination are determined directly.
This avoids reliance on iterative approaches that are required for non-linear
optimization problems (e.g., Levenberg-Marquardt) that are usually more
computationally time consuming and not guaranteed to result in the global
minimum. ChIMES models have been extended to include 4–body interactions in a
similar fashion in reactive MD simulations16, 19, though truncation of the
ChIMES total energy with the 3–body term has proven sufficient for
determination of the DFTB repulsive energy.

Our training set for the ChIMES force field was determined from DFT-MD
simulations of amorphous TiO2 run at temperatures of 2250 and 300 K using the
VASP code 20. For these calculations, we used the PBEsol functional21 with an
energy cutoff of 550 eV and $\Gamma$-point sampling of the Brillouin zone.
PBEsol was chosen due to its improved description of solids and their
surfaces. The amorphous phase allowed for improved sampling of a wide-range of
interatomic distances over the short-time scales of the simulations, which was
enhanced by including data from elevated temperature. Each MD simulation was
run on a system of 216 atoms for a total of 5 ps, with configurations taken
for $E_{\mathrm{rep}}$ training every 100 fs (in order to allow for decoupling
between training data), resulting in 50 training configurations taken from
each temperature. We also included additional 10 configurations from an MD
simulation run at 40 GPa in order to improve sampling of close interatomic
distances. This yielded a total of 110 training configurations.

The training set to determine $E_{\mathrm{rep}}$ was then computed by
subtracting the gradient of the DFTB electronic energy from our DFT reference
property, i.e., $\vec{F}_{\mathrm{TRAIN}}=\vec{F}_{DFT}-\vec{F}_{\mathrm{DFTB-
elec}}$. In this work, a spline fit of Ti-Ti dimer repulsive interactions was
included in the DFTB electronic calculations order to ensure that excessively
small inter-atomic distances were not approached during calculations. In
addition, O-O distances were poorly sampled in our training set and these
repulsive parameters were thus taken from the 3ob-0-1 parameter set and were
not a part of our fit. ChIMES parameters were determined through linear least-
squares fitting to the resulting ionic forces and the diagonal components of
the stress tensor for each of these MD snapshots. We set the ChIMES 2–body
polynomial order to 12 and the 3–body order to 8, similar to previous work22,
and solve for optimal coefficients using the Least-Angle Regression (LARS)23,
24 algorithm with a Least Absolute Shrinkage and Selection Operator (LASSO)25
regularization value of $10^{-4}$.

For the DFTB calculations, charges were converged with a tolerance of
$10^{-6}$ a.u. and forces with $10^{-4}$ a.u. Bulk calculations were performed
with an 8$\times$8$\times$8 Monkhorst-Pack 26 grid for the primitive unit-
cell. The rutile (110) surface was constructed as a 4$\times$2 surface unit
cell with 6 layers, the (100) surface as a 4$\times$2 unit cell with 6 layers
and the (001) surface as $4\times 4$ with 8 layers. The (001) anatase surface
had a 2$\times$2 unit cell with 8 layers, the (101) surface had a 8$\times$4
unit cell with 5 layers. All rutile structures were computed with a
2$\times$2$\times$1 k-point mesh, the anatase (001) surface with a
4$\times$4$\times$1 mesh and the (101) surface with $\Gamma$-only
calculations. Branching point energy calculations (discussed below) were
performed with the primitive unit cell of the anatase and rutile crystals
using an 8$\times$8$\times$8 k-point grid. For nanocrystals we performed
$\Gamma$-only calculations.

The nanocrystal interfaces were constructed using the JANUS code , which is
based on a quick hull algorithm for convex hulls 27. This algorithm allows to
rapidly identify convex hulls from a set of points, in this case the
coordinates of atoms from the nanocrystals in xyz format. The code then
combines simplices of the computed convex hull for all normals that form an
$\arccos$ of $0.99$ or larger. The largest facets of each individual
nanocrystal are then aligned to each other so that their surface normals are
parallel to each other and to the x-axis with a user specified distance
between them. The procedure is repeated with the nanocrystals being rotated
with respect to each other around the x-axis for angles between 0 and 180
degrees in steps of 15 degrees. Additionally, the nanocrystal gets displaced
in the yz-plane along the perimeter of a circle with a user specified radius.
This displacement occurs with angles between 0 and 360 degrees in steps of 30
degrees. This results in 156 interface geometries which are stored in an ASE
28 trajectory format. To create the interfaces we used the two smallest
relaxed nanocrystals of anatase and rutile. The nanocrystals were aligned with
a distance along x of 2 Å and a radius of 2 Å was used for the displacement in
the yz-plane.

Interface energies were computed as

$E_{\mathrm{inter}}=E_{\mathrm{tot}}-\left(E_{\mathrm{rut}}+E_{\mathrm{ana}}\right)$
(4)

where $E_{\mathrm{tot}}$ is the energy of the interface structure,
$E_{\mathrm{rut}}$ the energy of the isolated rutile nanocrystal contained in
the interface structure and $E_{\mathrm{ana}}$ accordingly the energy of the
anatase nanocrystal.

For comparative DFT calculations we used the Vienna Ab initio Simulation
Package 20 (VASP 5.4.4) with the PAW method to describe the cores. For
titanium the 3p electrons were treated as valence electrons. The plane wave
cut-off was chosen to be 420 eV, the augmentation cut-off was chosen as 840
eV. Geometry relaxations were performed until energy differences were
converged below an error of $10^{-3}$ eV. The k-point meshes were chosen in
accordance with the DFTB calculations.

As nanocrystal models we chose Wulff-type crystal structures for anatase
nanocrystals which expose (101) facets and quasi Wulff-type crystals for
rutile exposing (110), (100) and (101) facets. These nanocrystal structures
were chosen in accordance with Ref. 3 to allow for direct comparison.
Additionally, for reasons described below, we also considered another set of
quasi Wulff-type rutile nanocrystals by cutting particles along the (110) and
(101) planes with an even number of (110) layers across the waist of the
particle.

##  3 Results and discussion

###  3.1 Parameterization and validation of the new DFTB parameter set

The existing parameter sets for Ti-O interactions, the tiorg-0-1 7 and the
matsci-0-3 6 sets, were evaluated only for bulk structures and common surfaces
so far. By applying them to describe anatase Wulff-type nanocrystals, we faced
the issue that the HOMO-LUMO gap did not change with the particle size as
predicted by the theory of quantum confinement3 and as observed in analogous
PBE calculations. Instead of the expected linear function of $n^{-2/3}$ (with
$n$ being the number of TiO2 units) we observed surface states, which fall
into the band gap and give rise to HOMO-LUMO gap sizes lower than the bulk
value as shown in Fig. 1. In order to rule out geometrical effects, we also
performed single point calculations for those nanocrystals using their PBE
geometries, but found a similar wrong behavior with both sets. Efforts to
resolve the issue by changing the compression radii of the Ti-atom turned out
to be unsuccessful.

Figure 1: Calculated HOMO-LUMO gaps of anatase nanocrystals using PBE and DFTB
with the tiorg-0-1, matsci-0-3 and the newly developed tio2nano set. The PBE
results have been shifted to align with the values obtained by the new
parameter set for the smallest nanocrystal. The tiorg-SP values were obtained
by carrying out DFTB calculations with the tiorg-0-1 set at the PBE geometries
(without relaxing with DFTB).

As a next attempt, we tried to describe the Ti-O interaction by extending the
3ob-parameter set 10 with Ti. The electronic parameters of the Ti atom were
optimized by fitting on the theoretically expected linear behaviour of the
band gap with respect to $n^{-2/3}$. For this fit, anatase nanocrystals with
sizes of $n=84$, $165$ and $286$ and the bulk phase of anatase had been used
resulting in the compression radii reported in the Methods section.

Figure 2: Relaxed geometries of the Ti35O70 nanocrystal. The DFTB geometry
obtained using a two-body repulsive potential for the Ti-O interaction is
shown by balls and sticks, while the PBE geometry is represented by blue solid
lines. The geometries had been aligned at the top oxygen atom.

After fixing the electronic parameters, the two-center repulsive potentials of
the Ti-Ti and the Ti-O interactions had been derived using a symmetrized
titanum hcp crystal and a symmetrized rutile bulk crystal as reference,
respectively. The symmetrization ensures that the first neighbour Ti-Ti (in
hcp Ti) and Ti-O bonds (in rutile) have identical lengths, allowing for a
simple manual fitting procedure. When relaxing the nanocrystals with the
obtained repulsive potentials, the resulting geometries exhibited significant
deviations in the strongly undercoordinated tip regions as compared to PBE-
geometries (see Fig. 2 for an example). Unfortunately, those geometrical
changes lead to the appearance of spurious surface states in the gap. Our
attempts to enforce the correct geometry by tuning the two-center repulsive
functions were not successful. This is not surprising as the incorrect Ti-O-Ti
angle seems to be the main driving force behind the differences in the relaxed
geometries. To enforce correct nanocrystal geometries (and with that also the
correct electronic structure), we decided to represent the Ti-O and Ti-Ti
repulsive functions with the ChIMES force field, which allows the inclusion of
three center terms. The interface to the ChIMES software had been implemented
in a development version of the DFTB+ package. For the Ti-Ti repulsive
potential we kept the two-center one derived from the symmetrized titanium hcp
crystal, with ChIMES corrective terms. It is important to note, that the
performance of the Ti-Ti repulsive potential in describing Ti bulk phases was
not evaluated in this work.

The resulting parameterization performs well for anatase and rutile bulk as
well as for their surfaces and is able to describe systems for which the
previous parameter sets have failed. The bulk lattice parameters given in
Table 1 show a very good agreement with ab initio counterparts for rutile. For
anatase we observe a slight overestimation of the $a$ value, while $c$ is
underestimated.

| DFTBa | PBEa | PBEsola | BLYPb | Exp.c
---|---|---|---|---|---
rutile |  |  |  |  |
$a$ | 4.629 | 4.621 | 4.600 | 4.679 | 4.594
$c$ | 2.980 | 2.954 | 2.940 | 2.985 | 2.959
anatase |  |  |  |  |
$a$ | 3.887 | 3.792 | 3.780 | 3.828 | 3.784
$c$ | 9.293 | 9.640 | 9.580 | 9.781 | 9.515
Table 1: Calculated and experimental anatase and rutile bulk lattice
parameters (in Å). The superscripts $a$, $b$ and $c$ denote values obtained in
this work, in Ref. 29 and in Ref. 30, respectively.

The calculated surface energies are listed in Table 2. The current
parameterization slightly overestimates the surface energies for anatase (001)
and rutile (110) compared to PBE and PBEsol results. The anatase (101) and the
rutile (001) and (100) surfaces are underestimated compared to PBEsol results
and overestimated compared to PBE. The energetic ordering of the surfaces with
respect to each other is reproduced nicely within the tio2nano set. Overall
the current set shows a good agreement to ab initio methods.

| DFTBa | PBEa | PBEsola | PBE0b
---|---|---|---|---
anatase |  |  |  |
(001) | 1.36 | 0.95 | 1.23 | 1.25*
(101) | 0.49 | 0.45 | 0.80 | 0.58
rutile |  |  |  |
(100) | 1.07 | 0.62 | 1.11 | 0.83
(110) | 0.97 | 0.37 | 0.88 | 0.46
(001) | 1.47 | 1.33 | 1.79 | 1.59
Table 2: Calculated anatase and rutile surface energies in J/m2 for selected
surfaces. The index $a$ denotes values obtained in this work, while values
marked with $b$ were taken from Ref. 31, except the one indicated by *, which
was taken from Ref. 32.

Figure 3 shows the PBE gap values of anatase and rutile in comparison to the
ones obtained with the newly developed parameter set. The behavior of both
sets is now comparable to Ref. 3. The smallest anatase nanocrystal considered
shows with the ChIMES based repulsive function a relaxation which is very
close to the one obtained with the PBE calculation (see Fig. 4). Overall, the
new parameters reasonably describe the various forms and phases of titania.
The band alignments between the phases are described in the following section.

Figure 3: Gap energy behavior over particle size for anatase (left) and rutile
(right). The DFTB values were obtained with the new parameter set and are
displayed in addition to the PBE results. The PBE results were aligned with
the DFTB ones to match the value for the smallest nanocrystal structure. Note
that for anatase the smallest crystal size was not included in the linear fit
as that structure did not follow the quantum confinement behaviour for the gap
states. The straight solid lines represent a linear fit, while the dashed
lines serve only as guides to the eyes. Figure 4: PBE Ti35O70 relaxed
nanocrystal displayed by blue lines aligned at the top oxygen atom to the DFTB
relaxation with the new parameter set (shown by sticks and balls). Relaxation
at the tip is so close to the PBE counter part that the blue lines are almost
not visible in the representation.

###  3.2 Band alignment

The primary aim of the current paper is to establish a computationally
feasible approach that allows us to study band alignment across a TiO2
anatase-rutile nano-particle interface and to determine the role of the
interface in this alignment. The approach is based on the customized DFTB
parametrization described above. The quality of the new method is assessed by
comparison to calculations at the DFT-PBE level of theory.

First we present the bulk band alignment between anatase and rutile using the
branching point energy approach as it was done for TiO2 in Ref. 1 previously
and compare the PBE and the DFTB values. Afterwards, following Ref. 3, the
level alignment deduced from non-bonded nano-particle structures will be
discussed. We compare previous results to our DFT-PBE calculations as well as
to results obtained with DFTB. Finally, we investigate the influence of the
interface on the band alignments by DFTB-calculations on bonded nano-
particles.

####  3.2.1 Bulk band alignment

First we consider the band alignment between the bulk phases using the
branching point technique. The branching point energy $E_{\mathrm{BP}}$, or
charge neutrality level, can be computed from the average of the mid-level
states using the relation

$\displaystyle
E_{\mathrm{BP}}=\frac{1}{2N_{\mathbf{k}}}\sum_{\mathbf{k}}\left[\frac{1}{N_{\mathrm{CB}}}\sum_{i}^{N_{\mathrm{CB}}}\epsilon_{\mathrm{CB}}^{\mathbf{k},i}+\frac{1}{N_{\mathrm{VB}}}\sum_{i}^{N_{\mathrm{VB}}}\epsilon_{\mathrm{VB}}^{\mathbf{k},i}\right]\text{.}$
(5)

The index $\mathbf{k}$ runs over the $N_{\mathbf{k}}$ k-points being
considered. In every k-point $N_{\mathrm{VB}}$ valence band and
$N_{\mathrm{CB}}$ conduction band states with respective eigenvalues
$\epsilon_{\mathrm{CB}}^{\mathbf{k},i}$ and
$\epsilon_{\mathrm{VB}}^{\mathbf{k},i}$ are averaged. We chose
$N_{\mathrm{VB}}=1$ and $N_{\mathrm{CB}}=1$ independently of possible
degeneracies of the bands, as the alignment is known to depend only weakly on
the degeneracy and on the number of included bands 1. Figure 5 compares our
results for EBP to those in Ref. 1. The results obtained with the tio2nano set
compare well to the calculated PBE values as well as to the HSE06 values in
Ref. 1. All levels of theory predict that rutile and anatase act as a hole and
as an electron trap, respectively.

Figure 5: Anatase and rutile band alignments computed by the branching point
technique with DFTB (tio2nano), PBE and HSE06. The HSE06 values had been taken
from Ref. 1.

####  3.2.2 Band alignment in non-bonded nanocrystal systems

Following Ref. 3 we performed DFT-PBE and DFTB calculations for anatase and
rutile nanocrystals of different sizes. We have plotted the HOMO and LUMO
levels as well as their gap values as function of the number of TiO2 formula
units in the structure raised to the factor ${-\alpha/3}$ . The parameter
$\alpha$ had been set to 1.5 as this yielded a better fit to our data than the
value of 1.35 used in Ref. 3. Note, that we did not include the smallest
anatase crystal in the linear fit, as we found it to not follow the particle
in the box model (neither for PBE nor for DFTB). The PBE and DFTB results for
anatase are in good agreement with the PBEx data of Ref. 3. By extrapolating
the DFTB nanocrystal results to infinite particle sizes using a linear fit
with respect to $n^{-1.5/3}$, we obtained a band gap of $3.02$ eV as compared
to the bulk value of $2.84$ eV. A similar overestimation can be seen in the
extrapolation of the PBE data, which yields a gap of $2.25$ eV as compared to
the PBE bulk value of $2.09$ eV.

In the case of rutile nanocrystals, extrapolation of the band gap is
complicated by the presence of surface states. This leads to a non-linearity
of the band gap caused by a non-linear behavior of the VB-edge. As shown in
Fig. 6, although the highest lying occupied core-states clearly follow the
confinement linearly, an occupied surface state with non-linear behavior
appears in the gap. This leads to an initial narrowing of the gap with the
particle size before reaching the point where the surface and core states
cross and after which the gap follows a linear behaviour again. A similar
behavior can be found in the PBE-results as shown in Fig. 7. The observed
behavior of the surface states indicates that the simple text-book particle in
a box model can be applied only to core-states but not to surface states.

Figure 6: HOMO and LUMO levels for anatase (left) and rutile (right)
calculated with DFTB. For rutile both, the highest occupied core states (HOMO-
core) and the occupied surfaces states (HOMO-surf) are shown. The
corresponding levels for the second rutile crystal set are indicated with
(S2). The dashed connection lines serve only as guides to the eye. The
indicated gap sizes of the bulk phases are obtained from the linear fits of
the respective HOMO and LUMO states.

Figure 7: HOMO and LUMO levels for anatase (left) and rutile (right)
calculated with PBE. For rutile both, the highest occupied core states (HOMO-
core) and the occupied surfaces states (HOMO-surf) are shown. The dashed
connection lines serve only as guides to the eye. The indicated gap sizes of
the bulk phases are obtained from the linear fits of the respective HOMO and
LUMO states.

Comparing the extrapolated value of the band gap for rutile, excluding the
surface states, from the nanocrystals with the calculated band gap of the
bulk, we find that the extrapolated value slightly overestimates the bulk one
with DFTB (2.88 eV vs. 2.76 eV) and underestimates it with PBE (1.62 eV vs
1.79 eV). We also note that the rutile data indicate an inverse quantum
confinement trend even when excluding the surface states, in disagreement with
the PBEx data presented in Ref. 3. We are not entirely sure what the origin of
the discrepancy is. We note however, that both rutile slabs 33 as well as
nanowires 34 exclusively exposing (110) surfaces show an even-odd behaviour in
the valence and conduction band edges with number of layers. For odd numbers,
both display an inverse quantum confinement behaviour, while for even numbers
they display regular quantum confinement behaviour. The effect was found to be
more pronounced for the conduction band states.

Therefore, in addition, we also considered another set of quasi Wulff-type
rutile nanocrystals by cutting particles along the (110) and (101) planes with
an even number of (110) layers across the waist of the particle. In the
following we will refer to this new set as SET 2. These nanocrystals contain
76, 92, 100, 108, 116, 132 and 316 TiO2 formula units. In order to make the
crystals stoichiometric, we removed a number of Ti atoms at the (101) facets
as opposed to adding dangling oxygen ions as was done in Ref. 3. The particles
of this set are found to be more stable compared to those presented in Ref. 3,
see Fig. 8. We also tried building larger particles with even numbers of (110)
layers, but these particles became metallic and also became significantly less
stable. Using the stable and non-metallic particles of the new set we obtain
an extrapolated value (excluding surface states) for the band gap of 2.89 eV,
see Fig. 6. We note that the new set displays a regular quantum confinement
trend with a decreasing HOMO-core to LUMO gap with increasing crystal size.
While the HOMO-core states of the rutile nanocrystals of the new set are
located in the same energy region as those of the original set, the LUMO
states are located at higher energies. This behaviour is therefore consistent
with the even-odd behaviour with respect to the number of (110) layers for the
VB and CB positions for slabs and nanowires in Refs. 33 and 34, respectively.

Figure 8: Comparison of stability of rutile nanocrystal in the orginal set
from Ref. 3 (SET 1) and the new set (SET 2) as a function of size. Energies
are given per TiO2 formula units.

Given that particles of comparable stability may show contrasting behaviour in
terms of their HOMO and LUMO positions suggests that we need to be careful in
establishing scaling relations – explicit simulations are clearly required for
small rutile nanocrystals.

Using the HOMO and LUMO level positions obtained above, band level alignments
in combined nanocrystals, consisting of a rutile and an anatase nanocrystal,
can be predicted. Note, however, that such predictions would not take the
effect of the actual anatase-rutile interfaces into account which can have a
significant impact as will become clear from the presentation in the section
on Band alignment in bonded nanocrystals. We also note that much of such size
dependence in the alignment will be driven by the large variation in the HOMO
and LUMO position in the rutile nanocrystals. The results will also depend on
the type of rutile nanocrystals we use. Nevertheless, combining an anatase
nanocrystal with rutile nanocrystals of type presented in Ref. 3, one obtains
a so-called type II rutile alignment (using the nomenclature of Ref. 3) for
small rutile nanocrystals, where rutile acts as an electron trap and anatase
as a hole trap. Increasing the size of the rutile nanocrystal results in a
transition to a so-called type I alignment, where both the HOMO and the LUMO
states are located within the rutile nanocrystal. Finally, for even larger
rutile crystal sizes one obtains a so-called type II anatase alignment (as
also predicted by the branching point calculation with the bulk structures),
where rutile acts as a hole and anatase as an electron trap. The same
quantitative behaviour was derived for PBEx in Ref. 3.

On the other hand, combining an anatase nanocrystal with a rutile nanocrystal
from the new set (S2) one can obtain a number of different alignments. Most
notably, using the smallest rutile nanocrystal of the new set leads to an
alignment where both the HOMO and the LUMO states are located within the
anatase nanocrystal.

With the aforementioned difficulties in establishing robust trends in the HOMO
and LUMO levels in rutile nanocrystals, we conclude that we are unable to make
robust prediction in HOMO-LUMO level alignment between anatase and rutile
nanocrystals, at least for size range covered by our simulations, $\sim$ 1-4.5
nm, without resorting to explicit simulations.

####  3.2.3 Band alignment in bonded nanocrystals

Exploiting the efficiency of the DFTB method, we have investigated the
influence of the interface on the gap alignment by simulating explicit
interface structures between different orientations of the smallest anatase
and rutile nanocrystals. In the most stable structures we have obtained for
nanocrystals with a rutile (110) and anatase (101) interface, the nanocrystals
tend to align with an angle of around 30 degrees with respect to each other
along the rutile [001] and anatase [010] directions. This alignment, as shown
in Figure 9, maximizes the interaction area and the bond formation between the
two structures, minimizing also the interface energy which was computed
according to Equation 4.

Figure 9: Typical alignment of the most stable rutile-anatase nanocrystal
interface structures. The red and blue structures represent the rutile and
anatase nanocrystals, respectively.

The prediction using independent nanocrystals (neglecting the interface
effects) suggests a type I alignment with a mismatch of 0.07 eV for the HOMO
and 1.03 eV for the LUMO. This suggests a complete domination of the gap by
rutile states, due to the narrowing of its gap by surface states. The gap
alignments for all 156 investigated bonded arrangements, their average and the
value obtained from the non-interacting nanocrystals are displayed in Figure
10. The 5 configurations with the lowest energy has been also marked in the
figure. It is obvious that although the prediction using non-interacting
nanocrystals delivers the correct alignment type, the magnitude in the offset
are far off from both the average interface offset and the one obtained from
the most stable interface structures. This suggests, that interface effects
have a significant influence on the level alignment and must be accounted for
explicitly. One observes a slight correlation of the LUMO/HOMO alignments
along with the indicated diagonal line which would correspond to rigid shifts
of the bands edges without any change in band gap. This therefore represents
situations, where bands in both phases are shifted due to an interface dipole
that varies in magnitude depending on the relative orientation of the
nanocrystals, but where no new interface states are being created.

Figure 10: HOMO and LUMO alignment of TiO2 nanoparticles composed of rutile
and anatase nanocrystals. The alignment is positive if the respective edge
state in rutile is higher as in anatase. Blue dots indicate values obtained
from bonded nanocrystals, green crosses their average, and the red cross the
value obtained from the non-interacting nanocrystals. The five most stable
structures are marked with yellow triangles, but since two pairs have very
similar alignment values, only three triangles can be seen. The dashed
diagonal line serves a guide for the eye and moving parallel to this line
correspond to alignments with constant gap sizes for both nanocrystals. In the
line shown, the values for these gaps are 2.08 eV for rutile and 3.18 eV for
anatase and corespond to the HOMO-LUMO gaps of the isolated particles.

Some typical LUMO and HOMO wave functions of the interfaces are displayed in
Fig. 11. We found three different types of LUMO wave functions and one HOMO
wave function. The HOMO wave function is identical to the rutile HOMO-surface
state. The most common LUMO wave function is located exclusively in the rutile
nanoparticle, but differs considerably from the LUMO wave function of the
isolated rutile nanocrystal. The second most common is located directly at the
interface, while the least common one, which only occurs in four of the
investigated cases, resembles the original rutile LUMO wave function. The
interface in each of the five most stable structures is of the first type,
which is also the most common among all investigated nanoparticles. This
underlines the necessity to model nanocrystal interfaces explicitely, as the
LUMO edge states might differ considerably from the those obtained from the
independent particle model. This aspect would be important to consider not
least, when modeling catalytic reactions at the surface of nanocrystals.

Figure 11: Typical LUMO and HOMO wavefunction in the anatase-rutile
nanocrystals. Subfigure a) shows the least common LUMO wave function type
which is similar to the LUMO of an isolated rutile nanocrystal, b) shows the
most common interface wave function, where wave function is not present in the
interface region, c) shows the second most common interface type, where the
LUMO wave function is located at the interface. Subfigure d) displays the HOMO
wave function, which is similar in all interfaces and corresponds to the HOMO-
surface state of the isolated rutile nanoparticles.

##  4 Conclusions

We have developed a DFTB parameter set using the ChIMES force field with
three-body terms as repulsive potentials especially designed for rutile and
anatase nanocrystals. We found that the three-body term was crucial in order
to predict the correct relaxation around the low coordinated atoms at the apex
of anatase nanocrystals.

We used the new DFTB parametrization, the tio2nano set, to investigate
different anatase-rutile band alignment models in order to predict effective
charge carrier separation in mixed anatase-rutile systems. While the bulk band
alignment model predicts a type II anatase alignment where anatase acts as an
electron and rutile as a hole trap in accordance with the literature 1, 2, 5,
the nanocrystal model indicates a strong dependence of the band alignment type
on the crystal size in accordance to the predictions in Reference 3. The
detailed shape of small rutile nanocrystals also plays a crucial role. We
showed that two types of rutile nanocrystal which both predominantly expose
the (110) facets and which both have similar stability give rise to rather
different behaviours in terms of band alignment. This fact also makes it
difficult to establish robust rules of thumb when it comes to predicting band
alignment in a rutile/anatase nanocrystal mixture and underlines the necessity
for explicit simulations of those systems.

Using the efficiency of the DFTB method, we also investigated the effect of
the anatase/rutile nanocrystal interface in the band alignment. We optimized
the geometry of 156 anatase/rutile nanocrystal pairs where the mutual
orientations were systematically varied. While, all interface models show the
same type I alignment, the magnitude of the band-offset vary almost one
electronvolt with different orientations. This observation suggests that it is
important to consider interface effects in the band alignment and underpins
the importance of being able to extend the reach of electronic structure
simulations beyond the realms spanned by regular DFT.

The DFG grant RTG 2247 is acknowledged. JK would like to acknowledge the
Swedish National Strategic e-Science programme (eSSENCE). We would also like
to acknowledge Pavlin D. Mitev for his insightful comments. Prepared by LLNL
under Contract DE-AC52-07NA27344. Project 20-SI-004 with Brandon Wood as PI
was funded by the Laboratory Directed Research and Development Program at
LLNL.

## References

  * Deák et al. 2011 Deák, P.; Aradi, B.; Frauenheim, T. Band Lineup and Charge Carrier Separation in Mixed Rutile-Anatase Systems. _J. Phys. Chem. C_ 2011, _115_ , 3443––3446
  * Kullgren et al. 2014 Kullgren, J.; Huy, H. A.; Aradi, B.; Frauenheim, T.; Deák, P. Theoretical study of charge separation at the rutile–anatase interface. _Phys. Status Solidi Rapid Res. Lett._ 2014, _8_ , 566–570
  * Ko et al. 2017 Ko, K. C.; Bromley, S. T.; Lee, J. Y.; Illas, F. Size-Dependent Level Alignment between Rutile and Anatase TiO2 Nanoparticles: Implications for Photocatalysis. _J. Phys. Chem. Lett._ 2017, _8_ , 5593–5598
  * Lamiel-Garcia et al. 2017 Lamiel-Garcia, O.; Ko, K. C.; Lee, J. Y.; Bromley, S. T.; Illas, F. When Anatase Nanoparticles Become Bulklike: Properties of Realistic TiO2 Nanoparticles in the 1–6 nm Size Range from All Electron Relativistic Density Functional Theory Based Calculations. _J. Chem. Theory. Comput._ 2017, _13_ , 1785–1793
  * Hurum et al. 2003 Hurum, D. C.; Agrios, A. G.; Gray, K. A.; Rajh, T.; Thurnauer, M. C. Explaining the Enhanced Photocatalytic Activity of Degussa P25 Mixed-Phase TiO2 Using EPR. _J. Phys. Chem. B_ 2003, _107_ , 4545––4549
  * Luschtinetz et al. 2009 Luschtinetz, R.; Frenzel, J.; Milek, T.; Seifert, G. Adsorption of Phosphonic Acid at the TiO2 Anatase (101) and Rutile (110) Surfaces. _J. Phys. Chem. C_ 2009, _113_ , 5730–5740
  * Dolgonos et al. 2010 Dolgonos, G.; Aradi, B.; Moreira, N. H.; Frauenheim, T. An Improved Self-Consistent-Charge Density-Functional Tight-Binding (SCC-DFTB) Set of Parameters for Simulation of Bulk and Molecular Systems Involving Titanium. _J. Chem. Theory. Comput._ 2010, _6_ , 266–278
  * Balzaretti et al. 2021 Balzaretti, F.; Gupta, V.; Ciacchi, L. C.; Aradi, B.; Frauenheim, T.; öppen, S. Water reactions on reconstructed rutile TiO2: a DFT / DFTB approach. _in preparation_ 2021,
  * 9 The parameter set can be downloaded from http://dftb.org
  * Gaus et al. 2013 Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. _J. Chem. Theory. Comput._ 2013, _9_ , 338–354
  * Koziol et al. 2017 Koziol, L.; Fried, L. E.; Goldman, N. Using Force Matching To Determine Reactive Force Fields for Water under Extreme Thermodynamic Conditions. _J. Chem. Theory Comput._ 2017, _13_ , 135–146
  * Lindsey et al. 2017 Lindsey, R. K.; Fried, L. E.; Goldman, N. ChIMES: A Force Matched Potential with Explicit Three-Body Interactions for Molten Carbon. _J. Chem. Theory Comput._ 2017, _13_ , 6222–6229
  * Goldman et al. 2018 Goldman, N.; Aradi, B.; Lindsey, R. K.; Fried, L. E. Development of a Multicenter Density Functional Tight Binding Model for Plutonium Surface Hydriding. _J. Chem. Theory. Comput._ 2018, _14_ , 2652–2660
  * Kim et al. 2019 Kim, B.-H.; Kullgren, J.; Wolf, M. J.; Hermansson, K.; Broqvist, P. Multiscale Modeling of Agglomerated Ceria Nanoparticles: Interface Stability and Oxygen Vacancy Formation. _Frontiers in Chemistry_ 2019, _7_ , 203
  * Hourahine et al. 2020 Hourahine, B. et al. DFTB+, a software package for efficient approximate density functional theory based atomistic simulations. _J. Chem. Phys._ 2020, _152_ , 124101
  * Pham et al. 2021 Pham, C. H.; Lindsey, R. K.; Fried, L. E.; Goldman, N. Calculation of the detonation state of HN3 with quantum accuracy. _J. Chem. Phys._ 2021, _153_ , 224102
  * Armstrong et al. 2020 Armstrong, M. R.; Lindsey, R. K.; Goldman, N.; Nielsen, M. H.; Stavrou, E.; Fried, L. E.; Zaug, J. M.; Bastea, S. Ultrafast shock synthesis of nanocarbon from a liquid precursor. _Nature Communications_ 2020, _11_ , 353
  * Lindsey et al. 2020 Lindsey, R. K.; Goldman, N.; Fried, L. E.; Bastea, S. Many-body reactive force field development for carbon condensation in C/O systems under extreme conditions. _J. Chem. Phys._ 2020, _153_ , 054103
  * Lindsey et al. 2020 Lindsey, R. K.; Fried, L. E.; Goldman, N.; Bastea, S. Active learning for robust, high-complexityreactive atomistic simulations. _J. Chem. Phys._ 2020, _153_ , 134117
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_
  * Perdew et al. 2008 Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. _Phys. Rev. Lett._ 2008, _100_ , 136406
  * Dantanarayana et al. 2020 Dantanarayana, V.; Nematiaram, T.; Vong, D.; Anthony, J. E.; Troisi, A.; Cong, K. N.; Goldman, N.; Faller, R.; Moulé, A. J. Predictive Model of Charge Mobilities in Organic Semiconductor Small Molecules with Force-Matched Potentials. _J. Chem. Theory Comput._ 2020, _16_ , 3494–3503
  * Efron et al. 2004 Efron, B.; Hastie, T.; Johnstone, I.; Tibshirani, R. Least angle regression. _Ann. Stat._ 2004, _32_ , 407–499
  * Friedman et al. 2010 Friedman, J.; Hastie, T.; Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. _J. Stat. Softw._ 2010, _33_ , 1
  * Tibshirani 1996 Tibshirani, R. Regression shrinkage and selection via the lasso. _J. R. Statist. Soc. B_ 1996, _58_ , 267–288
  * Monkhorst and Pack 1976 Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. _Phys. Rev. B_ 1976, _13_ , 5188
  * Barber et al. 1996 Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. The quickhull algorithm for convex hulls. _ACM Transactions on Mathematical Software_ 1996, _22_ , 469
  * Larsen et al. 2017 Larsen, A. H. et al. The Atomic Simulation Environment—A Python library for working with atoms. _J. Phys.: Condens. Matter._ 2017, _29_ , 273002
  * Lazzeri et al. 2001 Lazzeri, M.; Vittadini, A.; Selloni, A. Structure and energetics of stoichiometric TiO2 anatase surfaces. _Phys. Rev. B_ 2001, _63_ , 155409
  * Abrahams and Bernstein 1971 Abrahams, S. C.; Bernstein, J. L. Rutile: Normal Probability Plot Analysis and Accurate Measurement of Crystal Structure. _J. Chem. Phys._ 1971, _55_ , 3206
  * Labat et al. 2008 Labat, F.; Baranek, P.; Adamo, C. Structural and Electronic Properties of Selected Rutile and Anatase TiO2 Surfaces: An ab Initio Investigation. _J. Chem. Theory. Comput._ 2008, _4_ , 341–352
  * Araujo-Lopez et al. 2016 Araujo-Lopez, E.; Varilla, L. A.; Seriani, N.; Montoya, J. A. TiO 2 anatase’s bulk and (001) surface, structural and electronic properties: A DFT study on the importance of Hubbard and van der Waals contributions. _Surface Science_ 2016, _653_ , 187–196
  * Bredow et al. 2004 Bredow, T.; Giordano, L.; Cinquini, F.; Pacchioni, G. Electronic properties of rutile $\mathrm{Ti}{\mathrm{O}}_{2}$ ultrathin films: Odd-even oscillations with the number of layers. _Phys. Rev. B_ 2004, _70_ , 035419
  * Hmiel and Xue 2012 Hmiel, A.; Xue, Y. Quantum confinement and surface relaxation effects in rutile TiO2 nanowires. _Phys. Rev. B_ 2012, _85_ , 235461