# Tunable orbital susceptibility in α-T 3 tight-binding models

F Piéchon, J-N Fuchs, A Raoux, G Montambaux

► To cite this version:

F Piéchon, J-N Fuchs, A Raoux, G Montambaux. Tunable orbital susceptibility in α-T 3 tight-binding models. Journal of Physics: Conference Series, 2015, 603, pp.012001 10.1088/1742-6596/603/1/012001. [hal-01265564](https://hal-fr.sorbonne-universite.fr/hal-01265564)

HAL Id: hal-01265564

https://hal.sorbonne-universite.fr/hal-01265564

Submitted on 1 Feb 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Tunable orbital susceptibility in α-$\mathcal{T}_3$ tight-binding models

F Piéchon<fnref n="1" />, J-N Fuchs<fnref n="1" />,<fnref n="2" />, A Raoux<fnref n="1" />,<fnref n="3" /> and G Montambaux<fnref n="1" />

E-mail: piechon@lps.u-psud.fr

**Abstract.** We study the importance of interband effects on the orbital susceptibility of three bands α-$\mathcal{T}_3$ tight-binding models. The particularity of these models is that the coupling between the three energy bands (which is encoded in the wavefunctions properties) can be tuned (by a parameter α) without any modification of the energy spectrum. Using the gauge-invariant perturbative formalism that we have recently developed [1], we obtain a generic formula of the orbital susceptibility of α-$\mathcal{T}_3$ tight-binding models. Considering then three characteristic examples that exhibit either Dirac, semi-Dirac or quadratic band touching, we show that by varying the parameter α and thus the wavefunctions interband couplings, it is possible to drive a transition from a diamagnetic to a paramagnetic peak of the orbital susceptibility at the band touching. In the presence of a gap separating the dispersive bands, we show that the susceptibility inside the gap exhibits a similar dia to paramagnetic transition.

## 1. Introduction

The orbital magnetic susceptibility of free electrons was computed long ago by Landau [2] and was found to be diamagnetic. Subsequently, Peierls extended Landau's result to the case of a single band tight-binding model. He found a formula for the orbital susceptibility that only depends on the zero-field band energy spectrum. This result already showed that the band structure can have a strong influence on the magnetic response of crystals [3]. For example, near a saddle point of the dispersion relation, the Peierls orbital susceptibility becomes paramagnetic [4]. However, one important effect was left out, namely the coupling between the bands in the case of several bands. A striking example is the possibility of having a finite orbital susceptibility inside the gap of a band insulator at zero temperature. This was understood as an inter-band effect by Fukuyama and Kubo on the example of bismuth [5]. Fukuyama also provided a compact and quite general linear-response formula that includes interband effects [6]. However, despite its many successes, this formula does not work for non-separable tight-binding models. The aim of this paper is to present orbital susceptibility results obtained from an exact linear response formula that we recently derived for tight-binding models [1]. In order to show the importance of inter-band (or band coupling) effects, we consider here a family of tunable three band tight-binding models constructed on a $\mathcal{T}_3$ lattice (or dice lattice) and that we call α-$\mathcal{T}_3$ [7]. The dice model was first introduced [8] as a simple 2-dimensional model that exhibits localized

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹ Laboratoire de Physique des Solides, CNRS UMR 8502, Université Paris-Sud, F-91405 Orsay Cedex, France</footnote>

<footnote n="2">² Laboratoire de Physique Théorique de la Matière Condensée, CNRS UMR 7600, Univ. Pierre et Marie Curie 4, place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="3">³ Département de Physique, École Normale Supérieure, 24 rue Lhomond, 75005 Paris, France</footnote>
and extended states as the same time. It also presents some peculiar properties under strong magnetic field [9]. These $\alpha$-$\mathcal{T}_3$ models depend on a real parameter $\alpha$ and have the important property that their zero-field energy spectrum –which is essentially that of graphene with an additional flat band– is independent of $\alpha$. However, the parameter $\alpha$ has a strong influence on the zero-field eigenstates and therefore on the Berry curvature. This influence is revealed by the orbital susceptibility that changes sign as a function of $\alpha$.

The paper is organized as follows. In §2 we present the three bands $\alpha$-$\mathcal{T}_3$ tight-binding models and characterize the key properties of their energy spectrum and wavefunctions. In §3, using the gauge-invariant perturbative formalism that we have recently developed [1], we provide a generic formula for the orbital susceptibility of $\alpha$-$\mathcal{T}_3$ tight-binding models. In §4 we apply this susceptibility formula to three characteristic examples of $\alpha$-$\mathcal{T}_3$ tight-binding models that exhibit respectively Dirac, semi-Dirac and quadratic band touching at low-energy. In §5 we summarize the main results of this study.

## 2. Tight-binding models on the $\mathcal{T}_3$ lattice

Starting from the honeycomb lattice with two sites $(A, B)$ per unit cell, the $\mathcal{T}_3$ or dice lattice is obtained by connecting additional $(C)$ sites at the center of each hexagon to the $B$ sites (see Fig. 1). The dice lattice is thus a triangular Bravais lattice with three sites $(A, B, C)$ per unit cell. We consider tight-binding models that consist of spinless electrons hopping on this lattice. In its simplest form, we allow for a constant onsite potential term $+\Delta$ on sites $A, C$ and $-\Delta$ on sites $B$ and an isotropic nearest-neighbors hopping with amplitude $c_\alpha t$ from $A$ to $B$ and $s_\alpha t$ from $C$ to $B$ with $c_\alpha = \frac{1}{\sqrt{1+\alpha^2}}$, $s_\alpha = \frac{\alpha}{\sqrt{1+\alpha^2}}$ such that $c_\alpha^2 + s_\alpha^2 = 1$. The real space representation of the corresponding Hamiltonian follows as

$$h = \sum_{\mathbf{r}_B} \left[ c_{\alpha} t (\delta_{\mathbf{r}_B - \mathbf{r}_A, \mathbf{r}_A} + \delta_{\mathbf{r}_B - \mathbf{r}_2, \mathbf{r}_A} + \delta_{\mathbf{r}_B - \mathbf{r}_3, \mathbf{r}_A}) |\mathbf{r}_B\rangle\langle\mathbf{r}_A| + s_{\alpha} t (\delta_{\mathbf{r}_B + \mathbf{r}_1, \mathbf{r}_C} + \delta_{\mathbf{r}_B + \mathbf{r}_2, \mathbf{r}_C} + \delta_{\mathbf{r}_B + \mathbf{r}_3, \mathbf{r}_C}) |\mathbf{r}_B\rangle\langle\mathbf{r}_C| \right] + \text{h.c.}, \quad (1)$$

where $\delta_1, \delta_2, \delta_3$ are the three vectors connecting nearest-neighbors sites. We assume that the localized orbital basis is orthogonal $(\langle\mathbf{r}_{j'}|\mathbf{r}_j\rangle = \delta_{\mathbf{r}_{j'}, \mathbf{r}_j})$ such that the position operator is purely diagonal $(\mathbf{r} = \sum_{\mathbf{r}_{j}=A,B,C} \mathbf{r}_j|\mathbf{r}_j\rangle\langle\mathbf{r}_j|)$. Introducing the Bloch states basis $|\mathbf{k}_{j=A,B,C}\rangle = \sum_{\mathbf{r}_j} e^{i\mathbf{k}\cdot\mathbf{r}_j}|\mathbf{r}_j\rangle$, for each wavevector $\mathbf{k} = (k_x, k_y)$, the Bloch Hamiltonian matrix associated to this $\alpha$-$\mathcal{T}_3$ tight-binding model reads

$$h_{\mathbf{k}} = \begin{pmatrix} \Delta & c_{\alpha} f_{\mathbf{k}} & 0 \\ c_{\alpha} f_{\mathbf{k}}^{*} & -\Delta & s_{\alpha} f_{\mathbf{k}} \\ 0 & s_{\alpha} f_{\mathbf{k}}^{*} & \Delta \end{pmatrix}, \quad (2)$$

where $f_{\mathbf{k}} = |f_{\mathbf{k}}|e^{-i\theta_{\mathbf{k}}} = t(e^{-i\mathbf{k}\cdot\delta_1} + e^{-i\mathbf{k}\cdot\delta_2} + e^{-i\mathbf{k}\cdot\delta_3})$. In the following, we will also consider generalized versions of this model that essentially consist in a modified $f_{\mathbf{k}}$. This class of models interpolates between the honeycomb ($\alpha = 0$) and the isotropic dice lattice ($\alpha = 1$).

The remarkable and interesting property of this class of models is that the energy band spectrum does not depend on $\alpha$ (see Fig. 1) whereas the eigenfunctions do. More quantitatively, the energy spectrum consists of two dispersive bands $\epsilon_{\pm, \mathbf{k}} = \pm\sqrt{\Delta^2 + |f_{\mathbf{k}}|^2}$ and a flat band at energy $\epsilon_0 = \Delta$. For $\Delta > 0$, the corresponding Bloch eigenfunctions $|s_{\alpha, \mathbf{k}}\rangle$ ($s = \pm, 0$) read

$$|+_{\alpha, \mathbf{k}}\rangle = \frac{\sqrt{1+c_{\mathbf{k}}}}{\sqrt{2}} \begin{pmatrix} c_{\alpha} e^{-i\theta_{\mathbf{k}}} \\ \frac{s_{\mathbf{k}}}{1+c_{\mathbf{k}}} \\ s_{\alpha} e^{i\theta_{\mathbf{k}}} \end{pmatrix}; \quad |-_{\alpha, \mathbf{k}}\rangle = \frac{s_{\mathbf{k}}}{\sqrt{2(1+c_{\mathbf{k}})}} \begin{pmatrix} c_{\alpha} e^{-i\theta_{\mathbf{k}}} \\ -\frac{1+c_{\mathbf{k}}}{s_{\mathbf{k}}} \\ s_{\alpha} e^{i\theta_{\mathbf{k}}} \end{pmatrix}; \quad |0_{\alpha, \mathbf{k}}\rangle = \begin{pmatrix} s_{\alpha} e^{-i\theta_{\mathbf{k}}} \\ 0 \\ -c_{\alpha} e^{i\theta_{\mathbf{k}}} \end{pmatrix} \quad (3)$$
<figure><img src="image_3.png" /><figcaption>Figure 1. (Color online). α-T$_{3}$ model: (a) The dice or T$_{3}$ lattice is constituted by two interpenetrating honeycomb lattices. Thick links, nearest-neighbor hoppings c$_{αt}$ from sites A to sites B (c$_{α}$ = 1/√1+α $^{2}$). Thin links, nearest-neighbor hoppings s$_{αt}$ from sites C to sites B (s$_{α}$ = α/√1+α $^{2}$). By varying α , the model interpolates between honeycomb ( α = 0) and dice ( α = 1). Potential term +Δ on sites A, C and − Δ on sites B . (b) Energy bands dispersion in k space, for Δ = 0: two dispersive bands and one flat band. This spectrum does not depend on α .</figcaption></figure>

where c$_{k}$ = Δ $_{√Δ}$2 + | f$_{k}$ | 2 and s$_{k}$ = | f$_{k}$ | $_{√Δ}$2 + | f$_{k}$ | 2 with c 2 k + s 2 k = 1. The Berry connection A s k = 〈 s$_{α,k}$ | i ∇$_{k}$ | s$_{α,k}$ 〉 of each band is then readily obtained as

$$
\mathcal{A}_{\mathbf{k}}^{0} = -\frac{1 - \alpha^2}{1 + \alpha^2} \nabla_{\mathbf{k}} \theta; \quad \mathcal{A}_{\mathbf{k}}^{+} = -\frac{1 + c_{\mathbf{k}}}{2} \mathcal{A}_{\mathbf{k}}^{0}; \quad \mathcal{A}_{\mathbf{k}}^{-} = -\frac{s_{\mathbf{k}}^2}{2(1 + c_{\mathbf{k}})} \mathcal{A}_{\mathbf{k}}^{0}, \qquad (4)
$$

such that quite generally A 0 $_{α,k}$+ A + α,k + A − α,k = 0. From this Berry connection it is possible to obtain, for each band, the so-called Berry curvature Ω s k = ∇$_{k}$ × A s $_{k}$. This Berry curvature usually presents some strongly peaked structure near specific k ∗ points where interband coupling is strong, even if there is a gap separating the bands. Moreover, by considering a k-space closed orbit C k ∗ around such a k ∗ point, it is possible to calculate a Berry phase that measures the winding of the phase θ$_{k}$ around such points k $^{∗}$. More precisely, considering closed orbits C k ∗ ( ε ) that corresponds to a constant energy ε = √ Δ 2 + | f$_{k}$ | $^{2}$, we deduce that for each band the Berry phase Φ s k ∗ = ∮ C k ∗ d k · A s k accumulated along such an orbit is given respectively by Φ 0 k ∗ = − Φ$_{B}$ W$_{C}$$_{k}$ ∗ with Φ ± k ∗ = − 1 $_{2}$(1 ± Δ $_{ϵ}$)Φ 0 k ∗ where Φ$_{B}$ = π 1 − α 2 1+ α 2 and W$_{C}$$_{k}$ ∗ = ∮ C k ∗ d k · ∇$_{k}$θ/ 2 π is the winding number of the orbit such that W$_{C}$ = 0 for a closed orbit around a regular k $^{∗}$-point and W$_{C}$$_{k}$ ∗ = ± 1 for an orbit encircling a Dirac point k ∗ [7, 10].

From this perspective, the main interest of the α -T$_{3}$ models described by Eq.(2) is that all these Berry quantities (connections, curvature and phase) are proportional to Φ$_{B}$ = π 1 − α 2 1+ α 2 (see Eq.(4)). As a consequence, the interband effects that are encoded by these Berry quantities are reduced upon increasing the value of α from 0 and totally vanish for α = 1. Anticipating on what follows, we point out that our results for the orbital susceptibility of α -T$_{3}$ models provide clear evidence that all interband effects are not only encoded by the Berry quantities.
### 3. Orbital susceptibility formula

For time reversal systems which we consider in the following, the orbital magnetization vanishes and the orbital susceptibility $\chi_{\text{orb}}(\mu, T)$ is the first measure of the sensitivity of the energy spectrum when it is placed in a uniform static perpendicular magnetic field $B = B\mathbf{u}_z$:

$$ \chi_{\text{orb}}(\mu, T) = -\frac{\mu_0}{S} \left. \frac{\partial^2 \Omega}{\partial B^2} \right|_{B=0} \qquad (5) $$

where $S$ is the area of the system, $T$ is the temperature, $\mu$ is the chemical potential, $\mu_0 = 4\pi \cdot 10^{-7} \text{S.I.}$ and $\Omega(T, \mu, B)$ is the grand canonical potential of the non-interacting electrons gas

$$ \Omega(T, \mu, B) = -T \int_{-\infty}^{+\infty} \ln \left( 1 + e^{-(\omega - \mu)/T} \right) \rho(\omega, B) \, d\omega, \qquad (6) $$

with $\rho(\omega, B)$ the field dependent density of states (Dos)

$$ \rho(\omega, B) = -\frac{1}{\pi} \Im \operatorname{Tr} G(\omega + i0^+, B) \qquad (7) $$

where $G(\omega, B)$ is the field dependent retarded Green's function (we use $k_B = 1$ and $\hbar = 1$).

In order to obtain the susceptibility, it is thus necessary to calculate the density of states $\rho(\omega, B)$ or the Green's function $G(\omega, B) = (\omega I - h(B))^{-1}$ associated to the field dependent Hamiltonian $h(B)$. For the tight-binding model (1), $h(B)$ is obtained by multiplying the real space hopping amplitudes by a gauge and position dependent Peierls phase:

$$ h(B) = \sum_{\mathbf{r}_B} \left[ c_{\alpha} t (\delta_{\mathbf{r}_B+\mathbf{r}_A, \mathbf{r}_A} + \delta_{\mathbf{r}_B+\mathbf{r}_A, \mathbf{r}_A} + \delta_{\mathbf{r}_B+\mathbf{r}_A, \mathbf{r}_A}) e^{-i\varphi_{\mathbf{r}_A, \mathbf{r}_B}} |\mathbf{r}_B\rangle\langle\mathbf{r}_A| + s_{\alpha} t (\delta_{\mathbf{r}_B-\mathbf{r}_A, \mathbf{r}_C} + \delta_{\mathbf{r}_B-\mathbf{r}_A, \mathbf{r}_C} + \delta_{\mathbf{r}_B-\mathbf{r}_A, \mathbf{r}_C}) e^{-i\varphi_{\mathbf{r}_C, \mathbf{r}_B}} |\mathbf{r}_B\rangle\langle\mathbf{r}_C| \right] + \text{h.c.}, \qquad (8) $$

where $\varphi_{\mathbf{r},\mathbf{r}'} = \frac{e}{\hbar} \int_{\mathbf{r}}^{\mathbf{r}'} \mathbf{A} \cdot d\mathbf{l}$ with $\mathbf{A}(\mathbf{r})$ the gauge dependent vector potential associated to the uniform magnetic field. Starting from the Hamiltonian $h(B)$ there are essentially two approaches to calculate $\rho(\omega, B)$.

The first approach consists in computing the exact magnetic field dependent energy spectrum $\epsilon_n(B)$ associated to $h(B)$; this leads to the so-called Hofstadter butterfly spectrum. Using such an approach, we have recently studied the orbital susceptibility of the $\alpha$-$\mathcal{T}_3$ model (2) for the case $\Delta = 0$ [7]. In particular we have shown that the orbital susceptibility $\chi_{\text{orb}}(\mu, T)$ strongly varies with the parameter $\alpha$. More precisely, at the Dirac point, $\chi_{\text{orb}}(\mu = 0, T)$ exhibits a continuous transition from a diamagnetic peak for $\alpha = 0$ (honeycomb-graphene) to a paramagnetic peak for $\alpha = 1$ (dice) (see Fig. 2). Away from the Dirac point, an opposite transition from paramagnetism to diamagnetism takes place such that a sum rule $\int d\mu \chi_{\text{orb}}(\mu, T) = 0$ is preserved for all $\alpha$ [7]. Despite these interesting results, such an approach suffers from being essentially numerical and thus it does not allow to fully understand how the dependence on the parameter $\alpha$ enters in the susceptibility.

In order to better understand the role of the parameter $\alpha$, we have developed a second approach which consists in calculating the susceptibility $\chi_{\text{orb}}(\mu, T)$ using our recently established gauge-invariant perturbative response formula [1]:

$$ \chi_{\text{orb}} = -\frac{\mu_0 e^2}{12\hbar^2} \frac{\Im m}{\pi S} \int_{-\infty}^{+\infty} d\omega \, n_F(\omega) \operatorname{Tr} \{gh^{xx}gh^{yy} - gh^{xy}gh^{xy} - 4(gh^xgh^xgh^ygh^y - gh^xgh^ygh^xgh^y)\}, \qquad (9) $$

where $n_F(\omega) = 1/(e^{\omega/\mu} + 1)$ is the Fermi function, $g(\omega) = (\omega I - h)^{-1}$ is the retarded Green's function associated to the zero-field Hamiltonian and $h^x = [x, h], h^{xx} = [x, [x, h]]$ are the single
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (Color online). Orbital susceptibility $\chi(\mu, T)$ obtained from the numerically computed Hofstadter spectrum [7]. $\chi$ (in units of the Landau band edge value $|\chi_L| = (\mu_0/16\pi)(e^2ta^2/\hbar)$) as a function of the chemical potential $\mu$ (in units of $t$) in the whole band for various $\alpha$ as indicated and for a temperature $T = 0.02t$. At $\mu = 0$ there is a transition from a diamagnetic peak for $\alpha = 0$ (red: graphene) to a paramagnetic peak for $\alpha = 1$ (blue: dice). Because of a sum rule, the orbital response at zero doping is systematically compensated by an opposite response at finite doping.</figcaption></figure>

and double commutators of the position operator with the zero-field Hamiltonian. As discussed in [1], the expression (9) is equivalent to a recent formula derived for graphene [11] but it differs from the well-known Fukuyama formula [6]. We stress however that when $h^{xy} \neq 0$, only formula (9) fully agrees with the susceptibility results obtained from direct numerical computation of the corresponding Hofstadter butterfly spectrum. Moreover, we remind that formula (9) is valid not only for Bloch electrons in infinite crystals but it also applies to disordered and finite systems as well.

In the present paper, we restrict to multiband Bloch electrons in infinite crystals. In that situation, the trace operator is rewritten $\text{Tr}(\bullet) = \sum_k \text{tr}(\bullet) = S \int \frac{d^2k}{4\pi^2} \text{tr}(\bullet)$ where the integration is performed over the first Brillouin zone (BZ) and $\text{tr}(\bullet)$ is the partial trace operator on the band index. Accordingly, the susceptibility formula follows as

$$ \chi_{\text{orb}} = -\frac{\mu_0 e^2 \Im m}{12 \pi} \int_{-\infty}^{+\infty} d\omega \int \frac{d^2 k}{4\pi^2} n_F(\omega) [U_k(\omega) - 4V_k(\omega)], \quad (10) $$

where

$$ \begin{aligned} U_k(\omega) &= \text{tr}\{(gh^{xx}gh^{yy} - gh^{xy}gh^{xy})_k\}, \\ V_k(\omega) &= \text{tr}\{(gh^xgh^xgh^ygh^y - gh^xgh^ygh^xgh^y)_k\}, \end{aligned} \quad (11) $$

with $g_k(\omega) = (\omega I - h_k)^{-1}$ the Green's function matrix associated to the zero-field Bloch Hamiltonian matrix and $h_k^j = \frac{\partial h_k}{\partial k_j}$, $h_k^{ij} = \frac{\partial^2 h_k}{\partial k_i \partial k_j}$ with $(i, j) \in (x, y)$. As detailed in the appendix, for the class of models described by Eq.(2), by defining the two components vectors $(\vec{x} = (x_1, x_2))$

$$ \vec{f}_k = (\Re f_k, \Im m f_k), \quad \vec{f}_k^i = \left( \frac{\partial \Re f_k}{\partial k_i}, \frac{\partial \Im m f_k}{\partial k_i} \right), \quad \vec{f}_k^{ij} = \left( \frac{\partial^2 \Re f_k}{\partial k_i \partial k_j}, \frac{\partial^2 \Im m f_k}{\partial k_i \partial k_j} \right), \quad (12) $$

and the quantities

$$ u_k^j = \frac{\vec{f}_k^j \cdot \vec{f}_k}{|f_k|}, \quad v_k^j = \frac{(\vec{f}_k^j \times \vec{f}_k)}{|f_k|}, \quad u_k^{ij} = \frac{\vec{f}_k^{ij} \cdot \vec{f}_k}{|f_k|}, \quad v_k^{ij} = \frac{(\vec{f}_k^{ij} \times \vec{f}_k)}{|f_k|}, \quad (13) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (color online) Honeycomb lattice model that determines the form $f_k$: $t$ (full lines) and $t_3$ (dashed lines) respectively first and third nearest-neighbor hopping amplitudes. The two sites in the unit cell are called A and B and are shown as blue and red dots. In a uniaxially compressed honeycomb lattice, there are two different values for the nearest-neighbor hopping amplitudes: $t$ for thin (non-vertical) lines and $t' \ge t$ for thick (vertical) lines.</figcaption></figure>

it is possible to obtain the following compact expressions for $U_k(\omega)$ and $V_k(\omega)$:

$$ 
\begin{aligned}
U_k(\omega) &= 2g_+g_-[(u_k^{xx}u_k^{yy} - u_k^{xy}u_k^{xy}) + (v_k^{xx}v_k^{yy} - v_k^{xy}v_k^{xy})] + 4g_+^2g_-^2(u_k^{xx}u_k^{yy} - u_k^{xy}u_k^{xy})|f_k|^2, \\
V_k(\omega) &= (u_k^xv_k^y - u_k^yv_k^x)^2 \left[ g_+^2g_-^2[1 + 3\left(\frac{1-\alpha^2}{1+\alpha^2}\right)^2] + 4g_+^3g_-^3|f_k|^2 \right].
\end{aligned} 
\quad (14) $$

where $g_{\pm}(\omega) = 1/(\omega - \epsilon_{\pm,k})$. Interestingly, these last expressions show that $U_k(\omega)$ and $V_k(\omega)$ have poles only on the two dispersive bands as if the flat band did not play any role. In fact, the implicit effect of the flat band is essentially encoded in the $\alpha$ dependent term appearing in $V_k(\omega)$. However it is quite remarkable that $U_k(\omega)$ is totally independent of $\alpha$.

## 4. Orbital susceptibility of low-energy $\alpha$-$\mathcal{T}_3$ models

In this part of the paper we use equations (13,14) to calculate the susceptibility of $\alpha$-$\mathcal{T}_3$ model for three different $f_k$. In order to allow for a full analytical calculation, for each example we consider the low-energy effective model Hamiltonian that correctly describes the energy spectrum near the minimum of the corresponding $|f_k|$. The three distinct forms of $f_k$ that we consider are defined on the honeycomb lattice model illustrated on Fig. 3, that corresponds to the $\alpha = 0$ limit in which the C sites are completely decoupled from A, B sites. The C sites are not shown on Fig. 3 for simplicity.

### 4.1. $\alpha$-$\mathcal{T}_3$ graphene model

We first consider the usual isotropic model of graphene which corresponds to $t' = t$ and $t_3 = 0$ in Fig. 3. In this situation,

$$ f_k = t(e^{-ik \cdot \delta_1} + e^{-ik \cdot \delta_2} + e^{-ik \cdot \delta_3}) = t(e^{-ik_y a} + 2e^{ik_y a/2} \cos(\sqrt{3}k_x a/2)) \quad (15) $$

For $\Delta = 0$, the corresponding $\alpha$-$\mathcal{T}_3$ energy spectrum (see Fig. 1(b)) exhibits linear band touching at the two inequivalent corners $\pm \mathbf{K}$ of the Brillouin zone. Introducing the valley index $\xi = \pm$, the effective low-energy model is obtained by expanding $f_k$ near $\xi \mathbf{K} + \mathbf{k} \approx v(\xi k_x - ik_y)$ with the velocity $v = \frac{3ta}{2}$ where a is the nearest-neighbor distance. Hereafter to simplify the
notations of most equations we define the pseudo wavevector $\kappa_{x,y} = vk_{x,y}$. For each valley, the low-energy $\alpha$-$\mathcal{T}_3$ Hamiltonian reads:

$$h_k = \begin{pmatrix} \Delta & c_\alpha(\xi\kappa_x - i\kappa_y) & 0 \\ c_\alpha(\xi\kappa_x + i\kappa_y) & -\Delta & s_\alpha(\xi\kappa_x - i\kappa_y) \\ 0 & s_\alpha(\xi\kappa_x + i\kappa_y) & \Delta \end{pmatrix}. \quad (16)$$

It is then immediate to obtain

$$\vec{f} = (\xi\kappa_x, -\kappa_y), \quad \vec{f}^x = (\xi, 0), \quad \vec{f}^y = (0, -1), \quad \vec{f}^{ij} = 0. \quad (17)$$

and

$$u^x = \xi v_y = \frac{\kappa_x}{\sqrt{\kappa_x^2 + \kappa_y^2}}, \quad u^y = -\xi v_x = -\frac{\kappa_y}{\sqrt{\kappa_x^2 + \kappa_y^2}}, \quad u^{ij} = v^{ij} = 0, \quad (18)$$

from which we deduce

$$\begin{gathered} U_k(\omega) = 0, \\ V_k(\omega) = \left[ \frac{1}{(\omega^2 - (\Delta^2 + \kappa_x^2 + \kappa_y^2))^2} \left[ 1 + 3 \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 \right] + \frac{4(\kappa_x^2 + \kappa_y^2)}{(\omega^2 - (\Delta^2 + \kappa_x^2 + \kappa_y^2))^3} \right]. \end{gathered} \quad (19)$$

Substituting these expressions in (10), summing over the two valleys and noting $\kappa_x = \varepsilon \cos \theta$, $\kappa_y = \varepsilon \sin \theta$, we obtain

$$\chi_{\text{orb}}(\mu) = \frac{2\mu_0 e^2}{12} \left(\frac{4}{4\pi^2}\right) \frac{3m}{\pi} \int_{-\infty}^{+\infty} d\omega \, n_F(\omega) \int_0^{2\pi} d\theta \\ \times \int_0^{\infty} d\varepsilon \, \varepsilon \left[ \frac{1}{(\omega^2 - (\Delta^2 + \varepsilon^2))^2} \left[ 1 + 3 \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 \right] + \frac{4\varepsilon^2}{(\omega^2 - (\Delta^2 + \varepsilon^2))^3} \right]. \quad (20)$$

Integrating over variables and following the procedure described in appendix B, we finally obtain

$$\chi_{\text{orb}}(\mu = 0) = \frac{\mu_0 e^2 v^2}{2\pi} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{1}{3} \right] \frac{n_F(\Delta) - n_F(-\Delta)}{2\Delta}. \quad (21)$$

which gives the correct behaviour for $\Delta = 0$. For $\alpha = 0$, Eq.(21) with $\Delta = 0$ recovers the famous diamagnetic peak originally found by McClure [12], whereas Eq.(21) with $\Delta \neq 0$ reproduces the in-gap diamagnetic plateau first obtained by Koshino and Ando [13]. Note however that these authors obtained their results by performing a low-field expansion of the grand potential calculated from direct summation over the effective Landau levels spectrum $\epsilon_n(B) = \pm \epsilon_B \sqrt{|n|}$ [12] or $\epsilon_n(B) = \pm \sqrt{\Delta^2 + \epsilon_B^2|n|}$ [13] where $\epsilon_B = v\sqrt{2eB}$. For $\alpha \neq 0$, Eq.(21) is also coherent with the results of such methods, albeit now with $\alpha$-dependent effective Landau levels spectrum $\epsilon_n(B) = \pm \epsilon_B \sqrt{|n+\gamma|}$ or $\epsilon_n(B) = \pm \sqrt{\Delta^2 + \epsilon_B^2|n+\gamma|}$ where $n \in \mathbb{Z}$ and $\gamma = \frac{\alpha^2}{1+\alpha^2}$ [7]. In this picture, using the relation $\gamma = \frac{1}{2} - \Phi_B/2\pi$ where $\Phi_B$ is the $\alpha$-dependent Berry phase of a Dirac cone [7, 10], the variation of the susceptibility with $\alpha$ is interpreted as a consequence of the variation of the corresponding Berry phase. As we already pointed out, for $\alpha = 1$ the Berry phase vanishes and so do interband effects encoded by the Berry curvature. Nevertheless, for $\alpha = 1$ there is still a paramagnetic susceptibility plateau inside the gap which is an indication that interband effects are still present. This result implies that some important interband effects are not encoded by the Berry curvature. More quantitatively we note that the $\alpha$ dependent prefactor of the susceptibility (21) changes sign at $\alpha_c = \frac{\sqrt{3}-1}{\sqrt{2}} \approx 0.518$. This signals a transition from a diamagnetic peak/plateau for $\alpha < \alpha_c$ to a paramagnetic peak/plateau for $\alpha > \alpha_c$. This value $\alpha_c = 0.518$ coincides with the numerical results obtained from computing the susceptibility with the Hofstadter spectrum [7]. As a final remark, we note that for the low-energy Hamiltonian (16) which is linear in $k$ and thus separable (i.e. $h_k^{xy} = 0$), the Fukuyama formula [6] would have given the same results.
## 4.2. $\alpha$-$\mathcal{T}_3$ semi-Dirac model

We now consider a model with $t' = 2t$ and $t_3 = 0$ (see Fig. 3). In that situation

$$f_k = t(2e^{-ik_y a} + 2e^{ik_y a/2} \cos(\sqrt{3}k_x a/2)) \quad (22)$$

For $\Delta = 0$, the corresponding $\alpha$-$\mathcal{T}_3$ energy spectrum (see Fig. 4) exhibits a semi-Dirac band touching at one of the $M$ points of the Brillouin zone [14]. The low-energy model obtained by expanding near this point reads $f_{M+k} \approx \frac{k_x^2}{2m_*} - ivk_y$ with the effective mass $m_* = \frac{1}{ta^2}$ velocity $v = \frac{ta}{h}$ and it features a linear-quadratic spectrum $\varepsilon_k = \pm\sqrt{(\frac{k_x^2}{2m_*})^2 + (vk_y)^2}$. As in previous example, to simplify the notations of most equations we define the pseudo wavevectors $\kappa_x = \frac{k_x}{\sqrt{2m_*}}$ and $\kappa_y = vk_y$. The low-energy $\alpha$-$\mathcal{T}_3$ Hamiltonian is:

$$h_k = \begin{pmatrix} \Delta & c_\alpha(\kappa_x^2 - i\kappa_y) & 0 \\ c_\alpha(\kappa_x^2 + i\kappa_y) & -\Delta & s_\alpha(\kappa_x^2 - i\kappa_y) \\ 0 & s_\alpha(\kappa_x^2 + i\kappa_y) & \Delta \end{pmatrix}, \quad (23)$$

such that

$$\vec{f} = (\kappa_x^2, -\kappa_y), \quad \vec{f}x = 2(\kappa_x, 0), \quad \vec{f}y = (0, -1) \quad (24)$$

and

$$u^x = \frac{2\kappa_x^3}{\sqrt{\kappa_x^4 + \kappa_y^2}}, \quad u^y = \frac{k_y}{\sqrt{\kappa_x^4 + \kappa_y^2}}, \quad v^x = \frac{-2\kappa_x\kappa_y}{\sqrt{\kappa_x^4 + \kappa_y^2}}, \quad v^y = \frac{k_x^2}{\sqrt{\kappa_x^4 + \kappa_y^2}}, \quad (25)$$

from which we deduce

$$\begin{align} U_k(\omega) &= 0, \nonumber \\ V_k(\omega) &= 4\kappa_x^2 \left[ \frac{1}{(\omega^2-(\Delta^2+\kappa_x^4+\kappa_y^2))^2} \left[1+3\left(\frac{1-\alpha^2}{1+\alpha^2}\right)^2\right] + \frac{4(\kappa_x^4+\kappa_y^2)}{(\omega^2-(\Delta^2+\kappa_x^4+\kappa_y^2))^3} \right]. \tag{26} \end{align}$$

Writing now $\kappa_x^2 = \varepsilon \cos \theta$ and $\kappa_y = \varepsilon \sin \theta$ with $\theta \in [0, \pi/2]$ we deduce

$$\chi_{\text{orb}}(\mu) = \frac{\mu_0 e^2}{12} \left(\frac{64}{4\pi^2}\right) \frac{3m}{\pi} \int_{-\infty}^{+\infty} d\omega \, n_F(\omega) \int_0^{\pi/2} d\theta \sqrt{\cos\theta} \times \int_0^{\infty} d\epsilon \, \epsilon^{3/2} \left[ \frac{1}{(\omega^2-(\Delta^2+\epsilon^2))^2} \left[ 1 + 3 \left( \frac{1-\epsilon^2}{1+\epsilon^2} \right)^2 \right] + \frac{4\epsilon^2}{(\omega^2-(\Delta^2+\epsilon^2))^3} \right] \quad (27)$$

Quite interestingly, apart from some prefactor, the main noticeable change between equations (20) and (27) is the effective density of states that appears in the integral: for the $\alpha$-$\mathcal{T}_3$ semi-Dirac model it is $\propto \varepsilon^{3/2}$ whereas it is $\propto \varepsilon$ in the $\alpha$-$\mathcal{T}_3$ graphene model. From formulae (B.4,B.9) given in appendix B, we finally obtain:

$$\chi_{\text{orb}}(\mu) = -\frac{\mu_0 e^2 v}{\sqrt{2m_*} \pi^2 \sqrt{2\pi}} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{1}{2} \right] \times \begin{cases} \frac{1}{\sqrt{|\mu|}}, & \Delta = 0, T = 0 \\ \left( \frac{4\Gamma(\frac{5}{4})^2}{\sqrt{\pi}} \right) \frac{1}{\sqrt{\Delta}}, & |\mu| < \Delta, T = 0. \end{cases} \quad (28)$$

For $\alpha = 0$ and $\Delta = 0$, the first expression coincides with the result recently obtained in [1] for a two band model at the semi-Dirac point. For any $\alpha$ it predicts a diamagnetic peak that scales like $1/\sqrt{\min(\mu, T)}$ in the gapless case $\Delta = 0$. In the presence of a gap, the second expression shows that the susceptibility plateau in the gap scales as $1/\sqrt{\Delta}$. As in the previous example, Eq.(28) predicts a diamagnetic to paramagnetic transition for $\alpha \ge \alpha_c$ at a critical value $\alpha_c = \sqrt{2} - 1 \approx 0.414$ which is smaller than in then previous example. We have verified that this value $\alpha_c = 0.414$ agrees with the numerical results obtained from computing the susceptibility
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (Color online). $\alpha$-$\mathcal{T}_3$ semi-Dirac model energy bands dispersion in $k$ space, for $\Delta = 0$. (a) Full spectrum with semi-Dirac band touching of the two dispersive bands. (b) low-energy spectrum.</figcaption></figure>

with the Hofstadter spectrum. We stress that for the $\alpha$-$\mathcal{T}_3$ semi-Dirac model, as yet, there is no exact analytical expression for the Landau levels $\epsilon_n(B)$ of the low-energy model. The Landau levels are related to the eigenvalues of a modified quartic oscillator [14]. However by computing the Hofstadter spectrum of this $\alpha$-$\mathcal{T}_3$ semi-Dirac model for small magnetic field and various $\alpha$, we have observed that the effective Landau levels $\epsilon_n(B)$ are well described by the approximate form $\epsilon_n \propto [(n + 1/2)B]^{2/3}$, first obtained in [14]. Interestingly, for $n \ge 3$ the numerical Landau levels $\epsilon_n(B)$ appear to be almost independent of $\alpha$. This observation implies that the dia- to paramagnetic transition is essentially driven by the variation of the $n = 0, 1, 2$ Landau levels. We also stress that for this semi-Dirac model, there is no Berry phase around the semi-Dirac point. However the Berry curvature is non-zero (for $\alpha \ne 1$) and it exhibits a double-peak structure near the semi-Dirac point. There is still a lack of a clear physical picture of the origin of the dia to paramagnetic transition for this $\alpha$-$\mathcal{T}_3$ semi-Dirac model. As a final remark, we note that since the low-energy Hamiltonian (23) is still separable (i.e. $h_k^{xy} = 0$), the Fukuyama formula [6] would have given the same result.

### 4.3. $\alpha$-$\mathcal{T}_3$ pseudo-bilayer model

The last model we consider has parameters $t' = t$ and $t_3 = 1/2$ in Fig. 5 [15]. In that situation,

$$f_k = t(e^{-ik_y a} + 2e^{ik_y a/2} \cos(\sqrt{3}k_x a/2)) + \frac{t}{2}(e^{2ik_y a} + 2e^{-ik_y a} \cos(\sqrt{3}k_x a)) \quad (29)$$

For $\Delta = 0$, the corresponding $\alpha$-$\mathcal{T}_3$ energy spectrum (see Fig. 5) exhibits quadratic band touching at $\pm \mathbf{K}$ similar to a bilayer graphene. The low-energy model obtained by expanding near $\xi \mathbf{K}$ reads $f_{\xi \mathbf{K}+\mathbf{k}} \approx \frac{1}{2m_*}(\xi k_x - i k_y)^2$ with the effective mass $m_* = \frac{2}{ta^2}$ and the energy spectrum $\varepsilon_k = \frac{k^2}{2m_*}$. As in previous sections, to simplify the notations of most equations, we define the pseudo wavevectors $\kappa_{x,y} = \frac{k_{x,y}}{\sqrt{2m_*}}$.

The low-energy $\alpha$-$\mathcal{T}_3$ Hamiltonian becomes:

$$h_k = \begin{pmatrix} \Delta & c_\alpha (\xi \kappa_x - i \kappa_y)^2 & 0 \\ c_\alpha (\xi \kappa_x + i \kappa_y)^2 & -\Delta & s_\alpha (\xi \kappa_x - i \kappa_y)^2 \\ 0 & s_\alpha (\xi \kappa_x + i \kappa_y)^2 & \Delta \end{pmatrix}. \quad (30)$$
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (Color online). $\alpha$-$\mathcal{T}_3$ pseudo-bilayer model energy bands dispersion in $\mathbf{k}$ space, for $\Delta = 0$. (a) Full energy spectrum with quadratic band touching. (b) low-energy spectrum.</figcaption></figure>

This effective Hamiltonian is now quadratic in $\kappa_{x,y}$ and it is not separable (i.e. $h_k^{x,y} \neq 0$). As a consequence, we do not expect the Fukuyama formula to give the correct result in that situation. Following similar steps as in previous examples we obtain

$$ \begin{aligned} \vec{f} &= (\kappa_x^2 - \kappa_y^2, -2\xi\kappa_x\kappa_y), & \vec{f}^x &= 2(\kappa_x, -\xi\kappa_y), & \vec{f}^y &= 2(-\kappa_y, -\xi\kappa_x), \\ \vec{f}^{xx} &= 2(1, 0), & \vec{f}^{yy} &= 2(-1, 0), & \vec{f}^{xy} &= 2(0, -\xi), \end{aligned} \quad (31) $$

such that by noting $\kappa_x = \sqrt{\varepsilon} \cos \theta$, $\kappa_y = \sqrt{\varepsilon} \sin \theta$ we deduce

$$ \begin{aligned} U_{\mathbf{k}}(\omega) &= -16 \left[ \frac{1}{(\omega^2 - (\Delta^2 + \varepsilon^2))} + \frac{\varepsilon^2}{(\omega^2 - (\Delta^2 + \varepsilon^2))^2} \right], \\ V_{\mathbf{k}}(\omega) &= 16\varepsilon^2 \left[ \frac{1}{(\omega^2 - (\Delta^2 + \varepsilon^2))^2} \left[ 1 + 3 \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 \right] + \frac{4\varepsilon^2}{(\omega^2 - (\Delta^2 + \varepsilon^2))^3} \right]. \end{aligned} \quad (32) $$

Substituting these expressions in (10), summing over the two valley, integrating over $\theta$ we find

$$ \begin{aligned} \chi_{\text{orb}}(\mu) &= \frac{\mu_0 e^2}{12} \left(\frac{8}{\pi}\right) \frac{3m}{\pi} \int_{-\infty}^{+\infty} d\omega \, n_F(\omega) \\ & \times \int_0^{\infty} d\varepsilon \left[ \frac{\varepsilon^2}{(\omega^2 - (\Delta^2 + \varepsilon^2))^2} \left[ 5 + 12 \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 \right] + \frac{16\varepsilon^4}{(\omega^2 - (\Delta^2 + \varepsilon^2))^3} + \frac{1}{(\omega^2 - (\Delta^2 + \varepsilon^2))^3} \right]. \end{aligned} \quad (33) $$

From Eqs.(B.4,B.9) given in appendix B, we finally obtain:

$$ \chi_{\text{orb}}(\mu) = -\frac{\mu_0 e^2}{\pi m_*} \times \begin{cases} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{3}{4} \right] \ln \frac{\epsilon_c}{|\mu|} - \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{11}{12} \right] & \Delta = 0, |\mu| < \epsilon_c, T = 0 \\ \\ \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{3}{4} \right] \ln \frac{2\epsilon_c}{\Delta} - \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{11}{12} \right] & |\mu| < \Delta, T = 0 \end{cases} \quad (34) $$

where it is necessary to introduce a cutoff scale $\epsilon_c$ [16, 17]. As in previous examples, there is a diamagnetic to paramagnetic transition of the logarithmic term at a critical value $\alpha_c = \sqrt{4-\sqrt{3}} = 0.263$, whereas there is a paramagnetic to diamagnetic transition of the second term at $\alpha_c = \sqrt{12-\sqrt{11}} = 0.147$. For $\alpha = 0$, these expressions coincide with previous results obtained in [16, 17]. In [17], the susceptibility was derived from using the analytical form of the effective Landau levels $\epsilon_n = \epsilon_B \sqrt{n(n-1)}$ ($\epsilon_B = \frac{eB}{m_*}$) that are associated to the low-energy model. A similar derivation for any $\alpha$ seems difficult to achieve even if the Landau level
sequence is known. In that perspective, we note that for the gapless $\alpha$-$\mathcal{T}_3$ pseudo-bilayer model it is possible to derive a topological Berry phase that varies with $\alpha$; this already gives some insight on the modified semiclassical Landau levels spectrum [10]. However as noted in [10], for the bilayer system the semiclassical Landau levels spectrum does not fully agree with the exact quantum Landau levels spectrum and therefore we do not expect the semiclassical Landau levels spectrum to be sufficient to describe correctly the susceptibility –especially the singular behaviour of the susceptibility near half filling $\mu = 0$.

## 5. Summary

In this work we have studied the importance of interband effects on the orbital susceptibility of three bands $\alpha$-$\mathcal{T}_3$ tight-binding models that were initiated in [7]. The particularity of the $\alpha$-$\mathcal{T}_3$ tight-binding models is that the coupling between the three energy bands (which is encoded in the wavefunctions properties) can be tuned (by a parameter $\alpha$) without any modification of the energy spectrum. To highlight the role of these interbands effects, we have examined the orbital susceptibility of these models. Using the gauge-invariant perturbative formalism that we have recently developed [1], we obtained the generic formula of the orbital susceptibility of the $\alpha$-$\mathcal{T}_3$ tight-binding models. More quantitatively we have calculated the orbital susceptibility of three distinct $\alpha$-$\mathcal{T}_3$ tight-binding models: $\alpha$-$\mathcal{T}_3$ graphene model, the $\alpha$-$\mathcal{T}_3$ semi-Dirac model and the $\alpha$-$\mathcal{T}_3$ pseudo-bilayer model. To obtain an analytical form of the susceptibility, we have only considered a low-energy Hamiltonian of these models; it correctly describes the energy spectrum near the band touching where it is expected that interband effects are the strongest. The main result of our study is that for each of these models, by varying the parameter $\alpha$ and thus the wavefunctions interband coupling, it is possible to drive a transition from a diamagnetic to a paramagnetic behavior, both in absence or presence of a gap separating the dispersive bands. In particular we emphasize that the in-gap susceptibility does not need to be diamagnetic. Moreover the existence of a finite in-gap (paramagnetic) susceptibility at $\alpha = 1$ (for each model) provides hints that some important interband effects are not encoded in the Berry curvature which vanishes for $\alpha = 1$.

## Appendix A. Determination of $U_k(\omega)$ and $V_k(\omega)$

The aim of this section is to give the key step to find the expressions $U_k(\omega)$ and $V_k(\omega)$ of Eq. (14). We first slightly rewrite their definitions Eq. (11):

$$ 
\begin{aligned}
U_k(\omega) &= \operatorname{tr}\{(gh^{xy}gh^{yy} - gh^{xy}gh^{xy})_k\}, \\
V_k(\omega) &= \operatorname{tr}\{(gh^xgh^y[gh^y, gh^x])_k\},
\end{aligned}
\qquad (A.1) 
$$

where $g_k(\omega) = (\omega I - h_k)^{-1}$ it the Green's function matrix associated to the zero-field Bloch Hamiltonian matrix Eq. (2), $h_k^i = \frac{\partial h_k}{\partial k_i}$ and $h_k^{ij} = \frac{\partial^2 h_k}{\partial k_i \partial k_j}$ with $(i, j) = (x, y)$. The calculation of $U_k, V_k$ is made simple by remarking that the matrices $h, g, h^j, h^{ij}$ can each be written in terms of the following three matrices only:

$$ S_{\pm} = \begin{pmatrix} 0 & \pm c_{\alpha} e^{-i\theta_k} & 0 \\ c_{\alpha} e^{i\theta_k} & 0 & \pm s_{\alpha} e^{-i\theta_k} \\ 0 & s_{\alpha} e^{i\theta_k} & 0 \end{pmatrix}; \quad S_0 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & -1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \qquad (A.2) $$

Note that the two matrices $S_{\pm}$ depend on both parameter $\alpha$ and wavevector $k$; moreover we stress that $S_{-}$ is an antihermitian matrix whereas $S_{+,0}$ are hermitian. With these definitions, it is immediate to verify the identities

$$ 
\begin{aligned}
h &= \Delta S_0 + |f_k|S_+, \\
h^j &= u_k^j S_+ + iv_k^j S_-, \\
h^{ij} &= u_k^{ij} S_+ + iv_k^{ij} S_-,
\end{aligned}
\qquad (A.3) 
$$
with

$$ 
\begin{aligned}
u_{\mathbf{k}}^{j} &= \frac{\vec{f}_{\mathbf{k}}^{j} \cdot \vec{f}_{\mathbf{k}}}{|f_{\mathbf{k}}|}, & v_{\mathbf{k}}^{j} &= \frac{(\vec{f}_{\mathbf{k}} \times \vec{f}_{\mathbf{k}}^{j})}{|f_{\mathbf{k}}|}, \\
u_{\mathbf{k}}^{ij} &= \frac{\vec{f}_{\mathbf{k}}^{ij} \cdot \vec{f}_{\mathbf{k}}}{|f_{\mathbf{k}}|}, & v_{\mathbf{k}}^{ij} &= \frac{(\vec{f}_{\mathbf{k}} \times \vec{f}_{\mathbf{k}}^{ij})}{|f_{\mathbf{k}}|},
\end{aligned} 
\quad (A.4) $$

where we have defined the two components vectors ($\vec{x} \equiv (x_1, x_2)$)

$$ 
\begin{aligned}
\vec{f}_{\mathbf{k}} &= (\Re e f_{\mathbf{k}}, \Im m f_{\mathbf{k}}), \\
\vec{f}_{\mathbf{k}}^{i} &= \left( \frac{\partial \Re e f_{\mathbf{k}}}{\partial k_i}, \frac{\partial \Im m f_{\mathbf{k}}}{\partial k_i} \right), \\
\vec{f}_{\mathbf{k}}^{ij} &= \left( \frac{\partial^2 \Re e f_{\mathbf{k}}}{\partial k_i \partial k_j}, \frac{\partial^2 \Im m f_{\mathbf{k}}}{\partial k_i \partial k_j} \right).
\end{aligned} 
\quad (A.5) $$

To obtain a simple form for $g$ we remark that the matrices $S_{\pm,0}$ appear to have very peculiar properties: $S_{\pm}^3 = \pm S_{\pm}$, $S_0^2 = I$, $[S_{\pm}, S_0]_+ = 0$. From these properties we deduce that $g = (\omega I - h_{\mathbf{k}})^{-1}$ can be written as $g = a_1 I + a_2 S_0 + a_3 S_+ + a_4 S_+^2$ where $a_{1,2,3,4}$ are determined from using the identity $(\omega - h)g = I$. We then obtain:

$$ g = g_+ g_- (\omega I - \Delta S_0 - |f_{\mathbf{k}}| S_+) - g_+ g_- g_0 |f_{\mathbf{k}}|^2 (I - S_+^2), \quad (A.6) $$

where $g_{\pm \mathbf{k}}(\omega) = \frac{1}{\omega - \epsilon_{\pm, \mathbf{k}}}$ and $g_0(\omega) = \frac{1}{\omega - \epsilon_0}$. Note that $g(\omega)$ has three poles corresponding to the three bands. Substituting identities Eqs.(A.3,A.6) in Eq.(A.1), we obtain the expressions Eq.(14):

$$ 
\begin{aligned}
U_{\mathbf{k}}(\omega) &= 2g_{+}g_{-}[(u_{\mathbf{k}}^{xx}u_{\mathbf{k}}^{yy} - u_{\mathbf{k}}^{xy}u_{\mathbf{k}}^{xy}) + (v_{\mathbf{k}}^{xx}v_{\mathbf{k}}^{yy} - v_{\mathbf{k}}^{xy}v_{\mathbf{k}}^{xy})] + 4g_{+}^2g_{-}^2(u_{\mathbf{k}}^{xx}u_{\mathbf{k}}^{yy} - u_{\mathbf{k}}^{xy}u_{\mathbf{k}}^{xy})|f_{\mathbf{k}}|^2, \\
V_{\mathbf{k}}(\omega) &= (u_{\mathbf{k}}^{x}v_{\mathbf{k}}^{y} - u_{\mathbf{k}}^{y}v_{\mathbf{k}}^{x})^2 \left[ g_{+}^2g_{-}^2[1 + 3\left(\frac{1-\alpha^2}{1+\alpha^2}\right)^2] + 4g_{+}^3g_{-}^3|f_{\mathbf{k}}|^2 \right].
\end{aligned} 
\quad (A.7) $$

## Appendix B. Pole decomposition and integration over variables $\epsilon$ and $\omega$

The generic form of the susceptibility of the three different $\alpha$-$\mathcal{T}_3$ models is:

$$ 
\begin{aligned}
\chi_{\text{orb}} &\propto \frac{\Im m}{\pi} \int d\omega \, n_{\text{F}}(\omega) \int_{0}^{\infty} d\epsilon \, \nu(\epsilon) \left[ \frac{A\epsilon^4}{(\omega^2-(\Delta^2+\epsilon^2))^3} + \frac{B\epsilon^2}{(\omega^2-(\Delta^2+\epsilon^2))^2} + \frac{C}{(\omega^2-(\Delta^2+\epsilon^2))} \right] \\
&\propto \frac{\Im m}{\pi} \int d\omega \, n_{\text{F}}(\omega) \int_{0}^{\infty} d\epsilon \, \nu(\epsilon) \left[ A\epsilon^4(g_+g_-)^3 + B\epsilon^2(g_+g_-)^2 + C(g_+g_-) \right]
\end{aligned} 
\quad (B.1) $$

where $g_s(\omega) = 1/(\omega - \epsilon_s)$ and with $\epsilon_s = s\sqrt{\Delta^2 + \epsilon^2}$ ($s = \pm$). The effective density of states $\nu(\epsilon)$ and the parameters $A, B, C$ of the three models ($\alpha$-graphene, $\alpha$-semi-Dirac, $\alpha$-bilayer) are summarized in table B1.

<table><caption>**Table B1.** Effective density of states $\nu(\epsilon)$ and parameters $A, B, C$ of the three considered models: $\alpha$-graphene, $\alpha$-semi-Dirac and $\alpha$-bilayer.</caption>
  <thead>
    <tr>
      <th></th>
      <th>α-graphene</th>
      <th>α-semi-Dirac</th>
      <th>α-bilayer</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>ν(ε)</th>
      <td>1/ε</td>
      <td>1/√ε</td>
      <td>1</td>
    </tr>
    <tr>
      <th>A</th>
      <td>4</td>
      <td>4</td>
      <td>16</td>
    </tr>
    <tr>
      <th>B</th>
      <td>1 + 3 (1-α<sup>2</sup>)/(1+α<sup>2</sup>)<sup>2</sup></td>
      <td>1 + 3 (1-α<sup>2</sup>)/(1+α<sup>2</sup>)<sup>2</sup></td>
      <td>5 + 12 (1-α<sup>2</sup>)/(1+α<sup>2</sup>)<sup>2</sup></td>
    </tr>
    <tr>
      <th>C</th>
      <td>0</td>
      <td>0</td>
      <td>1</td>
    </tr>
  </tbody>
</table>
We now describe the key steps to perform the explicit integration over variables $\epsilon$ and $\omega$. The first step consists to separate the poles at $\epsilon_s$:

$$ \chi_{\text{orb}} \propto \sum_{s=\pm} \frac{3m}{\pi} \int d\omega \, n_F(\omega) \int_0^\infty d\epsilon \, \nu(\epsilon) \left[ \frac{A}{8} \frac{\epsilon^4}{\epsilon_s^3} g_s^3 + \left( \frac{B}{4} \frac{\epsilon^2}{\epsilon_s^2} - 3 \frac{A}{16} \frac{\epsilon^4}{\epsilon_s^4} \right) g_s^2 + \left( \frac{C}{2\epsilon_s} - \frac{B}{4} \frac{\epsilon^2}{\epsilon_s^3} + 3 \frac{A}{16} \frac{\epsilon^4}{\epsilon_s^5} \right) g_s \right] \quad (B.2) $$

The next step consists to use the identity $\frac{3m}{\pi} \int d\omega \, n_F(\omega) g_s^n(\omega) = -\frac{1}{(n-1)!} n_F^{n-1}(\epsilon_s)$ where $n_F^{(n)}$ is the $n^{th}$ derivative; we can rewrite:

$$ \chi_{\text{orb}} \propto - \sum_{s=\pm} \int_0^\infty d\epsilon \, \nu(\epsilon) \left[ \frac{A}{16} \frac{\epsilon^4}{\epsilon_s^3} n_F''(\epsilon_s) + \left( \frac{B}{4} \frac{\epsilon^2}{\epsilon_s^2} - 3 \frac{A}{16} \frac{\epsilon^4}{\epsilon_s^4} \right) n_F'(\epsilon_s) + \left( \frac{C}{2\epsilon_s} - \frac{B}{4} \frac{\epsilon^2}{\epsilon_s^3} + 3 \frac{A}{16} \frac{\epsilon^4}{\epsilon_s^5} \right) n_F(\epsilon_s) \right] \quad (B.3) $$

From this point, we separate the study of case $\Delta = 0$ from case $\Delta \neq 0$.

### Gapless models, $\Delta = 0$;

For $\Delta = 0$, since $\epsilon_s = s\epsilon$ we first rewrite:

$$ \chi_{\text{orb}} \propto - \int_0^\infty d\epsilon \, \nu(\epsilon) \left[ \frac{A}{16}\epsilon(n_F''(\epsilon) - n_F''(-\epsilon)) + \left(\frac{B}{4} - 3\frac{A}{16}\right)(n_F'(\epsilon) + n_F'(-\epsilon)) + \left(\frac{C}{2} - \frac{B}{4} + 3\frac{A}{16}\right) \frac{n_F(\epsilon) - n_F(-\epsilon)}{\epsilon} \right] \quad (B.4) $$

We now note that for the three considered cases, $\nu(\epsilon) \sim \epsilon^p$ with $-1 \le p \le 0$; this property permits an integration by part of the term proportional to $n_F''$ and $n_F$:

$$ \begin{aligned} \int_0^\infty d\epsilon \epsilon^{p+1} (n_F''(\epsilon) - n_F''(-\epsilon)) &= \left[ \epsilon^{p+1} (n_F'(\epsilon) + n_F'(-\epsilon)) \right]_0^\infty - (p+1) \int_0^\infty d\epsilon \epsilon^p (n_F'(\epsilon) + n_F'(-\epsilon)), \\ \int_0^\infty d\epsilon \epsilon^{p-1} (n_F(\epsilon) - n_F(-\epsilon)) &= \left[ \frac{\epsilon^p}{p} (n_F(\epsilon) - n_F(-\epsilon)) \right]_0^\infty - \int_0^\infty d\epsilon \frac{\epsilon^p}{p} (n_F'(\epsilon) + n_F'(-\epsilon)), \end{aligned} \quad (B.5) $$

where the last line requires $p < 0$. Using these identities for the three considered cases ($p = -1, -1/2, 0$) and using the parameters of table B1 we deduce:

$$ \begin{aligned} \chi_{\text{orb}} &\propto \frac{3}{2} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{1}{3} \right] n_F'(0), & \alpha\text{-graphene} \\ \chi_{\text{orb}} &\propto \frac{3}{4} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{1}{2} \right] \int_0^\infty d\epsilon \frac{n_F'(\epsilon) + n_F'(-\epsilon)}{\sqrt{\epsilon}} & \alpha\text{-semi-Dirac} \\ \chi_{\text{orb}} &\propto 3 \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{3}{4} \right] \int_0^\infty d\epsilon \frac{n_F(\epsilon) - n_F(-\epsilon)}{\epsilon} + \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{11}{12} \right] & \alpha\text{-bilayer} \end{aligned} \quad (B.6) $$

### Gapped models, $\Delta \neq 0$:

For the gapped case, we can also perform some integration by part of terms proportional to $n_F''$ and $n_F$. More precisely using the identities $\frac{\partial n_F^n}{\partial \epsilon} = \frac{\epsilon}{\epsilon_s} n_F^{n+1}$ and $\frac{\partial \epsilon_s}{\partial \epsilon} = \frac{\epsilon}{\epsilon_s}$ we first obtain

$$ \begin{aligned} \int_0^\infty d\epsilon \frac{\epsilon^{p+4}}{\epsilon_s^3} n_F'' &= - \int_0^\infty d\epsilon \frac{\epsilon^{p+2}}{\epsilon_s^2} (p + 3 - 2 \frac{\epsilon^2}{\epsilon_s^2}) n_F', \\ \int_0^\infty d\epsilon \frac{\epsilon^{p+2}}{\epsilon_s^3} n_F &= \left[ - \frac{\epsilon^{p+1}}{\epsilon_s} n_F \right]_0^\infty + \int_0^\infty d\epsilon \left[ \frac{(p+1)\epsilon^p}{\epsilon_s} n_F + \frac{\epsilon^{p+2}}{\epsilon_s^2} n_F' \right], \\ \int_0^\infty d\epsilon \frac{\epsilon^{p+4}}{\epsilon_s^5} n_F &= \frac{1}{3} \left[ - \frac{\epsilon^{p+3}}{\epsilon_s^3} n_F \right]_0^\infty + \frac{1}{3} \int_0^\infty d\epsilon \left[ \frac{(p+3)\epsilon^{p+2}}{\epsilon_s^3} n_F + \frac{\epsilon^{p+4}}{\epsilon_s^4} n_F' \right], \\ &= \left[ - \frac{\epsilon^{p+3}}{\epsilon_s^3} n_F \right]_0^\infty + \frac{p+3}{3} \left[ - \frac{\epsilon^{p+1}}{\epsilon_s} n_F \right]_0^\infty + \frac{1}{3} \int_0^\infty d\epsilon \left[ \frac{(p+3)(p+1)\epsilon^p}{\epsilon_s} n_F + \frac{(p+3)\epsilon^{p+2}}{\epsilon_s^2} n_F' + \frac{\epsilon^{p+4}}{\epsilon_s^4} n_F' \right] \end{aligned} \quad (B.7) $$
from which we deduce:

$$ \chi_{\text{orb}} \propto -\sum_{s=\pm} \left[ \frac{A}{16} \left[ -\frac{\epsilon^{p+3}}{\epsilon_s^3} n_F \right]_0^\infty + \left[ \frac{B}{4} - \frac{A}{16}(p+3) \right] \left[ \frac{\epsilon^{p+1}}{\epsilon_s} n_F \right]_0^\infty + \left( \frac{C}{2} - (p+1) \left[ \frac{B}{4} - \frac{A}{16}(p+3) \right] \right) \int_0^\infty d\epsilon \frac{\epsilon^p}{\epsilon_s} n_F(\epsilon_s) \right]. \quad (B.8) $$

For the three considered cases ($p = -1, -1/2, 0$) with the parameters of table B1 we finally obtain:

$$ \begin{aligned} \chi_{\text{orb}} &\propto \frac{3}{4} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{1}{3} \right] \frac{n_F(\Delta)-n_F(-\Delta)}{\Delta} & \alpha\text{-graphene,} \\ \chi_{\text{orb}} &\propto \frac{3}{8} \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{1}{2} \right] \int_{\Delta}^{\infty} d\epsilon \frac{n_F(\sqrt{\Delta^2+\epsilon^2})-n_F(-\sqrt{\Delta^2+\epsilon^2})}{\sqrt{\epsilon}\sqrt{\Delta^2+\epsilon^2}} & \alpha\text{-semi-Dirac} \\ \chi_{\text{orb}} &\propto 3 \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{3}{4} \right] \int_{\Delta}^{\infty} d\epsilon \frac{n_F(\sqrt{\Delta^2+\epsilon^2})-n_F(-\sqrt{\Delta^2+\epsilon^2})}{\sqrt{\Delta^2+\epsilon^2}} - \left[ \left( \frac{1-\alpha^2}{1+\alpha^2} \right)^2 - \frac{11}{12} \right] & \alpha\text{-bilayer} \end{aligned} \quad (B.9) $$
[1] Raoux A Piéchon F Fuchs J-N and Montambaux G 2014 *arXiv:1411.5940*, accepted in Phys. Rev. B.

[2] Landau L D 1930 *Z. Phys.* **64**, 629.

[3] Peierls R 1933 *Z. Phys.* **80**, 763.

[4] Vignale G 1991 Phys. Rev. Lett **67**, 358.

[5] Fukuyama H and Kubo R 1970 *J. Phys. Soc. Jap.* **28**, 570.

[6] Fukuyama H 1970 *Phys. Lett.* **32A**, 111; Fukuyama H 1971 *Prog. Theor. Phys.* **45**, 704.

[7] Raoux A Morigi M Fuchs J-N Piéchon F and Montambaux G 2014 *Phys. Rev. Lett.* **112**, 026402.

[8] Sutherland B 1986 *Phys. Rev. B* **34**, 5208.

[9] Vidal J Mosseri R and Doucot B, 1998 *Phys. Rev. Lett.* **81**, 5888.

[10] Fuchs J-N Piéchon F Goerbig M-O and Montambaux G 2010 *Eur. Phys. J. B* **77** 351.

[11] Gómez-Santos G and Stauber T 2011 *Phys. Rev. Lett.* **106**, 045504.

[12] McClure J W 1956 *Phys. Rev.* **104**, 666.

[13] Koshino M and Ando T 2007 *Phys. Rev. B* **76**, 085425.

[14] Dietl P Piéchon F and Montambaux G 2008 *Phys. Rev. Lett.* **100**, 236405.

[15] Montambaux G 2012 *Eur. Phys. J. B* **85**, 375.

[16] Safran S A 1984 *Phys. Rev. B* **30**, 421.

[17] Koshino M 2011 *Orbital Magnetism of Graphenes*, in *Physics and Applications of Graphene - Theory*, Mikhailov S (Ed.), InTech.