# Sensitivity and spectral control of network lasers

Dhruv Saxena  The Blackett Laboratory, Department of Physics, Imperial College
London, London SW7 2BW, United Kingdom  Alexis Arnaudon  Department of
Mathematics, Imperial College London, London SW7 2AZ, United Kingdom  Oscar
Cipolato  The Blackett Laboratory, Department of Physics, Imperial College
London, London SW7 2BW, United Kingdom  Michele Gaio  The Blackett Laboratory,
Department of Physics, Imperial College London, London SW7 2BW, United Kingdom
Alain Quentel  The Blackett Laboratory, Department of Physics, Imperial
College London, London SW7 2BW, United Kingdom  Sophia Yaliraki  Department of
Chemistry, Imperial College London, London SW7 2AZ, United Kingdom  Dario
Pisignano  NEST, Istituto Nanoscienze-CNR and Scuola Normale Superiore, Piazza
San Silvestro 12, I-56127 Pisa, Italy  Dipartimento di Fisica, Università di
Pisa, Largo B. Pontecorvo 3, I-56127 Pisa, Italy  Andrea Camposeo  NEST,
Istituto Nanoscienze-CNR and Scuola Normale Superiore, Piazza San Silvestro
12, I-56127 Pisa, Italy  Mauricio Barahona  Department of Mathematics,
Imperial College London, London SW7 2AZ, United Kingdom  Riccardo Sapienza
The Blackett Laboratory, Department of Physics, Imperial College London,
London SW7 2BW, United Kingdom

###### Abstract

Recently, random lasing in complex networks Gaio  _et al._ (2019) has shown
efficient lasing over more than 50 localised modes, promoted by multiple
scattering over the underlying graph. If controlled, these network lasers can
lead to fast-switching multifunctional light sources with synthesised
spectrum. Here, we observe both in experiment and theory high sensitivity of
the network laser to the spatial shape of the pump profile, with mode
intensity variation of up to 280% for a non-homogeneous 7% pump decrease. We
solve the nonlinear equations within the steady state ab-initio laser theory
(SALT) approximation Ge  _et al._ (2010) over a graph and we show selective
lasing of around $90$% of the top modes, effectively programming the spectrum
of the lasing networks. In our experiments with polymer networks, this high
sensitivity enables control of the lasing spectrum through non-uniform pump
patterns. We propose the underlying complexity of the network modes as the key
element behind efficient spectral control opening the way for the development
of optical devices with wide impact for on-chip photonics for communication
Wang  _et al._ (2017), sensing Fan and Yun (2014) and computation Shen  _et
al._ (2017).

Lasers with a well defined emission frequency and direction have
revolutionised many fields, from material processing to biophysics and
communication, just to mention a few. Traditionally, the spectral properties
of the laser are inherited directly from the modes of the passive cavity,
which is usually designed to suppress multimode lasing and favour single-mode
operation. In contrast, random lasers are an unconventional lasing
architecture where light is amplified in a multimode scattering medium, thus
supporting many lasing modes at random frequencies Sapienza (2019); Luan  _et
al._ (2015); Wiersma (2008). The ensuing low-coherence, multi-frequency,
fluctuating laser radiation has applications in low-coherence imaging Cao  _et
al._ (2019) and super-resolution spectroscopy Boschetti  _et al._ (2020), but
is not suited for technologies that require fine control of the lasing
emission at specific frequencies, such as signal processing, spectroscopic
sensing, communication or optical computing. An experimental challenge is
therefore how to achieve spectral selection in a controlled manner from such
random lasing architectures. Indeed, spectral selection has been observed in
powder random lasers when the pump laser was elongated in one direction
Leonetti and Lopez (2013), in disordered toroidal cavities with varying
spatial patterns of the pump laser Liew  _et al._ (2014), and in a one-
dimensional opto-fluidic random laser excited by a structured pump profile
Bachelard  _et al._ (2014).

Recently, a novel type of random lasers called _network lasers_ was introduced
in Gaio  _et al._ (2019); Lepri  _et al._ (2017). Network lasers consist of
active single-mode waveguides connected according to a network topology. The
passive modes of such systems are captured by quantum graphs Gaio  _et al._
(2019) and scattering matrix models Lepri  _et al._ (2017). Yet to take into
consideration mode competition and nonlinear interactions, one must go beyond
such passive models and solve the Maxwell-Bloch equation Conti and Fratalocchi
(2008) or its steady state ab-initio laser theory (SALT) approximation Ge  _et
al._ (2010) on a graph. This leads to a problem in _nonlinear quantum graphs_
, recently studied in the context of the nonlinear Schrödinger equation
Gnutzmann and Waltner (2016); Besse  _et al._ (2021), but not yet considered
to formulate the spectral control of network lasers.

Beyond photonic systems, how to design network structures or their inputs to
produce specific dynamic behaviours is a central question in many areas, such
as in the haemodynamics of arterial networks Parker and Jones (1990), power
grids Schaub  _et al._ (2014), brain networks Gu  _et al._ (2015), or acoustic
waves in elastic networks Hu  _et al._ (2008). In conventional networks Newman
(2003), simple graph-theoretical measures are often sufficient to controllably
characterise and produce network outputs Schaub  _et al._ (2019). However,
such simple network measures are rendered unsatisfactory in nonlinear quantum
graphs due to the complex interplay between graph structure and dynamical
processes Berkolaiko and Kuchment (2013). Here we show that the underlying
complexity of the nonlinear quantum graphs associated with random lasing can
be harnessed to achieve a high degree of design control on the lasing
emissions. We demonstrate experimentally and numerically that the complex
emission spectrum of nanophotonic network lasers can be efficiently and
precisely controlled through optimisation of spatially non-uniform pump
patterns.

##  I Results and Discussion

Figure 1: Spectral sensitivity of network lasers with pump profile. a
Fluorescence image of a photonic network with interconnected dye-doped polymer
nanofibres. As a guide to the eye, we highlight the graph topology (edges as
blue lines, nodes as orange dots) over a few links as an example. Inset:
scanning electron microscopy image of a node of the network formed by
annealing two nanofibers (scale bar = 1 $\mu$m). b Lasing spectrum of the
network in a obtained with uniform illumination (blue) and with a slightly
modified pump pattern (yellow), both at pump fluence of $1.5$ mJ cm-2 pulse-1.
The insets show the respective pump patterns; the illumination area on sample
is $300\times 480~{}\mu$m and the modified pattern has pump removed from a
small area $60\times 60~{}\mu$m in the centre of the illuminated rectangle. c
A model planar photonic network modelled as a Buffon graph open at the
boundaries. We highlight edges with more than $50\%$ of the maximal amplitude
of the electric field for a delocalised mode (orange edges, mode 1) and a
localised mode (green edges, mode 2). d Numerical calculations of the lasing
spectrum from the Buffon network in c obtained with netSALT. Spectra at pump
power $D_{0}=0.01$ (SALT units) show $\sim$50 lasing modes within the gain
spectrum of dye (red dashed line). Note the suppression of mode 1 when
changing from uniform pumping (blue, pump profile in inset) to a pump missing
the two edges supporting the largest electric field amplitudes for mode 1
(yellow, pump profile in inset).

Network laser spectral sensitivity. The network lasers examined here are
planar and built from dye-doped polymer nanofibers physically joined together
at the nodes Gaio  _et al._ (2019), resulting in graph-like structures with an
average node degree of $4$ and edge lengths ranging from $10-100$ $\mu$m (Fig.
1a). Lasing is experimentally obtained by optical pumping using a custom-built
lasing microscope (see Methods). When uniformly pumped over a $300\times 480$
$\mu$m2 rectangular area, the networks lase from multiple modes, with narrow
linewidths ($\sim$50 pm), as shown in Fig. 1b. These modes are formed by
interference of light over multiple closed loops in the network and amplified
by optical gain in the network links. Typically, $30$ to $100$ lasing modes
are observed within the gain bandwidth of the dye.

The lasing spectrum is very sensitive to changes in the spatial profile of the
optical excitation. When the experimental pump pattern is modified so that a
small central area of $60\times 60$ $\mu$m2 is not excited (corresponding to a
$7$% reduction of the net pump energy delivered to the sample), we observe a
drastic change in the lasing spectrum, which is stable upon multiple
illumination (Fig. 1b). Some modes are amplified (up to $280$%) while others
are attenuated (down to $20$%), and even new modes (not lasing under the
uniform pump) lase (see SI Fig. S1).

To understand the sensitivity of the network laser to non-homogeneous pump
profiles, we developed netSALT, which solves the nonlinear interaction of the
optical waves on the network, modelling the lasing process within the SALT
approximation Ge  _et al._ (2010) (see Methods and SI for full details). The
netSALT model includes amplification/loss on graph edges and mode competition.
Under uniform pumping, the predicted spectrum in Fig. 1d is qualitatively
similar to the experimental one, with similar number of modes (see SI Fig.
S2c).

The high sensitivity arises because of a large number of modes competing for
gain, some delocalised and other localised, as shown in Fig. 1c. The network
modes are spatially coupled as they partially overlap on graph edges; in this
particular graph, there are 450 modes within a spectral range of 35 nm. If we
select the mode with highest modal amplitude (mode $1$ in Fig. 1d) and turn
off the pump illumination from the two edges supporting the largest electric
field amplitudes for this mode (pump profile shown in inset), mode $1$ does
not lase anymore and overall most lasing modes change amplitude (see Fig. 1d
and SI Fig. S2a).

Figure 2: Achieving single mode lasing through non-uniform pumping (theory).
a-b NetSALT calculations for lasing on a Buffon graph with non-uniform pumping
indicate various processes for mode selection. a Mode trajectories in the
complex $k$ plane (shown over large Re($k$) range, with red dots identifying
the passive modes and black dots identifying the modes that lase. The non-
uniform pump profile used to obtain these trajectories is shown in d, right-
hand panel inset. b Zoom of a over a small range of wavelengths with
trajectories shown in colour, where the colour scale indicates the pump
strength. Three modes (1, 99 and 134) are highlighted to exemplify the
different processes involved in mode selection. Mode trajectories with uniform
pumping are shown by solid grey lines. c Example mode profiles for two network
modes (1 and 134) before pumping (passive, $D_{0}=0$) and at threshold for
mode 1 ($D_{0}=0.004$) and at pump strength $D_{0}=0.01$ for mode 134. Minimal
change in profile is observed for mode 1, whereas mode 134 reshapes
significantly. d Modal intensities as a function of pump power for a uniform
pump and a pump chosen to maximise lasing of mode $1$ (see Methods), with pump
profiles represented as green edges in insets. e Heatmap of the modal
amplitudes of the first $50$ modes (along each row) under $50$ patterns
optimised for each mode (each column). The computations correspond to the top
50 modes of the Buffon graph ordered in descending order of $\mathcal{Q}$
factor. Optimisation of the pump profile leads to good mode selectivity.

Processes involved in mode selection. The high sensitivity of the network
laser to pump illumination can be used for designing the pump to either select
or suppress lasing from certain modes. The main underlying processes that
determine the lasing spectrum under a non-uniform pump are: a) efficient
pumping of a mode to reach threshold at lowest pump power compared to all
other modes; b) mode reshaping; and c) mode competition. To illustrate these
processes, we use netSALT to calculate the lasing modes of a complex network
when illuminated with a non-uniform pump profile with 50% fill fraction (shown
in inset of right panel in Fig. 2d), chosen to maximise the modal amplitude
for mode 1. As shown in Fig. 2a, the network has 454 passive modes (red dots)
in the region of the complex $k$ plane chosen for calculations, $k$ being the
complex wavenumber. Of these, only $208$ modes reach lasing threshold
(|Im($k$)|$=0$) when pumped with a pump strength $D_{0}=0.01$ with the given
non-uniform pump profile. Of these $208$ modes that can potentially lase, only
$15$ modes (black filled circles) actually lase due to strong mode competition
according to netSALT calculation.

Fig. 2b is a zoom-in of the complex $k$ plane, highlighting three modes,
labelled 1, 134 and 99 and their trajectories for uniform (grey line) and non-
uniform (coloured line) illumination (see also SI Fig. S3). These trajectories
are obtained by calculating and tracking the modes as we increase $D_{0}$,
with small increments ($1.4\times 10^{-4}$). Some modes (e.g. mode 1) move
directly and rapidly towards the lasing threshold (Im$(k)=0$) under the
increase of $D_{0}$, while others (e.g. mode 99) undergo nonlinear shifts in
resonance frequencies and thus reach lasing threshold at higher values of
$D_{0}$. In other cases, modes (e.g. mode $134$) can move away from the lasing
threshold and never reach it within our range of pump power. This behaviour is
due to the second process, i.e. mode reshaping, which changes the mode
amplitude on each edge and therefore modifies the condition for resonance. If
this reshaping moves a mode towards the boundaries of the network, such as for
mode $134$ (see Fig. 2c), the mode becomes so lossy that the increase of the
pump power is not enough to reach threshold; hence the trajectory of this mode
in the $k$ complex plane stalls. Lastly, mode competition, which is the
nonlinear interaction for gain above threshold due to spatial hole burning,
affects the modal intensity of competing modes as well as their effective
lasing thresholds, also called interacting lasing thresholds Ge  _et al._
(2010). Mode competition depends on many factors, including spatial overlap
between modes, mode frequencies with respect to the gain spectrum and pump
power required to reach threshold, and therefore is affected by the pumping
efficiency and mode reshaping.

Theoretical modal control. The complex modal interaction and the before-
mentioned three processes can be exploited to achieve mode selection by
adaptive pumping. We give one example in Fig. 2d, where we compare uniform
pumping against a pump optimised to lase mode $1$ (shown in Fig. 2b). The
light-in light-out (LL) curves (or modal amplitudes as a function of pump
power) show clear improvement in the mode suppression ratio for mode $1$ with
optimised pumping. Even if other modes lase at higher pump power $D_{0}$, the
intensity of this mode dominates across the power range. Notice that mode $1$
lases first, with a large gap of lasing threshold with the next lasing mode
(see SI Fig. S4).

In general, finding the right illumination pattern to achieve a desired lasing
spectrum, e.g. single mode operation, is not a trivial task (see Methods for a
description of the optimisation). Furthermore, to find pump patterns that are
physically relevant, experimental limitations on the pump spatial resolution,
pump power and optical gain have to be considered. We note that the naive
approach of pumping the edges where the target mode has a large electric field
does not always ensure single mode lasing, in particular for modes that are
spatially delocalised or have high losses (see SI Fig. S5). Instead, with
optimised pump profiles, we can lase $143$ out of the $200$ modes, with a
suppression ratio larger than one, and lase $102$ with a ratio larger than
two. The matrix in Fig. 2e shows in each row the modal amplitudes of the
optimised pumping of the first $50$ modes. These modes are arranged in the
matrix in order of $\mathcal{Q}$-factor, where mode 1 has the highest
$\mathcal{Q}$-factor. A large value on the diagonal corresponds to a good
performance single mode lasing and low values on the off-diagonal indicates
strong suppression of the unwanted lasing modes. We observe that control can
be achieved across a large frequency window, even far from the gain maximum,
as well as for relatively lossy modes (see SI Fig. S5). After pump
optimisation, 90% of the top 50 modes (and 70% of the top 200 modes) can be
controlled with amplitude more than double than any other mode. We remark that
most of the obtained optimal pump profiles have only a partial correspondence
to the target mode profiles.

Figure 3:  Experimental spectral control of network laser. a-b Lasing modes in
the network laser are controlled by selective illumination of the individual
network links. This is done by shaping the pump laser using a digital
micromirror device (DMD). The patterns are projected on the sample and cover a
rectangular area of $300\times 480$ $\mu$m2. c Plot of the quality function at
each iteration for the optimisation of mode at 585.9 nm (shown in red in d-f).
The pump pattern is optimised using a derivative-free, greedy iterative
algorithm (see Methods) that optimises the quality function Eqn. (1) to
improve the mode suppression ratio leading to a progressive suppression of the
unwanted modes while the selected one is maintained. d Lasing is controlled
from a multimode spectrum (grey) for homogeneous pumping to a single mode,
shown here for two examples at 583 nm (blue) and 586 nm (red), obtained when
illuminating with the patterns shown in insets. Bi-modal lasing, i.e.
enhancement of both modes at 583 and 586 nm (green) is also achieved under a
different illumination pattern. Fill fraction of the optimised patterns are
0.57, 0.46, 0.68, for the blue, red and green spectra, respectively. The
evolution of the emitted intensity as a function of pump fluence when pumping
with the uniform and optimised pump patterns are shown in e and f,
respectively.

Experimental demonstration of spectral control. The high mode selectivity of
network lasers predicted through numerical calculations is observed
experimentally. Following the approach of Bachelard  _et al._ (2014); Leonetti
and Lopez (2013); Liew  _et al._ (2015), we use a digital micromirror device
(DMD) to project different pump patterns on the sample (Fig. 3a-b),
discretised into binary intensity pixels. Limitations on the spatial
resolution of the pump, maximum amount of power available for pumping and the
amount of gain in the medium constrain the parameter space to find physically
relevant pump patterns. Patterns are optimised using a derivative-free, greedy
iterative algorithm (see Methods). The mode suppression ratio between the
target mode and other lasing modes is computed at each iteration to form a
quality function to be maximised (see Methods and Fig. 3c). The results of
such optimisations for the first and fourth largest modes under uniform
pumping (grey) are presented in the top and middle panels of Fig. 3d. The LL
curves under uniform and optimised pumps (Fig. 3e and f, respectively) shows a
successful suppression of undesired lasing modes, while maintaining the
intensity of the target mode. Additional results of single mode lasing
optimisation from different areas of the sample and at larger pump power are
given in SI (see SI Fig. S6 and S7). Furthermore, we experimentally
demonstrate that it is possible to optimise for concurrent lasing of several
modes, as shown for two modes in the bottom panels of Fig. 3d-f. We
numerically confirmed this result in SI Fig. S8, and assess the experimental
stability of the spectra by switching repeatedly between different pumps in SI
Fig. S7.

Discussion. Our procedure for single-mode optimisation converges in $\sim$100
steps, both in experiments and theory, while exploring a configuration space
of $2^{160}$ configurations (for a 16x10 pixels discretisation). This
remarkable efficiency indicates the existence of correlations emerging from
the underlying physical constrain of light waves on a network, whose modes are
not random and discontinuous, but localised and continuous.

Underlying all the physical processes that occur in a network laser when
optically pumped, the network and its specific topology dictates the form of
the lasing spectrum. The network topology defines the mode density, their
$\mathcal{Q}$ factor and spatial distribution, which affects controllability
of the lasing spectrum. For example, as an extreme case, simple graphs such as
rings (see SI Fig. S9) have fewer modes, with only one that can be controlled
in the centre of the gain. In contrast, very large planar networks will
contain many modes over a wide range of localisations with many possible
single lasing regimes. The large controllability of our network lasers stems
from its high structural complexity, with many cycles, producing multiple
scattering from the disordered graph. If complexity is reduced by adding
spatial correlations, as in a periodic network, the degree of lasing control
is strongly reduced (see SI Fig. S10).

In conclusion, we have shown experimentally and numerically that network
lasers inherently display a large spectral control, for over 90% of the top
$50$ modes, via the choice of the pump profile. The degree of control stems
from the network complexity, and could be increased for further flexibility or
decreased for improved resiliency. Further design of lasing networks may lead
to improved spectral and directional control Hisch  _et al._ (2017) and could
also be extended to other systems described by wave propagations on networks
Parker and Jones (1990); Peach  _et al._ (2021). In addition, network lasers
are naturally adapted for on-chip integration, and could be made out of
semiconductor materials to power next generation programmable light sources
Wang  _et al._ (2017), optical sensors Fan and Yun (2014) and neuromorphic
optical processors Shen  _et al._ (2017).

##  II Methods

Experiments on single mode lasing optimisation. Polymer nanofiber networks
were pumped using a $\lambda=532$ nm pulsed laser (TEEM Microchip, pulse width
500 ps) and the emission was detected using a grating spectrometer (Princeton
Instruments Isoplane-320) with 1800 gr mm-1 holographic grating ($0.05$ nm
resolution) and CCD camera (Princeton Instruments Pixis 400). A digital
micromirror device (DMD, Ajile AJD-4500) was used for beam shaping, resulting
in a rectangular illumination spot of 300 $\times$ 480 $\mu$m on the sample.

A derivative-free, greedy iterative algorithm was used to find the optimised
pump patterns. Firstly, a coarse grid (8 $\times$ 5 grid with each pixel
corresponding to 60 $\times$ 60 $\mu$m size on sample) was used. Starting from
the pixel closest to the centre of the grid, each pixel was switched off
consecutively and the change in the intensity of the selected lasing mode was
calculated using the recorded spectral counts. For a given lasing peak $p$, we
calculated the following quality function at each optimisation step $n$:

$\displaystyle\Phi_{n}=\frac{a_{n}}{a_{n-1}}-1\text{, with
}a_{n}=\dfrac{I_{p}(n)}{\frac{1}{M}\sum_{m=1;m\neq p}^{M}I_{m}(n)}$ (1)

where $a_{n}$ is the ratio of the intensity of the selected lasing peak $p$ to
the average intensity of the top $M$ strongest lasing peaks ($M=10$ in our
experiments), all under pump pattern $n$. If $\Phi_{n}>0$, the patch was kept
off the pump, otherwise it was switched back on, and the routine was iterated.
The final pattern from a first run was then fed as the initial pattern for a
subsequent re-run with a finer grid (patch sizes of 30 $\times$ 30 $\mu$m) for
further optimisation.

Numerical construction of Buffon graphs. Buffon graphs were generated by
drawing lines on a plane at random points with random slope. The intersections
of all the lines within a square region on the plane were obtained and the
length of the line segments between intersections calculated. If a segment
length was smaller than a minimum distance of $1$ $\mu$m, the intersection
points were merged together to the median point. The final set of intersection
points and line segments was then used to specify the graph vertices and
adjacency matrix. The Buffon graphs used for numerical calculations were
constructed to be similar to the polymer nanofiber networks, with 96 nodes,
131 edges, average degree 4, and mean edge length 23.8 $\mu$m.

Numerical model: SALT on networks (netSALT). Lasers are usually described with
two-level Maxwell-Bloch equations and numerically solved using finite
difference methods Ge  _et al._ (2008). An alternative, computationally
efficient approach is to approximate these equations assuming stationarity of
the population inversion and adopting the slowly-varying envelope
approximation, resulting in the so-called SALT model Türeci  _et al._ (2008);
Ge  _et al._ (2010). The SALT model can be solved for arbitrary geometries,
provided an efficient solver is available to compute the mode profiles in the
lasing cavity.

Here, our cavity has the structure of a complex network, which we approximate
as a quasi-1D system, where edges of the network are simple 1D cavities
coupled via the nodes of the graph. This assumes that most of the light
propagates in the direction of the edges, and that the complex scattering
processes at the nodes can be well approximated with Neumann boundary
conditionsKuchment (2004). These two approximations are fundamental for what
we call the netSALT model, i.e. SALT on networks. For full details on netSALT,
see the SI, but we give here a summary.

The SALT equation for a one-dimensional cavity is

$\displaystyle\partial_{x}^{2}u_{\mu}+\left(\epsilon_{ij}+\frac{D_{0}(\delta_{\mathrm{pump}})_{ij}\gamma_{\mu}}{1+\sum_{\nu}I_{\nu}\Gamma_{\nu}|u_{\nu}|^{2}}\right)k_{\mu}^{2}u_{\mu}=0\,,$
(2)

where $u_{\mu}$ is the normalised mode electric field and
$\delta_{\mathrm{pump}}$ is the pump profile (equal to $1$ on edges
illuminated by the pump and $0$ otherwise). The other parameters are:
$\epsilon_{ij}$ the dielectric constant on each edge; $D_{0}$ the pump
strength;
$\gamma_{\mu}=\gamma_{\perp}/((k_{\mu}-k_{\text{a}})^{2}+i\gamma_{\perp})$ the
Lorentzian gain curve; and $\Gamma_{\mu}=-\mathrm{Im}(\gamma_{\mu})$ the gain
linewidth. The electric field ($u_{\mu}$) and pump strength ($D_{0}$ ) in SALT
equations are dimensionless and can be converted to physical units (SI) of
electric field and inversion density, using
$E_{\mu}=u_{\mu}(\hbar\sqrt{\gamma_{\perp}\gamma_{\parallel}})/2g$
$\text{Vm}^{-1}$ and
$D=D_{0}(\epsilon_{0}\hbar\gamma_{\perp})/k_{\text{a}}^{2}g^{2}$
$\text{cm}^{-3}$ Ge  _et al._ (2010). The parameters $k_{\text{a}}$ and
$\gamma_{\perp}$ used in the netSALT calculations depend on the length units
chosen for the edges. In our calculations we used $k_{\text{a}}=10.68$
$\mu\text{m}^{-1}$ and $\gamma_{\perp}=0.5$ $\mu\text{m}^{-1}$.

To solve this equation, one needs the boundary conditions for each edge
matched at each node of the underlying network. We use the theory of quantum
graphs to derive a matrix equation for the electric field at the node. For
each edge, we have $\eta_{ij}(x)$ obeying

$\displaystyle\partial_{x}^{2}\eta_{ij}(x)+(n_{ij}k)^{2}\eta_{ij}(x)=0\qquad\forall\,(ij)\,,$
(3)

where $n_{ij}$ is the index of refraction of the edge $(ij)$. This has
solutions of the form

$\displaystyle\eta_{ij}(x)=\lambda^{+}_{ij}e^{ikn_{ij}x}+\lambda^{-}_{ij}e^{ikn_{ij}(l_{ij}-x)}\,,$
(4)

where $\lambda^{\pm}_{ij}$ are the wave amplitudes, one to one with the wave
amplitude $\eta_{i}$ at node $i$. One can recast the boundary conditions at
the nodes into a matrix $L(k)$ (see SI), such that the passive modes with
wavenumber $k_{\mu}$ satisfy

$\displaystyle L(k_{\mu})\,\bm{\eta}=0\,,$ (5)

where $\bm{\eta}$ is the vector containing the node wave amplitudes $\eta_{i}$
as components.

The wave equation Eqn. (2) (see SI) with nonlinear coupling between modes
cannot be solved directly, but we obtain an approximation in several steps.
First, we search for passive modes (without pump), i.e. with $D_{0}=0$. These
modes have a complex wavenumber $k_{\mu}$, whose imaginary part is related to
the loss of the mode via the standard $\mathcal{Q}$-factor

$\displaystyle\mathcal{Q}_{\mu}=\frac{\mathrm{Real}(k_{\mu})}{2\,|\mathrm{Im}(k_{\mu})|}\,.$
(6)

For each mode, we then search for the pump power $D_{0,\mu}$ for which
$\textrm{Im}(k_{\mu}(D_{0,\mu}))=0$ where $k_{\mu}(D_{0})$ solves Eqn. (2)
without the denominator in the nonlinear term. The wavenumber obtained is the
one of the so-called threshold lasing mode denoted here $u_{\mu}$. We then
assume that above lasing, these modes do not change their profile
significantly, thus the nonlinear coupling between the lasing modes due to the
spatial hole burning term can be approximated by a matrix equation (see SI).
The lasing modes obtained through the approximated solution to Eqn. (2) are
then given as $I_{\mu}u_{\mu}$, where $I_{\mu}$ is the mode amplitude computed
from this matrix equation.

Numerical individual mode lasing optimisation. To numerically optimise a pump
profile to single lase a specific mode, we would ideally maximise the ratio of
modal amplitude of the target mode over the largest next lasing mode. However,
as this quantity is numerically expensive to compute (due to the need to track
modes in the complex plane), we approximate it using the overlapping factor
Eqn. (17), as an indication of the change of lasing threshold, and write the
optimisation as a linear program (see SI).

##  III Acknowledgements

RS and DS acknowledge support from The Engineering and Physical Sciences
Research Council (EPSRC), grant number EP/T027258, and the European Community.
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 800410. M. Moffa and A. Portone are acknowledged for sample
preparation. AA, SNY and MB acknowledge funding from EPSRC grant EP/N014529/1
supporting the EPSRC Centre for Mathematics of Precision Healthcare.

##  IV Code availability

Code with examples available on github at https://github.com/arnaudon/netSALT.

##  V Author contributions

Author contributions are defined based on the CRediT (Contributor Roles
Taxonomy). Conceptualisation: RS, MB. Data Curation: DS, AA. Formal analysis:
DS, AA, OC, MB. Funding acquisition: DS, MB, RS. Investigation: DS, AA, OC,
MG, AQ, AC. Methodology: DS, AA, OC, MG, AQ, MB, RS. Project administration:
MB, RS. Software: AA. Supervision: SY, DP, AC, MB, RS. Validation: DS, AA.
Visualisation: DS, AA, OC, MB, RS. Writing – review and editing: All. DS and
AS contributed equally to this work.

##  VI Competing financial interests

The authors declare no competing financial interests.

## References

  * Gaio  _et al._ (2019) M. Gaio, D. Saxena, J. Bertolotti, D. Pisignano, A. Camposeo, and R. Sapienza, Nature Communications 10, 226 (2019).
  * Ge  _et al._ (2010) L. Ge, Y. Chong, and A. D. Stone, Physical Review A 82, 063824 (2010).
  * Wang  _et al._ (2017) Z. Wang, A. Abbasi, U. Dave, A. D. Groote, S. Kumari, B. Kunert, C. Merckling, M. Pantouvaki, Y. Shi, B. Tian, K. V. Gasse, J. Verbist, R. Wang, J. Z. Weiqiang Xie, Y. Zhu, J. Bauwelinck, X. Yin, Z. Hens, J. V. Campenhout, B. Kuyken, R. Baets, G. Morthier, D. V. Thourhout, and G. Roelkens, Lasers and Photonics Review (2017).
  * Fan and Yun (2014) X. Fan and S.-H. Yun, Nature Methods 11, 141 (2014).
  * Shen  _et al._ (2017) Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, and M. Soljačić, Nature Photonics 11, 441 (2017).
  * Sapienza (2019) R. Sapienza, Nature Review Physics 1, 690 (2019).
  * Luan  _et al._ (2015) F. Luan, B. Gu, A. S. Gomes, K.-T. Yong, S. Wen, and P. N. Prasad, NanoToday 10, 168 (2015).
  * Wiersma (2008) D. S. Wiersma, Nature Physics 4, 359 (2008).
  * Cao  _et al._ (2019) H. Cao, R. Chriki, S. Bittner, A. A. Friesem, and N. Davidson, Nature Review Physics 1, 156 (2019).
  * Boschetti  _et al._ (2020) A. Boschetti, A. Taschin, P. Bartolini, A. K. Tiwari, L. Pattelli, R. Torre, and D. S. Wiersma, Nature Photonics 14, 177 (2020).
  * Leonetti and Lopez (2013) M. Leonetti and C. Lopez, Applied Physics Letters 102, 071105 (2013).
  * Liew  _et al._ (2014) S. F. Liew, B. Redding, L. Ge, G. Solomon, and H. Cao, Appl. Phys. Lett. 104, 231108 (2014).
  * Bachelard  _et al._ (2014) N. Bachelard, S. Gigan, X. Noblin, and P. Sebbah, Nature physics 10, 426 (2014).
  * Lepri  _et al._ (2017) S. Lepri, C. Trono, and G. Giacomelli, Physical Review Letters 118, 123901 (2017).
  * Conti and Fratalocchi (2008) C. Conti and A. Fratalocchi, Nature Physics 4, 794 (2008).
  * Gnutzmann and Waltner (2016) S. Gnutzmann and D. Waltner, Physical Review E 93, 032204 (2016).
  * Besse  _et al._ (2021) C. Besse, R. Duboscq, and S. L. Coz, arXiv preprint arXiv:2103.09650 (2021).
  * Parker and Jones (1990) K. H. Parker and C. J. H. Jones, Journal of Biomechanical Engineering 112, 322 (1990).
  * Schaub  _et al._ (2014) M. T. Schaub, J. Lehmann, S. N. Yaliraki, and M. Barahona, Network Science 2, 66 (2014).
  * Gu  _et al._ (2015) S. Gu, F. Pasqualetti, M. Cieslak, Q. K. Telesford, A. B. Yu, A. E. Kahn, J. D. Medaglia, J. M. Vettel, M. B. Miller, S. T. Grafton, and D. S. Bassett, Nature Communications 6, 8414 (2015).
  * Hu  _et al._ (2008) H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Nature Physics 4, 945 (2008).
  * Newman (2003) M. E. Newman, SIAM review 45, 167 (2003).
  * Schaub  _et al._ (2019) M. T. Schaub, J.-C. Delvenne, R. Lambiotte, and M. Barahona, Phys. Rev. E 99, 062308 (2019).
  * Berkolaiko and Kuchment (2013) G. Berkolaiko and P. Kuchment,  _Introduction to quantum graphs_, 186 (American Mathematical Soc., 2013).
  * Liew  _et al._ (2015) S. F. Liew, L. Ge, B. Redding, G. S. Solomon, and H. Cao, Phys. Rev. A 91, 043828 (2015).
  * Hisch  _et al._ (2017) T. Hisch, M. Liertzer, D. Pogany, F. Mintert, and S. Rotter, Phys. Rev. Lett. 111, 023902 (2017).
  * Peach  _et al._ (2021) R. L. Peach, A. Arnaudon, and M. Barahona, “Relative, local and global dimension in complex networks,”  (2021).
  * Ge  _et al._ (2008) L. Ge, R. J. Tandy, A. D. Stone, and H. E. Türeci, Optics Express 16, 16895 (2008).
  * Türeci  _et al._ (2008) H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, Science 320, 643 (2008).
  * Kuchment (2004) P. Kuchment, Waves in Random Media 14, S107 (2004).
  * Gnutzmann and Smilansky (2006) S. Gnutzmann and U. Smilansky, Advances in Physics 55, 527 (2006).
  * Arnaudon and Barahona (2022) A. Arnaudon and M. Barahona, In preparation (2022).
  * Liertzer  _et al._ (2012) M. Liertzer, L. Ge, A. Cerjan, A. Stone, H. E. Türeci, and S. Rotter, Physical Review Letters 108, 173901 (2012).
  * Esterhazy  _et al._ (2014) S. Esterhazy, D. Liu, M. Liertzer, A. Cerjan, L. Ge, K. Makris, A. Stone, J. Melenk, S. Johnson, and S. Rotter, Physical Review A 90, 023816 (2014).
  * Newton and Chamoun (2007) P. K. Newton and G. Chamoun, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 463, 1525 (2007).

## Supplementary Information

### NetSALT: Extension of SALT theory to lasing networks

We derive here netSALT, the numerical model used in the simulations of lasing
networks in this work. The accompanying code is available at
https://github.com/arnaudon/netSALT.

#### Open quantum graphs

A quantum graph is a metric graph (i.e., a graph with lengths $l_{ij}$
associated to each edge $(ij)$ and accompanying length variable
$x\in[0,l_{ij}]$) with a function $\eta(x)$ defined on each edge (thus defined
on the entire graph) that satisfies the Helmholtz differential equation

$\displaystyle\partial_{x}^{2}\eta_{ij}(x)+(n_{ij}k)^{2}\eta_{ij}(x)=0\qquad\forall
ij\,,$ (7)

where the complex numbers $n_{ij}$ correspond to the index of refraction of
the edge $(ij)$. This equation being linear, it has solutions of the form

$\displaystyle\eta_{ij}(x)=\lambda^{+}_{ij}e^{ikn_{ij}x}+\lambda^{-}_{ij}e^{ikn_{ij}(l_{ij}-x)}\,,$
(8)

where the complex-valued numbers $\lambda^{\pm}_{ij}$ represent the left- and
right-propagating wave amplitudes. The continuity of $\eta(x)$ at each node is
ensured by considering the edge function $\eta_{ij}(x)$ evaluated on the
nodes, such that

$\displaystyle\eta_{i}=\eta_{ij}(0)\quad\text{and}\quad\eta_{j}=\eta_{ij}(l_{ij})\,.$

The conservation of energy at each node $i$ can be shown to be equivalent to

$\displaystyle\left(L(k)\,\bm{\eta}\right)_{i}=\sum_{j\sim
i}n_{ij}\frac{\eta_{i}\left(e^{ikn_{ij}l_{ij}}+e^{-ikn_{ij}l_{ij}}\right)-2\eta_{j}}{e^{ikn_{ij}l_{ij}}-e^{-ikn_{ij}l_{ij}}}=0\,\qquad\forall
i\,,$ (9)

where the sum is over the nodes adjacent to $i$, and the matrix $L(k)$,
dependent on the wavenumber, is a node matrix acting on the node vector
$\bm{\eta}$ with components $\eta_{i}$, the value of $\eta(x)$ at each node.
We refer to Gnutzmann and Smilansky (2006); Berkolaiko and Kuchment (2013) for
more details on the derivation of this equation. The matrix $L(k)$ can be
expressed in terms of an extension of the graph incidence matrix, which allows
the simplification of the calculations of various quantities, see Arnaudon and
Barahona (2022) and below. The condition (9) corresponds to an eigenvalue
problem

$L(k)\,\bm{\eta}=\bm{0},$

so, equivalently one can solve the corresponding scalar equation

$\displaystyle\mathrm{det}\left(L(k)\right)=0\,,$ (10)

for discrete wavenumbers indexed as $k_{\mu}$. Numerically, we solve this
equation using the smallest eigenvalue of $L(k)$, which is efficient to
compute with sparse matrices. Notice that the node representation of quantum
graphs is not usual, as it contains a denominator term that diverges when
$kn_{ij}l_{ij}\to n\pi$, with $n=1,2,\ldots$, causing instabilities in the
numerical solution. This scenario happens in rare cases, when an edge has its
length divided by $\pi$ exactly proportional to the wavenumber. We only
encountered this issue for graphs with several same length edges, which is
fixed by adding a small noise on the edge lengths (or node positions).

Each edge with one open end (node of degree $1$) is considered to be outside
of the cavity, and admits no incoming wave. This is simply written as a
projection of the matrix $L$, where elements corresponding to the outgoing
waves are projected out, thus allowing them to take any value (and not
enforced to be vanishing from the right hand side of (9)). This condition
makes the quantum graph open, or lossy, and any solution of (9) must have a
complex wavenumber $k_{\mu}$. In the sequel, we will make the distinction
between the inner edges corresponding to the lasing cavity, and the outer
edges, corresponding to the open boundary of the cavity, also known as the
last scattering surface in laser theory. For example, we will use the
shorthand notation
$\int_{\mathrm{in}}dx=\sum_{ij\in\mathrm{in}}\int_{0}^{l_{ij}}dx$ for
integration over the inner edges of the cavity.

For each passive mode $k_{\mu}$, the standard $\mathcal{Q}$-value is given as

$\displaystyle\mathcal{Q}_{\mu}=\frac{\mathrm{Real}(k_{\mu})}{2|\mathrm{Im}(k_{\mu})|}\,.$
(11)

#### The SALT equation

The SALT equation Ge  _et al._ (2010) describes the interaction of lasing
modes under non-uniform pumping. The pump is described by an edge unit vector
$\delta_{\mathrm{pump}}$ of the cavity and an amplitude $D_{0}$.

In addition, the lasing modes, defined as modes with $\mathrm{Im}(k_{\mu})=0$,
are taken of the form $\Phi_{\mu}(x)=\sqrt{I_{\mu}}u_{\mu}(x)$ where
$\mathrm{I}_{\mu}$ is the modal intensity, and $u_{\mu}$ is the mode profile,
normalised as

$\displaystyle\int_{in}\delta_{\mathrm{pump}}u_{\mu}^{2}dx=1\,.$ (12)

On a single edge $(ij)$, the SALT equation Ge  _et al._ (2010) is a nonlinear
extension of the Helmholtz equation (7) given by

$\displaystyle\partial_{x}^{2}u_{\mu,ij}+\left(n_{ij}^{2}+\frac{D_{0}\delta_{\mathrm{pump},ij}\gamma_{\mu}}{1+\sum_{\nu}I_{\nu}\Gamma_{\nu}|u_{\nu,ij}|^{2}}\right)k_{\mu}^{2}u_{\mu,ij}=0\,,$
(13)

where $\gamma_{\mu}=\frac{\gamma_{\perp}}{k_{\mu}-k_{a}+i\gamma_{\perp}}$ is
the Lorentzian gain curve and $\Gamma_{\mu}=-\mathrm{Im}(\gamma_{\mu})$ the
gain linewidth.

We will not re-derive this equation here from several approximations of the
Maxwell-Bloch equation, but refer to Ge  _et al._ (2010) for more details and
only mention that the one of the main assumption of the SALT model is the
steady state assumption, or stationary inversion approximation, where the
inversion population (denoted by $D(x,t)$ in Ge  _et al._ (2008)) is taken to
be constant in time. We refer to Ge  _et al._ (2008); Liertzer  _et al._
(2012); Esterhazy  _et al._ (2014) for more detailed studies on the validity
and generalisations of this approximation it.

#### Finding threshold lasing modes

Before computing the modal amplitudes $I_{\mu}$, we need the threshold lasing
modes, solutions of the linear equation

$\displaystyle\partial_{x}^{2}u_{\mu,ij}+\left(n_{ij}^{2}+D_{\mu,th}\delta_{\mathrm{pump}}\gamma_{\mu}\right)k_{\mu,th}^{2}u_{\mu,ij}=0\quad\mathrm{where}\quad\mathrm{Im}(k_{\mu,th})=0\,.$
(14)

This is an implicit equation for the lasing threshold $D_{\mu,th}$, the
threshold wavenumber $k_{\mu,th}$ and the threshold lasing mode profile
$u_{\mu,ij}$.

Solving this equation must involves an iterative algorithm on the value of
$D_{0}$ to reach the condition $\mathrm{Im}(k_{\mu,th})=0$, where the secular
equation (10) is solved at each step. When $D_{0}$ is updated, use the so-
called Brownian Ratchet algorithm Newton and Chamoun (2007) to search for the
corresponding $k_{\mu}(D_{0})$. This algorithm consists in proposing random
moves in the complex plane of wavenumbers, and accepting only the ones
decreasing the smallest eigenvalue of $L(k)$, and stop the search when a
certain threshold is reached. The size of the proposed moves is adjusted
according to how far we expect the mode to have moved.

To speed up the search of threshold lasing modes, we first estimate the
location of a mode with a different $D_{0}$ by assuming that the mode profiles
do not change with pump, i.e. $u_{\mu,ij}(D_{0})=\eta_{\mu,ij}$. First, recall
that $\eta_{\mu,ij}$ are the passive modes, solution of

$\displaystyle\partial_{x}^{2}\eta_{\mu,ij}+n_{ij}^{2}k_{\mu,0}^{2}\eta_{\mu,ij}=0\,.$
(15)

Multiplying (14) by $\eta_{\mu,ij}$ and integrating over the cavity, we obtain

$\displaystyle
k_{\mu}(D_{0})=\frac{k_{\mu}(0)}{\sqrt{1+D_{0}\gamma_{\mu}f_{\mu,\mathrm{pump}}}}\,.$
(16)

where $f_{\mu}$ is the pump overlapping factor of mode $\mu$, defined as

$\displaystyle
f_{\mu}(\delta_{\mathrm{pump}})=\frac{\int_{\mathrm{in}}\delta_{\mathrm{pump}}\eta_{\mu,ij}^{2}dx}{\int_{\mathrm{in}}n^{2}_{ij}\eta_{\mu,ij}^{2}dx}\,.$
(17)

To estimate the the pump strength at threshold, we use
$\mathrm{Im}(k_{\mu,th})=0$ in (16) to get

$\displaystyle
D_{\mu,th}(\delta_{\mathrm{pump}})\approx-\frac{1}{\mathcal{Q}_{\mu}\Gamma_{\mu}\mathrm{Real}(f_{\mu}(\delta_{\mathrm{pump}}))}\,,$
(18)

where we also used the fact that $\mathrm{Real}(\gamma_{\mu})$ is small for
high $\mathcal{Q}$ modes.

Similarly, to obtain an estimation of the complex wavenumber for an updated of
pump power $D_{0}^{\prime}=D_{0}+\delta D_{0}$ from a mode with pump power
$D_{0}$, i.e. $(u_{\mu}(D_{0}),k_{\mu}(D_{0}))$, we use, instead of (16), the
equation

$\displaystyle
k_{\mu}(D_{0}^{\prime})=k_{\mu}(D_{0})\sqrt{\frac{1+D_{0}\gamma_{\mu}f_{\mu,\mathrm{pump}}}{1+D_{0}^{\prime}\gamma_{\mu}f_{\mu,\mathrm{pump}}}}\,,$
(19)

where $\gamma_{\mu}$ is now evaluated at $k_{\mu}(D_{0})$. This equation is
obtained similarly equation (16), by replacing the passive mode with a pumped
mode.

Hence, to find the threshold lasing modes, we linearly increase $D_{0}$ with
small steps, use (19) as a starting point for the Brownian ratchet algorithm
to find the next partially pumped mode, until we reach $Im(k)=0$, then use a
binary search (together with Brownian ratchet) to locate the exact position
(with some search threshold) of the lasing threshold $D_{\mu,th}$, lasing mode
wavenumber $k_{\mu,th}$.

#### Interacting modal intensities

Once the threshold lasing modes are found, we can estimate their modal
intensities as a function of the pump power $D_{0}$. For this, we assume that
the mode profiles above threshold are the same as the mode profiles at
threshold, and the threshold wavenumbers $k_{\mu,th}$ remain the same above
threshold. With these approximation, corresponding to the single pole
approximation of Ge  _et al._ (2010), we can estimate the modal intensities of
each mode, given a pump profile $\delta_{\mathrm{pump}}$ and a pump strength
$D_{0}$.

From (13) and using the normalisation (12), we follow Ge  _et al._ (2010) to
arrive at the matrix equation

$\displaystyle\sum_{\nu}T_{\mu\nu}I_{\nu}=\frac{D_{0}}{D_{th,\mu}}-1\,,$ (20)

where the sum is over lasing modes only, and the interaction matrix $T$ has
elements defined as

$\displaystyle
T_{\mu\nu}=\Gamma_{\nu}\mathrm{Real}\left(\int_{\mathrm{in}}|u_{\nu}|^{2}u_{\mu}^{2}\delta_{\mathrm{pump}}(x)dx\right)\,.$
(21)

Note that this matrix does not have an explicit dependence on the dielectric
constant. Notice that the real part is an approximation, as this quantity has
small complex part in general. Given $D_{0}$, the modal intensities are simply
found as

$\displaystyle
I_{\mu}(D_{0})=\sum_{\nu}T^{-1}_{\mu\nu}\left(\frac{D_{0}}{D_{th,\nu}}-1\right)\,,$
(22)

if the set of lasing modes (indexed as $\nu$) are known. To find the lasing
mode, we follow again Ge  _et al._ (2010), and first compute the interacting
lasing thresholds $D_{int,\mu}$. For the first lasing mode, the interaction
threshold will be the lasing threshold, but for the next lasing modes,
interaction with the currently lasing modes will increase this value, until it
reaches $\infty$, and no more modes can lase (called gain clamping).

To compute a lasing threshold mode, we assume that we have lasing $N$ modes,
and we seek to compute the interacting threshold of the next mode, indexed
$\mu_{N+1}$. At exactly $D_{0}=D_{int,\mu_{N+1}}$, the mode $\mu_{N+1}$ will
not lase, so $I_{\mu_{N+1}}=0$, which, after some manipulation, gives

$\displaystyle
D_{int,\mu_{N+1}}=D_{th,\mu_{N+1}}\left(1+\sum_{\mu_{i}=0}^{\mu_{N}}T_{\mu_{N+1}\mu_{i}}I_{\mu_{i}}\left(D_{int,\mu_{N+1}}\right)\right)\,,$
(23)

an implicit equation for the interacting lasing threshold. Being linear, we
can simply rearrange terms to get

$\displaystyle
D_{int,\mu_{N+1}}=D_{th,\mu_{N+1}}\frac{1-\sum_{i=0,j=0}^{N}T_{\mu_{N+1}\mu_{i}}T^{-1}_{\mu_{i}\mu_{j}}}{1-\sum_{i=0,j=0}^{N}\frac{D_{th,\mu_{N+1}}}{D_{th,\mu_{j}}}T_{\mu_{N+1}\mu_{i}}T^{-1}_{\mu_{i}\mu_{j}}}\,.$
(24)

The next lasing mode is therefore the mode $\mu_{N+1}$ with the smallest value
of $D_{int,\mu_{N+1}}$. At some point, the denominator will become negative,
corresponding to gain clamping regime, where all other modes are suppressed by
currently lasing modes, see Ge  _et al._ (2010) for more details on that.
Sometimes, a lasing mode can stop lasing, due to a negative slope in (22), in
which case, this mode is removed from the list of lasing modes and will not
contribute anymore to this equation for the search of the next lasing mode.

The solution of this equation thus provides the so-called LL curves, with
modal intensities of all the modes as a function of the pump power $D_{0}$, as
piece-wise linear functions, or approximation of lasing spectra at a given
pump power, if some artificial lasing linewidth are added.

#### Pump optimisation in NetSALT with linear programming

To numerically optimise the pump profile in netSALT, we cannot evaluate the
modal intensities, as this will result in a costly and slow algorithm.
Instead, we use the linear approximation of lasing threshold (17) relying on
the pump overlapping factor (18). Being linear with the pump profile, the
overlapping factor can be written as a scalar product
$f_{\mu}(\delta_{\mathrm{pump}})=\sum_{ij}\delta_{\mathrm{pump},ij}f_{\mu,ij}=\boldsymbol{\delta}\cdot\boldsymbol{f}_{\mu}$,
where $f_{\mu,ij}$ is the overlapping factor for a pump only defined in the
edge $(ij)$. The optimal pump $\widehat{\boldsymbol{\delta}}_{\mu}$ is then a
result of the minimisation problem of the form

$\displaystyle\widehat{\boldsymbol{\delta}}_{\mu}=\mathrm{argmin}_{\boldsymbol{\delta}}\frac{\mathrm{max}_{\nu}\,\boldsymbol{\delta}\cdot
a_{\nu}}{\boldsymbol{\delta}\cdot a_{\mu}}\,,$ (25)

where $a_{\nu}=\boldsymbol{f}_{\nu}\mathcal{Q}_{\nu}\Gamma_{\nu}$ This cost
function may lead to small pump profiles, thus we modify it by adding an extra
term parametrised by an hyper-parameter, or regulariser $\epsilon>0$,

$\displaystyle\widehat{\boldsymbol{\delta}}^{\epsilon}_{\mu}=\mathrm{argmin}_{\boldsymbol{\delta}}\frac{\mathrm{max}_{\nu}\,\boldsymbol{\delta}\cdot
a_{\nu}+\epsilon}{\boldsymbol{\delta}\cdot a_{\mu}}\,,$ (26)

which results in a family of solution with various coverage of the network
surface area on the target mode profile.

To solve this integer problem, we relax the integer-valued vector
$\boldsymbol{\delta}$ to a real vector $0<\boldsymbol{x}<1$ a solution of

$\displaystyle\widehat{x}_{\mu}^{\epsilon}=\min_{x}\frac{\max_{\nu}{a_{\nu}^{T}x}+\epsilon}{a_{\mu}^{T}x}\,.$
(27)

From this solution, we propose a pump with all edges such that $x_{ij}>0$,
then remove edges which have a small impact on the cost, considered as noise
from the SALT approximation, which may potentially reduce the resulting modal
suppression ratio.

To solve the relaxed problem with $\boldsymbol{x}$, we rewrite it as a linear
program by adding an additional variable $m$ to represent the maximum in the
numerator and by using the Charnes-Cooper transformation

$\displaystyle y=\frac{x}{a_{\mu}^{T}x},\qquad t=\frac{1}{a_{\mu}^{T}x}\,.$
(28)

The corresponding linear program is

$\displaystyle\begin{split}\min_{y,m,t}&\,m+\epsilon t\\\ a_{\nu}^{T}y&\leq
m,\quad\forall\nu\\\ a_{\mu}^{T}y&=1\\\ 0<y_{i}&<t,\quad
i=0,\ldots,n\end{split}$ (29)

for which the solution of the original problem is given as $x=\frac{1}{t}m$.
We solve this linear problem using the public python software PuLP, available
at https://github.com/coin-or/pulp. The results of this optimisation on the
Buffon graph are illustrated in Fig. S5.

To optimise a pump for multi-mode lasing (see Fig. S8) in this linear
programming framework, we replace the denominator of the cost by the sum over
the $a_{\mu}\rightarrow\sum_{\mu}^{\prime}a_{\mu}^{\prime}$ of $\mu^{\prime}$
mode we which to lase together.

####  .0.1 Comparison with mode matching optimisation

A simpler strategy to optimise the pump profile for single lasing a particular
mode would be to assume that only pumping edges with large electric field of
the target mode will work.

We apply this method by selecting edges supporting the largest amplitudes of
the mode we want to single lase such that the cost function defined above is
minimised (taken with $\epsilon=0$). Then, as for the optimisation, we remove
edges from the pump which have a small impact on the cost. The results of this
optimisation on the Buffon graph are illustrated in Fig. S5, and shown to
produce small pump if a mode has only a few edges with most of its electric
field amplitude, or large one for highly delocalised modes. Globally, it is
outperformed by the optimised pump with linear programming, but sometimes
result in better suppression ratio, when the linear approximation used in the
cost function is not representative enough of the modal amplitudes.

### Classical laser geometries

Here we validate our netSALT calculations by modelling some simple laser
cavities.

#### 1D-cavity laser

We model a 1D-cavity laser with non-uniform index and non-uniform pump profile
from Ge  _et al._ (2010), as shown in Fig. S11a. Optical feedback due to
reflection at the two ends of the cavity (positions $x=0$ and $x=1$) is taken
into account by adding edges with unit index of refraction to the line graph,
and imposing open boundary conditions at the outer nodes (in red). The index
of refraction is set to $1.5$ on the left $1/4$ of the cavity and set to $3$
on the remaining inner edges. The pump is applied to the left half of the
cavity (on inner edges shown in green).

The remaining panels Fig. S11 reproduce the results of Ref. Ge  _et al._
(2010). The mode profiles of the first lasing mode (Fig. S11b) matches exactly
with Fig 3a and modal intensities (Fig. S11c) matches with Fig. 6 of Ref Ge
_et al._ (2010), respectively. The threshold lasing frequencies and the non-
interacting lasing thresholds (not shown) also match the values reported in
the paper. This example is available in the github repository.

#### Ring laser

In a ring with real index $n$ and length $L$, the modes lie on the real axis
and are given by $k_{m}=\frac{2m\pi}{nL}$, where $m$ is a positive integer. To
model a ring laser with a finite $Q$ factor, we require complex refractive
index $\tilde{n}$ on the edges. Let $\tilde{n}=n+i\kappa$, then the modes are
given by complex values:

$\displaystyle\text{Re}(k_{m})$
$\displaystyle=\frac{n}{n^{2}+\kappa^{2}}\frac{2m\pi}{L}\qquad\text{Im}(k_{m})=-\frac{\kappa}{n}\text{Re}(k)$

Loss is therefore defined via $\kappa$, or equivalently by the $\mathcal{Q}$
factor (Eqn. (11)). This example is reproducible in the github repository.
Fig. S9 shows example netSALT calculation of a uniformly pumped micro-ring
laser with cavity length $10$ $\mu$m and refractive index $1.5+0.005i$.

### Additional calculations in NetSALT

We collect here additional formulae of the netSALT model described above.

#### Matrix representation of $L(k)$

The quantum graph equation $L(k)\eta=0$ with the wavenumber dependent matrix
can be written in term of matrices with analogues in classical graph theory.
Indeed, this matrix can be interpreted as a quantum graph Laplacian of the
form

$\displaystyle L(k)\eta:=B^{T}(k)W^{-1}(k)B(k)\eta\,,$

with the matrix $B$ an extension of the incidence matrix with elements of the
form

$\displaystyle\begin{split}B_{i,ij}&=-1\\\
B_{j,ij}&=e^{ikl_{ij}}\,,\end{split}$

and the diagonal weight matrix $W$ is defined as

$\displaystyle W_{ij,ij}=W_{ji,ji}=e^{2ikl_{ij}}-1\,.$

We refer to Arnaudon and Barahona (2022) for the details on the derivation of
these equations.

#### Pump overlapping factor

The pump overlapping factor defined in (17) is explicitly given as

$\displaystyle f_{I}=\sum_{ij\in
I_{ij}}(\delta_{pump})_{ij}\int\eta_{\mu,ij(x)}^{2}dx=\eta^{T}(B^{d})^{T}W^{-1}\delta_{pump}ZW^{-1}B\eta\,,$

where the matrix $Z$ is

$\displaystyle
Z_{ij,ij}=\begin{pmatrix}\frac{e^{2ikl_{ij}}-1}{2ik}&l_{ij}e^{ikl_{ij}}\\\
l_{ij}e^{ikl_{ij}}&\frac{e^{2ikl_{ij}}-1}{2ik}\end{pmatrix}\,.$

#### Mode competition matrix

From the simplicity of this calculation, we consider $k_{\mu}$ to be complex
and contain the index of refraction and the pump term with $\gamma$, and we
drop all the edge indices $ij$. The matrix $\mathcal{T}_{\mu\nu}$

$\displaystyle
T_{\mu\nu}=\Gamma_{\nu}\mathrm{Real}\left(\frac{\int_{\mathrm{in}}|u_{\nu}|^{2}u_{\mu}^{2}\delta_{\mathrm{pump}}(x)dx}{\int_{\mathrm{in}}u_{\mu}^{2}\delta_{\mathrm{pump}}dx}\right)\,,$

has the following elements in its numerator

$\displaystyle\int_{\mathrm{in}}|u_{\nu}|^{2}u_{\mu}^{2}\delta_{\mathrm{pump}}(x)dx=\sum_{\mathbf{edges}\in\mathbf{pump}}\begin{pmatrix}|\lambda_{\nu}^{+}|^{2}\\\
\lambda_{\nu}^{+}\overline{\lambda}_{\nu}^{-}\\\
\overline{\lambda}_{\nu}^{+}\lambda_{\nu}^{-}\\\ |\lambda_{\nu}^{-}|^{2}\\\
\end{pmatrix}^{T}\begin{pmatrix}A&E&E&B\\\ C&F&F&D\\\ D&F&F&C\\\ B&E&E&A\\\
\end{pmatrix}\begin{pmatrix}(\lambda_{\mu}^{+})^{2}\\\
\lambda_{\mu}^{+}\lambda_{\mu}^{-}\\\ \lambda_{\mu}^{+}\lambda_{\mu}^{-}\\\
(\lambda_{\mu}^{-})^{2}\\\ \end{pmatrix}$

where

$\displaystyle A$
$\displaystyle=\frac{e^{i(k_{\nu}-\overline{k}_{\nu}+2k_{\mu})l}-1}{i(k_{\nu}-\overline{k}_{\nu}+2k_{\mu})}$
$\displaystyle B$
$\displaystyle=e^{2ik_{\mu}l}\frac{e^{i(k_{\nu}-\overline{k}_{\nu}-2k_{\mu})l}-1}{i(k_{\nu}-\overline{k}_{\nu}-2k_{\mu})}$
$\displaystyle C$
$\displaystyle=\frac{e^{i(k_{\nu}+2k_{\mu})l}-e^{-i\overline{k}_{\nu}l}}{i(k_{\nu}+\overline{k}_{\nu}+2k_{\mu})}$
$\displaystyle D$
$\displaystyle=\frac{e^{ik_{\nu}l}-e^{i(2k_{\mu}-\overline{k}_{\nu})l}}{i(k_{\nu}+\overline{k}_{\nu}-2k_{\mu})}$
$\displaystyle E$
$\displaystyle=e^{ik_{\mu}l}\frac{e^{i(k_{\nu}-\overline{k}_{\nu})l}-1}{i(k_{\nu}-\overline{k}_{\nu})}$
$\displaystyle F$
$\displaystyle=e^{ik_{\mu}l}\frac{e^{ik_{\nu}l}-e^{-i\overline{k}_{\nu}l}}{i(k_{\nu}+\overline{k}_{\nu})}\,,$

and the sum over the edges in the pump uses the dummy indices $ij$.

#### Edge mean of $|E|^{2}$ calculation

The mode solution has the form:
$\eta_{ij}(x)=\lambda^{+}_{ij}e^{ikn_{ij}x}+\lambda^{-}_{ij}e^{ikn_{ij}(l_{ij}-x)}$,
which has the value $\lambda^{-}_{ij}e^{ikn_{ij}l_{ij}}$ at $x=0$ and
$\lambda^{+}_{ij}e^{ikn_{ij}l_{ij}}$ at $x=l_{ij}$.

For brevity, we remove subscript and take the modulus squared:

$\displaystyle\eta(x)\overline{\eta(x)}=\lambda^{+}\overline{\lambda^{+}}e^{(ikn+\overline{ikn})x}$
$\displaystyle+\lambda^{+}\overline{\lambda^{-}}e^{\overline{ikn}l}e^{(ikn-\overline{ikn})x}$
$\displaystyle+\lambda^{-}\overline{\lambda^{+}}e^{iknl}e^{(\overline{ikn}-ikn)x}+\lambda^{-}\overline{\lambda^{-}}e^{(ikn+\overline{ikn})(l-x)}$

If we integrate this from $x=0$ to $l$, we get

$\displaystyle\langle|E|^{2}\rangle=\frac{1}{l}\int\eta(x)\overline{\eta(x)}dx=\frac{1}{l}\lambda^{+}\overline{\lambda^{+}}\frac{e^{(ikn+\overline{ikn})l}-1}{ikn+\overline{ikn}}$
$\displaystyle+\frac{1}{l}\lambda^{+}\overline{\lambda^{-}}\frac{e^{iknl}-e^{\overline{ikn}l}}{ikn-\overline{ikn}}$
$\displaystyle+\frac{1}{l}\lambda^{-}\overline{\lambda^{+}}\frac{e^{\overline{ikn}l}-e^{iknl}}{\overline{ikn}-ikn}+\frac{1}{l}\lambda^{-}\overline{\lambda^{-}}\frac{e^{(ikn+\overline{ikn})l}-1}{ikn+\overline{ikn}}$

This can be expressed in matrix form as

$\displaystyle\langle|E|^{2}\rangle=\frac{1}{l}\begin{pmatrix}\lambda^{+}\\\
\lambda^{-}\end{pmatrix}\begin{pmatrix}\frac{e^{(ikn+\overline{ikn})l}-1}{ikn+\overline{ikn}}&\frac{e^{iknl}-e^{\overline{ikn}l}}{ikn-\overline{ikn}}\\\
\frac{e^{iknl}-e^{\overline{ikn}l}}{ikn-\overline{ikn}}&\frac{e^{(ikn+\overline{ikn})l}-1}{ikn+\overline{ikn}}\end{pmatrix}\begin{pmatrix}\overline{\lambda^{+}}&\overline{\lambda^{-}}\end{pmatrix}$

which can be computed from the node solution.

#### Inverse participation ratio calculation

The inverse participation ration ($IPR$) provides a measure for the mode
spread over the graph, and is given by:

$\displaystyle
IPR_{\mu}=L_{tot}\frac{\sum_{ij}\int_{0}^{l_{ij}}|E_{\mu}|^{4}dx}{(\sum_{ij}\int_{0}^{l_{ij}}|E_{\mu}|^{2}dx)^{2}}\,,$
(30)

which can be evaluated analytically using the complex wave amplitudes on the
edges and the analytical solution of the electric field on each edge.

## Supplementary Figures

Figure S1: Sensitivity of network lasing spectrum. a Experimental spectrum
duplicated from Fig. 1b of main text. b Ratio of intensities calculated from
the spectra in a. When the pump pattern is modified by removing the pump from
small central area (which results in a reduction of the delivered pump power
by 7 %), some lasing peaks increase by a factor 2.8 while others are
attenuated down to 0.2 of the initial intensity.  Figure S2: Missing edges
experiment. a We show the modal amplitudes of uniform vs missing edge pump
profiles, the colour corresponds to the overlap with 2 missing edges, computed
as the real part of the $f$ factor with a pump localised on these two edges
only. Modes with a large overlap with these edges are suppressed, such as the
target mode $1$ in red, and others appear, such as the mode $6$ in green. b We
show that overlap with missing edges (same quantity as the colour in panel a
against the difference in lasing thresholds of all modes (black), lasing modes
(blue) and target mode $1$ (red) and new mode $6$. We observe a strong
correlation between these two quantities, indicating that the $f$ factor of
each mode on edges is a good indicator of resulting changes in non-interacting
lasing thresholds. c Number of lasing modes as a function of pump power for
uniform and missing edges pump profiles. d-e We display the profile of modes
$1$ (in d) and $6$ (in e with missing edges indicated by arrows. We see a
visual large overlap between these modes, but a small overlap of mode $6$ with
the missing edges (thus a small difference in lasing threshold, as shown in in
panel c.  Figure S3: Mode trajectories under uniform and optimised pumping.
a-b. Uniform pump profile and optimised pump profile for mode $1$, as also
shown in main manuscript Fig. 2d insets. Dark green segments indicate region
on edges that are pumped. c. Mode trajectories in the complex $k$ plane with
the uniform (top) and optimised (bottom) pump profile with increasing pump
power $D_{0}$, values shown by colour bar. The trajectories are linear for all
modes with uniform pumping. While gain in the network amplifies many modes to
reach threshold Im$(k)=0$, only a subset of modes (indicated by purple
diamonds) actually lase due to mode competition and gain saturation. The
passive modes (without gain) are shown by red dots and modes are labelled in
descending order of $Q$-factor.  Figure S4: Zoom in of LL curves shown in Fig
2c, where the lasing threshold of mode 1 becomes the smallest and the
difference between mode 1 threshold and the next lasing threshold becomes
large. Figure S5:  Mode control with optimisation and mode matching. a
Interacting lasing thresholds vs. lasing threshold, showing a large range of
amount of modal interactions. b Overlap with pump ($f$ factor) vs. lasing
threshold difference with respect to uniform pumping, showing that low lasing
thresholds modes have large overlap with the pump (each dot is a mode for $20$
optimisations). c Modal ratio vs $\mathcal{Q}$-factor for each of $100$
optimisations (blue) and mode matching (black) pumps. Crosses are for pumps
with small surface area (smaller than $2\%$ of the surface area). The
optimisation and mode matching were not limited in surface area profiles, and
mode matching works well for small pumps on localised modes (with not other
modes lasing) but is not realistic experimentally. Overall, the optimisation
provides pumps with larger areas and better ratios. With optimisation, 66/100
have ratio larger than 2, and 90/100 larger than 1, while for mode matching,
39/100 have larger than 2 and 59/100 larger than one. d Same as c, but with
$IPR$ instead of $Q$-factor, showing no dependence on mode localisation. e
Fraction of pump area for optimisation and mode matching, showing that mode
matching either under or over estimates the pump area. f Distribution of
suppression ratio for all pumps in dashed, and large pumps (larger than $2\%$
of surface area) in thick lines. g Full controllability matrix, where each row
contains the modal amplitudes of a linear programming optimisation. h Same as
g but for mode matching optimisation, showing larger off diagonal values, thus
worse single lasing regimes than with the optimisation.  Figure S6:
Optimisation of a high threshold mode. Optimisation performed on the
$11^{\textrm{th}}$ strongest intensity mode from the uniform pumped spectrum.
a. Optical image of the sample when illuminated with the optimised pattern.
Data is taken from the same area of the sample as in Fig. 3 of the main text.
b. Normalised spectrum at full power with uniform pump (grey) and optimised
pump pattern (purple). c-d. Light in-light out (LL) curve showing intensity of
lasing modes with pump power with uniform and optimised pump profile,
respectively.  Figure S7: Spectral modulation using optimised pump patterns.
Lasing spectrum from the network laser is stable and can be illuminated with
different optimised pump patterns to alternate between different lasing modes.
An example of alternating 10 times with two single mode lasing spectra in (a)
and between single mode and three mode lasing spectra in (b). Respective pump
patterns are shown on top right insets. Data in a and b was taken from a
different region of the sample compared to Fig. 3 of the main text.  Figure
S8: Multi-mode optimisation with netSALT a Mode profile of the four modes
targeted to lase, where mode $1$ and $7$ with the most overlap. b Optimised
pump profiles to lase 2, 3 or 4 of these modes. c Modal intensities as a
function of pump power for these three multi-mode optimisations, showing a
faster increase of the targeted mode, despite the modal competitions. d
Synthetic spectra of the three optimisations at $D_{0}=0.02$.  Figure S9:
Example of a ring laser. a. Ring graph with $L=10~{}\mu\textrm{m}$ and uniform
index $n=1.5+0.005i$. b. Electric field intensity profile of the lasing mode
(axial order $m=36$) at threshold $D_{0,th}=0.015$. c. Lasing spectra at
$D_{0}=0.05$ showing intensities of the 5 lasing modes. Gain spectrum (red
dashed line) has the following parameters: $k_{a}=15$ and $\gamma_{\perp}=3$.
Bottom panel shows the mode trajectories of the passive modes (open circles)
at different $D_{0}$ values: $D_{0}=0.005$, $0.01$, $0.015$ (black circles).
d. LL curve for the ring laser. Inset shows zoomed in view of the first three
lasing modes and their non-interacting thresholds (vertical dashed lines).
Figure S10: Periodic vs random network. a-b Controllability map for a periodic
network, a, in the form of a grid, which show poor control of the lasing
action as compared to a similar size network with a random topology, b, which
instead can be well controlled. Two representative modes are shown in figure c
and d, respectively.  Figure S11: Example of 1D laser cavity. a Schematic of
1D laser cavity with non-uniform index profile and non-uniform pump-profile.
Graph representation of the example is shown below. b Normalised electric
field intensity profile of the lasing mode at threshold (blue) and without
pumping (orange dashed). c Lasing intensity as a function of pump intensity
$D_{0}$ and spectrum at $D_{0}=1.2$ in the inset.