Effect of disorder on 2D topological merging transition from a Dirac semi-metal to a normal insulator

For $\bar {\Delta } \geqslant 0$, integrating over q in eq. (5a) we obtain

$$\begin{equation} X_0 = \varepsilon + \frac{ \sqrt{2} \alpha X_0 }{\sqrt{\bar{\Delta}-X_0}} \, K\left(2+\frac{2 \bar{\Delta} }{X_0-\bar{\Delta}}\right), \end{equation} \tag{ 8 }$$

where $K(z)=\int _0^{1} {{\mathrm { d}}t} [{{(1-t^2)(1-z t^2)}}]^{-1/2}$ is the complete elliptic integral of the first kind¹. Equation (8) has two complex conjugate solutions corresponding to the retarded and advanced functions. We always find for $\bar {\Delta } > \bar {\Delta }_c=\pi ^2 \alpha ^2/2$ the DOS of an insulator: there exists an energy gap $\varepsilon _c=\bar {\Delta } + O(\alpha )$ such that the solution of eq. (8) is purely real for 0 ⩽ ε ⩽ ε_c, and thus, the DOS ρ(ε) vanishes. The DOS computed using eq. (7) and the imaginary part of the numerical solution of eq. (8) is shown in fig. 3. For |ε − ε_c| ≫ 0 one can solve eq. (8) by iterations: this provides an expansion of the solution in powers of disorder strength. The first iteration derived by taking X₀ → ε + i0⁺ in the r.h.s. of eq. (8) gives a linear in disorder correction and thus the DOS of the clean system found in [13] up to the shift $\Delta \to \bar {\Delta }$. It is also shown in fig. 3 for comparison. To study the solution of the SCBA equation in the vicinity of the energy gap ε_c, we replace the complete elliptic integral in eq. (8) by its asymptotics that yields

$$\begin{equation} X_0 = \varepsilon +\frac{\alpha}{2} \sqrt{{X_0 }} \left[\ln \left(\frac{32 X_0}{X_0-\bar{\Delta}}\right)\pm i \pi \right], \end{equation} \tag{ 9 }$$

where " + /–" correspond to the retarded and advanced functions, respectively. Equation (9) can be written as X₀ = ε + f(X₀) so that the energy gap ε_c and the corresponding value X₀(ε_c) = X_c satisfy: X_c = ε_c + f(X_c) and f'(X_c) = 1. In the limit of weak disorder we obtain

$$\begin{equation} \varepsilon_c \!=\! \bar{\Delta} \left\{1 \!-\! \frac{\alpha}{2 \sqrt{\bar{\Delta}}}\left[ 1 \!-\! \ln \left(\frac{\alpha}{64 \sqrt{\bar{\Delta}}}\right) \right]\right\} \!+\! O({\alpha}^2 \ln {\alpha}), \end{equation} \tag{ 10 }$$

and the expansion of the DOS around this point for ε ≳ ε_c:

$$\begin{equation} \rho (\varepsilon) = \frac{\sqrt{m}}{\alpha c \pi^2 } \left[ \frac{2(\varepsilon-\varepsilon_c)}{f''(X_{c})} \right]^{1/2} + O((\varepsilon-\varepsilon_c)^{3/2}). \end{equation} \tag{ 11 }$$

Thus, the DOS of the disordered system in the insulating phase vanishes as a power law with the exponent $\frac 12$ in contrast to the jump in the DOS of the clean system (see the bottom plot in fig. 3).

Zoom In Zoom Out Reset image size

In the case $\bar {\Delta }<0$, the integration over q in eq. (5a) leads to the equation

$$\begin{equation} X_0 = \varepsilon + \frac{ \sqrt{2} \alpha X_0 }{\sqrt{\bar{\Delta}+X_0}} \left[ K\left(\frac{2X_0}{\bar{\Delta}+X_0}\right) + 2 i K\left(\frac{\bar{\Delta}-X_0}{\bar{\Delta}+X_0}\right) \right]. \end{equation} \tag{ 12 }$$

The imaginary part of the numerical solution of eq. (12) gives the DOS which is shown in the upper plot of fig. 3. Iterating eq. (12) one can obtain the solution for $|\varepsilon \!-\!\bar {\Delta }| \gg 0$. The DOS of the clean system found in [13] (up to the shift (6)) is recovered by taking the limit X₀ → ε + i0⁺ in the r.h.s. of eq. (12) and substituting the solution to eq. (7). This DOS diverges logarithmically at $\varepsilon =|\bar {\Delta }|$ while the DOS of the disordered system is continuous and smooth (see the upper plot in fig. 3).

Let us now focus on the region of small $\overline {\Delta }$. The DOS computed numerically using the SCBA eqs. (8) and (12) is shown in fig. 4 as a function of $\bar {\Delta }$ in the vicinity of the transition $\bar {\Delta }=0$. One can see that in the presence of disorder there is a region around the transition point $\bar {\Delta }=0$ characterized by a finite DOS at zero energy, ρ(ε = 0) > 0. The DOS at ε = 0 strictly vanishes for $\bar {\Delta }>\bar {\Delta }_c=\pi ^2 \alpha ^2/2$. Thus, in the presence of disorder the transition (semi-Dirac) point between the Dirac semi-metal and the insulator is replaced by a new intermediate regime which we call the disordered semi-Dirac regime. The disordered Dirac phase smoothly evolves into this intermediate regime, while the SCBA suggests a sharp continuous transition at $\bar {\Delta }=\bar {\Delta }_c$ between this disordered semi-Dirac phase and the insulating phase. Exactly at the transition we find an algebraic behavior of the DOS

$$\begin{equation} \rho(\varepsilon;\bar{\Delta}=\bar{\Delta}_c)= \left(\frac34\right)^{1/6} \frac{\sqrt{m}}{ c \pi } ( \pi \alpha )^{1/3} {\varepsilon }^{1/3}, \end{equation} \tag{ 13 }$$

which is displayed in fig. 4. The extrapolation towards the disordered semi-metallic regime is obtained for large negative $\bar {\Delta }<0$: in this regime we recover a DOS at zero energy decaying as

$$\begin{eqnarray} &&\rho(0)= \frac{8|\bar{\Delta}|\sqrt{m}}{\alpha c \pi^2} e^{- \sqrt{|\bar{\Delta}|} /(\alpha\sqrt{2})}. \end{eqnarray} \tag{ 14 }$$

Hence, instead of a transition the system undergoes a crossover to disordered Dirac fermions. Indeed, this behavior is fully consistent with the SCBA picture for the disordered Dirac fermions. For graphene it predicts a finite DOS ρ(ε) = 2Γ₀/(π²v²₀α_D) for ε ≪ Γ₀ where Γ₀ = Δ_De^−1/α_D is the energy scale generated by disorder. Here Δ_D is the bandwidth, v₀ —the Fermi velocity and α_D —the dimensionless strength of disorder [27, 28]. This would imply a finite DOS at zero energy for all considered types of disorder. However, this contradicts the results obtained for random gauge and random mass disorder using more reliable methods. For instance, for random gauge disorder one expects ρ(ε) = ε^2/z−1 with z = 1 + α_D in weak-disorder case (α_D < 2) [23] and z = (8α_D)^1/2 − 1 in strong-disorder case (α_D > 2) [25]. Similarly, with random mass disorder, the DOS remains linear in ε up to logarithmic corrections [30]. Only for scalar potential disorder the DOS saturates at a finite value in the vicinity of the Dirac point [27]. It was argued in [24] that the failure of SCBA in the vicinity of the Dirac point for ε < Γ₀ is due to importance of the diagrams with crossed disorder lines, neglected within the SCBA which takes into account only the noncrossed ones. Another possible reason for the failure of the SCBA is the divergence of the fermion wavelengths at the Dirac point rendering the SCBA uncontrollable, i.e., without a small parameter. Nevertheless one can still rely on the SCBA for energies ε ≫ Γ₀. The question of reliability of the SCBA in the vicinity of the topological transition is even more subtle.

Zoom In Zoom Out Reset image size

It is interesting to note that the SCBA equations can be analytically solved exactly at the point $\bar {\Delta }=0$, inside the "disordered semi-Dirac" regime. Indeed using that $K(2+0^+/z)= C(1-\mathrm {sign}[\mathrm {Im}\, z]i)/ \sqrt {2}$ where $C = { \Gamma \left ({1}/{4}\right )^2}/{(4 \sqrt {\pi } )}$ one can rewrite eq. (8) as $X_0=\varepsilon +(1\mp i) C {\alpha } \sqrt {-X_0}$, where " − / + " corresponds to the retarded/advanced Green's function. We obtain the DOS at the point $\bar {\Delta }=0$

$$\begin{eqnarray} &&\rho (\varepsilon)= \frac{\sqrt{m} C^2 \alpha}{c \pi^2} \left[ 1+ \frac1{\sqrt{2}} \sqrt{1+\sqrt{1+ \frac{4 \varepsilon^2}{C^4 {\alpha}^4} }} \right], \end{eqnarray} \tag{ 15 }$$

such that $\rho (0)= 2{\sqrt {m} C^2 \alpha }/{(c \pi ^2)}$. Note that the expansion of eq. (15) in powers of the disorder strength α is in fact an expansion in ${\alpha }/{\sqrt {\varepsilon }}$. Its truncation to a finite order in α appears to be pathological for ε ≪ α². We interpret this as a sign of the probable breakdown of the SCBA for a small enough energy. In this regime, a more refined technique is required [26].

Before we consider a weak-disorder renormalization group (RG) approach to the topological transition let us first discuss the consequences of our SCBA results for the behavior of specific heat and spin magnetic susceptibility [13, 36]. Assuming that the chemical potential is temperature independent and fixed at zero energy the specific heat can be written as eq. (18) in [13]. In the clean system the low-T specific heat interpolates from a T² behavior far in the semi-metal phase to an activated behavior in the insulating phase. Exactly at the transition it exhibits a T^3/2 behavior. In the presence of disorder the low-T specific heat at the critical point should behave as T^4/3. The magnetic susceptibility χ is related to the DOS by eq. (10.11) in [37]. In the clean system, χ is activated in the insulator, χ ∼ T in the semi-metal phase and χ ∼ T^1/2 at the transition. With disorder, the scaling behavior of the susceptibility at the transition is modified to χ ∼ T^1/3.

The standard way to generalize the model to d dimensions is to introduce the generators of the d-dimensional Clifford algebra. In even dimensions these generators can be represented by the matrices γ_i (i = 1,...,d) of dimension 2^d/2. Generally one also defines γ_d+1 ≡ γ_S = i^−d/2γ₁...γ_d. For the case of odd dimensions one constructs the matrices γ_i (i = 1,...,d − 1) similarly to the case of even dimension and adds γ_d ≡ γ_S = i^−(d−1)/2γ₁...γ_d−1. These matrices satisfy the anticommutation relations: γ_iγ_j + γ_jγ_i = 2δ_ijγ₀, i,j = 1,...,[[d + 1]], where $\gamma _0=\mathbb {I}$ is the identity matrix and [[d + 1]] stands for d + 1 if d is even, and d if it is odd. The matrices γ_i and γ_S replace the Pauli matrices σ_i (i = x,y) and σ_z in d dimensions. We define the two matrix vectors γ_∥: = (γ₁,...,γ_l) and γ_⊥: = (γ_l+1,...,γ_d) and replace in Hamiltonian (1) ∂_x → ∂_∥ which acts in l-dimensional subspace and ∂_y → ∂_⊥ which acts in (d − l)-dimensional subspace. Introducing n copies of the original system and using the replica trick to average over disorder we can write the effective replicated action as

$$\begin{eqnarray} & &\hspace{-10pt} \mathcal{S} = \int {\rm d}^d \mathbf{r}\left \{ \sum\limits_{a=1}^n \bar{\psi}_{a}\left[ \gamma_{\parallel} \left(\Delta - \frac{\partial_{\|}^2}{2m}\right) - i c \gamma_{\bot} \partial_{\bot} - i \varepsilon \right]\psi_{a} \right.\nonumber \\ & &\hspace{-10pt}\left. - \frac{\pi c}{\sqrt{m}} \sum\limits_{a,b=1}^n \sum\limits_{j=0}^{[[d+1]]} \alpha_{j} [\bar{\psi}_{a}(\mathbf{r})\gamma_{j} \psi_{a}(\mathbf{r})] [\bar{\psi}_{b}(\mathbf{r})\gamma_{j} \psi_{b}(\mathbf{r})] \vphantom{\sum\limits_{a=1}^n \bar{\psi}_{a}} \right\}, \qquad \end{eqnarray} \tag{ 16 }$$

where d^dr: = d^lr_∥ d^d−lr_⊥. The properties of the initial disordered model are obtained in the limit n → 0. Dimensional analysis of the action (16) yields the bare dimensions: [r_∥] = μ^−1/2, [r_⊥] = μ⁻¹, $[\psi ]= [\bar {\psi }] = \mu ^{(2d-l-2)/4}$ and [α_j] = μ^δ, where μ is an arbitrary momentum scale and δ = 2 + l/2 − d. Thus the upper critical dimension of the model (16) is d_c = 2 + l/2 and disorder is i) relevant for δ > 0, e.g., for 2D semi-Dirac fermions with δ = 1/2; ii) marginally relevant for 2D Dirac fermions with δ = 0; iii) irrelevant for 3D Weyl (δ = −1) and semi-Weyl fermions (δ = −1/2).

We focus here on the case of 2D semi-Dirac fermions. To renormalize the model (16) one has to eliminate UV divergences in the loop expansion by introducing counter terms. The one-loop counterterms to the propagator and disorder are given by the diagrams shown in fig. 5. To one-loop order the m and c do not get corrected while Δ gets an additive shift similar to that in eq. (6). We use dimensional regularization in which this shift vanishes [34]. The diagrams (c) and (d) in fig. 5 can generate extra terms with new algebraic structures like γ_iγ_jγ_k⊗γ_iγ_jγ_k and so on. In d = d_c − δ dimensions they have to be treated as independent new interactions reflecting the fact that within dimensional regularization this theory is not multiplicatively renormalizable [39]. For the Dirac fermions this problem shows up to orders higher than one loop and can be solved, for example, using a cutoff in strictly 2D so that the extra interactions do not arise. In contrast, for the semi-Dirac fermions the extra interactions infinitely proliferate in the renormalization already to the one-loop order. This cannot be cured by using a different regularization scheme since the upper critical dimension $\frac 52$ is noninteger. To avoid this obstacle we restrict ourselves to the presence of a single type of disorder. It is easy to see that the diagrams (c) and (d) in fig. 5 with any type of disorder always contribute to α₁ ≡ α_x. Thus we can derive a closed one-loop flow equation by considering a model with only the α_x type of disorder. The corresponding RG equations read

$$\begin{equation} -\mu \partial_\mu \ln \varepsilon = 1+\alpha_x,\quad -\mu \partial_\mu \alpha_x = \delta\, \alpha_x -\alpha_x^2, \end{equation} \tag{ 17 }$$

where we have included the value of the one-loop integral in redefinition of α_x. Equations (17) possess a fixed point α*_x = δ. Integrating them up to the scale corresponding to ε ≃ 1, we find ε = μ^z with z = 1 + α*_x = 1 + δ + O(δ²). The DOS behaves as $\rho \sim \langle \bar {\psi } \psi \rangle$, which leads to its scaling as ρ ∼ μ^3/2−z or equivalently in terms of energy as

$$\begin{equation} \rho(\varepsilon) \sim \varepsilon^{(3/2-z)/z}= \varepsilon^{(1/2-\delta)/(1+\delta)+O(\delta^2)}. \end{equation} \tag{ 18 }$$

Hence we find that for $\delta =\frac 12$ the DOS exponent vanishes to one-loop order suggesting a constant DOS at zero energy when only the α_x disorder is present. This result improves upon the critical behavior extracted from the SCBA analysis which was shown to be invalid at low energy. Whether the exponent vanishing is accidental or holds beyond the one-loop order remains an open question.

Zoom In Zoom Out Reset image size

To analyze the scaling behavior of the full disordered 2D semi-Dirac fermions, we find it more technically tractable to consider an alternative generalized model given by the effective action

$$\begin{eqnarray*} & & \mathcal{S} = \int {\rm d}^d \mathbf{r}\left \{ \sum\limits_{a=1}^n \bar{\psi}_{a}\left[ \sigma_x \left(\Delta - \frac{\nabla_{\|}^2}{2m}\right) - i c \sigma_y \partial_y - i \varepsilon \right]\psi_a \right. \\ & & \left. - \frac{\pi c}{\sqrt{m}} \sum\limits_{a,b=1}^n \sum\limits_{j=0,x,y,z} \alpha_{j} [\bar{\psi}_{a}(\mathbf{r})\sigma_{j} \psi_{a}(\mathbf{r})] [\bar{\psi}_{b}(\mathbf{r})\sigma_{j} \psi_{b}(\mathbf{r})] \vphantom{\sum\limits_{a=1}^n \bar{\psi}_{a}} \right\}, \end{eqnarray*} \tag{ 20 }$$

where d^dr: = d^d−1r_∥ dr_⊥. The advantage of this model is that it involves no infinite Clifford algebra. Dimensional analysis of this model suggests that its upper critical dimension is d_c = 3. We can construct a perturbative in disorder RG around this dimension using expansion in a small parameter δ = (3 − d)/2 similar to that done for the model (16). The corrections to energy and disorder strengths are given by the same diagrams shown in fig. 5. The one-loop RG equations read

$$\begin{eqnarray*} & & -\mu \partial_\mu \ln \varepsilon = 1+\alpha_0+\alpha_x +\alpha_y+\alpha_z, \\ & & -\mu \partial_\mu \alpha_0 = \delta\, \alpha_0 + 4 \alpha_0 (\alpha_0+2\alpha_x+\alpha_y+\alpha_z) + 4 \alpha_x \alpha_z, \nonumber \\ & & -\mu \partial_\mu \alpha_x = \delta\, \alpha_x + 2 \left(\alpha_0^2+\alpha_x^2+\alpha_y^2+\alpha_z^2\right) +4 \alpha_0 \alpha_z, \nonumber\\ & & -\mu \partial_\mu \alpha_y = \delta\, \alpha_y + 4 \alpha_x \alpha_y, \nonumber \\ & & -\mu \partial_\mu \alpha_z = \delta\, \alpha_z -4 \alpha_z (\alpha_0-2\alpha_x-\alpha_y+\alpha_z) + 4 \alpha_0 \alpha_x. \nonumber \end{eqnarray*}$$

Unfortunately, the solutions of these RG equations always flow towards a strong-disorder regime which describes Anderson localization, but is not within the reach of a weak-disorder RG analysis. In this generalized model, we do no find a critical perturbative fixed point, as opposed to the model with only α_x disorder.