Zero-field spin resonance in graphene
with proximity-induced spin-orbit coupling

We consider a monolayer graphene attached to a substrate made of, e.g, TMD, or heavy metal. Strong SOC in the substrate induces SOC in graphene, which can be generically of both Rashba and VZ types. For completeness, we also allow for an intrinsic KM term. Following Refs. Kane and Mele, 2005a; Rashba, 2009; Stauber and Schliemann, 2009; Wang et al., 2015; Cummings et al., 2017; Garcia et al., 2018, we adopt the following low-energy Hamiltonian:

$\displaystyle\hat{H}_{0}$

$\displaystyle=$

$\displaystyle v_{F}(\tau_{z}\hat{s}_{0}\hat{\sigma}_{x}k_{x}+\hat{s}_{0}\hat{\sigma}_{y}k_{y})+\Delta\hat{s}_{0}\hat{\sigma}_{z}+\frac{\lambda_{\rm KM}}{2}\tau_{z}\hat{s}_{z}\hat{\sigma}_{z}+\frac{\lambda_{\rm R}}{2}(\tau_{z}\hat{s}_{y}\hat{\sigma}_{x}-\hat{s}_{x}\hat{\sigma}_{y})+\frac{\lambda_{\rm Z}}{2}\tau_{z}\hat{s}_{z}\hat{\sigma}_{0},$

(1)

where $v_{F}$ is the Dirac velocity; k is the electron momentum measured either from the $K$ or $K^{\prime}$ point of the graphene Brillouin zone, $\lambda_{\rm KM,Z,R}$ are the coupling constants of the KM, VZ, and Rashba spin-orbit interactions, respectively, $\Delta$ is the gap due to substrate-induced asymmetry between the A and B sites of the honeycomb lattice, $\hat{\sigma}_{i}$ and $\hat{s}_{i}$ are the Pauli matrices in the sublattice (pseudospin) and spin spaces, respectively (with $\hat{\sigma}_{0}$ and $\hat{s}_{0}$ being the unity matrices in the corresponding spaces), and $\tau_{z}=\pm 1$ labels the $K$ and $K^{\prime}$ points.

The first three terms in Eq. (1) describe (massive) Dirac fermions near the $K$ and $K^{\prime}$ points, the rest of the terms describe SOC. The intrinsic SOC present locally on the sublattice sites of graphene gives rise to two terms in the Hamiltonian, which are symmetric and antisymmetric combinations of local SOCs on the A and B sites. Because free-standing graphene is invariant under sublattice exchange, only the symmetric combination, which was identified by Kane and Mele,Kane and Mele (2005a) survives. It couples spins, sublattices (pseudo-spins), and valleys, as indicated by the fourth term in Eq. (1). The presence of a substrate brings about new features. First, the breaking of $z\rightarrow-z$ inversion symmetry induces a Rashba-type SOC. At the $K$and $K^{\prime}$ points, the Rashba Hamiltonian contains the leading, momentum-independent term [the fifth term in Eq. (1)] and the subleading, linear-in-$k$ term.Kane and Mele (2005b); Rashba (2009) We assume that doping is low enough, i.e., $k_{F}\ll\lambda_{R}/\alpha$, where $\alpha$ is the Rashba parameter for the linear-in-$k$ coupling, such that the latter can be neglected, and Rashba SOC will be described by the fifth term in Eq. (1). Second, if the substrate also breaks the sublattice symmetry, the anti-symmetric combination of atomic SOCs gives rise to a VZ SOC–the last term in Eq. (1).

The Hamiltonian in Eq. (1) can be simplified further. First, it is well known that the KM coupling is much weaker than other types of SOC. While theoretical estimates place $\lambda_{\rm KM}$ in the range from 1 $\mu$eV Min et al. (2006); Yao et al. (2007) to 25-50 $\mu$eV,Boettger and Trickey (2007); Konschuh et al. (2010) a recent ESR experiment reports the value of $42.2$  $\mu$eV. Sichau et al. (2019) This is much smaller than typical values of 1-10 meV for $\lambda_{\rm R}$ and $\lambda_{\rm Z}$ for graphene on TMD substrates, Wang et al. (2015, 2016); Cummings et al. (2017); Wakamura et al. (2018); Zihlmann et al. (2018); Wakamura et al. (2019) and thus the KM term will be ignored in what follows. Next, substrate-induced sublattice asymmetry is expected to open a gap of magnitude $\Delta$ at the Dirac points. This gap endows graphene with Berry curvature which, in combination with SOC, leads to interesting consequences for EDSR, such as a Hall component of the induced current, as discussed recently in Ref. Raines et al., . In the regime of $\lambda_{\mathrm{R}},\lambda_{\mathrm{Z}}\ll\mu$ (where the FL theory of Secs. III.2-III.4 is valid), the presence of the gap gives rise to the Hall conductivity but does not qualitatively affect the collective modes, as it amounts only to changes in the Fermi velocity and SOC parameters. Given also that angular-resolved photoemission of graphene on TMD substrates has not detected gaps in the Dirac spectrum, Coy Diaz et al. (2015); Pierucci et al. (2016); Henck et al. (2018) we will neglect the asymmetry gap in this study. With these simplifications, the Hamiltonian is reduced to

$\displaystyle\hat{H}_{0}$

$\displaystyle=$

$\displaystyle v_{F}(\tau_{z}\hat{s}_{0}\hat{\sigma}_{x}k_{x}+\hat{s}_{0}\hat{\sigma}_{y}k_{y})+\frac{\lambda_{\rm Z}}{2}\tau_{z}\hat{s}_{z}\hat{\sigma}_{0}+\frac{\lambda_{\rm R}}{2}(\tau_{z}\hat{s}_{y}\hat{\sigma}_{x}-\hat{s}_{x}\hat{\sigma}_{y}).$

(2)

If only Rashba SOC is present ($\lambda_{\rm Z}=0$), the eigenvalues and eigenstates of the Hamiltonian (2) are given byRashba (2009); Stauber and Schliemann (2009)

$$\begin{split}\varepsilon_{\alpha\beta;\tau_{z}}({\bf k})&=\alpha\left(\sqrt{v_{F}^{2}k^{2}+\left(\frac{\lambda_{\rm R}}{2}\right)^{2}}+\beta\frac{\lambda_{\rm R}}{2}\right);\\ \ket{\alpha\beta;\tau_{z}}&=\frac{1}{\sqrt{2(1+(\epsilon_{\alpha\beta})^{2\tau_{z}})}}\left(\begin{array}[]{c}-i\alpha\beta\tau_{z}e^{-i(1+\tau_{z})\theta}\\ -i\alpha\beta e^{-i\theta}(\epsilon_{\alpha\beta})^{\tau_{z}}\\ \tau_{z}e^{-i\tau_{z}\theta}(\epsilon_{\alpha\beta})^{\tau_{z}}\\ 1\end{array}\right),~{}\epsilon_{\alpha\beta}=\varepsilon_{\alpha\beta;\tau_{z}}({\bf k})/v_{F}k,\end{split}$$

(3)

where $\alpha=\pm 1$ denotes the conduction/valence band, $\beta=\pm 1$ denotes the SOC-split chiral subbands, $\tau_{z}=\pm 1$ denotes the $K/K^{\prime}$ valleys, respectively, and $\theta$ is the azimuthal angle of k.

The energy spectrum for a realistic case of $\lambda_{\mathrm{R}}<\mu$ is shown in Fig. 1a, where we choose $\mu>0$ without loss of generality. The chiral nature of the bands is evident from the expectation values of the spin operators in each of the subbands:

$$\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{x}|\alpha\beta;\tau_{z}\rangle=-\frac{\beta\sin\theta}{\sqrt{1+r^{2}}},~{}~{}\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{y}|\alpha\beta;\tau_{z}\rangle=\frac{\beta\cos\theta}{\sqrt{1+r^{2}}},~{}~{}\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{z}|\alpha\beta;\tau_{z}\rangle=0,$$

(4)

where $\hat{\mathcal{S}}_{a}\equiv\hat{s}_{a}\hat{\sigma}_{0}$ are the components of the spin operator and $r\equiv\lambda_{\rm R}/2v_{F}k$. It is clear from Eq. (4) that ${\bf k}\cdot\hat{\boldsymbol{\mathcal{S}}}=0$, which means that the spin is perpendicular to the momentum in each subband, as shown in Fig. 1c. Note that $\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{i}|\alpha\beta;\tau_{z}\rangle$ is independent of $\tau_{z}$. This means that the chiral structure is the same at the $K$ and $K^{\prime}$ points.

To determine which transitions can be excited by an ac magnetic field in an ESR measurement, we also need the intravalley matrix elements of the spin operator for off-diagonal transitions between spin-split subbands. Using the eigenstates from Eq. (3), we obtain

	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{x}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle\frac{i\left(\alpha^{\prime}\beta^{\prime}e^{i\theta}-\alpha\beta e^{-i\theta}\right)\left(\epsilon_{\alpha\beta}+\epsilon_{\alpha^{\prime}\beta^{\prime}}\right)}{2\sqrt{\left(1+\epsilon^{2}_{\alpha\beta}\right)\left(1+\epsilon^{2}_{\alpha^{\prime}\beta^{\prime}}\right)}},$		(5a)
	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{y}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle\frac{\left(\alpha^{\prime}\beta^{\prime}e^{i\theta}+\alpha\beta e^{-i\theta}\right)\left(\epsilon_{\alpha\beta}+\epsilon_{\alpha^{\prime}\beta^{\prime}}\right)}{2\sqrt{\left(1+\epsilon^{2}_{\alpha\beta}\right)\left(1+\epsilon^{2}_{\alpha^{\prime}\beta^{\prime}}\right)}},$		(5b)
	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{z}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle\frac{\left(\alpha\beta\alpha^{\prime}\beta^{\prime}-1\right)\left(1+\epsilon_{\alpha\beta}\epsilon_{\alpha^{\prime}\beta^{\prime}}\right)}{2\sqrt{\left(1+\epsilon^{2}_{\alpha\beta}\right)\left(1+\epsilon^{2}_{\alpha^{\prime}\beta^{\prime}}\right)}}.$		(5c)

A response to an in-plane ac magnetic field is controlled by the matrix elements in Eqs. (5a) and (5b). These vanish if $\epsilon_{\alpha\beta}+\epsilon_{\alpha^{\prime}\beta^{\prime}}=0$, which can happen only if $\beta=\beta^{\prime}$ and $\alpha=-\alpha^{\prime}$. Thus transitions between conduction and valence subbands of the same chirality are forbidden. From Fig. 1a, we see that the frequencies of the allowed transitions for an in-plane field are $\Omega=\lambda_{\mathrm{R}}$ and $\Omega=2\mu\pm\lambda_{\mathrm{R}}$. A response to an out-of-plane magnetic field is controlled by the matrix element in Eq. (5c), which is non-zero if $\alpha\beta\alpha^{\prime}\beta^{\prime}\neq 1$. This condition allows for transitions between the states with opposite subband (conduction vs valence) or chirality indices but not both. As shown in Fig. 1a, the frequencies of the allowed transitions for an out-of-plane field are $\Omega=\lambda_{\mathrm{R}}$ and $\Omega=2\mu$.

In an ESR experiment, one measures the imaginary part of the spin susceptibility, which can be calculated using the Kubo formula, see Appendix A.3.1. The results are shown in Figs. 2a and 2b. In agreement with the selection rules, there is a resonance at $\Omega=\lambda_{\mathrm{R}}$ both in $\mathrm{Im}\chi_{xx}(\Omega)$ and $\mathrm{Im}\chi_{zz}(\Omega)$, and onsets of continua at $\Omega=2\mu\pm\lambda_{\mathrm{R}}$ in $\mathrm{Im}\chi_{xx}(\Omega)$ and at $\Omega=2\mu$ in $\mathrm{Im}\chi_{zz}(\Omega)$.

To understand the selection rules for transitions induced by an in-plane ac electric field in an EDSR measurement, we need the matrix elements of the velocity operator, which coincides with the velocity operator for Dirac fermions without SOC coupling:

$$\hat{{\bf v}}=\boldsymbol{\nabla}_{\bf k}\hat{H}_{0}=v_{F}\tau_{z}\hat{s}_{0}\hat{\sigma}_{x}\hat{x}+v_{F}\hat{s}_{0}\hat{\sigma}_{y}\hat{y}.$$

(6)

Using the eigenstates from Eq. (3), we obtain

	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{v}_{x}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle\tau_{z}v_{F}\frac{\epsilon_{\alpha\beta}e^{-i\theta\tau_{z}}+\epsilon_{\alpha^{\prime}\beta^{\prime}}e^{i\theta\tau_{z}}+\alpha\alpha^{\prime}\beta\beta^{\prime}\left(\epsilon_{\alpha\beta}e^{i\theta\tau_{z}}+\epsilon_{\alpha^{\prime}\beta^{\prime}}e^{-i\theta\tau_{z}}\right)}{2\sqrt{1+\epsilon_{\alpha\beta}^{2}}\sqrt{1+\epsilon_{\alpha^{\prime}\beta^{\prime}}^{2}}},$
	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{v}_{y}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle i\tau_{z}v_{F}\frac{\epsilon_{\alpha\beta}e^{-i\theta\tau_{z}}-\epsilon_{\alpha^{\prime}\beta^{\prime}}e^{i\theta\tau_{z}}-\alpha\alpha^{\prime}\beta\beta^{\prime}\left(\epsilon_{\alpha\beta}e^{i\theta\tau_{z}}-\epsilon_{\alpha^{\prime}\beta^{\prime}}e^{-i\theta\tau_{z}}\right)}{2\sqrt{1+\epsilon_{\alpha\beta}^{2}}\sqrt{1+\epsilon_{\alpha^{\prime}\beta^{\prime}}^{2}}}.$		(7)

In contrast to the case of magnetic driving, we see that the matrix elements of the velocity are finite for any combination of $\alpha$, $\beta$, $\alpha^{\prime}$, and $\beta^{\prime}$, and thus transitions between all the subbands can be excited by an ac electric field (except for those within the valence bands, which are forbidden by the Pauli principle). Therefore, we expect to see a resonance in the conductivity at $\Omega=\lambda_{\mathrm{R}}$ and onsets of continua at $\Omega=2\mu\pm\lambda_{\mathrm{R}}$ and $\Omega=2\mu$, as shown in Fig. 1b. In an EDSR experiment, one measures the real part of the optical conductivity, which can be also calculated using the Kubo formula (see. Appendix A.3.2). The conductivity shown in Fig. 3a (without the Drude part) indeed exhibits all the features following from the selection rules.

In the opposite limiting case, when only VZ SOC is present ($\lambda_{\mathrm{R}}=0$), the eigenvalues and eigenstates of the Hamlitonian (2) are given by

$$\begin{split}\varepsilon_{\alpha\beta;\tau_{z}}({\bf k})&=\alpha\left(v_{F}k+\beta\frac{\lambda_{\rm Z}}{2}\right),\\ \ket{\alpha\beta;\tau_{z}}&=\frac{1}{2\sqrt{2}}\left(\begin{array}[]{c}1+\tau_{z}\alpha\beta\\ (\tau_{z}\alpha+\beta)e^{i\tau_{z}\theta}\\ 1-\tau_{z}\alpha\beta\\ (\tau_{z}\alpha-\beta)e^{i\tau_{z}\theta}\end{array}\right)\end{split}$$

(8)

with the same notations as in Eq. (3). The corresponding spectrum is shown in Fig. 1b.

The expectation values of the spin operators are

$$\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{x}|\alpha\beta;\tau_{z}\rangle=0,~{}~{}\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{y}|\alpha\beta;\tau_{z}\rangle=0,~{}~{}\langle\alpha\beta;\tau_{z}|\hat{\mathcal{S}}_{z}|\alpha\beta;\tau_{z}\rangle=\tau_{z}\alpha\beta.$$

(9)

In contrast to the case of Rashba SOC, spins are no longer chiral but Ising-like, and are polarized in the opposite directions in the SOC-split bands, see Fig. 1d. Within a given electron or hole spin-split subband, the direction of polarization is also opposite at the $K$ and $K^{\prime}$ points. This is expected as the system preserves time-reversal symmetry.

To understand the selection rules for transitions probed by ESR, we need the off-diagonal matrix elements of the spin operator, which are given by

$$\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{x}|\alpha\beta;\tau_{z}\rangle=-\delta_{\alpha\alpha^{\prime}}\delta_{\beta,-\beta^{\prime}},~{}~{}\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{y}|\alpha\beta;\tau_{z}\rangle=i\tau_{z}\alpha\beta\delta_{\alpha\alpha^{\prime}}\delta_{\beta,-\beta^{\prime}},~{}~{}\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{z}|\alpha\beta;\tau_{z}\rangle=\tau_{z}\alpha\beta\delta_{\alpha\alpha^{\prime}}\delta_{\beta\beta^{\prime}}.$$

(10)

That $\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{\mathcal{S}}_{z}|\alpha\beta;\tau_{z}\rangle$ is purely diagonal implies that an out-of-plane magnetic field cannot induce any inter-band transitions; this is so because $\hat{\mathcal{S}}_{z}$ commutes with the Hamiltonian. Furthermore, an in-plane magnetic field can only induce a transition at $\Omega=\lambda_{\mathrm{Z}}$ between the spin-split branches of the conduction band, but there are no continua of spin-flip transitions, see Fig. 1d. Correspondingly, $\mathrm{Im}\chi_{zz}(\Omega)=0$ while $\mathrm{Im}\chi_{xx}(\Omega)$ exhibits a single peak at $\Omega=\lambda_{\mathrm{Z}}$, as shown in Fig. 2c.

The selections rules for transitions probed by an ac electric field follow from the matrix elements of velocity, given by

	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{v}_{x}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle\frac{\tau_{z}v_{F}}{4}\left[\alpha e^{i\theta\tau_{z}}\alpha^{\prime}+e^{-i\theta\tau_{z}}\right](1+\alpha\alpha^{\prime}\beta\beta^{\prime}),$
	$\displaystyle\langle\alpha^{\prime}\beta^{\prime};\tau_{z}|\hat{v}_{y}|\alpha\beta;\tau_{z}\rangle$	$\displaystyle=$	$\displaystyle\frac{\tau_{z}v_{F}}{4i}\left[\alpha e^{i\theta\tau_{z}}-\alpha^{\prime}e^{-i\theta\tau_{z}}\right](1+\alpha\alpha^{\prime}\beta\beta^{\prime}).$		(11)

It follows from the last equation that an in-plane ac electric field can induce transitions only between the valence and conduction band, with a simultaneous flip of chirality. As shown in Fig. 1e, the frequencies of these transitions are $\Omega=2\mu\pm\lambda_{\mathrm{Z}}$, but there is no resonance peak at $\Omega=\lambda_{\mathrm{Z}}$. Indeed, the corresponding conductivity in Fig. 3b shows no resonance but only a double step at $\Omega=2\mu\pm\lambda_{\mathrm{Z}}$, which is just a Pauli threshold at $\Omega=2\mu$ split by VZ SOC.

As we saw in the previous sections, there are two groups of transitions that can be excited by external electric and magnetic fields: i) those between the spin-split subbands of the conduction band at $\Omega=\lambda_{\mathrm{R}},\lambda_{\mathrm{Z}}$ and ii) those between the spin-split conduction and valence bands at $\Omega=2\mu,2\mu\pm\lambda_{\mathrm{R}},2\mu\pm\lambda_{\mathrm{Z}}$. Given that $\lambda_{\mathrm{R}},\lambda_{\mathrm{Z}}\ll\mu$ in real systems, the transition frequencies in the second group are near the direct absorption threshold at $\Omega\approx 2\mu$. If the Coulomb interaction between optically excited holes and conduction electrons is accounted for, absorption starts at the indirect threshold of $\Omega=\mu$, via the same mechanism as in Auger damping of photoexcited carriers in doped semiconductors.Gavoret et al. (1969); Pimenov et al. (2017); Goyal et al. In addition, the interaction between electrons in the conduction band also leads to absorption for $\Omega<2\mu$,Sharma et al. (2021) and this contribution is comparable to Auger’s one for $\Omega\sim\mu$. For $\Omega\gtrsim\mu$, the linewidth due to both types of damping is large, on the order of $g^{2}\mu$, where $g$ is the dimensionless coupling constant of the Coulomb interaction. Therefore, for moderate and strong interaction, both Auger and intraband damping are expected to smear the fine features near $2\mu$, induced by SOC. For this reason, we will ignore transitions at $\Omega\approx 2\mu$ and focus on the low-energy part of the spectrum at $\Omega\approx\lambda_{\mathrm{R}},\lambda_{\mathrm{Z}}\ll\mu$.

For this range of frequencies, it makes sense to derive an effective low-energy Hamiltonian for the spin-split conduction band. We start by transforming Eq. (2), written in the sub-lattice basis $\{\hat{a}_{\uparrow},\hat{b}_{\uparrow},\hat{a}_{\downarrow},\hat{b}_{\downarrow}\}^{T}$, to the subband basis $\{\hat{c}_{\uparrow},\hat{c}_{\downarrow},\hat{v}_{\uparrow},\hat{v}_{\downarrow}\}^{T}$, in which the Dirac part of $\hat{H}_{0}$ is diagonal. The transformation is effected via

$\displaystyle\hat{a}_{\varsigma}=\frac{\hat{c}_{\varsigma}+\hat{v}_{\varsigma}}{\sqrt{2}},\;b_{\varsigma}=\frac{\hat{c}_{\varsigma}-\hat{v}_{\varsigma}}{\sqrt{2}}\tau_{z}e^{i\tau_{z}\theta},\;\varsigma={\uparrow,\downarrow}.$

(12)

An ac magnetic field $\hat{\bf{b}}_{0}B_{0}e^{-i\Omega t}$ is accounted for by adding the Zeeman term $\hat{\bf s}\cdot\hat{\bf b}_{0}\hat{\sigma}_{0}\Delta_{\mathrm{Z}}/2$ to the Hamiltonian (2), where $\Delta_{\mathrm{Z}}=g\mu_{B}B_{0}e^{-i\Omega t}$, $g$ is the effective Landé-factor, and $\mu_{B}$ is the Bohr magneton. An ac electric field ${\bf E}_{0}e^{-i\Omega t}$ is accounted for by a gauge transformation ${\bf k}\to{\bf k}+(e/c){\bf A}$ in Eq. (2), where ${\bf A}=(c/i\Omega){\bf E}_{0}e^{-i\Omega t}$ is the vector-potential. In the new basis, and in the presence of both magnetic and electric fields, the Hamiltonian can be written as a $4\times 4$ block-matrix

$\displaystyle\hat{H}_{0}^{\rm band}=\left(\begin{array}[]{cc}\hat{H}_{cc}&\hat{H}_{cv}\\ \hat{H}_{vc}&H_{vv}\end{array}\right),$

(15)

where the $2\times 2$ blocks are given by

	$\displaystyle\hat{H}_{cc}$	$\displaystyle=$	$\displaystyle\hat{s}_{0}v_{F}k+\frac{\lambda_{\mathrm{R}}}{2}(\hat{\bf k}\times\hat{\bf s})\cdot\hat{\bf z}+\frac{\lambda_{\mathrm{Z}}}{2}\tau_{z}\hat{s}_{z}+\frac{\Delta_{\mathrm{Z}}}{2}\hat{\bf{s}}\cdot\hat{\bf{b}}_{0}+\hat{s}_{0}\frac{ev_{F}}{c}(\hat{\bf k}\cdot{\bf A}),$
	$\displaystyle\hat{H}_{vv}$	$\displaystyle=$	$\displaystyle-\hat{s}_{0}v_{F}k-\frac{\lambda_{\mathrm{R}}}{2}(\hat{\bf k}\times\hat{\bf s})\cdot\hat{\bf z}+\frac{\lambda_{\mathrm{Z}}}{2}\tau_{z}\hat{s}_{z}+\frac{\Delta_{\mathrm{Z}}}{2}\hat{\bf{s}}\cdot\hat{\bf{b}}_{0}-\hat{s}_{0}\frac{ev_{F}}{c}(\hat{\bf k}\cdot{\bf A}),$
	$\displaystyle\hat{H}_{cv}$	$\displaystyle=$	$\displaystyle i\tau_{z}\hat{s}_{0}\frac{ev_{F}}{c}(\hat{\bf k}\times{\bf A})\cdot\hat{\bf z}-i\frac{\lambda_{\mathrm{R}}}{2}\tau_{z}(\hat{\bf k}\cdot\hat{\bf s})$
	$\displaystyle\hat{H}_{vc}$	$\displaystyle=$	$\displaystyle-\hat{H}_{cv},$		(16)

with $\hat{\bf k}={\bf k}/k$ and $\hat{\bf z}$ being the unit vector normal the plane. Note that the coupling between the conduction and valence bands arises via the electric field and Rashba SOC. This is due to the selection rules discussed in Secs. II.2 and II.3, which say that, at $\Omega\ll\mu$, only Rashba SOC allows for transitions between the conduction and valence bands. Note also that the Rashba terms in $\hat{H}_{cc}$ and $\hat{H}_{vv}$ are almost the same as in a 2D electron gas (up to the dependence on the magnitude of ${\bf k}$).

Next, we project out the valence band via a standard downfolding procedure. Namely, we find the Green’s function of $\hat{H}_{0}^{\rm band}$ in Eq. (15), $G=(\epsilon-\hat{H}_{0}^{\rm band})^{-1}$, and read off its $cc$ element. This yields the effective low-energy Hamiltonian as $\hat{\mathcal{H}}^{\rm proj}_{cc}=\hat{H}_{cc}+\hat{H}_{cv}(\epsilon-\hat{H}_{vv})^{-1}\hat{H}_{vc}$. To leading order in the external fields and SOC, the eigenvalue $\epsilon$ can be replaced by $v_{F}k$, leading to

$$\hat{\mathcal{H}}^{\rm proj}_{cc}=\hat{s}_{0}v_{F}k+\frac{\lambda_{\mathrm{R}}}{2}(\hat{\bf k}\times\hat{\bf s})\cdot\hat{\bf z}+\frac{\lambda_{\mathrm{Z}}}{2}\tau_{z}\hat{s}_{z}+\frac{\Delta_{\mathrm{Z}}}{2}\hat{\bf{s}}\cdot\hat{\bf{b}}_{0}+\hat{s}_{0}\frac{ev_{F}}{c}(\hat{\bf k}\cdot{\bf A})+\frac{e\lambda_{\mathrm{R}}}{2ck}({\bf A}\times\hat{\bf k})\cdot\hat{\bf z}(\hat{\bf k}\cdot\hat{\bf s}).$$

(17)

Note that the electric field couples to spins only due to Rashba SOC, which is the origin of the EDSR effect. In the absence of external fields, SOC splits the conduction band into two subbands with energies

$$\varepsilon_{\pm}=v_{F}k\pm\frac{1}{2}\lambda_{\text{SOC}},\;\text{where}\;\lambda_{\text{SOC}}=\sqrt{\lambda_{\mathrm{R}}^{2}+\lambda_{\mathrm{Z}}^{2}},$$

(18)

and the resonance occurs at $\Omega=\lambda_{\text{SOC}}$.

In this section, we investigate the effect of electron-electron interaction on ESR and EDSR. As we shall demonstrate, the most prominent effect of the interaction is to split the resonances into two. This splitting is controlled by the coupling constants of the various interaction channels, which enables one to extract these important parameters from the measured spectra. Since we are interested only in energies much smaller than the Fermi energy, the effect of electron-electron interaction can be accounted for within a Fermi liquid (FL) theory for conduction electrons only, while the interaction with holes can be assumed to be absorbed into the coupling constants of the FL theory.

First, we discuss the structure of the FL theory for a two-valley system in the absence of SOC, developed recently in Ref. Raines et al., 2021. The interaction vertices are shown in Fig. 4. The solid and dashed lines depict electrons in the $K$ and $K^{\prime}$ valleys, respectively. Diagrams $a$ and $b$ describe intra-valley scattering. Each of these diagrams has an exchange partner (not shown), in which the outgoing states are swapped. Diagram $c$ describes an inter-valley scattering event, in which electrons stay in their respective valleys. The momentum transfer in such an event is less than or equal to $2k_{F}$. Diagram $d$ is an exchange partner to diagram $c$, in which electrons are swapped between the valleys. The momentum transfer in such an event is close to the distance between the $K$ and $K^{\prime}$ points, $|{\bf K}-{\bf K^{\prime}}|\sim 1/a$, where $a$ is the lattice constant. For $k_{F}a\ll 1$, the matrix element of the Coulomb interaction for diagram $d$ is much smaller than that for diagram $c$, and will be neglected in what follows. In this case, the valley index plays a role of conserved isospin, and we have an $SU(2)\times SU(2)$-invariant FL. The interaction between quasiparticles of such a FL is described by the Landau interaction function

	$\displaystyle\nu_{F}^{*}\hat{f}({\bf k},{\bf k}^{\prime})$	$\displaystyle=$	$\displaystyle F^{s}({\bf k},{\bf k}^{\prime})\delta_{\upsilon_{1}\upsilon_{3}}\delta_{\upsilon_{2}\upsilon_{4}}\delta_{\varsigma_{1}\varsigma_{3}}\delta_{\varsigma_{2}\varsigma_{4}}+F^{a}({\bf k},{\bf k}^{\prime})\delta_{\upsilon_{1}\upsilon_{3}}\delta_{\upsilon_{2}\upsilon_{4}}(\hat{\bf{s}}_{\varsigma_{1}\varsigma_{3}}\cdot\hat{\bf{s}}_{\varsigma_{2}\varsigma_{4}})$		(19)
		$\displaystyle+$	$\displaystyle G^{a}({\bf k},{\bf k}^{\prime})(\hat{\boldsymbol{\tau}}_{\upsilon_{1}\upsilon_{3}}\cdot\hat{\boldsymbol{\tau}}_{\upsilon_{2}\upsilon_{4}})\delta_{\varsigma_{1}\varsigma_{3}}\delta_{\varsigma_{2}\varsigma_{4}}+H({\bf k},{\bf k}^{\prime})(\hat{\boldsymbol{\tau}}_{\upsilon_{1}\upsilon_{3}}\cdot\hat{\boldsymbol{\tau}}_{\upsilon_{2}\upsilon_{4}})(\hat{\bf{s}}_{\varsigma_{1}\varsigma_{3}}\cdot\hat{\bf{s}}_{\varsigma_{2}\varsigma_{4}}),$

where $\nu_{F}^{*}$ is the renormalized density of states at the Fermi energy and $v$ ($\varsigma$) labels valley (spin). Components $F^{s}$ and $F^{a}$, which are present also in the single-valley case, describe direct and exchange interaction between fermions in the same valley, respectively. Component $G^{a}$ describes exchange interaction between different valleys, which would be present even for spinless fermions. Finally, component $H$ describes exchange interaction between both spins and valleys. The components of the Landau function are defined on the Fermi surface, i.e., for $k=k^{\prime}=k_{F}$, depend only on angle $\vartheta$ between ${\bf k}$ and ${\bf k}^{\prime}$, and can be characterized by angular harmonics, e.g., $F^{a}_{m}=\int d\vartheta F^{a}(\vartheta)e^{im\vartheta}/2\pi$, and the same for other components.

As long as SOC can be treated as a perturbation, i.e., for $\lambda_{\mathrm{R}},\lambda_{\mathrm{Z}}\ll\mu$, the Landau function in Eq. (19) can also be used to describe dynamics of spins in the presence of SOC, similar to how it was done in Refs. Shekhter et al., 2005; Ashrafi et al., 2013; Kumar and Maslov, 2017 for a 2D electron gas with Rashba and Dresselhaus SOC.

Accounting for the valley degree of freedom and in the absence of external fields, one can write the projected Hamiltonian (17) becomes

$$\hat{\mathcal{H}}^{\rm proj}_{cc}=v_{F}k\hat{\tau}_{0}\hat{s}_{0}+\frac{\lambda_{\mathrm{R}}}{2}\hat{\tau}_{0}(\hat{\bf k}\times\hat{\bf s})\cdot\hat{\bf z}+\frac{\lambda_{Z}}{2}\hat{\tau}_{z}\hat{s}_{z}.$$

(20)

Accordingly, the quasiparticle energy can be written as the sum of the equilibrium part (eq) and a correction due to the Landau functional (LF), whereas the equilibrium part is further separated into a spin-independent part and a correction due to SOC:

	$\displaystyle\hat{\varepsilon}({\bf k},t)$	$\displaystyle=$	$\displaystyle\hat{\varepsilon}_{\text{eq}}({\bf k})+\delta\hat{\varepsilon}_{\mathrm{LF}}({\bf k},t),$		(21a)
	$\displaystyle\hat{\varepsilon}_{\text{eq}}({\bf k})$	$\displaystyle=$	$\displaystyle\hat{\tau}_{0}\hat{s}_{0}v^{*}_{F}(k-k_{F})+\delta\hat{\varepsilon}_{\text{SO}}({\bf k}),$		(21b)
	$\displaystyle\delta\hat{\varepsilon}_{\text{SO}}({\bf k})$	$\displaystyle=$	$\displaystyle\frac{\lambda_{\mathrm{R}}^{*}}{2}\hat{\tau}_{0}(\hat{\bf k}\times\hat{\bf s})\cdot\hat{\bf z}+\frac{\lambda_{\mathrm{Z}}^{*}}{2}\hat{\tau}_{z}\hat{s}_{z},$		(21c)
	$\displaystyle\delta\hat{\varepsilon}_{\mathrm{LF}}({\bf k},t)$	$\displaystyle=$	$\displaystyle\text{Tr}^{\prime}\int\frac{d^{2}p^{\prime}}{(2\pi)^{2}}\hat{f}({\bf k},{\bf k}^{\prime})\hat{n}({\bf k}^{\prime},t).$		(21d)

Here, $\delta\hat{n}({\bf k},t)$ is the density matrix, $f({\bf k},{\bf k}^{\prime})$ is given by Eq. (19), ^′ indicates the valley/spin state of a quasiparticle with momentum ${\bf k}^{\prime}$, $v_{F}^{*}$ is the renormalized Fermi velocity, and $\lambda_{\mathrm{R}}^{*}=\lambda_{\mathrm{R}}/(1+F_{1}^{a})$,Shekhter et al. (2005) $\lambda_{\mathrm{Z}}^{*}=\lambda_{\mathrm{Z}}/(1+H_{0})$Raines et al. are the renormalized spin-orbit coupling constants. Equations (21a)-(21d) need to be solved self-consistently along with the kinetic equation for the density matrix

$$i\frac{\partial\hat{n}({\bf k},t)}{\partial t}=[\hat{\varepsilon}({\bf k},t),\hat{n}({\bf k},t)].$$

(22)

[Our kinetic equation does not contain the effects of Berry curvature because we neglected the asymmetry gap in the single-particle spectrum.] As in Refs. Shekhter et al., 2005; Kumar and Maslov, 2017, we introduce a set of rotated Pauli matrices in the spin space

$$\hat{\xi}_{0}=\hat{s}_{0},\,\,\,\,\,\,\hat{\xi}_{1}({\bf k})=-\hat{s}_{z},\,\,\,\,\,\,\hat{\xi}_{2}({\bf k})=\cos\theta\hat{s}_{x}+\sin\theta\hat{s}_{y},\,\,\,\,\,\,\hat{\xi}_{3}({\bf k})=\sin\theta\hat{s}_{x}-\cos\theta\hat{s}_{y},$$

(23)

and parametrize $\delta\hat{n}({\bf p,t})$ as the sum of the equilibrium (eq) and fluctuating (fl) parts: Raines et al. (2021, )

	$\displaystyle\hat{n}({\bf k},t)$	$\displaystyle=$	$\displaystyle\hat{n}_{\text{eq}}({\bf k})+\delta\hat{n}_{\text{fl}}({\bf k},t),$		(24a)
	$\displaystyle n_{\text{eq}}({\bf k})$	$\displaystyle=$	$\displaystyle\hat{\tau}_{0}\hat{\xi}_{0}n_{F}+n_{F}^{\prime}\hat{\varepsilon}_{\text{SO}}({\bf k}),$		(24b)
	$\displaystyle\delta\hat{n}_{\text{fl}}({\bf k},t)$	$\displaystyle=$	$\displaystyle n_{F}^{\prime}\hat{\tau}_{0}\hat{\xi}_{0}a({\bf k},t)+\delta\hat{n}_{\text{sv}}({\bf k},t),$		(24c)
	$\displaystyle\delta{\hat{n}_{\text{sv}}}({\bf k},t)$	$\displaystyle=$	$\displaystyle n_{F}^{\prime}\left[\hat{\tau}_{0}{\bf u}({\bf k},t)\cdot\hat{\boldsymbol{\xi}}({\bf k})+{\bf w}({\bf k},t)\cdot\hat{\boldsymbol{\tau}}\hat{\xi}_{0}+M_{\alpha\beta}({\bf k},t)\hat{\tau}_{\alpha}\hat{\xi}_{\beta}({\bf k})\right],$		(24d)

where $n_{F}$ is the Fermi function, $\alpha,\beta\in\{1,2,3\}$, and $n_{F}^{\prime}\equiv\partial_{\varepsilon}n_{F}(\varepsilon)$. Vector ${\bf u}$ and tensor $M_{\alpha\beta}$ describe oscillations of the uniform magnetization

$$\tilde{S}_{\alpha}=-\frac{g\mu_{B}}{2}\mathcal{S}_{\alpha}=-\frac{g\mu_{B}}{2}\int\frac{d^{2}k}{(2\pi)^{2}}\text{Tr}\left[\delta\hat{n}({\bf k},t)\,\hat{s}_{\alpha}\right]=\frac{g\mu_{B}\nu_{F}^{*}}{8}\int\frac{d\theta_{\bf k}}{2\pi}u_{\beta}({\bf k},t)\text{Tr}\left[\hat{\xi}_{\beta}\hat{s}_{\alpha}\right]$$

(25)

and valley-staggered magnetization

$\displaystyle\tilde{M}_{\alpha}=-\frac{g\mu_{B}}{2}\mathcal{M}_{\alpha}=-\frac{g\mu_{B}}{2}\int\frac{d^{2}k}{(2\pi)^{2}}\text{Tr}\left[\delta\hat{n}({\bf k},t)\,\tau_{z}\hat{s}_{\alpha}\right]=\frac{g\mu_{B}\nu_{F}^{*}}{8}\int\frac{d\theta_{\bf k}}{2\pi}\mathcal{N}_{\beta}({\bf k},t)\text{Tr}\left[\hat{\xi}_{\beta}\hat{s}_{\alpha}\right],$

(26)

respectively, where $\mathcal{N}_{\gamma}({\bf k},t)\equiv M_{3\gamma}({\bf k},t)$ with $\gamma=1\dots 3$. Vector ${\bf w}$ describes oscillations in the valley occupancy, which are decoupled from both magnetizations and will not be considered below. Expanding ${\bf u}$ and $\boldsymbol{\mathcal{N}}$ over the set of angular harmonics as ${\bf u}({\bf k},t)=\sum_{m}e^{im\theta}{\bf u}^{(m)}(t)$ and $\boldsymbol{\mathcal{N}}({\bf k},t)=\sum_{m}e^{im\theta}\boldsymbol{\mathcal{N}}^{(m)}(t)$, we obtain for the components of $\boldsymbol{\mathcal{S}}$

$\displaystyle\mathcal{S}_{x}=\nu_{F}^{*}\left[\frac{u_{2}^{(-1)}+u_{2}^{(+1)}}{2}+\frac{u_{3}^{(-1)}-u_{3}^{(+1)}}{2i}\right],\;\mathcal{S}_{y}=\nu_{F}^{*}\left[\frac{u_{2}^{(-1)}-u_{2}^{(+1)}}{2i}-\frac{u_{3}^{(-1)}+u_{3}^{(+1)}}{2}\right],\;\mathcal{S}_{z}=-\nu_{F}^{*}u_{1}^{(0)}.$

(27)

The expressions for $\mathcal{M}_{i}$ are obtained from the last equation by replacing $u_{\gamma}^{(\pm 1)}\to\mathcal{N}_{\gamma}^{(\pm 1)}$.

The equations of motion for $u^{(m)}_{\alpha}(t)$ and $\mathcal{N}^{(m)}_{\alpha}(t)$ are obtained by tracing out the corresponding components of Eq. (22). The full system of equations is presented in Appendix (B). To understand the dynamics of the system, it is instructive to consider first the case of $\lambda_{Z}^{*}=0$, when the equations for $u^{(m)}_{\alpha}(t)$ and $\mathcal{N}^{(m)}_{\alpha}(t)$ decouple. The $m=\pm 1$ harmonics, which enter the in-plane components of $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{M}}$, satisfy

	$\displaystyle\dot{u}_{1}^{(\pm 1)}$	$\displaystyle=$	$\displaystyle\lambda_{\mathrm{R}}^{*}\left[f_{+}u_{2}^{(\pm 1)}\pm if_{-}u_{3}^{(\pm 1)}\right],\;\dot{u}_{2}^{(\pm 1)}=-\lambda_{\mathrm{R}}^{*}fu_{1}^{(\pm 1)},\dot{u}_{3}^{(\pm 1)}=0;$		(28a)
	$\displaystyle\dot{\mathcal{N}}_{1}^{(\pm 1)}$	$\displaystyle=$	$\displaystyle\lambda_{\mathrm{R}}^{*}\left[h_{+}\mathcal{N}_{2}^{(\pm 1)}\pm ih_{-}\mathcal{N}_{3}^{(\pm 1)}\right],\;\dot{\mathcal{N}}_{2}^{(\pm 1)}=-\lambda_{\mathrm{R}}^{*}h\mathcal{N}_{1}^{(\pm 1)},\;\dot{\mathcal{N}}_{3}^{(\pm 1)}=0,$		(28b)

where

	$\displaystyle f$	$\displaystyle=$	$\displaystyle 1+F^{a}_{1},\;f_{+}=1+(F^{a}_{0}+F^{a}_{2})/2,\;f_{-}=(F^{a}_{0}-F^{a}_{2})/2,$
	$\displaystyle h$	$\displaystyle=$	$\displaystyle 1+H_{1},\;h_{+}=1+(H_{0}+H_{2})/2,\;h_{-}=(H_{0}-H_{2})/2.$		(29)

Equations (28a) and (28b) describe independent oscillations of the in-plane magnetization and valley-staggered magnetization with frequencies $\Omega_{-}=\lambda_{\mathrm{R}}^{*}\sqrt{ff_{+}}$ and $\Omega_{+}=\lambda_{\mathrm{R}}^{*}\sqrt{hh_{+}}$, respectively. From the structure of eigenvectors, one can deduce that both modes are linearly polarized.

The out-of-plane components of $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{M}}$ are expressed via the $m=0$ harmonics of ${\bf u}$ and $\boldsymbol{\mathcal{N}}$, which satisfy a set of equations similar to Eqs. (28a) and (28b) with different combinations of the Landau parameters [see Appendix (B)]. The mode frequencies in the $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{M}}$ sectors are $\Omega_{-,z}=\lambda_{\mathrm{R}}^{*}\sqrt{f(f_{+}+f_{-})}$ and $\Omega_{+,z}=\lambda_{\mathrm{R}}^{*}\sqrt{h(h_{+}+h_{-})}$, respectively.

Valley-Zeeman SOC mixes up the $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{M}}$ sectors. The frequencies of the in-plane modes (with $m=\pm 1$) for $\lambda_{\mathrm{Z}}^{*}\neq 0$ are given by

	$$\begin{split}\Omega_{\pm}^{2}=&\lambda_{\mathrm{R}}^{*2}\bigg{(}\frac{ff_{+}+hh_{+}}{2}\bigg{)}+\lambda_{\mathrm{Z}}^{*2}\bigg{(}f_{+}h_{+}+f_{-}h_{-}\bigg{)}\pm\Omega_{0}^{2},\end{split}$$		(30a)
where
	$$\begin{split}\Omega_{0}^{2}=&\left[\lambda_{\mathrm{R}}^{*4}\bigg{(}\frac{ff_{+}-hh_{+}}{2}\bigg{)}^{2}+\lambda_{\mathrm{Z}}^{*4}(f_{-}h_{+}+h_{-}f_{+})^{2}+\lambda_{\mathrm{R}}^{*2}\lambda_{\mathrm{Z}}^{*2}(ff_{-}+hh_{-})(f_{-}h_{+}+h_{-}f_{+})\right]^{1/2}.\end{split}$$		(30b)

In assigning the $\pm$ indices to the modes, we assumed that $ff_{+}<hh_{+}$ and $f_{-}h_{+}+h-_{-}f_{+}<0$. These two conditions are satisfied in the most realistic case, when i) both the intra- and intervalley interactions are attractive, i.e., $F^{a}_{m},H_{m}<0$; ii) the intravalley interaction is stronger than the intervalley one, i.e., $|F^{a}_{m}|>|H_{m}|$; and iii) both $|F^{a}_{m}|$ and $|H_{m}|$ decrease with $m$ monotonically. In this case, $\Omega_{+}>\Omega_{-}$ for any ratio of $\lambda_{\mathrm{R}}^{*}$ to $\lambda_{\mathrm{Z}}^{*}$. The structure of the eigenstates indicates that for both modes $\boldsymbol{\mathcal{S}}$ and $\boldsymbol{\mathcal{M}}$ are perpendicular to each other and phase lagged by $\pi/2$.