A Generalized Model for Linear-Periodically-Time-Variant Circulators

Magnetic-free non-reciprocity based on linear-periodically-time-variant (LPTV) circuits has received significant research and commercial attention since it could revolutionize wireless communications. LPTV circuits are formed by two main components: linear-time-invariant (LTI) networks and periodically-modulated switches. The modulated switches are the core elements to break the reciprocity of LTI networks. To understand and design LPTV circulators, a universal and intuitive analytical model is required. However, such model does not exist as it is extremely challenging to accurately model and fully understand the LPTV behaviour of energy storage networks. To address this limitation, this work introduces a novel analysis method, which is broadly applicable to any LPTV networks, and validates it experimentally. The novelty of this methodology comes from two main contributions: (1) modelling of the switch as a resistor in parallel with a current-controlled current source; (2) the decomposition of the LPTV network into the linear superposition of two LTI networks. We apply this technique to model the exact behaviour of an LPTV circulator in the frequency domain.

phase-shifted N-path filters [20][21][22][23] , sequentially-switched delayed lines (SSDL) 17,25,26,[28][29][30] , and distributedly modulated capacitors (DMC) 19,27 . Generally, AMB enables symmetric circulators and requires a modulation frequency that is a small fraction of the filter centre frequency when high quality factor resonators are employed, hence reducing the power consumption of the active modulation process 38 . Moreover, the LTI network modulation could be implemented through switches or varactors. Switches are more attractive because they offer simpler implementation, larger modulation ratio, and easier integration 15 . Therefore, the AMB topology is the focus of this work. Figure 1a shows the proposed general AMB circulator topology 15 . The proposed circulator consists of two identical LTI networks with 3-fold rotational symmetry parametrically modulated by a switch matrix. Figure 1b shows three digital pulse trains with the same period, but a phase difference of 120° with respect to each other. These pulses are the modulation signals used to drive the switches of the black branch of the LTI network. The modulation signals of the red branch complement that of the black one such that the radio frequency (RF) input from the antenna is commutated between the two branches and no power is lost when the switches are in the off state in one branch.
Although the magnet-less non-reciprocal networks have been actively investigated, the analysis of the circuit in Fig. 1 remains challenging due to the complexity of LPTV networks based on modulated switches. Switches are extremely common electronic components; however, it comes as a surprise that, over the last eight decades since sampling mixers came into use in 1930s, the interaction between switches and LTI networks has not been fully modelled. Instead, only a few simple LTI networks modulated by switches have been analysed in either time-or frequency-domain. For example, SSDLs have been qualitatively described in the time-domain by treating transmission lines as non-dispersive elements 17,25,26,30 . Alternatively, sampling mixers have been analysed in the frequency-domain using switch-resistor or switch-capacitor models 38 where the output and input signals are readily related. More complex models have been developed for N-path filters. E. Klumperink et al. 39 15,31 Varactors L-C tanks [5][6][7][8][9][10][11][12] MEMS Resonators 16 Sequentially-Switched Delay Lines (SSDL) Switches Transmission Lines 25,26,28,29 Acoustic Delay Lines 17,30 Phase Shifted N-Path Filters Switches Transmission Lines [20][21][22][23] Distributedly Modulated Capacitors (DMC) Varactors Transmission Lines 19,27 Table 1. Summary of non-reciprocal device topologies based on LPTV circuits. (b) There are 120° rotational phase relationships between modulation signals (square wave pulses). T 0 is the modulation period (1/frequency) and T p is the pulse width. Duty cycle, α, is defined by the ratio of T p to T 0 . a complete overview on the history of N-path filters. In the modelling of N-path filters, assumptions such as loss-less sampling 40,41 and non-overlapping between switch clock signals 42 are made to simplify the analysis. These assumptions may offer a good approximation of the circuit behaviour, yet they limit the scope of the methods and are not applicable to a generalized circulator model. The most common techniques to analyse N-path filters are based on polyphase "kernels" [43][44][45] and conversion matrices 46 . A kernel is a single-path sub-circuit. This method decomposes the LPTV circuits into multiple polyphase kernels, then each kernel is analysed by invoking LPTV theory, and finally the kernel states are combined to obtain the circuit states. The LPTV theory reveals that 45 , the output spectrum, V o , of an LPTV network, is a summation of an infinite number of frequency-shifted and filtered input spectrum, V i , which is where H n (ω − nω 0 ) is known as "harmonic transfer functions (HTFs)", n is the harmonic index, and ω 0 is the modulation frequency. Finding HTFs is the goal of the analysis of LPTV circuits. The analysis of circulators based on phase-shifted N-path filters has used this kernel method 20,43 . However, such analysis is neither flexible as it cannot handle arbitrary LPTV circuits and arbitrary overlap between clock signals of different kernels, nor intuitive since the results contain an infinite number of translated signals in frequency. Hence, it is quite challenging to analyse a generalized circulator architecture by using polyphase kernels. By contrast, the use of conversion matrices is a more systematic approach, which relies on matrix-form expressions of the HTFs for basic electronic components, such as resistors, capacitors, and inductors. With this method, the analysis of arbitrary LPTV circuits can be performed in a similar way to that of LTI circuits. Nonetheless, it requires to include all explicit component-by-component conversion matrices in the output expression and thus it is not applicable to a generalized LPTV network 46 . A variant of the conversion matrix method was developed to perform the analysis on generalized LPTV networks in ref. 29 . However, it requires to use a large number of square Floquent Scattering Matrices (FSM), which are non-trivial to obtain and make this method more complex. Furthermore, the above methods do not provide a standalone model of a switch by itself without which it remains difficult to fully understand LPTV behaviour. When varactors are used instead of switches to modulate the LTI networks, then the LPTV circulators can be more readily analysed by coupled-mode theory (CMT) such as in N. A. Estep et al. 10 and R. Fleury et al. 47 . This is possible since only first-order intermodulation products have to be considered to obtain an accurate description of the LPTV circulator. However, the proposed models are not applicable to switch based LPTV networks since higher-order mixing must be taken into consideration for accurate description of the circulators.
To address some of the limitations of the aforementioned methods, this work proposes a novel approach to analyse an LPTV network by decomposing it into a finite number of LTI networks and accurately describing its behaviour with semi-analytical equations. This approach is enabled by modelling a periodically modulated switch with a resistor in parallel with a current-controlled current source (CCCS). The resistor represents the on-impedance of the switch, whereas the CCCS current is controlled by the current flowing in the resistor. Such model is independent of the networks the switch connects to and does not impose any restrictions on duty cycles, operating frequencies, and clock overlaps. By using this model of the switch, we are able to separate the description of the switch modulation function from the LTI networks and derive an analytical expression of the proposed generalized circulator as shown in Fig. 1. It is interesting to note that, through this method, the analysis of the LPTV circuit is reduced to the analysis of the comprising LTI circuits in which the switches are effectively on. The circulator model is validated by applying it to the description of a circulator synthesized using off-the-shelf-components. It is worth noting that the same overarching method can be extended to the analysis of any other LPTV circuits with only minor modifications of the component description.
The rest of the paper is organized in the following manner: the switch model in frequency domain is first presented and the interaction between the switch and any arbitrary LTI network is discussed; then such model is used to simplify the description of LPTV networks and derive the semi-analytical S-parameters of the generalized circulator in Fig. 1; finally, the generalized circulator analytical model is validated by comparing its results to the experimental data.

Standalone Switch Model in Frequency Domain
This work proposes an innovative and illustrative model that describes the dynamic behaviour of a switch periodically toggled on and off using a resistor in parallel with a current-controlled current source (CCCS) (see Fig. 2a and refer to Methods Section for the derivation). The resistor models the frequency dependent (ω) on-impedance of the switch, z S (ω). The CCCS models the effect of periodic modulation. The current of the CCCS is controlled by the current flowing through the resistor. The quantitative relationship between the CCCS and the current through the resistor is presented in this section and is derived analytically in the Methods Section.
Given a carrier frequency, ω RF , and a modulation frequency, ω 0 (=2π/T 0 , where T 0 is the period of the modulation signals), the steady-state current through the switch (and the network it connects with) contains intermodulation products at ω RF ± nω 0 , where n is an integer. Termed as "harmonic index" in this work, n can theoretically span from negative to positive infinity. However, considering that most systems have a finite bandwidth (BW) and that the magnitude of the intermodulation products is proportional to 1/n, it is enough to include a finite number of intermodulation products in the analysis. If the maximum harmonic index is N, there are (2N + 1) frequency components, i.e. from (ω RF − Nω 0 ) to (ω RF + Nω 0 ). Roughly, N only needs to be greater than BW/ω 0 , which ensures that all effective intermodulation products fall within the band of the system regardless of the relative position of the carrier frequency with respect to the band centre frequency. More broadly, the larger the value of N, i.e. the number of higher-order intermodulation products included in the analysis, the more accurately the www.nature.com/scientificreports www.nature.com/scientificreports/ model describes the behaviour of the network. Conversely, when N = 0, the number of frequency components is 1 and only the carrier frequency is analysed, which models the network without modulation present (i.e., the switches are always on).
The LPTV currents generated by periodic switching are represented in the frequency domain by vectors comprised of complex current phasors (frequency information is omitted but implied by the positions and subscripts of phasors), as shown in Fig. 2a. As explained in the Methods Section, the correlation between I S0 , I Z0 , and I C0 can be expressed as where α is the duty cycle, θ is the phase delay of the modulation signal, Y is a diagonal spectral admittance matrix of the switch defined in Eq. (33), and C(α, θ) is a (2N + 1)-order matrix dictating the mapping relationship in frequency conversion defined in Eq. (34), whose elements are strongly related to the complex Fourier transform coefficients of the switch's periodic behaviour in time-domain. In the two extreme cases of no modulation, α = 0 or α = 1, C(0, θ) is a (2N + 1)-order negative identity matrix, while C(1, θ) is a (2N + 1)-order zero matrix, which respectively produce: z S (ω) is the spectral on-impedance of the switch, which is inversely related to its on-admittance, y S (ω). (b) The interaction between a switch and a general LTI network, and its corresponding equivalent circuit. The circuit is driven by a voltage source with internal impedance, Z 0 . The source impedance is called "Port". (c,d) The two LTI circuits and corresponding transfer functions used to analyse the LPTV circuit in (b). The italic lower-case "i" means current scalar, while the bold upper-case "I" represents current vector, with italic upper-case "I" referring to the current elements in the vector.
www.nature.com/scientificreports www.nature.com/scientificreports/ Intuitively, Eq. (4) represents the behaviour of an open switch for zero duty cycle, while Eq. (5) indicates that no modulation is present for 100% duty cycle, i.e. the switch is always closed. These two prosperities are independent of ω 0 and N.
When the periodically modulated switch interacts with an arbitrary LTI network driven by a voltage source, the switch can be replaced by its equivalent model shown in Fig. 2b. By applying linear superposition theory, the analysis of the circuit in Fig. 2b can be decomposed into the analysis of two LTI circuits as shown in Fig. 2c,d, which consider the contribution of a voltage source (corresponding to the driving source) and a current source (corresponding to the CCCS generated by the switch), respectively. The contribution to the currents at the locations of interest, for example, the port (see Fig. 2b) and the equivalent resistor in this work, can be characterized by four transfer functions, P 0 (ω), R 0 (ω), Q 0 (ω), and T 0 (ω) as defined in Fig. 2c and d. It is worth noting that the circuit in Fig. 2c is equivalent to the circuit in Fig. 2b without modulation. As proven in the Methods Section, the CCCS pumps currents into the equivalent resistor, which adds to the circuit current without modulation (Fig. 2c) to yield the following steady-state current: is the current through the equivalent resistor when the switch is always on, A T0 represents the absorption matrix that depicts the ability of the switch equivalent resistor to absorb current from the CCCS at different intermodulation frequencies, as defined in Eq. (37). Similarly, the modulation modifies the port current, I P0 , in a way that is the current through the port without modulation, A Q0 represents the absorption matrix that depicts the ability of the termination port to absorb current from the CCCS at different intermodulation frequencies, as defined in Eq. (38).
Looking at Eqs (6) and (7), it is interesting to note that the steady-state currents of the resulting LPTV circuits in which the switches are used (Fig. 2b) are composed of two parts: the initial current without modulation and the contribution from the CCCS due to modulation. Therefore, the analysis of the overall current in the LPTV network can be easily conducted by deriving the transfer functions of two LTI circuits as shown in Fig. 2c,d, i.e. by using P 0 (ω), R 0 (ω), Q 0 (ω), and T 0 (ω). In summary, the proposed switch model permits treating LPTV circuits (in Fig. 2b) as the linear superposition of two LTI circuits' states. This point will be illustrated further in the following section.

Generalized Analysis of the Circulator Circuit
In RF domain, S-parameters are commonly used to describe the network behavior 38 . Our ultimate goal is to derive the closed-form S-parameters of the proposed circulator in Fig. 1a. Due to the 3-fold rotational symmetry of the circulator topology, there are only three independent S-parameters: S 11 , S 21 , and S 31 . Among them, |S 21 | is referred to as insertion loss (IL), while |S 21 |/|S 31 | is called isolation, both being expressed in dB. These two parameters are the most important ones in characterizing the circulator performance. Generally speaking, one wants to minimize IL while maintaining high isolation. To compute the S-parameters, a single-tone voltage source (with carrier frequency, ω RF ) excites the circuit at Port 1 (Ant), as shown in Fig. 3a. From the above analysis of the switch, we can replace the switches with resistors and CCCSs (see Fig. 3b). This circuit can be solved by applying linear superposition theory, which corresponds to the analysis of the two core LTI circuits as shown in Fig. 3c,d.
For convenience, we first define p-and r-transfer functions for the circuit in Fig. 3c: as well as q-and t-transfer functions for the circuit in Fig. 3d: Intuitively, p-and q-transfer functions describe the fraction of current that each port can absorb from the voltage source and current source, respectively, while r-and t-transfer functions indicate the respective capability of the voltage source and current source to pump currents into the equivalent resistor. It is interesting to note that, P I (ω) and P IJ (ω) are related to the S-parameters (S S S , , ) of the circuit in Fig. 3a without modulation (namely, when switches are always on) in the following way: because of the symmetric circuit topology.
For brevity in the derivation of the circulator model, Table 2 shows the defined naming convention for currents at different locations in the circuit. By using the naming convention, we can re-write equations for each CCCS based on Eq. (1) as: i i , and θ i = (i − 1)2π/3, i = 1, 2, 3. Similar to Eqs (6) and (7), the current through equivalent resistors and ports are expressed as the sum of the initial current without modulation and the contribution from CCCSs due to modulation. With the assumption of 3-fold rotational symmetry of the LTI networks, these currents can be expressed as  where i, j, and k are the permutation of 1, 2, and 3, h and l are the permutation of 1 and 2, and A X represents the absorption matrix that depicts the ability of ports or switch equivalent resistors to absorb current from CCCS, and it is written as  Steady-state current of i-th port without modulation, defined as    I  I  I  I  I  I , ,0 where i = 1, 2, 3. All components are zero except its element at the carrier frequency. The current element at the frequency (ω RF + nω 0 ) in I Pi Steady-state current of ih-th switch's equivalent resistor without modulation, defined as  www.nature.com/scientificreports www.nature.com/scientificreports/ (23) The same sub-matrices in matrix A T ∼ and A Q ∼ have the same colours to help the reader recognize how these two matrices are arranged. Combining Eqs (17) and (18), I Z can be calculated as where Ĩ d is a 6(2N + 1)-order identity matrix. Substituting Eqs (24) and (17) In Eq. (25), I P (0) and I Z (0) represent the circuit response to the exciting voltage source without modulation. These matrices are easy to compute by LTI theory. The non-zero elements in I P (0) and I Z (0) are: The S-parameters for the proposed circulator can be obtained by where I Pi,0 (i = 1, 2, 3) is the carrier frequency steady-current current at i-th port as defined in Table 2. Equations (25)- (27) constitute the complete semi-analytical model of the proposed circulator as shown in Fig. 1. As evident from the derivation process, the model is agnostic of the LTI network and the switch performance and only needs linear matrices to describe the behaviour of these components. This is the most powerful aspect of this model, which dramatically simplifies the description of LPTV networks.

Experimental Verification of the Circulator Model
To verify the previously described circulator model, we set up a circulator circuit as shown in Fig. 4a and b by using MiniCitcuits discrete RF high isolation switches, ZFSWHA-1-20+, and bandpass filters (BPF), SBP-21.4+. Figure 5a,b show the frequency responses of the selected switch and filter, respectively. The filter centre frequency is 21.4 MHz and its 3dB-bandwidth (BW) is 7.4 MHz. The insertion loss (IL) of the filter and switch at 21.4 MHz are 0.81 dB and 0.61 dB, respectively. As shown in Fig. 5a, since the magnitudes of y 1 and y 2 are at least 100 times smaller than y S , it is reasonable to neglect their shunt components in the model of the switch and only take its series admittance, y S , into account. As for the modulation signals, three pulse trains were generated by two synchronized 2-output pulse generators, Agilent 81110 A, which are then fed to a hex inverter, 74HCT04, to generate 3 complementary pairs of square waves, as shown in Fig. 4c. The phase differences between these signals were monitored by a 4-channel oscilloscope, Agilent DSO6014, and are automatically adjusted to be 120° by a MATLAB program. By sweeping the modulation frequency and duty cycle of the signal applied to the switches, we obtained a family of responses for the circulator. Figure 6a shows an overlapped contour map between circulator isolation and IL for different modulation characteristics. The shaded area is the region where modulation parameters are chosen to offer both good isolation and IL. It can be concluded that trade-offs between IL and isolation exist because the best IL and the largest isolation occur for different modulation parameters. We are not interested in explaining the reasons behind these trade-offs, which can be found in C. Xu et al. 31 , but rather focus on the validation of the proposed theoretical model. We chose two example points, i.e. Point A and B in Fig. 6a, to compare the theoretical and experimental results, as shown in Fig. 6b and c, respectively. Point A is an example of the circulator operating in the shaded region, which offers IL of 5.7 dB and isolation of ~15 dB over a 3-dB bandwidth of ~1.5 MHz. Point B is another set of parameters providing slightly lower IL (6.4 dB) yet much larger isolation (>30 dB) for some particular frequencies. For both cases, S 11 is well below −10 dB. As a corner case, Fig. 6d shows the circulator response for α = 1, which effectively turns on only the normal set of switches (black ones in Fig. 1a) and thus blocks the complementary branch of LTI network. This configuration is reciprocal since the switches are kept either on or off and thus no modulation effect is imparted to the networks. The reciprocity is validated by S 21 = S 31 .
It is observed that the theoretical results overlap well with the experimental ones in all three cases (Fig. 6b-d). To guarantee such good agreement between theory and experiment for the first two cases (Fig. 6b,c), a maximum harmonic index N = 14 is used in this work, which is slightly greater than BW/ω 0 (equal to 9.3 and 13.7 for Point A and B, respectively). Minor discrepancies exist between the model and experiments. These originates from the deviation of the circulator from 3-fold rotational symmetry due to differences in the BPFs transfer function as evident from the measured S 21 and S 31 shown in Fig. 6d. Figure 7 shows how the choice of smaller values of N (=0, 1, 5, 10) affects the prediction accuracy of the proposed circulator model. In particular, N = 0 (Fig. 7a,e) means no intermodulation products are accounted for, which essentially represents no modulation of the circulator. Therefore, Fig. 7a,e have the same theoretical (2019) 9:8718 | https://doi.org/10.1038/s41598-019-45013-5 www.nature.com/scientificreports www.nature.com/scientificreports/ response, in which S 21 = S 31 , regardless of the modulation parameters. As expected, N = 1, i.e. only ω RF and ω RF ± ω 0 being taken into consideration, cannot accurately predict the circulator behaviours, as shown in Fig. 7b,f. It is also interesting to note that, N = 5 is large enough to predict the in-band performance of the circulator, because the magnitudes of the intermodulation products are proportional to 1/n. By comparing Fig. 7c,d and Fig. 7g-h (N = 5, 10) to Fig. 6b,c (N = 14), respectively, it is worth noting that increasing N improves the accuracy in predicting the out-of-band response of the circulator. Interestingly, for circulator operated at Point A, the prediction with N = 10 (Fig. 7c) has no significant difference from that with N = 14 (Fig. 6b). Hence, selecting N slightly greater than BW/ω 0 should be considered as a general rule of thumb to obtain an accurate description of the circulator network. As for the third case where α = 1 (Fig. 6d), any N ≥ 0 produces the same results from the circulator model as expected from Eqs (4) and (5) Further details of reproducing plots herein can be referred to Supplementary Information for "A Generalized Model for Linear-Periodically-Time-Variant Circulators".
It is worth noting that the model does not make any assumption on the nature of the control signals. For example, the clock signal of any switch overlaps with that of two others at any given time in both cases that were www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ analysed. Furthermore, even two different duty cycles (44% and 56%) were used for the case in Point B. This further proves the power of our generalized circulator model.
In conclusion, we have proven the validity of our circulator model, which can be used to predict circulator behaviours accurately, given specific switch performance, LTI response, modulation frequency, and duty cycle. It should be emphasized that no particular LTI network property is used in the derivation of the circulator model. Therefore, the model is broadly valid for all LTI networks with 3-fold symmetry. Moreover, it is easy to extend the model to account for general LTI networks and arbitrary switch configurations. It is also worth pointing out www.nature.com/scientificreports www.nature.com/scientificreports/ that a simplified model of the switch was used in this theoretical analysis for the experimental validation of the circulator model. A perfect open was used to model the switch in the off state and shunt parasitics were neglected in both the on and off states. Although beyond the scope of this work, the finite off-impedance and shunt parasitic can be separated from the switch and treated as part of LTI networks so that no modification to the form of the semi-analytical model is needed. Similarly, the switching speed was assumed to be infinite, but this is not a requirement for the proposed methodology. Any real switch signals in the time-domain can be described by complex Fourier transform coefficients (as explained in Methods section), hence easily incorporating the impact of non-ideal switching into the model.

Conclusions
This work rigorously derived and experimentally validated a generalized frequency-domain semi-analytical model for LPTV circulators. The overall model is enabled by describing the behaviour of the switch using a resistor in parallel with a CCCS. With the help of the switch model, the analysis of LPTV circuits can be reduced to the linear superposition of LTI circuits. Although applied to a particular circulator network, the proposed model can be used to describe any LPTV circulators.

Methods
The switch behaviour is described by building a model of its I-V relationship under modulation. In the off-state, the switch presents a small admittance when in the off-state and a large one, y S (ω), when in the on-state. Since a switch with non-zero off-admittance can be treated as equivalent to a zero off-admittance switch in parallel with a component that has the same admittance of the switch's off-admittance incorporated in the LTI network, we can consider the switch off-state admittance to be zero without imposing any particular restrictions on the validity of the model. Therefore, without loss of generality, we can assume that the switch admittance spectrum is equal to y S (ω) Siemens when in on-state and 0 Siemens when in off-state. When a voltage of u S0 (t) = U S0 exp(jω RF t) at the carrier frequency, ω RF , is applied across the switch (Fig. 8a), the corresponding excited steady-state current, i S0 (assuming i S0 exists), through the switch in the time domain is  (Fig. 9b), is used to describe y(t) in this work so that: where α is the duty cycle, and τ = T 0 •θ/(2π) is the time delay of the modulation signal. In this case, c n = (1 − e −j2αnπ )/(j2nπ) and particularly, c 0 = α. Combining Eqs (28) and (29) produces, , infinite intermodulation products (ω RF ± nω 0 ) are generated due to the mixing effect. The first term Ohm's law as if a voltage U S0 exp(jω RF t) was applied across an equivalent resistor of z S (ω) Ohms, while the rest can be treated as current controlled current sources (CCCSs) with a signal equal to C n RF 0, 0 dependent on i Z0,0 at the intermodulation frequencies, (ω RF + nω 0 ), in parallel with the equivalent resistance of the switch, as shown in Fig. 8b. Equation (31) also implies that once there is current Z n RF 0, 0 flowing through the equivalent resistor, CCCSs are generated accordingly. Upon the interaction of current sources with external circuits, the current sources inject currents into the equivalent resistor, then causing another iteration of mixing process. This process repeats until the circuit currents and voltages converge at steady state. In this state, there are infinite intermodulation frequency currents flowing through the equivalent resistor in parallel with infinite intermodulation frequency CCCSs, which are represented in a compact form as shown in Fig. 8c. For convenience, the vectors comprised of current phasors (frequency information is omitted but implied by the positions and subscripts of phasors) are used to represent the circuit states. In this work, we call them phasor vectors. The steady-state currents, I Z0 and I C0 , can be decomposed into the sum of infinite iterative currents, I Z Now, let us consider the interaction between a switch and a general LTI electrical network, which is driven by a voltage source, as shown in Fig. 10. The circuit states can be completely described by phasor vectors, I P0 , I Z0 , and www.nature.com/scientificreports www.nature.com/scientificreports/ , respectively. As explained in the main manuscript, it is not necessary to compute all infinite number of intermodulation terms in the circuit state vectors. Therefore, in Fig. 10, only (2N + 1) intermodulation frequencies, i.e. from (ω RF − Nω 0 ) to (ω RF + Nω 0 ), are taken into account. Fig. 10 also shows the correlation between iterative currents generated by the frequency mixing due to the switch modulation. To express the modulation effect, rewriting Eq. (31) in a matrix form produces where Y is a diagonal admittance matrix related to the switch's on-admittance, y S (ω) = 1/z S (ω), which is given by

Data Availability
Data presented in this work will be made available by the authors upon appropriate request.