# Supplementary materials for “Hyperbolic Spin Waves in Magnetic Polar Metals”

Abhishek Kumar<EMAIL_ADDRESS>Department of Physics and
Astronomy, Rutgers University, Piscataway, New Jersey 08854, USA  Département
de physique and Institut quantique, Université de Sherbrooke, Sherbrooke,
Québec, Canada J1K 2R1  Premala Chandra  Department of Physics and Astronomy,
Rutgers University, Piscataway, New Jersey 08854, USA  Pavel A. Volkov
<EMAIL_ADDRESS>Department of Physics and Astronomy, Rutgers
University, Piscataway, New Jersey 08854, USA  Department of Physics, Harvard
University, Cambridge, Massachusetts, 02138 USA  Department of Physics,
University of Connecticut, Storrs, Connecticut 06269, USA

###### Contents

  1. I Dirac-point-vicinity model
    1. I.1 Operator RPA: an Equation of motion approach
    2. I.2 Spin-exciton bound state
    3. I.3 Self-energy
      1. I.3.1 Zero Rashba case: Silin-Leggett mode
      2. I.3.2 Finite Rashba case
    4. I.4 Wavefunction of the collective mode frequency
  2. II Generation of spin-current by in-plane component of the tilted magnetic field for Dirac-point-vicinity model
    1. II.1 Dispersion of the collective mode frequency in the presence of tilted magnetic field
    2. II.2 Wavefunction of the collective mode frequency in the presence of tilted magnetic field
    3. II.3 Generation of spin-current
      1. II.3.1 Derivation of spin-expectation value
      2. II.3.2 Wavepacket representation
      3. II.3.3 Optical excitation of spin-excitons
      4. II.3.4 Charge current
    4. II.4 Derivation of X-wave solution of the hyperbolic Schrodinger equation
  3. III Spherical Fermi surface model
    1. III.1 Chiral spin resonances
    2. III.2 Chiral spin waves in the presence of tilted magnetic field
      1. III.2.1 ${\bf q}\parallel\hat{x}$
      2. III.2.2 ${\bf q}\parallel\hat{y}$
      3. III.2.3 ${\bf q}\parallel\hat{z}$

In this supplementary material (SM), we provide derivation of various
equations presented in the Main Text (MT). We provide technical details of the
calculation of collective mode and its dispersion for both the Dirac-point-
vicinity (DPV) model and the spherical Fermi surface (SFS) model, as discussed
in the MT. For the DPV model (Sec. I) we use equations of motion (EoM)
approach within operator RPA (random phase approximation) to calculate
collective modes, while for the SFS model (Sec. III) we use standard
diagrammatic RPA.

##  I Dirac-point-vicinity model

We consider full Hamiltonian of the system near Dirac points of the Fermi
surface as given by Eq. (1) of the MT. The eigenvalues and eigenstates of the
single-particle Hamiltonian (Eq. (1) of the MT) are given by

$\begin{split}\varepsilon_{\bf k}^{\gamma}&=\pm
v_{F}k_{\perp}+\frac{k_{\parallel}^{2}}{2m}+\frac{\gamma}{2}\epsilon_{\bf
k};\\\ \epsilon_{\bf
k}&\equiv\varepsilon_{+}-\varepsilon_{-}=\sqrt{\Delta_{Z}^{2}+4\alpha^{2}k_{\parallel}^{2}},\\\
\ket{f_{{\bf
k}\gamma}}&=\left(\begin{array}[]{c}\frac{i\gamma\left(\gamma\Delta_{Z}+\epsilon_{\bf
k}\right)^{1/2}}{\sqrt{2\epsilon_{\bf k}}}\frac{k_{-}}{|{\bf
k}_{\parallel}|}\\\ \frac{\left(-\gamma\Delta_{Z}+\epsilon_{\bf
k}\right)^{1/2}}{\sqrt{2\epsilon_{\bf k}}}\end{array}\right),\end{split}$ (1)

where $\gamma=\pm 1$ denotes the chirality of the spin-split subbands. The
electron-electron interaction is treated short ranged:

$\hat{H}_{U}=\frac{U}{2V}\sum_{{\bf q},{\bf p},{\bf
p}^{\prime}}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{{\bf p}+{\bf
q},\sigma}c^{\dagger}_{{\bf p}^{\prime}-{\bf q},\sigma^{\prime}}c_{{\bf
p}^{\prime},\sigma^{\prime}}c_{{\bf p},\sigma}.$ (2)

For the EoM approach, it is easier to work in the chiral basis. We use the
transformation [1]

$\begin{split}a_{{\bf
k}+}^{\dagger}&=\frac{1}{\sqrt{2}}\left(\sqrt{1+\frac{\Delta_{Z}}{\epsilon_{\bf
k}}}c_{{\bf k}\uparrow}^{\dagger}+e^{i\phi_{\bf
k}}\sqrt{1-\frac{\Delta_{Z}}{\epsilon_{\bf k}}}c_{{\bf
k}\downarrow}^{\dagger}\right)\\\ a_{{\bf
k}-}^{\dagger}&=\frac{1}{\sqrt{2}}\left(\sqrt{1-\frac{\Delta_{Z}}{\epsilon_{\bf
k}}}c_{{\bf k}\uparrow}^{\dagger}-e^{i\phi_{\bf
k}}\sqrt{1+\frac{\Delta_{Z}}{\epsilon_{\bf k}}}c_{{\bf
k}\downarrow}^{\dagger}\right)\end{split}$ (3)

and rewrite the full Hamiltonian ($\hat{H}=\hat{H}_{0}+\hat{H}_{U}$) in the
chiral basis as [1, 2]:

$\begin{split}H^{ch}=\sum_{{\bf p}\gamma}\varepsilon_{{\bf
p}\gamma}a^{\dagger}_{{\bf p}\gamma}a_{{\bf p}\gamma}+\frac{U}{2V}\sum_{{\bf
q}^{\prime}{\bf p}{\bf
p}^{\prime}}\sum_{\gamma_{1}\gamma_{2}}\sum_{\gamma_{1}^{\prime}\gamma_{2}^{\prime}}\langle
f_{{\bf p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}|f_{{\bf
p}\gamma_{1}}\rangle\langle f_{{\bf p}^{\prime}-{\bf
q}^{\prime},\gamma_{2}^{\prime}}|f_{{\bf p}^{\prime}\gamma_{2}}\rangle
a^{\dagger}_{{\bf p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}a^{\dagger}_{{\bf
p}^{\prime}-{\bf q}^{\prime},\gamma_{2}^{\prime}}a_{{\bf
p}^{\prime}\gamma_{2}}a_{{\bf p}\gamma_{1}},\end{split}$ (4)

where $\varepsilon_{{\bf k}\gamma}$ and $|f_{{\bf p}s}\rangle$ are the
eigenvalue and eigenstate (1) of the single-particle Hamiltonian (see Eq. (1)
of the MT), respectively.

###  I.1 Operator RPA: an Equation of motion approach

We use EoM approach [2] to calculate the collective mode frequency and the
corresponding wavefunction near the Dirac points of the Fermi surface (FS).
The creation operator ($Q_{\bf q}^{\dagger}$) for the exciton-like mode with
momentum ${\bf q}$ is defined as

$Q_{\bf q}^{\dagger}=\sum_{{\bf k}s_{1}s_{2}}C_{{\bf k}{\bf
q}}^{s_{1}s_{2}}a_{{\bf k}+{\bf q},s_{1}}^{\dagger}a_{{\bf k}s_{2}},$ (5)

which obeys the Heisenberg equation of motion

$[H^{ch},Q_{\bf q}^{\dagger}]=\Omega_{\bf q}Q_{\bf q}^{\dagger},$ (6)

where $\Omega_{\bf q}$ is the exciton frequency. For the case of Hubbard
interaction the EoM for the exciton creation operator can be written as

$\begin{split}[H^{ch},Q_{\bf q}^{\dagger}]=&\sum_{{\bf k}s_{1}s_{2}}C_{{\bf
k}{\bf q}}^{s_{1}s_{2}}(\varepsilon_{{\bf k}+{\bf q},s_{1}}-\varepsilon_{{\bf
k}s_{2}})a_{{\bf k}+{\bf q},s_{1}}^{\dagger}a_{{\bf k}s_{2}}\\\
&+\frac{U}{V}\sum_{{\bf k}s_{1}s_{2}}C_{{\bf k}{\bf q}}^{s_{1}s_{2}}\sum_{{\bf
p}{\bf
q}^{\prime}}\sum_{\gamma_{1}\gamma_{1}^{\prime}\gamma_{2}^{\prime}}\langle
f_{{\bf p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}|f_{{\bf
p}\gamma_{1}}\rangle\Big{(}\langle f_{{\bf k}+{\bf q}-{\bf
q}^{\prime},\gamma_{2}^{\prime}}|f_{{\bf k}+{\bf q},s_{1}}\rangle a_{{\bf
p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf k}+{\bf q}-{\bf
q}^{\prime},\gamma_{2}^{\prime}}^{\dagger}a_{{\bf k}s_{2}}a_{{\bf
p}\gamma_{1}}\\\ &\hskip 241.84842pt-\langle f_{{\bf k}s_{2}}|f_{{\bf k}+{\bf
q}^{\prime},\gamma_{2}^{\prime}}\rangle a_{{\bf p}+{\bf
q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf k}+{\bf
q},s_{1}}^{\dagger}a_{{\bf k}+{\bf q}^{\prime},\gamma_{2}^{\prime}}a_{{\bf
p}\gamma_{1}}\Big{)}.\end{split}$ (7)

The next step is to do operator RPA which reduces the product of four
fermionic operators to that of two [3]. In particular, we approximate the
product of four fermionic operators with a sum of two-fermion terms such that
their averages over the ground state (which is a Slater determinant for the
unperturbed system) coincide. In particular, we use:

$a^{\dagger}b^{\dagger}cd\to H_{appr}=a^{\dagger}d\langle
b^{\dagger}c\rangle+\langle a^{\dagger}d\rangle
b^{\dagger}c-a^{\dagger}c\langle b^{\dagger}d\rangle-\langle
a^{\dagger}c\rangle b^{\dagger}d-\langle a^{\dagger}d\rangle\langle
b^{\dagger}c\rangle+\langle a^{\dagger}c\rangle\langle b^{\dagger}d\rangle,$
(8)

where $a,b,c,d$ are fermionic annihilation operators. Let us illustrate that
principle

$\langle a^{\dagger}b^{\dagger}cd\rangle=\langle a^{\dagger}d\rangle\langle
b^{\dagger}c\rangle-\langle a^{\dagger}c\rangle\langle
b^{\dagger}d\rangle=\langle H_{appr}\rangle,$ (9)

assuming that the average is taken over a Slater determinant state.

For our case, the two-operator averages take the form:

$a_{{\bf p}^{\prime}\gamma^{\prime}}^{\dagger}a_{{\bf
p}\gamma}\rightarrow\langle 0|a_{{\bf
p}^{\prime}\gamma^{\prime}}^{\dagger}a_{{\bf p}\gamma}|0\rangle=\delta_{{\bf
p}{\bf p}^{\prime}}\delta_{\gamma\gamma^{\prime}}n_{{\bf p}\gamma},$ (10)

where $n_{{\bf p}\gamma}$ is the occupation number for electronic states with
momentum ${\bf p}$ from the band $\gamma$. Using Eq. (10), the four fermionic
terms in the RHS of Eq. (7) can be split into direct (Hartree) and exchange
(Fock) terms:

$\begin{split}a_{{\bf p}+{\bf
q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf k}+{\bf q}-{\bf
q}^{\prime},\gamma_{2}^{\prime}}^{\dagger}a_{{\bf k}s_{2}}a_{{\bf
p}\gamma_{1}}=&\,\,a_{{\bf p}+{\bf
q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf p}\gamma_{1}}\delta_{{\bf
q}{\bf q}^{\prime}}\delta_{\gamma_{2}^{\prime}s_{2}}n_{{\bf k}s_{2}}+a_{{\bf
k}+{\bf q}-{\bf q}^{\prime},\gamma_{2}^{\prime}}^{\dagger}a_{{\bf
k}s_{2}}\delta_{{\bf p}+{\bf q}^{\prime},{\bf
p}}\delta_{\gamma_{1}^{\prime}\gamma_{1}}n_{{\bf p}\gamma_{1}}\\\ &-a_{{\bf
p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf k}s_{2}}\delta_{{\bf
p},{\bf k}+{\bf q}-{\bf
q}^{\prime}}\delta_{\gamma_{1}\gamma_{2}^{\prime}}n_{{\bf
p}\gamma_{1}}-a_{{\bf k}+{\bf q}-{\bf
q}^{\prime},\gamma_{2}^{\prime}}^{\dagger}a_{{\bf p}\gamma_{1}}\delta_{{\bf
p}+{\bf q}^{\prime},{\bf k}}\delta_{\gamma_{1}^{\prime}s_{2}}n_{{\bf
k}s_{2}}\\\ a_{{\bf p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf
k}+{\bf q},s_{1}}^{\dagger}a_{{\bf k}+{\bf
q}^{\prime}\gamma_{2}^{\prime}}a_{{\bf p}\gamma_{1}}=&\,\,a_{{\bf p}+{\bf
q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf p}\gamma_{1}}\delta_{{\bf
q}{\bf q}^{\prime}}\delta_{\gamma_{2}^{\prime}s_{1}}n_{{\bf k}+{\bf
q},s_{1}}+a_{{\bf k}+{\bf q},s_{1}}^{\dagger}a_{{\bf k}+{\bf
q}^{\prime}\gamma_{2}^{\prime}}\delta_{{\bf p}+{\bf q}^{\prime},{\bf
p}}\delta_{\gamma_{1}^{\prime}\gamma_{1}}n_{{\bf p}\gamma_{1}}\\\ &-a_{{\bf
p}+{\bf q}^{\prime},\gamma_{1}^{\prime}}^{\dagger}a_{{\bf k}+{\bf
q}^{\prime}\gamma_{2}^{\prime}}\delta_{{\bf p},{\bf k}+{\bf
q}}\delta_{\gamma_{1}s_{1}}n_{{\bf k}+{\bf q},s_{1}}-a_{{\bf k}+{\bf
q},s_{1}}^{\dagger}a_{{\bf p}\gamma_{1}}\delta_{{\bf p}{\bf
k}}\delta_{\gamma_{1}^{\prime}\gamma_{2}^{\prime}}n_{{\bf k}+{\bf
q}^{\prime},\gamma_{2}^{\prime}}.\end{split}$ (11)

Note that the c-number parts appearing in Eq. (8) vanish identically here due
to nonzero momentum ${\bf q}^{\prime}$. Substituting Eq. (11) into Eq. (7),
and then using Eq. (6) we get

$\begin{split}\sum_{{\bf k}s_{1}s_{2}}C_{{\bf k}{\bf
q}}^{s_{1}s_{2}}(\Omega_{\bf q}-\varepsilon_{{\bf k}+{\bf
q},s_{1}}+\varepsilon_{{\bf k}s_{2}})a_{{\bf k}+{\bf
q},s_{1}}^{\dagger}a_{{\bf k}s_{2}}&=\frac{U}{V}\sum_{{\bf
k}s_{1}s_{2}}\langle f_{{\bf k}+{\bf q},s_{1}}|f_{{\bf k}s_{2}}\rangle a_{{\bf
k}+{\bf q},s_{1}}^{\dagger}a_{{\bf k}s_{2}}\sum_{{\bf
k}^{\prime}\gamma_{1}\gamma_{2}}C_{{\bf k}^{\prime}{\bf
q}}^{\gamma_{1}\gamma_{2}}\langle f_{{\bf k}^{\prime}\gamma_{2}}|f_{{\bf
k}^{\prime}+{\bf q},\gamma_{1}}\rangle\\\ &\hskip
184.9429pt\times\left(n_{{\bf k}^{\prime}\gamma_{2}}-n_{{\bf k}^{\prime}+{\bf
q},\gamma_{1}}\right)\\\ &+\frac{U}{V}\sum_{{\bf k}s_{1}s_{2}}C_{{\bf k}{\bf
q}}^{s_{1}s_{2}}\sum_{{\bf k}^{\prime}\gamma_{1}\gamma_{2}}\langle f_{{\bf
k}^{\prime}\gamma_{1}}|f_{{\bf k}^{\prime}\gamma_{1}}\rangle n_{{\bf
k}^{\prime}\gamma_{1}}\Big{(}\langle f_{{\bf k}+{\bf q},\gamma_{2}}|f_{{\bf
k}+{\bf q},s_{1}}\rangle a_{{\bf k}+{\bf q},\gamma_{2}}^{\dagger}a_{{\bf
k}s_{2}}\\\ &\hskip 170.71652pt-\langle f_{{\bf k}s_{2}}|f_{{\bf
k}\gamma_{2}}\rangle a_{{\bf k}+{\bf q},s_{1}}^{\dagger}a_{{\bf
k}\gamma_{2}}\Big{)}\\\ &+\frac{U}{V}\sum_{{\bf k}s_{1}s_{2}}a_{{\bf k}+{\bf
q},s_{1}}^{\dagger}a_{{\bf k}s_{2}}\sum_{{\bf
k}^{\prime}\gamma_{1}\gamma_{2}}C_{{\bf k}^{\prime}{\bf
q}}^{\gamma_{1}\gamma_{2}}\langle f_{{\bf k}^{\prime}\gamma_{2}}|f_{{\bf
k}s_{2}}\rangle\langle f_{{\bf k}+{\bf q},s_{1}}|f_{{\bf k}^{\prime}+{\bf
q},\gamma_{1}}\rangle\\\ &\hskip 184.9429pt\times\left(n_{{\bf
k}^{\prime}+{\bf q},\gamma_{1}}-n_{{\bf k}^{\prime}\gamma_{2}}\right)\\\
&+\frac{U}{V}\sum_{{\bf k}s_{1}s_{2}}C_{{\bf k}{\bf q}}^{s_{1}s_{2}}\sum_{{\bf
k}^{\prime}\gamma_{1}\gamma_{2}}\Big{(}\langle f_{{\bf
k}^{\prime}\gamma_{2}}|f_{{\bf k}\gamma_{1}}\rangle\langle f_{{\bf
k}s_{2}}|f_{{\bf k}^{\prime}\gamma_{2}}\rangle a_{{\bf k}+{\bf
q},s_{1}}^{\dagger}a_{{\bf k}\gamma_{1}}\\\ &\hskip 85.35826pt-\langle f_{{\bf
k}^{\prime}\gamma_{2}}|f_{{\bf k}+{\bf q},s_{1}}\rangle\langle f_{{\bf k}+{\bf
q},\gamma_{1}}|f_{{\bf k}^{\prime},\gamma_{2}}\rangle a_{{\bf k}+{\bf
q},\gamma_{1}}^{\dagger}a_{{\bf k}s_{2}}\Big{)}n_{{\bf
k}^{\prime},\gamma_{2}}.\end{split}$ (12)

The second and fourth terms in the RHS of Eq. (12) are identified as Hartree
(tadpole diagram) and Fock (exchange diagram) self-energy corrections,
respectively. Indeed, same terms are obtained if the interacting part of the
Hamiltonian (second term of Eq. (4)) is modified using the Hartree-Fock
approximation:

$\delta H_{HF}=\frac{U}{V}\sum_{{\bf p}{\bf
p}^{\prime}}\sum_{\gamma_{1}\gamma_{2}^{\prime}\gamma_{2}}\langle f_{{\bf
p}^{\prime}\gamma_{1}}|f_{{\bf p}^{\prime}\gamma_{1}}\rangle\langle f_{{\bf
p}\gamma_{2}^{\prime}}|f_{{\bf p}\gamma_{2}}\rangle n_{{\bf
p}^{\prime}\gamma_{1}}a^{\dagger}_{{\bf p}\gamma_{2}^{\prime}}a_{{\bf
p}\gamma_{2}}-\frac{U}{V}\sum_{{\bf p}{\bf
p}^{\prime}}\sum_{\gamma_{1}\gamma_{2}^{\prime}\gamma_{2}}\langle f_{{\bf
p}^{\prime}\gamma_{2}}|f_{{\bf p}\gamma_{1}}\rangle\langle f_{{\bf
p}\gamma_{2}^{\prime}}|f_{{\bf p}^{\prime}\gamma_{2}}\rangle n_{{\bf
p}^{\prime}\gamma_{2}}a^{\dagger}_{{\bf p}\gamma_{2}^{\prime}}a_{{\bf
p}\gamma_{1}}+const.$ (13)

Indeed, evaluating $[\delta H_{HF},Q_{\bf q}^{\dagger}]$ we obtain:

$\begin{split}[\delta H_{HF},Q_{\bf q}^{\dagger}]&=\frac{U}{V}\sum_{{\bf
k}s_{1}s_{2}}C_{{\bf k}{\bf q}}^{s_{1}s_{2}}\sum_{{\bf p}{\bf
p}^{\prime}}\sum_{\gamma_{1}\gamma_{2}}\langle f_{{\bf
p}^{\prime}\gamma_{1}}|f_{{\bf p}^{\prime}\gamma_{1}}\rangle n_{{\bf
p}^{\prime}\gamma_{1}}\left(\langle f_{{\bf k}+{\bf q},\gamma_{2}}|f_{{\bf
k}+{\bf q},s_{1}}\rangle a^{\dagger}_{{\bf k}+{\bf q},\gamma_{2}}a_{{\bf
k}s_{2}}-\langle f_{{\bf k}s_{2}}|f_{{\bf k}\gamma_{2}}\rangle
a^{\dagger}_{{\bf k}+{\bf q},s_{1}}a_{{\bf k}\gamma_{2}}\right)\\\
&-\frac{U}{V}\sum_{{\bf k}s_{1}s_{2}}C_{{\bf k}{\bf q}}^{s_{1}s_{2}}\sum_{{\bf
p}{\bf p}^{\prime}}\sum_{\gamma_{1}\gamma_{2}}\Big{(}\langle f_{{\bf
p}^{\prime}\gamma_{2}}|f_{{\bf k}+{\bf q},s_{1}}\rangle\langle f_{{\bf k}+{\bf
q},\gamma_{1}}|f_{{\bf p}^{\prime}\gamma_{2}}\rangle a^{\dagger}_{{\bf k}+{\bf
q},\gamma_{1}}a_{{\bf k}s_{2}}\\\ &\hskip 113.81102pt-\langle f_{{\bf
p}^{\prime}\gamma_{2}}|f_{{\bf k}+{\bf q},\gamma_{1}}\rangle\langle f_{{\bf
k}+{\bf q},s_{2}}|f_{{\bf p}^{\prime}\gamma_{2}}\rangle a^{\dagger}_{{\bf
k}+{\bf q},s_{1}}a_{{\bf k}\gamma_{2}}\Big{)}n_{{\bf
p}^{\prime}\gamma_{2}},\end{split}$ (14)

which is equivalent to second and fourth terms of Eq. (12). The physical
meaning of these terms is a renormalization of the parameters of the single-
particle Hamiltonian ($\hat{H}_{0}$, as defined in the first term of Eq. (4))
[3]. In the main text, we assume these correction to be already included in
the parameters of $H_{0}$ and do not write them out explicitly.

The first and third terms of Eq. (12) correspond to interaction energy of
collective particle-hole excitations beyond Hartree-Fock. Importantly, the
indices $\gamma_{1},\gamma_{2}$ in these terms do not need to correspond to
$s_{1},s_{2}$, indicating that different types of modes (inter- and intra-
band) can mix. To analyze these terms, we project the operator equations (12)
to the specific components, i.e., evaluate $\langle 0|a_{{\bf k}+{\bf
q},s_{1}}\eqref{eq}a_{{\bf k}s_{2}}^{\dagger}|0\rangle$ which results in the
following equation:

$\begin{gathered}C_{{\bf k}{\bf q}}^{s_{1}s_{2}}(\Omega_{\bf
q}-\tilde{\varepsilon}_{{\bf k}+{\bf q},s_{1}}+\tilde{\varepsilon}_{{\bf
k}s_{2}})=T_{1}^{{\bf k}{\bf q},s_{1},s_{2}}+T_{3}^{{\bf k}{\bf
q},s_{1},s_{2}},\\\ T_{1}^{{\bf k}{\bf q},s_{1},s_{2}}\equiv\frac{U}{V}\langle
f_{{\bf k}+{\bf q},s_{1}}|f_{{\bf k}s_{2}}\rangle\sum_{{\bf
k}^{\prime}\gamma_{1}\gamma_{2}}C_{{\bf k}^{\prime}{\bf
q}}^{\gamma_{1}\gamma_{2}}\langle f_{{\bf k}^{\prime}\gamma_{2}}|f_{{\bf
k}^{\prime}+{\bf q},\gamma_{1}}\rangle\left(n_{{\bf
k}^{\prime}\gamma_{2}}-n_{{\bf k}^{\prime}+{\bf q},\gamma_{1}}\right),\\\
T_{3}^{{\bf k}{\bf q},s_{1},s_{2}}=\frac{U}{V}\sum_{{\bf
k}^{\prime}\gamma_{1}\gamma_{2}}C_{{\bf k}^{\prime}{\bf
q}}^{\gamma_{1}\gamma_{2}}\langle f_{{\bf k}^{\prime}\gamma_{2}}|f_{{\bf
k}s_{2}}\rangle\langle f_{{\bf k}+{\bf q},s_{1}}|f_{{\bf k}^{\prime}+{\bf
q},\gamma_{1}}\rangle\left(n_{{\bf k}^{\prime}+{\bf q},\gamma_{1}}-n_{{\bf
k}^{\prime}\gamma_{2}}\right).\end{gathered}$ (15)

Here we assumed that the single-particle energies have been renormalized by
the Hartree-Fock correction as discussed above.

###  I.2 Spin-exciton bound state

We now discuss bound state in the spin sector, which will now be called as
“spin-exciton”. For the exciton case, we set $s_{1}=+$ and $s_{2}=-$, which
refers to the transition between spin-split subbands and in the presence of
interactions ($U$) result into a bound state in the spin sector. In Eq. (15),
we then identify three types of terms: (1) $\gamma_{1}=\gamma_{2}$, (2)
$\gamma_{1}=+,\gamma_{2}=-$ and (3) $\gamma_{1}=-,\gamma_{2}=+$. Explicitly
separating these terms we have:

$\begin{split}&\hskip 113.81102pt\hat{T}_{1}^{{\bf k},+,-}+\hat{T}_{3}^{{\bf
k},+,-}=(1)+(2)+(3),\\\ (1)&:\frac{U}{V}\sum_{{\bf
k}^{\prime},\gamma=\pm}C_{{\bf k}^{\prime}{\bf q}}^{\gamma\gamma}\left(n_{{\bf
k}^{\prime}\gamma}-n_{{\bf k}^{\prime}+{\bf q},\gamma}\right)\big{[}\langle
f_{{\bf k}+{\bf q},+}|f_{{\bf k}-}\rangle\langle f_{{\bf
k}^{\prime}\gamma}|f_{{\bf k}^{\prime}+{\bf q},\gamma}\rangle-\langle f_{{\bf
k}^{\prime}\gamma}|f_{{\bf k}-}\rangle\langle f_{{\bf k}+{\bf q},+}|f_{{\bf
k}^{\prime}+{\bf q},\gamma}\rangle\big{]}\\\ (2)&:-\frac{U}{V}\sum_{{\bf
k}^{\prime}}C_{{\bf k}^{\prime}{\bf q}}^{+-}\left(n_{{\bf
k}^{\prime}-}-n_{{\bf k}^{\prime}+{\bf q},+}\right)\big{[}\langle f_{{\bf
k}^{\prime}-}|f_{{\bf k}-}\rangle\langle f_{{\bf k}+{\bf q},+}|f_{{\bf
k}^{\prime}+{\bf q},+}\rangle-\langle f_{{\bf k}+{\bf q},+}|f_{{\bf
k}-}\rangle\langle f_{{\bf k}^{\prime}-}|f_{{\bf k}^{\prime}+{\bf
q},+}\rangle\big{]}\\\ (3)&:\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}{\bf q}}^{-+}\left(n_{{\bf k}^{\prime}+{\bf q},-}-n_{{\bf
k}^{\prime}+}\right)\big{[}\langle f_{{\bf k}^{\prime}+}|f_{{\bf
k}-}\rangle\langle f_{{\bf k}+{\bf q},+}|f_{{\bf k}^{\prime}+{\bf
q},-}\rangle-\langle f_{{\bf k}+{\bf q},+}|f_{{\bf k}-}\rangle\langle f_{{\bf
k}^{\prime}+}|f_{{\bf k}^{\prime}+{\bf q},-}\rangle\big{]}\end{split}$ (16)

To simplify the results, we use the form of eigenstates (1) for
$k_{\parallel}\ll\Delta_{z}/2\alpha$ (valid near gapped Dirac points $\pm
k_{F}\hat{z}$) and corresponding relations:

$\begin{gathered}|f_{{\bf
k}+}\rangle\approx\frac{1}{\sqrt{1+\alpha^{2}k_{\parallel}^{2}/\Delta_{Z}^{2}}}\begin{pmatrix}1\\\
-i\alpha k_{+}/\Delta_{Z}\end{pmatrix},\,\,|f_{{\bf
k}-}\rangle\approx\frac{1}{\sqrt{1+\alpha^{2}k_{\parallel}^{2}/\Delta_{Z}^{2}}}\begin{pmatrix}-i\alpha
k_{-}/\Delta_{Z}\\\ 1\end{pmatrix};\\\ \langle f_{{\bf k}\pm}|f_{{\bf
k}^{\prime}\pm}\rangle\approx 1+O([k,k^{\prime}]^{2}),\,\,\langle f_{{\bf
k}\pm}|f_{{\bf
k}^{\prime}\mp}\rangle\approx\frac{i\alpha(k_{\mp}-k^{\prime}_{\mp})}{\Delta_{Z}}+O([k,k^{\prime}]^{3}).\end{gathered}$
(17)

Note that the overall phase factor, $ik_{-}/|{\bf k}_{\parallel}|$, that one
would obtain in the expression of $|f_{{\bf k}+}\rangle$ (1) in the
$k_{\parallel}\ll\Delta_{Z}/2\alpha$ limit is omitted for transparency.
Indeed, the current form of eigenstates (17) suggests two spin flip bands in
the large $\Delta_{Z}$ limit. We observe that the contribution (1) in Eq. (16)
vanishes at ${\bf q}=0$ because the Fermi functions cancel each other and the
leading order $q$-dependent correction is $\propto\alpha q/\Delta_{Z}$ at
small $q$. Moreover, the wavefunction components $C^{\gamma\gamma}_{{\bf
k}^{\prime}{\bf q}}$ correspond to intraband excitations, where the collective
modes are plasmons [3]. To treat plasmons, we would need to take into account
the long-range character of the Coulomb interaction before screening. However,
since plasmon energies are different from exciton ones and typically large for
metals, the corrections due to this coupling would scale as $U(\alpha
q/\Delta_{Z})^{2}/(\omega_{pl}-\Omega_{ex})$ and therefore can be neglected
for weak coupling or small ${\bf q}$.

Contribution (2) in Eq. (16) is the most important one, leading to the
existence of a bound state. Assuming a weakly-coupled exciton, we write Eq.
(15), neglecting (3) [to be discussed below], at leading order in the relevant
momenta compared to $\Delta_{Z}/2\alpha$:

$\begin{split}C_{{\bf k}{\bf q}}^{+-}&=-\frac{U}{V}\frac{\sum_{{\bf
k}^{\prime}}C_{{\bf k}^{\prime}{\bf q}}^{+-}\left(n_{{\bf
k}^{\prime}-}-n_{{\bf k}^{\prime}+{\bf q},+}\right)}{\Omega_{\bf
q}-\tilde{\varepsilon}_{{\bf k}+{\bf q},+}+\tilde{\varepsilon}_{{\bf k}-}}\\\
\Rightarrow 1&=-\frac{U}{V}\sum_{{\bf k},a=\pm}\frac{\left(n_{{\bf
k}-}-n_{{\bf k}+{\bf q},+}\right)}{\Omega_{\bf
q}-\tilde{\Delta}_{Z}-\frac{({\bf k}_{\parallel}+{\bf
q}_{\parallel})^{2}}{2m_{+}}-\frac{k_{\parallel}^{2}}{2m_{-}}-av_{F}q_{z}}\\\
\Rightarrow 1&=-\frac{U}{V}\sum_{{\bf k},a=\pm}\frac{\left(n_{{\bf
k}-}-n_{{\bf k}+{\bf q},+}\right)}{\Omega_{\bf
q}-\tilde{\Delta}_{Z}-\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left({\bf
k}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}-\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}-av_{F}q_{z}},\end{split}$
(18)

where $m_{\pm}^{-1}=2\alpha^{2}/\tilde{\Delta}_{Z}\pm m^{-1}$,
$\tilde{\Delta}_{Z}$ is the self-energy renormalized Zeeman energy (to be
discussed in Sec. I.3) and $a=\pm 1$ refers to two Dirac points ($\pm
k_{F}\hat{z}$) of the Fermi surface. Let us now solve the last line of Eq.
(18). First, near Dirac points we write Fermi functions as $n_{{\bf
k}-}\approx 1$ and $n_{{\bf k}+}=0$. Next, we split the ${\bf k}$-integral
into its parallel (${\bf k}_{\parallel}$) and perpendicular ($k_{\perp}$)
components and make the substitution: ${\bf k}_{\parallel}\to{\bf
k}_{\parallel}-{\bf q}_{\parallel}(m_{-}/(m_{+}+m_{-}))$. Finally, we write
$\Omega_{\bf q}\approx\tilde{\Delta}_{Z}-E_{B}({\bf q})$, where $E_{B}({\bf
q})$ is the binding energy, and obtain

$\begin{split}1&=U\int_{-\tilde{\Delta}_{Z}/2v_{F}}^{\tilde{\Delta}_{Z}/2v_{F}}\frac{dk_{\perp}}{2\pi}\int\frac{d^{2}k_{\parallel}}{(2\pi)^{2}}\left[\frac{1}{E_{B}({\bf
q})+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}k_{\parallel}^{2}+\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}+v_{F}q_{z}}+\frac{1}{E_{B}({\bf
q})+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}k_{\parallel}^{2}+\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}-v_{F}q_{z}}\right]\\\
&=\frac{U\tilde{\Delta}_{Z}}{4\pi^{2}v_{F}}\int_{0}^{\infty}dk_{\parallel}k_{\parallel}\left[\frac{1}{E_{B}^{+}({\bf
q})+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}k_{\parallel}^{2}}+\frac{1}{E_{B}^{-}({\bf
q})+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}k_{\parallel}^{2}}\right];\,\,\,\,E_{B}^{\pm}({\bf
q})\equiv E_{B}({\bf q})+\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}\pm
v_{F}q_{z}\\\
&=\frac{U\tilde{\Delta}_{Z}}{4\pi^{2}v_{F}}\frac{m_{+}m_{-}}{m_{+}+m_{-}}\text{ln}\left[\left(\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\right)^{2}\frac{\Lambda^{4}}{E_{B}^{+}({\bf
q})E_{B}^{-}({\bf q})}\right],\end{split}$ (19)

where $\Lambda\gg\sqrt{2m_{+}m_{-}E_{B}^{\pm}({\bf q})/(m_{+}+m_{-})}$ is the
high momentum cutoff. Solving the above equation, we obtain for the collective
mode frequency as $\Omega_{\bf q}\approx\tilde{\Delta}_{Z}-E_{B}({\bf q})$,
where the binding energy $E_{B}({\bf q})$ is obtained as

$\begin{split}E_{B}^{+}({\bf q})E_{B}^{-}({\bf
q})&=\Lambda^{4}\left(\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\right)^{2}\text{Exp}\left[-\frac{4\pi^{2}v_{F}}{U\Delta_{Z}}\frac{m_{+}+m_{-}}{m_{+}m_{-}}\right]\\\
\Rightarrow E_{B}({\bf
q})&=\sqrt{\Lambda^{4}\left(\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\right)^{2}\text{Exp}\left[-\frac{4\pi^{2}v_{F}}{U\Delta_{Z}}\frac{m_{+}+m_{-}}{m_{+}m_{-}}\right]+v_{F}^{2}q_{z}^{2}}-\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}.\end{split}$
(20)

Finally, using $m_{\pm}^{-1}\equiv 2\alpha^{2}/\tilde{\Delta}_{Z}\pm m^{-1}$
and upon expanding for small $q_{z}$, we can write the collective mode
frequency and its dispersion:

$\Omega_{\bf
q}\approx\tilde{\Delta}_{Z}-E_{B}^{0}+\frac{q_{\parallel}^{2}}{2m}\frac{(4m^{2}\alpha^{4}-\tilde{\Delta}_{Z}^{2})}{4m\alpha^{2}\tilde{\Delta}_{Z}}-\frac{v_{F}^{2}q_{z}^{2}}{2E_{B}^{0}};\,\,\,\,E_{B}^{0}\equiv\frac{2\alpha^{2}\Lambda^{2}}{\tilde{\Delta}_{Z}}\text{Exp}\left[-\frac{8\pi^{2}\alpha^{2}v_{F}}{U\tilde{\Delta}_{Z}^{2}}\right],$
(21)

where the form of $\tilde{\Delta}_{Z}$ is derived in Sec. I.3 (see Eq. (28))
while discussing self-energy. As we can see, the mode dispersion is hyperbolic
when $2m\alpha^{2}\gg\tilde{\Delta}_{Z}$, while it is parabolic when
$2m\alpha^{2}\ll\tilde{\Delta}_{Z}$. Indeed, for the former, the dispersion in
the plane is opposite to that perpendicular to the plane. On the other hand,
for the latter, the dispersion is parabolically downward in all directions.

We can now discuss the contribution (3) of Eq. (16) and argue that it brings a
small correction to the contribution (2). Let us do this in the limit $q=0$
(corrections can be shown to be $O(\alpha^{2}q^{2}/\Delta_{Z}^{2})$) and in
the vicinity of Dirac points, i.e., $n_{-}=1,n_{+}=0$. We get coupled
equations:

$\begin{split}\left(\Omega-\tilde{\Delta}_{Z}-\frac{2\alpha^{2}{\bf
k}_{\parallel}^{2}}{\tilde{\Delta}_{Z}}\right)C_{{\bf
k}}^{+-}&=-\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}}^{+-}+\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}}^{-+}\langle f_{{\bf k}^{\prime}+}|f_{{\bf k}-}\rangle\langle
f_{{\bf k}+}|f_{{\bf k}^{\prime}-}\rangle,\\\
\left(\Omega+\tilde{\Delta}_{Z}+\frac{2\alpha^{2}{\bf
k}_{\parallel}^{2}}{\tilde{\Delta}_{Z}}\right)C_{{\bf
k}}^{-+}&=\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}}^{-+}-\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}}^{+-}\langle f_{{\bf k}^{\prime}-}|f_{{\bf k}+}\rangle\langle
f_{{\bf k}-}|f_{{\bf k}^{\prime}+}\rangle.\end{split}$ (22)

Let us treat the new terms perturbatively. In that case
$\Omega\approx\tilde{\Delta}_{Z}>0$ and therefore $C_{{\bf k}}^{-+}$ can be
approximated by

$C_{{\bf
k}}^{-+}\approx\frac{U}{V\left(2\tilde{\Delta}_{Z}+\frac{2\alpha^{2}{\bf
k}_{\parallel}^{2}}{\tilde{\Delta}_{Z}}\right)}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}}^{+-}\langle f_{{\bf k}^{\prime}-}|f_{{\bf k}+}\rangle\langle
f_{{\bf k},-}|f_{{\bf k}^{\prime},+}\rangle+O(U^{2}).$ (23)

This results in the correction to the equation for $C^{+-}_{\bf k}$:

$\begin{gathered}\left(\Omega-\tilde{\Delta}_{Z}-\frac{2\alpha^{2}k_{\parallel}^{2}}{\tilde{\Delta}_{Z}}\right)C_{{\bf
k}}^{+-}=-\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}}^{+-}+\delta,\\\ \delta\approx\frac{U}{V}\sum_{{\bf
k}^{\prime}}\frac{U}{V}\sum_{{\bf k}^{\prime\prime}}\frac{\langle f_{{\bf
k}^{\prime}+}|f_{{\bf k}-}\rangle\langle f_{{\bf k},+}|f_{{\bf
k}^{\prime},-}\rangle\langle f_{{\bf k}^{\prime\prime}-}|f_{{\bf
k}^{\prime}+}\rangle\langle f_{{\bf k}^{\prime},-}|f_{{\bf
k}^{\prime\prime},+}\rangle}{2\tilde{\Delta}_{Z}+\frac{2\alpha^{2}k_{\parallel}^{\prime
2}}{\tilde{\Delta}_{Z}}}C_{{\bf k}^{\prime\prime}}^{+-},\\\
|\delta|<\frac{U}{V}\sum_{{\bf k}^{\prime\prime}}\frac{U}{V}|C_{{\bf
k}^{\prime\prime}}^{+-}|\sum_{{\bf
k}^{\prime}}\frac{1}{2\tilde{\Delta}_{Z}+\frac{2\alpha^{2}k_{\parallel}^{\prime
2}}{\tilde{\Delta}_{Z}}}\sim\frac{U}{V}\frac{U\tilde{\Delta}_{Z}^{2}}{\alpha^{2}v_{F}}\log\frac{\alpha\Lambda}{\tilde{\Delta}_{Z}}\sum_{{\bf
k}^{\prime\prime}}|C_{{\bf k}^{\prime\prime}}^{+-}|\ll\frac{U}{V}\sum_{{\bf
k}^{\prime\prime}}|C_{{\bf k}^{\prime\prime}}^{+-}|.\end{gathered}$ (24)

Note that the upper bound in the last line is likely an overestimate, since
the matrix elements turn to zero when ${\bf k}^{\prime}={\bf k}$ or ${\bf
k}^{\prime\prime}={\bf k}^{\prime}$. For the last inequality, we argue that
$(U\tilde{\Delta}_{Z}^{2}/\alpha^{2}v_{F})\text{ln}(\alpha\Lambda/\tilde{\Delta}_{Z})\ll
1$. Indeed, in Eq. (19) the integral suggests that
$(U\tilde{\Delta}_{Z}^{2}/\alpha^{2}v_{F})\text{ln}(\alpha\Lambda/E_{B}^{0})\sim
1$. Since $E_{B}^{0}\ll\tilde{\Delta}_{Z}$, the last inequality in Eq. (24) is
justified. This clearly demonstrates that $\delta$ is small compared to the
other term in the RHS of the first line of Eq. (24). Hence, we can neglect the
effects of coupling between $C^{+-}$ and $C^{-+}$.

###  I.3 Self-energy

In this section, we discuss self-energy which renormalizes parameters of the
single-particle Hamiltonian $\hat{H}_{0}$. We consider second and fourth terms
of Eq. (12), also illustrated in Eq. (14), for the exciton case: $s_{1}=+$ and
$s_{2}=-$. The second term vanishes trivially, so it is only the fourth term
which is relevant. We rewrite the integral equation (15) as

$\begin{split}&C_{{\bf k}{\bf q}}^{+-}\left(\Omega_{\bf q}-\varepsilon_{{\bf
k}+{\bf q},+}+\varepsilon_{{\bf k}-}-\delta\Sigma_{{\bf k}{\bf
q}}\right)=T_{1}^{{\bf k}{\bf q},+,-}+T_{3}^{{\bf k}{\bf q},+,-},\\\
\text{where}\,\,\,&\delta\Sigma_{{\bf k}{\bf q}}=\frac{U}{V}\sum_{{\bf
k}^{\prime}}\Big{[}\Big{\\{}\big{|}\langle f_{{\bf k}^{\prime}-}|f_{{\bf
k}-}\rangle\big{|}^{2}-\big{|}\langle f_{{\bf k}^{\prime}-}|f_{{\bf k}+{\bf
q},+}\rangle\big{|}^{2}\Big{\\}}n_{{\bf k}^{\prime}-}+\Big{\\{}\big{|}\langle
f_{{\bf k}^{\prime}+}|f_{{\bf k}-}\rangle\big{|}^{2}-\big{|}\langle f_{{\bf
k}^{\prime}+}|f_{{\bf k}+{\bf q},+}\rangle\big{|}^{2}\Big{\\}}n_{{\bf
k}^{\prime}+}\Big{]}\end{split}$ (25)

####  I.3.1 Zero Rashba case: Silin-Leggett mode

Let us first discuss the limiting case: $q=0$ and $\alpha=0$ (no Rashba). This
case must give Silin-Leggett mode which is nothing but the Larmor frequency
$\Omega=\Delta_{Z}$ [4]. The self-energy is simply $\delta\Sigma_{\bf
k}=(U/V)\sum_{{\bf k}^{\prime}}(n_{{\bf k}^{\prime}-}-n_{{\bf k}^{\prime}+})$.
At $\alpha=0$, the Fermi surface is spherical so no Dirac point approximation
can be used here. Using $n_{{\bf k}\pm}\equiv
n_{F}(k^{2}/2m\pm\Delta_{Z}/2-\mu)$, where $n_{F}$ is the Fermi function and
$\mu$ is the chemical potential, and then doing ${\bf k}$-integral assuming
$\Delta_{Z}\ll\mu$, we get $\delta\Sigma_{\bf k}\approx U\nu_{F}\Delta_{Z}/2$,
where $\nu_{F}\equiv mk_{F}/\pi^{2}$ is the total 3D density of states. Hence,
the self-energy correction leads to the renormalization of Zeeman energy:
$\Delta_{Z}\to\tilde{\Delta}_{Z}\equiv\Delta_{Z}(1+U\nu_{F}/2)$.

We can now solve for the collective mode frequency (18) and obtain:
$\Omega=\tilde{\Delta}_{Z}(1-U\nu_{F}/2)=\Delta_{Z}(1-U^{2}\nu_{F}^{2}/4)$.
The Kohn’s theorem states that the Silin-Leggett mode is a bare Larmor
frequency, protected from renormalization due to interaction. This is indeed
the case here and the quadratic $U$ correction is an artifact of not
calculating the self-energy correction self-consistently. Indeed, the exact
self-consistent calculation of self-energy suggests
$\tilde{\Delta}_{Z}=\Delta_{Z}/(1-U\nu_{F}/2)$ [5]. The factor in the
denominator is cancelled exactly by the one appearing in the collective mode
frequency. So, the exact cancellation of self-energy and vertex corrections
leads to the collective mode as Larmor frequency, $\Omega=\Delta_{Z}$. Hence,
in our expression we simply ignore the $U^{2}$ correction.

####  I.3.2 Finite Rashba case

Let us now discuss self-energy correction for finite Rashba $(\alpha\neq 0)$.
We will still assume $2\alpha k_{\parallel}\ll\Delta_{Z}$, which is valid near
Dirac-points on the Fermi surface. Moreover, to be consistent with an earlier
work [1], we recover the additional phase factor, $ik_{-}/|{\bf
k}_{\parallel}|$, in the limiting form of $|f_{{\bf k}+}\rangle$ near Dirac
point. The eigenstate (17) then modifies to

$\begin{gathered}|f_{{\bf k}+}\rangle\approx\frac{ik_{-}/|{\bf
k}_{\parallel}|}{\sqrt{1+\alpha^{2}k_{\parallel}^{2}/\Delta_{Z}^{2}}}\begin{pmatrix}1\\\
-i\alpha k_{+}/\Delta_{Z}\end{pmatrix},\,\,|f_{{\bf
k}-}\rangle\approx\frac{1}{\sqrt{1+\alpha^{2}k_{\parallel}^{2}/\Delta_{Z}^{2}}}\begin{pmatrix}-i\alpha
k_{-}/\Delta_{Z}\\\ 1\end{pmatrix};\end{gathered}$ (26)

Using this, we obtain the form of $\delta\Sigma_{{\bf k}{\bf q}}$ (25) as

$\begin{split}\delta\Sigma_{{\bf k}{\bf
q}}&=\frac{U\Delta_{Z}^{2}}{2V}\left(\frac{1}{\epsilon_{\bf
k}}+\frac{1}{\epsilon_{{\bf k}+{\bf q}}}\right)\sum_{{\bf
k}^{\prime}}\frac{n_{{\bf k}^{\prime}-}-n_{{\bf k}^{\prime}+}}{\epsilon_{{\bf
k}^{\prime}}}\\\ &=\frac{U\Delta_{Z}^{2}}{2}\left(\frac{1}{\epsilon_{\bf
k}}+\frac{1}{\epsilon_{{\bf k}+{\bf
q}}}\right)\frac{\Delta_{Z}^{2}}{8\pi^{2}\alpha^{2}v_{F}}\text{ln}\frac{2\alpha^{2}\Lambda^{2}}{\Delta_{Z}^{2}},\end{split}$
(27)

where $\Lambda\gg\Delta_{Z}/\sqrt{2}\alpha$ is the same high-momentum cut off
as obtained in the bound state energy (21) above. Also, the Fermi function,
$n_{{\bf k}^{\prime}\pm}$, appearing above has been assumed to take the form
near Dirac points: $n_{{\bf k}^{\prime}-}=1$ and $n_{{\bf k}^{\prime}+}=0$.
Combining self-energy (27) with the single-particle dispersion as in
$(\Omega_{\bf q}-\varepsilon_{{\bf k}+{\bf q},+}+\varepsilon_{{\bf
k}-}-\delta\Sigma_{{\bf k}{\bf q}})$, we find that only $\Delta_{Z}$
renormalizes due to interaction:

$\Delta_{Z}\to\tilde{\Delta}_{Z}\approx\Delta_{Z}\left(1+\frac{U\Delta_{Z}^{2}}{8\pi^{2}\alpha^{2}v_{F}}\text{ln}\frac{2\alpha^{2}\Lambda^{2}}{\Delta_{Z}^{2}}\right)$
(28)

From the same logic as presented while arguing $C_{\bf k}^{-+}\ll C_{\bf
k}^{+-}$ (see text below Eq. (24)), we can say that the second term of Eq.
(28) is small compared to 1.

###  I.4 Wavefunction of the collective mode frequency

In this section we will provide technical details of the calculation of
wavefunction (more specifically $C_{{\bf k}{\bf q}}^{+-}$) of the spin
excitonic wave (21). Let us rewrite Eq. (15) for the exciton case (for brevity
we simply write $\tilde{\Delta}_{Z}$ as $\Delta_{Z}$):

$\begin{split}C_{{\bf k}{\bf q}}^{+-}(\Omega_{\bf q}-\varepsilon_{{\bf k}+{\bf
q},+}+\varepsilon_{{\bf k}-})=\frac{U}{V}\sum_{{\bf k}^{\prime}}C_{{\bf
k}^{\prime}{\bf q}}^{+-}\left(n_{{\bf k}^{\prime}+{\bf q},+}-n_{{\bf
k}^{\prime}-}\right)\big{[}\langle f_{{\bf k}^{\prime}-}|f_{{\bf
k}-}\rangle\langle f_{{\bf k}+{\bf q},+}|f_{{\bf k}^{\prime}+{\bf
q},+}\rangle-\langle f_{{\bf k}+{\bf q},+}|f_{{\bf k}-}\rangle\langle f_{{\bf
k}^{\prime}-}|f_{{\bf k}^{\prime}+{\bf q},+}\rangle\big{]}\end{split}$ (29)

Using Eq. (26), the matrix element in Eq. (29) within large-$\Delta_{Z}$
approximation can be written as

$\begin{split}\langle f_{{\bf k}^{\prime}-}|f_{{\bf k}-}\rangle\langle f_{{\bf
k}+{\bf q},+}|f_{{\bf k}^{\prime}+{\bf
q},+}\rangle&\approx\frac{(k_{+}+q_{+})}{|({\bf k}+{\bf
q})_{\parallel}|}\frac{(k^{\prime}_{-}+q_{-})}{|({\bf k}^{\prime}+{\bf
q})_{\parallel}|},\\\ \langle f_{{\bf k}+{\bf q},+}|f_{{\bf k}-}\rangle\langle
f_{{\bf k}^{\prime}-}|f_{{\bf k}^{\prime}+{\bf
q},+}\rangle&\approx\frac{\alpha^{2}q_{+}q_{-}}{\Delta_{Z}^{2}}\frac{(k_{+}+q_{+})}{|({\bf
k}+{\bf q})_{\parallel}|}\frac{(k^{\prime}_{-}+q_{-})}{|({\bf k}^{\prime}+{\bf
q})_{\parallel}|}.\end{split}$ (30)

Using Eq. (30), the Eq. (29) simplifies to

$C_{{\bf k}{\bf
q}}^{+-}=U\left(1-\frac{\alpha^{2}q_{\parallel}^{2}}{\Delta_{Z}^{2}}\right)\frac{(k_{+}+q_{+})}{|({\bf
k}+{\bf q})_{\parallel}|}\frac{N_{\bf q}}{\Omega_{\bf q}-\varepsilon_{{\bf
k}+{\bf q},+}+\varepsilon_{{\bf k}-}},$ (31)

where

$N_{\bf q}\equiv\frac{1}{V}\sum_{{\bf k}^{\prime}}C_{{\bf k}^{\prime}{\bf
q}}^{+-}\frac{k^{\prime}_{-}+q_{-}}{|({\bf k}^{\prime}+{\bf
q})_{\parallel}|}(n_{{\bf k}^{\prime}+{\bf q},+}-n_{{\bf k}^{\prime}-}).$ (32)

Note that an extra $q_{\parallel}^{2}$ correction term shows up in Eq. (31) on
the RHS, which is a result of interband matrix element at finite $q$ (second
term of Eq. (29)). This $q_{\parallel}^{2}$ correction should have been
considered in the dispersion of spin-exciton frequency (21) as well. However,
we argue that this correction is very small, so ignored in Eq. (21). Indeed,
if we notice, this factor merely dresses the interaction: $U\to
U(1-\alpha^{2}q_{\parallel}^{2}/\Delta_{Z}^{2})$ (31). If we plug this back in
Eq. (21), it leads to a positive correction in $q_{\parallel}^{2}$ which is
good for hyperbolicity of the dispersion. However, the correction is of order
$\mathcal{O}\left((E_{B}^{0}/\Delta_{Z})(\alpha^{2}v_{F}/U\Delta_{Z}^{2})\right)$,
which is small since $E_{B}^{0}\ll\Delta_{Z}$. So it was ignored.

We now calculate the normalization factor $N_{\bf q}$. As discussed in the MT,
$N_{\bf q}$ can be obtained by using the condition, $\langle 0|[Q_{\bf
q},Q_{{\bf q}^{\prime}}^{\dagger}]|0\rangle=\delta_{{\bf q}{\bf q}^{\prime}}$,
where $Q_{\bf q}^{\dagger}$ (5) is the general exciton creation operator.
Using the commutation relation

$\begin{split}[Q_{\bf q},Q_{{\bf q}^{\prime}}^{\dagger}]=\sum_{{\bf
k}s_{1}s_{2}}\bigg{\\{}\sum_{s_{2}^{\prime}}a_{{\bf k}s_{2}}^{\dagger}a_{{\bf
k}+{\bf q}-{\bf q}^{\prime},s_{2}^{\prime}}(C_{{\bf k}{\bf
q}}^{s_{1}s_{2}})^{\dagger}(C_{{\bf k}+{\bf q}-{\bf q}^{\prime},{\bf
q}^{\prime}}^{s_{1}s_{2}^{\prime}})-\sum_{s_{1}^{\prime}}a_{{\bf k}+{\bf
q}^{\prime}s_{1}^{\prime}}^{\dagger}a_{{\bf k}+{\bf q},s_{1}}(C_{{\bf k}{\bf
q}}^{s_{1}s_{2}})^{\dagger}(C_{{\bf k}{\bf
q}^{\prime}}^{s_{1}^{\prime}s_{2}})\bigg{\\}}.\end{split}$ (33)

and then averaging over the ground state, we get the relation $\sum_{{\bf
k}s_{1}s_{2}}|C_{{\bf k}{\bf q}}^{s_{1}s_{2}}|^{2}(n_{{\bf k}s_{2}}-n_{{\bf
k}+{\bf q}s_{1}})=1$. Choosing $\\{s_{1},s_{2}\\}\equiv\\{+,-\\}$ and then
using Eq. (31) in this, we get

$|N_{\bf
q}|^{2}=\frac{1}{U^{2}\left(1-\alpha^{2}q_{\parallel}^{2}/\Delta_{Z}^{2}\right)^{2}}\left[\sum_{\bf
k}\frac{n_{{\bf k}-}-n_{{\bf k}+{\bf q},+}}{(\Omega_{\bf q}-\varepsilon_{{\bf
k}+{\bf q},+}+\varepsilon_{{\bf k}-})^{2}}\right]^{-1},$ (34)

where we recall the form of $\varepsilon_{{\bf p}\gamma}$ in Eq. (1). The
wavefunction finally is given by

$|C_{{\bf k}{\bf q}}^{+-}|^{2}=\frac{1}{(\Omega_{\bf q}-\varepsilon_{{\bf
k}+{\bf q},+}+\varepsilon_{{\bf k}-})^{2}}\left[\sum_{{\bf
k}^{\prime}}\frac{n_{{\bf k}^{\prime}-}-n_{{\bf k}^{\prime}+{\bf
q},+}}{(\Omega_{\bf q}-\varepsilon_{{\bf k}^{\prime}+{\bf
q},+}+\varepsilon_{{\bf k}^{\prime}-})^{2}}\right]^{-1}.$ (35)

Now it is left to do ${\bf k}^{\prime}$ integral in the above equation (35).

To begin the discussion, we first assume ${\bf q}\parallel\hat{x},\hat{y}$.
The ${\bf k}$-integral in Eq. (35) can be done the same way as discussed in
Eq. (19). Writing $\Omega_{\bf
q}\approx\Delta_{Z}-\sqrt{(E_{B}^{0})^{2}+v_{F}^{2}q_{z}^{2}}+q_{\parallel}^{2}/2(m_{+}+m_{-})$
(using Eqs. (18) - (21)), where $m_{\pm}^{-1}=2\alpha^{2}/\Delta_{Z}\pm
m^{-1}$, we calculate:

$\begin{split}\sum_{\bf k}\frac{n_{{\bf k}-}-n_{{\bf k}+{\bf
q}_{\parallel},+}}{(\Omega_{{\bf q}_{\parallel}}-\varepsilon_{{\bf k}+{\bf
q}_{\parallel},+}+\varepsilon_{{\bf k}-})^{2}}&=\sum_{\bf k}\frac{n_{{\bf
k}-}-n_{{\bf k}+{\bf
q}_{\parallel},+}}{\left(\Delta_{Z}-\sqrt{(E_{B}^{0})^{2}+v_{F}^{2}q_{z}^{2}}+\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}-\varepsilon_{{\bf
k}+{\bf q}_{\parallel},+}+\varepsilon_{{\bf k}-}\right)^{2}}\\\
&\approx\frac{\Delta_{Z}V}{2\pi
v_{F}}\sum_{a=\pm}\int\frac{d^{2}k_{\parallel}}{(2\pi)^{2}}\frac{1}{\left(\sqrt{(E_{B}^{0})^{2}+v_{F}^{2}q_{z}^{2}}+av_{F}q_{z}+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left({\bf
k}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)\right)^{2}}\\\
&\approx\frac{\Delta_{Z}^{2}V}{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right).\\\
\text{Eq.~{}\eqref{norm}}\Rightarrow\,\,|N_{\bf
q}|^{2}&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{VU^{2}\Delta_{Z}^{2}}\left(1-\frac{\alpha^{2}q_{\parallel}^{2}}{\Delta_{Z}^{2}}\right)^{-1}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}.\end{split}$
(36)

Substituting (36) into (35), we get

$\begin{split}|C_{a{\bf k}{\bf
q}_{\parallel}}^{+-}|^{2}&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}V}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\frac{1}{\left(\Omega_{\bf
q}-\Delta_{Z}-av_{F}q_{z}-\frac{({\bf k}_{\parallel}+{\bf
q}_{\parallel})^{2}}{2m_{+}}-\frac{k_{\parallel}^{2}}{2m_{-}}\right)^{2}}\\\
&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}V}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\left(\sqrt{(E_{B}^{0})^{2}+v_{F}^{2}q_{z}^{2}}+av_{F}q_{z}+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left({\bf
k}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}\right)^{-2}.\end{split}$
(37)

As a conclusion, the form of the coefficient of the wavefunction near two
Dirac points can be written as

$\begin{split}|C_{(+k_{Fz},{\bf k}_{\parallel}),{\bf
q}}^{+-}|^{2}&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}V}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\left(E_{B}^{0}+\frac{v_{F}^{2}q_{z}^{2}}{2E_{B}^{0}}+v_{F}q_{z}+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left({\bf
k}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}\right)^{-2},\\\
|C_{(-k_{Fz},{\bf k}_{\parallel}),{\bf
q}}^{+-}|^{2}&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}V}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\left(E_{B}^{0}+\frac{v_{F}^{2}q_{z}^{2}}{2E_{B}^{0}}-v_{F}q_{z}+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left({\bf
k}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}\right)^{-2},\\\
\end{split}$ (38)

where $m_{\pm}^{-1}\equiv 2\alpha^{2}/\Delta_{Z}\pm m^{-1}$ and $E_{B}^{0}$ is
given in Eq. (21).

##  II Generation of spin-current by in-plane component of the tilted magnetic
field for Dirac-point-vicinity model

We now demonstrate that a finite spin-current can be generated by optically
exciting the exciton-like collective mode in the presence of a tilted magnetic
field. In this section, we will show that the in-plane component of the field
is crucial for the generation of spin-current. The single-particle Hamiltonian
can be written as

$\hat{H}=\pm
v_{F}k_{\perp}\hat{\sigma}_{0}+\frac{k_{\parallel}^{2}}{2m}\hat{\sigma}_{0}+\frac{\Delta_{Z}}{2}\hat{\sigma}_{z}+\frac{\Delta_{y}}{2}\hat{\sigma}_{y}+\alpha(k_{y}\hat{\sigma}_{x}-k_{x}\hat{\sigma}_{y})$
(39)

where $\Delta_{Z}=g\mu_{B}B\cos\theta_{\bf B}$ and
$\Delta_{y}=g\mu_{B}B\sin\theta_{\bf B}$, with $\theta_{\bf B}$ as the polar
angle of the magnetic field ${\bf B}$. We assume large weight of the field to
be along $\hat{z}$-direction, while a small weight along the plane, say
$\hat{y}$, such that $\Delta_{y}\ll\Delta_{Z}$. The eigenvalues and
eigenstates of the Hamiltonian (39) are given by

$\begin{split}\varepsilon_{\bf k}^{s}&=\pm
v_{F}k_{z}+\frac{k_{x}^{2}+k_{y}^{2}}{2m_{b}}+\frac{s}{2}\epsilon_{\bar{{\bf
k}}},~{}\epsilon_{\bar{{\bf
k}}}=\sqrt{4\alpha^{2}(\bar{k}_{x}^{2}+k_{y}^{2})+\Delta_{Z}^{2}};\\\
\ket{f_{{\bf
k}s}}&=\left(\begin{array}[]{c}\frac{is\big{(}s\Delta_{Z}+\epsilon_{\bar{{\bf
k}}}\big{)}^{1/2}}{\sqrt{2\epsilon_{\bar{{\bf
k}}}}}\frac{(\bar{k}_{x}-ik_{y})}{\sqrt{\bar{k}_{x}^{2}+k_{y}^{2}}}\\\
\frac{\big{(}-s\Delta_{Z}+\epsilon_{\bar{{\bf
k}}}\big{)}^{1/2}}{\sqrt{2\epsilon_{\bar{{\bf
k}}}}}\end{array}\right),\,\,\bar{k}_{x}=k_{x}-\frac{\Delta_{y}}{2\alpha},\end{split}$
(40)

where $s=\pm 1$ denotes the chirality of the spin-split subbands. To discuss
the spin-current, we first need to know how the dispersion of the collective
mode and its wavefunction modify due to the presence of in-plane magnetic
field. In order to do so, we repeat the same analysis as done in Secs. I.2 and
I.4, but now also in the presence of in-plane magnetic field.

###  II.1 Dispersion of the collective mode frequency in the presence of
tilted magnetic field

Near the Dirac point, the matrix element of the integral equation (15) for the
exciton case ($s_{1}=+$ and $s_{2}=-$) decouples again and reduces to the same
form as in Eq. (18), except that $k_{x}$ only in the particle-hole excitation
energy ($\epsilon_{\bf k}$) is replaced by $\bar{k}_{x}$, where $\bar{k}_{x}$
is defined in Eq. (40). Accounting this change ($k_{x}\to\bar{k}_{x}$), we
obtain the integral equation as

$1=-U\sum_{a=\pm}\int_{-\Delta_{Z}/2v_{F}}^{\Delta_{Z}/2v_{F}}\frac{dk_{\perp}}{2\pi}\int\frac{d^{2}k_{\parallel}}{(2\pi)^{2}}\frac{1}{\Omega_{\bf
q}-av_{F}q_{z}-\frac{q_{\parallel}^{2}+2{\bf k}_{\parallel}\cdot{\bf
q}_{\parallel}}{2m}-\frac{1}{2}\epsilon_{\bar{{\bf
k}}_{\parallel}}-\frac{1}{2}\epsilon_{\bar{{\bf k}}_{\parallel}+{\bf
q}_{\parallel}}},$ (41)

where

$\epsilon_{\bar{{\bf k}}_{\parallel}+{\bf
q}_{\parallel}}\approx\Delta_{Z}+\frac{4\alpha^{2}\left((\bar{k}_{x}+q_{x})^{2}+(k_{y}+q_{y})^{2}\right)}{2\Delta_{Z}}$
(42)

is the modified finite-$q$ particle-hole excitation energy in the presence of
an in-plane magnetic field ($\Delta_{y}$), with
$\bar{k}_{\parallel}=(\bar{k}_{x}^{2}+k_{y}^{2})^{1/2}$ and $\bar{k}_{x}$
given by Eq. (40). For simplicity, we do variable substitution,
$k_{x}\to\bar{k}_{x}+\Delta_{y}/2\alpha$, and transform the integral in Eq.
(41) on $k_{x}$ to that on $\bar{k}_{x}$ ($k_{y}$ remains the same):

$\begin{split}1=-U\sum_{a=\pm}\int_{-\Delta_{Z}/2v_{F}}^{\Delta_{Z}/2v_{F}}\frac{dk_{\perp}}{2\pi}\int\frac{d^{2}\bar{k}_{\parallel}}{(2\pi)^{2}}\frac{1}{\Omega_{\bf
q}-\Delta_{Z}-av_{F}q_{z}-\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left(\bar{{\bf
k}}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}-\frac{q_{\parallel}^{2}}{2(m_{+}+m_{-})}-\frac{\Delta_{y}q_{x}}{2m\alpha}}\\\
\end{split}$ (43)

As we can see, the in-plane magnetic field brings an extra last term in the
denominator of the above integrand, as compared to that in the absence of it
(18). We calculate the above integral (43) in the same way as done in Eq.
(19), and then solve for $\Omega_{\bf q}$:

$\Omega_{\bf
q}\approx\Delta_{Z}-\frac{2\alpha^{2}\Lambda^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{8\pi^{2}\alpha^{2}v_{F}}{U\Delta_{Z}^{2}}\right]-\frac{v_{F}^{2}q_{z}^{2}}{\frac{4\alpha^{2}\Lambda^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{8\pi^{2}\alpha^{2}v_{F}}{U\Delta_{Z}^{2}}\right]}+\frac{\Delta_{y}}{2m\alpha}q_{x}+\frac{(q_{x}^{2}+q_{y}^{2})}{2m}\frac{(4m^{2}\alpha^{4}-\Delta_{Z}^{2})}{4m\alpha^{2}\Delta_{Z}}.$
(44)

The in-plane magnetic field leads to a linear in $q_{x}$ term in the
dispersion, which tells us that at $q=0$ the exciton will move with velocity
$\Delta_{y}/2m\alpha$ in the $\hat{x}$-direction. This is useful in the
generation of spin-current which is the subject of discussion in Sec. II.3.

###  II.2 Wavefunction of the collective mode frequency in the presence of
tilted magnetic field

We now discuss the wavefunction. As we saw, the dispersion along $\hat{y}$ and
$\hat{z}$ remains unaffected by $\Delta_{y}$ (44). Hence, the $q_{y}$\- and
$q_{z}$-dependent parts of the wavefunction also will be unaffected by
$\Delta_{y}$ and takes the same form as in Eq. (45). However, the
$k_{\parallel}$ everywhere is replaced by $\bar{k}_{\parallel}$, which is
where the dependence of $\Delta_{y}$ shows up in the wavefunction. Indeed, the
explicit calculation, along the same lines as discussed in Secs. I.4 and II.1,
suggests

$\begin{split}|C_{(+k_{Fz},{\bf k}_{\parallel}),{\bf
q}}^{+-}|^{2}&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}V}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\left(E_{B}^{0}+\frac{v_{F}^{2}q_{z}^{2}}{2E_{B}^{0}}+v_{F}q_{z}+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left(\bar{{\bf
k}}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}\right)^{-2},\\\
|C_{(-k_{Fz},{\bf k}_{\parallel}),{\bf
q}}^{+-}|^{2}&\approx\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}V}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\left(E_{B}^{0}+\frac{v_{F}^{2}q_{z}^{2}}{2E_{B}^{0}}-v_{F}q_{z}+\frac{m_{+}+m_{-}}{2m_{+}m_{-}}\left(\bar{{\bf
k}}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}\right)^{-2},\\\
\end{split}$ (45)

where $|\bar{{\bf
k}}_{\parallel}|\equiv\bar{k}_{\parallel}=(\bar{k}_{x}^{2}+k_{y}^{2})^{1/2}$,
with $\bar{k}_{x}\equiv k_{x}-\Delta_{y}/2\alpha$ given in Eq. (40),
$m_{\pm}^{-1}\equiv 2\alpha^{2}/\Delta_{Z}\pm m^{-1}$ and $E_{B}^{0}$ is given
in Eq. (21). It can be checked that the wavefunction obeys the normalization
condition: $\sum_{\bf k}|C_{{\bf k}{\bf q}}^{+-}|^{2}=1$, or more explicitly
$\sum_{{\bf k},a=\pm}|C_{(ak_{Fz},{\bf k}_{\parallel}),{\bf q}}^{+-}|^{2}=1$.

###  II.3 Generation of spin-current

In this section, we will discuss the generation of spin-current by in-plane
magnetic field. More specifically, we will first calculate spin-expectation
value in the exciton basis in Sec. II.3.1 and find that the exciton carries a
spin almost 1. Next, we will demonstrate the wavepacket representation of the
exciton wavefunction in Sec. II.3.2. Finally, we will discuss the optical
excitation of spin-excitons in Sec. II.3.3 which is needed to provide an
estimate of the exciton generation rate. We will conclude this section by
discussing that these excitons carry pure spin-current and any accompanying
charge current is zero (Sec. II.3.4).

####  II.3.1 Derivation of spin-expectation value

In the second quantized notation, the total spin operator is given by

$\hat{\textbf{S}}=\sum_{{\bf k}\alpha\alpha^{\prime}}a_{{\bf
k}\alpha}^{\dagger}\textbf{S}_{\alpha\alpha^{\prime}}a_{{\bf
k}\alpha^{\prime}},\,\,\,\,\text{where}\,\,\,\textbf{S}_{\alpha\alpha^{\prime}}=\frac{1}{2}\langle
f_{{\bf k}\alpha}|\hat{\boldsymbol{\sigma}}|f_{{\bf
k}\alpha^{\prime}}\rangle,$ (46)

where $|f_{{\bf k}\alpha}\rangle$ is the single-particle eigenstate given by
Eq. (1) and $\hat{\boldsymbol{\sigma}}$ is a Pauli spin matrices. Using the
form of spin-exciton state as

$|xc({\bf q})\rangle\equiv C_{{\bf k}{\bf q}}^{+-}a_{{\bf k}+{\bf
q},+}^{\dagger}a_{{\bf k}-}|FG\rangle,$ (47)

the expectation value of the spin operator can be written as

$\begin{split}\langle xc({\bf q})|\hat{\textbf{S}}|xc({\bf
q})\rangle=\frac{1}{2}\sum_{{\bf k}_{1}{\bf k}_{2}{\bf
k}_{3}}\sum_{\gamma_{3}\gamma_{4}}\langle FG|\left(C_{{\bf k}_{1}{\bf
q}}^{+-}\right)^{*}a_{{\bf k}_{1}-}^{\dagger}a_{{\bf k}_{1}+{\bf q},+}&a_{{\bf
k}_{2}\gamma_{3}}^{\dagger}\langle f_{{\bf
k}_{2}\gamma_{3}}|\hat{\boldsymbol{\sigma}}|f_{{\bf k}_{2}\gamma_{4}}\rangle
a_{{\bf k}_{2}\gamma_{4}}\left(C_{{\bf k}_{3}{\bf q}}^{+-}\right)a_{{\bf
k}_{3}+{\bf q},+}^{\dagger}a_{{\bf k}_{3}-}|FG\rangle,\end{split}$ (48)

As we can see, there are total six fermion operators. We now adopt operator
RPA which amounts to using fermion anticommutation relation and collect only
terms which are reduced to two fermion operators [2]. Then we use the
property: $\langle 0|a_{{\bf p}\alpha}^{\dagger}a_{{\bf
p}\alpha^{\prime}}|0\rangle=n_{{\bf p}\alpha}\delta_{\alpha\alpha^{\prime}}$,
where $n_{{\bf k}\alpha}$ is the equilibrium Fermi function. Within RPA, we
finally get [2]

$\langle xc({\bf q})|\hat{\textbf{S}}|xc({\bf q})\rangle=\frac{1}{2}\sum_{\bf
k}|C_{{\bf k}_{1}{\bf q}}^{+-}|^{2}\Big{(}\langle f_{{\bf k}_{1}+{\bf
q},+}|\hat{\boldsymbol{\sigma}}|f_{{\bf k}_{1}+{\bf q},+}\rangle-\langle
f_{{\bf k}_{1}-}|\hat{\boldsymbol{\sigma}}|f_{{\bf
k}_{1}-}\rangle\Big{)}\left(n_{{\bf k}_{1}-}-n_{{\bf k}_{1}+{\bf
q},+}\right).$ (49)

To derive Eq. (49) we select creation-annihilation pairs in Eq. (48) and using
the fermion anticommutation relation contract the number of fermion operators.
We start by selecting the pair $a_{{\bf k}_{2}\gamma_{4}}a_{{\bf k}_{3}+{\bf
q},+}^{\dagger}$ and write this as $\left(\delta_{{\bf k}_{2},{\bf k}_{3}+{\bf
q}}\delta_{\gamma_{4},+}-a_{{\bf k}_{3}+{\bf q},+}^{\dagger}a_{{\bf
k}_{2}\gamma_{4}}\right)$. In the reduced four fermion operators, we select
another such pair to reduce the same to two fermions. The averaging of two
fermions then gives the first term of Eq. (49). In the process of fermion
reduction using its anticommutation relation, we will have two more terms with
four fermions and six fermions. The second term of Eq. (49) comes from
reducing the six fermions operator term to two fermions operator in the same
way as it is done for the first term. So, finally we arrive at Eq. (49) after
averaging over the ground state Fermi gas $|FG\rangle$. In the contracting
process, we have multiple other terms with four fermion and six fermion
operators. All these terms vanish upon averaging over the ground state.

Equation (49) is the main equation as all of our results for the spin-
expectation value follows from this equation only. Using the forms of
$|C_{{\bf k}{\bf q}}^{+-}|^{2}$ (45) and $|f_{{\bf p}\pm}\rangle$ (40), the
spin-expectation value in the exciton basis (49) up to a leading order
correction in $q$ can be calculated:

$\begin{split}\langle\hat{\textbf{S}}_{z}\rangle&=\frac{1}{2}\sum_{\bf
k}|C_{{\bf k}{\bf q}}^{+-}|^{2}\left(\langle f_{{\bf k}+{\bf
q},+}|\hat{\sigma}_{z}|f_{{\bf k}+{\bf q},+}\rangle-\langle f_{{\bf
k}-}|\hat{\sigma}_{z}|f_{{\bf k}-}\rangle\right)\\\
&=\frac{1}{2}\sum_{\bar{{\bf k}}}|C_{\bar{{\bf k}}{\bf
q}}^{+-}|^{2}\left(\frac{\Delta_{Z}}{\epsilon_{\bar{{\bf k}}+{\bf
q}_{\parallel}}}+\frac{\Delta_{Z}}{\epsilon_{\bar{{\bf
k}}}}\right)\approx\sum_{\bar{{\bf k}}}|C_{\bar{{\bf k}}{\bf
q}}^{+-}|^{2}\left(1-\frac{\alpha^{2}}{\Delta_{Z}^{2}}\left((\bar{{\bf
k}}_{\parallel}+{\bf q}_{\parallel})^{2}+\bar{{\bf
k}}_{\parallel}^{2}\right)\right)\\\ &\approx
1-\frac{\alpha^{2}(q_{x}^{2}+q_{y}^{2})}{2\Delta_{Z}^{2}}\left(1+\frac{\Delta_{Z}^{2}}{4m^{2}\alpha^{4}}\right)+\frac{v_{F}^{2}q_{z}^{2}}{2\Delta_{Z}E_{B}^{0}}\text{log}\frac{2\alpha^{2}\Lambda^{2}}{\Delta_{Z}E_{B}^{0}}+\mathcal{O}(q^{4}),\\\
\langle\hat{\textbf{S}}_{x}\rangle&=\frac{1}{2}\sum_{\bf k}|C_{{\bf k}{\bf
q}}^{+-}|^{2}\left(\langle f_{{\bf k}+{\bf q},+}|\hat{\sigma}_{x}|f_{{\bf
k}+{\bf q},+}\rangle-\langle f_{{\bf k}-}|\hat{\sigma}_{x}|f_{{\bf
k}-}\rangle\right)\\\ &=\frac{1}{2}\sum_{\bar{{\bf k}}}|C_{\bar{{\bf k}}{\bf
q}}^{+-}|^{2}\left(\frac{2\alpha(k_{y}+q_{y})}{\epsilon_{\bar{{\bf k}}+{\bf
q}_{\parallel}}}+\frac{2\alpha k_{y}}{\epsilon_{\bar{{\bf k}}}}\right)\\\
&\approx-\frac{q_{y}}{2m\alpha}+\mathcal{O}(q_{\parallel}^{3}),\\\
\langle\hat{\textbf{S}}_{y}\rangle&=\frac{1}{2}\sum_{\bf k}|C_{{\bf k}{\bf
q}}^{+-}|^{2}\left(\langle f_{{\bf k}+{\bf q},+}|\hat{\sigma}_{y}|f_{{\bf
k}+{\bf q},+}\rangle-\langle f_{{\bf k}-}|\hat{\sigma}_{y}|f_{{\bf
k}-}\rangle\right)\\\ &=-\frac{1}{2}\sum_{\bar{{\bf k}}}|C_{\bar{{\bf k}}{\bf
q}}^{+-}|^{2}\left(\frac{2\alpha(\bar{k}_{x}+q_{x})}{\epsilon_{\bar{{\bf
k}}+{\bf q}_{\parallel}}}+\frac{2\alpha\bar{k}_{x}}{\epsilon_{\bar{{\bf
k}}}}\right)\\\
&\approx\frac{q_{x}}{2m\alpha}+\mathcal{O}(q_{\parallel}^{3}).\end{split}$
(50)

As we can see, at $q=0$, $\langle\hat{\textbf{S}}_{z}\rangle\approx 1$, while
$\langle\hat{\textbf{S}}_{x,y}\rangle$ vanishes. There is no $q_{z}$-dependent
correction in $\langle\hat{S}_{x}\rangle$ and $\langle\hat{S}_{y}\rangle$
because the ${\bf k}$-integral vanishes when ${\bf q}\parallel\hat{z}$.
Indeed, since $\epsilon_{\bar{{\bf k}}+{\bf q}_{\parallel}}$ does not depend
on $q_{z}$, a simple observation suggests that the ${\bf k}$-integral is odd.

####  II.3.2 Wavepacket representation

It is required to calculate $\langle xc({\bf q})|\hat{\textbf{S}}({\bf
r})|xc({\bf q}^{\prime})\rangle$. Using Eq. (47) for the exciton wavefunction,
and the extension of Eq. (48), we can write

$\begin{split}\langle xc({\bf q})|\hat{\textbf{S}}({\bf r})|xc({\bf
q}^{\prime})\rangle=\frac{1}{2}\sum_{{\bf k}_{1}{\bf k}_{3}}\sum_{{\bf
k}_{2}{\bf k}_{2}^{\prime}}\sum_{\gamma_{3}\gamma_{4}}e^{-i({\bf k}_{2}-{\bf
k}_{2}^{\prime})\cdot{\bf r}}\langle FG|\left(C_{{\bf k}_{1}{\bf
q}}^{+-}\right)^{*}a_{{\bf k}_{1}-}^{\dagger}a_{{\bf k}_{1}+{\bf q},+}a_{{\bf
k}_{2}\gamma_{3}}^{\dagger}&\langle f_{{\bf
k}_{2}\gamma_{3}}|\hat{\boldsymbol{\sigma}}|f_{{\bf
k}_{2}^{\prime}\gamma_{4}}\rangle a_{{\bf k}_{2}^{\prime}\gamma_{4}}\\\
\times&\left(C_{{\bf k}_{3}{\bf q}^{\prime}}^{+-}\right)a_{{\bf k}_{3}+{\bf
q}^{\prime},+}^{\dagger}a_{{\bf k}_{3}-}|FG\rangle.\end{split}$ (51)

Doing operator RPA the same way as in Eq. (48) to obtain Eq. (49), we obtain

$\begin{split}\langle xc({\bf q})|\hat{\textbf{S}}({\bf r})|xc({\bf
q}^{\prime})\rangle=\frac{1}{2}e^{-i({\bf q}-{\bf q}^{\prime})\cdot{\bf
r}}\sum_{{\bf k}_{1}}\left(C_{{\bf k}_{1}{\bf
q}}^{+-}\right)^{*}\bigg{(}&\langle f_{{\bf k}_{1}+{\bf
q},+}|\hat{\boldsymbol{\sigma}}|f_{{\bf k}_{1}+{\bf
q}^{\prime},+}\rangle\left(C_{{\bf k}_{1}{\bf q}^{\prime}}^{+-}\right)\\\
&-\langle f_{{\bf k}_{1}+{\bf q}-{\bf
q}^{\prime},-}|\hat{\boldsymbol{\sigma}}|f_{{\bf k}_{1}-}\rangle\left(C_{{\bf
k}_{1}+{\bf q}-{\bf q}^{\prime},{\bf
q}^{\prime}}^{+-}\right)\bigg{)}\big{(}n_{{\bf k}_{1}-}-n_{{\bf k}_{1}+{\bf
q},+}\big{)},\end{split}$ (52)

which reduces to

$\langle xc({\bf q})|\hat{\textbf{S}}({\bf r})|xc({\bf
q}^{\prime})\rangle=\frac{1}{2}e^{-i({\bf q}-{\bf q}^{\prime})\cdot{\bf
r}}\sum_{{\bf k}_{1}}\left(C_{{\bf k}_{1}{\bf
q}}^{+-}\right)^{*}\bigg{(}\langle f_{{\bf k}_{1}+{\bf
q},+}|\hat{\boldsymbol{\sigma}}|f_{{\bf k}_{1}+{\bf
q}^{\prime},+}\rangle\left(C_{{\bf k}_{1}{\bf q}^{\prime}}^{+-}\right)-\langle
f_{{\bf k}_{1}+{\bf q}-{\bf q}^{\prime},-}|\hat{\boldsymbol{\sigma}}|f_{{\bf
k}_{1}-}\rangle\left(C_{{\bf k}_{1}+{\bf q}-{\bf q}^{\prime},{\bf
q}^{\prime}}^{+-}\right)\bigg{)},$ (53)

near the Dirac point where we considered empty conduction ($n_{{\bf k}+}=0$)
and filled valence ($n_{{\bf k}-}=1$) bands. Assuming small $q$ and
$q^{\prime}$, we get Eq. (11) of the MT.

####  II.3.3 Optical excitation of spin-excitons

In this section, we will provide technical details of the optical excitation
of spin-excitons. The excitons form inside the spin gap near Dirac points of
the Fermi surface. To excite them, one shines the light with frequency in
resonance with the optical gap. According to the Fermi’s golden rule, the
transition probability in the long wavelength regime ($q\to 0$) of this
excitation is given by $|\langle xc|\hat{\textbf{j}}\cdot{\bf
A}|FG\rangle|^{2}$, where $\hat{\textbf{j}}$ is current operator proportional
to the velocity operator ($\hat{{\bf v}}$)

$\begin{split}\hat{{\bf
v}}=\left(\frac{k_{x}}{m}-\alpha\hat{\sigma}_{y},\frac{k_{y}}{m}+\alpha\hat{\sigma}_{x},\pm
v_{F}\right),\end{split}$ (54)

and ${\bf A}$ is the vector potential. Since the main purpose of this step is
to know the number of excitons excited per unit time, we assume for simplicity
that $\Delta_{y}=0$. We calculate

$\begin{split}\langle xc|\hat{\textbf{j}}\cdot{\bf A}|FG\rangle&=\sum_{{\bf
k}{\bf k}^{\prime}}\sum_{\alpha\alpha^{\prime}}\left(C_{\bf
k}^{+-}\right)^{*}\langle FG|a_{{\bf k}-}^{\dagger}a_{{\bf k}+}a_{{\bf
k}^{\prime}\alpha}^{\dagger}\langle f_{{\bf
k}^{\prime}\alpha}|\hat{\textbf{j}}\cdot{\bf A}|f_{{\bf
k}^{\prime}\alpha^{\prime}}\rangle a_{{\bf
k}^{\prime}\alpha^{\prime}}|FG\rangle\\\ &=\sum_{{\bf k}\alpha}\left(C_{\bf
k}^{+-}\right)^{*}\langle FG|\left(\langle f_{{\bf
k}+}|\hat{\textbf{j}}\cdot{\bf A}|f_{{\bf k}\alpha}\rangle a_{{\bf
k}-}^{\dagger}a_{{\bf k}\alpha}-\langle f_{{\bf
k}\alpha}|\hat{\textbf{j}}\cdot{\bf A}|f_{{\bf k}-}\rangle a_{{\bf
k}\alpha}^{\dagger}a_{{\bf k}+}\right)|FG\rangle\\\ &=\sum_{\bf k}\left(C_{\bf
k}^{+-}\right)^{*}\langle f_{{\bf k}+}|\hat{\textbf{j}}\cdot{\bf A}|f_{{\bf
k}-}\rangle\left(n_{{\bf k}-}-n_{{\bf k}+}\right)\\\ &=eU\alpha N^{*}\sum_{\bf
k}\left\\{\frac{k_{-}(\Delta_{Z}k_{x}+i\epsilon_{\bf k}k_{y})}{\epsilon_{\bf
k}(\Omega_{0}-\epsilon_{\bf
k})(k_{x}^{2}+k_{y}^{2})}A_{x}+\frac{k_{-}(\Delta_{Z}k_{y}-i\epsilon_{\bf
k}k_{x})}{\epsilon_{\bf k}(\Omega_{0}-\epsilon_{\bf
k})(k_{x}^{2}+k_{y}^{2})}A_{y}\right\\}\left(n_{{\bf k}-}-n_{{\bf
k}+}\right)\\\ &=\frac{eU\alpha
N^{*}\Delta_{Z}V}{4\pi^{2}v_{F}}(A_{x}-iA_{y})\int_{0}^{\infty}dkk\frac{\Delta_{Z}+\epsilon_{\bf
k}}{\epsilon_{\bf k}(\Omega_{0}-\epsilon_{\bf k})}\\\ &\approx\frac{eU\alpha
N^{*}\Delta_{Z}V}{4\pi^{2}v_{F}}(A_{x}-iA_{y})\int_{0}^{\infty}dkk\frac{2\Delta_{Z}}{\Delta_{Z}(-E_{B}^{0}-2\alpha^{2}k^{2}/\Delta_{Z})}\\\
&\approx-\frac{eUN^{*}\Delta_{Z}^{2}V}{8\pi^{2}\alpha
v_{F}}(A_{x}-iA_{y})\text{log}\frac{2\alpha^{2}\Lambda^{2}}{E_{B}^{0}\Delta_{Z}},\,\,\,\,\text{where}\,\,\Lambda\gg\sqrt{\frac{E_{B}^{0}\Delta_{Z}}{2\alpha^{2}}}\end{split}$
(55)

where $V$ is volume and we used Eq. (31) to write $C_{\bf k}^{+-}$, with
$N^{*}\equiv N_{\bf q}^{*}(q=0)$ (34) as the normalization factor. Also, in
the above equation $\Lambda$ is the high-momentum cut off,
$\Omega_{0}\approx\Delta_{Z}-E_{B}^{0}$ (with $E_{B}^{0}\ll\Delta_{Z}$) is the
exciton frequency at $q=0$ and $\epsilon_{\bf k}$ is the particle-hole
excitation energy given by Eq. (42) at $q=0$ and $\Delta_{y}=0$ near Dirac
points ($\pm k_{F}\hat{z}$) of the Fermi surface. In the second last line of
Eq. (55), at leading order, the integrand is simplified as
$(\Delta_{Z}+\epsilon_{\bf k})/\epsilon_{\bf k}(\Omega_{0}-\epsilon_{\bf
k})\approx-2/(E_{B}^{0}+2\alpha^{2}k^{2}/\Delta_{Z})$ according to the
approximate form of $\epsilon_{\bf k}$ (42). Finally, evaluating $|\langle
xc|\hat{\textbf{j}}\cdot{\bf A}|FG\rangle|^{2}$, with
$|N|^{2}=8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}/VU^{2}\Delta_{Z}^{2}$ (36), we get

$\begin{split}\Gamma=\frac{2\pi}{\hbar}|\langle xc|\hat{\textbf{j}}\cdot{\bf
A}|FG\rangle|^{2}\delta(\Omega-\Omega_{0})=\frac{e^{2}\Delta_{Z}^{2}E_{B}^{0}V}{4\pi\hbar
v_{F}}\left(\text{log}\frac{2\alpha^{2}\Lambda^{2}}{E_{B}^{0}\Delta_{Z}}\right)^{2}\frac{E_{x}^{2}+E_{y}^{2}}{\Omega^{2}}\delta(\Omega-\Omega_{0}),\end{split}$
(56)

where $E_{x,y}=i\Omega A_{x,y}$ is the electric field component of
electromagnetic (EM) wave. To estimate the number of excitons per unit volume,
we approximate delta-function by Lorentzian assuming finite lifetime
($\Delta^{-1}$) of excitons:
$2\pi\delta(\Omega-\Omega_{0})\approx[2\Delta/((\Omega-\Omega_{0})^{2}+\Delta^{2})]_{\Delta\to
0}$.

We need to know the form of power ($P$) of the EM wave in terms of electric
field. The energy flux per unit area per unit time is given by ${\bf S}=({\bf
E}\times{\bf B})/\mu_{0}$, where $\mu_{0}$ is the permeability of free space.
The direction of ${\bf S}$ is given by the direction of propagation of the EM
wave. In free space, $B=E/c$, where $c$ is speed of light, we have
$S=E^{2}/c\mu_{0}=c\varepsilon_{0}E^{2}$, where we used
$\varepsilon_{0}\mu_{0}=1/c^{2}$. Upon time averaging, we get $\langle
S\rangle=c\varepsilon_{0}E^{2}/2$. Since $S$ is power per unit area, we can
write $P={\bf S}\cdot{\bf A}$, where $A$ is area. The average power is,
therefore, $\langle P\rangle=c\varepsilon_{0}AE_{\text{max}}^{2}/2$, or
$E_{\text{max}}^{2}=2\langle P\rangle/c\varepsilon_{0}A$. Upon recovering
correct factors of $\hbar$, the optical excitation rate per unit volume
($\Gamma/V$) near the resonance ($\Omega\approx\Omega_{0}\approx\Delta_{Z}$)
can be written as

$\frac{\Gamma}{V}=\frac{2e^{2}E_{B}^{0}}{4\pi^{2}\varepsilon_{0}c\hbar^{2}v_{F}\Delta}\frac{P}{A}\left(\text{log}\frac{2\alpha^{2}\Lambda^{2}}{E_{B}^{0}\Delta_{Z}}\right)^{2},$
(57)

where $E_{B}^{0}$ is the binding energy at $q=0$ given by Eq. (21) and $P/A$
is power per unit area. Here $\Lambda$ is the high-momentum cutoff introduced
earlier in the definition of binding energy $E_{B}^{0}$ (21).

####  II.3.4 Charge current

The charge current is defined as

$I_{c}^{v_{i}}\propto\langle\hat{v}_{j}\rangle=\langle xc({\bf
q})|\hat{v}_{i}|xc({\bf q})\rangle=\sum_{\bf k}|C_{{\bf k}{\bf
q}}^{+-}|^{2}\left(\langle f_{{\bf k}+{\bf q},+}|\hat{v}_{i}|f_{{\bf k}+{\bf
q},+}\rangle-\langle f_{{\bf k}-}|\hat{v}_{i}|f_{{\bf k}-}\rangle\right),$
(58)

where $\hat{{\bf v}}$ is the velocity operator (54). As we know, the state
$|f_{{\bf k}\alpha}\rangle$ is a normalized eigenstate: $\langle f_{{\bf
k}\alpha}|f_{{\bf k}\alpha}\rangle=1$. Moreover, near Dirac points since the
velocity along $\hat{z}$ direction is constant, $\hat{v}_{z}=\pm v_{F}$, the
average of $\hat{v}_{z}$ for two bands cancels each other, and there is no
charge current along $\hat{z}$. So, $I_{c}^{v_{z}}=0$.

Let us now discuss charge current in $\hat{x}$ and $\hat{y}$ directions. Using
Eqs. (40), (45) and (54), we can write the expression for charge current along
$\hat{x}$ direction as

$\begin{split}I_{c}^{v_{x}}&\propto\sum_{\bf k}|C_{{\bf k}{\bf
q}}^{+-}|^{2}\left[\langle f_{{\bf k}+{\bf
q},+}|\left(\frac{k_{x}+q_{x}}{m}\hat{\sigma}_{0}-\alpha\hat{\sigma}_{y}\right)|f_{{\bf
k}+{\bf q},+}\rangle-\langle f_{{\bf
k}-}|\left(\frac{k_{x}}{m}\hat{\sigma}_{0}-\alpha\hat{\sigma}_{y}\right)|f_{{\bf
k}-}\rangle\right]\\\ &=\frac{q_{x}}{m}+\alpha\sum_{\bf k}|C_{{\bf k}{\bf
q}}^{+-}|^{2}\Big{[}-\langle f_{{\bf k}+{\bf q},+}|\hat{\sigma}_{y}|f_{{\bf
k}+{\bf q},+}\rangle+\langle f_{{\bf k}-}|\hat{\sigma}_{y}|f_{{\bf
k}-}\rangle\Big{]}\,\,\,\,(\text{using normalization condition $\sum_{\bf
k}|C_{{\bf k}{\bf q}}^{+-}|^{2}=1$})\\\ &=\frac{q_{x}}{m}+\alpha\sum_{\bf
k}|C_{{\bf k}{\bf
q}}^{+-}|^{2}\left[\frac{2\alpha(\bar{k}_{x}+q_{x})}{\epsilon_{\bar{{\bf
k}}+{\bf q}}}+\frac{2\alpha\bar{k}_{x}}{\epsilon_{\bar{{\bf k}}}}\right]\\\
&\approx\frac{q_{x}}{m}+\frac{2\alpha^{2}}{\Delta_{Z}}\sum_{\bf k}|C_{{\bf
k}{\bf q}}^{+-}|^{2}\left(2\bar{k}_{x}+q_{x}\right)\,\,\,\,(\text{using
$\epsilon_{\bf k}\approx\Delta_{Z}$ \eqref{ph_in-plane} at leading order})\\\
&=\frac{q_{x}}{m_{+}}+\frac{8\pi^{2}\alpha^{2}v_{F}E_{B}^{0}}{\Delta_{Z}^{2}}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}\frac{4\alpha^{2}}{\Delta_{Z}}\sum_{{\bf
k},a=\pm}\frac{\bar{k}_{x}}{\left\\{E_{B}^{0}+\frac{v_{F}^{2}q_{z}^{2}}{2E_{B}^{0}}+av_{F}q_{z}+\frac{2m_{+}m_{-}}{m_{+}+m_{-}}\left(\bar{{\bf
k}}_{\parallel}+{\bf
q}_{\parallel}\frac{m_{-}}{m_{+}+m_{-}}\right)^{2}\right\\}^{2}}\\\ &\text{In
the second term, the small-$q$ expansion leads to the finite contribution only
from $\bar{k}_{x}q_{x}$ term}\\\ &\Rightarrow
I_{c}^{v_{x}}\propto\frac{q_{x}}{m_{+}}-\frac{q_{x}}{m_{+}}\left(1+\frac{v_{F}^{2}q_{z}^{2}}{2(E_{B}^{0})^{2}}\right)^{-1}=0.\end{split}$
(59)

In the last line of the above equation, the correction term is an artifact of
small-$q$ expansion which results into beyond-quadratic-dispersion
approximation term. Hence it is safely ignored. The calculation along same
lines for the charge current along $\hat{y}$-direction, it can be verified
that $I_{c}^{v_{y}}=0$ as well.

###  II.4 Derivation of X-wave solution of the hyperbolic Schrodinger equation

In this section, we will provide technical details of the derivation of
envelope $X$ wave solution (see Eq. (13) of the MT) of the hyperbolic
Schrodinger equation (HSE) given by Eq. (11) of the MT. In Eq. (11) of the MT,
we make the transformation, $x,y,z\to
x^{\prime}\sqrt{p_{1}},y^{\prime}\sqrt{p_{1}},z^{\prime}\sqrt{p_{2}}$, and
write the HSE as

$\left(i\frac{\partial}{\partial t}-\Omega_{0}\right)\Psi({\bf
r}^{\prime},t)=\left[-\left(\frac{\partial^{2}}{\partial x^{\prime
2}}+\frac{\partial^{2}}{\partial y^{\prime
2}}\right)+\frac{\partial^{2}}{\partial z^{\prime
2}}-\frac{ip_{3}}{\sqrt{p_{1}}}\frac{\partial}{\partial
x^{\prime}}\right]\Psi({\bf r}^{\prime},t).$ (60)

We use the following ansatz in the above equation: $\Psi({\bf
r}^{\prime},t)=e^{-i\Omega_{0}t}e^{-ip_{3}x^{\prime}/2\sqrt{p_{1}}}e^{ip_{3}^{2}t/4p_{1}}\phi(x^{\prime},y^{\prime},z^{\prime},t)$.
The HSE (60) reduces to

$i\frac{\partial\phi}{\partial t}=-\left(\frac{\partial^{2}}{\partial
x^{\prime 2}}+\frac{\partial^{2}}{\partial y^{\prime
2}}\right)\phi+\frac{\partial^{2}\phi}{\partial z^{\prime 2}}.$ (61)

To solve this differential equation, we use
$\phi(x^{\prime},y^{\prime},z^{\prime},t)=\varphi(\xi^{\prime},y^{\prime},z^{\prime})e^{iv^{\prime}x^{\prime}/2}e^{-itv^{\prime
2}/4}$, where $\xi^{\prime}\equiv x^{\prime}-v^{\prime}t$, with
$v^{\prime}\equiv v/\sqrt{p_{1}}$ as the wave propagating velocity, and obtain

$\left[\frac{\partial^{2}}{\partial\xi^{\prime
2}}+\frac{\partial^{2}}{\partial y^{\prime 2}}-\frac{\partial^{2}}{\partial
z^{\prime 2}}\right]\varphi(\xi^{\prime},y^{\prime},z^{\prime})=0.$ (62)

The non-spreading solution of this equation has been discussed in detail in
Ref. [6], but with the propagating velocity ($v$) assumed to be along
$\hat{z}$. Extending the analysis to our case, we can write the general
solution of Eq. (62):

$\varphi(\xi^{\prime},y^{\prime},z^{\prime})=\int_{0}^{\infty}d\alpha
f(\alpha)J_{0}(\alpha r_{\perp}^{\prime})e^{i\alpha
z^{\prime}},\,\,\,\text{where}\,\,\,r_{\perp}^{\prime}\equiv(\xi^{\prime
2}+y^{2})^{1/2}.$ (63)

Finally, the full solution of Eq. (60) can be written as

$\Psi({\bf
r}^{\prime},t)=e^{-i\Omega_{0}t}e^{-ip_{3}x^{\prime}/2\sqrt{p_{1}}}e^{ip_{3}^{2}t/4p_{1}}e^{iv^{\prime}x^{\prime}/2}e^{-itv^{\prime
2}/4}\int_{0}^{\infty}d\alpha f(\alpha)J_{0}(\alpha
r_{\perp}^{\prime})e^{i\alpha z^{\prime}},$ (64)

which, in the original cartesian basis $(x,y,z,t)$ can be written as

$\begin{split}\Psi({\bf
r},t)=e^{-ip_{3}(x-p_{3}t/2)/2p_{1}}e^{iv(x-vt/2)/2p_{1}}e^{-i\Omega_{0}t}\int_{0}^{\infty}d\alpha
f(\alpha)J_{0}(\alpha r_{\parallel}/\sqrt{p_{1}})e^{i\alpha
z/\sqrt{p_{2}}},\end{split}$ (65)

where $r_{\parallel}\equiv\big{(}(x-vt)^{2}+y^{2}\big{)}^{1/2}$. As we can
see, in the above solution (65) the velocity $v$ is completely arbitrary. For
an exponentially decaying spectrum $f(\alpha)=e^{-\alpha\gamma}$ [6, 7], where
$\gamma$ is some arbitrary width, if we assume
$v=p_{3}\equiv\Delta_{y}/2m\alpha$, we obtain a fundamental linear $X$ wave
solution:

$\Psi_{X}({\bf
r},t)=\sqrt{p_{1}p_{2}}e^{-i\Omega_{0}t}\left[p_{2}(x-p_{3}t)^{2}+p_{2}y^{2}+p_{1}(\sqrt{p_{2}}\gamma-
iz)^{2}\right]^{-1/2},$ (66)

as given by Eq. (11) of the MT. This is a shape preserving envelope $X$ wave
solution of the HSE (60) which propagates with velocity
$v_{xc}=p_{3}=\Delta_{y}/2m\alpha$ in the $\hat{x}$ direction. We note that
the above general $X$-wave solution (65) is obtained for the linear HSE. The
general $X$ wave solution for the non-linear HSE has been discussed in Refs.
[6, 8, 9, 10, 11, 12, 7], which is beyond the scope of our work.

##  III Spherical Fermi surface model

In the diagrammatic RPA approach we consider tilted magnetic field such that
its components along one of the planar directions, say $\hat{y}$, in addition
to the $\hat{z}$-direction as considered in the previous section (Sec. I), is
also finite. The single-particle Hamiltonian for the system can be written as

$\hat{H}=\frac{k^{2}}{2m}\hat{\sigma}_{0}+\frac{\Delta_{Z}}{2}\hat{\sigma}_{z}+\frac{\Delta_{y}}{2}\hat{\sigma}_{y}+\alpha(k_{y}\hat{\sigma}_{x}-k_{x}\hat{\sigma}_{y}),$
(67)

where $m$ is the band mass, $\sigma_{i}$ is the Pauli matrix in the spin
space, $k$ is the electron momenta, and $\Delta_{Z}=g\mu_{B}B\cos\theta_{B}$
and $\Delta_{y}=g\mu_{B}B\sin\theta_{B}\sin\phi_{B}$, with $g$ as the electron
Lande factor, $\mu_{B}$ as the Bohr magneton, $B$ as the magnitude of the
magnetic field and $\theta_{B}(\phi_{B})$ as the polar (azimuthal) angle of
magnetic field, are Zeeman energy splittings due to magnetic field components
along $z$ and $y$ directions, respectively. Note that the polar $(\theta_{B})$
and azimuthal $(\phi_{B})$ angles of magnetic field are different from those
of electronic momenta which will be denoted as $\theta_{\bf k}$ and $\phi_{\bf
k}$, respectively. To avoid further confusion, we will use $\Delta_{Z}$ and
$\Delta_{y}$ for perpendicular and in-plane magnetic fields, correspondingly,
from now onwards. We will treat the electron-electron interaction at short
ranges (2).

The eigenvalues and eigenstates of the single-particle Hamiltonian (67) are
given by

$\begin{split}\varepsilon_{\bf
k}^{s}&=\frac{k^{2}}{2m}+\frac{s}{2}\epsilon_{\bar{{\bf
k}}},~{}\epsilon_{\bar{{\bf
k}}}=\sqrt{4\alpha^{2}(\bar{k}_{x}^{2}+k_{y}^{2})+\Delta_{Z}^{2}};\\\
\ket{s}&=\left(\begin{array}[]{c}\frac{is\big{(}s\Delta_{Z}+\epsilon_{\bar{{\bf
k}}}\big{)}^{1/2}}{\sqrt{2\epsilon_{\bar{{\bf
k}}}}}\frac{(\bar{k}_{x}-ik_{y})}{\sqrt{\bar{k}_{x}^{2}+k_{y}^{2}}}\\\
\frac{\big{(}-s\Delta_{Z}+\epsilon_{\bar{{\bf
k}}}\big{)}^{1/2}}{\sqrt{2\epsilon_{\bar{{\bf
k}}}}}\end{array}\right),\bar{k}_{x}=k_{x}-\frac{\Delta_{y}}{2\alpha},\end{split}$
(68)

where $s=\pm 1$ denotes the chirality of the spin-split subbands. Notice that
the eigenstate for the spherical Fermi surface (SFS) model (68) is same as
that for the DPV model (40).

The collective modes in the spin sector of an interacting system show up as
poles of the full spin susceptibility defined as

$\chi_{ij}(i\Omega_{n},{\bf q})=-\Pi_{ij}^{U}(i\Omega_{n},{\bf q}),$ (69)

where

$\Pi_{ij}^{U}(Q)=\int_{K}\text{Tr}\big{[}\hat{\sigma}_{i}\hat{G}^{0}(K)\hat{\gamma}_{j}\hat{G}^{0}(K+Q)\big{]}$
(70)

is the full spin polarization bubble, with $Q=(i\Omega_{n},{\bf q})$,
$\int_{K}\equiv T\sum_{\omega_{m}}\int d^{3}k/(2\pi)^{3}$, $\epsilon_{m}$ as
the Matsubara frequency, $\hat{\gamma}_{i}$ as the $2\times 2$ spin-
interaction vertex,

$\begin{split}\hat{G}^{0}(i\epsilon_{m},&{\bf
k})=\sum_{s}\frac{\ket{s}\bra{s}}{i\epsilon_{m}+\mu-\varepsilon_{\bf
k}^{s}}=\sum_{s}\hat{D}_{s}({\bf k})g(i\epsilon_{m},{\bf k}),\\\
\hat{D}_{s}({\bf
k})&=\frac{1}{2}\bigg{[}\hat{\sigma}_{0}+s\frac{\Delta_{Z}}{\epsilon_{\bar{{\bf
k}}}}\hat{\sigma}_{z}+s\frac{2\alpha}{\epsilon_{\bar{{\bf
k}}}}(\hat{\sigma}_{x}k_{y}-\hat{\sigma}_{y}\bar{k}_{x})\bigg{]},\\\
g_{s}(i\epsilon_{m},&{\bf k})=\frac{1}{i\epsilon_{m}+\mu-\varepsilon_{\bf
k}^{s}}\end{split}$ (71)

as the single-particle Green’s function and $\mu$ as the chemical potential.
Here, we term $g_{s}(i\epsilon_{m},{\bf k})$ as a chiral Green’s function. The
eigenvalues ($\varepsilon_{\bf k}^{s}$), eigenstates ($\ket{s}$) and
$\epsilon_{\bf k}$ appearing in the definition of Green’s function (71) are
given in Eq. (68). To obtain the full spin susceptibility, we perform a random
phase approximation (RPA) which sums up polarization bubbles containing ladder
series of vertex corrections. For a doped 3D semiconductor, the interaction
line inside the bubble is supposed to be a screened Coulomb potential.
However, a closed form of the full polarization bubble (70) for a screened
Coulomb potential cannot be obtained analytically; therefore, we adopt an
approximation in which the interaction is replaced by a momentum-independent
constant (short-range Hubbard $U$) [13]. Within this assumption, the spin-
interaction vertex $\hat{\gamma}_{i}$ can be written as

$\hat{\gamma}_{j}=\sigma_{j}-U\int_{K}\hat{G}^{0}(K)\hat{\gamma}_{j}\hat{G}(K+Q).$
(72)

We expand $\hat{\gamma}_{j}$ over a completet set of Pauli matrices

$\hat{\gamma}_{j}=M_{j}^{a}\hat{\sigma}_{a},$ (73)

where $a\in 0,1,2,3$ and the coefficient $M_{j}^{a}$ forms a $4\times 4$
matrix. Substituting Eq. (73) back into Eq. (72), one can find

$\begin{split}\hat{\Pi}^{U}&=\hat{\Pi}^{0}\hat{M},\\\
\hat{M}&=\bigg{(}\hat{I}+\frac{U}{2}\hat{\Pi}^{0}\bigg{)}^{-1},\end{split}$
(74)

where $\hat{\Pi}_{0}$ is the non-interacting bubble and $\hat{I}$ is the
$4\times 4$ identity matrix. Collective modes are given by poles of
$\Pi_{ij}^{U}$ which are solutions of $\text{Det}(\hat{M}^{-1})=0$.

###  III.1 Chiral spin resonances

We will first discuss $q=0$ spin collective modes for the case when
$\Delta_{y}=0$. Later we will include small $\Delta_{y}$ such that
$\Delta_{y}\ll\Delta_{Z},\Delta_{R}\ll\mu$, where $\Delta_{R}=2\alpha k_{F}$
is the Rashba energy. We will see later that there are two Rashba energy
scales in the system, $\Delta_{R}$ and $m\alpha^{2}$ with
$m\alpha^{2}\ll\Delta_{R}$, and $\Delta_{y}$ will be assumed to be small
compared to both the energy scales. The finite $q$ case will be discussed
assuming $q^{2}/2m\ll v_{F}q\ll\Delta_{y}$, the main goal of which is to track
the dispersion of a particular chiral spin mode which is robust and remains
exponentially close to the lower edge of the spin-flip particle-hole
continuum.

As mentioned in Eq. (74), to know collective modes, we require the form of all
non-zero components of the non-interacting spin polarization bubble (assuming
spherical FS and $q=0$ and $\Delta_{y}=0$):

$\begin{split}\Pi_{xx}(i\Omega_{n})&=-\frac{\nu_{F}}{4}\Bigg{(}2+\frac{(-\Omega_{n}^{2}+\Delta_{Z}^{2})}{\Delta_{R}\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}}\text{Log}\Bigg{[}\frac{\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}+\Delta_{R}}{\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}-\Delta_{R}}\Bigg{]}\Bigg{)},\\\
\Pi_{yy}(i\Omega_{n})&=\Pi_{xx}(i\Omega_{n}),\\\
\Pi_{xy}(i\Omega_{n})&=\frac{\nu_{F}}{4}\frac{2\Omega_{n}\Delta_{Z}}{\Delta_{R}\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}}\text{Log}\Bigg{[}\frac{\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}+\Delta_{R}}{\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}-\Delta_{R}}\Bigg{]},\\\
\Pi_{yx}(i\Omega_{n})&=-\Pi_{xy}(i\Omega_{n}),\\\
\Pi_{zz}(i\Omega_{n})&=-\frac{\nu_{F}}{4}\Bigg{(}4-\frac{2(\Omega_{n}^{2}+\Delta_{Z}^{2})}{\Delta_{R}\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}}\text{Log}\Bigg{[}\frac{\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}+\Delta_{R}}{\sqrt{\Omega_{n}^{2}+\Delta_{R}^{2}+\Delta_{Z}^{2}}-\Delta_{R}}\Bigg{]}\Bigg{)},\end{split}$
(75)

where $\nu_{F}\equiv mk_{F}/\pi^{2}$ is the 3D density of states assuming
spherical FS. While obtaining above bubbles (75), we assumed the spin
splitting energy ($\epsilon_{\bar{{\bf k}}}$) to be small compared to the
Fermi energy such that the Fermi function, $n_{F}(\xi_{\bf
k}+s/2\epsilon_{\bar{{\bf k}}})$, which appears after frequency summation in
the form of the bubble (actually a difference of Fermi functions,
corresponding to transition between two spin-split subbands, shows up:
$n_{F}(\xi_{\bf k}+\epsilon_{\bar{{\bf k}}}/2)-n_{F}(\xi_{\bf
k}-\epsilon_{\bar{{\bf k}}}/2)$) can be expanded as $n_{F}(\xi_{\bf
k}+(s/2)\epsilon_{\bar{{\bf k}}})\approx n_{F}(\xi_{\bf
k})+(s/2)\epsilon_{\bar{{\bf k}}}n_{F}^{\prime}(\xi_{\bf k})$, where $\xi_{\bf
k}=k^{2}/2m-\mu$. At $T=0$, $n_{F}(\xi_{\bf k})\approx\Theta(-\xi_{\bf k})$
and $n_{F}^{\prime}(\xi_{\bf k})\approx-\delta(\xi_{\bf k})$. One can see that
all the bubbles above (75) show resonances at $\Omega=\Delta_{z}$ and
$\Omega=\sqrt{\Delta_{z}^{2}+\Delta_{R}^{2}}$, where $\Omega$ is the retarded
frequency, which are identified as lower and upper edges of the spin-flip
single-particle continuum. The form of the coefficient matrix $\hat{M}$
suggests decoupling of the in-plane, {x, y}, and out-of-plane, {z}, sectors:

$\text{Det}[\hat{M}^{-1}]=\bigg{(}1+\frac{U}{2}\Pi_{zz}\bigg{)}\bigg{[}\bigg{(}1+\frac{U}{2}\Pi_{xx}\bigg{)}^{2}+\frac{U^{2}}{4}\Pi_{xy}^{2}\bigg{]}=0.$
(76)

At weak coupling ($u\ll 1$, where $u\equiv U\nu_{F}/4>0$), there exists only
one mode, in the in-plane sector, which is robust and remains exponentially
close to the lower edge of the continuum, as given by Eq. (11) of the MT. The
mode in the out-of-plane sector, on the other hand, is Landau damped at weak
coupling; it is peeled off the continuum only at strong enough coupling which
is beyond the regime of validity of our theory. The exact analytical form of
collective modes from Eq. (76) cannot be calculated because of the
transcendental nature of the equation for $\Omega$ that one would encounter at
some point. However, near the lower edge of the continuum, i.e.,
$\Omega\approx\Delta_{Z}$, the analytical form of this mode can certainly be
obtained by assuming $\Omega=\Delta_{Z}-E_{B}$, with $E_{B}\ll\Delta_{Z}$ as
the binding energy. Expanding Eq. (75) for small $E_{B}$, the form of bubbles
(75) (up to Log singularity) can be written as

$\begin{split}\Pi_{xx}&\approx\frac{\nu_{F}}{4}\bigg{(}-2+\frac{2\Delta_{Z}^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)},\\\
\Pi_{xy}&\approx\frac{\nu_{F}}{4}\bigg{(}\frac{2i\Delta_{Z}^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)},\\\
\Pi_{zz}&\approx\frac{\nu_{F}}{4}\bigg{(}-4-\frac{4\Delta_{Z}E_{B}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}.\end{split}$
(77)

As one can see that for $E_{B}\to 0$, the Log in the above equation diverges.
One can immediately find from the formula of $\Pi_{zz}$ in Eq. (77) that the
Log-singular term is vanishingly small due to the multiplicative factor
$E_{B}$, unlike that in the planar components of the bubble. This tells us
that for the mode to be well lived in this sector, the coupling $u$ has to be
strong to reinforce the smallness imposed by the binding energy. But then at
strong coupling, the in-plane mode tends to become soft indicating
ferromagnetic instability (in our case of magnetized polar metals, this
instability could be ferroelectric-like), as discussed in the context of 2D
semiconductor heterostructures [14, 15]. So, we remain at weak coupling and
consider only the mode that exists in the in-plane sector below the continuum;
the out of plane mode exists inside the continuum at weak coupling and is
Landau damped. In the in-plane sector, essentially we are required to solve

$1+\frac{U}{2}\Pi_{xx}=\pm i\frac{U}{2}\Pi_{xy}$ (78)

according to Eq. (76). The one with the minus sign in Eq. (78) does not result
in any mode because the term containing leading Log-singularity cancels
identically. The one with the plus sign gives collective mode as written in
Eq. (11) of the MT. For further convenience, let’s rewrite the form of this
mode explicitly:

$\Omega\approx\Delta_{Z}-\frac{2\Delta_{R}^{2}}{\Delta_{Z}}\text{Exp}\Big{[}-\frac{\Delta_{R}^{2}(1-u)}{\Delta_{Z}^{2}2u}\Big{]}.$
(79)

Note that one cannot recover Silin-Leggett mode by substituting $\Delta_{R}=0$
in the above frequency (79): while writing (77), we assumed $E_{B}$ as the
smallest energy scale such that in addition to $E_{B}\ll\Delta_{Z}$ we also
have $E_{B}\ll\Delta_{R}$. Putting $\Delta_{R}=0$ in Eq. (79) will simply
violate our assumption, $E_{B}\ll\Delta_{Z},\Delta_{R}$.

###  III.2 Chiral spin waves in the presence of tilted magnetic field

We now calculate the dispersion of the collective mode (79) at small $q$.
Moreover, we now also turn on the in-plane component of the magnetic field,
$\Delta_{y}$, assuming $\Delta_{y}\ll\Delta_{Z},\Delta_{R}$. At finite $q$ and
$\Delta_{y}$, all the components of $\hat{\Pi}^{0}$, including the one in the
charge sector, are non-zero. Since now $\Delta_{y}$ is non-zero, the mirror
symmetry along $\hat{x}$ is broken. The dispersion of the collective mode
along $\hat{x}$\- and $\hat{y}$-directions, therefore, will not be the same,
same as what we saw for the DPV model (44). Since $\Delta_{y}$ gives a boost
along $\hat{x}$, see the single-particle Hamiltonian (67) and its eigenvalue
(40) where $\Delta_{y}$ dresses up $k_{x}$ to
$\bar{k}_{x}=k_{x}-\Delta_{y}/2\alpha$, the dispersion along $\hat{x}$ should
now start off linearly. In the $\hat{y}$-direction, however, it should still
start off quadratically only, as discussed in the context of DPV model (Sec.
II.1). The linear dispersion along $\hat{x}$ can be understood from symmetry:
the presence of $\Delta_{y}$ ensures broken mirror symmetry along the
$\hat{x}$ which explains the origin of boost along $\hat{x}$-direction in the
dispersion. Since the mirror symmetry along $\hat{y}$ is preserved, there is
no boost and, therefore, the dispersion should be quadratic. In the following,
we will calculate the dispersion of collective mode by first considering ${\bf
q}\parallel\hat{x}$ (Sec. III.2.1), followed by ${\bf q}\parallel\hat{y}$
(Sec. III.2.2) and then ${\bf q}\parallel\hat{z}$ (Sec. III.2.3).

####  III.2.1 ${\bf q}\parallel\hat{x}$

The bare polarization bubble $\hat{\Pi}^{0}$ at arbitrary $q$, assuming it to
be along $x$-direction, can be written as

$\begin{split}\Pi_{ij}({\bf
q},i\Omega_{n})&=\frac{1}{2}\sum_{s,\bar{s}}\int_{K}f_{ij}g_{\bar{s}}(i\omega_{m},{\bf
k})g_{s}\big{(}i(\omega_{m}+\Omega_{n}),{\bf k}+{\bf q}\big{)},\\\
&=\frac{1}{2}\sum_{s,\bar{s}}\int_{K}f_{ij}\frac{1}{i\omega_{m}-\xi_{\bf
k}^{\bar{s}}}\frac{1}{i(\omega_{m}+\Omega_{n})-\xi_{{\bf k}+{\bf
q}}^{s}},\end{split}$ (80)

where $\xi_{\bf k}^{s}=\varepsilon_{\bf k}^{s}-\mu$, with $\varepsilon_{\bf
k}^{s}$ as the eigenenergy (68), and $f_{ij}$ is the coherence factor
calculated as

$\begin{split}f_{00}&=1+s\bar{s}\frac{\Delta_{Z}^{2}+4\alpha^{2}\big{(}\bar{k}_{x}(\bar{k}_{x}+q)+k_{y}^{2}\big{)}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{0x}&=s\frac{2\alpha
k_{y}}{\epsilon_{\bar{k}_{x}+q,k_{y}}}+\bar{s}\frac{2\alpha
k_{y}}{\epsilon_{\bar{k}_{x},k_{y}}}-is\bar{s}\frac{\Delta_{Z}2\alpha
q}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{x0}&=f_{0x}(i\rightarrow-i),\\\
f_{0y}&=-s\frac{2\alpha(\bar{k}_{x}+q)}{\epsilon_{\bar{k}_{x}+q,k_{y}}}-\bar{s}\frac{2\alpha\bar{k}_{x}}{\epsilon_{\bar{k}_{x},k_{y}}},\\\
f_{y0}&=f_{0y},\\\
f_{0z}&=(s+\bar{s})\frac{\Delta_{Z}}{\epsilon_{\bar{k}_{x}+q,k_{y}}}+is\bar{s}\frac{4\alpha^{2}k_{y}q}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{z0}&=f_{0z}(i\rightarrow-i),\\\
f_{xx}&=1-s\bar{s}\frac{\Delta_{Z}^{2}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}}+s\bar{s}\frac{4\alpha^{2}\big{(}k_{y}^{2}-\bar{k}_{x}(\bar{k}_{x}+q)\big{)}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{yy}&=1-s\bar{s}\frac{\Delta_{Z}^{2}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}}-s\bar{s}\frac{4\alpha^{2}\big{(}k_{y}^{2}-\bar{k}_{x}(\bar{k}_{x}+q)\big{)}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{zz}&=1+s\bar{s}\frac{\Delta_{Z}^{2}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}}-s\bar{s}\frac{4\alpha^{2}\big{(}k_{y}^{2}+\bar{k}_{x}(\bar{k}_{x}+q)\big{)}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{xy}&=is\frac{\Delta_{Z}}{\epsilon_{\bar{k}_{x}+q,k_{y}}}-i\bar{s}\frac{\Delta_{Z}}{\epsilon_{\bar{k}_{x},k_{y}}}-s\bar{s}\frac{4\alpha^{2}\big{(}k_{y}(\bar{k}_{x}+q)+\bar{k}_{x}k_{y}\big{)}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{yx}&=f_{xy}(i\rightarrow-i),\\\
f_{xz}&=is\frac{2\alpha(\bar{k}_{x}+q)}{\epsilon_{\bar{k}_{x}+q,k_{y}}}-i\bar{s}\frac{2\alpha\bar{k}_{x}}{\epsilon_{\bar{k}_{x},k_{y}}}+s\bar{s}\frac{2\Delta_{Z}2\alpha
k_{y}}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{zx}&=f_{xz}(i\rightarrow-i),\\\ f_{yz}&=is\frac{2\alpha
k_{y}}{\epsilon_{\bar{k}_{x}+q,k_{y}}}-i\bar{s}\frac{2\alpha
k_{y}}{\epsilon_{\bar{k}_{x},k_{y}}}-s\bar{s}\frac{\Delta_{Z}2\alpha(2\bar{k}_{x}+q)}{\epsilon_{\bar{k}_{x},k_{y}}\epsilon_{\bar{k}_{x}+q,k_{y}}},\\\
f_{zy}&=f_{yz}(i\rightarrow-i).\end{split}$ (81)

In the regime $q^{2}/2m\ll v_{F}q\ll\Delta_{y}\ll\Delta_{Z},\Delta_{R}$, the
form of all the bubbles near the lower edge of the spin-flip continuum
($\Omega=\Delta_{Z}$ is still the location of the lower edge of the continuum
even at finite $\Delta_{y}$ because of our assumption
$\Delta_{y}\ll\Delta_{R}$ [16]) up to leading singularity can be written as

$\begin{split}\Pi_{00}&=\frac{\nu_{F}}{8}\Bigg{[}q^{2}\frac{4\alpha^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}+q^{2}\Delta_{y}^{2}\frac{4\alpha^{2}}{2\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\Bigg{]},\\\
\Pi_{0x}&=\frac{\nu_{F}}{8}\Bigg{[}iq\frac{4\alpha\Delta_{Z}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}+iq\Delta_{y}^{2}\frac{4\alpha\Delta_{Z}}{2\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\Bigg{]},\\\
\Pi_{x0}&=-\Pi_{0x},\\\
\Pi_{0y}&=\frac{\nu_{F}}{8}\Bigg{[}-q\frac{4\alpha\Delta_{Z}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}-q\Delta_{y}^{2}\frac{4\alpha\Delta_{Z}}{2\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\Bigg{]},\\\
\Pi_{y0}&=\Pi_{0y},\\\
\Pi_{0z}&=\frac{\nu_{F}}{8}\Bigg{[}-iq\frac{4\alpha}{\Delta_{R}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}-iq\Delta_{y}^{2}\frac{4\alpha}{2\Delta_{R}^{3}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\Bigg{]},\\\
\Pi_{z0}&=-\Pi_{0z},\\\
\Pi_{xx}&=\frac{\nu_{F}}{8}\Bigg{[}\bigg{(}-4+\frac{4\Delta_{Z}^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}+\Delta_{y}^{2}\bigg{(}\frac{2\Delta_{Z}^{2}}{\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}\\\
&\hskip
28.45274pt+\frac{v_{F}q\Delta_{y}}{E_{B}}\frac{4\Delta_{Z}^{2}}{\Delta_{R}^{3}}+\frac{2v_{F}^{2}q^{2}}{\Delta_{R}^{2}}\frac{\Delta_{Z}(4m^{2}\alpha^{4}-\Delta_{Z}^{2})}{E_{B}\Delta_{R}^{2}}+\mathcal{O}(q^{2}\Delta_{y}^{2})\Bigg{]},\\\
\Pi_{yy}&=\Pi_{xx},\\\
\Pi_{zz}&=\frac{\nu_{F}}{8}\Bigg{[}\bigg{(}-8-\frac{8\Delta_{Z}E_{B}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-\Delta_{y}^{2}\bigg{(}\frac{2}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}\\\
&\hskip 28.45274pt-
v_{F}q\Delta_{y}\bigg{(}\frac{8\Delta_{Z}}{\Delta_{R}^{3}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-\frac{4v_{F}^{2}q^{2}}{\Delta_{R}^{2}}\bigg{(}\frac{4m^{2}\alpha^{4}-2\Delta_{Z}^{2}}{\Delta_{R}^{2}}\bigg{)}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}+\mathcal{O}(q^{2}\Delta_{y}^{2})\Bigg{]},\\\
\Pi_{xy}&=\frac{\nu_{F}}{8}\Bigg{[}\bigg{(}i\frac{4\Delta_{Z}^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}+\Delta_{y}^{2}\bigg{(}i\frac{2\Delta_{Z}^{2}}{\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}\\\
&\hskip
28.45274pt+i\frac{v_{F}q\Delta_{y}}{E_{B}}\frac{4\Delta_{Z}^{2}}{\Delta_{R}^{3}}+i\frac{2v_{F}^{2}q^{2}}{\Delta_{R}^{2}}\frac{\Delta_{z}(4m^{2}\alpha^{4}-\Delta_{Z}^{2})}{E_{B}\Delta_{R}^{2}}+\mathcal{O}(q^{2}\Delta_{y}^{2})\Bigg{]},\\\
\Pi_{yx}&=-\Pi_{xy},\\\
\Pi_{xz}&=\frac{\nu_{F}}{8}\Bigg{[}i\frac{4\Delta_{y}\Delta_{Z}}{\Delta_{R}^{2}}\bigg{(}1-\frac{E_{B}\Delta_{Z}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-i\frac{4v_{F}q\Delta_{Z}^{2}}{\Delta_{R}^{3}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}+\mathcal{O}(q\Delta_{y}^{2})\Bigg{]},\\\
\Pi_{zx}&=-\Pi_{xz},\\\
\Pi_{yz}&=\frac{\nu_{F}}{8}\Bigg{[}\frac{4\Delta_{y}\Delta_{Z}}{\Delta_{R}^{2}}\bigg{(}1-\frac{E_{B}\Delta_{Z}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-\frac{4v_{F}q\Delta_{Z}^{2}}{\Delta_{R}^{3}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}+\mathcal{O}(q\Delta_{y}^{2})\Bigg{]},\\\
\Pi_{zy}&=\Pi_{yz},\end{split}$ (82)

where $E_{B}$ is the binding energy with the assumption
$E_{B}\ll\Delta_{Z},\Delta_{R}$. As we see, all the components of the bubble
are non-zero at finite $q$ and in the presence of $\Delta_{y}$, up to
$\mathcal{O}(q\Delta_{y})$ and $\mathcal{O}(q^{2})$. However, we argue that
near the lower edge of the continuum, i.e., $\Omega\approx\Delta_{Z}$, the
determinant equation can still be decoupled as in Eq. (76).

To begin the discussion, we first consider only $\Delta_{y}\neq 0$ and set
$q=0$. All the charge components of the bubble ($\Pi_{0i}$, $i\in(0...3)$) are
zero at $q=0$. The structure of the bubble matrix is now such that the in-
plane and the out-of-plane sectors are coupled; this is ensured by $\Pi_{xz}$
and $\Pi_{yz}$ being non-zero. However, it can be argued that at
$\Omega\approx\Delta_{Z}$, $\Pi_{xz}$ and $\Pi_{yz}$ do not affect the
determinant equation, $\text{Det}[\hat{M}^{-1}]=0$ (74), significantly.
Indeed, the determinant equation that is to be solved for $\Omega$ can be
written as

$\bigg{(}1+\frac{U}{2}\Pi_{zz}\bigg{)}\bigg{[}\bigg{(}1+\frac{U}{2}\Pi_{xx}\bigg{)}^{2}+\frac{U^{2}}{4}\Pi_{xy}^{2}\bigg{]}-\frac{U^{3}}{4}\Pi_{xy}\Pi_{yz}\Pi_{xz}+\frac{U^{2}}{4}\bigg{(}1+\frac{U}{2}\Pi_{zz}\bigg{)}\big{(}\Pi_{xz}^{2}-\Pi_{yz}^{2}\big{)}=0$
(83)

Clearly, for $E_{B}\to 0$, the first term of Eq. (83) diverges as
$\text{log}^{2}E_{B}$ while the second term as $\text{log}E_{B}$: log-
divergence in $\Pi_{iz}$, for $i\in(1...3)$, vanishes as
$E_{B}\text{log}E_{B}$, so their leading contribution is only up to some
finite value. This is the reason the third term of Eq. (83) has no divergence
and starts off with a finite value. Hence, compared to the first term, the
second term and the third term of Eq. (83) can be ignored, hence again
recovering the decoupled form of the determinant equation to be solved for the
collective mode frequency $\Omega$. At weak coupling ($u\ll 1$), the finite
$\Delta_{y}$ correction to the collective mode frequency (79), therefore,
mainly comes from the corresponding corrections in $\Pi_{xx}$ and $\Pi_{xy}$.

Now we consider finite $q$ and present the argument on decoupling of {x, y}-
and {z}-sectors once again. We can observe that all the charge components of
the bubble, $\Pi_{0i}$ for $i\in(0...3)$, diverges as
$q^{\alpha}\text{log}E_{B}$. In the spin components, the planar components
($\Pi_{xx}$ and $\Pi_{xy}$) diverge as $q^{\alpha}/E_{B}$, while the out-of-
plane components ($\Pi_{zz}$, $\Pi_{xz}$ and $\Pi_{yz}$) diverge as
$q^{\alpha}\text{log}E_{B}$. Clearly, $1/E_{B}$ divergence is stronger than
$\text{log}E_{B}$, so at finite $q$ the last two terms of Eq. (83) can be
again ignored compared to the first term.

Finally, at $\Omega\approx\Delta_{Z}$, the determinant equation (83) can be
written in a decoupled form, i.e., keeping only the first term at leading
order, and the dispersion of the collective mode (79) can be calculated for
tilted magnetic field:

$\Omega\approx\Delta_{Z}-\frac{2\Delta_{R}^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{\Delta_{R}^{2}}{\Delta_{Z}^{2}}\frac{1-u}{2u}\left(1+\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\right)^{-1}\right]+v_{F}q_{x}\bigg{(}\frac{\Delta_{y}}{\Delta_{R}}\bigg{)}+v_{F}^{2}q_{x}^{2}\bigg{(}\frac{4m^{2}\alpha^{4}-\Delta_{Z}^{2}}{2\Delta_{R}^{2}\Delta_{Z}}\bigg{)}+\mathcal{O}(q_{x}^{2}\Delta_{y}^{2}).$
(84)

As expected, the leading order $q$-dependent correction to the mode frequency
starts off linearly with the amplitude proportional to $\Delta_{y}$,
consistent with the results obtained for the DPV model (44).

####  III.2.2 ${\bf q}\parallel\hat{y}$

One can similarly calculate the dispersion along $\hat{y}$-direction. However,
instead of doing it explicitly, one can directly write down its form knowing
that there will not be any linear in $q_{y}$ (boost) term. Indeed, as
discussed in the beginning of this section, $\Delta_{y}$ breaks mirror
symmetry along $\hat{x}$-direction. This naturally explains the leading order
$q$-dependent correction to the collective mode frequency be linear in
$q_{x}$. However, since the mirror symmetry along $\hat{y}$ is preserved, the
dispersion should be quadratic in $q_{y}$. From the explicit calculation, we
find

$\Omega\approx\Delta_{Z}-\frac{2\Delta_{R}^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{\Delta_{R}^{2}}{\Delta_{Z}^{2}}\frac{1-u}{2u}\left(1+\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\right)^{-1}\right]+v_{F}^{2}q_{y}^{2}\bigg{(}\frac{4m^{2}\alpha^{4}-\Delta_{Z}^{2}}{2\Delta_{R}^{2}\Delta_{Z}}\bigg{)}+\mathcal{O}(q_{y}^{2}\Delta_{y}^{2}),$
(85)

which is same as that for ${\bf q}\parallel\hat{x}$, except for the linear in
$q_{x}$ term. We note that the term of order $q_{x,y}^{2}\Delta_{y}^{2}$ is
not necessarily going to be the same for $q_{x}$ and $q_{y}$.

The coefficient of $q_{x}^{2}$ (84) and $q_{y}^{2}$ (85) terms in the
collective mode has a turning point which has a crucial implication. If for
the time being we assume $\Delta_{y}=0$, then the dispersion along $\hat{x}$
and $\hat{y}$ is quadratic and the coefficient of it is positive when
$\Delta_{Z}<2m\alpha^{2}$. This is crucial in the sense that if the dispersion
along $\hat{z}$ is always monotonic, then the overall dispersion of the
collective mode is hyperbolic in nature as long as $\Delta_{Z}<2m\alpha^{2}$.
To see this explicitly, we calculate the dispersion along $\hat{z}$-direction
in the next section.

####  III.2.3 ${\bf q}\parallel\hat{z}$

The general equation of non-interacting polarization bubbles for ${\bf q}$
along $\hat{z}$-direction can be obtained by writing $k_{z}$ as $k_{z}+q_{z}$.
Since the coherence factor $f_{ij}$ (81) depends only on in-plane components
of ${\bf k}$, the $q$-dependence to the bubble will not arise from here. The
only source for $q$-dependence then is the chiral part of the Green’s
function, $g_{s}(i\epsilon_{m},{\bf k})$, as defined in Eq. (71). Near the
lower edge of the spin-flip continuum, $\Omega\approx\Delta_{Z}$, which is
possible for interband transitions only, the contribution from charge sector
vanishes identically which can be understood from Eq. (81) by putting $q=0$
and $s=-\bar{s}$. Therefore, only the spin contribution survives, the
corresponding coherence factors is same as that in Eq. (81) at $q=0$. The
polarization bubbles can be calculated up to leading singularity as

$\begin{split}\Pi_{xx}&=\frac{\nu_{F}}{8}\Bigg{[}\bigg{(}-4+\frac{4\Delta_{Z}^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}+\Delta_{y}^{2}\bigg{(}\frac{2\Delta_{Z}^{2}}{\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-\frac{2v_{F}^{2}q^{2}}{E_{B}^{2}}\frac{\Delta_{Z}^{2}}{\Delta_{R}^{2}}\bigg{(}1-\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\bigg{)}\Bigg{]},\\\
\Pi_{yy}&=\Pi_{xx},\\\
\Pi_{zz}&=\frac{\nu_{F}}{8}\Bigg{[}\bigg{(}-8-\frac{8\Delta_{Z}E_{B}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-\Delta_{y}^{2}\bigg{(}\frac{2}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}-\frac{4v_{F}^{2}q^{2}}{\Delta_{R}^{2}}\frac{\Delta_{Z}}{E_{B}}+\frac{v_{F}^{2}q^{2}}{\Delta_{R}^{2}}\frac{\Delta_{y}^{2}}{E_{B}^{2}}\Bigg{]},\\\
\Pi_{xy}&=\frac{\nu_{F}}{8}\Bigg{[}i\frac{4\Delta_{Z}^{2}}{\Delta_{R}^{2}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}+i\Delta_{y}^{2}\frac{2\Delta_{Z}^{2}}{\Delta_{R}^{4}}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}-i\frac{2v_{F}^{2}q^{2}}{E_{B}^{2}}\frac{\Delta_{Z}^{2}}{\Delta_{R}^{2}}\bigg{(}1-\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\bigg{)}\Bigg{]},\\\
\Pi_{yx}&=-\Pi_{xy},\\\
\Pi_{xz}&=\frac{\nu_{F}}{8}\Bigg{[}i\frac{4\Delta_{y}\Delta_{Z}}{\Delta_{R}^{4}}\bigg{(}\Delta_{R}^{2}-E_{B}\Delta_{Z}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}+i\Delta_{y}\frac{2v_{F}^{2}q^{2}\Delta_{Z}^{2}}{\Delta_{R}^{4}E_{B}}\Bigg{]},\\\
\Pi_{zx}&=-\Pi_{xz},\\\
\Pi_{yz}&=\frac{\nu_{F}}{8}\Bigg{[}\frac{4\Delta_{y}\Delta_{Z}}{\Delta_{R}^{4}}\bigg{(}\Delta_{R}^{2}-E_{B}\Delta_{Z}\text{Log}\frac{\Delta_{Z}E_{B}}{2\Delta_{R}^{2}}\bigg{)}+\Delta_{y}\frac{2v_{F}^{2}q^{2}\Delta_{Z}^{2}}{\Delta_{R}^{4}E_{B}}\Bigg{]},\\\
\Pi_{zy}&=\Pi_{yz}.\end{split}$ (86)

At small $q$ and $\Delta_{y}$, the determinant equation (83) (all the charge
components are zero, so this equation is valid also at finite $q_{z}$) can be
argued to decouple in the same way as discussed in Sec. III.2.1. Assuming
$q=0$ first, it can be observed in Eq. (86) that $\Pi_{zz}$, $\Pi_{xz}$ and
$\Pi_{yz}$ are finite as $E_{B}\to 0$. On the other hand, all the planar
components, such as $\Pi_{xx}$ and $\Pi_{xy}$ diverges logarithmically as
$\text{log}E_{B}$. This is the reason we can ignore $\Pi_{xz}$ and $\Pi_{yz}$
and consider the determinant equation (83) as decoupled. Now we consider
finite $q$. As contrary to that for $q_{x}$, the $\Pi_{xz}$ and $\Pi_{yz}$
components also have $q^{\alpha}/E_{B}$ divergence. However, now the planar
components $\Pi_{xx}$ and $\Pi_{xy}$ diverge as $q^{\alpha}/E_{B}^{2}$ which
is stronger that $q^{\alpha}/E_{B}$. This is the reason we can again consider
(also at finite $q_{z}$) the determinant equation (83) as decoupled. Using Eq.
(78), with the form of associated polarization bubbles from Eq. (86), the
dispersion of the collective mode (79) along $\hat{z}$-direction can be
calculated:

$\begin{split}\Omega\approx\Delta_{Z}-\frac{2\Delta_{R}^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{\Delta_{R}^{2}}{\Delta_{Z}^{2}}\frac{1-u}{2u}\left(1+\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\right)^{-1}\right]-\frac{v_{F}^{2}q_{z}^{2}}{\frac{4\Delta_{R}^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{\Delta_{R}^{2}}{\Delta_{Z}^{2}}\frac{1-u}{2u}\left(1+\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\right)^{-1}\right]}\left(1-\frac{\Delta_{y}^{2}}{\Delta_{R}^{2}}\right).\end{split}$
(87)

The dispersion of the mode (87) is monotonic and negative along
$\hat{z}$-direction, which, combined with the argument presented above for
${\bf q}$ along $\hat{x}$ (84) and $\hat{y}$ (85) indicate hyperbolicity
(assuming $\Delta_{y}=0$) only when $\Delta_{Z}<2m\alpha^{2}$. If
$\Delta_{Z}>2m\alpha^{2}$, then the mode dispersion is parabolic. This is also
consistent with our result for the DPV model as obtained in Eq. (44).

We emphasize that one cannot approach to the Silin-Leggett limit here by
putting $\Delta_{R}=0$ due to the same reason mentioned just below Eq. (79).
Moreover, the approach to the non-interacting limit by putting $u=0$ is also
not valid here because the $q$-dependent correction to the binding energy
($E_{B}$) is obtained perturbatively to the $q=0$ term which required
$E_{B}(q=0)\neq 0$. As we can see from the above equation that since
$E_{B}(q=0)\sim e^{-1/u}$, the $u=0$ limit makes $E_{B}(q=0)=0$, which would
violate the perturbation limit above.

Combining Eqs. (84), (85) and (87), the overall dispersion of the collective
mode for the spherical FS model can be written as

$\begin{split}\Omega&\approx\Delta_{z}-p_{1}-p_{2}\delta_{q_{y},0}\delta_{q_{z},0}+\frac{1}{p_{2}}\frac{\Delta_{y}^{2}}{4\Delta_{R}^{2}}\left[\left(v_{F}q_{x}+p_{2}\frac{2\Delta_{R}}{\Delta_{y}}\right)^{2}+v_{F}^{2}q_{y}^{2}\right]-v_{F}^{2}q_{z}^{2}\frac{1}{2p_{1}}\left(1-\frac{\Delta_{y}^{2}}{\Delta_{R}^{2}}\right),\\\
&\hskip
14.22636pt\text{where}\,\,p_{1}=\frac{2\Delta_{R}^{2}}{\Delta_{Z}}\text{Exp}\left[-\frac{\Delta_{R}^{2}}{\Delta_{Z}^{2}}\frac{1-u}{2u}\left(1+\frac{\Delta_{y}^{2}}{2\Delta_{R}^{2}}\right)^{-1}\right]\,\,\text{and}\,\,p_{2}=\frac{\Delta_{y}^{2}\Delta_{Z}}{2(4m^{2}\alpha^{4}-\Delta_{Z}^{2})}\\\
&\hskip 14.22636pt\text{with}\,\,\delta_{q_{i},0}\,\,\text{as Kronecker
delta}.\end{split}$ (88)

## References

  * Garate and Franz [2011] I. Garate and M. Franz, Phys. Rev. B 84, 045403 (2011).
  * Efimkin  _et al._ [2012] D. K. Efimkin, Y. E. Lozovik, and A. A. Sokolik, Nanoscale research letters 7, 163 (2012).
  * Egri [1985] I. Egri, Physics Reports 119, 363 (1985).
  * Silin [1957] V. P. Silin, J. Exptl. Theoret. Phys. (U.S.S.R.) 33, 1227 (1957).
  * Maiti  _et al._ [2016a] S. Maiti, M. Imran, and D. L. Maslov, Phys. Rev. B 93, 045134 (2016a).
  * Conti and Trillo [2004] C. Conti and S. Trillo, Phys. Rev. Lett. 92, 120404 (2004).
  * Xu  _et al._ [2020] H. Xu, C. Hang, and G. Huang, Phys. Rev. A 101, 053832 (2020).
  * Conti  _et al._ [2003] C. Conti, S. Trillo, P. Di Trapani, G. Valiulis, A. Piskarskas, O. Jedrkiewicz, and J. Trull, Phys. Rev. Lett. 90, 170406 (2003).
  * Efremidis  _et al._ [2009] N. K. Efremidis, G. A. Siviloglou, and D. N. Christodoulides, Physics Letters A 373, 4073 (2009).
  * Sedov  _et al._ [2015] E. S. Sedov, I. V. Iorsh, S. M. Arakelian, A. P. Alodjants, and A. Kavokin, Phys. Rev. Lett. 114, 237402 (2015).
  * Voronych  _et al._ [2016] O. Voronych, A. Buraczewski, M. Matuszewski, and M. Stobińska, Phys. Rev. B 93, 245310 (2016).
  * Gianfrate  _et al._ [2018] A. Gianfrate, L. Dominici, O. Voronych, M. Matuszewski, M. Stobińska, D. Ballarini, M. De Giorgi, G. Gigli, and D. Sanvitto, Light: Science & Applications 7, 17119 (2018).
  * Maiti  _et al._ [2015] S. Maiti, V. Zyuzin, and D. L. Maslov, Phys. Rev. B 91, 035106 (2015).
  * Maiti  _et al._ [2016b] S. Maiti, M. Imran, and D. L. Maslov, Phys. Rev. B 93, 045134 (2016b).
  * Kumar and Maslov [2017] A. Kumar and D. L. Maslov, Phys. Rev. B 95, 165140 (2017).
  * Kumar  _et al._ [2023] A. Kumar, P. Chandra, and P. A. Volkov, Phys. Rev. B 108, 075162 (2023).