Orbital magnetic susceptibility of type-I, II, and III massless Dirac fermions in two dimensions

In general, orbital magnetic susceptibility is given by [3]

where the summation is over the fermion Matsubara frequency, $\varepsilon_{n}=(2n+1)\pi k_{\rm B}T$ and the wave number $\bm{k}$. Note that we neglect the spin degrees of freedom for simplicity. In other words, the results shown in this paper is the orbital magnetic susceptibility per spin degrees of freedom. $V$ is the volume of the system (or the area for two dimensions). Tr is to take the trace over the Bloch bands and $\mathcal{G}$ is the abbreviation of the matrix form of thermal Green’s function ${\mathcal{G}}(\bm{k},i\varepsilon_{n})$ in the Bloch representation. The velocity $v^{i}_{\bm{k}}$ ($i=x,y$) is defined in the matrix form as $v^{i}_{\bm{k}}=\frac{\partial\mathcal{H}_{\bm{k}}}{\partial k_{i}}$, $e$ is the elementary charge, $\hbar$ is the reduced Planck constant, and $k_{\rm B}$ is the Boltzmann constant.

In the following, we consider two-band models whose Bloch Hamiltonians are generically given as

where $\bm{h}_{\bm{k}}=(h^{x}_{\bm{k}},h^{y}_{\bm{k}},h^{z}_{\bm{k}})$, and $\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ are the Pauli matrices. The dispersion relation of this Hamiltonian is given by $\varepsilon_{\bm{k},\pm}=h^{0}_{\bm{k}}\pm|\bm{h}_{\bm{k}}|$. Then the band touching occurs at which $|\bm{h}_{\bm{k}}|=0$ is satisfied. $h_{\bm{k}}^{0}$ gives tilting.

In this Hamiltonian, we obtain thermal Green’s function as

with $i\tilde{\varepsilon}_{n}:=i\varepsilon_{n}-h_{\bm{k}}^{0}+\mu$ ($\mu$ is the chemical potential) and $D=(i\tilde{\varepsilon}_{n})^{2}-|\bm{h}_{\bm{k}}|^{2}$. The velocities $v^{i}_{\bm{k}}$ ($i=x,y$) are given by

Before discussing the lattice model, we consider a single massless Dirac Hamiltonian in two-dimension with tilting. In this case, we use

where $v_{\rm F}$ is the Fermi velocity and the parameter $\alpha$ represents the degree of tilting in the $k_{x}$-direction. The dispersion relation is given by $\varepsilon_{\bm{k},\pm}=-\alpha v_{\rm F}\hbar k_{x}\pm v_{\rm F}\hbar k$ with $k=|\bm{k}|$. The type I (II) corresponds to $0\leq\alpha<1$ ($\alpha>1$) and The type III corresponds to $\alpha=1$ where the Dirac cone touches the $k_{x}$-$k_{y}$ plane. In the case without tilting ($\alpha=0$), it is well known that at $T=0$ the orbital magnetic susceptibility $\chi$ has a delta function peak with negative sign at $\mu=0$. At finite temperature $\chi$ is proportional to $-1/T$ [16, 17].

In the Hamiltonian (6), we obtain

Taking the trace in Eq. (1) and using the variables $x:=v_{\rm F}\hbar k_{x}$ and $y:=v_{\rm F}\hbar k_{y}$, we obtain the orbital magnetic susceptibility as

where $i\tilde{\varepsilon}_{n}:=i\varepsilon_{n}+\alpha x+\mu$ and $D=(i\tilde{\varepsilon}_{n})^{2}-x^{2}-y^{2}$. Then, after we carry out the integration by parts, $\chi$ becomes

From this, one may think that $\chi$ might change its sign when $\alpha$ becomes larger than 1 (i.e., the type II). However, it is not so simple because of the presence of the Fermi surface that extends to infinity in the $k_{x}$-$k_{y}$ plane in the type-II and type-III.

After a detailed calculation shown in Appendix A, we obtain

Rather surprisingly, this shows that $\chi$ is exactly zero in the type II Dirac Hamiltonian. [As we show in the following sections, $\chi$ has finite but small value in the type-II lattice Hamiltonian where the Fermi surface closes.] For $0\leq\alpha<1$, taking the summation over Matsubara frequency with use of $F(z)=1/(e^{\beta z}+1)$ with $\beta=1/(k_{\rm B}T)$, we obtain

where $f(\epsilon)=1/(e^{\beta(\epsilon-\mu)}+1)$ is the Fermi-Dirac distribution function and we have taken the residue at $z=\mu$. As $T$ goes to zero, $f^{\prime}(0)$ approaches $-\delta(\mu)$ as in the well-known result. On the other hand, when $\mu=0$ and at finite temperatures, $\chi$ is proportional to $-1/T$ as expected.

For the present models, we calculate the orbital magnetic susceptibility numerically. As for the summation over the Matsubara frequency, we use the standard method to replace it to the energy integral with Fermi distribution function as

with

Here $G_{\bm{k}}(\epsilon)=(\epsilon+i\Gamma-\mathcal{H}_{\bm{k}})^{-1}$ is the retarded Green’s function. Note that $\Gamma$ corresponds to the damping rate. Note also that $\frac{\partial^{2}\mathcal{H}_{\bm{k}}}{\partial k_{x}\partial k_{y}}=0$ holds for the models (i)-(iii). Several groups have shown that the orbital magnetic susceptibility in the lattice model has the form [28]

or [6]

instead of $\Theta(\epsilon)$ in Eq. (17), where $v^{ij}_{\bm{k}}=\frac{\partial^{2}\mathcal{H}_{\bm{k}}}{\partial k_{i}\partial k_{j}}$. In the present case where $\frac{\partial^{2}\mathcal{H}_{\bm{k}}}{\partial k_{x}\partial k_{y}}=0$ holds, it is apparent that the results using Eq. (LABEL:eq:Gomez_formula) is the same as those using Eq. (17). We also find that Eq. (23) is equivalent to Eq. (17) in the present case as far as the boundary contributions in the integration by parts vanish. (See the appendix B for details.) Therefore, the above three formulae for $\tilde{\Theta}(\epsilon)$ give the same results.

Figures 3(a), 3(b), and 3(c) show the $\mu$ dependence of the orbital magnetic susceptibility for the models (i), (ii), and (iii), respectively. Here we set $T=0$ and $\Gamma=0.06$. Focusing on the Dirac point ($\mu=0$), all three models exhibit qualitatively similar behavior. Specifically, we see that the type-I Dirac fermions exhibit divergingly large diamagnetic susceptibility, which coincides with the result in the previous section for the continuum model and the previous works [16, 17, 18, 27, 28, 6, 29]. Meanwhile, the type-II Dirac fermions do not show divergingly large susceptibility. It has a small positive value and the $\mu$ dependence around the Dirac point is very small, which is in sharp contrast to the type-I case where the steep drop of $\chi$ around the Dirac point is seen. Note that in the continuum model in Sec. II.2, the orbital magnetic susceptibility vanishes, meaning that the finite values of $\chi$ obtained here originate from the lattice model. As for the type-III case, we see the diamagnetic susceptibility, but its value is much smaller than that for the type-I. Comparing among three models, we see that the ratio between the orbital magnetic susceptibility of the type-III and the type-I becomes smaller as the number of the Dirac cones increases.

Away from the Dirac points, we see several peaks of $\chi$, all of which are positive-signed (i.e., paramagnetic). We elaborate on the origin of these peaks in the next section.

Figure 4 shows the temperature dependence of $\chi$ with $\mu=0$. Again, qualitative behavior is common among the three models. We see that only the type-I Dirac fermions exhibit steep temperature dependence, namely, the diamagnetic susceptibility becomes larger as $T\rightarrow 0$. In contrast, $\chi$ exhibits a small temperature dependence for the types-II and III Dirac fermions.

Summarizing these results, we find that only the type-I Dirac fermions show characteristic divergingly large diamagnetism at the Dirac points, which has a sharp dependence on $\mu$ and $T$, whereas the type-II Dirac fermions show small paramagnetism which is rather insensitive to $\mu$ and $T$. The type-III fermions show the diamagnetism which is moderately dependent on $\mu$ but rather insensitive to $T$.

Unlike the type-I case, in the type-II case, the Fermi surface extends the region away from the Dirac points. Thus, here we argue whether the Dirac points play essential roles for the orbital magnetic susceptibility when the chemical potential is at the Dirac point. To this end, we compute the momentum-resolved contributions to the orbital magnetic susceptibility ²²2Note that $\bar{\chi}_{\bm{k}}$ is not divided by the volume.:

with

First, we examine how the orbital magnetic susceptibility vanishes in the continuum model in the type-II case. We calculate $\bar{\chi}_{\bm{k}}$ using Eqs. (6) and (8) with $\alpha=2$ at $\mu=0$ with $\Gamma=0.06$. As shown in Fig. 5, we can see that there are delta-function-like positive peaks maximized at the Dirac point along $\alpha\cos\theta\pm 1=0$ [$\theta=\tan^{-1}(k_{y}/k_{x})$], which correspond to the Fermi surfaces of the continuum model. On the other hand, in the region where both $\alpha\cos\theta+1>0$ and $\alpha\cos\theta-1<0$ hold (i.e., the region where only the lower band is occupied), there are negative contributions. As pointed out in Eq. (11), the orbital magnetic susceptibility exactly vanishes after performing the integration over $\bm{k}$. Therefore, we expect that the cancellation occurs between the large positive contribution at the Fermi surface and the negative contributions from the other region in the continuum model.

Let us proceed to the lattice models. In Fig. 6, we plot $\bar{\chi}_{\bm{k}}$ for the type-II case at $\mu=0$. The upper panels are the map for the entire Brillouin zone. We see that the main contributions originate from the Fermi surface. The middle panels are the zoom-up near a Dirac point. We see that $\bar{\chi}_{\bm{k}}$ sharply oscillates along the direction perpendicular to the tilting of the cone. To be specific, $\bar{\chi}_{\bm{k}}$ has a large positive value right at the Dirac point, whereas it has a large positive value along the Fermi surfaces and a large negative value in the two areas out of four areas away from the Dirac point. For closer look at this, in the lower panels, we classify the momentum space with respect to the band occupation. There are classified into three types as indicated by the colors. Specifically, the orange (green) dots represent the region where both upper and the lower bands are occupied (unoccupied), and the purple dots represent the region where only the lower band is occupied. As mentioned above, we see that the positive contribution comes from the Fermi surface, maximized at the Dirac point, while the large negative value comes from the region represented by the purple dots. These behaviors are qualitatively the same behavior as in the continuum model shown in Fig. 5. However, the exact cancellation does not occur in the lattice models. As a result of this subtle cancellation, the susceptibility for the type-II case at $\mu=0$ is moderately paramagnetic [Fig. 3]. This is in contrast with the type-I Dirac fermions that contribute divergingly to the diamagnetic susceptibility.

In Sec. I, we have addressed the importance of the interband effect in the orbital magnetic susceptibility. In this regard, it is worth interpreting our results in terms of the intraband versus interband contributions. To extract the intraband contribution in the present results, we employ the Landau-Peierls formula [51, 52] to calculate the orbital magnetic susceptibility that purely comes from the intraband contribution:

We note that the damping rate $\Gamma$ is neglected in this formula, hence we cannot compare the results for $\chi^{\mathrm{LP}}$ with $\chi$ presented in Fig. 3 in a quantitative manner. Nevertheless, we can extract the features of the $\mu$ dependence of $\chi$, as we will show later.

Figure 7 shows the $\chi$ and $\chi^{\rm LP}$ as functions of $\mu$ with $\mu\geq 0$. Note that we set $T=0.06$ in calculating $\chi^{\rm LP}$, so that we can compare $\chi^{\rm LP}$ with $\chi$ at the finite $\Gamma$ in an ad hoc manner. We see that the paramagnetic peaks of $\chi$ can be accounted for by $\chi^{\rm LP}$. As mentioned in Sec. III, the peak positions correspond to the van Hove singularities. Note that the emergence of the paramagnetic peaks at the van Hove singularities coincides with the various two-dimensional tight-binding models such as the single-orbital square lattice [6, 53] and the honeycomb lattice [6, 28, 29].

We also see that $\chi^{\rm LP}$ (dashed lines in Fig. 7) are finite for $\mu\sim 0$, which does not reproduce the behavior of $\chi$ near $\mu=0$ for all three types of the Dirac cones, meaning that the interband contribution plays an essential role, not for the type-I but also for types-II and III.