Topological magnon amplification

TL;DR

This paper proposes a method to induce and detect topological magnon edge modes using electromagnetic driving, enabling applications like magnon amplification and lasers in spintronics.

Contribution

It introduces a novel electromagnetic driving technique to realize and observe topological magnon edge modes and demonstrates their potential for magnonic devices.

Findings

Electromagnetic driving causes edge mode instabilities.

Large non-equilibrium magnon edge currents are generated.

Potential for topological magnon amplifiers and lasers.

Abstract

Topology is quickly becoming a cornerstone in our understanding of electronic systems. Like their electronic counterparts, bosonic systems can exhibit a topological band structure, but in real materials it is difficult to ascertain their topological nature, as their ground state is a simple condensate or the vacuum, and one has to rely instead on excited states, for example a characteristic thermal Hall response. Here we propose driving a topological magnon insulator with an electromagnetic field and show that this causes edge mode instabilities and a large non-equilibrium steady-state magnon edge current. Building on this, we discuss several experimental signatures that unambiguously establish the presence of topological magnon edge modes. Furthermore, our amplification mechanism can be employed to power a topological travelling-wave magnon amplifier and topological magnon laser, with applications in magnon spintronics. This work thus represents a step toward functional topological magnetic materials.

Figures

Equations

$${\hat{H}}_{\mathrm{int}}=\sum _k\frac{g_k}{2}\left({\hat{a}}_{-k}^{†}{\hat{a}}_k^{†}\hat{b}+{\hat{b}}^{†}{\hat{a}}_k{\hat{a}}_{-k}\right).$$

$$i\frac{d}{dt}{\mathbf{A}}_k=\left(\begin{array}{cc}{\tilde{\omega }}_k-i\frac{\gamma }{2}& \mathcal{E}\\ -\mathcal{E}& -{\tilde{\omega }}_{-k}-i\frac{\gamma }{2}\end{array}\right){\mathbf{A}}_k,$$

$${\omega }_{k,\pm }=\frac{{\omega }_k-{\omega }_{-k}}{2}-\frac{i\gamma }{2}\pm \sqrt{\frac{({\tilde{\omega }}_k+{\tilde{\omega }}_{-k}{)}^2}{4}-\mathcal{E}}.$$

$${\hat{H}}_0=\sum _{⟨ij⟩}\left[-J{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j+{\mathbf{D}}_{ij}\cdot ({\hat{\mathbf{S}}}_i\times {\hat{\mathbf{S}}}_j)\right]-g_{\mathrm{L}}{\mu }_{\mathrm{B}}{\mathbf{H}}_0\cdot \sum _i{\hat{\mathbf{S}}}_i.$$

$${\hat{H}}_0=-\frac{1}{2}(J+i{D}_z)\sum _{⟨mn⟩}{\hat{a}}_m^{†}{\hat{a}}_n+\text{H.c.}+h\sum _m{\hat{a}}_m^{†}{\hat{a}}_m+K_0,$$

$${\hat{\mathbf{P}}}_{mn}$$

$$\equiv {\mathbf{Q}}_{mn}{\hat{a}}_m{\hat{a}}_n+\text{H.c.}+\cdots .$$

$$\hat{H}=\sum _{k,s}{\tilde{\omega }}_{k,s}{\hat{b}}_{k,s}^{†}{\hat{b}}_{k,s}-\frac{1}{2}\left[{\mathbf{E}}_0\cdot {\tilde{\mathbf{Q}}}_{s{s}^{\prime }}(k){\hat{b}}_{k,s}{\hat{b}}_{-k,s^{\prime }}+\text{H.c.}\right],$$

$$i\frac{d}{dt}{\mathbf{A}}_k=\left(\begin{array}{cc}{\tilde{\omega }}_k-i\frac{\gamma +\eta |{\alpha }_k{|}^2}{2}& \mathcal{E}\\ -\mathcal{E}& -{\tilde{\omega }}_{-k}-i\frac{\gamma +\eta |{\alpha }_k{|}^2}{2}\end{array}\right){\mathbf{A}}_k.$$

$$J_{\mathrm{SS}}={\int }_{-\Lambda }^{\Lambda }\frac{dq}{2\pi }|{\alpha }_{\pi +q}{|}^2{\nu }_{\pi +q},$$

$$J_{\mathrm{SS}}(\mathcal{E})\approx \frac{5{\omega }_{\pi }^{\prime }}{6\eta \pi }{\left(\frac{4{\mathcal{E}}^2-{\gamma }^2}{({\omega }_{\pi }^{\prime \prime }{)}^2}\right)}^{\frac{1}{4}}(2\mathcal{E}-\gamma ).$$

$$J_{\mathrm{net}}^x={\kappa }_{xy}({\mathcal{E}}_{\mathrm{avg}}){\Delta }_y\mathcal{E},{\kappa }_{xy}(\mathcal{E})\equiv \frac{d{J}_{\mathrm{SS}}(\mathcal{E})}{d\mathcal{E}}$$

$${\Delta }_{x}T=\frac{{\Delta }_y\mathcal{E}}{{\kappa }_{xx}}{\kappa }_{xy}(\mathcal{E}).$$

$$H_0=K_0-\frac{1}{2}\{(J+i{D}_z)\sum _{k,l_y}[e^{-ik/4}{a}_{1,k,l_y}^{†}{a}_{2,k,l_y}$$

$$+e^{-ik/4}{a}_{2,k,l_y}^{†}{a}_{3,k,l_y}+2\cos (k/2)a_{3,k,l_y}^{†}{a}_{1,k,l_y}$$

$$+e^{ik/4}{a}_{1,k,l_y}^{†}{a}_{2,k,l_y-1}+e^{ik/4}{a}_{2,k,l_y}^{†}{a}_{3,k,l_y+1}]+\text{H.c.}\}$$

$$H=H_0-H_{\text{rot}}+H_{\text{amp}}.$$

$$H=\sum _k\left[a_{k,s}^{†}{\mu }_{k,s{s}^{\prime }}{a}_{k,s^{\prime }}+\frac{1}{2}\left(a_{k,s}^{†}{\nu }_{k,s{s}^{\prime }}{a}_{-k,s^{\prime }}^{†}+\text{H.c.}\right)\right]$$

$$\frac{d}{dt}|a_k⟩=-i{\sigma }_z{h}_k|a_k⟩,$$

$$h_k=\left(\begin{array}{cc}{\mu }_k& {\nu }_k\\ {\nu }_k^{†}& {\mu }_{-k}^T\end{array}\right).$$

$$i\left(\begin{array}{c}{\dot{\alpha }}_k\\ {\dot{\alpha }}_{-k}^{\ast }\end{array}\right)=\left(\begin{array}{cc}{\tilde{\omega }}_k-i\frac{\gamma +\eta |{\alpha }_k{|}^2}{2}& \mathcal{E}\\ -\mathcal{E}& -{\tilde{\omega }}_{-k}-i\frac{\gamma +\eta |{\alpha }_k{|}^2}{2}\end{array}\right)\left(\begin{array}{c}{\alpha }_k\\ {\alpha }_{-k}^{\ast }\end{array}\right).$$

$$\frac{1}{2}{q}^2{\omega }_{\pi }^{\prime \prime }-\frac{i}{2}(\gamma +\eta |{\alpha }_k{|}^2)+z\mathcal{E}$$

$$-\frac{1}{2}{q}^2{\omega }_{\pi }^{\prime \prime }-\frac{i}{2}(\gamma +\eta |z{\alpha }_k{|}^2)-\frac{\mathcal{E}}{z}$$

$$|{\alpha }_k{|}^2=\frac{1}{\eta }\left(-2i\mathcal{E}z-i{q}^2{\omega }_{\pi }^{\prime \prime }-\gamma \right).$$

$$|{\alpha }_{\pi +q}{|}^2=\frac{1}{\eta }\left(\sqrt{4{\mathcal{E}}^2-q^4({\omega }_{\pi }^{\prime \prime }{)}^2}-\gamma \right).$$

$$J_{\mathrm{SS}}$$

$$=\frac{2{\omega }_{\pi }^{\prime }\sqrt{2\mathcal{E}}}{3\pi \eta \sqrt{{\omega }_{\pi }^{\prime \prime }}}$$

$$-\gamma \sqrt[4]{1-{\gamma }^2/(4{\mathcal{E}}^2)}\}.$$

$${\dot{n}}_n-i[H_0,n_n]={\dot{n}}_n-i\sum _m[h_m,n_n]=0,$$

$$H_0=\sum _n{h}_n.$$

$$H_{\mathrm{Hubbard}}=-\sum _{⟨ij⟩}\left[{\mathbf{c}}_i^{†}\left({\tau }_{ij}+{\mathbf{d}}_{ij}\cdot \mathbf{\sigma }\right){\mathbf{c}}_j+\text{H.c.}\right]+U\sum _j{n}_{j\uparrow }{n}_{j\downarrow },$$

$${\tilde{H}}_{\mathrm{eff}}=\sum _{⟨ij⟩}\frac{4t_{ij}^2}{U}{\mathbf{S}}_i^{†}{\mathcal{J}}_{ij}{\mathbf{S}}_j,$$

$${\mathcal{J}}_{ij}{\mathbf{S}}_j=\cos (2{\theta }_{ij}){\mathbf{S}}_j+\sin (2{\theta }_{ij})({\mathbf{S}}_j\times {\mathbf{n}}_{ij})+2{\sin }^2({\theta }_{ij}){\mathbf{n}}_{ij}({\mathbf{n}}_{ij}\cdot {\mathbf{S}}_j).$$

Full text

I Introduction

While fermionic topological insulators have a number of clear experimental signatures accessible through linear transport measurements Hasan and Kane (2010); Qi and Zhang (2011), noninteracting bosonic systems with topological band structure have a simple condensate or the vacuum as their ground state Vishwanath and Senthil (2013), making it more difficult to ascertain their topological nature. Their excited states, however, may carry signatures of the topology of the band structure, for example in form of a thermal Hall response Katsura et al. (2010); Onose et al. (2010); Hirschberger et al. (2015). There is great interest in certifying and exploiting topological edge modes in bosonic systems, as they are chiral and robust against disorder, making them a great resource to realize backscattering-free waveguides Haldane and Raghu (2008); Wang et al. (2009) and potentially topologically protected travelling-wave amplifiers Peano et al. (2016). It has been predicted that topological magnon insulators (TMI) are realized, e.g., in kagome planes of certain pyrochlore magnetic insulators as a result of Dzyaloshinskii-Moriya (DM) interaction Katsura et al. (2010); Zhang et al. (2013); Mook et al. (2014). To date, there exists only indirect experimental proof, via neutron scattering measurements of the bulk band structure in Cu[1,3-benzenedicarboxylate (bdc)] Chisnell et al. (2015), and observation of a thermal magnon Hall effect in Lu2V2O7 Onose et al. (2010) and Cu(1,3-bdc) Hirschberger et al. (2015). The main obstacle is that magnons are uncharged excitations and thus invisible to experimental tools like STM or ARPES with spatial resolution. An unambiguous experimental signature, such as the direct observation of an edge mode in the bulk gap is hampered by limitations in energy resolution (in resonant x-ray scattering) or signal strength (in neutron scattering) Chisnell et al. (2015).

Here, we propose driving a magnon edge mode to a parametric instability, which, when taking into account nonlinear damping, induces a non-equilibrium steady state with a large chiral edge mode population. In such a state the local polarization and magnetization associated with the edge mode are coherently enhanced, which could enable direct detection of edge modes via neutron scattering. Crucially, we show that selective amplification of edge modes can be achieved while preserving the stability of the bulk modes and thus the magnetic order. Another key experimental signature we predict is that applying a driving field gradient gives rise to a temperature gradient along the transverse direction, thus establishing what one might call a driven Hall effect (DHE). Topological magnon amplification has further uses in magnon spintronics Chumak et al. (2015), providing a way to amplify magnons and to build a topological magnon laser Harari et al. (2018); Bandres et al. (2018). Our work on driving topological edge modes in magnetic materials complements previous investigations in ultracold gases Galilo et al. (2015, 2017), photonic crystals Peano et al. (2016), and most recently arrays of semiconductor microresonators Harari et al. (2018); Bandres et al. (2018) and graphene Plank et al. (2018).

II Results

II.1 Edge mode parametric instability

Before discussing a microscopic model, we show qualitatively how a parametric instability may arise from anomalous magnon pairing terms in a chiral one-dimensional waveguide. We consider bosonic modes $\{\hat{a}_{k}\}$ with energies $\{\omega_{k}\}$ (for example the magnon edge mode between the first and second band, as in Fig. 1), labelled by momentum $k$, interacting with another bosonic mode (electromagnetic field mode) $\hat{b}$, as described by generic three-wave mixing ($\hbar=1$)

[TABLE]

Under strong coherent driving, the bosonic annihilation operator $\hat{b}$ can be replaced by its classical amplitude $\hat{b}\approx\beta\exp(-i\Omega_{0}t)\gg 1$, yielding an effective Hamiltonian $\hat{H}=\hat{H}_{0}+g_{k}\beta[\hat{a}_{-k}^{\dagger}\hat{a}_{k}^{\dagger}\exp(i\Omega_{0}t)+\hat{a}_{k}\hat{a}_{-k}\exp(-i\Omega_{0}t)]$. The second term produces magnon pairs with equal and opposite momentum. The time-dependence can be removed by passing to a rotating frame with respect to $\sum_{k}(\Omega_{0}/2)\hat{a}_{k}^{\dagger}\hat{a}_{k}$. From the Hamiltonian it is straightforward to derive the equations of motion, which couple particles at momentum $k$ with holes at $-k$. Neglecting fluctuations, we focus on the classical amplitudes of the fields $\alpha_{k}=\langle\hat{a}_{k}\rangle$ and include a phenomenological linear damping rate $\gamma_{k}$ to take into account the various damping processes present in such materials Chisnell et al. (2015); Chernyshev and Maksimov (2016). As we are interested in amplification around a small bandwidth, we neglect the momentum dependence of the coupling $g_{k}\simeq g$ and damping $\gamma_{k}\simeq\gamma$, arriving at

[TABLE]

where we have introduced the frequency relative to the rotating frame $\tilde{\omega}_{k}=\omega_{k}-\Omega_{0}/2$, the vector $\mathbf{A}_{k}=(\alpha_{k},\alpha_{-k}^{*})$, and the overall coupling strength $\mathcal{E}=g\beta$. The eigenvalues of the dynamical matrix Eq. 2 are the complex energies

[TABLE]

If the coupling $\mathcal{E}$ exceeds the energy difference between pump photons and magnon pair (the detuning) $\tilde{\omega}_{k}+\tilde{\omega}_{-k}=\omega_{k}+\omega_{-k}-\Omega_{0}$, the square root becomes imaginary. If further its magnitude exceeds $\gamma$, more particles are created than dissipated, causing an instability and exponential growth of the number of particles in this mode. Eventually the growth is limited by nonlinear effects, as discussed below. Despite its simplicity, Eq. 3 provides an accurate account of the fundamental instability mechanism in two-dimensional kagome TMIs, as is illustrated by the quantitative agreement Fig. 2. This forms the key ingredient for directly observing chiral magnon edge modes.

II.2 Microscopic model

Turning to a more realistic model, we consider spins on the vertices of an insulating kagome lattice ferromagnet that interact via Heisenberg and Dzyaloshinskii-Moriya (DM) interaction

[TABLE]

Here, $\mathbf{D}_{ij}$ is the DM vector that can in principle differ from bond to bond, but is heavily constrained by lattice symmetries. $\mathbf{H}_{0}$ is an externally applied magnetic field, $\mu_{\mathrm{B}}$ the Bohr magneton, $g_{\mathrm{L}}$ the Landé g-factor, and $J$ is the Heisenberg interaction strength. This model has been found to describe the thermal magnon Hall effect in Lu2V2O7 Onose et al. (2010), as well as the bulk magnon band structure of Cu(1,3-bdc) Chisnell et al. (2015).

The low energy excitations around the ferromagnetic order are magnons, whose bilinear Hamiltonian is obtained from a standard Holstein-Primakoff transformation to order $1/S$ along the direction of magnetization, i.e., $\hat{S}^{+}=\sqrt{2s}\,\hat{a},\hat{S}^{-}=\sqrt{2s}\,\hat{a}^{\dagger},\hat{S}^{z}=s-\hat{a}^{\dagger}\hat{a}$ Katsura et al. (2010); Zhang et al. (2013), yielding

[TABLE]

where $K_{0}$ is a constant, the sum ranges over bonds directed counterclockwise in each triangle, and we have chosen the magnetic field to point along $z$, introducing $h\equiv g_{\mathrm{L}}\mu_{\mathrm{B}}H_{0}^{z}$.

To second order, the Hamiltonian only contains the component of $\mathbf{D}_{ij}$ along $z$ ($D_{z}$), which is the same for all bonds due to symmetry. We take the unit cell to be one upright triangle (red in Fig. 1), with sites $\bm{\rho}_{1}=(0,0)$, $\bm{\rho}_{2}=(-1,\sqrt{3})/4$, $\bm{\rho}_{3}=(-1/2,0)$. The unit cells form a triangular Bravais lattice generated by the lattice vectors $\bm{\delta}_{1}=(1,0)$, $\bm{\delta}_{2}=(1,\sqrt{3})/2$, $\bm{\delta}_{3}=\bm{\delta}_{2}-\bm{\delta}_{1}=(-1,\sqrt{3})/2$. For nonzero $D_{z}$, the bands in this model are topological Zhang et al. (2013); Chisnell et al. (2015) causing exponentially localized edge modes to appear within the band gaps.

The effect of an oscillating electric field on magnons in a TMI is characterized by the polarization operator, which can be expanded as a sum of single-spin terms, products of two spins, three spins, etc. Moriya (1968) Lattice symmetries restrict which terms may appear in the polarization tensor Moriya (1968). In the pyrochlore lattice, the polarization due to single spins (linear Stark effect) vanishes, as each lattice site is a centre of inversion, such that the leading term contains two spin operators. The associated tensor can be decomposed into the isotropic (trace) part $\bm{\pi}$, as well as the anisotropic traceless symmetric and antisymmetric parts $\bm{\Gamma}$ and $\mathbf{D}$, viz. $\hat{\mathbf{P}}_{jl}=(\bm{\pi}_{jl}\delta^{\beta\gamma}+\mathbf{\Gamma}_{jl}^{(\beta\gamma)}+\mathbf{D}_{jl}^{[\beta\gamma]})\hat{S}_{j}^{\beta}\hat{S}_{l}^{\gamma}$ (sum over $\beta,\gamma$ implied). Kagome TMIs generically have a nonzero anisotropic symmetric part, which implies the presence of anomalous magnon pairing terms in the spin-wave picture

[TABLE]

The polarization enters the Hamiltonian via coupling to the amplitude of the electric field, $\hat{H}(t)=\hat{H}_{0}-\mathbf{E}(t)\cdot\hat{\mathbf{P}}$, thus introducing terms that create a pair of magnons while absorbing a photon. Pair production of magnons is a generic feature of antiferromagnets (via $\bm{\pi}$), Moriya (1968) but since in ferromagnets it relies on anisotropy, it is expected to be considerably weaker. A microscopic calculation based on a third-order hopping process in the Fermi-Hubbard model at half filling reveals that $|\mathbf{Q}|=ae(t/U)^{3}/2$, where $a$ is the lattice vector, $e$ the elementary charge, $t$ the hopping amplitude, and $U$ the on-site repulsion (see Supplementary Note 1).

As in the chiral waveguide model, assume an oscillating electric field $\mathbf{E}(t)=\mathbf{E}_{0}\cos(\Omega_{0}t)$. We consider an infinite strip with $W$ unit cells along $y$, but remove the lowest row of sites to obtain a manifestly inversion-symmetric model. Diagonalizing the undriven Hamiltonian $\hat{H}_{0}$ (4), we label the eigenstates $b_{k,s}$ by their momentum along $x$ and an index $s\in\{1,2,\cdots,3W-2\}$. After performing the rotating-wave approximation, the full Hamiltonian reads

[TABLE]

where we have introduced $\tilde{\omega}_{k,s}=\omega_{k,s}-\Omega_{0}/2$, and $\mathbf{\tilde{Q}}_{ss^{\prime}}(k)$, which characterizes the strength of the anomalous coupling between two modes. It is obtained from $\mathbf{Q}_{mn}$ through Fourier transform and rotation into the energy eigenbasis (cf. Supplementary Note 3).

As in the one-dimensional waveguide model, a pair of modes is rendered unstable if their detuning $\Delta_{k,ss^{\prime}}=\omega_{k,s}+\omega_{-k,s^{\prime}}-\Omega_{0}$ is smaller than the anomalous coupling between them. The detuning $\Delta_{k,ss^{\prime}}$ varies quickly as a function of $k$ except at points where the slopes of $\omega_{k,s}$ and $\omega_{-k,s^{\prime}}$ coincide to first order, which happens at $k=0,\pi$ when $s=s^{\prime}$. At those values of $k$, the energy matching condition is fulfilled for a broader range of wavevectors, which leads to a larger amplification bandwidth. However, the edge modes are only localized to the edge around $k=\pi$, such that driving around $k=\pi$ is most efficient, which we consider here (cf. Fig. 1). Expanding the dispersion to second order around this point, we find $\omega_{\pi+q}\simeq\omega_{\pi}+q\omega_{\pi}^{\prime}+(q^{2}/2)\omega_{\pi}^{\prime\prime}+\mathcal{O}(q^{3})$, yielding $\Delta_{\pi+q}=2\omega_{\pi}-\Omega_{0}+q^{2}\omega_{\pi}^{\prime\prime}+\mathcal{O}(q^{4})$. Placing the pump at $\Omega_{0}=2\omega_{\pi}$ thus makes magnon pairs around $k=\pi$ resonant, on a bandwidth of order $\sqrt{\mathcal{E}/\omega_{\pi}^{\prime\prime}}$. For weak driving, where the bandwidth is low, higher-order terms in the dispersion relation can be neglected, and this simple calculation captures the amplification behaviour extremely well, as we illustrate in Fig. 2a. We calculate the band structure and find the unstable modes numerically (see Methods), with $\tilde{\omega}_{k,s}$ and $\mathbf{\tilde{Q}}_{ss^{\prime}}$ obtained from a microscopic model detailed in Supplementary Note 2, and plot the resulting band structure with instabilities in Fig. 1.

As we have seen above, an instability requires the anomalous terms to overcome the linear damping and the effective detuning. Linear damping, which we include as a uniform phenomenological parameter $\gamma$, has important consequences, as it sets a lower bound for the amplitude of the electrical field required to drive the system to an instability. It also ensures bulk stability. We have seen that there are three conditions for a parametric instability. First, there has to exist a pair of modes whose lattice momenta add to 0 (or $2\pi$); Second, the sum of their energy has to match the pump frequency; and Third, the strength of their anomalous interaction has to overcome both their detuning and their damping. While momentum and energy matching is by design fulfilled by the edge mode, there is a large number of bulk mode pairs that also fulfil it. We show in Supplementary Note 4 that choosing the polarization to lie along $y$ increases the anomalous coupling for modes with wavevector close to $\pi$ and that the coupling is small for almost all bulk mode pairs. The reason for this is that the bulk modes are approximately standing waves along $y$, and most bulk mode pairs have differing numbers of nodes, such that their overlap averages to zero. The remaining modes with appreciable anomalous coupling a far detuned in energy. This way, robust edge state instability can be achieved without any bulk instabilities, as demonstrated in Fig. 1. Bulk stability is crucial for the validity of the following discussion.

In the presence of an instability, the linear theory predicts exponential growth of edge magnon population. In a real system, the exponential growth is limited by nonlinear damping, for which we introduce another uniform parameter $\eta$, in the same spirit as Gilbert damping in nonlinear Landau-Lifshitz-Gilbert equations Rückriegel et al. (2018), such that Eq. 2 becomes

[TABLE]

Microscopically, such damping arises from the next order in the spin-wave expansion that allows four-wave mixing. While the linear theory only predicts the instability, Eq. 8 predicts a steady-state magnon occupation given through $|\alpha_{\pi+q}|^{2}=\eta^{-1}(\sqrt{4\mathcal{E}^{2}-q^{4}(\omega_{\pi}^{\prime\prime})^{2}}-\gamma)$ (cf. Methods), which we show in Fig. 2b.

II.3 Experimental signatures

TMIs exhibit a magnonic thermal Hall effect at low temperatures Katsura et al. (2010); Onose et al. (2010). A similar effect occurs when the magnon population is not thermal, but a consequence of coherent driving, realizing a driven Hall effect (DHE).

We can calculate the steady-state edge magnon current from the occupation calculated above,

[TABLE]

where $\nu_{\pi+q}\simeq\omega_{\pi}^{\prime}+q\omega_{\pi}^{\prime\prime}$ is the group velocity and $\Lambda=\sqrt[4]{(4\mathcal{E}^{2}-\gamma^{2})/(\omega_{\pi}^{\prime\prime})^{2}}$ is the range over which the steady-state population is finite (which coincides with the range over which the modes become unstable). While the integral can be done exactly (cf. Methods) an approximation within $\pm 5\%$ is given through

[TABLE]

For $2\mathcal{E}\gg\gamma$, a characteristic scaling of steady-state current with driving strength appears, $J_{\mathrm{SS}}\propto\mathcal{E}^{3/2}$, distinct from the linear dependence one would expect for standard heating.

In Fig. 3b,c, we demonstrate that the steady-state edge current depends on the drive strength in a fashion that is well described by Eq. 10. The order-of-magnitude equilibration time can be estimated from the solution to $\dot{\alpha}=(1/2)(\mathcal{E}-\eta|\alpha|^{2})\alpha$, and for $\eta/\mathcal{E}\gg 1$ it evaluates to $t_{\mathrm{eq}}\sim\mathcal{E}^{-1}\log(\mathcal{E}/\eta)\sim 10^{4}J^{-1}$ for our chosen values of $\mathcal{E}$ and $\eta$.

A DHE arises when a rectangular slab of size $L_{x}\times L_{y}$ is driven by a field with a gradient along $y$, as sketched in Fig. 3a. If $L_{x},L_{y}\gg\nu_{\pi}t_{\mathrm{eq}}$, the edges equilibrate to a steady-state magnon population governed by Eq. 10. The difference between the steady-state magnon currents on top and bottom edge corresponds to a net energy current $J_{\mathrm{net}}^{x}$ along $x$, which to first order in the drive strength difference $\Delta_{y}\mathcal{E}$ can be written Vinkler-Aviv and Rosch (2018)

[TABLE]

where one should note that in this non-equilibrium setting $\kappa_{xy}$ is not a proper conductivity as in conventional linear response. The net edge current causes one side of the system to heat up faster, resulting in a temperature difference transverse to the gradient. As the edge magnons decay along the edge, the reverse heat current is carried by bulk modes. For small temperature differences the heat current follows the temperature gradient linearly and thus the averaged temperature difference $\Delta_{x}T=\int dy\,[T(L_{x},y)-T(0,y)]/L_{y}=J_{\mathrm{net}}^{x}/\kappa_{xx}$. The temperature difference can thus be written in terms of the applied field strength difference

[TABLE]

As a word of caution, we note that this relation relies on several key assumptions. To begin with, temperature is in fact not well defined along the edge, as there is a non-equilibrium magnon occupation. Edge magnons decay at a certain rate into phonons, which can be modelled as heating of the phonon bath. If the equilibration time scale of the latter is fast compared to the heating rate through magnon decay, one can at least associate a local temperature to the phonons. Similarly, the bulk magnon modes can be viewed as a fast bath for the magnon edge mode and similar considerations apply. Even if these assumptions are justified, the two baths do not need to have the same temperature. Next, the heat conductivity associated to magnons and phonons differ in general, such that the $\kappa_{xx}$ appearing in Eq. 12 can only be associated with the bulk heat conductivity if the temperatures of the two baths are equal. Some of these complications have been recognized to also play an important role in measurements of the magnon thermal Hall effect Vinkler-Aviv and Rosch (2018).

While the above-mentioned concerns make quantitative predictions difficult, the DHE is easily distinguishable from the thermal Hall effect, due to the strong dependence of the temperature difference $\Delta_{x}T$ on drive frequency and polarization, as well as the fact that below the cutoff $2\mathcal{E}=\gamma$ no instability occurs and that $J_{\mathrm{net}}\propto\sqrt{\mathcal{E}_{\mathrm{avg}}}$ for $2\mathcal{E}_{\mathrm{avg}}\gg\gamma$, rather than the linear dependence one would expect from standard heating. In certain materials such as Cu(1,3-bdc), the appearance or disappearance of the topological edge modes can be tuned with an applied magnetic field Chisnell et al. (2015), a property that could be used to further corroborate the results of such an experiment.

A number of other experimental probes might be used to certify a large edge magnon current and thus the presence of edge states. On the one hand, with a large coherent magnon population in a given mode, the local magnetic field and electric polarization associated to that mode will be enhanced. In particular techniques that directly probe local magnetic or electric fields, such as neutron scattering Chisnell et al. (2015); Yao et al. (2018) or x-ray scattering, which to date are not powerful enough to resolve edge modes Chisnell et al. (2015), would thus have a coherently enhanced signal, for example by almost two orders of magnitude when taking the conservative parameters in Fig. 2. On the other hand, heterostructures provide a way to couple the magnons out of the edge mode into another material Rückriegel et al. (2018), for example one with a strong spin Hall effect, in which they can be detected more easily. In this setup, again the fact that the edge magnons have a large coherent population should make their signal easily distinguishable from thermal noise.

II.4 Material realizations

The model of a kagome lattice ferromagnet with DM interaction has been found to describe the thermal magnon Hall effect in Lu2V2O7 Onose et al. (2010), as well as the bulk magnon band structure of Cu(1,3-bdc) Nytko et al. (2008); Chisnell et al. (2015). These materials are in fact 3D pyrochlore lattices, which can be pictured as alternating kagome and triangular lattices along the [111] direction. However, their topological properties can be captured by considering only the kagome planes Zhang et al. (2013); Mook et al. (2014); Chisnell et al. (2015) (shown in Fig. 1), thus neglecting the coupling between kagome and triangular planes. It has been suggested that the effect of the interaction may be subsumed into new effective interaction strengths Mook et al. (2014) or into an effective on-site potential Zhang et al. (2013). Typical values for strength of the DM and Heisenberg interactions lie between $|\mathbf{D}|/J\approx 0.18$ Chisnell et al. (2015), $J\approx 0.6\pm 0.1\,\text{meV}\simeq 150\pm 30\,\text{GHz}$ Chisnell et al. (2015) in Cu(1,3-bdc) and $|\mathbf{D}|/J\approx 0.32$ Onose et al. (2010), $J\approx 3.4\,$meV$\,\approx 0.82\,$THz Zhang et al. (2013) in Lu2V2O7. The energy of the edge states close to $k=\pi$ is approximately $J$, such that the applied drive needs to be at a frequency $\omega_{0}/2\pi=0.3$–$1.6\,$THz. While experimentally challenging, low THz driving down to 0.6 THz has recently been achieved Karch et al. (2011); Plank et al. (2018). Furthermore, the magnon energy can be tuned by applied magnetic fields.

An instability requires $E_{0}ae(t/U)^{3}\gtrsim\gamma$. With $a\simeq 10\,$Å Zhang et al. (2013), $J\simeq 1\,$meV, and assuming $t/U\simeq 0.1$, we can estimate the minimum field strength required to overcome damping $\gamma_{k}\simeq 10^{-4}J$ to be $E_{0}\simeq 10^{5}\,$V/m, although for quantitative estimates one would require both accurate values for the damping of the edge modes (at zero temperature) and $t/U$. This is accessible in pulsed operation Karch et al. (2011); Plank et al. (2018), and perhaps in continuous operation through the assistance of a cavity.

Since the qualitative behaviour we describe can be derived from general and phenomenological considerations, we expect it to be robust and present in a range of systems, as long as they allow for anisotropy, i.e., if bonds are not centres of inversion. We thus expect that topological magnon amplification is also possible in recently discovered topological honeycomb ferromagnet CrI3 Chen et al. (2018).

III Discussion

We have shown that appropriate electromagnetic driving can render topological magnon edge modes unstable, while leaving the bulk modes stable. The resulting non-equilibrium steady state has a macroscopic edge magnon population. We present several strategies to certify the topological nature of the band structure, namely, implementing a driven Hall effect (DHE), direct detection with neutron scattering, or by coupling the magnons into a material with a spin Hall effect.

Our work paves the way for a number of future studies. As we have pointed out, edge mode damping plays an important role here. One might expect their damping to be smaller than that of generic bulk modes as due to their localization they have a smaller overlap to bulk modes. This suppression should be compounded by the effect of disorder Rückriegel et al. (2018), which may further enhance the feasibility of our proposed experiments. On the other hand, rough edges will have an influence over the matrix element between drive and edge modes, leading to variations in the anomalous coupling strength. Phonons in the material are crucial for robust thermal Hall measurements Vinkler-Aviv and Rosch (2018) and could possibly mix with the chiral magnon mode Thingstad et al. (2019), which motivates full microscopic calculations.

An exciting prospect is to use topological magnon amplification in magnon spintronics. There have already been theoretical efforts studying how magnons can be injected into topological edge modes with the inverse spin Hall effect Rückriegel et al. (2018). Given an efficient mechanism to couple magnons into and out of the edge modes, our amplification mechanism may enable chiral travelling-wave magnon amplifiers, initially proposed in photonic crystals Peano et al. (2016). Even when simply seeded by thermal or quantum fluctuations, the large coherent magnon steady state could power topological magnon lasers Harari et al. (2018), with tremendous promise for future application in spintronics. In the near future, we hope that topological magnon amplification can be used for an unambiguous discovery of topological magnon edge modes.

Acknowledgements.

We are grateful to Ryan Barnett, Derek Lee, Rubén Otxoa, Pierre Roy, and Koji Usami for insightful discussions and helpful comments. DM acknowledges support by the Horizon 2020 ERC Advanced Grant QUENOCOBA (grant agreement 742102). AN holds a University Research Fellowship from the Royal Society and acknowledges support from the Winton Programme for the Physics of Sustainability and the European Union’s Horizon 2020 research and innovation programme under grant agreement No 732894 (FET Proactive HOT).

Methods

Numerical Calculation

For the numerical calculation, we choose a manifestly inversion-symmetric system obtained by deleting the lowest row of sites, a situation that is depicted in Fig. 1, where the tip of the lowest blue triangle is part of a unit cell whose other sites are not included. For example, repeating the star shown in Fig. 1 along $x$ would result in an inversion-symmetric strip with $W=3$. A Fourier transform of Eq. 5 along $x$ yields a $3W-2$ by $3W-2$ Hamiltonian matrix for each momentum $k$

[TABLE]

Note that we take $\hbar=1$. Diagonalizing this matrix yields single-particle energy eigenstates with annihilation operator $b_{k,s}$, and a Hamiltonian $H_{0}=\sum_{k,s}\omega_{k,s}b_{k,s}^{\dagger}b_{k,s}$. The resulting band structure is shown in Fig. 1. In our convention, the lowest bulk band has Chern number $\operatorname{sgn}D_{z}$, the middle bulk band [math] and the top bulk band $-\operatorname{sgn}D_{z}$ (calculated, e.g., through the method described in Ref. Fukui et al., 2005). Accordingly, there is one pair of edge modes in each of the bulk gaps, one right-moving localized at the lower edge and one left-moving at the upper.

Including the anomalous terms obtained from a calculation based on the Fermi-Hubbard model at half filling yields the full Hamiltonian Eq. 7. By means of a Bogoliubov transformation we obtain the magnon band structure and the unstable states Blaizot and Ripka (1986), which form the basis for Fig. 1 and Fig. 2. The inclusion of nonlinear damping yields Eq. 8, which has been used to calculate Fig. 3. In the end, we calculate the current by evaluating the expectation value of the particle current or energy current operator, which can be obtained for a given bond from the continuity equation Bernevig and Hughes (2013). The current across a certain cut of the system is obtained by summing the current operators for all the bonds that cross it. As the system we study is inversion symmetric, the total current in the $x$ direction vanishes. In order to specifically find the edge current, we thus define a cut through half of the system, for example from the top edge to the middle.

Unstable modes in Bogoliubov-de Gennes equation

We consider the full Hamiltonian

[TABLE]

$H_{0}-H_{\text{rot}}$ gives rise to the band structure shown in Fig. 1 above, while $H_{\text{amp}}$ contains the anomalous terms. The idea of this section is to calculate which modes in Eq. 14 are unstable. Ideally, those should be the relevant edge modes, and only those. It turns out that this is possible in presence of linear damping.

Following Ref. Peano et al., 2016, we define the vector $|a_{k}\rangle=(a_{k,1},a_{k,2},\dots,a_{k,N},a_{k_{0}-k,1}^{\dagger},\dots,a_{k_{0}-k,N}^{\dagger})^{T}$, where the index combines the label $l_{y}$ and the site label in the unit cell and therefore runs from $1$ to $N=3W-2$. The Hamiltonian can generically be written

[TABLE]

where $\mu_{k}$ originates from $H_{0}-H_{\text{rot}}$ and $\nu_{k}$ from $H_{\text{amp}}$. This form makes it evident that $\mu_{k}=\mu_{k}^{\dagger}$ and $\nu_{k}=\nu_{k}^{T}$. The equation of motion for this vector can be found from the Hamiltonian above and is

[TABLE]

with $\sigma_{z}=\operatorname{diag}(1,1,\dots,-1,-1,\dots)$ ($N$ “$+1$”s and $N$ “$-1$”s), with

[TABLE]

We can then solve the eigenvalue problem and find stable and unstable modes. Furthermore, we can find the time-evolution for operators in the Heisenberg picture from Eq. 16. It is simply $|a_{k}(t)\rangle=e^{-i\sigma_{z}h_{k}t}|a_{k}(0)\rangle$.

Steady state of nonlinear equations of motion

We start from the equations of motion (8) given in the main text, repeated here for convenience

[TABLE]

In the steady state, $|\alpha_{k}|^{2}=\text{const.}$, so we use the ansatz $\alpha_{k}=\exp(i\Delta t)\bar{\alpha}$, and $\alpha_{-k}^{*}=\exp(i\Delta t)z\bar{\alpha}$ for some complex numbers $z$, $\bar{\alpha}$ and real frequency $\Delta$. As we are only interested in a narrow range of momenta, we expand the dispersion relation to second order, as in the main text. The pump frequency is set to match the edge mode at $k=\pi$, i.e., $\Omega_{0}=2\omega_{\pi}$. As a consequence, $\tilde{\omega}_{\pi+q}=q\omega_{\pi}^{\prime}+q^{2}\omega_{\pi}^{\prime\prime}/2+\mathcal{O}(q^{3})$.

If there is an instability, the solution $\bar{\alpha}=0$ is unstable. Assuming $\bar{\alpha}\neq 0$ (thus $z\neq 0$), and for $\Delta=q\omega_{\pi}^{\prime}$, we find the set of equations

[TABLE]

Multiplying the second equation by $|z|^{2}$, and subtracting the complex conjugate of the resulting equation from the first equation, one can show that $|z|^{2}=1$. With this condition Eqs. 19 and 20 coincide, such that we can solve them for the intensity

[TABLE]

This equation has solutions if and only if $4\mathcal{E}^{2}\geq q^{4}(\omega_{\pi}^{\prime\prime})^{2}+\gamma^{2}$, which coincides with the condition for the instability. If this condition is fulfilled, we have

[TABLE]

The steady-state edge magnon current

[TABLE]

where $F(k,m)$ is the elliptic integral of the first kind.

Particle current operator

The particle current operator is obtained from the continuity equation for the number of magnons. We have

[TABLE]

where $h_{n}$ are local Hamiltonians defined through

[TABLE]

The second term in Eq. 24 can be interpreted as a sum of the particle currents from $n$ to the neighbouring sites $m$.

Supplementary Note 1: Effective spin Hamiltonian from Fermi-Hubbard model

In order to support our qualitative analysis above, we derive the polarization tensor in a kagome TMI from a Fermi-Hubbard model at half filling, with an on-site Coulomb repulsion $U$ much larger than the hopping $t$, $t/U\ll 1$. As is well known, the low-energy physics can be described by perturbing around the Mott insulator state MacDonald et al. (1988). A contribution to the polarization arises to third order in the hopping Bulaevskii et al. (2008) (hopping around a triangle), which is derived in detail below.

Following Zhu et al. Zhu et al. (2014), we consider a one-band Hubbard model with SOC

[TABLE]

with $\mathbf{c}_{i}=(c_{j\uparrow},c_{j\downarrow})^{T}$, and SOC vector $\mathbf{d}_{ij}=i\mathbf{n}_{ij}t_{ij}\sin(\theta_{ij})$. The parameterization of $\mathbf{d}$ in terms of a unit vector $\mathbf{n}_{ij}$ and an angle $\theta_{ij}$ will become useful later.

Zhu et al. Zhu et al. (2014) derive the effective low-energy spin Hamiltonian to second order in the hopping, which is

[TABLE]

where $S_{j}^{\nu}=\tfrac{1}{2}\mathbf{c}_{j}^{\dagger}\cdot\sigma^{\nu}\cdot\mathbf{c}_{j}$, and $\mathcal{J}_{ij}$ is the exchange tensor pertaining to bond $\langle ij\rangle$ and can be written

[TABLE]

The three terms give rise to the isotropic Heisenberg interaction, to asymmetric, and to symmetric exchange anisotropy, respectively.

Supplementary Note 2: Polarization tensor

The direction in which the polarization may point is constrained in the same way as the DM vector associated to each bond [cf. Eq. 4]. The reflection symmetry around the plane orthogonal to each bond constrains the vector to lie in this symmetry plane. In addition, the lattice is three-fold rotation symmetric, as well as inversion symmetric around lattice sites, such that the direction of one vector determines that of all others. The precise direction and magnitude of the vector may be obtained, for example, from perturbation theory in the Fermi-Hubbard model at half filling, as we demonstrate below. The results in the main text in principle only require that the anisotropic part of the polarization tensor is nonzero, which is allowed whenever bonds are not centres of inversion.

Microscopically, the anisotropy is due to spin-orbit coupling (SOC). We follow Zhu et al. Zhu et al. (2014), who derive the electric polarization as a third-order hopping process, which is the lowest-order relevant contribution. It is given through Zhu et al. (2014)

[TABLE]

where $k$ is the third site in the loop and with

[TABLE]

The vector $\mathbf{p}_{0,ij}$ points into the triangle, orthogonal to the bond $\langle ij\rangle$ and in the plane of the triangle. Importantly, this means that, when following bonds along a straight line, their polarization changes sign from bond to bond.

The angle $\theta$ parametrizes the relative strength of the SOC. To first order in $\theta$, the only scalar quantity one can construct with one vector are of the form $\mathbf{n}\cdot(\mathbf{S}\times\mathbf{S})$, which does not have the form we are interested in. Hence we expand to second order

[TABLE]

with

[TABLE]

Physically, in the original Hubbard Hamiltonian, only the SOC term can generate spin flips, which is why we need to go to second order in the SOC to obtain anomalous pairing terms that generate two magnons from one photon. In the spin wave picture, terms such as $(\mathbf{S}_{i}\cdot\mathbf{n}_{ij})(\mathbf{S}_{j}\cdot\mathbf{n}_{ij})$ and $(\mathbf{S}_{i}\cdot\mathbf{n}_{ij})(\mathbf{S}_{j}\cdot\mathbf{n}_{ijk})$ can lead to anomalous terms, which can lead to instabilities and thus amplification. In order to make progress, we need to apply the general model Eq. 27 to our particular problem. Note that while Eq. 27 predicts a positive $J$, experiments show that $J$ is in fact negative. This is a result of other contributions, such as exchange. Thus, the measured $J$ cannot be used to determine the angle $\theta_{ij}$ in Eq. 27. Instead, the angle needs to be fitted independently, or determined from a measurement of the spin-orbit effect. Comparing Eq. 27 to Eq. 4, we can identify

[TABLE]

where $\mathbf{D}_{ij}$ is the vector in Eq. 4, $\tan(2\theta_{ij})=-|\mathbf{D}_{ij}|/J_{\mathrm{SE},ij}$ quantifies the strength of the DM interaction relative to the Heisenberg coupling from superexchange $J_{\mathrm{SE}}$, and we have used the subscript SE to denote the superexchange contribution. In principle, all these quantities can differ from site to site, but here we study a translation-invariant Hamiltonian, which simplifies the description considerably.

By lattice symmetry, $\mathbf{D}$ has to lie in the plane orthogonal to the bonds (since that is a symmetry plane). In the pyrochlore lattice, each bond is part of two triangles. The net DM interaction is the sum of the contribution from each triangle. If we consider the corner-sharing cube that surrounds the tetrahedron, the DM vector lies in the plane of the cube face that also encompasses the bond, as derived for instance in Ref. Elhajal et al., 2005. If we choose the upright triangles in Fig. 1 to be part of tetrahedra pointing into the plane (and thus the upside-down triangles are part of tetrahedra pointing out of the plane), and consider the bond lying along $x$ in an upright triangle, we have $\mathbf{n}_{12}=-(\sqrt{2/3})\mathbf{\hat{z}}-1/\sqrt{3}\mathbf{\hat{y}}$ ($\mathbf{\hat{z}}$ points out of the plane, i.e., our coordinate system is right handed with Fig. 1 being the $xy$-plane, with $x$ being horizontal and $y$ vertical). The DM vectors for the other bonds in the triangle can be obtained through rotation by $2\pi/3$ around $z$. The DM vectors in upside-down triangle then follow from reversing the vectors in the upright triangle ($\mathbf{v}\to-\mathbf{v}$). This argument assumes ordered bonds (here counterclockwise in all triangles).

This determines $\mathbf{n}_{ij}$ and $\theta_{ij}\equiv\theta=(1/2)\tan^{-1}(D/J_{\mathrm{SE}})$. The spin-orbit contribution is assumed to be weak, such that $\theta$ is small. In analogy to a charged particle picking up a U(1) phase when hopping in a loop penetrated by a magnetic field, $\theta_{ijk}$ and $\mathbf{n}_{ijk}$ parameterize the SU(2)-phase that is picked up by the electron spin when hopping around the triangle Zhu et al. (2014). Writing the hopping part of the original Hubbard Hamiltonian

[TABLE]

we can identify $\mathcal{A}_{ij}\equiv\exp(i\theta_{ij}\mathbf{n}_{ij}\cdot\bm{\sigma}).$ The lowest order contribution to the polarization comes from a third order hopping process around a triangle, during which an electric spin picks up the total rotation

[TABLE]

This defines $\theta_{ijk}$ and $\mathbf{n}_{ijk}$. Due to translation and rotation symmetries, $\theta_{ijk}$ is the same for all bonds, and given through

[TABLE]

The sign is ambiguous, and we have chosen $\theta_{ijk}>0$ in the second equality. The vector $\mathbf{n}_{ijk}$ depends on the bond we consider. For the bond connecting site 1 and 2 in the same unit cell (i.e., the lower edge in an upright triangle), we have

[TABLE]

The vectors $\mathbf{n}_{231}$, $\mathbf{n}_{312}$ can be obtained from $\mathbf{n}_{123}$ through rotation by $2\pi/3$ and $4\pi/3$ around $z$.

Terms such as $\mathbf{S}_{i}\cdot\mathbf{S}_{j}$, $S_{i}^{z}S_{j}^{z}$, $\mathbf{S}_{i}\times\mathbf{S}_{j}$ cannot change the angular momentum along $z$ and thus do not lead to anomalous terms. In the second order (in $\theta$) contribution to the polarization [cf. Eq. 32], we have two promising terms. The second, however, yields

[TABLE]

and thus does not contribute to second order. The remaining term is

[TABLE]

where we have subtracted the component of the vector $\mathbf{n}_{ij}$ along $z$, because it does not lead to anomalous terms.

Supplementary None 3: Amplification Hamiltonian

Independent of whether we justify the existence of an anomalous term via symmetry considerations (Section II.2) or the microscopic derivation (Supplementary Note 2), the amplification Hamiltonian takes same form, due to symmetry constraints. Note that due to symmetry, $\mathbf{P}_{ij}$ is constrained to lie in the plane perpendicular to the bond. As the $z$ component is irrelevant, the resulting amplification Hamiltonian depends only on the modulus of the in-plane component of the polarization $\mathbf{P}$. Taking only the relevant term, the amplification Hamiltonian is written as (note the slightly odd convention where by addition of the indices, such as $m+n$ or $1+2$, we mean that we take the vectors to those sites and add them)

[TABLE]

where $\mathbf{Q}_{mn}=\mathbf{n}_{mn}-\hat{\mathbf{z}}n_{mn}$ is perpendicular to the bond and points outside for upright triangles and inside for upside-down triangles. In fact, it is irrelevant whether $\mathbf{Q}_{mn}$ points in or out, since $\mathbf{Q}_{mn}\to-\mathbf{Q}_{mn}$ leaves Eq. 40 unchanged. In the spin-wave picture, we have $(\mathbf{S}_{m}\cdot\mathbf{Q}_{mn})(\mathbf{S}_{n}\cdot\mathbf{Q}_{mn})=S^{-}_{m}S^{-}_{n}(Q_{mn}^{+})^{2}/4+\text{H.c.}+\cdots$, where $Q^{\pm}\equiv Q^{x}\pm iQ^{y}$. Since $\mathbf{Q}_{\mathrm{mn}}$ points out of the upward facing triangles (they are parallel to $\mathbf{n}_{mn}$), we have $Q^{+}_{12}=e^{i\pi/6}/\sqrt{3}$, $Q^{+}_{23}=e^{5i\pi/6}/\sqrt{3}$, and $Q^{+}_{31}=-i/\sqrt{3}$. We end up with

[TABLE]

where we have used $(Q_{12}^{+})^{2}=\tfrac{1}{3}e^{i\pi/3}$, $(Q_{23}^{+})^{2}=\tfrac{1}{3}e^{-i\pi/3}$, $(Q_{31}^{+})^{2}=-\tfrac{1}{3}$. The minus sign between the two terms in the round and square brackets above stems from the fact that the induced polarization switches sign going from a bond to an adjacent one. $\Omega_{0}$ is the frequency of the incoming radiation, $\mathbf{E}_{0}=\hat{\mathbf{e}}E_{0}$ its polarization and amplitude. If the radiation is polarized along $z$, at least to this order in perturbation theory, it has no effect on the TMI, thus we choose it to lay in the plane.

Recall $\mathbf{p}_{0,ij}=8ea\frac{t_{ij}t_{jk}t_{ki}}{U^{3}}(2\bm{\rho}_{k}-\bm{\rho}_{i}-\bm{\rho}_{j})$ (where $i,j,k\in\{1,2,3\}$ and are all distinct). Then $2\bm{\rho}_{1}-\bm{\rho}_{2}-\bm{\rho}_{3}=(3,-\sqrt{3})/4$, $2\bm{\rho}_{2}-\bm{\rho}_{1}-\bm{\rho}_{3}=(0,\sqrt{3})/2$ and $2\bm{\rho}_{3}-\bm{\rho}_{1}-\bm{\rho}_{2}=-(3,\sqrt{3})/4$. We further define

[TABLE]

proportional to the strength of the electric field, and go into the rotating frame with respect to the Hamiltonian $H_{\text{rot}}=\frac{\Omega_{0}}{2}\sum_{\alpha,l_{y},p_{x}}a_{\alpha,p_{x},l_{y}}^{\dagger}a_{\alpha,p_{x},l_{y}}.$ We arrive at

[TABLE]

From this form it is clear that the terms at $\pm k$ couple, so that it is we should combine negative and positive momenta. Finally, choosing $\mathbf{\hat{e}}=\mathbf{\hat{y}}$, this leads to

[TABLE]

Supplementary Note 4: Amplification matrix element

In the section above we have derived the amplification Hamiltonian

[TABLE]

where here the generic indices $s,s^{\prime}$ contain both the site label $\alpha$ and the unit cell label $l_{y}$.

Diagonalizing the bilinear undriven Hamiltonian $H_{0}=\sum_{k}\mathbf{a}_{k}^{\dagger}\mathsf{\mu}_{k}\mathbf{a}_{k}=\sum_{k}\mathbf{b}_{k}^{\dagger}\mathsf{\omega}_{k}\mathbf{b}_{k}$, where $\mathbf{b}_{k}=\mathsf{U}_{k}^{\dagger}\mathbf{a}_{k}$ is a vector containing the annihilation operators of the energy eigenmodes and $\omega_{k}$ is a diagonal matrix. Writing Eq. 45 in terms of energy eigenstates, we obtain

[TABLE]

We can investigate the coupling strength between the various modes numerically, as is done in Fig. 4. The first conclusion, when considering the coupling matrix in the energy eigenbasis for wavevectors close to $\pi$ is that the anomalous coupling in between the edge modes is among the largest. Comparable coupling strength is only achieved in between modes in differing bulk bands, as is seen from the diagonal lines in the off-diagonal blocks. This can be appreciated by thinking about the form of the bulk wavefunctions along $y$, which are approximately standing waves with [math] to $N_{y}-1$ nodes. Since the matrix element between two bulk modes contains their product (with a constant applied field), summed over $y$, bulk modes with differing numbers of nodes approximately sum to zero. In between bands, the number of nodes within a unit cell changes, such that a full cancellation no longer occurs.

We next plot the maximum coupling strength between any of the modes as a function of wavevector. From this plot we conclude that the anomalous coupling is most efficient around $k\approx\pi$. This result can be understood to some degree by looking at the form of the amplification Hamiltonian Eq. 43. Choosing the polarization of the applied field to lie along $y$, the first term coupling sites $1$ and $3$ is dominant. In Fourier space this term has the functional form $\sin(k)$, such that it is largest around $\pi$, which roughly matches the shape in Fig. 4. This conclusion is strongly dependent on the polarization we choose for the applied field. We can plot the same quantities for a polarization along $x$, which turns off the coupling between sites 1 and 3. In this case the maximum coupling strength no longer lies around $k=\pi$, which is plotted in Fig. 4. Finally, this demonstrates one of the reasons why the agreement between the chiral waveguide model and the microscopic two-dimensional model is so good, namely that the matrix element is near unity (in units of $2\mathcal{E}$).

As we emphasize in the main text, the anomalous coupling strength is only one of the factors that influence whether a mode pair would become unstable under driving. For example, all bulk mode pairs close to $k=\pi$ are far detuned in energy and thus cannot become unstable, regardless of the strength of their anomalous coupling.