# Wave Analysis and Homogenization of Spatiotemporally Modulated Wire Medium

Michael Kreiczer and Yakir Hadad  The school of Electrical Engineering, Ramat-
Aviv, Tel-Aviv University, Israel, 69978<EMAIL_ADDRESS>
###### Abstract

In this paper we develop homogenization theory for spatiotemporally modulated
wire medium. We first solve for the modal waves that are supported by this
composite medium, we show peculiar properties such as extraordinary waves that
propagate at frequencies below the cut-off frequency of the corresponding
stationary medium. We explain how these unique solutions give rise to an
extreme Fresnel drag that exists already with weak and slow spatiotemporal
modulation. Next, we turn to derive the effective material permittivity that
corresponds to each of the first few supported modes, and write the average
fields and Poynting’s vector. Nonlocality, nonreciprocity, and anisotropy due
to the spatiotemporal modulation direction, are three inherent properties of
this medium, and are clearly seen in the effective material parameters. As a
figure of merit, we also derive the effective permittivity of a plasma medium
with spatiotemporally modulated plasma frequency. This comparison is
interesting since the plasma medium can be considered as the effective medium
that is obtained by a stationary wire medium. We validate that homogenization
and spatiotemporal variation are not necessarily interchangeable operations.
And indeed, in certain parameter regimes the homogenization should be
performed directly on spatiotemporally modulated composite medium, rather than
first homogenize the stationary medium and then phenomenologically introduce
the effect of the space-time modulation.

††preprint: APS/123-QED

##  I Introduction

The study electromagnetic wave dynamics in time-varying media goes back
several decades ago [1, 2]. Recently, there is a new blooming of research in
this direction especially in the context of new developments toward the next
generation in the research of metamaterials. Such a progress is accompanied by
new-technological abilities that make the actual realization of engineered
time-varying media plausible. The benefits from removing the restriction
caused by time invariance, and adding another dimension for engineering - the
time, has opened a whole new area of technological possibilities. Linear time
variant (LTV) metamaterials have been recently used for various applications,
such as: inverse prism, frequency conversions, temporal band gap, magnetless
nonreciprocity, time reversal and effective permittivity realization by
temporal switching [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19]. Moreover, since many of the known physical bounds on wave engineering are
based on the LTI assumption, by enabling the wave network to be time-varying
we may expect overcoming these physical LTI bounds [20, 21]. See for example
[22] for the Bode-Fano bound for impedance matching for LTV microwave circuit,
and [23] for overcoming the bandwidth bound in matching small antennas.

On a parallel route. Effective medium theory which is classical in
electromagnetics and dates back to the 1904 work of Maxwell Garnett [24] has
become a cornerstone in the context of metamaterials. It is an essential tool
to model and analyze effective media with various types of meta-atoms [25,
26]. The effective medium approach, also termed homogenization, allows us to
transform from the composite material, and its microscopic structure to a new
equivalent material having, macroscopic, smoothed properties. A particularly
interesting example of effective medium is the wire medium of its various
types. In its simplest form, the wire medium consists of a two dimensional
lattice of infinitely long perfectly conducting wires surrounded by an
homogeneous dielectric material. Despite its apparent simplicity that enables
a relatively easy to follow analytical modelling, interestingly, the wire
medium gives rise to various peculiar properties such as strong anisotropy,
effective plasma-like dispersion, as well as strong-nonlocality [27, 28, 29,
30, 31, 32, 33, 34, 35]. Various extensions to the basic wire medium have been
studied, including loaded-wire medium that is relatively straightforward to
analyze with various loading types [36].

In this paper we study electromagnetic waves that are propagating in a
capacitively-loaded wire medium that undergo a spatiotemporal modulation. See
illustration in Fig. 1. Our goal in this work is two fold. First, taking into
account the complete interaction between the spatiotemporally modulated loaded
wires, we explore the various wave phenomena that can be supported by this
medium. Specifically, we study ($i$) non-reciprocal and anisotropic
propagation due to the spatiotemporal modulation, ($ii$) the existence of
extraordinary waves that propagate at frequencies that are below the so-called
plasma frequency of the stationary wire medium (that is, in the absence of
modulation), even with negligible modulation index, and thus are akin to the
Whistler mode in magnetized plasma [37], and consequently ($iii$) the
emergence of Fresnel drag [38] already with weak and slow modulation. The
second goal, is to derive effective medium properties when the meta-atoms (the
loaded wires) are time-modulated. In particular we explore the effect of
taking into account the time-modulation already in the homogenization process,
as opposed to first homogenize the stationary medium and later introducing the
time-modulation into the effective properties. In other words, we try to
address what is the effect of the spatiotemporal modulation on the
homogenization. Related to that, is the issue of dispersion in the presence of
time modulation, this has been recently tackled in several works such as [39,
40]. Here, we provide a study for the specific wire-medium case using a quasi
frequency-domain approach, and take the complete inter-wire interaction that
comprises the composite medium, as discussed below. We pay a particular
attention to the wave dynamics near the plasma frequency, there, due to the
sensitivity of the zero crossing point of the effective transverse
permittivity of the stationary wire medium, interesting phenomena pop up as
soon as we introduce weak spatiotemporal modulation.

The paper is arranged as follow. For the sake of its self-continency, in Sec.
II, we open the paper with a brief review on stationary capacitively loaded
wire media. Then, in Sec. III we introduce spatiotemporal modulation to the
capacitive loads and and provide a through mathematical formulation to solve
the modal problem. In Sec. IV, focusing on narrow frequency range near the
plasma-frequency of the corresponding stationary lattice, we derive analytic
expressions for the dispersion relations of the first lowest order modes
supported in the lattice, and study the wave dynamics in this parameters
domain. Later, in Sec. V, we derive the effective medium parameters of several
supported modes, in addition we derive expressions for the averaged electric
and magnetic fields, and the corresponding Poynting’s vector in this case.
Finally, in Sec. VI, we compare between the wave phenomena in the time-
modulated wire medium, and the wave phenomena predicted when the time-
modulation is introduced into the effective medium properties of the
stationary wire-medium, i.e., after homogenization.

##  II Stationary capacitively loaded wire medium

In the following sections we turn to analyze wave phenomena in
spatiotemporally modulated loaded-wire media. Unavoidably, we will frequently
refer to the derivation and the main results of the corresponding, stationary,
case. Therefore, to make the paper self-contained, in this section we briefly
overview _stationary_ , capacitively loaded wire medium. We mainly follow the
analysis in [36] for general wire loading, and stress the more relevant
aspects for the spatiotemporal problem to follow.

Figure 1: (a) Capacitively loaded wire. (b) The wire medium, on a rectangular
lattice with unit-cell dimensions $a$ along $x$ and $b$ along $y$. In this
work we focus on wave that propagate transverse to the wires with
$\vec{q}=q_{x}\hat{x}+q_{y}\hat{y}$.

Assume that a $\hat{z}$ polarized electromagnetic wave is propagating inside
an infinite loaded wire medium as shown in Fig. 1. We assume that the wires
are surrounded by vacuum with permittivity and permeability $\epsilon_{0}$ and
$\mu_{0}$, respectively. The lattice points are given by:
${\vec{R}}_{m,l}=ma\hat{x}+lb\hat{y}$. The induced current on each of the
wires is given by $\vec{I}=\alpha\vec{E}^{\mbox{\small loc}}$ where
$\vec{E}^{\mbox{\small loc}}$ is the local electric field, namely, the
electric field at the wire location but in the absence of the wire itself, and
$\alpha$ is the wireś susceptibility. If the wires are periodically loaded by
lumped impedance $Z_{L}$, with periodicity, $\Delta\ll\lambda$, the inverse
susceptibility satisfies [36],

${\alpha^{-1}}\left(\omega\right)={\alpha_{0}^{-1}}\left(\omega\right)+\frac{Z_{L}}{\Delta},$
(1)

where ${\alpha_{0}^{-1}}\left(\omega\right)=\frac{{\eta
k}}{4}H_{0}^{\left(2\right)}\left({k{r_{0}}}\right)$. For serial loading of
lumped capacitors $Z_{L}/\Delta=1/(j\omega\tilde{C}_{0})$ and
$\tilde{C}_{0}=C_{0}\Delta$. Here $\eta=120\pi\Omega$, and $k=\omega/c$ are
the free space impedance and wavenumber, respectively, $\omega$ is the radial
frequency, $c$ is the speed of light in vacuum, $r_{0}$ is the wire radius,
and $H_{0}^{(2)}$ denotes the zero order Hankel function of the second type.
In light of the lattice periodicity the Floquet-Bloch wave solution is
expected,

${I_{m,l}}={I_{0}}{e^{-j\left({{q_{x}}am+{q_{y}}bl}\right)}}$ (2)

where $q_{x},q_{y}$ stand for the transverse (to $z$) components of the
wavenumber in the wire medium. Note that throughout this paper we assume for
simplicity that $q_{z}=0$ which is in accord with our initial assumption that
the electric field is polarized solely along $z$. Under the assumptions above,
the wavenumber dispersion is found by solving,

$E_{z}^{\mbox{\small loc}}\left(I_{0}\right)=\alpha^{-1}(\omega)I_{0}$ (3)

where $E_{z}^{\mbox{\small loc}}$ denotes the local field acting on the wire
at $m=l=0$. Following [36] it is given by,

$E_{z}^{\mbox{\small loc}}=-I_{0}\frac{\eta k}{2}\left(S_{0}+\sum_{m\neq
0}S_{m}e^{-jq_{x}am}\right),$ (4)

with,

$\displaystyle{S_{0}}=\sum\limits_{l=1}^{\infty}{H_{0}^{\left(2\right)}\left({k{R_{0l}}}\right)\cos{q_{y}}bl}$
(5a)
$\displaystyle{S_{m}}=\frac{1}{2}\sum\limits_{l=-\infty}^{\infty}{H_{0}^{\left(2\right)}\left({k{R_{ml}}}\right){e^{-j{q_{y}}bl}}}.$
(5b)

and $R_{ml}=\sqrt{(ma)^{2}+(lb)^{2}}$. The infinite slowly converging series
in Eqs. (5a) and (5b) may be converted into fast converging series using e.g.,
Poisson summation [36, 41]. For a dense grid, $ka,kb\ll 1$, the dispersion
relation may be approximated by

$\omega\tilde{L}_{w}-\frac{1}{\omega\tilde{C}_{0}}=\frac{k\eta}{ab(k^{2}-q^{2})}$
(6)

where $q^{2}=q_{x}^{2}+q_{y}^{2}$ and

$\tilde{L}_{w}=\frac{\eta}{2c}\left(\frac{1}{\pi}\ln\frac{b}{2\pi
r_{0}}\\!+\\!\sum_{l=1}^{\infty}\frac{\coth(\pi al/b)-1}{\pi
l}\\!+\\!\frac{a}{6b}\right)$ (7)

denotes the per unit length intrinsic wire inductance. We note that for the
dense grid case, the dispersion equation depends only on
$q=\sqrt{q_{x}^{2}+q_{y}^{2}}$. Then, an equation for the effective
permittivity along the longitudinal (to $z$) axis is obtained by
$\varepsilon_{z}=\frac{q^{2}}{k^{2}}$. This yields

${\varepsilon_{z}}\left(\omega\right)={\varepsilon_{0}}\left[{1-\frac{{k_{p}^{2}{\left(k\right)}}}{{{k^{2}}}}}\right].$
(8)

with

$k_{p}^{2}(k)=\frac{\eta k}{ab(\omega\tilde{L}_{w}-1/\omega\tilde{C}_{0})}.$
(9)

where $k=\omega/c$, and $k_{p}$ is termed here and henceforth, the plasma
frequency. Nevertheless, note that unlike the unloaded wire media, for which
effectively $\tilde{C_{0}}\rightarrow\infty$ and consequently
$k_{p}^{2}=(\eta/c)/ab\tilde{L}_{w}=(1/\tilde{C}_{b}\tilde{L}_{w})/c^{2}$ with
$\tilde{C}_{b}=ab\epsilon_{0}$ is frequency independent, when the wires are
loaded the behaviour is not of a typical plasma since $k_{p}$ depends on $k$.
In the following we shall term $\tilde{C}_{b}$ \- the per-unit-length
background capacitance.

##  III Wire media with spatiotemporally modulated capacitive loading

Our first goal in this paper is to study wave propagation in a
spatiotemporally modulated wire media. To that end, we generalize the
conventional approach for stationary LTI wire media that was described above.
Our analysis is based on the concept of harmonic balance that can be regarded
as a _quasi_ -frequency domain technique. In this sense, our method is less
general then, e.g., [42] that suggests, in principle, a direct time-domain
venue to that problem. However, since we focus on the particular case of small
time-harmonic capacitance perturbation, we find this methodology, as discussed
below, to be more direct and therefore a better fit for this problem.

###  III.1 The response of a single wire

Assume that the loading capacitance on each of the wires is modulated as,

$C\left(t\right)={C_{0}}+\delta C\cos\left({\Omega
t-\varphi}\right)\,\,,\,\,m=\frac{{\delta C}}{{{C_{0}}}}\ll 1.$ (10)

Here and henceforth, $\Omega$ denotes the modulation frequency and $m$ the
modulation depth. In the absence of time modulation, the time-domain
counterpart of Eq. (3) is given by $\widehat{{E^{\mbox{\small
loc}}}}\left(t\right)=\widehat{\alpha^{-1}}(t)*\widehat{I}\left(t\right)$
where ‘wide-hat’ $\widehat{X}(t)$ represents the time-domain, inverse-Fourier
transform, counterpart of $X(\omega)$, and $X$ stands for $E^{\mbox{\small
loc}}$, $\alpha^{-1}$, and $I$. This time-domain perspective enables us to
naturally introduce the effect of the time-modulation on the capacitors (see
appendix VIII.2)

$\widehat{E^{\mbox{\small
loc}}}\left(t\right)=\widehat{\alpha_{0}^{-1}}(t)*\widehat{I}\left(t\right)+\frac{1}{{\tilde{C}\left(t\right)}}\int\limits_{-\infty}^{t}{I\left(\tau\right)d\tau}.$
(11)

with $\tilde{C}(t)=C(t)\Delta$. By using ${1}/(1+x)\approx 1-x$ for $x\ll 1$,
inserting Eq. (10) in (11), and transforming to the frequency domain we get,

$\begin{array}[]{l}{E^{\mbox{\small
loc}}}\left(\omega\right)=\,\left[{{\alpha_{0}^{-1}}\left(\omega\right)+\frac{1}{{j\omega{\tilde{C}_{0}}}}}\right]I\left(\omega\right)\\\
\\\
\\!-\frac{{m}}{{2}}\left[{\frac{{{e^{-j\varphi}}}}{{j\left({\omega-\Omega}\right){\tilde{C}_{0}}}}I\left({\omega-\Omega}\right)+\frac{{e^{j\varphi}}}{{j\left({\omega+\Omega}\right){\tilde{C}_{0}}}}}{I\left({\omega+\Omega}\right)}\right].\end{array}$
(12)

Consider now a particular nominal frequency $\omega_{0}$, which may be, for
instance, the excitation frequency of a monochromatic impinging wave, or the
excitation frequency of a localized source within the bulk. Then, in light of
the harmonic time-modulation the local field has the following frequency
dependence,

$E^{\mbox{\small
loc}}\left(\omega\right)=\sum\limits_{n=-\infty}^{\infty}E_{n}^{\mbox{\small
loc}}\delta\left({\omega-{\omega_{n}}}\right),$ (13)

with $\omega_{n}=\omega_{0}+n\Omega$, and $\delta\left(\cdot\right)$ denotes
Dirac’s delta. Then, for the current we have,

$I\left(\omega\right)=\sum\limits_{n=-\infty}^{\infty}{{I_{n}}\delta\left({\omega-{\omega_{n}}}\right).}$
(14)

Here, ${E_{n}^{\mbox{\small loc}}}$ and $I_{n}$ denote the $n$’th harmonic
amplitude of the local field and the induced current, respectively. By
plugging Eq. (13) and Eq. (14) into Eq. (12), and by balancing the
coefficients of equal harmonics we find,

$\begin{array}[]{l}E_{n}^{\mbox{\small
loc}}=\left[{{\alpha_{0}^{-1}}\left({{\omega_{n}}}\right)+\frac{1}{{j{{\omega_{n}}}{\tilde{C}_{0}}}}}\right]{I_{n}}\\\
\\\
-\frac{m}{2}\left[\frac{e^{-j\varphi}}{j\omega_{n-1}\tilde{C}_{0}}I_{n-1}+\frac{e^{j\varphi}}{j\omega_{n+1}\tilde{C}_{0}}I_{n+1}\right].\\\
\end{array}$ (15)

###  III.2 The time-modulated lattice

Our goal is to analyze a spatiotemporally modulated wire media. While the
temporal modulation is introduced by the modulation frequency parameter,
$\Omega$, the space modulation is introduced by setting the phase $\varphi$ in
Eq. (10) for each of the wires. To that end, we introduce the vector
$\vec{\zeta}=\hat{x}\zeta\cos\xi+\hat{y}\zeta\sin\xi$ in the $xy$ plane, that
represents the direction and magnitude at which the capacitors phase $\varphi$
is accumulated. Thus, in order to achieve the effect of synthetic motion, we
choose the phase of the capacitance of the wire with indexes $(m,l)$ to be
${\varphi_{m,l}}=\vec{\zeta}\cdot{\vec{R}}_{m,l}=ma\zeta\cos\xi+bl\zeta\sin\xi$.
By controlling the phase modulation direction, spatial dispersion is obtained
as shown below. Spatial dispersion has been demonstrated by spatiotemporally
modulated metasurface has been shown in [43, 44]. Substantial spatial
dispersion may be achieved also in different means e.g. in [45] where the
authors utilized dielectric rods with different diameter instead of a constant
loaded wire.

Due to Floquet-Bloch theorem, the currents take the form,

${I_{m,l}}=\sum\limits_{n=-\infty}^{\infty}{{A_{n}}{e^{-j{{\vec{q}}_{n}}\cdot{\vec{R}}_{m,l}}}\delta\left({\omega-{\omega_{n}}}\right)}$
(16)

where $n$ represents the temporal harmonic number, and the corresponding wave
vector is given by (see appendix VIII.3 for derivation),

${{\vec{q}}_{n}}={{\vec{q}}_{0}}+n\vec{\zeta}$ (17)

with
${{\vec{q}}_{0}}=\hat{x}{q_{0,x}}+\hat{y}{q_{0,y}}=\hat{x}{q_{0}}\cos{\theta_{0}}+\hat{y}{q_{0}}\sin{\theta_{0}}$.
Note that ${{\vec{q}}_{0}}$ represents the corresponding wave vector at the
nominal frequency $\omega_{0}$. We can now use the expression for the local
field, and derive the dispersion relation and the equations for the
eigenmodes. We follow the approach applied for stationary medium in Sec. (II),
and perform the necessary modifications due to the spatiotemporal modulation
of the medium. Assuming an infinite lattice, with currents given by Eq. (16),
the local field on the $(m,l)=(0,0)$ wire reads,

$\begin{array}[]{l}{E^{\mbox{\small
loc}}}=\sum\limits_{n=-\infty}^{\infty}{E_{n}^{\mbox{\small
loc}}}\delta\left({\omega-{\omega_{n}}}\right)\\\
\,=-\sum\limits_{n=-\infty}^{\infty}{\sum\limits_{m,l\neq
0,0}{{A_{n}}G_{nml}{e^{-j{{\vec{q}}_{n}}\cdot{\vec{R}}_{m,l}}}\delta\left({\omega-{\omega_{n}}}\right)}}\end{array}$
(18)

where
$G_{nml}=({{\eta{k_{n}}}}/{4})H_{0}^{\left(2\right)}\left({{k_{n}}{R_{ml}}}\right)$
with $R_{ml}=\sqrt{(ma)^{2}+(lb)^{2}}$ and $k_{n}=\omega_{n}/c$, denoting the
two dimensional Green’s function at the n’th harmonic for a source located at
$\vec{R}_{ml}$ and observer at $\vec{R}_{00}$. By applying Eq. (18) with Eq.
(12), and balancing between equal harmonics at the two sides of the resulting
in equation, we find a recursive relation between the different temporal
harmonics of the induced current on the $(m,l)=(0,0)$ wire,

$a_{n}A_{n-1}+b_{n}A_{n}+c_{n}A_{n+1}=0,$ (19)

where,

$\displaystyle{a_{n}}=\frac{m}{2}\frac{1}{{{\omega_{n-1}\tilde{C}_{0}}}},$
(20a) $\displaystyle\begin{array}[]{l}{b_{n}}=\frac{\eta
k_{n}}{2}\left[\frac{1}{\pi}\ln\frac{b}{2\pi
r_{0}}+\frac{1}{b\beta_{n,x}^{(0)}}\frac{\sin\beta_{n,x}^{(0)}a}{\cos\beta_{n,x}^{(0)}-\cos
q_{n,x}a}+\right.\\\ \left.\sum_{l\neq
0}{\left(\frac{1}{b\beta_{n,x}^{(l)}}\frac{\sin\beta_{n,x}^{(l)}a}{\cos\beta_{n,x}^{(l)}-\cos
q_{n,x}a}-\frac{1}{2\pi|l|}\right)}-\frac{1}{\omega_{n}\tilde{C}_{0}}\right],\end{array}$
(20d) $\displaystyle
c_{n}=\frac{m}{2}\frac{1}{{{\omega_{n+1}\tilde{C}_{0}}}},$ (20e)

with $\beta_{n,x}^{(l)}=-j\sqrt{\left(q_{n,y}+\frac{2\pi
l}{b}\right)^{2}-k_{n}^{2}},\,\mbox{Re}\left\\{\sqrt{\cdot}\right\\}\\!>\\!0$.
For a dense grid, as in Sec. II, we approximate $b_{n}$ in Eq. (20d),

$b_{n}=\omega_{n}\tilde{L}-\frac{1}{\omega_{n}\tilde{C}_{0}}-\frac{\eta
k_{n}}{ab(k_{n}^{2}-q_{n}^{2})}$ (21)

with,

$q_{n}^{2}={\left({{{\vec{q}}_{0}}+n\vec{\zeta}}\right)^{2}}=q_{0}^{2}+2n\zeta{q_{0}}\cos\gamma+{n^{2}}{\zeta^{2}}.$
(22)

Here, $(\vec{x})^{2}=\vec{x}\cdot\vec{x}$, and $\gamma=\theta_{0}-\xi$. In Eq.
(21), $\tilde{C}_{0}=C_{0}\Delta$ as before and

$\tilde{L}=\tilde{L}_{w}+\tilde{L}_{0}\mbox{ with
}\tilde{L}_{0}=\frac{L_{0}}{\Delta}$ (23)

where $L_{0}$ is an inductance that may be connected in series with the
loading capacitor $C(t)$. Effectively it simply increases the intrinsic wire
inductance. Note that as in the case of stationary wire medium with dense
grid, also here the approximation for $b_{n}$ in Eq. (21) depends only on
$q_{n}^{2}$, whereas the exact expression for $b_{n}$ in Eq. (20d) depends
also on $q_{n,x}$.

Figure 2: (a) The time-modulated loaded wire. $L_{0}$ increases the self
inductance of the wire, and $C(t)$ introduces the time-modulated capacitance.
These lumped elements are periodically loaded on the wire with periodicity
$\Delta$, thus creating the effective wire time-modulated susceptibility. (b)
The wire medium. Each wire is modulated with the same temporal frequency
$\Omega$, but with a different phase that is determined by the modulation wave
vector $\vec{\zeta}$. The propagation takes place with a fundamental wave
vector $\vec{q}_{0}$.

As opposed to the stationary loaded wire medium case where the effective
longitudinal permitivity may be analytically expressed using Eq. (8) with Eq.
(9), for the spatiotemporally modulated case, the interaction between the
different temporal frequency harmonics should be included. To that end, Eq.
(19) is used. The latter can be represented using a tridiagonal matrix of
infinite rank. We assume that for large enough $N$, the currents
$A_{-\left(N+1\right)},A_{N+1}$ are negligible. In this case the infinite
matrix can be approximated by a finite square matrix

${\underline{\underline{A}}^{N}}=\left[{\begin{array}[]{*{20}{c}}\ddots&\ddots&0&0&{}\hfil&{}\hfil&{}\hfil\\\
\ddots&{{b_{-2}}}&{{c_{-2}}}&0&0&{}\hfil&{}\hfil\\\
0&{{a_{-1}}}&{{b_{-1}}}&{{c_{-1}}}&0&0&{}\hfil\\\
0&0&{{a_{0}}}&{{b_{0}}}&{{c_{0}}}&0&0\\\
{}\hfil&0&0&{{a_{1}}}&{{b_{1}}}&{{c_{1}}}&0\\\
{}\hfil&{}\hfil&0&0&{{a_{2}}}&{{b_{2}}}&\ddots\\\
{}\hfil&{}\hfil&{}\hfil&0&0&\ddots&\ddots\end{array}}\right].$ (24)

For nontrivial wave solutions we look for the dispersion relation
$\vec{q}_{0}(\omega)$ that nullifies the determinant,

$|\underline{\underline{A}}^{N}|=0.$ (25)

In the following we explore the wave dynamics in the bulk as obtained from
this dispersion relation for a particularly interesting case near the plasma
frequency of the corresponding stationary wire medium.

##  IV Near the plasma frequency of the stationary medium

In the following we explore the wave dynamics of the spatiotemporally
modulated medium with parameters that are near the plasma frequency of the
corresponding stationary medium. Namely, with
$\omega=\omega^{0}+\delta\omega$, where $\omega^{0}$ is the frequency that
nullifies $\epsilon_{z}(\omega)$ in Eq. (8). Thus

$\omega^{0}=\sqrt{\frac{1}{\tilde{L}}\left(\frac{1}{\tilde{C}_{0}}+\frac{1}{\tilde{C}_{b}}\right)}=\sqrt{\frac{1}{{\tilde{L}\tilde{C}_{0}}}}\sqrt{1+\psi}$
(26)

where $C_{b}$ is defined after Eq. (9), and

$\psi=\frac{\tilde{C}_{b}}{\tilde{C}_{0}}.$ (27)

In the following we will also use $k^{0}=\omega^{0}/c$.

###  IV.1 Analytic approximation for the lower order solutions with weak
modulation

Assuming weak modulation, $m\ll 1$, and as a result the dominant harmonics
beside the fundamental one will be $n=\pm 1$. Note that this assumption
excludes the possibility that higher order harmonic will be strongly excited
in this case due to the presence of some resonance mechanism at some other
frequency. Under this assumption the infinite rank matrix
$\underline{\underline{A}}$ will be replaced by a tridiagonal $3\times 3$
matrix $\underline{\underline{A}}^{1}$ ($N=1$ in Eq. (25)). The dispersion
relation in this case is found by nullifying its determinant, i.e.,

$\frac{{{a_{0}}{c_{-1}}}}{{{b_{-1}}}}+\frac{{{a_{1}}{c_{0}}}}{{{b_{1}}}}={b_{0}}.$
(28)

By plugging the expressions for $a_{n}$ and $c_{n}$ that are given in Eq.
(20), with the dense grid approximation Eq. (21), together with Eq. (27), and
with $q_{1},q_{-1}$ given in Eq. (17), the non-trivial solution requirement in
Eq. (28) turns into a sixth order polynomial equation for $q_{0}$,

$\frac{m^{2}\psi^{2}}{4}\left[\frac{1}{X_{-1}(\omega,q_{0})}+\frac{1}{X_{1}(\omega,q_{0})}\right]=X_{0}(\omega,q_{0})$
(29a) with
$X_{n}(\omega,q_{0})=\frac{k_{n}^{2}}{\left(k^{0}\right)^{2}}\left(1+\psi\right)-\psi-\frac{k_{n}^{2}}{k_{n}^{2}-q_{n}^{2}}.$
(29b)

Let us the denote the $n$’th root of Eq. (29a) by $q_{0}^{\\#n}$ where
$n=1..6$. Two of the six solutions of this polynomial equation will represent
_a perturbation over the two solutions of the corresponding stationary wire
medium_. Since we focus here on the dispersion near the plasma frequency of
the corresponding stationary medium, i.e., in the vicinity of $\omega^{0}$, it
is reasonable to assume that two of the solutions satisfy
${{q_{0}^{\\#1}}},\,{{q_{0}^{\\#2}}}\ll k^{0}$. In that case, in the left-
hand-side of Eq. (29a) we will use $\omega=\omega^{0}$, and $q_{0}=0$, and
thus $X_{1},X_{-1}\sim 1$. Since $m\ll 1$, the plasma frequency of the
spatiotemporally modulated medium, namely, the frequency at which
$q_{0}^{\\#1,2}$ experience the transition from being purely imaginary to
purely real shifts by $\delta\hat{\omega}$ with respect to $\omega^{0}$ of the
stationary medium where,

$\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}\approx\left[\frac{1}{X_{1}(\omega^{0},0)}+\frac{1}{X_{-1}(\omega^{0},0)}\right]\,\frac{{{m^{2}}{\psi^{2}}}}{{8\left({1+\psi}\right)}}.$
(30)

Then, the first two solutions $q_{0}^{\\#1,2}$ for the fundamental modes, up
to the shift in $\delta\hat{\omega}$ are obtained by solving
$X_{0}(\omega,q_{0})=0$ and are given by,

$\begin{array}[]{l}q_{0}^{\\#1,2}=\pm{k^{0}}\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}}\right)\times\\\
\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\sqrt{1-\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}}\right)}^{2}}\left({1+\psi}\right)-\psi}}}.\end{array}$
(31)

We can see that under the dense grid approximation, solutions ${\\#1,\\#2}$
are isotropic, independent of the relative angle $\gamma$ between
spatiotemporal modulation axis and the direction of the wave propagation.

To find the other four roots of Eq. (29a) we shall assume that these solutions
are not small compare to $k^{0}$ near $\omega^{0}$. As a result, the right-
hand-side of Eq. (29a) $X_{0}(\omega^{0},q_{0})\sim 1$. On the other hand, in
light of the small factor $m^{2}$ in the left-hand-side, the only way to
balance the two sides is if either $X_{1}$ or $X_{-1}$ are nearly vanishing.
Therefore the remaining four roots are approximated by solving
$X_{1}(\omega,q_{0})=0$ and $X_{-1}(\omega,q_{0})=0$. Each one of the two
equations yields two solutions,

$\displaystyle\begin{array}[]{l}q_{0}^{\\#3,4}=-\zeta\cos\gamma\pm\\\ \\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,{k_{1}}\sqrt{-{{\left({\frac{\zeta}{{{k_{1}}}}}\right)}^{2}}{{\sin}^{2}}\gamma+1-\frac{1}{{\frac{{k_{1}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}\left({1+\psi}\right)-\psi}}}\end{array}$
(32d) $\displaystyle\begin{array}[]{l}q_{0}^{\\#5,6}=\zeta\cos\gamma\pm\\\ \\\
\,\,\,\,\,\,\,\,\,\,\,\,\,{k_{-1}}\sqrt{-{{\left({\frac{\zeta}{{{k_{-1}}}}}\right)}^{2}}{{\sin}^{2}}\gamma+1-\frac{1}{{\frac{{k_{-1}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}\left({1+\psi}\right)-\psi}}}\end{array}$
(32h)

If imaginary, the square roots in Eqs. (31) and (32) should be chosen to
guarantee physical solutions, namely, the wave must decay at infinity. In Fig.
3 we show the dispersion of the six solutions as obtained by the analytical
approximated relations in Eqs. (31) and (32), the dispersion is shown for the
following structure parameters
$a=0.07\lambda^{0},b=0.1\lambda^{0},r_{0}=0.0001\lambda^{0}$, and
$\theta_{0}=0$. As opposed to the structural parameters, the _modulation_
parameters differ between four typical cases that are summarized in the table
below.

| $m=\delta C/C_{0}$ | $\Omega/\omega^{0}$ | $\zeta/k^{0}$ | $\tilde{L}_{0}/\tilde{L}_{w}$ | $\xi$ | $\psi$
---|---|---|---|---|---|---
Case I | 0.1 | 0.15 | 0.1 | 4.3 | $-\pi/3$ | 0.27
Case II | 0.35 | 0.2 | 1.0 | 4.0 | 0 | 0.2
Case III | 0.06 | 0.25 | 1.2 | 7.6 | $\pi/20$ | 1.07
Case IV | 0.25 | 0.3 | 1.0 | 9.3 | $-\pi$ | 1.47
Table 1: Four parameter cases used in the numerical calculations in this
section.

The complex dispersion diagrams, real and imaginary $q_{0}$, are shown as a
function of $\delta\omega$ which is the detuning frequency from $\omega^{0}$
\- the plasma frequency of the stationary lattice. In all four cases the lines
are color-coded, blue for $q_{0}^{\\#1,2}$ \- the perturbation on the solution
of the stationary medium, whereas green and red lines are used for the new
solutions due to the spatiotemporal modulation, $q_{0}^{\\#3,4}$ and
$q_{0}^{\\#5,6}$, respectively. The colored lines in Fig. 3 are obtained by
using the approximated solutions that are given in Eqs. (31) and (32d,32h), to
show the analytic approximation validity, as a comparison in the figures we
also draw by black circles the numerical solution of Eq. (29a).

Figure 3: Dispersion relation as obtained by solutions #1-#6 in Eqs.
(31)-(32). The fundamental wavenumber $q_{0}$ of each one of the solutions is
shown as a function of the frequency $\delta\omega/\omega_{0}$ for the four
parameter cases that are tabulated in Table 1. (a) Case I, (b) Case II, (c)
Case III, and (d) Case IV.

The good agreement between the analytical approximations, and the numerical
calculations, allow us to use our analytical results for synthesis. For
example, an interesting characteristic for the roots of
$q_{0}^{\\#3},q_{0}^{\\#4},q_{0}^{\\#5}$ and $q_{0}^{\\#6}$, is that we can
create a different response for different angle $\gamma$ and in particular
reach effective Fresnel drag dispersion accompanied with strong anisotropy for
waves propagating in the wire medium bulk by using _slow_ spatiotemporal
modulation.

###  IV.2 Effective Fresnel drag with weak and slow spatiotemporal modulation

The spatiotemporal modulation enables to reach a dispersion that shares
similarities with the dispersion of waves in a moving medium. However,
typically this effect at the fundamental harmonic frequency is rather weak and
requires time modulation much faster then the guided signal frequency [38].
This is not the case however when we spatiotemporally modulating the loaded
wire medium. In this section we show that it is the _interplay_ between the
_spatiotemporal modulation_ and the _cutoff of the fundamental harmonic_ below
$\omega^{0}$ \- the plasma frequency of the corresponding stationary medium -
that enables to obtain effectively Fresnel drag with no actual motion and
using slow and weak modulation. To see this effect, we use Eq. (31), and Eqs.
(32d,32h) in order to extract the relation between $\delta\omega$ and the
transverse wave number components $q_{0x}=q_{0}\cos\theta_{0}$ and
$q_{0y}=q_{0}\sin\theta_{0}$ for a given set of medium and modulation
parameters. By solving Eq. (31) for $\delta\omega$ we obtain,

${\frac{{\delta\omega}}{{{\omega^{0}}}}^{\\#1,2}}=\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}+\frac{{\frac{{q_{0}^{2}}}{{2{{\left({{k^{0}}}\right)}^{2}}}}}}{{\left[{1+\psi\left({1-\frac{{q_{0}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}}\right)\,}\right]}}$
(33)

whereas by solving Eqs. (32d) and (32h) for $\delta\omega$ we get,

${\frac{{\delta\omega}}{{{\omega^{0}}}}^{\pm
1}}=\frac{{1+\left({1-\frac{\psi}{{1+\psi}}}\right)\frac{{q_{\pm
1}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}-{{\left({1\pm\frac{\Omega}{{{\omega^{0}}}}}\right)}^{2}}}}{{4{{\left({1\pm\frac{\Omega}{{{\omega^{0}}}}}\right)}^{2}}-2\left[{1+\frac{{q_{\pm
1}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}}\right]}}$ (34)

where the $+1$ [$-1$] corresponds to solutions #3 and #4 [#5 and #6] in Eqs.
(32), and $q_{\pm 1}$ is given using Eq. (22). In the latter we note that
$q_{0}^{2}=q_{x}^{2}+q_{y}^{2}$, $\gamma=\theta_{0}-\xi$ where $\xi$ is the
direction of the modulation and $\tan\theta_{0}=q_{y}/q_{x}$.

Figure 4: Dispersion diagram of the fundamental solutions #1 and #2 and high
order solutions #3-#4. (a) With $\xi=0$ the fundamental harmonic dispersion is
centered at the origin of the ($q_{x},q_{y}$) plane, however, its cutoff
frequency is $\delta\omega/\omega_{0}\approx 0$. (b) As opposed to that, the
high order solutions are propagating already at negative frequencies. These
solutions are essentially nonreciprocal in light of their asymmetric
dispersion. (c) and (d) are as (a) and (b) but with $\xi=\pi/4$.

These dispersion relations are shown in Fig. 4 below. In (a) with $\xi=0$, a
color map that shows $\delta\omega/\omega_{0}$ as a function of the transverse
wave vector $(q_{x},q_{y})$ is presented, showing that the fundamental
solutions are propagating only at positive frequencies (to be more precise
$\frac{{\delta\omega}}{{{\omega^{0}}}}^{\\#1,2}=\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}\approx
0$). However, the $\pm 1$ harmonics (see Eq. (34)), with strongly _asymmetric_
dispersion, are propagating at _negative_ frequencies, as shown in Fig. 4(b).
As a result, in the spatiotemporally modulated wire medium effective motion
may be exhibited at the first propagating harmonic and with slow modulation.
This behaviour may be controlled by adjustments of the direction of
modulation, as shown in Fig. 4(c) for the fundamental harmonic, and in (d) for
the $\pm 1$ harmonic, this time with spatiotemporal modulation along the
$\xi=45\deg$ direction.

##  V Homogenization

###  V.1 Effective modal permittivity

Using the wavenumber dispersion of the modal solutions given in Eqs. (31)-(32)
it is possible to define effective permittivity for each one of the modal
solutions found above. To that end we first write free space wavenumber as

${{k^{2}}}={\left({{k^{0}}}\right)^{2}}{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)^{2}}.$
(35)

Then, the homogenized effective permittivity of the modal solution #$i$ is
given by

$\varepsilon_{r:i}={\left({q_{0}^{\\#i}}\right)^{2}}/{k^{2}}.$ (36)

Note that this effective relative permittivity connects the free space wave
number with the _fundamental_ guided mode wavenumber of each one of the modal
solutions. It is important to stress that each one of these modes, with
$q_{0}^{\\#i}$, consists of infinite space-time harmonics as dictated by Eq.
(16). Once $\varepsilon_{r:i}$ is known, $q_{0}^{\\#i}$ can be calculated.
Then, using Eq. (17), the wavenumber of each one of the infinite harmonics
that consists the $i$ modal solution is immediately found. Specifically, for
the first six fundamental modal solutions in the spatioteporally modulated
wire medium we have,

$\displaystyle{\varepsilon_{r:1}}={\varepsilon_{r:2}}={\left({\frac{{1+\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}}}{{1+\frac{{\delta\omega}}{{{\omega^{0}}}}}}}\right)^{2}}\left[{1-\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}}\right)}^{2}}\left({1+\psi}\right)-\psi}}}\right]$
(37a)
$\displaystyle{\varepsilon_{r:3,4}}=\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}}}{\left[{-\frac{\zeta}{{{k^{0}}}}\cos\gamma\pm\frac{{{k_{1}}}}{{{k^{0}}}}\sqrt{-{{\left({\frac{\zeta}{{{k_{1}}}}}\right)}^{2}}{{\sin}^{2}}\gamma+1-\frac{1}{{\frac{{k_{1}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}\left({1+\psi}\right)-\psi}}}}\right]^{2}}$
(37b)
$\displaystyle{\varepsilon_{r:5,6}}=\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}}}{\left[{\frac{\zeta}{{{k^{0}}}}\cos\gamma\pm\frac{{{k_{-1}}}}{{{k^{0}}}}\sqrt{-{{\left({\frac{\zeta}{{{k_{-1}}}}}\right)}^{2}}{{\sin}^{2}}\gamma+1-\frac{1}{{\frac{{k_{-1}^{2}}}{{{{\left({{k^{0}}}\right)}^{2}}}}\left({1+\psi}\right)-\psi}}}}\right]^{2}}$
(37c)

As evident by Eq. (37a), clearly the effective permittivities that correspond
with the first two modal solutions exhibit a reciprocal behaviour and are
independent of the direction of propagation, and moreover they become positive
only above $\delta\omega/\omega_{0}\approx 0$. This is what we expect from the
a _stationary_ wire medium as given in Eq. (8). This is reasonable since these
two solutions in the spatiotemporally modulated wire media are essentially a
weak _perturbation_ over the solutions of the conventional stationary wire
media. Note however that by comparing the terms in Eq. (8) and Eqs. (37a), it
is clear that the modulation parameters affect to some extent the values of
the permittivity, albeit not its general properties. Moreover, while Eq. (37a)
exhibits a reciprocal, stationary-medium-like, behaviour, there are additional
higher order harmonics (see Eq. (17)) that break the reciprocity. As opposed
to $\varepsilon_{r:1,2}$, the additional effective permittivities, i.e.,
$\varepsilon_{r:3}$ and on, demonstrate a substantial nonreciprocity already
at the fundamental harmonic. Even more important is the fact that these
solutions are propagating already with negative frequencies
$\delta\omega/\omega_{0}<0$ in the region where the first two solutions are
evanescent. This implies, as we already discussed in the previous section, in
the context of Fresnel drag, that a substantial nonreciprocity can be observed
in the spatiotemporally modulated wire media already with weak and slow
modulation parameters. We would like also to stress that these modes are
propagating bellow the plasma frequency _for any modulation depth_ $m$, _even
if extremely small_. Therefore, these modes share similarities with the
extraordinary Whistler mode in magneto-plasma that also propagates below the
plasma frequency with any magnetic biasing [37]. This behaviour is shown in
Fig. 5 that shows for the parameters that are given case II of Table 1 a plot
of the effective permittivity for the first six fundamental modal solutions,
as a function of $\delta\omega/\omega^{0}$ and for two main situations. In
(a,b), real and imaginary parts, where the propagation is parallel to the
modulation axis, $\vec{q}_{0}||\vec{\zeta}$, and in (c,d) where
$\vec{q}_{0}\perp\vec{\zeta}$. The effective permittivity for the stationary
(LTI medium) is shown in black circles, while the solutions for the
spatiotemporally modulated lattice are shown in the colored lines. In Fig. 5
the high anisotropy is evident due to the directional preference dictated by
the spatiotemporal modulation vector $\vec{\zeta}$. In addition, the two
fundamental solutions #1,2 are essentially identical to the two counter-
propagating solutions of the corresponding stationary (LTI) medium (shown in
black circles), and as such they propagate only with $\delta\omega>0$. In
contrast, the higher order solutions may be highly nonreciprocal, except for
the case where the propagation is transverse to the modulation vector (c,d).
As an example, see the counter propagating solutions #3 and #4 in Fig. 5(a)
that are associated with different effective permittivities. Moreover, it is
seen that the spatiotemporal modulation yields a magnetized plasma-like
extraordinary wave propagation in the sense that propagation becomes alowed
with any modulation index below the so-called plasma-frequency $\omega^{0}$.
As opposed to Fig. 5 that uses $\psi=0.2\ll 1$, in Fig. 6 we show similar
results but with larger capacitive loading, so that $\psi=1.07\sim 1$. Here,
the effect of the spatiotemporal modulation becomes stronger, as evident by a
comparison between Fig. 5 and Fig. 6. Also here, the fundamental solutions #1
and #2 exhibit permittivity that is essentially identical to that of the
corresponding unmodulated lattice. However, as opposed to Fig. 5, here, the
extraordinary modes that are supported below the plasma frequency are much
more dispersive, and the non-reciprocity becomes more dominant, as can be
viewed by comparing the counter propagating solution pairs #3,4 and #5,6. We
will show in Sec. VI that as opposed to the effective permittivities that are
shown in Fig. 6 for parameters case III, the effective permittivities shown in
Fig. 5 for parameters case II can be adequately reproduced in a model in which
the spatiotemporal modulation is introduced after the homogenization.

Figure 5: Effective permittivity of spatiotemporally modulated wire medium.
Parameters of case II in Table 1. (a) effective permittivity of the modal
solutions when the propagation is collinear with the modulation. Circles: no
modulation - stationary wire media. Cont. lines: blue-fundamental modal
solutions #1,2 are nearly identical to these of the stationary medium, and are
reciprocal. The other higher order solutions #3,4 and #5,6 exhibit deviation
with respect to the stationary medium solutions and substantial
nonreciprocally. (b) as (a) but for propagation normal to the direction of
modulation. The solutions are reciprocal in this case, and thus the strong
anisotropy due to the space-time modulation clearly emerges by comparing the
two figures (a) and (b). Figure 6: Effective permittivity of the modal
solutions for case III in Table 1. As in Fig. 5 except that the fact that
$\psi\sim 1$ implies much stronger effect by the modulation as evident by
comparing with Fig. 5.

###  V.2 Averaged field equations and Poyinting’s vector

For each one of the the modal solutions, with $q_{0}=q_{0}^{\\#i}$, we write
the governing electrodynamic equations for the _averaged_ fields. We begin by
calculating the average fields in the unit-cell at the origin. To that end, we
use the expression for the currents in Eq. (2). For the fundamental harmonic
of each one of the modal solutions, the electric field can be written in the
same way as in Eq. (4), but with the modification for the entire unit cell,
except on the wire itself,

$\begin{array}[]{l}\vec{E}\left({\vec{r}}\right)=-\frac{{\hat{z}\eta
k}}{4}\sum\limits_{m,l}{H_{0}^{\left(2\right)}\left({k\left|{{{\vec{R}}_{ml}}-\vec{r}}\right|}\right){I_{0}}{e^{-j\left({{q_{x}}am+{q_{y}}bl}\right)}}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,=-\frac{{\hat{z}\eta
k{I_{0}}}}{2}\sum\limits_{m=-\infty}^{\infty}{{S_{m}}\left({\vec{r}}\right){e^{-j{q_{x}}am}}}\end{array}$
(38)

with,

$\begin{array}[]{l}{S_{m}}\left({\vec{r}}\right)=\frac{1}{2}\sum\limits_{l=-\infty}^{\infty}{H_{0}^{\left(2\right)}\left({k\left|{{{\vec{R}}_{ml}}-\vec{r}}\right|}\right){e^{-j{q_{y}}bl}}}\\\
\,\,\,\,\,\,\,\,=-\frac{1}{b}\sum\limits_{l=-\infty}^{\infty}{\frac{{{e^{jy\left({{q_{y}}+\frac{{2\pi
l}}{b}}\right)}}{e^{-ja\left|{m-\frac{x}{a}}\right|{\beta_{x,l}}}}}}{{{\beta_{x,l}}\,}}}.\end{array}$
(39)

The last summation in Eq. (39) was obtained with the aid of Poisson summation.
This summation is rapidity converging and van be expressed in closed form
using a geometrical series. This finally leads to an expression for the
electric field, with separation between the first propagating harmonic and the
remaining Floquet harmonics that are all evanescent in the dense grid
approximation and as long as the modulation frequency is low enough so that
the array can be considered as dense also at the highest temporal harmonic
$\omega_{n}$ that cannot be neglected (see Appendix VIII.4),

$\vec{E}\left({\vec{r}}\right)\,\approx\,\,\,-\frac{{\hat{z}\eta
k{I_{0}}}}{{2b}}\left\\{\begin{array}[]{l}\frac{{\left({1+jy{q_{y}}}\right)}}{{\sqrt{{k^{2}}-q_{y}^{2}}\,}}\left[{1+\frac{{1-j\left({a+x}\right)\sqrt{{k^{2}}-q_{y}^{2}}+j{q_{x}}a}}{{ja\sqrt{{k^{2}}-q_{y}^{2}}-j{q_{x}}a}}+\frac{{1-j\left({a-x}\right)\sqrt{{k^{2}}-q_{y}^{2}}-j{q_{x}}a}}{{ja\sqrt{{k^{2}}-q_{y}^{2}}+j{q_{x}}a}}}\right]\\\
-2\sum\limits_{l=1}^{\infty}{\frac{{{e^{jy\left({\frac{{2\pi
l}}{b}}\right)}}\left({1+j{q_{y}}y}\right)}}{{\left({\frac{{2\pi
l}}{b}+{q_{y}}}\right)\,}}\left[\begin{array}[]{l}{e^{-\left|x\right|\frac{{2\pi
l}}{b}}}\left({1-{q_{y}}\left|x\right|}\right)+\frac{{{e^{-2\pi
l\frac{{a+x}}{b}}}\left({1-{q_{y}}a+j{q_{x}}a-{q_{y}}x}\right)}}{{1-{e^{-2\pi
l\frac{a}{b}}}\left({1-{q_{y}}a+j{q_{x}}a}\right)}}\\\ +\frac{{{e^{-2\pi
l\frac{{a-x}}{b}}}\left({1-{q_{y}}a+j{q_{x}}a+{q_{y}}x}\right)}}{{1-{e^{-2\pi
l\frac{a}{b}}}\left({1-{q_{y}}a-j{q_{x}}a}\right)}}\end{array}\right]}\end{array}\right\\}.$
(40)

In order to simplify the calculation of the average fields inside the unit
cell, since the wire is assumed to be thin $r_{0}\ll a,b$ we will neglect the
fact that the field inside the PEC wire is zero. This simplification enables
us to integrate inside the unit cell without complicated boundaries. This
assumptions is supported by the fact that,

$\begin{array}[]{l}\frac{{2\pi\int\limits_{0}^{{r_{0}}}{{E^{w}}rdr}}}{{ab}}=2\pi\frac{{\eta
k{I_{0}}}}{{4ab}}\int\limits_{0}^{{r_{0}}}{H_{0}^{\left(2\right)}\left({{k_{0}}r}\right)rdr}\\\
\,\,\,\,\,\,\sim
O\left({\frac{{r_{0}^{2}}}{{ab}}\ln{k_{0}}{r_{0}}}\right)\end{array}$ (41)

which implies that the resulting error by this approximation is
$O(r_{0}^{2}/ab)$. The average field is then calculated by,

$\left\langle{\vec{E}}\right\rangle\approx\frac{1}{{ab}}\int\limits_{-\frac{a}{2}}^{\frac{a}{2}}{\int\limits_{-\frac{b}{2}}^{\frac{b}{2}}{\vec{E}\left({\vec{r}}\right)dxdy}}.$
(42)

Thus,

$\left\langle{\vec{E}}\right\rangle\,\,\,\approx\,\,\,-\frac{{\hat{z}\eta{k_{0}}{I_{0}}}}{2}\left\\{\begin{array}[]{l}\frac{2}{{jab\left({k_{0}^{2}-q_{0}^{2}}\right)}}-\frac{1}{{b\sqrt{k_{0}^{2}-q_{y}^{2}}\,}}\\\
-2\sum\limits_{l=1}^{\infty}{\frac{{{q_{y}}b{{\left({-1}\right)}^{l}}}}{{2l\pi
a\left({\frac{{2\pi
l}}{b}+{q_{y}}}\right)\,}}\left[\begin{array}[]{l}\frac{{-{q_{y}}b+2l\pi+{e^{-l\pi\frac{a}{b}}}\left({{q_{y}}al\pi+{q_{y}}b-2l\pi}\right)}}{{2{l^{2}}{\pi^{2}}}}+\\\
\left({\frac{{2{e^{-\left({2\pi
l\frac{a}{b}+{q_{y}}a-j{q_{x}}a}\right)}}\sinh\frac{a}{2}\left({\frac{{2\pi
l}}{b}+{q_{y}}}\right)}}{{b\left({\frac{{2\pi
l}}{b}+{q_{y}}}\right)}}}\right)\left[\begin{array}[]{l}\frac{1}{{1-{e^{-2\pi
l\frac{a}{b}}}\left({1-{q_{y}}a+j{q_{x}}a}\right)}}+\\\ \frac{1}{{1-{e^{-2\pi
l\frac{a}{b}}}\left({1-{q_{y}}a-j{q_{x}}a}\right)}}\end{array}\right]\end{array}\right]}\end{array}\right\\}$
(43)

For a dense array, the infinite summation in Eq. (43) can be neglected, and
thus the average field may be simplified as

$\left\langle{{{\vec{E}}}}\right\rangle\approx-\hat{z}\eta{k_{0}}{I_{0}}\left[{\frac{1}{{jab\left({k_{0}^{2}-q_{0}^{2}}\right)}}-\frac{1}{{2b\sqrt{k_{0}^{2}-q_{y}^{2}}\,}}}\right]$
(44)

The result in Eq. (44) has a good agreement with [27] (see appendix VIII.5).
Following [35], we can write the relation between the macroscopic (averaged)
fields,

$\displaystyle\vec{\nabla}\times{{\vec{E}}_{macro}}=-j\omega{\mu_{0}}{{\vec{H}}_{macro}}$
(45a)
$\displaystyle\vec{\nabla}\times{{\vec{H}}_{macro}}=\frac{I}{{ab}}\hat{z}+j\omega{\varepsilon_{0}}{{\vec{E}}_{macro}}.$
(45b)

Therefore, the macroscopic magnetic field,

${{\vec{H}}_{macro}}=-\frac{1}{{jk\eta}}\vec{\nabla}\times{{\vec{E}}_{macro}}=-\frac{1}{{jk\eta}}j{{\vec{q}}_{0}}\times\left\langle{\vec{E}}\right\rangle.$
(46)

The average Poynting’s vector reads [35],

${S_{avg}}=\frac{1}{2}{\mathop{\rm
Re}\nolimits}\left\\{{{{\vec{E}}_{macro}}\times{{\vec{H}}^{*}}_{macro}+\frac{{{\varphi_{w}}{I^{*}}}}{{ab}}\hat{z}}\right\\}.$
(47)

And since in our case the currents are independent of $z$, $\varphi_{w}=0$,
and we may write,

${S_{avg}}=\frac{1}{2}{\mathop{\rm
Re}\nolimits}\left\\{{{{\vec{E}}_{macro}}\times{{\vec{H}}^{*}}_{macro}}\right\\}=\frac{{{q_{0}}}}{{2k\eta}}{\left|{\left\langle
E\right\rangle}\right|^{2}}.$ (48)

##  VI Perturbed continues media near the plasma region

In the previous sections we have homogenized and derived the effective
permittivity of the various modal solutions that are supported by the
spatiotemporally modulated wire medium. However, there is another way to go.
We could, first homogenize the stationary wire medium. Leading to its famous
plasma-like behaviour, and then, introduce the space-time modulation on the
effective properties of the homogenized material. Evidently, the latter
approach is expected to be less accurate since it does not take into account
the microscopic properties of the modulation, however, it is surely easier to
perform. In this section we compare the two approaches. Specifically, near its
resonance frequency, the effective permittivity of a stationary wire medium
takes the plasma-like form,

${\varepsilon_{r}}=1-\frac{{\omega_{p}^{2}}}{{{\omega^{2}}}}.$ (49)

The spatiotemporal modulation of the structure may now be introduced in the
effective permittivity parameter $\omega_{p}$.

In order to be able to perform a fair comparison it is essential to define
what it the analogue of the capacitance modulation discussed in the previous
sections. Clearly, in the homogenized model, the only quantity we may modulate
is $\omega_{p}$, but then, we need to ask for a given capacitance modulation
$\delta C$ in Eq. (10), what is the equivalent measure of modulation in
$\omega_{p}$?

To that end, we may estimate the behaviour for a small change in $\omega_{p}$
and the frequency. Recall that in Sec. IV we have defined the frequency that
nullifies the effective permittivity of the loaded, but stationary, wire
medium as $\omega^{0}$ (see Eq. (26)). In the continuous model that we
consider here in Eq. (49), $\omega_{p}$ takes the role of $\omega^{0}$. In
order to estimate the small variation behaviour of Eq. (49) we denote
$\omega_{p}=\omega^{0}+\delta\omega_{p}$, and
$\omega=\omega^{0}+\delta\omega$. Then, for frequencies near the plasma
frequency, and for small variations in the plasma frequency, Eq. (49) becomes,

$\begin{array}[]{l}{\varepsilon_{r}}\approx
2\left({\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta{\omega_{p}}}}{{{\omega^{0}}}}}\right).\end{array}$
(50)

By performing a similar procedure for the effective permittivity of the
stationary wire medium that is given in Eq. (8), and comparing the
permittivity near the plasma frequency, we obtain that a variation in
$\omega_{p}$ is related to a variation in $C$ via,

$\frac{{\delta{\omega_{p}}}}{{{\omega^{0}}}}=-\frac{\psi}{2}\,\frac{{\delta
C}}{{{C_{0}}}}$ (51)

when $\psi\ll 1$. If $\psi\gtrapprox 1$ then the simple local plasma model in
Eq. (49) is inadequate to describe the effective model of the loaded wire
medium. Let us now time-modulate the plasma frequency,

${\omega_{p}}\left({\vec{r},t}\right)={\omega^{0}}+\delta{\omega_{p}}\cos\left({\Omega
t-\varphi}\right),$ (52)

with $\varphi={\vec{\zeta}}\cdot{\vec{R}}$ and ${\vec{R}}=x\hat{x}+y\hat{y}$.
Then, we approximate,

$\omega_{p}^{2}\left({\vec{r},t}\right)\approx{\left({{\omega^{0}}}\right)^{2}}\left[{1+2\frac{{\delta{\omega_{p}}}}{{{\omega^{0}}}}\cos\left({\Omega
t-\varphi}\right)}\right].$ (53)

We assume a $\hat{z}$ polarized electric field as in the original wire system.
Then, write Faraday’s law in the frequency domain
$\vec{\nabla}\times\vec{H}\left(\omega\right)=j\omega\hat{z}\varepsilon\left(\omega\right)E\left(\omega\right)$,
and convert to the time domain ($1/j\omega\mapsto$ integration),

$\vec{\nabla}\times\vec{H}\left(t\right)=\hat{z}{\varepsilon_{0}}\frac{{\partial
E\left(t\right)}}{{\partial
t}}+\hat{z}{\varepsilon_{0}}\omega_{p}^{2}\left(t\right)\int\limits_{-\infty}^{t}{E\left(\tau\right)d\tau}.$
(54)

We can now return back to the frequency domain,

$\begin{array}[]{l}\vec{\nabla}\times\vec{H}\left(\omega\right)=j\hat{z}\omega{\varepsilon_{0}}\vec{E}\left(\omega\right)+\\\
\,\,\,\,\hat{z}{\varepsilon_{0}}{\left({{\omega^{0}}}\right)^{2}}\left\\{\begin{array}[]{l}\delta\left(\omega\right)+\\\
\frac{{\delta{\omega_{p}}}}{{{\omega^{0}}}}\left[\begin{array}[]{l}{e^{j\varphi}}\delta\left({\omega-\Omega}\right)+\\\
{e^{-j\varphi}}\delta\left({\omega+\Omega}\right)\end{array}\right]\end{array}\right\\}*\frac{{\vec{E}\left(\omega\right)}}{{j\omega}}\end{array}$
(55)

where $*$ denotes convolution. We will derive the equation for the $n$th
frequency harmonics, and also use the known identity for plane waves
$\vec{\nabla}=-j\vec{q}$

$\begin{array}[]{l}-\frac{1}{{{\varepsilon_{0}}}}j\vec{q}\times{{\vec{H}}_{n}}=\left[{{\omega_{n}}-\frac{{{{\left({{\omega^{0}}}\right)}^{2}}}}{{{\omega_{n}}}}}\right]\hat{z}j{E_{n}}+\\\
\,\,\,\,\,\,\,\,\,\,\,\,\hat{z}\frac{{{\omega^{0}}\delta{\omega_{p}}{e^{j\varphi}}{E_{n-1}}}}{{j{\omega_{n-1}}}}+\,\hat{z}\frac{{{\omega^{0}}\delta{\omega_{p}}{e^{-j\varphi}}{E_{n+1}}}}{{j{\omega_{n+1}}}}\end{array}.$
(56)

Now using in addition Faraday’s equation, we obtain,

${q^{2}}{E_{z}}=-j{\omega_{n}}\mu{\varepsilon_{0}}\frac{{\left({-j\vec{q}}\right)}}{{{\varepsilon_{0}}}}\times\vec{H}.$
(57)

We will now get an equation for the electric field, which can be written as an
infinite tri-diagonal matrix. Let us now focus in the case of linear phase
$\varphi=\vec{\zeta}\cdot\vec{R}$.

$a_{n-1}e^{j\varphi}E_{n-1}-b_{n}E_{n}+a_{n+1}e^{-j\varphi}E_{n+1}=0,$ (58)

with,

$\displaystyle{a_{n}}=\frac{{{k_{0}}}}{{{k_{n}}}}\frac{{\delta{\omega_{p}}}}{{{\omega^{0}}}}$
(59a)
$\displaystyle{b_{n}}=\frac{{{k_{n}}}}{{{k^{0}}}}-\frac{{{k^{0}}}}{{{k_{n}}}}-\frac{{q_{n}^{2}}}{{{k^{0}}{k_{n}}}}$
(59b)

By assuming, as in previous sections, that only the three fundamental
harmonics are dominant, the recursive relation in Eq. (58) may be truncated.
Yielding a $3\times 3$ homogenous system for $E_{n}$, with $n=-1,0,1$. Then,
to get a nontrivial solution we require that its determinant vanishes.
Yielding,

$\left(\frac{\delta\omega_{p}}{\omega^{0}}\right)\left[\frac{a_{-1}}{b_{-1}}+\frac{a_{1}}{b_{1}}\right]=b_{0}$
(60)

in a similar way to which Eq. (28) was derived. Later, we use Eq. (17) and get
the dispersion equation for the various wave solutions supported by the
system. For the sake of brevity we omit the expressions, nevertheless, these
solutions are then used to calculate the effective permittivity of the various
modes, as done in Sec. V. These read,

$\displaystyle{\varepsilon_{r:1}}={\varepsilon_{r:2}}=\frac{{\left({2+\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}}\right)\left({\frac{{\delta\omega}}{{{\omega^{0}}}}-\frac{{\delta\hat{\omega}}}{{{\omega^{0}}}}}\right)}}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}}}$
(61a)
$\displaystyle{\varepsilon_{r:3,4}}=\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}}}{\left[{-\frac{\zeta}{{{k^{0}}}}\cos\gamma\pm\frac{{{k_{1}}}}{{{k^{0}}}}\sqrt{-{{\left({\frac{\zeta}{{{k_{1}}}}}\right)}^{2}}{{\sin}^{2}}\gamma+1-\frac{{{k^{0}}}}{{{k_{1}}}}}}\right]^{2}}$
(61b)
$\displaystyle{\varepsilon_{r:5,6}}=\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}}}{\left[{\frac{\zeta}{{{k^{0}}}}\cos\gamma\pm\frac{{{k_{1}}}}{{{k^{0}}}}\sqrt{-{{\left({\frac{\zeta}{{{k_{1}}}}}\right)}^{2}}{{\sin}^{2}}\gamma+1-\frac{{{k^{0}}}}{{{k_{1}}}}}}\right]^{2}}$
(61c)

In order to compare with the effective primitivity that is obtained by
directly homogenizing the spatiotemporally modulated wire medium we use the
transform
$\frac{{\delta{\omega_{p}}}}{{{\omega^{0}}}}=-\frac{\psi}{2}\,\frac{{\delta
C}}{{{C_{0}}}}$ for the cases listed in Table 1, and thus get a new set of
corresponding four cases that are listed in Table 2.

| $m=\delta\omega_{p}/\omega^{0}$ | $\Omega/\omega^{0}$ | $\zeta/k^{0}$ | $\theta_{0}$ | $\xi$
---|---|---|---|---|---
Case I | 0.014 | 0.15 | 0.1 | $\pi/3$ | $0$
Case II | 0.035 | 0.2 | 1.0 | $0$ | $0$
Case III | 0.03 | 0.25 | 1.2 | $\pi/5$ | $\pi/4$
Case IV | 0.18 | 0.3 | 1.0 | $\pi$ | $0$
Table 2: Four parameter cases used in the numerical examples for the effective
permittivity in a continuous plasma model with spatiotemporally modulated
plasma-frequency. The values here correspond to these in Table 1 for the wire
medium.

where here the modulation depth $m=\delta\omega_{0}/\omega^{0}$. Specifically,
we calculate the effective permittivities for the modal solutions that
correspond to cases II and III. These are shown in Fig. 7 and Fig. 8, and
should be compared with the wire-medium counterparts in Fig. 5 and Fig. 6
above.

Figure 7: Effective permittivity of the modal solutions in a spatiotemporally
modulated continuous plasma-like model. Results for case II in Table 2. This
figure should be compared with Fig. 5 that was calculated by brute force
homogenization of the spatiotemporally modulated wire medium. (a) effective
permittivity when the propagation is collinear with the modulation. Circles:
no modulation - stationary plasma medium. Cont. lines: blue-fundamental modal
solutions #1,2 are nearly identical to these of the stationary medium, and are
reciprocal. The other higher order solutions #3,4 and #5,6 exhibit deviation
with respect to the stationary medium solutions and nonreciprocally. (b) as
(a) but for propagation normal to the direction of modulation. The solutions
are reciprocal in this case, and thus the strong anisotropy due to the space-
time modulation clearly emerges by comparing the two figures (a) and (b). A
comparison with Fig. 5 demonstrates strong similarities between the results of
the two homogenization approaches. Figure 8: Effective permittivity,
spatiotemporally modulated continuous plasma-like model. As in Fig. 7 but
calculated for case III in Table 2. These results should be compared with Fig.
6, exhibiting large differences between the two homogenization approaches in
the case $\psi\sim 1$.

Specifically, it is seen that the two homogenization approaches are nearly
equivalent, with minor differences for cases with $\psi\ll 1$, as seen for
example by comparing Fig. 7 and Fig. 5. This behaviour changes for cases with
larger capacitive loading $\psi\gtrapprox 1$ as seen for example by comparing
Fig. 8 and Fig. 6. In this case since the capacitive loading is not
negligible, the plasma model in Eq. (49) becomes less accurate already for the
stationary medium case [36]. Nevertheless, while the effect of the
spatiotemporal modulation on the fundamental solutions #1,2 is moderate, the
effect on the higher order solutions, that dominate below $\omega^{0}$, is
substantial, despite the fact that the modulation index in this case $m=0.03$
is small.

##  VII Conclusion

In this paper we have developed homogenization theory for spatiotemporally
modulated wire medium. Our analysis takes into account the complete
interaction between the time-modulated wires, both in space and time. The
dispersion relations that we derive for the modes that are supported by the
lattice demonstrate peculiar wave phenomena such as guidance below the cutoff
frequency (so called plasma-frequency) of the stationary medium, with
modulation depth as weak as desired. This behaviour shares similarities with
the extraordinary mode that is guided in magnetized plasma, parallel to the
magnetization direction, the so called Whistler mode. Furthermore, in light of
the ability to guide below the plasma frequency, we show that this system
provides a mean to achieve substantial effective Fresnel drag with weak and
slow modulation. In addition to these wave phenomena we calculate,
analytically, the effective permittivity of the first low order solutions, and
moreover derive expressions for the averaged (macroscopic) fields, and
Poynting’s vector. Lastly, we compare between two homogenization approaches of
spatiotemporally modulated wire medium, first, that includes the modulation in
the homogenization, and second that introduces the modulation into the
effective parameters of the homogenized corresponding stationary (LTI) medium.
We show that under certain conditions these approaches provide similar
results, however, they may substantially deviate in other parameters regimes.

## Acknowledgment

This research was supported by the Israel Science Foundation (grant No.
1353/19).

##  VIII Appendixes

###  VIII.1 Evaluation of the dispersion of LTI wire medium

This calculation can be found entirely in textbooks, such as in [36],
nevertheless, for the sake of completeness and self-continency, we briefly
provide it below. The local field at the reference point $m=l=0$,

$\begin{array}[]{l}E_{z}^{loc}=-\frac{{\eta k{I_{0}}}}{4}\sum\limits_{m,l\neq
0,0}{H_{0}^{\left(2\right)}\left({k{R_{ml}}}\right){e^{-j\left({{q_{x}}am+{q_{y}}bl}\right)}}}\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,=-\frac{{\eta
k{I_{0}}}}{2}\left({{S_{0}}+\sum\limits_{m\neq
0}{{S_{m}}{e^{-j{q_{x}}am}}}}\right).\end{array}$ (62)

Where $S_{m}$ denotes the summation for fixed $m$, over the indices $l$. For
$m=0$ we have,

$\begin{array}[]{l}{S_{0}}=\sum\limits_{l=1}^{\infty}{H_{0}^{\left(2\right)}\left({k{R_{0l}}}\right)\cos{q_{y}}bl}\\\
\,\,\,\,\,\,=\frac{1}{{b\sqrt{{k^{2}}-q_{y}^{2}}}}-\frac{1}{2}+\\\
\,\,\,\,\,\,\,\,\,\,\,\,\,\,\frac{j}{\pi}\left\\{{\ln\frac{{bk}}{{4\pi}}+\gamma+\frac{1}{2}\sum\limits_{l\neq
0}{\left[{\frac{{-2\pi
j}}{{b{\beta_{x,l}}}}-\frac{1}{{\left|l\right|}}}\right]}}\right\\}.\end{array}$
(63)

Here $\gamma=0.5772$ is the Euler constant. For $m\neq 0$ we use the Poisson
summation which becomes highly efficient for dense grid, and get,

$\begin{array}[]{l}{S_{m}}=\frac{1}{2}\sum\limits_{l=-\infty}^{\infty}{H_{0}^{\left(2\right)}\left({k{R_{ml}}}\right){e^{-j{q_{y}}bl}}}\\\
\,\,\,\,\,\,\,\,=-\frac{1}{b}\sum\limits_{l=-\infty}^{\infty}{\frac{{{e^{-ja\left|m\right|{\beta_{x,l}}}}}}{{{\beta_{x,l}}\,}}}.\end{array}$
(64)

By substituting Eq. (5b) to Eq. (4), and changing the order of summations, and
using,

$\sum\limits_{m\neq
0}{{e^{-j\left({{\beta_{x,l}}a\left|m\right|+{q_{x}}am}\right)}}}=\frac{{j\sin{\beta_{x,l}}a}}{{\cos{\beta_{x,l}}a-\cos{q_{x}}a}}-1.$
(65)

We can rewrite Eq. (3),

$\begin{array}[]{l}\frac{{\eta
k}}{4}H_{0}^{\left(2\right)}\left({k{r_{0}}}\right)+\frac{1}{{j\omega\tilde{C}_{0}}}\,=\\\
\,\,\,\,\,\,-\frac{{\eta
k}}{2}\left\\{{{S_{0}}+\sum\limits_{l=-\infty}^{\infty}{\frac{1}{{{\beta_{x,l}}}}\left[{\frac{{j\sin{\beta_{x,l}}a}}{{\cos{\beta_{x,l}}a-\cos{q_{x}}a}}-1}\right]}}\right\\}\end{array}$
(66)

Note that we assumed small losses in he background medium, so that convergence
is secured. We will assume that the wire radius is very small comparing to the
separation between the wires, and use the approximation for the Hankel
function [46],

$H_{0}^{\left(2\right)}\left({k{r_{0}}}\right)\approx
1-j\frac{2}{\pi}\left({\ln\frac{{k{r_{0}}}}{2}+\gamma}\right).$ (67)

This leads to the dispersion relation,

$\begin{array}[]{l}\frac{j}{\pi}\ln\frac{{k{r_{0}}}}{2}-\frac{2}{{jk\eta\omega\tilde{C}_{0}}}\,=\frac{1}{{b\sqrt{{k^{2}}-q_{y}^{2}}}}+\\\
\,\,\,\,\,\,\frac{j}{\pi}\left[{\ln\frac{{bk}}{{4\pi}}+\frac{1}{2}\sum\limits_{l\neq
0}{\left({\frac{{-2\pi
j}}{{b{\beta_{x,l}}}}-\frac{1}{{\left|l\right|}}}\right)}}\right]+\\\
\,\,\,\,\,\sum\limits_{l=-\infty}^{\infty}{\frac{1}{{{\beta_{x,l}}}}\left[{\frac{{j\sin{\beta_{x,l}}a}}{{\cos{\beta_{x,l}}a-\cos{q_{x}}a}}-1}\right]}.\end{array}$
(68)

And for a dense grid ($ka,kb\ll 1$), with the approximations for trigonometric
functions with small arguments we can get,

$\begin{array}[]{l}\ln\frac{b}{{2\pi{r_{0}}}}+\sum\limits_{l=1}^{\infty}{\frac{{\coth\frac{{\pi
al}}{b}-1}}{l}}-\frac{{2\pi}}{{\eta v{k^{2}}\tilde{C}_{0}}}\\\
\,\,\,\,\,\,\,\,=\frac{{2\pi}}{{ab\left[{{k^{2}}-\left({q_{x}^{2}+q_{y}^{2}}\right)}\right]}}-\frac{{\pi
a}}{{6b}}.\end{array}$ (69)

###  VIII.2 Derivation of Eq. (11)

We begin by transforming the expression for the frequency domain
susceptibility of a loaded wire, i.e., Eq. (1), into the time-domain. It easy
to notice that $\alpha^{-1}$ behaves as an impedance to unit length, and
$E^{loc}$ behaves as an external “voltage” source (obviously, with dimensions
[V/m]). Therefore, we can consider our physical model to be equivalent to the
simple serial circuit that is shown in Fig. 9 with stationary impedance
$\alpha_{0}^{-1}$ per unit length (involves resistance as well as reactance)
and time modulated capacitance.

Figure 9: The equivalent circuit of a wire illuminated by an external field

The current on the capacitor, reads

$I=\frac{d}{{dt}}\left({\tilde{C}\left(t\right){E_{C}}\left(t\right)}\right)$
(70)

where $E_{C}$ denotes the “voltage” on the capacitance. Then, Eq. (11) is
straightforward.

###  VIII.3 Derivation of Eq. (17)

The wire media behaves as an infinite periodic crystal, and therefore we can
use Floquet-Bloch theorem, with the expression in Eq. (16). By substituting
Eq. (16) into Eq. (15), and using the explicit expression for the phase of the
capacitors
${\varphi_{m,l}}=\vec{\zeta}\cdot\vec{r}=ma\zeta\cos\xi+bl\zeta\sin\xi$. We
get for the $m,l$ wire, and for n-th mode:

$\begin{array}[]{l}E_{n}^{inc}{e^{-j{{\vec{q}}_{n}}\cdot\vec{R}_{ml}}}=\left[{{\gamma_{0}}{{\omega_{n}}}+\frac{1}{{j{{\omega_{n}}}{\tilde{C}_{0}}}}}\right]{A_{n}}{e^{-j{{\vec{q}}_{n}}\cdot\vec{R}_{ml}}}-\\\
A_{n-1}\frac{m}{2}\frac{1}{j\omega_{n-1}\tilde{C}_{0}}e^{-j(\vec{\zeta}+\vec{q}_{n-1})\cdot\vec{R}_{ml}}-\\\
A_{n+1}\frac{m}{2}\frac{1}{j\omega_{n+1}\tilde{C}_{0}}e^{-j(-\vec{\zeta}+\vec{q}_{n+1})\cdot\vec{R}_{ml}}.\end{array}$
(71)

The expression in (71) must be valid for all the wires $m,l$. Therefore, we
must require that

$-\vec{\zeta}-{{\vec{q}}_{n-1}}=-{{\vec{q}}_{n}},\mbox{ and
}\vec{\zeta}-{{\vec{q}}_{n+1}}=-{{\vec{q}}_{n}}$ (72)

Therefore, we get ${{\vec{q}}_{n}}={{\vec{q}}_{0}}+n\vec{\zeta}$ as in Eq.
(17).

###  VIII.4 Calculating the field in unit cell

In order to evaluate the electric field in the unit cell at the origin, we use
Eq. (38),

$\begin{array}[]{l}\vec{E}\left({\vec{r}}\right)=\\\
\,\,\,\,-\frac{{\hat{z}\eta
k{I_{0}}}}{{2b}}\sum\limits_{m=-\infty}^{\infty}{\sum\limits_{l=-\infty}^{\infty}{\frac{{{e^{jy\left({{q_{y}}+\frac{{2\pi
l}}{b}}\right)}}{e^{-ja\left|{m-\frac{x}{a}}\right|{\beta_{x,l}}}}{e^{-j{q_{x}}am}}}}{{{\beta_{x,l}}\,}}}}\end{array}$
(73)

Then we approximate ${\beta_{x,l}}=-j\sqrt{{{\left({{q_{y}}+\frac{{2\pi
l}}{b}}\right)}^{2}}-{k^{2}}}\approx-j\left({\frac{{2\pi
l}}{b}+{q_{y}}}\right)$. This leads to,

$\begin{array}[]{l}\vec{E}\left({\vec{r}}\right)=-\frac{{\hat{z}\eta
k{I_{0}}}}{{2b}}\cdot\\\
\,\,\,\,\,\,\sum\limits_{l=-\infty}^{\infty}{\frac{{{e^{jy\left({{q_{y}}+\frac{{2\pi
l}}{b}}\right)}}}}{{{\beta_{x,l}}}}\left[\begin{array}[]{l}\frac{{{e^{-j\left({a+x}\right){\beta_{x,l}}}}{e^{j{q_{x}}a}}}}{{1-{e^{-ja{\beta_{x,l}}}}{e^{j{q_{x}}a}}}}+{e^{-j\left|x\right|{\beta_{x,l}}}}\\\
+\frac{{{e^{-j\left({a-x}\right){\beta_{x,l}}}}{e^{-j{q_{x}}a}}}}{{1-{e^{-ja{\beta_{x,l}}}}{e^{-j{q_{x}}a}}}}\end{array}\right]}\end{array}$
(74)

For a dense grid ${\beta_{x,l}}$ may be further approximated,

${\beta_{x,l}}\approx\left\\{{\begin{array}[]{*{20}{c}}{\sqrt{{k^{2}}-q_{y}^{2}}}&{l=0}\\\
{-j\left({\frac{{2\pi\left|l\right|}}{b}+{q_{y}}}\right)}&{l\neq
0}\end{array}}\right.$ (75)

We can see that we get ${\beta_{x,l}}\approx{\beta_{x,-l}}$. These eventually
lead to Eq. (40).

###  VIII.5 Heavily dense grid

We compare Eq. (44) with the average over the $z$-component of the electric
field as calculated in [27],

$\left\langle{E_{z}}\right\rangle=\left({j\omega\tilde{L}+{Z_{w}}}\right)I+\frac{{\partial\varphi}}{{\partial
z}}$ (76)

Where ${\tilde{L}}$ is the effective inductance per unit length of the wire,
and ${Z_{w}}$ is the impedance loading on the wire. In our case, for a PEC
wire with capacitance load, we have
${Z_{w}}=\frac{1}{{j\omega\tilde{C}_{0}}}$. Also, since we are focusing on the
case of a propagating wave in a direction normal to the wires, there is no
variation along the $z$-axis and thus we can nullify the potential derivative
$\frac{{\partial\varphi}}{{\partial z}}=0$. Then we get (76),

$\begin{array}[]{l}\left\langle{{E_{z}}}\right\rangle=jI\left[{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right){\omega^{0}}\tilde{L}-\frac{1}{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right){\omega^{0}}\tilde{C}_{0}}}}\right]\\\
\,\,\,\,\,\,\,\,\,\approx
jI\frac{\eta}{{ab{k^{0}}}}\left[{1+\left({1+2\psi}\right)\left({\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}\right]\end{array}$
(77)

The expression in (77) is valid for a very dense grid in a stationary medium.
Now to show that our averaged field in Eq. (44) is consistent with that. Note
that the second term in Eq. (44) should be neglected for a heavily dense grid.
Then, using the Eq. (31) for the wave number in the media (with zero
perturbation $\delta C=0$), we get,

$\begin{array}[]{l}\left\langle{\vec{E}}\right\rangle\,\approx-\frac{{\hat{z}\eta{k_{0}}{I_{0}}}}{{ab}}\left\\{{\frac{1}{{j\left({k_{0}^{2}-q_{0}^{2}}\right)}}}\right\\}\\\
\,\approx-\frac{{\hat{z}\eta{k_{0}}{I_{0}}}}{{ab}}\left\\{{\frac{1}{{j\left({k_{0}^{2}-{{\left({{k^{0}}}\right)}^{2}}{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}\left[{1-\frac{1}{{{{\left({1+\frac{{\delta\omega}}{{{\omega^{0}}}}}\right)}^{2}}\left({1+\psi}\right)-\psi}}}\right]}\right)}}}\right\\}\\\
\,\approx
j\frac{{\hat{z}\eta{I_{0}}}}{{ab{k^{0}}}}\left\\{{\frac{{1+\frac{{\delta\omega}}{{{\omega^{0}}}}}}{{1+2\frac{{\delta\omega}}{{{\omega^{0}}}}-2\left({1+\psi}\right)\frac{{\delta\omega}}{{{\omega^{0}}}}}}}\right\\}\\\
\,\approx
j\frac{{\hat{z}\eta{I_{0}}}}{{ab{k^{0}}}}\left\\{{1+\left({1+2\psi}\right)\frac{{\delta\omega}}{{{\omega^{0}}}}}\right\\}.\end{array}$
(78)

Identical to Eq. (77).

## References

  * [1] R. L. Fante, “Transmission of electromagnetic waves into time-varying media,” _IEEE Trans. Ant. Prop._ , 193 pp. 417-424. (1971).
  * [2] F. R. Morgenthaler “Velocity modulation of electromagnetic waves,” _IRE Trans. Microw. Theory Tech._ 6 pp. 167-172 (1958).
  * [3] A. Akbarzadeh, N. Chamanara, and C. Caloz, “Inverse prism based on temporal discontinuity and spatial dispersion,” _Opt. lett._ 43 (14) pp. 3297-3300 (2018).
  * [4] S. F. Preble, Q. Xu, and M. Lipson, “Changing the colour of light in a silicon resonator,” _Nat. Photon._ 1, pp. 293–296 (2007).
  * [5] Daniel Torrent, Olivier Poncelet and Jean-Chirstophe Batsale, “Nonreciprocal Thermal Material by Spatiotemporal Modulation,” _Phys. Rev. Lett_. 120, 125501 (2018).
  * [6] H. Nassar, X. C. Xu, A. N. Norris and G. L. Huang. “Modulated phononic crystals: Non-reciprocal wave propagation and Willis materials,” _J. Mech. Phys. Sol._ 101, 2017, pp. 10-29 (2017).
  * [7] N. A. Estep, D. L. Sounas, J. C. Soric, and A. Alú, “Magnetic-free non-reciprocity and isolation based on parametrically modulated coupled-resonator loops,” _Nat. Phys._ 10 (12) pp. 923-927 (2014).
  * [8] J. R. Zurita-Sánchez, P. Halevi, J. C. Cervantes-González, “Reflection and transmission of a wave incident on a slab with a time-periodic dielectric function,” _Phys. Rev. A_ , 79, 053821 (2009).
  * [9] J. S. Martínez-Romero, O. M. Becerra-Fuentes, and P. Halevi, “Temporal photonic crystals with modulations of both permittivity and permeability,”. _Phys. Rev. A_ , 93, 063813 (2016).
  * [10] V. Bacot, M. Labousse, A. Eddi, M. Fink, and E. Fort, “Time reversal and holography with spacetime transformations,” _Nat. Phys._ 12, pp. 972–977 (2016).
  * [11] S. Vezzoli, V. Bruno, C. DeVault, T. Roger, V. M. Shalaev, A. Boltasseva, M. Ferrera, M. Clerici, A. Dubietis, and D. Faccio, “Optical time reversal from time-dependent Epsilon-Near-Zero media,” _Phys. Rev. Lett._ 120, 043902 (2018).
  * [12] M. A. Salem, C. Caloz, “Wave propagation in periodic temporal slabs,” _9th European Conference on Antennas and Propagation (EuCAP)_ pp. 2164-3342 (2015).
  * [13] M. A. Salem, and C. Caloz, “Temporal photonic crystals: Causality versus periodicity,” _International Conference on Electromagnetics in Advanced Applications (ICEAA)_ (2015).
  * [14] G. Lerosey, J. De Rosny, A. Tourin, A. Derode, G. Montaldo, and M. Fink, “Time reversal of electromagnetic waves,” _Phys. Rev. Lett._ 92, 193904 (2004).
  * [15] V. Pacheco-Peña, “Effective medium concept in temporal metamaterials,” _Nanophotonics_ , 9, (2), pp. 379–391 (2020).
  * [16] D. M. Solis and N. Engheta, “A Generalization of the Kramers-Kronig Relations for Linear Time-Varying Media,” _Phys. Rev. B_ 103, 144303 (2021).
  * [17] Y. Hadad, D. L. Sounas, and A. Alú, “Space-time gradient metasurfaces,” _Phys. Rev. B_ 92, 100304(R) (2015).
  * [18] Y. Hadad, J. C. Soric, A. Alú, “Breaking temporal symmetries for emission and absorption,” _Proc. Nat. Acad. Sci._ 113 (13), pp. 3471-3475 (2016).
  * [19] Y. Hadad, D. L. Sounas, “Space-Time Modulated Loaded-Wire Metagratings for Magnetless Nonreciprocity and Near-Complete Frequency Conversion,” _arXiv_ :1906.00215 (2019).
  * [20] R. E. Collin, _Foundations for Microwave Engineering_ , McGraw-Hill (1966).
  * [21] S. A. Tretyakov, “Electromagnetic Metamaterials: Past, Present and Future,” (2016).
  * [22] A. Shlivinski and Y. Hadad, “Beyond the Bode-Fano Bound: Wideband Impedance Matching for Short Pulses Using Temporal Switching of Transmission-Line Parameters,” _Phys. Rev. Lett._ 121, 204301 (2018).
  * [23] H. Li, and A. Alú, “Temporal Switching to Extend the Bandwidth of Thin Absorbers,” _Optica_ 8, (1), pp. 24-29 (2021)
  * [24] J. C. M. Garnett, “Colours in metal glasses and in metallic films,” _Phil. Trans. Royal Soc. of London Series A_ , 203, pp. 385-420 (1904).
  * [25] W. Rotman, “Plasma imulation by artificial dielectrics and parallel plate media,” _IRE Trans. Ant. Prop._ 10 (1962).
  * [26] K. E. Golden, “A Study of Artificial Dielectrics,” Technical Report No. TR-269 (4280-10)-4 (1964).
  * [27] S. I. Maslovski, and M. G. Silveirinha, “Nonlocal permittivity from a quasistatic model for a class of wire media,” _Phys. Rev. B_ 80, 245101 (2009).
  * [28] P. A. Belov, S. A. Tretyakov and A. J. Viitanen, “Dispersion and Reflection Properties of Artificial Media Formed By Regular Lattices of Ideally Conducting Wires,” _J. Electromag. Waves and App._ 16 pp. 1153-1170 (2012).
  * [29] N. A. Nicorovichi, R. C. McPhedran and L. C. Botten, “Photonic bad gaps for array of perfectly conducting cylinders,” _Phys. Rev. E_ 52, 1135 (1995).
  * [30] S. K. Chin, N. A. Nicorovici and R. C. McPhedran, “Green’s function and lattice sums for electromagnetic scattering by a square array of cylinders,” _Phys. Rev. E_ 49, 4590 (1994).
  * [31] V. Kuzmiak, A. A. Maradudin, and F. Pincemin, “Photonic band structures of two dimensional systems containing metallic components,” _Phys. Rev. B_ 50 (23) pp. 16835-16844 (1994).
  * [32] V. Kuzmiak, A. A. Maradudin and A. R. Mcgurn, “Photonic band structures of two dimentional systems fabricated from rods of a cubic polar crystal,” _Phys. Rev. B_ 55 (7), pp. 4298 (1997).
  * [33] K. Sakoda, “Photonic bands of metallic systems I. Principle of calculation and accuracy,” _Phys. Rev. B_ 64, 045116 (2001).
  * [34] M. M. Sigalas, “Metallic Photonic band gap materials,” _Phys. Rev. B_ 52, 11744 (1995).
  * [35] M. G. Silveirinha and S. I. Maslovski, “Radiation from elementary sources in a uniaxial wire medium,” _Phys. Rev. B_ 85, 155125 (2012).
  * [36] S. Tretyakov, _Analytical Modeling in Applied Electromagnetics_ , Artech House, Boston (2003).
  * [37] A. Ishimaru, _Electromagnetic wave propagation, radiation, and scattering_ , Prentice-Hall, New-Jersy (1991).
  * [38] P. A. Huidobroa, E. Galiffib, S. Guenneau, R. V. Crasterc, and J. B. Pendryb, “Fresnel drag in space–time-modulated metamaterials,” _Proc. Nat. Acad. Sci._ 116 (50) pp. 24943-24948 (2019).
  * [39] D. M. Solís, R. Kastner, N. Engheta, “Time-Varying Materials in Presence of Dispersion: Plane-Wave Propagation in a Lorentzian Medium with Temporal Discontinuity,” _arXiv:_ 2103.06142, (2021)
  * [40] D. M. Solís and N. Engheta, “Functional analysis of the polarization response in linear time-varying media: A generalization of the Kramers-Kronig relations,” _Phys. Rev. B_ 103, 144303 (2021).
  * [41] R. E. Collin, _Field theory of guided waves_ , Wiely (1991).
  * [42] M. S. Mirmoosa, T. T. Koutserimpas, G. A. Ptitcyn, S. A. Tretyakov, and R. Fleury, “Dipole polarizability of time-varying particles,” _arXiv:2002.12297v2_ (2020).
  * [43] Y. Mazor and A. Alú, “Nonreciprocal hyperbolic propagation over moving metasurfaces,” _Phys. Rev. B_ 99, 045407 (2019).
  * [44] Y. Mazor and A. Alú, “One Way Hyperbolic Metasurfaces Based on Synthetic Motion,” _IEEE Trans. Ant. Prop._ 68, (3) pp. 1739 - 1747 (2020).
  * [45] T. Boyd, J. Gratus, P. Kinsler,R. Letizia, and R. Seviour, “Mode Profile Shaping in Wire Media: Towards An Experimental Verification,” _Appl. Sci._ 8 (8), p. 1276, (2018).
  * [46] F. W. J. Olver, and L. C. Maximon, _Bessel functions_ in _NIST Handbook of Mathematical Functions_ , Cambridge University Press (2010).