PRECISE Photonic hybRid EleCtromagnetIc SolvEr

The Photonic hybRid EleCtromagnetic SolvEr (PRECISE) is a Matlab based library to model large and complex photonics integrated circuits. Each circuit is modularly described in terms of waveguide segments connected through multiport nodes. Linear, nonlinear, and dynamical phenomena are simulated by solving the system of differential equations describing the effect to be considered. By exploiting the steady state approximation of the electromagnetic field within each node device, the library can handle large and complex circuits even on desktop PC. We show that the steady state assumption is fulfilled in a broad number of applications and we compare its accuracy with analytical model (coupled mode theory) and experimental results. PRECISE is highly modular and easily extensible to handle equations different from those already implemented and is, thus, a flexible tool to model the increasingly complex photonic circuits.


I. INTRODUCTION
T HE continuous evolution of silicon photonic produces Photonic Integrated Circuits (PIC) of endlessly increasing complexity.To model the optical response of large PICs requires intensive computational effort if approached by directly solving Maxwell equations with either FDTD or FEM methods.Several numerical tools are already available to accurately model basic building blocks and simple PICs (e.g., MEEP [1], FEM as Lumerical [2] and COMSOL [3], etc.).Other opensource tools like PhotonTorch [4] and Symphony [5] permit the analysis of linear systems only but lack the possibility to handle Non-Linear (NL) effects.To our best knowledge, only CAPHE [6] is able to handle nonlinear photonics devices.Yet the analysis of complex circuits would greatly benefit from proper coarse-grain modeling, considering only the most representative aspects of PIC nodes and simplifying the circuit description to reduce the overall computational effort.A significant simplification of the PIC descriptions will permit to increase the size of the circuit to be modeled as well as to include nonlinearities that are often at the heart of many photonic functions and nonlinearities are even increasing their interest in the emerging field of neuromorphic photonics [7]- [11].The description of NL optical phenomena requires solving systems of differential equations (ODEs) that are inherently serial problems and cannot efficiently exploit parallelization techniques (albeit recent libraries propose improvements [12]).
Here we propose the Photonic hybRid EleCtromagnetIc SolvEr (PRECISE) Toolkit to ease the modeling of linear and nonlinear spectra and the dynamics of complex PICs by assuming steady-state conditions of the electromagnetic fields within each fundamental circuit block.We will show that the steady-state approximation is fulfilled in a large number of devices designs and for state-of-the-art applications and we confirm the excellent accuracy of our model by comparing its results with different analytical schemes and experimental results.The modular structure of PRECISE does not rely on a specific assumption about the PIC building blocks and it handles any type of passive devices found in integrated photonics.This fact enables the study and design of both known and novel photonic structures and to sweep over broad space parameters.

II. THEORETICAL BACKGROUND AND PICS DESCRIPTION
PRECISE describes PICs as waveguides (WGs) connected by nodes and interacting with the surrounding via I/O ports.The idea of our code is to evaluate complex PICs without resorting to the direct solution of ODE systems at the sub ps time-scale.For purely linear PICs, the system is described in terms of symbolic equations that are inverted and solved to calculate the parameters of interest.In the case of nonlinear PICs, PRECISE assumes that the field intensity instantaneously achieves its steady-state, within each node.Such approximation (similar to [13]) is satisfied whenever the optical field stationary state is reached on timescales smaller than the characteristic dynamics of the physical processes defining the nonlinear behavior (for example, the free carriers or the material's temperature).As we will show, this approximation is often fulfilled and greatly simplifies the solution of ODE systems.

A. Description of PIC
PICs are described programmatically by decomposing them into basic components: any PIC is essentially formed of WGs segments that are connected together at nodes.To completely describe a PIC we need to define three fundamentals elements: 1) Nodes are logical blocks that groups a number of ports together.A node can have as many input/output ports as necessary and its generic connectivity is described by the following matrix equation: where E in and E out are 1 × n column vectors and C is an n × n matrix.In its simplest form coupling only includes constant coupling parameters, but a frequency dependent coupling can also be included, thanks to the symbolic definition of the nodes themselves.For example, a 2 × 2 directional coupler is represented by a node with four ports (shown in Figure 1a).Its coupling matrix is The user can add a node add_node and define the coupling coefficients of its ports using add_link, as it will be explained better later.2) WG segments connect the output port (i) of a node to an input port (j) of the same or different node, as described in Figure 1b.We consider the solutions of the wave equation as a monochromatic plane wave propagating along the WG (z): whose transverse amplitude profile u(x, y) and propagation constant β belong to one of the guided modes.
In principle a WG may support many modes, however, they are often designed to support only the fundamental mode and our toolkit restricts to simulations of this case.
The nonlinearities of PICs are introduced with a small variation of the effective index in In particular, when studying silicon-on-insulator PICs, the effects usually considered are the third-order nonlinear optical effects, the free carrier effects, and the thermo-optical effect, and so n ef f,N L is a function of the field amplitude |A| 2 (optical intensity), the free carrier density N , and the temperature of the material T .
The fields propagate symmetrically through the WGs following E(t) is the complex amplitude of the electromagnetic (EM) field, ν is its optical frequency, c 0 is the speed of light in vacuum, L is the WG length, and n ef f and n ef f,N L are, respectively, the linear and nonlinear part of the effective index of the mode propagating within the WG.The user can define a new WG between existing ports using add_guide.
3) The last important piece is the definition of one or more input ports to the whole system.In order to designate an input port, the user should call the function add_input with the port index as an argument (see later).The input field of ports that are not connected with WGs or defined as inputs is assumed to be zero.
When these three elements are correctly defined, PRECISE automatically creates the system of parametric equations representing the structure.For example, by combining two 2 × 2 couplers and 2 WG segments, a microring (MR) resonator in the add-drop filter configuration is obtained, as depicted in Figure 1c.The set of equations derived from this optical system reads: where E in/out j are the complex field amplitudes, k 1 and k 2 are the coupling constants, and the phases φ 1 and φ 2 are defined in Equation (5).Such description of electromagnetic fields is of general validity and allows both counter-propagating fields as well as refractive index variations due to nonlinear effects.The direct solution of this symbolic system of equations is inefficient and becomes quickly intractable with the increasing number of PIC elements.Thus, the system of symbolic equations is transformed in its matrix form and the numerical values of all parameters are substituted before its inversion.To do this, PRECISE has a built-in function that takes as input the parameter's values and returns the numerical matrices describing the system.These matrices are then fed to the linear system solver algorithm, which solves them by inversion and finds the complex amplitude of the fields at the output port of each and every connection node.For example, in the case of the add-drop filter configuration with one input port (see Figure 1c), the resulting function becomes: and its output is where ν is the optical frequency, in 1 is the complex field amplitude of the input signal, k are the (two) coupling constants of the coupling regions, L are the (two) lengths of each WG segment, the n ef f and n ef f,N L are the linear and nonlinear effective refractive indices of the (two) WG respectively, and E out 1 , . . ., E out 8 are the complex amplitudes of EM fields at the output ports.The resulting spectrum for this case, i.e. a single MR resonator in the add-drop filter configuration, is shown in Figure 2.Each line is obtained taking the square modulus of the given field and then is normalized for the square modulus of the input field, eventually obtaining the transmission: 191 191.5 192 192.5 193 193.5 194 194 The general procedure to create an optical system (all ports are indexed with positive integers) is: % initialize the class object S = optical_system( ); % add nodes and connect them with links % and waveguides ; n is the number of ports .

S.add_node( n );
% add direct links between the ports , such as those % found in a 2 x2 coupler .It follows the equation % ao =K * bi , with a ( output ) and b ( input ) port % vectors and K a coefficient matrix .

S.add_link( a, b, K );
% connects the output port p1 to the input port p2 % and the output port p2 to the input port p1 % with a waveguide .

S.add_guide( p1, p2 );
% define the port p0 as inputs of the structure S.add_input( p0 ); % generate the function that evaluates the fields f = S.generate_pfe( ); The last command, generate_pfe, generates the function that parametrically evaluates the output fields, as reported in Equation ( 7).This function is called parametric field evaluation (PFE) and is generated automatically by the toolkit.To simplify the modeling of PICs further, our toolkit enables a modular definition of new structures.Hence, for example, a MR resonator structure in the add-drop filter configuration is easily defined with:

B. Temporal solution and nonlinear effects
In principle, the temporal evolution of the EM field in a PIC is the solution of a system of differential equations describing the relevant parameters: the optical field complex amplitude at each port, the temperature and the free carrier concentration in each WG, etc.Such a method is computationally intensive as it requires taking time steps small enough to track the complex amplitude dynamic (∼ ps).On the other hand, if we assume that the optical field reaches the stationary state at each time step before any relevant change happens to either temperature or free carriers concentration, then we can consider an instantaneous transmission of the input optical field throughout the system, which is easily evaluated with the linear description obtained in the previous section.To correctly describe the system in its steady-state, the n ef f,N L is evaluated at each time step with the current values of temperature and free carriers.
The ODE solver algorithm used is MATLAB's own 'ode23' algorithm, a non-stiff differential equations solver implementing an adaptive step size.The different steps that the ODE solver carries out at each iteration are: 1) evaluation of refractive index variation n N L 0 due to nonlinear effects involving temperature and free carriers [15]- [17]: The effective refractive index change due to nonlinear effects is then obtained using the relation where Γ is the optical mode confinement factor, defined as in [17].2) evaluation of the optical fields using the parametric field evaluation function as seen in Equation (7).Here the input signal at time t becomes the input field amplitude in.
3) evaluation of the differential change for the free carrier population and the temperature in each WG, adapted from [18]: I = 1 2 c 0 0 E 2 is the irradiance of the field circulating in the WG, where 0 is the vacuum permittivity, and U = Ln g E 2 /c 0 its energy, evaluated from the fields E, while ∆T and ∆N are the deviations of the temperature and the Free Carrier (FC) density from their equilibrium value, respectively.Moreover, the power generating heat is (adapted from [18]): The other parameters are described in Table I.
Notice that if the first two terms in Equation ( 9) are negligible, points 1 and 2 are evaluated only once.Otherwise, their evaluation should continue until self the consistency condition is found.In this case, however, the computational complexity becomes similar to that of the CMT and so does the simulation time.In essence, the dynamic of effects faster than roughly ten times the cavity lifetime (see the more in-depth discussion in the Section IV) cannot be accurately described by PRECISE (i.e. as the Kerr effect and the Two Photon Absorption (TPA) contributions) since they are neglected during the updates of the optical fields in Equation ( 9).Yet, for silicon photonics, these effects are often orders of magnitude smaller than free carriers and thermal dispersion ones.However, TPA is still considered in Equation ( 13) for the evaluation of the total absorbed power.

C. Parametric solutions
An important goal of this toolkit is to allow the analysis of systems by exploring their parameter space, which may have several dimensions even for small systems (i.e. for a single MR we might consider: MR radius's, wavelength, losses, coupling strengths, etc.).For this reason, in addition to the fact that the functions evaluating the linear field propagation The values whose sources are marked with an asterisk * should be considered of the same order of magnitude of the given sources.The left half of the table contains values typical of silicon photonics structure, whereas the right half contains quantities whose values have been adapted to describe the structure studied in Section III-B.
and its nonlinear evolution are parametric, the toolkit helps the user to explore the parameter space methodically.A specific class (dataloader) stores all the values that any parameter assumes, such that any desired combination can be retrieved via the indexing of a single integer number, going from 1 to the total number of combinations N C .This set of parameters is then given as arguments to the function representing the system that, in turn, is evaluated by the ODE solver.The benefits of this technique are twofold.The first advantage is that it is straightforward to run parallel evaluations of many configurations, using Matlab's parfor loop, thus increasing the overall computational efficiency.The second one is that each temporal solution is analyzed right after its evaluation and, once the interesting figures-of-merit are extracted, the data of those time series can be eliminated.A given solution can be re-generated at any time by simply looking at its scalar index.Section II-D will explain this more in detail.These three characteristics, together, contribute to a reduced memory footprint and a faster overall execution time, particularly for large parameter spaces.

D. Data analysis and visualization
For each combination of parameters, the time evolution of several quantities is evaluated.PRECISE calculates the complex field amplitude at each port and the deviation of the temperature and the FC density from their equilibrium values.For example, in the case of a MR, the toolkit evaluates the evolution of twelve time series: 8 fields, 2 temperatures, and 2 free carriers, since each MR is described by two WG segments.The output dataset would then be three-dimensional, with size N C ×N S ×N t , where N C is the number of combinations, N S is the number of time series per combination, and N t is the number of samples in each time series.When long traces are evaluated or when the N C is large, the amount of data that is generated becomes quickly intractable.
Since in almost any case N t > N C > N S , it is preferable to synthesize the results by extracting a (scalar) figure-ofmerit from each time series.The resulting dataset is bidimensional and has size N C × N S .Apart from the most straightforward (i.e.minimum, maximum, and average values, from which the signal visibility is derived), custom functions are defined to evaluate more complex figures of merit, such as the period of the most significant FFT com-ponent (periodFFT), the number of the extremal points (extremals), and the unique values of the extremal points (non-scalar, extremals_counts).
Even though the figure of merit condenses the temporal dimension into a scalar value, the dimension along which the results of different combinations are stacked represents the parameter space explored by all the combinations.This is usually n-dimensional and hence requires ad-hoc visualizations.Our toolkit provides several functions and classes to properly visualize this information.A first class (dashboard) offers a way to observe the temperature and the FC density average, maximum, and minimum values for all the combinations evaluated.It also plots their time evolution as well as the absolute value of the complex field amplitude for a specific combination of parameters selected by the user.A separate class (array_plot) adds the possibility to observe any generic figure of merit for all the combinations related to two parameters at a time, e.g., the frequency against the input power or the coupling coefficient against the MR's radius.All 'views' can be updated with different parameters and/or for different combinations.

III. USE CASES
In this section, we present several use cases to demonstrate the capabilities of our toolkit.The first two examples rely on linear spectra calculation of devices with peculiar structures.The third case describes the temporal dynamics of a nonlinear resonator.

A. Linear spectrum
The simplest way of using PRECISE is to simulate linear transmission spectra of integrated optical structures.The outcome of this simulation is a set of complex field amplitudes, from which the transmission and enhancement factor values are easily extracted.We propose two examples of resonant devices with peculiar structures and compare the modeling results with experimental data.
1) Coupled Resonators Optical Waveguide structure: A CROW is composed of side-coupled multiple MRs, sandwiched between two bus WGs in "add-drop" configuration [20].We show the result of a simulation that considers a CROW made of four rings with the input signal entering from both "input" and "add" ports with a π/2 phase difference.This particular configuration is known as "interferometric band interleaver" [21], [22] and its scheme is reported in Figure 1e.
In Figure 3 one can observe that the model correctly predicts the free spectral range of the resonances and their complex shape due to the interaction between the rings.We coupled the PRECISE solver with an optimization routine base on Matlab's 'particleswarm' to find the optimal set of parameters describing the CROW spectral response, considering all radii and all coupling coefficients equal.The outcome of the optimization correctly describes the free spectral range of the resonance as well as their complex shape due to the interaction between the rings, as reported in [14], [22].The parameters are reported in Table II  2) Taiji microring resonator structure: This example describes a Taiji MR resonator composed of bus WGs and reflective facets.A complete in-deep description of this device can be found in [23].Briefly, the Taiji is a particular MR configuration described by a node with two ports designed such that the reflection of either port depends on which port the input field enters the device itself (Figure 1d).The device is reciprocal as the total transmission of the structure is the same for both propagation directions, but it shows an asymmetric reflection response.Our model correctly describes all the Taiji spectral features: the asymmetric reflection, the symmetric transmission, as well as the Fabry-Perot effect of the facets.Figure 4 shows the back-reflection (Rx) and the transmission (Tx) for the backward (bwd) mode and the transmission only (Tx) for the forward mode (f wd), for both the numerical model and the experimental data.In the Taiji the total reflection is due to counter-propagating fields, thus the good agreement with the experimental data confirms that those fields are correctly handled by the toolkit.Similar to the CROW case, we employed an optimization routine to find the set of parameters that best approximate the spectral response.Table III reports the simulation parameters

B. Temporal dynamic of nonlinear microring resonator
Considering a single MR, for high enough fields intensity, we expect to observe a region around the MR's resonance exhibiting the so-called "self pulsing" behavior [13].The "self pulsing" arises as a dynamical modulation of the optical field at the output of the MR as a response to a continuous wave input signal.The specific shape and period of the output signal depend non trivially on several parameters, e.g., device materials and geometry, injected optical power, or detuning from the resonance.For the case reported in Figure 5, we consider a single ring resonator, whose radius is 7 µm and whose coupling coefficients are k 1 = k 2 = 0.17, probed with a signal of constant optical power.The values of all the other parameters are the same reported in Table I.The red line of Figure 5 shows the NL response of the "drop" port (||E out 6 || 2 ) for three different combinations of input power and frequency detuning.In the same graphs, the temporal evolution calculated with a CMT model for the same input parameters is shown with the blue line.Our code correctly describes both the frequency and the amplitudes of the signal modulation, with only small differences compared to the CMT model.These are due to the different approximations considered by CMT and by our code.One of the most important is that while the resonance frequency emerges automatically in our code, it is imposed as a parameter in the CMT model, thus resulting in a possible misalignment of the frequencies evaluated in PRECISE and CMT.Another important difference is that in PRECISE the microring cavity is composed of separated waveguide segments connecting at the coupling regions, while in the CMT it is a unique waveguide.This enables PRECISE to have different values for the physical quantities in each section of the microring.Moreover, the dispersion of n ef f and the waveguide losses are approximated differently in the two codes.Overall the matching of the two codes is good  ).The region where the microring does not self-pulse is shown as having a 0 ps period.The horizontal white lines highlight the microring resonance frequency, at 0GHz detuning.(a) Shows the values for PRECISE, whereas (b) shows the results for the corresponding CMT implementation.The highlighted points refer to the curves shown in Figure 5. and it is also confirmed by the results described in the next paragraph.
Figure 6 shows the stability maps obtained analyzing the complex field amplitude of the "drop" port ( E out 6 2 ) in a MR resonator as a function of the input power and of the frequency detuning.The maps report the principal non-zero Fourier component of the spectra of the MR in the self-pulsing regime: panel (a) refers to our code, while panel (b) reports the results for the CMT one.The three points highlighted in the map refer to the spectra reported in Figure 5.
As expected, the system shows self-pulsing behavior in a limited space region around the resonance frequency and at input power above a certain threshold.The actual shape of this region is heavily dependent on the details of the system, as it is defined by the interplay between the optical energy in the cavity, the free carriers population, and the temperature of the material.The shape of the self-pulsing region is very similar between the two implementations, with only minor differences in the actual period of the oscillations.

IV. BENCHMARKS
We perform two different benchmarks to test both the accuracy and the speed of the library.
The accuracy of the model has been tested by comparing the propagation of optical pulses of different widths sent through a single MR in add-drop configuration near its resonance frequency (see Figure 7).The four panels describe the variation of some of the fundamental system parameters: (top-left) the peak power in the T-port, (top-right) peak power in the Dport, (bottom left) MR temperature increase, and (bottomright) free carrier concentration increase.For all simulations the maximum peak intensity has been kept constant, thus the overall total pulse power increases proportionally to its duration.The main difference between PRECISE and the CMT is visible in the dynamics of the MR loading, where the CMT correctly describes the dynamics of the energy pile up in the cavity, while our model assumes a steady-state condition.On the other hand, both ∆T and ∆n are correctly described.The time scale that PRECISE correctly handles can be inferred by using the cavity lifetime as a scaling factor.The vertical black lines indicate the cavity lifetime and it is clear that PRECISE accurately describes the system dynamics for timescale one order of magnitude greater than the cavity lifetime.For dynamics with characteristic times slower than τ , the difference between the two models is around 2%.
Finally, we run a benchmark to show the performance of PRECISE by running a suite of simulations on CROW structures of increasing complexity.We choose the CROW as it is composed of multiple rings and so by a number of elements (WGs and ports) proportional to the number of the rings.All the simulations have been run with Matlab 2020a on a laptop PC equipped with an AMD Ryzen 5 2500U (15W, 4 cores, 8 threads) and 8 GB RAM.Performance metrics are shown in Figure 8a, where all graphs are plotted in a log-log scale and the values are obtained with a statistic of 10 consecutive runs.Figure 8a shows the time required to initialize the symbolic system and to obtain the parametric function that evaluates the EM fields at each port as a function of the number of elements in the system.The system generation time shows a polynomial trend which sees the increase of time needed by one order of magnitude for an increase in the number of elements by almost two orders of magnitude.Figure 8b shows the time required for a single evaluation of the parametric function, where the actual value has been obtained from a batch evaluation of 1001 points.In this case, the linear evaluation time shows initially a polynomial trend which then seems to approach a power-law trend for higher number of system elements.Finally, Figure 8c shows the scaling of the evaluation time on the number of processes, with three lines representing the parallel evaluation of the system in 51, 101, and 201 frequency points, respectively.With just four workers we see a reduction of the computation time of 38% for 51 points, 47% for 101 points, and 51% for 201 points.In a real case scenario, the user would probably simulate several hundreds or thousands of combinations.In this case, we can expect the performance improvement obtained with the parallel evaluation to be much more beneficial than in the small benchmarks presented here.
V. CONCLUSION PRECISE is a highly modular toolkit that allows simulating time-resolved linear and nonlinear regimes of large PICs.The software is based on a set of general differential equations that describe any passive PICs with the only approximation of steady-state condition of the electric field within each block at each time step.The engine is highly expandable to enable the description of different systems of equations (i.e. in [24] a double decay time has been introduced to correctly describe the free carriers' decay dynamics).Our results demonstrate that the steady-state approximation is rather general and enables the modeling of large circuits without sacrificing the accuracy of the simulations.Furthermore, the toolkit allows for efficient parallelization of workloads and enables to sweep over broad parameter spaces.The validity of the library is demonstrated by the comparison with another analytical model (CMT) and, more importantly, with experimental results with peculiar PIC devices such as the CROW or the Taiji MR.The dynamic of nonlinear effects is also demonstrated with the analysis of the self-pulsing regime.

APPENDIX FIT PARAMETERS OF THE LINEAR SPECTRA
The following tables contain the fit parameters for the models described in Section III-A1 and Section III-A2.The data can be easily reproduced by feeding the set of parameters to the 'linear_crow' and 'linear_taiji' functions found in 'precise.examples'.

8 Fig. 2 :
Fig. 2: Microring resonator linear response.Transmission (power) associated with the field at the through port (red), and at the drop port (black), and (power) enhancement factors associated with the fields circulating in the right half (blue) and left half (yellow) ring respectively.The transmission correlated to the counter-propagating fields E out 1 , E out 4 , E out 5 , and E out 8 are evaluated too, but are identically zero.

Fig. 3 :
Fig. 3: CROW structure (interferometric band interleaver) linear response.Transmission (power) of the fields at the through port (cyan) and at the drop port (red).The experimental values (black circles) of the transmission in the drop port are shown as well.

Fig. 5 :
Fig.5: Microring nonlinear response in the drop port (||E out 6 || 2 ) for three different combinations of input power and frequency.The curves shown above are obtained from longer simulations and the time interval is chosen so that the response can be considered as a stationary state.

Fig. 6 :
Fig. 6: Instability region colored with the period of the most important Fourier-component observed on the drop port ( E out 6 2).The region where the microring does not self-pulse is shown as having a 0 ps period.The horizontal white lines highlight the microring resonance frequency, at 0GHz detuning.(a) Shows the values for PRECISE, whereas (b) shows the results for the corresponding CMT implementation.The highlighted points refer to the curves shown in Figure5.

Fig. 7 :
Fig. 7: Comparison between PRECISE and CMT models for the propagation of Gaussian pulses of varying width through an add-drop filter near its resonance frequency.(c) MR temperature difference, (d) Free carrier concentration, (a) peak power in the Through port, and (b) peak power in the Drop port.Vertical black lines point to three times (3τ ) the cavity lifetime.Continuous blue lines indicate the results of PRECISE, while dashed green lines are the CMT model (left axis).The red, solid line reports the residuals (right axis).

Fig. 8 :
Fig.8: Performance analysis for the CROW structure, carried out by taking the average and the standard deviation on ten consecutive runs.Each additional ring in the CROW adds two elements to the system.(a) Time necessary to initialize the system and generate the corresponding symbolic equations.(b) Time necessary for a single linear evaluation of the EM fields.Obtained from a batch evaluation of 1001 points.(c) Time necessary to simulate the nonlinear evolution of a CROW with two rings for 51, 101, and 201 frequency points as a function of the number of parallel workers (processes).Matlab's own parpool initialization time is neglected.Time may vary with the task.
Each system is firstly defined by a set of nodes (small squares), which group a number of ports (black dots).The WGs may then connect the ports and are represented by thick gray lines.The dashed arrows illustrate the input/output fields at each port and are identified by integer numbers.Input and output of the whole system are instead depicted as red and blue arrows respectively.(a) Example of a node with four ports, showing both input and output fields.
(b) Example of WG connecting two ports belonging to different nodes.(c) Scheme of a MR resonator built with two 2 × 2 nodes, for a total of eight ports, and two WGs.There are two input fields, input and add, and two output fields, through and drop.(d) Scheme of a Taiji MR.This structure is composed of three 2 × 2 nodes and four 1 × 1 nodes, for a total of 20 ports, and eight WGs.(e) Scheme of a CROW structure with four rings.The structure is composed of five 2 × 2 nodes, for a total of 20 ports, and eight WGs.The input/output fields are the same as the MR, but, in this case, both the input and add fields are used simultaneously.

TABLE I LIST
OF THE PARAMETERS THAT APPEAR IN EQUATION (12) WITH THEIR DESCRIPTION AND DEFINITION. .
L 2πR + 10 µm k are field coupling coefficients and their indexes identify the ports of the node in which they are used.L are the WG lengths and the indexes identify the port at the ends of the WG.All indexes refer to Figure1e.
are field coupling coefficients and their indexes identify the ports of the node in which they are used.L are the WG lengths and the indexes identify the port at the ends of the WG.All indexes refer to Figure1d.