Spatially Dispersive Inhomogeneous Electromagnetic Media with Periodic Structure

Spatially dispersive (also known as non-local) electromagnetic media are considered where the parameters defining the permittivity relation vary periodically. Maxwell's equations give rise to a difference equation corresponding to the Floquet modes. A complete set of approximate solutions is calculated which are valid when the inhomogeneity is small. This is applied to inhomogeneous wire media. A new feature arises when considering spatially dispersive media, that is the existence of coupled modes.


Introduction
The standard method of calculating the Bloch modes for a given lattice is to prescribe the constitutive relations for each point within the unit cell. The corresponding dispersion relations for the bulk material are then calculated either theoretically or numerically. One may interpret this dispersion relation in terms of an anisotropic permittivity and permeability Bulk and µ Bulk for a homogeneous medium. These, in general, are temporally and spatially dispersive. However, usually the permittivity and permeability and µ of the materials used to construct the unit cell are considered to be only temporally dispersive, that is, they depend on the frequency ω §. Up to now very little theoretical research has considered periodic media where the materials inside the unit cell itself are spatially dispersive, i.e. (ω, k) and µ(ω, k), where k is the wave number. A spatially dispersive medium responds not only to a signal at a particular point but to signals in the neighbourhood of that point. For this reason it is also known as a non-local medium.
The use of a spatially dispersive inhomogeneous medium enables one to shape the propagation modes, extending the possibilities from the usual sinusoidal shape. This may have a number of advantages. For example in an accelerating drift tube, one could flatten the mode shape to enable higher acceleration for a given peak field. By contrast in signal transmission one may desire a higher peak for a given total energy.
In order to construct an inhomogeneous unit cell it is necessary to consider a medium where the permittivity and permeability depend on both the wave vector k and the position x within the unit cell, i.e. (ω, k, x) and µ(ω, k, x). In order to make sense of the arguments of (ω, k, x) we no longer interpret (ω, k, x) as the Fourier transform of the permittivity response function (t, x) but in terms of a differential equation [1]. This is given in equation (3) below. An alternative interpretation of (ω, k, x) is in terms of a susceptibility kernel [2,3].
This approach is in contrast with most theoretical articles which consider inhomogeneous spatially dispersive media. In [4][5][6][7][8][9][10][11][12] the inhomogeneity is restricted to looking at two homogeneous media which are connected and investigating the corresponding additional boundary conditions. For this article we will assume that there is a linear constitutive relationship between the polarization field P (t, x) = D(t, x)− 0 E(t, x) and electric field E(t, x). To simplify the analysis we assume: there is no magnetization so that H = µ −1 0 B and that all fields are functions of time t and one spatial coordinate x = x 1 , thus independent of x 2 , x 3 . In frequency domain, k 2 = k 3 = 0 and we set k = k 1 . The electric and polarization fields are either linearly polarized transverse modes or longitudinal modes. For transverse modes E(t, x) = E(t, x)e 2 , P (t, x) = P (t, x)e 2 and B(t, x) = B(t, x)e 3 in the (e 1 , e 2 , e 3 ) frame. These assumptions automatically satisfies the two non-dynamic source free Maxwell's equations ∇ · B = 0 and ∇ · D = 0 and the remaining Maxwell's § Note that in this paper we use ω to denote the temporal frequency rather than the angular frequency for which it is commonly used. equations can be combined to give Alternatively we assume that the electric and polarization fields longitudinal, and the magnetic field vanishes. I.e. E = E(t, x)e 1 , P = P (t, x)e 1 and B = 0. Maxwell's equations are automatically satisfied if The permittivity relationship between P (t, x) and E(t, x) analysed here is where β(ω, x) is the finite propagation speed of the polarisation wave multiplied by a constant, giving it dimensions of length, and L(ω, x) is a dimensionless constitutive quantity. The Fourier transform of P (t, x) with respect to t is given bŷ If L(ω, x) = L(ω) and β(ω, x) = β(ω) are independent of x then we can Fourier transform (3) to give is the Fourier transform of P (t, x). Such a constitutive relation can be made using diagonally crossed wires [1]. If L(ω) = ((ω − iλ 0 ) 2 + ω 2 P )/ω 2 P then (5) is the generalisation of a single-resonance Lorentz model of permittivity [13] to include spatial dispersion. For longitudinal modes, one can also use a wire medium [14,15]. Using the longitudinal component of given by [15, eqn. (4b)] where, the grid spacing b in the same in both y and z and For transverse modes, we consider a periodic structure where L satisfies L(ω, x + a) = L(ω, x) (9) Comparing notation: 2πω (here) = k (Song) and 2πk (here) = q x (Song) . We set a (Song) = b (Song) = b (here) . and β(ω, x) = β(ω) is independent of x. We further assume that the amplitude of the inhomogeneity is dominated by the first mode, that is In the case of longitudinal modes, it is the plasma frequency k p given in (8) which is periodic in x, caused by a periodic variation in r.
In both the transverse and longitudinal cases, the goal in this article is to find solutions to the source free Maxwell equations, i.e. the dispersion relations. However since the medium is inhomogeneous it is not possible to find single mode solutions of the form e 2πi(ωt+kx) . In section 2 we combine the polarisation differential equation (3) with Maxwell's equation to get a single 4th order ODE. Since this has periodic coefficients, this is an example of Floquet's equation. Hence we look for Floquet's modes of the form where κ, 0 ≤ κ < 1 is the phase. This gives rise to a difference equation for P q Λ(ω)P q−1 + f q (ω)P q + Λ(ω)P q+1 = 0 (12) where in the transverse case This is shown in the appendix. In the Longitudinal case, for a wire medium, where k 2 0 = k 2 p | r=r 0 is a constant. We can trivially solve (12) by arbitrarily fixing P 0 and P 1 and then applying (12) iteratively to find all subsequent P q . However, in general, this will lead to the coefficients P q diverging, i.e. |P q | → ∞ as q → ±∞. The Fourier transform of this would therefore be a non-smooth wave which may not even be continuous. Therefore, for physical solutions, we demand that In section 3 we make the assumptions that the inhomogeneity Λ(ω) is a small constant. This enables us to give a dispersion relation for ω and κ, as well as the harmonics P q . In the unperturbed case Λ = 0, then (12) becomes f q (ω n )P n q = 0, which we solve by setting ω n = Ω n and P n q = δ n q where f n (Ω n ) = 0. We find that, analogous to non spatially dispersive media, there exist an infinite set of modes, (ω n ,P n (x)), where n ∈ Z. Thus we write (11) as where We "normalise" (16) by setting The relationship between ω n and κ which satisfy (18) corrected by the perturbation (17) is a dispersion relation. An example of the dispersion relations, both of the unperturbed system and the perturbed system, is given in figure 1 for the transverse case and figure  4 for the longitudinal case. We do not solve (12) exactly. Instead we find approximate solution to some order, that is we solve for some order O(Λ r ). In this article we solve (20) at least to order O(Λ 3 ), which is where the new features appear. However in many cases we can do much better. For example, with the uncoupled modes given by (33) below, P n q = O(Λ |n−q| ) so that we can trivially solve (20) to order O(Λ |n−q|−1 ), although in fact we solve (20) to order O(Λ |n−q|+2 ).
There are two classes of solutions which we call uncoupled modes and coupled modes. The uncoupled modes have the property that and are given by Solving (12), (13), (17) and (18), gives rise to a dispersion relation between κ and ω n . An example is given in figure 1. The existence of coupled modes is a new feature that arises with spatially dispersive media. These occur when (21) is not satisfied. That is there exist integers n and m and a frequency Ω n ∈ C such that n = m and In this case the modes corresponding to P n and P m , (16), are coupled by the perturbation. These are explored in section 4. As stated we wish to satisfy (20) to order O(Λ 3 ). This can easily be achieved if n − m ≥ 4, where we say that coupled modes have decoupled to order O(Λ 3 ). However if n − m ≤ 3 then new solutions appear. These require taking higher order expansions for ω n and P n q . Furthermore if n − m ≤ 2 then the frequency is shifted. In the case when κ = 0 one can use the symmetry in order to construct the odd and even modes. An example of the lowest two modes is given in figures 2 and 3.
We see in figure 3 that the spatial profile of the lowest modes are significantly altered by the presence of the spatial dispersion. This suggests that with higher order harmonics cos (4πx/a) etc. in (10), one will be able to choose the electric field profile.  Figure 2. The odd P o 1 (red) and even P e 1 (black) transverse modes for L 0 ≡ 1, Λ ≡ 0.75, a ≡ 1, β ≡ 1. In this case κ = 0, ω o 1 = 0.399 and ω e 1 = 0.753i.
x P (x)  P (x) x In section 5 we present an alternative numerical method for finding the permitted frequencies andP (x). Finally in the concluding section 6 we discuss the application of inhomogeneous periodic structures.

Floquet's Modes
For transverse modes, taking the Fourier transform of (1) with respect to t gives Here prime is the partial differentiation with respect to x. In most cases, in the following, we will not explicitly write the ω argument inÊ,P , L and β.
Taking the Fourier transforms of (3) and (24) with respect to x we get Spatially Dispersive Inhomogeneous Electromagnetic Media with Periodic Structure 8 Combining these into a single convolution equation gives or equivalently hence (11). In the appendix we show that from (10), (26) and (28) the coefficients P q satisfy the difference equation (12) with f q given by (13).
Observe that having higher harmonics, such as cos (4πx/a) in (10), will result in (12) being replaced by a higher order difference equation. Likewise letting β depend periodically on x will also increase the order of the difference equation.

Longitudinal modes in wire media with periodic variation
The case for longitudinal modes is very similar. Using (7), (8) with the radius of the wire varying periodically, see figure 6 r where k 2 0 = k 2 p | r=r 0 . Substituting into (8) we see We may consider simply setting r(x) = r 0 +r 1 cos(2πx/a), which will produce terms such as cos(2πx/a) 2 and will increase the degree of the difference equation (12). For small enough r 1 then we can ignore the higher terms cos(2πx/a) 2 , then Λ ≈ (r 1 k 4 p b 2 )/(2πr 0 ). Combining (25), (2) and (28) we once again get the difference equation (12) where f q (ω) is given by (14). The dispersion relation for ω and κ is plotted in figure 4. Observe the the perturbed modes are dramatically different near n = 0 and κ = 1 2 . This is near the point where the modes couple. The detailed analysis of this behaviour is future work. An example distorted mode is given in figure 5. The result is a three dimensional periodic structure, with periods (a, b, b). This could also be analysed using numerical techniques, however this would be more time consuming.

Approximate Analytic solutions
As stated, we assume here that the inhomogeneity of the medium given by Λ in (10) is small and give approximate analytic solutions to (12) which are valid to at least O(Λ 3 ).
Let Ω n be a solution to (18), f n (Ω n ) = 0. For each Ω n we can find the corresponding P n q such that (12) is satisfied to some order given below.
In the case of transverse modes, from (13) it is clear that where Q = n + κ and hence κ = Q( mod 1) and n = Q , the largest integer n ≤ Q. Thus giving the shape of figure 1. The perturbation given in (32) below also has this property. In the longitudinal case F (Q, ω) is replaced by F (Q, ω) = ω 2 − k 2 0 − Q 2 /a 2 . In this section we consider the generic case, i.e. that there are no other integer m such that f m (Ω n ) = 0. The special cases when there exists n = m ∈ Z and Ω n such that f n (Ω n ) = f m (Ω n ) = 0 are considered below in section 4. Let and where These (ω n , P n q ) solve the difference equation (12) to an order of Λ which depends on n and q as follows and ΛP n−1 n + f n (ω n )P n n + ΛP n+1 Spatially The proof of (35) and (36) is given in the appendix. As in standard perturbation theory, improvement in the accuracy of (32) and (33) can be made by higher order corrections.

Coupled Modes
In this section we consider the case when there are two distinct integers n and m and Ω such that f n (Ω) = f m (Ω) = 0. By writing (13) as a quadratic for (q + κ) 2 we see that where One of two situation can occur. First, n and m correspond to the same root C ± (Ω), so that This implies 2κ = −n − m and so either κ = 0 or κ = 1 2 . These will be considered in sections 4.5 and 4.6 below.
In the second case set Then eliminating κ we have either Since, for most of ω, C + (ω) ± C − (ω) are continuous function these situations will also occur in general.
In the Longitudinal case, from (14) we see that only the case (39) can occur. Since two modes (ω n , P n q ) and (ω m , P m q ) are coupled, we find that, in all cases, there are two new modes.
It may be possible for specially chosen L 0 (ω) and β(ω) that there may exist three or four roots integer roots to (13). However these will not be considered. Similarly one can construct L 0 (ω) and β(ω) such that f n (Ω) = 0 or f m (Ω) = 0. In which case formulae such as (45) will not be valid. These will also not be considered.
In the special cases the dispersion curves become isolated points relating ω and κ. For example when κ = 0 and f (Ω) = 0 .
Recall our goal is to solve (20) to at least order O(Λ 3 ). We define the right hand side of (20) as Q n q , i.e. Q n q = ΛP n q−1 + f q (ω n )P n q + ΛP n q+1 (42) so that we require Q n q = O(Λ 3 ).
Spatially Dispersive Inhomogeneous Electromagnetic Media with Periodic Structure 11

Decoupled coupled modes.
In the case when n − m ≥ 4 then the coupled modes, decouple to order O(Λ 3 ), then the two solutions have the same ω n given by (32) and (18). The corresponding modes are given by The second solution, P m q , is given by In this case n − m = 1 we have two solutions given by where ω n = ±(F n F n−1 ) −1/2 The corresponding modes are given by Equations (45)  Of all the cases the case when n − m = 2 is the most complicated. The two modes are coupled at order O(Λ 2 ). Therefore it is necessary to take the second order perturbation to see the coupling. The solution for ω n is given by where α is the solution to the quadratic The mode is described by We see that the Λ coefficients of P n−2 and the Λ 2 coefficients of P n−3 and P n−1 are related by the quantity α 2 . However it would require perturbing to order O(Λ 4 ) to determine this quantity.

Coupled modes, n − m = 3
In the case n − m = 3 one solution is given by and the mode given by where α 3 cannot be determined order at order O(Λ 3 ). The other solution is given by exchanging n and m, i.e. by (50) and (51) where m − n = 3.
The remaining case is when n = 1. In this case some simplification arises and the solutions may be interpreted as odd and even modes. The frequencies are and the corresponding P n q P 1,e and where R 1 1 = 1 and When κ = 1 2 then from (13) we see that f q (ω) = f −1−q (ω) for all ω. Hence F q = F −1−q and F q = F −1−q . Thus since f n (Ω n ) = 0 then f −1−n (Ω n ) = 0.
Again if n ≥ 3 then the modes decouple. However if n = 0 or n = 1 then the modes remain coupled and we can use the results of n − m = 1 or n − m = 3 respectively. Unlike the case when κ = 0 no significant simplification takes place. However we note that (45) simplifies to

Numerical Approaches
A numerical approximation scheme, which is valid if Λ L 0 , is as follows: Choose an integer N ≥ 2. Then assume that P m ≈ 0 for |m| > N thus truncating the infinite set of equation given by (12) to a set of 2N + 1 linear equations for {P −N , . . . , P N }. Write this in matrix language M b = 0 where M is a (2N + 1) × (2N + 1) matrix with X X X X X y 9 Inhomogeneous periodic cross wire structure to create inhomogeneous periodic spatially dispersive medium in cavity.
Protons Drift tubes Optimised electric field profile

Conclusion
In this article we describe a model for spatially dispersive media with a periodic structure. Two types of media are considered, one based on a single-resonance Lorentz model of permittivity, the other on a wire medium. Approximate analytic solutions to Maxwell's equations are found for such media for small magnitudes in the inhomogeneities. We demonstrate the existence of coupled modes. Higher order correction should be possible by repeating the procedure. If one were to solve (20) to order O(Λ 4 ) then one would expect coupled modes where n − m = 4 would remain coupled. By combining the solutions above with additional boundary conditions [4-9,16] one will be able to predict the behaviour or electromagnetic radiation as it passes though a slab of periodic spatially dispersive medium.
It should be possible to construct a meta-material with such periodic spatially dispersive properties, especially the wire medium. These can be used to test that the packets given above are indeed supported. In addition, numerical simulations using, for example, FDTD codes in a periodic cell, can be made.
As stated in the introduction, one of the key advantages for using inhomogeneous spatially dispersive periodic structures, is that one can choose the shape to the propagating wave. For example, one may consider adding a spatially dispersive media to a drift tube in order to optimise the fields the particles experience. See figure 7. Another example, in data transmission, it to shape the wave to have a higher peak value. For example for longitudinal modes ifÊ(x) = −P (x) = 3 −1/2 ( cos(2πx/a) + cos(6πx/a) + cos(10πx/a)), then the peak field increases by √ 3 but the root mean square per length ofÊ is unchanged when compared to the sinusoidal case.