Properties of translation operator and the solution of the eigenvalue and boundary value problems of arbitrary space–time periodic metamaterials

There is a recent interest in understanding and exploiting the intriguing properties of space–time metamaterials. In the current manuscript, the time periodic circuit theory is exploited to introduce an appropriate translation operator that fully describes arbitrary space–time metamaterials. It is shown that the underlying mathematical machinery is identical to the one used in the analysis of linear time invariant periodic structures, where time and space eigen-decompositions are successively employed. We prove some useful properties the translation operator exhibits. The wave propagation inside the space time periodic metamaterial and the terminal characteristics can be rigorously determined via the expansion in the operators eigenvectors (space–time Bloch waves). Two examples are provided that demonstrate how to apply the framework. In the first, a space time modulated composite right left handed transmission line is studied and results are verified via time domain computations. Furthermore, we apply the theory to explain the non-reciprocal behaviour observed on a nonlinear transmission line manufactured in our lab. Bloch-waves are computed from the extracted circuit parameters. Results predicted using the developed machinery agree with both measurements and time domain analysis. Although the analysis was carried out for electric circuits, the approach is valid for different domains such as acoustic and elastic media.


SYE, 0000-0002-6527-9861
There is a recent interest in understanding and exploiting the intriguing properties of space-time metamaterials. In the current manuscript, the time periodic circuit theory is exploited to introduce an appropriate translation operator that fully describes arbitrary space-time metamaterials. It is shown that the underlying mathematical machinery is identical to the one used in the analysis of linear time invariant periodic structures, where time and space eigen-decompositions are successively employed. We prove some useful properties the translation operator exhibits. The wave propagation inside the space time periodic metamaterial and the terminal characteristics can be rigorously determined via the expansion in the operators eigenvectors (space-time Bloch waves). Two examples are provided that demonstrate how to apply the framework. In the first, a space time modulated composite right left handed transmission line is studied and results are verified via time domain computations. Furthermore, we apply the theory to explain the non-reciprocal behaviour observed on a nonlinear transmission line manufactured in our lab. Bloch-waves are computed from the extracted circuit parameters. Results predicted using the developed machinery agree with both measurements and time domain analysis. Although the analysis was carried out for electric circuits, the approach is valid for different domains such as acoustic and elastic media.

Introduction
Recently, there has been a surge of interest in studying the properties of space-time modulated systems that arise from the intrinsic asymmetric interaction of the space-time harmonics [1,2]. Such asymmetry breaks the principle of reciprocity, hence enabling the design of novel non-reciprocal devices such as nonreciprocal antenna [3][4][5], magnet-less circulators [6][7][8], one-way beam splitters [9], isolators [10,11] and space-time modulated metasurfaces [12,13]. On top of that, systems that possess time and/or space-time periodic elements may not be constrained by the physical limitations of linear time invariant (LTI) systems. For instance, a time modulated reactance does not necessarily result in a total reflection of an incident wave impinging the structure, hence enabling the accumulation of energy [14]. Additionally, a time switched transmission line was demonstrated to have a broadband matching capability not limited by the Bode-Fano criteria that impose a return loss/bandwidth trade-off [15,16]. Taravati & Caloz [17] provide an excellent review on the developments, applications and methods of analysis of space-time media.
The interest in space-time modulated media dates back to the middle of the last century in the context of understanding the properties of distributed parametric amplifiers [18][19][20][21][22][23]. The salient properties of such media emerge from Bloch-Floquet theory, which states that the eigensolutions are the sum of infinite space time harmonics. Unlike LTI systems, space-time modulation leads to an asymmetric dispersion relation, where at a given frequency ω, the forward and backward wavenumbers are not necessarily equal [1,2].
In general, space-time modulated systems are built from nonlinear lumped elements. The interaction of the nonlinearity and a strong pump results in a spatio-temporal modulation of one or more system parameters [24,25]. In a conventional treatment, it is generally assumed that the system is distributed such that the modulation wavelength is much longer than the length of the unit cell. Such an assumption permits describing the system via partial differential equations (PDEs). Upon the application of the Bloch-Floquet condition, the PDEs reduce to an infinite system of homogeneous equations, where the eigenmodes can be calculated from its non-trivial solutions. Very recently, we have developed a circuit based framework that extends the theory of linear time periodic (LTP) circuits and systems, developed in [26][27][28], to space-time structures [29]. The dispersion relation emerges from the application of the Bloch condition to a unit cell that links its input and output time harmonics. Therefore, it is valid for both electrically long and electrically short systems. Furthermore, the framework enables the exploration of various structures such as Composite Right Left Handed (CRLH) transmission lines (TLs) and non-sinusoidal periodic modulation [29]. The governing equations reduce to generalized telegraphist's equations when the unit cell is infinitesimally small.
In the current manuscript, we exploit the circuit approach [29] to develop a translation operator and explore its basic properties. We show that the approach is, in principle, equivalent to how periodic structures are described via the diagonalization of the translation operator, where the eigenmodes are nothing but the Bloch waves. Additionally, for a generic unit cell and an arbitrary periodic modulation, the boundary value problem is solved via the expansion of the solution inside the structure in terms of the eigenmodes. Section 2 starts with a brief review of how the immittance matrix, a generalization of the immittance circuit parameter in LTI systems, emerges from Bloch-Floquet theorem. We then proceed by showing how elements in cascade combine and how the ABCD parameters change between unit cells. In §3, we focus on the eigenvalue problem that describes the system modal behaviour. The invariance of eigenvalues and eigenvectors resulting from the transformation of the system translation operator is presented. We also show that for a generic space-time circuit, the eigensolutions along the modulation line are equivalent. Furthermore, the section discusses how the LTI and LTP systems are formally equivalent in the sense of eigen-decompositions in time and spatial domains. Section 4 demonstrates how the driven modal solution is expanded in terms of the system eigenmodes. Additionally, expressions of the transmission and reflection coefficients are derived. In §5, two systems are studied. In §5.1, a CRLH TL is fully described using the developed machinery and results are compared to time domain simulation. Dispersion relations, eigenvalues, eigenvectors, waveforms and transmission coefficient are computed for both the right-hand (RH) and left-hand (LH) regimes. Subsection 5.2 applies the framework to a nonlinear RH TL that has been fabricated in our lab. The modulation is achieved via a strong pump and hence, the TL operates in the sonic regime [22]. Dispersion relations, eigenvalues and waveforms are computed and compared to simulation. Furthermore, the scattering parameters are calculated and compared to measurements.

Translational properties of immittance and ABCD matrices
In LTP circuits, the Bloch Floquet theorem allows the voltage and current harmonics to be related via immittance matrices [26,29]. An immittance matrix is best visualized as a generalization of the concept of immittance, a scalar complex quantity, in LTI systems. Before proceeding with the detailed description, it is worth noting that we represent the (m, n) element in matrix A, using the notation A n m , i.e. the subscript (superscript) represents the row (column). Without loss of generality, consider a time modulated capacitanceCðtÞ. The eigensolutions of LTP systems are in the form of p(t)exp (iωt), p(t + T ) = p(t) [29,30], hence the voltage-current relation can be described compactly by the matrix equation ð2:1Þ The entries of an arbitrary kth row can be determined from the zeroth row, wherẽ Consider now a structure where the above capacitance is modulated via a travelling wave with speed n, i.e.C ðt, xÞ ¼C t À x n , and x is a multiple of the underlying spatial lattice distance p (i.e. x = np, n = −∞, …, − 2, − 1, 0, 1, 2, …, ∞). Again, expandingC in its Fourier components, and noting that the modulation frequencyv and wavenumberb are related byv ¼nb, imply that This means that the elements in a given row of an admittance matrixỸðxÞ are those in the same row ofỸð0Þ, but multiplied by a phasor that rotates in the counter clockwise direction as we go from left to right. Additionally for a fixed column, the elements from top to bottom are multiplied by a clockwise rotating phasor.
In the subsequent development, the ABCD parameters of the unit cell will be shown to play a significant role. They are formed via the multiplication of LTP impedance and admittance matrices. Therefore, it is crucial to understand the properties of the product of matrices representing space-time modulated elements. Consider for instance the series Z and shunt Y of a lumped right-handed transmission line. It can be shown that the (r, r − s) element of ðZYÞ rÀs r at position x ðZYÞ rÀs r ðxÞ ¼ ðZYÞ rÀs r ð0Þ e Àisbx : Therefore, we have the following important property: Property 2.1. For a space-time periodic structure consisting of a cascade of space-time periodic unit cells, the ABCD parameters X = A, B, C and D for a unit cell x away from the origin are related to the ones at the origin by where G qÀp W exp ðÀi½q À pbpÞ. Hence, the ABCD parameters follow the same transformation of immittance matrices (equation (2.3)). Note that if the modulation wavelengthl W 2p=b is a multiple of p, and when x is a multiple ofl, G becomes the identity matrix and X p q ðxÞ ¼ X p q ð0Þ as expected.

Eigenvalue problem and dispersion relation
The harmonics at the terminals of the nth unit cell are related by the ABCD transfer matrix or in the more convenient form where C r W ðV½r, I ½rÞ t is an infinite-dimensional vector that stores the amplitude of all time harmonics at x = rp; and P n is the ABCD matrix of the nth unit cell. We seek solutions of the form and G rr ¼ G r ¼ expðÀirbpÞ and zero otherwise. The condition (3.2) is equivalent to seeking a travelling wave solution of the form P 1 r¼À1 C 0r e i½vrtÀb r np , where b r W b þ rb. Therefore, (3.1) and (3.2) can be combined to give the eigenvalue problem (EVP) ð3:3Þ Note that T n W P n L acts as a translation operator, where exp (iβp) and C n are its eigenvalue and eigenvector, respectively. Furthermore, T n is a function of the operating frequency ω. If T n is indeed a translation operator, we should expect that the eigen-solutions are independent of the unit cell. The next property shows that T n is indeed independent of the unit cell (i.e. independent of n).
is a solution of The proof of property (3.1) is presented in the electronic supplementary material. Equations (3.4) and (3.5) imply that regardless of the unit cell used, the EVP will always result in a unique propagation constant β. Additionally, the kth component of the eigenvector changes in a way that is equivalent to the phase delay of the kth harmonic along a unit cell, which is equal to kbp. Therefore, the solution of the EVP (3.3) is independent of the unit cell.
The above result can be generalized to It is worth digressing here and discussing how the above formulation has the same mathematical foundations used to describe LTI systems. The discussion provides a deeper insight into the mathematical structure of the framework and highlights the physical similarities and differences between LTI and LTP systems. We refer to figure 1 that depicts a conceptual diagram of the applied flow. Figure 1a illustrates the logical flow one usually applies to describe LTI circuits. The basic relations of a two port circuit are given via differential and algebraic equations, leading to an LTI state space representation of the system. This in turn allows the output parameters to be related to the input ones via the convolution with the 2 × 2 impulse response matrix, shown in the left side of the figure by the A(t), B(t), C(t) and D(t) functions. The convolution operation with the impulse response can be easily diagonalized, where the set of complex functions exp (iωt) form the eigenfunctions of the convolution operator P Ã . Hence, it allows the system to be represented by simple royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 matrix multiplications at each frequency ω as shown in the middle panel of figure 1a. In other words, the frequency domain representation is merely an expansion of the input and output parameters in the eigen-functions basis of P Ã . For arbitrary time varying circuits, the output and input parameters can be related via the aid of the state transition matrix, rendering the analysis of such general systems a bit challenging. Fortunately for LTP systems, Floquet theorem tells us that the eigenfunctions are in the form P Q r expðiv r tÞ. Therefore, the input and output voltages and currents become infinitely long complex vectors that store the coefficients Q r and are connected via the multiplication with LTP ABCD matrix, as equation (3.1) emphasizes. In another words, the LTP ABCD matrix is the operator that connects the output and input vectors expanded in the LTP temporal eigen-functions.
When two port LTI networks are cloned to form a periodic structure, the network ABCD matrix does not change as one moves from one unit cell to the other. Such ABCD matrix maps the input applied to one of its terminals to the output, hence it represents the translation operator T (i.e. T = P). T is a linear operator in the two-dimensional complex space C 2 . Diagonalization of T is equivalent to finding its eigen-functions, which are nothing but the Bloch waves. The eigenvalues are conveniently represented as exp (iβp), where p is the unit cell length. Therefore for each ω, representing a temporal eigenvalue, there are two spatial eigenvalues ±β that diagonalize the spatial operator. The ordered pair (ω, ± β) describes the propagation of waves in a periodic structure at ω. On the other hand, property 2.1 implies that the ABCD matrix P n in a space-time modulated structure changes from one point to the other (i.e. depends on n). Nevertheless, property 3.1, along with equation (3.3) shows that T n does indeed act as a translation operator. Since the space of time harmonics is infinite, the diagonalization of T n (i.e. enforcing the Bloch condition) results in an infinite number of eigen-functions as the last panel of figure 1b illustrates. Note that if the modulation is removed L reduces to the identify matrix and T n in equation (3.3) becomes the ABCD matrix P n at frequencies ω r . Although the mathematical spaces for the LTI and LTP systems are different, the circuit based formalism suggests that space-time structures are formally equivalent to the well-known periodic structures. The infinite dimension of the LTP space and its translation operator highlights the richness in the spectrum (i.e. the eigen-functions) of space-time periodic structures, which eventually may lead to the possibility of breaking the physical limitations of LTI circuits. Ultimately, the properties of the LTI translation operator are dictated from the constraints of its basic circuit elements (stability, passivity, reciprocity and physical realizability), the properties of LTP operators hinge upon the fundamental characteristics of LTP elements, an area yet to be explored. Figure 1. Eigen-decomposition of (a) LTI and (b) LTP systems in the frequency and subsequently in the spatial domains.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 For any given frequency ω and in the limit of infinitesimal modulation strength, the A, B, C and D matrices become diagonal and the eigenmodes become the solution of the LTI system at frequencies ω r . The propagation constants β r satisfy the well-known dispersion relation [31,32] b r p ¼ cos À1 A r r ðvÞ þ D r r ðvÞ 2 : Furthermore, the eigenvectors become the Bloch waves. Figure  For the subsequent discussion, it is useful to introduce the harmonic shift operator S U .
. . t n that has the following effect i.e. S U shifts the vector C up by one position. Similarly S D W S U À1 shifts the vector down by one position.
If (ω, β) is a solution of (3.3), then the following is a general property of arbitrary space time modulated structures (Proof in electronic supplementary material). The previous property should not be surprising. The change v ! v þ lv and b ! b þ lb is equivalent to re-numbering the harmonics. Such property can be exploited to understand how the different modes at a given frequency behave by restricting the analysis to the zero-order branches only. Consider for example the dispersion relation shown in figure 2b. Suppose one is interested in the modal behaviour at v ¼ 1:5 a:u. Three modes are highlighted in the figure. Property 3.4 shows that the (2) eigenvector is a shifted copy of the (2') one. The (2') mode is at the intersection point of the 0 + and 1 − branches. As is already established when the modulation strength is relatively large, this point corresponds to the centre of a bandgap, where there is a strong interaction with the −1 harmonic [21,22]. Therefore, the voltage of this mode will contain non-zero entries only in the −1th and 0th royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 locations. Property 3.4 implies that at v ¼ 1:5 a:u:, the (2) mode has a zero entry in the 0th location and hence does not couple to the excitation at ω.
For the subsequent analysis, the kth eigenvalue and the corresponding eigenvector will be identified by the (k) superscript. In general, the mode voltage V ðkÞ and current I ðkÞ are related via the Bloch The general solution is the superposition of all modes v½n ¼

Boundary value problem
Consider the structure in figure 3 that represents a generic space-time modulated structure that is connected to a source and load. The incident (v (inc) , i (inc) ) and reflected (v (ref ) , i (ref ) ) waves appear on a transmission line of characteristic impedance Z 0 (usually a 50 V microstrip or coaxial TL) that connects the structure to a voltage source v s (t). At the input side, v s (t) = 2v (inc) (t). At the input of the first unit cell (x = 0), the boundary conditions are imposed by requiring the voltage and current to be continuous at the input and output ports. Therefore, for a sinusoidal excitation v s (t) = V s cos (ωt + ϕ) and a load impedance Z L , one gets where V inc1,k r W ðV ðkÞ r þ Z 0 I ðkÞ r Þ is the contribution of the kth mode to the wave incident on P1 and V inc = V s e iϕ /2. Equation (4.1) shows that the coefficients a k are such that the net effect of the branches is balanced with the excitation at frequency ω and they destructively interfere at any other harmonics. At the load side, where V inc2,k r W ðV ðkÞ r À Z L I ðkÞ r Þ=2 is the wave reflected from the load Z L when the output port is terminated in Z L . For the remainder of the manuscript, it is assumed that Z L = Z 0 (i.e. the structure is terminated in the reference impedance Z 0 ). Therefore, (4.2) implies that the a k coefficients are the ones that result in a null reflection from the load at all harmonics.
In practical applications and away from the luminal regime, only a limited number N H of harmonics are significant. For convenience, we consider N H to be an odd number 2N s þ 1, where N s ¼ 0, 1, 2, . . .. This selection allows the symmetric inclusion of harmonics from ÀN s to N s . Furthermore, we consider the number of branches to be 2N H to account for forward and backward waves. For each harmonic, the truncated versions of (4.1) and (4.2) provide two equations in the 2N H v s (t) royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 7 coefficients a k . Taking all N H harmonics into account, we end up with a system of 2N H equations in 2N H unknowns that can be written as From the solution of (4.3), the transmission coefficient of the rth harmonic S ðr,0Þ 21 can be calculated.

Results and discussion
To verify and demonstrate the use of the machinery developed in § §3 and 4, we will apply the framework to analyse two main structures. The first is a space-time periodic CRLH TL. Such an idealistic model allows a thorough analysis of the propagation behaviour that can be compared with state space time domain simulations. Next, we use the framework to reproduce and give insight into the nonreciprocial behaviour observed on a nonlinear right-handed transmission line (NL RH TL) that has been manufactured in our lab. Although the analysis carried out below is not meant to be exhaustive, it provides a useful and a systematic procedure to describe complex space-time periodic structures.

Composite right-left handed space-time modulated TL
The CRLH consists of N = 40 unit cells such as the one shown in figure 4, where the right-handed capacitance C R is space-time modulated. The first unit cell is connected to a source of impedance 50 V. The load is also assumed to be 50 V. KCL and KVL along with the current and voltage relations in the time domain are used to derive a state space model (SSM) of the circuit that can be written as where x is an N × 1 vector that stores the state variables (current in inductors and voltages across capacitors), A is a N × N matrix, B is a N × 1 vector that connects the input excitation to the states. The unit cell shown in figure 4 can be divided into three sub-units: (1) the LTI series impedance Z se ≜iωL R − i/ωC L , (2) shunt admittance 1/ωL L , and (3) the LTP admittanceỸ R . Hence, the ABCD matrix of the unit cell can be formed by cascading its three sub-units. Therefore, the different eigenvalues e ib ðkÞ p and the corresponding eigenvectors ðV ðkÞ , I ðkÞ Þ t are determined through the use of (3.3). Subsequently, when the TL is excited by a sinusoidal source of frequency ω, the boundary value problem (4.3) is solved and the modes coefficients a k are computed. The LTP dispersion relation when C R is sinusoidally modulated, i.e. C R ¼ C R0 ½1 þ M cosðvt ÀbnpÞ is obtained from the eigenvalues as depicted in figure 5a. The LTI dispersion relation (when M = 0) is superimposed to highlight the right-hand (RH) and left-hand (LH) regions. The LH region is in the low-frequency range, frequencies below 1 a.u., where the phase and group velocities are opposite [33,34]. For a balanced operation, the series and shunt resonances were both set to unity [34]. Figure 5b,c shows a close up view of the dispersion relation in the RH and LH regions, respectively. Different modes are highlighted and labelled. It is worth noting that whenever two branches appear to intersect, two complex conjugate propagations constants (γ = −iβ) are generated. The point of intersection represents the centre of a bandgap, where strong coupling with time harmonics may be significant. For instance consider figure 5b, where points 7 and 8 represents two eigenvectors that have two complex conjugate propagation constants.
The eigenvalues, when the frequency is at the centre of the RH BG (f ¼ 1:5 a:u.), are depicted in figure 6a. Unlike the first four, the higher eigenvalues result in complex conjugate propagation constants. This is not royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 surprising since they correspond to points inside BGs as figure 5b shows. For any given eigenvalue, the magnitude of the components of the corresponding eigenvector is plotted in figure 6b. For a given eigenvector (mode), the y-axis represents the strength of the rth harmonic. According to (3.6), the waveform inside the space-time periodic structure is the linear superposition of the different eigenvectors. Figure 6d plots the magnitude of the expansion coefficient a k . Clearly, the wave behaviour is dominated by the 8th eigenvector, which corresponds to one of the modes inside the BG of the main branch as illustrated in figure 5b. For such modes, figure 6b shows that the 0th and −1th harmonics are dominant and of the same order of magnitude, inferring a strong interaction between the fundamental and its −1th harmonic. Additionally, there is a small contribution coming from the 6th and 9th modes. The figure also shows that the 6th and 9th modes have significant components at the 0, + 1 and 0, − 1 enteries, respectively. Note that the behaviour of the two modes can also be deduced using property 3.4. Indeed, the equivalent mode of the 6th (9th) one on the main branch is inside the first BG. Hence the equivalent mode has non-zero values at the -1 and 0 entries. Since the eigenvector of the 6th (9th) mode is a down (up) shifted copy, it has non-zero values at the +1, 0 (−1, 0) entries, in agreement with figure 6b.  Figure 4. Unit cell of a space-time modulated CRLH TL. The right-handed capacitance C R is modulated as a travelling wave royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 To assess how accurately the LTP approach can predict the wave behaviour inside the structure, the waveform at the middle of the RH BG, at the three frequencies ω, v Àv and v þv are calculated using (3.6) and compared with the solution of the SSM. The time domain data obtained from the SSM simulation is transformed to the frequency domain, where the frequencies of interest are isolated. Figure 6c reports the amplitude of the three harmonics. As shown, there is an excellent agreement between LTP and SSM. Additionally, the amplitude of the main component at ω rapidly decreases as the wave penetrates into the structure, where it is scattered (mainly) in the −1 harmonic back to the source. Furthermore, there is a non vanishing contribution, coming from the +1 harmonic, as a result of the excitation of the 6th mode.
The same procedure is repeated but for ω at the centre of the LH BG (figure 5c). Unlike the RH BG, the incident and modulating waves are contra-directional. This is due to the left handedness of the CRLH in this regime. Therefore, the incident wave scatters in the v þv (blue shifted) as figure 7c highlights. The scattering, however, is not as strong as in the RH BG case. This is due to the smaller magnitude of the real part of the eigenvalue (figure 7a) and witnessed by the slight reduction of the amplitude of the fundamental component (figure 7c). It is expected that a modulation of L L or C L will result in a wider LH BG. Note also that unlike the dominant mode in the RH BG (mode 8, figure 6a), the dominant mode in the LH BG has a negative phase velocity as witnessed by the imaginary part of βp (mode 7, figure 7a).
Finally, the transmission coefficient is calculated over a wide frequency range that includes both the RH and LH BGs. Figure 8a,b presents the results, when the incident and modulation waves are co-and contra-directional, respectively. The figures show that LTP-based calculations are in a very good agreement with SSM. Furthermore, space-time modulation has the effect of attenuating the transmitted signal (−15 dB for the LH BG and −30 dB for the RH BG), compared to almost 0 dB when modulation is absent.  The capacitance of the varactor C v depends on the voltage across its terminals u(t) = u m (t) + u s (t), where u m (t) and u s (t) are the voltages due to the modulating input and signal, respectively. The current through the varactor is given by iðtÞ ¼ C n ðuÞ du dt : Since C v ðuÞ ¼ C v ðu m þ u s Þ % C n ðu m Þ þ u s dC n =dtj um , it can be shown that the current i s (t) due to the signal excitation is where C n (u m ) is the capacitance evaluated at u m (t), which is periodic with frequencyv. Hence, the system is linearized about the limit cycle steady state [35]. Figure 9c shows the varactor's equivalent circuit. R s models the ohmic losses in the semiconductor bulk, contact and bondwires, L s accounts for the inductance of the bondwires, and C represents the varactor capacitance. The different circuit parameters were extracted from measuring the S parameters at different bias voltage and fitting the response via the use of the vector fitting technique [36]. At low frequency, the microstrip line can be described by lumped circuits as in figure 9b. The p/2 microstrip line section is modelled as a lumped LC network, such that L = τ d Z c and C = τ d /Zc, where τ d and Z c are the delay and characteristic impedance of the microstrip, respectively. The lumped circuit approximation allows the convenient representation of the NL RH TL in an SSM form. In this case, the biasing circuit and blocking capacitances can be included as in figure 9b.

Dispersion relation
As a first step, the LTI dispersion relation of the structure is extracted from measuring the small signal S parameters for different bias voltages and compared to the circuit models. Figure 10 shows the dispersion relation curves of four bias voltages. Clearly, both the lumped and distributed circuit models are in agreement with measurement, confirming the validity of the lumped circuit model. In the presence of the modulating signal with frequencyv, the varactor capacitance C n becomes periodic with a period of 2p=v. Therefore, it can be expanded in Fourier series Since the amplitude of modulation is large (u m ≫ u s ), the DC capacitanceC 0 may deviate from the small signal value. In the subsequent analysis, the DC and first harmonic only (C 0 andC 1 ) will be considered. They are calculated from a time domain simulation of NLRHTL. The system differential equations are solved to compute the voltages u m across the different varactors. Consequently,C 0 and C 1 are calculated from the Fourier transform of C(u m ). Figure 11a shows the computed spectrum of C v across the 10th varactor. The DC capacitanceC 0 has increased from approximately 1.2 pF to 1.35 pF. The modulation strength M WC 1 =C 0 % 0:2 when the excitation is approximately 10 À 15 dBm. Figure 11b demonstrates howC 0 andC 1 change from one unit cell to the other. Although not constant, we will assume that bothC 0 andC 1 are constants and fixed to their average values. This assumption allows the application of the LTP formalism and captures the main essence of the system, as will be shown below.
The unit cell can be represented by the block diagram in figure 12. Hence the ABCD matrix of the unit cell is T ¼ T LTI T LTP T LTI . For the microstrip, the ABCD parameters are identical to the LTI counterpart, but calculated at each harmonic frequency ω r .
The LTP block T LTP represents the ABCD parameters of the shunt varactor, which is modelled by a shunt time periodic admittanceỸ sh ¼ ðZ se þỸ À1 Þ À1 ¼Ỹðe þ Z seỸ Þ À1 . The second term in the last expression is the inverse of a tridiagonal matrix and can be computed using closed-form expressions as in [37].
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 13 The speed of modulationn is determined from the phase ϕ ofC 1 , wherê which is expected to slightly deviate from the TL LTI speed. For a given modulation frequency f m and strength M, the dispersion relation can be determined from the solution of (3.3). Figure 13 depicts the dispersion relation for f m ¼ 1 GHz, and when the modulation propagates in forward (figure 13a) and backward (figure 13b) directions. As shown,n is very close to the LTI speed, suggesting that the LTP system is in the sonic regime [21,38].
To explore the interaction between the different modes and how they contribute to the overall propagation, consider the situation where the modulation and signal are co-directional. Using 14 modes, the eigenvectors are calculated as reported in figure 14b. The TL is excited with a sinusoidal signal of frequency f ¼ 2:82 GHz. As will be shown later, at this frequency, maximum non-reciprocity is observed. The plot shows the magnitude of the components of each eigenvector normalized to its maximum value. Modes of interest are the ones that strongly couple with the input excitation; hence they have significant components at ω (or the 0th harmonic as highlighted in figure 14b) and can potentially be excited. Additionally, the BVP (4.3) is invoked to compute the different a k values that in turn determine the strength of the excited modes as figure 14d shows. The waveforms at different frequencies are the superposition of the corresponding eigenvectors as presented in figure 14c and confirmed with SSM in figure 14a. Furthermore, figure 14e demonstrates that the waveforms can be approximated by the dominant eigenmodes (i.e. the ones that couple with the input excitation such that their expansion coefficients a k are non-vanishing). The signal at ω is significantly reduced at the output due to the interaction with its harmonics, mainly the −1 harmonic. The absence of bandgaps in the dispersion plots in figure 13 suggests that this type of interaction is passive in nature (i.e. β is   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367 The strong interaction between the fundamental and its −1 harmonic is apparent from the measured output spectrum ( figure 14f ). Here, the input port was fed by an RF source that was swept over a frequency range around 2.82 GHz and the output of the spectrum was measured by a spectrum analyzer. The spectrum shows that once the modulation is turned on, the interaction is mainly with the −1 harmonic. Note that modulation and its higher harmonics (1, 2 and 3 GHz) appear as spikes in the measured spectrum. When the modulation and signal are contra-directional, as figure 16a demonstrates, the eigenvalues are generally different from those calculated above. The dispersion relation shows an increase in the separation between the forward branches as in figure 16c. Hence, the incident wave is expected to strongly couple to the main branch, labelled by the mode number 14. Note that other higher modes, for instance mode 15, are wrapped back to the negative side once βp exceeds π. It is worth noting from the computed eigenvectors (figure 16d) that the 9th and 10th modes have a significant component at the 0th harmonic. However, due to the increased separation between the branches in the forward direction such modes are not excited. Therefore, one may conclude that when the signal and modulation are contra-directional the propagation is basically that of the LTI system. Indeed, the calculated a k coefficients (figure 16b) show that coupling is mainly with the 14th mode. Therefore, the mode couples with the forward main branch and the structure appears to be transparent in this mode of operation.

Transmission coefficient
Finally, the modes are superimposed and the transmission coefficient is calculated for both the co-and contra-directional modes of operations. The transmission coefficient at the fundamental frequency is calculated using SSM and the process is repeated over the 0-4 GHz range. Figure 17 shows that as a consequence of space-time modulation and the asymmetric interaction between harmonics in the forward and backward directions, strong non-reciprocity between the forward and backward propagation arises. As has been shown, this is due to the passive interaction between the fundamental mode and its lower harmonic when the modulation and signal are co-directional. For an input frequency of 2.82 GHz, the coherent length is 20 unit cells long. Hence, maximum energy is transferred to the −1th harmonic, reducing the signal at the output port. In the opposite direction, however, the distances between the forward branches are widened and the effect of modulation is negligible. Figure 17c,f reveals that such non-reciprocal behaviour demonstrates itself in the measured scattering parameters. The baseline, however, is reduced by approximately 30 dB due to the presence of the directional coupler. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210367

Conclusion
The time periodic circuit theory was exploited to derive some of the properties of the infinite-dimensional spatial translation operator of space-time modulated circuits. The modal behaviour of a generic spacetime periodic structure can be explained from the solution of the system eigenvalue problem. Additionally, we showed that the translation operator guarantees that solutions are invariant under spatial translation. Furthermore, it was shown that all points in the (β, ω) plane parallel to the modulation velocityn are equivalent in the sense that the eigenvectors are related by a shift operator. The waveforms inside the space-time periodic circuit and the time periodic scattering parameters were determined through the expansion of the total solution in terms of the eigenmodes, and after imposing the suitable boundary conditions. Two examples were discussed. In the first, a space-time modulated CRLH TL was studied using the developed approach and compared with time domain simulation. In the second example, the non-reciprocal behaviour observed on a nonlinear TL was explained. This was made possible via the extraction of circuit parameters from measurements that were then used to predict the wave behaviour inside the TL and its effect on the terminal properties. It was shown that the passive interaction between different harmonics resulted in an observed nonreciprocal behaviour, where the difference between forward and backward transmission coefficients S 21 − S 12 can be significant. The frequencies at which non-reciprocity occurred and its strength agree with time domain simulation and measurements.