Anderson localization of two-dimensional massless pseudospin-1 Dirac particles in a correlated random one-dimensional scalar potential

TL;DR

This paper investigates how correlated disorder affects Anderson localization of two-dimensional massless pseudospin-1 Dirac particles, deriving analytical and numerical results for localization length across various regimes and parameters.

Contribution

It provides the first detailed analysis of disorder correlation effects on localization length for pseudospin-1 Dirac particles, including analytical formulas and scaling regimes.

Findings

Localization length diverges as incident angle approaches zero.

Localization length depends non-monotonically on disorder correlation length.

Multiple scaling regimes are identified depending on disorder strength and energy.

Abstract

We study theoretically Anderson localization of two-dimensional massless pseudospin-1 Dirac particles in a random one-dimensional scalar potential. We focus explicitly on the effect of disorder correlations, considering a short-range correlated dichotomous random potential at all strengths of disorder. We also consider a $\delta$-function correlated random potential at weak disorder. Using the invariant imbedding method, we calculate the localization length in a numerically precise way and analyze its dependencies on incident angle, disorder correlation length, disorder strength, energy, wavelength and average potential over a wide range of parameter values. In addition, we derive analytical formulas for the localization length, which are very accurate in the weak and strong disorder regimes. From the Dirac equation, we obtain an expression for the effective wave impedance, using which we explain several conditions for delocalization. We also deduce a condition under which the localization length vanishes. For all cases considered, the localization length depends non-monotonically on the disorder correlation length and diverges as $\theta^{-4}$ as the incident angle $\theta$ goes to zero. As the disorder strength is varied from zero to infinity, we find that there appear three different scaling regimes. As the energy or wavelength is varied from zero to infinity, there appear three or four different scaling regimes with different exponents, depending on the value of the average potential. The crossovers between different scaling regimes are explained in terms of the disorder correlation effect.

Equations

$$\mathcal{H}=v_F\left(S_x{p}_x+S_y{p}_y\right)+UI,$$

$$S_x=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right),S_y=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0& -i& 0\\ i& 0& -i\\ 0& i& 0\end{array}\right).$$

$$p_x=\frac{\hslash }{i}\frac{d}{dx},p_y=\hslash k_y,$$

$$\mathcal{H}\psi =E\psi ,$$

$${\psi }_1$$

$${\psi }_3$$

$$\frac{d}{dx}\left(\frac{\hslash v_F}{E-U}\frac{d{\psi }_2}{dx}\right)+\left(\frac{E-U}{\hslash v_F}-\frac{\hslash v_F{k_y}^2}{E-U}\right){\psi }_2=0.$$

$$k=\frac{E}{\hslash v_F}.$$

$$k_x=k\cos \theta ,k_y=k\sin \theta ,$$

$$ϵ=1-u,$$

$$\frac{d}{dx}\left(\frac{1}{ϵ}\frac{d{\psi }_2}{dx}\right)+{k_x}^2ϵ{\eta }^2{\psi }_2=0,$$

$${\eta }^2$$

$$U=U_0+\delta U(x),$$

$$⟨\delta U(x)\delta U(x^{\prime })⟩=G\delta (x-x^{\prime }),⟨\delta U(x)⟩=0,$$

$$⟨\delta U(x)\delta U(x^{\prime })⟩=S^2\exp \left(-|x-x^{\prime }|/l_c\right),⟨\delta U(x)⟩=0.$$

$$\displaystyle\psi_{2}\left(x,L\right)=\left\{\begin{array}[]{ll}e^{-ik_{x}\left(x-L\right)}+r(L)e^{ik_{x}\left(x-L\right)},&x>L\\ t(L)e^{-ik_{x}x},&x<0\end{array},\right.$$

$$\frac{d}{dx}\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right)=A\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right),$$

$$f_1={\psi }_2,f_2=\frac{\hslash v_F}{E-U}\frac{d{\psi }_2}{dx},$$

$$A=\left(\begin{array}{cc}0& \frac{E-U}{\hslash v_F}\\ \frac{\hslash v_F{k_y}^2}{E-U}-\frac{E-U}{\hslash v_F}& 0\end{array}\right).$$

$$\frac{1}{k}\frac{dr}{dl}=-\frac{i\cos \theta }{2}ϵ{\left(r-1\right)}^2$$

$$+\frac{i}{2\cos \theta }\left(ϵ-\frac{{\sin }^2\theta }{ϵ}\right){\left(r+1\right)}^2,$$

$$\frac{1}{k}\frac{dt}{dl}=-\frac{i\cos \theta }{2}ϵ\left(r-1\right)t$$

$$+\frac{i}{2\cos \theta }\left(ϵ-\frac{{\sin }^2\theta }{ϵ}\right)\left(r+1\right)t,$$

$$\xi =-\underset{L\to \infty }{\lim }\left(\frac{L}{⟨\ln T⟩}\right).$$

$$\displaystyle\bigg{|}\frac{\delta u}{\epsilon_{0}}\bigg{|}\ll 1,$$

$${ϵ}_0=1-u_0,u_0=\frac{U_0}{E},\delta u=\frac{\delta U}{E}.$$

$$\frac{1}{ϵ}=\frac{1}{{ϵ}_0-\delta u}\approx \frac{1}{{ϵ}_0}\left(1+\frac{\delta u}{{ϵ}_0}\right).$$

$$-\frac{1}{k}\frac{d⟨\ln T⟩}{dl}=\frac{g{D_1}^2}{{\cos }^2\theta }+\frac{{ϵ}_0{C}_1}{\cos \theta }\mathrm{Im}(Z_1)$$

$$+\frac{2g{D}_0{D}_1}{{\cos }^2\theta }\mathrm{Re}(Z_1)+\frac{g{D_1}^2}{{\cos }^2\theta }\mathrm{Re}(Z_2),$$

$$Z_n=⟨r^n⟩,g=\frac{G}{4\hslash v_{F}E},C_1=\left(1-\frac{1}{{{ϵ}_0}^2}\right){\sin }^2\theta ,$$

$$D_0=1+{\cos }^2\theta +\frac{{\sin }^2\theta }{{{ϵ}_0}^2},D_1=\left(1+\frac{1}{{{ϵ}_0}^2}\right){\sin }^2\theta .$$

$$\frac{1}{k\xi }=\frac{g{D_1}^2}{{\cos }^2\theta }+\frac{{ϵ}_0{C}_1}{\cos \theta }\mathrm{Im}[Z_1(l\to \infty )]$$

$$+\frac{2g{D}_0{D}_1}{{\cos }^2\theta }\mathrm{Re}[Z_1(l\to \infty )]$$

$$+\frac{g{D_1}^2}{{\cos }^2\theta }\mathrm{Re}[Z_2(l\to \infty )].$$

Full text

I introduction

In Dirac materials, the quasiparticles obey an effective Dirac-type equation and their dynamics resembles that of relativistic particles. Since the experimental isolation of single-layer graphene in 2004, the interest in Dirac materials has increased explosively motivated by their promising potential applications in various devices and many interesting physical properties [1, 2]. In addition to single-layer and bilayer graphene, the study of Dirac materials has expanded to other pseudo-relativistic materials such as topological insulators, pseudospin-$N$ Dirac materials, Weyl and Dirac semimetals and Kane semiconductors and the list keeps growing [3, 4, 5, 6, 7, 8, 9]. Since Dirac materials and metamaterials can be realized in condensed-matter systems, photonic-crystal structures and cold-atom optical lattices, researches in pseudo-relativistic systems have become very popular in many disciplines of physics including condensed matter physics, optics and atomic physics [10, 11, 12, 13, 14].

If we ignore effects such as intervalley scattering, electron spin and electron-electron and electron-phonon interactions, the quasiparticles in graphene are described by a two-dimensional (2D) pseudospin-1/2 Dirac equation for massless particles, where the pseudospin represents the two different sublattices of the underlying honeycomb lattice [3]. More recently, the study of pseudospin-1/2 systems has been extended to those with general pseudospin $N$ ($=1$, 3/2, 2, $\cdots$) [15, 16]. In the massless case, all these pseudospin-$N$ Dirac systems display the Klein tunneling effect, which is manifested as a total transmission of normally-incident particles through an arbitrary scalar potential barrier [17, 18, 19]. In the case of pseudospin-1 systems, an omnidirectional total transmission phenomenon termed super-Klein tunneling, which occurs when the particle energy is precisely one half of the constant potential barrier, has attracted much interest [20, 21, 22, 23, 24]. In addition to these, there have been many recent theoretical studies exploring various properties of pseudospin-1 systems [25, 26].

The unique transport properties of Dirac materials also influence the nature of Anderson localization in random potentials [27, 28, 29]. Differences in the types of wave equations, the material properties and the nature of disorder can influence Anderson localization strongly [30, 31, 32, 33, 34, 35, 36, 37]. Since the wave equation obtained from the Dirac equation is of a substantially different type from the Schrödinger equation and the electromagnetic wave equation, we expect conceptually new localization phenomena to arise in Dirac materials. Anderson localization in pseudospin-1/2 systems in a random one-dimensional (1D) scalar potential has been studied theoretically by several authors [38, 39, 40, 41, 42, 43]. It has been shown that localization is destroyed at normal incidence due to the Klein tunneling effect. Close to the normal incidence, the localization length $\xi$ has been shown to scale as $\xi\propto\theta^{-2}$, where $\theta$ is the incident angle. As the disorder strength increases from zero to infinity, the localization length shows a non-monotonic behavior, such that it initially decreases, attains a minimum value and then increases to infinity. This type of counterintuitive delocalization effect induced by strong disorder has been interpreted in terms of effective impedance [41].

Anderson localization of pseudospin-1 Dirac particles in a random 1D scalar potential was studied theoretically in two recent papers [42, 43]. In Ref. [42], the localization length for a random multilayer structure, where all layers had the same thickness and the potential in each layer was a random variable distributed uniformly in the range $[-W,W]$, was calculated numerically using the transfer matrix method by averaging over 4000 random configurations. In addition, some analytical results on the localization length were obtained using the surface Green function method. The authors reported a very peculiar transition behavior such that the localization length diverged as $\xi\propto\sin^{-4}\theta$ at normal incidence when $W$ was smaller than the particle energy $E$, while it diverged as $\xi\propto\sin^{-2}\theta$ when $W$ was larger than $E$.

In Ref. [43], the localization length was calculated as a function of the wavelength in the long wavelength limit using the same method as in Ref. [42] for three different types of random multilayer models. In the first model, the potential in each layer was randomly selected from the two values $V_{0}+W$ and $V_{0}-W$. In this binary case, the authors showed that the localization length scaled as $\xi\propto\lambda^{6}$ if $V_{0}=0$ and as $\xi\propto\lambda^{4}$ if $V_{0}\neq 0$ in the long wavelength limit. When $V_{0}\neq 0$, they also reported that there appeared a sharp peak of the localization length at the value of $\lambda$ corresponding to $E=V_{0}$, signifying the onset of delocalization. In the second model, the potential in the $i$-th layer was randomly selected from the two values $W(1+\delta_{i})$ and $-W(1+\delta_{i})$, where $\delta_{i}$ is a random number distributed uniformly in a small interval $[-Q,Q]$. In this case, the authors showed that $\xi$ scaled as $\xi\propto\lambda^{4}$ in the long wavelength limit. In the third model, layers with zero potential were alternated periodically with those with a random potential and the scaling of $\xi$ was found to be $\xi\propto\lambda^{2}$ in the long wavelength limit. It was argued that this nonuniversal dependence of the scaling exponent on the specific type of the random potential was a characteristic of pseudospin-1 Dirac systems.

We point out that the random multilayer structures considered in Refs. [42] and [43] have short-range disorder correlations imbedded in them. However, the influence of these correlations on localization phenomena was not studied explicitly. In Ref. [42], the incident angle dependence was investigated carefully, but the scaling dependence on other parameters was not analyzed in detail especially in the strong disorder regime. On the other hand, the results reported in Ref. [43] were limited to the long wavelength limit. Therefore it is highly desirable to study the scaling behavior of a single model in the entire wavelength (or energy) range from zero to infinity and the effect of disorder correlations on the the crossover between different scaling behaviors.

Recently, there has been much interest in the roles of short-range and long-range disorder correlations in localization phenomena [33, 44, 45, 46]. In this paper, we extend this approach to the localization in Dirac systems. We consider the continuum Dirac equation in 2D for massless pseudospin-1 particles in a random 1D scalar potential and calculate the localization length using the invariant imbedding method (IIM) [47, 48, 49, 50, 51, 52, 53, 54]. In order to study the scaling behavior of the localization length in the entire energy range and the influence of disorder correlations on it, we consider a short-range correlated dichotomous (or binary) random potential which is explicitly characterized by a correlation length $l_{c}$ as well as its strength. In addition, we consider a $\delta$-correlated (or uncorrelated) random potential in the weak disorder regime and compare the results with those from the short-range correlated model. One of the main advantages of our IIM is that we can perform the disorder averaging analytically in an exact way and convert the original random problem into an equivalent nonrandom one. Therefore we can avoid repeating calculations for a large number of random configurations and then averaging over the results. Another advantage is that it is possible to derive analytical formulas for the localization length which are extremely accurate in the weak and strong disorder regimes.

Using our method, we investigate the dependencies of the localization length on incident angle, disorder correlation length, disorder strength, energy, wavelength and average potential over a wide range of parameter values. We find different scaling behaviors in the different regions of the parameter space and explain the crossovers between them in terms of the disorder correlation effect. We also derive an expression for the effective wave impedance, using which we interpret several delocalization phenomena.

Massless pseudospin-1 systems are characterized by two Dirac cones intersected by a flat band at the Dirac point. It has been known that they can be realized in 2D lattices such as the Lieb, dice and kagome lattices [55]. Recently, there has been much progress in the experimental realization of pseudospin-1 systems. Several candidate 2D materials, which may be described by the pseudospin-1 Dirac equation, have been discovered [56, 57]. There also have been many attempts to construct artificial 2D lattice structures displaying the properties of pseudospin-1 systems, based on cold-atom optical lattices, photonic crystal structures and artificial electronic lattices [58, 59, 60]. Due to this rapid development, we expect that it will be possible to explore the unique properties of these systems experimentally in the near future.

The rest of this paper is organized as follows. In Sec. II, we introduce the pseudospin-1 Dirac equation and the two different kinds of random potentials used in this study. We also derive an analytical expression for the wave impedance. In Sec. III, the IIM for the calculation of the localization length is described and the invariant imbedding equations for the two random models are derived. In Sec. IV, we apply the perturbation expansion method to the invariant imbedding equations and derive analytical formulas for the localization length in the weak and strong disorder regimes. In Sec. V, we present detailed numerical results obtained using the IIM and discuss the dependencies of the localization length on incident angle, disorder correlation length, disorder strength, energy, wavelength and average potential. We conclude the paper in Sec. VI.

II Model

The effective Hamiltonian that describes massless pseudospin-1 Dirac-like particles moving in the 2D $xy$ plane in a 1D scalar potential $U=U(x)$ takes the form

[TABLE]

where $v_{F}$ is the Fermi velocity and $I$ is the $3\times 3$ unity matrix. The $x$ and $y$ components of the pseudospin-1 operator, $S_{x}$ and $S_{y}$, are represented by

[TABLE]

In the case where the potential $U$ depends only on $x$, the $x$ and $y$ components of the momentum operator, $p_{x}$ and $p_{y}$, are given by

[TABLE]

where $k_{y}$ is the constant of motion.

The time-independent Dirac equation in 2D for the three-component vector wave function $\psi$ [$=\left(\psi_{1},\psi_{2},\psi_{3}\right)^{\rm T}$] is

[TABLE]

which is a set of three coupled first-order differential equations for $\psi_{1}$, $\psi_{2}$ and $\psi_{3}$. We can eliminate $\psi_{1}$ and $\psi_{3}$ using the equations

[TABLE]

and obtain a single wave equation for $\psi_{2}$ of the form

[TABLE]

We assume that a plane wave described by $\psi_{2}$ is incident obliquely from the free region $x>L$ where $U=0$ onto the nonuniform region in $0\leq x\leq L$ where $U=U(x)$, and then transmitted to the free region $x<0$ where $U=0$. The wave number in the free regions, $k$, is related to the particle energy $E$ by

[TABLE]

The negative $x$ component of the wave vector in the incident and transmitted regions, $k_{x}$, and the $y$ component of the wave vector, $k_{y}$, are given by

[TABLE]

where $\theta$ is the incident angle.

We introduce a dimensionless quantity $\epsilon=\epsilon(x)$ defined by

[TABLE]

where $u=U/E$. In the incident and transmitted regions where $U=0$, we have $\epsilon=1$. In terms of $\epsilon$, the wave equation, Eq. (6), can be written concisely as

[TABLE]

where $\eta$ is defined by

[TABLE]

We notice that the wave equation of this form looks identical to that for $p$-polarized electromagnetic waves propagating normally in a medium with the wave impedance function given by $\eta(x)$. In the free regions where $\epsilon=1$, $\eta^{2}$ is identically equal to 1. Therefore, if $\eta^{2}$ is unity in the inhomogeneous region $0\leq x\leq L$ as well, the impedance is matched throughout the space and there will be no wave scattering.

There are two cases where the impedance matching can be achieved. If the incident angle $\theta$ is zero, $\eta^{2}$ is identically equal to 1 regardless of the magnitude and the functional form of $U(x)$. This is the well-known Klein tunneling phenomenon occurring at normal incidence. The other case occurs when $\epsilon=-1$ corresponding to a constant potential barrier of height $U=2E$. In this case, $\eta^{2}$ is unity for all $\theta$ and therefore an omnidirectional total transmission arises. This phenomenon has been termed super-Klein tunneling [20, 21, 22, 23, 24].

In this paper, we are interested in the localization of pseudospin-1 Dirac-like particles in a random 1D scalar potential. We assume that in the region $0\leq x\leq L$, $U$ is a random function of $x$ given by

[TABLE]

where $U_{0}$ is the disorder-averaged value of $U$ and $\delta U(x)$ is the randomly-fluctuating part of $U$ with zero mean. We will consider two different types of $\delta U(x)$. In the first case, it is a Gaussian random function satisfying

[TABLE]

where the notation $\langle\cdots\rangle$ denotes averaging over disorder and $G$ is a parameter characterizing the strength of disorder. This case corresponds to uncorrelated white noise. In the second case, we assume that $\delta U(x)$ is a short-range correlated dichotomous (or binary) Gaussian random function, which fluctuates randomly between the two values $S$ and $-S$ and satisfies

[TABLE]

In this case, $S^{2}$ measures the strength of disorder and the parameter $l_{c}$ is the correlation length of disorder. By short-range correlation, we mean that the correlation function $\langle\delta U(x)\delta U(x^{\prime})\rangle$ decays exponentially versus $|x-x^{\prime}|$. If it decays as a power law, the random potential is called long-range correlated.

Our method is based on deriving a set of equivalent nonrandom differential equations starting from the original random wave equation, where the effect of randomness is taken care of analytically. Later, we will show that in the case of $\delta$-correlated randomness, this method can be applied only when the disorder is sufficiently weak, while in the case of dichotomous randomness, it can be applied to arbitrary strengths of disorder. Therefore our numerical calculations will be mainly focused on the case of short-range correlated dichotomous randomness.

III Invariant imbedding method

We can solve the wave equation, Eq. (6), using the IIM. In this method, we first calculate the reflection and transmission coefficients $r$ and $t$ defined by the wave functions in the incident and transmitted regions:

[TABLE]

where $r$ and $t$ are regarded as functions of $L$. The first step in the derivation of the invariant imbedding equations for $r$ and $t$ is to rewrite Eq. (6) in the form

[TABLE]

where $A$ ia a $2\times 2$ matrix. With the definitions

[TABLE]

we can easily obtain

[TABLE]

Then we follow the procedure described in Ref. [54] and derive

[TABLE]

where $r$, $t$ and $\epsilon$ are functions of $l$. Here, the variable $l$ represents the thickness of the system in the direction of the inhomogeneity and is called the imbedding parameter. The quantity $r(l)$ denotes the reflection coefficient of a hypothetical system of thickness $l$ and $\epsilon(l)$ has the same functional form as $\epsilon(x)$.

For any functional form of $U$ and for any values of $kL$ and $\theta$, we can integrate these equations from $l=0$ to $l=L$ using the initial conditions $r(0)=0$ and $t(0)=1$ and obtain $r(L)$ and $t(L)$. The reflectance $R$ and the transmittance $T$ are obtained using $R=|r|^{2}$ and $T=|t|^{2}$. In the absence of dissipation, the law of energy conservation requires that $R+T=1$. In this paper, we are mainly interested in calculating the localization length $\xi$ defined by

[TABLE]

The differential equation satisfied by $\ln T$ can be easily derived from the second of Eq. (22).

The parameter $\epsilon$ in Eq. (22) is a random function of $l$. Therefore these equations are stochastic differential equations with random coefficients. The nonrandom differential equations satisfied by the disorder-averaged quantities such as $\langle R\rangle$, $\langle T\rangle$ and $\langle\ln T\rangle$ can be derived from Eq. (22) using standard methods of stochastic differential equations. We will use Furutsu-Novikov formula [61, 62] in the case of $\delta$-correlated randomness and the formula of differentiation of Shapiro and Loginov [63] in the case of short-range correlated dichotomous randomness.

III.1 $\delta$-correlated random potential

In the case of a $\delta$-correlated random potential, we have difficulty in handling the random function $\epsilon$ appearing in the denominators of some coefficients in Eq. (22). It is possible to raise it to the numerator if the randomness is sufficiently weak such that

[TABLE]

where

[TABLE]

Then we can use the approximation

[TABLE]

It is now straightforward to derive the equation for $\langle\ln T\rangle$ using Furutsu-Novikov formula, which takes the form

[TABLE]

where $Z_{n}$ ($n=1,2,\cdots$) and the parameters $g$, $C_{1}$, $D_{0}$ and $D_{1}$ are defined by

[TABLE]

From Eqs. (23) and (27), we find that the localization length $\xi$ is given by

[TABLE]

Using the first of Eq. (22) and Furutsu-Novikov formula, we can also derive an infinite number of coupled differential equations for $Z_{n}$:

[TABLE]

where $C_{0}$ is defined by

[TABLE]

These equations are supplemented with the initial conditions $Z_{0}=1$ and $Z_{n}=0$ for $n\geq 1$. In the large $l$ limit corresponding to the case where the thickness of the random region diverges, all $Z_{n}$’s should approach constants independent of $l$, if the disorder-averaged value of $\epsilon$ is a constant independent of $x$. Then we can set the left-hand side of Eq. (30) to zero and obtain an infinite number of coupled algebraic equations. We solve these equations numerically by a systematic truncation method [48] and obtain $Z_{1}(l\to\infty)$ and $Z_{2}(l\to\infty)$, using which we calculate the localization length.

III.2 Short-range correlated dichotomous random potential

In the case of a short-range correlated dichotomous random potential, the parameter $\epsilon$ takes only two values, $\epsilon_{0}+\sigma$ and $\epsilon_{0}-\sigma$, where $\sigma$ is defined by $\sigma=S/E$. Then $1/\epsilon$ fluctuates between $1/(\epsilon_{0}+\sigma)$ and $1/(\epsilon_{0}-\sigma)$, the average of which is $\epsilon_{0}/({\epsilon_{0}}^{2}-\sigma^{2})$. From this, we obtain the identity

[TABLE]

In contrast to the $\delta$-correlated case, there is no approximation involved here. The nonrandom differential equation satisfied by $\left\langle\ln{T}\right\rangle$ takes the form

[TABLE]

where $W_{n}$ ($n=1,2,\cdots$), $\tilde{C}_{1}$ and $\tilde{D}_{1}$ are defined by

[TABLE]

The localization length is given by

[TABLE]

In order to obtain $Z_{1}$ and $W_{1}$ in the large $l$ limit, we need to derive an infinite number of coupled nonrandom differential equations satisfied by $Z_{n}$’s and $W_{n}$’s using the formula of differentiation of Shapiro and Loginov, which takes the form

[TABLE]

where $j$ is an arbitrary nonnegative integer and the function $f$ satisfies an ordinary differential equation with random coefficients expressed in terms of $\zeta$ [63]. The dichotomous Gaussian random function $\zeta$ satisfies the condition

[TABLE]

By substituting $f$ and $\zeta$ with $r^{n}$ and $\delta u$ respectively and taking $j=0,1$, we obtain

[TABLE]

where

[TABLE]

They are supplemented with the initial conditions $Z_{0}=1$, $Z_{n}=0$ for $n\geq 1$ and $W_{n}=0$ for all $n$. In the large $l$ limit, all $Z_{n}$’s and $W_{n}$’s are constants independent of $l$. Then we can set the left-hand sides of Eq. (38) to zero and obtain an infinite number of coupled algebraic equations. We solve these equations numerically by a systematic truncation method and obtain $Z_{1}(l\to\infty)$ and $W_{1}(l\to\infty)$, using which we calculate the localization length.

IV Analytical expressions for the localization length in the weak and strong disorder regimes

Though we can solve the invariant imbedding equations, Eqs. (27), (30), (33) and (38), numerically for general parameter values, it is highly instructive to apply the perturbation theory to them to derive analytical expressions for the localization length in some limiting cases.

IV.1 Weak disorder regime in a $\delta$-correlated random potential

In the case of a $\delta$-correlated Gaussian random potential, the imbedding equations, Eqs. (27) and (30), have been derived assuming that the disorder is sufficiently weak. Therefore we apply the perturbation theory to those equations only in the weak disorder regime expressed by Eq. (24). We write the reflection coefficient in the large $l$ limit, $r$, as $r=r_{0}+\delta r$, where $r_{0}$ is the reflection coefficient from an interface between free space and a half-space nonrandom medium with the parameter $\epsilon_{0}$. The expression for $r_{0}$ takes the form

[TABLE]

where $\tilde{p}$ is defined by

[TABLE]

We expand $Z_{n}$ as

[TABLE]

Using this and the infinite number of algebraic equations obtained from Eq. (30), we get an infinite number of coupled equations for $\langle\left(\delta r\right)^{n}\rangle$ for all integers $n$. We expand these averages in terms of the small perturbation parameter $g$. From analytical considerations and numerical calculations, we can show that the leading terms for $\langle\delta r\rangle$ and $\langle\left(\delta r\right)^{2}\rangle$ are of the first order in $g$, whereas that for $\langle\left(\delta r\right)^{3}\rangle$ is of the second order in $g$, except at incident angles close to the critical angle of total reflection $\theta_{c}$ satisfying ${\epsilon_{0}}^{2}=\sin^{2}\theta_{c}$. Therefore, we substitute

[TABLE]

into Eq. (30) when $n=1,2$ in the large $l$ limit and obtain two coupled equations for $\langle\delta r\rangle$ and $\langle\left(\delta r\right)^{2}\rangle$. We solve these equations analytically and substitute the resulting expressions into Eq. (29) to the leading order in $g$. From this, we obtain the localization length of the form

[TABLE]

where $\Theta$ is the step function, $\Theta(x)=1$ for $x>0$ and 0 for $x<0$. We remind that the weak disorder expansion cannot be applied if $\epsilon_{0}\approx 0$ (that is, $E\approx U_{0}$) or ${\epsilon_{0}}^{2}\approx\sin^{2}\theta$.

From Eq. (44), it follows that there is a symmetry under the sign change of $\epsilon_{0}$ and $\theta$. We also notice that in the total reflection (or tunneling) regime where $|\sin\theta|>|\epsilon_{0}|$ and when the disorder parameter $g$ is sufficiently small, $\xi$ increases as $g$ increases from zero. This is an example of the well-known disorder-enhanced tunneling phenomenon [41, 64, 65, 66, 67, 68]. On the other hand, if $|\epsilon_{0}|$ is greater than $|\sin\theta|$, we have $\xi\propto g^{-1}$, which simply means that localization is enhanced by weak disorder.

When $\theta$ is very close to zero, we find that $k\xi\approx[{\epsilon_{0}}^{4}/(4g)]\theta^{-4}$, which diverges at $\theta=0$ due to the Klein tunneling effect. This $\theta^{-4}$ dependence of the localization length near $\theta=0$ is a unique characteristic of pseudospin-1 systems and is distinct from the $\theta^{-2}$ dependence occurring in pseudospin-1/2 systems. If the average potential $U_{0}$ is zero, we obtain

[TABLE]

This result can be compared with the dependence $\xi\propto\sin^{-4}\theta$ reported in Ref. [42], which has been derived using the transfer matrix method and the surface Green function method for a random multilayer model, when the average potential is zero and the disorder strength is smaller than a certain critical value. In that model, all layers have the same thickness and the potential in each layer is a random variable distributed uniformly in the range $[-W,W]$. In Ref. [42], it has also been reported that $\xi\propto\sin^{-2}\theta$ if the disorder strength $W$ is larger than the critical value given by $W_{c}=E$. Our method for the $\delta$-correlated random potential cannot be applied to the strong disorder regime. However, we can solve exactly the case of a short-range correlated dichotomous random potential for arbitrary strengths of disorder. In Secs. IV.2 and IV.3, we will show that the $\theta^{-4}$ dependence near $\theta=0$ is valid in that model regardless of the strength of disorder and the size of the correlation length.

Next, we consider the dependence of the localization length on the energy of the incident particle, or equivalently, on the frequency and the wavelength of the incident wave. For that purpose, it is more convenient to normalize $\xi$ by the wave number associated with the random potential, $k_{u}$, defined by

[TABLE]

which is independent of energy. For simplicity, let us consider only the case where the average potential $U_{0}$ is zero. Then it is trivial to show that

[TABLE]

which is independent of energy (and, equivalently, of frequency and wavelength). This result is applied in the weak disorder limit, which corresponds to the large energy (or frequency) or short wavelength limit. However, we will show in the next subsection that the above result is an artifact of the $\delta$-correlated random model, the spectral density of which has no ultraviolet cutoff, and more realistic short-range correlated random models exhibit different behavior in the limit of asymptotically large energy.

IV.2 Weak disorder regime in a short-range correlated dichotomous random potential

In the case of a short-range correlated dichotomous random potential, the wave number in the medium with dichotomous randomness where $\epsilon=\epsilon_{0}\pm\sigma$ is given by $\epsilon_{0}\left(1\pm\sigma/\epsilon_{0}\right)k$. Therefore it is natural to represent the strength of disorder by the parameter $s$ defined by

[TABLE]

where $\epsilon_{0}\neq 0$. In the weak disorder regime where $s$ is sufficiently small, we can rewrite $\tilde{C}_{0}$, $\tilde{C}_{1}$ , $\tilde{D}_{0}$ and $\tilde{D}_{1}$ in Eqs. (34) and (39) as

[TABLE]

Similarly to Sec. IV.1, we write $r=r_{0}+\delta r$, where $r_{0}$ is given by Eq. (40), and expand $Z_{n}$ as in Eq. (42) and $W_{n}$ as

[TABLE]

From analytical considerations and numerical calculations, we can show that both $|\langle(\delta r)^{j}\rangle|$ and $|\langle(\delta r)^{j}\delta u\rangle|$ decrease rapidly as $j$ increases. We can also demonstrate that the leading terms for $\langle\delta r\rangle$, $\langle(\delta r)^{2}\rangle$ and $\langle\delta r\delta u\rangle$ are of the first order in $s$, while those for $\langle(\delta r)^{3}\rangle$ and $\langle(\delta r)^{2}\delta u\rangle$ are of the second order in $s$. From these considerations, we derive the analytical expression for the localization length in the weak disorder regime of the form

[TABLE]

From this, we find that there is a symmetry under the sign change of $\epsilon_{0}$ and $\theta$ and the localization length diverges as $\theta^{-4}$ at $\theta=0$. In the total reflection regime where $|\sin\theta|>|\epsilon_{0}|$, the localization length increases as the disorder strength $\sigma$ increases, only if

[TABLE]

Therefore, in contrast to the case of $\delta$-correlated randomness, the disorder-enhanced tunneling phenomenon in the present case occurs only when the correlation length is sufficiently large and $|\epsilon_{0}|$ is not too small.

We also observe that when $|\epsilon_{0}|$ is larger than $|\sin\theta|$, the expression given by Eq. (51b) reduces to Eq. (44) in the limit where $kl_{c}\to 0$ and $kl_{c}\sigma^{2}\to 2g$. Therefore, except in the total reflection regime, the case of $\delta$-correlated randomness can be considered as that of short-range correlated randomness in the limit where the correlation length vanishes.

When $|\epsilon_{0}|$ is larger than $|\sin\theta|$, we have the dependence $\xi\propto\sigma^{-2}$ regardless of the values of other parameters including the correlation length. In addition, we have a non-monotonic dependence of $\xi$ on the correlation length $l_{c}$. When $l_{c}$ is sufficiently small, $\xi$ is inversely proportional to $l_{c}$, whereas when $l_{c}$ is sufficiently large, it is proportional to $l_{c}$, as we can see clearly from the approximate expressions derived from Eq. (51b):

[TABLE]

where $l_{c}^{\rm min}$ is the value of $l_{c}$ at which the localization length takes a minimum value and is given by

[TABLE]

When the average potential is zero, these expressions are simplified to

[TABLE]

Since $k\sqrt{{\epsilon_{0}}^{2}-\sin^{2}\theta}$ is the $x$ component of the wave vector in the random region in an averaged sense, we can rewrite Eq. (54) as

[TABLE]

where we have defined an effective wavelength $\lambda_{\rm eff}$ by

[TABLE]

Therefore the two different regimes described by Eqs. (53a) and (53b) are distinguished by the relative size of the correlation length and the effective wavelength. In the regime where the effective wavelength is much larger than the correlation length, the model is effectively uncorrelated and the correlation effect can be ignored. In the opposite regime where the effective wavelength is much smaller than the correlation length, the correlation effect should be important.

We next consider the dependence of the localization length on the energy and the wavelength. In the present model, we normalize $\xi$ by the wave number associated with the random potential, $k_{d}$, defined by

[TABLE]

which is independent of energy. For simplicity, we consider only the case where the average potential is zero such that $\epsilon_{0}=1$. Then we obtain

[TABLE]

which can be approximated in two forms

[TABLE]

depending on the relative size of $k/k_{d}$ ($=E/S$). This result is applied in the weak disorder limit, corresponding to the large energy (or frequency) or short wavelength limit. We notice that Eq. (60a), which is independent of energy, is the same as the corresponding expression in the $\delta$-correlated case, Eq. (47), since the parameter $k_{u}$ is equal to $2{k_{d}}^{2}l_{c}$ if we identify $kl_{c}\sigma^{2}$ with $2g$. In the asymptotically large energy limit, however, we have to apply Eq. (60b), which shows that $\xi$ has the dependence $\xi\propto k^{2}$, $E^{2}$, $\omega^{2}$, $\lambda^{-2}$, where $\omega$ is the frequency and $\lambda$ is the wavelength. It is straightforward to verify that similar scaling behaviors are obtained when the average potential is nonzero.

The boundary between the two scaling regions is given by the condition $kl_{c}\cos\theta=1/2$. We find that in the region described by Eq. (60a), the disorder correlation effect is unimportant, while in the asymptotically large energy region described by Eq. (60b), correlations play a crucial role. In a previous paper on the localization of pseudospin-1/2 Dirac electrons in a $\delta$-correlated random scalar potential, we have reported that the localization length is independent of energy in the large energy limit. This result was due to the use of a $\delta$-correlated random potential and the true behavior in the asymptotically large energy limit for more realistic short-range correlated random models has to be the same $\xi\propto E^{2}$ behavior as obtained here. If the correlation length is very small, however, a wide regime where $\xi$ is independent of energy and wavelength should precede the asymptotic regime.

IV.3 Strong disorder regime in a short-range correlated dichotomous random potential

In the strong disorder regime where $s\gg 1$, we can rewrite $\tilde{C}_{0}$, $\tilde{C}_{1}$ , $\tilde{D}_{0}$ and $\tilde{D}_{1}$ in Eqs. (34) and (39) as

[TABLE]

We note that in this perturbation regime, the parameter $\epsilon_{0}\approx 0$ is allowed. We write $r=r_{\infty}+\delta r$, where $r_{\infty}$ is the reflection coefficient from an interface between free space and an infinitely-disordered half-space medium with the parameter $\epsilon_{0}$ and is given by

[TABLE]

Similarly to Sec. IV.2, we expand $Z_{n}$ and $W_{n}$ as

[TABLE]

We can show that both $|\langle(\delta r)^{j}\rangle|$ and $|\langle(\delta r)^{j}\delta u\rangle|$ decrease rapidly as $j$ increases. We can also demonstrate the scaling relationships

[TABLE]

From these considerations, we can derive an analytical expression of the localization length in the strong disorder regime of the form

[TABLE]

Again, there is a symmetry under the sign change of $\epsilon_{0}$ and $\theta$. The localization length depends on $\theta$ as $\xi\propto\sin^{-4}\theta$ for all parameter values. It also diverges as $\epsilon_{0}\to 0$, $l_{c}\to\infty$ or $\sigma\to\infty$. The divergence of $\xi$ in the strong disorder limit, which also occurs in the pseudospin-1/2 case, is counterintuitive and can be understood from the form of the wave impedance, Eq. (11). We find that in the strong disorder limit where $\delta u$ is statistically much larger than 1, the impedance approaches to a constant given by $\eta\approx|\cos\theta|^{-1}$. Therefore, as the disorder strength approaches to infinity, the system becomes effectively nonrandom and the localization length should diverge for all $\theta$.

We can also understand the delocalization occurring when $\epsilon_{0}$ vanishes from the form of the wave impedance. In that case, $\epsilon^{2}$ is always equal to $\sigma^{2}$ and nonrandom. Therefore the impedance is nonrandom for any incident angle and complete delocalization arises. However, we stress that this type of delocalization is strictly limited to the dichotomous randomness and will not arise in more generic models.

It is easy to understand the delocalization occurring when $l_{c}\to\infty$. In a short-range correlated model, the potential can be considered nonrandom within the correlation length $l_{c}$. Therefore in the limit where $l_{c}$ diverges, the potential is effectively nonrandom and localization is destroyed.

Similarly to the weak disorder case, we have a non-monotonic dependence of $\xi$ on the correlation length $l_{c}$, as we can see from the approximate expressions derived from Eq. (65):

[TABLE]

where $\tilde{l}_{c}^{\rm min}$ is the value of $l_{c}$ at which the localization length takes a minimum value and is given by

[TABLE]

This condition can also be interpreted in terms of the effective wavelength. If we take the absolute value of the geometric average of the wave numbers in the dichotomous random potential, $(\epsilon_{0}+\sigma)k$ and $(\epsilon_{0}-\sigma)k$, as the effective wave number $k_{\rm eff}$, we obtain $k_{\rm eff}=k|\sigma^{2}-{\epsilon_{0}}^{2}|/{|\epsilon_{0}|}$. From this, we can define the effective wavelength in the strong disorder regime as $\tilde{\lambda}_{\rm eff}=2\pi|\epsilon_{0}|/(k\sigma^{2})$. Then Eq. (67) can be written similarly to Eq. (56).

We notice that in the regime where $l_{c}$ is sufficiently small, $\xi$ is proportional to $\sigma^{2}$, while in the regime where $l_{c}$ is sufficiently large, $\xi$ is proportional to $\sigma^{6}$. The $\xi\propto\sigma^{6}$ scaling behavior and the crossover between the two scaling regimes have never been obtained before. We stress that these contrasting scaling behaviors are universal and occur regardless of the value of the average potential $U_{0}$, unless $U_{0}=E$. We also point out that $\tilde{l}_{c}^{\rm min}$ is inversely proportional to $\sigma^{2}$, and therefore in the strong disorder regime, the minimum of the localization length occurs at a rapidly-decreasing value of $kl_{c}$ as $\sigma$ increases.

We next consider the dependence of the localization length on the energy and the wavelength. We normalize $\xi$ by the wave number associated $k_{d}$ defined by Eq. (58) and rewrite Eq. (65) as

[TABLE]

where $\tilde{u}_{0}$ is the normalized average potential defined by $\tilde{u}_{0}=U_{0}/S$. In the strong disorder regime, it is necessary to distinguish the cases where the average potential is zero or nonzero. We first consider the case where $\tilde{u}_{0}=0$. Then $k_{d}\xi$ can be approximated in two forms

[TABLE]

depending on the relative size of $k/k_{d}$. This result is applied in the strong disorder limit, corresponding to the small energy (or frequency) or long wavelength limit. The boundary between the two scaling regions is given by the condition $k=2{k_{d}}^{2}l_{c}$, which is equivalent to $\tilde{\lambda}_{\rm eff}=4\pi l_{c}$. We find that in the region described by Eq. (69b), the disorder correlation effect is unimportant, while in the asymptotically small energy region described by Eq. (69a), correlations play a crucial role and $\xi$ has the dependence $\xi\propto k^{-6}$, $E^{-6}$, $\omega^{-6}$, $\lambda^{6}$. If the correlation length is extremely small, then a wide regime given by Eq. (69b), where $\xi\propto k^{-4}$, $E^{-4}$, $\omega^{-4}$, $\lambda^{4}$, should precede the asymptotic regime.

On the other hand, if the average potential $\tilde{u}_{0}$ is nonzero, we can approximate Eq. (68) as

[TABLE]

which applies if $k/k_{d}\ll|\tilde{u}_{0}|$. Therefore, in this case, we have the dependence $\xi\propto k^{-4}$, $E^{-4}$, $\omega^{-4}$, $\lambda^{4}$ in the asymptotically small energy or long wavelength limit. In Ref. [43], from a numerical study of a binary random multilayer model, the $\xi\propto\lambda^{6}$ ($\xi\propto\lambda^{4}$) scaling behavior was obtained in the long wavelength limit, when the average potential was zero (nonzero). These results agree with ours in the asymptotic long wavelength limit. However, the crossover between the two behaviors, which would occur when the average potential was zero, was not investigated in Ref. [43]. An important advantage of our result is that different scaling behaviors for different values of the average potential are incorporated in a single formula, Eq. (68).

V Numerical results

V.1 Incident angle dependence

In this section, we present the results of our comprehensive numerical calculations obtained using the IIM and discuss the dependencies of the localization length on various parameters such as incident angle, disorder correlation length, disorder strength, energy and wavelength. We first consider the incident angle dependance.

In Fig. 1, we consider a short-range correlated dichotomous random potential and plot the normalized localization length $k\xi$ as a function of the incident angle $\theta$ in the weak disorder regime with $s=0.01$ and the strong disorder regime with $s=100$, when the average value of the potential $U_{0}$ is zero and the normalized correlation length $kl_{c}$ is equal to 0.01 and 100. The numerical results obtained using the IIM are compared with the analytical formulas in the weak and strong disorder regimes, Eqs. (51b) and (65). The agreements are perfect. All curves show that the localization length diverges rapidly as $\theta$ approaches to zero and complete delocalization occurs at normal incidence.

In Fig. 2, we plot the normalized localization length $k\xi$ versus $|\theta|$ in a log-log plot in the weak and strong disorder regimes. The curves shown were obtained for different values of the parameters $u_{0}$, $kl_{c}$ and $\sigma^{2}$. We find that all curves show the same divergent scaling behavior $k\xi\propto\theta^{-4}$ near $\theta=0$ regardless of the parameter values.

We have also performed a numerical calculation for a $\delta$-correlated random potential in the weak disorder regime and found the similar $\theta^{-4}$ scaling behavior. In Fig. 3, we plot $k\xi$ versus $\theta$ when the average potential $U_{0}$ is zero and the disorder parameter in the $\delta$-correlated case $g$ is equal to 0.00005 and compare it with the result obtained for the short-range correlated case with $kl_{c}=0.01$ and $\sigma^{2}=0.01$. For the parameters chosen to satisfy the condition $kl_{c}\sigma^{2}=2g$, we find that the agreement between the two results is perfect.

V.2 Disorder correlation length dependence

It is well-known that disorder correlations play a significant role in localization phenomena. In Fig. 4, we plot the normalized localization length $k\xi$ versus normalized correlation length $kl_{c}$ in a log-log plot, when the incident angle $\theta$ is fixed to $30^{\circ}$ and $\epsilon_{0}$ ($=1-u_{0}$) is 1 or 2. We consider the weak disorder regime where $s=0.001$, the strong disorder regime with $s=1000$ and the intermediate disorder regime where $s=5$. In all cases, the localization length shows a non-monotonic dependence on $l_{c}$. As $l_{c}$ increases from zero to infinity, the localization length initially decreases as $\xi\propto{l_{c}}^{-1}$, attains a minimum value at $l_{c}^{\rm min}$ or $\tilde{l}_{c}^{\rm min}$ given by Eqs. (54) and (67), and then increases as $\xi\propto{l_{c}}$. In the weak and strong disorder regimes, the numerical results obtained using the IIM agree very well with the analytical formulas, Eqs. (51b) and (65), respectively. We have checked that in the weak disorder regime, the value of $l_{c}$ at which $\xi$ takes a minimum, $l_{c}^{\rm min}$, is independent of the disorder strength $s$ (or $\sigma$). As $s$ increases, however, this value decreases gradually and, in the strong disorder regime, $\tilde{l}_{c}^{\rm min}$ behaves as $\tilde{l}_{c}^{\rm min}\propto\sigma^{-2}$ as given in Eq. (67).

V.3 Disorder strength dependence

We now consider the dependence of the localization length on the strength of disorder $\sigma^{2}$ in the short-range correlated model. In Fig. 5, we plot the normalized localization length $k\xi$ versus $\sigma^{2}$ in a log-log plot, when $kl_{c}=0.001$, $U_{0}=0$ and $\theta=30^{\circ}$. A very wide range of $\sigma^{2}$ from $10^{-4}$ to $10^{6}$ is considered. First, we notice that the localization length has an extremely sharp dip at $\sigma^{2}=1$. We have carefully checked that the localization length actually vanishes at this value, implying an extreme localization. This unique phenomenon in pseudospin-1 systems arises due to the flat band located at $E=U$ and is directly related to the singularity of the wave equation, Eq. (6), at $E=U$, or equivalently, to that of the invariant imbedding equation, Eq. (22), at $\epsilon=0$. For our dichotomous random potential, the singularity condition corresponds to $E=U_{0}\pm S$, which is equivalent to $\epsilon_{0}=\pm\sigma$. When $U_{0}$ is zero, this condition becomes $\sigma^{2}=1$.

We notice that there are three distinct scaling regions. In the weak disorder region where $\sigma^{2}\ll 1$, the localization length decreases as $\xi\propto\sigma^{-2}$ as $\sigma$ increases from zero. This scaling behavior is independent of the size of the correlation length as can be seen from Eqs. (53a) and (53b) and agrees very well with the analytical formula, Eq. (51b). In the strong disorder region where $\sigma^{2}\gg 1$, we notice that $\xi$ increases monotonically as the disorder strength increases, which implies that localization is destroyed by infinitely strong disorder, similarly to the pseudospin-1/2 case [41]. There are two different scaling behaviors given by Eqs. (66a) and (66b). The crossover between them occurs at the value of $\sigma^{2}$ determined by Eq. (67), namely, $\sigma^{2}=\epsilon_{0}/(2kl_{c})=500$ for the parameter values used here. As we have discussed in Sec. IV.3, in the region where $\xi\propto\sigma^{2}$, the effective wavelength is sufficiently larger than the correlation length and the correlation effect is negligible. In contrast, in the asymptotically strong disorder region where $\xi\propto\sigma^{6}$, the effective wavelength is much smaller than the correlation length and the correlation effect becomes highly relevant. The occurrence of the $\sigma^{6}$ scaling behavior and the crossover between the two scaling regions have not been obtained before.

Next, we consider the case where $U_{0}$ is nonzero. We first point out that the three scaling regions where $\xi$ is proportional to $\sigma^{-2}$, $\sigma^{2}$ or $\sigma^{6}$ also occur in this case, except in the total reflection regime where $|\epsilon_{0}|<|\sin\theta|$. In Fig. 6, we plot the normalized localization length $k\xi$ versus $\sigma^{2}$ in a log-log plot, when $kl_{c}=0.5$, $\theta=30^{\circ}$ and $\epsilon_{0}=0.1$, 0.49, 0.5 and 0.6. We verify that the sharp dips of the curves occur precisely at the $\sigma^{2}$ values equal to ${\epsilon_{0}}^{2}$ as is expected. For the incident angle $\theta=30^{\circ}$, the critical value of $|\epsilon_{0}|$ is equal to 0.5. We find that the localization length shows a nontrivial scaling behavior of the form $\xi\propto\sigma^{-2/3}$ in the weak disorder region at $\epsilon_{0}=0.5$, whereas it scales as $\xi\propto\sigma^{-2}$ if $\epsilon_{0}>0.5$. Contrasting scaling behaviors of this kind are observed in all systems showing the disorder-enhanced tunneling phenomenon [41, 68]. We have performed a similar calculation for the $\delta$-correlated case and found that at the critical value, $\xi$ scales as $\xi\propto g^{-1/3}$, which is equivalent to the result in the short-range correlated case. Precisely the same scaling behavior was previously obtained in the pseudospin-1/2 case [41].

In contrast to the $\delta$-correlated case, the disorder-enhanced tunneling phenomenon in the present case does not occur in the whole total reflection regime, but is limited to the parameter region where both Eq. (52) and the condition $|\epsilon_{0}|<|\sin\theta|$ are satisfied. For the parameter values used here, this gives the bound $0.455<\epsilon_{0}<0.5$. For $\epsilon_{0}=0.49$, we find that $\xi$ initially increases to a maximum and then decreases as $\sigma$ increases, while for $\epsilon_{0}=0.1$, $\xi$ decreases monotonically until the dip at $\sigma^{2}={\epsilon_{0}}^{2}$.

Finally, in Fig. 7, we compare the scaling behaviors in the weak disorder limit between the $\delta$-correlated case and the short-range correlated case, when $U_{0}=0$ and $\theta=30^{\circ}$. We find that the agreement between the two is perfect in the weak disorder limit, if we identify $g=kl_{c}\sigma^{2}/2$. The sharp dips in the short-range correlated case occur precisely at $kl_{c}\sigma^{2}/2=kl_{c}/2$. We observe that the result from the $\delta$-correlated model does not show a sharp dip because our method in this case does not capture the singularity effect due to the flat band.

V.4 Energy and wavelength dependence

In Secs. IV.2 and IV.3, we have discussed the energy and wavelength dependence of the localization length in the weak and strong disorder regimes in the short-range correlated case. We have found that when the average potential $U_{0}$ is zero, there should appear four different scaling regions depending on energy or wavelength. In Fig. 8(a), we show the result of the IIM calculation when $k_{d}l_{c}=0.001$, $U_{0}=0$ and $\theta=30^{\circ}$. The sharp dip occurs at the expected position $k/k_{d}=E/S=1$, which is obtained from the condition $E=U_{0}\pm S=\pm S$. We see clearly that there are four scaling regions, where $\xi$ is proportional to $E^{-6}$, $E^{-4}$, $E^{0}$ or $E^{2}$, or equivalently, to $\lambda^{6}$, $\lambda^{4}$, $\lambda^{0}$ or $\lambda^{-2}$. The IIM result is compared with the analytical formulas, Eqs. (69a), (69b), (60a) and (60b) and the agreements are quite good.

The crossover between different scaling behaviors is predicted to occur at $k/k_{d}=2k_{d}l_{c}=0.002$ in the small energy region and at $k/k_{d}=1/(2k_{d}l_{c}\cos\theta)\approx 577.35$ in the large energy region and the curve shows a good agreement with the predictions. As we have shown in Secs. IV.2 and IV.3, both crossovers occur when the relative magnitudes of the disorder correlation length and the effective wavelength change in the weak and strong disorder regimes respectively. In the $E^{-4}$ and $E^{0}$ scaling regions, the effective wavelength is much larger than the correlation length and the correlation effect is negligible, while in the $E^{-6}$ and $E^{2}$ scaling regions, the effective wavelength is much smaller than the correlation length and the correlation effect is important.

In Fig. 8(b), we compare the curves obtained for different values of $k_{d}l_{c}$. When $k_{d}l_{c}$ is 0.00001, the crossovers should occur at $k/k_{d}=0.00002$ and $k/k_{d}\approx 57735$. Since these values are outside of the range shown here, we should observe only the $E^{-4}$ and $E^{0}$ scaling behaviors, as can be verified from the figure. On the other hand, when $k_{d}l_{c}$ is 0.1, the crossovers occur at $k/k_{d}=0.2$ and $k/k_{d}\approx 5.77$.

In Sec. IV.3, we have proved that when the average potential $U_{0}$ is nonzero, the scaling behavior in the asymptotically small energy or long wavelength limit should be $\xi\propto E^{-4}$, $\lambda^{4}$ instead of $\xi\propto E^{-6}$, $\lambda^{6}$. In Fig. 9, we plot the normalized localization length $k_{d}\xi$ versus $k/k_{d}$ when $\theta=30^{\circ}$, $k_{d}l_{c}=0.1$ and the parameter $\tilde{u}_{0}$ ($=U_{0}/S$) takes the values 0.001, 2 and 1. We expect that sharp dips should occur at $k/k_{d}=\tilde{u}_{0}\pm 1$. This condition gives $k/k_{d}=1.001$ in Fig. 9(a), $k/k_{d}=1,3$ in Fig. 9(b) and $k/k_{d}=2$ in Fig. 9(c), which are confirmed in the figures. In addition, we observe that sharp delocalization peaks appear when the condition $E=U_{0}$ is satisfied. This condition is equivalent to $k/k_{d}=\tilde{u}_{0}$, which is also confirmed in the figures. We have checked carefully that the localization length diverges precisely at $k/k_{d}=\tilde{u}_{0}$.

Except when $\tilde{u}_{0}$ is either 0 or 1, we confirm that the scaling behavior in the asymptotically small energy limit is indeed given by $\xi\propto E^{-4}$. However, if $\tilde{u}_{0}$ is very close to zero as in Fig. 9(a), $\xi$ follows the scaling behavior $\xi\propto E^{-6}$ first and then crosses over to $\xi\propto E^{-4}$ as the energy decreases to zero. The $\xi\propto E^{-2}$ scaling behavior shown in Fig. 9(c) for $\tilde{u}_{0}=1$ is peculiar and is different from the other cases. Our perturbation theory in the weak and strong disorder regimes given in Secs. IV.2 and IV.3 cannot be applied to this case, since the perturbation parameter $s$ is given by $s=\sigma^{2}/(1-\sigma)^{2}$, which approaches to 1 in the small energy limit where $\sigma$ diverges. Therefore $s$ is neither large nor small and the perturbation expansion does not work. Since $U_{0}$ is equal to $S$ when $\tilde{u}_{0}=1$, our dichotomous random potential fluctuates randomly between 0 and $2S$. In Ref. [43], it was reported that in a random superlattice structure where layers with zero potential and those with a random potential were alternated periodically, the localization length scaled as $\xi\propto\lambda^{2}$ in the long wavelength limit. Though this case and our case with $\tilde{u}_{0}=1$ are similar and show the same scaling behavior, there is a subtle difference between them. In our case, the potential fluctuates randomly between 0 and a constant value $2S$, while in Ref. [43], it alternates periodically between 0 and a random value. The common feature is the occurrence of the regions where the potential is identically equal to zero.

We elucidate the scaling behavior in the $\tilde{u}_{0}=1$ case further by calculating the localization length for $\tilde{u}_{0}$ slightly different from 1. In Fig. 10, we plot the normalized localization length $k_{d}\xi$ versus $k/k_{d}$ when $\theta=30^{\circ}$, $k_{d}l_{c}=0.1$ and $\tilde{u}_{0}$ is very close to 1. As is expected, sharp dips occur at $\tilde{u}_{0}\pm 1$ and sharp peaks occur at $\tilde{u}_{0}$. We find that as the energy decreases to zero, the scaling follows $\xi\propto E^{-2}$ first and then crosses over to $\xi\propto E^{-4}$. Therefore the $E^{-2}$ scaling behavior strictly occurs only when $\tilde{u}_{0}=1$.

Finally, in Fig. 11, we show how the localization length diverges when $E$ approaches $U_{0}$, or equivalently, when $\epsilon_{0}$ approaches zero for several parameter values. We find that $\xi$ always diverges as $\xi\propto{\epsilon_{0}}^{-2}$. The sharp delocalization peaks appearing in Figs. 9 and 10 obey this behavior. A similar result was reported in Ref. [43].

VI Conclusion

In this paper, we have studied the Anderson localization of 2D massless pseudospin-1 Dirac particles in a random 1D scalar potential theoretically. We have explored the effect of disorder correlations by solving the Dirac equation with short-range correlated dichotomous random potential for all strengths of disorder. Using the invariant imbedding method, we have calculated the localization length in a numerically exact manner and analyzed its dependencies on incident angle, disorder correlation length, disorder strength, energy, wavelength and average potential extensively over a wide range of parameter values. We have also derived concise analytical expressions for the localization length, which are extremely accurate in the weak and strong disorder regimes. Using the effective wave impedance derived from the pseudospin-1 Dirac equation, we have obtained several delocalization conditions for which the localization length diverges. We have also obtained a condition for which the localization length vanishes. For all cases studied in this paper, we have found that the localization length depends non-monotonically on the correlation length and diverges as $\theta^{-4}$ at normal incidence. As the disorder strength and the energy (or wavelength) vary from zero to infinity, we have found that there appear several different types of scaling behaviors with different values of the scaling exponent. We have explained the crossover between different scaling behaviors in terms of the relative magnitude of the correlation length and the effective wavelength.

We hope our results will stimulate future experiments on localization in pseudospin-1 systems. Our approach can be easily adapted to the case where both the random scalar and vector potentials are present. It is also straightforward to apply our method to other pseudospin-$N$ systems such as pseudospin-3/2 and pseudospin-2 systems. These directions of research will be pursued in the future.

Acknowledgements.

This research was supported by the Basic Science Research Program through a National Research Foundation of Korea Grant (NRF-2019R1F1A1059024) funded by the Ministry of Education.