# Phonon Hall Viscosity in Magnetic Insulators

Mengxing Ye  Kavli Institute for Theoretical Physics, University of
California, Santa Barbara, CA 93106, USA  Lucile Savary  Université de Lyon,
École Normale Supérieure de Lyon, Université Claude Bernard Lyon I, CNRS,
Laboratoire de physique, 46, allée d’Italie, 69007 Lyon, France  Leon Balents
Kavli Institute for Theoretical Physics, University of California, Santa
Barbara, CA 93106, USA

###### Abstract

The Phonon Hall Viscosity is the leading term evincing time-reversal symmetry
breaking in the low energy description of lattice phonons. It may generate
phonon Berry curvature, and can be observed experimentally through the
acoustic Faraday effect and thermal Hall transport. We present a systematic
procedure to obtain the phonon Hall viscosity induced by phonon-magnon
interactions in magnetic insulators under an external magnetic field. We
obtain a general symmetry criterion that leads to non-zero Faraday rotation
and Hall conductivity, and clarify the interplay between lattice symmetry,
spin-orbit-coupling, external magnetic field and magnetic ordering. The
symmetry analysis is verified through a microscopic calculation. By
constructing the general symmetry-allowed effective action that describes the
spin dynamics and spin-lattice coupling, and then integrating out the spin
fluctuations, the leading order time-reversal breaking term in the phonon
effective action, i.e. the phonon Hall viscosity, can be obtained. The
analysis of the square lattice antiferromagnet for a cuprate Mott insulator,
Sr2CuO2Cl2, is presented explicitly, and the procedure described here can be
readily generalized to other magnetic insulators.

Introduction.— In a material without time-reversal symmetry, the elastic
stress may contain a non-dissipative term proportional to the time derivative
of the elastic strain, the coefficient of which defines the phonon Hall
viscosity (PHV) Barkeshli  _et al._ (2012); Tuegel and Hughes (2017);
Shapourian  _et al._ (2015). The PHV is analogous to the Hall viscosity in a
fluid, which has been heavily studied in quantum Hall states Avron  _et al._
(1995); Avron (1998); Read (2009); Read and Rezayi (2011), and
superconductors/superfluidsRead (2009); Read and Rezayi (2011); Bradlyn  _et
al._ (2012). The PHV can be represented as a momentum space gauge field for
phonons, and as such encodes the Berry curvature of phonon eigenstates.
Multiple mechanisms can generate a PHV Sytcheva  _et al._ (2010); Qin  _et
al._ (2012); Vinkler-Aviv and Rosch (2018); Zhang  _et al._ (2019); Chen  _et
al._ (2020); Ye  _et al._ (2020), including subtle band effects Saito  _et
al._ (2019). In this paper, we consider the PHV generated by the
magnetoelastic coupling of strain and spins Sytcheva  _et al._ (2010); Ye  _et
al._ (2020), taking into account exchange and the Zeeman interaction of the
spins with an external magnetic field.

The phonon Hall viscosity may be probed in experiments such as acoustic
Faraday rotation and thermal Hall transport. The first is a measure of broken
degeneracy between right and left circularly polarized acoustic phonons; the
Faraday rotation angle is proportional to the PHV Boiteux  _et al._ (1971);
Sytcheva  _et al._ (2010); Tuegel and Hughes (2017). The second describes a
transverse heat current in response to a thermal gradient, and the part of
such heat current carried by phonons – the phonon Hall effect – is
proportional to the PHV (in the ballistic phonon regime). Recent thermal Hall
measurements in cuprate compounds Grissonnanche  _et al._ (2019); Boulanger
_et al._ (2020); Grissonnanche  _et al._ (2020) suggest that the phonon Hall
effect may dominate the thermal Hall signal in both the Mott insulating and
pseudogap phases of cuprates.

In this paper, we present a systematic procedure to obtain the PHV in magnetic
insulators, based on the coupling of spins to elastic strains. We concentrate
mainly on the PHV terms that contribute to the acoustic Faraday rotation and
phonon Hall conductivity. We first introduce the effective field theory and
the expression for the viscosity coefficients. We follow with a thorough
symmetry analysis in different phases and regimes, arriving at a symmetry
criterion for a nonzero PHV, which we verify through a microscopic
calculation. We conclude with a discussion of experimental implications, spin-
orbit coupling, and open questions.

Phonon effective action and magnetoelastic coupling—The PHV term in the
effective linear response action is Avron (1998); Barkeshli  _et al._ (2012),

$\displaystyle\mathcal{S}_{\rm
PHV}=\frac{1}{2}\sum_{{\bm{q}},\Gamma,\Gamma^{\prime}}\int\mathop{}\\!\mathrm{d}t\,\,\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})(\mathcal{E}_{\Gamma,{\bm{q}}}\dot{\mathcal{E}}_{\Gamma^{\prime},-{\bm{q}}}-\dot{\mathcal{E}}_{\Gamma,{\bm{q}}}\mathcal{E}_{\Gamma^{\prime},-{\bm{q}}}),$
(1)

where $\mathcal{E}_{ij}=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i})$ is
the strain tensor, $i,j$ label Euclidean coordinates, and ${\bm{u}}$ is the
lattice displacement field. Here $\dot{\mathcal{E}}=\partial_{t}\mathcal{E}$,
and $\Gamma,\Gamma^{\prime}$ are the irreducible representations (irreps) of a
group whose precise definition we make clear below. The Hall viscosity
$\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})$ describes the stress
$\varSigma_{\Gamma}$ resulting from a time-dependent strain
$\mathcal{E}_{\Gamma^{\prime}}$, i.e.
$\varSigma_{\Gamma}\sim\delta\mathcal{S}_{\rm
PHV}/\delta{\mathcal{E}_{\Gamma}}\sim\eta_{\Gamma\Gamma^{\prime}}\dot{\mathcal{E}}_{\Gamma^{\prime}}$
Read and Rezayi (2011). As we will discuss at length below, symmetry
constrains greatly which $\Gamma\Gamma^{\prime}$ lead to a nonvanishing
$\eta^{H}_{\Gamma\Gamma^{\prime}}$. We allow for a momentum dependence of the
Hall viscosity
$\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})=-\eta^{H}_{\Gamma^{\prime}\Gamma}(-{\bm{q}})$.111
Symmetry allows additional PHV terms which are not expressed solely in terms
of the strain but also in terms of the time derivative of the rotation
$\mathcal{M}_{ij}=\frac{1}{2}(\partial_{i}u_{j}-\partial_{j}u_{i})$. Such
terms, however, cannot be induced by conventional magnetoelastic coupling
which involves strain only.

The Lagrangian density for the magnetoelastic coupling can be expressed as

$\displaystyle\mathcal{L}_{sl}=\sum_{\Gamma,a}\lambda_{\Gamma,a}\,\mathcal{E}_{\Gamma}({\bm{x}},\tau)\cdot\mathcal{O}_{\Gamma,a}({\bm{x}},\tau),$
(2)

where $\lambda_{\Gamma,a}$ denotes the magnetoelastic coupling strength in
units of energy to the $\Gamma$th component of the strain field
$\mathcal{E}_{\Gamma}$ of (each independent copy $a$ of) the spin operator
composite $\mathcal{O}_{\Gamma,a}({\bm{x}},\tau)$ which transforms under the
same $\Gamma$ irrep.

The Hall viscosity itself is then obtained by extracting terms linear in the
phonon energy $\omega$ (reflecting the fact that the PHV is time-reversal
odd). More precisely, the Hall viscosity “coefficient”
$\eta^{H}_{\Gamma\Gamma^{\prime}}$ is obtained from the correlation of
$\mathcal{O}_{\Gamma}$, $\mathcal{O}_{\Gamma^{\prime}}$ as Ye  _et al._
(2020):

$\frac{i\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})}{\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}=\lim_{\omega\rightarrow
0}\partial_{\omega}\sum_{a,a^{\prime}}\langle\mathcal{O}_{\Gamma
a}(-{\bm{q}},\tau)\cdot\mathcal{O}_{\Gamma^{\prime}a^{\prime}}({\bm{q}},0)\rangle_{-i\omega+0^{+}},$
(3)

where
$\langle\mathcal{A}\rangle_{\omega_{n}}\equiv\int_{0}^{\beta}\mathop{}\\!\mathrm{d}\tau
e^{i\omega_{n}\tau}\langle T_{\tau}\mathcal{A}(\tau)\rangle$ denotes the
Fourier transform, $\tau$ is imaginary time, $\omega_{n}=2\pi n/\beta$ is the
$n$th bosonic Matsubara frequency and $T_{\tau}$ denotes imaginary-time
ordering.

Magnetic system.—We now specify the theory to calculate the correlations of
the $\mathcal{O}$ operators. We study a two-dimensional spin system which
orders antiferromagnetically in zero field. Because we will compute
$\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})$ for small momentum ${\bm{q}}$ at
low temperature (compared to the Debye temperature and magnon bandwidth), it
is enough to consider low-energy spin and lattice fluctuations. We therefore
employ a path integral approach, and describe the low energy dynamics in terms
of continuous fields: the irrep components of the strain tensor
($\mathcal{E}_{\Gamma}$) and the staggered (${\bm{n}}$) and uniform
(${\bm{m}}$) magnetizations. The latter two fields are described by the well-
known nonlinear sigma model (see SM), with the constraints $|\bm{n}|^{2}=1$
and $\bm{n}\cdot\bm{m}=0$.

The $\mathcal{O}$ correlations are obtained from the nonlinear sigma model as
follows. We will consider a weak applied field along the $z$-axis,
$\bm{h}=h\bm{\hat{z}}$. We assume that in the zero-field limit the uniform
magnetization vanishes, $\lim_{h\rightarrow 0}\langle\bm{m}\rangle=0$, and the
Néel vector $\bm{n}$ lies along the $x$-axis, $\lim_{h\rightarrow
0}\langle\bm{n}\rangle=n_{0}\bm{\hat{x}}$. Then constraints are solved via
$\bm{n}=(n_{0},\bm{n}_{\scriptscriptstyle\perp})$, with
$n_{0}=\sqrt{1-\bm{n}_{\scriptscriptstyle\perp}^{2}}$,
$\bm{m}=(m_{x},\bm{m}_{\scriptscriptstyle\perp})$ and
$m_{x}=-\bm{n}_{\scriptscriptstyle\perp}\cdot\bm{m}_{\scriptscriptstyle\perp}/n_{0}$.
We also define ${\bf{m}}$ through $\bm{m}=\chi h\bm{\hat{z}}+{\bf{m}}$,
highlighting the field-induced uniform magnetization $\chi h\bm{\hat{z}}$
($\chi$ is the magnetic susceptibility along the $z$-axis).
$\bm{n}_{\scriptscriptstyle\perp},{\bf{m}}_{\scriptscriptstyle\perp}$ form two
canonically conjugate pairs of Gaussian fields and may be considered “spin-
wave variables”. We also define ${\upomega}_{\alpha,\bm{k}}$ to be the
dispersion for the magnon branch $\alpha\in\\{y,z\\}$. In the case of a
Heisenberg model with exchange $J$, plus weak anisotropies of the order of
$\Delta_{\alpha}$, the dispersion relation takes the form
$\upomega_{\alpha,\bm{k}}=v_{m}\sqrt{k_{x}^{2}+k_{y}^{2}+\delta_{\alpha}^{2}}$,
where $v_{m}$ is the magnon velocity and $\delta_{\alpha}$ are the magnon gaps
in units of inverse length: $\delta_{y}=\Delta_{y}/v_{m}$,
$\delta_{z}=\sqrt{\Delta_{z}^{2}+h^{2}}/v_{m}$.

From this, we find that correlators of the form $\langle nn\rangle,\langle{\rm
mm}\rangle$ are even functions of frequency and given in the SM, and only the
“mixed” correlators $\langle n{\rm m}\rangle$ are odd in frequency:

$\displaystyle\langle n_{\alpha,\bm{k}}{\rm
m}_{\bar{\alpha},-\bm{k}}\rangle_{\omega_{n}}=\epsilon_{x\alpha\bar{\alpha}}\,\frac{\omega_{n}a_{0}^{2}}{S}\mathcal{D}_{\alpha}(\bm{k},\omega_{n}),$
(4)

where $\epsilon$ is the Levi-Civita tensor, $\overline{\alpha}\neq\alpha,x$,
$a_{0}$ is the lattice constant of the magnetic layer, and
$\mathcal{D}_{\alpha}^{-1}(\bm{k},\omega_{n})=\omega_{n}^{2}+\upomega_{\alpha,\bm{k}}^{2}$.

The magnetic operators $\mathcal{O}_{\Gamma,a}$ may be expressed as
polynomials in the staggered and uniform magnetizations $\bm{n},\bm{m}$. In
order to understand their form, we must proceed to a thorough symmetry
analysis.

|  zero field | ${\bm{h}}=h\bm{\hat{z}}$
---|---|---
|  lattice and spin |  lattice effective
paramagnet | $\mathsf{G}=P4/mmm1^{\prime}$ | $\mathsf{G}(\bm{0},h\bm{\hat{z}})=P4/mm^{\prime}m^{\prime}$ | $\mathsf{G}^{\rm eff}(\bm{0},h\bm{\hat{z}})=4/mm^{\prime}m^{\prime}$
high sym. AFM | $\mathsf{G}(\bm{\hat{x}},\bm{0})=\langle i,\mathcal{T}X,\mathcal{T}Y,C_{2x},\mathcal{T}C_{2z}\rangle$ | $\mathsf{G}(\bm{\hat{x}},h\bm{\hat{z}})=\langle i,XY,\mathcal{T}C_{2y},XC_{2z}\rangle$ | $\mathsf{G}^{\rm eff}(\bm{\hat{x}},h\bm{\hat{z}})=\langle i,\mathcal{T}C_{2y},C_{2z}\rangle$
low sym. AFM | $\mathsf{G}(\bm{\hat{e}},\bm{0})=\langle i,\mathcal{T}X,\mathcal{T}Y,\mathcal{T}C_{2z}\rangle$ | $\mathsf{G}(\bm{\hat{e}},h\bm{\hat{z}})=\langle i,XY,XC_{2z}\rangle$ | $\mathsf{G}^{\rm eff}(\bm{\hat{e}},h\bm{\hat{z}})=\langle i,C_{2z}\rangle$
Table 1:  Summary of the magnetic groups (in the Hermann-Mauguin notation) for
crystallographic layer group $G=P4/mmm$. Rows denote different values for
$\bm{n}$: without magnetic order, and for the two distinct orientations of the
staggered magnetization in the $xy$ plane (here $\bm{\hat{e}}$ is a generic
vector not along a high-symmetry axis in the plane), and columns specify the
zero and finite field cases, and the effective magnetic group for the
effective lattice theory. Here $i$ is inversion,
$P4/mmm1^{\prime}=P4/mmm\times({\rm Id}+\mathcal{T})$ is a grey group,
$P4/mm^{\prime}m^{\prime}=\left({\rm C}_{4h}+\mathcal{T}\times({\rm
D}_{4h}-{\rm C}_{4h})\right)\ltimes P$ is a black-white group in which
$m^{\prime}$ denotes the vertical mirror composed with $\mathcal{T}$. The
symbol $\langle\cdot\rangle$ indicates the group generated by the “$\cdot$”
operations.

Symmetry group considerations.— We call $G$ the lattice space group of the
crystal in the paramagnetic phase. We also define $\mathsf{G}$ to be its
associated (paramagnetic) magnetic space group, $\mathsf{G}=G\times\\{{\rm
Id},\mathcal{T}\\}$. Spontaneous ordering with staggered magnetization
$\bm{n}$ and/or an external field $\bm{h}$ reduce the system’s symmetry to
$\mathsf{G}(\bm{n},\bm{h})$. For example, the subgroup of $\mathsf{G}$ which
preserves the AFM order defines the magnetic space group of the AFM,
$\mathsf{G}(\bm{n},\bm{0})$, and the subgroup which preserves the Zeeman
interaction with the magnetic field $\bm{h}$ defines the magnetic space group
of the paramagnet in a field, $\mathsf{G}(\bm{0},\bm{h})$. In simple cases,
the magnetic space group of the AFM in an external field,
$\mathsf{G}(\bm{n},\bm{h})$, can be found as the intersection set of
$\mathsf{G}(\bm{n},\bm{0})$ and $\mathsf{G}(\bm{0},\bm{h})$.

Ultimately the PHV is defined in the effective action, $\mathcal{S}_{\rm PHV}$
from Eq. (1), purely in the phonon space. It results, in the framework
described here, from integrating out the magnetic degrees of freedom
$\bm{m},\bm{n}$, i.e. from carrying out the integral in Eq. (A2) of the SM.
One may then observe a simplifying feature of the representation of the
symmetry operations on the lattice terms. In particular, the translation
operation acts as the identity on the lattice strain field. Indeed, in
$\mathcal{L}_{sl}$, $\mathcal{O}_{\Gamma}$ and $\mathcal{E}_{\Gamma}$ are
independently translationally invariant. As a consequence, the appropriate
representation of $\mathsf{G}(\bm{n},\bm{h})$ acting on the terms in the
effective phonon action (in the presence of AFM order and in a field) is that
of the magnetic point group $\mathsf{G}^{\rm eff}(\bm{n},\bm{h})$ obtained
through the group morphism $\Pi:(W,t)\mapsto(W,0)$, where
$(W,t)\in\mathsf{G}(\bm{n},\bm{h})$ and $W,t$, are respectively the
(anti)linear and translational parts of $(W,t)$ and $0$ denotes here the zero
translation. The morphism theorem applies and

$\mathsf{G}^{\rm
eff}(\bm{n},\bm{h})\cong\mathsf{G}(\bm{n},\bm{h})/\mathsf{T}(\bm{n},\bm{h}),$
(5)

i.e. $\mathsf{G}^{\rm eff}(\bm{n},\bm{h})$ is isomorphic to the factor group
$\mathsf{G}(\bm{n},\bm{h})/\mathsf{T}(\bm{n},\bm{h})$, where
$\mathsf{T}(\bm{n},\bm{h})$ is the pure translation subgroup of
$\mathsf{G}(\bm{n},\bm{h})$. The group $\mathsf{G}^{\rm eff}(\bm{n},\bm{h})$
may now be used to analyze the PHV. In particular, distinct irreps
$\Gamma,\Gamma^{\prime}$ under ${\sf G}$ may collapse to the same irrep under
$\mathsf{G}^{\rm eff}(\bm{n},\bm{h})$, allowing the corresponding strain
components $\mathcal{E}_{\Gamma},\mathcal{E}_{\Gamma^{\prime}}$ to couple in
the PHV term.

Application.— As an example, we now consider a 3D crystal composed of finite
thickness regions whose symmetry is given by the layer group $G=P4/mmm$ (layer
group number 61, see e.g. Ref. Aroyo  _et al._ (2006)), whose point group is
${\rm D}_{4h}$ in the sense that $G$ is the semi-direct product of the point
group $4/mmm$ (${\rm D}_{4h}$) and the 2D translation group $P$ of a square
Bravais lattice. The layer group is the natural framework to describe those
properties of quasi-2d solids which are independent of the stacking structure
(however, some stackings may reduce the point group symmetries, and thereby
relax restrictions on the PHV). The cuprate Mott insulator Sr2CuO2Cl2 with
space group $I4/mmm$ (number 139) is such an example. The translation group
$P$ is generated by translations $X,\,Y$ within the $xy$ plane, i.e. by the
Bravais lattice vectors $\bm{\hat{x}},\bm{\hat{y}}$, respectively. To
characterize the lattice strain field, only the _point_ group is relevant, and
the decomposition of the strain tensor $\mathcal{E}$ into irreducible
representations (irreps) of ${\rm D}_{4h}$ reads:

$\displaystyle\mathcal{E}_{A_{1}}$
$\displaystyle=\mathcal{E}_{xx}+\mathcal{E}_{yy},\,\mathcal{E}_{B_{1}}=\mathcal{E}_{xx}-\mathcal{E}_{yy},\,\mathcal{E}_{B_{2}}=\mathcal{E}_{xy},$
$\displaystyle\mathcal{E}_{E}$
$\displaystyle=\\{\mathcal{E}_{Ex},\mathcal{E}_{Ey}\\}=\\{\mathcal{E}_{xz},\mathcal{E}_{yz}\\}.$
(6)

In the quasi-2D limit, we ignore the inter-layer spin-lattice coupling, so
$\mathcal{E}_{zz}$ in the $A_{1}$ irrep does not enter the Hall viscosity.

Consider next the determination of $\mathsf{G}(\bm{n},\bm{h})$ and
$\mathsf{G}^{\rm eff}(\bm{n},\bm{h})$ for this system. For concreteness, we
choose the field along the $z$-axis, ${\bm{h}}=h\bm{\hat{z}}$, and place
$\bm{n}$ in the $xy$ plane (which is energetically preferred for this field
orientation). For our purpose, we distinguish two cases: the “high symmetry”
AFMs for which $\bm{n}$ is aligned with a high-symmetry direction, i.e. to the
$x$ or $y$ axes or at 45 degrees between them, and the “low symmetry” AFMs
where $\bm{n}$ takes any other in-plane orientation. The paramagnetic and
magnetic groups for zero and nonzero fields with the above provisos are given
in Table 1.

For this field orientation, we are interested in $\kappa_{xy}$ and the
acoustic Faraday effect for sound waves propagating along $\bm{\hat{z}}$. Both
of these are odd under time-reversal and under vertical mirrors, and hence
even under their combination. This is compatible with all three cases of
$\mathsf{G}^{\rm eff}(\bm{n},\bm{h})$ in a non-zero applied field (see the
final column of Table 1), and hence both effects should be non-zero in these
situations. To determine them, we will need the PHV terms with the same
symmetries, i.e. those which are odd under break time-reversal $\mathcal{T}$
and vertical mirrors, and which are invariant under $\mathsf{G}^{\rm
eff}(\bm{n},\bm{h})$: the only compatible PHV coefficients are
$\eta^{H}_{B_{1},B_{2}}$ and $\eta^{H}_{E_{x},E_{y}}$. Other field
configurations can be analyzed similarly – see SM.

Contributing magnetic operators.—We are now in a position to ask about the
allowed forms of the $\mathcal{O}$ operators which appear in
$\mathcal{L}_{sl}$ and which contribute to a nonzero $\eta^{H}$, as well as
their correlations. The magnetic operators $\mathcal{O}$ are polynomials of
the Gaussian fields
$\bm{n}_{\scriptscriptstyle\perp},{\bf{m}}_{\scriptscriptstyle\perp}$, so that
we may calculate correlators of the $\mathcal{O}$’s using Wick’s theorem. We
keep only those which effectively yield “mixed” correlators of ${\rm m}$ and
$n$.

In the absence of a magnetic field, the PHV vanishes. Indeed, because $\bm{m}$
is odd under the $\mathcal{T}X$ and $\mathcal{T}Y$ symmetries of the AFM
state, all allowed polynomials are even in $\bm{m}$. In an external magnetic
field however, the uniform magnetization ${\bm{m}}$ acquires a static
component ${\bm{m}}_{0}=\chi h\bm{\hat{z}}$ so that ${\rm m}_{x}=-\chi hn_{z}$
to linear order in $h$. One can thereby obtain an odd-in-frequency
contribution to the two-point $\mathcal{O}$ correlators: this occurs in the
cases for which one magnetic operator $\mathcal{O}$ is even and the other is
odd in ${\bf m}_{\scriptscriptstyle\perp}$. To lowest order, the important
terms are $\mathcal{O}_{B_{1}}\sim h{\rm m}_{z}$, $\mathcal{O}_{B_{2}}\sim
n_{y}$, $\mathcal{O}_{E_{x}}\sim n_{z}$, $\mathcal{O}_{E_{y}}\sim h{\rm
m}_{y}$ (see SM Sec.B3) —contributions from higher powers of magnon operators
are parametrically small by a factor of $k_{B}T/J$ or $v_{m}\delta_{\alpha}/J$
(see SM Sec. D). Appropriately combined, the latter give

$\displaystyle\eta_{B_{1},B_{2}}^{H}\sim-h\langle{\rm
m}_{z}n_{y}\rangle,\quad\eta_{E_{x},E_{y}}^{H}$ $\displaystyle\sim h\langle
n_{z}{\rm m}_{y}\rangle.$ (7)

A careful analysis that restores the units (see SM Sec. A1) gives

$\displaystyle\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})=\frac{\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}S^{3}\chi
h}{d_{z}a_{0}^{2}}\mathcal{D}_{\alpha(\Gamma\Gamma^{\prime})}({\bm{q}}_{\perp},\omega=0)=\frac{\gamma_{\Gamma\Gamma^{\prime}}}{q_{x}^{2}+q_{y}^{2}+\delta_{\alpha}^{2}},$
(8)

where
$\gamma_{\Gamma\Gamma^{\prime}}=\frac{hS^{2}\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}{v^{3}_{m}gd_{z}a_{0}}$,
$\alpha(B_{1},B_{2})=y$, $\alpha(E_{x},E_{y})=z$, and ${\bm{q}}_{\perp}$
projects the momentum ${\bm{q}}$ to the xy plane. The factor of $d_{z}$, the
inter-layer spacing in the $z=c$ direction, is necessary to convert to bulk
three-dimensional elasticity.

_Acoustic Faraday effect.—_ This effect can be understood as arising from
splitting of the degeneracy between right- and left-circularly polarized sound
waves. This allows a transverse linearly-polarized wave propagating along
$\bm{\hat{z}}$ at a frequency $\omega_{\rm ph}$ to undergo a Faraday rotation.
The Faraday rotation angle $\Phi$ per unit length $L$ is simply given by the
difference in wavenumber between right- and left-circularly polarized waves at
a given frequency $\omega_{\rm ph}$. By calculating the dispersion relation
including the PHV, we find $\Phi/L$ is related to the Hall viscosity
coefficient $\eta^{H}_{E_{x}E_{y}}$ through (see the SM)

$\displaystyle\frac{\Phi}{L}=\frac{\eta^{H}_{E_{x}E_{y}}\omega_{\rm
ph}^{2}}{v_{T}^{3}\rho}=\frac{\gamma_{E_{x}E_{y}}\omega_{\rm
ph}^{2}}{v_{T}^{3}\rho\delta_{z}^{2}}+{O}(\omega_{\rm ph}^{3}).$ (9)

where $\rho$ is the mass density of the lattice, $v_{T}$ is the asymptotic
long wavelength transverse sound wave velocity for this propagation direction,
and recall $\delta_{z}$ is the out-of-plane magnon gap in units of inverse
length.

Figure 1: Intrinsic phonon thermal Hall conductivity, $\kappa_{xy}/T$ v.s. $T$
(rescaled in units of $v_{\rm ph}\delta_{\alpha}$). Only the value of
$F_{\alpha}(\tilde{T})$ in Eq. (12) is plotted, and the overall coefficient
not shown in the plot is set by the RHS of Eq. (13).

_Thermal Hall conductivity.—_ Following Ref. Qin  _et al._ (2012), the
intrinsic (by which we mean independent of impurities) non-collisional thermal
Hall conductivity $\kappa_{jk}$ is determined by the Berry curvature
$\varOmega^{i}_{{\bm{q}},\sigma}$ and the phonon dispersions
$\upomega_{{\bm{q}},\sigma}$:

$\displaystyle\kappa_{jk}=-\frac{1}{T}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}E\,E^{2}\sigma_{jk}(E)\frac{\mathop{}\\!\mathrm{d}n^{\text{eq}}_{B}(E)}{\mathop{}\\!\mathrm{d}E},$
(10)

where we defined the “conductivity” $\sigma_{jk}$:

$\displaystyle\sigma_{jk}(E)=-\sum_{\sigma}\int\frac{\mathop{}\\!\mathrm{d}^{3}q}{(2\pi)^{3}}\epsilon_{ijk}\varOmega^{i}_{{\bm{q}},\sigma}\Theta(E-\upomega_{{\bm{q}},\sigma}),$
(11)

where $\epsilon$ is the Levi-Civita tensor, the subscript $\sigma$ labels the
phonon branch, and $n_{B}^{\rm eq}$ is the Bose-Einstein distribution.

A non-zero Berry curvature $\bm{\varOmega}$ is obtained by finding the phonon
eigenfunctions including the PHV, which induces a vector potential in the
phonon Hamiltonian (see Sec. C2 of the SM). The leading nonzero terms are
linear in $\eta^{H}$, and to this order contributions to $\kappa_{H}$ from
$\eta^{H}_{\Gamma\Gamma^{\prime}}$ add–each is labeled
$\kappa_{H}^{\Gamma\Gamma^{\prime}}$ below. Defining a dimensionless
temperature $\tilde{T}=T/(v_{\rm ph}\delta_{\alpha})$, where $v_{\rm ph}$ is
the mean sound velocity, we find

$\displaystyle\frac{\kappa^{{\Gamma\Gamma^{\prime}}}_{xy}}{T}\sim\frac{\gamma_{\Gamma\Gamma^{\prime}}}{\rho
v_{\rm ph}}F_{\alpha(\Gamma\Gamma^{\prime})}(\tilde{T}),$ (12)

where $F_{\alpha(\Gamma\Gamma^{\prime})}$ is a scaling function. The detailed
forms of $F_{\alpha(\Gamma\Gamma^{\prime})}(\tilde{T})$ for different
$\Gamma\Gamma^{\prime}$ are given in the SM, and are numerically evaluated and
plotted in Fig. 1. Importantly, we find that the temperature scaling of
$\kappa_{H}/T$ is sensitive to the ${\bm{q}}$ dependence of $\eta^{H}$. If
$\eta^{H}$ is a constant, straightforward power counting shows that
$\kappa_{H}/T\sim T^{d-1}$, where $d$ is the spatial dimension. For $\eta^{H}$
with the ${\bm{q}}$ dependence of Eq. (8), the spin gap $\sim\delta_{\alpha}$
introduces an additional scale that controls the temperature dependence, and
$\kappa_{H}/T$ is non-monotonic. $\kappa_{H}/T$ increases from $T=0$ as
$T^{2}$, reaches a maximum at $T_{\rm max}\sim v_{\rm ph}\delta_{\alpha}$. At
$T/v_{\rm ph}\delta_{\alpha}\gg 1$, $\kappa_{H}/T$ decreases as $1/T^{\zeta}$,
where $\zeta\sim 1$ (see Fig. 1).

_Summary and discussion.—_ In this work, we presented a symmetry analysis of
the phonon Hall viscosity and resulting phonon Berry curvature in a
magnetoelastically coupled system. The procedure includes accounting for
separate scales for anisotropy-induced magnon gaps and the applied magnetic
field, and follows from a magnetic space group symmetry analysis. The symmetry
predictions were checked for a low-energy magnon model by an explicit
calculation in terms of spin correlation functions. We found that the phonon
Hall viscosity $\eta^{H}$ retains non-trivial dependence on the the scaling
variable ${\bm{q}}\delta$, even when the phonon momentum is small
${\bm{q}}a_{0}\ll 1$, due to the small spin gap $\sim\delta$.

As an example, we modeled the cuprate Mott insulator Sr2CuO2Cl2 with
tetragonal symmetry. We showed that PHV induces thermal Hall conductivity via
both Hall viscosities $\eta^{H}_{B_{1}B_{2}}$ and $\eta^{H}_{E_{x}E_{y}}$,
while the acoustic Faraday effect arises only from $\eta^{H}_{E_{x}E_{y}}$.

The treatment above implicitly uses spin-orbit coupling (SOC) throughout, in
the symmetry analysis and through the forms of the spin-lattice couplings.
Without SOC, one instead requires that the microscopic Hamiltonian (before any
spontaneous symmetry breaking) in the absence of an applied magnetic field has
a global SO(3) spin-rotation symmetry $SO(3)_{s}$. The latter acts only on
spin indices, and is distinct from and independent of the space group and
time-reversal symmetries, which act on coordinates, spatial derivatives, and
strain indices. Under the assumption that the applied magnetic field couples
only via the Zeeman interaction, then the symmetry constraints on the PHV are
significantly more stringent in the absence of SOC. We forgo a general
discussion here, but give a simple argument that whenever the applied magnetic
field and all local ordered spin moments lie in a single plane, the PHV
vanishes. Indeed, under those assumptions, all Zeeman and exchange fields are
invariant under the operation $\mathcal{T}\left[C_{2\bm{\hat{w}}}\right]_{s}$,
where the second operation is a $C_{2}$ rotation in spin space about an axis
$\bm{\hat{w}}$ normal to the plane containing the spins and field. Because
spin rotations do not act on the strain, this operation is indistinguishable
from time-reversal symmetry in the lattice effective action, and thus PHV is
prohibited. In a real material with weak but nonzero SOC, there will be an
additional smallness of the PHV due to weak SOC.

We now turn to an estimate of the magntitude of the phonon Hall effect induced
by the mechanism in this paper. The characteristic energy $v_{\rm
ph}\delta_{\alpha}$ and the maximum of $\kappa_{xy}/T$ can be related to
parameters that may be obtained from experiments/ab-initio calculations. The
thermal conductivity is most conveniently expressed as $\kappa_{xy}/T$ per
layer and in terms of the ratio of the phonon to the magnon velocity
$\Upsilon=v_{\rm ph}/v_{m}$. In units of thermal conductance, we have:

$\displaystyle\frac{d_{z}\kappa^{\Gamma\Gamma^{\prime}}_{xy}}{T}\sim\frac{\gamma_{\Gamma\Gamma^{\prime}}d_{z}}{\rho
v_{\rm
ph}}\frac{k_{B}^{2}}{\hbar}=\Upsilon\frac{\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}{J\,v_{m}/d_{z}}\frac{h}{\rho
a_{0}^{2}d_{z}v_{\rm ph}^{2}}\frac{k_{B}^{2}}{\hbar}\times{O}(1).$ (13)

Here $d_{z}$ is the inter-layer distance, and the numerical coefficient
(${O}(1)$) depends on the specific microscopic model. Because the
magnetoelastic couplings may be considered spatial derivatives of the
(anisotropic) exchange, and the magnon velocity is also set by exchange, we
expect that
$\frac{\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}{J\,v_{m}/d_{z}}\lesssim 1$
(smallness due to weak SOC also enters here). It could be more precisely
evaluated through ab initio calculations. The remaining factors behave as
$\frac{h}{\rho a_{0}^{2}d_{z}v_{\rm ph}^{2}}\sim h(C_{\rm
3D}a_{0}^{2}d_{z})^{-1}$, where $C_{\rm 3D}$ is the three-dimensional bulk
modulus. The latter is generally an eV energy so that for a magnetic field
$B=15$ T, i.e. $h=g_{s}\mu_{B}B\sim 30$ K, $h(C_{\rm
3D}a_{0}^{2}d_{z})^{-1}\sim 10^{-4}$. For cuprates, $\Upsilon<1$ due to the
large exchange energy. Thus, the thermal Hall conductivity due to the PHV
mechanism is smaller by a factor of at least $10^{-4}$ than measured values
for cuprates.

This leaves open other mechanisms for the phonon Hall effect in the cuprates.
It is possible that the effects of the PHV are enhanced by impurity scattering
of phonons, not included here (see a recent preprint Guo and Sachdev (2021)
for discussion on related issues). There might also be alternative mechanisms
to generate larger PHV, e.g. from charged impurities. Finally, phonons may
aquire chirality via skew scattering of phonons from spins. This mechanism,
which is independent of the PHV, requires a true non-equilibrium transport
treatment.

A partial means to distinguish between these possibilities is to compare
directly the acoustic Faraday and phonon Hall effect measurements to see if
they can be consistently related to comparable PHVs. Beyond the cuprates, the
above discussion can be easily extended to any antiferromagnet, and may be
used to guide a search for large PHVs. At a general level, it is clear that
large spin-orbit coupling or non-coplanar magnetic order (which evades the
weak SOC smallness), and large magnetoelastic couplings are beneficial for
enhancing the PHV. We expect our results will be useful to guide future
experiments and computations in this active area.

###### Acknowledgements.

We acknowledge Léo Mangeolle for a collaboration on a related project and
helpful discussions. M.Y. also benefited from collaborations with Natalia
Perkins and Rafael Fernandes on a past work on phonon Hall viscosity. L.B. was
supported by the DOE, Office of Science, Basic Energy Sciences under Award No.
DE-FG02-08ER46524. M.Y. is supported in part by the Gordon and Betty Moore
Foundation through Grant GBMF8690 to UCSB and by the National Science
Foundation under Grant No. NSF PHY-1748958. L.S. acknowledges funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (Grant agreement No. 853116, “TRANSPORT”) as
well as from the ToRe Idex-Lyon breakthrough program.

## References

  * Barkeshli  _et al._ (2012) Maissam Barkeshli, Suk Bum Chung, and Xiao-Liang Qi, “Dissipationless phonon hall viscosity,” Phys. Rev. B 85, 245107 (2012).
  * Tuegel and Hughes (2017) Thomas I. Tuegel and Taylor L. Hughes, “Hall viscosity and the acoustic faraday effect,” Phys. Rev. B 96, 174524 (2017).
  * Shapourian  _et al._ (2015) Hassan Shapourian, Taylor L. Hughes, and Shinsei Ryu, “Viscoelastic response of topological tight-binding models in two and three dimensions,” Phys. Rev. B 92, 165131 (2015).
  * Avron  _et al._ (1995) J. E. Avron, R. Seiler, and P. G. Zograf, “Viscosity of quantum hall fluids,” Phys. Rev. Lett. 75, 697–700 (1995).
  * Avron (1998) J. E. Avron, “Odd viscosity,” Journal of Statistical Physics 92, 543–557 (1998).
  * Read (2009) N. Read, “Non-abelian adiabatic statistics and hall viscosity in quantum hall states and ${p}_{x}+i{p}_{y}$ paired superfluids,” Phys. Rev. B 79, 045308 (2009).
  * Read and Rezayi (2011) N. Read and E. H. Rezayi, “Hall viscosity, orbital spin, and geometry: Paired superfluids and quantum hall systems,” Phys. Rev. B 84, 085316 (2011).
  * Bradlyn  _et al._ (2012) Barry Bradlyn, Moshe Goldstein, and N. Read, “Kubo formulas for viscosity: Hall viscosity, ward identities, and the relation with conductivity,” Phys. Rev. B 86, 245309 (2012).
  * Sytcheva  _et al._ (2010) A. Sytcheva, U. Löw, S. Yasin, J. Wosnitza, S. Zherlitsyn, P. Thalmeier, T. Goto, P. Wyder, and B. Lüthi, “Acoustic faraday effect in ${\text{tb}}_{3}{\text{ga}}_{5}{\text{o}}_{12}$,” Phys. Rev. B 81, 214415 (2010).
  * Qin  _et al._ (2012) Tao Qin, Jianhui Zhou, and Junren Shi, “Berry curvature and the phonon hall effect,” Phys. Rev. B 86, 104305 (2012).
  * Vinkler-Aviv and Rosch (2018) Yuval Vinkler-Aviv and Achim Rosch, “Approximately quantized thermal hall effect of chiral liquids coupled to phonons,” Phys. Rev. X 8, 031032 (2018).
  * Zhang  _et al._ (2019) Xiaoou Zhang, Yinhan Zhang, Satoshi Okamoto, and Di Xiao, “Thermal hall effect induced by magnon-phonon interactions,” Phys. Rev. Lett. 123, 167202 (2019).
  * Chen  _et al._ (2020) Jing-Yuan Chen, Steven A. Kivelson, and Xiao-Qi Sun, “Enhanced thermal hall effect in nearly ferroelectric insulators,” Phys. Rev. Lett. 124, 167601 (2020).
  * Ye  _et al._ (2020) Mengxing Ye, Rafael M. Fernandes, and Natalia B. Perkins, “Phonon dynamics in the kitaev spin liquid,” Phys. Rev. Research 2, 033180 (2020).
  * Saito  _et al._ (2019) Takuma Saito, Kou Misaki, Hiroaki Ishizuka, and Naoto Nagaosa, “Berry phase of phonons and thermal hall effect in nonmagnetic insulators,” Phys. Rev. Lett. 123, 255901 (2019).
  * Boiteux  _et al._ (1971) M. Boiteux, P. Doussineau, B. Ferry, J. Joffrin, and A. Levelut, “Acoustical faraday effect in antiferromagnetic ${\mathrm{cr}}_{2}$${\mathrm{o}}_{3}$,” Phys. Rev. B 4, 3077–3088 (1971).
  * Grissonnanche  _et al._ (2019) Gaël Grissonnanche, Anaëlle Legros, Sven Badoux, Etienne Lefrançois, Victor Zatko, Maude Lizaire, Francis Laliberté, Adrien Gourgout, Jianshi Zhou, Sunseng Pyon, Tomohiro Takayama, Hidenori Takagi, Shimpei Ono, Nicolas Doiron-Leyraud, and Louis Taillefer, “Giant thermal hall conductivity in the pseudogap phase of cuprate superconductors,” Nature 571, 376–380 (2019).
  * Boulanger  _et al._ (2020) Marie-Eve Boulanger, Gaël Grissonnanche, Sven Badoux, Andréanne Allaire, Étienne Lefrançois, Anaëlle Legros, Adrien Gourgout, Maxime Dion, C. H. Wang, X. H. Chen, R. Liang, W. N. Hardy, D. A. Bonn, and Louis Taillefer, “Thermal hall conductivity in the cuprate mott insulators nd2cuo4 and sr2cuo2cl2,” Nature Communications 11, 5325 (2020).
  * Grissonnanche  _et al._ (2020) G. Grissonnanche, S. Thériault, A. Gourgout, M. E. Boulanger, E. Lefrançois, A. Ataei, F. Laliberté, M. Dion, J. S. Zhou, S. Pyon, T. Takayama, H. Takagi, N. Doiron-Leyraud, and L. Taillefer, “Chiral phonons in the pseudogap phase of cuprates,” Nature Physics (2020), 10.1038/s41567-020-0965-y.
  * Note (1) Symmetry allows additional PHV terms which are not expressed solely in terms of the strain but also in terms of the time derivative of the rotation $\mathcal{M}_{ij}=\frac{1}{2}(\partial_{i}u_{j}-\partial_{j}u_{i})$. Such terms, however, cannot be induced by conventional magnetoelastic coupling which involves strain only.
  * Aroyo  _et al._ (2006) Mois I. Aroyo, Asen Kirov, Cesar Capillas, J. M. Perez-Mato, and Hans Wondratschek, “Bilbao crystallographic server. ii. representations of crystallographic point groups and space groups,” Acta Crystallographica Section A Foundations of Crystallography 62, 115–128 (2006).
  * Guo and Sachdev (2021) Haoyu Guo and Subir Sachdev, “Extrinsic phonon thermal Hall transport from Hall viscosity,” arXiv e-prints , arXiv:2103.02614 (2021), arXiv:2103.02614 [cond-mat.str-el] .
  * Sachdev (2009) Subir Sachdev, “Quantum phase transitions,”  (2009), 10.1017/cbo9780511973765.
  * Koster (1963) George F. Koster,  _Properties of the thirty-two point groups_ (M.I.T. Press, Cambridge, Mass., 1963) p. 104 p.
  * Yildirim  _et al._ (1995) T. Yildirim, A. B. Harris, Amnon Aharony, and O. Entin-Wohlman, “Anisotropic spin hamiltonians due to spin-orbit and coulomb exchange interactions,” Phys. Rev. B 52, 10239–10267 (1995).
  * Benfatto  _et al._ (2006) L. Benfatto, M. B. Silva Neto, A. Gozar, B. S. Dennis, G. Blumberg, L. L. Miller, Seiki Komiya, and Yoichi Ando, “Field dependence of the magnetic spectrum in anisotropic and dzyaloshinskii-moriya antiferromagnets. ii. raman spectroscopy,” Phys. Rev. B 74, 024416 (2006).

## Supplemental Material

In the supplemental material (SM), we present additional details of the
modeling and computations.

Sec. A reviews the detailed derivation of the correlation function of
staggered and ferromagnetic field in Eq. (4) of the main text based on the
non-linear sigma model formulation. Summaries of the point group symmetry
operations and symmetry allowed magnetic operators in the spin-lattice
coupling [Eq. (2) in main the text] are given in Sec. B. In Sec. C, we show
details to analyze the acoustic Faraday and phonon thermal Hall effects.

For completeness, a few additional considerations are presented. We first show
that the next order in $1/S$ contribution to $\eta^{H}$ from two magnon
fluctuations are small in magnitude at low temperature, and can thus be
ignored (Sec. D). Next, we show the microscopic Hamiltonian allowed by crystal
symmetry in Sec. E. A comparison between the microscopic Hamiltonian and the
low energy effective action analysis is helpful to infer the spin gap and
reveal the microscopic origin of the relevant spin-lattice coupling for a
material. Finally, the phonon Hall viscosity with in-plane magnetic field is
discussed (Sec. F).

##  A Non-linear sigma model formulation of spin dynamics

In this section, we derive the low energy field theory and correlations
function for the antiferromagnetic state. Before doing so, we address briefly
the validity of the low energy approach. Since
$\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})$ arises from virtual rather than
on-shell magnetic excitations, there is no kinetic constraint, and spin
fluctuations at all energies may contribute to it. However, high energy
magnetic excitations are still polynomially suppressed by their energies, and
indeed we have checked that, within a full spin wave calculation,
contributions from high-energy magnons are suppressed by additional powers of
$k_{B}T/J$ with respect to the low energy ones. This justifies the low energy
approach.

We proceed with the standard non-linear sigma model (NLSM) formulation for
collinear antiferromagnets and obtain the two-point correlation functions for
a staggered field (${\bm{n}}$) and a ferromagnetic field ($\bm{m}$) in an
external magnetic field ($\bm{h}$) in a two-dimensional (2d) spin system. It
gives Eq. (4) in the main text.

Following Ref. Sachdev (2009), the coherent state path integral for a
Heisenberg spin with spin value $S$ at site ${\bf{r}}$ can be obtained in the
basis of the (unit) vector field ${\bm{e}}_{\bf{r}}$, which is defined through
$\bm{\hat{S}}_{\bf{r}}\ket{{\bm{e}}_{\bf{r}}}=S{\bm{e}}_{\bf{r}}\ket{{\bm{e}}_{\bf{r}}}$
and $|{\bm{e}}_{\bf{r}}|^{2}=1$. For a collinear antiferromagnet, the ground
state has spins oriented in opposite directions on the two sublattices
(defined as $A$, $B$ sublattice hereafter). The low energy spin dynamics can
be described by a set of continuous fields, which include the staggered
${\bm{n}}\sim{\bm{e}}_{A}-{\bm{e}}_{B}$ and uniform
$\bm{m}\sim{\bm{e}}_{A}+{\bm{e}}_{B}$ magnetization fields, where the
subscript $A,\,B$ label the A, B sublattice. To be accurate, we use

$\displaystyle{\bm{e}}_{\bf{r}}=(-1)^{\bf{r}}{\bm{n}}_{{\bf{r}}}\sqrt{1-(\bm{m}_{{\bf{r}}})^{2}}+\bm{m}_{{\bf{r}}},$
(14)

where $\bm{m}$ is the uniform magnetization per site in units of the
saturation magnetization ($=S$ semiclassically), and $(-1)^{\bf{r}}$ is a sign
equal to $+1$ on the A sublattice and $-1$ on the B sublattice. The effective
spin action (for the isotropic Heisenberg model) is

$\displaystyle\mathcal{Z}_{s}$
$\displaystyle=\int\mathop{}\\!\mathcal{D}{\bm{n}}\mathop{}\\!\mathcal{D}\bm{m}\delta({\bm{n}}^{2}-1)\delta({\bm{n}}\cdot\bm{m})\exp{(-\mathcal{S}_{s})}$
$\displaystyle\mathcal{S}_{s}$
$\displaystyle=\frac{1}{2}\int\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y\mathop{}\\!\mathrm{d}\tau\\{\frac{v_{m}}{ga_{0}}[({\bm{\nabla}}_{x}{\bm{n}})^{2}+({\bm{\nabla}}_{y}{\bm{n}})^{2}]+\frac{v_{m}gS^{2}}{a^{3}_{0}}\,\bm{m}^{2}-\frac{2iS}{a_{0}^{2}}\bm{m}\cdot\left({\bm{n}}\times\frac{\partial{\bm{n}}}{\partial\tau}-i{\bm{h}}\right)\\},$
(15)

In the second equation, we introduced the spin wave (magnon) velocity
$v_{m}\sim Ja_{0}S$, with $a_{0}$ the lattice constant, and the coefficient of
${\bm{m}}^{2}$ defines the coupling $g$ with $g\sim S^{-1}$. We ignore the
spin wave velocity along $\bm{\hat{z}}$ because of the much weaker interlayer
spin exchange. In the standard procedure the Gaussian field $\bm{m}$ is then
integrated out. As $\bm{m}$ does not have its own dynamics, the action can be
obtained by replacing $\bm{m}$ with its saddle-point solution,

$\displaystyle\bm{m}=\chi\left(i{\bm{n}}\times\frac{\partial{\bm{n}}}{\partial\tau}+{\bm{h}}-{\bm{n}}({\bm{n}}\cdot{\bm{h}})\right),$
(16)

which defines the susceptibility $\chi=a_{0}/v_{m}gS$.

Below, we consider an antiferromagnet with the Néel vector along the $x$-axis,
i.e. $\lim_{h\rightarrow 0}\langle\bm{n}\rangle=n_{0}\bm{\hat{x}}$, and the
magnetic field perpendicular to the Néel vector. As the spin-lattice coupling
is most transparently expressed in terms of ${\bm{n}}$, $\bm{m}$ fields to
reveal the symmetries, our goal below is to obtain the correlators of
${\bm{n}}$ and $\bm{m}$. The staggered field can be parameterized by
${\bm{n}}=\\{n_{0},n_{y},n_{z}\\}=(n_{0},{\bm{n}}_{\scriptscriptstyle\perp})$,
where $n_{0}=\sqrt{1-{\bm{n}}_{\scriptscriptstyle\perp}^{2}}$ is the order
parameter, and $n_{y,z}$ are transverse fluctuations (spin waves). Note from
Eq. (16) that $\bm{m}$ can be decomposed into components of zeroth order in
${\bm{h}}$ and first order in ${\bm{h}}$. The first-order term includes the
static part, i.e. the field induced uniform magnetization $\chi{\bm{h}}$ as
well as $-\chi{\bm{n}}({\bm{n}}\cdot{\bm{h}})$ to satisfy the constraint
$\bm{n}\cdot\bm{m}=0$ at first order in ${\bm{h}}$. We then define $\bf{m}$
through $\bm{m}=\chi{\bm{h}}+\bf{m}$, and within linear spin wave theory,
${\rm m}_{x}=-\chi n_{x}({\bm{n}}\cdot{\bm{h}})\rightarrow-\chi
n_{0}(n_{y}h_{y}+n_{z}h_{z})$. Plugging $\chi{\bm{h}}$ and ${\rm m}_{x}$ to
$\mathcal{S}_{s}$ into Eq. (15), we obtain the Lagrangian density in terms of
${\bm{n}}_{\scriptscriptstyle\perp}$ and $\bf{m}_{\scriptscriptstyle\perp}$.

#### With an external field ${\bm{h}}=h\bm{\hat{z}}$.

$\displaystyle\mathcal{S}_{s}=\frac{1}{2}\sum_{{\bm{k}},n}\begin{pmatrix}n_{y}&{\rm
m}_{y}&n_{z}&{\rm
m}_{z}\end{pmatrix}_{{\bm{k}},n}\begin{pmatrix}\frac{1}{gv_{m}a_{0}}(v_{m}^{2}\bm{k}^{2}+\Delta_{y}^{2})&&&-\frac{S\omega_{n}}{a_{0}^{2}}\\\
&\frac{v_{m}gS^{2}}{a_{0}^{3}}&-\frac{S\omega_{n}}{a_{0}^{2}}&\\\
&\frac{S\omega_{n}}{a_{0}^{2}}&\frac{1}{gv_{m}a_{0}}(v_{m}^{2}\bm{k}^{2}+h^{2}+\Delta_{z}^{2})&\\\
\frac{S\omega_{n}}{a_{0}^{2}}&&&\frac{v_{m}gS^{2}}{a_{0}^{3}}\\\
\end{pmatrix}\begin{pmatrix}n_{y}\\\ {\rm m}_{y}\\\ n_{z}\\\ {\rm m}_{z}\\\
\end{pmatrix}_{-{\bm{k}},-n}$ (17)

Here, the Fourier transform follows the convention
$f_{{\bm{k}}}=\frac{1}{L}\int\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y\,e^{-i{\bm{k}}\cdot{\bm{x}}}f({\bm{x}})$,
where $L$ is the linear size of the sample ($\int\\!dxdy\,1=L^{2}$). Inverting
the matrix, we obtain the correlator:

$\displaystyle\langle T_{\tau}\begin{pmatrix}n_{\alpha}\\\ {\rm
m}_{\bar{\alpha}}\end{pmatrix}_{\bm{k}}\begin{pmatrix}n_{\alpha}&{\rm
m}_{\bar{\alpha}}\end{pmatrix}_{-\bm{k}}\rangle_{\omega_{n}}=\mathcal{D}_{\alpha}(\bm{k},\omega_{n})\begin{pmatrix}v_{m}ga_{0}&\frac{\omega_{n}a_{0}^{2}}{S}\,\epsilon_{x\alpha\bar{\alpha}}\\\
-\frac{\omega_{n}a_{0}^{2}}{S}\,\epsilon_{x\alpha\bar{\alpha}}&\frac{\upomega^{2}_{\alpha,\bm{k}}a_{0}^{3}}{v_{m}gS^{2}}\end{pmatrix}.$
(18)

Here $a_{0}$ is the lattice constant,
$\mathcal{D}_{\alpha}^{-1}(\bm{k},\omega_{n})=\omega_{n}^{2}+\upomega_{\alpha,\bm{k}}^{2}$,
$\upomega_{\alpha,\bm{k}}$ is the dispersion for the magnon branch
$\alpha\in\\{y,z\\}$, such that
$\upomega_{\alpha,\bm{k}}=v_{m}\sqrt{k_{x}^{2}+k_{y}^{2}+\delta_{\alpha}^{2}}$,
the spin gap is determined by the XYZ anisotropy ($\sim\Delta_{\alpha}$) and
external field strength $h$ by $\delta_{y}=\Delta_{y}/v_{m}$,
$\delta_{z}=\sqrt{\Delta_{z}^{2}+h^{2}}/v_{m}$.

#### With an external field ${\bm{h}}=h\bm{\hat{y}}$.

Because here also ${\bm{h}}\perp\langle{\bm{n}}\rangle$, the situation is
similar to that when ${\bm{h}}=h\bm{\hat{z}}$. The effective spin action
simply changes according to $\Delta_{y}^{2}\rightarrow(\Delta_{y}^{2}+h^{2})$,
$\Delta_{z}^{2}+h^{2}\rightarrow\Delta_{z}^{2}$ in Eq. (17). Consequently,
$\delta_{y}=\sqrt{\Delta_{y}^{2}+h^{2}}/v_{m}$ and
$\delta_{z}=\Delta_{z}/v_{m}$.

###  A.1 Normalization of operators, Fourier conventions, etc.

Here we discuss the conventions used to obtain the Hall viscosity form given
in the main text, Eq. (8). We begin with a consideration of units. Hall
viscosity, Eq. (1), is defined as a coefficient in a three-dimensional elastic
theory. We employ Fourier conventions for the strain which are appropriate to
a three-dimensional system, so that

$\mathcal{E}_{\Gamma,\bm{q}}=\frac{1}{\sqrt{V}}\int\\!d^{3}\bm{x}\,\mathcal{E}_{\Gamma}(\bm{x})e^{-i\bm{q}\cdot\bm{x}},\qquad\mathcal{E}_{\Gamma}(\bm{x})=\frac{1}{\sqrt{V}}\sum_{\bm{q}}e^{i\bm{q}\cdot\bm{x}}\mathcal{E}_{\Gamma,\bm{q}},$
(19)

where $V$ is the volume of the system. With this convention, since the real
space strain $\mathcal{E}_{\Gamma}(\bm{x})$ is dimensionless,
$\mathcal{E}_{\Gamma,\bm{q}}$ has dimensions of $1/\sqrt{V}$. Then Eq. (1)
implies that, because the action is dimensionless, $\eta^{H}$ has units of
inverse volume.

Next consider the spin-lattice coupling in Eq. (2). The Lagrange density has
units of energy density, so that the combination
$\lambda_{\Gamma}\mathcal{O}_{\Gamma}$ must have units of energy density. We
assigned energy units to $\lambda$, which requires $\mathcal{O}_{\Gamma}$ to
scale as a number density. How this is precisely realized depends upon our
treatment of the third dimension. We will proceed here with the treatment as a
discrete layered system, so that the corresponding action is
$\mathcal{S}_{sl}=\sum_{z}\int\\!dxdy\,\mathcal{L}_{sl}$. Then
$\mathcal{O}_{\Gamma}$ has units of inverse length squared. Let us see how
these factors appear in the derivation of Eq. (8).

All the contributions to the spin-lattice coupling arise microscopically from
expressions (see Sec. E) of the form
$\sum_{\bf{r}}\lambda_{\Gamma}\mathcal{E}_{\Gamma}(SS)_{\Gamma,{\bf{r}}}$,
where $(SS)_{\Gamma,{\bf{r}}}$ represents a sum of spin bilinears in the
vicinity of site ${\bf{r}}$. We convert this to continuum fields using Eq.
(14) and $\sum_{\bf{r}}\rightarrow a_{0}^{-2}\sum_{z}\int\\!dxdy\,$, which
gives

$\mathcal{O}_{B_{1}}=\frac{S^{2}\chi h{\rm
m}_{z}}{a_{0}^{2}},\qquad\mathcal{O}_{B_{2}}=\frac{S^{2}n_{0}n_{y}}{a_{0}^{2}},\qquad\mathcal{O}_{E_{x}}=\frac{S^{2}n_{0}n_{z}}{a_{0}^{2}},\qquad\mathcal{O}_{E_{y}}=\frac{S^{2}\chi
h{\rm m}_{y}}{a_{0}^{2}}.$ (20)

The factors of $S^{2},\chi,1/a_{0}^{2},n_{0}$ are subsumed in the $\sim$ in
the discussion of the main text. Now consider evaluating Eq. (3). We must take
care due to the combination of the three-dimensional Fourier transform
convention for elasticity with our two-dimensional magnetism theory. Consider
the $B_{1}-B_{2}$ contribution. What actually arises in the effective action
is

$\mathcal{S}_{sl}^{B_{1}B_{2}}=\sum_{z,z^{\prime}}\int\\!dxdyd\tau\,\int\\!dx^{\prime}dy^{\prime}d\tau^{\prime}\,\lambda_{B_{1}}\lambda_{B_{2}}\mathcal{E}_{B_{1}}(x,y,z,\tau)\mathcal{E}_{B_{2}}(x^{\prime},y^{\prime},z^{\prime},\tau^{\prime})\langle\mathcal{O}_{B_{1}}(x,y,z,\tau)\mathcal{O}_{B_{2}}(x^{\prime},y^{\prime},z^{\prime},\tau^{\prime})\rangle.$
(21)

Since we assume no spin correlations between layers, the summand is non-zero
only for $z^{\prime}=z$. Inserting the three-dimensional Fourier expression of
Eq. (19) for the strains gives

$\mathcal{S}_{sl}^{B_{1}B_{2}}=\frac{1}{V}\sum_{z}\sum_{\bm{q},\bm{q}^{\prime}}\int\\!dxdyd\tau\,\int\\!dx^{\prime}dy^{\prime}d\tau^{\prime}\,\lambda_{B_{1}}\lambda_{B_{2}}\mathcal{E}_{B_{1},\bm{q}}(\tau)\mathcal{E}_{B_{2},\bm{q}^{\prime}}(\tau^{\prime})e^{i(\bm{q}_{\perp}\cdot\bm{x}_{\perp}+\bm{q}^{\prime}_{\perp}\cdot\bm{x}^{\prime}_{\perp})}e^{i(q_{z}+q^{\prime}_{z})z}\langle\mathcal{O}_{B_{1}}(x,y,\tau)\mathcal{O}_{B_{2}}(x^{\prime},y^{\prime},\tau^{\prime})\rangle.$
(22)

The sum over $z$ gives $N_{z}\delta_{q_{z}+q^{\prime}_{z},0}$, where
$N_{z}=L_{z}/d_{z}$ is the number of layers. Combining this with the $1/V$
prefactor gives $\frac{L_{z}}{d_{z}}\frac{1}{V}=\frac{1}{d_{z}L^{2}}$. We can
now use the $1/L^{2}$ to form the prefactors of the two-dimensional Fourier
transform for each of the two $\mathcal{O}$ operators. Hence

$\mathcal{S}_{sl}^{B_{1}B_{2}}=\sum_{\bm{q}_{\perp},\bm{q}^{\prime}_{\perp}}\sum_{q_{z}}\frac{\lambda_{B_{1}}\lambda_{B_{2}}}{d_{z}}\int\\!d\tau
d\tau^{\prime}\,\mathcal{E}_{B_{1},\bm{q}}(\tau)\mathcal{E}_{B_{2},\bm{q}^{\prime}}(\tau^{\prime})\langle\mathcal{O}_{B_{1},-\bm{q}_{\perp}}(\tau)\mathcal{O}_{B_{2},-\bm{q}^{\prime}_{\perp}}(\tau^{\prime})\rangle$
(23)

Due to momentum conservation, $\bm{q}^{\prime}_{\perp}=-\bm{q}_{\perp}$ is the
only non-zero correlator, and we obtain using Eq. (20) finally

$\eta^{H}_{B_{1},B_{2}}(\bm{q})=\frac{\lambda_{B_{1}}\lambda_{B_{2}}}{d_{z}}\frac{S^{2}\chi
h}{a_{0}^{2}}\frac{S^{2}n_{0}}{a_{0}^{2}}\times\left[-i\partial_{\omega}\langle{\rm
m}_{z,-\bm{q}_{\perp}}n_{y,\bm{q}_{\perp}}\rangle_{\omega_{n}\rightarrow-i\omega+0^{+}}\right].$
(24)

Inserting Eq. (4) for the ${\rm m}_{z}-n_{y}$ correlation function inside the
square bracket, one obtains the result in Eq. (8) for $\Gamma=B_{1}$,
$\Gamma^{\prime}=B_{2}$.

We took some pains to present this in great detail for clarity, but the result
can also be understood schematically on dimensional grounds: the three-
dimensional Fourier tranform differs from the two-dimensional one by a factor
of the inverse of the square root of a length in the $z$ direction. Converting
the 2d to 3d Fourier conventions for the two $\mathcal{O}$ operators appearing
at second order in the spin-lattice coupling, one obtains an overall factor of
$1/d_{z}$. This is the factor in the first term in Eq. (24). The remaining
factors were explained previously as arising from the conversion from the
lattice to the 2d continuum theory.

We note furthermore that it would have been possible to formulate the magnetic
correlations in three dimensions as well, which is in a sense more general,
and would also avoid some of this confusion. We opted for the present
formulation in order to emphasize, as discussed in the main text, that the
results apply to any three-dimensional structure composed of such 2d layers,
and that no three-dimensional magnetic correlations are required to induce the
desired PHV terms, even those which involve the inherently three-dimensional
$E_{x},E_{y}$ strains.

##  B Group theory analysis

Here, the relevant symmetry groups and symmetry operations on the strain and
spin fields are listed for reference. We consider the example discussed in the
main text, i.e. a crystal with layer group symmetry $P4/mmm$ (number 61 of
layer group). The symmetry allowed magnetic operators that appear in the spin-
lattice coupling [Eq. (2) in the main text] are also presented in Tabs. 3 and
4.

###  B.1 List of point group symmetry operations

The generators of the point group symmetry for the underlying crystal, ${\rm
D}_{4h}$, are ${\rm D}_{4h}=\langle C_{4z},\sigma_{h},\sigma_{v}\rangle$. They
give 16 point group elements: ${\rm D}_{4h}=\\{{\rm
Id},2C_{4z},C_{2z},2C_{2}^{\prime},2C_{2}^{\prime\prime},i,\sigma_{h},2\sigma_{v},2\sigma_{d},2S_{4}\\}$
Koster (1963). Here ${\rm Id}$ denotes the identity. $2C_{2}^{\prime}$ denotes
the two $\pi$ rotation transformations, around the $\bm{\hat{x}}$ or
$\bm{\hat{y}}$ axes, respectively. $2C_{2}^{\prime\prime}$ denotes the two
$\pi$ rotations around the diagonal axes
$\bm{\hat{a}}=\frac{1}{\sqrt{2}}(\bm{\hat{x}}+\bm{\hat{y}})$ and
$\bm{\hat{b}}=\frac{1}{\sqrt{2}}(-\bm{\hat{x}}+\bm{\hat{y}})$, respectively.
$\sigma_{v}$ and $\sigma_{d}$ are mirror transformations with a mirror plane
perpendicular to the $xy$ plane, whose normal direction is along
$\bm{\hat{x}},\bm{\hat{y}}$ and diagonal axis $\bm{\hat{a}},\bm{\hat{b}}$,
respectively. $\sigma_{h}$ is a horizontal mirror reflection. The convention
is shown in Fig. 2, such that the inversion transformation $i$ can be obtained
from $i=C_{2\alpha}\sigma_{\alpha}$, where $C_{2\alpha}=C_{2x,2y,2z,2a,2b}$,
$\sigma_{\alpha}=\sigma_{vx,vy,h,da,db}$. Here, $\alpha=x,y,z,a,b$ denotes the
axis of two-fold rotation for $C_{2\alpha}$, and the axis is also the normal
vector of mirror plane for the respective $\sigma_{\alpha}$, with
$\alpha=vx,vy,h,da,db$.

An external field ${\bm{h}}=h\bm{\hat{z}}$ breaks the point group ${\rm
D}_{4h}$ down to ${\rm C}_{4h}=\langle
C_{4z},\sigma_{h}\rangle=\\{E,2C_{4z},C_{2z},i,\sigma_{h},2S_{4}\\}$ through
the Zeeman term $H_{Z}={\bm{h}}\cdot\sum_{\bf{r}}{\bm{S}}_{\bf{r}}$. In terms
of the antiunitary symmetries, $\mathcal{T}$ is broken, while the Zeeman term
is invariant under the antiunitary symmetries
$2\mathcal{T}{C_{2}^{\prime}}\oplus 2\mathcal{T}{C_{2}^{\prime\prime}}\oplus
2\mathcal{T}{\sigma_{v}}\oplus 2\mathcal{T}{\sigma_{d}}$. This gives the
black-white magnetic point group $4/mm^{\prime}m^{\prime}=\left({\rm
C}_{4h}+\mathcal{T}\times({\rm D}_{4h}-{\rm C}_{4h})\right)$ as listed in Tab.
I.

An external field ${\bm{h}}=h\bm{\hat{y}}$ breaks ${\rm D}_{4h}$ down to ${\rm
C}^{\prime}_{2h}=\langle C_{2y},i\rangle=\\{E,C_{2y},i,\sigma_{vy}\\}$, whose
abstract group structure is the same as that of ${\rm C}_{2h}$. It also
preserves the antiunitary symmetries
$\mathcal{T}{C_{2z}}\oplus\mathcal{T}{\sigma_{h}}\oplus\mathcal{T}{\sigma_{vx}}\oplus\mathcal{T}{C_{2x}}$,
which can be generated by e.g. $\mathcal{T}{C_{2z}}$ and ${\rm
C}^{\prime}_{2h}$.

Figure 2: Point group symmetry transformations for ${\rm D}_{4h}$. The label
$\otimes$ denotes horizontal mirror symmetry $\sigma_{h}$ that relates
$\times$ and $\circ$.

###  B.2 Symmetry operations

The point group symmetry operations acting on the strain field
$\mathcal{E}_{\Gamma}$, the spin vector ${\bm{S}}$, the lattice coordinate
$\\{x,y,z\\}$, the continuous spin fields ${\bm{n}},\bm{m}$ are listed in
Table. 2. For brevity, only the point group symmetry generators listed in Tab.
1 and Tab. 5 are shown.

| $C_{4z}$ | $\sigma_{vx}$ | $\sigma_{h}$ | $i$ | $C_{2z}$ | $C_{2y}$ | $C_{2x}$
---|---|---|---|---|---|---|---
$\mathcal{E}_{xx}+\mathcal{E}_{yy}$ ($\mathcal{E}_{A_{1}}$) |  |  |  |  |  |  |
$\mathcal{E}_{xx}-\mathcal{E}_{yy}$ ($\mathcal{E}_{B_{1}}$) | $-\mathcal{E}_{xx}+\mathcal{E}_{yy}$ |  |  |  |  |  |
$\mathcal{E}_{xy}$ ($\mathcal{E}_{B_{2}}$) | $-\mathcal{E}_{xy}$ | $-\mathcal{E}_{xy}$ |  |  |  | $-\mathcal{E}_{xy}$ | $-\mathcal{E}_{xy}$
$\\{\mathcal{E}_{xz},\,\mathcal{E}_{yz}\\}$ ($\mathcal{E}_{E}$) | $\\{-\mathcal{E}_{yz},\,\mathcal{E}_{xz}\\}$ | $\\{-\mathcal{E}_{xz},\,\mathcal{E}_{yz}\\}$ | $\\{-\mathcal{E}_{xz},\,-\mathcal{E}_{yz}\\}$ |  | $\\{-\mathcal{E}_{xz},\,-\mathcal{E}_{yz}\\}$ | $\\{\mathcal{E}_{xz},\,-\mathcal{E}_{yz}\\}$ | $\\{-\mathcal{E}_{xz},\,\mathcal{E}_{yz}\\}$
$\\{S_{x},S_{y}\\}$ | $\\{-S_{y},S_{x}\\}$ | $\\{S_{x},-S_{y}\\}$ | $\\{-S_{x},-S_{y}\\}$ |  | $\\{-S_{x},-S_{y}\\}$ | $\\{-S_{x},S_{y}\\}$ | $\\{S_{x},-S_{y}\\}$
$S_{z}$ |  | $-S_{z}$ |  |  |  | $-S_{z}$ | $-S_{z}$
$\\{x,y\\}$ | $\\{-y,x\\}$ | $\\{-x,y\\}$ |  | $\\{-x,-y\\}$ | $\\{-x,-y\\}$ | $\\{-x,y\\}$ | $\\{x,-y\\}$
$z$ |  |  | $-z$ | $-z$ |  | $-z$ | $-z$
$n_{x}/m_{x}$ | $-n_{y}/-m_{y}$ |  | $-n_{x}/-m_{x}$ |  | $-n_{x}/-m_{x}$ | $-n_{x}/-m_{x}$ |
$n_{y}/m_{y}$ | $n_{x}/m_{x}$ | $-n_{y}/-m_{y}$ | $-n_{y}/-m_{y}$ |  | $-n_{y}/-m_{y}$ |  | $-n_{y}/-m_{y}$
$n_{z}/m_{z}$ |  | $-n_{z}/-m_{z}$ |  |  |  | $-n_{z}/-m_{z}$ | $-n_{z}/-m_{z}$
Table 2:  Symmetry transformation of the elastic strain tensor, spin vector,
lattice coordinates, and continuous spin fields. The first row lists the
important point group transformations to find the magnetic space group
generators discussed above. The blank space in the table denotes the variable
is invariant under the corresponding transformation. Note that from the
relations $iC_{2\alpha}=\sigma_{\alpha}$ and
$C_{2\alpha}C_{2\beta}=C_{2\gamma},C_{2\alpha}\sigma_{\beta}=\sigma_{\gamma}$,
where $C_{2\alpha}=C_{2x,2y,2z,2a,2b}$,
$\sigma_{\alpha}=\sigma_{vx,vy,h,da,db}$ and $\alpha,\beta,\gamma$ are
mutually orthogonal basis, other symmetry transformations not listed can be
generated.

###  B.3 Symmetry allowed spin-lattice coupling

Based on Table 2, it is straightforward to classify polynomials of
$n_{\mu},m_{\mu}$ that couples to the strain field by irreps. This is given in
Table 3. Note that we restricted our list to terms without spatial
derivatives, as these suffer additional suppression by temperature factors
well below the Debye temperature. Furthermore, we kept terms only to linear
order in $m_{x}$, because in the spin wave expansion, $m_{x}$, being
longitudinal, is already quadratic in the low energy transverse fields at
zeroth order in $h$.

| $n$ | $nn$ | $mm$ | $nmm$
---|---|---|---|---
$\mathcal{O}_{B_{1}}$ |  | $n_{y}n_{y},n_{z}n_{z}$ | $m_{y}m_{y},m_{z}m_{z}$ |  $n_{y}m_{x}m_{y}$,$n_{z}m_{x}m_{z}$
$\mathcal{O}_{B_{2}}$ | $n_{y}$ |  | $m_{x}m_{y}$ | $n_{y}m_{y/z}m_{y/z},n_{z}m_{y}m_{z}$
$\mathcal{O}_{E_{x}}$ | $n_{z}$ |  | $m_{x}m_{z}$ | $n_{z}m_{y/z}m_{y/z},n_{y}m_{y}m_{z}$
$\mathcal{O}_{E_{y}}$ |  | $n_{y}n_{z}$ | $m_{y}m_{z}$ |  $n_{y}m_{x}m_{z}$,$n_{z}m_{x}m_{y}$
Table 3:  Operators arranged by irrep for the high symmetry antiferromagnet
such that Eq. (2) in the main text is invariant under
$\mathsf{G}(\bm{\hat{x}},{\bm{0}})$. Other factors, e.g.
$S^{2},\chi,1/a_{0}^{2},n_{0}$, have been omitted in the table.

As described in the text, the operators can be expressed in terms of
$\bm{n}_{\scriptscriptstyle\perp}$ and $\bf{m}_{\scriptscriptstyle\perp}$
using the NLSM constraints, after taking into account an external Zeeman field
$\bm{h}$. This leads to the forms in Table 4 by replacing $m_{z}=\chi h+{\rm
m}_{z}$, $m_{x}={\rm m}_{x}=-\chi n_{z}h$.

| $n$ | $nn$ | $mm$ | $nmm$
---|---|---|---|---
$\mathcal{O}_{B_{1}}$ |  | $n_{y}n_{y},n_{z}n_{z}$ | $h\,{\rm m}_{z}$ |
$\mathcal{O}_{B_{2}}$ | $n_{y}$ |  | $-h\,n_{z}{\rm m}_{y}$ | $h\,n_{y}{\rm m}_{z},h\,n_{z}{\rm m}_{y}$
$\mathcal{O}_{E_{x}}$ | $n_{z}$ |  | $-h\,n_{z}{\rm m}_{z}$ | $h\,n_{z}{\rm m}_{z},h\,n_{y}{\rm m}_{y}$
$\mathcal{O}_{E_{y}}$ |  | $n_{y}n_{z}$ | $h\,{\rm m}_{y}$ |
Table 4: Magnetic operators up to quadratic order in the transverse spin wave
fluctuations in the high symmetry AFM in the presence of a small field along
the $z$-axis. Other factors, e.g. $S^{2},\chi,1/a_{0}^{2},n_{0}$, have been
omitted in the table. See Eq. (20) for the complete expression of
$\mathcal{O}$ at linear order in $n_{\mu},{\rm m}_{\mu}$.

##  C Experimental Implications

Our starting point is the effective phonon Lagrangian. In Fourier space, it
reads
$\mathcal{L}_{ph}=\sum_{\bm{q}}\mathcal{L}^{(0)}_{ph}({\bm{q}})+\mathcal{L}_{\rm
PHV}({\bm{q}})$, where

$\displaystyle\mathcal{L}_{ph}({\bm{q}})=$
$\displaystyle\frac{1}{2}\rho\dot{{\bm{u}}}^{T}_{-{\bm{q}}}\dot{{\bm{u}}}_{{\bm{q}}}-\frac{1}{2}{\bm{u}}^{T}_{-{\bm{q}}}\mathsf{M}_{\bm{q}}{\bm{u}}_{{\bm{q}}}+\dot{{\bm{u}}}^{T}_{-{\bm{q}}}\mathsf{A}_{\bm{q}}{\bm{u}}_{{\bm{q}}}.$
(25)

Here, $\rho$ is the lattice mass density, the sans serif font denotes a matrix
in Euclidean space, the first and second terms are the harmonic acoustic
phonon Lagrangian $\mathcal{L}^{(0)}_{ph}$, and the third term comes from the
phonon Hall viscosity and $\mathsf{A}_{\bm{q}}=-\mathsf{A}_{\bm{q}}^{T}$. For
simplicity, we will ignore the anisotropy in $\mathcal{L}^{(0)}_{ph}$ in the
evaluation of the phonon Berry curvature. $\mathsf{M}_{\bm{q}}$ for a 3D
isotropic elastic medium is

$\displaystyle\mathsf{M}_{\bm{q}}=\begin{pmatrix}c_{1}q_{x}^{2}+c_{2}q^{2}&c_{1}q_{x}q_{y}&c_{1}q_{x}q_{z}\\\
c_{1}q_{x}q_{y}&c_{1}q_{y}^{2}+c_{2}q^{2}&c_{1}q_{y}q_{z}\\\
c_{1}q_{x}q_{z}&c_{1}q_{y}q_{z}&c_{1}q_{z}^{2}+c_{2}q^{2}\\\ \end{pmatrix}.$
(26)

Here $c_{1,2}$ are the elastic modulus tensor coefficients. The eigenmodes
include two degenerate transverse acoustic waves with sound wave velocity
$v_{T}=\sqrt{c_{2}/\rho}$, and one longitudinal wave with
$v_{L}=\sqrt{(c_{1}+c_{2})/\rho}$. However, the anisotropy of a crystal with
lower symmetry breaks the degeneracy and mixes transverse and longitudinal
waves at a generic momentum.

The antisymmetric $\mathsf{A}$ matrices for $\eta^{H}_{B_{1}B_{2}}$,
$\eta^{H}_{E_{x}E_{y}}$ (relevant when the field is along z-axis) are

$\displaystyle\mathsf{A}_{\bm{q}}^{B_{1}B_{2}}=\eta^{H}_{B_{1}B_{2}}({\bm{q}})\begin{pmatrix}0&(q_{x}^{2}+q_{y}^{2})&0\\\
-(q_{x}^{2}+q_{y}^{2})&0&0\\\ 0&0&0\\\
\end{pmatrix},\qquad\mathsf{A}_{\bm{q}}^{E_{x}E_{y}}=\eta^{H}_{E_{x}E_{y}}({\bm{q}})\begin{pmatrix}0&q_{z}^{2}&q_{y}q_{z}\\\
-q_{z}^{2}&0&-q_{x}q_{z}\\\ -q_{y}q_{z}&q_{x}q_{z}&0\\\ \end{pmatrix},$ (27)

where $\eta^{H}_{\Gamma\Gamma^{\prime}}({\bm{q}})$ have been obtained in the
main text. We reproduce the result here:
$\eta^{H}_{B_{1}B_{2}}({\bm{q}})=\frac{\gamma_{B_{1}B_{2}}}{(q_{x}^{2}+q_{y}^{2})+\delta_{y}^{2}}$,
$\eta^{H}_{E_{x}E_{y}}({\bm{q}})=\frac{\gamma_{E_{x}E_{y}}}{(q_{x}^{2}+q_{y}^{2})+\delta_{z}^{2}}$,
where $\delta_{y}=\Delta_{y}/v_{m}$,
$\delta_{z}=\sqrt{\Delta_{z}^{2}+h^{2}}/v_{m}$, and
$\gamma_{\Gamma\Gamma^{\prime}}=\frac{hS^{2}\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}{v^{3}_{m}gd_{z}a_{0}}$.

###  C.1 Acoustic Faraday effect

The acoustic Faraday effect can be observed in an anisotropic medium only when
the acoustic wave is propagating along a high symmetry direction, such that
the transverse waves remain degenerate and do not mix with the longitudinal
one. In our case, for a tetragonal lattice crystal with ${\rm D}_{4h}$ point
group symmetry, there are indeed two degenerate transverse modes and one
longitudinal mode at ${\bm{q}}=q_{z}\bm{\hat{z}}$, due to the $C_{4z}$
symmetry in the little group of the high symmetry line. Upon applying a
magnetic field along $\bm{\hat{z}}$, the magnetoelastic-coupling-induced
phonon Hall viscosity in the $\eta^{H}_{E_{x}E_{y}}$ channel lifts the
degeneracy, and the left/right circularly polarized components are eigenmodes
and non-degenerate.

To be specific, the Lagrangian of a wave with frequency $\omega$ propagating
along $\bm{\hat{z}}$ with ${\bm{q}}=q\bm{\hat{z}}$ can be expressed as:

$\displaystyle\mathcal{L}_{ph}(q\bm{\hat{z}})=\frac{1}{2}{\bm{u}}^{T}_{-{\bm{q}}}\begin{pmatrix}\rho\omega^{2}-c_{2}q^{2}&i2\eta^{H}_{E_{x}E_{y}}q^{2}\omega&\\\
-i2\eta^{H}_{E_{x}E_{y}}q^{2}\omega&\rho\omega^{2}-c_{2}q^{2}&\\\
&&\rho\omega^{2}-(c_{1}+c_{2})q^{2}\\\ \end{pmatrix}{\bm{u}}_{{\bm{q}}}.$ (28)

While Eq. (28) is obtained from Eq. (25), which is only valid for an isotropic
medium, it also describes the acoustic phonon along $\bm{\hat{z}}$ for a
tetragonal group. The transverse modes can be diagonalized as left
($\bm{\hat{e}}_{+}=\frac{1}{\sqrt{2}}\\{1,i,0\\}^{T}$) and right
($\bm{\hat{e}}_{-}=\frac{1}{\sqrt{2}}\\{1,-i,0\\}^{T}$) circularly-polarized
waves, with $\rho\omega^{2}-c_{2}q_{\pm}^{2}\mp
2\eta^{H}_{E_{x}E_{y}}q_{\pm}^{2}\omega=0$. Consequently, the transverse
linearly-polarized wave along $\bm{\hat{z}}$ at frequency $\omega_{\rm ph}$
undergoes a Faraday rotation. The Faraday rotation angle per unit length,
which is given by the difference of the left and right wave numbers, $q_{+}$
and $q_{-}$ respectively, is

$\displaystyle\frac{\Phi}{L}=\frac{1}{2}(q_{+}-q_{-})=\frac{\omega_{\rm
ph}^{2}\eta^{H}_{E_{x}E_{y}}(q\bm{\hat{z}})}{v_{T}^{3}\rho}+{O}(\omega_{\rm
ph}^{3}),$ (29)

where $v_{T}=\sqrt{c_{2}/\rho}$ is the transverse sound wave velocity along
$\bm{\hat{z}}$. Note that the magnon spectrum is dispersionless along
$q\bm{\hat{z}}$ due to the weak interlayer spin interactions, so
$\eta^{H}_{E_{x}E_{y}}(q\bm{\hat{z}})=\frac{\gamma_{E_{x}E_{y}}}{\delta_{z}^{2}}$
from Eq. (8). This gives Eq. (9) in the main text.

###  C.2 Phonon thermal Hall conductivity

In this section, we show more details of the calculation of $\kappa^{H}(T)$.
To compute the Berry curvature of the phonon bands, we first construct the
effective phonon Hamiltonian $H_{ph}^{\text{eff}}$ from the Legendre
transformation of the Lagrangian $\mathcal{L}_{ph}$. The equation of motion
for the pair of canonical variables
$\psi_{\bm{q}}=\\{{\bm{u}}_{\bm{q}},{\bm{p}}_{\bm{q}}\\}^{T}$ can be obtained
from $\dot{\psi}_{\bm{q}}=-i[\psi_{\bm{q}},H_{ph}^{\text{eff}}]$ ($\hbar=1$).
This gives the band Hamiltonian $\mathcal{H}_{\bm{q}}$ in the basis of
$\psi_{\bm{q}}$. The phonon band energies/eigenstates can then be computed as
eigenvalues/eigenvectors of $\mathcal{H}_{\bm{q}}$.

To be specific, we find

$\displaystyle H_{ph}^{\text{eff}}$
$\displaystyle=\sum_{\bm{q}}\frac{\left({\bm{p}}_{-{\bm{q}}}-{\bm{A}}_{-{\bm{q}}}\right)\left({\bm{p}}_{{\bm{q}}}-{\bm{A}}_{{\bm{q}}}\right)}{2\rho}+\frac{1}{2}{\bm{u}}^{T}_{-{\bm{q}}}\mathsf{M}_{\bm{q}}{\bm{u}}_{{\bm{q}}}$
(30)

where the canonical momentum is defined as
${\bm{p}}_{{\bm{q}}}=\frac{\partial\mathcal{L}_{ph}}{\partial\dot{{\bm{u}}}_{-{\bm{q}}}^{i}}=\rho\dot{{\bm{u}}}_{{\bm{q}}}+\mathsf{A}_{\bm{q}}{\bm{u}}_{\bm{q}}=\rho\dot{{\bm{u}}}_{{\bm{q}}}+{\bm{A}}_{\bm{q}}$.
We then obtain the equation of motion in matrix form as
$\dot{\psi}_{\bm{q}}=-i\mathcal{H}_{\bm{q}}\psi_{\bm{q}}$, where

$\displaystyle\mathcal{H}_{\bm{q}}=\frac{i}{\rho}\begin{pmatrix}-\mathsf{A}_{\bm{q}}&\mathbb{I}\\\
-\rho\mathsf{M}_{\bm{q}}&-\mathsf{A}_{\bm{q}}\end{pmatrix}.$ (31)

Defining $\xi^{L(R)}_{{\bm{q}},\sigma}$ as the left/right eigenvector of
$\mathcal{H}_{\bm{q}}$ for the phonon branch $\sigma$ with eigenvalue/spectrum
$\upomega_{{\bm{q}},\sigma}$, the Berry curvature is

$\displaystyle\varOmega_{{\bm{q}},\sigma}^{i}=-\operatorname{Im}\left[\epsilon_{ijk}\partial_{q_{j}}\xi^{L}_{{\bm{q}},\sigma}\partial_{q_{k}}\xi^{R}_{{\bm{q}},\sigma}\right].$
(32)

As $\mathcal{H}_{\bm{q}}$ is not Hermitian, $\xi^{L(R)}_{{\bm{q}},\sigma}$ are
not Hermitian conjugate with each other, and they are normalized with
$\xi^{L}_{{\bm{q}},\sigma^{\prime}}\xi^{R}_{{\bm{q}},\sigma}=\delta_{\sigma\sigma^{\prime}}$.

The relation between the phonon Berry curvature and the intrinsic thermal Hall
conductivity has been obtained in Ref. Qin  _et al._ (2012).
$\kappa^{\text{ph}}_{jk}$ in terms of $\varOmega^{i}_{{\bm{q}},\sigma}$ and
$\upomega_{{\bm{q}},\sigma}$ is

$\displaystyle\frac{\kappa^{\text{ph}}_{jk}}{T}=-\frac{1}{T^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}E\,E^{2}\sigma_{jk}(E)\frac{\mathop{}\\!\mathrm{d}n_{B}^{\text{eq}}(E)}{\mathop{}\\!\mathrm{d}E}=\int_{0}^{\infty}\mathop{}\\!\mathrm{d}xx^{2}\frac{e^{x}}{(e^{x}-1)^{2}}\sigma_{jk}(xT)$
(33)

where

$\displaystyle\sigma_{jk}(E)=-\sum_{\sigma}\int\frac{\mathop{}\\!\mathrm{d}^{3}q}{(2\pi)^{3}}\epsilon_{ijk}\varOmega^{i}_{{\bm{q}},\sigma}\Theta(E-\upomega_{{\bm{q}},\sigma})=-\frac{1}{4\pi^{2}}\int_{0}^{\epsilon}\mathop{}\\!\mathrm{d}E^{\prime}\,E^{\prime
2}\sum_{\sigma}\int\frac{\mathop{}\\!\mathrm{d}\phi_{\bm{q}}}{2\pi}\mathop{}\\!\mathrm{d}\theta_{\bm{q}}\sin\theta_{\bm{q}}\frac{1}{v_{\sigma}^{3}}\varOmega^{i}_{{\bm{q}},\sigma}\Theta(E-E^{\prime}).$
(34)

$n_{B}^{\text{eq}}(E)=1/(e^{\beta E}-1)$ is the equilibrium Bose distribution
function. $\Sigma_{\sigma}$ sums over all eigenstates with positive energy.
$\theta_{{\bm{q}}},\phi_{\bm{q}}$ are the polar coordinates for ${\bm{q}}$,
defined through
${\bm{q}}=\\{q\sin\theta_{\bm{q}}\cos\phi_{\bm{q}},q\sin\theta_{\bm{q}}\sin\phi_{\bm{q}},q\cos\theta_{\bm{q}}\\}$.

Importantly, the behavior of $\kappa^{H}/T$ v.s. temperature is determined by
the functional form of $\sigma_{jk}(E)$ as we analyze below. As an example, we
compute $\kappa^{H}_{xy}(T)/T$ when the field is ${\bm{h}}=h\bm{\hat{z}}$
perturbatively in the Hall viscosity coefficient, keeping the first order in
$\eta^{H}$. The antisymmetric $\mathsf{A}$ matrices for
$\eta^{H}_{B_{1}B_{2}}$, $\eta^{H}_{E_{x}E_{y}}$ are presented in Eq. (27).

To linear order in $\eta^{H}$, the contributions from different Hall viscosity
terms simply add up. We find the Berry curvature
$\varOmega^{z,{\Gamma\Gamma^{\prime}}}_{{\bm{q}},\sigma}$ from
$\eta^{H}_{{\Gamma\Gamma^{\prime}}}$ of the two degenerate transverse phonon
bands $\sigma=T_{1},T_{2}$ and one longitudinal phonon band $\sigma=L$ as

$\displaystyle\varOmega^{z,B_{1}B_{2}}_{{\bm{q}},T_{1}}+\varOmega^{z,B_{1}B_{2}}_{{\bm{q}},T_{2}}=\frac{\gamma_{B_{1}B_{2}}v_{T}^{2}\sin^{2}\theta_{{\bm{q}}}\left(v_{L}^{2}+3v_{T}^{2}\right)\left[-4\cos
2\theta_{{\bm{q}}}\left(3\delta_{y}^{2}v_{T}^{2}+(v_{T}q)^{2}\right)-20\delta_{y}^{2}v_{T}^{2}+3(v_{T}q)^{2}\cos
4\theta_{{\bm{q}}}+(v_{T}q)^{2}\right]}{2\rho
qv_{T}(v_{T}-v_{L})(v_{L}+v_{T})\left(2\delta_{y}^{2}v_{L}^{2}-(v_{T}q)^{2}\cos
2\theta_{{\bm{q}}}+(v_{T}q)^{2}\right)^{2}}$
$\displaystyle\varOmega^{z,B_{1}B_{2}}_{{\bm{q}},L}=\frac{\gamma_{B_{1}B_{2}}v_{L}^{2}\sin^{2}\theta_{{\bm{q}}}\left(3v_{L}^{2}+v_{T}^{2}\right)\left[-4\cos
2\theta_{{\bm{q}}}\left(3\delta_{y}^{2}v_{L}^{2}+(v_{L}q)^{2}\right)-20\delta_{y}^{2}v_{L}^{2}+3(v_{L}q)^{2}\cos
4\theta_{{\bm{q}}}+(v_{L}q)^{2}\right]}{2\rho
v_{L}q(v_{L}-v_{T})(v_{L}+v_{T})\left(2\delta_{z}^{2}v_{L}^{2}-(v_{L}q)^{2}\cos
2\theta_{{\bm{q}}}+(v_{L}q)^{2}\right)^{2}}$
$\displaystyle\varOmega^{z,E_{x}E_{y}}_{{\bm{q}},T_{1}}+\varOmega^{z,E_{x}E_{y}}_{{\bm{q}},T_{2}}=\frac{\gamma_{E_{x}E_{y}}v_{T}^{2}\cos^{2}\theta_{{\bm{q}}}\left(v_{L}^{2}+3v_{T}^{2}\right)\left[-12\cos
2\theta_{{\bm{q}}}\left(\delta_{z}^{2}v_{T}^{2}+(v_{T}q)^{2}\right)-4\delta_{z}^{2}v_{T}^{2}+3(v_{T}q)^{2}\cos
4\theta_{{\bm{q}}}+9(v_{T}q)^{2}\right]}{\rho
v_{T}q(v_{T}-v_{L})(v_{L}+v_{T})\left(2\delta_{z}^{2}v_{T}^{2}-(v_{T}q)^{2}\cos
2\theta_{{\bm{q}}}+(v_{T}q)^{2}\right)^{2}}$
$\displaystyle\varOmega^{z,E_{x}E_{y}}_{{\bm{q}},L}=\frac{\gamma_{E_{x}E_{y}}v_{L}^{2}\cos^{2}\theta_{{\bm{q}}}\left(3v_{L}^{2}+v_{T}^{2}\right)\left[-12\cos
2\theta_{{\bm{q}}}\left(\delta_{z}^{2}v_{L}^{2}+(v_{L}q)^{2}\right)-4\delta_{z}^{2}v_{L}^{2}+3(v_{L}q)^{2}\cos
4\theta_{{\bm{q}}}+9(v_{L}q)^{2}\right]}{\rho
v_{L}q(v_{L}-v_{T})(v_{L}+v_{T})\left(2\delta_{z}^{2}v_{L}^{2}-(v_{L}q)^{2}\cos
2\theta_{{\bm{q}}}+(v_{L}q)^{2}\right)^{2}},$ (35)

where $q=|{\bm{q}}|$. From Eqs. (35) and (34), we find

$\displaystyle\sigma_{xy}^{B_{1}B_{2}}(v_{\rm ph}\delta_{y}\,t)$
$\displaystyle=-\frac{4\gamma_{B_{1}B_{2}}}{v_{\rm
ph}\rho}\left(-\frac{1}{1+t^{2}}-\frac{\left(3+4t^{2}\right)\log\left(\sqrt{1+t^{2}}+t\right)}{t^{3}\left(1+t^{2}\right)^{3/2}}+\frac{3}{t^{2}}\right)=-\frac{4\gamma_{B_{1}B_{2}}}{v_{\rm
ph}\rho}f_{y}(t)$ $\displaystyle\sigma_{xy}^{E_{x}E_{y}}(v_{\rm
ph}\delta_{z}\,t)$ $\displaystyle=-\frac{8\gamma_{E_{x}E_{y}}}{v_{\rm
ph}\rho}\left(\frac{\left(3+2t^{2}\right)\log\left(\sqrt{1+t^{2}}+t\right)}{t^{3}\left(1+t^{2}\right)^{3/2}}-\frac{3}{t^{2}}\right)=-\frac{8\gamma_{E_{x}E_{y}}}{v_{\rm
ph}\rho}f_{z}(t).$ (36)

Corrections at order $(v_{L}-v_{T})/v_{L}$ are subleading and nonsingular, and
so they are ignored in the above expressions. $v_{\rm ph}$ denotes the
averaged sound velocity. Note that $\sigma_{xy}^{\Gamma\Gamma^{\prime}}(E)$ is
only a function of $t=E/(v_{\rm ph}\delta_{\alpha(\Gamma\Gamma^{\prime})})$,
where $\alpha(B_{1},B_{2})=y,\alpha(E_{x},E_{y})=z$. For both
$\sigma_{xy}^{B_{1}B_{2}}$, $\sigma_{xy}^{E_{x}E_{y}}$, the function in the
parentheses, defined as $f_{y,z}(t)$, increases as $\sim t^{2}$ for $t\ll 1$,
and decreases to zero as $\sim\frac{1}{t^{2}}$ for $f_{y}(t)$, and as
$\sim\frac{\ln t}{t^{2}}$ for $f_{z}(t)$ at $t\gg 1$ (see Fig. 3).

Figure 3: Behavior of $f_{y,z}(t)$ defined in Eq. (36).

From Eq. (33) and (36), defining a dimensionless temperature
$\tilde{T}=T/(v_{\rm ph}\delta_{\alpha})$,

$\displaystyle\frac{\kappa_{xy}^{\Gamma\Gamma^{\prime}}}{T}\sim
F_{\alpha(\Gamma\Gamma^{\prime})}\left(\tilde{T}\right),\qquad
F_{\alpha(\Gamma\Gamma^{\prime})}(\tilde{T})=\int_{0}^{\infty}\mathop{}\\!\mathrm{d}xx^{2}\frac{e^{x}}{(e^{x}-1)^{2}}f_{\alpha}(x\tilde{T}).$
(37)

It is straightforward to see that the behavior of $\kappa^{H}/T$ is only a
function of $\tilde{T}=T/(v_{\rm ph}\delta_{\alpha})$. Also, noting that the
Boltzmann distribution decays exponentially at high temperature,
$\frac{\kappa_{xy}^{\Gamma\Gamma^{\prime}}}{T}$ at $\tilde{T}\ll 1$($T\ll
v_{\rm ph}\delta_{\alpha}$) and $\tilde{T}\gg 1$($T\gg v_{\rm
ph}\delta_{\alpha}$) may be inferred from $f_{\alpha}(t)$ analytically. When
$\tilde{T}\ll 1$, it mainly comes from $f_{\alpha}(t)\sim t^{2}$ at $t\ll 1$,
so we have $\frac{\kappa_{xy}^{\Gamma\Gamma^{\prime}}}{T}\sim T^{2}$. When
$\tilde{T}\gg 1$, the scaling of
$\frac{\kappa_{xy}^{\Gamma\Gamma^{\prime}}}{T}$ is more complicated as it
depends on $f_{\alpha}(t)$ in the whole range. From the decaying behavior of
$f_{\alpha}(t)$ when $t\gtrsim 1$, we find that
$\frac{\kappa_{xy}^{\Gamma\Gamma^{\prime}}}{T}\sim 1/T^{\zeta}$ at $T\gg
v_{\rm ph}\delta_{\alpha}$, with $\zeta\sim 1$ up to a logarithmic factor. In
the intermediate regime, $\frac{\kappa_{xy}^{\Gamma\Gamma^{\prime}}}{T}$ can
be fit to an exponential empirically as
$e^{-T/T_{0}^{{\Gamma\Gamma^{\prime}}}}+\text{constant}$. From the numerical
fit, we find $T_{0}^{B_{1}B_{2}}\sim 2.3v_{\rm ph}\delta_{y}$,
$T_{0}^{E_{x}E_{y}}\sim 3.3v_{\rm ph}\delta_{z}$.

##  D Contribution to $\eta^{H}$ from four-point correlation functions

In this section, we discuss contributions to $\eta^{H}$ from four-point
correlation functions of the spin fields. We argue that contributions at this
order are subleading at low temperature, and can thus be ignored.

We consider the contribution to $\eta^{H}_{B_{1}B_{2}}$
($\eta^{H}_{E_{x}E_{y}}$ is similar). From Eq. (3) and Tab. 3, we find

$\displaystyle\langle\mathcal{O}_{B_{1}}(-{\bm{q}},\tau)\mathcal{O}_{B_{2}}({\bm{q}},0)\rangle_{\omega_{n}}$
$\displaystyle\sim$
$\displaystyle\frac{1}{L^{2}}\sum_{{\bm{k}},{\bm{k}}^{\prime}}\left(\varsigma_{1}\langle
n_{y,-\bm{k}-{\bm{q}}_{\scriptscriptstyle\perp}}(\tau)n_{y,\bm{k}}(\tau)n_{y,\bm{k}^{\prime}+{\bm{q}}_{\scriptscriptstyle\perp}}(0){\rm
m}_{z,-\bm{k}^{\prime}}(0)\rangle_{\omega_{n}}+\varsigma_{2}\langle
n_{z,-\bm{k}-{\bm{q}}_{\scriptscriptstyle\perp}}(\tau)n_{z,\bm{k}}(\tau)n_{z,\bm{k}^{\prime}+{\bm{q}}_{\scriptscriptstyle\perp}}(0){\rm
m}_{y,-\bm{k}^{\prime}}(0)\rangle_{\omega_{n}}\right)$ $\displaystyle\sim$
$\displaystyle T\frac{1}{L^{2}}\sum_{{\bm{k}},m}\left(\varsigma_{1}\langle
n_{y,-\bm{k}-{\bm{q}}_{\scriptscriptstyle\perp}}{\rm
m}_{z,\bm{k}+{\bm{q}}_{\scriptscriptstyle\perp}}\rangle_{\omega_{m}}\langle
n_{y,\bm{k}}(\tau)n_{y,-\bm{k}}\rangle_{-\omega_{m}+\omega_{n}}+\varsigma_{2}\langle
n_{z,-\bm{k}-{\bm{q}}_{\scriptscriptstyle\perp}}{\rm
m}_{y,\bm{k}+{\bm{q}}_{\scriptscriptstyle\perp}}\rangle_{\omega_{m}}\langle
n_{z,\bm{k}}(\tau)n_{z,-\bm{k}}\rangle_{-\omega_{m}+\omega_{n}}\right).$ (38)

where $\varsigma_{1},\varsigma_{2}$ are dimensionless coefficients that
determine the magnitude and sign of the two types of four-point correlations.
Following the analysis in Sec. A.1, an overall factor
$\frac{\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}{d_{z}}\left(\frac{S^{2}}{a_{0}^{2}}\right)^{2}\chi
h$ is absorbed into $\sim$.

For simplicity, we consider the limit when $q=|{\bm{q}}|\rightarrow 0$. After
summing over the Matsubara frequencies $\omega_{m}$ and analytical continuing
to real frequencies, $\eta^{H}_{B_{1}B_{2}}$ from Eq. (38) can be expressed as

$\displaystyle\left[\eta_{B_{1},B_{2}}^{H}({\bm{q}})\right]^{(2)}\sim$
$\displaystyle\frac{\lambda_{\Gamma}\lambda_{\Gamma^{\prime}}}{d_{z}}\left(\frac{S^{2}}{a_{0}^{2}}\right)^{2}\chi
h\,\frac{a_{0}^{2}}{S}\,v_{m}ga_{0}$
$\displaystyle\frac{1}{a_{0}^{2}}\frac{1}{N_{x}N_{y}}\sum_{\bm{k}}\left\\{\frac{\varsigma_{1}}{(\upomega_{y,{\bm{k}}}+\upomega_{y,{\bm{k}}-{\bm{q}}_{\scriptscriptstyle\perp}})^{2}}\left(\frac{1}{\upomega_{y,{\bm{k}}}}\coth{\frac{\beta\upomega_{y,{\bm{k}}}}{2}}+\frac{1}{\upomega_{y,{\bm{k}}-{\bm{q}}_{\scriptscriptstyle\perp}}}\coth{\frac{\beta\upomega_{y,{\bm{k}}-{\bm{q}}_{\scriptscriptstyle\perp}}}{2}}\right)\right.$
$\displaystyle\qquad\qquad\quad\,\,\left.-\frac{\varsigma_{2}}{(\upomega_{z,{\bm{k}}}+\upomega_{z,{\bm{k}}-{\bm{q}}_{\scriptscriptstyle\perp}})^{2}}\left(\frac{1}{\upomega_{z,{\bm{k}}}}\coth{\frac{\beta\upomega_{z,{\bm{k}}}}{2}}+\frac{1}{\upomega_{z,{\bm{k}}-{\bm{q}}_{\scriptscriptstyle\perp}}}\coth{\frac{\beta\upomega_{z,{\bm{k}}-{\bm{q}}_{\scriptscriptstyle\perp}}}{2}}\right)\right\\}$
$\displaystyle\xrightarrow{q\rightarrow 0}$
$\displaystyle\gamma_{B_{1}B_{2}}v_{m}^{3}ga_{0}\frac{1}{a_{0}^{2}}\int_{-\frac{\pi}{a_{0}}}^{\frac{\pi}{a_{0}}}\int_{-\frac{\pi}{a_{0}}}^{\frac{\pi}{a_{0}}}\frac{\mathop{}\\!\mathrm{d}^{2}k}{(2\pi/a_{0})^{2}}\left\\{\frac{\varsigma_{1}}{\upomega_{y,{\bm{k}}}^{3}}\coth{\frac{\beta\upomega_{y,{\bm{k}}}}{2}}-\frac{\varsigma_{2}}{\upomega_{z,{\bm{k}}}^{3}}\coth{\frac{\beta\upomega_{z,{\bm{k}}}}{2}}\right\\}$
$\displaystyle\sim$
$\displaystyle\gamma_{B_{1}B_{2}}v_{m}^{3}ga_{0}\frac{1}{Tv_{m}^{2}}\left[\varsigma_{1}\mathcal{F}\left(\frac{v_{m}\delta_{y}}{2T}\right)-\varsigma_{2}\mathcal{F}\left(\frac{v_{m}\delta_{z}}{2T}\right)\right]$
$\displaystyle\sim$
$\displaystyle\gamma_{B_{1}B_{2}}v_{m}ga_{0}\begin{cases}\frac{1}{v_{m}\delta_{y}}&T\ll
v_{m}\delta_{y}\\\ \frac{T}{(v_{m}\delta_{y})^{2}}&T\gg
v_{m}\delta_{y}\end{cases}$ (39)

where
$\mathcal{F}(x)=\int_{x}^{\infty}\mathop{}\\!\mathrm{d}z\frac{1}{z^{2}}\coth
z$, and $\delta_{y}\ll\delta_{z}$ has been applied to obtain the last line. As
our interest here is the relative strength of four-point contribution compared
to two-point contribution, numerical coefficients are subsumed in the
“$\sim$”. The superscript $``(2)"$ denotes that it is a second order
contribution. Note that the contribution from 2 magnons is much smaller than
that from 1 magnon by a factor
$\frac{gT}{v_{m}a_{0}^{-1}}\sim\frac{T}{JS^{2}}\ll 1$. For this reason, only
contributions from 1 magnon terms are considered further in the main text.

##  E Microscopic Hamiltonian

In this section, we present the microscopic spin Hamiltonian and spin-lattice
coupling Hamiltonian. The microscopic analysis is helpful to infer the spin
gap and relevant spin-lattice coupling in the low energy effective action
derived in the main text from a symmetry analysis.

The microscopic Hamiltonian must be invariant under the grey magnetic space
group $\mathsf{G}$ (see Table 1 ) of the crystal. It is useful to decompose
the part involving the spin as:

$\displaystyle H$ $\displaystyle=H_{\rm s}+H_{\rm sl}$ $\displaystyle=H_{\rm
s}^{(0)}+\lambda_{SOC}H_{\rm s}^{(1)}+H_{\rm sl}^{(0)}+\lambda_{SOC}H_{\rm
sl}^{(1)},$ (40)

where ${\rm s}$ refers to “spin” and ${\rm sl}$ to “spin-lattice coupling”. We
assume weak SOC (as e.g. appropriate for cuprates), and accordingly separate
terms of 0th order in SOC ($H_{\rm s}^{(0)}+H_{\rm sl}^{(0)}$) and those that
require SOC ($H_{\rm s}^{(1)}+H_{\rm sl}^{(1)}$). This is defined by symmetry:
$H_{\rm s}^{(0)}+H_{\rm sl}^{(0)}$ also has a global $SO(3)$ spin-rotation
symmetry, so is also invariant under $SO(3)_{s}$, while $H_{\rm
s}^{(1)}+H_{\rm sl}^{(1)}$ must break $SO(3)_{s}$.

_$H_{\rm s}$_ –The spin Hamiltonian reads

$\displaystyle H_{s}^{(0)}$
$\displaystyle=\sum_{i}J({\bm{S}}_{\bf{r}}\cdot{\bm{S}}_{{\bf{r}}\pm\hat{{\bm{x}}}}+{\bm{S}}_{\bf{r}}\cdot{\bm{S}}_{{\bf{r}}\pm\hat{{\bm{y}}}})$
$\displaystyle H_{s}^{(1)}$
$\displaystyle=\sum_{\bf{r}}J^{\prime}_{1}(S^{x}_{\bf{r}}S^{x}_{{\bf{r}}\pm\bm{\hat{x}}}-S^{y}_{\bf{r}}S^{y}_{{\bf{r}}\pm\bm{\hat{x}}}-S^{x}_{\bf{r}}S^{x}_{{\bf{r}}\pm\bm{\hat{y}}}+S^{y}_{\bf{r}}S^{y}_{{\bf{r}}\pm\bm{\hat{y}}})+J^{\prime}_{z}(S^{z}_{\bf{r}}S^{z}_{{\bf{r}}+\bm{\hat{x}}}+S^{z}_{\bf{r}}S^{z}_{{\bf{r}}+\bm{\hat{y}}})$
$\displaystyle\quad+J^{\prime}_{2}(S^{x}_{\bf{r}}S^{y}_{{\bf{r}}\pm\bm{\hat{a}}}+S^{y}_{\bf{r}}S^{x}_{{\bf{r}}\pm\bm{\hat{a}}}-S^{x}_{\bf{r}}S^{y}_{{\bf{r}}\pm\bm{\hat{b}}}-S^{y}_{\bf{r}}S^{x}_{{\bf{r}}\pm\bm{\hat{b}}})$
(41)

where, as before,
$\bm{\hat{a},\hat{b}}=\frac{1}{\sqrt{2}}(\pm\bm{\hat{x}}+\bm{\hat{y}})$ are
second neighbor vectors. The nearest neighbor exchange dominates and $J\sim
0.1$ eV. From Refs. Yildirim  _et al._ (1995); Benfatto  _et al._ (2006), the
anisotropy gap $\Delta_{y,z}$ can be estimated for Sr2CuO2Cl2 from the spin
exchange interactions that require SOC. The XXZ anisotropic exchange
$J^{\prime}_{z}$ is positive, which stabilizes the staggered order in the $xy$
plane. $J^{\prime}_{z}$ introduces a gap for out-of-plane fluctuations
$\Delta_{z}=4SJ^{\prime}_{z}\sim 3$ meV, where the factor of $4$ comes from
the coordination number of square lattice. While $J^{\prime}_{1,2}\ll J$,
$J^{\prime}_{1}$ is still important to introduce an effective mass gap to the
in-plane Goldstone mode at order $1/S$. In Ref. Yildirim  _et al._ (1995), by
studying the dependence of the quantum zero point energy on the angle $\phi$
of the Néel order relative to $\bm{\hat{x}}$, it is shown that
$J^{\prime}_{1}$ stabilize Néel order with $\phi=\pi/4\text{ mod }\pi/2$, and
thus the in-plane Goldstone mode acquires an effective gap, which is estimated
to be $\Delta_{y}=0.05$ meV when $S=1/2$.

_$H_{\rm sl}$_ – In this work, we assume the spin-lattice coupling comes from
the modification to the exchange coupling when the direction and length of the
bond change due to an elastic deformation. The microscopic spin-lattice
Hamiltonian can be obtained by finding all the terms composed of
$\mathcal{E}_{\Gamma}$ and $\mathcal{O}_{\Gamma}$ invariant under
$\mathsf{G}$, where $\mathcal{E}_{\Gamma}$ is the strain tensor decomposed
into irreps of ${\rm D}_{4h}$ (see Eq. (6) in the main text), and on the
lattice we take $\mathcal{O}_{\Gamma}$ to be composed of bilinears of spin
operators on nearby sites, chosen to transform under ${\rm D}_{4h}$ with the
same irrep $\Gamma$. In real space, we find

$\displaystyle H_{\rm sl}^{(0)}=$
$\displaystyle\lambda_{A_{1}}\sum_{{\bf{r}}}(\mathcal{E}_{A_{1}})_{{\bf{r}}}({\bm{S}}_{{\bf{r}}}\cdot{\bm{S}}_{{\bf{r}}\pm\bm{\hat{x}}}+{\bm{S}}_{{\bf{r}}}\cdot{\bm{S}}_{{\bf{r}}\pm\bm{\hat{y}}})+\frac{\lambda_{B1}}{2}\sum_{{\bf{r}}}(\mathcal{E}_{B_{1}})_{{\bf{r}}}\left({\bm{S}}_{{\bf{r}}}\cdot{\bm{S}}_{{\bf{r}}+\bm{\hat{x}}}-{\bm{S}}_{{\bf{r}}}\cdot{\bm{S}}_{{\bf{r}}+\bm{\hat{y}}}\right)+$
$\displaystyle\frac{\lambda_{B2}}{2}\sum_{{\bf{r}}}(\mathcal{E}_{B_{2}})_{{\bf{r}}}\left({\bm{S}}_{{\bf{r}}}\cdot{\bm{S}}_{{\bf{r}}+\bm{\hat{a}}}-{\bm{S}}_{{\bf{r}}}\cdot{\bm{S}}_{{\bf{r}}+\bm{\hat{b}}}\right).$
(42)

Here we denote by $(\mathcal{E}_{\Gamma})_{\bm{x}}$ the strain at position
$\bm{x}$ – note that since the strain is slowly varying, shifts of this
coordinate by order one displacements do not modify the results. The XXZ
anisotropy may be ignored here in the spin-lattice coupling as it is
negligibly small and not essential (it is, however, important in the spin
Hamiltonian itself). The anisotropic spin-lattice coupling terms, keeping spin
operators on the nearest possible pairs of spins in each channel, give $H_{\rm
sl}^{(1)}$ in the form

$\displaystyle H_{\rm sl}^{(1)}=$
$\displaystyle\frac{\lambda^{\prime}_{B1}}{2}\sum_{{\bf{r}}}(\mathcal{E}_{B_{1}})_{{\bf{r}}}\left(S^{x}_{{\bf{r}}}S^{x}_{{\bf{r}}\pm{\bm{\hat{x}}}}-S^{y}_{r}S^{y}_{{\bf{r}}\pm{\bm{\hat{y}}}}+S^{x}_{{\bf{r}}}S^{x}_{{\bf{r}}\pm{\bm{\hat{y}}}}-S^{y}_{{\bf{r}}}S^{y}_{{\bf{r}}\pm{\bm{\hat{x}}}}\right)+$
$\displaystyle\frac{\lambda^{\prime}_{B2}}{2}\sum_{{\bf{r}}}(\mathcal{E}_{B_{2}})_{{\bf{r}}}\left(S^{x}_{{\bf{r}}}S^{y}_{{\bf{r}}\pm\bm{\hat{a}}}+S^{y}_{{\bf{r}}}S^{x}_{{\bf{r}}\pm\bm{\hat{a}}}+S^{x}_{{\bf{r}}}S^{y}_{{\bf{r}}\pm\bm{\hat{b}}}+S^{y}_{{\bf{r}}}S^{x}_{{\bf{r}}\pm\bm{\hat{b}}}\right)+$
$\displaystyle\frac{\lambda^{\prime}_{E}}{2}\sum_{{\bf{r}}}(\mathcal{E}_{E_{x}})_{{\bf{r}}}\left[(1+\beta)\left(S^{x}_{{\bf{r}}}S^{z}_{{\bf{r}}\pm{\bm{\hat{x}}}}+S^{z}_{{\bf{r}}}S^{x}_{{{\bf{r}}}\pm{\bm{\hat{x}}}}\right)+(1-\beta)\left(S^{x}_{{\bf{r}}}S^{z}_{{{\bf{r}}}\pm{\bm{\hat{y}}}}+S^{z}_{{\bf{r}}}S^{x}_{{{\bf{r}}}\pm{\bm{\hat{y}}}}\right)\right]+$
$\displaystyle\frac{\lambda^{\prime}_{E}}{2}\sum_{{\bf{r}}}(\mathcal{E}_{E_{y}})_{{\bf{r}}}\left[(1-\beta)\left(S^{y}_{{\bf{r}}}S^{z}_{{{\bf{r}}}\pm{\bm{\hat{x}}}}+S^{z}_{{\bf{r}}}S^{y}_{{{\bf{r}}}\pm{\bm{\hat{x}}}}\right)+(1+\beta)\left(S^{y}_{{\bf{r}}}S^{z}_{{{\bf{r}}}\pm{\bm{\hat{y}}}}+S^{z}_{{\bf{r}}}S^{y}_{{{\bf{r}}}\pm{\bm{\hat{y}}}}\right)\right].$
(43)

Note that in the $\Gamma=E$ channel, there are two symmetry allowed spin-
lattice coupling terms, of strength $\lambda^{\prime}_{E}$ and
$\beta\lambda^{\prime}_{E}$, where $\beta$ is a dimensionless parameter
defined by Eq. (43).

The above microscopic spin-lattice couplings give rise to continuum spin-
lattice couplings of the form in Eq. (2) of the main text by taking a
continuum limit. In particular, one uses the NLSM decomposition with
$\bm{S}_{\bm{r}}=S\bm{e}_{\bm{r}}$ with $\bm{e}_{\bm{r}}$ in Eq. (14), and
expresses the Hamiltonian thereby in terms of $\bm{n}$ and $\bm{m}$, which are
presumed to be slowly varying functions of position. To zeroth order in the
gradient expansion of these fields, derivatives
$\partial_{\mu}\bm{n},\partial_{\mu}\bm{m}$ are neglected. The result then
takes the schematic form

$\displaystyle H_{\rm
sl}^{(loc)}\sim\sum_{z}\frac{1}{a_{0}^{2}}\int\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$
$\displaystyle\left[\lambda_{A_{1}}\mathcal{E}_{A_{1}}\bm{s}\cdot\bm{s}+\lambda^{\prime}_{B1}\mathcal{E}_{B_{1}}\left(s^{x}s^{x}-s^{y}s^{y}\right)+\lambda^{\prime}_{B2}\mathcal{E}_{B_{2}}\left(s^{x}s^{y}+s^{y}s^{x}\right)+\right.$
$\displaystyle\left.\lambda^{\prime}_{E}\mathcal{E}_{E_{x}}\left(s^{x}s^{z}+s^{z}s^{x}\right)+\lambda^{\prime}_{E}\mathcal{E}_{E_{y}}\left(s^{y}s^{z}+s^{z}s^{y}\right)\right].$
(44)

Here each term represents a sum of two contributions, one with
$\bm{s}=S\bm{m}$, and another with $\bm{s}=S\bm{n}$, and some sign differences
may appear between these terms. Here, we summed over 2d magnetic layers
($\sum_{z}$), where the layer index is $z$. In the limit of decoupled magnetic
layers, we considered the continuous limit as
$\sum_{\bf{r}}\rightarrow\sum_{z}\frac{1}{a_{0}^{2}}\int\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}y$
(see the discussion of Eq. (20) for more details). From Eq. (44), we see that
$H_{\rm sl}$ in all channels $\Gamma\neq A_{1}$ requires breaking $SO(3)_{s}$,
and thus requires SOC microscopically. Eq. (44) give results consistent with
the general symmetry analysis in terms of continuous spin fields ${\bm{n}}$
and $\bm{m}$ (c.f. Table 3).

##  F Phonon Hall viscosity for in-plane magnetic field

For completeness, we have also analyzed the phonon Hall viscosity for in-plane
magnetic fields. For concreteness, the field is applied along the $y$-axis.
Such a scenario can be relevant to the phonon Hall effect when the heat
current is applied perpendicular to the CuO2 planes, i.e. upon studying
$\kappa_{xz}$. The analysis follows that for the out-of-plane field presented
in the main text. Here, we summarize the main results.

The magnetic space group for a paramagnet and AFM with staggered order in the
$xy$ plane at an arbitrary $\phi\neq 0\text{ mod }\pi/4$ and $\phi=0$ relative
to $\bm{\hat{x}}$, when the field is along $\bm{\hat{y}}$, is summarized in
Table 5.

|  zero field | ${\bm{h}}=h\bm{\hat{y}}$
---|---|---
|  lattice and spin |  lattice effective
paramagnet | $\mathsf{G}=P4/mmm1^{\prime}$ | $\mathsf{G}(\bm{0},h\bm{\hat{y}})=\langle X,Y,C_{2y},i,\mathcal{T}{C_{2z}}\rangle$ | $\mathsf{G}^{\rm eff}(\bm{0},h\bm{\hat{y}})=\langle C_{2y},i,\mathcal{T}{C_{2z}}\rangle$
high sym. AFM | $\mathsf{G}(\bm{\hat{x}},\bm{0})=\langle i,\mathcal{T}X,\mathcal{T}Y,C_{2x},\mathcal{T}C_{2z}\rangle$ | $\mathsf{G}(\bm{\hat{x}},h\bm{\hat{y}})=\langle i,XY,\mathcal{T}C_{2z}\rangle$ | $\mathsf{G}^{\rm eff}(\bm{\hat{x}},h\bm{\hat{y}})=\langle i,\mathcal{T}C_{2z}\rangle$
low sym. AFM | $\mathsf{G}(\bm{\hat{e}},\bm{0})=\langle i,\mathcal{T}X,\mathcal{T}Y,\mathcal{T}C_{2z}\rangle$ | $\mathsf{G}(\bm{\hat{e}},h\bm{\hat{y}})=\langle i,XY,\mathcal{T}C_{2z}\rangle$ | $\mathsf{G}^{\rm eff}(\bm{\hat{e}},h\bm{\hat{y}})=\langle i,\mathcal{T}C_{2z}\rangle$
Table 5:  Summary of the magnetic space group (in the Hermann-Mauguin
notation) for layer group $G=P4/mmm$ with ${\bm{h}}=h\bm{\hat{y}}$. Same to
Tab. 1 in the main text, rows denote different values for $\bm{n}$: without
magnetic order, and for the two distinct orientations of the staggered
magnetization in the $xy$ plane (here $\bm{\hat{e}}$ is a generic vector not
along a high-symmetry axis in the plane), and columns specify the zero and
finite field cases, and the effective magnetic group for the effective lattice
theory. The symbol $\langle\cdot\rangle$ indicates the group generated by the
“$\cdot$” operations.

Note that $\eta^{H}_{E_{x}E_{y}}=0$ due to the $\mathcal{T}{C_{2z}}$ symetry,
so there is no AFE. To generate a nonzero phonon Berry curvature
$\varOmega^{y}$, both time-reversal and mirror symmetries perpendicular to the
$xz$ plane ($\sigma_{h}$ and $\sigma_{vx}$ in the ${\rm D}_{4h}$ group) should
be broken. The symmetry-allowed Hall viscosity that generates a nonzero phonon
Berry curvature reads

$\displaystyle\mathcal{S}_{{\rm
PHV},h\bm{\hat{y}}}=\int\mathop{}\\!\mathrm{d}^{3}x\mathop{}\\!\mathrm{d}\tau\,\\{\eta^{H}_{A_{1}E_{x}}(\mathcal{E}_{A_{1}}\dot{\mathcal{E}}_{E_{x}}-\dot{\mathcal{E}}_{A_{1}}\mathcal{E}_{E_{x}})+\eta^{H}_{B_{1}E_{x}}(\mathcal{E}_{B_{1}}\dot{\mathcal{E}}_{E_{x}}-\dot{\mathcal{E}}_{B_{1}}\mathcal{E}_{E_{x}})+\eta^{H}_{B_{2}E_{y}}(\mathcal{E}_{B_{2}}\dot{\mathcal{E}}_{E_{y}}-\dot{\mathcal{E}}_{B_{2}}\mathcal{E}_{E_{y}})\\}.$
(45)

To obtain the Hall viscosity induced by spin-lattice coupling, we find
$\mathcal{O}_{\Gamma}({\bm{x}},\tau)$ (see Eq. (2)) in terms of the spin
fields as tabulated in Table 6.

| $\epsilon\,n$ | $\epsilon\,nn$ | $\epsilon\,mm$ | $\epsilon\,nmm$
---|---|---|---|---
$\mathcal{O}_{A_{1}}$ |  | $n_{y}n_{y},n_{z}n_{z}$ | $h\,{\rm m}_{y}$ |
$\mathcal{O}_{B_{1}}$ |  | $n_{y}n_{y},n_{z}n_{z}$ | $h\,{\rm m}_{y}$ |
$\mathcal{O}_{B_{2}}$ | $n_{y}$ |  | $-h\,n_{y}{\rm m}_{y}$ | $h\,n_{y}{\rm m}_{y},h\,n_{z}{\rm m}_{z}$
$\mathcal{O}_{E_{x}}$ | $n_{z}$ |  | $-h\,n_{y}{\rm m}_{z}$ | $h\,n_{z}{\rm m}_{y},h\,n_{y}{\rm m}_{z}$
$\mathcal{O}_{E_{y}}$ |  | $n_{y}n_{z}$ | $h\,{\rm m}_{z}$ |

.

Table 6: Magnetic operators in the high symmetry AFM in the presence of a
small field along the $y$-axis. Other factors, e.g.
$S^{2},\chi,1/a_{0}^{2},n_{0}$, have been omitted in the table.

At leading order in $1/S$, this gives

$\displaystyle\eta_{A_{1},E_{x}}^{H}\sim h\langle{\rm
m}_{y}n_{z}\rangle,\quad\eta_{B_{1},E_{x}}^{H}\sim h\langle{\rm
m}_{y}n_{z}\rangle,\quad\eta_{B_{2},E_{y}}^{H}\sim h\langle{\rm
m}_{z}n_{y}\rangle.$ (46)

Note that at leading order in $1/S$, the Hall viscosity coefficients are
determined by the same sets of spin correlators as when
${\bm{h}}=h\bm{\hat{z}}$. However, this is restricted to when the in-plane
field is perpendicular to the staggered order. If the external field also has
components parallel to the staggered order, we find that $\langle
n_{y}n_{z}\rangle_{\omega}\sim ih\omega$ also contribute to a nonzero Hall
viscosity.