## Abstract

We show how spatial dispersion can be used as a mechanism to customize the longitudinal profiles of electric fields inside modulated wire media, using a fast and remarkably accurate 1D inhomogeneous model. This customization gives fine control of the sub-wavelength behaviour of the field, as has been achieved recently for transverse fields in simpler slotted-slab media. Our scheme avoids any necessity to run a long series of computationally intensive 3D simulations of specific structures, in order to iteratively converge (or brute-force search) to an empirical ‘best-performance’ design according to an abstract figure-of-merit. Instead, if supplied with an ‘ideal waveform’ profile, we could now calculate how to construct it directly. Notably, and unlike most work on photonic crystal structures, our focus is specifically on the field profiles because of their potential utility, rather than other issues such as band-gap control, and/or transmission and reflection coefficients.

Published by The Optical Society under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

## 1. Introduction

Detailed control over electromagnetic field profiles (the field strength variation through space) and other electromagnetic properties is an invaluable ability. In existing accelerator technology, crab cavities [1] are used to provide a tailored field profile for control of the particle bunch; and photonic crystals are used for band gap, transmission, and dispersion control [2,3]. More recently, the emerging multidisciplinary field of metaphotonics [4] attempts to further such aims with subwavelength (meta)material design, using both electric and magnetic interactions and their cross-coupling. Such technology promises a remarkable degree of control, enabling zero (or negative) electromagnetic response, chirality control, and cloaking.

Our interest here is to be able to design any electric field profile on a wavelength scale; that is, to sculpt the carrier wave of a propagating field itself. This idea has previously been investigated in the context of nonlinearity-induced carrier shocking [5] and by harmonic synthesis [6,7]. More recent work has outlined the potential for mode profiling using modulated dielectric slabs [8,9], but only for transverse field modes. However, here we aim to control the *longitudinal* field component, not the transverse one. For this we use an all-dielectric wire medium, made of high permittivity rods surrounded by air in a rectangular lattice as in Fig. 1. This choice is crucial to our success, because it is spatially dispersive and natively supports modes in which the longitudinal electric field component can dominate the transverse ones.

The method outlined in the following sections allows for the rapid and accurate generation of any chosen field profile in the wire media structure; further, once some uniform-radius reference simulations are done in 3D, all further design need only be based on a simple 1D model. This high level of control could be used to tailor field profiles for specific applications such as enhancing ionization in high harmonic generation [5]. One could also imagine the creation of electric fields with high gradients without the usually associated high peak field values (i.e a flattened mode), which can lead to unwanted nonlinear effects. Alternatively, an electromagnetic profile with a more pronounced peak would be useful in signal processing as it would make for easier detection amongst any background noise.

However, one significant area of our interest is that of accelerator applications such as the shaping of electron bunches or uses in the plasma ionisation involved in laser wakefield accelerators [10, 11]. This is because of our new ability to gain detailed control over the field profiles and gradients, and not just aim at a general redistribution of field intensity as a means to moderate losses or heating. For example, in accelerator applications, an RF pillbox cavity is used as a way to increase the available accelerating voltage with a minimum power loss. Typical structural features involve the use of a “nose-cone” profile to concentrate the electric field near the particle beam passage; but the design involves an iterative process of numerically optimizing the cavity figures of merit based on full 3D simulations [12,13]. This is fundamentally different from our method, which requires only a 1D specification of the field profile, and then – without iteration – gives the radius variation needed.

Of course, to be truly useful in accelerator applications, we would prefer the field profiles to be co-moving with the particle bunches, or to have control over the time domain profiles. In what we present here the profiles have stationary spatial profiles and ordinary sinusoidal behavior in time; nevertheless our results show the first important step towards our eventual goal.

We will show here that field modes with a *longitudinal* electric field exist in dielectric wire media. Starting with results based on uniform-radius wires, we show that by implementing a design where it includes a variation of the radius (see Fig. 1), we can generate a modulation of the mode profile. In particular, we focus on periodic variations of the wire radius, where the electric field can then be shown to obey Mathieu’s equation [14], i.e.

*q*and

*a*control the properties of the solution; sample solutions which we will use as target mode profiles are shown on Fig. 2. In particular, for a given

*q*there always exists a discrete spectrum of

*a*values such that the solutions

*y*(

*σ*) are periodic over 2

*π*[14]. In the trivial case where

*q*= 0, we find that

*a*∈ {0, 1, 4, 9, …}. In this work the focus will be purely on periodic solutions of Mathieu’s equation, although we expect that our method can be applied to treat non-periodic solutions, and could even be generalized beyond the Mathieu equation to generate a range of non-periodic field profiles.

As can be seen in Fig. 2 some of the periodic solutions of Mathieu’s equation fit the description of the potentially desirable profiles discussed earlier. The flatter profile which could minimise the harmful impact of nonlinear effects, whereas the profile with more pronounced peaks could allow easier detection above backgound noise. Either might be useful to control the field experienced by particles located within a larger bunch in an accelerator, ensuring either a more uniform acceleration or an enhanced chirp. Our challenge here is to show how modulations in the wire radii can be mapped onto a sufficiently simple electric field equation, which then also mimics Mathieu’s equation; this is therefore similar in principle to the rather more straightforward mapping from the physics to Mathieu’s equation possible in the context of Josephson junctions [15]. This is achieved in section 4, which shows how longitudinal modes obeying Mathieu’s equation are supported, thus allowing use to generate Mathieu-like mode profiles in 3D structures using an efficient 1D design process. In addition to our theoretical modelling, we have validated our procedure for a full 3D wire medium by means of simulations using CST Microwave Studio [16].

## 2. Inhomogeneous and spatially dispersive 1D model

Wire media are a well known class of metamaterials consisting of a lattice of parallel wires, with the radius *r* of the wires being small compared to their spacing [17]. It should be noted that wire media act as metamaterials when the Bragg scattering wavelength is less than the Mie scattering wavelength (i.e when low permittivity rods are tightly spaced or rods with large spacing have high permittivity). As the ratio between permittivity and lattice spacing decreases wire media start to enter the regime of photonic crystals [18]. One of the unique features of these structures is that they support purely electric longitudinal modes which have a plasma-like dispersion relation.

Wire media are complex three-dimensional inhomogeneous media, which makes them difficult to study without introducing a simpler model for the system. The model we will use approximates the lattice of wires as a one-dimensional inhomogeneous media which exhibits spatial dispersion. Using this model will greatly simplify our efforts to study the behaviour and manipulation of the profile of longitudinal electric modes in wire media. If we can then make predictions about these modes which are validated in full 3D simulations, it will also validate the model upon which those predictions were based. We will see later that this 1D model works remarkably well.

In our model for longitudinal modes in a 1D medium the desired fields are given by

*z*axis, our chosen longitudinal orientation. Maxwell’s equation is automatically satisfied and the requirement that

**=**

*D***implies that**

*0**ε*= 0, an ENZ material [19]. Our choice of wire media for our studies was motivated by the existence in the literature of a homogenized 1D model for the uniform-radius case, which exhibits spatial dispersion. Such a medium would allow us to make the simplifications detailed above, as well as possess the required ENZ nature. However, in order to find a way of controlling the mode profiles we need an

*inhomogeneous*1D model: it needs to be both inhomogeneous (i.e. have a position

*z*dependence), as well as spatially dispersive (i.e. have a wavevector

*k*dependence), which is challenging as

*k*and

*z*are Fourier conjugate variables [17, 20, 21]. We start with the undamped hydrodynamic Lorentz model which is spatially dispersive and homogenous. The constitutive relation [22–24] for this model in the frequency-wavevector representation is usually written as

*f*=

*ω*/2

*π*, and likewise for wavevector k we will use the alternative (inverse wavelength)

*κ*=

*k*/2

*π*. It follows that in order for

**P**= −

*ε*

_{0}

**E**(i.e.

**D**= 0) then: This is a plasma like dispersion relation where

*cβ*is the polariton velocity,

*f*

_{0}is the polariton resonance frequency, and

*f*

_{p}is the plasma frequency parameter. If we rearrange Eq. (4) and Fourier transform back into the frequency-physical space representation then we have

*λ*has units of frequency, and is

Since Eq. (7) is written in physical space representation, we can now allow for an inhomogeneous permittivity, following [17]. This is achieved by allowing for a position dependent modulation of the the cut-off, which we denote by replacing *λ*^{2} by Λ^{2}(*z*). Although we might also let the relative polariton velocity *β* vary, we do not, since as we will show later, for our chosen parameters it stays nearly constant. This means that (7) becomes

*z*in Λ

^{2}(

*z*) is a result of

*f*

_{0}or

*f*

_{p}having a

*z*dependence, or of both. We can now see that (8) can be made equivalent to Mathieu’s euqation (1), simply by rescaling the spatial coordinate with

*σ*= 2

*πz*/

*L*, and setting the cut-off modulation Λ

^{2}(

*z*) to be

By doing this we have reconciled the equations of our wire media with Mathieu’s equation, a process which required making either the plasma frequency *f*_{p} or the resonance frequency *f*_{0} dependent on *z*. This can be achieved by realizing that these frequencies depend on the radius *r* of the wires in the wire media. Our goal therefore is to find the appropriate variation in the wire radius *r*, which we will denote *R*(*z*), to implement in our 3D simulations. In section 3 below we use set of full 3D simulations with uniform wires to calculate the dependence of dispersion relation (5) on the radius. Thus for some specified wire radius *r*, the cut-off value is obtained by numerical simulation; where of course this simulated cut-off, denoted *λ*_{s}(*r*), is in turn dependent on how the frequencies *f*_{p} and *f*_{0} vary with the wire radius *r*. Hence (5) can be rewritten as

*R*(

*z*) we simply solve where ${\lambda}_{\mathrm{s}}^{2}\left(r\right)$ is given by (10) Λ

^{2}(

*z*) is given by (9). We will see in the next section that

*λ*

_{s}(

*r*) follows a simple and invertible fitting formula (12), meaning that we can solve for (11) analytically.

We now need to combine these results with a 3D simulation of our structure. The flowchart shown in Fig. 3 shows our method for designing 3D structures on the basis of our 1D predictions. It requires that we find longitudinal plasma-like modes in uniform-radius 3D simulations which persist for a range of wire radii. This gives us the cut-off values as a function of the wire radius ${\lambda}_{\mathrm{s}}^{2}\left(r\right)$, and we can match these up with the 1D model’s required modulation in the cut-off. The result will be that we can find the necessary wire radius as a function of position for a full 3D mode-profiling simulation. The next section will describe the uniform-radius 3D simulations that give us ${\lambda}_{\mathrm{s}}^{2}\left(r\right)$.

## 3. 3D numerical calculations of mode dispersion: ${\lambda}_{\mathrm{s}}^{2}\left(r\right)$

Here we find longitudinal electric modes with a plasma-like dispersion relation in uniform-radius 3D simulations which persist for a range of wire radii. These CST simulations were based on a single unit cell as shown on Fig. 4, containing a dielectric rod of radius *r*, and with periodic boundary conditions in all directions, thus representing an infinite wire media. The calculations in this section all incorporate cells of this type, but with differing wire radii *r*, so that we can find *λ*^{2}(*r*) for the range of 3D structures we are interested in.

Our standard (reference) simulation was for an infinite dielectric wire medium, where the wires had a permittivity of *ε* = 1600*ε*_{0} (e.g. ceramics [25]) and a radius of 0.2mm, with lattice constants of 15mm and 13.06mm in the *x* and y (transverse) directions respectively. This use of a non-square unit cell cross-section ensures there is no ambiguity due to degeneracy between *x* and *y* orientations in the simulation output. Note that in our 3D description, *κ* is retained (solely) as the Fourier conjugate variable to the longitudinal spatial axis *z*. Transverse wavevectors do play a role in that they help determine the mode shapes, and in particular whether or not any given mode has our desired longitudinal property. As such they are dependent on 1/*L _{x}* and 1/

*L*, but their specific values play no role in our calculations here.

_{y}These simulations, as well as demonstrating the expected predominance of transverse modes, confirmed the existence of longitudinal electric modes in wire media structures. The longitudinal modes found for the parameters above have a high band index, as determined by counting the bands found by CST from lowest frequency to highest. However, note that band index determination in numerical results is often complicated or made ambiguous due to band crossings and the choice of unit cell sizing. Their cut-off frequency is about *λ*_{s} = 6GHz, as can be seen on Fig. 5. However, note that the index of the longitudinal mode differs for different simulation parameters, since the changing geometry has an effect on the details of the bandstructure.

These longitudinal modes have the significant features that the field between the wires was both strongly longitudinal and anti-parallel to the field in the wires, as shown in Fig. 6. This is ideal for applications where the resulting field is intended to accelerate or decelerate particles along their path in a free space region.

The dispersion properties of these longitudinal modes are shown in Fig. 7, where we see that they obey the plasma-like dispersion relation predicted by theoretical analysis. This was confirmed by the plotting of *f*^{2} vs. *κ*^{2} graphs which show a straight line, from which the cut-off ${\lambda}_{\mathrm{s}}^{2}$, and polariton velocity *β*, can be straightforwardly obtained, although *f* vs. *κ* graphs are often plotted instead. In our simulations, we found that the dominant effect of radius variation is on ${\lambda}_{\mathrm{s}}^{2}$, with *β* being approximately constant for all simulations at *β* = 0.96.

We now use our simulations to find the dispersion relations for a useful range of wire radii, and thus find (by fitting) the corresponding cut-offs. On Fig. 8 we can see how these change with wire radius, and by fitting this data, it is the possible to find a relation for the radius dependence of the cut-off:

Here*C*,

*A*and

*ρ*are fitting constants that depend on the wire permittivity, unit cell, and mode index. Although this relation is only approximate, the agreement should be sufficient so long as we stay within the range of data taken for the radius, between 0.1mm and 0.5mm. As can be seen on Fig. 8, these parameter fits are remarkably accurate, with

*exceptionally*small deviations from the model behaviour. For example, we can see on Fig. 7 the small discrepancy between the dispersion relation predicted by this fit, and that from the numerical data on which it was originally based. This difference is negligible on the scale of the results on Fig. 8, but increases for larger radii.

We repeated this process for thicker rods with lower permittivities, an important consideration when we wish to move towards real-world fabrication possibilities. We verified that mode profiling was still viable at lower relative permittivities, i.e. at both *ε* = 50, *ε* = 100, and *ε* = 400. As can be seen in Fig. 8, for a given frequency, the radius required for longitudinal modes to appear becomes larger for a lower relative permittivity. This means that the mode profiling structures for *ε* = 50 involve larger wire radii in the order of millimetres. Such lower index materials are not only easier to find, but may be very useful in practical applications where such thin rods may be either difficult to fabricate or too fragile. For the results shown on Fig. 8, the fitting parameters for (12) are shown on table 1.

## 4. The wire shape needed to match a Mathieu equation: *R*(*z*)

We now need to describe how to match the properties of our wire media, and the fields they can support, to the mathematical form of Mathieu’s equation. This comparison will tell us how to turn a mathematical Mathieu profile into the wire radius variation that will support it. Inserting the definition (11) into (9) allows us to see that we simply require

*β*,

*a*,

*q*and the parameters which define ${\lambda}_{\mathrm{s}}^{2}\left(r\right)$ are determined either by the simulations or by the choice of profile shape. We see that to determine

*R*(

*z*) from (13) we need only specify the frequency

*f*and the repetition length

*L*; but although we have a lot of freedom in our choice of

*L*and

*f*, there are some constraints. Notably, by taking the maximum and minimum values of the cosine function we obtain the concomittant maximum and minimum values of the cut-off, being From figure 8 we see that ${\lambda}_{\mathrm{s}}^{2}\left(r\right)$ is bounded so we require ${\lambda}_{-}^{2}>C$, where

*C*is given by (12). Further, the requirements resulting from actually constructing the wire medium device may also demand a value for the minimum radius which will then in turn constrain ${\lambda}_{+}^{2}$.

From (14) we see there is an inverse relationship between the range of *λ*^{2} and the length.

*L*and a large cut-off modulation, or a long

*L*and a small modulation.

For the work presented here we choose a frequency near that of the cut-off, which means that our field profiles are in effect standing waves. However, there is the possibility of investigating whether we can construct propagating field profiles, i.e. ones which have a finite phase velocity, by working higher up the dispersion curves. Since *β* < 0, then for small *κ*, the waves are superluminal; whereas for large *κ* the waves are instead subluminal. In figure 9 we see that setting *f* = 20GHz then one obtains a phase velocity of about *c*. Although we would then expect the profiles to also become time dependent, such an achievement would bring us closer to our goal of matching the wave profile speed to the particle speed, enhancing both their interaction and the utility of that interaction. However, due to the extra complications and more demanding computations involved, this is beyond the the parameter ranges we have examined to date, and such high modes may not have the desired hyperbolic response., Note however that since our modes are not sinusoidal, and are therefore consist of several spatial harmonics, one can no longer easily assign a single phase velocity.

Now that we have a method for determining the remaining two free parameters, *f* and *L*, it is now possible to solve (13) for *R*(*z*). This has the form

*ε*= 1600 wire medium, the necessary Γ

*parameters are given by table 2. Implementing this*

_{i}*R*(

*z*) function in a dielectric wire medium should then satisfy (13), allowing for the support of longitudinal modes with the desired mode profiles, based on the chosen Mathieu function.

## 5. Results and discussion

We followed the methodology outlined above to create varying wire unit cells for our three different wave profiles. The results of these simulations are shown on figs. 10, 11, 12, and they confirm that our method accurately reproduces the desired longitudinal field profiles given the correctly specified wire media structures. As *R*(*z*) will be a periodic function it is fairly simple to implement this profile in a unit cell with periodic boundary conditions in CST, although it should be noted that a feature of Mathieu’s equation means that the period (wavelength) of the Mathieu mode will be *L*, which is double the period of the wire variation. Though an exact replication of the radius variation is not possible in CST it is possible to make an approximation of the variation with a large number of connected conical frustums.

We have seen above, notably in Fig. 8, that despite our initial preference for very thin high permittivity wires, our approach still has currency even for thicker wires with lower permittivity. We are currently checking how far this process can be pushed, i.e. to what base wire thickness and what low permittivity.

Note that even materials that would be entirely impractical for wires if used on their own – such as BaTiO3-based ceramics [25], which can have *ε* ~ 1000+ in the GHz – might still be used if encapsulated in an outer sheath of more tractable material. Simulations have shown that adding a low-index sheath over the high-index core with varying radius does not change the outcome of our (uniform wire) simulations significantly. We can see in the results on Fig. 13 that the dispersion properties of the wire media can be sufficiently unaffected by the sheath, so that as a consequence the profile shaping will also still persist.

We are currently working on extending our numerical results into the time domain, either with external driving or internal excitation. However, obtaining such results can be problematic, due to the greatly increased computational requirements. Such simulations also require a specification for how to drive the mode profiles, as indeed is also needed for eventual experiments. Preliminary indications are that driving our mode profiling waveguides externally might be difficult due to poor input coupling. One resolution of this is to place a dipole antenna *internally*, and so drive from inside the medium, as certainly worked well for our slotted medium mode profiling investigations [8, 9]. Indeed, it is very convenient that our goal of achieving free-space mode profiling makes such a step so simple.

Beyond these theoretical extensions, two more steps need to be taken to make our concept application ready. First, we need to consider the actual structure of a device, which will most likely consist of a finite wire array housed in a metallic box. Since the field patterns present here have strong components – albeit of opposite sign – both in the wire and in free space, the best positioning for the box walls with respect to the wires needs to be investigated. Second, we have currently focussed on standing wave field profiling, rather than the travelling wave profiles needed for our aim of finding accelerator applications. However, now we have established that field profiling is possible with a high degree of accuracy, we are testing configurations aimed at building the wire media using a mode basis that is not centred around near cut-off frequencies and hence low *κ* values in the dispersion curves, but at a median *κ* values which have a faster group velocity.

Our next step will be to add time domain shaping to our toolbox, either in concert with the spatial profiling shown here, or by itself. This would use multiple frequency components, and greatly enhance the possibilities for the bunch control of charged particles in accelerator applications.

## 6. Conclusion

In this paper we have extended the theoretical analysis for metal wire media to dielectric wire media; with both reducing to a 1D homogeneous model incorporating spatial dispersion. We have shown that these dielectric media support modes with longitudinal electric fields (i.e. field parallel to the wires) with a plasma-like dispersion relation; and how the plasma frequency *f*_{p} and mode cut-off *λ*^{2} depends on the chosen radius for the wires.

This then formed a basis for our method incorporating a varying wire radius, and allowed us to show how the resulting dispersion relation for the modes could be manipulated into corresponding to Mathieu’s equation. We could then create mode profiles which replicated solutions to Mathieu’s equation; thus providing an effective approach for implementing specific Mathieu functions (solutions) as physical electric field profiles. As part of the process we required a fit for the dependence of the cut-off on the wire radius, and we found that a simple exponential model gave remarkable agreement, greatly simplifying the prediction of the wire radial variations needed. The exceptional accuracy of this fitting function suggests that a theoretical investigation should yield a useful and accurate radius-dependent simplified model of spatial dispersion, making it possible to dispense with our 3D uniform-radius unit-cell CST simulations.

We applied our method to a selection of chosen mode profiles, namely a flat-top, peaked, and triangular waveforms. The results of the frequency-domain simulations for these examples showed very good agreement, not only for straightforward single unit-cell calculations, but also for wires stabilised in an external cladding. We also discussed routes for future work, including time domain simulations, investigation of options for exciting/driving the mode profiles, device-like simulations involving multiple wires in a waveguide, and the possibilities for group velocity engineering of the modes.

## Funding

Science and Technology Facilities Council (STFC) (the Cockcroft Institute ST/G008248/1 and ST/P002056/1); Engineering and Physical Sciences Research Council (EPSRC) (the Alpha-X project EP/J018171/1 and EP/N028694/1).

## Acknowledgments

Paul Kinsler would like acknowledge the hospitality of Imperial College London.

## References and links

**1. **S. Verdu-Andres, S. Belomestnykh, I. Ben-Zvi, R. Calaga, Q. Wu, and B. Xiao, “Crab cavities for colliders: past, present and future,” Nuclear and Particle Physics Proceedings **273–275**, 193–197 (2016).

**2. **P. S. Russell and T. A. Birks, “Hamiltonian optics of nonuniform photonic crystals,” J. Lightwave Techn. **17**(11), 1982–1988 (1999). [CrossRef]

**3. **P. S. J. Russell, T. A. Birks, and F. D. Lloyd-Lucas, “Confined electrons and photons: New physics and applications,” in NATO Advanced Study Institute on Confined Electrons and Photons - New Physics and Applications, 340, 585–633 (1995).

**4. **A. Baev, P. N. Prasad, H. Agren, M. Samoc, and M. Wegener, “Metaphotonics: An emerging field with opportunities and challenges,” Phys. Rep. **594**, 1–60 (2015). [CrossRef]

**5. **S. B. P. Radnor, L. E. Chipperfield, P. Kinsler, and G. H. C. New, “Carrier-wave self-steepening and applications to high-order harmonic generation,” Phys. Rev. A **77**(3), 033806 (2008). [CrossRef]

**6. **H.-S. Chan, Z.-M. Hsieh, W.-H. Liang, A. H. Kung, C.-K. Lee, C.-J. Lai, R.-P. Pan, and L.-H. Peng, “Synthesis and measurement of ultrafast waveforms from five discrete optical harmonics,” Science **331**(6021), 1165–1168 (2011). [CrossRef] [PubMed]

**7. **J. A. Cox, W. P. Putnam, A. Sell, A. Leitenstorfer, and F. X. Kärtner, “Pulse synthesis in the single-cycle regime from independent mode-locked lasers using attosecond-precision feedback,” Opt. Lett. **37**(17), 3579–3581 (2012). [CrossRef] [PubMed]

**8. **J. Gratus, P. Kinsler, R. Letizia, and T. Boyd, “Electromagnetic mode profile shaping in waveguides,” Appl. Phys. A **123**(1), 108 (2017). [CrossRef]

**9. **J. Gratus, P. Kinsler, R. Letizia, and T. Boyd, “Subwavelength mode profile customisation using functional materials,” J. Phys. Commun. **1**(2), 025003 (2017). [CrossRef]

**10. **P. Piot, Y. E. Sun, J. G. Power, and M. Rihaoui, “Generation of relativistic electron bunches with arbitrary current distribution via transverse-to-longitudinal phase space exchange,” Phys. Rev. ST Accel. Beams **14**(3), 022801 (2011). [CrossRef]

**11. **F. Albert, A. G. R. Thomas, S. P. D. Mangles, S. Banerjee, S. Corde, A. Flacco, M. Litos, D. Neely, J. Vieira, and Z. Najmudin, “Laser wakefield accelerator based light sources: potential applications and requirements,” Plasma Physics and Controlled Fusion **56**(8), 084015 (2014). [CrossRef]

**12. **M. Larrañaga, R. Enparantza, and C. Plostinar, “Design optimisation of the re-bunching cavities for the front end test stand at RAL,” in Proceedings of the Linear Accelerator Conference (LINAC), (2010), paper MOP080.

**13. **S. W. Shin and J. S. Chai, “Optimization of the RF cavity of the medical purpose electron linac by using genetic algorithm,” in Proceedings of the Linear Accelerator Conference (LINAC), (2014), paper TUPP133.

**14. ** MathWorld, “Mathieu functions,” http://mathworld.wolfram.com/MathieuFunction.html.

**15. **S. A. Wilkinson, N. Vogt, D. S. Golubev, and J. H. Cole, “Approximate solutions to Mathieu’s equation,” http://arxiv.org/abs/1710.00657

**16. ** CST AG, “CST Studio Suite,” http://www.cst.com/.

**17. **J. Gratus and M. McCormack, “Spatially dispersive inhomogeneous electromagnetic media with periodic structure,” J. Opt. **17**(2), 025105 (2015). [CrossRef]

**18. **M. V. Rybin, D. S. Filonov, K. B. Samusev, P. A. Belov, Y. S. Kivshar, and M. F. Limonov, “Phase diagram for the transition from photonic crystals to dielectric metamaterials,” Nature Commun. **6**, 10102 (2015). [CrossRef]

**19. **D. Ramaccia, F. Scattone, F. Bilotti, and A. Toscano, “Broadband compact horn antennas by using EPS-ENZ metamaterial lens,” IEEE Trans. Antennas Propag. **61**(6), 2929–2937 (2013). [CrossRef]

**20. **J. Gratus and R. Tucker, “Covariant constitutive relations and relativistic inhomogeneous plasmas,” J. Math. Phys. **54**(4), 042901 (2011). [CrossRef]

**21. **J. Gratus, P. Kinsler, M. McCall, and R. T. Thompson, “On spacetime transformation optics: temporal and spatial dispersion,” New J. Phys. **18**(12), 123010 (2016). [CrossRef]

**22. **P. A. Belov, R. Marqués, S. I. Maslovski, I. S. Nefedov, M. Silveirinha, C. R. Simovski, and S. A. Tretyakov, “Strong spatial dispersion in wire media in the very large wavelength limit,” Phys. Rev. B **67**(11), 113103 (2003). [CrossRef]

**23. **W. Song, Z. Yang, X. Q. Sheng, and Y. Hao, “Accurate modeling of high order spatial dispersion of wire medium,” Opt. Express **21**(24), 29836–29846 (2013). [CrossRef]

**24. **C. Ciraci, J. B. Pendry, and D. R. Smith, “Hydrodynamic model for plasmonics: A macroscopic approach to a microscopic problem,” ChemPhysChem **14**(6), 1109–1116 (2013). [CrossRef] [PubMed]

**25. **J. Li, H. Kakemoto, S. Wada, and T. Tsurumi, “Dielectric properties of BaTiO_{3}-based ceramics measured up to GHZ region,” J. Electroceram. **21**(1), 427–430 (2008). [CrossRef]