Exact localization length for s-polarized electromagnetic waves incident at the critical angle on a randomly-stratified dielectric medium

The interplay between Anderson localization and total internal reflection of electromagnetic waves incident near the critical angle on randomly-stratified dielectric media is investigated theoretically. Using an exact analytical formula for the localization length for the Schr\"odinger equation with a Gaussian $\delta$-correlated random potential in one dimension, we show that when the incident angle is equal to the critical angle, the localization length for an incident $s$ wave of wavelength $\lambda$ is directly proportional to $\lambda^{4/3}$ throughout the entire range of the wavelength, for any value of the disorder strength. This result is different from that of a recent study reporting that the localization length at the critical incident angle for a binary multilayer system with random thickness variations is proportional to $\lambda$ in the large $\lambda$ region. We also discuss the characteristic behaviors of the localization length or the tunneling decay length for all other incident angles. Our results are confirmed by an independent numerical calculation based on the invariant imbedding method.


Introduction
Despite of the intensive studies for more than five decades, Anderson localization of quantum particles and classical waves continues to fascinate many researchers [1][2][3][4]. Recent studies have identified new systems displaying Anderson localization and revealed new aspects of the phenomenon [5][6][7][8]. Localization of noninteracting quantum particles in a random potential in one dimension is an old problem that has been studied extensively for a long time [9,10]. Among many theoretical results, we revisit an exact analytical formula for the localization length in the presence of a Gaussian δ-correlated random potential derived long ago by several authors [11][12][13][14]. This formula can be applied to the cases where both the strength of the random potential and the ratio between the average value of the potential and the particle energy are arbitrary.
Due to the similarity between the Schrödinger equation and the electromagnetic wave equation for s-polarized waves, the analytical formula for the localization length mentioned above can also be applied to s waves incident obliquely on randomly-stratified (or randomly-layered) dielectric media. When the disorder-averaged value of the dielectric permittivity is smaller than the permittivity in the incident region, a modified total internal reflection phenomenon occurs near and above the critical angle [15,16]. In this paper, we are especially interested in the interplay between Anderson localization and total internal reflection near the critical angle, which has drawn some renewed theoretical and experimental interest recently [17,18]. In particular, we show that when the incident angle is precisely equal to the critical angle, the localization length for an incident s wave of wavelength λ is directly proportional to λ 4/3 throughout the entire range of the wavelength. This result is different from that of a recent theoretical study based on the transfer matrix method reporting that the localization length at the critical incident angle for a binary multilayer system with random layer thicknesses of deeply subwavelength scale drawn from a uniform distribution, which has short-range disorder correlations embedded in it, is proportional to λ in the large λ region [17]. In addition, we discuss the agreement and the disagreement between our results and those of two previous works [16,19]. The possible effects of disorder correlations on localization phenomena are also discussed. Our results are confirmed by an independent numerical calculation based on the invariant imbedding method.

Model
We are interested in the propagation and the Anderson localization of an s-polarized plane electromagnetic wave of frequency ω and vacuum wave number k 0 (= ω/c) in isotropic random dielectric media. The wave is incident obliquely on a stratified random medium, where the dielectric permittivity ǫ varies randomly only in the z direction. We assume that the random medium lies in 0 ≤ z ≤ L and the wave propagates in the xz plane. In the s wave case, then, the complex amplitude of the electric field, E, satisfies where q is the x component of the wave vector. We assume that the wave is incident from the region where ǫ = ǫ 1 and z > L and transmitted to the region where ǫ = ǫ 1 and z < 0. When θ is the angle of incidence, the quantity q is equal to k sin θ, where k = √ ǫ 1 k 0 . In the inhomogeneous slab located in 0 ≤ z ≤ L, ǫ is given by where ǫ is the disorder-averaged value of ǫ and δǫ(z) is a δ-correlated Gaussian random function with zero average. We assume that ǫ is independent of z. It is convenient to introduce a normalized variableǫ defined bỹ where The random variable δǫ satisfies where g 0 is a parameter measuring the strength of disorder and has the dimension of a length. With some simple manipulation, we can rewrite the wave equation in the form This equation has precisely the same form as the Schrödinger equation for a quantum-mechanical particle in a random potential in one dimension where E (= 2 κ 2 /2m 0 ) is the energy of the incoming particle of mass m 0 and V(z) is the potential. We assume that the potential is equal to 0 for z < 0 and z > L and V(z) = V 0 + δV(z) for 0 ≤ z ≤ L, where V 0 is a real constant and δV(z) is a Gaussian random function with zero mean and a white-noise spectrum: The correspondences between Eq. (6) and Eq. (7) are summarized in Table 1. s wave equation Schrödinger equation The tunneling of a quantum-mechanical particle occurs when the particle energy E (> 0) is smaller than V 0 . This condition corresponds to −a > cos 2 θ, or equivalently, sin θ > ǫ /ǫ 1 , which gives the usual expression for the critical angle, sin θ c = ǫ /ǫ 1 , when 0 < ǫ < ǫ 1 .

Lyapunov exponent
We consider the transmission of a quantum particle through the random potential V(z). We assume that the particle, which is described by a plane wave of unit magnitude ψ(z) = exp [iκ(L − z)], is incident onto the barrier from the region where z > L and transmitted to the region where z < 0. The complex reflection and transmission coefficients, r = r(L) and t = t(L), are defined by the wave functions outside the random region: We are interested in the Lyapunov exponent γ, which can be defined by in terms of the transmittance T (= |t| 2 ). When E is greater than V 0 , γ can be related to the the localization length ξ by In the tunneling situation where E < V 0 , ξ is more appropriately called as the tunneling decay length. The exact analytical expression for γ has been derived previously by several authors and takes the form [14] γ =D 1/3 F(X), and the function F(X) is expressed in terms of the Airy functions Ai and Bi and their derivatives (denoted by primes) as Using this and the correspondences shown in Table 1, we derive the analytical expression of ξ for s waves propagating obliquely in randomly-stratified dielectric media, which takes the form When ǫ is smaller than ǫ 1 , there exists a critical angle θ c , at which the argument of F in Eq. (15) vanishes. Then we find that for s waves incident at the critical angle, ξ c [= ξ(θ = θ c )] is precisely given by in the entire range of frequencies. We can rewrite this expression in terms of the wavelength in the incident region, λ (= 2π/k), as where we have inserted the numerical value of F(0) (≈ 0.364506). In Sec. 5, we will introduce a dimensionless disorder parameter g (= kg 0 /2), in terms of which we have for all values of g (> 0). We emphasize that not only the exponent but also the proportionality constant is universal in this equation.
The pure power law of Eq. (16) is due to the use of the δ-correlated model. If the disorder has a finite correlation length, then the validity of this power law will be restricted to a low-frequency or weak-disorder regime. In the context of the Schrödinger equation, Eq. (7), the analogue of this law has been shown to hold over the entire energy range for a δ-correlated random potential [12]. In the tight-binding model, where the disorder has a finite correlation length, it has been shown that the power law only describes the weak-disorder regime [13].
When θ is smaller than θ c and g is sufficiently small, the argument of F in Eq. (15) approaches to large negative values. Then, from the property of F such that we obtain and in the weak disorder or small frequency limit. If we include the prefactor, this relation can be written as In Ref. 16, the localization of s waves in a short-range-correlated multilayer system with a Gaussian random distribution of the layer refractive indices has been studied using the transfer matrix method. It has been shown that the localization length at the critical angle, ξ c , is related to the width of the distribution of the refractive index, ζ, by in the long-wavelength region. It has also been shown that, when the wave is incident normally, ξ satisfies in the long-wavelength region. This result can be extended to all incident angles below θ c . The randomly fluctuating part of the refractive index, δn, is proportional to that of the dielectric permittivity, δǫ, and therefore we have δn(z)δn(z ′ ) ∝ δǫ (z)δǫ(z ′ ) . Since the coefficient of the δ function, 2g 0 , in Eq. (5) is proportional to the width of the corresponding Gaussian probability distribution, we conclude that our disorder parameter g 0 is proportional to ζ 2 . Then the scaling relations, Eqs. (23) and (24), are consistent with Eqs. (17) and (22), from which it follows that ξ c ∝ g −1/3 0 and ξ ∝ g −1 0 respectively. Since g 0 should appear as the dimensionless combination kg 0 and ξ has the dimension of a length, we also obtain ξ c ∝ λ 4/3 when θ = θ c and ξ ∝ λ 2 when θ < θ c . In the long-wavelength region, these results agree precisely with ours obtained using the δ-correlated model.

Influence of the correlation of disorder
The δ-correlated model studied here has no correlations, while the multilayer models used in Refs. 16 and 17 have short-range correlations embedded in them. The influence of correlations on localization has been a major topic of research in recent years [6,20]. It has been shown that long-range correlations generally have dramatic effects on localization [21-23], while shortrange correlations, except for in some special cases [24,25], have weaker effects on it.
The multilayer model of Ref. 17, which is a kind of periodic-on-average structure, has been designed specifically to study the influence of short-range correlations on localization experimentally [18], and its application to other wave systems is limited. On the other hand, the δ-correlated model is more generic and is valuable in that it can be applied to all wave systems to study their long-wavelength properties.
Properties of Gaussian random systems can be characterized by the spectral density, S(k), which is the spatial Fourier transform of the disorder correlation function such as δǫ (z)δǫ(z ′ ) . For the δ-correlated model given by Eq. (5), we obtain which is a constant independent of k. The feature that S(k) is finite in the k → ∞ limit and ∫ S(k)dk diverges is unphysical, and therefore the prediction of this model in the shortwavelength region should not be trusted.
The simplest kind of a short-range-correlated model, which can be treated analytically, is given by where l c is the disorder correlation length. In the limit where σ → ∞, l c → 0 and σ 2 l c → g 0 , this correlation function reduces precisely to Eq. (5). Its spectral density is given by In the long-wavelength limit (k → 0), this function approaches to the constant 2σ 2 l c . The power-law dependences of the localization length on λ in Eqs. (17) and (22)  Under the RG transformation [26], it can be shown that the k-dependent term in Eq. (27) is irrelevant in the macroscopic limit and does not contribute to determining the critical exponent. Therefore, in the long-wavelength regime, the critical exponents for the two models given by Eqs. (5) and (26) will be the same, while, in the short-wavelength regime, these models will behave differently.
Somewhat more general short-range-correlated models are expected to have the spectral density of the series form S(k) = ∞ n=0 a n k 2n .
Under the RG transformation, all terms with n > 0 can be shown to flow to zero and only the leading constant term a 0 is relevant. Therefore the critical exponents for this class of models will also be the same as those for the δ-correlated model. This may be the reason why the exponents for the model in Ref. 16 are the same as our exponents for both θ = θ c and θ < θ c . Using the RG theory terminology, these two parameter regions are described by two different RG fixed points. We note that δ-correlated models have been routinely used in a large number of studies on static critical phenomena in random systems [28,29] and also in those on dynamic critical phenomena driven by δ-correlated noise [30]. There can be exceptions to the argument given above. They are usually related to the presence of a hidden symmetry such as the gauge symmetry. In the literature, there are examples of low-dimensional wave systems, which include dimer type models [24] and multilayer models containing alternate positive index and negative index layers [25], where the long-wavelength scaling behavior is nontrivially affected by short-range correlations. These models may have a symmetry that constrains to change the long-wavelength behavior of the spectral density. More specifically, a possible scenario to allow a system with short-range correlations to behave differently from the δ-correlated model is that it has a special symmetry that forces the constant term a 0 to remain strictly zero under the RG transformation and S(k) contains a nonanalytic leading term of the form |k| β with 0 < β < 2. An argument of a similar kind has been used previously in an RG study of a system with a gauge symmetry [31]. In Ref. 17, the authors have reported the results of analytical and numerical analysis showing that ξ c ∝ λ at θ = θ c . Does this result suggest that the model in Ref. 17 has a relevant perturbation or a hidden symmetry that makes its universality class different from that of Ref. 16? We consider this to be a very interesting question which deserves further investigation.
Long-range correlations have much more dramatic effects on critical phenomena. For instance, the model in Ref. 21 has the spectral density scaling as S(k) ∝ k −p in the small k limit (p > 0). Under the RG transformation, this term is more relevant than the constant term and dominates the macroscopic critical behavior. The critical exponents will depend on the value of p continuously and a nontrivial phase transition can occur in such systems. It is a very interesting problem to study the influence of long-range correlations on the scaling behavior considered in this paper.

Invariant imbedding method
The reflection and transmission coefficients for s waves are defined similarly to Eq. (9). For an s wave of unit magnitudeẼ(x, z) = E(z) exp (iqx) = exp [ip(L − z) + iqx], where p = k cos θ, incident on the medium, they are defined bỹ Using the invariant imbedding method [32-36], we derive exact differential equations satisfied by r and t: which are integrated from l = 0 to l = L using the initial conditions r(0) = 0 and t(0) = 1. We can use Eq. (30) in calculating the disorder averages of various physical quantities consisting of r and t. We are mainly interested in calculating the localization length ξ of the wave defined by In order to obtain ξ, we need to compute ln T (l) in the l → ∞ limit. The nonrandom differential equation satisfied by ln T is obtained using Eq. (30) and Novikov's formula [37] and takes the form where we have used the definitions g = kg 0 /2 and Z n = r n (n = 1, 2). From this, we obtain Therefore the problem of calculating ξ reduces to that of calculating Z 1 and Z 2 in the l → ∞ limit.
It turns out that in order to obtain Z 1 and Z 2 , we need to compute Z n for all positive integers n. An infinite number of coupled nonrandom differential equations satisfied by Z n 's are obtained using Eq. (30) and Novikov's formula and have the form where Z n satisfies the initial condition Z n (l = 0) = 0. The magnitude of Z n decays (more rapidly for the larger value of g) as n increases. Based on this observation, we solve the infinite number of coupled equations using a simple truncation method [38]. In the l → ∞ limit, the left-hand side of Eq. (34) vanishes and we get a set of algebraic equations. This set of algebraic equations takes the form of a five-term linear recursion relation and can provide another way to derive Eq. (15). In Fig. 1, we show the inverse of the normalized localization length (or tunneling decay length if θ > θ c ) as a function of the dimensionless disorder parameter g in a log-log plot for various values of the incident angle θ, when ǫ /ǫ 1 = 0.75 and θ c = 60 • . The curves obtained from the analytical formula, Eq. (15), are compared with the numerical results obtained using the invariant imbedding method, which are designated by square dots. The agreements are perfect in all cases.

Numerical results
When θ is equal to θ c , we note that the curve (designated by the red color) is a straight line in the entire range of g. This is obvious from Eq. (18), from which we have As g (= kg 0 /2) increases to large numbers, or equivalently, as the disorder strength or the frequency increases to large values, the argument X of the function F in Eq. (15) approaches to zero. Then it follows that all curves converge to the straight line obtained when θ = θ c as g becomes large, as demonstrated clearly in Fig. 1. In Ref. 19, it has been claimed that the localization length for electromagnetic waves in one dimension in the high frequency limit is proportional to ω −2/3 , while it is proportional to ω −2 in the low frequency limit. Though the low frequency result is correct, the high frequency result differs clearly from our result that ξ ∝ ω −4/3 .
When θ is smaller than θ c and g is sufficiently small, all curves approach straight lines with the slope equal to one in the logarithmic variables. This well-known behavior is obtained from Eq. (20), which leads to the ξ ∝ λ 2 scaling in Eq. (22). We note that similar results were obtained for the short-range-correlated multilayer model studied in Ref. 16.
When θ is above the critical angle, the curves show a nonmonotonic behavior, namely the tunneling decay length first increases, reaches a maximum, and then decreases, as g increases from zero. This can be understood easily from the asymptotic expansion of F [14]: which leads to The leading effect of a weak disorder is to enhance the tunneling decay length. The related phenomenon of disorder-enhanced tunneling transmission has been studied previously by several authors [14,[39][40][41].   Fig. 2(b) of Ref. 17, which was obtained for λ = 1 µm. For the numerical comparison, we have chosen the disorder parameter g 0 = 0.00225 µm. With this choice of g 0 , we find that the agreement is excellent. We also find that the curve is smooth and shows no sharp feature across θ = θ c . This smooth behavior is a consequence of the fact that the scaling function F(X) in Eq. (15) is regular for all finite values of X. This implies that the distinction between the localization length at θ < θ c and the tunneling decay length at θ > θ c is blurred in the presence of disorder. In Fig. 3, we plot ξ versus wavelength for several values of θ near θ c and for a fixed value of g 0 and compare the result with Fig. 2(a) of Ref. 17. We could not find a value of g 0 that gives good agreements for all considered values of θ. In Fig. 3(a), we compare our results obtained using g 0 = 0.0003 µm with Fig. 2(a) of Ref. 17. When θ is equal to 58 • and 59 • , the agreement is pretty good and when θ is 59.5 • , it is not too bad. For θ = 60 • , however, our result differs greatly from that of Ref. 17. In Fig. 3(b), we compare our results obtained using g 0 = 0.000072 µm with Fig. 2(a)  The fact that we have obtained very good agreements in some curves is interesting and may suggest that the two models have close similarities in these wavelength scales despite of the obvious differences on a short-range scale. The reason why we do not obtain good agreements in all curves using the same disorder parameter is not clear at this stage and we believe it is worthwhile to elucidate it. In Ref. 17, it has been proposed that the data points for θ = θ c = 60 • should follow a straight line. However, it is interesting to notice that, in this wavelength scale, our result shown in Fig. 3(b) demonstrates that the data can be represented very well by the λ 4/3 dependence. It will be very helpful to compare the results of this paper with the calculations of Ref. 17 using equivalent parameters on an equal ground in a future work.

Conclusion
In this paper, we have studied the interplay between Anderson localization and total internal reflection of electromagnetic waves incident near the critical angle on randomly-stratified dielectric media. In particular, we have shown that when the incident angle is equal to the critical angle, the localization length for an incident s wave of wavelength λ is directly proportional to λ 4/3 throughout the entire range of the wavelength. A further theoretical investigation of the scaling behavior for various kinds of disorder correlations and experimental studies of the results reported here on randomly-layered dielectric structures are highly desired.

Funding
National Research Foundation of Korea Grant funded by the Korean Government (NRF-2015R1A2A2A01003494).