Customized longitudinal electric field profiles: using spatial dispersion in dielectric wire arrays

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.


I. 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 dielectric wire medium, a choice 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 * https://orcid.org/0000-0003-1597-6084; j.gratus@lancaster.ac.uk † https://orcid.org/0000-0001-5744-8146 ‡ https://orcid.org/0000-0002-  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 FIG. 1: A diagram showing a standard uniform-radius wire medium (left) compared with a varying-radius wire medium (right) for wave profile shaping. By understanding and characterizing the properties of a uniform medium, we can use them as a basis to predict the behaviour of varying-radius medium, and so perform sub-wavelength mode profile shaping. The important feature of such wire media is that they naturally support modes with longitudinal electric fields, i.e. fields parallel to the wires themselves. 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.
This is a well established differential equation with a set of known solutions, where the free parameters 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 IV, 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].

II. 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 where e z is the unit vector along the z axis, our chosen longitudinal orientation. Maxwell's equation is automatically satisfied and the requirement that D = 0 implies that ε = 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 uniformradius 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][23][24] for this model in the frequencywavevector representation is usually written aŝ However, in what follows we will use ordinary (not angular) frequencies 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 (4) and Fourier transform back into the frequency-physical space representation then we have where the cut-off λ has units of frequency, and is Since (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 It is not necessary to consider whether the dependence on 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 equation (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 III 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 In order to find the the radius variation R(z) we simply solve where λ 2 s (r) is given by (10) and Λ 2 (z) is given by (9) 1 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 λ 2 s (r), 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 λ 2 s (r).
Here we find longitudinal electric modes with a plasmalike dispersion relation in uniform-radius 3D simulations Aim: Validate a chosen mode for wire media and develop a method for mode profile shaping Choose a simple model. We use a 1D inhomogeneous spatially dispersive model.
Choose a desired profile which you want to replicate. We use Mathieu functions Use the model to solve for the cut-off Λ 2 (z) which satisfies your profile.
Find a combination of lattice parameters, wire radius and electrical permittivity which supports longitudinal electric modes.
Vary the radius of the system and use the dispersion relations for each radius to find the radius dependence of the cut-off λ 2 s (r) Combine and solve for R(z) λ 2 s (R(z)) = Λ 2 (z).
Implement varying wires with R(z) into a CST model. 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 de- scription, κ 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 y , but their specific values play no role in our calculations here.
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 pa- We found that that for wires of radius r = 0.20mm (black, circles) that the cut-off is λs 6.64GHz; that for r = 0.30mm (blue, squares) that the cut-off is λs 6.05GHz; and r = 0.40mm (red, diamonds) that the cut-off is λs 5.92GHz. The solid lines indicate the parabolic fits to the dispersion data; whereas the dashed lines show the predicted dispersion made on the basis of the later phenomenological fit of (12). rameters above have a high band index 2 and a cut-off frequency of 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 λ 2 s , 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 λ 2 s , with β being approximately constant for all simulations at β = 0.96.
We now use our simulations to find the dispersion re-2 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.
lations 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 I.   (12).

IV. 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  FIG. 8: The dependence of λs on the radius r of a uniform wire. Note that the vertical axis does not extend down to λs = 0. Data is shown for four cases, with different permittivity values (i.e. ε ∈ {50, 100, 400, 1600}); where data points for for smaller radii than those shown significantly affected the quality of the otherwise remarkably good fit, so we excluded such points from our analysis. The dashed lines indicate the large radius (asymptotic) value for λs, which gives the value of C in (12).
profile into the wire radius variation that will support it. Inserting the definition (11) into (9) allows us to see that we simply require The factors β, a, q and the parameters which define λ 2 s (r) 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 λ 2 s (r) is bounded so we require λ 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 λ 2 + . From (14) we see there is an inverse relationship between the range of λ 2 and the length.
Thus the same overall mode shape can be designed to have different periodicities, and be made with (e.g.) either a short 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 3 , enhancing both their interaction and the utility of that interaction. 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 To generate the radius variation needed for each of our selected mode profiles, and for our standard ε = 1600 wire medium, the necessary Γ i parameters are given by table II. Implementing this R(z) function in a dielectric These were based on a cell length of L z,peak = 691mm a frequency of f peak = 7.3GHz, and a cut-off variation with λ± = 6.0, 9.9GHz. Note that this profile shape requires much longer cell lengths than the others, as a result of its a, q Mathieu parameters.   (20). For each of the three profiles, Γ1 = 0.5406647836mm for all three shapes.
wire medium should then satisfy (13), allowing for the support of longitudinal modes with the desired mode profiles, based on the chosen Mathieu function.

V. 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 FIG. 13: The flat-top field profile obtained using a simulation with a coated dielectric wire with L = 317.4mm and permittivity ε = 50, whose profile was based on wire data for an unclad wire, and varied between 0.69mm and 0.71mm. The cladding covers over the dielectric wire out to an outer radius of 3mm, and has a relative permittivity of ε = 2. 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 BaTiO3based 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.

VI. 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 dis-cussed 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.