Artificial dispersion via high-order homogenization: magnetoelectric coupling and magnetism from dielectric layers

We investigate a high-order homogenization (HOH) algorithm for periodic multi-layered stacks. The mathematical tool of choice is a transfer matrix method. Expressions for effective permeability, permittivity and magnetoelectric coupling are explored by frequency power expansions. On the physical side, this HOH uncovers a magnetoelectric coupling effect (odd-order approximation) and artificial magnetism (even-order approximation) in moderate contrast photonic crystals. Comparing the effective parameters' expressions of a stack with three layers against that of a stack with two layers, we note that the magnetoelectric coupling effect vanishes while the artificial magnetism can still be achieved in a centre-symmetric periodic structure. Furthermore, we numerically check the effective parameters through the dispersion law and transmission property of a stack with two dielectric layers against that of an effective bianisotropic medium: they are in good agreement throughout the low-frequency (acoustic) band until the first stop band, where the analyticity of the logarithm function of the transfer matrix () breaks down.


Introduction
There is a vast amount of literature on the homogenization of moderate contrast periodic structures with the classical effect of artificial anisotropy. A less well-known effect is artificial magnetism and low-frequency stop bands through averaging processes in high-contrast periodic structures [1][2][3][4][5].  In layman terms, O'Brien and Pendry observed in 2002 that rather than using the LC-resonance of conducting split ring resonators to achieve a negative permeability as an average quantity [1,6], one can design periodic structures displaying some Mie resonance, e.g. by considering a cubic array of dielectric spheres of high-refractive index [2]. These resonances give rise to the heavy-photon bands in photonic crystals that are responsible for a low-frequency stop band corresponding to a range of frequencies wherein the effective permeability takes extreme values [4]. The so-called high-contrast homogenization [5] predicts this effect, which is associated with the lack of a lower bound for a frequency-dependent effective parameter deduced from a spectral problem reminiscent of Helmholtz resonators in mechanics.
In this paper, we achieve such a magnetic activity without resorting to a high-contrast material. The route we propose is based upon a homogenization approach for high frequencies, i.e. when the period of a multi-layered structure approaches the wavelength of the applied field. The extension of classical homogenization theory [7][8][9] to high frequencies is of pressing importance for physicists working in the emerging field of photonic crystals and metamaterials [10], but applied mathematicians also show a keen interest in this topic [3,11,12], where the periodic structures at sub-wavelength scales (λ/10 to λ/6) [2,6] can clearly be regarded as almost homogeneous.
The tool of choice for our one-dimensional model is the transfer matrix method, which allows for analytical formulae as shown in §2, a high-order homogenization (HOH) method is proposed wherein the Baker-Campbell-Hausdorff formula (the BCH formula, an extension of Sophus Lie theorem) was implemented; we stress that ideas contained therein can be extended to two-dimensional periodic structures like lamellar gratings, see [13], because our model extends those of Rytov [14]. Importantly, we not only achieve magnetic activity in moderate contrast dielectric structures, but also unveil some artificial magnetoelectric coupling (a particular case of bianisotropy) in §3, where a multi-layered stack with an alternation of two dielectric layers is considered. As proposed by Pendry, the bianisotropy is yet another route towards negative refraction [15]. However, it is noted that the artificial magnetoelectric coupling vanishes in a periodic stack with centre symmetry, where an extension of the HOH method applied to a stack with m (≥ 3) alternative layers is explored in §4. Furthermore, a correction factor is investigated in §5 to estimate the asymptotic error, which is proportional to 1/n p with p the approximation order. We then numerically check the effective parameters through the dispersion law and transmission property between the multi-layers and the effective medium in §6: a good agreement throughout the low frequency band up to the first stop band illustrates the equivalence between the structured dielectric medium and its effective bianisotropic counterpart. We finally consider a frequency power expansion of the transfer matrix, which is analytical in the whole complex plane in §7 and draw some concluding remarks in §8. Figure 1a shows a schematic of periodic multi-layered stack with a unit cell made of two homogeneous dielectric layers L 1 and L 2 of respective thicknesses h 1 and h 2 (h = h 1 + h 2 =h 1 /n + h 2 /n), where n stands for the number of the unit cells in the stack. The permittivities and the permeabilities of the two layers are denoted by ε 1 , ε 2 and μ 1 = μ 2 = μ 0 (with μ 0 the permeability of vacuum), respectively. It should be noted that the whole thickness of the stackh (= nh) is a constant and it is comparable with the wavelength λ of the applied field ash/λ ≈ 1. Starting from the case when n = 1, the stack is the most simple structure consisting of only two dielectric layers, in which case the effective medium theory [16][17][18] cannot be applied when h/λ ≈ 1. However, if we increase n to a large enough constant, then the system contains n times smaller unit cells, with a thickness much smaller than the wavelength of the applied field (e.g. h/λ 1). Hence, a homogeneous medium with permittivity ε eff , permeability μ eff and magnetoelectric coupling K eff shown in figure 1b can be assumed to behave as an effective medium for such a multi-layer. The approximation of the multi-layer by the effective medium will be all the more accurate than n is large, a fact which will be proved in the following sections.

Mathematical set-up of the problem
and the constitutive relations for non-magnetic isotropic dielectric media, with ε m the permittivity in the mth homogeneous layer located in the domain L m of R 3 . Then, a Fourier decomposition is introduced for both electric and magnetic fields aŝ with k 1 and k 2 the projections of wavevector k on x 1 and x 2 axes, respectively; the wavevector for an oblique plane wave is denoted by k = (k 1 , k 2 , k 3 ). Applying the decomposition (2.3) to (2.1) and (2.2), we derive an ordinary differential equation (involving a 4 × 4-matrix and a four-component column vector) [19] ∂F ∂x 3 (ω, k 1 , k 2 , whereF is a column vector containing the tangential components of the Fourier-transformed electric fieldÊ and magnetic fieldĤ, i.e. their components along x 1 and x 2 axes. Here, in order to simplify subsequent calculations, we define a new set of coordinates (x , x ⊥ , x 3 ): The component x is along the direction of wave vector k = (k 1 , k 2 ), x ⊥ is along k = (−k 2 , k 1 ), which is perpendicular to x . In other words, the new set of coordinates is a rotation of the previous coordinates around the x 3 -axis. The change of coordinates from (x 1 , x 2 , x 3 ) to (x , x ⊥ , x 3 ) can be expressed as Note that, thanks to the symmetry of the geometry, the parameters of the multi-layer are invariant under this transformation [20]. Hence, (2.4) can be recast in the new coordinate system as with the column vectorsF whereÊ ,Ĥ andÊ ⊥ ,Ĥ ⊥ are the components of the electric and magnetic fields along the x -axis and the x ⊥ -axis, respectively. Correspondingly, M m is a 4 × 4 matrix, the components of which can be expressed as where k 2 = k · k and k is the two-dimensional component of the three-dimensional wavevector k after projection in the plane (x 1 , x 2 ). The matrix M m (ω, k) is independent of x 3 in each homogeneous layer, the solution of (2.6) in the layer L m is simplyF (2.9) The exponential above is well defined as a power series of matrix M m (ω, k) and defines the transfer matrix in the mth layer of thickness h m , which relates the total tangential components of electric and magnetic fields at the two ends of the slab. Note that this definition only differs by a change of basis from an equivalent definition of the transfer matrix which relates the upward and the downward propagation fields at the two ends of the slab. As this power series has an infinite radius of convergence, the transfer matrix is analytical with respect to the three independent variables ω, k 1 and k 2 . For an arbitrary permittivity profile (with the classical assumption of upper and lower bounds of permittivity greater than ε 0 (the permittivity of vacuum) uniformly in position x and the number of layers n), analyticity is proved using a Dyson expansion [21].

(b) Main homogenization result
From the knowledge of the transfer matrix in multi-layered stacks, we deduce that the most general expression of the transfer matrix [22] can be obtained from the following homogenized constitutive equations and where ε eff and μ eff are tensors of rank 2, which represent, respectively, the (anisotropic) effective permittivity and permeability as follows: of the effective medium, the matrix J corresponds to 90 • rotation around the x 3 -axis, and K eff is the bianisotropic parameter measuring the magnetoelectric coupling effect Now, applying (2.3) to (2.1) and (2.11), we obtain with matrix where the 2 × 2 blocks are defined by These parameters are all unknowns at this stage which we derive from a homogenization algorithm. The transfer matrix is correspondingly, (2.18)

High-order homogenization algorithm for multi-layered stack
As we have derived the transfer matrix for the multi-layered stack and postulated its structure for the effective medium in the previous section, it follows that the description of the homogenization procedure shown in figure 1 can be expressed as This means that the two structures should exhibit the same transmission property, where M eff is the unknown to be calculated. Note that the left-hand side of (3.1) is a product of two exponential functions, which can be approximated by introducing the BCH formula (an extension of the Sophus Lie theorem, see [23] Furthermore, the expressions for the effective parameters in (2.12) and (2.13) can be derived by comparing the two matrices on the left-and right-hand sides of (3.4). First, we consider the zeroth-order approximation in (3.4), it yields M eff ≈ M 1 f 1 + M 2 f 2 with the filling fractions f 1 (= h 1 /h) and f 2 (= h 2 /h), respectively, and the effective parameters are They are identical to the effective parameters obtained in [17,18,[24][25][26] by classical homogenization: the effective permittivity and permeability are equal to their average within two dielectric layers, while the magnetoelectric coupling is zero. If we go further by taking the first-order approximation, we obtain As both M 1 and M 2 are off-diagonal matrices, then their commutator leads to a diagonal matrix, the components of which correspond to those of M eff in (2.15), i.e.
provides the first-order correction to the leading order approximation (classical homogenization) in (3.5). This first-order correction is encompassed in the following magnetoelectric coupling coefficients: Note that the tensor K eff is not only frequency-dependent but also exhibits spatial dispersion. The latter leads to K = K ⊥ when k = 0. Furthermore, if we consider the second-order correction, a term with 'double commutator' will appear in the asymptotic expansion with commutator of M 1 and M 2 given in (3.7), and double commutator According to the definitions of σ m in (2.8), σ K and σ K in (2.17), we have Thus, (3.11) reduces to ]. Hence, the terms arising from 'double commutator' lead to which is again an off-diagonal matrix. Substituting (3.14) into (3.10) and comparing with the form of matrix M eff in (2.15), we find and Finally, the expressions of ε eff , μ eff are as follows: All the effective parameters are frequency-dependent and with spatial dispersion. Note that, the expressions above remain well-defined in the limit ω → 0 since the ratio k/ω is fixed by the angle of incidence. These expressions turn out to be equivalent to the ones reported in [14,25,27], where the effective refractive index n eff is expanded using a power series of period-to-wavelength ratio Λ/λ. Taking eqn (5) (s-polarized incidence is considered) in the paper of Gu & Yeh [27] as an illustration, the dispersion relation for a two-component layered medium is approximated by taking the fourth order of O[(Λ/λ) 2 ] as where K and β are the z and x components of the Bloch wavevector, a and b are the thicknesses of the alternating layers, Λ = a + b is the period, n 1 and n 2 are the indices of refraction of the corresponding layers, c the velocity of light in vacuum (c −2 = ε 0 μ 0 ), and Comparing with the notation in our formulae, we have which contains the terms of order ω 2 and ω 4 . On the other hand, for the effective medium in our HOH process, the dispersion relation of k 3 versus ω is Let us substitute the expressions of effective parameters (3.16) into (3.21) and collect the terms up to order ω 4 in the calculation process, we obtain the same formula as (3.20). Similar calculation applied to a p-polarized incident wave shows that again our homogenization provides exactly the same effective index. Hence, it is stressed that the effective parameters ε eff , μ eff and K eff identified using the HOH algorithm contain more information than the single effective index parameter in [14,25,27], e.g. artificial magnetism and magnetoelectric coupling from periodic dielectrics that is not captured in the derivation of the refraction index.
In the asymptotic process, we have noticed that magnetoelectric coupling comes from the odd-order approximation, while artificial magnetism and high-order corrections to permittivity emerge from the even-order approximation in (3.3). This can be explained in the following way: the matrix M m (m = 1, 2) for the dielectric layer is off-diagonal, the terms of the odd-order correction usually contain odd commutators, hence a diagonal matrix will result, the components of which correspond to K and K ⊥ . However, the terms of the even-order correction contain even commutators, hence the resulting matrix is always off-diagonal. Physically, this introduces artificial magnetism and high-order corrections to permittivity.
Moreover, these results are fully consistent with descriptions in terms of spatial dispersion [28,29] where, expanding the permittivity in power series of the wavevector, first order yields optical activity and second order magnetic response. The equivalence of these two descriptions (frequency and wavevector power series) is confirmed by considering a unit cell with a centre of symmetry, for example, a stack of three homogeneous layers (permittivity ε m and thickness , it is found that K eff = 0, and thus it is retrieved that both magnetoelectric coupling and optical activity vanish in a medium with a centre of symmetry [28]. The present expansion in power series of frequency provides a new explanation for artificial magnetism and magnetoelectric coupling. The analytical expressions (3.16) of effective parameters can be used to analyse artificial properties. In particular, we show from (3.16) that: artificial magnetism, previously proposed with high contrast [2,3,12], can be obtained with arbitrary low contrast; and magnetoelectric coupling, previously achieved in Ω-composites [30], can be present in simple one-dimensional multi-layers. Note that, one can obtain more accurate asymptotic expressions for the effective parameters with more terms in (3.9) and (3.16), by taking higher order approximations in (3.3).

Extension of Baker-Campbell-Hausdorff formula for m layers
In §3, we have introduced our HOH algorithm for a periodic multi-layered stack consisting of an alternation of two dielectric layers, wherein the BCH formula is implemented. In this section, we investigate the extension of HOH to a stack with m layers in a unit cell, therefore a new form of BCH formula should be explored. We start with m = 3, which means a multi-layered stack with an alternation of three dielectric layers is considered and the thickness of each layer in a unit cell is h m with m = 1, 2, 3. Then, the transfer matrix of one unit cell can be expressed as or in a more compact form which defines a product of matrix exponentials. Obviously, it can be developed using an iteration of the BCH formula. First, we suppose and A can be derived through (3.3) The term A (0) is the zeroth-order approximation, while A (i) represents the ith (i ≥ 1) order correction and more precisely Here, we assume Z = Z (0) + Z (1) + Z (2) + · · · , where Z (m) is the mth (m ≥ 1) order correction for Z. The zeroth-order approximation Z (0) is simply the sum of A 1 , A 2 and A 3 , The first-order correction including a single commutator of these three matrices A 1 , A 2 and A 3 is . (4.9) and the second order including a double commutator is A similar algorithm holds for the third and higher orders, but it will not be further explored here.
As the BCH formula for (4.2) has been derived, one can easily perform the homogenization of a multi-layered stack with an alternation of three dielectric layers. Here, we assume that the third layer of the unit cell is identical to the first layer, i.e. ε 3 = ε 1 and h 3 = h 1 , as well as M 3 = M 1 ; taking (4.9) with A 1 = A 3 , we have Z (1) = 0. It should be noted that all the odd-order corrections in the HOH asymptotics vanish, which can be attributed to the centre-symmetric property of the structure [22]. By contrast, even orders rule the approximation process in that case. Applying the formulae (4.8)-(4.10) to (4.1), one deduces the expressions for the effective parameters at the second-order approximation as follows: and The effective magnetoelectric coupling tensor K eff is equal to zero because Z (2p+1) = 0, and only artificial magnetism and high order corrections to the permittivity persist. A similar calculation can be applied to higher order approximation, e.g. the expressions of these parameters at the fourth-order approximation under a normal incidence are discussed in [33]. So far, we have discussed the HOH asymptotic for a multi-layered stack consisting of an alternation of two layers as well as three layers; and the BCH formula has also been amended correspondingly. If we extend this asymptotic procedure to a more general case, i.e. we consider a stack with an alternation of m (≥ 3) layers, then the transfer matrix of a unit cell becomes Once again, tedious iteration of BCH in (4.12) can produce all the formulae for different orders of approximation. Here, we just list the formulae from the zeroth-order approximation up to second-order correction: These formulae can be checked by taking m = 3 and then compared with equations (4.8)-(4.10).
Apart from an iteration of the BCH formula, another method to obtain the approximation for Z would be to expand each exponential function in the form of a Taylor series and collect the terms with the same order.

Corrector for high-order homogenization asymptotics
In this section, we introduce the corrector for the asymptotic error in the HOH algorithm, where a structure with thickness constant with respect to the frequency or the wavelength should be considered. Considering the multi-layered stack and its effective medium shown in figure 1 Straightforward upper bounds for S n and T n are: Then developing the exponential functions in (5.5) as a series, In periodic classical homogenization [7,34], only the leading order term is kept in the asymptotic procedure, hence the corrector is of order 1/n. Here, we find that: Let us emphasize that our iterative procedure amounts to adding more and more terms in the asymptotic expansion of classical homogenization and thus it improves the corrector of classical homogenization. Similar ideas have been implemented in the high-frequency homogenization recently developed by Craster et al. [11], however with no correctors being derived therein. It would also be interesting to see how randomness would affect our correctors: at order zero, the corrector is known to vary between n −1/2 and n −1 [35]. Furthermore, using the same proof process, we can also obtain the second-order estimate for the ansatz (5.2). Upper bounds for S n and T n are once again straightforward:

Proposition 5.2. If A 1 and A 2 in (5.2) are bounded by the number A /2, then
Moreover, developing the exponential function in (5.15) as a series, It is noted that as the approximation order increases, the speed of the convergence defined by the difference between the transfer matrices of the multi-layers and the effective medium increases by a factor 1/n, hence it seems natural to conjecture that for the higher order approximation, the estimate between the transfer matrices of the multi-layers and the effective medium will be much more accurate with an error of 1/n p , with p the order taken in HOH approximation process.
Similarly, the asymptotic corrector can be derived for the stack with m layers. Here, we explore the corrector for HOH approximation in a multi-layered stack with three layers, where the permittivities are ε 1 , ε 2 , ε 3 and thicknesses areh 1 /n,h 2 /n,h 3 /n in a unit cell. Transfer matrices of the multi-layered stack approach the effective medium as follows: Then developing the exponential function as a series, The proof indicates that the corrector is of order n −1 when taking the classical homogenization (zeroth-order approximation) for a multi-layered stack with an alternation of three layers, this can be adopted for the m layers case. Correctors for higher order approximations as well as for a multi-layered stack consisting of an alternation of m layers can be obtained by the same algorithm.

Numerical calculations: dispersion law and transmission property
In this section, we numerically investigate the asymptotic estimate between the multi-layered stack and its effective medium obtained from HOH algorithm. We consider both their spectral (dispersion law) and scattering (transmission curve) properties. According to the previous analysis, the transfer matrix defined by the exponential function of matrix M is analytical, and it can be expanded as a Taylor series, taking the transfer matrix of the effective medium as an example, Considering an s-polarized incident wave, the column vector in (2.7) is defined bŷ The matrix M eff in (2.15) is a 2 × 2 matrix, and we obtain where I is the 2 × 2 identity matrix. Similarly, the transfer matrix of the dielectric layer L m is The transfer matrix T of the unit cell consisting of two dielectric layers is derived from the above expression as follows:

(a) Dispersion law
A general expression for the dispersion law in a periodic structure is defined by the trace of the transfer matrix T of a single period [36][37][38]. As the eigenvalues and eigenvectors of T (and thus any power of T) are the Bloch phase factors and Bloch states of the periodic structure (in the limit of infinite n), respectively, further physical insight can be achieved in the single period matrix T. For a multi-layered stack with two layers, we have [39] tr(T) 2 = cos(β 1 h 1 ) cos(β 2 h 2 ) − 1 2 where β m is defined in (6.7), while for the effective medium, tr(T eff ) 2 = cos(k eff h). (6.10) Here, k eff defined by (6.4) can be obtained by substituting the expressions of the effective permittivity, permeability and magnetoelectric coupling, which are derived from the HOH algorithm.
Considering a normal incident plane wave (k = 0), the s-polarization and p-polarization coincide, as ε = ε ⊥ , μ = μ ⊥ , K = K ⊥ . Hence, we take an s-polarized incident wave as an example, and assume the two dielectric layers of the stack are Glass and Silicon, respectively; the relative permittivities are ε 1 /ε 0 = 2 and ε 2 /ε 0 = 12 and the filling fractions are f 1 = 0.8 and f 2 = 0.2. For the sake of illustration, in the HOH algorithm, we take the third-, seventh-and 19th-order approximations for the effective medium. The curves of the effective permittivity, permeability and magnetoelectric coupling at the 19th-order approximation versus normalized frequency are depicted in figure 2a, whereω = ωh/c 0 , the expressions of these effective parameters are omitted here to save space. It is observed that all three curves increase with the frequency, wherein the effective permeability (dashed line) has values greater than 1, and effective magnetoelectric coupling (dotted-dashed line) is non-vanishing. In other words, artificial magnetism and magnetoelectric coupling can indeed be achieved from dielectrics through HOH asymptotics, as has been predicted in the theoretical analysis of §3. Figure 2b shows the dispersion law of the stack as well as that of the effective medium at different orders approximation, which is obtained by substituting the expressions of the effective parameters into k eff of (6.4) and then (6.10). The lower edge of the first stop band is denoted byω g = 0.202, where tr(T)/2 = −1 [40]. It should be noted that a good agreement between the dispersion laws of the stack and the effective medium can be observed in the lower frequency band, and the asymptotics of these two curves can be improved by taking higher order approximation, e.g. the 19th-order approximation generates a curve (dotted-dashed line), which nearly coincides with that of the stack (solid line) from zero frequency (quasi-static limit) up to a normalized frequency around 0.18.
Moreover, obliquely incident plane wave in s-polarization as well as in p-polarization are also analysed, wherein the same parameters of the stack are taken and the incident angle is θ i = 30 • , the dispersion laws of the multi-layered stack and the effective medium are shown in figure 3. The third-and seventh-order approximations are considered in the HOH algorithm, x 3 Figure 4. Schematic of the reflection and transmission for an incident wave on a stack of n identical unit cells, the thickness is nh, the upper and lower regions being vacuum.
similarly, asymptotics improve in conjunction with higher order approximation for both s-and p-polarized waves.

(b) Transmission property
Apart from the dispersion law, asymptotics between the transmission curves of the multi-layered stack and the effective medium is another important feature to be checked. Figure 4 shows a schematic of the reflection and transmission for an incident wave U i on a stack of n identical unit cells, the thickness of which is denoted by nh; the upper and lower spaces of the stack are supposed to be vacuum with permittivity ε 0 and permeability μ 0 . We assume the incident wave in form of Fourier decomposition of (2.3) iŝ where k 2 3 = ω 2 ε 0 μ 0 − k 2 1 − k 2 2 , and U 0 is the amplitude of the incident electromagnetic field. The reflected and transmitted waves arê Assume that T nh is the transfer matrix of the stack with thickness nh, we havê  Similarly, the transmission curves of the stack and the effective medium under an oblique incidence with θ i = 30 • in s-polarization as well as p-polarization are shown in figure 6. Once more, the asymptotic approximation between the two transmission curves of the stack and the effective medium is quite good throughout the low frequency band, and improves with higher order approximation as it transpires from the dispersion laws of figure 3.

(c) On logarithm of transfer matrix and analyticity
From figures 3 to 6, it is noted that the asymptotics break down around the lower edge of the first stop bandω g between the dispersion laws and the transmission curves of the multilayered stack and its effective medium, no matter how high the order of approximation is. This invalidity is owing to the power-series expansion of X = iM eff h in (3.3) which diverges atω g . Indeed, we choose BCH formula to obtain the approximation for matrix M eff and extract all the information required for effective permittivity, permeability and magnetoelectric coupling, by taking X = log{exp[A] exp[B]} in (3.3). In complex analysis, a branch of log(z) is a continuous function L(z) defined on a connected open subset G in the complex plane, such that L(z) is a logarithm of z for each z in G [42]. An open subset G is chosen as the set C \ R − obtained by removing the branch cut (thick solid line) along the negative real axis and the branch point where Q is a square matrix consisting of the eigenvectors of T, and Λ is a diagonal matrix whose components are the eigenvalues (denoted by λ ± ) and λ ± = a ± i 1 − a 2 (6.23) with a = tr(T)/2 the half trace of transfer matrix T. Furthermore, in order to derive the effective matrix M eff for those effective parameters, we take the logarithm of T  g (denoted by filled point). In the lower frequency band, we choose λ + and λ − (starting from +1) that approach −1 via the upper half-circle path and the lower half-circle path, respectively, where the paths lie in the open subset G: log(λ ± ) is always analytical and unique. However, at the edge of the stop band where λ + = λ − = a = −1 on the negative real axis, the logarithm is no longer analytical when its arguments meet at the branch cut of the logarithm: this implies that an expression of M eff as a power series of the frequency ω has its radius of convergence bounded byω g . In other words, the effective parameters lose their efficiency for the approximation at frequencies higher than that of the first stop band, but they work just fine throughout the lower pass band. In order to achieve all frequency homogenization for a periodic structure, a new set of effective parameters (e.g. effective refractive index and surface impedance) should be introduced [33], where the analytical property of the transfer matrix in the complex plane is ensured.

Frequency power expansion of the transfer matrix
Although the function log{exp[A] exp[B]} is no longer analytical at the lower edge of the first stop band, the transfer matrix T = exp[A] exp[B] is analytical in the whole complex plane and can be approached by a power series at any frequency, a fact that will be numerically checked in this section. In mathematics, an exponential function can be approximated by a Taylor 3 3! + · · · 1 + ωN 1 + (ωN 1 ) 2 2! + (ωN 1 ) 3 3! + · · · (7.2) with notation ωN i = iM i h i . Let us expand and organize (7.2) in powers of ω, better approximate T eff , one needs to push the approximation to the next even power p (i.e. p = 8, the thin solid line), which changes the curvature and gives a sharper estimate in the stop band region, although its intersection with the horizontal axis defines an approximate position for the upper edge of the stop band. This can be improved by taking a larger p in (7.4). Altogether, the larger the p in (7.4), the more accurate the approximation between the dispersion curves of the effective medium and that of the multi-layered stack. Moreover, according to the expression of the transmission coefficient in (6.20), we take p = 8 in (7.4) for T eff , the transmission curves of the stack (solid line) and the effective medium (dashed line) are compared in figure 8b. A good agreement between these two curves can be observed up to the first stop band, as predicted in figure 8a. Similar calculation can be applied to an oblique incidence. This demonstrates that the transfer matrix of the effective medium can be approached as a frequency power series at any frequency.

Concluding remarks
We provide a rigorous HOH algorithm for one-dimensional moderate contrast photonic crystals, where the period of the structure approaches the wavelength of the applied field. From an expression of transfer matrices in terms of exponential functions, Sophus Lie and BCH formulae are applied to produce HOH asymptotics. The analytical expressions of the effective parameters are derived for a stack with two layers in §3, where the artificial magnetism and magnetoelectric coupling effect are achieved in a moderate contrast periodic structure. In §4, we explore the extension of HOH algorithm to a stack with an alternation of m dielectric layers, and derive the expressions of the effective parameters for a centre symmetric stack: the magnetoelectric coupling vanishes while the artificial magnetism can be achieved with non-resonant periodic structures. Furthermore, the corrector for the approximation of a finite stack by its effective medium has been discussed in §5: the error estimate is of order 1/n p with p the order of the approximation. Finally, based on the expressions of the effective parameters, we numerically validate our approximation method by comparing both the dispersion law and the transmission property of the stack and its effective medium in §6. The good agreement between these curves demonstrates that the HOH approximation is efficient throughout the lower pass band, while at the edge of the first stop band the logarithm function is no longer analytical. Finally, we investigate the approximation for the transfer matrix instead of the matrix M eff of the effective medium by frequency power expansions. The dispersion laws as well as the transmission curves of both the multi-layered stack and the effective medium are explored in §7, and the excellent numerical agreement confirms that the transfer matrix of the effective medium can be approached by a power series at any frequency.