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 na and nb. Attention is paid to the model of the so-called quarter stack with perfectly matched layers (the same unperturbed by disorder impedances, Za = Zb, and optical path lengths, nada = |nb|db, with da and db 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 Lloc is enormously large in comparison to 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 Lloc for weak disorder and any wave frequency ω. In the limit ω → 0 one gets a quite specific dependence, L−1loc∝σ4ω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 na and nb, the thicknesses da and db slightly fluctuate as well (positional disorder). We show how conventional localization recovers with the addition of positional disorder.


Introduction
Due to the remarkable progress in manufacturing one-dimensional (1D) systems with given transport characteristics, interest in a rigorous analysis of various 1D models has greatly increased in recent decades (see e.g. book [1] and the references therein). Of particular interest are bilayer structures in optics [2] and electromagnetics [3], semiconductor superlattices [4], in devices with alternating quantum wells and barriers in electronics, etc.
Unlike 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 the famous effect of a 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] has shown that one has to use specific techniques in order to derive correct results. The origin of the discovered effect was found to be related to a kind of resonance emerging for weak disorder. As was later understood, such an anomaly is due to the non-uniform 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 the band edges of energy spectra (see e.g. [9] and the 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 system emerges 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 and 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, the distribution of the wave phase turns out to be highly non-uniform in the vicinity of band edges or resonances, a fact that leads to the 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, underlying correlations either of the 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 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, not only in systems with intentionally included correlations; one such example will be discussed in this paper.
New perspectives in creating 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 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 the 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 an enormously fast divergence of the localization length, L −1 loc ∝ ω 6 for ω → 0. A 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-onaverage potentials.
As was found in [23], for the model discussed in [14,15] the phase of the wave propagating through the RH-LH array with fluctuating refractive indices is described by a highly nonuniform distribution. This effect is similar to that arising in the tight-binding Anderson model for the energy close to the band edges. Further analysis [23] has led to a remarkable conclusion that the localization length diverges in the quadratic approximation of disorder strength. In the next study [24] a new method was suggested allowing one to resolve the above problem of nonconventional 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 dependences with κ ≈ 6 have to be attributed to the not large enough size N of samples used in numerical calculations.
The aim of this paper is two-fold. Firstly, 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. Secondly, 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, when changing 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 layer 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 the mixed stack where the a -layers contain the RH material while b -layers are composed of the 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 positive. In what follows, we consider the situation when the disorder originates 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 nth unit (a, b) cell; 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 , σ 2 η 1, 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 in 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 1D 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 nth unit (a, b) cell can be presented as a superposition of two standing waves E a (x) = E a (x a n ) cos k an (x − x a n ) + k −1 an E a (x a n ) sin k an (x − x a n ) (2.7a) inside a n layer, where x a n x x b n , The coordinates x a n and x b n 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 b n ) = E a (x a n ) cos ϕ an + k −1 an E a (x a n ) sin ϕ an , E a (x b n ) = −k an E a (x a n ) sin ϕ an + E a (x a n ) cos ϕ an , (2.9a) (2.9b) 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) and (2.9b) with the boundary conditions (2.6) at the interfaces x i = x b n and x a n+1 , yields the recurrent relations describing the wave transfer through the nth 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 the left-hand edge of the nth unit (a, b) cell a E a (x a n ), P n = (c/ω)Z 1/2 a E a (x a n ). (2.13) 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 and 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 , (2.14) 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 .
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 and 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 time-dependent 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 , Q n = R n cos θ n , P n = R n sin θ n . (2.16) 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 1D 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 aand 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 of the Bloch phase γ dependent 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 a sufficiently large number N of unit 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 a 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, 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 alayers are positive, whereas in b -layers the parameters ε 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 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 the properties of the localization length, therefore, in the 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: Q n+1 = Q n cos γ + P n sin γ , P n+1 = −Q n sin γ + P n cos γ . (3.6) 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 length of optical path 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 nonzero, γ = 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 RH-LH cells has to be invisible for an observer. For such a kind of system, the unperturbed trajectory in the phase space (Q, P) degenerates into a single stationary point, a fact which is drastically different when compared with the homogeneous RH-RH stack.
Thus, when using 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 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 the 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 , then the localization length caused by the positional disorder diverges. Indeed, the Anderson localization originates 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) and (2.2b). Consequently, the matching condition (3.3) remains valid even in the presence of positional disorder. Therefore, there is no reflection from the interfaces, even if these interfaces 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 this fact into account, as a first step it is reasonable to analyze 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 the 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 and 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 the definition (2.18) of 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 [8,[25][26][27]. Specifically, we rewrite the map (4.3) in the continuum limit where random variables η a (t) dt and η 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] dW a,b = 0, dW a dW b = 2σ 2 η δ ab dt. (4.7) 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 ρ(θ), 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).

Right-handed (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 allows one 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 θ n and substituting the result into (2.18), with a further expansion of the logarithm within the quadratic approximation in disorder. Then, 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] d/L loc ≡ λ = σ 2 η sin 2 ϕ, ϕ = ωn a d a /c.

RH-left-handed (LH) matched quarter stack
A principally different situation emerges for the mixed RH-LH matched quarter stack for which In this case the Bloch phase vanishes independent 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 non-uniform phase distribution The Lyapunov exponent can be derived according to its definition (2.18) and quadratic approximation (4.3) of the Hamiltonian θ -map, taking into account the conditions (6.1). In the second order of the perturbation theory in the weak compositional disorder η a,b (n) one gets 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 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 an evaluation of high order terms in ρ(θ ) is not possible with the method [8, 25-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 that 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 gives a hint of what can be done analytically. Specifically, for the RH-RH stack, the trajectory in the phase space shown in figure 2(a) looks like it is created by a combination of two processes. The first one is a random increase of the radius 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 2(c), the phase distribution is uniform, a fact which has been already used for the analytical evaluation of the Lyapunov exponent. The right panel in figure 2 demonstrates a 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 in order to pass properly from variables Q, P to new ones by both rotating the axes and rescaling the ellipse into a circle, the trajectory in new variables will be of the same kind as those shown in figure 2(a). 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 the 'old' variables, Q n and P n , to the 'new' ones,Q n andP n , Q n = α 1/2 Q n cos τ + α 1/2 P n sin τ, P n = −α −1/2 Q n sin τ + α −1/2 P n cos τ, (6.5) 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Ã 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 τ, (6.6) 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 the 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 Note that the new Fokker-Plank equation (6.8) differs from the old one (6.2) only in the function 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 thẽ θ-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 a fluctuating radius. These data correspond to our expectation that the distribution of a new phase can be considered as uniform.
The results (6.6) and (6.11) allow us to calculate the Lyapunov exponent λ for the mixed RH-LH matched quarter stack with a 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 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 a uniform distribution function, we arrive at the 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 [14,15], cannot 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) and (6.16) for the RH-LH quarter stack with ζ = −6/5, as well as (5.4) for the RH-RH structure. For N = 10 6 and 10 8 , the data are obtained with an additional average over 100 realizations of disorder and 10 realizations for N = 10 10 , while for N = 10 12 only one realization is used.

RH-LH matched quarter stack: effect of positional disorder
In previous sections we presented a full analytical description of a non-trivial 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, compared with (4.2), σ 2 η 1, (σ η ϕ) 2 1 and (σ ϕ) 2 1. (7.1) Note that impedances (2.2a) and (2.2b) depend only on the compositional disorder η a,b (n), whereas phase shifts (2.10a) and (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) and (2.1b). An 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 the 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 ρ(θ ) → 1/π for σ 2 → 0.
This fact is in 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) yields 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 non-zero, σ 2 = 0. This fact is a direct consequence of the layer matching (3.3) and is in 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 a simple, however, non-trivial result This asymptotic expression 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 a fundamental role, changing the abnormal octal frequency dependence (6.16), L −1 loc ∝ σ 4 η ω 8 , into a conventional 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

Conclusion
In this paper we have studied two 1D layered models which consist of both RH and 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], a quite unexpected Anderson-type localization emerges for the RH-LH model, which 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-uniform distribution of the phase of the wave propagating along the bilayer structure. Unexpectedly, it was discovered [23] that with weak fluctuations of dielectric constants ε an and ε bn , the Lyapunov exponent vanishes in the second order of the 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 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; a brief communication was published in [24]. Above, we have presented an analysis of the problem with many details that are important in order 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 assumed that ε a = µ a = 1 and ε 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 this 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 one to evaluate the Lyapunov exponent take a quite simple form; 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 the 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 (the unperturbed impedances are equal, Z a = Z b ) and derive the localization length caused by compositional disorder. 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 RH-LH structure represents a quarter stack (optical pass lengths are the same, n a d a = |n b |d b ). 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 . Note that the effective refractive index n is zero for the RH-LH quarter stack, see (3.4). If n = (n a d a − |n b |d b )/(d a + d b ) is non-zero, 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 n is sufficiently small, n σ 2 η . On the contrary, when n is not small, the RH-LH structure no longer belongs to a quarter-stack class (n a d a = |n b |d b ). Then, the Bloch phase is of finite value and the Lyapunov exponent is described by the universal expression valid for both RH-LH and RH-RH systems. This expression has the standard quadratic frequency dependence at the bottom of the spectrum, L loc ∝ σ −2 η ω −2 [9,23,29]. As was mentioned in the review [9], one of the important and quite specific features of the mixed RH-LH quarter stack is that it can propagate waves only when the layers are perfectly matched. Otherwise, when the basic a -and b -layers are mismatched (ϕ a = −ϕ b = ϕ, but Z a = Z b ), in accordance with the dispersion relation (3.2), the spectral bands disappear; therefore, the transmission is absent, apart from a discrete set of resonant frequencies (5.5) for which ϕ = sπ and γ = 0. A numerical study of the Anderson localization in such a situation was suggested in [19]. However, the corresponding analysis seems to be highly non-trivial, therefore the problem remains open. Thus, the obtained results show that the non-conventional localization observed numerically in [14][15][16][17] is very fragile with respect to the choice of model parameters and can be hardly observed experimentally. The fragility of the non-conventional localization was stated for the first time in [15].
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 and ϕ bn , the other is absorbed by the impedances Z an and Z bn . As one can see that even if the compositional disorder is of a white-noise type, it 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 has never been discussed in literature.
In our study special attention 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 in destroying strong correlations between the disordered phase shifts and impedances, therefore, for recovering the conventional dependence L −1 loc ∝ σ 2 η ω 2 when ω → 0. The positional disorder in both thicknesses d an and d bn has been analyzed 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 intercorrelations 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 a RH-RH or RH-LH array. The main result is that for such structures the phase distribution is uniform; therefore, the analytical expressions for the Lyapunov exponent can be obtained within the standard perturbation theory. As we have mentioned above, the same statement stems from the analysis of RH-RH and RH-LH arrays with fluctuating refractive indices [9,23,29], however, only when n a d a = |n b |d b in the latter case. Moreover, in the standard case of mismatched no-quarter-stack bilayer array the contributions of two uncorrelated (positional and compositional) disorders are universal and additive for both RH-RH and RH-LH stacks. Therefore, the frequency dependence has to be universal, L −1 loc ∝ ω 2 when ω → 0. These results indicate once more that the abnormal localization emerges in a very specific situation and is not only due to the inclusion of LH materials into the structure.
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-uniform 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 a vanishing Bloch 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.
In our study we did not ask an important question about the role of absorption or amplification that can be present in experiments. As was shown in [30,31], the addition of a weak absorption in the standard tight-binding Anderson model results in another type of nonconventional Anderson localization [32], in comparison with that discussed above. The analysis of the question whether (or how) the absorption or amplification (or both) can influence the localization in quarter stacks with metamaterials, needs special treatments.