Non-conventional Anderson localization in a matched quarter stack with metamaterials

We study the problem of non-conventional Anderson localization emerging in bilayer periodic-on-average structures with alternating layers of materials with positive and negative refraction indices $n_a$ and $n_b$. Main attention is paid to the model of the so-called quarter stack with perfectly matched layers (the same unperturbed by disorder impedances, $Z_a=Z_b$, and optical path lengths, $n_ad_a=|n_b| d_b$, with $d_a$, $d_b$ being the thicknesses of basic layers). As was recently numerically discovered, in such structures with weak fluctuations of refractive indices (compositional disorder) the localization length $L_{loc}$ is enormously large in comparison with the conventional localization occurring in the structures with positive refraction indices only. In this paper we develop a new approach which allows us to derive the expression for $L_{loc}$ for weak disorder and any wave frequency $\omega$. In the limit $\omega \rightarrow 0$ one gets a quite specific dependence, $L^{-1}_{loc}\propto\sigma^4\omega^8$ which is obtained within the fourth order of perturbation theory. We also analyze the interplay between two types of disorder, when in addition to the fluctuations of $n_a$, $n_b$ the thicknesses $d_a$, $d_b$ slightly fluctuate as well (positional disorder). We show how the conventional localization recovers with an addition of positional disorder.


Introduction
Due to a remarkable progress in manufacturing one-dimensional (1D) systems with given transport characteristics the interest to a rigorous analysis of various 1D models has greatly increased in recent decades (see, e.g., book [1] and references therein). Of a particular interest are bilayer structures in optics [2] and electromagnetics [3], semiconductor superlattices [4], in the devices with alternating quantum wells and barriers in electronics, etc.
Unlike the disordered models with continuous potentials for which the theory is fully developed, the analysis of the Kroing-Penney type models meets serious theoretical difficulties. One of the first non-trivial problems was a famous effect of band center anomaly in the standard tight-binding Anderson model. As was found numerically [5], the value of the localization length at the band center did not support simple analytical predictions based on the conventional perturbation theory [6]. More careful analysis performed in [7,8] have shown that one has to use specific technics in order to derive correct results. The origin of the discovered effect was found to be related to a kind of resonances emerging for weak disorder. As was later understood, such an anomaly is due to a non-homogeneous distribution of the phase of a wave function, generated in the process of wave propagation along a structure. Similar effects are also known to occur in the vicinity of band edges of energy spectra (see, for example, [9] and references therein).
To date, it is understood that the effect of resonance transmission of a wave through any periodic structure poses a similar problem of a correct description of various transport characteristics for the models that are periodic-on-average. The latter term refers to various kinds of weak disorder added to a strictly periodic potential. This type of systems emerge naturally since in practice on the top of periodic structures such as photonic lattices or electron superlattices, weak random variations of parameters are experimentally unavoidable. Such variations can occur for the width of barriers or wells, electromagnetic or optical characteristics of materials, effective masses of electrons in superlattices, etc. For the frequency of waves far from the resonances the phase distribution of propagating waves is typically constant. In this case the analytical treatment of the transmission coefficient or localization length is relatively easy. However, in the vicinity of band edges or resonances the distribution of wave phase turns out to be highly non-homogeneous, the fact that leads to a dramatic complication of a rigorous analysis (see review [9]).
Another effect which needs specific analytical tools in its analysis is the influence of various correlations in a disorder. As is already understood, the underlying correlations, either of short or long-range type, can significantly suppress or enhance the localization length, therefore, strongly influence the global characteristics of the transmission [9]. In particular, by imposing specific long-range correlations one can practically create the devices with given transmission/reflection characteristics. This effect of enhancement or suppression of the Anderson localization was recently observed experimentally in the waveguides with both the bulk [10] or surface [11] correlated disorder. It should be stressed that specific long-range correlations can emerge naturally in physical systems, and not only in the systems with intentionally included correlations. One of such examples will be discussed in this paper.
New perspectives in creating the devices with unusual transport characteristics are related to specific optic properties of metamaterials embedded in periodic structures [12,13,14,15,16,17,18,19,20,21,22]. A particular system of an increasing interest is an array of two alternating a -and b -layers with equal optical path lengths, n a d a = |n b |d b . Here the a -layers are made of right-handed (RH) material with positive refractive index, n a > 0 and thickness d a . In contrast, the b -layers refer to the lefthanded (LH) material with negative index, n b < 0, and thickness d b . In such a model the phase shift of the wave gained in the a -layer is fully compensated by the subsequent shift in the next b -layer. As a result, the total phase shift vanishes after passing any of N units constructed by (a, b) cells. If, in addition, the layers are matched (Z a = Z b ), this results in a kind of invisibility of the structure for an observer.
Recent numerical data [14] obtained for the array of two matched alternating RH and LH layers with d a = d b and weakly disordered refractive indices n a ≈ |n b | ≈ 1, have demonstrated enormously fast divergence of the localization length, L −1 loc ∝ ω 6 for ω → 0. More detailed numerical study [15] has shown that the power κ in dependence L −1 loc ∝ ω κ increases with the sample length, and approaches the value κ ≈ 8.78. This result is in contrast with the conventional dependence L −1 loc ∝ ω 2 known to occur for many models, both with continuous and periodic-on-average potentials.
As was found in article [23], for the model discussed in [14,15] the phase of wave propagating through the RH-LH array with fluctuating refractive indices is described by a highly non-uniform distribution. This effect is similar to that arising in the tightbinding Anderson model for the energy close to band edges. Further analysis [23] has led to a remarkable conclusion that the localization length diverges in the quadratic approximation in disorder strength. In the next study [24] a new method was suggested allowing one to resolve the above problem of non-conventional Anderson localization. It was shown that the inverse localization length within the fourth-order perturbation theory is described by the asymptotics L −1 loc ∝ σ 4 ω 8 for small ω. Therefore, the found dependencies with κ ≈ 6 have to be attributed to a not large enough size N of samples used in numerical calculations.
The aim of this paper is two-fold. First, we present a full description of the method used in [24] when analyzing the abnormal localization, numerically observed in [14,15]. We generalize the model studied in [14,15,23,24] and give many details which are important for understanding our method, together with a comprehensive discussion of the obtained results. Second, we present new analytical results for the model with two kinds of disorder, namely, with the positional and compositional disorders, and show the interplay between two mechanisms of localization. We show how the abnormal localization reduces to the conventional one in dependence on the relative strengths of the two disorders.

Model formulation: Hamiltonian map
We consider the propagation of an electromagnetic wave of frequency ω in an infinite dielectric array (stack) of two alternating a -and b -layers. Every kind of layers is respectively specified by their thickness d a,b , dielectric permittivity ε a,b , magnetic permeability µ a,b , refractive index n a,b = √ ε a,b µ a,b , impedance Z a,b = µ a,b /n a,b and wave number k a,b = ωn a,b /c. We shall analyze two systems: the homogeneous stack in which both a -and b -layers are made of RH materials, and mixed stack where the a -layers contain the RH material while b -layers are composed of the left-handed (LH) material.
As is known, in the RH medium all optic parameters are positive. On the contrary, in the LH medium the permittivity, permeability and, consequently, the refractive index are negative, however, the impedance remains to be positive. In what follows, we consider the situation when the disorder is originated from random variations of both layer thicknesses d a,b (positional disorder ), as well as from the fluctuations of ε a,b (compositional disorder ). Specifically, we assume a weakness of both types of disorder: and ε an = ε a [1 + η a (n)] 2 , n an = n a [1 + η a (n)], Here the index n enumerates the n-th unit (a, b) cell, and we assume that the magnetic permeabilities µ a and µ b are disorder-independent. The random sequences ̺ a,b (n) and η a,b (n) imposing, respectively, the positional and compositional disorder are specified by white-noise entries with the zero average and small variances σ 2 Hereinafter, the angular brackets . . . stand for the averaging over different realizations of a random structure (ensemble averaging) or along its single realization (spatial averaging), which is assumed to be equivalent due to ergodicity. Numerically, when generating random sequences ̺ a,b (n) and η a,b (n) we use the flat distribution on a finite interval. However, our analytical expressions are valid for any distribution with finite variance.
Within every a -or b -layer the electric field of the wave, obeys the one-dimensional Helmholtz equation with two boundary conditions at the interfaces between neighboring layers, The x-axis is directed along the array of bilayers with x = x i standing for the interface coordinate. The prime implies the derivative with respect to x. The general solution of the wave equation (2.5) inside the n-th unit (a, b) cell can be presented as a superposition of two standing waves, inside a n layer, where x an x x bn ; The coordinates x an and x bn refer to the left-hand edges of successive a n and b n layers, respectively. The thickness of individual layers is defined as The solutions (2.7a) and (2.7b) give rise to useful relations between the electric field E a,b and its derivative E ′ a,b at opposite boundaries of a n and b n layers, E a (x bn ) = E a (x an ) cos ϕ an + k −1 an E ′ a (x an ) sin ϕ an , E ′ a (x bn ) = −k an E a (x an ) sin ϕ an + E ′ a (x an ) cos ϕ an ; (2.9a) The disordered phase shifts ϕ an and ϕ bn depend on the cell index n due to randomized refractive indices n an and n bn , as well as via random thicknesses d an and d bn of the layers, Here the unperturbed phase shifts ϕ a and ϕ b are defined by the expressions: Then, the combination of relations (2.9a), (2.9b) with the boundary conditions (2.6) at the interfaces x i = x bn and x i = x a n+1 , yields the recurrent relations describing the wave transfer through the n-th unit (a, b) cell, Q n+1 = A n Q n + B n P n , P n+1 = −C n Q n + D n P n . (2.12) Here Q n and P n refer to the normalized electric field and its derivative, respectively, taken at left-hand edge of the n-th unit (a, b) cell, The normalization factors in (2.13) contain the unperturbed impedance Z a of a -layers, see definitions (2.2a) and (2.2b). The randomized factors A n , B n , C n , D n read A n = cos ϕ an cos ϕ bn − Z −1 an Z bn sin ϕ an sin ϕ bn , B n Z a = Z an sin ϕ an cos ϕ bn + Z bn cos ϕ an sin ϕ bn , C n Z −1 a = Z −1 an sin ϕ an cos ϕ bn + Z −1 bn cos ϕ an sin ϕ bn , D n = cos ϕ an cos ϕ bn − Z an Z −1 bn sin ϕ an sin ϕ bn . (2.14) As one can see, for a non-zero disorder the coefficients (2.14) depend on the cell index n due to randomized phase shifts ϕ an , ϕ bn and impedances Z an , Z bn . It should be stressed that the phase shifts are influenced by both the compositional and positional disorders, however, the impedances depend on the compositional disorder only, see relations (2.1a), (2.1b), (2.2a), (2.2b) and (2.10a), (2.10b). Note also that the recurrent relations (2.12) belong to the class of area-preserving maps whose determinant equals one, Remarkably, the relations (2.12) can be treated as the Hamiltonian map describing the evolution of trajectories in discrete time n for a linear oscillator subjected to a timedependent parametric force. In such a representation Q n and P n can be treated as the classical coordinate and momentum, respectively. Therefore, the problem of quantum localization can be formally reduced to the analysis of the properties of trajectories in the phase space (Q, P ). For the analytical treatment it is convenient to pass to polar coordinates, namely, to the radius R n and angle θ n , It can be shown that the Hamiltonian map in the radius-angle presentation gets the form, Then the localization length L loc can be derived as the inverse of the Lyapunov exponent λ in accordance with the definition [9,23], Here the average . . . is performed along the "time" n, with a consequent average over different realizations of disorder in order to reduce the fluctuations. Such a representation is very useful since instead of the two-dimensional map for R n and θ n the analysis of the Lyapunov exponent can be done in terms of the one-dimensional map for the phase θ n . It should be, however, noted that the expression (2.18) is not valid in the vicinity of band edges, where special methods have to be used (see [9] and references therein).

Unperturbed periodic counterpart: Matched quarter stack
Without disorder (̺ a,b = 0 and η a,b = 0) all unit cells are identical since the parameters of a -and b -layers do not depend on the cell index n. Therefore, the array of bilayers is periodic with the period d which is the unperturbed size of a unit (a, b) cell, As is known, the transmission through a periodic bilayer stack-structure is governed by the following dispersion relation [1], which defines the Bloch wave number κ ≡ γ/d. We remind that ϕ a and ϕ b are the unperturbed phase shifts given by (2.11). The dispersion relation (3.2) specifies the band-structure in the dependence of the Bloch phase γ on the wave frequency ω. Within the spectral bands where | cos γ| < 1, the solution γ(ω) is real, therefore, electromagnetic waves can propagate through the bilayer stack. Otherwise, γ is purely imaginary inside the spectral gaps where | cos γ| > 1; here the waves are the evanescent Bloch states attenuated on the scale of the order of |γ| −1 . As a result, inside spectral gaps the transmission is exponentially small for the samples with sufficiently large number N of unit (a, b) cells, N ≫ |γ| −1 .
In what follows, we assume that without disorder the basic a -and b -layers are perfectly matched which means that their unperturbed impedances are equal, In this case there are no waves reflected from the interfaces between matched layers. Therefore, in accordance with the dispersion relation (3.2), such a stack-structure is equivalent to a homogeneous medium with the linear spectrum and mean refractive index n, Note that for this structure there are no gaps in the spectrum (3.4), and the Bloch phase is simply the total phase shift of the wave after passing any unit (a, b) cell, Figure 1. (Color online) Unperturbed Hamiltonian map (3.6), (3.7) of periodic bilayer stack-structure.
As we noted above, for a homogeneous RH-RH array all optic parameters are positive quantities. On the other hand, for a mixed RH-LH stack the parameters ε a , µ a , and n a of a -layers are positive, whereas in b -layers ε b , µ b , and n b are negative. Therefore, the phase shift in the LH b -layers is negative, ϕ b = −ω|n b |d b /c, however, due to their definitions both impedances, Z a and Z b , remain to be positive. As one can see, the only difference in the dispersion relation (3.2) is due to the sign before the second term. Specifically, for the RH-LH stack the sign "minus" in (3.2) has to be replaced by "plus", as compared with the RH-RH stack. Correspondingly, in the expressions (3.4) and (3.5) for the RH-LH stack one has to make the substitution, n b → −|n b |. As we show below, the change of the sign in the dispersion relation results in a strong change of properties of the localization length, therefore, in basic characteristics of the transmission.
Let us now address the Hamiltonian map (2.12) -(2.17). For a periodic (without disorder) bilayer stack the coefficients (2.14) are independent of the cell index n. Therefore, in line with (3.3) and (3.5) the Hamiltonian map (2.12) takes the following form Correspondingly, the map (2.17) in the radius-angle presentation transforms into the relations In spite of apparently simple forms of the unperturbed maps (3.6), (3.7) their trajectories can be highly non-trivial. To show this, let us consider the matched quarter stack for which two basic a -and b -layers not only have the equal impedances (3.3), but also the same lengths of optical paths, For the homogeneous RH-RH array the refractive index of any b -layer is positive, n b > 0, therefore, the phase shifts are exactly the same, ϕ a = ϕ b , and the Bloch phase is always non-zero, γ = 2ϕ a . In this case the unperturbed map generates a circle in the phase space (Q, P ) with a fixed radius R n . As for the angle θ n , it changes by the Bloch phase γ in one step of the discrete time n, see Figure 1.
On the contrary, for the mixed RH-LH stack, the refractive index of b -layer is negative, n b < 0, therefore, ϕ b = −ϕ a . This means that the phase shift ϕ a of the wave, gained in any RH layer, is canceled by the subsequent negative shift ϕ b in the next LH layer. Therefore, the phase shift γ vanishes after passing every unit (a, b) cell, γ = 0. Together with the perfect transmission, this means that the structure consisting of a finite set of the RH-LH cells, has to be invisible for an observer. For such a kind of systems, the unperturbed trajectory in the phase space (Q, P ) degenerates into a single stationary point, the fact which is drastically different as compared with the homogeneous RH-RH stack.
Thus, when using the perturbation methods one has to take into account that in the zero order of perturbation the properties of trajectories of the Hamiltonian map remarkably depend on the type of a bilayer structure (fully conventional or having the metamaterial parts). Evidently, this fact should be reflected by the peculiarities of the Anderson localization. In particular, one can expect that the localization length for the RH-LH structure has a non-standard dependence on the disorder, as well as on control parameters of the underlying unperturbed system.

Weak compositional disorder: Phase distribution
Now we are in a position to turn on the disorder. As was shown in [9,28], if the unperturbed impedances of two basic layers are equal, Z a = Z b , the localization length caused by the positional disorder, diverges. Indeed, the Anderson localization is originated from the effect of multiple wave reflections. However, the impedances are independent of the disorder in the layer thicknesses d an , d bn , see definitions (2.1a), (2.1b), (2.2a), (2.2b). Consequently, the matching condition (3.3) remains valid even in the presence of positional disorder. Therefore, there is no any reflection from the interfaces even if these interfaces are randomly appear in the structure due to the positional disorder. Thus, for the structure with matched layers the positional disorder has an impact on transport characteristics only being accompanied by the compositional disorder, since the latter destroys the matching of layers.
Taking into account this fact, as a first step it is reasonable to analyse the localization length contributed by the compositional disorder alone (̺ a,b = 0). In this case all layers have constant thicknesses d a or d b , so that the size d = d a + d b of a unit (a, b) cell is also constant. The effect of positional disorder will be analyzed in the last part of the paper. In order to proceed, it is convenient to rewrite the definitions (3.3), (3.8) for the unperturbed matched quarter stack in the form, Hereafter, the upper sign corresponds to the homogeneous RH-RH array, while the lower sign is associated with the mixed RH-LH stack. In order to develop a proper perturbation theory, we assume the compositional disorder to be weak, In doing so, we substitute (4.1) into (2.14) with ̺ a,b = 0, and expand the coefficients A n , B n , C n , D n up to the second order in the compositional perturbation η a,b (n) ≪ 1.
Using the uncorrelated nature (2.3) of the disorder, from the exact perturbed θ-map (2.17) one can obtain its quadratic approximation, Here the functions U(θ) and W (θ) are defined by According to definition (2.18) of the localization length, one has to know the phase distribution that emerges as a result of iterations of the Hamiltonian map. When this distribution is uniform, the average in (2.18) can be performed relatively easy. However, in specific cases like the standard Anderson model with the energy at the band center or close to the band edges, the analytical derivation of the phase probability density ρ(θ) is a highly non-trivial task. The same kind of problem occurs in our model with RH-LH stacks. In order to obtain the phase distribution, one has to derive the stationary Fokker-Plank equation for ρ(θ). This can be done in the way described in [25,26,8,27]. Specifically, we rewrite the map (4.3) in the continuum limit, where random variables η a (t)dt, η b (t)dt are replaced, respectively, with the Wiener variables dW a = η a (t)dt and dW b = η b (t)dt. As a result, expression (4.5) transforms to the stochastic Itô equation, In accordance with (2.3), the Wiener variables have the white-noise properties [25]: Following the theory of stochastic differential equations [25], one can readily associate the Itô equation (4.6) with the stationary Fokker-Plank equation for the probability density ρ(θ), This equation should be complemented by the condition of periodicity and by the normalization condition, Note that the integration in (4.9) is performed within the interval π due to the periodicity of ρ(θ) with the period π. From the above equation (4.8) one can see that its solution is sensitive to whether the Bloch phase γ is non-zero (RH-RH array) or vanishes (RH-LH array).

RH-RH matched quarter stack
In such a structure the unperturbed optic path lengths of two basic layers are equal, and the Bloch phase (3.5) is non-zero, In this case the term in the Fokker-Plank equation (4.8) containing γ/σ 2 η prevails over the others for any, even arbitrarily small, value of the phase shift ϕ. Indeed, under weak disorder conditions (4.2) and small ϕ ≪ 1 one can get the relations σ 2 η W (θ)/γ ∼ σ 2 η ≪ 1 and σ 2 η U 2 (θ)/γ ∼ σ η (σ η ϕ) ≪ 1, which allow to neglect all the terms with W and U in (4.8). As a result, the phase distribution ρ(θ) is uniform within the first order of perturbation theory, The dimensionless inverse localization length d/L loc is derived by differentiating the map (4.3) with respect to θ and substituting the result into (2.18), with a further expansion of the logarithm within the quadratic approximation in disorder. Then, the subsequent averaging can be readily performed over random entries η a,b (n) and over the angle θ with the uniform distribution function (5.2). After some algebra, one gets [23,24], Note that in the θ-map (4.3) only linear terms in the compositional perturbations η a,b (n) contribute to the Lyapunov exponent, since the last quadratic term vanishes after averaging over θ. The final expression (5.3) is in a complete correspondence with the results previously obtained in the papers [29,23,24]. When the phase shift ϕ is small, the result (5.3) yields the asymptotics, This gives rise to the standard quadratic ω-dependence λ ∝ σ 2 η ω 2 in the limit ω → 0. It should be noted that the localization length L loc (ω) defined by (5.3), exhibits the Fabry-Perot resonances emerging when the phase shift ϕ of the wave passing any layer is multiple to π, i.e., when the wave frequency ω meets the condition ω/c = sπ/n a d a s = 1, 2, 3, . . . .

RH-LH matched quarter stack
The principally different situation emerges for the mixed RH-LH matched quarter stack for which In this case the Bloch phase vanishes independently of the value of the phase shift ϕ. As a result, the functions (4.4) turn out to be related to each other, W (θ) = −U(θ)U ′ (θ), and the Fokker-Plank equation (4.8) reduces to d dθ Solving this equation with the conditions (4.9) of periodicity and normalization, we get a highly nonuniform phase distribution, ρ(θ) = 1 π ϕ 2 − sin 2 ϕ U(θ). The averaging in (6.4) is performed over the angle θ with the distribution function ρ(θ) determined by (6.3). Since the denominator U(θ) in (6.3) is the same as the coefficient in (6.4), we come to a very unexpected result: the Lyapunov exponent vanishes within the second order perturbation theory for any value of the phase shift ϕ [23]. This means that in order to derive a non-vanishing Lyapunov exponent, one has to go beyond the second order perturbation theory. For this, one has to obtain the expressions for both the θ-map and the phase distribution ρ(θ) in the next (fourth order) approximation, which is not a simple task. The problem is that the evaluation of high order terms in ρ(θ) is not possible with the method [25,26,8,27] applied above. The reason is that the θ-map written within the fourth order approximation, cannot be reduced to the Fokker-Plank equation within the linear theory of differential stochastic equations. The only result which here can be drawn is that the Lyapunov exponent for the RH-LH matched quarter stack has to be proportional to σ 4 η , which is in contrast with the conventional quadratic dependence emerging in various disordered models.
In order to proceed further, here we suggest another method that turns out to be very effective for deriving the expression for the Lyapunov exponent. The numerical data in Figure 2 give a hint of what can be done analytically. Specifically, for the RH-RH stack the trajectory in the phase space shown in Figure 2a, looks like it is created by the combination of two processes. The first one is a random increase of the radius for the RH-LH array from map (2.12) (histogram), and from expression (6.3) (solid curve). Due to periodicity, ρ(θ + π) = ρ(θ), only the range 0 θ < π is shown in (c) and (d) (after [24]). R n defined by (2.17). The second process is a uniform-like filling of a circle that is confirmed by the phase distribution ρ(θ) generated by the exact Hamiltonian map. As one can see in Figure 2c, the phase distribution is uniform, the fact which has been already used for the analytical evaluation of the Lyapunov exponent. The right panel in Figure 2 demonstrates completely different effect: for the RH-LH stack the points fill an ellipse characterized by some aspect ratio, and rotated by some angle τ with respect to the axes. This ellipse randomly increases in time, apparently keeping both the aspect ratio and the angle τ . Thus, one can suppose that if to pass properly from the variables Q, P to new ones by both rotating the axes and rescaling ellipse into the circle, the trajectory in new variables will be of the same kind as those shown in Figure 2a. In such a way one can expect that the distribution of a new phase θ n will be uniform, at least, approximately.
Following this idea, we make the linear transformation from "old" variables Q n , P n to the "new" ones Q n , P n , Q n = α 1/2 Q n cos τ + α 1/2 P n sin τ, where the rotating angle τ and rescaling α are the parameters to be specified. Note that in the new variables the expressions (2.12) and (2.16) -(2.18) conserve their forms, however, with the new factors, A n = A n cos 2 τ + (B n − C n ) sin τ cos τ + D n sin 2 τ, B n α −1 = B n cos 2 τ − (A n − D n ) sin τ cos τ + C n sin 2 τ, C n α = C n cos 2 τ + (A n − D n ) sin τ cos τ + B n sin 2 τ, D n = D n cos 2 τ − (B n − C n ) sin τ cos τ + A n sin 2 τ, which replace the old ones (2.14). Now the distribution function ρ( θ) for new phase θ can be found starting from the quadratic expansion of the exact θ-map (2.17) with new coefficients (6.6). Taking into account the relations (6.1), we come to the following equation, The stationary Fokker-Plank equation corresponding to this new θ-map can be obtained in the same way as described before. Then, one gets d Note that the new Fokker-Plank equation (6.8) differs from the old one (6.2) only in the [ϕ + sin ϕ cos(2τ − ϕ)](cos 2 θ + 1). (6.9) One can see that V ( θ) = −U(θ) at α = 1 and τ = 0, when new and old phases coincide. From equation (6.8) with additional conditions similar to (4.9), one readily gets that the phase distribution is, indeed, uniform, ρ( θ) = 1/π, when the function V ( θ) is actually independent of the angle θ, i.e., when its derivative with respect to θ vanishes, With this condition and definition (6.9), we can specify the values of τ and α which determine the transformation from old to new phase-space variables. Also, one can obtain explicitly the θ-independent expression for the function V ( θ), The data presented in Figure 3 confirm the success of our approach: in the new variables the trajectory looks like a circle with the fluctuating radius. These data correspond to our expectation that the distribution of a new phase can be considered as uniform.
The results (6.6), (6.11) allow us to calculate the Lyapunov exponent λ for the mixed RH-LH matched quarter stack with the compositional disorder. For a weak disorder (4.2), the corresponding asymptotics of the θ-map reads Here the function Y ( θ) is defined as Surprisingly, the fourth-order terms in the expansion of the θ-map do not contribute in the fourth-order approximation for the Lyapunov exponent. This happens as a result of the averaging over new phase θ with the uniform distribution. Such an average is in the spirit of the perturbation theory: the next fourth-order terms can be approximated by their average over the phase, with the phase distribution obtained in the second order approximation. Now we substitute the expression (6.12) into definition (2.18) for the Lyapunov exponent. Taking into account that λ vanishes within the quadratic approximation in disorder, we expand the logarithm up to fourth-order terms in the perturbation. After, using the correlation properties (2.3) and averaging over θ with uniform distribution function, we arrive at final result [24], (6.14) Here is the so-called excess kurtosis, the parameter which is specified by the form of the distribution of random variable η a.b (n). In general, ζ can have any value within the range −2 ζ < ∞. In particular, for the Gaussian and flat distributions it equals ζ = 0 and ζ = −6/5, respectively. Expression (6.14) determines the asymptotics for small phase shift ϕ, As one can see, a quite surprising frequency dependence of the Lyapunov exponent, λ ∝ σ 4 η ω 8 , emerges when ω → 0. Thus, the dependence λ ∝ ω 6 numerically obtained for small values of ω in Refs. [14,15], can not be considered as the asymptotical result. From the analysis of the correct expression (6.14) it follows that the dependence λ ∝ ω 6 should be regarded as the intermediate one, apparently emerging due to not sufficiently large lengths N over which the average of λ was performed. Figure 4 demonstrates an excellent agreement between the numerical data for the localization length and analytical results derived above. The data have been obtained by the iteration of the exact map (2.12), with the use of definition (2.18), and compared with the analytical expressions (6.14), (6.16) for the RH-LH quarter stack with ζ = −6/5, as well as (5.3), (5.4) for the RH-RH structure. For sequences of unit-cell numbers N = 10 5 , 10 7 and 10 9 the data are obtained with the ensemble averaging over 100 realizations of disorder, while for N = 10 12 only one realization is used.

RH-LH matched quarter stack: Effect of positional disorder
In the previous sections we presented a full analytical description of a nontrivial effect of compositional disorder in the mixed matched quarter stack with the LH material. On the other hand, as is noted in Section 4, the positional disorder itself in a bilayer stack with matched layers (Z a = Z b ) does not result in the localization [9,28]. Below we show that in the presence of an additional compositional disorder, one has to take into account its non-trivial influence on the localization length. Specifically, we consider an interplay between both disorders, and show how the conventional frequency dependence of the localization length, L −1 loc ∝ ω 2 , (for small ω) emerges due to the influence of the positional disorder.
Since we consider the mixed RH-LH matched quarter stack (4.1), (6.1), it is reasonable to treat from the beginning the transformed Hamiltonian map, which is described by equations (2.12), (2.16) and (2.17) with factors (6.6) instead of (2.14). The angle τ and rescaling parameter α are specified by expressions (6.11).
As before, in order to develop the perturbation theory we assume both compositional and positional disorders to be weak, compare with (4.2), Note that the impedances (2.2a), (2.2b), depend only on the compositional disorder η a,b (n), whereas the phase shifts (2.10a), (2.10b) are randomized due to both compositional η a,b (n) and positional ̺ a,b (n) disorders. The latter is additionally incorporated into the layer thicknesses (2.1a), (2.1b).
The expansion of the coefficients (6.6) and the recurrent relation (2.17) within the quadratic approximation in the disorders η a,b (n), ̺ a,b (n) yields the following θ-map, This expression differs from (6.12) in the second and fourth terms that take into account the effect of positional disorder. Functions V (ϕ) and Y ( θ) are respectively determined by (6.11) and (6.13). As for the function G( θ), it is given by 3) The stationary Fokker-Plank equation corresponding to the θ-map (7.2) reads As one can see, in this equation the first term containing the variance σ 2 η , is responsible for the compositional disorder, and the second one (with σ 2 ̺ ) arises due to positional disorder. As both disorders do not correlate with each other, the equations (7.2) and (7.4) do not include the cross-correlated terms with the product σ η σ ̺ .
The solution to (7.4) satisfying the conditions of periodicity and normalization (4.9) is with the normalization function I(ϕ) determined as Remarkably, when the positional disorder vanishes the distribution function ρ( θ) reduces to the uniform one, This fact is in a complete accordance with the result obtained in the previous section. The dimensionless inverse localization length d/L loc is obtained with the use of the map (7.2) in the same way as described above. Taking into account the white-noise properties (2.3), the averaging over θ with the probability density (7.5), (7.6) From Figures 5 and 6 one can see that expression (7.8) displays a very good agreement with the numerical data for the Lyapunov exponent obtained from the exact map (2.12) for the mixed RH-LH matched quarter stack. It is important to emphasize that this expression is applicable for any ratio between the variances σ 2 η and σ 2 ̺ . Also, one should note that it is valid within a wide frequency range restricted to the requirement (7.1) of a weak disorder only. Remarkably, in the absence of compositional disorder, σ 2 η = 0, the Lyapunov exponent (7.8) vanishes even if the positional disorder is nonzero, σ 2 ̺ = 0. This fact is a direct consequence of the layer matching (3.3), and is in the accordance with the previously obtained results [9,28]. The first term in (7.8) prevails over the second one and correctly describes the localization length when the positional disorder predominates, σ 2 η ≪ σ 2 ̺ , or when both disorders are of the same order of magnitude, σ 2 η ∼ σ 2 ̺ , see Figure 5. The most important and interesting consequence of expression (7.8) is that at any finite but sufficiently small strength σ 2 η of the compositional disorder, its second term is negligible. For this case we arrive at simple, however, non-trivial result, This asymptotics coincides with the Lyapunov exponent (5.3) for the RH-RH matched quarter stack and provides the conventional quadratic frequency dependence, λ ∝ σ 2 η ω 2 when ω → 0, see Figure 5. One should emphasize that the result (7.9) is valid when the positional disorder predominates. However, due to the layer matching (3.3), it does not contribute to the Lyapunov exponent. Thus, expression (7.9) turns out to be provided by small compositional disorder, instead of a large positional one. In spite of this fact, the presence of positional disorder plays the fundamental role, changing the abnormal octal frequency dependence (6.16), L −1 loc ∝ σ 4 η ω 8 , into the convetional quadratic one, L −1 loc ∝ σ 2 η ω 2 , when ω → 0. The second term in (7.8) coincides with the Lyapunov exponent (6.14) caused by the compositional disorder only. However, when the positional disorder is also present, this term can significantly contribute to the Lyapunov exponent if the compositional disorder prevails over the positional one, σ 2 ̺ ≪ σ 2 η . At first glance, in this case one can neglect the terms with σ 2 ̺ in the radicand of the first summand in (7.8) and in the radicand of the normalization function (7.6). In other words, one can suggest that the distribution function (7.5) is uniform (7.7). However, a more careful analysis shows that this is not true at relatively low frequencies ω, when the phase shift ϕ is sufficiently small, ϕ 6 σ 4 η ≪ σ 2 ̺ ≪ σ 2 η . For such frequencies the Lyapunov exponent again obeys the standard quadratic ω-dependence originated from the first term in (7.8) and is described by expression (7.9). Nevertheless, when the condition σ 2 ̺ ≪ σ 2 η is met, the second term does contribute to the localization length in the regime of high and intermediate frequencies.
This fact is clearly displayed in Figure 6.

Conclusion
In this paper we have studied two 1D layered models which consist of both right-handed (RH) and left-handed (LH) materials for alternative a -and b -layers of thicknesses d a and d b , respectively. The problem considered above is an analytical evaluation of the Lyapunov exponent whose inverse value is the localization length. As was found in a set of numerical studies [14,15], for the RH-LH model a quite unexpected Andersontype localization emerges, that has led to an intensive discussion in literature. In previous papers [23,24] we have shown that the origin of this abnormal localization can be associated with a highly non-homogenous distribution of phase of the wave propagating along the bilayer structure. Unexpectedly, it was discovered [23] that with weak fluctuations of dielectric constants ε an , ε bn the Lyapunov exponent vanishes in the second order of perturbation theory. This fact has shed light on the origin of the observed non-conventional localization. On the other hand, it was understood that the analytical treatment of this localization is of a very difficult task since one has to develop a proper perturbation theory beyond the second-order approximation.
The problem of this unusual localization was rigorously solved with the use of a new method, and brief communication was published in [24]. Above, we have presented an analysis of the problem with many details that are important to understand our approach, together with a discussion concerning the mechanism of the observed phenomenon.
We have to emphasize that originally the effect of abnormal localization was observed in [14,15] for a quite specific bilayer model for which without disorder the RH-LH array represents an ideal mixed stack. Specifically, it was assume that ε a = µ a = 1, ε b = µ b = −1, hence the impedances are trivially equal, Z a = Z b = 1. The additional assumption was that the layer thicknesses are also equal, d a = d b . In the present study we consider a more general model of the so-called matched quarter stack for which the condition of equal optical pass lengths, n a d a = |n b |d b , and the equality of the impedances, Z a = Z b , are independent of each other. Thus, our results can be applied to a more general class of physical systems.
In our approach we use the reduction of the underlying problem of a quantum localization to the study of the classical Hamiltonian map of a quite transparent structure. In this way the evaluation of the localization length can be performed in terms of local instability of trajectories in the phase space of the corresponding classical model. Then, with the transformation to radius-angle variables in the phase space, the equations allowing to evaluate the Lyapunov exponent take a quite simple form, and the problem is reduced to the analytical treatment of the angle distribution which emerges after a sufficiently long time of iteration of the classical map. It should be noted that the angle in this classical map is nothing but the phase of wave function in the quantum problem. Therefore, even with the reduction of the quantum model to the classical one the physical correspondence between two representations remains transparent.
Our approach allows one to resolve the problem of the Anderson localization for bilayer RH-LH arrays with the matched layers, and derive the localization length. As we show, the peculiarity of this model is entirely due to the zero value of the unperturbed Bloch phase that happens when the unperturbed impedances are equal. We rigorously derive the expression for the Lyapunov exponent which appears to be defined by the fourth order of perturbation theory, and which is valid for any value of frequency ω. Our results prove that for small ω the localization length is enormously large, L loc ∝ σ −4 η ω −8 . The generalization of our result is due to the θ-map (4.3), according to which the same dependence for L loc is expected to occur when the effective unperturbed refractive index n = (n a d a − |n b |d b )/(d a + d b ), see (3.4), is sufficiently small, n ≪ σ 2 η . Thus, the obtained results show that the non-conventional localization observed numerically in Refs. [14,15,16,17] is very fragile with respect to the choice of model parameters, and can be hardly observed experimentally.
It is interesting that the discussed effect of anomalous localization can be considered as a kind of strongly correlated disorder, one part of which is embedded into the randomized phase shifts ϕ an , ϕ bn and the other is absorbed by the impedances Z an , Z bn . As one can see, the compositional disorder, even if it is of a white-noise type, results in strong correlations between fluctuations of phase shifts and impedances. To the best of our knowledge, this phenomenon of a natural emergence of correlated disorder from the white-noise disorder was never discussed in literature.
A special attention in our study was paid to the model with two kinds of disorder. Specifically, we have considered the case when in addition to the compositional disorder with the fluctuating values of dielectric permittivities ε an , ε bn , the positional disorder is also included (therefore, the thicknesses d an and/or d bn also slightly fluctuate). It is extremely important that the positional disorder randomly affects the phase shifts only and do not contribute to the impedances. As one can see, the positional disorder plays a fundamental role for the destroy of strong correlations between the disordered phase shifts and impedances, therefore, for the recovering of the conventional dependence L −1 loc ∝ σ 2 η ω 2 when ω → 0. The positional disorder in both thicknesses d an and d bn has been earlier analysed in [28]. Specifically, the general case of correlated disorder was analytically studied by assuming any kind of statistical correlations in thicknesses of a -and b -layers, as well as the inter-correlations between the fluctuations of two thicknesses. It was shown that for any ratio between d a and d b the localization length is governed by the unique expression, no matter whether the structure represents RH-RH or RH-LH array. The main result is that for such structures the phase distribution is flat, therefore, the analytical expressions for the Lyapunov exponent can be obtained within the standard perturbation theory. The same statement stems from the analysis of RH-RH and RH-LH arrays with fluctuating refractive indices [29], however, only when n a d a = |n b |d b in the latter case. These results also indicate that the abnormal localization emerges in a very specific situation, and not only due to the inclusion of left-handed materials into the structure.
Finally, we would like to note that the effect of vanishing of the unperturbed Bloch phase shift γ in the θ-map (4.3), is a typical effect occurring at band edges in the tight-binding Anderson model as well as in the Kronig-Penney models. Indeed, in these models after one period of perturbation the unperturbed Bloch phase γ vanishes when approaching the band edges. This results in a very specific non-homogeneous distribution of the wave phase, therefore, to a non-standard expression for the Lyapunov exponent. However, in the considered model of the mixed RH-LH matched quarter stack, the peculiarity is that vanishing Block phase emerges independently of the value of frequency ω. This is a principal difference as compared with the localization occurring in the Anderson and Kronig-Penney models for which the zero Bloch phase emerges at band edges, therefore, only for specific values of frequency.