De Haas-van Alphen effect and quantum oscillations as a function of temperature in correlated insulators

TL;DR

This paper investigates the behavior of excitonic insulators, revealing how hybridization, magnetic field, and temperature influence quantum oscillations and magnetization, with analytical insights into their characteristics.

Contribution

It introduces a self-consistent study of hybridization effects in excitonic insulators, showing how magnetic field and temperature modulate quantum oscillations and magnetization.

Findings

Hybridization gap decreases with magnetic field, indicating diamagnetism.

De Haas-van Alphen oscillation amplitude varies periodically with inverse temperature.

Analytical expressions for the position and height of dominant oscillation peaks.

Abstract

We theoretically study a model of excitonic insulators which show de Haas-van Alphen oscillations as well as periodic dependence of the magnetization on inverse temperature. The insulating behavior is due to the Coulomb interaction driven hybridization of fermions at the crossing point of their energy bands. We study this hybridization self-consistently and recover known results by A. Allocca and N. Cooper, SciPost Phys. {\bf 12}, 123 (2022). In addition to known results, we show that the hybridization gap decreases with the magnetic field which corresponds to the diamagnetism. Furthermore, we show that the amplitude of the de Haas-van Alphen is oscillating with inverse temperature with a period defined by a combination of the hybridization gap and magnetic field. We analytically obtain the position and the height of the first and dominant peak of these oscillations.

Equations

$$M_{\mathrm{osc}}=\frac{2\mu }{\pi {\theta }_0}{\left(\frac{2\pi {\theta }_0}{B}\right)}^{\frac{3}{2}}\left(1-\frac{B}{2\pi {\theta }_0}\right)\sin \left(\frac{2\pi \mu }{B}-\pi \right)e^{-4\pi \frac{{\theta }_0}{B}},$$

$$\displaystyle H_{0}=\int_{\bf k}\bar{\psi}\left[\begin{array}[]{cc}\xi_{\bf k}&0\\ 0&\epsilon_{0}\end{array}\right]\psi\equiv\int_{\bf k}\bar{\psi}h\psi,$$

$$H_{\mathrm{int}}=U{\int }_x{\bar{\varphi }}_{\mathrm{d}}(x){\varphi }_{\mathrm{d}}(x){\bar{\varphi }}_{\mathrm{f}}(x){\varphi }_{\mathrm{f}}(x),$$

$$\displaystyle H_{\mathrm{int};1}=\int_{\bf k}\bar{\psi}\left[\begin{array}[]{cc}0&\theta\\ \bar{\theta}&0\end{array}\right]\psi\equiv\int_{\bf k}\bar{\psi}h_{\mathrm{int};1}\psi.$$

$${ϵ}_{\mathbf{k},\pm }=\frac{{\xi }_{\mathbf{k}}+{ϵ}_0}{2}\pm \sqrt{{\left(\frac{{\xi }_{\mathbf{k}}-{ϵ}_0}{2}\right)}^2+\bar{\theta }\theta },$$

$$\theta =U\theta {\int }_{\mathbf{k}}\frac{{\mathcal{F}}_{{ϵ}_{\mathbf{k},+}}-{\mathcal{F}}_{{ϵ}_{\mathbf{k},-}}}{{ϵ}_{\mathbf{k},+}-{ϵ}_{\mathbf{k},-}},$$

$$1=UB\sum _n\frac{{\mathcal{F}}_{{ϵ}_{n,+}}-{\mathcal{F}}_{{ϵ}_{n,-}}}{\sqrt{(Bn+\frac{B}{2}-\mu {)}^2+4\bar{\theta }\theta }},$$

$$\sum _{n=0}g(n)={\int }_0^{\infty }g(x)dx+\sum _{p\ne 0}{\int }_0^{\infty }{e}^{i2\pi px}g(x)dx,$$

$$\theta ={\theta }_0+2{\theta }_0\cos \left(2\pi \frac{\mu }{B}-\pi \right)e^{-\frac{4\pi {\theta }_0}{B}}\sqrt{\frac{B}{2\pi {\theta }_0}},$$

$$\delta F_{\mathrm{osc}}$$

$$\times \cos \left(\frac{2\pi \mu }{B}-\pi \right)e^{-\frac{4\pi {\theta }_0}{B}},$$

$$\theta ={\theta }_0(T)+2{\theta }_0(T)\cos \left(2\pi \frac{\mu }{B}-\pi \right)R(T),$$

$$R_n(T)={\int }_0^{+\infty }\frac{e^{i2\pi nx}}{x+\frac{1}{2}}{e}^{-\frac{2\pi {\theta }_0}{B}\left[\frac{1}{2\pi \frac{T}{{\theta }_0}\left(x+\frac{1}{2}\right)}+2\pi \frac{T}{{\theta }_0}\left(x+\frac{1}{2}\right)\right]}dx,$$

$$R(T)\approx R_0(T)\left[1+2\sum _{n=1}\cos \left(\frac{{\theta }_0}{T}\bar{z}n-\pi n\right)e^{-\frac{2{\pi }^2{n}^2}{|g(\bar{z})|}}\right],$$

$$\bar{z}=\frac{4\pi {\theta }_0}{B+\sqrt{B^2+(4\pi {\theta }_0{)}^2}},$$

$$T_{\mathrm{peak}}=\frac{4{\theta }_0^2}{B+\sqrt{B^2+{\left(4\pi {\theta }_0\right)}^2}}.$$

$$T_{\mathrm{peak}}=\frac{{\theta }_0}{\pi }.$$

$$M_{\mathrm{osc}}=$$

$$\displaystyle H_{0}=\int_{\bf k}\bar{\psi}\left[\begin{array}[]{cc}\xi_{\bf k}&0\\ 0&\epsilon_{0}\end{array}\right]\psi\equiv\int_{\bf k}\bar{\psi}h\psi,$$

$$\displaystyle\bar{\psi}=\left[\bar{\phi}_{\mathrm{d}},~{}\bar{\phi}_{\mathrm{f}}\right],~{}~{}~{}\psi=\left[\begin{array}[]{c}\phi_{\mathrm{d}}\\ \phi_{\mathrm{f}}\end{array}\right].$$

$$H_{\mathrm{int}}=U{\int }_x{\bar{\varphi }}_{\mathrm{d}}(x){\varphi }_{\mathrm{d}}(x){\bar{\varphi }}_{\mathrm{f}}(x){\varphi }_{\mathrm{f}}(x)=-U{\int }_x{\bar{\varphi }}_{\mathrm{d}}(x){\varphi }_{\mathrm{f}}(x){\bar{\varphi }}_{\mathrm{f}}(x){\varphi }_{\mathrm{d}}(x),$$

$$-\left(\bar{\theta }+{\bar{\varphi }}_{\mathrm{d}}{\varphi }_{\mathrm{f}}U\right)U^{-1}\left(\theta +U{\bar{\varphi }}_{\mathrm{f}}{\varphi }_{\mathrm{d}}\right)+U{\bar{\varphi }}_{\mathrm{d}}{\varphi }_{\mathrm{f}}{\bar{\varphi }}_{\mathrm{f}}{\varphi }_{\mathrm{d}}=-\bar{\theta }{U}^{-1}\theta -{\bar{\varphi }}_{\mathrm{d}}{\varphi }_{\mathrm{f}}\theta -\bar{\theta }{\bar{\varphi }}_{\mathrm{f}}{\varphi }_{\mathrm{d}}.$$

$$i{S}_0=i{\int }_{C}dt\bar{\psi }(x)(ϵ-h)\psi (x)$$

$$\equiv i{\int }_{\mathbf{r}}{\int }_{-\infty }^{+\infty }dt\hat{\bar{\Psi }}(x)(ϵ-{\hat{H}}_0)\hat{\Psi }(x),$$

$$\hat{\bar{\Psi }}=\hat{\bar{\psi }}{\hat{L}}^{-1},\hat{\Psi }=\hat{L}{\hat{\sigma }}_3\hat{\psi },$$

$$\displaystyle\hat{L}=\frac{1}{\sqrt{2}}\left[\begin{array}[]{cc}1&-1\\ 1&1\end{array}\right],~{}~{}~{}\hat{\sigma}_{3}=\left[\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right]$$

$$i{S}_{\mathrm{int};1}=-i{\int }_{C}dt{H}_{\mathrm{int};1}$$

$$=-i{\int }_{\mathbf{r}}{\int }_{-\infty }^{+\infty }dt\hat{\bar{\Psi }}(x)\left[\hat{\theta }{\tau }^{+}+\hat{\bar{\theta }}{\tau }^{-}\right]\hat{\Psi }(x)$$

$$\displaystyle\tau^{+}=\left[\begin{array}[]{cc}0&1\\ 0&0\end{array}\right],~{}~{}~{}\tau^{-}=\left[\begin{array}[]{cc}0&0\\ 1&0\end{array}\right]$$

$$\displaystyle\hat{\sigma}_{0}=\left[\begin{array}[]{cc}1&0\\ 0&1\end{array}\right],~{}~{}~{}\hat{\sigma}_{1}=\left[\begin{array}[]{cc}0&1\\ 1&0\end{array}\right]$$

$$i{S}_{\mathrm{int};2}=-\frac{i}{2}{\int }_x{\mathrm{Tr}}_{\mathrm{K}}\left[\hat{\bar{\theta }}(x)U^{-1}{\hat{\sigma }}_1\hat{\theta }(x)\right],$$

$$\displaystyle iS_{0}+iS_{\mathrm{int};1}=i\int_{\bf r}\int_{-\infty}^{+\infty}dt\left\{\hat{\bar{\Psi}}(x)\hat{\sigma}_{0}\left[\begin{array}[]{cc}\epsilon-\hat{\xi}_{\bf k}&-\theta\\ -\bar{\theta}&\epsilon-\epsilon_{0}\end{array}\right]\hat{\Psi}(x)+\hat{\bar{\Psi}}(x)\hat{\sigma}_{1}\left[\begin{array}[]{cc}0&-\theta_{\mathrm{q}}\\ -\bar{\theta}_{\mathrm{q}}&0\end{array}\right]\hat{\Psi}(x)\right\},$$

Full text

We are motivated by recent experiments where quantum oscillations in Kondo insulator SmB6 were observed LiScience2014 ; SuchitraScience2015 ; SuchitraNatPhys2017 ; SuchitraScience2020 . This material becomes insulating below about $T_{\mathrm{c}}=40\mathrm{K}$ due to the gapping of the Fermi surface. Observed quantum oscillations occur only in the magnetization, i.e. de Haas - van Alphen (dHvA) effect, while electric resistivity doesn’t show any oscillations (absence of Shubnikov - de Haas effect). The latter is due to after all insulating behavior, while occurance of dHvA oscillations in insulators poses a mystery. Not to mention that the observed frequency of dHvA oscillations of the insulator is that of the metallic phase of the system, before it turned insulating via a gapping mechanism. In addition, the frequency is similar to the one of the LaB6 metallic material, which never turns insulating and which has a similar band structure with the metallic phase of the SmB6.

There are theories which propose emergence of neutral quasiparticles with Fermi surface that somehow nevertheless couple to the magnetic field to show dHvA effect but not the electric field, or of fermions that fractionalize to either Majoranas or other exotic structures BaskaranArxiv ; SodemannChowdhurySenthil ; VarmaPRB . Whether these theories can explain the experiments LiScience2014 ; SuchitraScience2015 ; SuchitraNatPhys2017 ; SuchitraScience2020 is currently under debate. There are theories which obtain dHvA quantum oscillations in insulators based on the non-interacting fermions KnolleCooperPRL15 ; PalArxiv2022 . SmB6 insulator is known to be strongly correlated (see for example Hewson ), and it is rather odd to expect that the non-interacting picture will fully explain its physical properties.

In our theory the two fermions hybridize with each other via the repulsive Coulomb interaction. As a result, a gap opens up at their energy band crossing points resulting in an insulator when the Fermi level is at the crossing point. Free energy gets minimized as a result of the hybridization. The magnetic field quantizes the energy spectrum of fermions forming Landau levels. Then an intuitive picture of the effect of quantium oscillations in insulators is that when the Landau levels of two fermion branches cross each other, the hybridization is maximised, while when the Landau levels distance apart, the hybridization is minimized. To see this effect analytically, the hybridization must be treated self-consistently and obtained from the non-linear equation. Finally, since the free energy gain is proportional to the second power of the hybridization, the magnetization and specific heat are expected to be oscillating functions of the magnetic field. It turns out that the frequency of oscillations equals to the energy of the energy bands crossing point, and is related to the Fermi energy of the metallic phase before the hybridization (this fact was known since KnolleCooperPRL15 ). The response of the system due to the described mechanism is found to be diamagnetic (Landau diamagnetism) and the non-oscillating part of the magnetization is proportional to the magnetic field $B$. The oscillating part of the magnetization $M$ at small magnetic fields and at zero temperature can be best described with

[TABLE]

where $\mu$ is the Fermi energy of the metallic phase and $2\theta_{0}$ is the non-oscillating part of the hybridization gap, and where we keep $B$ in units of cyclotron frequency. Smallness of the magnetic field is given by the $8\pi\theta_{0}>B$ condition, which falls in in to the experimental regime. Despite of the exponential suppression, in reality of the experiments in which $2\theta_{0}\approx 40K$ with the effective mass picked to be $m=0.2m_{\mathrm{e}}$ (where $m_{\mathrm{e}}$ is bare electron’s mass), Eq. (1) gives sizable contribution at small magnetic fields. We plot oscillating part of the torque in Fig. (1). Having presented our main results for zero temperature, in the rest of the paper we outline analytical arguments which lead to them and expand the list of our results to the non-zero temperature. Each step of the calculations is outlined in details in the Supplemental Material (SM) SM .

We consider a model Hamiltonian keeping in mind a SmB6 system.

[TABLE]

where $\xi_{\bf k}=\frac{{\bf k}^{2}}{2m}-\mu$, notation $\int_{\bf k}(..)=\int\frac{d{\bf k}}{(2\pi)^{2}}(..)$ was used, then $\epsilon_{0}=E_{0}-\mu=\mathrm{const}$, which we will set $\epsilon_{0}=0$ without losing generality, where $\mu$ is the chemical potential, and $\bar{\psi}=\left[\bar{\phi}_{\mathrm{d}},~{}\bar{\phi}_{\mathrm{f}}\right],$ and the same for $\psi$. In order to present analytic arguments we assume that the system is two-dimensional. Three-dimensional case is slightly more analytically involved, but will also be described by the presented below arguments. We consider spinless fermions for simplicity, and note that all of the results below will not change if the $f-$ fermion’s spectrum had a dispersion. Interaction between the two types of fermions is

[TABLE]

where $U>0$ corresponds to repulsion. We decouple interaction Eq. (5) using the Hubbard-Stratonovich transformation, and keeping in mind the transformation is made with the action, we write the result, $-\left(\bar{\theta}+\bar{\phi}_{\mathrm{d}}\phi_{\mathrm{f}}U\right)U^{-1}\left(\theta+U\bar{\phi}_{\mathrm{f}}\phi_{\mathrm{d}}\right)+U\bar{\phi}_{\mathrm{d}}\phi_{\mathrm{f}}\bar{\phi}_{\mathrm{f}}\phi_{\mathrm{d}}=-\bar{\theta}U^{-1}\theta-\bar{\phi}_{\mathrm{d}}\phi_{\mathrm{f}}\theta-\bar{\theta}\bar{\phi}_{\mathrm{f}}\phi_{\mathrm{d}}.$ Bilinear in fermion fields term introduces hybridization between $d-$ and $f-$ fermions, described by the Hamiltonian

[TABLE]

Assuming that $\theta$ and $\bar{\theta}$ are constants in space (mean field approximation), the spectrum of fermions described by Hamiltonian $h+h_{\mathrm{int};1}$ becomes

[TABLE]

which are schematically plotted in Fig. (2). There is a gap opening at the cross-section of the energy bands of $d-$ and $f-$ fermions due to the hybridization.

Next, we need to find the structure of the hybridization. For that, as usual, we integrate our fermions and minimize the action with respect to the hybridization, and get self-consistent equation for the hybridization,

[TABLE]

where ${\cal F}_{\epsilon}=\tanh\left(\frac{\epsilon}{2T}\right)$ is the Fermi-Dirac distribution function where $T$ is the temperature.

A magnetic field perpendicular to the plane of the system is added. It results in the Landau quantization of the energy spectrum, such that $\xi_{\bf k}\rightarrow\xi_{n}=B(n+\frac{1}{2})-\mu$, where $n>0$ is an integer denoting corresponding Landau level. We note that here and everywhere below we count $B$ in terms of cyclotron energy $\omega_{B}=\frac{eB}{mc}$, where $m$ is the mass of $d-$ fermions. With that we rewrite Eq. (10) for the system in the applied magnetic field,

[TABLE]

where now summation is over the Landau levels. To sum over the Landau levels, we utilize Poisson summation formula

[TABLE]

where $g(n)$ is some function. We will first obtain all of the results for $T=0$. With that when the Fermi energy is in the gap, the upper band is deoccupied meaning ${\cal F}_{\epsilon_{n,+}}=1$, and the lower band is completely filled, ${\cal F}_{\epsilon_{n,-}}=-1$. As a result, the non-linear equation is solved with

[TABLE]

in the limit $8\pi\frac{\theta_{0}}{B}>1$, which allowed us to pick only $p=\pm 1$ in the Poisson summation formula because harmonics with $n>1$ decay with a $e^{-\frac{4\pi\theta_{0}}{B}n}$ factor. Here $\theta_{0}=\sqrt{\mu^{2}-\frac{B^{2}}{4}}e^{-\frac{1}{4U}}$ is the non-oscillating part of the solution. The free energy of the hybridized system reads $F=-\theta^{2}+F_{\Sigma}$, where $F_{\Sigma}$ is the part of the free energy corresponding to the oscillation in magnetic field of the non-interacting system, i.e. to the system with constant hybridization $\theta_{0}$. We note that it is not $-\frac{\theta^{2}}{2U}$ that enters the free energy as it gets cancelled troughout the derivation, but just like in the BCS theory it is rather $-\theta^{2}$. Collecting oscillation terms from both parts of free energy, i.e. from $-\theta^{2}$ and $F_{\Sigma}$, we get

[TABLE]

here $1$ in the round brackets is the first main result of the present paper, it originated from the oscillations of the hybridization and is due to the interactions. The term $-\frac{B}{2\pi\theta_{0}}$ in the round brackets is non-interacting, came from oscillating part of $F_{\Sigma}$ and corresponds to the results of Refs. KnolleCooperPRL15 ; PalArxiv2022 . Essentially, it corresonds to the oscillations in magnetic field of the fully filled fermion band $\epsilon_{{\bf k};-}$. It is clear that the latter term is smaller than the former at small magnetic fields defined by the $2\pi\frac{\theta_{0}}{B}>1$ condition. When $\frac{B}{2\pi\theta_{0}}=1$ the amplitude vanishes, and at larger fields the non-interacting part will become dominant. We note that vanishing AokiJPSJ2014 and a decrease Li_PRX2022 of the amplitude was observed in strongly correlated systems AokiJPSJ2014 ; Li_PRX2022 and attributed to the field driven metamagnetic transition. Proposed here mechanism might be relevant to the experiments AokiJPSJ2014 ; Li_PRX2022 , but more details are needed. We differentiate the oscillating part of minus the free energy Eq. (De Haas-van Alphen effect and quantum oscillations as a function of temperature in correlated insulators) with respect to the magnetic field and obtain Eq. (1) by keeping only the term $\propto\mu$.

We observe that the response of the system due to the hybridization mechanism is diamagnetic. This is precisely due to the magnetic field dependence of the non-oscillating part of the hybrdidization $\theta_{0}$ defined after Eq. (13). Indeed, the non-oscillating part of the magnetization is $M\propto-\partial_{B}(-\theta_{0}^{2})=-\frac{B}{2}e^{-\frac{1}{2U}}$. Physics of obtained diamagnetism is that the system wants to expel the magnetic field in order to minimize the energy associated with the Coulomb repulsion. Another point is that the calculations suggest that we could have had any quadratic in magnetic field structure under the square root in the definition of $\theta_{0}$, for example $\sqrt{\mu^{2}-\frac{B^{2}}{4}+aB}$, where $a$ is some energy, then there would be spontaneous magnetic order in the system characterized by a finite magnetization. This is experimentally not the case LiScience2014 , hence the choice of $a=0$. Moreover, we wish to point out that $-\frac{B^{2}}{4}$ in $\theta_{0}=\sqrt{\mu^{2}-\frac{B^{2}}{4}}e^{-\frac{1}{4U}}$ is due to the lowest Landau level.

Let us now discuss temperature dependence of the amplitude of the dHvA oscillations. The non-linear equation at finite temperature is given in Eq. (10), where the temperature is in the Fermi-Dirac distribution functions. With all the technical steps presented in the SM SM to the paper, we here breifly outline the steps. Eq. (10) is modified to the presence of the magnetic field, summation over the Landau levels is performed again with the help of the Poisson summation formula Eq. (12). The non-oscillating in magnetic field terms are barely changed at finite temperatures. However, the terms which result in dHvA effect are drastically modified by the non-zero temperature. The distribution functions introduce residues, and it turns out it is convenient to integrate over the Poisson variable by introducing a contour which embraces certain residues which ensure convergence, such that the integral over the Poisson variable is reduced to the sum over the residues. Then the sum over the residues is calculated by using the Poisson summation formula Eq. (12) for the second time. As a result of the summation an oscillation with the inverse temperature of the hybrdization appears. Then the obtained solution for the hybridization Eq. (13) at zero temperature gets updated in accord with the mentioned steps as

[TABLE]

where $R(T)=\sum_{n=-\infty}^{+\infty}R_{n}(T)$ is a result of the Poisson summation formula with elements

[TABLE]

where $x+\frac{1}{2}$ is reminiscent of the Matsubara frequency. We plot $R_{0}(T)$ in Fig. 3 to suggest that except for the shoulder it is consistent with standard Lishits-Kosevich formula.

It is now instructive to perform the steepest descent approximation in evaluating the $R_{0}(T)$, we get the following limiting expressions $R_{0}(T=0)=e^{-\frac{4\pi\theta_{0}}{B}}\sqrt{\frac{B}{2\theta_{0}}},$ and $R_{0}(T\rightarrow\infty)\approx\frac{B}{4\pi^{2}T}e^{-\frac{4\pi\theta_{0}}{B}}e^{-\pi\frac{\theta_{0}}{2H}\left(\frac{2\pi T}{\theta_{0}}\right)^{2}}.$ Indeed, the $T=0$ limit is consistent with the Eq. (13), and the $T\rightarrow\infty$ decay is consistent with the Lishits-Kosevich expression.

Having understood $R_{0}(T)$ let us now add higher harmonics. We work in the regime when the solution for the gap Eq. (15) is valid. We plot $R(T)$ as a function of temperature in Figs. (4) and (5). to describe shown dependence we can approximate the integral with the steepest descent method and obtain in the vicinity of the peak,

[TABLE]

where saddle-point is at

[TABLE]

and $g(\bar{z})=\left(\frac{2\pi T}{\theta_{0}}\right)^{2}\frac{1}{\bar{z}^{3}}\left(\bar{z}-\frac{4\pi\theta_{0}}{B}\right)$, and a first maximum in the temperature dependence of Eq. (17) occurs at $\frac{\theta_{0}}{T_{\mathrm{peak}}}\bar{z}=\pi$ which reads

[TABLE]

Let us first discuss the case of relatively small magnetic fields, $4\pi\theta_{0}>B$. In particular, in Fig. (4) we plot $R(T)$ for $2\theta_{0}=B$. The position of the peak defined in Eq. (19) can be approximated to be

[TABLE]

To justife our approximation, we plot $\propto\cos\left(\frac{\theta_{0}}{T}\right)$ together with $R(T)$ in Fig. (4). The approximation Eq. (17) works well in the vicinity of the peak, away from it one has to include more terms to make sure the function is positive-definite. As a proof of principle, when the parameter $\frac{2\theta_{0}}{B}$ increases, the amplitude of Eq. (17) exponentially decreases due to the decaying $e^{-4\pi\frac{\theta_{0}}{B}}$ term, however, with that more oscillation cycles appear. For example, see Fig. (5) for the plot for $\frac{2\theta_{0}}{B}=3$, where the position of the peak is consistent with analytical predictions Eq. (20). In addition, we analytically found that $R(T_{\mathrm{peak}})\approx 2e^{-\frac{4\pi\theta_{0}}{B}}$, which is indeed plotted in Fig. (4) and in the left plot in Fig. (5). Non-zero temperature expression for the magnetization reads as

[TABLE]

where $R(T)$ is given by primarily Eq. (16), and with Eqs. (17) approximating the behavior in the vicinity of the peak.

Let us now discuss relatively large magnetic fields $B\sim 4\pi\theta_{0}$ and $B>4\pi\theta_{0}$. We plot $R(T)$ in Fig. (5) right for $\frac{2\theta_{0}}{B}=0.2$. Unfortunately, the steepest descent method fails to work well in this regime. We, nevertheless, still can predict the existence of the peak in the amplitude due to the oscillations. From the right plot of Fig. (5) we observe that for $\frac{2\theta_{0}}{B}=0.2$ the position of the peak is at $\frac{2\pi T^{(\mathrm{plot})}_{\mathrm{peak}}}{\theta_{0}}\approx 1.32$ which should be compared with the corresponding quantity from the Eq. (19), which is $\frac{2\pi T^{(\mathrm{analytics})}_{\mathrm{peak}}}{\theta_{0}}\approx 0.96$. It turns out that for whatever reason a better approximation is $T^{*}_{\mathrm{peak}}=\frac{2\theta_{0}^{2}}{B}$ which is a high field limit of Fig. (5) and which gives $\frac{2\pi T^{*}_{\mathrm{peak}}}{\theta_{0}}\approx 1.25$. A peak in the amplitude of the dHvA oscillations was experimentally observed at small temperatures in SuchitraScience2015 ; SuchitraNatPhys2017 ; SuchitraScience2020 . In paricular of interest is Fig. 7c in SuchitraScience2020 , where a mismatch of Lifshits-Kosevich formula with the observed temperature dependence of the amplitude was observed. We think that the peak structure in $R(T)$ plotted in the right figure of Fig. (5) might explain the mismatch.

In addition to the discussed above, we find a new regime of oscillations occuring in rather strong magnetic fields, $B>8\pi\theta_{0}$. There in the course of the oscillations the hybridization gap vanishes at the minimum of the cycle and then emerges back at the peak. We speculate that such an oscillation can undergo with a phase separation to metallic and insulating domains. Unfortunately, it can’t be experimentally found in the SmB6 material due to the physical parameters of the system, however, InAs/GaSb quantum wells where value of the Coulomb interaction can be controlled and where quantum oscillations in the insulating regime were observed ExpInAs/GaSb1 ; ExpInAs/GaSb2 , are very promising candidates.

We note that in previous works KnolleCooperPRL15 ; PalArxiv2022 it has been numerically shown that non-interacting insulating fermion system shows a deviation from the Lifshits-Kosevich temperature dependence of the dHvA oscillation amplitude at small temperatures, and a peak in the amplitude has been numerically plotted in these papers. However, it is only in the present paper it was recognized that the deviation from Lifshits-Kosevich formula is actually an oscillation of the amplitude with inverse temperature, and the peak can be associated with the first cycle of the oscillation. Note that there are not that many instances of oscillations of physical observables with temperature. We are aware of one theoretical prediction of occurance of mesoscopic oscillations of conductance as a function of temperature in the disordered conductors due to the mesoscopic fluctuations SpivakZyuzinCobden . These oscillations are different from the ones predicted in this paper.

In regards of topology, here we were working with ordinary insulator. There are theoretical claims that SmB6 may be a topological insulator DzeroSunGalitskiColeman . Then, if instead we were to think of a topological insulator, we can fit our model to Volkov-Pankratov model of topological insulators SovietTI by promoting the hybridization gap to momentum and spin dependences as $\theta\rightarrow\theta({\bf k}\cdot{\bm{\sigma}})$, as well as adding $k-$ dependence to the mass in SovietTI . Although, the hybridization gap will now be proportional to the Landau level index, at the Fermi level it can be approximated as a constant, i.e. $\propto\theta Bn_{\mathrm{F}}$, and all the predictions of this paper hold. The only complication is to find such an interaction that will result in momentum and spin dependences of the hybridization, and solve the self-consistency equation with this interaction. If necessary, we leave the details to future research.

Finally, we became aware of a similar work AlloccaCooper which studied the hybridization self-consistently, just like it is done in the present paper. We chose to leave our scientific feedback on this work in the SM SM rather than here.

To conclude, we have studied de Haas-van Alphen oscillations in correlated insulators. We studied a mechanism when the osillations originate from the magnetic field dependence of the hybridization. The hybridization is due to the Coulomb repulsion and is treated self-consistently. It turned out that this mechanism is the dominant at low magnetic fields as compared to studied non-interacting mechanism KnolleCooperPRL15 ; PalArxiv2022 . The amplitude vanishes at a field $B=2\pi\theta_{0}$. In experiments AokiJPSJ2014 ; Li_PRX2022 a transition in the dHvA oscillations marked by vanishing in AokiJPSJ2014 and reduction in Li_PRX2022 of the amplitude was observed. Moreover, we showed that the amplitude of dHvA oscillations is itself a periodic function of inverse temperature. The position of the peak is analytically estimated and from where one can extract the value of the hybridization gap. In experiments SuchitraScience2015 ; SuchitraNatPhys2017 ; SuchitraScience2020 deviation from the Lifshits-Kosevich expression was observed at small temperatures. Our work might be relevant to these experiments.

Acknowledgements. The author thanks A.M. Finkel’stein, M.M. Glazov, A. Kamenev, and A.Yu. Zyuzin for helpful discussions. VAZ is grateful to Pirinem School of Theoretical Physics and to Weizmann Institute of Science for hospitality during the Summer of 2022, and in particular to Y. Oreg and A. Stern for support. VAZ is supported by the Foundation for the Advancement of Theoretical Physics and Mathematics BASIS.

Supplemental Material to ”De Haas-van Alphen effect and quantum oscillations as a function of temperature in correlated insulators”

.1 Model

[TABLE]

where $\xi_{\bf k}=\frac{{\bf k}^{2}}{2m}-\mu$ and $\epsilon_{0}=E_{0}-\mu=\mathrm{const}$ (we will later assume $\epsilon_{0}=0$), where $\mu$ is the chemical potential. Arguments presented below will also be valid for the inverted band structure, where $\epsilon_{0}\rightarrow-\xi_{\bf k}$, as well as for $\epsilon_{0}\rightarrow\frac{{\bf k}^{2}}{2m_{2}}-\mu_{2}$ spectrum. Spinors are defined as

[TABLE]

Interaction between the two types of fermions is

[TABLE]

where $U>0$ corresponds to repulsion. Hubbard-Stratonovich transformation of the action reads

[TABLE]

In order to obtain self-consistent equation for the hybridization we employ the Keldysh technique (for example see Ref. [1SM]), although, the task at hand can be formulated in Matsubara space. It is our personal choice to work in Keldysh technique. In the Main Text in order not to burden the reader we don’t even mention that we formulated our theory in the Keldysh space.

[TABLE]

where we have performed Larkin-Ovchinnikov rotation defining

[TABLE]

where $\hat{\bar{\psi}}=(\bar{\psi}_{+},~{}\bar{\psi}_{-})$, and

[TABLE]

are matrices acting in Keldysh space.

[TABLE]

where

[TABLE]

are matrices acting in d-f space, and

[TABLE]

are Pauli matrices acting in Keldysh space. We have defined classical and quantum combinations of the interaction fields, namely $\theta_{\mathrm{cl}/\mathrm{q}}=\frac{1}{2}(\theta_{+}\pm\theta_{-})$ and $\bar{\theta}_{\mathrm{cl}/\mathrm{q}}=\frac{1}{2}(\bar{\theta}_{+}\pm\bar{\theta}_{-})$. For convenience we will be calling $\theta_{\mathrm{cl}}\equiv\theta$ and $\bar{\theta}_{\mathrm{cl}}\equiv\bar{\theta}$.

[TABLE]

where $\mathrm{Tr}_{\mathrm{K}}$ is the trace over the Keldysh components of the matrices.

Action describing the fermions is

[TABLE]

where $\hat{\xi}_{\bf k}$ stands for operator acting on the coordinate. Whenever, Fourier transformation is performed, $\hat{\xi}_{\bf k}\rightarrow\xi_{\bf k}=\frac{{\bf k}^{2}}{2m}-\mu$. Integrating fermions out, we get

[TABLE]

where $\hat{H}_{\mathrm{int};1;q}$ contains quantum componet of the interaction with the Hubbard-Stratonovich field,

[TABLE]

We recall,

[TABLE]

where

[TABLE]

and

[TABLE]

where ${\cal F}_{\epsilon}=\tanh\left(\frac{\epsilon}{2T}\right)$ is the Fermi-Dirac distribution function with $T$ being the temperature. Keldysh part of the Green function is

[TABLE]

Details of the Green function are

[TABLE]

where

[TABLE]

is the spectrum. Then, for example,

[TABLE]

such that

[TABLE]

From where we derive self-consistent equation for the hybridization. For that we vary the action over the $\bar{\theta}_{\mathrm{q}}$ and set it to zero,

[TABLE]

and we get

[TABLE]

where we have

[TABLE]

Simple check shows ${\cal F}_{\epsilon_{{\bf k},+}}-{\cal F}_{\epsilon_{{\bf k},-}}>0$ and $\epsilon_{{\bf k},+}-\epsilon_{{\bf k},-}>0$ such that the equation can be satisfied with a non-trivial $\theta$.

.2 Zero temperature

In magnetic field and at $T=0$ (also recall we are working in two-dimensions) the equation reads

[TABLE]

where factor of $2$ is due to the sum of the distribution functions.

We use Poisson summation formula

[TABLE]

in order to estimate the sum. Non-oscillating term is

[TABLE]

where $\Lambda$ is the upper cut-off, $f=\frac{\mu}{B}-\frac{1}{2}$, and we assumed $b<\Lambda-f$ and $b<f$. The other integral requires some care. We wish to study the lowest harmonic of oscillation and pick $p=\pm 1$ from the sum. Here is the step-by-step calculation,

[TABLE]

Consider now

[TABLE]

over the contour $C_{+}$ shown in the Fig. (1).

Then

[TABLE]

and the integrals over the arcs vanished due to their exponential decay as the radius goes to infinity. In the expression above notation $\pm 0$ should be understood from which side of the real axis, positive or negative, the variable $z$ approaches zero. This is needed to pick the correct list of the multivalued square root. Therefore,

[TABLE]

Therefore,

[TABLE]

[TABLE]

where we remind that $b=\frac{2\theta}{B}$.

.2.1 $8\pi\frac{\theta}{B}>1$ case

In this case we expand the logarithm as

[TABLE]

This is now a small correction to the non-linear equation Eq. (74),

[TABLE]

and can be treated perturbatively. Assume we found

[TABLE]

which is

[TABLE]

where $fB=\mu-\frac{B}{2}$, and where we set $\Lambda-\mu=\mu$ in order to make sure that the magnetization is absent when the magnetic field is zero, i.e. there is no spontaneous magnetic order in the system. We comment on that in more details in the subsection C.

In order to find oscillating part of the hybridization we perturb $\theta=\theta_{0}+\delta\theta$ and get

[TABLE]

and the solution reads

[TABLE]

.2.2 $8\pi\frac{\theta}{B}\sim 1$ case

For example, we keep in mind a situation when $b=0.1$, then $4\pi b=8\pi\frac{\theta}{B}\approx 1.2$. Then, we estimate last term in the brackets in the right-hand side of Eq. (74) to be

[TABLE]

which is still a small correction as compared to the first term in the brackets of the right-hand side of Eq. (74). We just can not expand the logarithm anymore. We get the solution

[TABLE]

.2.3 $8\pi\frac{\theta}{B}<1$ case

The non-linear equation reads

[TABLE]

where, recall $b=\frac{2\theta}{B}$ and $fB=\mu-\frac{B}{2}$, and $B\Lambda=2\mu$ to make sure there is no spontaneous magnetization in the system. We then derive a solution for the hybridization,

[TABLE]

When $\cos\left(2\pi\frac{\mu}{B}-\pi\right)=0$, the gap is

[TABLE]

Moreover, it is clear that the solution is on the boarderline of its existance. Namely, when the sign of the argument of the natural exponent changes sign to positive, the solution no longer obeys $8\pi\frac{\theta}{B}<1$ and $\theta<\mu$ smallness conditions which were assumed in the course of the derivation. Even next harmonics of the oscillation which were omitted in the calculation can’t help the argument to not change the sign to positive. Hence, large magnetic fields causes a phase transition by either closing the gap or changing the solution to another limit, i.e. to $8\pi\frac{\theta}{B}>1$ case, when the oscillations are rather suppressed and the argument changes sign back to negative. It is quite likely this process is accompanied with a phase separation of the two possibilities. We anticipate that the oscillations are going to be antisymmetric with respect to highest and lowest phases of the oscillation, besides the existence of the higher harmonics which become less exponentially suppressed.

.3 Free energy and magnetization

Here we wish to derive the free energy $F$ within the Keldysh technique. We will be following Ref. [2SM].

[TABLE]

where $S=\int dT\frac{1}{T}\frac{dU}{dT}$ is entropy, and the internal energy $U$ is typically $U=\int_{\epsilon}{\cal F}_{\epsilon}A_{\epsilon}$, where $A_{\epsilon}$ is a function of $\epsilon$ and ${\cal F}_{\epsilon}=\tanh\left(\frac{\epsilon}{2T}\right)$ is the distribution function, then

[TABLE]

and therefore,

[TABLE]

Here $1$ in the square brackets should be dropped in order to define a proper free energy, hence the right arrow. Then we can define the free energy by figuring out the internal energy, which can be done within the Keldysh technique. For that, following the lines of Ref. [2SM], we introduce heat density sources in to the action.

[TABLE]

where $\hat{\eta}=\eta_{1}\hat{\sigma}_{0}+\eta_{2}\hat{\sigma}_{1}$ are the heat density sources, and $\hat{\xi}_{\bf k}$ stands for operator acting on the coordinate. Heat density is given by varying the partition functino with respect to the source,

[TABLE]

Internal energy $U$ is defined by this average. Let us find out what $A_{\epsilon}$ defined in Eq. (88) equals to in our model. After fermions are integrated out

[TABLE]

where $\hat{G}$ and $\hat{h}_{\mathrm{int};q}$ were defined in the subsection .1. Then

[TABLE]

where in the last equality sign we utilized self-consistent gap equation. Another term reads

[TABLE]

Summing the two we get

[TABLE]

which gives the correct expression for the free energy upon calculating entropy.

We study the free energy at $T=0$, it is given by

[TABLE]

The sum over the Landau levels in the second term will be made again with the help of Poisson summation formula Eq. (61) .

[TABLE]

We calculate

[TABLE]

where, recall, $\Lambda B=2\mu$ and $fB=\mu-\frac{B}{2}$. This calculation of the free energy is in line with the BCS theory, i.e the energy gain is equal to $-\theta^{2}$. We note that the response from the non-interacting term (when hybridization is treated as a constant rather than self-consistently) is diamagnetic, although it is not immediately clear from Eq. (112).

Now the sum over the harmonics. We again pick only $p=\pm 1$ from the sum. Then we are left with the integral, which in the small magnetic field regime, i.e. $\frac{1}{2\pi b}<1$, is estimated to be

[TABLE]

where by $\rightarrow$ we extracted only oscillating part, and we remind $b=\frac{2\theta_{0}}{B}$. We used integration by parts once in order to get rid of divergence at infinity. After that we utilized the same contour integration as was used in Eq. (66). Therefore, for $4\pi b>1$ we get for the oscillating part of free energy,

[TABLE]

Oscillating term from $-\theta^{2}$ is obtained by plugging in Eq. (80) there,

[TABLE]

Here is why we have been choosing $\Lambda B=2\mu$ in all of our calculations. The non-oscillating part of the free energy obtained from Eq. (115) reads $-\theta_{0}^{2}=-\left(\mu^{2}-\frac{B^{2}}{4}\right)e^{-\frac{1}{2U}}$, and hence there is no term here propotional to the magnetic field in free energy. If there was a term propotional to the magnetic field, the system would have had a spontaneous magnetization, which is experimentally not the case. We also conclude from this argument that the response of the system is diamagnetic, because non-oscillating part of the magnetization is $M_{0}=-\partial_{B}(-\theta_{0}^{2})=-\frac{B}{2}e^{-\frac{1}{2U}}$. The system wants to expel the magnetic field such that the Coulomb energy of the system is minimized.

We compare now the two oscillating terms in the free energy, the one came from Eq. (112), which is due to the Coulomb interaction, and the other from the non-interacting picture given by Eq. (114),

[TABLE]

Oscillating part of the magnetization reads

[TABLE]

where only the $\propto\mu$ term was kept because it is the largest parameter here. Clearly, the non-interacting term which is second in the brackets is subleading in the case of small magnetic fields, $\frac{B}{2\pi\theta_{0}}<1$. Recall, that we have assumed $\frac{B}{8\pi\theta_{0}}<1$ in our derivations. The amplitude vanishes at $\frac{B}{2\pi\theta_{0}}=1$, and at higher fields the non-interacting term will dominate. It is not surprising as in this limit the field becomes larger than the hybridization gap and the system can be thought of as metallic. Even so, in this limit, $\frac{B}{8\pi\theta_{0}}>1$, we have predicted formation of the orbital domains due to the Coulomb interaction.

.3.1 Comment on the Ref. [3SM] A. Allocca and N. Cooper, SciPost Phys. 12, 123 (2022)

Although self-consistent equation on the hybdridization was utilized in Ref. [3SM] and there is a correspondence of their intermediate steps with the steps made in the present paper till the Eq. (13), however, further predictions and conclusions reached in Ref. [3SM] are different from the ones obtained in the present paper. Importantly, contrary to what is found here, the contribution to the dHvA effect from the magnetic field dependence of the hybridization is claimed in Ref. [3SM] to have the same field dependence as the non-interacting part found in Ref. KnolleCooperPRL15 of the Main Text. Whereas in the present paper it is found to be the dominant contribution at small magnetic fields as compared to the non-interacting part. Moreover, present paper predicts that the amplitude of the dHvA oscillations vanishes at $B=2\pi\theta_{0}$ (similar to beating). Ref. [3SM] does not have any discussion on the temperature dependence of the effect. In addition, the higher harmonics of dHvA osillation are found here to be exponentially suppressed as expected, contrary to what is claimed in Ref. [3SM].

.4 Temperature dependence

Here we study temperature dependence of oscillations obtained in previous subsections. For that we update the non-linear equation to account for the temperature,

[TABLE]

We follow the same steps as were in $T=0$ case. In case of finite temperature the integral Eq. (62) gets modified as

[TABLE]

This is because ${\cal F}_{x}=1-\frac{2}{e^{\frac{x}{T}}+1}$. Here, again $b=\frac{2\theta_{0}^{2}}{B}$ and $f=\frac{\mu}{B}-\frac{1}{2}$. We plot temperature dependence of the last two terms in the square brackets in the right hand side of Eq. (119) in Fig. In addition, let us do some analytics

[TABLE]

where we kept only Eq. (121) after integrating by parts, as the other is less singular at $z=0$, and assumed $\Lambda-f\approx f$ in the distribution function when taking the upper and lower limits. This assumption is valid at small temperatures, when the derivative over the distribution function is almost a delta function. The upper and lower limits are fixed and are never equal to infinity.

Now, the most interesting part. The integral corresponding to oscillations at zero temperature Eq. (66) gets modified as

[TABLE]

Let us work with the temperature dependent parts only. We define the $+1$ harmonic to be

[TABLE]

where we made $z=b\sinh(y)$ change of variable, which resulted in $dz=b\cosh(y)dy$, $\sqrt{z^{2}+b^{2}}=b\cosh(y)$ and $-z+\sqrt{z^{2}+b^{2}}=be^{-y}$ while $z+\sqrt{z^{2}+b^{2}}=be^{y}$, and $\frac{Bb}{2T}=\frac{\theta_{0}}{T}$. Harmonic with $-1$ is obtained from Eq. (126) by complex conjugation, i.e. ${\cal R}^{(-1)}_{\mathrm{osc}}(T)=\left[{\cal R}^{(+1)}_{\mathrm{osc}}(T)\right]^{*}$.

We introduce contour integration shown in Fig. (3) and obtain the expression for $W_{-}$ defined in Eq. (127),

[TABLE]

We found it better to come back to $\tanh$ in order to take care of the divergence at infinity. The chosen contour of integration shown in Fig. (3) is such that the integral over the arc converges because of the exponent in the numerator of the integrand in Eq. (129). The residues of the integrand of $W_{-}$ are at $y=y_{n}$ and defined as

[TABLE]

and according to the contour Fig. (3) they are with $n<0$, such that $\mathrm{Im}[y_{n}]>0$ and real part can be of any sign. Then

[TABLE]

Then

[TABLE]

It is not clear how to analytically sum up the series above, however, in order to gain some insight, we can approximate the series by applying the Poisson summation formula and performing the steepest descent method. We can then numerically plot the result of summation and compare it with the approximate analytics. Before doing it, as a remark, let us contrast with the standard Lifshits-Kosevich expression for the temperature dependence.

.4.1 Poisson summation and steepest descent method

To sum the series Eq. (137) we again employ Poisson summation formula Eq. (61).

[TABLE]

The non-oscillating part of the formula reads

[TABLE]

we numerically estimate it and plot it as a function of temperature in Fig. (4).

We can estimate the integral using steepest descent method,

[TABLE]

where we defined

[TABLE]

In Fig. (5) we plot justification for the steepest descent method. The original curve is plotted in blue, while the curve corresponding to the steepest descent approximation is shown in red. This figure suggests that the position of the extremum is captured correctly, but for large magnetic fields (right plot in Fig. (5)) the approximation only qualitatively describes the integral. We will show in the following that the position of the extremum defines period of quantum oscillations with inverse temperature.

(i) We first consider $2\pi b>1$ case. The steepest descent approximation is shown in the left plot in Fig. (5). Then $\bar{z}\approx 1$ and $|g(1)|=2\pi b\left(\frac{2\pi T}{\theta_{0}}\right)^{2}$, and we get

[TABLE]

which we plot in Fig. (4). At $T=0$ the integral over $x$ can be approximated as Gaussian, i.e. from $-\infty$ to $+\infty$ in our case by shifting $x\rightarrow x+\frac{1}{2}-\frac{\theta_{0}}{2\pi T}$ and setting $-\infty$ in the lower limit, then

[TABLE]

At $T\rightarrow\infty$

[TABLE]

where is the complimentary error function, whose expansion at large values of argument we used.

Now we pick the $\pm n$ harmonics in the formula Eq. (61),

[TABLE]

Since steepest descent method gives $\bar{z}=1$, we predict that there is a peak in $R(T)$ which occurs at

[TABLE]

In Fig. (6) we plot $R(T)$ for $b=1$, $b=3$ and $b=4$. There the peak indeed occurs at a value determined in Eq. (152).

Let us estimate the height of the peak. Numerics suggest, see Fig. (6), that the peak has a value of almost $2e^{-2\pi b}$, let us show it is indeed the case. We found it convenient to not perform the steepest descent for this task, but rather do the following transformations after application of the Poisson summation formula,

[TABLE]

now goes our approximation, because $\cosh(z)$ has a minimum at $z=0$ (reminiscent of steepest descent approximation made above), $e^{z}\sim 1+z$ in the vicinity of $z=0$. Therefore,

[TABLE]

For $b=1$ we have $4e^{-2\pi b*0.5}\approx 0.17$ and the other terms in the series are practically zero. For $b=3$ we have $4e^{-2\pi b*0.5}\approx 3*10^{-4}$ which is also practically zero. Therefore, we conclude that

[TABLE]

which is indeed numerically plotted in Fig. (6), recalling that $b=\frac{2\theta_{0}}{B}$. This result, we remind, is valid for $2\pi b=\frac{4\pi\theta_{0}}{B}>1$.

(ii) Let us now consider $2\pi b\sim 1$ and $2\pi b<1$ cases.

In particular we want to first focus on $b=0.2$ case as we think it might be relevant to the experiment. For this value the steepest descent saddle point is defined by the most general expression,

[TABLE]

Calculations and Fig. (5b) suggest we have to be careful with the steepest descent and take in to account higher orders of approximation. In particular, in Fig. (5b) we see that the position of the saddle point is captured well with the approximation, however, the tail isn’t. The tail is going to contribute to the amplitude of the integral.

[TABLE]

the peak occurs at

[TABLE]

which for small values of $\frac{4\pi\theta_{0}}{B}$ becomes $T_{\mathrm{peak}}=\frac{2\theta_{0}^{2}}{B}$. Let us check that our approximation works. From right plots of Fig. (7) we observe that for $\frac{2\theta_{0}}{B}=0.2$ the position of the peak is at $\frac{2\pi T^{(\mathrm{plot})}_{\mathrm{peak}}}{\theta_{0}}\approx 1.32$ which should be compared with the corresponding quantity from the Eq. (164), which is $\frac{2\pi T^{(\mathrm{analytics})}_{\mathrm{peak}}}{\theta_{0}}\approx 0.96$. It turns out that for whatever reason better approximation is $T^{*}_{\mathrm{peak}}=\frac{2\theta_{0}^{2}}{B}$, which gives $\frac{2\pi T^{*}_{\mathrm{peak}}}{\theta_{0}}=2\pi b\approx 1.25$. For $\frac{2\theta_{0}}{B}=0.1$ we have $\frac{2\pi T^{(\mathrm{plot})}_{\mathrm{peak}}}{\theta_{0}}\approx 0.78$ which is compared with $\frac{2\pi T^{(\mathrm{analytics})}_{\mathrm{peak}}}{\theta_{0}}\approx 0.57$ and $\frac{2\pi T^{*}_{\mathrm{peak}}}{\theta_{0}}=2\pi b\approx 0.62$. We conclude that the steepest descent approximation is a poor approximation for $\frac{4\pi\theta_{0}}{B}<1$ case.

.4.2 Conventional Lifshits-Kosevich expression for the temperature dependence

As a remark, let us contrast our calculation with the standard Lifshits-Kosevich expression for the temperature dependence of dHvA oscillation of a metal. There integral defining the temperature dependence of the dHvA oscillation amplitude instead of Eqs. (127), (129), and (137) equals (for example see Ref. [4SM], E.M. Lifshits and L.P. Pitaevskii “Statistical Physics, Part 2: Course of Theoretical Physics - Vol. 9”) to

[TABLE]

where $\alpha>0$. The contour of integration is the same as in Fig. (3), and the residues are at $z_{n}=2\pi i\left(n+\frac{1}{2}\right)$ with $n>0$. Then,

[TABLE]

where we used

[TABLE]

We see that there is no quantum oscillations as a function of temperature in metallic case.

.5 References in Supplemental Material (SM)

A. Kamenev, Field theory of non-equilibrium systems (Cambridge, University Press, 2012). 2. 2.

G. Schwiete and A.M. Finkel’stein PRB 90, 155441 (2014). 3. 3.

A. Allocca and N. Cooper, SciPost Phys. 12, 123 (2022). 4. 4.

E.M. Lifshits and L.P. Pitaevskii “Statistical Physics, Part 2: Course of Theoretical Physics - Vol. 9” Elsevier 2014.