The Input–Output Equivalent Sources Method for Fast Simulations of Distributed Nonlinearities in Bulk Acoustic Wave Resonators and Filters

This work presents a new method to analyze weak distributed nonlinear (NL) effects, with a focus on the generation of harmonics (H) and intermodulation products (IMD) in bulk acoustic wave (BAW) resonators and filters composed of them. The method consists of finding equivalent current sources [input–output equivalent sources (IOES)] at the H or IMD frequencies of interest that are applied to the boundary nodes of any layer that can contribute to the nonlinearities according to its local NL constitutive equations. The new methodology is compared with the harmonic balance (HB) analysis, by means of a commercial tool, of a discretized NL Mason model, which is the most used model for NL BAW resonators. While the computation time is drastically reduced, the results are fully identical. For the simulation of a seventh-order filter, the IOES method is around 700 times faster than the HB simulations.


I. INTRODUCTION
HE analysis of nonlinearities occurring in passive devices, although being weak, are payed a growing attention, due to their effects in the high demanded performance of nowadays receivers [1]. Passive intermodulation (PIM) that can produce receiver blocking, might be caused by components placed close to the antenna if this is shared by the transmitter and the receiver. PIM might arise by different causes [2], such as degradation of materials, bad RF connections, thermal effects or intrinsic nonlinearities arising inside of the components: Manuscript received Oct. 2020. This work was supported in part by the Catalan Gov. through grant 2017 SGR 813, and by the Spanish Gov through grants TEC2017-84817-C2-2-R, TEC2017-88343-C4-2-R. C. Collado, J. Mateu, M. González-Rodrigues, R. Perea-Robles are with the Signal Theory and Communications Dept. from Universitat Politècnica de Catalunya (UPC), Barcelona 08034, Spain (e-mail: collado@tsc.upc.edu). R. Aigner is with Qorvo, Inc., Orlando, FL antennas, isolators, connectors, filters, etc. All of them must meet low PIM specifications becoming a challenging task for the RF designers of such devices.
Proper nonlinear models and simulation techniques help to reduce the time-to-market if those models can predict, before manufacturing, the PIM performance of these devices, which, in a general sense, are size comparable to the wavelength. A proven good approach is to work with distributed nonlinear models in which the nonlinear effect is locally described, allowing to find shape and size-independent nonlinear material parameters that can be later used for the design of any device. The main drawback of this approach is that the device must be discretized in many unit-cells describing the nonlinearities locally, and it leads to a high computation time, especially if the device has a lot of potential distributed nonlinear components, such as multiplexers, which includes many resonators.
In comparison with highly nonlinear components, such as diodes or transistors, passive components are weakly nonlinear. For weak nonlinearities we understand that the generated H and IMD in a nonlinear circuit are significant, but there is not significant saturation or detuning at the fundamental signals. Under these premises, nonlinear analyses can be simplified if we do not constrain the analysis to the strict accomplishment of power conservation [3], [4].
This document describes the development of a new method of analysis for the simulation of H/IMD generation in passive distributed circuits. The method is called the Input-Output Equivalent Sources (IOES) method and it allows to perform extremely fast and robust simulations of distributed weak nonlinear circuits. But the simplification of weak nonlinearities, this method is mathematically rigorous, it is not restricted to narrowband scenarios, and remix effects are also considered.
Though it could be applicable to any technology, such as superconducting or ferroelectric devices, this article describes the application of this method to electroacoustic devices, particularly, to Bulk Acoustic Wave (BAW) resonators and filters, due to its impact on current communication devices.
Electroacoustic technology has become the main driving technology for the development of nowadays portable communication devices, with very complex multisystem RF stages comprising many stringent filtering components sharing a single antenna [5]. Despite this technology is by far the best current solution to locate many high-performance filters in a very reduced area, it suffers from intrinsic nonlinearities that may set limits to the performance of the receivers.
In the last two decades, a lot of progress has been made to model and simulate the nonlinear behavior of the electroacoustic devices [6]- [8]. The most complete description of a nonlinear model of BAW resonators is found in [8]. The models are based on the nonlinear constitutive equations of the piezoelectricity, which are time-domain locally defined and then applied at any point of the circuit. This model has been extensively used for finding the derivative material constants of the piezoelectric material appearing in those equations. This characterization process is usually performed by analyzing measurements of one or more resonators, through Harmonic Balance (HB) simulations using commercial circuit simulators.
Recently, it has been demonstrated that peculiar nonlinear effects observed in some resonators, are due to further materials beyond the piezoelectric, such as the Silicon Dioxide (SiO2) commonly used in the acoustic reflector of the Solidly Mounted Resonators (SMR) [9], [10]. That means that, in some cases, a comprehensive model must discretize, not only the piezoelectric layers, but also some other potential layers, which leads to nonlinear problems with a lot of nonlinear unit-cells.
When simulating a single resonator, the simulations are time consuming, but the problem is enough small to keep the simulation time into a few minutes (of course depending on the frequency points to be analyzed). Although convergence problems may appear depending on the input power and the separation between fundamental frequencies in a two-tone experiment, this analysis method has been proven its efficiency and usefulness. Nevertheless, when simulating more complex devices, such as filters or duplexers, the problem size is too big to keep the computation time into reasonable values.
Several alternatives have been developed during the last years to reduce the computation time significantly. Under certain assumptions the number of nonlinear sources of a distributed model can be reduced going to the most simplified model consisting on having only two nonlinear sources in an equivalent nonlinear Butterworth Van Dyke (BVD) model [11], [12]. These models are much faster than those distributed, and they have been satisfactorily used in the design of filtering devices to evaluate or minimize the nonlinear performance in the early stages of the design process. Nevertheless, those lumped models are accurate only in a limited bandwidth or for certain nonlinear scenarios that are usually required by the customers, they do not have the capability of predicting unexpected nonlinear effects in a broadband scenario due to the lack of the distributed nature of the problem.
The analysis of distributed models with weak nonlinearities has been previously proposed in [13]. The method described in [13] assumes weak nonlinearities and it is based on a frequency domain description of the nonlinear equations for each H/IMD of interest, which allows to solve the whole discretized circuit avoiding HB techniques. As a result, the computation time is significantly improved, and the convergence problems are circumvented. However, this analysis involves the inversion of matrices that are still too big for complex devices like filters or multiplexers.
The aim of the IOES method is to avoid the inversion of big matrices to linearly solve the distributed circuit at the fundamental frequencies, which allows us to find the nonlinear sources at a given harmonic or IMD, and even to solve the whole circuit at these new frequencies. This can be achieved without losing mathematical rigor and without simplifications.
Section II briefly reviews the nonlinear constitutive equation and the well-known distributed Mason model to frame the complex distributed problem and state the framework we have used to compare the new IOES method with conventional HB techniques. Section III describes an alternative method that avoid HB techniques to simulate weak nonlinear circuits. This technique, based on a Volterra series analysis [14], was presented in [13] and it is briefly described in this article as an intermediate step towards the new IOES method, which is useful to unveil some concepts and equations that are later used. Section IV is the core of this work. It describes the IOES method including a new description of the 4-port ABCD matrix of a nonlinear Mason circuit and all the algebraic manipulations that are required to figure out the equivalent sources. These equivalent sources allow to solve big distributed weak nonlinear circuits much faster than the conventional HB techniques and the method described in [13]. Section V illustrates some examples that are used to validate this new technique. The first example sets a comparison between the IOES method and HB techniques demonstrating that both techniques provide exactly the same response, being the IOES method more than two orders of magnitude faster. The second example shows previously published measurements in [13] and the corresponding IOES simulations just to briefly illustrate that the new method is perfectly flexible to include external components, such as the measurement system previously characterized, as it is always required in real situations. The third example shows new simulations of a 7 th filter to outline the main advantage of this technique, that is, the bigger the nonlinear problem, the faster the IOES method is in comparison with HB and [13].

II. CONSTITUTIVE EQUATIONS AND NONLINEAR MODEL
For completeness of the article we describe very briefly the constitutive equations and the distributed nonlinear Mason model, which had been extensively detailed in [8], [15]- [17].

A. Constitutive Time-Domain equations
The constitutive electroacoustic equations of the piezoelectricity [8] describe the interaction between stress T, strain S, electric field E and electric displacement D using the constants c E , e and  S , stiffness, piezoelectric and permittivity respectively: where c D is the stiffened elasticity defined as = + 2 , and = + , = . (4) The equations for the non-piezoelectric layers are much simpler where , 2 and 3 are the elastic constant of a given material and its nonlinear derivatives.

B. Distributed Nonlinear Mason Model
The circuit model is the same as the one used in [8]- [10], [15]- [17]. Equations (1)-(4) are implemented using the distributed Mason model where the voltage represents the force, and the current represents the particle velocity. The equivalent circuit discretizes the piezoelectric layer into many unit-cells (slabs of thickness z) accounting for electro-acoustic interactions and wave propagation in the thickness direction, which is modeled as an acoustic transmission line (ATL) with distributed parameters [18]: where , , are the resonators' area, mass density and viscosity. Figure 1 shows the unit-cell [8]- [10], [15]- [17] modelling a thin z slab including the nonlinear sources Vc and Tc of (4).
Note that Tc is scaled by the area Fc=-A·Tc since the voltage represents the force. The whole layer is modelled cascading a finite number of unit-cells. The number of unit-cells will depend on the highest frequency of interest, being 100 unit-cells a good compromise between broadband reliability and computational time speed.
The acoustic wave propagation through the nonpiezoelectric layers can be modelled with a conventional T-network circuit vA being the phase velocity and = 1 2 2 √ 2 the attenuation constant [18]. If the layer is considered nonlinear, the acoustic field magnitudes must be calculated at each point along its thickness and a discretized model, with at least 100 unit-cells per layer like the one in Fig 2b, is used with nonlinear voltage sources Tc and the T-circuit corresponding to a small z thickness.
The equivalent circuit of the whole resonator is depicted in Fig.3. The piezoelectric, in this example Aluminum Nitride (AlN) and some layers that may affect the nonlinear behavior, for example the silicon dioxide (SiO2) layers, are divided into many unit-cells, whereas the other layers, considered as linear, are simulated using the circuit of Fig.2a It has been demonstrated [13] than the circuit of Fig. 3 is very useful to find the nonlinear material parameters that characterize the nonlinearities of SMR-BAW resonators. Those parameters are later used to predict the behavior of new devices with other areas, layer thicknesses or stack configuration. In [13], the circuit is solved by HB techniques using Symbolically Defined Device (SDD) components of ADS [19] to implement the time-domain nonlinear equations (1)- (5). However, the required discretization degree, around 400 nonlinear unit-cells for a single resonator, makes this analysis method not very agile for simulating filters or multiplexers.

NONLINEARITIES
An alternative for solving this kind of big distributed problems is to take advantage of the weak nonlinear behavior of these passive devices. This section describes the approach outlined in [13] and introduces some useful concepts that will be later used in the IOES method, described in Section IV.

A. Distributed Problem
As a distributed effect, the H's and IMD's are generated along the stack depending on the field magnitudes (independent variables in (1)-(5)) at any point and evaluated at the frequencies that can generate a given H/IMD. Those frequencies are the frequencies of the external sources that feed a given device, but also could be H's or IMD's that could contribute to other H or IMD due to the so-called remix effects.
At any point of the stack, the nonlinear terms of (1)-(5) generate new frequency components and the final field distribution at these frequencies will depend on the boundary conditions that the stack of materials and the external loads (for example source and load impedance in a 2-port device) impose. Finally, the output power signal that couples to the external load (or loads) is the required "useful" magnitude from a practical point of view. Therefore, two fundamental aspects are required:  The field magnitudes that can generate a given H/IMD must be found at any point of the stack to calculate the distributed nonlinear sources.
 The whole circuit must be solved at the new frequencies to find out how the generated signal along the stack couples to the output load.

B. Weak Nonlinearities and Frequency-Domain Equations
When analyzing passive devices with intrinsic nonlinearities caused by inherent physical mechanisms, as the piezoelectricity is in this case, the generated harmonics or IMD's are typically low enough to keep unperturbed the response of the device to the fundamental frequencies. That means that no saturation or detuning effects are significant at moderate power levels. Despite of that, the generated spurious signals are perfectly measurable and even can cause receiver blocking on certain scenarios.
Under this assumption, frequency defined models can be developed and even the used HB algorithms might not be necessary. For example, it can be easily demonstrated [1] that the nonlinear sources at the second harmonic (H2) according to (2) will be: where the field magnitudes ∆ 2 1 , ∆ 2 1 are phasors at the H2 frequency and 1 , 1 phasors at the fundamental frequency

1.
In a two-tone experiment with 1 and 2 as fundamental frequencies, the nonlinear sources of the IMD3 at 21-2 for example, without considering remix effects, will be [1]: If remixing effects are considered (Appendix I), the frequency-defined equations become a little more complicated since several combinations of generated frequencies may contribute to a given IMD.

C. Solution of the Frequency-Defined Nonlinear Circuit
The analysis is based on the admittance matrix description of the discretized circuit of Fig. 3 and the concept is quite simple. The whole circuit is linearly solved at the fundamental frequencies and the field magnitudes (equivalent current and voltages) are found at any node of the circuit. Then, the nonlinear sources at a given H/IMD are calculated and the circuit is solved again at that desired frequency. That provides the output power at the load.
This concept, described in [13], was implemented in Matlab ® . In terms of computation resources, the time required to analyze a one-tone experiment in a single resonator was around four times lower than the time required by a HB analysis using ADS, whereas the time reduction with the new methodology presented in this paper is by several orders of magnitude. Note that in [13], the size of the matrices to be inverted are quite big. Since the circuit is fully discretized, many nodes are required to get directly the strain S and the electric field E at each position that is required to calculate the nonlinear sources. Just as an example, the required matrices size of the resonator described in [13] was around 2000x2000 and several equation systems had to be solved (inversion of matrices) for calculating the H/IMD generated in a two-tone experiment.
For bigger problems like filters, though this analysis becomes faster than ADS, it still consumes too much time since the computation time increases exponentially with the matrices size. That makes it not feasible for optimizing complex devices in the early design stages considering the nonlinearities.
The IOES method described in the next section improves drastically the computational time.

IV. INPUT-OUTPUT EQUIVALENT SOURCE METHOD
The main objective of this new method is to reduce the size of the admittance matrices, which must be inverted, to the minimum possible. The minimum required nodes are the boundaries of any nonlinear layer. The voltages at these boundaries allow us to mathematically calculate the distribution of the field magnitudes at any point of the stack, which are needed to calculate the nonlinear sources.
As an example, Fig. 4 shows the boundary nodes (9 nodes in this example), marked in black circles, of a one-port SMR BAW resonator if only the AlN and SiO2 layers are considered nonlinear. The corresponding 9x9 admittance matrix at a given frequency i will be denoted as Yr,i. The admittance matrices of the 4-port piezoelectric (AlN block) and each 2-port section of transmission line (W, SiO, etc.), which are used to construct Yr,i, are given in Appendix II. For the sake of clarity, we will henceforth detail each step of the analysis particularizing to the H2 (2·1) generation. The new analysis is performed following the next steps: 1. The non-discretized circuit of Fig. 4 is solved at 1 to get the voltages at the boundary nodes using Yr,1 2. The standing wave field distributions S(z) and E(z) at 1 at any potential nonlinear layer are mathematically calculated instead of solving the whole discretized circuit.

Knowing S(z) and E(z), the distributed nonlinear sources
Fc(z) and Vc(z) at H2 are calculated using (4), (5) and (8). 4. A mathematical process is developed to find equivalent current sources at H2, which would produce the same boundary node voltages than the full distributed model. This process can be mathematically done without inverting big matrices and we call these sources Input-Output Equivalent Sources, as they will be applied at the boundary nodes of a nonlinear guided-wave medium (layers in this case). 5. IOES applied to the same boundary node are added to form the vector of source currents 2 1 and the H2 node voltages are calculated using [ 2 1 ] = [Y r,2 1 ]\[ 2 1 ]. 6. The H2 output power is calculated.
Next subsections describe in further detail the previous steps.

A. Standing wave field distributions
By knowing the boundary node voltages obtained when solving the circuit at the fundamental frequencies, conventional microwave analysis can be used to get the field distributions as a function of the position z inside a given nonlinear layer.
1) Nonpiezoelectric layer The voltage distribution along the layer, or ATL, will be: where V1 and V2 are the voltages at the input (z=z1) and output (z=z2) of a given layer that have been found in step 1 of the analysis, and  is the propagation constant.
V(z) allows to calculate the strain S(z) at each point z of the nonpiezoelectric layer by inspection of the circuit of a unit-cell ( Fig.2a) with z0 and the definition of S (relative displacement of particles S=du/dz ): where VT(z) is the voltage at the junction of the T-network, A is the resonator area and c is the elastic constant. VT(z) can be calculated from the input and output voltage of a unit-cell using 2) Piezoelectric layer For the piezoelectric layer, we know its four boundary node voltages from step 1 of the process. Two of them V1 and V2 correspond to the ATL, and the other two (V3 and V4) correspond to the electrical terminals.
According to the Mason model of Fig. 1, the electrical displacement D at a given frequency remains constant along the discretized model since Iin= -j··D·A and the voltage drop at the acoustic side of the transformer (see Fig.1) will be: where TR is the transformer ratio of the Mason model (TR =-A·e/z). Therefore, if is constant along the stack, (10) can be easily adapted for a piezoelectric layer to get V(z): Now, the remaining step is to find Iin as a function of the known boundary voltages, which is straightforward using the Y-matrix of the piezoelectric layer (see Appendix II): = 31 1 + 32 2 + 33 3 + 34 4 . (15) Using (13)-(15), we will obtain V(z) along the ATL of the piezoelectric layer, (12) gives VT(z) and (11) the strain S(z).
It can be easily demonstrated that the electric field E(z) can be found from the voltage drop between the electrical ports using Iin and VT(z) Once the field distributions S(z) and E(z) are known at the fundamental frequencies, we describe how to get the IOES in the next subsection.

B. IO equivalent sources.
The fourth step of the six-step process previously described is the core of this new method of analysis. The objective is to find the equivalent nonlinear current sources to be applied to the reduced Y-matrix and in turn to speed up the simulations without simplifications and with no penalty on the computational time. The key point is that these equivalent nonlinear (NL) sources cannot lose the distributed nature of the problem and must rigorously account for the local generation of H/IMD at any point of a nonlinear layer. This can be done following the next procedure, which is based on the conventional ABCD matrix of a nonpiezoelectric layer, and a new 4-port ABCD matrix description of the Mason's unit-cell including the nonlinear sources.

1) IOES of a Nonpiezoelectric Layer
Following with the example of finding the H2 of a shortcircuited resonator (Fig. 4), the input signal at the fundamental frequency 1 is applied at the input electrical port of the ALN layer (node 3 in Fig. 5) and the linear circuit is solved providing the voltages at the nine nodes of the circuit.
For simplicity we detail first the procedure of finding the IOES for the non-piezoelectric layers. Let us illustrate the procedure with the SiO2 layer, discretized in N unit-cells, just beneath the bottom electrode as indicated in Fig  Knowing the voltages and currents at the input and output nodes (4 and 5 in Fig. 5), the discrete voltages V1,V2…VN+1 (we denote V(z) for simplicity) along this SiO2 layer are calculated using (10). Then (11), (12) are used to calculate the strain at any point 1 ( ). This is the independent variable of the distributed NL sources at H2, ,2 1 (z), which can be calculated with the counterpart of (8) for a nonpiezoelectric layer: The ABCD matrix of a nonpiezoelectric unit-cell of length z as the one in Fig. 6 is where = 0 ℎ ( ) and = 0 ℎ( · ) . We can rewrite (18) as where If we cascade two identical unit-cells (differing only in the NL source Fc) as it is illustrated in Fig. 7 and we can rewrite last equation as where It is clear then that the input current-voltage pair is related to the output current-voltage through the linear ABCD matrix of the whole section, plus two additional sources, that in this case correspond to a voltage and a current sources applied at the input port of the whole section, as it is illustrated in Fig. 8a. Note that in this figure zs and zp correspond to the whole layer (ATL of length d=N·z).
Appendix III shows some tips to consider for fast calculations of the equivalent sources of (23) and note that, although these equivalent sources are finally applied to a discrete point of the circuit, the distributed effect of the NL sources is mathematically guaranteed with (23). However, as we work with the admittance matrix, we are interested to convert those equivalent sources to only current equivalent sources. It can be demonstrated that the circuit of Fig.8a is completely equivalent to the circuit in Fig.8b   The currents Ie1 and Ie2 are the IOES of this layer and they can be applied to the reduced admittance matrix Y r,2 1 to solve the circuit at the desired frequency as illustrated in Fig. 9. The only thing we must do is, to place the IOES into the appropriate rows of the current source vector and, if no other NL layers were considered, solve [ 2 1 ] = [Y r,2 1 ]\[ 2 1 ] according to the step 5 of the process. The voltage node 3 (see Fig.9) will then be used to calculate the dissipated H2 power.
Note that this method does not imply the inversion of big matrices since the number of nodes is extremely reduced (nine in this case) in comparison with the whole distributed model of Section III.C.
For the other two layers of SiO2 between nodes 6-7 and 8-9 of Fig 5, the process will be the same. Note that the overall size of the matrices to be inverted will not increase.

1) IOES of a Piezoelectric Layer
In this case, the V and I at the four nodes (three nodes for a 1-port resonator) of the ALN layer (See Fig. 4), which have been obtained from a linear analysis at the fundamental frequencies, allows to mathematically find the voltage and current magnitudes at any point of the Mason model, either in its acoustic part or in its electrical part. That will allow to find the independent variables S(z) and E(z) required to calculate the NL sources (4) at the targeted frequencies (8), (9) or the corresponding equations for other H/IMD. If we intend to apply the same procedure as it was done for the simpler non-piezoelectric layers, a 4-port ABCD matrix (we will denote 4-ABCD) must be defined for the Mason model of Fig. 10. This concept resembles to the P-matrix description but relating total magnitudes instead of traveling waves. As far as we know the 4-port description of an ABCD matrix for electroacoustic devices, had not yet previously reported, so no external references are included.
The 4-ABCD matrix of a piezoelectric unit-cell of length z as the one in Fig. 10 is where [4 ] is the lineal 4-ABCD matrix of a differential z section, and the right terms are added terms that affect the  Fig. 11. Cascading 4-ABCD matrices of a piezoelectric layer.
As it was done before, we can cascade many unit cells as it is illustrated in Fig.11 In this case the equivalent sources are as shown in Fig. 12a, and again, for practical issues and use only the Y-matrix, we are interested in finding equivalent current sources as it is shown in Fig. 12b  This can be done using the equation  Fig.13, and we just need to place them at the appropriate rows of the 2 1 vector. Note that in the example of one-port resonator Ieq4 is shortcircuited.
Once all the IOES of any nonlinear layer are calculated we will solve for [ 2 1 ] = [Y r,2 1 ]\[ 2 1 ] to get the node voltages and the power of H2 dissipated at the load.  Fig.13 IOES of a piezoelectric layer are applied to the resonator.

C. Remix effects
In BAW resonators, second-order nonlinearities are more significant than third-order nonlinearities. It has been previously reported that the contribution of remix effects to the generation of IMD3 and H3 may not be negligible in comparison with the direct generation [10].
Remix effects generate third order H/IMD from the second order nonlinear equations that directly generate 2 nd order H/IMD. That means that for example, the generated H2 at 2·1 can be again remixed with 1 to generate H3 at 3·1. In the case of IMD3, several combinations of 2 nd order IMD with the fundamentals can contribute to a given 3 rd order IMD. For example: 2·1 can be remixed with 2 to provide 2·1-2 but also 2-1 may remix with 1. All the involved frequencies must be considered into the equations and thus, the corresponding standing wave pattern must be calculated to obtain the nonlinear sources at a given IMD/H. This fact, slows down the simulations and is not always required. However, we have also considered remixing effects for completeness of the article and to make later fare comparisons with conventional HB analysis of the whole distributed model. In this case, the sixth-step process described in Section IV must be extended.
For example, if we want to calculate the contribution to the H3 due to remix effects, we need to calculate the field distribution of the H2 at any point of the stack once we have completed step 1 to 5 for the H2. This can be done using the ABCD and 4-ABCD matrices of a unit-cell of nonpiezoelectric and piezoelectric layers, respectively. See Appendix I for more details. Once the IOES of H3 due to remix effects are calculated, these current sources must be just added to the IOES of H3 due to direct effects (step 5) and execute the step 6 for the H3.

V. EXAMPLES
This section outlines two examples and compare the computational time required by a commercial circuit simulator using HB techniques and the proposed IOES method.

A. Two-tone Experiment in a single resonator
The first example simulates a two-tone experiment that allows us to find out H's and IMD's generated in short circuited SMR BAW resonator. The resonator was designed to operate at the LTE B30 frequency band and the SiO2 layers of the acoustic reflector were designed for proper confinement of the acoustic wave and for compensating changes in temperature. To take care of the temperature fluctuations, thicker SiO2 layers are required, whose contribution to the overall generation of H's and IMD are quite significant [10]. Details of this resonator can be found in [10].
We have simulated a two-tone experiment with two 10 dBm input sources whose frequencies are separated 10 MHz, directly connected to the input electrical port of the resonator. The central frequency of both signals is swept from 2.2 GHz to 2.5 GHz with 1 MHz steps resulting in 301 frequency points.
The circuit was simulated with HB limiting the maximum order of each fundamental frequency (and maximum mixing order) to 3 and is further compared with the IOES method. For a fare comparison between both methods, all the equations for H and IMD that were simulated by HB: 1, 2, 21, 22, 2-1, 1+2, 31, 32, 21-2, 22-1, 21+2, and 22+1, were included in the IOES analysis and remixing effects were also considered. Fig. 14 shows a comparison between the HB simulations and the proposed method. As it can be seen, HB (dashed) and IOES (continuous) traces perfectly overlap. The computational time for the HB simulator was 366 s, the method described in [13] requires 210 s, and the IOES spent 1.4 s, which means that the IOES analysis was more than 250 times faster than the HB simulation and 150 times faster than [13]. Note that the IOES method is valid for any separation of the frequency of the two tones as the conventional HB simulations of the full distributed model are.
As an example of practical application, Fig. 15 and Fig. 16, compares IOES simulations with real measurements of a B30 resonator that were already presented and deeply discussed in [10]. The effects of the measurement system [20], characterized by a 4-port network, whose S-parameter were measured and conveniently stored, can be easily integrated in the IOES procedure using its Y matrix since the IOES method is based on the admittance description of the whole circuit. There is no penalty in terms of the computational time for including the measurement system and the corresponding interpolation process at all the frequencies appearing in the IOES method. As it can be seen in Fig. 15 and Fig. 16, the smooth ripple appearing in the experimental data due to the measurement system is very well modeled into the simulations.  The measured harmonics and IMD3, corresponding to the experiment outlined above, exhibit several peaks along the frequency range. Those peaks are due to the nonlinearities of the SiO2 layers [13]. The residual H3 and IMD3 of the measurement system by itself was around -90 dBm and -80 dBm for the H3 and IMD3 respectively.

B. Filter
The second example corresponds to a simulated 7-order ladder filter, only designed to illustrate the IOES methodology proposed in this paper.
The filter configuration is outlined in Fig. 17 and it is created by means of B30 resonators as the one in previous example. As mentioned above, the Y-matrix description of the problem makes easier the interconnection of resonators.
We performed a two-tone (input power of 10 dBm) test with the same number of analyzed H's/IMD's that in previous example, both for the HB and the IOES simulations.  Fig. 18 compares the simulated output power at the fundamental frequencies, IMD3, H2 and H3 using both methods. As it can be seen, both methods provide identical results. The HB simulation took 3083s long to simulate 101 frequency points, [13] took 2604 s, whereas IOES simulation took only 4.8 s. That means that IOES method is around 700 times faster than the commercial simulator. This is a very significant improvement in the speed of simulations, which allows us to simulate broadband scenarios. For example, Fig.  19 shows a simulation in a span of 2 GHz with 1001 frequency points that was performed in 49 s. As it can be seen, there are a couple of hot-frequency bands around 3 GHz with very high H2, H3 and IMD3. This kind of broadband simulations could anticipate possible potential interfering harmful scenarios in a real environment.

VI. CONCLUSIONS
The new IOES proposed method has demonstrated to be very useful for the analysis of big distributed nonlinear circuits as it is the case of SMR BAW resonator and filters. In the first case, the computational time is reduced 200 times in comparison with conventional analysis based on HB techniques or the method proposed in [13]. For the analysis of 7 th order filter the reduction is 700 times. This improvement will be more important for even larger nonlinear problems, such as higher order filters or multiplexers because the IOES method does not involve the inversion of large matrices, as occurs with HB techniques when solving the whole discretized circuit. The larger the nonlinear problem, the better the advantages of the IOES method. Fig. 19. Broadband IOES simulations of the output power in a 7 th order filter at the fundamental (red), H2 (blue), IMD3 (green) and H3 (black) frequencies.
The method uses an admittance matrix description of the circuit, which makes easy the interconnection of components or external elements, and it is mathematically rigorous without losing the distributed effect of the problem, although it is valid only for weak nonlinearities. Note that this important limitation does not apply to the more generalized HB techniques. The main drawback of the method is that the nonlinear equations are written in the frequency domain for each H/IMD to be analyzed. This implies a loss of flexibility and for example in a 3-tone experiment, additional equations should be written for each potentially harmful H or IMD. This drawback might also be an advantage if the user is only interested in a given H or IMD because the procedure could be restricted to perform only the needed calculations for the targeted frequency, reducing even more the computational time.

A. Standing wave patterns at a targeted frequency for remix purposes (No-piezoelectric layers)
To illustrate the required steps to evaluate the remix effects, let us assume we want to calculate the H3 generated by remix effects. Then, we will need to find the standing wave pattern at the frequency H2, which cannot be done as it was done before for the fundamental signals since the nonlinear sources at H2 are distributed along the layer and therefore these sources affect to the voltage and current distribution at H2.
The procedure we have followed starts solving the circuit for H2 as it has been described in Section IV. Then we get the input and output voltages V1,2ω1 and V2,2ω1 at the beginning and at the end of the transmission line and be using the Y parameters (Appendix II) of the total layer, we can obtain the current I1 at the input of the layer as (see Fig 9) where Ie1 is the input equivalent source at H2 that was previously calculated. Once the input I-V pair is known, we can apply (20) to get the pair I-V at the output of the first unit-cell since we already know the values of Fc1 of H2. Equation (20) is then sequentially repeated for the N unit-cells to get the distribution V(z) at the H2 that allows us to find S(z).

B. Standing wave patterns at a targeted frequency for remix purposes (piezoelectric layers)
The procedure is like the one described for the nonpiezoelectric layer, but in this case, we will use the 4ABCD matrix of the piezoelectric layer. Once the circuit is solved for the H2, we know the voltages V1, V3 (right-side voltages of the 4-port network) and the right-side currents I1, I3 are calculated with the Y matrix (Appendix II) of the equivalent Mason model (Fig. 13) of the whole layer at the frequency H2. 1 = 11 · 1 + 12 · 2 + 13 · 3 + 14 · 4 − 1 3 = 31 · 1 + 32 · 2 + 33 · 3 + 34 · 4 − 3 . (32) Now we can apply (27) to find the next vector V2, I2, V4, I4 at each unit-cell, and repeat this procedure along the layer.
Finally, we obtain S(z) and E(z) as A similar procedure is followed to find the IMD3 due to remix effects. The main difference is that the remix effect of the IMD3 2f1-f2 for example might come from the mix of H2 at 2f1 with f2, and from the mix of the IMD2 f2-f1 with f1. Consequently, the magnitude distributions at 2f1 and at f2-f1 must be calculated.

C.
Admittance matrices The admittance matrix of a section of transmission line like the one in Fig A2. a, is whereas the admittance matrix of the 4-port Mason model ( Fig. A2. b) is , where the parameters A, B, C, and D are calculated accordingly to their subscript that indicates the length i·z, where i=1…N. A similar matrix operation is done for the piezoelectric layer.