# Torus bifurcations of large-scale swarms having range dependent
communication delay

Ira B. Schwartz<EMAIL_ADDRESS>U.S. Naval Research Laboratory,
Code 6792, Plasma Physics Division, Washington, DC 20375, USA  Victoria
Edwards  U.S. Naval Research Laboratory, Code 5514, Navy Center for Applied
Research in Artificial Intelligence, Washington, DC 20375, USA  Sayomi
Kamimoto  Department of Mathematics, George Mason University, Fairfax
Virginia, 22030, USA  Klimka Kasraie  Aerospace, Transportation and Advanced
Systems Laboratory of the Georgia Tech Research Institute, Atlanta, GA 30332
M. Ani Hsieh  Mechanical Engineering and Applied Mechanics University of
Pennsylvania, Philadelphia, PA 19104 USA  Ioana Triandaf1 Jason Hindes1

###### Abstract

Dynamical emergent patterns of swarms are now fairly well established in
nature, and include flocking and rotational states. Recently, there has been
great interest in engineering and physics to create artificial self-propelled
agents that communicate over a network and operate with simple rules, with the
goal of creating emergent self-organizing swarm patterns. In this paper, we
show that when communicating networks have range dependent delays, rotational
states which are typically periodic, undergo a bifurcation and create swarm
dynamics on a torus. The observed bifurcation yields additional frequencies
into the dynamics, which may lead to quasi-periodic behavior of the swarm.

Swarming behavior occurs when a large number of self-propelled agents interact
using simple rules. Natural swarms of biological systems have been observed at
a range of length scales forming complex emergent patterns. Engineers have
drawn inspiration from these natural systems, resulting in the translation of
swarm theory to communicating robotic systems. Example applications of
artificial swarms include: exploration and mapping, search and rescue, and
distributed sensing and estimation. Through continued development, an
additional parameter of delay in communication between artificial agents has
become important to consider. Specifically, it was previously discovered, that
communication delay will create new rotational patterns which are not observed
without delay, both theoretically and experimentally. Here we extend the
understanding of communication delays to reveal the effects of range dependent
delay, where the communication between agents depends on the distance between
agents. The results of the research show that by including range dependent
delay, new rotational states are introduced. We show how these new states
emerge, discuss their stability, and discuss how they may be realized in large
scale robotic systems. In improving our theoretical understanding of predicted
swarm behavior modeled in simulation we can better anticipate what will happen
experimentally. Additionally, it is possible to leverage the predicted
autonomous behaviors to try and force different swarm behavior.

##  I Introduction

Swarming behavior, which we define as the emergence of spatio-temporal group
behaviors from simple local interactions between pairs of agents, is
widespread and observed over a range of application domains. Examples can be
found in biological systems over a range of length scales, from aggregates of
bacterial cells and dynamics of skin cells in wound healing Budrene and Berg
(1995); Polezhaev  _et al._ (2006); Lee  _et al._ (2013) to dynamic patterns
of fish, birds, bats, and even humans Tunstrø m  _et al._ (2013); Helbing and
Molnar (1995); Giuggioli, McKetterick, and Holderied (2015a); Lee (2006).
These systems are particularly interesting because they allow simple
individual agents to achieve complex tasks in ways that are scalable,
extensible, and robust to failures of individual agents. In addition, these
swarming behaviors are able to form and persist in spite of complicating
factors such as delayed actuation, latent communication, localized number of
neighbors each agent is able to interact with, heterogeneity in agent
dynamics, and environmental noise. These factors have been the focus of
previous theoretical research in describing the bifurcating spatial-temporal
patterns in swarms, as seen for example in Refs. Topaz and Bertozzi (2004);
Szwaykowska, Romero, and Schwartz (2015); Mier-y-Teran Romero, Forgoston, and
Schwartz (2011); Hindes, Szwaykowska, and Schwartz (2016). Likewise, the
application of swarms have been experimentally realized in areas, such as
mappingRamachandran, Elamvazhuthi, and Berman (2018), leader-followingMorgan
and Schwartz (2005); Wiech, Eremeyev, and Giorgio (2018), and density
controlLi  _et al._ (2017). To guarantee swarming behavior experimentally,
control is typically employed Tanner, Jadbabaie, and Pappas (2007); Gazi
(2005); Jadbabaie, Jie Lin, and Morse (2003); Viragh  _et al._ (2014); Desai,
Ostrowski, and Kumar (2001) to prove convergence to a given state by relying
on strict assumptions to guarantee the desired behavior. However, by relaxing
certain assumptions, a number of studies show that even with simple
interaction protocols, swarms of agents are able to converge to organized,
coherent behaviors in a self-emergent manner; i.e. autonomously without
control. Different mathematical approaches yielded a wide selection of both
agent-based Helbing and Molnar (1995); Lee (2006); Vicsek  _et al._ (2006);
Tunstrø m  _et al._ (2013) and continuum models that predict swarming
dynamics. Edelstein-Keshet, Grunbaum, and Watmough (1998); Topaz and Bertozzi
(2004); Polezhaev  _et al._ (2006). In almost all models, since the agents
have just a few simple rules, there exists only a relatively small number of
controllable parameters. The parameter set usually consists of a self-
propulsion force, a potential function governing attracting and repelling
forces between agents, and a communicating radius governing the local
neighborhood at which the agents can sense and interact with each other.

In both robotic and biological swarms, an additional parameter appears as a
delay between the time information is perceived and the actuation (reaction)
time of an agent. Such delays have now been measured in swarms of bats, birds,
fish, and crowds of peopleGiuggioli, McKetterick, and Holderied (2015b); Nagy
_et al._ (2010); Fehrenbach  _et al._ (2014). The measured delays are longer
than the typical relaxation times of the agents, and may be space and time
dependent. Robotic swarms experience communication delays which provide
similar effects to the delay experienced in natural swarms. Incorporating
stationary delays along with a minimal set of parameters in swarm models
results in multi-stability of rotational patterns in space Mier-y-Teran
Romero, Forgoston, and Schwartz (2012); Szwaykowska  _et al._ (2016); Edwards
_et al._ (2020); Hindes and Schwartz (2018); Szwaykowska, Schwartz, and Carr
(2018). In particular, for delays that equal and fixed, one observes three
basic swarming states or modes: Flocking, which is a translating center of
mass, Ring state, where the agents are splayed out on a ring in in phase about
a stationary center of mass, and a Rotating state, where the center of mass
itself rotates.

Synthetic robotic swarms have communication delays that naturally occur over
wireless networks, as a result of low bandwidthKomareji, Shang, and Bouffanais
(2018) resulting in delayed communication and multi-hop communicationOliveira,
Almeida, and Lima (2015). In cases where the delays are fixed and equal, and
the communication occurs on a homogeneous network, it is known that delays
create new rotational patterns , as has been verified both theoretically and
experimentally Szwaykowska  _et al._ (2016); Edwards  _et al._ (2020).
However, in situations with robots, even simple communication models are based
on the distance between agents ying Ani Hsieh  _et al._ (2004); Hsieh  _et
al._ (2007). Following from these models, if one assumes that the delays are
range dependent, the problem becomes one of studying state dependent delays
where delays depend implicitly on the relative positions between agents.

When placing swarms in realistic complex environments, delays are not
necessarily a continuous function of range, but rather it is the increasing
probability of delays increasing stochastically when agents move further away
from one another beyond a certain radius Fink, Ribeiro, and Kumar (2012,
2013). That is, the rate of communication becomes spatially dependent, whereby
near agents see a signal with a fast rate of communication, but due to shading
and fading of signals, communication rates are slowed and complex outside a
given radius. Underwater communication is an excellent swarm example where
delays outside a significant radius impart rates of communication of one to
two orders of magnitude greater than local communication rates Arrichiello
_et al._ (2009).

The swarm model that follows takes a globally coupled swarm, and explicitly
relaxes the fixed delay assumption, by including range dependent delay based
on a fixed communication radius. We show that when range dependent delays are
included, new frequencies are introduced and generate bifurcations to a torus.
The result is a milling type of swarm that depends on just a few parameters.
The results here are important for robotic swarming where one of the goals is
to produce desired patterns autonomously, without external controls. The
pattern formations predicted here show how delayed information, whether coming
from communication, actuations, or both, impacts the stability of swarm
states, such as ring and/or rotating states. By revealing those parameter
regions where patterns are destabilized, we provide a comprehensive
characterization of the autonomously accessible swarm states in the presence
of range-dependent delay.

##  II The swarm model

Consider a swarm of delay-coupled agents in $\bm{{R}}^{2}$. Each agent is
indexed by $i\in\\{1,\ldots,N\\}$. We use a simple but general model for
swarming motion. Each agent has a self-propulsion force that strives to
maintain motion at a preferred speed and a coupling force that governs its
interaction with other agents in the swarm. The interaction force is defined
as the negative gradient of a pairwise interaction potential $U(\cdot,\cdot)$.
All agents follow the same rules of motion; however, mechanical differences
between agents may lead to heterogeneous dynamics; this effect is captured by
assigning different acceleration factors (denoted $\kappa_{i}$) to the agents.
In this paper, we assume $\kappa_{i}=1$ for all $i$. For the effect of
heterogeneity on the swarm bifurcations, see Szwaykowska, Romero, and Schwartz
(2015).

Agent-to-agent interactions occur along a graph
$\mathcal{{G}}=\\{\mathcal{{V}},\mathcal{{E}}\\}$, where $\mathcal{{V}}$ is
the set of vertexes $v_{i}$ in the graph and $\mathcal{{E}}$ is the set of
edges $e_{ij}$. The vertices correspond to individual swarm agents, and edges
represent communication links; that is, agents $i$ and $j$ communicate with
each other if and only if $e_{ij}\in\mathcal{{E}}$. All communications links
are assumed to be bi-directional, and all communications occur with a time
delay $\tau$. That is, range dependence is not included. Let
${\bm{{r}}}_{i}\in\mathbb{R}^{2}$ denote the position of the agent $i$ and let
$\mathcal{{N}}_{i}=\\{v_{j}\in\mathcal{{V}}\mathrel{\mathop{\mathchar
58\relax}}e_{ij}\in\mathcal{{E}}\\}$ denote its set of neighbors of agent $i$.
The motion of agent $i$ is governed by the following equation:

$\ddot{\mathbf{r}}_{i}=\kappa_{i}(1-\mathinner{\\!\left\lVert\dot{\mathbf{r}}_{i}\right\rVert}^{2})\dot{\mathbf{r}}_{i}-\kappa_{i}\sum_{j\in\mathcal{{N}}_{i}}\nabla_{x}U({\mathbf{r}}_{i}(t),{\mathbf{r}}_{j}^{\tau}(t)),$
(1)

where superscript $\tau$ is used to denote time delay, so that
${\mathbf{r}}_{j}^{\tau}(t)={\mathbf{r}}_{j}(t-\tau)$,
$\mathinner{\\!\left\lVert\cdot\right\rVert}$ denotes the Euclidean norm, and
$\nabla_{x}$ denotes the gradient with respect to the first argument of $U$.
The first term in Eq. 1 governs self-propulsion, where the speed has been
normalized to unity. That is, without coupling the agents always asymptote to
unit speed.

To analyze the dynamics of a large scale swarm, we use a harmonic interaction
potential with short-range repulsion.

$U(\mathbf{r}_{i},\mathbf{r}_{j}^{\tau})=c_{r}e^{-\frac{\mathinner{\\!\left\lVert\mathbf{r}_{i}-\mathbf{r}_{j}\right\rVert}}{l_{r}}}+\frac{a}{2N}\mathinner{\\!\left\lVert\mathbf{r}_{i}-\mathbf{r}_{j}^{\tau}\right\rVert}^{2}.$
(2)

In Eq. 1, it is assumed that the communication delay, $\tau$, is independent
of the distance, or range, between any pair of agents. (Notice that the
exponent of the repulsion term is independent of the delay since the repulsion
force is local.) With the addition of delays in the network, it was shown in
homogeneous communication networks that in addition to the usual dynamical
translating and milling (or ring) states, for sufficiently large $\tau$, new
rotational states emerge Szwaykowska  _et al._ (2016). In particular, for a a
given attractive coupling strength, there is a delay that destabilizes the
periodic ring state into a rotating state, in which the agents coalesce to a
small group and move around a fixed center of rotation; this behavior is quite
different from the ring state where agents are spread out in a splay state
phase. The rotating state is only observed with delay introduced in the
communication network, and it appears through a Hopf bifurcation.

However, in real-world robotic swarms, communication delays are not uniform
between all pairs of agents; delays may be stochastic or even state-dependent.
For example, if agents are communicating over a multi-hop network, the delay
will increase with the number of hops required to send a message from one
agent to the other, and in general will scale with the separation between
them. In order to handle range dependent delays, we will make an approximation
that depends on a communication range radius.

###  II.1 Approximating range dependent delayed coupling

For the coupling term, we are interested in introducing an approximation to
range based coupling delay. Since all communicating agents send signals with
some delay, we compute relative distances defined as

$\displaystyle D^{\tau}_{i,j}$
$\displaystyle\equiv||{\mathbf{r}}_{i}-{\mathbf{r}}_{j}^{\tau}||.$ (3)

We define a Heaviside function, $H(x)$, that is zero when $x\leq 0$ and 1
otherwise, and we employ global coupling based on a spring potential. For our
range dependent metric, we let $\epsilon\geq 0$ denote the range radius.
Suppose that when the separation between two agents is small, that is less
than $\epsilon$, then sensing between two agents is almost immediate. In
practice, the time needed for sensing depends on several factors, such as
actuation times, and so distances in practice are computed with delay.
Therefore, we model the coupling term for the $i^{th}$ agent as

$\displaystyle
C_{i}({\mathbf{r}}_{i},{\mathbf{r}}_{j},{\mathbf{r}}_{j}^{\tau},\epsilon)=-\frac{a}{N}(\nabla_{x}U({\mathbf{r}}_{i}(t),{\mathbf{r}}_{j}^{\tau}(t)))H(D^{\tau}_{i,j}-\epsilon)$
$\displaystyle-\frac{a}{N}(\nabla_{x}U({\mathbf{r}}_{i}(t),{\mathbf{r}}_{j}(t)))(1-H(D^{\tau}_{i,j}-\epsilon)),$
(4)

where the first coupling term has delay turned on since the distance is
outside a ball of radius $\epsilon$, while the second term has no delay since
the distance is within the $\epsilon$ ball. The resulting swarm model with
range dependence from Eq. II.1 is now

$\ddot{\mathbf{r}}_{i}=\kappa_{i}(1-\mathinner{\\!\left\lVert\dot{\mathbf{r}}_{i}\right\rVert}^{2})\dot{\mathbf{r}}_{i}-\kappa_{i}\sum_{j\in\mathcal{{N}}_{i}}C_{i}({\mathbf{r}}_{i},{\mathbf{r}}_{j},{\mathbf{r}}_{j}^{\tau},\epsilon).$
(5)

If the delayed distance is within an $\epsilon$ ball, then we evaluate the
coupling without delay. Otherwise the coupling is delayed. Thus the coupling
function takes into account when delay is active or not between pairs of
communicating agents, and depends on the range radius, $\epsilon$.

The Heaviside function of the right hand side of Eq. 9 renders the
differential delay equation derivatives discontinuous, and as such poses a
numerical integration problem. To mollify the lack of smoothness, we
approximate $H(x)$ by letting
$H(x)\approx\frac{1}{\pi}\arctan(kx)+\frac{1}{2}$, where $k\gg 1$ and
constant, and limits on the Heaviside function as $k\rightarrow\infty$.

Using only the delayed distance to compute a range dependent coupling assumes
that any measurement is not instantaneous. If one were to be able to compute
the ideal situation where delay would not be a sensing factor, then certain
issues would need to be resolved, which we do not consider here.

###  II.2  Numerical simulations of full swarms

Examples of simulations using the swarm model with the range dependent
coupling are shown below. Here the number of agents, $N=150$, and the coupling
strength, $a=2.0$. For the remainder of the analysis, we set $c_{r}=0$, and
note that the attractors persist when the repulsive amplitude is sufficiently
small Szwaykowska  _et al._ (2016). (See supplementary material for a video of
the dynamics with small repulsion.)

Figure 1: Three snapshots of swarm state in space for
$\epsilon=0.01,a=2.0,\tau=1.75$. Sample times
$t_{0},t_{1}=t_{0}+20,t_{2}=t_{0}+40$ .

Note that even when $\epsilon$ is very small, as shown in Fig. 1, we observe a
mix of clustered states which are a combination of pure ring and rotation
states. The agents tend to cluster into local groups, and the clusters move in
clockwise and counter-clockwise directions as in the ring state. Here,
however, the phase differences between agents are non-uniform. When examining
a single random agent, as shown in Fig 2, it is periodic with a sharp
frequency of rotation, and the relative positions of all individual agents are
phase locked. When considering the center of mass of the positions over all
agents, $\bm{{R}}\equiv\frac{1}{N}\sum_{i}\bm{{r}}_{i}$, the center of mass
does small amplitude oscillations about a fixed point (not shown).

Figure 2: Swarm ring state for $\epsilon=0.01,a=2.0,\tau=1.75$. (a)Time series
of the x-component of a single agent. (b)The power spectrum showing a sharp
frequency. (c)A phase portrait of the orbit of a single agent. The red point
denotes the center of mass.

As the radius $\epsilon$ increases, instability of the periodic mixed state
occurs, giving rise to more complicated behavior, as seen in Fig. 3. New
frequencies are introduced, causing the ring state to appear as a quasi-
periodic attractor. Moreover, the dynamics of the center of mass has its own
non-trivial dynamics which includes the effects of new frequencies. By
examining the Poincare map of the attractors, the instability gives rise to
dynamics which we conjecture is motion on a torus. Letting $(M_{x},M_{y})$
denote the time averaged center of mass over all agents, we compute the
sequence ${x(t_{i}),i=1..M}$ when $y(t_{i})=0$ and $x(t_{i})>M_{x}$. The
result is shown in the two panels in Fig. 4. Panel (a) shows a complicated
toroidal motion after transients are removed of the center of mass in Fig. 3c.
For a single frequency, the dynamics of the center of mass would be a single
fixed point. The addition of new frequencies is revealed in the Poincare map
as complicated motion on a torus. For larger values of $\epsilon$, the motion
on the torus converges to a periodic attractor in panel (b).

Figure 3: Swam instability $\epsilon=0.25,a=2.0,\tau=1.75$. (a)Time series of
the x-component of a single agent. (b)The Power spectrum showing a slight
broadening and birth of a new frequency. (c)A phase portrait of the orbit of a
single agent.

Figure 4: Poincare map of Eqs. 1-II.1 for (a) $\epsilon=0.25$, (b)
$\epsilon=0.5$. Other parameters are fixed: $a=2.0,\tau=1.75$. See text for
details.

##  III Mean-Field Equation of Range Dependent Delay Coupled Swarm

In order to shed some light on the origin of the bifurcation to dynamics on a
torus, we examine the full swarm model from a mean-field perspective. The mean
field is much lower dimensional, and a full bifurcation analysis may be done.
We consider the case of all-to-all communication. Let

$\bm{R}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{r}_{i}$

and

$\mathbf{r}_{i}=\bm{R}+\delta\mathbf{r}_{i},$

where $\delta\mathbf{r}_{i}$ is a fluctuation term with the identity, and

$\sum_{i=1}^{N}\delta\mathbf{r}_{i}=0.$ (6)

Then we can write Eq. 5 as

$\displaystyle\ddot{\bm{R}}+\delta\ddot{\mathbf{r}}_{i}$
$\displaystyle=(1-|\dot{\bm{R}}+\delta\dot{\mathbf{r}}_{i}|^{2})(\dot{\bm{R}}+\delta\dot{\mathbf{r}}_{i})$
$\displaystyle-\frac{a}{N}\sum_{j=1,j\neq
i}^{N}((\bm{R}+\delta\mathbf{r}_{i})-(\bm{R}^{\tau}+\delta\mathbf{r}_{j}^{\tau}))C_{1,i}$
$\displaystyle-\frac{a}{N}\sum_{j=1,j\neq
i}^{N}((\bm{R}+\delta\mathbf{r}_{i})-(\bm{R}+\delta\mathbf{r}_{j}))C_{2,i},$
(7)

where

$\displaystyle C_{1,i}$
$\displaystyle=H(\lVert{\mathbf{r}_{i}-\mathbf{r}_{j}^{\tau}\rVert}-\epsilon)$
$\displaystyle=H(\lVert{(\bm{R}+\delta\mathbf{r}_{i})-(\bm{R}^{\tau}+\delta\mathbf{r}_{j}^{\tau})\rVert}-\epsilon)$
$\displaystyle=H(\lVert{\bm{R}-\bm{R}^{\tau}+\delta\mathbf{r}_{i}-\delta\mathbf{r}_{j}^{\tau}\rVert}-\epsilon)$

and

$C_{2,i}=1-C_{1,i}.$

We use the following to reduce the equations of motion to the mean field: From
Eq. 6, we note

$\displaystyle\sum_{i=1}^{N}\delta\mathbf{r}_{i}^{\tau}$
$\displaystyle=\sum_{j=1,j\neq
i}^{N}\delta\mathbf{r}_{j}^{\tau}+\delta\mathbf{r}_{i}^{\tau}=0\iff$
$\displaystyle-\sum_{j=1,j\neq
i}^{N}\delta\mathbf{r}_{j}^{\tau}=\delta\mathbf{r}_{i}^{\tau}.$ (8)

We further assume that all perturbations from the mean,
$\delta\mathbf{r_{i}}$, are all negligible. (This is always true if the
coupling amplitude is sufficiently large.) In addition, we use the fact that
$\displaystyle{\frac{a(N-1)}{N}}$ limits to $a$, as $N\to\infty$. Therefore,
we obtain mean field approximation for the center of mass of range dependent
coupled delay case:

$\ddot{\bm{R}}=(1-\lvert\dot{\bm{R}}\rvert^{2})\cdot\dot{\bm{R}}-a(\bm{R}-\bm{R}^{\tau})\cdot
H(\lVert{\bm{R}-\bm{R}^{\tau}\rVert}-\epsilon)$ (9)

##  IV Numerical Analysis of the mean field equation

###  IV.1 Examples of rotational attractors

As in the case for the full multi-agent system, we see the existence of
periodic behavior for $\tau$ sufficiently below an instability threshold, as
shown in the time series of Fig. 5. As we increase $\tau$, we expect the
periodic orbit to lose stability, resulting in a new attractor. In particular,
one notices the emergence of a new frequency in addition to the existing
dominant one, as shown in Fig. 6 The additional frequency usually implies a
bifurcation to dynamics on a torus, or a higher dimensional torus.

Figure 5: Periodic motion of the mean field Eq. 9 for
$\epsilon=0.01,a=0.64,\tau=1.6$. (a) Time series of the x-component of the
mean field. (b) Power spectra of the time series.

We now investigate this transition by tracking the stability via monitoring
the Floquet exponents corresponding to the periodic orbit. For a general
differential delay equation given by
$\dot{\bm{x}}(t)=\bm{F}(\bm{x}(t),\bm{x}(t-\tau))$, if
$\bm{\phi}(t)=\bm{\phi}(t+T)$ for all $t\geq 0$, then stability is determined
by examining the linearized equation along $\bm{\phi}(t)$:

$\displaystyle\dot{\bm{X}}(t)$
$\displaystyle=\frac{\partial\bm{F}}{\partial\bm{x}(t)}(\bm{\phi}(t),\bm{\phi}(t-\tau))\bm{X}(t)$
$\displaystyle+\frac{\partial\bm{F}}{\partial\bm{x}(t-\tau)}(\bm{\phi}(t),\bm{\phi}(t-\tau))\bm{X}(t-\tau).$
(10)

The stability of the periodic solution is determined by the spectrum of the
time integration operator $U(T,0)$ which integrates Eq. IV.1 around $\phi(t)$
from time t = 0 to t = T. This operator is called the monodromy operator and
its (infinite number of) eigenvalues, which are independent of the initial
state, are called the Floquet multipliers Hale (1977). For autonomous systems,
it is necessary and sufficient there exists a trivial Floquet multiplier at 1,
corresponding to a perturbation along the periodic solution Hartung  _et al._
(2006); Hale and Lunel (1993). The periodic solution is stable provided all
multipliers (except the trivial one) have modulus smaller than 1; it is
unstable if there exists a multiplier with modulus larger than 1. Bifurcations
occur whenever Floquet multipliers move into or out of the unit circle.
Generically three types of bifurcations occur in a one parameter continuation
of periodic solutions: a turning point, a period doubling, and a torus
bifurcation where a branch of quasi-periodic solutions originates and where a
complex pair of multipliers crosses the unit circle Hale (1977).

Figure 6: Quasi-periodic motion of the mean field Eq. 9. (a) Time series of
the x-component of the mean field. Solid (red) line denotes period length of
dominant spectral peak. Dashed line denotes period length of secondary peak.
(b) Power spectra of the time series.

We have tracked a set of stable periodic orbits for various radii of
$\epsilon$, and located the change in stability by computing the Floquet
multipliers. The results plotted in Fig. 7 show that for a range of radii
$\epsilon$, there exists a bifurcation to a torus at some delay. Notice that
as $\epsilon$ increases, there results an increase in the size of the orbits,
which qualitatively agrees with our full agent based simulations.

Figure 7: Bifurcation plot showing the norm of the periodic orbits as a
function of delay $\tau$. Parameter a=0.68. Red (blue) markers denote unstable
(stable) orbits.Cyan symbols denote the change in stability where a pair of
complex eigenvalues cross the imaginary axis.

Since there exists a range of delays which destabilize periodic swarm dynamics
for each $\epsilon$, we summarize the onset of torus bifurcations by plotting
the locus of points at which stability changes as a a function of coupling
amplitude and delay. The results are plotted in Fig. 8.

Figure 8 is revealing in that it shows a functional relationship of the
bifurcation onset that is similar over a range of $\epsilon$. For larger
values of $\epsilon$, it is clear that lower values of delay and coupling are
required to generate bifurcations. This holds true over two orders of
$\epsilon$. For a fixed value of $\epsilon$, we also see monotonic
relationship between delay and coupling strength, so that it is easier for
smaller delays to destabilize periodic motion for larger coupling strengths.

Figure 8: Plotted is the locus of points at which torus bifurcations emerge as
a function of coupling amplitude $a$, delay $\tau$ for various range radii
$\epsilon$ for the mean field Eq. 9.

##  V Conclusions

We considered a new model of a swarm with delay coupled communication network,
where the delay is considered to be range dependent. That is, given a range
radius, delay is on if two agents are outside the radius, and zero otherwise.
The implication is that small delays do not matter if the agents are close to
each other.

The additional range dependence creates a new set of bifurcations not
previously seen. For general swarms without delay, the usual states consist of
flocking (translation) or ring / rotational state (milling), with agents
spread in phase. With the addition of a fixed delay, a rotational state
bifurcates that has all agents in phase and rotate together Hindes  _et al._
(2020). Range dependence introduces a new rotational bifurcating state that
exhibits behavior observed as a new mixed state combining dynamics of both
ring and rotating states.

The radius parameter $\epsilon$, was used to quantify the bifurcation of the
rotational mixed state. For small $\epsilon$, we see dynamics for the full
swarm shows clustered counter-rotational behavior that is periodic. This
agrees for small radius values in the mean field description as well. As the
radius increases, the mixed periodic state generates new frequencies in the
full model, which are manifested as torus bifurcations in the mean field. Mean
field analysis was done by tracking Floquet multipliers that cross the
imaginary axis as complex pairs. Frequency analysis explicitly shows the
additional frequencies in the mean field.

Finally, we tracked the locus of coupling amplitudes and delay for various
values of $\epsilon$ locating the parameters at which torus bifurcation occur.
The results reveal that as $\epsilon$ increases, torus bifurcations onset at
lower values of coupling amplitude and delay. The implications are that more
complicated behavior than periodic motion has a greater probability of being
observed in both theory and experiment if range dependence of delay is
included.

##  VI Supplementary Material

The videos show the attractor of a swarm consisting of N=300 agents. Fixed
parameters for the three videos are $a=2.0,\tau=1.75$ The parameters for zero
radius (delay is on all the time) are $\epsilon=0.0,c_{r}=0.05$, and
$l_{r}=0.05$ for a baseline, are shown in Video1_eps_0p0.mp4.

The parameters corresponding to Fig. 2 are $\epsilon=0.01,c_{r}=0.01$, and
$l_{r}=0.05$ are shown in Video2_eps_0p01.mp4. The video shows that the
attractor persists when repulsive forces are local and weak. Similar behavior
is observed when N=150, which is used in Fig. 1 without repulsion; i.e.,
$c_{r}=0$. The parameters for corresponding to Fig. 3 are
$\epsilon=0.25,c_{r}=0.05$, and $l_{r}=0.05$, shown in Video3_eps_0p25.mp4.

###### Acknowledgements.

IBS, JH, IT and KK gratefully acknowledge ONR for their support under
N0001412WX20083, N0001420WX00034, and the NRL Base Research Program
N0001420WX00410. VE is supported under the NRL Karles Fellowship Program, JON
55-N2Q4-09. SK was supported through the GMU Provost PhD award as part of the
Industrial Immersion Program. MAH is supported by ONR No. N00014-18-1-2580 and
ARL DCIST CRA W911NF-17-2-0181.

Data sharing is not applicable to this article as no new data were created in
this study.

## References

  * Budrene and Berg (1995) E. O. Budrene and H. C. Berg, “Dynamics of formation of symmetrical patterns by chemotactic bacteria,” Nature 376, 49–53 (1995).
  * Polezhaev  _et al._ (2006) A. A. Polezhaev, R. A. Pashkov, A. I. Lobanov, and I. B. Petrov, “Spatial patterns formed by chemotactic bacteria Escherichia coli,” The International Journal of Developmental Biology 50, 309–314 (2006).
  * Lee  _et al._ (2013) R. M. Lee, D. H. Kelley, K. N. Nordstrom, N. T. Ouellette, and W. Losert, “Quantifying stretching and rearrangement in epithelial sheet migration,” New Journal of Physics 15 (2013), 10.1088/1367-2630/15/2/025036.
  * Tunstrø m  _et al._ (2013) K. r. Tunstrø m, Y. Katz, C. C. Ioannou, C. Huepe, M. J. Lutz, and I. D. Couzin, “Collective states, multistability and transitional behavior in schooling fish,” PLoS computational biology 9, e1002915 (2013).
  * Helbing and Molnar (1995) D. Helbing and P. Molnar, “Social force model for pedestrian dynamics,” Physical Review E 51, 4282–4286 (1995).
  * Giuggioli, McKetterick, and Holderied (2015a) L. Giuggioli, T. J. McKetterick, and M. Holderied, “Delayed Response and Biosonar Perception Explain Movement Coordination in Trawling Bats,” PLOS Computational Biology 11, e1004089 (2015a).
  * Lee (2006) S.-H. Lee, “Predator’s attack-induced phase-like transition in prey flock,” Physics Letters A 357, 270–274 (2006).
  * Topaz and Bertozzi (2004) C. M. Topaz and A. L. Bertozzi, “Swarming Patterns in a Two-Dimensional Kinematic Model for Biological Groups,” SIAM Journal on Applied Mathematics 65, 152–174 (2004).
  * Szwaykowska, Romero, and Schwartz (2015) K. Szwaykowska, L. M.-y.-T. Romero, and I. B. Schwartz, “Collective Motions of Heterogeneous Swarms,” IEEE Transactions on Automation Science and Engineering 12, 810–818 (2015), arXiv:1409.1042 .
  * Mier-y-Teran Romero, Forgoston, and Schwartz (2011) L. Mier-y-Teran Romero, E. Forgoston, and I. B. Schwartz, “Noise, Bifurcations, and Modeling of Interacting Particle Systems,” in  _Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems_ (2011) pp. 3905–3910.
  * Hindes, Szwaykowska, and Schwartz (2016) J. Hindes, K. Szwaykowska, and I. B. Schwartz, “Hybrid dynamics in delay-coupled swarms with mothership networks,” Phys. Rev. E 94, 032306 (2016).
  * Ramachandran, Elamvazhuthi, and Berman (2018) R. K. Ramachandran, K. Elamvazhuthi, and S. Berman, “An optimal control approach to mapping gps-denied environments using a stochastic robotic swarm,” in  _Robotics Research: Volume 1_, edited by A. Bicchi and W. Burgard (Springer International Publishing, Cham, 2018) pp. 477–493.
  * Morgan and Schwartz (2005) D. S. Morgan and I. B. Schwartz, “Dynamic coordinated control laws in multiple agent models,” Physics Letters A 340, 121 – 131 (2005).
  * Wiech, Eremeyev, and Giorgio (2018) J. Wiech, V. A. Eremeyev, and I. Giorgio, “Virtual spring damper method for nonholonomic robotic swarm self-organization and leader following,” Continuum Mechanics and Thermodynamics 30, 1091–1102 (2018).
  * Li  _et al._ (2017) H. Li, C. Feng, H. Ehrhard, Y. Shen, B. Cobos, F. Zhang, K. Elamvazhuthi, S. Berman, M. Haberland, and A. L. Bertozzi, “Decentralized stochastic control of robotic swarm density: Theory, simulation, and experiment,” in  _2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ (2017) pp. 4341–4347.
  * Tanner, Jadbabaie, and Pappas (2007) H. G. Tanner, A. Jadbabaie, and G. J. Pappas, “Flocking in fixed and switching networks,” IEEE Transactions on Automatic Control 52, 863–868 (2007).
  * Gazi (2005) V. Gazi, “Swarm aggregations using artificial potentials and sliding-mode control,” IEEE Transactions on Robotics 21, 1208–1214 (2005).
  * Jadbabaie, Jie Lin, and Morse (2003) A. Jadbabaie, Jie Lin, and A. S. Morse, “Coordination of groups of mobile autonomous agents using nearest neighbor rules,” IEEE Transactions on Automatic Control 48, 988–1001 (2003).
  * Viragh  _et al._ (2014) C. Viragh, G. Vasarhelyi, N. Tarcai, Szorenyi, and et al, “Flocking algorithm for autonomous flying robots,” Bioinspiration & biomimetics 9, 025012 (2014).
  * Desai, Ostrowski, and Kumar (2001) J. P. Desai, J. P. Ostrowski, and V. Kumar, “Modeling and control of formations of nonholonomic mobile robots,” in  _IEEE Transactions on Robotics and Automation_ , Vol. 17(6) (2001) pp. 905–908.
  * Vicsek  _et al._ (2006) T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel type of phase transition in a system of self-driven particles,” (2006), arXiv:0611743v1 [arXiv:cond-mat] .
  * Edelstein-Keshet, Grunbaum, and Watmough (1998) L. Edelstein-Keshet, D. Grunbaum, and J. Watmough, “Do travelling band solutions describe cohesive swarms? An investigation for migratory locusts,” Journal of Mathematical Biology 36, 515–549 (1998).
  * Giuggioli, McKetterick, and Holderied (2015b) L. Giuggioli, T. McKetterick, and M. Holderied, “Delayed response and biosonar perception explain movement coordination in trawling bats,” PLoS Comput Biol 11 (2015b).
  * Nagy  _et al._ (2010) N. Nagy, Z. Akos, D. Biro, and T. Vicsek, “Hierarchical group dynamics in pigeon flocks,” Nature 464, 890–893 (2010).
  * Fehrenbach  _et al._ (2014) J. Fehrenbach, J. Narski, J. Hua, S. Lemercier, A. Jelic, C. Appert-Rolland, S. Donikian, J. Pettré, and P. Degond, “Time-delayed follow-the-leader model for pedestrians walking in line,”  (2014), 10.3934/nhm.2015.10.579, arXiv:1412.7537 .
  * Mier-y-Teran Romero, Forgoston, and Schwartz (2012) L. Mier-y-Teran Romero, E. Forgoston, and I. B. Schwartz, “Coherent Pattern Prediction in Swarms of Delay-Coupled Agents,” IEEE Transactions on Robotics 28, 1034–1044 (2012), arXiv:arXiv:1205.0195v1 .
  * Szwaykowska  _et al._ (2016) K. Szwaykowska, I. B. Schwartz, L. Mier-y Teran Romero, C. R. Heckman, D. Mox, and M. A. Hsieh, “Collective motion patterns of swarms with delay coupling: Theory and experiment,” Phys. Rev. E 93, 032307 (2016).
  * Edwards  _et al._ (2020) V. Edwards, P. deZonia, M. A. Hsieh, J. Hindes, I. Triandaf, and I. B. Schwartz, “Delay-induced swarm pattern bifurcations in mixed-reality experiments,” (2020).
  * Hindes and Schwartz (2018) J. Hindes and I. B. Schwartz, “Rare slips in fluctuating synchronized oscillator networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 071106 (2018), https://doi.org/10.1063/1.5041377 .
  * Szwaykowska, Schwartz, and Carr (2018) K. Szwaykowska, I. B. Schwartz, and T. W. Carr, “State transitions in generic systems with asymmetric noise and communication delay,” in  _11th International Symposium on Mechatronics and its Applications (ISMA)_ (2018) pp. 1–6.
  * Komareji, Shang, and Bouffanais (2018) M. Komareji, Y. Shang, and R. Bouffanais, “Consensus in topologically interacting swarms under communication constraints and time-delays,” Nonlinear Dynamics 93, 1287–1300 (2018).
  * Oliveira, Almeida, and Lima (2015) L. Oliveira, L. Almeida, and P. Lima, “Multi-hop routing within tdma slots for teams of cooperating robots,” in  _2015 IEEE World Conference on Factory Communication Systems (WFCS)_ (2015) pp. 1–8.
  * ying Ani Hsieh  _et al._ (2004) M. ying Ani Hsieh, P. Srivastava, V. Kumar, and C. J. Taylor, “Composable communication constraint-based control,” in  _Mobile Robots XVII_, Vol. 5609, edited by D. W. Gage, International Society for Optics and Photonics (SPIE, 2004) pp. 192 – 200.
  * Hsieh  _et al._ (2007) M. A. Hsieh, A. Cowley, J. F. Keller, L. Chaimowicz, B. Grocholsky, V. Kumar, C. J. Taylor, Y. Endo, R. C. Arkin, B. Jung, D. F. Wolf, G. S. Sukhatme, and D. C. MacKenzie, “Adaptive teams of autonomous aerial and ground robots for situational awareness,” Journal of Field Robotics 24, 991–1014 (2007), https://onlinelibrary.wiley.com/doi/pdf/10.1002/rob.20222 .
  * Fink, Ribeiro, and Kumar (2012) J. Fink, A. Ribeiro, and V. Kumar, “Robust control for mobility and wireless communication in cyber–physical systems with application to robot teams,” Proceedings of the IEEE 100, 164–178 (2012).
  * Fink, Ribeiro, and Kumar (2013) J. Fink, A. Ribeiro, and V. Kumar, “Robust control of mobility and communications in autonomous robot teams,” IEEE Access 1, 290–309 (2013).
  * Arrichiello  _et al._ (2009) F. Arrichiello, D. Liu, S. Yerramall, A. Pereira, J. Das, U. Mitra, and G. Sukhatme, “Effects of underwater communication constraints on the control of marine robot teams,” (2009).
  * Hale (1977) J. K. Hale,  _Theory of Functional Differential Equations_ , Applied Mathematical Sciences (Springer-Verlag, New York, 1977).
  * Hartung  _et al._ (2006) F. Hartung, T. Krisztin, H. Walther, and J. Wu, “Chapter 5 functional differential equations with state-dependent delays: Theory and applications,” in  _Handbook of Differential Equations: Ordinary Differential Equations_, Handbook of Differential Equations: Ordinary Differential Equations, Vol. 3 (2006) pp. 435–545.
  * Hale and Lunel (1993) J. K. Hale and S. M. V. Lunel,  _Introduction to Functional Differential Equations_ (Springer, New York, 1993).
  * Hindes  _et al._ (2020) J. Hindes, V. Edwards, S. Kamimoto, I. Triandaf, and I. B. Schwartz, “Unstable oscillations and bistability in delay-coupled swarms,” (2020), arXiv:2002.12420 [nlin.AO] .


# Evidence for singular-phonon-induced nematic superconductivity in a
topological superconductor candidate Sr0.1Bi2Se3

Jinghui Wang  Kejing Ran  Shichao Li  Zhen Ma  Song Bao  Zhengwei Cai  Youtian
Zhang  National Laboratory of Solid State Microstructures and Department of
Physics, Nanjing University, Nanjing 210093, China  Kenji Nakajima  Seiko
Ohira-Kawamura  J-PARC Center, Japan Atomic Energy Agency, Tokai, Ibaraki
319-1195, Japan  P. Čermák  Jülich Centre for Neutron Science (JCNS) at Heinz
Maier-Leibnitz Zentrum (MLZ), Forschungszentrum Jülich GmbH, Lichtenbergstr.
1, 85748 Garching, Germany  Charles University, Faculty of Mathematics and
Physics, Department of Condensed Matter Physics, Ke Karlovu 5, 121 16, Praha,
Czech Republic  A. Schneidewind  Jülich Centre for Neutron Science (JCNS) at
Heinz Maier-Leibnitz Zentrum (MLZ), Forschungszentrum Jülich GmbH,
Lichtenbergstr. 1, 85748 Garching, Germany  Sergey Y. Savrasov  Department of
Physics, University of California, Davis, California 95616, USA  Xiangang Wan
<EMAIL_ADDRESS>Jinsheng Wen<EMAIL_ADDRESS>National Laboratory of Solid
State Microstructures and Department of Physics, Nanjing University, Nanjing
210093, China  Collaborative Innovation Center of Advanced Microstructures,
Nanjing University, Nanjing 210093, China

###### Abstract

Superconductivity mediated by phonons is typically conventional, exhibiting a
momentum-independent $s$-wave pairing function, due to the isotropic
interactions between electrons and phonons along different crystalline
directions. Here, by performing inelastic neutron scattering measurements on a
superconducting single crystal of Sr0.1Bi2Se3, a prime candidate for realizing
topological superconductivity by doping the topological insulator Bi2Se3, we
find that there exist highly anisotropic phonons, with the linewidths of the
acoustic phonons increasing substantially at long wavelengths, but only for
those along the [001] direction. This observation indicates a large and
singular electron-phonon coupling at small momenta, which we propose to give
rise to the exotic $p$-wave nematic superconducting pairing in the MxBi2Se3
(M=Cu, Sr, Nb) superconductor family. Therefore, we show these superconductors
to be example systems where electron-phonon interaction can induce more exotic
superconducting pairing than the $s$-wave, consistent with the topological
superconductivity.

Introduction
For superconductors, the central issue is what drives the otherwise repulsive
electrons into bound pairs, which collectively condense below the
superconducting transition temperature $T_{\rm{c}}$. In conventional
superconductors, it is known that the elementary excitations of lattice,
phonons, couple with electrons, and act as attractive force that pairs the
electrons, resulting in an isotropic $s$-wave superconducting pairing
symmetry.BCS1957theory  The reason behind this is that the electron-phonon
interaction is often nearly momentum independent. There are some cases where
phonons may play some role in the unconventional superconductivity, for
example, in YBa2Cu3O7-x, but how the electron-phonon coupling is related to
the presumable $d$-wave pairing is not clear at the moment.PhysRevB.38.4992 ;
PhysRevLett.65.915 ; prl359 ; PhysRevLett.70.1457  Therefore, one natural
approach to have a non-$s$-wave superconducting pairing, $d$\- or $p$-wave for
instance, is to resort to other interactions such as magnetic
couplings.RevModPhys.72.969 ; RevModPhys.75.657  Another possible route is to
seek for strongly momentum dependent electron-phonon
coupling,PhysRevLett.88.117005 ; PhysRevB.74.184503 ; PhysRevB.77.094502 ;
0953-8984-20-4-045227  but extensive research in three-dimensional materials
with weak spin-orbit coupling (SOC) has not been quite successful so far.

Figure 1: Schematic crystal structure and the experimental scheme of the
neutron scattering measurements. a, Schematic of the hexagonal crystal
structure of Sr0.1Bi2Se3. The dopant Sr ions are believed to be intercalated
in the van der Waals gaps between the quintuple layers consisting of Bi and Se
ions.du2017drive  The dashed lines denote a unit cell, while the shades denote
the octahedra formed by the Se ions, with Bi ions in the centre. b, Brillouin
zones in the ($H0L$) plane defined by two orthogonal axes [100] and [001] in
the reciprocal space, constructed based on the hexagonal notation illustrated
in a. The two long arrows point to (0, 0, -15) and (1, 0, -5), the two
strongest Bragg peaks available around which we measure the phonons. Arrows on
(0, 0, -15) and (1, 0, -5) represent the directions along which we plot the
phonon dispersions in Fig. 2. The letters a-d correspond to panels a-d in Fig.
2. Dashed lines illustrate the Brillouin zone boundaries, and the shades
represent the two measured Brillouin zones. Dots represent high-symmetry
points, with $\Gamma$, Z, F, and L in the measured Brillouin zones being
labelled. The wave vector Q is expressed as ($H,\,K,\,L$) reciprocal lattice
unit (rlu) of
$(a^{*},\,b^{*},\,c^{*})=(4\pi/\sqrt{3}a,\,4\pi/\sqrt{3}b,\,2\pi/c)$, with
$a=b=4.14$ Å, and $c=28.6$ Å. Figure 2: Phonon dispersions around two Bragg
peaks (0, 0, -15) and (1, 0, -5) measured at $T=17$ K. a and b, phonon
dispersing from (0, 0, -15) along the [100] and [001] directions,
respectively. c and d, phonon dispersing from (1, 0, -5) along the [100] and
[001] directions, respectively. The data plotted in a and b were collected on
a time-of-flight (TOF) spectrometer AMATERAS, and those in c and d were
collected on a triple-axis spectrometer (TAS) PANDA. For the TOF data plotted
in a and b, which are the phonon dispersions along one direction, one needs to
sum over intensities for a certain thickness along the other two directions to
improve the statistics. The data plotted against the $H$ direction in a were
obtained by integrating the intensities in ( $q$, 0+$K$, -15+$L$) with a
thickness of $K$ and $L$ ranging from -0.1 to 0.1 rlu, and -16 to -14 rlu,
respectively. Those in b were integrated with $H=K=[-0.05,0.05]$ rlu. $E$ and
$q$ represent the phonon energy and momentum, respectively. TA, LA, and TO
represent transverse acoustic, longitudinal acoustic and transverse optic
modes, respectively. Solid lines are the dispersions obtained by fitting the
energy scans at a series of $q$s as shown in Fig. 3. The white streaks along
the diagonal direction in b were due to the lack of detector coverage wherein.
Figure 3: Energy cuts on the phonon dispersions. a-d correspond to the energy
cuts on the dispersions plotted in Fig. 2a-d at different $q$ values. The cuts
are offset so that each cut can be visualized clearly. In c, due to the large
intensity difference of phonons at small and large $q$s, we plot the square
root of the intensities as a function of energy so that the cut profiles at
large $q$s can be visible. Solid lines are fits with Lorentzian functions
convoluted the instrumental resolutions, as indicated by the horizontal bars.
From the fits, we obtain the phonon dispersions, which are plotted as solid
lines in Fig. 2. Errors represent one standard deviation throughout the paper.
Figure 4: Evolution of the phonon linewidth with $q$. a and b, Phonons along
the [100] and [001] directions, respectively. Squares and triangles represent
results extracted from data collected on a TOF spectrometer AMATERAS and TAS
PANDA, respectively. For comparison, the linewidths of the Bi2Se3 phonons
measured on AMATERAS are also plotted. Lines through data are guides to eye to
illustrate the evolution of the phonon linewidth as a function of $q$.

Recently, the debate on the possibility of realizing topological-insulator-
derived topological superconductivity where the electron-electron correlation
is weak and the electron-phonon interaction dominates the superconducting
pairing, has made this fundamental question more
outstanding.RevModPhys.83.1057 ; RevModPhys.82.3045  In this regard, a
promising topological superconductor candidate Sr0.1Bi2Se3 with quasi-two-
dimensional structure and strong SOC has come to our attention. Sr0.1Bi2Se3 is
a member of the MxBi2Se3 (M=Cu, Sr, Nb) family, which become superconducting
by doping the topological insulator Bi2Se3.np5_398 ; np5_438 ;
PhysRevLett.104.057001 ; PhysRevB.84.054513 ; PhysRevLett.105.097001 ; np6_855
; PhysRevLett.107.217001 ; doi:10.1021/jacs.5b06815 ; PhysRevB.92.020506 ;
arXiv:1512.03519 ; :/content/aip/journal/apl/107/17/10.1063/1.4934590 ;
PhysRevMaterials.2.014201 ; np12_852 ; np13_123 ; PhysRevB.98.184509 ;
PhysRevB.94.180510 ; sr6_28632 ; Du2017 ; PhysRevX.7.011009 ;
PhysRevLett.105.097001 ; PhysRevB.90.100509 ; PhysRevB.94.180504 ;
PhysRevLett.113.046401 ; nc5_5144 ; PhysRevB.90.184512 ; PhysRevB.96.144504
Bi2Se3 crystallises in a layered hexagonal structure with the space group
$R\bar{3}m$ (Fig. 1a). Superconductivity in this family was initially
discovered in CuxBi2Se3 by intercalating Cu into the gaps between the Se
layers of Bi2Se3,PhysRevLett.104.057001 ; PhysRevB.84.054513  and subsequently
in SrxBi2Se3 and NbxBi2Se3, with the latter two having larger superconducting
volume fractions.doi:10.1021/jacs.5b06815 ; PhysRevB.92.020506 ;
arXiv:1512.03519  Experimental studies from nuclear magnetic
resonance,np12_852  specific heat,np13_123 ; PhysRevB.98.184509  penetration
depth,PhysRevB.94.180510  transport,sr6_28632 ; Du2017  and magnetic
torque,PhysRevX.7.011009  all indicate that the superconducting order
parameter breaks the in-plane three-fold crystalline rotation symmetry, and
exhibits a two-fold nematic pattern consistent with that of a $p$-wave
pairing.

A number of theoretical works have also proposed the exotic $p$-wave
superconductivity in this family.PhysRevLett.105.097001 ; nc5_5144 ;
PhysRevB.90.100509 ; PhysRevB.90.184512 ; PhysRevLett.113.046401 ;
PhysRevB.94.180504 ; PhysRevB.96.144504  Based on a short-range electron
density-density interaction in a two-orbital model, Fu and Berg proposed a
spin-triplet pairing with odd parity in Cu-doped Bi2Se3.PhysRevLett.105.097001
The presence of strong SOC and hexagonal warping in Bi2Se3 had been suggested
to play an important role in rendering the unconventional superconducting
behaviours.PhysRevB.90.100509 ; PhysRevB.94.180504  By performing symmetry
analyses, Yang et al. also pointed out the possibility of realizing a $p$-wave
superconductivity in CuxBi2Se3.PhysRevLett.113.046401  For such a system with
$p$ electrons where electronic correlations are weak, it is most likely that
the electron-phonon interaction will be the main driving force for the
superconductivity. Theoretically, it is indeed proposed that such an
non-$s$-wave superconducting pairing in CuxBi2Se3 is driven by the electron-
phonon interaction, which is usually expected to result in a conventional
$s$-wave superconductivity.nc5_5144 ; PhysRevB.90.184512  In particular, based
on first-principles calculations, a singular electron-phonon interaction at
long wavelength had been suggested in CuxBi2Se3.nc5_5144  This unusual
electron-phonon interaction was shown to make different pairing channels have
similar strength, and the $A_{\rm 2u}$ symmetric $p_{z}$ pairing should win if
the effect of the Coulomb pseudopotential $\mu^{*}$ was taken into
account.nc5_5144  Similarly, Brydon et al. also proposed an odd-parity
superconductivity in the same material by considering the combined effect of
the electron-phonon interaction and $\mu^{*}$.PhysRevB.90.184512  Based on an
effective model for Bi2Se3, the on-site repulsive interaction had been
carefully studied, and the results showed that the Coulomb interaction
generated repulsive interaction in both the s-wave and $A_{\rm 2u}$ paring
channels, but not in the $E_{\rm u}$ ($p_{x}$ or $p_{y}$ pairing) channel,
making the latter win over the former two.PhysRevB.96.144504  Therefore,
examining the electron-phonon interaction experimentally is the key to the
understanding of the superconducting pairing mechanism in the MxBi2Se3 (M=Cu,
Sr, Nb) superconductor family.

Here, using inelastic neutron scattering (INS), we measure the phonons on a
superconducting single crystal of Sr0.1Bi2Se3, isostructural to CuxBi2Se3 but
with a larger superconducting volume fraction.doi:10.1021/jacs.5b06815 ;
PhysRevB.92.020506  We show that the phonon linewidths of the acoustic mode
measured along the [001] direction increase significantly when approaching the
Brillouin zone centre. This observation reflects the singular electron-phonon
pairing interaction which diverges at small momenta ($q$s), and may be
responsible for the exotic $p$-wave superconducting pairing symmetry discussed
extensively.PhysRevB.94.180510 ; sr6_28632 ; np12_852 ; np13_123 ; Du2017 ;
PhysRevX.7.011009 ; PhysRevLett.105.097001 ; nc5_5144 ; PhysRevB.90.100509 ;
PhysRevB.90.184512 ; PhysRevLett.113.046401 ; PhysRevB.94.180504 ;
PhysRevB.96.144504

Results
Phonon dispersions. The Sr0.1Bi2Se3 single crystals with a $T_{\rm{c}}=3.2$ K
used in this work have been well characterised as described in the
Supplementary Fig. 1 and in ref. du2017drive, . We have carried out INS
measurements on the single crystal using an experimental scheme sketched in
Fig. 1b. We focus on two accessibly strongest Bragg peaks (0, 0, -15) and (1,
0, -5) and map out the phonons at low energies around these peaks. At (0, 0,
-15), as the wave vector Q is parallel to the [001] direction, phonons
propagating along [001] and [100] directions are purely longitudinal and
transverse modes, respectively. On the other hand, at (1, 0 -5), along both
[001] and [100] directions, there is a mix of the longitudinal and transverse
phonon modes, as these directions are in between the longitudinal and
transverse directions, as illustrated in Fig. 1b.

In Fig. 2, we show the phonon dispersions around (0, 0, -15) and (1, 0, -5)
for Sr0.1Bi2Se3 measured at $T=17$ K. These low-energy branches correspond to
the motions of the heavier Bi atoms.nc5_5144  Around (0, 0, -15), we observe
one dispersing transverse acoustic (TA) mode up to $\sim$5 meV along the [100]
direction (Fig. 2a); along the [001] direction, a softer longitudinal acoustic
(LA) mode is observed with a band top of $\sim$3 meV (Fig. 2b). Around (1, 0,
-5), we find one TA and one LA mode along the [100] direction (Fig. 2c); along
the [001] direction, besides the softer TA and LA modes, we also observe an
almost dispersionless transverse optic (TO) mode (Fig. 2d). The TO mode at
$q=0.4$ rlu has an energy of $\sim$4.85 meV, which is very close to the zone
centre optic mode at 4.82 meV obtained from Raman
scattering.PhysRevB.84.195118  Overall, these data are in reasonable agreement
with previous numerical resultsnc5_5144 ; PhysRevB.84.195118 ; apl100_082109 ;
PhysRevB.83.094301 , with phonons somewhat softer than those calculated in
ref. nc5_5144, .

We have also examined the phonons at $T=0.8$ K, below the $T_{\rm{c}}$ of 3.2
K, but do not observe any essential changes on the phonon spectra across the
$T_{\rm{c}}$. Therefore, in the rest of the paper, we will discuss our results
obtained at $T=17$ K. We note that for some superconductors with a higher
$T_{\rm{c}}$, there have been some reports showing changes for certain phonon
modes matching twice of the superconducting gap $2\Delta(T)$ across
$T_{\rm{c}}$,PhysRevLett.30.214 ; PhysRevB.12.4899 ; PhysRevLett.101.237002 ;
PhysRevB.82.024509  likely due to the participation in the electron-phonon
coupling for these modes.PhysRevB.56.5552

Linewidth broadenings near the zone centre. We perform energy cuts at a series
of $q$s on the phonon dispersions shown in Fig. 2, and the results of the cuts
are plotted in Fig. 3. From the cuts, we see that the phonon linewidth at
different $q$s along the [100] direction for the two TA modes around ($q$, 0,
-15) and (1+$q$, 0, -5), shown in Fig. 3a and c, respectively, do not show
noticeable changes. Similarly, as shown in Fig. 3d, around (1, 0, -5), along
the [001] direction, the linewidth of the LA and TO modes has no $q$
dependence. On the other hand, for the LA and TA phonons along the [001]
direction, shown in Fig. 3b and d, respectively, we do clearly see that they
become much broader when approaching the zone centre. For instance, in Fig.
3b, the phonon is sharp and almost resolution limited at large $q$s such as at
$q=1.0$-$1.3$ rlu. When $q=0.5$ rlu, the phonon becomes very broad. As can be
clearly seen from Fig. 3d, the TA mode becomes significantly broader for
$q\lesssim 1.0$ rlu, well before the energy of the LA mode getting too close
to be resolved from the TA mode. In other words, the linewidth increase
observed in Fig. 3d is not due to the merge of the TA and LA modes.

To better characterise the phonon linewidth, we have fitted the energy scans
in Fig. 3 using Lorentzian functions convoluted with the instrumental
resolution, and plotted the fitted width as a function of $q$ in Fig. 4. The
results show that along the [100] direction, the linewidth remains almost
constant for the TA modes around (0, 0, -15) and (1, 0, -5) (Fig. 4a); by
contrast, as shown in Fig. 4b, along the [001] direction, the width increases
from 0.15 meV at $q=1.3$ rlu to 1.0 meV at $q=0.5$ rlu. For comparison
purpose, we have performed similar measurements on the phonons of the undoped
compound Bi2Se3, and the obtained phonon dispersions are shown in
Supplementary Fig. 2. In Fig. 4, we also plot the phonon linewidth as a
function of $q$ for Bi2Se3. The results are very similar to those of
Sr0.1Bi2Se3, where the low-energy phonons only become broader along the [001]
direction at small $q$s. However, compared to Sr0.1Bi2Se3, the linewidth
increase is much reduced, indicating that Sr doping can enhance the broadening
effect at small $q$s. The evolution of the phonon linewidth has actually been
predicted in an isostructural compound CuxBi2Se3 in ref. nc5_5144, , where it
is shown that the acoustic phonons remain well defined at large $q$s along the
[001] direction or in the whole $q$ range along other directions. In other
words, according to the calculations,nc5_5144  the acoustic phonons will only
show broadenings at small $q$s along the [001] direction, which is exactly
what we show in Fig. 4.

Discussions
By far, we have demonstrated that the acoustic phonons become broad only along
the [001] direction at small $q$s, consistent with the theoretical
calculations.nc5_5144  What is the origin of such singular phonons and what is
the consequence? Since Bi2Se3 has both time-reversal and spatial-inversion
symmetry, every energy band is at least double degenerate. With the strong
SOC, the phonon displacement along the [001] direction at small $q$s which
breaks the spatial inversion symmetry efficiently lifts the double degeneracy,
resulting in a large and singular electron-phonon coupling matrix element
along this direction near the zone centre. Furthermore, because of an open-
cylinder-like electron pocket along the $\Gamma$-$Z$ direction centering at
the $\Gamma$ point,PhysRevB.88.195107  it is shown that there is a Fermi
surface nesting along this direction, and the nesting function has the largest
value as $q\rightarrow 0$.nc5_5144  The combining effect of the large and
singular electron-phonon interaction as well as the strong Fermi surface
nesting gives rise to the large linewidth for the phonons along the [001]
direction at small $q$s. The contribution of each phonon mode to the electron-
phonon coupling constant is proportional to the phonon linewidth divided by
the square of the phonon energy. Therefore, the broad acoustic phonons at
small $q$s observed here should dominate the electron-phonon coupling. As
shown in Fig. 4b, the linewidth increase is larger for Sr0.1Bi2Se3 than that
for the Bi2Se3 sample. In this sense, the Sr doping not only brings in the
electrons,du2017drive ; PhysRevLett.104.057001  but also enhances the
electron-phonon coupling, both of which drive the superconductivity in
Sr0.1Bi2Se3.

Based on a symmetry analysis, one can expect that this highly unusual
electron-phonon interaction results in a remarkable proximity of all pairing
channels. The strengths of the $s$\- and $p$-wave pairing interaction are
indeed shown to be comparable in this case.nc5_5144  Naively one may expect
that the Coulomb interaction $\mu^{\ast}$ only suppresses the $s$-wave pairing
and a moderate $\mu^{\ast}$ should result in an $A_{\rm 2u}$ ($p_{z}$)
symmetric pairing.nc5_5144  However, a realistic calculation reveals that the
Coulomb interaction in Bi2Se3 generates repulsive interaction in both the
$s$-wave and $A_{\rm 2u}$ pairing channels but not in the $E_{\rm u}$ ($p_{x}$
or $p_{y}$) channel.PhysRevB.96.144504  Considering that the strength of the
$E_{\rm u}$ pairing is comparable to that of the $s$-wave pairing from first-
principles linear response calculations, one can thus expect a two-fold
nematic $E_{\rm u}$ pairing which breaks the in-plane three-fold crystalline
rotation symmetry in this system.nc5_5144 ; PhysRevB.96.144504  Alternatively,
Fu and Berg suggest that the $E_{\rm u}$ pairing states can also be realized
as a consequence of the inter-orbital pairing, and the broadening of the [001]
phonons observed here may assist such inter-orbital
pairing.PhysRevLett.105.097001

As we show in Supplementary Fig. 3, for the superconducting sample Sr0.1Bi2Se3
we study here, both the resistance and upper critical field exhibit a clear
two-fold pattern, consistent with previous experimental reports on the nematic
$p$-wave superconductivity in the MxBi2Se3 (M=Cu, Sr, Nb) superconductor
family.np12_852 ; np13_123 ; PhysRevB.94.180510 ; sr6_28632 ; Du2017 ;
PhysRevX.7.011009  Our work thus provides concrete evidence that an odd-parity
$p$-wave pairing compatible with that of the topological superconductivity can
be induced by strongly anisotropic electron-phonon interaction. Very recently,
thousands of topological materials based on weakly-correlated systems have
been proposed theoretically.nature566_475 ; nature566_480 ; nature566_486
Since we show here that topological superconductivity can exist in such
systems, the search for topological superconductors should be much more
promising.

Methods
Single-crystal growth and characterisations. High-quality single crystals of
Bi2Se3 and Sr0.1Bi2Se3 were grown using the horizontal Bridgman method. The
stoichiometric raw materials (99.99% Sr, 99.99% Bi and 99.99% Se powders) of
nominal composition were mixed well and sealed in an ampoule in vacuum. The
ampoule was then placed in a tube furnace with a temperature gradient, heated
up to 850 ∘C, stayed for 2 days, and cooled from 850 ∘C to 610 ∘C at 1 ∘C per
hour for the melt to crystallise. The resistance and magnetization were
measured in a physical property measurement system (PPMS-9T, Quantum Design).
The PPMS was equipped with a sample rotator used to measure the angular
dependence of the electrical properties.

INS experiments. INS measurements were performed on the cold triple-axis
spectrometer (TAS) PANDA located at FRM-II, Germany, and cold-neutron time-of-
flight (TOF) disk-chopper spectrometer AMATERAS located at J-PARC, Tokai,
Ibaraki, Japan. On the TAS, we used a fixed final energy $E_{\rm{f}}$ mode.
For the experiment on PANDA, one single crystal of Sr0.1Bi2Se3 with a mass of
1.1 g was mounted in the ($H,\,0,\,L$) scattering plane. The sample mosaic was
about 5 degrees (FWHM, full width at half maximum) scanned around the Brag
peak (1, 0, -5). Double-focusing mode was used for both the monochromator and
analyser. The collimation is 80′-80′-sample-80′-80′. With these
configurations, the energy resolutions were about 0.09 meV for
$E_{\rm{f}}=5.1$ meV and 0.05 meV for $E_{\rm{f}}=3.5$ meV. To reduce higher-
order neutrons, one Be filter was placed after the sample. A closed-cycle
refrigerator (CCR) equipped with a 3He insert was used to reach the base
temperature of 0.8 K. In the TAS experiments, we measured the phonon
dispersion along the direction that matched the resolution ellipse to improve
the resolution. On the TOF spectrometer AMATERAS, a multiple-$E_{\rm{i}}$
(incident energy) mode with main $E_{\rm{i}}=10.5$ meV and resolution $\Delta
E=0.33$ meV at $E=0$ meV was used. For Sr0.1Bi2Se3, the 1.1-g sample measured
on PANDA was used on AMATERAS. We also measured a Bi2Se3 sample with a larger
mass of 1.6 g for comparison. Both samples were mounted with $c$ axis along
the beam direction and $a$ axis in the horizontal plane. The rotation axis was
perpendicular to $a$-$c$ plane, and the step was 0.5 degree around the Bragg
peaks (0, 0, -15) and (1, 0, -5). A 4He CCR with a temperature range of 5-300
K was used on AMATERAS. The wave vector Q is expressed as ($H,\,K,\,L$)
reciprocal lattice unit (rlu) of
$(a^{*},\,b^{*},\,c^{*})=(4\pi/\sqrt{3}a,\,4\pi/\sqrt{3}b,\,2\pi/c)$, with
$a=b=4.14$ Å and $c=28.6$ Å for both Bi2Se3 and Sr0.1Bi2Se3.

Data availability: Data are available from J.S.W. (Email<EMAIL_ADDRESS>upon reasonable request. The source data underlying Figs 2-4 and Supplementary
Figs 1-3 are provided as a Source Data file.

## References

  * (1) Bardeen, J., Cooper, L. N. & Schrieffer, J. R.  Theory of superconductivity.  _Phys. Rev._ 108, 1175–1204 (1957).
  * (2) Krantz, M., Rosen, H. J., Macfarlane, R. M. & Lee, V. Y.  Effect of oxygen stoichiometry on softening of Raman active lattice modes in YBa2Cu3Ox.  _Phys. Rev. B_ 38, 4992–4995 (1988).
  * (3) Friedl, B., Thomsen, C. & Cardona, M.  Determination of the superconducting gap in RBa2Cu3O7-δ.  _Phys. Rev. Lett._ 65, 915–918 (1990).
  * (4) Arai, M. _et al._ Anomaly of phonon state of superconducting YBa2Cu3O7 studied by inelastic neutron scattering.  _Phys. Rev. Lett._ 69, 359 (1992).
  * (5) Pyka, N. _et al._ Superconductivity-induced phonon softening in YBa2Cu3O7 observed by inelastic neutron scattering.  _Phys. Rev. Lett._ 70, 1457–1460 (1993).
  * (6) Tsuei, C. C. & Kirtley, J. R.  Pairing symmetry in cuprate superconductors.  _Rev. Mod. Phys._ 72, 969–1016 (2000).
  * (7) Mackenzie, A. P. & Maeno, Y.  The superconductivity of ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$ and the physics of spin-triplet pairing.  _Rev. Mod. Phys._ 75, 657–712 (2003).
  * (8) Varelogiannis, G.  Small- $\mathit{q}$ phonon-mediated superconductivity in organic $\kappa$-BEDT-TTF compounds.  _Phys. Rev. Lett._ 88, 117005 (2002).
  * (9) Schnell, I., Mazin, I. I. & Liu, A. Y.  Unconventional superconducting pairing symmetry induced by phonons.  _Phys. Rev. B_ 74, 184503 (2006).
  * (10) Alexandrov, A. S.  Unconventional pairing symmetry of layered superconductors caused by acoustic phonons.  _Phys. Rev. B_ 77, 094502 (2008).
  * (11) Ye, X.-S., Yao, Z.-J. & Li, J.-X.  Phonon-mediated unconventional superconductivity in NaxCoO${}_{2}\centerdot y$H2O compounds.  _J. Phys. Conden. Matter_ 20, 045227 (2008).
  * (12) Du, G. _et al._ Drive the Dirac electrons into Cooper pairs in SrxBi2Se3.  _Nat. Commun._ 8, 14466 (2017).
  * (13) Qi, X.-L. & Zhang, S.-C.  Topological insulators and superconductors.  _Rev. Mod. Phys._ 83, 1057–1110 (2011).
  * (14) Hasan, M. Z. & Kane, C. L.  Colloquium: Topological insulators.  _Rev. Mod. Phys._ 82, 3045–3067 (2010).
  * (15) Xia, Y. _et al._ Observation of a large-gap topological-insulator class with a single Dirac cone on the surface.  _Nat. Phys._ 5, 398 (2009).
  * (16) Zhang, H. _et al._ Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface.  _Nat. Phys._ 5, 438–442 (2009).
  * (17) Hor, Y. S. _et al._ Superconductivity in ${\mathrm{Cu}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$ and its implications for pairing in the undoped topological insulator.  _Phys. Rev. Lett._ 104, 057001 (2010).
  * (18) Kriener, M. _et al._ Electrochemical synthesis and superconducting phase diagram of CuxBi2Se3.  _Phys. Rev. B_ 84, 054513 (2011).
  * (19) Fu, L. & Berg, E.  Odd-parity topological superconductors: Theory and application to ${\mathrm{Cu}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. Lett._ 105, 097001 (2010).
  * (20) Wray, L. A. _et al._ Observation of topological order in a superconducting doped topological insulator.  _Nat. Phys._ 6, 855 (2010).
  * (21) Sasaki, S. _et al._ Topological superconductivity in ${\mathrm{Cu}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. Lett._ 107, 217001 (2011).
  * (22) Liu, Z. _et al._ Superconductivity with topological surface state in SrxBi2Se3.  _J. Am. Chem. Soc._ 137, 10512–10515 (2015).
  * (23) Shruti, Maurya, V. K., Neha, P., Srivastava, P. & Patnaik, S.  Superconductivity by Sr intercalation in the layered topological insulator ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. B_ 92, 020506 (2015).
  * (24) Qiu, Y. _et al._ Time reversal symmetry breaking superconductivity in topological materials.  _Preprint at https://arxiv.org/abs/1512.03519_ (2015).
  * (25) Han, C. Q. _et al._ Electronic structure of a superconducting topological insulator Sr-doped Bi2Se3.  _Appl. Phys. Lett._ 107 (2015).
  * (26) Li, Z. _et al._ Possible structural origin of superconductivity in Sr-doped $\mathrm{B}{\mathrm{i}}_{2}\mathrm{S}{\mathrm{e}}_{3}$.  _Phys. Rev. Mater._ 2, 014201 (2018).
  * (27) Matano, K., Kriener, M., Segawa, K., Ando, Y. & Zheng, G.-q.  Spin-rotation symmetry breaking in the superconducting state of CuxBi2Se3.  _Nat. Phys._ 12, 852–854 (2016).
  * (28) Yonezawa, S. _et al._ Thermodynamic evidence for nematic superconductivity in CuxBi2Se3.  _Nat. Phys._ 13, 123–126 (2017).
  * (29) Willa, K. _et al._ Nanocalorimetric evidence for nematic superconductivity in the doped topological insulator ${\mathrm{Sr}}_{0.1}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. B_ 98, 184509 (2018).
  * (30) Smylie, M. P. _et al._ Evidence of nodes in the order parameter of the superconducting doped topological insulator ${\mathrm{Nb}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$ via penetration depth measurements.  _Phys. Rev. B_ 94, 180510 (2016).
  * (31) Pan, Y. _et al._ Rotational symmetry breaking in the topological superconductor SrxBi2Se3 probed by upper-critical field experiments.  _Sci. Rep._ 6, 28632 (2016).
  * (32) Du, G. _et al._ Superconductivity with two-fold symmetry in topological superconductor SrxBi2Se3.  _Sci. Chin. Phys. Mech. & Astro._ 60, 037411 (2017).
  * (33) Asaba, T. _et al._ Rotational symmetry breaking in a trigonal superconductor Nb-doped ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. X_ 7, 011009 (2017).
  * (34) Fu, L.  Odd-parity topological superconductor with nematic order: Application to ${\mathrm{Cu}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. B_ 90, 100509 (2014).
  * (35) Venderbos, J. W. F., Kozii, V. & Fu, L.  Odd-parity superconductors with two-component order parameters: Nematic and chiral, full gap, and Majorana node.  _Phys. Rev. B_ 94, 180504 (2016).
  * (36) Yang, S. A., Pan, H. & Zhang, F.  Dirac and Weyl superconductors in three dimensions.  _Phys. Rev. Lett._ 113, 046401 (2014).
  * (37) Wan, X. & Savrasov, S. Y.  Turning a band insulator into an exotic superconductor.  _Nat. Commun._ 5, 4144 (2014).
  * (38) Brydon, P. M. R., Das Sarma, S., Hui, H.-Y. & Sau, J. D.  Odd-parity superconductivity from phonon-mediated pairing: Application to ${\mathrm{Cu}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$.  _Phys. Rev. B_ 90, 184512 (2014).
  * (39) Wu, F. & Martin, I.  Nematic and chiral superconductivity induced by odd-parity fluctuations.  _Phys. Rev. B_ 96, 144504 (2017).
  * (40) Gnezdilov, V. _et al._ Helical fluctuations in the Raman response of the topological insulator Bi2Se3.  _Phys. Rev. B_ 84, 195118 (2011).
  * (41) Wang, B.-T. & Zhang, P.  Phonon spectrum and bonding properties of Bi2Se3: Role of strong spin-orbit interaction.  _Appl. Phys. Lett._ 100, 082109 (2012).
  * (42) Cheng, W. & Ren, S.-F.  Phonons of single quintuple Bi2Te3 and Bi2Se3 films and bulk materials.  _Phys. Rev. B_ 83, 094301 (2011).
  * (43) Axe, J. D. & Shirane, G.  Influence of the superconducting energy gap on phonon linewidths in Nb3Sn.  _Phys. Rev. Lett._ 30, 214–216 (1973).
  * (44) Shapiro, S. M., Shirane, G. & Axe, J. D.  Measurements of the electron-phonon interaction in Nb by inelastic neutron scattering.  _Phys. Rev. B_ 12, 4899–4908 (1975).
  * (45) Weber, F. _et al._ Direct observation of the superconducting gap in phonon spectra.  _Phys. Rev. Lett._ 101, 237002 (2008).
  * (46) Weber, F. & Pintschovius, L.  Superconductivity-induced distortions of phonon lineshapes in niobium.  _Phys. Rev. B_ 82, 024509 (2010).
  * (47) Allen, P. B., Kostur, V. N., Takesue, N. & Shirane, G.  Neutron-scattering profile of $q\neq 0$ phonons in BCS superconductors.  _Phys. Rev. B_ 56, 5552–5558 (1997).
  * (48) Lahoud, E. _et al._ Evolution of the Fermi surface of a doped topological insulator with carrier concentration.  _Phys. Rev. B_ 88, 195107 (2013).
  * (49) Zhang, T. _et al._ Catalogue of topological electronic materials.  _Nature_ 566, 475–479 (2019).
  * (50) Vergniory, M. G. _et al._ A complete catalogue of high-quality topological materials.  _Nature_ 566, 480–485 (2019).
  * (51) Tang, F., Po, H. C., Vishwanath, A. & Wan, X.  Comprehensive search for topological materials using symmetry indicators.  _Nature_ 566, 486–489 (2019).

Acknowledgements
The work was supported by the National Natural Science Foundation of China
with Grant Nos 11822405, 11674157, 11525417, and 11834006, NSF Jiangsu
province in China with Grant No. BK20180006, Fundamental Research Funds for
the Central Universities with Grant No. 020414380117, and the Office of
International Cooperation and Exchanges of Nanjing University. S.Y.S. was
supported by NSF DMR with Grant No. 1832728. We thank Shengyuan Yang, Fan
Zhang, Haijun Zhang, Hai-Hu Wen, and Qiang-Hua Wang for stimulating
discussions.

Author contributions
J.S.W. and X.G.W. conceived the project. J.H.W. grew and characterised the
single crystals with assistance from K.J.R., S.C.L., Z.M., S.B., Z.W.C., and
Y.T.Z.. J.S.W. and J.H.W. carried out the neutron scattering experiments with
help from K.N., S.O.K., P.C., and A.S.. J.H.W. and J.S.W. analysed the data.
J.S.W., X.G.W., and J.H.W. wrote the paper with inputs and comments from
S.Y.S. and all other authors.

Correspondence
Correspondence and request for materials should be addressed to J.S.W. (Email:
<EMAIL_ADDRESS>and X.G.W. (Email: xgwan@nju.edu.cn).

Competing Interests
The authors declare no competing interests.