# Topological phonons in oxide perovskites controlled by light

Bo Peng,1 Yuchen Hu,1,2 Shuichi Murakami3,4,
Tiantian Zhang,3,4∗ Bartomeu Monserrat1,5†

1Cavendish Laboratory, University of Cambridge,
J. J. Thomson Avenue, Cambridge CB3 0HE, United Kingdom
2Department of Chemistry, University of Cambridge,
Lensfield Road,, Cambridge CB2 1EW, United Kingdom
3Department of Physics, Tokyo Institute of Technology,
Ookayama, Meguro-ku, Tokyo 152-8551, Japan
4Tokodai Institute for Element Strategy, Tokyo Institute of Technology,
Nagatsuta, Midori-ku, Yokohama, Kanagawa 226-8503, Japan
5Department of Materials Science and Metallurgy, University of Cambridge,
27 Charles Babbage Road, Cambridge CB3 0FS, United Kingdom

∗†To whom correspondence should be addressed;
E-mail: ∗<EMAIL_ADDRESS>†<EMAIL_ADDRESS>
> Perovskite oxides exhibit a rich variety of structural phases hosting
> different physical phenomena that generate multiple technological
> applications. We find that topological phonons – nodal rings, nodal lines,
> and Weyl points – are ubiquitous in oxide perovskites in terms of structures
> (tetragonal, orthorhombic, and rhombohedral), compounds (BaTiO3, PbTiO3, and
> SrTiO3), and external conditions (photoexcitation, strain, and temperature).
> In particular, in the tetragonal phase of these compounds all types of
> topological phonons can simultaneously emerge when stabilized by
> photoexcitation, whereas the tetragonal phase stabilized by thermal
> fluctuations only hosts a more limited set of topological phonon states.
> Additionally, we find that the photoexcited carrier concentration can be
> used to tune the topological phonon states and induce topological
> transitions even without associated structural phase changes. Overall, we
> propose oxide perovskites as a versatile platform in which to study
> topological phonons and their manipulation with light.

## Introduction

In the family of topological materials, semimetals have an important place in
the study of topological order because they exhibit topologically protected
surface states (?), an anomalous bulk transport phenomenon known as the
“quantum anomaly” (?, ?), diversified classifications like nodal rings and
lines (?, ?, ?) and Weyl and Dirac points (?, ?, ?, ?, ?), and they serve as a
platform for obtaining various other topological states, such as topological
(crystalline) insulators (?) and Chern insulators (?). Beyond electronic
structures, topological degeneracies in the spectra of other quasiparticles
such as excitons, magnons, and phonons, have drawn wide attention in the past
decade (?, ?, ?, ?, ?). One advantage of the latter group is that excitation
of these quasiparticles depends only on the energy of the external probe,
rather than being restricted to the bands near the Fermi level like in the
case of electrons. Among these quasiparticles, phonons are of particular
interest as basic emergent bosonic excitations associated with lattice
vibrations (?, ?, ?, ?, ?, ?), and contribute to many physical processes, such
as conventional superconductivity (?, ?) and the thermal Hall effect (?, ?, ?,
?). Thus, studies of topological phonons, and especially on the control of
topological phonons, are becoming an important field in condensed matter
physics.

Oxide perovskites have gained popularity as a material system in the past
decades due to hosting multiple competing phases including ferroelectric (?),
magnetic (?), or superconducting (?), and also exotic excitations such as
skyrmions (?, ?). Additionally, these competing phases are particularly
sensitive to external stimuli, including temperature (?), strain (?), pressure
(?), and composition (?). More recently, theoretical and experimental evidence
has demonstrated that light is yet another means by which it is possible to
control the crystal symmetry of perovskites (?, ?, ?, ?, ?). As topological
order is intimately related to symmetry, the versatility of the perovskite
family makes these structures promising platforms to also explore topological
properties. However, most oxide perovskites are large gap insulators (for
example BaTiO3, SrTiO3, and PbTiO3), and as a consequence they exhibit no
electronic topological states.

In this work, we show that the phonon spectrum of multiple noncentrosymmetric
perovskites can host three types of topological states: topological nodal
rings, nodal lines, and Weyl points; suggesting that topological phonons are
pervasive in different structural phases of oxide perovskites. Using the
tetragonal BaTiO3 phase as a prototype, we show that all these types of
topological phonon emerge simultaneously when this phase is stabilized by
photoexcitation. Additionally, topological order can be controlled with the
photoexcited carrier density, driving topological transitions without any
associated structural phase change, including the creation and annihilation of
Weyl phonons and switching between nodal-ring and nodal-line phonons. In
contrast, when the tetragonal phase of BaTiO3 is stabilized by thermal
fluctuations, it only exhibits a more limited number of topological states.
This is because the long-range Coulomb interaction leads to a large energy
splitting between the longitudinal and transverse optical phonons (LO-TO
splitting), reducing the possibility of topological degeneracies between these
bands. The photoexcitation route eliminates this problem because photoexcited
carriers can screen long-range interactions.

## Results

### Crystal structures and structural phase transitions in BaTiO3

As one of the most studied perovskite materials, BaTiO3 has a cubic ABO3-type
crystal structure at temperatures above $393$ K, with the A cation (Ba2+) at
the corners and the B cation (Ti4+) at the center of an octahedral cage of
oxygen atoms. Upon cooling, BaTiO3 undergoes an inversion symmetry-breaking
phase transition, giving rise to a tetragonal phase with a deformation along
the [001] direction (?). As shown in Fig. 1(a), the phase transition breaks
inversion symmetry by the polar displacements of Ti and O atoms along the $z$
direction, which changes the space group from $Pm\bar{3}m$ (No. 221) to $P4mm$
(No. 99). With further cooling, an orthorhombic $Amm2$ phase appears below
$278$ K and a rhombohedral $R3m$ phase follows below $183$ K. As a result of
these temperature-driven phase transitions, the room temperature tetragonal
phase exhibits imaginary harmonic phonon modes, as shown in Fig. 1(b). The
inclusion of anharmonic vibrations (?) stabilizes the tetragonal structure at
$300$ K, also shown in Fig. 1(b).

Figure 1: Crystal structure and phonon dispersion of tetragonal BaTiO3. (a)
Crystal structures of the high-temperature cubic phase (No. 221, $Pm\bar{3}m$)
and of the room-temperature tetragonal phase (No. 99, $P4mm$) of BaTiO3. (b)
Phonon dispersion of tetragonal BaTiO3 at the harmonic level and at $300$ K
including anharmonic vibrations. (c) Ground-state and excited-state potential
energy surfaces along the $\Gamma$ point phonon mode that is imaginary within
the harmonic approximation. (d) Phonon dispersion of tetragonal BaTiO3 at
$n=0.06$ $e$/f.u. (e) Bulk Brillouin zone at $n=0.06$ $e$/f.u. with two nodal
rings on the $q_{x}=0$ and $q_{y}=0$ mirror planes (orange circles, formed by
the 10th and 11th bands), one nodal line along the $\Gamma$-Z high-symmetry
line (brown line, formed by the 10th and 11th bands), one pair of triple Weyl
points at the intersection points between the nodal rings and nodal line
(violet cones, formed by the 10th, 11th and 12th bands), and four pairs of
Weyl points on the $q_{z}$ = $\pi$ plane (purple and yellow cones, formed by
the 10th and 11th bands).

In addition to temperature, multiple strategies have traditionally been
proposed to control structural phase transitions in perovskites, including
strain (?), pressure (?, ?), and composition (?). More recently, it has been
shown that photoexcitation provides an alternative route to stabilizing
multiple perovskite phases (?, ?, ?, ?, ?). Photoexcitation is relatively
economical, and can be easier to control compared to other strategies such as
heating to change temperature, material synthesis to change composition or
strain, or external pressure. Inspired by these recent discoveries, we find
that photoexcitation can also stabilize the imaginary phonon modes of the
tetragonal phase of BaTiO3: Fig. 1(c) shows that the double well potential
energy curve along the imaginary phonon mode of amplitude $Q$ at the $\Gamma$
point becomes a single well upon illumination, indicating the stabilization of
the crystal structure. The underlying physics is that changing the carrier
concentration $n$ induced by photoexcitation leads to changes in the potential
energy experienced by the ions, which in turn changes the interatomic force
constants and the phonon dispersion. For $0.04<n<0.10$ $e$/f.u., the $P4mm$
phase is not only dynamically stable, but it is also thermodynamically more
stable than the cubic $Pm\bar{3}m$ phase; whereas for $n>0.10$ $e$/f.u., the
$P4mm$ phase relaxes to the cubic phase as the latter becomes
thermodynamically more stable (?). The relative energy between the two phases
and their corresponding lattice constants under illumination are shown in the
Supplementary Materials. Figure 1(d) shows the phonon dispersion with a
photoexcited carrier density of $n=0.06$ $e$/f.u., confirming the dynamical
stabilization of the tetragonal phase.

### Topological phonons in BaTiO3

As a large gap insulator, BaTiO3 displays trivial topology in its electronic
band structure [Fig. 2(a)]. However, there are three types of topological
states in its phonon dispersion, which are nodal rings, nodal lines, and Weyl
points. Since the unit cell of tetragonal BaTiO3 has 5 atoms, there are 15
branches in the phonon spectrum. Band inversion between the 10th and 11th
bands (labelled by increasing energy) along the X-$\Gamma$ high-symmetry line
is protected by $M_{x}$ symmetry, which restricts the band crossing to a 1D
continuous ring/line on the mirror-invariant plane (?). Therefore a nodal ring
is formed on the $q_{x}=0$ plane. Because of the existence of an additional
$C_{4z}$ symmetry, there is another nodal ring located on the $q_{y}=0$ plane,
related to the one on the $q_{x}=0$ plane by $C_{4z}$ symmetry [orange circles
in Fig. 1(e)]. The combination of $C_{4z}$ and mirror symmetries also brings
out another type of topological phonon in BaTiO3: an endless nodal line along
the $\Gamma$-Z direction [brown line in Fig. 1(e)]. Intersection points
between these two nodal rings and the nodal line form two triple Weyl points
located along the Z’(0,0,-$\pi$)-$\Gamma$-Z(0,0,$\pi$) line. Another
topological feature of the band crossings between the 10th and 11th bands are
eight Weyl points on the $q_{z}=\pi$ plane [purple and yellow cones in Fig.
1(e)]. These eight Weyl points, located at generic momenta on the $q_{z}=\pi$
plane, are not protected by any crystalline symmetry due to the breaking of
the $M_{z}$ symmetry in the tetragonal lattice, so we need to apply more
advanced symmetry-based indicator theory to diagnose the topology.

To gain a better understanding of the topological nature of the eight Weyl
points present in photoexcited tetragonal BaTiO3, we explore the existence of
Weyl points from the perspective of topology as shown in Fig. 3(a), diagnosing
topological degeneracies by using symmetry-based indicators from subgroups (?,
?, ?). After obtaining the symmetry data at $\Gamma$, M, and X, we note that
the symmetry-indicator group for space group No. 99 is a trivial one, and the
symmetry data at high-symmetry momenta does not satisfy the compatibility
condition. Thus, we need to find a subgroup that has a nontrivial symmetry-
based indicator group and satisfies the compatibility condition (?, ?). Space
group No. 99 has five subgroups with the same size of unit cell, namely No. 75
($P4$), No. 25 ($Pmm2$), No. 6 ($Pm$), No. 3 ($P2$), and No. 1 ($P1$). Among
them, No. 3 is the maximal subgroup that has a nontrivial symmetry-based
indicator group $\mathbb{Z}_{2}$ and satisfies the compatibility condition.
Thus, subgroup No. 3 is used for further diagnosis. Space group No. 3 only has
a generator of $C_{2z}$ rotation symmetry, so the symmetry-based indicator
group $\mathbb{Z}_{2}$ corresponds to the $z_{2}$ Berry phase of the loop
enclosing half of the $q_{z}=0$
[X($\pi$,0,0)-M($\pi$,$\pi$,0)-M’(-$\pi$,$\pi$,0)-X’(-$\pi$,0,0)-X($\pi$,0,0)]
or $q_{z}=\pi$
[R($\pi$,0,$\pi$)-A($\pi$,$\pi$,$\pi$)-A’(-$\pi$,$\pi$,$\pi$)-R’(-$\pi$,0,$\pi$)-R($\pi$,0,$\pi$)]
plane. If $z_{2}=0$, there will be 0 $mod$ 4 Weyl points on the $q_{z}=0$ and
$q_{z}=\pi$ planes. Our calculations deliver eight Weyl points in total,
corresponding to the case $\mathbb{Z}_{2}=0$. Furthermore, considering the
additional mirror symmetries in space group No. 99, the eight Weyl points on
the $q_{z}=\pi$ plane are related by $M_{x}$ and $M_{y}$ symmetries: four of
them with left-hand chirality, while the other four with right-hand chirality.

Figure 2: Nodal lines and nodal rings in tetragonal BaTiO3 controlled by
light. (a) Schematic of the photoexcitation process in tetragonal BaTiO3.
(b)-(d) Nodal lines and nodal rings formed by the 10th and 11th phonon
branches in the bulk Brillouin zone for photoexcited carrier concentrations of
0.04 $e$/f.u., 0.07 $e$/f.u., and 0.12 $e$/f.u. The corresponding phonon
dispersions are shown in (e)-(g).

### Controlling topological phonons by light

As discussed above, both photoexcitation with $0.04<n<0.10$ $e$/f.u. and
temperatures between $278$-$393$ K can stabilize the tetragonal phase of
BaTiO3, as shown in Figs. 1(b) and (d). However, a main difference is that
without photoexcitation, the long-range Coulomb interaction leads to a large
LO-TO splitting (?), which tends to lift band degeneracies, as shown in Fig.
1(b). Consequently, the eight conventional Weyl points and the two triple Weyl
points that we find in illuminated BaTiO3 are not observed in temperature-
stabilized tetragonal BaTiO3 at 300 K (although the presence of nodal-line and
nodal-ring phonons is confirmed in the Supplementary Materials). In contrast,
under illumination, the photoexcited carriers lead to a strong free-carrier
screening that suppresses the LO-TO energy splitting and facilitates the
appearance of Weyl phonons in tetragonal BaTiO3. Another advantage of the
photoexcitation route is that, as discussed next, it can drive the phononic
system between different topological quantum states, including switching
between nodal-ring and nodal-line phonons, controlling the position of Weyl
points in momentum space, and creating and annihilating these topological
states.

As discussed above, changing the photoexcited carrier density brings about not
only the transition between nodal rings and nodal lines, but also the creation
of new nodal rings. As shown in Fig. 2(b)-(c), increasing the photoexcited
carrier density from $n=0.04$ $e$/f.u. to $n=0.07$ $e$/f.u. drives the two
nodal rings to become two nodal lines, and the corresponding phonon spectra
are shown in Fig. 2(e)-(f). Further increasing $n$ to $0.12$ $e$/f.u. (under
which conditions the tetragonal structure becomes cubic) leads to the
formation of two nodal rings around the R point on the $q_{x}=0$ and $q_{y}=0$
planes, as shown in Fig. 2(d) and (g). Thus, by changing $n$, nodal rings in
the phonon spectra can be created or transformed into nodal lines.

Figure 3: Weyl points in tetragonal BaTiO3 controlled by light. (a) Diagnosis
process for the eight Weyl points in BaTiO3 using symmetry-based indicator
theory for the lowest 10 bands. (b) Evolution of the Weyl points on the
$q_{z}$ = $\pi$ plane with increasing $n$ from $0.050$ to $0.085$ $e$/f.u. (c)
Phonon surface arcs for a phonon frequency of 11.42 THz at $n=0.07$ $e$/f.u.

In addition to nodal lines and nodal rings, Weyl phonons are also sensitive to
the photoexcited carrier concentration $n$. Fig. 3(b) shows the evolution of
eight pairs of Weyl points on the $q_{z}=\pi$ plane when increasing $n$ from
$0.050$ to $0.085$ $e$/f.u. (four extra pairs of Weyl points emerge when
$0.0695<n<0.0825$ $e$/f.u.). Below $0.050$ $e$/f.u., BaTiO3 exhibits no Weyl
points between the 10th and 11th branches in its phonon spectrum. At about
$n=0.050$ $e$/f.u., four pairs of Weyl points are created around the
$q_{x}=q_{y}$ and $q_{x}=-q_{y}$ planes [orange and maroon dots in Fig. 3(b)].
As $n$ increases, the Weyl points in each pair with different chirality move
away from each other and head to the $q_{x}=0$ and $q_{y}=0$ planes [yellow
and purple arrows in Fig. 3(b)]. Once $n$ reaches $0.0815$ $e$/f.u, Weyl
points meet on the $q_{x}=0$ and $q_{y}=0$ planes and annihilate in pairs
[light yellow and light magenta dots in Fig. 3(b)]. In addition, another four
pairs of Weyl points of opposite chirality are created around the
$q_{x}=q_{y}$ and $q_{x}=-q_{y}$ planes at $n=0.0695$ $e$/f.u., and annihilate
in pairs on the $q_{x}=\pi$ and $q_{y}=\pi$ planes when $n$ reaches $0.0825$
$e$/f.u. [blue and green dots in Fig. 3(b)]. Since no band inversion happens
at any high-symmetry points in this process, the Berry phase of the loop
R-A-A’-R’-R on the $q_{z}=\pi$ plane remains zero, which is consistent with 0
$mod$ 4 Weyl points on the $q_{z}=\pi$ plane with $\mathbb{Z}_{2}=0$ [Fig.
3(a)].

To have a better understanding of the Weyl phonons in tetragonal BaTiO3, we
calculate the surface local density of states (LDOS) from the imaginary part
of the surface Green’s function (?). The surface LDOS in Fig. 3(c) is
calculated along the (001) direction at $11.42$ THz and $n=0.07$ $e$/f.u. The
four pairs of Weyl points connect via surface arcs crossing the Brillouin zone
boundaries, so that each surface arc starts from one Weyl point with a
positive monopole charge and ends at another with a negative one. Overall, by
modulating $n$ we can control the creation and annihilation of Weyl points,
and the length of surface arcs in tetragonal BaTiO3.

### Topological phonons in PbTiO3

Multiple types of topological phonons can also be found in other perovskites
with the same $P4mm$ space group, and here we study PbTiO3 and SrTiO3 as
additional examples.

Figure 4: Topological phonons in tetragonal PbTiO3. (a) Phonon dispersion of
tetragonal $P4mm$ PbTiO3 at $n=0.025$ $e$/f.u. (b) Bulk Brillouin zone at
$n=0.025$ $e$/f.u. with four nodal rings on the $q_{x}=q_{y}$ and
$q_{x}=-q_{y}$ planes (orange circles), two nodal lines along the $\Gamma$-Z
and M-A high-symmetry lines (brown line), and eight Weyl points located on
$q_{z}=\pi$ plane (purple and yellow cones). All these three types of
topological phonons are composed by the 10th and 11th bands. (c) Phonon
surface arcs for a phonon frequency of 9.63 THz at $n=0.025$ $e$/f.u.

PbTiO3 has a tetragonal $P4mm$ structure for $T<600$ K, and a cubic structure
above that temperature. Athough the $P4mm$ phase has the lowest energy in a
wide photoexcited density range $n<0.125$ $e$/f.u. (?), it becomes dynamically
unstable when $n>0.05$ $e$/f.u. (shown in the Supplementary Materials). We
calculate the excited-state phonon dispersion at $n=0.025$ $e$/f.u., as shown
in Fig. 4(a). Similar to BaTiO3, PbTiO3 also has nodal rings, nodal lines, and
Weyl points between the 10th and 11th bands. But the triple Weyl points no
longer exist as the nodal rings and the nodal line do not touch. The unique
band crossings along the $\Gamma$-M high-symmetry line bring about four nodal
rings on the $q_{x}=q_{y}$ and $q_{x}=-q_{y}$ planes protected by $M_{xy}$ and
$M_{x\bar{y}}$ symmetries. The endless nodal lines along the $\Gamma$-Z and
M-A directions are robust as well, as they are protected by the $C_{4z}$
symmetry.

The surface LDOS in Fig. 4(c) is calculated along the (001) direction at
$9.63$ THz with $n=0.025$ $e$/f.u. Four pairs of Weyl points, created at
$n=0.005$ $e$/f.u., connect with each other via surface arcs crossing the
$q_{x}=q_{y}$ and $q_{x}=-q_{y}$ planes. We note that the location of these
four pairs of Weyl points can also be modulated by changing the
photoexcitation density until the lattice becomes dynamically unstable for
$n>0.05$ $e$/f.u.

### Topological phonons in SrTiO3

SrTiO3 has a different phase diagram compared to BaTiO3 and PbTiO3. At low
temperatures, it exhibits a centrosymmetric tetragonal phase of space group
$I4/mcm$, and undergoes a phase transition to a cubic phase above $106$ K. The
$P4mm$ phase can only be accessed with small negative strains (?, ?). The in-
plane lattice constants $a$ and $b$ of $P4mm$ SrTiO3 increase upon
photoexcitation. Therefore by fixing $a$ and $b$ to the ones in the dark, we
can induce in-plane negative strain in illuminated SrTiO3. This can be
realized experimentally by growing SrTiO3 on an appropriate substrate. As
shown in Fig. 5(a), after fixing the in-plane lattice constants $a$ and $b$
and relaxing only the out-of-plane lattice constant $c$ in the $P4mm$ phase,
we find that the tetragonal phase can have a lower energy than the cubic
$Pm\bar{3}m$ phase.

Figure 5: Topological phonons in tetragonal SrTiO3. (a) Relative energy
difference between the $P4mm$ and $Pm\bar{3}m$ phases with fixed in-plane
lattice constants $a$ and $b$, with the energy of the fully relaxed
$Pm\bar{3}m$ phase set to be zero. (b) Phonon dispersion of tetragonal $P4mm$
SrTiO3 at $n=0.09$ $e$/f.u. (c) Bulk Brillouin zone at $n=0.09$ $e$/f.u. with
two nodal lines on the $q_{x}=0$ and $q_{y}=0$ planes (orange lines), one
nodal line along the $\Gamma$-Z high-symmetry line (brown line), one pair of
triple Weyl points at the intersection points between the three nodal lines
(violet cones), and four pairs of Weyl points on the $q_{z}$ = $0$ plane
(purple and yellow dots).

$P4mm$ SrTiO3 exhibits imaginary phonon modes in the dark, while adding a
photoexcited charge density of $n=0.09$ $e$/f.u. stabilizes the structure, as
shown in Fig. 5(b). Similar to BaTiO3, the $P4mm$ phase of SrTiO3 has four
nodal lines in total on the $q_{x}=0$ and $q_{y}=0$ mirror planes and one
nodal line along the $\Gamma$-Z high-symmetry line in the bulk Brillouin zone.
The intersection points of the three nodal lines form one pair of triple Weyl
points. Different from BaTiO3, the four pairs of Weyl points are located on
the $q_{z}=0$ plane rather than the $q_{z}=\pi$ plane. Weyl phonons in SrTiO3
can also be manipulated by photoexcitation, and the critical $n$ for the
annihilation of the Weyl points on the $q_{x}=0$ and $q_{y}=0$ planes is 0.105
$e$/f.u., which is similar to that in BaTiO3.

## Ubiquity of topological phonons in perovskite oxides

The existence and coexistence of different types of topological phonons,
together with their tunability by light, occur in a variety of other
structural phases of perovskites beyond the tetragonal phases investigated
above. Topological phonons in the orthorhombic $Amm2$ and the rhombohedral
$R3m$ phases of BaTiO3 are shown in the Supplementary Materials. Furthermore,
topological phonons in oxide perovskites can also be manipulated by other
tuning parameters beyond photoexcitation, such as strain and temperature (for
details, see the Supplementary Materials). Putting these insights together
suggests that topological phonons are ubiquitous in the family of perovskite
oxides, which, given their widespread use and versatility, provide a promising
platform for exploring topological phonon physics and its interplay with other
phases.

Given the proposed ubiquity of topological phonons in perovskites, it should
be possible to directly identify them experimentally by looking for “isolated”
band crossings, and the prototype materials detailed above are promising
starting candidates. Topological bulk phonons can be measured by inelastic
neutron scattering (?, ?, ?) and inelastic X-ray scattering (?), while the
accompanying topological surface states can be detected by high resolution
electron energy loss spectroscopy (?). For light-induced topological phonons,
although it may be difficult to illuminate the sample to maintain a constant
photoexcited carrier density, transient photoexcitation could be used to
observe the evolution of nodal rings, nodal lines and Weyl points as a
function of the time-dependent photoexcited carrier density.

In terms of phenomena, topological phonons with non-zero Berry curvature may
contribute to the recently observed phonon thermal Hall effect in SrTiO3 (?,
?). In addition, novel phonon-phonon and electron-phonon scattering mechanisms
could arise from topological phonons, which may provide new insights on the
enhanced superconductivity of SrTiO3 (?, ?). The emergence of Weyl phonons may
be accompanied by other unique physical properties such as a pseudogauge field
with a one-way propagating bulk mode (?, ?, ?), topological negative
refraction (?), and nonlinear acoustic/optical responses (?, ?), which can
offer new routes for designing novel technologies like light-controlled
neuromorphic computing in phononic systems (?, ?). Our work shows that oxide
perovskites provide a promising platform to explore all of these phenomena and
applications.

## Conclusion

We find that topological phonons are ubiquitous in oxide perovskites and that
photoexcitation provides a promising route for their manipulation. As
examples, the noncentrosymmetric tetragonal phases of three oxide perovskites
(BaTiO3, PbTiO3, and SrTiO3) exhibit topological nodal rings, nodal lines, and
Weyl points in their phonon spectra. Remarkably, we find that photoexcitation
is the only way to obtain Weyl phonons in tetragonal BaTiO3, since the
thermally-stabilized tetragonal phase has a large LO-TO energy splittting in
the phonon spectrum that prevents band crossings. By contrast, photoexcited
carriers screen long-range interactions, suppressing the LO-TO energy
splitting and facilitating the band crossings. We also find that the
photoexcited carrier density can be used to tune the creation and annihilation
of Weyl points and nodal rings/lines without any associated structural phase
transitions. Topological phonons in oxide perovskites provide a promising
platform to study physical phenomena ranging from the phonon Hall effect to
superconductivity, and may also offer new technological opportunities such as
the realization of controllable topological quantum states for neuromorphic
computing.

## Methods

Density functional theory (DFT) calculations are performed using the Vienna ab
initio simulation package (vasp) with the projector-augmented-wave potential
method (?, ?). We use the generalized gradient approximation (GGA) with the
Perdew-Burke-Ernzerhof parameterization revised for solids (PBEsol) as the
exchange-correlation functional (?). A plane-wave basis set is employed with a
kinetic energy cutoff of $800$ eV and a $7\times 7\times 7$ k-mesh during
structural relaxation, which is stopped when forces are below 10-3 eV/Å. The
band structure of BaTiO3 is calculated using the HSE06 hybrid functional in
the presence of spin-orbit coupling (?). Photoexcited carriers are simulated
by promoting electrons from high-energy valence band states to low-energy
conduction band states. This $\Delta$ self-consistent field ($\Delta$SCF)
method introduces non-interacting electron-hole pairs by changing the
occupation numbers of the Kohn-Sham orbitals (?, ?, ?, ?), and is
computationally less demanding compared to other approaches like constrained
density functional theory (?) and excited-state force calculations (?).
Nevertheless, it gives consistent phonon spectra compared with those obtained
with constrained DFT (?), as shown in the Supplementary Materials. The
occupancies are fixed with a smearing of $0.01$ eV.

The crystal structures in the dark and under photoexcitation are fully relaxed
before the phonon calculations. We calculate the force constants with density
functional perturbation theory (DFPT) (?) in a $2\times 2\times 2$ supercell
with a $5\times 5\times 5$ k-mesh using vasp. The phonon dispersion is then
obtained using phonopy (?). We perform convergence tests on the supercell size
between $2\times 2\times 2$, $3\times 3\times 3$, and $4\times 4\times 2$, all
confirming the existence of Weyl points in illuminated BaTiO3. We calculate
the chirality of Weyl points by employing the Wilson-loop method to calculate
the Wannier charge center flow, $i.e.$, the monopole charge of a Weyl point
(?, ?). The phonon surface states are obtained using surface Green’s functions
as implemented in WannierTools (?). The finite-temperature phonon frequencies
including anharmonic contributions are calculated using a self-consistent ab
initio lattice dynamical method (?, ?, ?) in a $2\times 2\times 2$ supercell.
The self-consistent cycle is terminated after 240 iterations when the
difference in free energy is less than 0.2 meV, and the space group symmetry
is enforced on the resulting force constants. The Born effective charges are
calculated using DFPT to obtain the LO-TO splitting in the phonon dispersion
in the dark (?). Under illumination, the photoexcited electrons screen long-
range interactions and no LO-TO splitting occurs.

## References

  * 1. T. T. Heikkilä, N. B. Kopnin, G. E. Volovik, Flat bands in topological media, JETP Letters 94, 233– (2011).
  * 2. D. T. Son, B. Z. Spivak, Chiral anomaly and classical negative magnetoresistance of weyl metals, Phys. Rev. B 88, 104412 (2013).
  * 3. J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, N. P. Ong, Evidence for the chiral anomaly in the Dirac semimetal Na3Bi, Science 350, 413 (2015).
  * 4. H. Weng, Y. Liang, Q. Xu, R. Yu, Z. Fang, X. Dai, Y. Kawazoe, Topological node-line semimetal in three-dimensional graphene networks, Phys. Rev. B 92, 045108 (2015).
  * 5. A. Bernevig, H. Weng, Z. Fang, X. Dai, Recent Progress in the Study of Topological Semimetals, J. Phys. Soc. Jpn. 87, 041001 (2018).
  * 6. C. Fang, Y. Chen, H.-Y. Kee, L. Fu, Topological nodal line semimetals with and without spin-orbital coupling, Phys. Rev. B 92, 081201 (2015).
  * 7. Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, Z. Fang, Dirac semimetal and topological phase transitions in ${A}_{3}$Bi ($A=\text{Na}$, K, Rb), Phys. Rev. B 85, 195320 (2012).
  * 8. N. P. Armitage, E. J. Mele, A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * 9. S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, M. Z. Hasan, Discovery of a weyl fermion semimetal and topological fermi arcs, Science 349, 613–617 (2015).
  * 10. H. Weng, C. Fang, Z. Fang, B. A. Bernevig, X. Dai, Weyl semimetal phase in noncentrosymmetric transition-metal monophosphides, Phys. Rev. X 5, 011029 (2015).
  * 11. L. X. Yang, Z. K. Liu, Y. Sun, H. Peng, H. F. Yang, T. Zhang, B. Zhou, Y. Zhang, Y. F. Guo, M. Rahn, D. Prabhakaran, Z. Hussain, S.-K. Mo, C. Felser, B. Yan, Y. L. Chen, Weyl semimetal phase in the non-centrosymmetric compound taas, Nat Phys 11, 728–732 (2015).
  * 12. S.-Y. Yang, H. Yang, E. Derunova, S. S. P. Parkin, B. Yan, M. N. Ali, Symmetry demanded topological nodal-line materials, Advances in Physics: X 3, 1414631– (2018).
  * 13. J. He, X. Li, P. Lyu, P. Nachtigall, Near-room-temperature chern insulator and dirac spin-gapless semiconductor: nickel chloride monolayer, Nanoscale 9, 2246–2252 (2017).
  * 14. F. Wu, T. Lovorn, A. H. MacDonald, Topological exciton bands in moiré heterojunctions, Phys. Rev. Lett. 118, 147401 (2017).
  * 15. F.-Y. Li, Y.-D. Li, Y. B. Kim, L. Balents, Y. Yu, G. Chen, Weyl magnons in breathing pyrochlore antiferromagnets, Nature Communications 7, 12691– (2016).
  * 16. O. Stenull, C. L. Kane, T. C. Lubensky, Topological phonons and weyl lines in three dimensions, Phys. Rev. Lett. 117, 068001 (2016).
  * 17. Y. Liu, Y. Xu, S.-C. Zhang, W. Duan, Model for topological phononics and phonon diode, Phys. Rev. B 96, 064106 (2017).
  * 18. H. He, C. Qiu, L. Ye, X. Cai, X. Fan, M. Ke, F. Zhang, Z. Liu, Topological negative refraction of surface acoustic waves in a weyl phononic crystal, Nature 560, 61–64 (2018).
  * 19. T. Zhang, Z. Song, A. Alexandradinata, H. Weng, C. Fang, L. Lu, Z. Fang, Double-weyl phonons in transition-metal monosilicides, Phys. Rev. Lett. 120, 016401 (2018).
  * 20. H. Miao, T. T. Zhang, L. Wang, D. Meyers, A. H. Said, Y. L. Wang, Y. G. Shi, H. M. Weng, Z. Fang, M. P. M. Dean, Observation of double weyl phonons in parity-breaking fesi, Phys. Rev. Lett. 121, 035302 (2018).
  * 21. J. Li, Q. Xie, S. Ullah, R. Li, H. Ma, D. Li, Y. Li, X.-Q. Chen, Coexistent three-component and two-component weyl phonons in tis, zrse, and hfte, Phys. Rev. B 97, 054305 (2018).
  * 22. B. W. Xia, R. Wang, Z. J. Chen, Y. J. Zhao, H. Xu, Symmetry-protected ideal type-ii weyl phonons in cdte, Phys. Rev. Lett. 123, 065501 (2019).
  * 23. T. T. Zhang, H. Miao, Q. Wang, J. Q. Lin, Y. Cao, G. Fabbris, A. H. Said, X. Liu, H. C. Lei, Z. Fang, H. M. Weng, M. P. M. Dean, Phononic helical nodal lines with $\mathcal{PT}$ protection in ${\mathrm{mob}}_{2}$, Phys. Rev. Lett. 123, 245302 (2019).
  * 24. Y. Liu, X. Chen, Y. Xu, Topological phononics: From fundamental models to real materials, Adv. Funct. Mater. 30, 1904784– (2020).
  * 25. A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, S. I. Shylin, Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system, Nature 525, 73–76 (2015).
  * 26. K. Ahadi, L. Galletti, Y. Li, S. Salmani-Rezaie, W. Wu, S. Stemmer, Enhancing superconductivity in srtio<sub>3</sub> films with strain, Sci Adv 5, eaaw0120– (2019).
  * 27. S. Murakami, A. Okamoto, Thermal hall effect of magnons, J. Phys. Soc. Jpn. 86, 011010– (2016).
  * 28. J.-Y. Chen, S. A. Kivelson, X.-Q. Sun, Enhanced thermal hall effect in nearly ferroelectric insulators, Phys. Rev. Lett. 124, 167601 (2020).
  * 29. M. Hamada, S. Murakami, Phonon rotoelectric effect, Phys. Rev. B 101, 144306 (2020).
  * 30. X. Li, B. Fauqué, Z. Zhu, K. Behnia, Phonon thermal hall effect in strontium titanate, Phys. Rev. Lett. 124, 105901 (2020).
  * 31. R. E. Cohen, Origin of ferroelectricity in perovskite oxides, Nature 358, 136–138 (1992).
  * 32. Y. Moritomo, H. Kuwahara, Y. Tomioka, Y. Tokura, Pressure effects on charge-ordering transitions in perovskite manganites, Phys. Rev. B 55, 7549–7556 (1997).
  * 33. Y. Nahas, S. Prokhorenko, L. Louis, Z. Gui, I. Kornev, L. Bellaiche, Discovery of stable skyrmionic state in ferroelectric nanocomposites, Nature Communications 6, 8542– (2015).
  * 34. S. Das, Y. L. Tang, Z. Hong, M. A. P. Goncalves, M. R. McCarter, C. Klewe, K. X. Nguyen, F. Gómez-Ortiz, P. Shafer, E. Arenholz, V. A. Stoica, S.-L. Hsu, B. Wang, C. Ophus, J. F. Liu, C. T. Nelson, S. Saremi, B. Prasad, A. B. Mei, D. G. Schlom, J. Íniguez, P. García-Fernández, D. A. Muller, L. Q. Chen, J. Junquera, L. W. Martin, R. Ramesh, Observation of room-temperature polar skyrmions, Nature 568, 368–372 (2019).
  * 35. A. Mahmoud, A. Erba, K. E. El-Kelany, M. Rérat, R. Orlando, Low-temperature phase of batio3: Piezoelectric, dielectric, elastic, and photoelastic properties from ab initio simulations, Phys. Rev. B 89, 045103 (2014).
  * 36. K. H. Ahn, T. Lookman, A. R. Bishop, Strain-induced metal-insulator phase coexistence in perovskite manganites, Nature 428, 401–404 (2004).
  * 37. Y. Shirako, H. Kojitani, M. Akaogi, K. Yamaura, E. Takayama-Muromachi, High-pressure phase transitions of carho3 perovskite, Physics and Chemistry of Minerals 36, 455– (2009).
  * 38. L. L. Rusevich, G. Zvejnieks, E. A. Kotomin, M. M. Krzmanc, A. Meden, S. Kunej, I. D. Vlaicu, Theoretical and experimental study of (ba,sr)tio3 perovskite solid solutions and batio3/srtio3 heterostructures, J. Phys. Chem. C 123, 2031–2036 (2019).
  * 39. C. Paillard, E. Torun, L. Wirtz, J. Íñiguez, L. Bellaiche, Photoinduced phase transitions in ferroelectrics, Phys. Rev. Lett. 123, 087601 (2019).
  * 40. T. F. Nova, A. S. Disa, M. Fechner, A. Cavalleri, Metastable ferroelectricity in optically strained SrTiO3, Science 364, 1075– (2019).
  * 41. Y. Ahn, A. Pateras, S. D. Marks, H. Xu, T. Zhou, Z. Luo, Z. Chen, L. Chen, X. Zhang, A. D. DiChiara, H. Wen, P. G. Evans, Nanosecond optically induced phase transformation in compressively strained ${\mathrm{bifeo}}_{3}$ on ${\mathrm{laalo}}_{3}$, Phys. Rev. Lett. 123, 045703 (2019).
  * 42. M. Porer, M. Fechner, M. Kubli, M. J. Neugebauer, S. Parchenko, V. Esposito, A. Narayan, N. A. Spaldin, R. Huber, M. Radovic, E. M. Bothschafter, J. M. Glownia, T. Sato, S. Song, S. L. Johnson, U. Staub, Ultrafast transient increase of oxygen octahedral rotations in a perovskite, Phys. Rev. Research 1, 012005 (2019).
  * 43. D. M. Juraschek, M. Fechner, N. A. Spaldin, Ultrafast structure switching through nonlinear phononics, Phys. Rev. Lett. 118, 054101 (2017).
  * 44. H. F. KAY, R. G. RHODES, Barium titanate crystals, Nature 160, 126–127 (1947).
  * 45. P. Souvatzis, O. Eriksson, M. I. Katsnelson, S. P. Rudin, Entropy driven stabilization of energetically unstable crystal structures explained from first principles theory, Phys. Rev. Lett. 100, 095901 (2008).
  * 46. C. Fang, H. Weng, X. Dai, Z. Fang, Topological nodal line semimetals, Chinese Physics B 25, 117106– (2016).
  * 47. H. C. Po, A. Vishwanath, H. Watanabe, Symmetry-based indicators of band topology in the 230 space groups, Nature Communications 8, 50– (2017).
  * 48. Z. Song, T. Zhang, C. Fang, Diagnosis for nonmagnetic topological semimetals in the absence of spin-orbital coupling, Phys. Rev. X 8, 031069 (2018).
  * 49. J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, R.-J. Slager, Topological classification of crystalline insulators through band structure combinatorics, Phys. Rev. X 7, 041069 (2017).
  * 50. T. Zhang, L. Lu, S. Murakami, Z. Fang, H. Weng, C. Fang, Diagnosis scheme for topological degeneracies crossing high-symmetry lines, Phys. Rev. Research 2, 022066 (2020).
  * 51. W. Zhong, R. D. King-Smith, D. Vanderbilt, Giant lo-to splittings in perovskite ferroelectrics, Phys. Rev. Lett. 72, 3618–3621 (1994).
  * 52. Q. Wu, S. Zhang, H.-F. Song, M. Troyer, A. A. Soluyanov, WannierTools: An open-source software package for novel topological materials, Computer Physics Communications 224, 405–416 (2018).
  * 53. O. Diéguez, K. M. Rabe, D. Vanderbilt, First-principles study of epitaxial strain in perovskites, Phys. Rev. B 72, 144101 (2005).
  * 54. L. Ni, Y. Liu, C. Song, W. Wang, G. Han, Y. Ge, First-principle study of strain-driven phase transition in incipient ferroelectric srtio3, Physica B: Condensed Matter 406, 4145–4149 (2011).
  * 55. X. He, D. Bansal, B. Winn, S. Chi, L. Boatner, O. Delaire, Anharmonic eigenvectors and acoustic phonon disappearance in quantum paraelectric ${\mathrm{srtio}}_{3}$, Phys. Rev. Lett. 124, 145901 (2020).
  * 56. N. Choudhury, E. J. Walter, A. I. Kolesnikov, C.-K. Loong, Large phonon band gap in $\mathrm{Sr}\mathrm{Ti}{\mathrm{o}}_{3}$ and the vibrational signatures of ferroelectricity in $a\mathrm{Ti}{\mathrm{o}}_{3}$ perovskites: First-principles lattice dynamics and inelastic neutron scattering, Phys. Rev. B 77, 134111 (2008).
  * 57. X. Jia, S. Zhang, R. Sankar, F.-C. Chou, W. Wang, K. Kempa, E. W. Plummer, J. Zhang, X. Zhu, J. Guo, Anomalous acoustic plasmon mode from topologically protected states, Phys. Rev. Lett. 119, 136805 (2017).
  * 58. Y.-f. Yang, G.-M. Zhang, F.-C. Zhang, Universal behavior of the thermal hall conductivity, Phys. Rev. Lett. 124, 186602 (2020).
  * 59. J. J. Lee, F. T. Schmitt, R. G. Moore, S. Johnston, Y.-T. Cui, W. Li, M. Yi, Z. K. Liu, M. Hashimoto, Y. Zhang, D. H. Lu, T. P. Devereaux, D.-H. Lee, Z.-X. Shen, Interfacial mode coupling as the origin of the enhancement of tc in fese films on srtio3, Nature 515, 245–248 (2014).
  * 60. S. Roy, M. Kolodrubetz, N. Goldman, A. G. Grushin, Tunable axial gauge fields in engineered weyl semimetals: semiclassical analysis and optical lattice implementations, 2D Materials 5, 024001– (2018).
  * 61. H. Jia, R. Zhang, W. Gao, Q. Guo, B. Yang, J. Hu, Y. Bi, Y. Xiang, C. Liu, S. Zhang, Observation of chiral zero mode in inhomogeneous three-dimensional weyl metamaterials, Science 363, 148– (2019).
  * 62. V. Peri, M. Serra-Garcia, R. Ilan, S. D. Huber, Axial-field-induced chiral channels in an acoustic weyl system, Nature Physics 15, 357–361 (2019).
  * 63. J. Kim, K.-W. Kim, D. Shin, S.-H. Lee, J. Sinova, N. Park, H. Jin, Prediction of ferroelectricity-driven berry curvature enabling charge- and spin-controllable photocurrent in tin telluride monolayers, Nature Communications 10, 3965– (2019).
  * 64. P. O. Sukhachov, H. Rostami, Acoustogalvanic effect in dirac and weyl semimetals, Phys. Rev. Lett. 124, 126602 (2020).
  * 65. N. Farmakidis, N. Youngblood, X. Li, J. Tan, J. L. Swett, Z. Cheng, C. D. Wright, W. H. P. Pernice, H. Bhaskaran, Plasmonic nanogap enhanced phase-change devices with dual electrical-optical functionality, Sci Adv 5, eaaw2687– (2019).
  * 66. J. Feldmann, N. Youngblood, C. D. Wright, H. Bhaskaran, W. H. P. Pernice, All-optical spiking neurosynaptic networks with self-learning capabilities, Nature 569, 208–214 (2019).
  * 67. G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169–11186 (1996).
  * 68. G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 - 50 (1996).
  * 69. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * 70. J. Heyd, G. E. Scuseria, M. Ernzerhof, Hybrid functionals based on a screened Coulomb potential, J. Chem. Phys. 118, 8207 (2003).
  * 71. R. O. Jones, O. Gunnarsson, The density functional formalism, its applications and prospects, Rev. Mod. Phys. 61, 689–746 (1989).
  * 72. A. Görling, Density-functional theory for excited states, Phys. Rev. A 54, 3912–3915 (1996).
  * 73. A. Hellman, B. Razaznejad, B. I. Lundqvist, Potential-energy surfaces for excited states in extended systems, J. Chem. Phys. 120, 4593–4602 (2004).
  * 74. B. Peng, H. Zhang, W. Chen, B. Hou, Z.-J. Qiu, H. Shao, H. Zhu, B. Monserrat, D. Fu, H. Weng, C. M. Soukoulis, Sub-picosecond photo-induced displacive phase transition in two-dimensional mote2, npj 2D Materials and Applications 4, 14– (2020).
  * 75. F. Mauri, R. Car, First-principles study of excitonic self-trapping in diamond, Phys. Rev. Lett. 75, 3166–3169 (1995).
  * 76. S. Ismail-Beigi, S. G. Louie, Excited-state forces within a first-principles green’s function formalism, Phys. Rev. Lett. 90, 076401 (2003).
  * 77. S. Baroni, S. de Gironcoli, A. Dal Corso, P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515-562 (2001).
  * 78. A. Togo, I. Tanaka, First principles phonon calculations in materials science, Scripta Materialia 108, 1-5 (2015).
  * 79. A. A. Soluyanov, D. Vanderbilt, Computing topological invariants without inversion symmetry, Phys. Rev. B 83, 235401 (2011).
  * 80. R. Yu, X. L. Qi, A. Bernevig, Z. Fang, X. Dai, Equivalent expression of ${\mathbb{z}}_{2}$ topological invariant for band insulators using the non-abelian berry connection, Phys. Rev. B 84, 075119 (2011).
  * 81. D. Alfè, Phon: A program to calculate phonons using the small displacement method, Computer Physics Communications 180, 2622–2633 (2009).
  * 82. B. Cai, X. Chen, M. Xie, S. Zhang, X. Liu, J. Yang, W. Zhou, S. Guo, H. Zeng, A class of pb-free double perovskite halide semiconductors with intrinsic ferromagnetism, large spin splitting and high curie temperature, Mater. Horiz. 5, 961–968 (2018).
  * 83. M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, F. Bechstedt, Linear optical properties in the projector-augmented wave methodology, Phys. Rev. B 73, 045112 (2006).

## Acknowledgments

The authors gratefully acknowledge helpful discussions with Danny Bennett
(University of Cambridge) on the phase diagram of oxide perovskites, Prof.
Guntars Zvejnieks (University of Latvia) on the crystal structures of BaTiO3,
Bowen Hou and Prof. Hao Zhang (Fudan University) on the $\Delta$SCF method, Dr
Bo Cai (Nanjing University of Science and Technology) on finite-temperature
lattice dynamics, and Prof. Stefano Baroni (SISSA) on the screening effect of
photoexcited carriers. B.P. and B.M. acknowledge support from the Winton
Programme for the Physics of Sustainability, and B.M. also acknowledges
support from the Gianna Angelopoulos Programme for Science, Technology, and
Innovation. Y.H. acknowledges support from EPSRC grant EP/R512461/1 and
Trinity College Henry-Barlow Scholarship. S.M. and T.T.Z. acknowledge support
from the Tokodai Institute for Element Strategy (TIES) funded by MEXT Elements
Strategy Initiative to Form Core Research Center, and S.M. also acknowledges
support by JSPS KAKENHI Grant Number JP18H03678. The calculations were
performed using resources provided by the Cambridge Tier-2 system operated by
the University of Cambridge Research Computing Service
(http://www.hpc.cam.ac.uk) and funded by EPSRC Tier-2 capital grant
EP/P020259/1, and also with computational support from the UK Materials and
Molecular Modelling Hub, which is partially funded by EPSRC (EP/P020194), for
which access was obtained via the UKCP consortium and funded by EPSRC grant
ref EP/P022561/1.

Competing Interests: The authors declare that they have no competing
interests.

Author Contributions: B.M., T.T.Z. and B.P. devised the project idea. B.P.,
Y.H. and T.T.Z. performed the calculations. B.P, T.T.Z. and B.M. prepared the
main part of the manuscript. B.P., Y.H., S.M., T.T.Z. and B.M. discussed the
results and the ideas for their analysis and edited the manuscript.

Data Availability: All data needed to evaluate the conclusions in the paper
are present in the paper and/or the Supplementary Materials. Additional data
related to this paper may be requested from the authors.

## Supplementary materials

Supplementary Text
Figs. S1 to S9


# Ranking labs-of-origin for genetically engineered DNA using Metric Learning

Fernando H. F. Camargo<EMAIL_ADDRESS>AmalgamP.O. Box
1212AustinTexasUSA43017-6221 ,  Igor Muniz Soares<EMAIL_ADDRESS>AmalgamP.O.
Box 1212AustinTexasUSA43017-6221 and  Adriano Marques<EMAIL_ADDRESS>Amalgam1 Thørväld CircleAustinTexasUSA

###### Abstract.

With the constant advancements of genetic engineering, a common concern is to
be able to identify the lab-of-origin of genetically engineered DNA sequences.
For that reason, AltLabs has hosted the Genetic Engineering Attribution
Challenge to gather many teams to propose new tools to solve this problem.
Here we show our proposed method to rank the most likely labs-of-origin and
generate embeddings for DNA sequences and labs. These embeddings can also
perform various other tasks, like clustering both DNA sequences and labs and
using them as features for Machine Learning models applied to solve other
problems. This work demonstrates that our method outperforms the classic
training method for this task while generating other helpful information.

lab-of-origin, genetic engineering attribution, metric learning, triplet
network, deep learning, dna, rna, plasmid

††ccs: Information systems Similarity measures††ccs: Information systems
Learning to rank

##  1\. Introduction

Genetic engineering or genetic modification is one of the fastest-growing
fields in biotechnology. It is responsible for manipulating an organism’s
genes. It is a well-known technology with applications in numerous fields,
bringing several benefits to science, medicine, and industry, which is why it
is receiving increasing investments. Changing DNA makes it possible to create
drugs and vaccines while reducing costs and time to market in health and
industrial applications. However, alongside the benefits, there are possible
situations of improper usage. Modification of genes can lead to the accidental
or intentional release of new and dangerous biological agents.

Detecting genetic modifications and attributing them to a specific origin is a
recent and open problem. Tools have been developed for detecting genetic
engineering through genetic analysis of the organism genome. This is the case
of FELIX (Finding Engineering-Linked Indicators) (IARPA, [n.d.]), a program
developed by IARPA that can rapidly detect indicators of genome engineering
and has recently demonstrated its power by demystifying the rumors that the
COVID-19 virus was created in a laboratory. However, assigning the origin of a
specific engineered DNA to a lab is a more complex challenge.

AltLabs, a not-for-profit organization, hosted the Genetic Engineering
Attribution Challenge (com, [n.d.]), providing the participants with a dataset
of 63,017 DNA sequences (with their phenotype characteristics) designed by
1,314 different labs. The goal was to design machine learning models to, given
a DNA sequence, identify the most likely lab-of-origin. Since it is tough to
predict the correct lab because they might use similar design techniques, they
decided to evaluate the solutions using top ten accuracy (Alley et al., 2020).
It means that the model needs to place the correct lab-of-origin within the
ten most likely labs, according to the scores. It turns this problem into a
ranking problem instead of a simple classification, which might ask for
different techniques.

During the competition, we decided to create two branches. The first is a
traditional approach using classification models (with a softmax to output
probabilities for each lab). The second is the use of ranking models through
Metric Learning (Kulis, 2013) (more specifically, Triplet Networks (Hoffer and
Ailon, 2018)). This one has the goal of learning how to extract
embeddings111Embedding is a vectorial representation of some entity composed
of latent features. Humans cannot easily understand these features, but they
are beneficial for machine learning models. When we train a model to extract
embeddings, the goal is to place similar entities closer together in the
latent space. (also known as feature vectors) from DNA sequences and, at the
same time, learn the embeddings of the labs. A similarity measure (cosine
similarity) is then applied to generate the score between a DNA sequence and a
lab. This score, instead of trying to mimic a probability, describes how
similar they are.

Despite the difference in how we train those two kinds of models, we designed
them very similarly. Both of them shared the same preprocessing steps, and the
first layers were almost the same. They diverged a little bit because some
techniques worked better for one than the other. Nevertheless, the overall
structure of these models was mostly the same. Even though they were very
similar, our triplet networks consistently outperformed our classification
models, indicating they are more suitable for ranking tasks.

Furthermore, it is not only about raw accuracy. Suppose a new sequence from an
unknown lab goes through a regular classification model. In that case, it will
probably attribute it to one of the known labs while giving no clue about its
uncertainty. It is a known problem and has encouraged research in an area
called Out-Of-Distribution Detection (Chen et al., 2020)(Ren et al., 2019). On
the other hand, our triplet network is more robust to this problem. The
embedding extracted from such a sequence might differ significantly from all
the labs’ embeddings, generating low similarity scores. This way, we could set
a threshold and decide that it comes from an unknown lab if the score is below
it. There is also an embedding learned for an unseen lab. We train the model
to push it away from all the known labs. So, suppose the new sequence’s
embedding is more similar to this one. In that case, we can also assume an
out-of-distribution sequence. In summary, our proposed solution has a high
accuracy while also dealing with new labs and providing embeddings for
clustering and other tasks.

Figure 1. We used t-SNE (Hinton and Roweis, 2002) to project 200D into 2D and
produce this Figure. Each circle represents a DNA sequence, and the color
shows the lab-of-origin of such sequence. The lab is also presenting using a
text with its ID. We can also see that there are many sequences very similar
in the middle, even though they come from different labs. We theorize that
they are from small labs that use similar design techniques. It is also worth
noting that dimensionality reduction causes some losses, which might also
explain it.

##  2\. Proposed Method

###  2.1. Implementation Details

Our proposed method consists of different preprocessing techniques proved to
be efficient for DNA sequences in our experiments, along with a Neural
Network. As a preprocessing step, we used Byte Pair Encoding (BPE) algorithm
(Gage, 1994) to compress the DNA sequences. It works by looking for common
patterns and unifying them into tokens, increasing the vocabulary while
reducing the sequences’ size. It turned our vocabulary of 4 DNA bases into
1001 different tokens.

To deal with the fact that DNA sequences are circular, we used a circular
shift data augmentation. This augmentation was used during training and
inference, doing what is known as Test Time Augmentation. It means we feed the
model with multiple versions of the same sequence (randomly circular shifted).
We then take the average of the outputs. Further, to decrease the model’s
complexity, we also limited the sequence to 1000 tokens. This operation runs
after the circular shift, so the model receives different 1000 tokens
sequences per execution and combines the results.

Figure 2. (A) In the beginning, the negative (another lab) might be closer to
the anchor (DNA) than the positive (lab-of-origin). During training, we pull
the anchor and lab-of-origin towards each other while pushing the negative
away. This way, the positive will be closer to the anchor, and the negative
will be farther. (B) In the end, labs and their DNA sequences will be closer
to each other, forming groups. We can also expect both labs and DNA sequences
to be closer to other similar ones.

Both types of models are composed of a Convolutional Neural Network with
multiple kernels of different sizes, as proposed by Kim. We use it to extract
features from the sequence and concatenate them with the binary features
provided in the dataset. The difference between the classification and triplet
network models are after this part.

The classification model takes these features and pass them through a dense
layer and then the output layer. This output layer is a dense layer with
softmax as the activation function, resulting in probabilities for each lab.
On the other hand, the triplet network passes these features through a dense
layer that will generate our sequence embedding. In parallel, we have an
embedding layer that will learn the lab embeddings. For last, we use a cosine
similarity measure between the embeddings to give us the final output.

While we train the classification model using regular Supervised Learning, we
train the Triplet Network differently. We create triplets of anchors,
positives, and negatives and use them for training the model. In a Face
Recognition task, the anchor would be one picture of one person, the positive
would be another picture of the same person, and the negative would be any
picture of another person. The model learns to pull the anchor and positive
together (making pictures of the same person produce similar embeddings) while
pushing the negative away. In our case, the DNA sequence is the anchor, while
the positive is the source lab, and the negative is another lab. Our goal is
to generate embeddings in which the DNA sequences are close to their source
lab and far from other labs and their sequences. Figure 2 summarizes this
process.

To generate these triplets, we use the dataset to provide us with the anchor
(sequence) and positive (source lab). We then use a technique known as Hard
Negative Mining (Hermans et al., 2017) to select the negative (another lab).
It means that instead of selecting a random lab as a negative example, we
select the most challenging example. It would be the lab that is currently
closer to our sequence in the latent space.

###  2.2. Triplet Network Advantages

After training, we can use the Triplet Network to extract embeddings from DNA
sequences. It also contains a table of embeddings for each lab present in the
dataset and one embedding for an ”unseen” lab. These embeddings can be helpful
beyond the expected usage (ranking possible labs-of-origin given a DNA
sequence). We can also use them to compare and cluster DNA sequences and labs
in terms of design style, as shown in Figure 1. Speaking of which, we applied
K-Means (Lloyd, 1982) to cluster the lab embeddings to find out what is a good
number of groups for the 1,314 labs available. Our result, using the Elbow
Method, is that we can comfortably put them into nine groups. Such clustering
capabilities of our model could be instrumental in investigating the
relationship between these labs.

We also expect this model to be very robust for unseen labs. As pointed out in
Section 1, we can set a threshold and use the unseen embedding to detect that
a DNA sequence comes from an unknown lab. Besides that, suppose that an
accidental release happened and we need to find out the responsible. Even
though some labs are not part of our training dataset, we could ask for some
samples of sequences designed and extract their embeddings. We could average
these embeddings to generate an embedding for the lab. With this embedding at
hand, we could compute this lab’s similarity with the DNA sequence that was
accidentally released. It is worth noting that we can do that without
requiring the model to be retrained.

For last, these embeddings are also very feature-rich, which means we can use
them as input for other machine learning models, tackling other problems. For
example, we know that it is common to design DNA sequences derived from
multiple different sources. Suppose we want to identify the composition of a
given DNA. We could prepare a dataset and use our embeddings to serve as
inputs for another model. It would drastically accelerate the creation of such
a model.

Input: $lab\\_indices$ (B,), $anchor\\_embeddings$ (B, L),
$positive\\_embeddings$ (B, E)

Output: Batch of $negatives$ of shape (B, E)

1 $all\\_lab_{i}ndices\leftarrow$ indices of all labs repeated by B // (B, L)

2 $negative\\_labs\\_mask\leftarrow$ matrix of boolean values with $True$ for
every lab // (B, L)

3 $negative\\_labs\\_mask[:,lab\\_indices]\leftarrow False$ // (B, L)

4 // At this point, we have a mask with True for the negative labs and False
otherwise

5 $all\\_negative\\_lab\\_indices\leftarrow$
$all\\_labs[negative\\_labs\\_mask]$ reshaped // (B, L-1)

6 $all\\_negative\\_labs\\_embeddings\leftarrow
lab\\_embeddings[all\\_negative\\_labs]$ L2 normalized // (B, L-1, E)

7 $anchor\\_similarities\leftarrow$ dot product between $anchor\\_embeddings$
and $all\\_negative\\_labs\\_embeddings$ // (B, L-1)

8 $hardest\\_negative\\_lab_{i}ndices\leftarrow$ argmax of
$anchor\\_similarities$ // (B,)

9 $negative\\_embeddings\leftarrow
all\\_negative\\_labs\\_embeddings[:,hardest\\_negative\\_lab\\_indices]$ //
(B, E)

10 return _$negative\\_embeddings$_

Algorithm 1 Algorithm for Hard Negative Mining using tensors. The shape of
each tensor is at the end of each line as a comment, being B batch size, E the
embedding dimension, and L is the number of labs. It is worth noting that we
L2 normalize all the embeddings.

##  3\. Results

During the competition, we tuned the classification model and triplet network
to improve their performance. As a result, our best triplet network got 90.39%
top ten accuracy in the test dataset. In contrast, our best classification
model got 89.36%. Thus, it shows that the triplet network displays many
advantages and performs better than a classification model in such a ranking
task. It is worth noting that we are comparing both of them using a very
similar architecture. At the end of the competition, other competitors created
more accurate classifiers. However, we theorize that if we trained their
architecture with our scheme, they would be even more precise.

For the competition, we decided to combine these models into an ensemble to
boost our accuracy further. However, since they give very different outputs,
this approach posed several challenges. We first tried to compute each lab’s
average rank position known as Borda voting rule (Borda, 1781). However, the
results were worse than the individual models. We then decided to try the
Copeland voting rule (Saari and Merlin, 1996), which performed much better.
Our final submission using this ensemble got 91.67% top ten accuracy and the
10th position in the competition. We also got the 2nd position in the
Innovation Track, meaning that the organizers recognized our solution as
innovative.

##  4\. Conclusions and Future Work

Metric Learning is quite common in other areas (like Face Recognition and
Recommender Systems). However, to the best of our knowledge, it has never been
used in this context. We showed an innovative solution with promising results,
outperforming a classification model trained with Supervised Learning.
Moreover, it can also generate embeddings as a by-product, which are very
useful.

We believe that more advanced techniques could be applied to extract features
from the DNA sequences to improve them even further. New techniques like
Transformers and Graph Convolutional Networks would be good candidates for
this task. With better feature extraction, the embeddings would have higher
quality, improving overall accuracy.

## References

  * (1)
  * com ([n.d.]) [n.d.].  Genetic Engineering Attribution Challenge .  https://www.drivendata.org/competitions/63/genetic-engineering-attribution/page/164/.  https://www.drivendata.org/competitions/63/genetic-engineering-attribution/page/164/ Accessed: 2020-10-22.
  * Alley et al. (2020) Ethan C. Alley, Miles Turpin, Andrew Bo Liu, Taylor Kulp-McDowall, Jacob Swett, Rey Edison, Stephen E. Von Stetina, George M. Church, and Kevin M. Esvelt. 2020.  Attribution of genetic engineering: A practical and accurate machine-learning toolkit for biosecurity.  _bioRxiv_ (2020).  https://doi.org/10.1101/2020.08.22.262576 arXiv:https://www.biorxiv.org/content/early/2020/08/22/2020.08.22.262576.full.pdf
  * Borda (1781) J-C de Borda. 1781\.  Mémoire sur les élections au scrutin: Histoire de l’Académie Royale des Sciences.  _Paris, France_ 12 (1781).
  * Chen et al. (2020) Jiefeng Chen, Yixuan Li, Xi Wu, Yingyu Liang, and Somesh Jha. 2020.  Robust Out-of-distribution Detection for Neural Networks.  arXiv:2003.09711 [cs.LG]
  * Gage (1994) Philip Gage. 1994\.  A New Algorithm for Data Compression.  _C Users J._ 12, 2 (Feb. 1994), 23–38.
  * Hermans et al. (2017) Alexander Hermans, Lucas Beyer, and Bastian Leibe. 2017\.  In Defense of the Triplet Loss for Person Re-Identification.  arXiv:1703.07737 [cs.CV]
  * Hinton and Roweis (2002) Geoffrey Hinton and Sam Roweis. 2002.  Stochastic Neighbor Embedding. In _Proceedings of the 15th International Conference on Neural Information Processing Systems_ _(NIPS’02)_. MIT Press, Cambridge, MA, USA, 857–864.
  * Hoffer and Ailon (2018) Elad Hoffer and Nir Ailon. 2018.  Deep metric learning using Triplet network.  arXiv:1412.6622 [cs.LG]
  * IARPA ([n.d.]) IARPA. [n.d.].  FELIX.  https://www.iarpa.gov/index.php/research-programs/felix
  * Kim (2014) Yoon Kim. 2014.  Convolutional Neural Networks for Sentence Classification.  arXiv:1408.5882 [cs.CL]
  * Kulis (2013) Brian Kulis. 2013\.  Metric Learning: A Survey.  _Foundations and Trends® in Machine Learning_ 5, 4 (2013), 287–364.  https://doi.org/10.1561/2200000019
  * Lloyd (1982) S. Lloyd. 1982\.  Least squares quantization in PCM.  _IEEE Transactions on Information Theory_ 28, 2 (1982), 129–137.  https://doi.org/10.1109/TIT.1982.1056489
  * Ren et al. (2019) Jie Ren, Peter J. Liu, Emily Fertig, Jasper Snoek, Ryan Poplin, Mark A. DePristo, Joshua V. Dillon, and Balaji Lakshminarayanan. 2019.  Likelihood Ratios for Out-of-Distribution Detection.  arXiv:1906.02845 [stat.ML]
  * Saari and Merlin (1996) Donald G. Saari and Vincent R. Merlin. 1996.  The Copeland method.  _Economic Theory_ 8, 1 (01 Feb 1996), 51–76.  https://doi.org/10.1007/BF01212012