Low-frequency Quantum Oscillations from Interactions in Layered Metals

Metals composed of weakly-coupled, stacked layers possess a Fermi surface that slightly varies in size along the stacking direction. This appears in de Haas-van Alphen (dHvA) oscillations of the magnetisation with magnetic field as two close frequencies, corresponding to the two extremal Fermi surface cross-sectional areas. We show that, for layered materials of sufficiently high mobility, Coulomb interactions can have a dramatic effect on the form of the dHvA oscillations: there is also generically an oscillation at the small difference of the two large frequencies. We determine the size and form of this effect, and show that it probes the short-range part of the Coulomb interactions within the layered material. We argue that this interaction effect may explain recent experimental observations of anomalous low-frequency dHvA oscillations in the ultrapure delafossites.

Introduction -Modern materials science has produced a wide variety of new sorts of solid state systems, such as those characterized by non-trivial topology or strong interactions. The novel properties of these materials have led to reexaminations of previously well-understood phenomena. Quantum oscillations (QOs) of the magnetization as a function of (inverse) magnetic field [1], long interpreted in terms of the geometry of closed Fermi surfaces, and the standard Lifshitz-Kosevich (LK) theory [2], have been necessarily reexplored following the discovery of Weyl and nodal semimetals [3][4][5][6][7], quasiperiodic systems [8,9], and the observation of QOs in insulators [10][11][12][13][14][15][16][17], systems all featuring unusual or absent Fermi surfaces.
Following in this vein, recent de Haas-van Alphen (dHvA) studies on delafossites, a class of layered materials featuring strong interactions and high in-plane mobility, also display anomalous behavior-large, lowfrequency QOs at the difference of the two natural high frequencies related to extremal Fermi surface areas [18,19]. These observations are not readily explained by existing theory based on magnetic interactions [20], and require a re-examination of the theory of the dHvA effect in these new materials settings.
In this Letter, we show that Coulomb interactions can lead to difference-frequency oscillations of the magnetization, of a size that can account for these experimental observations. The mechanism that we identify is overlooked in long-standing theories of interaction effects on the dHvA effect [20][21][22], as it requires the retention of terms that are smaller than those needed to account for the high-frequency oscillations. However, in high-mobility systems such as the delafossites, the dominant dHvA oscillation at the difference frequency in both bilayer and multi-layer systems can arise from the Coulomb interaction mechanism we identify. In these circumstances our theory shows that the measurement of a differencefrequency oscillation is a signature of strong non-local interactions. Our results indicate that, in general, new considerations are warranted when analyzing dHvA oscillations in high-mobility metals.
To trace the origin of a difference frequency in QOs, first consider a non-interacting Fermi gas whose Fermi surface has two extremal orbits of similar areas. As the magnetic field varies, the Landau quantized orbits sweep through these two extremal areas, causing thermodynamic quantities to oscillate at the two frequencies f ± . Within a model of non-interacting electrons at fixed chemical potential these oscillations are independent. However, any nonlinearities that couple these two oscillations at f ± , can lead to a difference frequency δf = |f + −f − | (as well as side-bands). One source of nonlinearity is the magnetic interaction, by which the magnetic field B = µ 0 (H +M ) acquires an oscillatory component through the oscillating magnetization M (H). This leads to a difference frequency component of size related to the geometry-dependent demagnetization field [20]. However, this effect appears insufficient to explain the difference-frequency oscillations of magnetization seen in recent experiments [18,19].
As we will show, a much larger difference frequency oscillation can arise from Coulomb interactions. One consequence of the Coulomb interaction is that its long-range component forces overall charge neutrality in bulk 3D materials, requiring one to work at fixed electron number density n rather than fixed chemical potential. This introduces a non-linearity, mediated by oscillations in the chemical potential [23,24], that can readily produce a large difference frequency oscillation in the conductivity (Shubnikov-de Hass effect) [25]. While we find that this can also produce a difference frequency oscillation in the magnetization (i.e. in the dHvA effect), this effect by itself is very small. In fact, as we will show, the dominant effect arises from a full consideration of the Coulomb interactions that also takes account of the (screened) shortrange interactions that couple local fluctuations of the charge density. The combination of both short-and longrange components of the Coulomb interaction produces a much larger difference frequency oscillation than either can independently.
Bilayer (2D) -To illustrate the essential physics underlying this mechanism, we first consider a bilayer model. (We neglect possible strongly correlated phases that can arise from strong interactions in 2D Landau quantized systems.) We consider a pair of identical, parallel 2D gases of spinless electrons, with parabolic inplane dispersion of effective mass m * . Interlayer hopping, of amplitude ∆/2, splits the energy eigenstates into symmetric and antisymmetric subbands, and a perpendicular magnetic field B reorganizes the in-plane states into Landau levels. The spectrum is then l,s = ω c (l + 1/2) − s∆/2, with ω c = eB/m * the cyclotron frequency, l = 0, 1, 2, . . . indexing Landau levels, and s = ± marking even and odd subbands. For a chemical potential µ, the densities of the two subbands at B = 0 are n ± (0) = (µ ± ∆/2)m * /(2π 2 ) and there are two distinct Fermi surfaces. At non-zero B these give rise to QOs in the magnetization at frequencies f ± = en ± (0)/h = m * (µ ± ∆/2)/e . Within LK theory these are the only two frequencies in the magnetization oscillations.
At fixed chemical potential µ, the occupations of the two subbands, n ± (B), both oscillate with B. For ω c ∆ µ, so that many Landau levels are occupied and n ± have similar average values, these oscillations are on the scale of n Φ = B/Φ 0 , the degeneracy of each Landau level per unit area, with Φ 0 = h/e the magnetic flux quantum. Keeping only the first harmonic, we write where η 1 is introduced phenomenologically to account for any disorder and temperature effects that suppress oscillations. The oscillations in density Eq. (1) lead to oscillations in the Coulomb interactions that can also have components at the difference frequency. Taking the ionic background charge density to be n I in each well, and computing the expectation value of the electron-electron interactions in the non-interacting ground state leads to the interaction energy withn ≡ (n + + n − )/2 and δn ≡ n + − n − . The first term in (2) is the Hartree energy, related to the overall electrostatic energy of fluctuations of the total charge density. The second term is the Fock term, describing the suppression of short-range repulsion due to exchange [26].
Since the subband densities n ± oscillate at the frequencies f ± , the cross term n + n − produces a small component that oscillates at the difference frequency δf = |f + − f − |, with amplitude of order V η 2 n 2 Φ .
To account for the effects of the long-range part of the Coulomb interaction, we enforce charge neutrality, requiringn(B) ≡ n I . The chemical potential µ becomes an oscillating function of B, which to first order in η is With this oscillating µ(B) inserted in Eq. (1) one finds that the density n ± (B) oscillates not just at f ± but also at the difference frequency δf (and at other sidebands). The restriction to fixed n has a large effect on the interaction energy. The Hartree term vanishes exactly, leaving just the Fock exchange energy, −V [δn/2] 2 , which contains an oscillatory term of size This gives rise to oscillations at the difference frequency δf , with amplitude of order V η 2 n Φ δn(0). This is larger than the δf oscillation we identified at fixed chemical potential by a factor of δn(0)/n Φ = ∆/( ω c ), which can be large at low fields.
Multi-layer (3D) -We now turn to a model for a multi-layer metal. As above, each layer is described by a parabolic dispersion with effective mass m * , electrons may hop between adjacent layers with an amplitude that we now denote t ⊥ , and a perpendicular field B breaks inplane states into Landau levels. At nonzero temperature T the system is described by the action S 0 = n ,l,ky kzψ l,ky (k z , n ) (−i n + ξ l (k z )) ψ l,ky (k z , n ), (5) where ψ,ψ are the electron field operators in the energy eigenbasis, n = (2n + 1)πk B T is the Matsubara frequency, and ξ l (k z ) = ω c (l + 1/2) − 2t ⊥ cos(k z a ⊥ ) − µ is the single-particle energy measured from the Fermi level. The discrete eigenstate index s of the bilayer is replaced by the continuous quasi-momentum k z ∈ (−π/a ⊥ , π/a ⊥ ] describing dispersion along the c-axis, and we use the notation kz = π/a ⊥ −π/a ⊥ dk z /2π, where a ⊥ is the interlayer spacing. With our choice of Landau gauge, k y indexes the degenerate states in each Landau level with total number ky = An Φ , with A the sample area. The frequencies determined by the two extremal cross-sectional areas of the Fermi surface along k z are f ± = m * (µ ± 2t ⊥ )/e [27]. As before we consider spinless electrons, as this already shows the new effect.
For simplicity, in considering the interlayer interaction between local charge density fluctuations we approximate the (screened) Coulomb potential as only acting between pairs of nearest points on neighboring layers, , where ε is the permittivity and λ 2 is the area of the "patch" on each layer that participates in the interaction. The interaction term in the action is (6) where x = (τ, r). We also acquire a term S I = βAL z V λ 2 n I (n I /a ⊥ − 2n), where L z is the extent of the system in z, which depends only on the total electron density n and accounts for the interaction between electrons and lattice ions. Note that here n represents a number per unit 3D volume, whereas n I is a number per unit 2D area in-plane.
To analyze the interaction effects we proceed as follows. The interlayer interaction can be included to first order with the Hartree-Fock self-energy Σ, inserted into the full electron Green's function, It is important for our analysis to keep both the constant, zero-field part Σ 0 and oscillatory part Σ(B) of the self-energy; previous studies of the dHvA effect including interactions have used the general relation Σ Σ 0 to discard Σ entirely [20][21][22], but we find that this term produces the leading contribution to difference frequency oscillations in this sort of layered system.
From the Green's function G we calculate the grand potential Ω via standard field theoretic methods, expanding up to first order in the interaction constant V . We then obtain the free energy F as the Legendre transform of the grand potential, F (B, n) = Ω(B, µ(B, n)) + nµ(B, n), where n is the fixed electron density, and µ(B, n) is the oscillatory chemical potential needed to fix the electron density. To linear order in V , it is sufficient to fix the density using just the non-interacting part of the theory, i.e. defining µ(B, n) through n = −∂Ω/∂µ| V =0 ; corrections to µ that depend on V only lead to terms in F that are of order V 2 and higher [28]. The result of this procedure-the free energy F (B, n)-is the relevant thermodynamic potential for a system with fixed electron density, including all first order interlayer interaction effects.
We calculate the Hartree-Fock self-energy in the approximation that the interaction is independent of Landau-level index, This gives the exact form of Ω up to first order in V . Here we define the 3D number density of occupied states at momentum k z , n(k z , B) = n Φ d l n F ( )A( − ξ l (k z )), where A is the spectral density. With the 3D electron density fixed to n(B) = n = n I /a ⊥ we determine the oscillatory chemical potential µ(B, n) [28], then obtain the free energy F 3D = F (0) As in the 2D case, the interaction part F 3D is given solely by the Fock energy.
We analytically evaluate the oscillatory part of the free energy assuming the hierarchy of energy scales ω c 2t ⊥ µ. Equivalently this sets a hierarchy of 2D densities n Φ n ⊥ n I , where we define n ⊥ ≡ t ⊥ m * /π 2 . The dominant contributions to oscillations at f ± and δf are found to be [28] from the kinetic part, and from the interacting part, where are the Dingle factor and LK temperature factor, accounting for finite quasiparticle lifetime τ qp and nonzero temperature T respectively. Note that in both Eqs. (11) and (12) the δf terms are found to depend on the square of these factors, while the f ± terms only depend on a single power of each. We verify our analytic calculations by analyzing the clean, T → 0 limit of the system numerically. The chemical potential to fix n is evaluated numerically, the full F 3D are evaluated on a set of points evenly spaced in 1/B, a polynomial background is subtracted off, and the spectral content of the resulting oscillatory data is analysed via discrete Fourier transform. The resulting Fourier spectra are presented in Fig. 1. The amplitudes at f ± and δf given by this analysis are found to closely match the τ qp → ∞, T → 0 limit of the analytic results, Eqs. (11) and (12).
Discussion -In the above analysis we have determined the oscillations of the free energy of the 3D system. The oscillatory part of the magnetization can be  obtained as M = −∂ F 3D /∂B, and naturally separates into interacting and non-interacting parts in the same way as F 3D , M = M (0) + M (1) . The largest contributions arise from the derivative acting on the oscillating factors themselves, not the preceding amplitudes, as long as k B T, /τ qp t ⊥ . Thus, the oscillation amplitudes of M can be acquired from those of F 3D in Eqs. (11) and (12) by simply multiplying each term with a factor of 2πf /B 2 , where f is the oscillation frequency of that term.
To compare the predictions of this theory with the experimental observations on PtCoO 2 , we use the parameters reported in Ref. 19 [29]. Their investigation into short-range interactions in this material suggests an onsite Hubbard-like repulsion with U 6 eV. We take this to suggest short-range Coulomb interaction energies in general, including V , are on the scale of eV, and use λ 2 on the scale of the in-plane area of the PtCoO 2 unit cell. With these parameters, we find that the difference frequency amplitude of M (0) is less than that of M (1) for all fields B < 4πV λ 2 n ⊥ m * / e ∼ 600 T. The experimentally relevant low-field regime is far below this threshold, so within our theory dHvA oscillations at δf are dominated by the interaction contribution, M (1) . This amplitude is given by | M (1) | δf = V λ 2 δf R 2 D,1 R 2 T,1 /(π 2 Φ 2 0 a ⊥ ). Note that this depends on magnetic field strength only through the squares of the Dingle and LK temperature factors, while the corresponding amplitude of the noninteracting term M (0) depends on an additional factor of B itself. This form of | M (1) | δf qualitatively matches the experimental results of the field-dependence of the lowfrequency oscillations in Ref. 19. Furthermore, adding a factor of 2 due to spin, the magnitude of this term is consistent with the measurements [30]. Future measurements, focusing on the quantitative size of this lowfrequency component, have the potential to allow a better understanding of the strength and nature of interlayer interactions.
The careful reader may note that the inequality δf f ± , inherited from t ⊥ µ, suggests that the dHvA oscillations our theory predicts at δf should always be much smaller than those at f ± . This is in contrast with what is seen in experiment. However, this apparent conflict may be reconciled by considering small-angle scattering off long-range disorder, e.g. spacial inhomogeneity of µ. It is well understood that long-range disorder suppresses QO amplitudes, and that the mean free path determined from Dingle measurements can be smaller than what is found in transport, which is less sensitive to small-angle scattering. Indeed this is seen in experiment [19], where the mean free path determined from Dingle measurements is about ten times smaller than that found in transport. This significant small-angle scattering strongly suppresses main frequency oscillations, however the difference frequency effect we find, being independent of µ, is insensitive to this sort of disorder effect, allowing the two to be of a comparable scale. The starkly different dependence on small-angle scattering suggests that some measure of the strength of long-range disorder may be obtained by comparing Dingle measurements of difference and main frequency oscillations.
In summary, by considering oscillations of the selfenergy that are neglected in standard theories [20][21][22], we have shown that Coulomb interactions can give rise to a new form of low-frequency dHvA oscillations in high mobility layered materials. These difference frequency oscillations are highly sensitive to the short-range disorder, their amplitude being suppressed by the square of the Dingle factor, but are insensitive to long-range disorder. They can therefore become a dominant feature of oscillations in very high-mobility materials such as the delafossites. Our theory shows that the size of the difference frequency oscillation is a measure of the strength and form of the interactions beyond the on-site Hubbard U . Although we have focused on simple layered metals, our results suggest that dHvA studies may be of use to probe properties beyond just Fermi surface geometry in more general bilayer and multi-layer systems.
We thank Elena Hassinger and Andrew Mackenzie for helpful discussions, for sharing experimental data, and for comments on an earlier draft. This work is supported by EPSRC Grant No. EP/P034616/1 and by a Simons Investigator Award.

Grand Potential
The method of calculating the grand potential of a non-interacting multi-layer system is presented in Ref 23. Ignoring the effect of spin and assuming ω c 2t ⊥ we have where ν 0 = m * /(2π 2 ) is the 2D density of states per spin at zero field. The contribution to the grand potential at first order in the interlayer interaction can be simply calculated as well. We use the part of the Hartree-Fock self-energy that is diagonal in Landau level index, Because we are calculating the contributions to Ω only up to first order in V this approximation gives the correct result; directly computing the diagrams for the Hartree and Fock contributions to Ω gives the same result as we find below, though doing so obscures the point of departure from Ref. 21. The first order interaction term is then 16) where V = AL z is the volume of the system. In writing this result we have incorporated the contribution from the interaction with the ionic lattice, which contributes only at first order in V . It is clear from this form that if the electron density n is set equal to n I /a ⊥ then the first term in the final expression exactly vanishes, with the Hartree term cancelling against contribution from the ionic lattice.

Fixing Electron Density
When holding the particle density n constant, the chemical potential necessarily becomes a function of the applied field, µ = µ(B). We can get some insight into the form of µ(B) by first examining how the density changes for fixed µ in a non-interacting single layer. The degeneracy of each Landau level is N Φ = AB/Φ 0 , so the number of electrons in each level increases with B. Therefore the total particle number for low temperature, k B T ω c , and fixed µ suddenly decreases by N Φ every time a Landau level crosses above the Fermi energy, then rises again as the degeneracy of the still-filled states continues to increase. The overall result is for n(B) to oscillate around n = n(B = 0) in a sawtooth pattern.
To keep the density n constant instead, the chemical potential must "stick" to a single Landau level as it moves in energy-the filling fraction of the highest level decreases to compensate for increasing degeneracy. When the top level is completely emptied, µ(B) then jumps to the level next highest in energy. We see that the chemical potential must then oscillate in a sawtooth pattern around its value at B = 0, with the value of µ(B) equal to the energy of some particular Landau level.

The Impact of Interactions on the Chemical Potential
The above intuition ignores interactions, which in general will non-uniformly redistribute electrons between the levels of the system and therefore affect how the chemical potential must oscillate in order to keep the total density fixed. Consider now an interacting system. The grand potential may be written as where Ω V contains all interactions, so that Ω V =0 = Ω (0) . The electron density in the system is Define µ(B, n) to be the oscillatory chemical potential needed to exactly fix the electron density to the value n, i.e. n(B, µ(B, n)) ≡ n. We can separate µ(B, n) into oscillatory and non-oscillatory parts as where µ 0 is the constant value of the chemical potential for B = 0, and µ is the full oscillatory component, with µ V =0 ≡ µ (0) , so that µ V contains all oscillations induced by the interaction. Similarly we note n (0) (B, µ 0 + µ) V =0 = n (0) (B, µ 0 + µ (0) ) and define n (0) ) to be the entire part of n (0) (B, µ 0 + µ) that depends on interactions. The fixed-density condition then becomes with all effects of the interaction concentrated into the last two terms, marked with subscript V . The value of n is entirely independent of the interaction, equal to the ionic charge density in order for the system to be charge neutral overall. For the fixed-density condition to be satisfied we must then have We now examine the free energy of this system, related to the grand potential by a Legendre transform. Expanding Ω (0) around the non-interacting parts of the chemical potential, we have F (B,n) = Ω(B, µ(B, n)) + n µ(B, n) = Ω (0) (B, µ(B, n)) + n µ(B, n) + Ω V (B, µ(B, n)) where in the second line we use Eq. (A.21) to relate the derivative of Ω (0) evaluated at the non-interacting chemical potential to the fixed density n. We see that the interaction appears only within Ω V and in terms second order and higher in µ V . Because Ω V and µ V are themselves at least first order in the interaction, the only term in this expression that may occur at first order overall is obtained from the non-interacting part of the chemical potential inside Ω V . Therefore, if we discard contributions above first order, the free energy is where Ω (1) is the part of Ω V first order in interactions. Crucially, we see that µ V does not appear. We conclude, then, that if we only consider the interaction up to first order it is unnecessary to consider how interactions affect the fixed-density condition, and the only relevant oscillatory part of the chemical potential is determined through Eq. (A.21).

Oscillatory Chemical Potential
We follow Champel and Mineev [23] to derive an oscillatory chemical potential using the grand potential Ω 3D . As argued above, it suffices for our purposes to consider only the density coming from the non-interacting part of the grand potential. This is which we demand to be equal to the density n = n I /a ⊥ by letting µ = µ(B, n). Note that we arrive at exactly the same expression if we instead calculate n(B) from the spectral density, n(B) = n Φ kz d n F ( ) ∞ l=0 A( − ξ l (k z )). This then gives where µ 0 (n) = a ⊥ n/ν 0 = n I /ν 0 in our system, and R D,p and R T,p are the Dingle factor and Lifshitz-Kosevich temperature factor defined in the main text. In principle this equation determines µ(B, n) self-consistently. However, we note that the oscillatory part is much smaller than the leading constant term because ω c 2t ⊥ µ 0 and R D , R T ≤ 1, so better and better approximations to the exact form of µ(B, n) can be determined by iteration, starting with the replacement µ(B, n) → µ 0 in the right-hand side of the equation. For our purposes this starting point is a sufficient approximation, giving us Indeed, if we were to consider the next iteration by inserting this form of µ(B, n) in place of µ 0 within the sine, the additional terms we would acquire would be even smaller still and would only contribute very small corrections to our main results.

Free Energy Oscillations
Free energy is related to the grand potential through a Legendre transformation, Because µ(B, n) is independent of interactions, the free energy separates cleanly into non-interacting and interacting parts related to the corresponding parts of the grand potential, which we now analyze in turn focusing specifically on the oscillatory components at frequencies f ± and δf .
We start with the non-interacting part. We drop the terms of Ω (0) that are independent of µ since they will not yield any oscillations after substituting µ → µ(B, n) = µ 0 + µ(B). We have (A.29) The dominant contributions to oscillations at the frequencies of interest arise from the p = 1 term of this sum since the size of the summand decreases quickly with growing p. The same is true of the sum in the definition of µ itself, allowing the same approximation. Because ω c 2t ⊥ and R D , R T ≤ 1 we can expand the factors involving sine and cosine of µ, keeping just the single largest term in each case, Trigonometric identities let us combine the product of sines in the last term, giving one term oscillating at the sum of f α and f β and one term oscillating at their difference, which is either constant or equal to δf . Discarding all constant terms and all terms with an oscillation frequency other than f ± or δf , we are left with just F 3D as in Eq. (11), Turning now to the interacting part of the free energy we have, Employing the same techniques as in the calculation of Ω we obtain The constant term is much larger than the oscillatory term here, so for the free energy, in which this term appears squared, the dominant oscillatory contribution will result from the cross term. Keeping just this part, we have As above, the largest oscillatory contributions come from the p = 1 term of this sum, and also for the sum within µ.
With these approximations and expanding the sine and cosine of µ we have The first term gives oscillations at f ± . The trigonometric functions of the second term can be combined to give oscillations at the sum or difference of f α and f β . Keeping just the terms oscillating at f ± or δf we then obtain F As a point of comparison, we provide the oscillation amplitudes at the frequencies f ± and δf in both the case of fixed chemical potential Ω 3D and fixed electron density F 3D , further subdivided into interacting and non-interacting contributions. The expressions are written in terms of the densities n Φ = B/Φ 0 = ν 0 ω c and n ⊥ = 2ν 0 t ⊥ , as well as the Dingle and LK temperature factors. We see that the the restriction to fixed n does not affect the amplitude of oscillations at f ± , but it is responsible for the generation of large difference frequency oscillations. For fixed chemical potential there are no oscillations at δf at all in the non-interacting case, and in the interacting case they are very small, n Φ /n ⊥ = ω c /2t ⊥ 1 smaller than in the case of fixed density.
Amplitude at δf Amplitude at f±