Graph model for multiple scattering in lithium niobate on insulator integrated photonic networks

We present a graph-based model for multiple scattering of light in integrated lithium niobate on insulator (LNOI) networks, which describes an open network of single-mode integrated waveguides with tunable scattering at the network nodes. We first validate the model at small scale with experimental LNOI resonator devices and show consistent agreement between simulated and measured spectral data. Then, the model is used to demonstrate a novel platform for on-chip multiple scattering in large-scale optical networks up to few hundred nodes, with tunable scattering behaviour and tailored disorder. Combining our simple graph-based model with material properties of LNOI, this platform creates new opportunities to control randomness in large optical networks.


Introduction
Graphs, or networks, are commonly used structures for modelling physical systems.The utility of their geometry is twofold: on the one hand, complex dynamics can be modelled by a large network of simple elements and their pairwise connections.Such ideas have been applied to numerous examples ranging from quantum information transmittance [1], improving power grid robustness [2], and community structures in biological and social systems [3].In the context of optics, networks describe optical circuits, which can perform elaborate operations using connections between simple optical elements [4][5][6][7].On the other hand, a graph often serves as an appropriate model of dimension reduction and discretisation of a real physical medium [8], for example for computing the band structure of graphene [9] or carbon nanotubes [10].
The second case above often falls within the theory of quantum graphs [11].Quantum graph models may not necessarily describe systems which are quantum mechanical in nature; for our purposes, they describe waves (where interference comes into play) on metric graphs (where both the topology of connections and geometry of distances play a role).Classical realisations for quantum graphs have been achieved in the microwave domains, using networks of coaxial cables [12,13].For optics, such models have been applied to random lasers: in these examples, the graph edges are either composed of dye doped polymer nanofibres spun on chip [14,15] or optical fibres with coherent optical amplifying sections as gain [16].
In this work, we apply network-based models to linear photonic integrated circuits (PICs), as a platform for studying on-chip multiple scattering with fully controlled and pre-designed disorder.To this end, several modifications are made compared to existing quantum graph models for optical systems.For one, focus is put on the ability to fully specify the scattering behaviour at each scattering site in the network, which allows for greater flexibility and better agreement with integrated optical networks compared to Neumann-type continuity conditions commonly used in the quantum graph literature [11,17] and in previous work on integrated random laser [14,15].Secondly, our model focuses on open networks, i,e. the linear transport properties by looking at the continuous spectra, rather than discrete resonance modes.While closely related, the coupling to open ports modifies the resonances of a network, and the open model matches more closely with the physical measurements.It is also worth noting that this approach differs from traditional optical circuits [4][5][6][7], since backscattering is utilised as a source for generating interference.
Small scale graphs, which consist of a few to several dozens nodes, can be constructed to model traditional integrated optical devices such as ring resonators, and automatise the calculation for circuits composed of several coupled elements.We demonstrate that by taking into account the appropriate material properties, quantitative agreement is reached between simulated and measured spectra for integrated lithium niobate on insulator (LNOI) devices.LNOI is chosen as the physical platform due to its electro-optic (EO) property, allowing tuning of network geometry by applying voltage across waveguides [18,19].EO tuning comes into play when we look at scaled up network models, where complex random networks serve as a model for multiple light scattering [20].In such cases, tuning allows for control over scattering and propagation parameters on individual optical elements.
By implementing a scattering matrix model which is significantly simplified compared to numerical techniques such as finite-difference time domain (FDTD) and finite element method (FEM), we are able to explore integrated photonic networks with few hundred nodes as a new platform for studying multiple scattering.The nature of top-down design and fabrication allows one to control the exact type of randomness introduced into the network, as well as behaviour of individual scatterers.

The scattering matrix model
We start with the description of the scattering matrix model for integrated networks, distinguishing between two cases: 1) the closed network, which is the common model adopted for fibre lasing networks [16] and quantum graphs [11,17], and 2) the open network, which more realistically accounts for the coupling geometry in a measurement setup.

Closed network model and the secular equation
The integrated network is represented by a mathematical graph  = (, ), defined by a set of edges  and a set of vertices  (Fig. 1).For our purposes, the graph is assumed to be connected (i.e.there exists no isolated components) and undirected (i.e. the edges have no inherent orientation).The edges, or links, correspond to single-mode waveguides along which light propagates.The vertices, or nodes, correspond to the scattering elements in the network: they can be sites of reflection and power splitting, for example Bragg gratings and beam splitters.Note that we are restricting our system to 2D, where only in-plane scattering is considered; out-of-plane scattering is incorporated as loss.For this paper, we will use nodes to emphasise the fact that scattering sites are scattering elements with internal structure, and edges to highlight the metric properties beyond connectivity.Light propagation in the network is modelled as a superposition of two counter-propagating waves with amplitudes Φ + and Φ − , which are defined locally on each edge as where  ∈ [0, ] is the local spatial coordinate defined independently on each edge,  being the length of the edge.The set of amplitudes on all edges fully characterise the field distribution in a graph, and are computed using matrix equations as will be shown below.Edges are joined together by nodes, which introduce further phase difference and local interference effects that are encapsulated in the local scattering matrix for each node.A solution, or mode, on a graph is an admissible set of amplitudes, which are consistent with the scattering behaviour and geometry of the network.It is obtained by treating scattering (nodes) and propagation (edges) in conjunction through the matrix equations described below.
To construct the relevant matrices, we define two sets of field amplitudes: , the set of all incoming amplitudes at all the nodes, and , the set of all outgoing amplitudes.Each set is related to Φ + and Φ − by a simple phase factor, illustrated in Fig. 1a.The local scattering matrix at each node is defined by the optical element that it corresponds to; it is a  ×  matrix that maps  to , where the degree  is the number of edges attached to the node.The propagation matrix maps  to  along each edge, and has the form for an edge of length  between nodes  and .The overall scattering matrix  and propagation matrix  for an entire network is constructed by embedding local matrix elements at the correct index.It is clear in the edge-centred picture that either  or  alone is sufficient for describing a solution (i.e. a full set of amplitudes) in the network, with a set of phase differences between them.One small but important difference between these solutions and Φ (as specified in Eq. 1) is the direction: propagating along an edge from  to  always picks up a phase in the same sign    , while Φ + and Φ − pick up    and  −  , respectively, corresponding to a coordinate change  →  −  implicitly performed for the Φ − amplitudes.A solution on the graph is then an eigenvector  with eigenvalue 1 when mapped back onto itself via : (or alternatively, ( − 1) = 0).In analogy to the propagator in the quantum graph context [21], we define  () := (), where we make explicit that the k-dependence comes from the propagation along the edges, while the scattering at the node can be treated as wavelength independent. () is unitary when both  and  are lossless.The modes of the network then correspond to the nontrivial solutions of the secular equation and may only exist for a discrete set of   .

Open network model and the transmission matrix
The distinction between open and closed networks can be expressed as follows: networks with open edges, i.e. edges connected to nodes of degree 1, are classified as open; otherwise they are closed (Fig. 1c-d).The reason for this distinction is that the above characteristic equation formalism poses boundary conditions that are too restrictive for a physical device.Consider an external node with degree 1: the scattering matrix for this node is simply complex scalar, characterising the reflection back into the network from the boundary.In particular, if we want external edges to correspond to physical waveguides with no back-reflection, this scalar would be zero, resulting in a solution with zero field on the edge.The constraint matrix would then have the block matrix form where 1  is the  ×  identity matrix,  being the number of external edges;   and   denote the left and right sections of the remaining 2  −  rows,   being the total number of edges.In this basis, the solution vector is arranged to have the form  = (   ,   ,   )  , where   are the amplitudes on external edges travelling into the network,   are those on external edges coming out of the network, and   are amplitudes on internal edges.It is clear that no matrix of this form would satisfy Eq. 2 with nonzero   .Instead, we modify Eq. 2 to allow for nonzero input:  = ()  + (   , 0, 0)  , and therefore  = −( − 1) −1 (   , 0, 0)  .Using standard block matrix inversion formulas, we can pick out the relevant submatrix in −( − 1) −1 : The matrix (− −1    ) takes the input amplitudes into the entire network, and maps to all other amplitudes constituting a solution.In particular, we define the matrix   as the first  rows of (− −1    ), characterising the input-output relation in the external waveguides of the network via   =     .Under this construction one obtains a continuous spectrum in , the resonances of which correspond to the real part of the discrete eigenvalues of  () for the closed network [22].
One final caveat is that   differs from the true experimental input   by the phase picked up from travelling along external edges.In Fig. 1a, if   is the external node, this is the difference between   and    .We therefore redefine   →     , where   is the propagation submatrix applied to the external edges (if one starts by solving for , a similar correction is needed to go from   to   ).This ensures the final   matrix is reciprocal in both amplitude and phase.

Experimental verification on integrated lithium niobate resonators
To verify the model's capability to describe real integrated photonic networks, simulation results are compared with measurement data on LNOI components fabricated on chip.We demonstrate agreement between the simulations and measurements when material characteristics are accounted for in the model.In particular, we look at two classes of resonators: ring resonators and interferometric resonators.They correspond to two simple types of networks in our formalism.Indeed, a complex network can be interpreted as a complex set of coupled resonators, and these simple resonators can be viewed as the building blocks of larger systems.
The devices were fabricated on a 600 nm X-cut LNOI chip patterned using electron-beam lithography and argon ion milling.To expose the end-facets of the waveguides, the chip was diced.Polarisation maintaining lensed fibers were used to couple to the fundamental transverse electric (TE) mode of the waveguide.The spectral measurement on ring resonators use a tunable laser source and a high dynamic range power meter, while that on interferometric resonators use a broadband superluminescent diode source and an optical spectrum analyser.Details on the fabrication and characterisation techniques are presented in previous work [23].

Ring resonators
A single ring resonator, coupled to one bus waveguide via evanescent coupling, is a notch-type filter commonly found in integrated photonic devices [24].Figure 2a shows one such device, with 2b being the corresponding network.Note the ring can be modelled with one curved edge, but is here shown as a square, demonstrating that a desired balance between model size (number of edges) and abstraction (all straight edges, ease of plotting) can be achieved by properly choosing the scattering matrices.The top center node, corresponding to the evanescent coupler, .We assume | | 2 + || 2 = 1 for a lossless scattering node.
One can then follow the procedure described in Section 2.2, and recover the symmetric transmission between the two open ports (up to an additional phase picked up by traveling along the external bus waveguide) to be  16 = (1 −  − )/( −  − ),  being the circumference of the ring resonator, in agreement with the analytical solution for notch-type filters [24].
In order to reach agreement between simulation results and physical measurements, we define  eff = / eff (, ), scaling the vacuum wave number  by the effective refractive index  eff (, ). eff (, ) accounts for the waveguide geometry (as opposed to bulk lithium niobate), as well as its dispersion (-dependence), and is calculated using Lumerical MODE software.Birefringence of LNOI material is included by the parameter , the angle between a waveguide and the ordinary crystal axis on an x-cut lithium niobate thin film.The -dependence of  eff is calculated using finite element method, from which the index ellipsoid is obtained.
For the ring resonator, we assume  eff is uniform across the ring, and the value that it takes is the angle-average The free spectral range (FSR) is used to quantify the transmission spectra of ring resonators for comparison.In Fig. 2c, the FSR simulated using the scattering matrix model is compared against measurement data from several physical devices, as well as the theoretical value (FSR=  2 /  ,   the group index).Within the 1.49-1.63m wavelength window, the angle-averaged  eff is fitted by the polynomial   () = 0.0006 2 − 0.2652 + 2.3363, with wavelengths in  units.
Excellent agreement is reached between simulation results and several devices with different coupling strengths between the ring resonator and the bus waveguide (Fig. 2c), which verifies the resonance behaviour to be independent of coupling strength, and demonstrates good fabrication quality consistent across devices.Agreement with the theoretical value (where   = − / is computed with the  =  avg ()) further verifies that by encoding material properties in a single parameter  eff (, ), the model can capture quantitative features in the transmission spectra of integrated LNOI devices.

Interferometric resonator
An interferometric resonator is composed of a racetrack-shaped loop and two bus waveguides, each coupling to the loop at two points (Fig. 3a).It serves as an example device with more complexity than the single ring resonator [25], and where different orientations of waveguide components become important.The geometry is approximated by the network schematics shown in Fig. 3b, where the bent sections are approximated by straight waveguides of the same length, and the  eff () along the edge is taken to be the same as a straight waveguide placed at the average angle with respect to the thin film.For these devices in the given wavelength range, we have  eff (, 0) = 0.0327 2 − 0.3953 + 2.4392,  eff (, /2) = −0.0224 2 − 0.2353 + 2.3687, and interpolated for the tilted edges  eff (, 0.6622) = 0.0127 2 − 0.3375 + 2.4141, with  in  units.The simulated and measured throughput spectra for such a device are shown in Fig. 3c.The simulated spectrum (top) replicates the key features in the measured one, namely the higher frequency resonances of the central loop, as well as the modulating envelop resulting from the unbalanced Mach-Zehnder interferometers formed by the edge of the racetrack and the bus waveguides.A fast Fourier transform on the simulation and measurement data yield a peak at 0.175  −1 and 0.198  −1 respectively, corresponding to 7 and 7.92 envelop cycles in the 40nm window shown in Fig. 3c.The discrepancy may be attributed to the measured spectra being less evenly sampled and with more noise present, as well as to the approximation of bent sections by straight edges.Nonetheless, this example demonstrates that even with realistic simplifications to the device geometry, the simulation captures key characteristics in the transmission spectrum of a physical device.

Agreement with analytical solutions
Besides existing physical samples, there are also photonic networks one may wish to characterise before fabrication, but do not have simple solutions like the single ring resonator.A small-scale network with increased complexity is given by a system of N serially coupled ring resonators.We consider the case of N=2 coupled resonators and compare them with the available analytical solution [26], as shown in Fig. 4. The transmitted spectral shape varies significantly for different ring sizes ( 1 ,  2 ), as well as different  and  parameters at the degree-4 coupling nodes, all of which are captured exactly by simulation.This demonstrates that for predicting the transmission properties of composite structures, the model performs just as well as traditional analytical methods, without the need for extensive algebraic calculations for each geometry.Together, the three examples in this section show the scattering matrix model to be a reliable predictor for the transport properties in an integrated LNOI network, while being a drastic simplification from numerical methods such as FDTD and FEM.The combined accuracy and computational tractability thus allow us to apply the scattering matrix model to scaled up optical networks, in particular as a discretised multiple-scattering system.

Developing the model for larger networks
Having verified the scattering matrix model for small-scale systems, we consider larger networks characterised by random disorder, which are used as planar, on-chip platforms for multiple scattering of photons.Differently from light transport in 2D random scattering media [27,28], in random photonic networks scattering and propagation are totally decoupled.Moreover, output channels are limited to few open edges and no further spatial discretisation is needed.When all these features are utilised, advanced control on light transport and deposition throughout the network can be achieved.A key feature of this model is the ability to introduce arbitrary scattering functions at the nodes of a photonic network.We present numerical experiments that demonstrate how tuning the scattering anisotropy at individual nodes modifies overall transport properties.This is a useful tool for many applications requiring control of the dwell time of light within the network without affecting the network's size .

Randomised grid-like networks
The geometry that we select as a model to investigate light scattering in networks is inspired by unidirectional forward-scattering optical circuits, which include beam splitters (realised via evanescent couplers) and phase shifters to couple and manipulate light modes travelling in waveguides: for example, to achieve universal unitary gates in quantum optical circuits [29], and in particular LNOI quantum photonic processors [30].This is translated into grid-like networks with degree-4 nodes which, under the appropriate scattering matrix, can act exactly as forward-scattering beam splitters (Fig. 5a).The size of the network is characterised by two parameters: the width, or the number of ports on each side   , and the depth, or the number of scattering layers   .Starting from regular grids, we transition into the multiple scattering regime by including backscattering at the nodes, and introduce disorder by adding a random displacement to each node, resulting in Gaussian distribution in the edge lengths (Fig. 5b).This induces a random phase that is picked up by light between two scattering events.The relevant geometric parameters are detailed in the caption of Fig. 5.The scattering matrix at the degree-4 nodes are selected to be of the form with ,  ∈ [0, 1] being the splitting and backreflection parameters, respectively. 4 is unitary and symmetric, and for  = 0 reduces to a unidirectional forward-scattering beam splitter.Such matrices are not fully isotropic due to the presence of zero entries; they are nevertheless selected for the ease of parametrisation while satisfying the constraints of an ideal physical device (lossless and reciprocal).Such a matrix is realised experimentally by attaching two Bragg reflectors with reflectivity  to a beam splitter with splitting ratio  (see inset of Fig. 5a).
Due to the random nature of the network, the transmission at the output edges creates discretised speckle-like patterns which we use to monitor transport properties.The model can simulate solutions independently at a range of different wavelengths, which together generate transmission spectra between arbitrary input and output ports.

Transmission and internal field
Given the geometry and scattering behaviour described above, a network is specified by parameters   ,   ,  and .In the following studies, we keep the network width   = 30 and node splitting ratio  = 0.5 fixed, and vary the network depth   and node reflectivity .
For each network, the transmission spectrum between any two open ports can be computed by selecting the relevant element in the overall scattering matrix   .This is a much larger set of possible measurements compared to smaller devices described in the previous section.A few example spectra of varying network sizes   = {9, 25, 41} and backscattering strengths  = {0.1,0.2, 0.3} are shown in Fig. 6.Increasing   and  produces sharper resonances, as well as more complex spectra with more resonance peaks.
Since the graph model contains solutions on the entire network, one can go beyond inputoutput relations and look at internal field distributions.This information is inaccessible in experiment without introducing additional out of plane losses, and the simulation therefore provides additional insight into the scattering process for different network geometries.Fig. 7 shows the corresponding field distribution at maximum and minimum transmission wavelengths for a subset of the spectra shown in Fig. 6.For the   = 9 network, the field profiles show a distinct difference between high and low transmission between specified input and output ports: in the minimum case (Fig. 7b), the field is split and sent to either side of the target output, whereas for maximum transmission (Fig. 7a) it is focused on the same port.For the larger network (Fig. 7c), high transmission is generated by exciting a series of coupled resonances between the input and output, akin to the necklace states observed in other multiple scattering systems [31].The excited field has a higher intensity near the centre of the network, in contrast to the linear decay with network depth in the diffusive regime.The individual transmission channels provide a convenient experimental analogue to smaller devices, but for large networks it is instructive to look at the entire scattering matrix   , an example of which is shown in Fig. 8a (amplitude) and 8b (phase).Since the open ports are divided into the ones on the left and the right of a network (Fig. 5a), the   matrix is divided into four quadrants: two diagonal blocks of reflection back to the same side, and two off-diagonal blocks of transmission to the opposite side.We write this as  = , where  :=   =    due to reciprocity.The separation between  and  is clearly visible in the amplitude matrix, whereas the phase matrix shows fully mixed phases between − and .Furthermore, each column in   is the speckle pattern resulting from a single-port input, with a reflected component (in   or   ) and a transmitted component (in   or   ).The wavelength-dependent speckle generated by the random geometry is a resource that can be utilised, for example in a random spectrometer [32].The spectral resolution is provided by decorrelation between the speckle patterns at different wavelengths.Fig. 8c shows the transmitted speckle spectral correlation for networks with different   and , averaged over all input ports and wavelengths.A higher spectral resolution is observed for networks with larger sizes, as well as stronger backreflections at the nodes.The observed behaviour can be interpreted in the complementary view of delay time in multiple scattering processes: the resolution of an interferometer is improved when longer optical paths come into play.Therefore, spectral decorrelation results confirm the capability of this platform to increase light circulation via multiple scattering, and show that average properties can be tuned by the parameters  and   .One can also go beyond the speckle pattern from single-port inputs, and study the transmission matrix  as a whole by decomposing them into transmission eigenchannels.Fig. 8d shows the eigenvalues of  † , averaged over 500 random realisations, of networks with  = 0.1 and 0.3, and at varying depths.Here a realisation means one instance of a network generated with the specified type of randomness, i.e. uniform random displacement of the nodes.Ergodicity is the equivalence between averaging over wavelengths and realisations, and the reason why networks satisfy this condition will become clear in section 4.3.The abundance of open channels, i.e. an increase in density of near-unity eigenvalues is observed, in accordance with the bimodal distribution expected for continuous random scattering [33].This effect is diminished for higher , suggesting excessive backreflection at the nodes results in transmission deviating from the diffusion regime.This deviation can also been seen in the scaling of transmission with network size, discussed below.

Transport properties on a discrete length scale
To statistically study transport properties on network ensembles with different design parameters, we look at the scalar value    , i.e. the total transmitted power from unit input.Scaling laws, which relate    to system size, are a useful tool for extracting information from multiple scattering media, for example to measure the transport mean free path  * and the localisation length .Here we study both  * and  on random networks, with the added caveat that network size is better characterised by discrete measures (  ) rather than the physical size of the system.If one extends all edges in the network by a certain factor, a different random phase will be traversed between any pair of connected nodes, but the overall nature of the random scattering remains unchanged (see Fig. 5b: the edge length distribution will remain Gaussian, with sufficient variance to accommodate 2 phase difference between edges).Therefore, stretching or shrinking the size of a network is equivalent to generating another random realisation of a network of the same size.Since stretching the edge lengths is complementary to reducing wavelength, this edge-length-invariance (combined with wavelength-independent scattering at the nodes) further implies that the network scattering matrices are ergodic, meaning it is statistically equivalent to average over wavelengths and random realisations.The number of layers   is therefore is chosen to be the discrete measure of size in our grid-like networks, instead of the physical scale.As a result, the computed  and  * only obtain a physical length scale by relating them to relevant dimensions in the network.
For localisation length, we define a size-invariant version   such that  =     , where   is the average distance between layers in the grid network.This value is obtained via a linear fit to the scaling of ⟨ln(   )⟩ with   , where the slope is taken to be −1/  [31,34].The linear fit is performed for   ≥ 25, where localisation has begun to set in.The scaling is shown for various backscattering strengths in Fig. 9a, with   = 34.43,29.32, 24.28, 20.10 respectively for  = 0.1, 0.15, 0.2, 0.25.As expected, increased level of backscattering leads to stronger localisation.
Fig. 9b shows the scaling of ⟨1/   ⟩ with   .A linear trend following optical Ohm's law is observed for weaker backscattering strengths, whereas larger  deviate from this diffusive regime.Two sets of simulation are shown, one for  = 1500  (crosses) and one for  = 750  (circles), and both overlap.This is an illustration of the edge-length-invariance and ergodic property discussed above: since no wavelength-dependence is introduced at the scattering nodes, the statistical transmission properties over a random ensemble do not change with wavelength.
In analogy with scattering in continuous random media, one would compute  * using the slope of the linear fit following optical Ohm's law.However, such a law is derived by solving the diffusion equation in finite size continuous media [35], and it is not obvious how the method should be adapted for the discretised case on a planar network.Instead, we construct an embedded version of  * by relating it to the scattering mean free path   , which can be taken to be the average network edge length, and compare the results with the scaling data. * and   are related via the similarity relation  * =   /(1 − ), where  = ⟨cos⟩ = ()cos is the anisotropy factor, averaging over all scattering angles weighted by their scattering probability () [36,37].
Fig. 9c shows two simple schematics used to compute  for the nodes described in Eq. 6 and shown in Fig. 5a.The first one treats a node simply as a cross, characterised by the angle .In this case   () = − + (1 − ) + cos (1 − ) (1 − ).Similarly, one may take into account the symmetry of a beamsplitter node and construct a node symmetrically with respect to the input port, characterised by the angle , which gives   () = − + cos(1 − ).Setting  = 0.5 as before, Fig. 9d plots   () and   () for several different angles.In both cases the node approaches isotropic scattering ( = 0) for larger angles, or when  is increased.Comparing with the scaling of ⟨1/   ⟩ in Fig. 9a, we see that the onset of localisation occurs at  ≈ 0.2, which corresponds to roughly  ≈ 0.5 in the two node models.
The anisotropy is then plugged into the similarity relation to give  *  and  *  , shown in Fig. 9e as a function of  for various angles.These values are compared with the simulation results for the scaling of ⟨1/   ⟩: even though the exact form of Ohm's law is unknown, we assume it takes the general form (in analogy to finite size continuum media [38]) where   are constant parameters, and   is used as the measure for network size instead of .
Then for each , the inverse slope of the scaling shown in Fig. 9b (obtained via a linear fit) is linearly related to  * : By adjusting the fitting parameters  1 and  2 , the resulting  * () follows  *  and  *  for different scattering angles, for small  (Fig. 9b).The agreement between the two methods for computing  * fails for  ≥ 0.2, where the transmission scaling deviates from optical Ohm's law, and is badly fitted by a linear trendline.In this sense, calculating  * using simple node models and the similarity relation captures the the effect of tuning  on transmission scaling in the networks in the diffusive regime.Furthermore, if a network design has well-defined scattering angles at the nodes, this method provides a way to measure the phenomenological coefficients of optical Ohm's law, in a discretised planar setting.
Due to their discretised nature and the size-invariant scattering behaviours, extra care needs to be taken if one wishes to treat the networks as direct analogues of other 2D multiple scattering systems.However, this feature also allows for flexibility in network design, since larger photonic devices can be used as scattering elements without altering the scattering behaviour of the overall network.

Conclusion
By modelling light propagating in integrated waveguides as a scalar wave, and treating the propagation and scattering separately, we have constructed a simple numerical model for multiple scattering of light in large scale photonic networks.In particular, the application to the LNOI platform has been experimentally verified to reach consistent agreement with multiple physical devices, providing a computationally inexpensive tool for simulating photonic networks and circuits.Furthermore, the model has been applied to a class of randomised grid-like networks, where we study multiple scattering phenomena in systems with well-defined type of randomness and scattering behaviour at each site.Such networks produce spectral transmission behaviour similar to that in other random media, while being tunable by the backscattering strengths at individual nodes.The networks are also unique in the discrete and statistically edge-lengthinvariant nature, which allows for more flexibility in terms of network designs and on-chip fabrication.By combining the computational simplicity of a graph-based network model and the control of top-down LNOI fabrication, these photonic networks provide a new on-chip platform for studying 2D multiple scattering with user-specified disorder and control.

Fig. 1 .
Fig. 1.Network components in the scattering matrix model.(a) Amplitudes along an edge are related by the propagation matrix  = .Note that although they are equal in magnitude, a coordinate change exists between the propatation for   and Φ   − , as evidenced by their index orders.(b) Amplitudes at a node are related by the scattering matrix  =  .(c) An example for a closed network, with nodes of degrees 2, 3, and 4. (d) The same network becomes an open network after attaching leads/external nodes.External nodes, shown with hollow markers, are where one sets boundary conditions corresponding to experimental input and measures the output.

Fig. 2 .
Fig. 2. Experimental verification on a ring resonator.(a) SEM image of a notch type ring resonator structure fabricated on LNOI thin film.The box indicates the coupling region.(b) Graph used to simulate the ring resonator with the model described in Section 2. The ring could be represented with one curved edge, but is here shown as square to demonstrate the model does not have to match exact geometry of physical networks, as long as appropriate dimensions and scattering matrices are used.(c) Free spectral range (FSR) for the single ring resonator, theory (red dotted line) vs. measurement vs. simulation.Inset shows cross section at the coupling region labeled in (a): 300nm etched on 600nm thin film, top width 1000nm,   the distance between the bus waveguide and the ring structure.The resonances are unaffected by the evanescent coupling strength as expected, and robustness across devices indicates consistent fabrication quality.For the sake of clarity in presentation, the FSR for one in every four resonances is depicted.

Fig. 3 .
Fig. 3. Experimental verification on an interferometric resonator.(a) Optical microscope image of an interferometric resonator device.(b) Network geometry used to simulate the design, the arrow indicating the input and output ports for throughput transmission.(c) Throughput port transmission spectrum of the device, measurement (top) vs. simulation (bottom).The measurement data is normalised by fitting and subtracting a second order polynomial, to remove the background spectrum of the superluminescentlight emitting diode source.

Fig. 4 .
Fig. 4. Theoretical verification on serially coupled 2-ring resonators.(a) Graph representation of the device in the model.The throughput configuration, indicated by arrows for input and output ports, is used for computing the spectra.For a lossless and reflection-free system, this captures all the nontrivial transmission information.(b) Transmission spectra of two rings in series, with different ring sizes and coupling strengths at the degree-4 nodes. is the evanescent coupling parameter as in  4 defined for the ring resonator at the beginning of section 3.1, where  1 = cos(20) and  2 = cos(5).Theoretical results follow the analytical solutions given in [26].

Fig. 5 .
Fig. 5. (a) A randomised grid-like network (node markers omitted for legibility): light can enter and exit through the open ports on left and right, propagate along waveguides denoted by the lines, and scatter at nodes where the lines meet.The size of the network is defined by   , the number of open ports on one side of the network, and can be thought of as the width, and   , the number of internal layers.Inset: the node design corresponding to the chosen scattering matrix (Eq.6) consists of a beamsplitter section (blue underlay) and Bragg reflector section (orange underlay).(b) Edge length distribution of the network on the left, fitted with a Gaussian distribution.The average vertical distance between nodes is 3, horizontal 4, random displacement of each node uniformly sampled from [0, 1.5]  in each direction.For wavelengths around 1500, this ensures a fully randomised phase between two scattering events.

Fig. 6 .
Fig. 6.Transmission spectra between ports 10 and 40, for networks with   = 30 and (a)   = 9, (b)   = 25, and (c)   = 41.From top to bottom each geometry had local scattering matrices with reflectivities  = 0.1,0.2,0.3,respectively.As   increases, the peaks become narrower, similar to the transmission peak sizes in other multiple scattering systems.Furthermore, we see an increase in the number of peaks as  is increased.The triangles mark the maximum and minimum transmission points, for which the corresponding field distribution are shown in Fig. 7.

Fig. 7 .
Fig. 7. Internal field distributions for   = 30,  = 0.2 and (a) (b)   = 9, or (c)   = 41.Ports 10 and 40 are labelled with blue markers, and unit-power input is injected at port 10.(a) and (b) show the minimum and maximum transmission wavelengths respectively, whereas (c) shows the maximum transmission for a larger network.

Fig. 8 .
Fig. 8.The   matrix and transmission behaviour.(a) Amplitude and (b) phase components of the   matrix.Ports 1-30 are on the left side of the network, 31-60 on the right.Each column is the speckle pattern resulting from a single-port input at a single wavelength.(c) Spectral correlation between speckle patterns, averaged over all inputs and 500 random realisations.Sharper spectral resolution is observed for larger   or higher .(d) Transmission eigenspectra for networks with  = 0.1 and 0.3, with varying depths, averaged over 500 random realisations.The increase in the distribution at higher transmittance is observed, but the effect is diminished for higher .

Fig. 9 .
Fig. 9. Scaling of transmission with size of the network, characterised by number of layers   rather than physical size due to discreteness of networks.Averaged over 500 random realisations.(a) Linear scaling of ⟨ln(   )⟩ gives a discrete version of localisation length,   .(b) Scaling of ⟨1/   ⟩.For smaller  (<0.2), a linear trend is observed, following optical Ohm's law in the diffusive regime.For increased  (≥0.2), deviation from the linear trend indicate onset of localisation.Simulation results for  = 1500 (crosses) and  = 750 (circles) overlap due to the edge-length-invariant nature of scattering in the network system.(c) Two simple models for a degree-four node described by Eq. 6, with the corresponding (d) scattering anisotropy  and (e) transport mean free path  * shown for various scattering angles  and , calculated using the similarity relation.Black markers in (e) are  * values computed from numerical scaling results as in (b), scaled by different parameters  1 and  2 , defined in Eq. 8.They deviate from the similarity relation results for  ≥ 0.2, i.e. when optical Ohm's law no longer holds.