Homogenization of Maxwell’s equations in a layered system beyond the static approximation

The propagation of electromagnetic waves through a disordered layered system is considered in the paradigm of the homogenization of Maxwell’s equations. Although the accuracy of the effective dielectric permittivity and/or magnetic permeability is still unclear outside the static approximation, we show that the effective wave vector can be correctly introduced even in high-frequency cases. It is demonstrated that both the real and imaginary parts of the effective wave vector are self-averaging quantities connected by the Kramers–Kronig relations. We provide a unified approach to describe the propagation and localization of electromagnetic waves in terms of the effective wave vector. We show that the effective wave vector plays the same role in describing composite materials in electrodynamics as the effective dielectric permittivity does in statics.

The electrodynamic properties of artificial composite materials depend on their structure and the component compounds and can differ substantially from the properties of homogeneous materials 1 . Nevertheless, the description of metamaterials as homogeneous by means of their effective parameters (dielectric permittivity and magnetic permeability, chirality coefficients, etc.) remains a relevant problem in electrodynamics. The problem of replacing an inhomogeneous composite system with a homogeneous material with effective parameters is called the homogenization problem (detailed formulation of the problem is given in the literature 2 ). It is assumed that the effective parameters of a homogeneous system provide the same scattered field as that of an inhomogeneous system (see Fig. 1). That is, the effective parameters allow the description of light scattering by a composite system without considering the inhomogeneous structure of the system.
The accuracy of the homogenization problem is under active study in connection with the characteristic scales 3 -system length L, inhomogeneity size ξ and wavelength . As a rule, it is assumed that the linear size L of the composite system is sufficiently large and the characteristic size of the inhomogeneities is sufficiently small that the composite can be considered a homogeneous material in the macroscopic sense. In addition, if we have a procedure that provides the effective parameters for a finite sample of size L, then as L increases, these parameters should tend to some constant values (this property is called "self-averaging").
Notably, the homogenization problem in electrostatics differs substantially from the one in electrodynamics. In the static case (the rigorous theory of the homogenization problem is constructed only for the static case), Maxwell's equations are divided into two independent groups of equations-magnetic and electric. Therefore, the description of inhomogeneous media by means of effective dielectric permittivity ε eff and magnetic permeability µ eff is quite natural. The microscopic description of each of the fields (electric and magnetic) reduces to a dimensionless Laplace equation, and a scaling approach can be applied 4 . Currently, the greatest progress has been made in homogenizing Maxwell's equations for the static case. Namely, the scaling theory of G-convergence 4 and spectral theory [5][6][7][8] have been developed.
In contrast to statics, the electric and magnetic fields are related to each other in the case of dynamics. In particular, the amplitudes of the electric and magnetic fields of the plane wave are connected by admittance Y. In addition, the homogenization problem in dynamics is multiscale-there are additional scales associated with the equations, namely, the set of local wavelengths / √ εµ , where = 2π/k 0 is the wavelength in vacuum and k 0 = ω/c is the wavenumber. The presence of multiple scales leads to difficulties in the transition from a finite system to an infinite system 9 , and the results of the G-convergence theory cannot be applied in the general case. www.nature.com/scientificreports/ Thus, two questions arise when the theory of homogenization for dynamics is being constructed: which parameters should be chosen as effective ( ε eff , µ eff or n eff , Y eff or some others), and how should the homogenization procedure be implemented 2,10-13 ? For clarity, let us distinguish two separate approximations widely used in homogenization theory: the quasistatic approximation and the long wavelength approximation. The quasistatic approximation considers the homogenization of Maxwell's equations in electrostatics. However, the static response (e.g., polarization) of each inhomogeneity is replaced by the response at non-zero frequency. For example, the quasistatic method adopts the solution of the static homogenization problem for a one-dimensional system, which has been known for more than a century. However, the approach is invalid with respect to metamaterial description: resonant responses of inhomogeneities (Mie resonance, split-ring resonances, etc.) exist in the case of the temporal differential equation only. That is, the radiation illuminated by inhomogeneities should be calculated in the framework of a time-dependent problem. The second approximation is the long wavelength approximation, which considers the solution of Maxwell's equations in terms of ξ/ expansion, where ξ is the inhomogeneity size and is the wavelength. However, ab initio simulations of light scattering in layered systems 14 indicate a disagreement between the numeric experiment and the predictions of long wavelength approximation theory in the first order of ξ/ .
One of the approaches to the homogenization problem in the long wavelength approximation was proposed for periodic systems in the works of Wood et al. 15,16 and was further developed in [17][18][19] . The basic idea is to determine the effective dielectric permittivity of an inhomogeneous periodic dielectric medium by means of a Bloch wave vector k Bl = k 0 √ ε eff . It is implied that in the absence of magnetic components, we can assume k Bl = k 0 √ ε eff . In essence, this approach is equivalent to the fact that we transfer the independence of the electric and magnetic fields from statics to dynamics (that is, the effective permittivity and permeability are introduced independently by the homogenization procedure). Specifically, the impedance is implied to be equal to the inverse refractive index. This is a common case in the optics of natural media. However, this approach inaccurately describes metamaterials with artificial magnetism.
An alternative approach proposed by Rytov and Levin 20-26 for layered systems was generalized to two-dimensional and three-dimensional cases [27][28][29][30][31][32][33] . The basic idea is the averaging of the fields over the unit cell (or over its surface 13 ) of the periodic medium.
Two quantities are defined separately for the periodic system within Rytov's approach: the effective wave vector k Ryt eff = ε The effective wave vector is set equal to the Bloch wave vector k Bl . The effective admittance is obtained by field averaging over the unit cell volume: Y Ryt eff = �H�/�E� . However, this approach leads to a somewhat dubious result: either ε eff or µ eff has a negative imaginary part. Moreover, the introduced effective parameters fail validation by direct calculation of wave propagation through the finite system. The parameters calculated from scattering data by means of the R-T retrieval approach 34-36 oscillate with the system length, i.e., ε eff and µ eff are not self-averaging quantities 14,37 . The introduction of the parameters ε eff and µ eff might be incorrect for layered systems outside the static approximation [37][38][39] .
In summary, the theory of the homogenization problem beyond the static case is still under intense investigation.
Recently, we have shown 40 that the effective refractive index of a periodic system self-averages not only in the long wavelength approximation but also in high-frequency cases. However, periodic systems differ from disordered systems. Specifically, Anderson localization of light occurs in the latter.
In this paper we consider the normal incidence of light on the layered structure. The effective wave vector has only one non-zero component, which is parallel to the system axis. We refer to that component as the wave vector. We focus our attention on the properties of the effective wave vector k eff in the general case of disordered systems. We are interested in the self-averaging property of the effective wave vector defined by the transmission approach 40 . The self-averaging property of k eff is critical for the homogenization problem. The studies are performed by the numerical computation in different frequency domains and even in the high-frequency case. In addition, we generalize the Kramers-Kronig-like relations known for the periodic systems 40 to the case of the disordered systems, which is important for establishing the analytical relation between the localization length and the real part of the effective wave vector.

Results
Self-averaging of the effective wave vector. We study wave propagation through a disordered system composed of dielectric layers with different dielectric permittivities. For simplicity, we assume that all the layers have the same thickness.
The standard approach to the localization and propagation of light considers a disordered system immersed in a vacuum. In this case, the transmission and reflection coefficients and other scattering data are affected by www.nature.com/scientificreports/ both the Anderson localization 41 and the effect of reflection from boundaries between the vacuum and the disordered system. To separate the localization effect from the scattering at the boundaries, we immerse the system into a medium with an averaged dielectric permittivity (the averaging is taken over the permittivity probability distribution). The transmission coefficient of a one-dimensional disordered system can be represented in the following form: The propagation of an electromagnetic wave through the inhomogeneous layer can be viewed as a three-stage process: transmission of the external medium-layer boundary, propagation through the volume of the layer, and transmission of the layer-external medium boundary. The homogenization problem reduces to two independent problems from that perspective: determining the effective wave vector and considering the boundary conditions. To solve the problem of boundary conditions, regularization of the effective impedance by means of additional surface currents 42 or an additional effective layer at the boundary 43 was previously proposed. The problem of boundary conditions is also considered in the profound variation of Rytov's approach 13 , where the finite samples are studied. In the latter case, the field averaging is conducted over the unit cell, which is far from the system boundaries, to avoid boundary effects. In this paper, we shift our attention away from the boundary problem and focus on the effective wave vector.
Notably, the contribution of the boundaries of a system to the transmission coefficient does not increase (or increases more slowly than exponentially) as the system thickness increases. At the same time, propagation through the layer is characterized by an effective wave vector, whose contribution to the phase or exponential decay (associated with absorption or scattering of light) of the transmission coefficient increases linearly with the thickness L of the system. Thus, the reasonable approximation used in our discussion neglects the surface effects for large L (that is, the phase is much larger than π , i.e., L ≫ ). Therefore, the imaginary part of the logarithm of the transmission coefficient iIm(lnt) is a phase accumulated along the thickness L of the system. The immersion of the sample in a medium with an averaged permittivity reduces the reflections from the boundaries of the sample. It is helpful in the case of long wavelengths where the reduction of boundary contributions allows consideration of smaller systems (compared to the vacuum case).
Note that the phase φ of the transmission coefficient (1) is a multivalued function. The single branch of the function should be chosen for later discussion. We consider two procedures of phase restoration: "by frequency" and "by propagation". The phase restoration by frequency (by propagation) is based on the assumption that the phase φ is a continuous and monotonically increasing function of the frequency (coordinate). Both these procedures lead to the same value of the phase and are used in this paper. For convenience, we refer to the phase restored by frequency (or by propagation) as phase, and the phase reduced to the region [0, 2π] , as reduced phase.
The real part of ln|t|/L describes the decay of the wave amplitude, which is determined by the imaginary part of the wave vector. Thus, the effective wave vector for a thick system can be defined as follows 40 : These formulas can be combined into The definition of k eff by means of the transmission coefficient is valid for periodic systems as confirmed by numerical calculations 14,40 .
By definition (3), the imaginary part of the wave vector converges to the inverse length of the localization: where L loc is the localization length. The inverse length of localization γ shows the rate of exponential decay of localized state in the disordered system. The quantity γ is a self-averaging quantity called the Lyapunov exponent in the theory of Anderson localization 44 . That is, the imaginary part of the wave vector self-averages as the system length increases. That is also confirmed by the numerical calculations (see Fig. 2). The existence of the Anderson localization scale-the Anderson localization length-is equivalent to the self-averaging of Imk eff . Let us now consider the real part of the effective wave vector. Calculations show that the real part demonstrates similar behavior (Fig. 3): its distribution becomes narrower as the system length increases.
As seen from the graphs in Figs. 2 and 3 for the real and imaginary parts of the effective wave vector, selfaveraging occurs as the system thickness increases, i.e., the distribution function becomes narrower. The rate of self-averaging can be estimated. The variance of the distribution function as a thickness function is shown in Fig. 4. The dependency of the variance on the system length follows a power law σ k eff ∼ L −α . The exponent obtained from the calculations is α ≈ 0.5 . Thus, the effective wave vector is expected to reach a certain value for a sufficiently long system. In other words, both the imaginary and the real parts of the effective wave vector self-average.

(2)
Rek eff = Im(lnt) L , www.nature.com/scientificreports/ Thus, we have shown numerically that k eff determined by formula (4) self-averages as the thickness of the system increases not only in the long wavelength approximation but even outside of this approximation. Moreover, convergence occurs for the high-frequency cases when the wavelength is on the order of the layer thickness or even significantly smaller. At the same time, the effective wave vector acquires an imaginary part due to wave localization in the system. Notably, the number of layers is important for self-averaging (quantity averaging over the few layers is unlikely). However, self-averaging occurs even for high-frequency cases for sufficiently large (composed of many layers) systems.   www.nature.com/scientificreports/ Now, let us compare definition (4) with other definitions of the effective wave vector. Numerical calculations for an ensemble of disordered systems show that in the long wavelength approximation the value of the effective wave vector calculated by formula (2) coincides with that obtained from R-T retrieval 34 . The value of k eff averaged over the ensemble in the long wavelength limit for both definitions are close to the function k st = √ �ε�k 0 (which corresponds to effective permittivity in the static case ε st = �ε� ), where the dielectric permittivity is averaged over the distribution. Therefore, we assumed in the calculations that the admittance of the external medium is Y av = √ �ε� . Utilizing an external medium with this value of admittance makes it possible to reduce reflections from the boundaries of the system and, thus, reduce the dependency of the transmission coefficient on the boundaries. At the same time, the introduction of an admittance differing from that of vacuum is not essential if the size of the system is substantially larger than the wavelength.
The definition of the effective wave vector was introduced earlier in the literature 45 . The multilayered system is considered as a unit cell of an infinite photonic crystal. According to the approach 45 the effective wave vector is equal to wave vector k Bl of the Bloch wave in the photonic crystal. The definition was applied 45 to analyze the Anderson localization in layered systems. Specifically, it was shown that k Bl gains the non-zero imaginary part for almost all frequencies if the system length increases. Moreover, the imaginary part of k Bl tends to the Lyapunov exponent (5) if the system length increases, i.e., �Imk Bl � −−−→ L→∞ �Imk t � for almost all frequencies, where k t is the effective wave vector (4). Let us now consider the convergence �Rek Bl � −−−→ L→∞ �Rek t �. Numerical calculations show that the difference (k t − k Bl ) between the two definitions of the effective wave vector decreases as the system size increases. Notably, the real part of k Bl frequently changes with the thickness of the system by jumps of size ∼ π/L . This occurs because the Bloch wave vector corresponds to the band gap for almost all frequencies for sufficiently long system 45 . Although the phase shift on the period of the considered photonic crystal changes abruptly, k Bl and k t tend to the same value.
Thus, the effective wave vector introduced by (4) corresponds to the definitions given in literature 34,45 . Qualitative arguments on the self-averaging of the effective wave vector. Following Sheng 46 , let us consider two large (consisting of many layers and with a total thickness considerably larger than the wavelength) pieces of one disordered system referred to as pieces A and B. The transmission coefficients of A and B are t A and t B , respectively. The system composed of these two pieces has transmission coefficient t AB , which equals (in order of magnitude) t AB ≈ t A t B . Now, adding one more piece, referred to as C, we assume that t ABC ≈ t AB t C ≈ t A t B t C , with logarithmic accuracy This relation is valid for both modules and phases but only if we choose the phases reconstructed by propagation, that is, the phase of the wave is supposed to be an increasing function of system thickness.
Since ln t i are independent and identically distributed, 1 N N i=1 ln t i is asymptotically Gaussian distributed in accordance with the law of large numbers. The variance of the distribution is σ ∼ 1/ √ N , i.e., D = σ 2 ∼ 1/L when all the layers of the system have the same thickness. Therefore, the quantities φ/L and γ are self-averaging; moreover, D γ ∼ D k ∼ 1/L , which corresponds to the numerical experiment (see Fig. 4).

Frequency dependence of the real and imaginary parts on the effective wave vector.
Consider the dispersion properties of the introduced effective wave vector. The graphs of the dependencies are shown in Fig. 5 for layered systems composed of two types of layers. The effective wave vector is close to the value k st = √ �ε�k 0 in the long wavelength limit. However, as the frequency increases, the dependency changes and becomes k eff ≈ � √ ε�k 0 . The imaginary part of k eff becomes 0 at some points because one of the layers becomes transparent at these frequencies, i.e., the transmission coefficient of the layer equals ±1 . The dispersion dependency of the imaginary part of k eff behaves as Imk eff ∼ −2 in the long wavelength regime (see the inset in Fig. 5b), which is in accordance with the known behavior of the Lyapunov exponent.
Kramers-Kronig-like relations for the effective wave vector. Since both k = Rek eff and γ = Imk eff are self-averaging quantities, a possible connection between them and their physical meaning are of interest.
We examine this question by means of the transfer matrix method 47 (see "Methods" section). Clearly, the transmission coefficient of the layered system is an analytic function of the frequency 40 . The T-matrix of the whole system can be written in terms of the transmission and reflection coefficients as follows where t is the transmission coefficient (the transmission coefficient is independent of the direction), r L and r R are reflection coefficients for the left and right sides of the system, respectively. On the other hand, the T-matrix of the whole system is the product of the T-matrices of the homogeneous layers forming the system. All elements of these T-matrices are analytic functions of the frequency. Thus, the transmission coefficient, being a fractionally rational function of these elements, is also an analytic function of the frequency. Consequently, the function should also be analytical. Therefore, the following relations are valid: www.nature.com/scientificreports/ The Kramers-Kronig-like relations are well-known properties of optical parameters connected to the causality principle 48 . The Kramers-Kronig-like relations (9) and (10) show that the effective wave vector being the effective parameter of the disordered system still satisfies the essential principle of causality. The relations (9) and (10) are the important ones to reveal the physical meaning of the introduced effective wave vector. With respect to the definition of the Lyapunov exponent γ (ω) , the second relation can be rewritten as (we make the notation k = Rek eff ) Taking into account parity k(−ω) = −k(ω) (hence, dk/dω is an even function of frequency) and utilizing the relation between the wave vector and the density of states ( ρ(ω) = dk/dω ), one can obtain Performing integration over the frequency, one obtains the Jones-Herbert-Thouless formula 49,50 (9) dRek eff (ω) www.nature.com/scientificreports/ Thus, the frequency derivative of the real part of the effective wave vector corresponds to the density of states, which confirms that k eff has physical grounds to be considered an effective wave vector. phase randomization. It was shown in the preceding sections that the effective wave vector self-averages.
Let us now consider in detail the behavior of phase φ (1), which is directly related to the effective wave vector for each system realization as φ = −i ln (t/|t|) = Im ln t = k eff L . Although k eff self-averages and is consequently a deterministic quantity, the phase randomizes; more precisely, the reduced phase, which is a phase reduced to an interval [0, 2π] , randomizes. Phase randomization is used as the basis for developing a random T-matrix approach [51][52][53][54] . The apparent contradiction can easily be explained by noting that the variance of the wave vector tends to zero as D k = σ 2 k ∼ 1/L and the variance of the phase increases as D φ = L 2 σ 2 k ∼ L . Since the distribution of k eff is similar to the Gaussian distribution (see Fig. 3a), the distribution for the total phase is expected to be similar to the Gaussian one, but as the system thickness increases, the distribution shifts and becomes wider.
The numerical calculations confirm this behavior (Fig. 6). The distribution shifts to the right proportionally to the system size and broadens. Therefore, after the reduction to the interval [0, 2π] , the distribution of the reduced phase should be close to a constant for a sufficiently large system (Fig. 6b).
However, even for a thick system, the calculations show that the distribution function differs from a constant due to the presence of an additional correlation scale-the width of the layer.

Discussion
In this paper, we considered the propagation of an electromagnetic wave through a disordered layered system. The majority of approaches to the homogenization problem propose introducing corrections (which are proportional to the relation ξ/ between the inhomogeneity size and the wavelength) to the static effective dielectric permittivity and/or magnetic permeability. For example, the second-order theory describes chirality, magnetic dipoles www.nature.com/scientificreports/ and electric quadrupoles 55 . However, the theory experiences breakdown as the parameter k 0 d increases, and its validity depends on the system length 55 . These approaches have a limited field of application and are inaccurate even in the long wavelength approximation [37][38][39] : the effective parameters ε eff and µ eff vary with the number of layers (that is, they oscillate), which leads to breakdown of the approaches. However, in this work, we showed that the effective wave vector k eff tends to a constant value for different relations as the number of layers increases. In other words, it self-averages and is thus a macroscopic quantity. The self-averaging of k eff is observed not only in the long wavelength approximation ξ/ ≪ 1 but even in the high-frequency cases of ξ/ ∼ 1 and ξ/ ≫ 1.
The introduction of disorder into periodic systems substantially changes the spectral characteristics, such as the density of states and the transmission and reflection coefficients 56,57 . Moreover, the disorder in photonic structures leads to Anderson localization 41 . Nevertheless, the complex optical properties of disordered systems fit the effective wave vector description. The imaginary part of the proposed effective wave vector equals the Lyapunov exponent (inverse localization length) in the limit of infinite system length L → ∞.
To clarify the physical meaning of the real part of the introduced effective wave vector, we consider the frequency derivative of the effective wave vector. We show that this quantity is an analytical function of frequency (because it can be directly expressed in terms of the propagation coefficient) and thus satisfies the Kramers-Kronig-like relations. In this work, it was proven that these relations lead to the Jones-Herbert-Thouless relation, connecting the spectral density and the localization length of eigenstates, i.e., the real part of the introduced effective wave vector gives us the correct value of the density of states.
We show that k eff acts as the self-averaging parameter in electrodynamics instead of ε eff in statics. Homogenization of the effective wave vector enables a unified description of electromagnetic wave propagation and localization.

Methods
The numerical results given in the article are obtained by computation of transmission coefficient of wave through the multilayered system. The propagation of plane waves through layered system is studied via the transfer matrix (T-matrix) method 47 . The T-matrix T of a layered system relates the amplitudes of waves propagating in different directions. Specifically, the following relation exists for the case considered in the article where r and t are the reflection and transmission coefficients, correspondingly.
The T-matrix T j of a single layer placed into an external medium (to unify the description, it is convenient to place the layer between two adjacent imaginary layers of zero width) is as follows: where Y ext is the external medium admittance, ρ j = n j k 0 d is the optical thickness and Y j = n j is the admittnace in the case of normal incidence. The T-matrix of a system composed of N layers is obtained by multiplying the T-matrices of the layers in reverse order Equation (14) enables calculation of the transmission coefficient. Notably, the calculation of the transmission coefficient phase is performed according to the "by propagation" phase recovery procedure discussed in the text. The straightforward application of Eq. (4) yields the value of the effective wave vector.
The method allows for a mathematically precise solution of Maxwell's equations in one-dimensional geometry. Thus, the only inaccuracies in the calculations are the numerical errors due to the number precision limit. The smallness of the inaccuracies was monitored in the calculations by evaluating T-matrix invariants: |T| = 1 , and in the absence of absorption, T 12 = T * 21 , T 11 = T * 22 . The transfer matrix approach was implemented in C++ programming language.