Orbital magnetization of interacting Dirac fermions in graphene

TL;DR

This paper develops a formalism to compute the orbital magnetization of interacting Dirac fermions in graphene, addressing divergence issues and incorporating interactions via mean-field theory, providing insights into magnetic properties of the system.

Contribution

A new formalism for calculating orbital magnetization of interacting Dirac fermions that respects particle-hole symmetry and accounts for interactions through mean-field theory.

Findings

Orbital magnetization differs significantly between interacting and non-interacting cases.

The formalism successfully overcomes divergence issues in magnetization calculations.

Interactions influence charge, spin orderings, and Landau level exchange interactions.

Abstract

We present a formalism to calculate the orbital magnetization of interacting Dirac fermions under a magnetic field. In this approach, the divergence difficulty is overcome with a special limit of the derivative of the thermodynamic potential with respect to the magnetic field. The formalism satisfies the particle-hole symmetry of the Dirac fermions system. We apply the formalism to the interacting Dirac fermions in graphene. The charge and spin orderings and the exchange interactions between all the Landau levels are taken into account by the mean-field theory. The results for the orbital magnetization of interacting Dirac fermions are compared with that of non-interacting cases.

Equations

$$M=-(\partial \Omega /\partial B{)}_{T,\mu }$$

$$H=-t\sum _{⟨ij⟩s}{c}_{is}^{†}{c}_{js}+U\sum _j\delta n_{j\uparrow }\delta n_{j\downarrow }+\frac{1}{2}\sum _{i\ne j}{v}_{ij}\delta n_i\delta n_j$$

$$v(r)=\frac{e^2}{r}[1-\exp (-q_{0}r)],$$

$$\Omega$$

$${\Sigma }_{l{l}^{\prime }}(k,i\omega )=\delta \Phi /\delta G_{l^{\prime }l}(k,i\omega ),$$

$$M=-\frac{{\Omega }_{N-1}(B+\Delta B)-{\Omega }_N(B)}{\Delta B}.$$

$${\Omega }_N=B{S}_N(B)$$

$$M$$

$$+(N+1/2)[y_N(B)+B{y}_N^{\prime }(B)/(N-1/2)]$$

$$H_v(p)=s_v{p}_x{\sigma }_1+p_y{\sigma }_2,$$

$$H_{vn}=\sqrt{2Bn}{\sigma }_1,$$

$$\Omega$$

$$S_N(B)$$

$$M$$

$$-c_0\sqrt{N}(N+1/2)(1+\frac{1}{2N-1})$$

$$S_N(B)=\frac{2}{\pi }[\sum _{n=1}^N(-\sqrt{2Bn}-\mu )+\sum _{n=1}^{N_F}(\sqrt{2Bn}-\mu )-\mu ]$$

$$S_N(B)$$

$$M=\frac{2}{\pi }(N_F+1/2)\{\mu -\sqrt{2B(N_F+1/2)}[1+O(N_F^{-2})]\}.$$

$$\Phi$$

$$F(k)$$

$${\Sigma }_{H,l{l}^{\prime }}$$

$${\Sigma }_{X,l{l}^{\prime }}(k)$$

$$V_{l{l}^{\prime }}^x$$

$$G(k,i\omega )=\sum _{\lambda }\frac{{\psi }_{\lambda }(k){\psi }_{\lambda }^{†}(k)}{i\omega +\mu -{ϵ}_{\lambda }(k)}$$

$$[\sqrt{2Bn}{\sigma }_1+\Sigma (k)]{\psi }_{\lambda }(k)={ϵ}_{\lambda }(k){\psi }_{\lambda }(k).$$

$${ϵ}_{\lambda }(k)$$

$${\psi }_{+}(k)$$

$${\psi }_{-}(k)$$

$${ϵ}_0(0Ks)$$

$${ϵ}_0(0K^{\prime }s)$$

$${\rho }_a$$

$$m_l$$

$$F_{ab}(k)=[f_{+}(k)-f_{-}(k)]\frac{{ϵ}_1(k)+{\Sigma }_1(k)}{2E(k)},$$

$${\Sigma }_1(k)$$

Full text

I intrduction

The study of the properties of interacting Dirac (or Weyl) fermions in (topological) semimetals under a magnetic field is a fundamental subject of the condensed matter physics Castro ; Goerbig . One of the physical themes is to investigate the orbital magnetization (OM) of the Dirac fermions (DFs) with Coulomb interactions. The OM of an electron system is usually defined as Shi

[TABLE]

where $\Omega=\Omega(T,\mu,B)$, as a function of the temperature $T$ and the chemical potential $\mu$ and the magnetic field $B$, is the thermal dynamic potential. Equation (1) is equivalent to a statistical average of the OM operator Hirst . However, for Dirac (or Weyl) fermions, Eq. (1) is ill defined because the occupation of the Landau levels in the lower band leads to divergence of $\Omega$ and thereby $M$. For noninteracting DFs in graphene, $\Omega$ can be evaluated with a special method Gusynin ; Hesse ; SGB ; Slizovskiy by which the field $B$ dependent part of $\Omega$ is separated out. The effects of finite-temperature occupations and the impurity broadening of the Landau levels on the OM of the noninteracting DFs have been studied SGB ; Slizovskiy ; Koshino . Nonetheless, for interacting DFs, it is not easy to separate the $B$-dependent part of $\Omega$ from that of the independent part. Study of the OM of Dirac fermions with Coulomb interactions is lacking. How to calculate the OM of interacting DFs is still an open question. In this paper, we are developing a general approach for solving this problem and use it to calculate the OM of interacting Dirac fermions in graphene.

II formalism

The electrons in graphene are moving on a honeycomb lattice of carbon atoms. The Hamiltonian of the electrons with a neutralizing background is

[TABLE]

where $c^{\dagger}_{is}$ ($c_{is}$) creates (annihilates) an electron of spin $s$ in site $i$, $\langle ij\rangle$ sums over the nearest-neighbor (NN) sites, $t\approx$ 3 eV is the NN hopping energy, $\delta n_{is}=n_{is}-n_{s}$ is the number deviation of electrons of spin $s$ at site $i$ from the average occupation $n_{s}$, and $U$ and $v_{ij}$ are the Coulomb interactions between electrons. In real space, $v_{ij}=v(r_{ij})$ with $r_{ij}$ the distance between sites $i$ and $j$ is given by

[TABLE]

where $q_{0}$ is a parameter taking into account the wavefunction spreading effect in the short-range interactions between electrons. Here we take $q_{0}=0.5/a_{0}$ with $a_{0}\approx 2.46$ Å  as the lattice constant of graphene. For carrier concentration close to the charge neutrality point (CNP), one usually adopts the simplified continuum model. With the continuum model and using the mean-field theory (MFT, or the self-consistent Hartree-Fock approximation), we have recently studied the Landau quantization of the interacting electrons taking into account the charge and spin orderings and the exchange interactions between all the levels Yan .

According to the many-particle theory Luttinger , the thermodynamical potential $\Omega$ per unit volume of an electron system under a magnetic field $B$ is given by

[TABLE]

where $\Phi$ is the ‘free energy’ functional of the Green’s function $G$, $\Sigma$ is the self-energy, $k$ is the state index, $\omega$ is the fermionic Matsubara frequency, and $\eta$ is an infinitesimal small positive quantity. For Dirac fermions in graphene, $G$ and $\Sigma$ are $2\times 2$ matrices in the space of sublattices $a$ and $b$, and $k$ stands for $(n,v,s)$ with $n,v,s$ respectively the indexes of the Landau level (LL) and valley and spin Yan . The self-energy matrix element $\Sigma_{ll^{\prime}}(k,i\omega)$ with $l(l^{\prime})=a$ or $b$ is related with $\Phi$ by

[TABLE]

which ensures the microscopic conservation law being satisfied Baym . The point here is, after the summation over the Matsubara frequency, $\Omega$ can be expressed as the sum over the LLs from $n=0$ to $\infty$. We will use the units in which $\hbar=e=c=a_{0}=1$, the energy unit $\epsilon_{0}=\hbar v_{0}/a_{0}=1$ (with $v_{0}$ the Fermi velocity of electrons in graphene), and the unit of magnetic field $B_{0}=\hbar c/ea_{0}^{2}=1$.

To get rid of the divergence difficulty, we consider a system in momentum space containing finite LLs as shown in Fig. 1. The cutoff momentum is given by $k_{c}=\sqrt{(2N+1)B}$ where $N$ is the highest Landau index at the field $B$. The thermodynamic potential of this finite system is then given by $\Omega_{N}(B)$ (suppressing the $T$ and $\mu$ dependence for brevity). The number $N$ changes with $B$ varying for fixed $k_{c}$. When the magnetic field $B$ varies from $B=k_{c}^{2}/(2N+1)$ to $B+\Delta B$ with $\Delta B=2B/(2N-1)$, the index of the highest LL changes to $N-1$. We then define the OM of the finite system as

[TABLE]

The ratio given by Eq. (4) with $k_{c}\to\infty$ can be considered as the special limit of the derivative in Eq. (1). For sufficient large cutoff $k_{c}$, this definition should give rise to the result of the entire system. For low carrier concentration close to the charge neutrality point (CNP), the cutoff can be taken as $k_{c}=1$. Here, we should remark that our finite system of $N$ LLs is part of the whole system of infinite LLs. It does not mean we can consider an isolated system of only $N$ LLs from the beginning. For the Dirac fermions, the cutoff for such an isolated system leads to unphysical results. The consideration of such an isolated Dirac system is equivalent to thinking only the top $N$ LLs being occupied with the rest lower LLs as empty in the lower band. This is apparently unphysical.

Now that $\Omega_{N}(B)$ contains $N$ terms, we define

[TABLE]

and suppose each term in $S_{N}(B)$ be an analytical function of $B$. Write $S_{N-1}(B+\Delta B)=S_{N}(B+\Delta B)-y_{N}(B+\Delta B)$ with $y_{N}$ the $N$th term in the sum $S_{N}$. Then, by expanding $S_{N}(B+\Delta B)$ to order $(\Delta B)^{2}$ and $y_{N}(B+\Delta B)$ to order $\Delta B$, the OM can be expressed as

[TABLE]

where the primes mean the derivatives with respect to $B$.

III OM of noninteracting Dirac fermions

As an example, here, we consider the free Dirac fermions in graphene at zero temperature. The Hamiltonian of a single Dirac fermion is

[TABLE]

where $s_{v}=1$ (-1) for particle in valley $v=K$ ($K^{\prime}$), the momentum $\vec{p}$ in each valley is measured from the Dirac point, and $\sigma^{\prime}s$ are the Pauli matrices operating in the sublattice ($a,b$) space. Under a magnetic field $B$ applied perpendicularly to the system plane, the states of the Dirac fermions are given by the Landau quantization. In the LL representation, the Hamiltonian (7) reads

[TABLE]

where $n$ is the LL index. The LLs are obtained as $\epsilon_{\lambda}(n)=\lambda\sqrt{2Bn}$ with $\lambda=\pm$ for $n\neq 0$, and $\epsilon_{0}=0$ for $n=0$. At CNP and $T=0$, the LLs in the lower band are fully occupied while the LLs in the upper band are completely empty. The thermodynamic potential reads (see Appendix)

[TABLE]

where the $k$ sum in the first line is understood over the Landau index $n$ and the valley $v$ and the spin $s$. The sum $S_{N}(B)$ is then obtained as ram

[TABLE]

with $c_{0}=2\sqrt{2B}/\pi$ and $\zeta(-1/2)=-0.207886225$. We then have $2BS_{N}^{\prime}(B)=-4B^{2}S_{N}^{\prime\prime}(B)=S_{N}(B)$ and $By_{N}^{\prime}(B)=y_{N}/2=-c_{0}\sqrt{N}/2$. $M$ is calculated as

[TABLE]

Due to the expansion of $S_{N}(B+\Delta B)$ to second order in $\Delta B$ and $y_{N}(B+\Delta B)$ to linear $\Delta B$, there are precise cancellations from $O(N^{3/2})$ to $O(1)$ in Eq. (6). The result given by Eq. (10) is consistent with the existing one SGB ; Slizovskiy ; Ghosal .

At finite doping with chemical potential $\mu>0$ and $T=0$, the sum $S_{N}(B)$ is given by

[TABLE]

where the sums in the square brackets are, respectively, from the lower and upper bands with $N_{F}$ the index of highest LL below the chemical potential, and the last term $-\mu$ comes from the zero LL. Suppose $\mu<<k_{c}$, we then have $N_{F}<<N$. In the limit $B\to 0$, because of $N_{F}>>1$, we obtain

[TABLE]

The OM is given by

[TABLE]

As $B\to 0$, $M$ oscillates rapidly between $-\mu/2\pi$ and $\mu/2\pi$ with period $\Delta B=2B^{2}/\mu^{2}$. This is the de Haas-van Alphen oscillation. The average of $M$ vanishes, which leads to the vanishing orbital magnetic susceptibility $\chi=0$ (defined as the derivative of $M$ with respect to $B$ at $B=0$) at finite doping. On the other hand, at the CNP, because $M\propto\sqrt{B}$ as given by Eq. (10), $\chi$ diverges at $B=0$. This is consistent with the existing result Koshino ; JWM ; Safran , $\chi=-(2/3\pi)\delta(\mu)$, which is obtained by the response of uniform Dirac fermions to the magnetic field without considering the Landau quantization Safran ; Fukuyama .

IV MFT for interacting Dirac fermions in graphene

As in our previous work, we use the MFT to deal with the interactions between the electrons Yan ; Cote . By the MFT, the ‘free energy’ functional $\Phi$ is approximated as shown in Fig. 2(a). The self-energy is then obtained as in Fig. 2(b), which is independent of the Matsubara frequency. In terms of $G$, $\Phi$ is given by

[TABLE]

with $\Sigma(k)=\Sigma_{H}+\Sigma_{X}(k)$ and

[TABLE]

where $\beta=1/k_{B}T$, and $\Sigma_{H,ll^{\prime}}$ has been written in terms of the charge $\rho_{l}$ and the spin $m_{l}$ order parameters with $v_{c}$ and $U$ the corresponding interaction parameters. The first sum in the first equal of Eq. (11) is due to the direct Coulomb interaction, while the second sum comes from the exchange interaction. Here, $v_{\mu\nu}(0)$ and $v^{x}_{\mu\nu}(k,k^{\prime})$ are the interaction elements in the LL representation; they are dependent on the magnetic field $B$ Yan . The appearance of the extra term -1/2 in addition to the diagonal distribution function $F_{ll}(k^{\prime})$ in Eq. (12) originates from the interaction form of the system of DFs with a neutralizing background given in terms of the density-density multiplication instead of the normal order of the fermion operators. Corresponding to this term, there is a shift $V^{x}/2$ from the self-energy as shown in Eq. (11); this shift is not drawn in the diagrams in Fig. 2. Because of this shift, the particle-hole symmetry of the system is reflected by the invariance under the transform $\mu\to-\mu$ with $\mu=0$ at the CNP Yan3 .

In the LL’s picture, the Green’s function is given by

[TABLE]

where $\psi_{\lambda}(k)$ is the $\lambda$th eigen-wave function with eigen-energy $\epsilon_{\lambda}(k)$. The LLs $\epsilon_{\lambda}(k)$ and the wavefunctions $\psi_{\lambda}(k)$ are determined by Yan

[TABLE]

Express the self-energy matrix as $\Sigma(k)=\Sigma_{0}(k)\sigma_{0}+\Sigma_{1}(k)\sigma_{1}+\Sigma_{3}(k)\sigma_{3}$. The energy levels for $n\neq 0$ are obtained as

[TABLE]

and the corresponding wavefunctions are

[TABLE]

where $R_{\pm}(k)=\sqrt{1\pm\Sigma_{3}(k)/E(k)}/\sqrt{2}$. For $n=0$, the eigenstates are given by

[TABLE]

in valleys $v=K$ and $v=K^{\prime}$, respectively. The charge and spin orders are calculated by

[TABLE]

where $s_{0}=\sqrt{3}/2$ is the area of the unit cell, $B/2\pi$ is the spatial degeneracy of the Landau state, $\psi_{l\lambda}(k)$ is the $l$th component of $\psi_{\lambda}(k)$ and $s_{l}$ = 1 (-1) for $l=a$ ($b$), $s=1$ (-1) for spin-up (down), and $f_{\lambda}(k)=f(\xi_{\lambda})=1/[\exp(\beta\xi_{\lambda})+1]$ with $\xi_{\lambda}(k)=\epsilon_{\lambda}(k)-\mu$ is the Fermi distribution function.

Here, we need to pay special attention to the equation for the self-energy element $\Sigma_{ab}(k)~{}[=\Sigma_{ba}(k)]$ or $\Sigma_{1}(k)$ given by Eq. (12). Using the wavefunctions given by Eq. (21), we have

[TABLE]

with $\epsilon_{1}(k)=\sqrt{2Bn}$. Note that $F_{ab}(k)$ goes to $-1/2$ in the limit $n\to\infty$. Equation for $\Sigma_{1}(k)$ can be written as

[TABLE]

with $v_{ab}^{xv}(n,n^{\prime})=v_{ab}^{x}(k,k^{\prime})$ and

[TABLE]

By so doing, the sum over $n^{\prime}$ in Eq. (LABEL:s1) converges fast. For $\Sigma_{X,ll}(k)$, Eq. (12) is the proper form since $F_{ll}(k)-1/2$ goes to zero in the limit $n\to\infty$ and therefore the sum over $n^{\prime}$ converges quickly. Usually, the self-energy given by Eq. (12) is evaluated with a cutoff $k_{c}=1$ Barlas ; Hwang ; Kusminskiy . By the similar treatment, we have solved Eq. (12) with cutoff $k_{c}=1$ for DFs in a magnetic field in our previous work Yan . This cutoff has little effect on the low energy levels close to zero. However, it influences substantially the high levels. In particular, the LLs at the cutoff are strongly modified. As indicated in Sec. II, we should solve the equations of the self-energy for the LLs in the whole range $0\leq n<\infty$. Therefore, the revision given by Eq. (LABEL:s1) is necessary. The big task now is to calculate $V^{x}_{1}(n)$.

To calculate $V^{x}_{1}(n)$, we first consider the case of $B=0$ and look for an approximation scheme from the result. By the transform $T(\phi_{v})=Diag[1,\exp(i\phi_{v})]$ with $\phi_{v}$ the angle of momentum $(s_{v}k_{x},k_{y})$, the effective mean-field Hamiltonian reads

[TABLE]

which is independent on the angle $\phi_{v}$. Here, $k$ is understood as the momentum. The self-energy element $\Sigma^{vs}_{ab}(k)$ reads

[TABLE]

where $V$ is the volume (area) of the two-dimensional system, $\theta$ is the angle between $\vec{k}$ and $\vec{k}^{\prime}$, and $F^{vs}_{ab}(k)$ has the same form as given by Eq. (30) provided the Landau energy $\epsilon_{1}(k)$ is replaced with the energy $k$. We can revise Eq. (34) to get a similar form as Eq. (LABEL:s1) and obtain the corresponding $V_{1}(k)$ as

[TABLE]

where $f(k,q)=[(k-q)K(\alpha)+(k+q)E(\alpha)]/\pi k$ with $K(\alpha)$ and $E(\alpha)$ the elliptic integrals and $\alpha=2\sqrt{kq}/(k+q)$. Now, for the quantized interaction $V^{x}_{1}(n)$, a reasonable approximation is to replace the continuous momentum $k$ with the quantized one $k_{n}=\sqrt{2Bn}$ in $V_{1}(k)$,

[TABLE]

V thermodynamic potential

Using the result (see Appendix)

[TABLE]

we obtain $\Omega(B)$ under the MFT as

[TABLE]

From Eq. (37), we may get $\Omega_{N}(B)$. However, to maintain the particle-hole symmetry in $\Omega_{N}(B)$, we must revise the form.

We need to write the equations for the self-energy $\Sigma_{0,3}(k)$ more clearly

[TABLE]

where $\rho=\rho_{a},m_{0,3}=(m_{a}\pm m_{b})/2,g_{\lambda}(k)=f_{\lambda}(k)-1/2$. Note that $v^{K}(n,n^{\prime})=\sigma_{1}v^{K^{\prime}}(n,n^{\prime})\sigma_{1}$, we then obtain

[TABLE]

where we have used the relation

[TABLE]

which is independent on $n$. Using Eq. (39), we rewrite Eq. (37) in the form

[TABLE]

where $F_{0,1,3}$ are distribution functions defined as

[TABLE]

$f_{0}=f_{0}(0vs)$ is the Fermi distribution function of level $n=0$, and ${\sum_{k}}^{\prime}$ means $n\neq 0$.

Under the transform $\mu\to-\mu$, the self-energy components change as $\Sigma_{0}(nvs)=-\Sigma_{0}(nvs)$ and $\Sigma_{1}(nvs)=\Sigma_{1}(nvs)$ and $\Sigma_{3}(nvs)=-\Sigma_{3}(n\bar{v}s)$ or $\Sigma_{1,3}(nvs)=\Sigma_{1,3}(n\bar{v}s)$ (with $\bar{v}$ means $\bar{K}=K^{\prime}$ and $\bar{K}^{\prime}=K$). Note that the constant terms $\mu-V^{x}/4$ and $\mu/2-V^{x}/8$ in Eq. (41) will disappear in the final formula for $M$ because a cancellation between the terms $-S_{N}(B)$ and $(N+1/2)y_{N}(B)$ as indicated by Eq. (6). We can then conclude that $M$ is symmetric under the particle-hole transform. The function $S_{N}(B)$ can now be extracted from Eq. (41).

For calculating $S_{N}^{\prime}(B)$ and $S_{N}^{\prime\prime}(B)$, we need to derive the equations of the self-energy elements with respect to $B$ and solve them. The derivation is elementary but tedious. For brevity of the paper, we will not express these equations here.

VI orbital magnetization

We have numerically solved the equations for the self-energy $\Sigma(k),\partial\Sigma(k)/\partial B$, and $\partial^{2}\Sigma(k)/\partial B^{2}$. In the present calculation, the on-site interaction is set as $U/\epsilon_{0}=2$. The coupling constant of the interaction is $e^{2}/a_{0}\epsilon_{0}=2.2$. With the results for the self-energy and its derivatives, we calculate the OM at CNP and at finite carrier concentration.

Shown in Fig. 3 are the numerical results for the interacting and free DFs at CNP and at $T=0$. It is seen that the magnitude of the OM of interacting DFs (blue solid circles) is smaller than that of the free DFs (red circles). At $T=0$, there exists antiferromagnetic spin ordering in the interacting DFs catalyzed by the magnetic field as investigated in many works Yan ; Gusynin1 ; Herbut ; Kharitonov ; Lado ; Roy ; Khveshchenko ; Alicea ; Gusynin2 ; Jung ; Lukose . This spin ordering results in the splitting of the zero Landau levels. In the low energy zero-LL states, the spin-up and down electrons move in the sublattices $a$ and $b$, respectively. The spin ordering also modifies the electron distributions in the two sublattices at other LLs. Overall, in the presence of the spin ordering, the electrons cannot move freely in the whole lattice. Since the antiferromagnetic spin ordering acts as the obstacle for the orbital circumnutation, the OM is therefore weakened.

The black line in Fig. 3 represents the analytical formula Eq. (10) for $N\to\infty$. In the numerical calculation, the cutoff $N$ is finite given by $N=k_{c}^{2}/2B-1/2$ with $k_{c}=1$ (and $B$ in units of $B_{0}=1.1\times 10^{4}$ T). At small $B$ close to zero, since $N$ is sufficiently large, the numerical result (red circles) for the free DFs is in very good agreement with the analytical formula. The difference between them increases with increasing $B$. For $B\sim 8$ T, the numerical result seems still good.

In Fig. 4, we present the results at $T/\epsilon_{0}=0.01$ and at CNP. Since $T$ is high, there are many LLs within the temperature range. As a result, the OM of free DFs varies linearly with $B$ consistent with the existing result Slizovskiy . While for the interacting DFs, the OM is not linear in $B$ and its magnitude is larger than that of free DFs. At this high temperature, the spin ordering vanishes but the LLs of the DFs are strongly changed by the interactions through the self-energy $\Sigma_{1}$. $\Sigma_{1}$ gives rise to an enhancement of the velocity Yan2 ; Sarma ; Menezes , leading to fast orbital circumnutations. The nonlinear behavior of $M$ with $B$ implies the renormalized velocity varies with momentum. Because of the vanishing of spin ordering and the enhancement of the velocity, the OM of the interacting DFs is stronger than that of the free DFs.

Figure 5 exhibits the result for the OM of the interacting DFs at finite carrier concentration with $\mu/\epsilon_{0}=0.02$ and at $T/\epsilon_{0}=0.001$. The chemical potential for the free DFs is set as $\mu_{0}/\epsilon_{0}=0.00167$ so that the first LL in the upper band for both interacting DFs and the free DFs has almost the same position $B$. At the finite carrier concentration and temperature, the charge and spin orderings disappear and all the LLs are degenerated with degeneracy 4. The index of the first LL in the upper band is $n=1$. With the field $B$ varying, when the LLs pass cross the Fermi level, the OM shows the de Haas-van Alphen oscillations. The LL of $n=1$ is at about $B\approx 0.97$ T. Above this field, there are no LLs below the Fermi level in the upper band and the OM decreases monotonically with $B$.

For finite doping at very small $B$ and $T=0$, there are rapid de Haas-van Alphen oscillations similarly as that indicated in Sec. III for noninteracting DFs. For interacting DFs, however, the average of the oscillations should not be vanishing at small $B$. According to the perturbation theory, the system shows orbital paramagnetism at very small $B$ Principi ; with the first order perturbation calculation for Thomas-Fermi screened Coulomb interactions, it has been shown that the orbital magnetic susceptibility $\chi$ is positive for DFs in doped graphene. Therefore, the average $M$ should increase from $M=0$ with increasing the field $B$. At finite $T$, the oscillations are smeared by temperature. The average $M$ should be weakened by the thermal fluctuations. (We did not perform the calculation at very small $B$ because for which the cutoff number $N$ is so large that the accuracy requirement for the numerical calculation exceeds the ability of our computer.)

VII Remark

In the present approach, the OM is calculated by expanding the sum $S_{N}(B+\Delta B)$ [and the term $y_{N}(B+\Delta B)$] to second (first) order in $\Delta B$ as shown in Eq. (6). The formalism works only for the system of eigen-energy being linear in momentum $k$. For a Dirac or Weyl system of $\epsilon_{\lambda}(k)\to\lambda k^{\nu}$ as $k\to\infty$, the sum $S_{N}$ is order $N^{\nu/2+1}$. We need to expand $S_{N}(B+\Delta B)$ [$y_{N}(B+\Delta B)$] to $m$th ($m-1$th) order in $\Delta B$ with $m=[\nu/2]+2$. Here $[\nu/2]$ means the integer part of the number $\nu/2$. For example, for an $L$-layered graphene, since it has $[L/2]$ bilayer bands and $L$ mod 2 monolayer bands Koshino1 , we need to expand $S_{N}(B+\Delta B)$ to $(\Delta B)^{3}$ and $y_{N}(B+\Delta B)$ to $(\Delta B)^{2}$. By so doing, the unphysical part will be eliminated due to the precise cancellations between these expanded terms.

Though the system of infinitive LLs is considered, the contribution to the total OM comes mostly from the LLs below $k_{c}=1$ as reflected by the result for free DFs shown in Fig. 3. In graphene, the Dirac cone approximation to the energy bands of electrons is valid within the circle of radius $k_{c}=1$ in the momentum space. Therefore, the present result for the OM is a fairly good measure of that of electrons in graphene. However, as already stressed, we cannot isolate the LLs below $k_{c}$ from the entire system. The reason is that the high LLs (especially the LLs close to the cutoff) are strongly modified by the isolation. We have performed the numerical calculation for the isolated system. The consequence of the isolated system is that the magnitude of the OM is several orders larger than the result presented here; it becomes bigger and bigger as $B\to 0$ even not vanishing at $B=0$.

VIII conclusion

We have developed the approach for calculating the orbital magnetization of Dirac fermions. The main points in the formalism are: (1) To overcome the divergence difficulty due to the occupation in the lower band, the orbital magnetization is defined as the special limit for the derivative of the thermodynamic potential with respect to the magnetic field. (2) The equations for the self-energy and its derivatives with respect to the magnetic field need to be solved. (3) The particle-hole symmetry should be ensured in the partly sum of the thermodynamic potential. (4) The system with finite LLs is part of the entire system but not isolated from the rest of the entire system.

With the formalism, we have calculated the OM for interacting DFs in graphene and compared the results with that of the free DFs. At very low carrier concentration close to CNP, when the antiferromagnetic spin ordering catalyzed by the magnetic field exists, the OM is weakened. Without the spin and charge orderings, the OM is enhanced due to the velocity renormalization by interactions. At low temperature and finite carrier concentration, the de Haas-van Alphen oscillation appears in the OM as a function of magnetic field.

The present approach may be extended to study the OM of Weyl fermions in the topological semimetals as well.

Acknowledgements.

This work was supported by the National Basic Research 973 Program of China under Grant No. 2016YFA0202300 and the Robert A. Welch Foundation under Grant No. E-1146.

Appendix A

1. Matsubara-frequency sum. For calculating $\Omega$, we need the following sum,

[TABLE]

which is usually performed with the loop integral in the complex $z$ plane,

[TABLE]

where $f(\omega)=1/[\exp(\beta\omega)+1]$ and $\omega^{+}=\omega+0^{+}$. At $T=0$, it reduces to

[TABLE]

2. Exchange interaction. The exchange interaction should contain the screening effect due to the electron density fluctuations. For qualitatively reflecting the screening, we adopt the Thomas-Fermi form for exchange interaction given as

[TABLE]

where $v(q)$ is the Fourier transform of the interaction $v(r)$. In the continuum model, $v(q)$ is given by

[TABLE]

By using the long-wavelength limit of the density-density response function of free Dirac fermions with chemical potential $\mu$ at temperature $T$ Yan1 , the Thomas-Fermi wave number $q_{TF}$ is obtained as

[TABLE]