Anderson localization and Brewster anomaly of electromagnetic waves in randomly-stratified anisotropic media

Anderson localization of p-polarized waves and the Brewster anomaly phenomenon, which is the delocalization of p-polarized waves at a special incident angle, in randomly-stratified anisotropic media are studied theoretically for two different random models. In the first model, the random parts of the transverse and longitudinal components of the dielectric tensor, between which the longitudinal component is the one in the stratification direction, are assumed to be uncorrelated, while, in the second model, they are proportional to each other. We calculate the localization length in a precise way using the invariant imbedding method. From analytical considerations, we provide an interpretation of the Brewster anomaly as a phenomenon arising when the wave impedance is effectively uniform. Similarly, the ordinary Brewster effect is interpreted as an impedance matching phenomenon. We derive the existence condition for the Brewster anomaly and concise analytical expressions for the localization length, which are accurate in the weak disorder regime. We find that the Brewster anomaly can arise only when disorder is sufficiently weak and only in the second model with a positive ratio of the random parts. The incident angle at which the anomaly occurs depends sensitively on the ratio of the random parts and the average values of the tensor components. In the cases where the critical angle of total reflection exists, the angle at which the anomaly occurs can be either bigger or smaller than the critical angle. When the transverse and longitudinal components are uncorrelated, localization is dominated by the the transverse component at small incident angles. When only the longitudinal component is random, the localization length diverges as θ−4 as the incident angle θ goes to zero and is also argued to diverge for all θ in the strong disorder limit.

Since the discovery of the BA by Sipe et al [9], many authors have discussed different aspects of this phenomenon. Jordan et al have studied the BA occurring in randomly-layered anisotropic media consisting of alternating isotropic-uniaxial media using numerical calculations based on the transfer matrix method [17]. They have found that in the cases called mixed stacks, the BA is suppressed and does not occur, while it can occur in the cases called binary stacks. A similar model based on alternating metamaterial-uniaxial randomly layered stacks has been studied by del Barco et al using the transfer matrix method, where the BA is again found to be suppressed [19].
In this paper, we will present a unique perspective on the BA that it is a phenomenon arising when the effective wave impedance is uniform and non-random, which is made possible only in weakly-disordered media. By a similar argument, we will also argue that the ordinary Brewster effect arises when the wave impedance is completely matched throughout the space. Using an analytical method based on the invariant imbedding theory [24,25], we will derive precise conditions for the occurrence of the BA in randomly-stratified anisotropic media and derive concise analytical expressions for the localization length in the weak disorder regime for two different random models. These results will be compared with more accurate numerical results obtained using the invariant imbedding method and also with the previous results obtained for randomlylayered anisotropic media [17,19]. In addition, we will derive some interesting properties of localization in anisotropic media from general analytical considerations.

Model
We consider a random uniaxial medium, the dielectric permittivity tensor of which is diagonalized in the coordinate system (x, y, z) and is written as The medium is stratified along the z axis and the transverse and longitudinal tensor components, ò ⊥ and ò P , are random functions of z only. Plane electromagnetic waves of frequency ω and vacuum wave number k 0 (=ω/c) are assumed to propagate in the xz plane. Then the wave equations for the s-and p-polarized waves are completely decoupled. In this paper, we are only interested in the propagation of p waves, for which the y component of the magnetic field satisfies where q is the x component of the wave vector and a prime denotes a differentiation with respect to z. We assume that an inhomogeneous anisotropic medium of thickness L lies in 0zL and the waves are incident obliquely from a uniform dielectric region (z>L) and transmitted to another uniform dielectric region (z<0). The incident and transmitted regions are filled with ordinary isotropic media of the same kind, where ò (=ò 1 ) is a scalar quantity. When θ is the angle of incidence, q is equal to k sin q, where k k 0 1  = . From now on, we will assume that ò ⊥ and ò P are always normalized by ò 1 to simplify the notations, unless otherwise explicitly stated.
We consider two different random models. In Model I, we assume that ò ⊥ and ò P are independent random functions of z and satisfy The notation á ñ  denotes averaging over disorder and g^and g  are independent parameters characterizing the strength of disorder. On the other hand, in Model II, we consider the situation where the random components δò ⊥ (z) and δò P (z) are not independent, but proportional to each other such that where f is a real constant.

Invariant imbedding method
Since the BA occurs only for p waves, we focus on that case here. We consider a p wave of unit magnitude incident on the anisotropic medium. Using the invariant imbedding method and starting from equation (2), we derive exact differential equations satisfied by the reflection and transmission coefficients, r and t: where p (=kcosθ) is the negative z component of the wave vector in the incident and transmitted regions. We use equation (6) for the precise numerical calculation of the localization length ξ defined by where T is the transmittance given by T t 2 = | | .
The invariant imbedding equations for r and t, equation (6), are stochastic differential equations with random coefficients. In order to deal with the random term ò P appearing in the denominators of the coefficients in equation (6) by using known methods, we assume that the disorder in ò P is sufficiently weak so that We point out that this is the only approximation used in the present work. In contrast, the disorder in ò ⊥ can be of arbitrary strength. From the general considerations presented in section 4, we will show that the BA can occur only when disorder is sufficiently weak. Therefore, the condition given in equation (8) is one of the necessary conditions for the existence of the BA, rather than an approximation.

Model I
In order to obtain the localization length, we need to compute the average T L ln á ñ ( ) in the L  ¥ limit. The nonrandom differential equation satisfied by T ln á ñ can be obtained using the second of equation (6), Novikov's formula [26] and takes the form where Z n (n=1, 2) is equal to r n á ñ and the parameters C 0 , C 1 and C 2 are defined by The dimensionless disorder parameters g ⊥ and g P are given by In Model I, the random terms δò ⊥ and δò P are uncorrelated. This fact has played an important role in deriving equation (9). In the l  ¥ limit, the left-hand side of equation (9) approaches asymptotically to a constant equal to (kξ) −1 .
To calculate Z 1 and Z 2 for use in equation (9), we derive an infinite number of coupled nonrandom differential equations satisfied by Z n , where n is an arbitrary nonnegative integer, using the first of equation (6) and Novikov's formula. These equations turn out to take the form The initial conditions for Z n ʼs are Z 0 =1 and Z n (l=0)=0 for n>0. In the l  ¥ limit, the left-hand sides of these equations vanish and we obtain an infinite number of coupled algebraic equations, which are much easier to solve numerically than the coupled differential equations. The moments Z n with n>0 are coupled to one another and their magnitudes decrease rapidly as n increases. Based on this observation, we solve these equations numerically by a systematic truncation method [27].

Model II
In Model II, δò ⊥ (z) and δò P (z) are not independent, but proportional to each other. This condition leads to completely different equations for Z n and T ln á ñ for p waves. The equation for Z n in this case is written as where the parameters D 1 and D 2 are defined by tan cos , 1 tan cos . 14 1 2 The equation for the localization length takes the form

General considerations on the existence condition of the Brewster anomaly and the properties of localization 4.1. Argument based on the impedance matching condition
There is a very simple interpretation of the BA phenomenon, which has never been, to our knowledge, advocated before. Based on this interpretation, it is possible to explain both the Brewster effect and the BA phenomenon in a unified way. Furthermore, we can deduce some interesting properties of localization in anisotropic media. We begin by rewriting the wave equation, equation (2), in the following equivalent form: In these expressions, we remind again that ò ⊥ and ò P are quantities normalized by ò 1 . Therefore, in the incident and transmitted regions where ò ⊥ =ò P =1, η is equal to 1 for all θ. We notice that the wave equation written in the above form looks identical to that for p waves propagating normally in a medium with the wave impedance given by η(z). Before discussing the BA, it is instructive to examine the ordinary Brewster effect from the viewpoint of impedance matching. It is well-known that if the entire medium has a uniform impedance, waves are completely transmitted without any backward reflection. In our case, the uniform impedance condition requires η to be equal to 1 in the entire slab. From equation (17), it is straightforward to derive the incident angle θ b , which is nothing but the ordinary Brewster angle, for total transmission of p waves. We obtain tan 1 Obviously, the right-hand side of the above equation has to be positive for b q to exist. The same result has been obtained long ago by other authors [28,29]. In isotropic media, we have ò ⊥ =ò P (=ò). Then we reduce equation (18) to the well-known expression for the Brewster angle, The BA is a delocalization phenomenon arising at a special incident angle, θ B , when ò ⊥ and ò P are random functions of z. In order for delocalization to occur, the impedance η needs to be either a real constant or a realvalued nonrandom function of z. From the functional form of equation (17), we find that this cannot be realized if ò ⊥ and ò P are random functions of arbitrary strength of disorder. For sufficiently weak disorder, however, we substitute equation (3) to the first order in δò ⊥ and δò P . The only nontrivial possibility for η to be nonrandom in equation (19) is when a bb sin sin , 20 while both δò ⊥ and δò P are nonzero. If δò P is zero and b is equal to sin 2 θ, η becomes zero and the wave does not propagate. Therefore δò ⊥ and δò P have to be proportional to each other, as in our Model II, where δò P =fδò ⊥ .
We finally obtain b b af sin , 21 where (b+af ) has to be positive and bigger than b 2 to have a solution for θ B . In the case of isotropic media with f=1 and b=a, this reduces to the well-known result, a sin 2 B q = , derived originally by Sipe et al [9]. By substituting equation (21) into equation (19), we obtain the effective wave impedance when a p wave is incident at θ B given by In order to have a propagating wave, the wave impedance needs to be real, which gives an additional constraint such that f has to be positive. In other words, the random functions δò ⊥ (z) and δò P (z) need to be always of the same sign. When the constraints b+af>b 2 and f>0 are satisfied, the angle θ B is well-defined and the impedance η B is a positive real constant, which is not generally equal to 1. Since η B is not matched to that of the incident region in general, the wave incident at θ B on a randomly-stratified slab of finite thickness is partially reflected and the disorder-averaged transmittance is smaller than 1 and depends on the thickness. From the simple form of η in equation (17), we can also deduce several additional properties of localization in anisotropic media. If the incident angle θ is zero, then the dependence on ò P disappears in equation (17), which reduces to 2 1  h =^-. This implies that if ò ⊥ is nonrandom, Anderson localization does not occur at θ=0 for any random function ò P and the localization length ξ diverges. In addition, when θ is sufficiently close to zero, we find from equation (19) that the random term  d  has a coefficient proportional to θ 2 , which suggests that the strength of the   disorder, g P , will always appear as multiplied by θ 4 in this regime. In later sections, we will present an analytical formula and numerical results showing that, in the presence of only the ò P disorder, ξ is indeed proportional to g 4 1 q - ( ) when θ is sufficiently small. Another observation we make about the θ=0 case is that if ò ⊥ is equal to 1, that is, if the transverse tensor component is matched to the permittivity of the incident region, then the impedance is 1 in all regions of space and therefore the transmission has to be perfect regardless of the form of   . We notice that this behavior shows close similarity to the Klein tunneling of massless Dirac electrons entering a random scalar potential barrier normally [8,[30][31][32], where the scalar potential plays a similar role as ò P .
Next, we consider the situation where the longitudinal component ò P is very strongly disordered, while ò ⊥ is nonrandom. Then, in the numerator of equation (17), ò P dominates the sin 2 θ term with high probability and we get cos ) , which is nonrandom. Therefore the localization length has to diverge for all θ as g P goes to infinity, if ò ⊥ is positive and nonrandom. This implies that the dependence of ξ on g P is non-monotonic: as g P increases from zero to infinity, ξ initially decreases, then increases to infinity. This behavior is again similar to that obtained for massless Dirac electrons in a one-dimensional random scalar potential [8,33]. Our numerical method described in the previous section relies on the assumption that the disorder in   is sufficiently weak, and therefore it cannot be used to study the limit g  ¥  . However, we can use a method based on the formula of differentiation derived by Shapiro and Loginov [34] to study Anderson localization for arbitrarily strong disorder. This approach is beyond the scope of this paper and will be presented in a future publication.
Finally, we consider the case where ò P is equal to 1. Then the expression for the impedance reduces to 2 1  h =^-, which is independent of θ for any functional form of ò ⊥ . In this case, if there is only one scattering interface, the transmission is independent of the incident angle. However, if there are more than one interfaces such as in a uniform slab of finite thickness or in the case with inhomogeneous z ^( ), then the interference of multiply scattered waves will occur. This effect depends on p ( k cos q = ) and ò ⊥ (z), therefore the transmission and other characteristics depend on θ in general. As an example, we show in figure 1 the localization length as a function of θ calculated using the invariant imbedding method when ò P =1, a=2 and g ⊥ =0.01, 1, 100. We remind that when only ò ⊥ is random, our method can be applied to any large value of g ⊥ and the results shown here are exact. When g ⊥ is smaller than 1, ξ is very accurately given by which is a special case of equation (29) to be derived in the next section. We notice that the localization length has a strong θ dependence and diverges as θ approaches 90°. This divergence was pointed out previously by Jordan el al, who studied an alternating isotropic-uniaxial random layered medium using the transfer matrix method [17]. However, the behavior of their data shown in figure 4(a) of [17] is markedly different from ours in that, in their case, as θ increases from zero, ξ remains almost constant up to 60 q~, and then increases sharply to infinity as θ approaches 90°. Whether this difference is due to the difference in the models used or some other reason remains to be investigated.

Argument based on the Fresnel formula
Equivalently, we can derive the existence condition of the BA using the Fresnel formula. We consider our medium as consisting of a large number of very thin layers. The reflection coefficient between two neighboring layers is written as = in uniaxial media. We suppose that the wave is delocalized at an incident angle θ B . In order for delocalization to occur, the random variation of a and b should not cause any reflection, and therefore we have the no-reflection condition, p a p a = ¢ ¢. We write a¢ and b¢ as a a a d ¢ = + and b b b d ¢ = + , with δ a and δ b as small quantities. Substituting these into p a p a = ¢ ¢ and using the Taylor expansion, we obtain which implies that δ b has to be proportional to δ a . Therefore, only Model II can show the BA. If we define δ b=fδ a, the condition for the BA becomes identical to equation (21). The same conclusions can be deduced from the expressions for the localization length, equations (9) and (15). In one dimension, waves are localized in the presence of even an infinitesimally weak randomness, except for in some special cases. The fact that a p wave is delocalized at θ=θ B implies that disorder does not play any role in the wave propagation process and the reflection coefficient r is the same as the value in the absence of disorder, r 0 , given by   1 t a n , 2 7 from which we conclude that only Model II with f>0 can display the BA.

Analytical expressions for the localization length in the weak disorder regime
Starting from equations (9), (12), (13) and (15), it is possible to derive analytical expressions for the localization length in the weak disorder limit. We write r as r=r 0 +δ r. From numerical calculations, we have verified that r d á ñ and r 2 d á ñ ( ) are of the first order in disorder, while r 3 d á ñ ( ) is of the second order, except at incident angles close to the critical angle for total internal reflection. From this consideration, we substitute Z r r Z r r r r Z r r r r r , 2 , into equation (12) in the l  ¥ limit when n=1 and 2 and obtain two coupled equations for r d á ñ and r 2 d á ñ ( ) . We solve them analytically and substitute the results into equation (9) to the leading order in the disorder parameters. The final expression for the localization length for Model I is 3 2 and Θ is the step function, for x>0 and 0 for x<0. Similarly, we obtain the localization length for Model II as We have found numerically that both of these equations are quite accurate when the disorder parameters are sufficiently small, except near the region where w=0. In the isotropic case with f=1 and b=a, the second term of equation (

Numerical results
In figure 2, we show the normalized localization length, kξ, as a function of the incident angle for Model I, when a=2 and b=±1.5. We note that the case with a>0 and b<0 corresponds to a type I hyperbolic medium [35]. We consider three cases, where only ò ⊥ is random, only ò P is random and both ò ⊥ and ò P are random. In the first case, ξ increases (decreases) monotonically as θ increases when b=1.5 (b=−1.5), while, in the second case, it diverges at θ=0 and decreases monotonically as θ increases for both b=±1.5. The third case is a combination of the first two cases. These behaviors can be readily understood from the form of the function U defined by which, in the equivalent Schrödinger equation, plays the role of V(z)/E, where V(z) is the potential and E is the energy of an incident quantum particle. In the case of figure 2(a), we find that the strength of the δò ⊥ term decreases (increases) monotonically as θ increases when b=1.5 (b=−1.5), in consistence with the behavior of ξ. We notice that if b=1, the δò ⊥ term will vanish and ξ will diverge, as θ approaches 90°. This case corresponds to that shown in figure 4(a) of [17], where the longitudinal component of the refractive index is uniform and matched to that of the surrounding medium. In the case of figure 2(b), the strength of the δò P term increases from zero monotonically as θ increases from zero, regardless of the sign of b, which is again in consistence with the behavior of ξ. When only ò P is random, all normally incident waves are delocalized. We find that localization is dominated by the randomness of ò ⊥ at small incident angles. The nonmonotonic behavior of ξ shown in figure 2(c) when a=2 and b=1.5 is qualitatively similar to that shown in figure 2(b) of [17]. We point out that the system called mixed stack in [17] corresponds to Model I and will not show the BA, while that called binary stack can show it.
In figure 3, we plot kξ versus θ for Model II, when a=2, b=±1.5 and g ⊥ =0.01, for various values of f. we get the BA if f>0.375 ( f>1.875), in a perfect agreement with the numerical results. The entire curves as well as the values of θ B agree precisely with equation (31). The BA has also been observed in figure 1(a) of [17], where it is easy to see that the random variations of ò ⊥ and ò P are directly proportional to each other with a positive ratio, in consistence with our results. Next, we consider the situation where 0<b<1. There exists a critical angle of total reflection, θ c , given by b sin c q = . Then the BA can occur for both a>1 and a<−1 cases. In figure 4, we plot kξ versus θ for Model II, when a=±2, b=0.5 and g ⊥ =0.01, for various values of f. In the case of figure 4(a), the BA is possible for Figure 2. Normalized localization length, kξ, versus incident angle, θ, for p waves in Model I, when a=2, b=±1.5 and (a) g ⊥ =0.01, g P =0, (b) g ⊥ =0, g P =0.01, and (c) g ⊥ =g P =0.01. The numerical results obtained using the invariant imbedding method are compared with those obtained from the analytical formula, equation (29), which are designated by dots. any value of f>0. In the case of figure 4(b), it is possible only if 0<f<0.125. We note that θ B <θ c =45°if a>0, while θ B >θ c if a<0. Interestingly, in the corresponding non-disordered case with g ⊥ =0, the ordinary Brewster angle, θ b , which is given by b a ab sin 1 1 , exists only when a is negative. When a=−2 and b=0.5, θ b is equal to 60°(> θ c ) and has no direct relationship to θ B . When a=2 and b=0.5, no Brewster effect occurs in the clean case, still the BA can occur at an angle smaller than θ c in the random case.
In table 1, we make a comparison between the results of this work and those of [17]. We remind that our model with δ-correlated disorder is substantially different from the multilayer model of [17] and only qualitative comparisons can be made. One of the biggest differences is that the ordinary Brewster angle θ b is the same as the angle θ B where the BA would occur in [17], while those two angles are unrelated in our work.

Conclusion
In conclusion, we have studied Anderson localization and the BA of electromagnetic waves in random anisotropic media theoretically. We have presented a unique perspective on the BA that it is a phenomenon occurring when the effective wave impedance is uniform and non-random, which is possible only in weaklydisordered media. We have also argued that the Brewster effect occurs when the wave impedance is completely matched throughout the space. We have derived the existence condition for the BA and analytical expressions for the localization length and elucidated several interesting physical aspects. Our results can provide valuable insights in understanding the unique properties of some biological reflectors and designing novel photonic devices based on anisotropic media [36].