# Topological phonons in arrays of ultracold dipolar particles

Marco Di Liberto<EMAIL_ADDRESS>Andreas Kruckenhauser  Peter
Zoller  Mikhail A. Baranov  Institute for Quantum Optics and Quantum
Information of the Austrian Academy of Sciences, 6020 Innsbruck, Austria
Institute for Theoretical Physics, University of Innsbruck, 6020 Innsbruck,
Austria

###### Abstract

The notion of topology in physical systems is associated with the existence of
a nonlocal ordering that is insensitive to a large class of perturbations.
This brings robustness to the behaviour of the system and can serve as a
ground for developing new fault-tolerant applications. We discuss how to
design and study a large variety of topology-related phenomena for phonon-like
collective modes in arrays of ultracold polarized dipolar particles. These
modes are coherently propagating vibrational excitations, corresponding to
oscillations of particles around their equilibrium positions, which exist in
the regime where long-range interactions dominate over single-particle motion.
We demonstrate that such systems offer a distinct and versatile tool to
investigatea wide range of topological effects in a _single experimental
setup_ with a chosen underlying crystal structure by simply controlling the
anisotropy of the interactions via the orientation of the external polarizing
field. Our results show that arrays of dipolar particles provide a promising
unifying platform to investigate topological phenomena with phononic modes.

##  1 Introduction

Many of the intriguing phenomena in condensed matter systems have been
associated and tied to the notion of topology and to the existence of
topological states characterized by discrete topological indices (see, for
example, Refs. [1, 2, 3, 4]). The most prominent manifestations of the
underlying nontrivial topology in a many-body system are robust edge states
and quantized responses to external fields, as for the celebrated Quantum Hall
effect [5, 6]. The family of physical systems exhibiting topological phases
and effects is wide and includes topological insulators and superconductors
[7], crystalline insulators [8], ultracold atoms and molecules [9, 10],
photonic systems [11], electric circuits [12], and also acoustic and
mechanical systems [13, 14].

Ultracold dipolar particles with long-range and anisotropic interparticle
interactions provide a versatile and powerful platform for simulating many-
body quantum systems, including topological ones [15, 16, 17, 18]. Recent
progress in creating systems of dipolar particles in periodic arrays, as for
the case of magnetic atoms [19, 20], polar molecules [21] or Rydberg atoms
[22], offers a remarkable opportunity to experimentally implement and study
various topological lattice models in this setting. From the theory side, it
has been suggested to use dipolar particles to implement topological bands
[23, 24], Hopf insulators [25], chiral excitations [26], Weyl nodes [27],
topological bound states [28], or strongly-correlated many-body systems [29,
30, 31]. From the experimental side, the realizations of artificial gauge
fields [32] and a symmetry-protected topological phase [33] were recently
reported.

Figure 1: Model. (a) Honeycomb lattice with lattice spacing $a$ and lattice
vectors $\boldsymbol{\delta}_{1,2}=(a/2,\pm\sqrt{3}a/2)$,
$\boldsymbol{\delta}_{3}=(-a,0)$. The unit cell containing the two
inequivalent sites $A$ and $B$ is also highlighted. (b) Schematic picture of
two particles interacting through a dipolar interactions $V_{\textrm{int}}$.
The particles are locally trapped at the lattice sites by a deep harmonic
potential and their polarization axis is aligned by a strong external field
that is tilted with respect to the $\hat{z}$ axis by an angle $\Theta$. The
atoms equilibrium position is indicated by dashed circles and the relative
position is indicated as $\mathbf{r}_{n,\nu}$ according to the convention used
in the text.

The previously mentioned dipolar systems are characterized by very different
interaction energies that are naturally to be compared with the level spacing
of the confining potential. The interaction energy is smaller (and even much
smaller) for magnetic atoms and polar molecules and much larger for Rydberg
atoms. In this paper, we focus on the former situation which is typically
described in terms of Hubbard-like models characterized by hopping particles
and their interactions. In the presence of long-range dipole-dipole
interactions, one can access physical regimes where the individual particle
hopping is frozen by a deep optical lattice but excitations of individual
particles can propagate via an interaction-mediated resonant exchange
mechanism. The entire dynamics of the system in this regime is, therefore,
fully controlled by the interactions. A common choice is to consider the
_internal_ atomic excitations corresponding to transitions between, for
example, Zeeman sublevels for magnetic atoms or rotational states for polar
molecules, which has led to a number of theoretical proposals and experiments,
see for example Refs. [31, 19, 21].

Following a remarkable experimental progress in creating systems of dipolar
ultracold atoms and molecules in optical lattices, we address the problem of
excitation dynamics from a new perspective, and consider the propagation of
_vibrational_ excitations of particles, which are naturally present in these
systems. Such propagating excitations are analogous to optical phonons in
crystals in the sense that they have a finite energy gap. In our case, these
vibrational excitations correspond to transitions between $s$ and $p$ orbitals
in the local harmonic oscillator basis of each well and the gap corresponds to
the energy difference between these two levels. A large energy gap relaxes
experimental requirements on temperature and technical noise, giving more
precise control over the number of excitations, as compared to gapless
acoustic-like modes. We stress that the scenario considered in our work is
completely different from the one described by $p$-band Hubbard models, see
Ref. [34] and references therein. In those systems, the $p$-state atoms are
mobile and free to move across the lattice whereas in our setup they are
localized. An additional advantage of freezing the atomic motion is the
complete elimination of onsite particle collisions, which provide the dominant
decay mechanism in $p$-band Hubbard systems [35, 36, 37].

To be more specific, we consider polarized dipolar particles (magnetic atoms
or polar molecules) in their lowest energy internal state, which are placed
with unit filling in a deep optical lattice, thus suppressing single-particle
tunneling. The considered dynamical degrees of freedom are therefore encoded
in the collective vibrational modes - phonons - originating from the long-
range dipole-dipole interactions. We consider the specific case of the
honeycomb lattice and show how one can access a broad variety of topological
effects and states by controlling the polarization direction, thus exploiting
the _anisotropy_ of the dipole-dipole interaction. These include the
instability of topological band degeneracies with respect to crystalline
deformations, Chern and anomalous Floquet topological bands, a higher-charge
monopole with a synthetic dimension, and relativistic Landau levels generated
by a strain field.

We stress that, in our setup, the implementation of the topology-related
effects and states mentioned above only requires the control of the dipole
orientation, without changing the underlying optical lattice geometry.
Furthermore, the orbital character of the phonons provides an additional
degree of freedom that enriches the band structures. The combination of these
two factors offers definite advantages over the use of ensembles of mobile
atoms in the lowest band of optical lattices, where the observation of some of
these phenomena requires ad-hoc design and distinct experimental techniques.
In contrast, our proposal provides a unifying platform, a phononic topological
toolbox, to explore a large number of topological phenomena in a _single
experimental setup_ , which have a counterpart in various classical and
condensed matter phononic systems [38, 39, 40, 41, 42, 43, 44, 45, 46, 47].

The paper is organized as follows: in Sec. 2, we introduce the model and
derive the phonon dispersion relation; in Sec. 3, we analyze the topological
properties of the phonon bands and how to control them; in Sec. 4, we discuss
several aspects related to the probing methods for phonons; in Sec. 5, we
analyze the experimental platform and the energy scales for different dipolar
atoms and polar molecules; in Sec. 6, we draw our conclusions.

##  2 Phonon modes

###  2.1 Model

Let us consider a system of polarized dipolar particles (either atoms or
molecules) on a honeycomb lattice with one atom per site, as depicted in Fig.
1, where a deep local potential allows us to neglect single-particle tunneling
processes. The collective dynamics of the ensemble is determined by the
dipole-dipole interactions and it is governed by the Hamiltonian

$H=H_{0}+V_{\textrm{int}}\,,$ (1)

where

$H_{0}=\sum_{n\in\cal A,\cal
B}\frac{\mathbf{p}_{n}^{2}}{2m}+\frac{1}{2}m\omega^{2}\sum_{n\in\cal A,\cal
B}(\mathbf{r}_{n}-\mathbf{R}_{n})^{2}\,,$ (2)

describes the single-particle Hamiltonian, with
$\mathbf{r}_{n}\\!\equiv\\!\mathbf{r}(\mathbf{R}_{n})$ and
$\mathbf{p}_{n}\\!\equiv\\!\mathbf{p}(\mathbf{R}_{n})$ the coordinate and
momentum of the particle at the site $n$ of position $\mathbf{R}_{n}$, $m$ is
the particle mass, and $\omega$ the local harmonic oscillator frequency of the
potential minima, which we take as isotropic for simplicity of notation. We
assume the low-temperature condition $k_{B}T\ll\hbar\omega$ such that we
neglect the thermal occupation of excited states in each harmonic well. The
nearest-neighbor dipole-dipole interaction takes the form

$V_{\textrm{int}}=V_{\textrm{dd}}\sum_{n\in\mathcal{B},\nu}\frac{|\mathbf{r}_{n,\nu}|^{2}-3(\hat{\mathbf{d}}\cdot\mathbf{r}_{n,\nu})^{2}}{|\mathbf{r}_{n,\nu}|^{5}}\,,$
(3)

where we defined the inter-particle distance
$\mathbf{r}_{n,\nu}=\mathbf{r}(\mathbf{R}_{n})-\mathbf{r}(\mathbf{R}_{n}+\boldsymbol{\delta}_{\nu})$,
$V_{\textrm{dd}}$ is the interaction strength, $\boldsymbol{\delta}_{\nu}$ are
the lattice vectors shown in Fig. 1(a) and the anisotropy of dipolar
interactions is described by the orientation of the dipole moment
$\hat{\mathbf{d}}=(\sin\Theta\cos\Phi,\sin\Theta\sin\Phi,\cos\Theta)$, where
the angles are chosen as in Fig. 1(b). In particular, for $\Theta\\!=\\!0$ the
dipoles are oriented orthogonal to the honeycomb lattice plane and the
interactions are therefore isotropic.

###  2.2 Dispersion relation

In order to find the collective vibrational modes, we define the small
displacements $\mathbf{u}^{A}(\mathbf{R}_{n}+\boldsymbol{\delta}_{\nu})$,
$\mathbf{u}^{B}(\mathbf{R}_{n})$, $n\in\cal B$, with respect to the
equilibrium position

$\displaystyle\mathbf{r}(\mathbf{R}_{n})$
$\displaystyle=\mathbf{R}_{n}+\mathbf{u}^{B}(\mathbf{R}_{n})\,,$
$\displaystyle\mathbf{r}(\mathbf{R}_{n}+\boldsymbol{\delta}_{\nu})$
$\displaystyle=\mathbf{R}_{n}+\boldsymbol{\delta}_{\nu}+\mathbf{u}^{A}(\mathbf{R}_{n}+\boldsymbol{\delta}_{\nu})\,,$
(4)

and we obtain the dynamical matrix

$D_{ij}^{\alpha\beta}(\mathbf{R}^{\prime},\mathbf{R}^{\prime\prime})=m\omega_{i}^{2}\delta_{ij}\delta_{\alpha\beta}+\frac{\partial^{2}V_{\textrm{int}}}{\partial
u_{i}^{\alpha}(\mathbf{R}^{\prime})\partial
u_{j}^{\beta}(\mathbf{R}^{\prime\prime})}\,,$ (5)

with $i,j=x,y,z$ the orthogonal directions of oscillation. Notice that we have
allowed for different harmonic frequencies $\omega_{i}$ along the three
spatial dimensions at each lattice site. The Hamiltonian can therefore be
recast in momentum space as

$H=\sum_{i,\alpha,\mathbf{k}}\frac{p_{i,\mathbf{k}}^{\alpha}\,p_{i,-\mathbf{k}}^{\alpha}}{2m}+\frac{1}{2}\sum_{i,j,\alpha,\beta,\mathbf{k}}u_{i,\mathbf{k}}^{\alpha}D_{ij}^{\alpha\beta}(\mathbf{k})u_{j,-\mathbf{k}}^{\beta}\,,$
(6)

where we only kept quadratic terms in $u_{i,\mathbf{k}}^{\alpha}$ and the
explicit form of $D_{ij}^{\alpha\beta}(\mathbf{k})$ is given in Appendix A.
The phonon spectrum $\hbar\omega_{\mathbf{k}}$ and the corresponding
eigenmodes are then obtained by solving the eigenvalue problem [48]

$\mathbf{D}(\mathbf{k})\cdot\mathbf{u}_{\mathbf{k}}=m\omega_{\mathbf{k}}^{2}\,\mathbf{u}_{\mathbf{k}}\,.$
(7)

Let us now consider a large isotropic trapping frequency $\hbar\omega\gg
V_{\text{dd}}\ell^{2}/a^{5}$ where $\ell=(\hbar/m\omega)^{1/2}$ is the
harmonic oscillator length at each site. In this case, we can write the
eigenfrequencies as

$\omega_{\mathbf{k}}\approx\omega+\frac{\epsilon_{\mathbf{k}}}{2m\omega}\,,$
(8)

where $\epsilon_{\mathbf{k}}$ are the eigenvalues of
$\mathbf{D}^{0}(\mathbf{k})\equiv\mathbf{D}(\mathbf{k})|_{\omega=0}$. In the
rest of this work, we will show the corresponding phonon bands of dipolar
atoms on the honeycomb lattice obtained by diagonalizing
$\mathbf{D}^{0}(\mathbf{k})$ and we will use the quantity
$E_{\text{int}}=V_{\text{dd}}\ell^{2}/2a^{5}$ as energy unit, where the factor
two in the denominator is taken as in Eq. (8).

The phonon dispersion obtained from $\mathbf{D}^{0}(\mathbf{k})$ displays six
bands in total. The in-plane modes $u_{x,\mathbf{k}}$ and $u_{y,\mathbf{k}}$
are however decoupled from the out-of-plane modes $u_{z,\mathbf{k}}$ when
$\Theta=0$, as it can be directly seen by inspecting the explicit form of the
dynamical matrix. However, in the rest of this work we will consider the in-
plane and out-of-plane modes as independent also when $\Theta\neq 0$ by
assuming $|\omega_{z}-\omega_{x,y}|\gg\Delta E_{B}/\hbar$, where $\Delta
E_{B}$ is the typical bandwidth of the phonon modes.

##  3 Controlling topological effects

###  3.1 Band crossing points for $\Theta=0$

The full dispersion relation is shown in Fig. 2(a) in the regime of isotropic
interactions, namely for $\Theta=0$ and
$\omega_{x}\\!=\\!\omega_{y}\\!=\\!\omega_{z}\\!=\\!\omega$, which displays
several band crossing points [49, 50]. The in-plane modes display two types of
band crossings. There are two quadratic band touching points (QBTPs) at
$\boldsymbol{\Gamma}$ at different energies, whereas there are two Dirac
points with a relativistic linear dispersion at $\mathbf{K}$ and
$\mathbf{K}^{\prime}=-\mathbf{K}$. The dispersion of the in-plane phonon modes
shares many similarities with the single-particle tight-binding $p$-band
spectrum on the honeycomb lattice [51], which has been observed with
polaritons [52] and electronic lattices [53], except for the fact that here
the upper and lower bands are dispersive instead of flat. The out-of-plane
modes, which resemble the $p_{z}$-orbital modes of graphene [54], also display
Dirac cones at $\mathbf{K}$ and $\mathbf{K}^{\prime}$. Notice the existence of
a triple point degeneracy at the $\boldsymbol{\Gamma}$ point near zero energy
for $\Theta=0$, which originates from the translational symmetry of
interactions. However, this band degeneracy has no interesting topological
properties as the in-plane modes are decoupled from the out-of-plane modes.
Triple band degeneracies at the $\boldsymbol{\Gamma}$ point are common in
phononic systems and they can be nevertheless characterized by unconventional
topological indices [55].

Figure 2: Phonon dispersion for $\Theta=0$. (a) Dispersion of the in-plane
modes $u_{x,\mathbf{k}}$, $u_{y,\mathbf{k}}$ (solid line) and out-of-plane
modes $u_{z,\mathbf{k}}$ (dashed line) for $\Theta=0$, namely with the
polarization orthogonal to the honeycomb lattice plane. (b) Dirac dispersion
at $\mathbf{K}=(0,4\pi/3\sqrt{3}a)$ for the in-plane modes and (c)
corresponding winding of the $\mathbf{d}_{\text{2D}}$ vector displaying a
winding number $w=1$. (d) Upper quadratic band-touching point at
$\boldsymbol{\Gamma}=(0,0)$ and (e) corresponding winding of the
$\mathbf{d}_{\text{2D}}$ vector displaying a winding number $w=2$. The
different letters $U$, $L$, $D$ indicate the band crossing points analyzed in
the main text.

The topological properties of two-band crossings in the vicinity of a momentum
$\mathbf{k}_{0}$ are universally described by an effective continuum model for
the dynamical matrix that reads

$\mathbf{D}^{0}_{\text{eff}}(\mathbf{k}_{0},\mathbf{q})=d_{0}(\mathbf{q})\mathcal{I}_{2\times
2}+\sum_{i=x,y,z}d_{i}(\mathbf{q})\sigma_{i}\,,$ (9)

where $\mathbf{q}=\mathbf{k}-\mathbf{k}_{0}$ with $|\mathbf{q}|a\ll 1$,
$\sigma_{i}$ are the Pauli matrices and $\mathcal{I}_{2\times 2}$ is the
identity matrix. The derivation of the specific form of the
$\mathbf{D}^{0}_{\text{eff}}(\mathbf{k}_{0},\mathbf{q})$ for the different
cases discussed in this work is based on Lödwin’s perturbative method, which
we outline in Appendix B.

Let us focus for the moment on the in-plane modes $u_{x,\mathbf{k}}$ and
$u_{y,\mathbf{k}}$. At the $\mathbf{K}$ point, the leading term of the
effective dynamical matrix is linear in momentum (see Fig. 2(b)) and we can
therefore neglect quadratic terms originating from higher-order perturbation
theory contributions. The resulting effective model reads

$\mathbf{D}^{0}_{\text{eff},D}(\mathbf{K},\mathbf{q})=v_{F}(q_{y}\sigma_{x}-q_{x}\sigma_{y})\,,$
(10)

with $v_{F}=45V_{\text{dd}}/4a^{4}$, which describes a massless Dirac particle
in two dimensions. We have omitted a diagonal term
$d_{0}=27V_{\text{dd}}/2a^{5}$ that shifts the energies but does not change
the topological properties carried by the eigenvectors. As a result of the
effective theory description (10), the phonon modes experience a Berry phase
$\pi$ when adiabatically encircling the $\mathbf{K}$ point, corresponding to a
winding $w=1$ of the vector $\mathbf{d}_{\text{2D}}=(d_{x},d_{y})$, see Fig.
2(c). Unless time-reversal or inversion symmetries are broken, the theory
remains gapless and the winding number or the Berry phase is a well defined
topological invariant [56]. Similarly to what happens in graphene, the winding
at $\mathbf{K}^{\prime}=-\mathbf{K}$ is opposite to the one at $\mathbf{K}$.

In order to derive the effective models describing the QBTPs at the
$\mathbf{\Gamma}$ point, it is necessary to perform second order perturbation
theory, as they contribute with quadratic terms in momentum, which are the
leading terms of a QBTP. The dispersion near the upper QBTP, see Fig. 2(d),
namely near $\epsilon_{\boldsymbol{\Gamma}}=27V_{\text{dd}}/a^{5}$, therefore
reads

$\mathbf{D}_{\textrm{eff},U}^{0}(\boldsymbol{\Gamma},\mathbf{q})=\begin{pmatrix}t_{2}(q_{x}^{2}-q_{y}^{2})&t_{1}q_{x}q_{y}\\\
t_{1}q_{x}q_{y}&-t_{2}(q_{x}^{2}-q_{y}^{2})\end{pmatrix}\,,$ (11)

with $t_{1}=-45V_{\text{dd}}/8a^{3}$ and $t_{2}=t_{1}/2$. As for the Dirac
point, we did not include an overall diagonal contribution that reads
$d_{0}=3(144+7|\mathbf{q}|^{2}a^{2})V_{\text{dd}}/16a^{5}$. A similar
expression can be obtained for the lower QBTP around
$\epsilon_{\boldsymbol{\Gamma}}=0$. The effective theory is gapless and the
winding of $\mathbf{d}_{\text{2D}}$ around the $\boldsymbol{\Gamma}$ point is
$w\\!=\\!2$, as shown in Fig. 2(e), unless time-reversal or the discrete
$C_{3}$ rotational symmetries are broken [57]. We will explore how to break
the symmetry protection of the degeneracy points in the next sections.

###  3.2 Tilting the polarization axis: $C_{3}$ crystal symmetry breaking

The anisotropy of dipolar interactions provides a knob to manipulate the
phonon band structure by breaking certain crystal symmetries that are present
when the polarizing field is orthogonal to the lattice plane. In particular,
let us consider a tilt of the polarization axis by an angle $\Theta$, while
keeping $\Phi=0$, which breaks the discrete $C_{3}$ rotational symmetry. For
small $\Theta\ll 1$, the dynamical matrix for the upper QBTP, which we
discussed in the previous section, has a correction

$\Delta\mathbf{D}_{\textrm{eff},U}^{0}=\Theta^{2}\left(\Delta
d_{0}\,\mathcal{I}_{2\times 2}+\Delta d_{z}\sigma_{z}\right)\,,$ (12)

where $\Delta d_{0}=-81V_{\text{dd}}/2a^{5}$ and $\Delta
d_{z}\\!=\\!27V_{\text{dd}}/4a^{5}$. As shown in Fig. 3(a), the QBTP then
splits into two Dirac cones as a consequence of the momentum independent term
$\Delta d_{z}\sigma_{z}$ [58]. The winding of $\mathbf{d}_{\text{2D}}$ for
$\Theta\neq 0$ is shown in Fig. 3(b), which clearly displays a winding $w=1$
around each Dirac cone such that the winding $w=2$ of the QBTP is conserved.

The lower QBTP at $\epsilon_{\boldsymbol{\Gamma}}=0$ behaves differently. A
direct calculation shows that a small angle $\Theta$ does not affect the
effective theory and the touching point is therefore protected for a small
anisotropy of the dipolar interaction. However, for larger values of $\Theta$,
the bands are more severely affected. Among the four eigenvalues at
$\boldsymbol{\Gamma}$, two of them are always zero,
$\epsilon_{\boldsymbol{\Gamma}}^{1,2}=0$. The other two eigenvalues read
$\epsilon_{\boldsymbol{\Gamma}}^{3}=27V_{\text{dd}}(3+5\cos 2\Theta)/8a^{5}$
and $\epsilon_{\boldsymbol{\Gamma}}^{4}=27V_{\text{dd}}(1+7\cos
2\Theta)/8a^{5}$. This means that there are two tilt angles that provide a
triple degeneracy at $\epsilon_{\boldsymbol{\Gamma}}=0$. In particular, we
focus on the angle $\Theta_{c}=\arccos(-3/5)/2\approx 0.352\pi$. The spectrum
is shown in Fig. 3(c) and it is captured by the effective theory

$\mathbf{D}_{\textrm{eff},T}^{0}(\boldsymbol{\Gamma},\mathbf{q})=\begin{pmatrix}0&-iv_{x}q_{x}&iv_{y}q_{y}\\\
iv_{x}q_{x}&0&-iv_{y}q_{y}\\\ -iv_{y}q_{y}&iv_{y}q_{y}&0\end{pmatrix}\,,$ (13)

where $v_{x}\\!=\\!27V_{\text{dd}}/2a^{4}$ and $v_{y}\\!=\\!v_{x}/\sqrt{2}$,
which describes a two-dimensional (pseudo)spin-1 relativistic particle.
Indeed, the model above can be rewritten as
$\mathbf{D}_{\textrm{eff},T}^{0}(\boldsymbol{\Gamma},\mathbf{q})$ $=$
$\sum_{i=x,y,z}d_{i}(\mathbf{q})S_{i}=\mathbf{d}_{i}(\mathbf{q})\cdot\mathbf{S}$,
where $\mathbf{S}=(S_{x},S_{y},S_{z})=(\lambda_{2},\lambda_{5},\lambda_{7})$
is a spin-1 representation of the SU(2) algebra satisfying
$[S_{i},S_{j}]=i\,\epsilon_{ijk}S_{k}$, $\lambda_{i}$ are Gell-Mann matrices,
whereas $d_{x}=v_{x}q_{x}$, $d_{y}=-v_{y}q_{y}$ and $d_{z}=v_{y}q_{y}$. A
discussion of the topological characterization of a general spin-1 model in 2D
has been carried out in Ref. [59], where it is applied to the Lieb lattice
case. For the time-reversal invariant regime that we are considering here, a
topological invariant associated to each band is the Berry phase around the
degeneracy point. A direct calculation, however, shows that the Berry phase of
each band is trivially zero. In Sec. 3.4, we will however show that this
triple band-crossing point can be exploited to construct a higher-dimensional
higher-charge monopole by using the angle $\Phi$ and the harmonic trapping as
additional degrees of freedom.

We conclude this section by studying the effect of the polarization tilt on
the out-of-plane phonons $u_{z,\mathbf{k}}$. These modes are described by a
graphene-like Hamiltonian that reads

$\mathbf{D}^{0}_{z}(\mathbf{k})=\begin{pmatrix}-V(0)&V(\mathbf{k})\\\
V(\mathbf{k})^{\dagger}&-V(0)\end{pmatrix}\,,$ (14)

where
$V(\mathbf{k})=-J_{1}e^{i\mathbf{k}\cdot\boldsymbol{\delta}_{1}}-J_{2}e^{i\mathbf{k}\cdot\boldsymbol{\delta}_{2}}-J_{3}e^{i\mathbf{k}\cdot\boldsymbol{\delta}_{3}}$
and

$\displaystyle J_{1,2}$
$\displaystyle=-\frac{3V_{\text{dd}}}{8a^{5}}\left(11+13\cos
2\Theta\right)\,,$ $\displaystyle J_{3}$
$\displaystyle=-\frac{3V_{\text{dd}}}{2a^{5}}\left(-1+7\cos 2\Theta\right)\,.$
(15)

The corresponding dispersion relation is shown in Fig. 3(d) for $\Theta=0$
(solid line). Applying a tilt to the polarization along one of the symmetry
axis of the lattice, namely changing $\Theta$ while keeping $\Phi=0$, creates
an anisotropy among the couplings, $J_{3}\neq J_{1,2}$. This effect shifts the
Dirac cones position (and the overall energy of the band) until the Dirac
cones merge at the $\boldsymbol{\Gamma}$ point for
$\Theta=\Theta_{c1}=-\arccos(-1/3)/2\approx 0.304\pi$. For
$\Theta>\Theta_{1c}$, the spectrum gaps out since the two Dirac cones have
opposite winding [60, 58], as shown in Fig. 3(d). Such a merging transition
has been observed with ultracold atoms [61], photonic [62, 63, 64] and
electronic systems [65].

Figure 3: Effects of interactions anisotropy. (a) Splitting of the upper QBTP
shown in Fig. 2(d) for $\Theta\\!=\\!0.03\pi$ into two Dirac cones and (b) the
corresponding winding of the $\mathbf{d}_{\text{2D}}$ vector. (c) Dispersion
of the in-plane modes for $k_{x}=0$ at the critical angle
$\Theta_{c}=\arccos(-3/5)/2$. Inset shows the triple band-crossing point near
the $\boldsymbol{\Gamma}$ point. (d) Dispersion of out-of-plane modes for
$k_{x}=0$ as a function of $\Theta$: (solid line) $\Theta=0$, (dashed line)
$\Theta=0.28\pi$ and (dotted line) $\Theta=0.32\pi$ showing that Dirac cones
merge and then gap out. (e) Disconnected dimers regime for
$\Theta=\Theta_{2c}$ for the out-of-plane modes.

Furthermore, at a larger value of the tilt angle
$\Theta_{c2}=\arccos(-11/13)/2\approx 0.410\pi$, we instead find a complete
suppression of the couplings, $J_{1}(\Theta_{c2})=J_{2}(\Theta_{c2})=0$, thus
leaving only disconnected dimers coupled by $J_{3}(\Theta_{c2})\neq 0$ (see
Fig. 3(e)). This regime opens the opportunity to realize an anomalous Floquet
topological system [66], by tilting the dipolar axis along the three crystal
symmetry axes in a clockwise or counterclockwise sequence: at each time step,
the $u_{z,\mathbf{k}}$ modes are only coupled as dimers, like in Fig. 3(d),
and the bulk dynamics is fully localized. This out-of-equilibrium phase has
vanishing Chern number in the bulk bands but nevertheless topologically
protected edge modes. The ideal protocol where at each time step only
disconnected dimers are present and the corresponding Floquet quasienergy
bands are flat, as proposed by Kitagawa et al. [66], has been realized with
optical waveguides [67, 68]. The anomalous Floquet regime has also been
accessed with ultracold atoms in optical lattices by appropriately deforming
the relative position of the optical lattice potential minima [69]. However,
the ideal flat band case remains a more challenging scenario to be realized
with ultracold atoms [70]. The degree of control of the nearest-neighbor
couplings for the phonon modes of dipolar atoms therefore paves the way to
achieve the flat-band anomalous Floquet regime, where the effect of
interactions can be enhanced and thus play a nontrivial role.

###  3.3 High-frequency rotation of the polarization axis: time-reversal
symmetry breaking

Besides breaking crystalline symmetries, the control over the tilt angle can
also be exploited to break time-reversal symmetry by applying a fast rotation
of the polarization axis with respect to the $\hat{z}$ axis. As we shall show
below, this method can be used to mimic the effect of a synthetic magnetic
field on the phonon modes, as realized for neutral atoms by circularly shaking
the lattice [10].

Let us consider a polarization axis tilted by a fixed small angle $\Theta\ll
1$ and apply a time-dependence on the angle $\Phi=\Omega t$, where $\Omega$ is
the rotation frequency, which generates a fast precession of the polarization
axis. This method has already been used to control the dipole-dipole
interaction strength and sign [71, 72], under the adiabatic assumption that
the rotation frequency is small as compared, for example, to the Zeeman
splitting in the case of dipolar magnetic atoms. We assume the condition
$\hbar\Omega\gg\Lambda$ to be valid, with $\Lambda$ playing the role of an
energy cut-off for the time-dependent theory that allows us to apply a high-
frequency Floquet expansion. Let us focus on the upper QBTP at the
$\boldsymbol{\Gamma}$ point, since the lower QBTP is not affected by small
tilting angles, as discussed above. By performing a high-frequency expansion,
we obtain the effective Floquet (time-independent) dynamical matrix, see
Appendix C for further details,

$\mathbf{D}_{\textrm{eff},U}^{F}\approx\mathbf{D}_{\textrm{eff},U}^{(0)}+\Delta\mathbf{D}_{\textrm{eff},U}^{F}\,,$
(16)

where $\mathbf{D}_{\textrm{eff},U}^{(0)}$ is the lowest order term when
$\Omega\rightarrow\infty$ and $\Delta\mathbf{D}_{\textrm{eff},U}^{F}$ is the
leading correction in $1/\Omega$. The first term,
$\mathbf{D}_{\textrm{eff},U}^{(0)}$, has the same form as Eq. (11) with
parameters: $t_{1}\\!=\\!-45V_{\text{dd}}/8a^{5}$, $t_{2}\\!=\\!t_{1}/2$ and
$d_{0}\\!=\\!3\left[72(2\\!-\\!3\Theta^{2}\\!+\\!\Theta^{4})+7|\mathbf{q}|^{2}a^{2}\right]V_{\text{dd}}/16a^{5}$.
We have kept terms up to $\Theta^{4}$ as they will appear in the next-to-
leading term and assumed that they are of the same magnitude as the terms
$q^{2}a^{2}$ that characterize the QBTP. Therefore, we have neglected terms
$q^{2}a^{2}\Theta^{2}$, which are perturbatively negligible under the previous
assumption. Time-reversal symmetry is broken if
$\Delta\mathbf{D}_{\textrm{eff},U}^{F}$ contains a term proportional to
$\sigma_{y}$ that is even in $\mathbf{q}$. Indeed, a direct calculation yields

$\Delta\mathbf{D}_{\textrm{eff},U}^{F}=\sigma_{y}\frac{729\,\Theta^{4}}{32}\frac{V_{\text{dd}}^{2}}{a^{10}}\frac{(2m\omega)^{-1}}{2\hbar\Omega}\,,$
(17)

which confirms that the rotation of the dipolar axis has the effect of
breaking time-reversal symmetry and thus opens a gap at the
$\boldsymbol{\Gamma}$ point. This result is validated numerically in Fig.
4(a), where we compare the effective Floquet theory derived above with the
numerical quasienergy spectrum obtained by the full time-dependent dynamical
matrix $\mathbf{D}^{0}(\mathbf{k},t)$. The corresponding Berry curvature
$\Omega_{xy}(\mathbf{k})$ calculated from the lower band eigenvectors is shown
in Fig. 4(b) and has the typical annular shape expected when gapping a QBTP.
After integrating over momentum, one finds a Chern number
$\nu=\int\text{d}^{2}\mathbf{k}\,\Omega_{xy}(\mathbf{k})/2\pi=-1$, where the
integration area has been extended to a full 2D plane in momentum space
because the Berry curvature is localized in the vicinity of the
$\boldsymbol{\Gamma}$ point. A similar analysis can be performed for the
$\mathbf{K}$ and $\mathbf{K}^{\prime}$ points, yielding a theory of the form
(10) for $\mathbf{D}^{F}_{\text{eff}}(\mathbf{K},\mathbf{q})$. The
perturbative correction from Floquet theory is in this case

$\Delta\mathbf{D}^{F}_{\text{eff}}(\mathbf{K},\mathbf{q})=\sigma_{z}\frac{99225\,\Theta^{4}}{64}\frac{V_{\text{dd}}^{2}}{a^{10}}\frac{(2m\omega)^{-1}}{2\hbar\Omega}\,,$
(18)

which also breaks time-reversal symmetry for the Dirac model and opens a gap.
A similar analysis can be carried out at $\mathbf{K}^{\prime}$ and the
resulting effective theory has the same form, so the Berry curvature of the
two valleys adds up to $\nu=-1$ for the upper band. Notice that the upper band
contribution at the Dirac points and the lowest band contribution at the
$\boldsymbol{\Gamma}$ point actually belong to the same band (see Fig. 2).
Therefore, the total Chern number is $\nu=-2$.

Let us now summarize the Chern numbers of the resulting bands. In total, we
have three set of bands: an upper band with $\nu=+1$, with Berry curvature
around the $\boldsymbol{\Gamma}$ point; a middle band with $\nu=-2$ with Berry
curvature around the $\boldsymbol{\Gamma}$ and the $\boldsymbol{K}$,
$\boldsymbol{K}^{\prime}$ points; a lower set of two bands with $\nu=+1$ with
Berry curvature around the $\boldsymbol{K}$, $\boldsymbol{K}^{\prime}$ points.
Recall that the lower QBTP is not affected by the dipolar axis rotation and
therefore the two lowest bands remain degenerate at $\boldsymbol{\Gamma}$.

Figure 4: Time-reversal symmetry breaking for the QBTP. (a) Quasienergy
spectrum near the upper QBTP for $\Theta=0.127\pi$ and
$\Omega=4\,E_{\text{int}}/\hbar$ computed by diagonalizing the full time-
evolution operator (solid line) and the effective Floquet theory obtained in
perturbation theory (empty circles). (b) Berry curvature of the lowest band
with the typical annular shape of a QBTP.  Figure 5: Higher-charge monopole
using a synthetic dimension. (a) Exact quasienergy spectrum near the triple
point for $\Theta=\Theta_{c}=\arccos(-3/5)/2$, $\Phi_{0}=0.5$, $\alpha=\pi/2$,
$W_{0}=0.5\,E_{\text{int}}a^{-2}$, $\Omega=3\,E_{\text{int}}/\hbar$ (solid
line) and effective perturbative Floquet theory (dashed lines). (b) Berry
curvature of the lowest band on a sphere enclosing the singularity point
$q_{1}=q_{2}=q_{3}=0$. The corresponding monopole charge and Chern number is
$\nu=2$.

###  3.4 Modulation of the harmonic frequency: Higher-charge monopole using a
synthetic dimension

In Sec. 3.2, we found that the in-plane phonon dispersion relation displays a
peculiar behaviour at a critical angle $\Theta_{c}$ where it describes
relativistic excitations of pseudo-spin $S\\!=\\!1$. As anticipated, the band
topological properties are trivial, as a direct calculation of the Berry phase
of each band shows. However, the model (13) provides the starting point to
engineer a topological monopole in 3D, as we discuss here below. In order to
obtain this result, we will exploit the concept of synthetic dimension [73,
74].

Let us consider the correction to the dynamical matrix near the triple point
obtained when changing the angle $\Phi$, with $\Phi\ll\pi$, and the trapping
frequency along the $x$ direction
$\omega_{x}\rightarrow\omega+\delta\omega_{x}$, with
$\delta\omega_{x}\ll\omega$, while keeping $\omega_{y}=\omega$ fixed. We
obtain

$\displaystyle\Delta\mathbf{D}_{\textrm{eff},T}^{0}$
$\displaystyle=\begin{pmatrix}0&0&v\,\Phi\\\ 0&0&v\,\Phi\\\
v\,\Phi&v\,\Phi&2m\,\omega\,\delta\omega_{x}\end{pmatrix}=\frac{W}{3}\mathcal{I}-\frac{W}{\sqrt{3}}\lambda_{8}+v\,\Phi(\lambda_{4}+\lambda_{6})\,,$
(19)

where in the second equality we wrote $\Delta\mathbf{D}_{\textrm{eff},T}^{0}$
in terms of Gell-Mann matrices and we defined
$v=27V_{\textrm{dd}}/5\sqrt{2}a^{5}$ and $W=2m\,\omega\,\delta\omega_{x}$. By
noticing that the following commutation relation holds,
$i[\lambda_{8},\lambda_{4}+\lambda_{6}]=-\sqrt{3}(\lambda_{5}+\lambda_{7})$,
one realizes that perturbative Floquet theory will preserve the SU(2)
algebraic structure of the time-independent model
$\mathbf{D}_{\textrm{eff},T}^{0}$, Eq. (13), which is described by
$(\lambda_{2},\lambda_{5},\lambda_{7})$, as discussed in Sec. 3.2. Another
crucial observation is that the commutator calculated above yields
$\lambda_{5}$ and $\lambda_{7}$ with same sign, whereas they appear with
opposite sign in $\mathbf{D}_{\textrm{eff},T}^{0}$ (see Eq. (13)).

Let us then apply a time modulation of the form $\Phi=\Phi_{0}\cos(\Omega
t+\alpha)$ and $W=W_{0}\cos\Omega t$. We thus find an effective Floquet theory

$\mathbf{D}_{\textrm{eff},T}^{F}\approx\mathbf{D}_{\textrm{eff},T}^{0}-\frac{(2m\omega)^{-1}W_{0}\Phi_{0}v}{2\hbar\Omega}\sin\alpha\,(\lambda_{5}+\lambda_{7})\,.$
(20)

We can now define $q_{1}=v_{x}q_{x}$, $q_{2}=v_{y}q_{y}$ and
$q_{3}=(2m\omega)^{-1}\sin\alpha\,W_{0}\Phi_{0}v/2\hbar\Omega$ and obtain

$\mathbf{D}_{\textrm{eff},T}^{F}=\begin{pmatrix}0&-iq_{1}&i(q_{2}+q_{3})\\\
iq_{1}&0&-i(q_{2}-q_{3})\\\ -i(q_{2}+q_{3})&i(q_{2}-q_{3})&0\end{pmatrix}\,,$
(21)

which corresponds to the Hamiltonian of a relativistic particle of pseudo-spin
$S\\!=\\!1$ in the three-dimensional parameter space $(q_{1},q_{2},q_{3})$
111One can actually redefine $q_{2}+q_{3}\rightarrow-q_{2}$ and
$q_{2}-q_{3}\rightarrow q_{3}$ to obtain a more symmetric structure of the
matrix that corresponds to the model introduced in [75], but here we have
decided to keep the current parametrization to avoid further manipulation of
the model since each parameter $q_{i}$ can be controlled by a single
experimental parameter., as introduced in Ref. [75]. Similar types of
dispersions [76, 77, 78, 79, 80] can describe nontrivial modes beyond the
conventional massless Dirac/Weyl representation [81]. In particular, here we
have a monopole in momentum space with charge $\nu\\!=\\!2$, which provides a
Berry curvature flux on an arbitrary surface enclosing the triple band
degeneracy at $q_{1}=q_{2}=q_{3}=0$. This can be easily verified by going to
polar coordinates, $q_{1}=r\sin\theta\cos\phi$, $q_{2}=r\sin\theta\sin\phi$,
$q_{3}=r\cos\theta$. Let us then calculate the Berry curvature on the sphere
enclosing the singularity for the lowest band, which in polar coordinates
reads [82]

$\Omega_{\theta\phi}=i\sum_{n>0}\frac{\langle
u_{0}|\partial_{\theta}\mathbf{D}_{\textrm{eff},T}^{F}|u_{n}\rangle\langle
u_{n}|\partial_{\phi}\mathbf{D}_{\textrm{eff},T}^{F}|u_{0}\rangle-(\theta\leftrightarrow\phi)}{(\epsilon_{n}-\epsilon_{0})^{2}}\,,$
(22)

where we have indicated the eigenvalues of $\mathbf{D}_{\textrm{eff},T}^{F}$
as $\epsilon_{0,2}=\mp\sqrt{q_{1}^{2}+(q_{2}+q_{3})^{2}+(q_{2}-q_{3})^{2}}$
and $\epsilon_{1}=0$, whereas $|u_{m}\rangle$ indicate the corresponding
eigenvectors. The Chern number, namely the charge that is the source of Berry
flux, can then be straightforwardly obtained as
$\nu=(1/2\pi)\int_{S^{2}}\mathrm{d}\theta\mathrm{d}\phi\,\Omega_{\theta\phi}=2$,
whereas for the upper band one finds $\nu=-2$ and a vanishing Chern number for
the middle band [75].

We have verified these results by first comparing the exact Floquet phonon
quasienergies with the effective model (21), which is shown in Fig. 5(a).
Furthermore, we have numerically computed the Berry curvature
$\Omega_{\theta\phi}$ from the Floquet eigenstates of the full time-dependent
problem and the result is shown in Fig. 5(b). Numerical integration gives
$\nu=2$, in agreement with the result from the effective theory (21). A few
experiments have realized synthetic monopoles in dimensions $d\\!>\\!2$ with
ultracold atoms, for example a Weyl point in 3D [83] or a Yang monopole in 5D
[84], by exploiting internal atomic levels as parameter space. However, to the
best of our knowledge no realization of a 3D monopole of higher charge
($\nu=2$), as the one discussed here by combining real ($q_{1}$, $q_{2}$) and
synthetic ($q_{3}$) dimensions, has been achieved so far.

###  3.5 Gradient of the polarizing external field: strain and relativistic
Landau levels

In solid state systems as graphene, the possibility to stretch a material
along certain directions by applying strain offers the opportunity to observe
interesting physical phenomena. The application of strain has the effect of
deforming the crystal structure and it therefore changes the lattice spacings,
which become space dependent. As a result, under specific strain
configurations the effective theory near the Dirac cones resembles the one of
a two-dimensional relativistic particle in a magnetic field, which displays
gapped relativistic Landau levels [85]. Similar physics has been observed, for
instance, in optical systems [86, 87] where the lattice spacing can be
engineered at will in the fabrication process. The simulation of strain with
ultracold atoms instead has not been realized so far, although theoretical
proposals are already available [88, 89]. The purpose of this section is to
show how the physics of strain can be explored with dipolar particles. In
particular, we show how the gradient of the external polarizing (electric or
magnetic) field can be used to simulate a strain field for the out-of-plane
modes $u_{z,\mathbf{k}}$ and observe the corresponding (pseudo)-Landau levels.

For simplicity, we exploit the angle dependence of the couplings in Eq. (3.2)
by assuming that the orientation of the polarizing external field is space
dependent as $\Theta=\Theta_{0}+\tau\,x/a$, where $\Theta_{0}=\pi/4$ and
$\tau$ is a dimensionless parameter quantifying the angle gradient. In this
case, the couplings take the approximate form

$\displaystyle J_{1,2}$
$\displaystyle=-\left(\frac{33}{8}-\frac{39\tau}{4a}x\right)\frac{V_{\textrm{dd}}}{a^{5}}\,,$
$\displaystyle J_{3}$
$\displaystyle=-\left(-\frac{3}{2}-\frac{21\tau}{a}x\right)\frac{V_{\textrm{dd}}}{a^{5}}\,.$
(23)

The space independent anisotropic contribution provides a shift of the Dirac
cones, as we have already shown in Fig. 3(d), which are now located at
$\mathbf{K}_{s}=(0,\pm 2\arccos(2/11)/\sqrt{3}a)$. The space dependent part
instead generates the strain field, which translates into a gauge potential
near the Dirac points. Indeed, by expanding near $\mathbf{K}_{s}$, we obtain
the low-energy theory

$\displaystyle\mathbf{D}^{0}_{\text{eff},z}(\mathbf{K}_{s},\mathbf{q})=$
$\displaystyle\,-v_{F}^{x}(q_{x}+\mathcal{A}_{x})\sigma_{y}-v_{F}^{y}(q_{y}+\mathcal{A}_{y})\sigma_{x}-(27V_{\text{dd}}/4a^{5}-\mathcal{E}_{x}x)\,\mathcal{I}_{2\times
2}\,,$ (24)

where $\mathbf{A}=(0,270\,\tau x/11a^{2})$ is an emergent vector potential in
the Landau gauge corresponding to a pseudo-magnetic field
$\mathcal{B}=\partial_{x}\mathcal{A}_{y}-\partial_{y}\mathcal{A}_{x}=270\tau/11a^{2}$,
$\mathcal{E}_{x}=81\tau V_{\text{dd}}/2a^{6}$ is an emergent (pseudo)electric
field, $v_{F}^{x}=9V_{\text{dd}}/4a^{4}$ and
$v_{F}^{y}=9\sqrt{39}V_{\text{dd}}/8a^{4}$ are anisotropic velocities. Note
that in deriving the low-energy theory we have neglected the spatial
dependences of the Fermi velocities [90]. We have therefore found that near
the Dirac cone at $\mathbf{K}_{s}$ the model describes a particle in the
presence of a (pseudo)magnetic and a (pseudo)electric field. Near the other
cone at $-\mathbf{K}_{s}$, the opposite (pseudo)magnetic field will appear, as
time reversal is not broken [91].

The (pseudo)electric field $\mathcal{E}_{x}$ is quite strong. Indeed, by
direct inspection one obtains the value $\mathcal{E}_{x}/v_{F}\mathcal{B}>1$,
where $v_{F}=\sqrt{v^{x}_{F}v^{y}_{F}}$ is the Fermi velocity geometric mean.
As a result, the drift velocity $v_{D}=\mathcal{E}/\mathcal{B}$ is larger than
the mean Fermi velocity $v_{F}$, thus destroying the Landau level structure of
the spectrum originating from the (pseudo)magnetic field. An overview of the
different regimes is provided in Ref. [91]. The most convenient strategy to
observe Landau levels physics in our setup is therefore to compensate for
$\mathcal{E}_{x}$, _e. g._ by applying an external linear force of opposite
sign to the phonons. We discuss in the next section how such a force can be
engineered. When $\mathcal{E}_{x}=0$, the spectrum of the low-energy theory
(24) takes the form of relativistic Landau levels that reads

$\epsilon^{\text{LL}}_{n}=\pm\sqrt{v^{x}_{F}v^{y}_{F}\mathcal{B}|n|}\,,$ (25)

where we have omitted an overall energy shift $-27V_{\text{dd}}/4a^{5}$. In
Fig. 6(a), we show the exact spectrum for a system described by the hopping
amplitudes (3.5) with $\tau=0.001$ on a lattice that is infinite in the $y$
direction and has $N_{x}=50$ unit cells in the $x$ direction. In order to
compare with the effective Landau level picture discussed above, we subtracted
the space-dependent contribution $V(0)=-(J_{1}+J_{2}+J_{3})$ (see Eq. (14)),
which corresponds to taking $\mathcal{E}_{x}=0$ in the low-energy theory. In
Fig. 6(b), we compare the dispersion with the analytical result (25). From the
inspection of the spectrum, we see that the energy levels display a momentum
dependence, which originates from a spatial dependence of the Fermi velocities
that we dropped in the derivation of the effective low-energy theory [90].

Figure 6: Simulating relativistic Landau levels from strain. (a) Spectrum of
out-of-plane modes with $\Theta_{0}=\pi/4$ and $\tau=0.001$ on a lattice that
is infinite in the $y$ direction and has $N_{x}=50$ unit cells in the $x$
direction. We have subtracted the diagonal term $V(0)$. (b) Zoom of the
spectrum shown in panel (a) near one of the Dirac points showing (solid line)
the exact result, where the couplings are parametrized as in Eq. (3.2), and
(dashed line) the Landau levels $\epsilon^{\textrm{LL}}_{n}$ from Eq. (25) for
$n=1,2,3$.

The previous results are obtained by considering a small gradient of the
magnetic field around the tilt angle $\Theta_{0}=\pi/4$, which was chosen to
minimizes the corrections beyond linear order in Eqs. (3.5). However, large
tilt angles reduce the overall bandwidth and bring the Dirac points close to
each other near the gapping transition, as shown in Fig. 3(d). To circumvent
these possible limitations, it could therefore be convenient to reduce the
tilt angle in the experimental setting in order to increase the overall
bandwidth and the Landau levels energy spacing. In that case, however, the
corrections to Eqs. (3.5) will be of quadratic order in $x$ instead of cubic
order, which occurs for the case $\Theta_{0}=\pi/4$ discussed here. The strong
anisotropy, may also affect the wavepacket dynamics, as we shall see in the
next section. A distinct possibility to reach a more ideal regime where to
observe Landau level physics could be to start from a highly anisotropic
honeycomb lattice from the beginning with
$|\boldsymbol{\delta}_{3}|>|\boldsymbol{\delta}_{1,2}|$ and restore the
$C_{3}$ crystal symmetry through the polarization tilt. By following this
strategy, one could therefore optimize the bandwidth, suppress unwanted strong
anisotropic effects and have a sufficiently homogeneous pseudo(magnetic)
field.

Figure 7: Hall deflection of a phonon wavepacket. (a) Semiclassical wavepacket
dynamics for the out-of-plane modes with $\mathbf{F}=(0,F_{y})$,
$F_{y}=0.63\,E_{\text{int}}/a$ and width $\sigma=30a$. An onsite energy
difference $\Delta=4.5\,E_{\text{int}}/a$ has been applied to open the gap at
$\mathbf{K}$, thus generating a non-vanishing Berry curvature that is
responsible for the Hall drift along $x$. The small oscillations are caused by
non-adiabatic processes originating from transfer to the upper band. (b)
Initially prepared gaussian wavepacket, where the dashed vertical line
indicates the initial $\langle x\rangle$ position. (c) Final state of the
wavepacket showing a horizontal displacement with respect to the initial
$\langle x\rangle$ position. In panels $(b)$ and $(c)$, the normalized
wavepacket density $n$ is shown.

##  4 Probing methods

In this Section, we discuss several aspects related to the control of a phonon
wavepacket (see the end of the Section for its preparation) in order to
explore the topological properties of the phonon bands, as the Berry curvature
or the Berry phase. A strategy that is well established with ultracold atoms
requires a tilted potential, which plays the role of an electric field for
charged particles, to investigate the transport properties of wavepackets in
individual or multiple bands. As a result, the wavepacket satisfies
semiclassical equations, which for a non-degenerate band read

$\dot{\mathbf{r}}_{c}=\nabla_{\mathbf{k}}\epsilon(\mathbf{k}_{c})-\mathbf{F}\times\mathbf{\Omega}_{xy}(\mathbf{k}_{c})\,,\quad\dot{\mathbf{k}}_{c}=\mathbf{F}\,,$
(26)

where $\mathbf{r}_{c}$ and $\mathbf{k}_{c}$ are the center-of-mass position
and momentum of the wavepacket, $\mathbf{F}$ is the applied force,
$\epsilon(\mathbf{k}_{c})$ the band dispersion and
$\mathbf{\Omega}_{xy}(\mathbf{k}_{c})$ the Berry curvature [10]. The force
$\mathbf{F}$, which moves an atomic wavepacket in real space, is typically
obtained by accelerating the lattice or by applying a magnetic field gradient.
However, these techniques do not produce a force on a phonon wavepacket and a
different approach, which we detail below, must be followed.

In order to correctly reproduce the effects of an applied force on the phonon
modes, we find it useful to recast the formalism that we developed in Sec. 2.2
into a tight-binding model in real space, as described in Ref. [39] for the
phonon modes of trapped ion arrays. As outlined in Appendix D, the Hamiltonian
reads

$\displaystyle H$
$\displaystyle=\sum_{n}\hbar\tilde{\omega}_{i,n}\,a^{\dagger}_{i,n}a_{i,n}\\!+\\!\sum_{\begin{subarray}{c}i,j\\\
\langle n,n^{\prime}\rangle\end{subarray}}\frac{\hbar
D_{ij}^{nn^{\prime}}}{2m\sqrt{\tilde{\omega}_{i,n}\tilde{\omega}_{j,n^{\prime}}}}(a_{i,n}^{\dagger}a_{j,n^{\prime}}\\!+\\!\text{h.c.})$
$\displaystyle\equiv\sum_{i,n}\mu_{i,n}\,a^{\dagger}_{i,n}a_{i,n}+\sum_{\begin{subarray}{c}i,j\\\
\langle
n,n^{\prime}\rangle\end{subarray}}J_{ij}^{nn^{\prime}}(a_{i,n}^{\dagger}a_{j,n^{\prime}}+\text{h.c.})\,,$
(27)

where $\tilde{\omega}^{2}_{i,n}=\omega^{2}_{i,n}+D_{ii}^{nn}/m$ and we
neglected phonon non-conserving terms like $a_{i,n}a_{j,n^{\prime}}$ or
$a^{{\dagger}}_{i,n}a^{\dagger}_{j,n^{\prime}}$ under the rotating wave
approximation, which is valid when $\hbar\omega_{i,n}\gg
V_{\text{dd}}\ell^{2}/a^{5}$. Notice that we have recovered the result in Eq.
(8) for the spectrum of the model and at the same time we have recast the
Hamiltonian into an effective tight-binding description with $\mu_{i,n}$ the
onsite energies and $J_{ij}^{nn^{\prime}}$ the hopping amplitudes. The first
term in Eq. (4), namely $\mu_{i,n}$, clearly shows that a space-dependent
harmonic frequency, _e. g._ $\omega_{i,n}=\omega_{i,0}+\Delta\omega\,X_{n}/a$
corresponds to an onsite potential $\mu_{i,n}$ for the phonon modes with a
nonvanishing space gradient in the $x$ direction, which can therefore be used
to generate the force $\mathbf{F}$ acting on a phonon wavepacket. This
approach would naturally also lead to space-dependent couplings, as evident
from the second term of Eq. (4). However, adiabaticity requires that the
magnitude of the force must be sufficiently small to prevent Landau-Zener
transitions when the phonon wavepacket approaches a band gap. This would
reflect into a slow variation of the couplings $J_{ij}^{nn^{\prime}}$, which
needs to be minimized in order to avoid detrimental effects on the topological
gap.

In Fig. 7, we show an example of a gaussian wavepacket dynamics on the
honeycomb lattice for the $u_{z,\mathbf{k}}$ modes. We assume a uniform and
isotropic case, with $J_{1}\\!=\\!J_{2}\\!=\\!J_{3}\\!=\\!J$ and we break
inversion symmetry by changing the relative onsite energy between the $A$ and
$B$ sites by imposing an offset $\Delta$, which can be obtained by making the
harmonic frequencies on the two sublattices different, _i.e._
$\omega^{A}\neq\omega^{B}$. We prepare a gaussian wavepacket at
$\boldsymbol{\Gamma}$ and apply a linear potential corresponding to a force
$\mathbf{F}=(0,F_{y})$. The offset $\Delta$ opens a gap at the $\mathbf{K}$
and $\mathbf{K}^{\prime}$ points with a nonvanishing (but opposite at the two
valleys) Berry curvature. As a result, the wavepacket experiences a Hall
deflection after crossing the $\mathbf{K}$ point under the applied
$\mathbf{F}$ force, as expected from Eq. (26). This strategy can be applied to
control a phonon wavepacket in order to probe the Berry curvature generated by
the driving schemes discussed before or to measure the Berry phase around the
band degeneracy points using interferometric methods [92].

Figure 8: Cyclotron orbit motion in the presence of strain. (a) A gaussian
wavepacket of out-of-plane modes, width $\sigma=20a$, is prepared at the
$\boldsymbol{\Gamma}$ and rapidly accelerated to the $\textbf{K}_{s}$ point,
where it is left to evolve in the absence of any force. The polarization tilt
parameters are chosen as $\Theta_{0}=\pi/4$ and $\tau=0.001$. The large
difference in the $x$ and $y$ motion reflects the strong anisotropy of the
Dirac velocity. (b) Center of mass motion up to $t=12\,\hbar
E_{\text{int}}^{-1}$ after preparing the wavepacket at the two valleys
$\pm\textbf{K}_{s}$. The trajectory shows opposite chiral motion, as expected
from the underlying time-reversal symmetry of the model.

For the anisotropic and nonuniform case with space dependent couplings
$J_{\nu}(x)$, as in Eq. (3.5), we saw in Sec. 3.5 that a strong pseudo-
electric field contribution $\mathcal{E}_{x}$ must be compensated in order to
access the Landau levels regime. In this case, after noticing that the onsite
energies read $\tilde{\omega}=\omega\\!-\\!J_{T}(x)/2m\omega$, where
$J_{T}(x)\equiv(J_{1}(x)\\!+\\!J_{2}(x)\\!+\\!J_{3}(x))$, and therefore
$\mathcal{E}_{x}\propto\partial_{x}J_{T}(x)$, the strong electric field can be
suppressed by an appropriate space-dependent trapping frequency
$\omega(x)=\omega_{0}+\Delta\omega\,x/a$ by taking
$\Delta\omega/a=\partial_{x}J_{T}/2m\omega_{0}$ such that $\tilde{\omega}$ is
approximately a constant. If $\tilde{\omega}$ has no space space dependence or
a residual small one, the Landau level structure is preserved and the related
physics can therefore be observed. We prepared a wavepacket with momentum
$\mathbf{K}_{s}$ that we let evolve in the presence of strain. We show in Fig.
8(a) the expected dynamics, where the wavepacket performs a chiral motion
corresponding to cyclotron orbits, originating from the nonvanishing pseudo-
magnetic field $\mathcal{B}$ induced by strain. In Fig. 8(b), we verify that
the motion displays opposite chirality at $-\mathbf{K}_{s}$, as expected from
time-reversal symmetry, namely the fact that the pseudo-magnetic field
$\mathcal{B}$ changes sign from one valley to the other. The strong anisotropy
reflects into substantially different Fermi velocities in the two spatial
directions, which cause a larger displacement of the wavepacket along the $y$
direction as compared to the $x$ direction.

In order to resolve the phonons occupation and monitor the position of the
phonon wavepacket, as required to observe the Hall deflection and more
generally transport properties, we may take advantage of single-site
microscopy. A phonon excitation locally corresponds to a superposition of the
atomic (or molecular) wavefunction including the lowest $s$-orbital of each
site and the higher orbitals, dominantly the $p$-orbitals. A measurement of
the phonon wavepacket density could then proceed in two steps: in the first
step, a $\pi$ microwave pulse resonantly couples the $s$ orbitals to a
different atomic state that is not trapped by the laser confinement. As a
result, the ground state occupation can be removed, thus projecting the
wavefunction onto the excited orbitals, which map the phonon modes occupation.
In a second step, single site microscopy is then employed to image the
position over each lattice site, thus resolving the phonon wavepacket real-
space distribution.

While the discussion on the probing methods has mainly focused on the
transport properties of a phonon wavepacket to extract the topological
properties of the bands, as Berry phases or Berry curvatures, a distinct,
complementary and obvious possibility is to perform spectroscopic measurements
of the phonon bands. In this case, one can access the phonon dispersion and
directly map out the band singularities, their splitting or the gap opening
events. Raman [36] or Bragg [93] beams can therefore be exploited to perform
band spectroscopy of the phonon dispersion, while being at the same time
employed for the state preparation of a phonon wavepacket.

##  5 Experimental realizations

Let us now discuss potential experimental realizations of our phonon system
and present the corresponding estimates of the relevant energy scales. The
most natural candidates for dipolar particles are magnetic atoms (and also
diatomic molecules made of magnetic atoms) and polar molecules placed in a
constant magnetic or electric field, respectively.

More precisely, we consider magnetic atoms or diatomic molecules (Er, Dy atoms
or $\mathrm{Er_{2}}$, $\mathrm{Dy_{2}}$ diatomic molecules, for example) in
their lowest energy state in a magnetic field of the order of $0.1\div 1$
Gauss which corresponds to energy gaps of the order of $1\div
50\,h\,\mathrm{MHz}$ between Zeemann levels. In the case of polar molecules
(KRb, RbCs, NaK, NaRb, for example), we consider them in their electronic and
vibrational ground state in an external constant electric field $\mathbf{E}$
which induces an electric dipole moment. For our estimates, we take the
achievable value of $12\,\mathrm{kV}/\mathrm{cm}$. In Table 1, we present the
values of masses, magnetic/electric dipole moments, as well as the
corresponding values of the dipole-dipole interaction parameter
$V_{\mathrm{dd}}=\mu_{0}\mu^{2}/(4\pi)$ or
$V_{\mathrm{dd}}=d^{2}/(4\pi\epsilon_{0})$ for two magnetic ($\mu$) or
electric ($d$) dipole moments, respectively, which correspond to their
interaction at the distance $1\mu m$, for several magnetic/electric dipolar
particles.

| ${m}/{\mathrm{u}}$ | $\mu/\mu_{0}$ | $V_{\mathrm{dd}}/(h\,\mathrm{kHz}\leavevmode\nobreak\ \mu\mathrm{m}^{3})$
---|---|---|---
Er | $164$ | $7$ | $0.00064$
Dy | $161$ | $10$ | $0.0013$
Er2 | $328$ | $\lesssim 14$ | $0.0025$
Dy2 | $322$ | $\lesssim 20$ | $0.0052$
|  | ${\mathcal{D}}/{\mathrm{D}}$ | $d/{\mathrm{D}}$ |
KRb | $127$ | $0.57$ | $0.34$ | $0.017$
RbCs | $220$ | $1.17$ | 0.95 | $0.136$
NaK | $63$ | $2.72$ | $1.90$ | $0.547$
NaRb | $110$ | $3.2$ | 2.50 | $0.916$
LiRb | $91$ | $3.99$ | 2.35 | $0.836$
Table 1: Parameters for a variety of magnetic atoms/molecules and polar
molecules. For all polar molecular species the offset electric field is chosen
to be $|\mathbf{E}|=12\,\mathrm{kV}/\mathrm{cm}$. Parameters for different
magnetic atoms can be found in Ref. [94] and for polar molecules in Ref. [95].

To get an estimate of the phonon bandwidth, we consider a 1D situation where
the dipolar particles are placed in a deep optical lattice with a lattice
spacing $a=\lambda/2$ and a depth $V_{0}E_{R}$, where $E_{R}$ is the recoil
energy $E_{R}=\pi^{2}\hbar^{2}/(2ma^{2})$ and $V_{0}\gg 1$. In this case, the
hopping amplitude $t$ between different sites of the lattice is exponentially
small. For example, $t/E_{R}=(4/\sqrt{\pi})V_{0}^{3/4}\exp(-2\sqrt{V_{0}})$
[96] for the nearest-neighbor tunneling of particles between the lowest
Wannier states of each site, such that the individual particles tunneling
dynamics is frozen. When considering one particle per site, the spatial
degrees of freedom reduce to local oscillations/vibrations which, in the
harmonic approximations, can be characterized by the oscillator frequency
$\omega\\!=\\!\sqrt{2V_{0}}E_{R}/\hbar$ and by the oscillator length
$\ell\\!=\\!\sqrt{\hbar/(m\omega)}\\!=\\!a\pi^{-1}(2/V_{0})^{1/4}$. To be more
precise, the harmonic approximation, resulting in an equidistant spectrum, is
only valid for a few lowest energy states in a very deep optical lattices. For
experimentally realistic situations, the level spacing decreases with
increasing the energy, and $\hbar\omega$ gives the difference between the
ground and the first excited states ($s$ and $p$ orbitals in the 3D language).
By using the expression $12V_{\mathrm{dd}}/a^{5}$ for the off-diagonal element
of the dynamical matrix [see Eq. (A) adapted to the 1D case discussed here],
the phonon bandwidth $\Delta E_{B}$ is given by $\Delta E_{B}$ $=$
$24(\ell/a)^{2}V_{\mathrm{dd}}/a^{3}=24(\sqrt{2/V_{0}}/\pi^{2})V_{\mathrm{dd}}/a^{3}$,
such that $\Delta
E_{B}/(\hbar\omega)=(24/\pi^{2})(V_{0}E_{R})^{-1}V_{\mathrm{dd}}/a^{3}\sim
a^{-1}$.

In Fig. 9, we present the values of $\Delta E_{B}$ for a fixed lattice spacing
$a=0.266\mathrm{\mu m}$ and different lattice depths $V_{0}.$ We see that the
values of $\Delta E_{B}$ range from $10\div 300\,h\,\mathrm{Hz}$ for magnetic
systems to $1\div 50\,h\,\mathrm{kHz}$ for molecular ones. We mention also
that $\Delta E_{B}$ can reach values of more than $100\,h\,\mathrm{kHz}$ for
polar molecules with larger dipole moment (NaCs and LiCs have a dipole moment
$d=3.8\,\mathrm{D}$ and $d=3.67\,\mathrm{D}$, respectively, with an electric
field as in Table 1, for example) but this situation is not covered by our
approximations as the phonon bandwidth becomes larger than the harmonic
frequency.

Finally, we mention another possible experimental implementation of optical
phonons which is based on Rydberg atoms in tweezer arrays. This case, however,
requires more careful analysis because of much smaller values of the parameter
$\ell/a$, as compared to optical lattices. As a result, phonon bands of
similar widths would require much stronger values of the interparticle
interaction and, therefore, the creation of stable structures can be
problematic. Similar arguments also apply to the case of molecules with large
dipole moments.

Figure 9: Phonon bandwidth estimate for dipolar atoms and polar molecules. (a)
Phonon bandwidth in units of the harmonic trapping frequency as a function of
the lattice depth $V_{0}$ for a lattice spacing $a=0.266\,\mu\mathrm{m}$. (b)
Phonon bandwidth in energy units as a function of the lattice depth for the
same lattice spacing used in panel (a).

##  6 Discussion and perspective

In this work, we have shown that arrays of ultracold dipolar particles offer a
versatile platform to explore the physics of topological states of matter by
studying the vibrational phonon modes in a deep optical lattice. The current
experimental progress provides all the necessary ingredients to pursue this
approach. We have shown how the sole static and driven control of the dipole-
dipole interactions anisotropy can be instrumental to access a variety of
topological effects. While the results of this work have been derived for the
honeycomb lattice case, our methods are quite general. They can be applied to
other two-dimensional lattices, _e. g._ the Lieb or the Kagome, as well as to
one- and three-dimensional lattices, where new topological phenomena can be
accessed. Our analysis has been simplified under the condition that the
dynamics of in-plane modes is independent of the out-of-plane modes one, and
viceversa. However, relaxing this condition in three dimensions can offer an
additional nontrivial set of possibilities that is left to a future
investigation [55].

An interesting perspective to be investigated is the fate of the edge modes
that are expected in the gaps of the topological phonon bands. This discussion
requires a comprehensive analysis of the boundary details. The region of the
system where the density drops to zero provides a natural boundary for the
existence of edge modes. If the external confinement is a box potential, this
corresponds to a hard boundary and the edge modes will therefore be present.
Less clear is the existence of edge modes in a different scenario, namely in
the presence of a confining harmonic potential. This case corresponds to a
soft boundary because the particle density smoothly decreases, thus breaking
the condition of unit filling with unclear consequences for the dynamics of
phonon excitations in this region. Note also that our setup offers a distinct
and rather simple opportunity to create and study topological modes localized
at defects in the bulk of the system. Such defects with arbitrary shape can
indeed be created by removing dipolar particles with single-site addressing
techniques.

In our analysis, we assumed a low phonon density, which allowed us to neglect
phonon-phonon interactions. These will provide nontrivial scattering processes
for phonons as they will couple the different directions of oscillations,
namely their orbital degrees of freedom. Collisions therefore play an
important role to establish intra- and inter-band scattering events that can
lead, for example, to the “thermalization” of the phonon system. In this
sense, a comprehensive analysis of the regime with a large phonon density
offers the opportunity to explore remarkable orbital-like many-body phenomena
[37].

We mention also that the strategies considered in this paper, which exploit
the anisotropy of the dipole-dipole interaction to manipulate the band
structure of the phonon excitations, can also be applied to other types of
excitations whose transport is governed by a resonant dipole-exchange
mechanism. Besides the case of spin-like excitations with magnetic atoms and
polar molecules, an interesting possibility is offered by the high-energy
orbital excitations in arrays of Rydberg atoms [16, 24, 33, 32].

##  7 Acknowledgements

We acknowledge discussions with F. Ferlaino, N. Goldman, M. Mark, M. Norcia,
G. Palumbo, G. Salerno, J. Yu. We also acknowledge support by the European
Union program Horizon 2020 under Grant No. 817482 (PASQuanS), by the QuantERA
grant MAQS via the Austrian Science Fund FWF No I4391-N, by the Simons
Collaboration on Ultra-Quantum Matter, which is a grant from the Simons
Foundation (651440, PZ), and by a joint-project grant from the FWF (Grant No.
I 4426, RSF/Russia 2019).

##  Appendix A Dynamical matrix

The phonon spectrum is obtained by diagonalizing the dynamical matrix, which
for the system considered in this work reads

$\displaystyle D_{ij}^{\alpha\alpha}(\mathbf{k})$
$\displaystyle=m\omega_{i}^{2}\delta_{ij}+\sum_{\nu}f_{ij}^{\nu}\,,$
$\displaystyle D_{ij}^{AB}(\mathbf{k})$
$\displaystyle=-\sum_{\nu}f_{ij}^{\nu}\,e^{i\mathbf{k}\cdot\boldsymbol{\delta}_{\nu}}\,,$
$\displaystyle D_{ij}^{BA}(\mathbf{k})$
$\displaystyle=D_{ij}^{AB}(-\mathbf{k})\,,$ (28)

with

$\displaystyle\left(\frac{V_{\text{dd}}}{a^{5}}\right)^{-1}f_{ij}^{\nu}=$
$\displaystyle\,3\,\delta_{ij}\left[-1+5a^{-2}(\hat{\mathbf{d}}\cdot\boldsymbol{\delta}_{\nu})^{2}\right]-6\,\hat{d}_{i}\hat{d}_{j}$
$\displaystyle+30\,a^{-2}(\hat{\mathbf{d}}\cdot\boldsymbol{\delta}_{\nu})\left[\hat{d}_{i}\,\delta_{\nu}^{j}+\hat{d}_{j}\,\delta_{\nu}^{i}\right]$
$\displaystyle+15\,\delta_{\nu}^{i}\delta_{\nu}^{j}\left[1-7a^{-2}(\hat{\mathbf{d}}\cdot\boldsymbol{\delta}_{\nu})^{2}\right]\,.$
(29)

##  Appendix B Effective continuum theory

Near a band touching point occurring at $\mathbf{k}_{0}$ and for small momenta
$|\mathbf{q}|a\equiv|\mathbf{k}-\mathbf{k}_{0}|a\ll 1$, the effective two-band
theory can generally be cast in the form

$\mathbf{D}^{0}_{\text{eff}}(\mathbf{k}_{0},\mathbf{q})=d_{0}(\mathbf{q})\mathcal{I}_{2\times
2}+\sum_{i=x,y,z}d_{i}(\mathbf{q})\sigma_{i}\,,$ (30)

where $\sigma_{i}$ are the Pauli matrices and $\mathcal{I}_{2\times 2}$ is the
identity matrix. In order to obtain the form of the effective theory, we use
Lödwin’s method [97]. We therefore split the modes into a low-energy sector,
which lies in energy window near the degeneracy point, and into a high-energy
sector, which is treated in perturbation theory. We therefore construct the
low-energy projector
$P=\sum_{\rho,\sigma}^{\prime}|\epsilon_{\mathbf{k}_{0}}^{\rho}\rangle\langle\epsilon_{\mathbf{k}_{0}}^{\sigma}|$
and the high-energy projector
$Q=\sum_{\rho,\sigma}^{\prime\prime}|\epsilon_{\mathbf{k}_{0}}^{\rho}\rangle\langle\epsilon_{\mathbf{k}_{0}}^{\sigma}|$,
where $|\epsilon_{\mathbf{k}_{0}}^{\rho}\rangle$ represent the eigenvectors at
$\mathbf{k}_{0}$ and the symbols $\sum^{\prime}$ and $\sum^{\prime\prime}$
indicate that the summation is restricted to the low-energy and high-energy
manifolds, respectively. We obtain the effective low-energy theory

$\mathbf{D}^{0}_{\textrm{eff}}=P\mathbf{D}^{0}(\mathbf{k}_{0}\\!+\\!\mathbf{q})P+P\mathbf{\tilde{D}}Q(\epsilon_{\mathbf{k}_{0}}^{\prime}-Q\mathbf{\tilde{D}}Q)^{-1}(Q\mathbf{\tilde{D}}P)\,,$
(31)

where
$\mathbf{\tilde{D}}(\mathbf{q})\\!=\\!\mathbf{D}^{0}(\mathbf{k}_{0}\\!+\\!\mathbf{q})-\mathbf{D}^{0}(\mathbf{k}_{0})$
and $\epsilon_{\mathbf{k}_{0}}^{\prime}$ is the characteristic energy of the
low-energy sector.

##  Appendix C Floquet perturbation theory

Let us consider a theory for the phonon modes subject to a periodic drive
described by a dynamical matrix $\mathbf{D}^{0}(t)$. We can therefore
decompose the matrix as follows

$\mathbf{D}^{0}(t)=\sum_{l}e^{il\Omega t}\mathbf{D}^{(l)}\,,$ (32)

where $\mathbf{D}^{(l)}$ denotes the $l$-th Fourier component and $\Omega$ the
driving frequency. This problem can be described through Floquet theory [98]
by constructing a time-independent Floquet dynamical matrix $\mathbf{D}^{F}$
that governs the stroboscopic time evolution, which reads

$\frac{\mathbf{D}^{F}}{2m\omega}=\frac{i\hbar}{T}\log\mathbf{U}(T)\,,$ (33)

where
$\mathbf{U}(t)=\mathcal{T}\exp(-i\hbar^{-1}\int_{0}^{t}\text{d}t\,\mathbf{D}(t)/2m\omega)$
is the time-evolution operator and $T=2\pi/\Omega$.

By performing a high-frequency expansion, the Floquet effective (time-
independent) dynamical matrix can be computed as

$\mathbf{D}^{F}\approx\mathbf{D}^{(0)}+\frac{(2m\omega)^{-1}}{\hbar\Omega}\sum_{l>0}\frac{1}{l}[\mathbf{D}^{(+l)},\mathbf{D}^{(-l)}]\,,$
(34)

where we have indicated only the leading order in $1/\Omega$.

##  Appendix D Tight-binding description

The phonon modes description that we have used throughout the manuscript is
based on the construction of a dynamical matrix describing the small
vibrations of the ensemble. However, under the circumstances described here
below, the description can be recast in a tight-binding model picture. The
Hamiltonian in momentum space reads

$H=\sum_{i,\alpha,\mathbf{k}}\frac{p_{i,\mathbf{k}}^{\alpha}\,p_{i,-\mathbf{k}}^{\alpha}}{2m}+\frac{1}{2}\sum_{i,j,\alpha,\beta,\mathbf{k}}u_{i,\mathbf{k}}^{\alpha}D_{ij}^{\alpha\beta}(\mathbf{k})u_{j,-\mathbf{k}}^{\beta}\,,$
(35)

which can be cast in real space as

$\displaystyle H$
$\displaystyle=\sum_{i,n}\frac{p_{i,n}^{2}}{2m}+\frac{1}{2}\sum_{i,j,n,n^{\prime}}u_{i,n}D_{ij}^{nn^{\prime}}u_{j,n^{\prime}}\,$
(36)
$\displaystyle=\sum_{i,n}\frac{p_{i,n}^{2}}{2m}+\frac{1}{2}m\sum_{i,n}\tilde{\omega}^{2}_{i,n}u_{i,n}^{2}+\frac{1}{2}\sum_{\begin{subarray}{c}i,j\\\
\langle
n,n^{\prime}\rangle\end{subarray}}u_{i,n}D_{ij}^{nn^{\prime}}u_{j,n^{\prime}}\,,$

where $n,n^{\prime}$ indicate the lattice sites (irrespective of the
sublattice degree of freedom), $i,j=x,y,z$ and
$\tilde{\omega}^{2}_{i,n}=\omega_{i,n}^{2}+D_{ii}^{nn}/m$. We now introduce
the local harmonic oscillator creation and annihilation operators as follows

$\displaystyle u_{i,n}$
$\displaystyle=\sqrt{\frac{\hbar}{2\,m\,\tilde{\omega}_{i,n}}}(a_{i,n}+a_{i,n}^{\dagger})\,,$
(37) $\displaystyle p_{i,n}$
$\displaystyle=i\sqrt{\frac{m\,\hbar\,\tilde{\omega}_{i,n}}{2}}(a_{i,n}^{\dagger}-a_{i,n})\,,$
(38)

and rewrite the Hamiltonian as

$\displaystyle H$
$\displaystyle=\sum_{n}\hbar\tilde{\omega}_{i,n}\,a^{\dagger}_{i,n}a_{i,n}\\!+\\!\sum_{\begin{subarray}{c}i,j\\\
\langle n,n^{\prime}\rangle\end{subarray}}\frac{\hbar
D_{ij}^{nn^{\prime}}}{2m\sqrt{\tilde{\omega}_{i,n}\tilde{\omega}_{j,n^{\prime}}}}(a_{i,n}^{\dagger}a_{j,n^{\prime}}\\!+\\!\text{h.c.})$
$\displaystyle\equiv\sum_{i,n}\mu_{i,n}\,a^{\dagger}_{i,n}a_{i,n}+\sum_{\begin{subarray}{c}i,j\\\
\langle
n,n^{\prime}\rangle\end{subarray}}J_{ij}^{nn^{\prime}}(a_{i,n}^{\dagger}a_{j,n^{\prime}}+\text{h.c.})\,,$
(39)

where we neglected phonon non-conserving terms like $a_{i,n}a_{j,n^{\prime}}$
or $a^{{\dagger}}_{i,n}a^{\dagger}_{j,n^{\prime}}$ under the rotating wave
approximation (RWA), which is valid when $\hbar\omega_{i,n}\gg
V_{\text{dd}}\ell^{2}/a^{5}$. As a result, the RWA allows us to reinterpret
the phonon dynamics in terms of particles hopping on a lattice. This picture
become extremely convenient in order, for instance, to understand how to
manipulate a phonon wavepacket, as discussed in the main text.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Chen  _et al._ [2012] X. Chen, Z.-C. Gu, Z.-X. Liu, and X.-G. Wen, Science 338, 1604 (2012).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Rachel [2018] S. Rachel, Reports on Progress in Physics 81, 116501 (2018).
  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * Laughlin [1981] R. B. Laughlin, Physical Review B 23, 5632 (1981).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Fu [2011] L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Goldman  _et al._ [2016] N. Goldman, J. C. Budich, and P. Zoller, Nature Physics 12, 639 (2016).
  * Cooper  _et al._ [2019] N. R. Cooper, J. Dalibard, and I. B. Spielman, Rev. Mod. Phys. 91, 015005 (2019).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Rev. Mod. Phys. 91, 015006 (2019).
  * Lee  _et al._ [2018] C. H. Lee, S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, and R. Thomale, Communications Physics 1, 39 (2018).
  * Huber [2016] S. D. Huber, Nature Physics 12, 621 (2016).
  * Ma  _et al._ [2019] G. Ma, M. Xiao, and C. T. Chan, Nature Reviews Physics 1, 281 (2019).
  * Lahaye  _et al._ [2009] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Reports on Progress in Physics 72, 126401 (2009).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nature Physics 6, 382 (2010).
  * Baranov  _et al._ [2012] M. A. Baranov, M. Dalmonte, G. Pupillo, and P. Zoller, Chemical Reviews 112, 5012 (2012).
  * Gross and Bloch [2017] C. Gross and I. Bloch, Science 357, 995 (2017).
  * de Paz  _et al._ [2013] A. de Paz, A. Sharma, A. Chotia, E. Maréchal, J. H. Huckans, P. Pedri, L. Santos, O. Gorceix, L. Vernac, and B. Laburthe-Tolra, Phys. Rev. Lett. 111, 185305 (2013).
  * Baier  _et al._ [2016] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Z. Cai, M. Baranov, P. Zoller, and F. Ferlaino, Science 352, 201 (2016).
  * Yan  _et al._ [2013] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Nature 501, 521 (2013).
  * Labuhn  _et al._ [2016] H. Labuhn, D. Barredo, S. Ravets, S. de Léséleuc, T. Macrì, T. Lahaye, and A. Browaeys, Nature 534, 667 (2016).
  * Peter  _et al._ [2015] D. Peter, N. Y. Yao, N. Lang, S. D. Huber, M. D. Lukin, and H. P. Büchler, Phys. Rev. A 91, 053617 (2015).
  * Weber  _et al._ [2018] S. Weber, S. de Léséleuc, V. Lienhard, D. Barredo, T. Lahaye, A. Browaeys, and H. P. Büchler, Quantum Science and Technology 3, 044001 (2018).
  * Schuster  _et al._ [2021] T. Schuster, F. Flicker, M. Li, S. Kotochigova, J. E. Moore, J. Ye, and N. Y. Yao, Phys. Rev. Lett. 127, 015301 (2021).
  * Syzranov  _et al._ [2014] S. V. Syzranov, M. L. Wall, V. Gurarie, and A. M. Rey, Nature Communications 5, 5391 (2014).
  * Syzranov  _et al._ [2016] S. V. Syzranov, M. L. Wall, B. Zhu, V. Gurarie, and A. M. Rey, Nature Communications 7, 13543 (2016).
  * Salerno  _et al._ [2020] G. Salerno, G. Palumbo, N. Goldman, and M. Di Liberto, Phys. Rev. Research 2, 013348 (2020).
  * Manmana  _et al._ [2013] S. R. Manmana, E. M. Stoudenmire, K. R. A. Hazzard, A. M. Rey, and A. V. Gorshkov, Phys. Rev. B 87, 081106 (2013).
  * Yao  _et al._ [2013] N. Y. Yao, A. V. Gorshkov, C. R. Laumann, A. M. Läuchli, J. Ye, and M. D. Lukin, Phys. Rev. Lett. 110, 185302 (2013).
  * Yao  _et al._ [2018] N. Y. Yao, M. P. Zaletel, D. M. Stamper-Kurn, and A. Vishwanath, Nature Physics 14, 405 (2018).
  * Lienhard  _et al._ [2020] V. Lienhard, P. Scholl, S. Weber, D. Barredo, S. de Léséleuc, R. Bai, N. Lang, M. Fleischhauer, H. P. Büchler, T. Lahaye, and A. Browaeys, Phys. Rev. X 10, 021031 (2020).
  * de Léséleuc  _et al._ [2019] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Science 365, 775 (2019).
  * Li and Liu [2016] X. Li and W. V. Liu, Reports on Progress in Physics 79, 116401 (2016).
  * Isacsson and Girvin [2005] A. Isacsson and S. M. Girvin, Phys. Rev. A 72, 053604 (2005).
  * Müller  _et al._ [2007] T. Müller, S. Fölling, A. Widera, and I. Bloch, Phys. Rev. Lett. 99, 200405 (2007).
  * Wirth  _et al._ [2011] G. Wirth, M. Ölschläger, and A. Hemmerich, Nature Physics 7, 147 (2011).
  * Prodan and Prodan [2009] E. Prodan and C. Prodan, Phys. Rev. Lett. 103, 248101 (2009).
  * Bermudez  _et al._ [2011] A. Bermudez, T. Schaetz, and D. Porras, Phys. Rev. Lett. 107, 150501 (2011).
  * Kane and Lubensky [2014] C. L. Kane and T. C. Lubensky, Nature Physics 10, 39 (2014).
  * Süsstrunk and Huber [2015] R. Süsstrunk and S. D. Huber, Science 349, 47 (2015).
  * Stenull  _et al._ [2016] O. Stenull, C. L. Kane, and T. C. Lubensky, Phys. Rev. Lett. 117, 068001 (2016).
  * Salerno  _et al._ [2016] G. Salerno, T. Ozawa, H. M. Price, and I. Carusotto, Phys. Rev. B 93, 085105 (2016).
  * Li  _et al._ [2020] J. Li, L. Wang, J. Liu, R. Li, Z. Zhang, and X.-Q. Chen, Phys. Rev. B 101, 081403 (2020).
  * Wei  _et al._ [2021] Q. Wei, X. Zhang, W. Deng, J. Lu, X. Huang, M. Yan, G. Chen, Z. Liu, and S. Jia, Nature Materials 20, 812 (2021).
  * Luo  _et al._ [2021] L. Luo, H.-X. Wang, Z.-K. Lin, B. Jiang, Y. Wu, F. Li, and J.-H. Jiang, Nature Materials 20, 794 (2021).
  * Lange  _et al._ [2022] G. F. Lange, A. Bouhon, B. Monserrat, and R.-J. Slager, Phys. Rev. B 105, 064301 (2022).
  * Bruus and Flensberg [2004] H. Bruus and K. Flensberg,  _Many-Body Quantum Theory in Condensed Matter Physics: An Introduction_ (Oxford University Press, 2004).
  * Schnyder [2018] A. P. Schnyder, Topological Matter School (2018).
  * Feng  _et al._ [2021] X. Feng, J. Zhu, W. Wu, and S. A. Yang, Chinese Physics B 30, 107304 (2021).
  * Wu and Das Sarma [2008] C. Wu and S. Das Sarma, Phys. Rev. B 77, 235107 (2008).
  * Jacqmin  _et al._ [2014] T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Phys. Rev. Lett. 112, 116402 (2014).
  * Gardenier  _et al._ [2020] T. S. Gardenier, J. J. van den Broeke, J. R. Moes, I. Swart, C. Delerue, M. R. Slot, C. M. Smith, and D. Vanmaekelbergh,  _ACS Nano_ , ACS Nano 14, 13638 (2020).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Park  _et al._ [2021] S. Park, Y. Hwang, H. C. Choi, and B.-J. Yang, Nature Communications 12, 6781 (2021).
  * Bernevig and Hughes [2013] B. A. Bernevig and T. L. Hughes,  _Topological Insulators and Topological Superconductors_ (Princeton University Press, 2013).
  * Sun  _et al._ [2009] K. Sun, H. Yao, E. Fradkin, and S. A. Kivelson, Phys. Rev. Lett. 103, 046811 (2009).
  * Montambaux  _et al._ [2018] G. Montambaux, L.-K. Lim, J.-N. Fuchs, and F. Piéchon, Phys. Rev. Lett. 121, 256402 (2018).
  * He  _et al._ [2012] Y. He, J. Moore, and C. M. Varma, Phys. Rev. B 85, 155106 (2012).
  * Montambaux  _et al._ [2009] G. Montambaux, F. Piéchon, J.-N. Fuchs, and M. O. Goerbig, Phys. Rev. B 80, 153412 (2009).
  * Tarruell  _et al._ [2012] L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature 483, 302 (2012).
  * Bellec  _et al._ [2013] M. Bellec, U. Kuhl, G. Montambaux, and F. Mortessagne, Phys. Rev. Lett. 110, 033902 (2013).
  * Rechtsman  _et al._ [2013a] M. C. Rechtsman, Y. Plotnik, J. M. Zeuner, D. Song, Z. Chen, A. Szameit, and M. Segev, Phys. Rev. Lett. 111, 103901 (2013a).
  * Milićević  _et al._ [2019] M. Milićević, G. Montambaux, T. Ozawa, O. Jamadi, B. Real, I. Sagnes, A. Lemaître, L. Le Gratiet, A. Harouri, J. Bloch, and A. Amo, Phys. Rev. X 9, 031010 (2019).
  * Kim  _et al._ [2015] J. Kim, S. S. Baik, S. H. Ryu, Y. Sohn, S. Park, B.-G. Park, J. Denlinger, Y. Yi, H. J. Choi, and K. S. Kim, Science 349, 723 (2015).
  * Kitagawa  _et al._ [2010] T. Kitagawa, E. Berg, M. Rudner, and E. Demler, Phys. Rev. B 82, 235114 (2010).
  * Mukherjee  _et al._ [2017] S. Mukherjee, A. Spracklen, M. Valiente, E. Andersson, P. Öhberg, N. Goldman, and R. R. Thomson, Nature Communications 8, 13918 (2017).
  * Maczewsky  _et al._ [2017] L. J. Maczewsky, J. M. Zeuner, S. Nolte, and A. Szameit, Nature Communications 8, 13756 (2017).
  * Wintersperger  _et al._ [2020] K. Wintersperger, C. Braun, F. N. Ünal, A. Eckardt, M. Di Liberto, N. Goldman, I. Bloch, and M. Aidelsburger, Nature Physics 16, 1058 (2020).
  * Quelle  _et al._ [2017] A. Quelle, C. Weitenberg, K. Sengstock, and C. Morais Smith, New Journal of Physics 19, 113010 (2017).
  * Giovanazzi  _et al._ [2002] S. Giovanazzi, A. Görlitz, and T. Pfau, Phys. Rev. Lett. 89, 130401 (2002).
  * Tang  _et al._ [2018] Y. Tang, W. Kao, K.-Y. Li, and B. L. Lev, Phys. Rev. Lett. 120, 230401 (2018).
  * Celi  _et al._ [2014] A. Celi, P. Massignan, J. Ruseckas, N. Goldman, I. B. Spielman, G. Juzeliūnas, and M. Lewenstein, Phys. Rev. Lett. 112, 043001 (2014).
  * Ozawa and Price [2019] T. Ozawa and H. M. Price, Nature Reviews Physics 1, 349 (2019).
  * Bradlyn  _et al._ [2016] B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353, aaf5037 (2016).
  * Green  _et al._ [2010] D. Green, L. Santos, and C. Chamon, Phys. Rev. B 82, 075104 (2010).
  * Fulga and Stern [2017] I. C. Fulga and A. Stern, Phys. Rev. B 95, 241116 (2017).
  * Zhu  _et al._ [2017] Y.-Q. Zhu, D.-W. Zhang, H. Yan, D.-Y. Xing, and S.-L. Zhu, Phys. Rev. A 96, 033634 (2017).
  * Tan  _et al._ [2018] X. Tan, D.-W. Zhang, Q. Liu, G. Xue, H.-F. Yu, Y.-Q. Zhu, H. Yan, S.-L. Zhu, and Y. Yu, Phys. Rev. Lett. 120, 130503 (2018).
  * Hu and Zhang [2018] H. Hu and C. Zhang, Phys. Rev. A 98, 013627 (2018).
  * Armitage  _et al._ [2018] N. Armitage, E. Mele, and A. Vishwanath, Reviews of Modern Physics 90, 015001 (2018).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Wang  _et al._ [2021] Z.-Y. Wang, X.-C. Cheng, B.-Z. Wang, J.-Y. Zhang, Y.-H. Lu, C.-R. Yi, S. Niu, Y. Deng, X.-J. Liu, S. Chen, and J.-W. Pan, Science 372, 271 (2021).
  * Sugawa  _et al._ [2018] S. Sugawa, F. Salces-Carcoba, A. R. Perry, Y. Yue, and I. B. Spielman, Science 360, 1429 (2018).
  * Guinea  _et al._ [2010] F. Guinea, M. I. Katsnelson, and A. K. Geim, Nature Physics 6, 30 (2010).
  * Rechtsman  _et al._ [2013b] M. C. Rechtsman, J. M. Zeuner, A. Tünnermann, S. Nolte, M. Segev, and A. Szameit, Nature Photonics 7, 153 (2013b).
  * Jamadi  _et al._ [2020] O. Jamadi, E. Rozas, G. Salerno, M. Milićević, T. Ozawa, I. Sagnes, A. Lemaître, L. Le Gratiet, A. Harouri, I. Carusotto, J. Bloch, and A. Amo, Light: Science & Applications 9, 144 (2020).
  * Tian  _et al._ [2015] B. Tian, M. Endres, and D. Pekker, Phys. Rev. Lett. 115, 236803 (2015).
  * Jamotte  _et al._ [2022] M. Jamotte, N. Goldman, and M. Di Liberto, Communications Physics 5, 30 (2022).
  * Salerno  _et al._ [2015] G. Salerno, T. Ozawa, H. M. Price, and I. Carusotto, 2D Materials 2, 034015 (2015).
  * Goerbig [2011] M. O. Goerbig, Rev. Mod. Phys. 83, 1193 (2011).
  * Duca  _et al._ [2015] L. Duca, T. Li, M. Reitter, I. Bloch, M. Schleier-Smith, and U. Schneider, Science 347, 288 (2015).
  * Ernst  _et al._ [2010] P. T. Ernst, S. Götze, J. S. Krauser, K. Pyka, D.-S. Lühmann, D. Pfannkuche, and K. Sengstock, Nature Physics 6, 56 (2010).
  * Frisch [2014] A. Frisch,  _Dipolar Quantum Gases of Erbium_ , Ph.D. thesis (2014).
  * Aymar and Dulieu [2005] M. Aymar and O. Dulieu,  _The Journal of Chemical Physics_ , The Journal of Chemical Physics 122, 204302 (2005).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Lim  _et al._ [2020] L.-K. Lim, J.-N. Fuchs, F. Piéchon, and G. Montambaux, Phys. Rev. B 101, 045131 (2020).
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Phys. Rev. X 4, 031027 (2014).