1 Introduction↩︎

Magnetic quantum oscillations, such as the Shubnikov-de Haas (SdH) and de Haas-van Alphen (dHvA) effects, are among the most powerful probes of electronic structure in metals and semimetals [1]–[10]. They arise from the successive crossing of Landau levels through the Fermi energy as the magnetic field is varied. These oscillations are typically periodic in inverse magnetic field, with a frequency determined by the extremal Fermi-surface orbit area [4], [5], [11], [12]. As a result, quantum oscillations form a cornerstone of modern Fermiology [8]. While the conventional picture involves a single oscillation frequency, several mechanisms can modify this behavior. Examples include aperiodic oscillations [13]–[19], interaction-induced anomalous oscillations [20]–[23], and amplitude beating between nearby frequencies [24]–[27]. Beating occurs when two oscillatory components with slightly different frequencies or phases interfere, producing a slowly varying envelope with nodes where the oscillation amplitude vanishes. Since the frequency is directly related to the Fermi-surface area, the beating pattern contains information about the origin of the frequency mismatch.

In graphene-based systems, several distinct mechanisms can generate such a mismatch in quantum oscillation frequencies. A strain-induced pseudomagnetic field makes the two valleys experience different effective magnetic fields, and hence different oscillation frequencies [28]–[36]. Alternatively, the two oscillatory channels can acquire different Fermi-surface areas through unequal valley populations, valley-dependent energy shifts, spin-orbit coupling-induced band splittings, or Kekulé distortions [37]–[44]. Although these mechanisms have different microscopic origins, they can produce qualitatively similar beating envelopes, making it difficult to identify the microscopic mechanism giving rise to the beating from the oscillation pattern alone.

This raises a natural question: how can the SdH beating spectrum be used to distinguish between different microscopic mechanisms that give rise to the beating pattern? A recent study showed that in strained monolayer graphene, the carrier density at the beating nodes scales quadratically with magnetic field, enabling extraction of the pseudomagnetic field from bulk transport measurements [31]. This demonstrates that the node positions contain information beyond the visual appearance of the beating envelope and suggests a more general diagnostic framework based on the scaling of the node trajectories.

In this work, we develop such a framework using Onsager’s quantization relation. We derive general scaling relations for the critical carrier density \(N_{c,j}\) and magnetic field \(B_{c,j}\) corresponding to the \(j\)th beating node, where \(j=0,1,2,\ldots\) labels successive nodes. Figure 1 shows two representative cases in the carrier density and magnetic field plane. A strain-induced pseudomagnetic field gives \(N_{c,j}\propto (2j+1)B_{c,j}^{2}\), whereas a density-independent valley imbalance gives \(N_{c,j}\propto (2j+1)B_{c,j}\). A constant valley- or spin-dependent energy splitting also produces a quadratic magnetic field dependence, but with a distinct node-index dependence, \(N_{c,j}\propto (2j+1)^2B_{c,j}^{2}\). Thus, the magnetic field dependence provides the first diagnostic, while the node-index dependence separates mechanisms that otherwise have the same field scaling. Table 1 summarizes the resulting diagnostics for pseudomagnetic fields, valley imbalance, intrinsic valley splitting, valley-Zeeman spin-orbit coupling, and Y-type Kekulé distortion. Our results establish that, when the beating is governed by two oscillatory channels, the beating-node trajectories provide a quantitative probe of pseudomagnetic fields and valley- or spin-dependent electronic structure in graphene-based systems.

2 Scaling relations from Onsager’s quantization condition↩︎

Since beating arises from the interference between two closely spaced quantum oscillation channels, Onsager’s quantization condition provides a natural starting point for relating the node positions to carrier density, magnetic field, and node index. For a nondegenerate oscillatory channel with fixed spin and valley indices, the semiclassical Onsager quantization condition reads \[\left( n+\gamma \right)\frac{eB}{h} =\frac{S_0(\varepsilon)}{4\pi^2} \equiv N_0(\varepsilon)~. \label{onsager}\tag{1}\] Here, \(n\) is the Landau level (LL) index, \(S_0(\varepsilon)\) is the extremal area of the constant-energy orbit, \(\gamma\) is the phase offset, and \(N_0(\varepsilon)\) is the zero-field carrier density for one spin and one valley.

For a magneto-oscillatory quantity \(X\), each contribution has the generic form \(\delta X \propto \cos \left(2\pi n + \phi_0\right)\), where \(n\) is fixed by Eq. 1 and \(\phi_0\) is a phase offset [4], [5], [7], [19]. Beating in physical responses arises from two nearby oscillatory contributions characterized by slightly different Landau level indices, \(n_1\) and \(n_2\). Assuming comparable amplitudes and no additional relative phase varying between the two channels, the total response is \[\begin{align} \delta X_{\rm tot} &\propto& \cos \left(2\pi n_1+\phi_0\right) + \cos \left(2\pi n_2+\phi_0\right)\nonumber\\ &=& 2\cos\left(\pi\delta n\right) \cos\left(\pi\bar n+\phi_0\right), \label{beat95form} \end{align}\tag{2}\] where \(\delta n=n_1-n_2\) and \(\bar n=n_1+n_2\). Since \(\delta n \ll \bar n\), the oscillation is modulated by the slowly varying envelope \(\cos(\pi\delta n)\). A beating node occurs when this envelope vanishes, \[|\delta n|=j+\frac{1}{2}, \qquad j=0,1,2,\cdots .\] We denote the applied field and total carrier density at the \(j\)-th node by \(B_{c,j}\) and \(N_{c,j}\), respectively. The corresponding Landau level filling factor is given by \(\nu_{c,j}=hN_{c,j}/(eB_{c,j})\).

The phase mismatch can arise from two physically distinct sources: a field-like asymmetry or a Fermi-surface asymmetry. To make this distinction explicit, consider graphene with two valleys, \(\xi=\pm\), and allow the valleys to have different carrier densities. We denote the spin-degenerate zero-field carrier density in valley \(\xi\) by \(N_\xi\), so that \(N_\xi=2N_0\) and the total carrier density is \(N=N_+ + N_-\). We will use this notation for total density, and spin summed density for each valley, throughout the manuscript. We also allow the valleys to experience different effective magnetic fields, \[B_\xi=B-\xi B_{\rm pm},\] where \(B_{\rm pm}\) is a valley-contrasting pseudomagnetic field (PMF). Such PMFs can arise from nonuniform strain in hexagonal two-dimensional systems [29], [30], [45].

The Onsager condition for the spin-summed density in valley \(\xi\) is \[n_\xi+\gamma =\frac{hN_\xi}{2e(B-\xi B_{\rm pm})} \approx \frac{hN_\xi}{2eB}\left(1+\xi\frac{B_{\rm pm}}{B}\right), \label{gen95onsager}\tag{3}\] valid to leading order in \(B_{\rm pm}/B\), where \(n_\xi\) denotes the LL index for valley \(\xi\). We have assumed that the two channels have the same phase offset \(\gamma\). A constant relative offset shifts the node indexing but does not change the leading field or density powers. Subtracting the two valley quantization conditions gives \[\label{gen95scaling} \delta n =\frac{h}{2eB} \left[ \delta N +N\frac{B_{\rm pm}}{B} \right],\tag{4}\] where \(\delta N=N_+-N_-\) is the valley carrier density imbalance. The two terms in Eq. 4 have different physical origins. The term \(N B_{\rm pm}/B\) is field-like: it appears when the two valleys experience different effective magnetic fields. The term \(\delta N\) is Fermi-surface-like: the two channels have different carrier densities and therefore different SdH frequencies even without a field-like asymmetry.

In this illustration, the Fermi-surface asymmetry comes from differential valley population. More generally, the same density-imbalance term can represent valley-dependent band shifts, spin-orbit-induced spin splitting, velocity imbalance, or multiple bands contributing as distinct oscillatory channels. These realizations are organized in Sec. 5. We now focus on the beating-node scaling when only one asymmetry, either field-like or Fermi-surface-like, is present.

2.1 Field-like asymmetry↩︎

Long-wavelength nonuniform strain generates a valley-contrasting pseudomagnetic field in graphene-based hexagonal systems, which preserves global time-reversal symmetry. As a result, carriers in the two valleys experience effective magnetic fields \(B\pm B_{\rm pm}\). In these systems, there is no zero-field Fermi-surface asymmetry from valley or spin splitting, so \(\delta N=0\). The phase mismatch responsible for beating therefore originates solely from field-like asymmetry. Applying the node condition to Eq. 4 gives \[\label{quadratic95B95scaling} N^{\rm PMF}_{c,j} =\frac{e}{h|B_{\rm pm}|}(2j+1)B_{c,j}^2~.\tag{5}\] Thus, a PMF-induced field-like asymmetry produces a quadratic magnetic field dependence and a linear node-index dependence in the critical density, \(N^{\rm PMF}_{c,j}\propto (2j+1)B_{c,j}^2\). This scaling provides a direct way to extract the magnitude of the pseudomagnetic field from beating-node positions in non-uniformly strained graphene [31].

2.2 Fermi-surface asymmetry↩︎

In the absence of PMF, beating in SdH oscillations originates only from unequal Fermi-surface areas. Then Eq. 4 reduces to \[|\delta N(N_{c,j})| =\frac{e}{h}(2j+1)B_{c,j}~. \label{FS95node95condition}\tag{6}\] The field dependence of the node trajectory is therefore controlled by how the density imbalance depends on the total carrier density. This density, or Fermi-surface, asymmetry can arise from SOC or from valley splitting. The two channels experience the same applied field, but their extremal orbit areas differ. In several cases, the density imbalance itself depends on the total density. For generality, we write \(|\delta N(N)|=A N^r\), where \(A\) and \(r\) are material- and mechanism-specific parameters. Combining this form with Eq. 6 , we obtain \[N_{c,j} =\left[\frac{e}{hA}(2j+1)B_{c,j}\right]^{1/r}. \label{FS95general95scaling}\tag{7}\] Contrasting this density-imbalance case with Eq. 5 , we see that within this power-law class the beating-node density scales with the node index and the magnetic field with the same power.

A useful case is a density imbalance proportional to the total density. This corresponds to an imbalance in the valley populations described by a polarization parameter \(p_v\in[-1,1]\), such that \(N_\xi=(1-\xi p_v)N/2\). The magnitude of the valley population imbalance is then \(|\delta N|=|p_v|N\), corresponding to \(r=1\). If spin degeneracy is preserved, Eq. 7 gives \[\label{linear95B95scaling} N^{\rm VP}_{c,j} =\frac{e}{h|p_v|}(2j+1)B_{c,j}~.\tag{8}\] A concrete example of this constant-polarization scaling is Y-type Kekulé distortion in graphene [31], [46], [47], which produces a Fermi-velocity imbalance between the two inequivalent valleys. See Sec. 5.1 for details.

In contrast, a constant energy splitting of two Dirac valleys gives \(|\delta N|\propto N^{1/2}\). For example, if the two Dirac valleys are shifted by energies \(\pm \Delta_0\), the resulting carrier density imbalance is \(|\delta N|\simeq \dfrac{2|\Delta_0|}{\sqrt{\pi}\hbar v_F}N^{1/2}\), as discussed in Sec. 5.2. Substituting \(r=1/2\) in Eq. 7 yields \[N_{c,j}\propto (2j+1)^2B_{c,j}^2~.\] Although both the PMF-induced and valley-contrasting energy-splitting trajectories are quadratic in \(B_{c,j}\), their node-index factors are different. A PMF gives \((2j+1)B_{c,j}^2\), whereas a Dirac energy splitting gives \((2j+1)^2B_{c,j}^2\). This distinction offers a direct way to identify the microscopic origin of the observed beating in experiments.

2.3 Universal scaling ratio↩︎

The same two-channel interference picture also gives a simple mechanism-independent relation. The rapid SdH oscillations are controlled by the average phase of the two channels, whereas the beating envelope is controlled by their phase difference; see Eq. 2 . The number of fast oscillations between neighboring beating nodes is therefore set by the ratio of the average LL index to the LL-index mismatch.

At the \(j\)-th beating node, the envelope condition is \(|\delta n|=j+\frac{1}{2}\). The local number of full rapid oscillations between neighboring beating nodes is \[\mathcal{N}_{{\rm osc},j} = \frac{\bar n}{2|\delta n|} = \frac{n_1+n_2}{2j+1}. \label{Nosc95universal}\tag{9}\] The factor of two in the denominator appears because the interval between two neighboring envelope nodes corresponds to a half-period of the beat envelope.

For the two-channel beating considered in this work, the filling factor counts the total number of occupied LLs from both channels. With spin degeneracy included, this gives \[\nu_{c,j}=2(n_1+n_2).\] Combining this with Eq. 9 , we obtain the universal ratio \[\frac{\nu_{c,j}}{\mathcal{N}_{{\rm osc},j}} = 2(2j+1). \label{universal95ratio}\tag{10}\] This ratio is universal in the sense that it does not depend on whether the beating originates from a PMF, valley polarization, intrinsic valley splitting, SOC, or Kekulé distortion. The microscopic mechanism determines how \(N_{c,j}\) and \(B_{c,j}\) scale, but these details cancel in Eq. 10 . The ratio, therefore, does not identify the mechanism by itself. Instead, it provides a consistency check that the observed modulation is governed by two nearby oscillatory channels with weak splitting and reliable node indexing in oscillation data. Our microscopic calculations support the same relation, as shown in Table 1.

We can use the above analysis as a reliable diagnostic tool under the following conditions: two dominant oscillatory channels with comparable amplitudes, weak splitting, reliable node indexing, slowly varying damping factors and smooth prefactors, and no large density inhomogeneity. Outside this regime, the node trajectory can still indicate the source of beating, but it should be interpreted more carefully.

3 Quantum oscillations in graphene↩︎

We first develop the reference oscillatory response in monolayer graphene in the absence of any field-induced or Fermi-surface asymmetry. This will serve as a baseline for comparing with the cases with beating nodes later. For a single valley of gapped monolayer graphene, the low-energy Hamiltonian is \[\label{Ham95dirac} {\cal H}^\xi=\hbar v_F(\xi k_x\sigma_x+k_y\sigma_y)+\Delta\sigma_z~,\tag{11}\] where \(\xi=\pm\) labels the valleys [48], \(v_F\) is the Fermi velocity, and \(\Delta\) is a staggered sublattice potential. The finite \(\Delta\) allows us to keep the valley-dependent zeroth Landau level explicit, while the main scaling results below are evaluated in the gapless limit when appropriate. In a uniform magnetic field, we substitute \(\hbar{\boldsymbol{k}}\rightarrow\hbar{\boldsymbol{k}}+e{\boldsymbol{A}}\), with \({\boldsymbol{B}}=\boldsymbol{\nabla}\times{\boldsymbol{A}}\). The Landau levels are \[\begin{align} \label{LL95disp} \varepsilon_n = \begin{cases} \lambda \sqrt{\Delta^2+2n\hbar^2\omega_c^2}~, & n\neq 0~,\\ -\xi\Delta~, & n=0~, \end{cases} \end{align}\tag{12}\] where \(\lambda=\pm1\), \(\omega_c=v_F/l_B\), and \(l_B=\sqrt{\hbar/eB}\). In the gapless limit used for the main scaling discussion, \(\Delta=0\), the two valleys have identical spectra and identical Fermi-surface areas; see Fig. 2 (a).

Disorder broadens the Landau levels and produces a smooth oscillatory density of states. Keeping the first harmonic, the spin-summed DOS in valley \(\xi\) can be written as \[D_\xi(\varepsilon)=D_{\rm LL0}+D_0(\varepsilon) \left[1+2\Omega_D(\varepsilon,B) \cos\left(\frac{2\pi F(\varepsilon)}{B}\right)\right]~. \label{DOS95osc}\tag{13}\] The detailed derivation is given in Appendix 9. Here, \(D_{\rm LL0}=\dfrac{1}{\pi^2 l_B^2 \Gamma_0^2} \sqrt{\Gamma_0^2-\dfrac{(\varepsilon+\Delta)^2}{4}}~\Theta\left(2\Gamma_0-|\varepsilon+\Delta|\right)\) denotes the zeroth-LL contribution. Here, \(\Theta(x)=1\) for \(x>0\) and \(0\) otherwise is the step function, and \(\Gamma_0\) is the impurity-induced LL broadening parameter. The zero-field DOS per spin-degenerate valley is \(D_0(\varepsilon)= |\varepsilon|/( \pi \hbar^2 v_F^2)\). The factor \(\Omega_D(\varepsilon,B)=\exp [- 2 \pi^2 ( \Gamma_0\varepsilon)^2/(\hbar \omega_c)^4 ]\) accounts for damping due to LL broadening. Equation 13 contains a zeroth-LL contribution and the oscillatory contribution from the nonzero Landau levels. The oscillatory contributions are periodic in \(1/B\), with frequency \(F(\varepsilon)=(\varepsilon^2-\Delta^2)/(2e\hbar^2v_F^2)\). Figure 2 (b) shows the corresponding Landau level quantum oscillations in the DOS as a function of energy.

The longitudinal magnetoconductivity has the same oscillatory phase. Starting from the collisional conductivity and including finite-temperature damping, as summarized in Appendix 10, the valley-resolved result is \[\label{sigma95osc95sdh} \begin{align} \sigma_{xx}^\xi(N_\xi,B_\xi) &\approx \sigma_0(N_\xi,B_\xi) \\ &\quad\times \left[ 1+2\Omega_T\Omega_D \cos\left(\frac{2\pi F(N_\xi)}{B_\xi}\right) \right]~, \end{align}\tag{14}\] where \[\begin{align} \sigma_0(N_\xi,B_\xi) &= \frac{e^2}{h}\frac{n_{\rm im}U_0^2}{2k_s^2\Gamma_0} \frac{2N_\xi}{e\hbar v_F^2B_\xi} \\ &\quad\times \frac{\pi\hbar^2v_F^2N_\xi+2\Delta^2}{\sqrt{2\pi\hbar^2v_F^2N_\xi+\Delta^2}}. \end{align}\] Here, \(N_\xi=(\mu^2-\Delta^2)/(2\pi\hbar^2v_F^2)\) is the spin-summed carrier density in valley \(\xi\), and therefore \(F(N_\xi)=hN_\xi/(2e)\). The parameters \(U_0 = e^2/(2 \epsilon_0 \epsilon_r)\), \(n_{\rm im}\), \(k_s\), \(\epsilon_0\), and \(\epsilon_r\) denote the screened impurity strength, impurity density, screening wave vector, vacuum permittivity, and relative permittivity, respectively. The factor \(\Omega_T=(T/T_D)/\sinh(T/T_D)\) is the thermal damping factor, while \(\Omega_D\) is the disorder-broadening factor.

In pristine graphene, \(N_+=N_-=N/2\) and \(B_+=B_-=B\). In this case, the two valley contributions have the same phase \(2\pi F(N/2)/B\) and add without producing a beating pattern, \[\sigma_{xx}(N,B)=\sum_{\xi=\pm}\sigma_{xx}^\xi(N/2,B)~.\] Figure 2 (c) shows the oscillatory part of the longitudinal magnetoconductivity, \(\delta \sigma_{xx}=\sigma_{xx}-2\sigma_0(N/2,B)\) in the \(B\)-\(N\) plane. The alternating red and blue stripes denote successive SdH oscillations. Since the two valley channels have the same oscillation frequency in this reference case, no beating nodes appear. This provides the baseline against which the beating patterns generated by field-like or Fermi-surface asymmetries are compared below.

4 Field-like asymmetry: pseudomagnetic field↩︎

Long-wavelength strain in graphene acts as a valley-dependent gauge field and generates a pseudomagnetic field [28]. For a honeycomb lattice, the strain-induced vector potential is \[\boldsymbol{A}_{\rm pm} =\frac{\hbar\beta\kappa}{2a_0}(u_{xx}-u_{yy},-2u_{xy})~,\] where \(\beta\approx2\)–\(3\) is the Grüneisen parameter [49], [50], \(a_0\) is the lattice constant, and \(\kappa\approx1/3\) describes the elastic response. The corresponding pseudomagnetic field \({\boldsymbol{B}}_{\rm pm}=\boldsymbol{\nabla}\times\boldsymbol{A}_{\rm pm}\) couples with opposite sign to the two valleys. It therefore preserves global time-reversal symmetry, while making the effective magnetic field valley dependent.

For smooth strain profiles whose correlation length exceeds the cyclotron radius, the PMF can be treated as locally uniform along an orbit [31]. For such cases, strain enters the oscillatory response through valley-dependent magnetic field \(B_\xi=B-\xi B_{\rm pm}\), implying valley-dependent Landau quantization; see Fig. 3 (a). Consequently, the low-energy Hamiltonian can be written in terms of the valley-dependent kinetic momentum \(\boldsymbol{\Pi}_\xi=\hbar{\boldsymbol{k}}+e{\boldsymbol{A}}-e\xi{\boldsymbol{A}}_{\rm pm}\), where \(\boldsymbol{A}\) is the vector potential for real magnetic field.

This is a field-like asymmetry: the two valleys have equal carrier densities, \(N_\xi=N/2\), but they are quantized by different effective fields. Thus, the microscopic perturbation enters through \(B_\xi\), not through a zero-field density imbalance. Substituting \(N_\xi=N/2\) and \(B_\xi=B-\xi B_{\rm pm}\) into Eq. 14 , and neglecting the weak \(B_{\rm pm}\) dependence of the damping factors and prefactors, we obtain, \[\begin{align} \sigma_{xx}(N,B) &\approx& 2\sigma_0 \bigg[ 1+\Omega_T\Omega_D \bigg\{ \cos\left(\frac{2\pi F(N/2)}{B+B_{\rm pm}}\right)\nonumber\\ &&+ \cos\left(\frac{2\pi F(N/2)}{B-B_{\rm pm}}\right) \bigg\} \bigg]~. \label{sigma95beat95quad} \end{align}\tag{15}\] Here, \(\sigma_0\), \(\Omega_T\), and \(\Omega_D\) are evaluated at the common density \(N/2\) and applied field \(B\). In the weak-PMF regime, \(|B_{\rm pm}|\ll B\), this becomes \[\begin{align} \sigma_{xx}(N,B) &\approx& 2\sigma_0 \bigg[ 1+2\Omega_T\Omega_D \cos\left(\frac{2\pi F(N/2)}{B}\right)\nonumber\\ &&\times \cos\left(\frac{2\pi F(N/2)B_{\rm pm}}{B^2}\right) \bigg]~. \label{sigma95beat95quad952} \end{align}\tag{16}\] The first cosine is the rapid SdH oscillation, while the second cosine is the PMF-induced envelope, which produces the beating nodes in the SdH oscillations. Fig. 3 (b) shows the resulting nodes in \(\delta\sigma_{xx}\). Since its phase scales as \(1/B^2\), the envelope is periodic in \(1/B^2\), and the nodes are not equally spaced in \(1/B\); see Fig. 3 (b).

The envelope vanishes when \[\frac{2\pi F(N_{c,j}/2)|B_{\rm pm}|}{B_{c,j}^2} =\left(j+\frac{1}{2}\right)\pi .\] Using \(F(N/2)=hN/(4e)\) gives a node trajectory controlled by the magnitude of the PMF, \[\label{n95c95pmf} N^{\rm PMF}_{c,j} =\frac{e}{h|B_{\rm pm}|}(2j+1)B_{c,j}^2~.\tag{17}\] Note that the above expression exactly matches that obtained from Onsager’s quantization relation, see Eq. 5 . Figure 3 (c) shows \(\delta \sigma_{xx}\) as a function of \(B\) and \(N\). Alternating blue and red stripes denote SdH oscillations, while prominent white stripes mark beating nodes with different \(j\) values. The nodes follow \(N_c \propto B_j^2\) in the \(N\)-\(B\) plane. Tracking this trajectory in oscillatory magnetoresistance enables bulk extraction of \(B_{\rm pm}\) and has recently revealed millitesla-scale pseudomagnetic fields in strained graphene [31].

The associated critical filling factor and local oscillation count are \[\nu^{\rm PMF}_{c,j} =\frac{hN^{\rm PMF}_{c,j}}{eB_{c,j}} =\frac{(2j+1)B_{c,j}}{|B_{\rm pm}|}, \qquad \mathcal{N}^{\rm PMF}_{{\rm osc},j} =\frac{B_{c,j}}{2|B_{\rm pm}|}.\] The oscillation count increases linearly with \(B\), consistent with Fig. 3 (b). Hence \(\nu_{c,j}^{\rm PMF}/\mathcal{N}^{\rm PMF}_{{\rm osc},j}=2(2j+1)\), as required by the two-channel ratio in Sec. 2.3. The PMF node condition carries a single power of \((2j+1)\) multiplying \(B_{c,j}^2\). This node-index factor will differ from the Dirac energy-splitting case below, even though both trajectories are quadratic in \(B_{c,j}\).

The same PMF scaling extends to bilayer graphene. Although monolayer graphene has \(N\propto\mu^2\) and bilayer graphene has \(N\propto\mu\), the beating condition is expressed in carrier density, which is always proportional to the Fermi-surface area in two dimensions. Equation 17 therefore remains valid for the bilayer model treated in Appendix 11 to leading order in the PMF.

5 Fermi-surface asymmetry driven beating↩︎

We now turn to beating generated by Fermi-surface asymmetry. In this case, the two oscillatory channels experience the same applied magnetic field, but their extremal Fermi-surface areas differ. Equation 4 then reduces to \[|\delta N(N_{c,j})| = \frac{e}{h}(2j+1)B_{c,j}~. \label{FS95node95condition95main}\tag{18}\] The node trajectory is therefore determined by the density dependence of the channel density imbalance \(|\delta N|\). This section implements the classification developed in Sec. 2.2: mechanisms with \(|\delta N|\propto N\) give a linear \(N_{c,j}\)-\(B_{c,j}\) trajectory, whereas mechanisms with \(|\delta N|\propto N^{1/2}\) give a quadratic trajectory with an extra power of \((2j+1)\).

5.1 Linear-density imbalance↩︎

5.1.1 Constant fractional imbalance↩︎

A constant fractional valley imbalance gives the simplest \(r=1\) Fermi-surface asymmetry. Such an imbalance may be produced by optical pumping, such as circularly polarized excitation, that populates the valleys unequally while leaving the applied magnetic field common to both valleys. We write the spin-summed valley densities as \(N_\xi=(1-\xi p_v)N/2\), where \(p_v\in[-1,1]\) and \(N=(N_+ + N_-)\). This gives \(|\delta N|=|N_+-N_-|=|p_v|N\). In this scenario, the polarization parameter \(p_v\) is independent of carrier density, as appropriate for an externally imposed nonequilibrium valley imbalance; see Fig. 4 (a). Substituting \(N_\xi\) into Eq. 14 , and evaluating the smooth prefactors and damping factors at the average density, gives \[\begin{align} \label{pop95osc} \sigma_{xx}(N,B) &\approx& 2\sigma_0 \bigg[ 1+2\Omega_T\Omega_D \cos\left(\frac{2\pi F(N/2)}{B}\right)\nonumber\\ &&\times \cos\left(\frac{2\pi F(N|p_v|/2)}{B}\right) \bigg]~. \end{align}\tag{19}\] The first cosine is the rapid SdH oscillation and the second is the Fermi-surface envelope. Unlike the PMF envelope in Fig. 3 (b), this envelope in Fig. 4 (b) is periodic in \(1/B\).

The node condition gives \[\label{vp95linear95B95scaling} N^{\rm VP}_{c,j} = \frac{e}{h|p_v|}(2j+1)B_{c,j}~.\tag{20}\] Figure 4 (c) shows the corresponding beating pattern and the linear node trajectory in the \(N\)-\(B\) plane. The associated filling factor and local rapid-oscillation count are \[\nu^{\rm VP}_{c,j} = \frac{2j+1}{|p_v|}, \qquad \mathcal{N}^{\rm VP}_{{\rm osc},j} = \frac{1}{2|p_v|}.\] Both quantities are independent of the magnetic field, while the field-independent oscillation count is also evident from Fig. 4 (b). This is a useful experimental distinction from PMF-driven beating, where \(\nu_{c,j}^{\rm PMF}\) and \(\mathcal{N}^{\rm PMF}_{{\rm osc},j}\) grow linearly with \(B_{c,j}\).

5.1.2 Velocity imbalance↩︎

A Y-type Kekulé distortion gives a second realization of the \(r=1\) case. It changes the Fermi velocity in a valley-dependent way rather than shifting the two Dirac points in energy [51]. Such distortions can arise in strained graphene [46] or graphene grown on copper substrates [47]. The spin-summed valley carrier densities are \(N^{\xi}= \mu^2/(2 \pi \hbar^2 (v_F^{\xi})^2)\), where \(v_F^{\xi}\) is the valley-resolved Fermi velocity. For a small velocity imbalance \(\delta v_F=v_F^- -v_F^+\) with \(|\delta v_F|\ll v_F\), the resulting density imbalance is \(|\delta N| \simeq N|\delta v_F|/(v_F)\). Therefore, the Y-type Kekulé case can be obtained from the constant-polarization result by the replacement \(|p_v|\to|\delta v_F|/v_F\), giving \[N^{\rm YK}_{c,j} = \frac{e v_F}{h|\delta v_F|}(2j+1)B_{c,j}~.\] Thus, constant fractional imbalance and Y-type Kekulé distortion should be diagnosed by the same linear \(B_{c,j}\) and \((2j+1)\) scaling, even though their microscopic origins are different.

5.2 Square-root-density imbalance↩︎

5.2.1 Energy splitting of Dirac Fermi surfaces↩︎

A constant energy splitting of two Dirac cones induced by a magnetic proximity effect gives a different density power. The splitting is fixed in energy, so the induced density imbalance grows only as the Fermi momentum. A minimal valley-dependent model is \[\label{Ham95VP95int} \mathcal{H}^\xi = \hbar v_F(\xi k_x\sigma_x+k_y\sigma_y) + \xi\Delta_v\sigma_0~,\tag{21}\] where \(\Delta_v\) is the valley-dependent energy shift [52], [53]. The spin-summed carrier density in valley \(\xi\) is \(N_\xi=(\mu-\xi\Delta_v)^2/(2\pi\hbar^2v_F^2)\). In the limit of weak valley splitting, \(|\Delta_v|\ll \mu\), a linear expansion yields a carrier density imbalance \[|\delta N| = |N_+-N_-| \simeq \frac{2|\Delta_v|}{\sqrt{\pi}\hbar v_F}\sqrt{N}.\] Thus, intrinsic valley splitting belongs to the \(r=1/2\) Fermi-surface class. Its node trajectory is \[N^{\rm int}_{c,j} = \frac{(2j+1)^2}{16\pi} \frac{e^2v_F^2}{\Delta_v^2} B_{c,j}^2~.\] This trajectory is quadratic in \(B_{c,j}\), like the PMF trajectory, but the node-index dependence is different. A field-like PMF gives \((2j+1)B_{c,j}^2\), whereas a Dirac energy splitting gives \((2j+1)^2B_{c,j}^2\).

5.2.2 Valley-Zeeman spin-orbit coupling↩︎

Valley-Zeeman spin-orbit coupling in graphene gives the same \(r=1/2\) class when \(s_z\) remains an approximately good quantum number. This regime is relevant for graphene proximitized by transition-metal dichalcogenides [43], [54], where the dominant proximity SOC often acts as a valley-dependent Zeeman field [43], [55]. The low-energy Hamiltonian is \[\label{Ham95dirac95SOC} {\cal H}^{\xi,s_z} = \hbar v_F(\xi k_x\sigma_x+k_y\sigma_y) + \xi s_z\Delta_{\rm VZ}\sigma_0~.\tag{22}\] The energy shift depends on \(\eta=\xi s_z=\pm1\). The relevant two oscillatory channels are therefore the two \(\eta\) sectors, not individual spin or valley branches. The \(\eta=+1\) sector contains \((\xi,s_z)=(+,+)\) and \((-,-)\), while \(\eta=-1\) contains \((+,-)\) and \((-,+)\). Each sector remains twofold degenerate, but the two sectors have different Fermi-surface areas.

Including this twofold degeneracy, the density in sector \(\eta\) is \(N^{\rm VZ}_{\eta} = (\mu-\eta\Delta_{\rm VZ})^2/(2\pi\hbar^2v_F^2)\). This has the same structure as constant energy splitting, with \(\Delta_v\) replaced by \(\Delta_{\rm VZ}\). Hence, we have \[N^{\rm VZ}_{c,j} = \frac{(2j+1)^2}{16\pi} \frac{e^2v_F^2}{\Delta_{\rm VZ}^2} B_{c,j}^2~.\] This mapping assumes that spin remains a good quantum number, i.e., spin mixing is negligible. Strong Rashba SOC mixes spin and sublattice degrees of freedom and can modify the effective splitting, while Kane-Mele SOC mainly opens a gap and is typically weaker than the valley-Zeeman term in graphene- devices [43], [54], [55]. Depending on the relative SOC strengths, the beating pattern can acquire corrections or be suppressed [56]. Beating from spin-split bands has been observed in graphene- heterostructures [43].

6 Universality of scaling relations↩︎

The preceding sections show that a beating-node trajectory contains two distinct pieces of diagnostic information. The first is the magnetic field power. A field-like PMF produces a quadratic trajectory, \(N_{c,j}\propto B_{c,j}^2\), while a Fermi-surface imbalance with \(|\delta N|\propto N\) produces a linear trajectory, \(N_{c,j}\propto B_{c,j}\). A Fermi-surface energy splitting, for which \(|\delta N|\propto N^{1/2}\), also produces a quadratic trajectory. Thus, the field dependence alone does not always identify the microscopic mechanism.

The second diagnostic is the node-index dependence. It separates the two quadratic cases, \[N^{\rm PMF}_{c,j}\propto (2j+1)B_{c,j}^2, \qquad N^{\rm split}_{c,j}\propto (2j+1)^2B_{c,j}^2~.\] Here, “split” denotes intrinsic valley splitting or valley-Zeeman spin splitting in the weak-spin-mixing limit. This distinction is important: both mechanisms give \(B^2\) scaling, but they differ in the node-index dependence.

Beyond these mechanism-dependent node trajectories, the two-channel interference picture also gives a mechanism-independent consistency check: \[\frac{N_{c,j}}{\mathcal{N}_{{\rm osc},j}} = 2(2j+1)\frac{e}{h}B_{c,j}, \qquad \frac{\nu_{c,j}}{\mathcal{N}_{{\rm osc},j}} = 2(2j+1).\] This ratio tests whether the observed beating is governed by two nearby oscillatory channels with weak splitting and reliable node indexing.

For data analysis, the diagnostic can be applied in three steps. First, extract the beating-node positions in the \(N\)-\(B\) plane. Second, fit both the magnetic field power and the node-index dependence of \(N_{c,j}(B_{c,j})\). Third, use \(\nu_{c,j}/\mathcal{N}_{{\rm osc},j}=2(2j+1)\) as a two-channel consistency check and not as a mechanism identifier. The mechanism is constrained by the node trajectory. Table 1 summarizes this diagnostic hierarchy. The table is organized according to the underlying asymmetry class and its microscopic realization. The PMF belongs to the field-like class, whereas constant fractional valley polarization, Y-type Kekulé distortion, intrinsic valley splitting, and valley-Zeeman SOC belong to Fermi-surface classes distinguished by the density power in \(|\delta N|=A N^r\).

Within the two-channel oscillation regime, this formulation is not limited to the monolayer Dirac dispersion. The scaling relations are expressed in terms of carrier density and Fermi-surface area to leading order, so they remain valid whenever Onsager quantization applies, and the oscillatory response is dominated by two nearby frequencies. This is why the PMF scaling extends to bilayer graphene, as verified explicitly in Appendix 11.

Several effects can complicate the direct use of these diagnostics. More than one splitting mechanism may contribute simultaneously, a strong Rashba SOC can invalidate the simple valley-Zeeman mapping, interactions can renormalize the oscillatory spectrum, and disorder or density inhomogeneity can mimic envelope suppression. In such cases, the combined diagnostic still constrains the origin of beating, but it should not be read as a unique mechanism assignment. Extending this analysis to nonlinear quantum oscillations [9], [57], [58] is left for future work.

7 Conclusion↩︎

We have developed a unified framework for identifying the microscopic origin of quantum oscillation beating in graphene-based systems. Starting from Onsager’s quantization relation, we showed that the trajectories of beating nodes in the carrier density vs. magnetic field plane encode direct information about the underlying source of the phase mismatch between oscillatory channels. A central result of this work is that the magnetic field and node-index dependence of the beating nodes provides complementary diagnostics of the underlying mechanism. In particular, strain-induced pseudomagnetic fields, constant fractional valley imbalance, and Dirac-band energy splittings generate distinct scaling trajectories, allowing mechanisms that produce visually similar beating envelopes to be distinguished quantitatively. The resulting classification is summarized in Table 1 and is independent of the specific transport coefficient used to detect the oscillations.

Since the analysis is formulated in terms of carrier density and Onsager quantization, it applies when the oscillatory response is governed by two nearby frequencies with well-defined node indexing. The framework, therefore, provides a practical route for constraining valley- and spin-dependent electronic structure from quantum oscillation measurements. More generally, beating-node spectroscopy can serve as a simple probe of symmetry-breaking fields and band splittings in graphene and related two-dimensional quantum materials.

8 Acknowledgments↩︎

A. Adhikary is supported by the Institute Fellowship, IIT Kanpur. S. Das acknowledges the Ministry of Education, Government of India, for funding support through the Prime Minister’s Research Fellowship. A. Agarwal acknowledges funding from the Core Research Grant by ANRF (Sanction No. CRG/2023/007003), Department of Science and Technology, India.

9 Derivation of oscillatory DOS of SLG↩︎

In this Appendix, we calculate the oscillatory DOS of monolayer graphene using the Green’s function formalism, with disorder-induced broadening encoded through the self-energy. The retarded self-energy in the presence of impurity scattering satisfies \[\begin{align} \label{self95energy} \Sigma^-(\varepsilon) = \Gamma_0^2 \sum_n \frac{1}{\varepsilon - \varepsilon_n - \Sigma^-(\varepsilon)}, \end{align}\tag{23}\] where \(\Gamma_0^2\) is the impurity scattering strength, and the sum is over all energy levels. The resulting DOS from the imaginary part of the self-energy can be evaluated as [59]

\[\label{DOS95self95energy} D(\varepsilon)=\frac{2}{2\pi l_B^2}\,\mathrm{Im}\left[ \frac{\Sigma^{-}(\varepsilon)}{\pi \Gamma_0^2}\right].\tag{24}\] Here, the prefactor \(\dfrac{2}{2\pi l_B^2} = 2\dfrac{eB}{h}\) accounts for the Landau level degeneracy (including spin) per valley. For the zeroth Landau level, taking \(\varepsilon_n = -\Delta\) in Eq. 23 , the self-energy equation reduces to a quadratic form, yielding \(\Sigma^{-}(\varepsilon)=\dfrac{1}{2}\left[(\varepsilon+\Delta)\pm\sqrt{(\varepsilon+\Delta)^2 - 4\Gamma_0^2}\right]\). Substitution into Eq. 24 gives \[D_{\rm LL0}=\dfrac{1}{\pi^2 l_B^2 \Gamma_0^2} \sqrt{\Gamma_0^2-\frac{(\varepsilon+\Delta)^2}{4}}~\Theta(2\Gamma_0-|\varepsilon+\Delta|)~~\text{with}~\Theta(x)=\begin{cases} 1~~\text{for}~~x>0\\ 0~~\text{otherwise} \end{cases}.\] The zeroth LL does not contribute to the oscillatory DOS. Including contributions from non-zero Landau levels, we substitute \(\varepsilon_n = \sqrt{2n(\hbar\omega_c)^2 + \Delta^2}\) into Eq. 23 , which yields a self-consistent equation for self-energy as \(\Sigma^-(\varepsilon)=\Gamma_0^2\sum_{n} f(n)\) where \(f(n) =\big(\varepsilon - \varepsilon_n - \Sigma^-(\varepsilon)\big)^{-1}\). To evaluate the sum over \(n\), we use the residue theorem to write the complex integral as \[\begin{align} \sum_{n=0}^\infty f(n) = -\frac{1}{2\pi i} \oint_C f(z) \pi \cot(\pi z) \, dz= - {\rm Res}[f(z)\pi \cot(\pi z)]=- \left. \frac{\pi \cot(\pi z)}{\frac{d}{dz} (\varepsilon -\varepsilon_z - \Sigma^-(\varepsilon) )} \right|_{z_0}, \end{align}\] where the poles of \(f(z)\) are given by \(z_0 = \frac{(\varepsilon -\Sigma^-(\varepsilon))^2 -\Delta^2}{2 (\hbar \omega_c)^2}\) yields \[\begin{align} \sum_{n=0}^\infty f(n)&=& \frac{\pi (\varepsilon -\Sigma^-(\varepsilon))}{(\hbar \omega_c)^2} \cot\left(\pi \frac{[(\varepsilon -\Sigma^-(\varepsilon))^2 -\Delta^2]}{2 (\hbar \omega_c)^2} \right). \end{align}\] Since \(\Sigma^-(\varepsilon)\) is a complex quantity, we can express it as \(\Sigma^-(\varepsilon) = \Sigma_r + i \Sigma_i\). Hence, the self-consistent expression for the self-energy is given by (assuming \(\varepsilon > \Sigma^-(\varepsilon)\)) \[\Sigma^-(\varepsilon) \approx \Gamma_0^2 \frac{ \pi \varepsilon}{(\hbar \omega_c)^2} \left( \cot(\pi n_0) \right),\] with \(n_0=u-iv\). The quantities \(u\) and \(v\) are given by \(u = \dfrac{(\varepsilon - \Sigma_r)^2 - \Sigma_i^2 - \Delta^2}{2 (\hbar \omega_c)^2} \approx \dfrac{\varepsilon ^2 - \Delta^2}{2 (\hbar \omega_c)^2}\) and \(v = \dfrac{ (\varepsilon - \Sigma_r) \Sigma_i}{ (\hbar \omega_c)^2} \approx \dfrac{ \varepsilon \Sigma_i}{ (\hbar \omega_c)^2}\). The identity \(\cot(\pi(u-i v)) = \dfrac{\sin(2\pi u) + i \sinh(2\pi v)}{ \cosh(2\pi v)-\cos(2\pi u) }\) [60] gives the imaginary part of the self-energy as \[\begin{align} {\rm Im}[\Sigma^-(\varepsilon)] & \approx & \Gamma_0^2 \frac{ \pi \varepsilon}{(\hbar \omega_c)^2} \frac{\sinh(2\pi v)}{ \cosh(2\pi v)-\cos(2\pi u)} =\Gamma_0^2 \frac{\pi \varepsilon}{(\hbar \omega_c)^2} \left[1+ 2 \sum_{p=1}^{\infty} e^{-p 2\pi v} \cos(2\pi p u) \right]. \label{img95self} \end{align}\tag{25}\] Here, we have used the relation \(\dfrac{\sinh 2\pi v}{ \cosh 2\pi v-\cos 2\pi u } = 1 + 2 \sum_{p=1}^\infty e^{-p 2\pi v} \cos (p 2\pi u)\) [60]. We retain only the first harmonic, i.e., \(p=1\). In the limit \(\varepsilon \Sigma_i \ll (\hbar \omega_c)^2\), \(\Sigma_i\) can be obtained iteratively. The first iteration gives \(\Sigma_i = \Gamma_0^2 \pi \varepsilon/(\hbar \omega_c)^2\). Substituting Eq. 25 into Eq. 24 and adding the zeroth Landau level contribution, we obtain the total DOS as follows \[\label{DOS95osc3} D(\varepsilon)=D_{\rm LL0}+D_0(\varepsilon) \left[1+2 \Omega_D(\varepsilon) \cos \left( \frac{2 \pi F(\varepsilon)}{B} \right)\right].\tag{26}\] where \(D_0(\varepsilon)=|\varepsilon|/(\pi \hbar^2 v_F^2)\) is the zero-field DOS, \(\Omega_D(\varepsilon)\) denotes disorder-induced damping, and the oscillation frequency is given by \(F(\varepsilon) = (\varepsilon^2 - \Delta^2)/(2e \hbar v_F^2)\). Equation 26 gives the oscillatory density of states in a high magnetic field for a particular valley, including impurity broadening.

10 Derivation of finite-temperature oscillatory magnetoconductivity↩︎

This Appendix derives the finite-temperature magnetoconductivity in three steps: start from the zero-temperature collisional conductivity, rewrite the Landau-level sum as an energy integral using the oscillatory DOS, and then evaluate the thermal average that generates the Lifshitz-Kosevich damping factor.

Starting from the low-energy Hamiltonian in Eq. 11 for monolayer gapped graphene, the high-field zero-temperature collisional magnetoconductivity is [61] \[\begin{align} \label{sigma95n1} \sigma_{xx}^{\xi} = \frac{e^2}{h}\frac{n_{\rm im} U_0^2}{2\pi k_s^2 l_B^2 \Gamma_0} \sum_{n} \left[ \frac{\varepsilon_n^2-\Delta^2}{\hbar^2 \omega_c^2}\left(1+3\frac{ \Delta^2 }{\varepsilon_n^2}\right) - \xi \frac{4\Delta}{\varepsilon_n} \right] \delta(\varepsilon_n- \mu). \end{align}\tag{27}\] Equation 27 is the oscillatory magnetoconductivity from the \(n\neq 0\) Landau levels. The terms proportional to \(n\) produce the even-in-\(B\) oscillations that survive even without valley polarization. By contrast, the term linear in \(\Delta\) is tied to the valley index \(\xi\) and yields an odd-in-\(B\) contribution. It vanishes after summing over both valleys, but can survive in valley-polarized, time-reversal-symmetry-broken graphene, consistent with the observed odd-parity magnetoresistance in proximitized samples [45]. Since this odd-in-\(B\) term is much smaller than the even-in-\(B\) contribution [45], [61], we neglect it below and focus on the part that controls the SdH oscillations.

For the \(K\) valley, we write the finite-temperature oscillatory magnetoconductivity as

\[\sigma_{xx} = \frac{\sigma_0}{\pi l_B^2 } \sum_{n} A(\varepsilon_n)\, (-\partial_{\varepsilon} f)~~\text{with}~ ~\sigma_0 = \frac{e^2}{h} \frac{n_{\rm im} U_0^2}{2 k_s^2 \Gamma_0} ~~\text{and} ~~A(\varepsilon_n) = \frac{\varepsilon_n^2-\Delta^2}{\hbar^2 \omega_c^2}\left(1+3\frac{ \Delta^2 }{\varepsilon_n^2}\right).\] The next step is to replace the discrete Landau-level sum by an energy integral weighted by the DOS. Using \(\frac{2}{2\pi l_B^2} \int dn\, g(\varepsilon_n) = \int g(\varepsilon)\, D(\varepsilon)\, d\varepsilon\) [60], we obtain \[\sigma_{xx} = \sigma_0 \int A(\varepsilon)\, D_0(\varepsilon)\, \left[1+2 \Omega_D(\varepsilon) \cos\!\left( \frac{2 \pi F(\varepsilon)}{B} \right)\right] (-\partial_\varepsilon f)\, d\varepsilon.\] To isolate the temperature dependence, we assume that \(A(\varepsilon)\), \(D_0(\varepsilon)\), and \(\Omega_D(\varepsilon)\) vary slowly on the scale \(k_B T\) compared to the oscillatory phase. We therefore evaluate them at \(\varepsilon=\mu\) and pull them out of the integral,

\[\sigma_{xx} \approx \sigma_0\, A(\mu)\, D_0(\mu) \left[1+2 \Omega_D(\mu) I(T)\right],\] where all temperature dependence is contained in \[I(T)=\int_{-\infty}^{\infty} d\varepsilon\, \left(-\partial_\varepsilon f\right) \cos\!\left( \phi \right), \qquad \phi=\frac{2 \pi F(\varepsilon)}{B}.\] Expanding the phase around the chemical potential, \(\phi(\varepsilon)=\phi(\mu)+(\varepsilon-\mu)\, \phi'(\mu)\), and dropping the sine term, we obtain \[\begin{align} I(T)&=&\cos \left[\phi(\mu)\right] \int_{-\infty}^{\infty} \frac{\beta_T \cos[(\varepsilon-\mu)\phi'(\mu)]}{4\cosh^2 \left(\beta_T(\varepsilon-\mu)/2\right)}\, d\varepsilon, \end{align}\] where \(\beta_T = 1/(k_B T)\) is the inverse temperature. Using the substitution \(y=\frac{\beta_T(\varepsilon-\mu)}{2}\), with \(d\varepsilon=\frac{2}{\beta_T}dy\), we have \[\begin{align} I(T)&=& 2\cos \left[\phi(\mu)\right] \int_{-\infty}^{\infty} \frac{\cos \left[\frac{2\phi'(\mu)}{\beta_T}y\right]}{4\cosh^2 y}\, dy = \frac{T/T_D}{\sinh(T/T_D)}\, \cos \left[\phi(\mu)\right]. \end{align}\] Here, we used the standard integral \(\int_{-\infty}^{\infty} \frac{\cos(ay)}{\cosh^2 y}\, dy = \frac{\pi a}{\sinh\left(\pi a/2\right)}\). This step produces the familiar thermal damping factor, with Dingle temperature \(T_D = (\hbar\omega_c)^2/(2\pi^2 k_B \mu)\). The finite-temperature conductivity for valley \(\xi=+1\) then becomes \[\sigma_{xx} \approx \sigma_0\, A(\mu)\, D_0(\mu) \left[1+2 \frac{T/T_D}{\sinh(T/T_D)} \Omega_D(\mu) \cos \left[\phi(\mu)\right]\, \right].\] The magnetoconductivity for valley \(\xi\) is then \[\begin{align} \sigma_{xx}^\xi \approx \frac{e^2}{h}\frac{n_{\rm im} U_0^2}{2k_s^2 \Gamma_0} \frac{\mu}{\pi (\hbar v_F)^2} \left[ \frac{\mu^2 -\Delta^2}{2(\hbar \omega_c)^2}\left(1+3\frac{\Delta^2}{\mu^2}\right) \right] \left[1+ 2 \Omega_T \Omega_D(\mu) \cos\left( \frac{2\pi F(\mu)}{ B}\right) \right] . \label{sigma95osc95mu} \end{align}\tag{28}\] Here, \(\Omega_T = \dfrac{T/T_D}{\sinh(T/T_D)}\) is the temperature-dependent broadening factor. This is the Lifshitz-Kosevich form quoted in the main text. Substituting \(N_\xi=(\mu^2-\Delta^2)/(2 \pi \hbar^2 v_F^2)\), \(F(N_\xi)=h N_\xi/(2 e)\), and allowing the field to be valley dependent, \(B\to B_\xi\), gives the single-valley expression

\[\begin{align} \sigma_{xx}^\xi \approx \frac{e^2}{h}\frac{n_{\rm{im}} U_0^2}{2 k_s^2 \Gamma_0} \frac{2N_\xi}{ e \hbar v_F^2 B_\xi} \frac{\pi \hbar^2 v_F^2 N_\xi + 2\Delta^2}{\sqrt{2 \pi \hbar^2 v_F^2 N_\xi + \Delta^2}}\bigg[1+ 2 \Omega_T \Omega_D \cos\left( \frac{2\pi F(N_\xi)}{ B_\xi}\right) \bigg] . \label{sigma95osc} \end{align}\tag{29}\] This expression is the characteristic Lifshitz-Kosevich form of quantum oscillations in monolayer graphene and is the starting point for the beating analysis in the main text.

11 Derivation of scaling relation in BLG↩︎

This Appendix extends the strain-induced PMF scaling analysis to BLG. We first review the Landau levels in the low-energy Hamiltonian and then derive the oscillatory DOS, magnetoconductivity, and resulting scaling behavior.

Starting from the full four-band tight-binding description of Bernal-stacked bilayer graphene [62], projection onto the two low-energy bands near charge neutrality gives [18], [63] \[{\cal H} = -\frac{1}{2m}\left[\left(p_x^2 - p_y^2\right)\sigma_x+ \xi\,(p_y p_x + p_x p_y)\,\sigma_y\right]+ \Delta\,\sigma_z ,\] where \(m=\gamma_1/(2 v_F^2)\) is the effective mass, \(\gamma_1\) denotes the interlayer dimer hopping parameter, \(v_F\) is the Fermi velocity of monolayer graphene, and \(\Delta\) denotes a sublattice potential. In the presence of a magnetic field (\(\boldsymbol{B}=B\hat{z}\)), the momentum \(\boldsymbol{p}\) is replaced by \((\boldsymbol{p}+e\boldsymbol{A})\). Under the Landau gauge \(\boldsymbol{A} = (0,xB,0)\), the energy eigenvalues are given by \[\varepsilon_n = \begin{cases} \lambda \sqrt{\Delta^2 + \hbar^2 \omega_c^2 n(n-1)}, & n \neq 0,1, \\ -\xi\Delta, & n = 0,1. \end{cases}\] where \(\omega_c=eB/m\) is the cyclotron frequency. Substituting \(\cos{\alpha}=\Delta/\varepsilon_n\), the wavefunction for \(n \neq 0,1\) can be written as \[\Phi_{n,k_y}^{\xi=1}(X) = \frac{e^{ik_y y}}{\sqrt{L_y}} \begin{pmatrix} \cos{(\alpha/2)} \phi_{n-2}(X) \\ \sin{(\alpha/2)} \phi_n(X) \end{pmatrix}, \Phi_{n,k_y}^{\xi=-1}(X) = \frac{e^{ik_y y}}{\sqrt{L_y}} \begin{pmatrix} \cos{(\alpha/2)} \phi_{n}(X) \\ \sin{(\alpha/2)} \phi_{n-2}(X) \end{pmatrix}. \label{wavefun}\tag{30}\] For the beating pattern, we work in the high-Landau-level regime, where many distinguishable LLs are occupied. The energy spectrum then simplifies to \(\varepsilon_n \approx \sqrt{\Delta^2 + n^2 \hbar^2 \omega_c^2 }\), as reported in Ref. [[29]].

Following the monolayer analysis, we evaluate the oscillatory DOS in BLG using the self-energy formalism. The corresponding self-energy follows from Eq. 23 , where poles of \(f(z)\) are given by \(z_0 = \sqrt{(\varepsilon -\Sigma^-(\varepsilon))^2 -\Delta^2}/(\hbar \omega_c)\). For simplicity, we take \(\Delta=0\), giving \[\label{self95energy952} \Sigma^-(\varepsilon) \approx \Gamma_0^2 \frac{ \pi }{\hbar \omega_c} \left( \cot(\pi n_0) \right),\tag{31}\] where \(n_0=u-iv\) with \(u = (\varepsilon - \Sigma_r)/ \hbar \omega_c\) and \(v = \Sigma_i/\hbar \omega_c\). Substituting Eq. 31 into Eq. 24 , we obtain the density of states \[\begin{align} D(\varepsilon) &=& D_0 \left[1+ 2 \Omega_D\cos\left( \frac{2\pi F(\varepsilon)}{B}\right) \right] . \label{DOS95osc2} \end{align}\tag{32}\] where \(D_0 = m/(\pi \hbar^2)\) denotes the zero-field DOS, \(\Omega_D=\exp\left(-\dfrac{2\pi^2\Gamma_0^2}{\hbar^2\omega_c^2}\right)\) is the damping factor arising from Landau-level broadening, and \(F(\varepsilon)=m\varepsilon/(e\hbar)\) represents the oscillation frequency. Equation 32 gives the oscillatory density of states in a high magnetic field for a particular valley.

To obtain the finite-temperature magnetoconductivity and beating behavior, we start from the zero-temperature oscillatory magnetoconductivity for BLG. The longitudinal conductivity arises mainly from charged-impurity scattering of cyclotron orbits, giving the collisional conductivity \[\label{final} \sigma_{xx}=C(B) \sum_{n} \left[(4n-2)\left(1+\frac{\Delta^2}{\varepsilon_n^2}\right) - \xi\frac{8\Delta}{\varepsilon_n} \right]\delta(\varepsilon_n - \mu).\tag{33}\] As in monolayer graphene, we can extend the analysis to the finite-temperature limit. For simplicity, we take \(\Delta=0\). In the high-Landau-level regime, substituting the carrier density per valley, \(N_\xi=m\mu/(\pi \hbar^2)\), the finite-temperature magnetoconductivity for valley \(\xi\) becomes \[\begin{align} \label{sigma95blg} \sigma_{xx}^\xi \approx && \sigma_0\left[1+ 2 \Omega_T \Omega_D \cos\left( \frac{2\pi F(N_\xi)}{B_\xi}\right) \right]~~\text{with} ~~\sigma_0=\frac{e^2}{h}\frac{n_{\rm im} U_0^2}{2k_s^2 \Gamma_0} \frac{4 m N_\xi}{\hbar e B_\xi}. \end{align}\tag{34}\] where \(\sigma_0\) is the non-oscillatory contribution and \(F(N_\xi)=hN_\xi/2e\) is the oscillation frequency. Here, \(\Omega_T = \dfrac{T/T_D}{\sinh(T/T_D)}\) is the temperature-dependent broadening factor, with Dingle temperature \(T_D = \hbar\omega_c/(2\pi^2 k_B)\). This is the Lifshitz-Kosevich form for oscillatory magnetoconductivity in graphene.

In the presence of a valley-contrasting pseudomagnetic field, we substitute \(B_\xi =B-\xi B_{\rm{pm}}\) and \(N_\xi=N/2\) (\(N\) is the total carrier density) into Eq. 34 to obtain the conductivity in each valley. Neglecting the weak \(B_{\rm pm}\) dependence of the prefactors and damping factors, summing over the two valleys, \(\sigma_{xx}(N, B) = \sum_\xi \sigma^{\xi}_{xx} (N/2, B_\xi)\), gives the total magnetoconductivity \[\begin{align} \sigma_{xx}(N,B) \approx 2\sigma_0 \bigg[1+ 2 \Omega_T \Omega_D \cos\left( \frac{2\pi F(N/2)}{ B}\right)\cos\left( \frac{2\pi F(N B_{\rm{pm}} /2B ) }{B} \right) \bigg]~. \end{align}\] The \(j\)-th beating node occurs when \(\cos\left( \dfrac{2\pi F(N B_{\rm{pm}} /2B ) }{B} \right) =0\), giving \(h N |B_{\rm pm}|/(2 e B_{c,j}^2) = (2j + 1)/2\). In terms of total carrier density, we have \[N^{\rm{PMF}}_{c,j} = \frac{e B_{c,j}^2}{ h |B_{\rm pm}|} (2 j+1)~,\] This is the same scaling obtained from Onsager quantization and from the microscopic single-layer calculation.

References↩︎