# Augmenting On-Chip Microresonator through Photonic Inverse Design

Geun Ho Ahn1, Ki Youl Yang1, Rahul Trivedi1,2, Alexander D. White1, Logan Su1,
Jinhie Skarda1, Jelena Vučković1,†
1E. L. Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA.
2Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Str. 1, 85748 Garching,
Germany.
†Corresponding author<EMAIL_ADDRESS>
###### Abstract

Recent advances in the design and fabrication of on-chip optical
microresonators has greatly expanded their applications in photonics, enabling
metrology, communications, and on-chip lasers. Designs for these applications
require fine control of dispersion, bandwidth and high optical quality
factors. Co-engineering these figures of merit remains a significant
technological challenge due to design strategies being largely limited to
analytical tuning of cross-sectional geometry. Here, we show that photonic
inverse-design facilitates and expands the functionality of on-chip
microresonators; theoretically and experimentally demonstrating flexible
dispersion engineering, quality factor beyond 2 million on the silicon-on-
insulator platform with single mode operation, and selective wavelength-band
operation.

On-chip optical microresonators are critical to a wide range of applications,
including but not limited to lasers, optical communication, nonlinear optics,
and quantum opticsVahala:2003:Nature ; Lipson:2005:Nature ; Englund:2005:PRL .
While relying on the same basic operational principles, resonator designs for
different technological applications require vastly different design
considerations. In particular, several important metrics, including
dispersionOkawachi:2014:OL ; Yang:2016:NaturePhotonics ;
Weiner:2017:NatureCommunications ; Papp:2019:ACSPhotonics ;
Bowers:2020:PhotonicsResearch , quality factorBowers:2021:NaturePhotonics ;
Blumenthal:2021:NatureCommunications ; Yang:2018:NaturePhotonics ;
Kippenberg:2018:Optica  and operation bandwidthTang:2018:Optica ;
Loncar:2019:Optica ; Lipson:2015:NatureCommunications , need to be co-
optimized with the specific application dictating their relative importance.
Current microresonator design methodologies mostly rely on tuning a small set
of design parameters, such as the resonator geometry or the material stack,
guided by an analytical or numerical modelling of the device performance.
Consequently, it is a challenge to hit multiple design targets.

An alternative strategy is to use _inverse design_ — here, the designer first
constructs a cost function, expressed in terms of the electric fields
generated by the device under different excitation conditions, which captures
the required design targets. A gradient based optimization algorithm can then
be used to efficiently explore a large space of devices to optimize this cost
function. Application of this methodology to designing traditional integrated
optical devices, such as wavelength demultiplexers
Piggott:2015:NaturePhotonics , mode converters inv_mode_2016 , grating
couplers Su:18  has resulted in devices that achieve high performance, high
bandwidth and low cross-talk within a much smaller device footprint when
compared to their traditionally designed counterparts. Here, we consider the
application of inverse design for on-chip microresonator design — we
numerically and experimentally demonstrate this with three applications (i)
engineering the microresonator resonances (dispersion engineering), (ii)
achieving high quality (Q) factors in multi-mode waveguide-based resonators
and (iii) designing wavelength-band selective resonances in microresonators.

The first design problem that we address is dispersion engineering — several
applications in nonlinear optics rely on finely tuned microresonator
dispersions. For instance, the dispersion of the microresonator needs to be
finely controlled for frequency comb applications, both to phase match for the
initiation of optical parametric oscillation Kippenberg:2018:Science , as well
as to control the comb’s spectral shape and bandwidth Kippenberg:2018:Science
; Loncar:2019:Nature ; Afltouni:2020:CLEO . Traditional dispersion engineering
approaches attempt to design the frequency-dependent phase accumulated by an
optical field propagating in the microresonator by either tuning the resonator
cross-sectional geometryOkawachi:2014:OL ; Yang:2016:NaturePhotonics ,
carefully selecting a material stackKippenberg:2012:OE  or utilizing grating
structuresPapp:2019:ACSPhotonics ; Bowers:2020:PhotonicsResearch . While these
designs are able to initiate the nonlinear wave mixing processes, it is
difficult for them to achieve finer spectral control over the nonlinear
process Bell:Natpho . Recently, more sophisticated design methods that couple
multiple resonator modes have been employed to expand this design space —
these include photonic molecules Loncar:2019:NaturePhotonics ;
Weiner:2017:NatureCommunications ;
Kippenberg:2018:NatureCommunications:visiblecomb ;
Vahala:2017:NatureCommun:visiblecomb  and mode hybridization effects at
photonic crystal edges Yu:2021:NaturePhotonics ; Yu:2020:PhotonicResearch . By
posing the design of either a frequency-dependent phase or the coupling matrix
between multiple resonator modes as an inverse design problem, we numerically
and experimentally demonstrate a wide-range of dispersion profiles that are
relevant for nonlinear optics applications.

Figure 1: Inverse designed intracavity element for microring resonator
frequency engineering (a) SEM of fabricated Si microring resonator with the
partially reflective element in the microring. The microring resonator is
designed with a ring radius of 100 $\mu m$ and the width of the waveguide is
500 $nm$. The scale bar here is 20$\mu m$ and 500$nm$ respectively. (b)
Transmission spectrum of the fabricated device with the inverse designed
intracavity element. (c) simulated reflectance spectrum of the corresponding
inverse designed intracavity element with increasing reflectance as frequency
increases. (d) Measured integrated dispersion of the split fundamental modes,
indicated by blue and red points for each side of the split modes.

Along with achieving fine dispersion control, it is also critical in many
applications for microresonators to have a high quality (Q) factor – the next
design challenge we address. A high Q factor is crucial for the resonant
enhancement of lasingBowers:2021:NaturePhotonics , nonlinear wave mixing
Srinivasan:2019:NaturePhotonics , light-matter interactionDirk_PhysRevLett
and narrow band filters for on-chip communicationWatts:2012:OL ;
Watts:2014:NatureCommunications ; Ram:2018:Nature . Apart from material
absorption, the Q-factor of the resonator modes are limited by the
microresonator surface roughness. Consequently, one of the key strategies for
improving the quality factor is to optimize the fabrication process to improve
surface roughness to reduce the scattering
lossesLipson:2015:NatureCommunications ; Watts:2012:OL . Alternatively, the
interaction of the resonator mode with the rough resonator surface can be
reduced by using wider-width waveguides to form the
microresonatorWatts:2012:OL . However, wider-width waveguides in general
support multiple propagating modes, which can couple to the fundamental
resonator mode and adversely impact its Q-factor. An adiabatic taper can be
used to suppress higher order modes Kordts:16 , but while this is an effective
way of achieving high quality factor single mode operation, the adiabatic
taper makes the resonator very long, limiting applications to those that can
tolerate a small free-spectral range (FSR). In this paper, we demonstrate that
using inverse design, the coupling of the fundamental resonator mode to higher
order waveguide modes can be prohibited within a much shorter device area,
allowing for flexible engineering of the microresonator FSR and Q-factor.

Finally, we consider the problem of designing microresonators that support
resonant modes at two distinct and well-separated resonant bands. Such
wavelength selective band operation is critical for nonlinear interaction
between optical wavelengths that are far apart while preventing leakage into
undesired bands mckenna2021ultralowpower  which could significantly reduce the
efficiency of the nonlinear process. Traditional resonator designs find it
difficult to implement such frequency-selective resonance formations while
also maintaining high Q-factor and dispersion control. We resolve this problem
by posing the frequency selectivity of the resonator as an inverse design
problem, and numerically and experimentally demonstrate wavelength selective
band operation.

Inverse Design of Microresonator Dispersion.

Figure 2: Inverse designed dispersion engineered Fabry-Perot resonators (a)
Schematic of inverse designed FP resonator and conceptual depiction of
dispersion compensation from the inverse designed reflectors, analogous to
chirped mirrors. (b) Measured integrated dispersion of three inverse designed
FP resonators with cavity length of 200 $\mu m$ with inverse designed
reflectors with dispersion targets. Blue, red, and black lines correspond to
fit lines for the three FP resonators, showing anomalous, near-zero, and
normal dispersion. (c) SEM images of the three different inverse designed
reflectors used for the dispersion engineered resonators. (d) OPO simulation
of the FP resonator with the experimental dispersion coefficients $D_{2}$ and
$D_{3}$ derived from dispersion measurements of the three FP resonators with
cavity length of 480 $\mu m$

Due to their cylindrical symmetry, whispering gallery mode resonators possess
degenerate resonant modes which propagate in the clockwise and counter-
clockwise directions. By systematically coupling these two modes, we can
control the positions of their resonances and engineer the desired dispersion
profile. Here, we introduce a computational method to achieve this through
inverse designed intracavity partially reflective elements
(PREs)Yang:2020:NaturePhotonics .

When a PRE is introduced in the resonator, it couples the counter propagating
resonant modes, lifting the mode degeneracy. The coupling strength, which is
dictated by the reflectivity of the element, in turn determines the position
of resonances of the resulting hybridized modes. More specifically, as we show
in supplementary section I, the resonance splitting between the hybridized
modes is given by $\Gamma=\theta/\pi\times\text{FSR}$, where FSR is the free-
spectral range of the resonator and
$\theta=\arctan\left({|r|}/{\sqrt{1-|r|^{2}}}\right)$, with $|r|$ being the
reflectance of the PRE at the resonant frequency of the unhybridized modes.
Consequently, a desired microring dispersion, specified as the spectral
position of resonant modes, can be translated to a target frequency-dependent
reflectance of the PRE, which can then be inverse designed. As a
demonstration, we inverse design a PRE that has higher reflectance at higher
frequency and lower reflectance at lower frequency. FIG 1a shows the scanning
electron microscope (SEM) image of the inverse design PRE embedded in a
microring resonator fabricated on a SOI stack. FIG 1b shows the transmission
spectrum of the microring resonator when excited through the coupled bus
waveguide – we indeed observe that the splitting between the hybridized
microring resonator modes, and hence their resonant position, is dictated by
the engineered reflectance spectrum of the PRE (FIG 1c). The measured
integrated microring resonator dispersion ($D_{int}$) for the red-shifted and
blue-shifted hybridized modes is depicted in FIG 1d — we clearly observe a
correlation between the extent of mode splitting and the reflectance spectrum.

Figure 3: Inverse design of high Q FP resonator with single mode operation.
(a) Schematic of the inverse designed FP resonator with inverse designed
multimode reflectors allowing for a wide width waveguide for the resonator and
single mode waveguide outputs. Conceptual picture showing the functionality of
reflector, as it reflects fundamental mode well, and can couple out
fundamental mode, while other higher order modes are not supported. (b) SEM
image of the single mode inverse designed reflector with wider width waveguide
(1.5 $\mu m$) for reflection of fundamental mode. Single mode waveguide of 0.5
$\mu m$ width is used on the other side to couple out the resonances from the
resonator. The scale bar is 2 $\mu m$ (c) Simulated reflectance of the
depicted inverse design reflector when excited with different TE modes,
showing its selective reflection of the fundamental mode only. (d) Zoomed in
transmission spectrum through the FP resonator, showing single mode operation
separated by a FSR of 19.8 GHz. (e) Q-factor fit of a resonance of the inverse
designed FP resonator, showing a loaded Q-factor of 2.03M. (f) Q factor
statistics of resonance distribution from three different fabrication runs of
the identical FP resonator design.

A similar design technique can be applied to control the dispersion of on-chip
fabry-perot (FP) resonators, which are composed of two reflectors with a
straight waveguide section in the middle. This length of waveguide $L_{D}$,
together with the group index $n_{g}$ of the propagating waveguide mode
controls the FSR of the microresonator, as FSR = $c/2n_{g}L_{D}$.
Traditionally, the dispersion of the FP resonator is controlled by engineering
the waveguide cross-section and the length of the
waveguidePapp:2019:ACSPhotonics , or by employing chirped mirrors Chirped1994
. To achieve finer spectral control over the dispersion of the resonator, we
can alternatively use inverse-design to engineer the frequency-dependent phase
imparted by the mirrors at the ends of the FP resonator (FIG 2a). Such
resonator with inverse designed reflectors is schematically shown in Fig. 2a.
Notably, this concept of dispersion engineering using optimized reflectors
with different dispersion targets is a generalization of chirped mirrors used
in pulse engineering in lasers, where the distribution of multilayer coatings
of thin films is optimized to provide engineered group delay dispersion
Chirped1994 . As demonstrative examples, we use inverse designed reflectors to
make the FP resonator dispersion anomalous, where group velocity dispersion is
positive, or be normal, where group velocity dispersion is negative. In
addition, the optimized reflectors can be very compact, with their width and
length being only 0.5 $\mu m$ and 3.0 $\mu m$ respectively. The optimized
dispersion engineered reflectors work by compensating phase upon the
reflection. Notably, the inverse design technique finds the optimal structure
that provides the desired wavelength dependent phase compensation in a compact
footprint manner. Fig. 2b illustrates the experimental integrated dispersion
profiles of the microresonators formed by three different inverse designed
reflectors with different dispersion target, as shown in the SEM images of the
reflectors in Fig. 2c. The cavity length used is 200 $\mu m$. We demonstrate
tunable dispersion in equal length cavities ranging from anomalous to near-
zero to normal, with $D_{2}$ values ranging from 18 MHz to -168 MHz. In
addition, the length of the cavities can be varied, achieving even more
flexible dispersion tuning, as shown in FIG S2. To explore the potential
implications of such dispersion engineered microresonators, we simulate
optical parametric oscillation (OPO) based on the experimentally derived
dispersion coefficients. Detailed implementation of the OPO simulation of a
Fabry-Perot resonator with two photon absorption of silicon is shown in
supplementary section II. Fig. 2d shows the result of OPO simulations with
three different dispersion profiles derived from the inverse designed
resonators with a 480 $\mu m$ cavity length, as shown in FIG S2. Importantly,
when the dispersion of the microresonator is engineered to be normal, we do
not see any OPO side band generation, as expected. When the dispersion of the
microresonator is engineered to be anomalous at different strengths, the
microresonators show OPOs generated at different resonance modes.

High Quality Factor Single Mode Resonator.

Figure 4: Inverse design of selective wavelength band resonator (a) A
graphical representation of parametric coupling in a typical resonator and in
a selective band resonator, formed by selective reflection bands. (b)
Schematic of utilized inverse designed selective wavelength reflector and its
simulated electric field profile at the chosen wavelength of operation ($1300$
$nm$ and $1550$ $nm$). (c) Simulated reflectance profile of the inverse
designed selective reflection band reflector, showing high reflectance at
O-band and C-band with low reflectance in other wavelength bands. (d) Measured
transmission spectrum of the selective wavelength band resonator, showing
sharp resonances in both O-band and C-band, while the resonances outside these
bands are broad.

To achieve high microresonator quality factors, the key is reducing the
interaction of the optical fields with the microresonator surface, thereby
reducing scattering losses due to surface roughness. As previously mentioned,
this can often be achieved by using wide-width waveguides which tend to
support fundamental waveguide modes that are well confined within the
waveguide. However, wider waveguides often support multiple waveguide modes,
which can couple to the fundamental mode and introduce loss.

To demonstrate how inverse design can be used to support this goal, we design
a low loss FP resonator built out of a wide-width waveguide that supports only
a single resonant mode. This is achieved by inverse designing reflectors that
selectively reflect the fundamental mode of the wide waveguide and couple the
transmitted field directly to single mode waveguides, suppressing higher-order
modes in the process (FIG 3a). FIG 3b shows the SEM of the inverse designed
reflector fabricated on a SOI stack, and FIG 3c shows its simulated
reflectance when excited with different waveguide modes. We clearly see that
the reflector only acts as a near-perfect mirror for the fundamental waveguide
mode, while the remaining waveguide modes are allowed to transmit through and
hence do not interfere with the fundamental mode operation. Furthermore, this
functionality is achieved within a footprint of 3.5 $\mu m$ by 1.5 $\mu m$,
which is significantly smaller than the conventional tapers that can be used
to accomplish the same task. FIG 3d shows the transmission spectrum of the FP
resonator and FIG 3e shows the zoomed-in spectrum of a single resonator mode —
we experimentally achieve a loaded Q-factor of 2.03 million at
telecommunication wavelengths in the foundry-compatible SOI system with our
design. Furthermore, the design is robust to fabrication variation, which is
demonstrated in FIG 3f, where we measure the Q-factor of different fabrication
runs of the same design and consistently achieve a high Q factor.

Selective Multiple Wavelength Band Operation.

The last design problem that we address is that of a microresonator with
resonant frequencies in different frequency bands. One possible application of
such a device would be to reduce pump leakage into undesired modes during a
non-linear wave mixing process and decrease conversion efficiency from the
mode competition process mckenna2021ultralowpower ; stone2021conversion . This
is schematically depicted in FIG 4a, wherein we contrast between parametric
down-conversion in a regular resonator, where coupling to many undesired
resonant modes is possible, and in a band-selective resonator, where such
couplings are absent. To achieve this with FP resonators, we can utilize
photonic inverse design to optimize for high reflectance at chosen wavelength
bands, and suppress the reflection at undesired wavelength bands. As a
demonstration, we optimize a silicon reflector that has high reflectance at
($1300$ $nm$ and $1550$ $nm$) wavelength bands, which correspond to O-band and
C-band, and reflectance suppression in the middle band ($1410$ $nm$). The
resulting inverse-designed device, with a footprint of $8\mu$m $\times$ 0.5
$\mu$m is shown in FIG 4b along with its simulated electric field profiles at
the selected wavelengths. As shown in the electric field profiles, the single
mirror reflects two well separated wavelengths ($1300$ $nm$ and $1550$ $nm$)
effectively; its simulated reflectance spectrum is shown in FIG 4c. The
schematic of the FP formed by the reflector is shown in FIG S3. To optically
characterize the resonator, we side coupling the resonator to a waveguide-
cavity directional coupler. We note that due to the difficulty in coupling
equally to a resonator across distant wavelength bands
Srinivasan:2019:NaturePhotonics  the resonances in O-band is weakly coupled,
while the resonances in the C-band is close to critically coupled. This can be
mediated by engineering coupler designs further
Srinivasan:2019:NaturePhotonics . FIG 4d depicts the measured transmission
spectrum of the selective wavelength band microresonator. This resonator can
achieve moderate loaded Q-factor in both of the selected bands (Q=1.5 $\sim$
7K), while the suppressed wavelength bands show loaded quality factor as low
as 600, which provides more than an order of magnitude difference between the
selected and suppressed wavelength bands.

Conclusion

In summary, we demonstrate photonic inverse design approaches for on-chip
microresonator designs that enable flexible engineering of important figure of
merits such as dispersion, high Q-factor operation and selective wavelength
band operation. As the proposed techniques can be combined together to
engineer for multiple desirable metrics at the same time, we foresee the
combined optimized techniques achieving nearly arbitrary resonator designs.
Furthermore, given this is the demonstration of a general approach not limited
to the SOI system used here as a model platform, it can easily extend to other
material platforms including those with transparency in the visible spectrum
and those with enhanced nonlinear properties like second order nonlinearity.
For instance, selective band operation can be extended to platforms such as
lithium niobateLoncar:2018:Nature  and silicon
carbideLukin:2020:NaturePhotonics , where we can target f-2f nonlinear
operations in visible and telecom wavelength bands. Even in silicon, achieving
a near-zero dispersion resonator opens up the possibility of resonator
enhanced EOM combs Afltouni:2020:CLEO ; Loncar:2019:Nature  and mid-IR
nonlinear wave mixingBell:Natpho  in commercial foundries, which is critical
for on-chip communications and spectroscopy.

## References

  * (1) Vahala, K. J.  Optical microcavities.  _Nature_ 424, 839–846 (2003).
  * (2) Xu, Q., Schmidt, B., Pradhan, S. & Lipson, M.  Micrometre-scale silicon electro-optic modulator.  _Nature_ 435, 325–327 (2005).
  * (3) Englund, D. _et al._ Controlling the spontaneous emission rate of single quantum dots in a two-dimensional photonic crystal.  _Physical Review Letters_ 95, 013904 (2005).
  * (4) Okawachi, Y. _et al._ Bandwidth shaping of microresonator-based frequency combs via dispersion engineering.  _Optics Letters_ 39, 3535–3538 (2014).
  * (5) Yang, K. Y. _et al._ Broadband dispersion-engineered microresonator on a chip.  _Nature Photonics_ 10, 316–320 (2016).
  * (6) Kim, S. _et al._ Dispersion engineering and frequency comb generation in thin silicon nitride concentric microresonators.  _Nature Communications_ 8, 1–8 (2017).
  * (7) Yu, S.-P., Jung, H., Briles, T. C., Srinivasan, K. & Papp, S. B.  Photonic-crystal-reflector nanoresonators for Kerr-frequency combs.  _ACS Photonics_ 6, 2083–2089 (2019).
  * (8) Zhang, Z., Norman, J. C., Liu, S., Malik, A. & Bowers, J. E.  Integrated dispersion compensated mode-locked quantum dot laser.  _Photonics Research_ 8, 1428–1434 (2020).
  * (9) Jin, W. _et al._ Hertz-linewidth semiconductor lasers using cmos-ready ultra-high-Q microresonators.  _Nature Photonics_ 1–8 (2021).
  * (10) Puckett, M. W. _et al._ 422 million intrinsic quality factor planar integrated all-waveguide resonator with sub-MHz linewidth.  _Nature Communications_ 12, 1–8 (2021).
  * (11) Yang, K. Y. _et al._ Bridging ultrahigh-Q devices and photonic circuits.  _Nature Photonics_ 12, 297–302 (2018).
  * (12) Pfeiffer, M. H. _et al._ Ultra-smooth silicon nitride waveguides based on the damascene reflow process: fabrication and loss origins.  _Optica_ 5, 884–892 (2018).
  * (13) Liu, X. _et al._ Ultra-high-Q UV microring resonators based on a single-crystalline AlN platform.  _Optica_ 5, 1279–1282 (2018).
  * (14) Desiatov, B., Shams-Ansari, A., Zhang, M., Wang, C. & Lončar, M.  Ultra-low-loss integrated visible photonics using thin-film lithium niobate.  _Optica_ 6, 380–384 (2019).
  * (15) Griffith, A. G. _et al._ Silicon-chip mid-infrared frequency comb generation.  _Nature Communications_ 6, 1–5 (2015).
  * (16) Piggott, A. Y. _et al._ Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer.  _Nature Photonics_ 9, 374–377 (2015).
  * (17) Callewaert, F., Butun, S., Li, Z. & Aydin, K.  Inverse design of an ultra-compact broadband optical diode based on asymmetric spatial mode conversion.  _Scientific Reports_ 6, 32577 (2016).
  * (18) Su, L. _et al._ Fully-automated optimization of grating couplers.  _Opt. Express_ 26, 4023–4034 (2018).
  * (19) Kippenberg, T. J., Gaeta, A. L., Lipson, M. & Gorodetsky, M. L.  Dissipative kerr solitons in optical microresonators.  _Science_ 361 (2018).
  * (20) Zhang, M. _et al._ Broadband electro-optic frequency comb generation in a lithium niobate microring resonator.  _Nature_ 568, 373–377 (2019).
  * (21) Idjadi, M. H., Arab, S. & Aflatouni, F.  Optical frequency comb generation in silicon by recursive electro-optic modulation.  In _CLEO: Science and Innovations_ , SF3O–5 (Optical Society of America, 2020).
  * (22) Riemensberger, J. _et al._ Dispersion engineering of thick high-q silicon nitride ring-resonators via atomic layer deposition.  _Optics express_ 20, 27661–27669 (2012).
  * (23) Runge, A. F. J., Hudson, D. D., Tam, K. K. K., de Sterke, C. M. & Blanco-Redondo, A.  The pure-quartic soliton laser.  _Nature Photonics_ 14, 492–497 (2020).
  * (24) Zhang, M. _et al._ Electronically programmable photonic molecule.  _Nature Photonics_ 13, 36–40 (2019).
  * (25) Karpov, M., Pfeiffer, M. H., Liu, J., Lukashchuk, A. & Kippenberg, T. J.  Photonic chip-based soliton frequency combs covering the biological imaging window.  _Nature communications_ 9, 1–8 (2018).
  * (26) Lee, S. H. _et al._ Towards visible soliton microcomb generation.  _Nature communications_ 8, 1–8 (2017).
  * (27) Yu, S.-P. _et al._ Spontaneous pulse formation in edgeless photonic crystal resonators.  _Nature Photonics_ 15, 461–467 (2021).
  * (28) Yu, X., Rao, A., Moille, G., Westly, A. D. & Srinivasan, K.  Universal frequency engineering tool for microcavity nonlinear optics: multiple selective mode splitting of whispering-gallery resonances.  _Photonics Research_ 8, 1676–1686 (2020).
  * (29) Lu, X. _et al._ Efficient telecom-to-visible spectral translation through ultralow power nonlinear nanophotonics.  _Nature Photonics_ 13, 593–601 (2019).
  * (30) Englund, D. _et al._ Controlling the spontaneous emission rate of single quantum dots in a two-dimensional photonic crystal.  _Phys. Rev. Lett._ 95, 013904 (2005).
  * (31) Biberman, A., Shaw, M. J., Timurdogan, E., Wright, J. B. & Watts, M. R.  Ultralow-loss silicon ring resonators.  _Optics Letters_ 37, 4236–4238 (2012).
  * (32) Timurdogan, E. _et al._ An ultralow power athermal silicon modulator.  _Nature Communications_ 5, 1–11 (2014).
  * (33) Atabaki, A. H. _et al._ Integrating photonics with silicon nanoelectronics for the next generation of systems on a chip.  _Nature_ 556, 349–354 (2018).
  * (34) Kordts, A., Pfeiffer, M. H. P., Guo, H., Brasch, V. & Kippenberg, T. J.  Higher order mode suppression in high-q anomalous dispersion sin microresonators for temporal dissipative kerr soliton formation.  _Opt. Lett._ 41, 452–455 (2016).
  * (35) McKenna, T. P. _et al._ Ultra-low-power second-order nonlinear optics on a chip (2021).  eprint 2102.05617.
  * (36) Yang, K. Y. _et al._ Inverse-designed non-reciprocal pulse router for chip-based LiDAR.  _Nature Photonics_ 14, 369–374 (2020).
  * (37) Szipöcs, R., Ferencz, K., Spielmann, C. & Krausz, F.  Chirped multilayer coatings for broadband dispersion control in femtosecond lasers.  _Optics letters_ 8, 1676–1686 (2020).
  * (38) Stone, J. R., Moille, G., Lu, X. & Srinivasan, K.  Conversion efficiency in kerr microresonator optical parametric oscillators: From three modes to many modes (2021).  eprint 2109.04910.
  * (39) Wang, C. _et al._ Integrated lithium niobate electro-optic modulators operating at cmos-compatible voltages.  _Nature_ 562, 101–104 (2018).
  * (40) Lukin, D. M. _et al._ 4H-silicon-carbide-on-insulator for integrated quantum and nonlinear photonics.  _Nature Photonics_ 14, 330–334 (2020).

Methods

Photonic Inverse Design
Stanford Photonics Inverse DesignSoftware (SPINS) was utilized to optimize the
various reflectors, and intracavity element shown through out the work. SPINS
can provide design of various on-chip elements as previously demonstrated.

Data availability The data that support the plots within this paper and other
findings of this study are available from the corresponding author upon
reasonable request.

Code availability An open source of the photonic optimization software used in
this paper is available at https://github.com/stanfordnqp/spins-b.

Acknowledgment We acknowledge insightful discussion with F.Aflatouni, S.Papp,
and K.Srinivasan. This work is funded by the DARPA under the PIPES and LUMOS
program. We thank G.Keeler and the programme management teams for discussions
throughout the project. G. H. A. acknowledges support from STMicroelectronics
Stanford Graduate Fellowship (SGF) and Kwanjeong Educational Foundation. R. T.
acknowledges Max Planck Harvard research center for quantum optics (MPHQ)
postdoctoral fellowship. A. W. acknowledges the Herb and Jane Dwight Stanford
Graduate Fellowship. J. S. acknowledges support from Cisco Systems Stanford
Graduate Fellowship (SGF)

Competing interests. The authors declare they have no competing financial
interests.

## Supplementary Information to
Augmenting On-Chip Microresonator through Photonic Inverse Design

Geun Ho Ahn1, Ki Youl Yang1, Rahul Trivedi1,2, Alexander D. White1, Logan Su1,
Jinhie Skarda1, Jelena Vučković1,†

1E. L. Ginzton Laboratory, Stanford University, Stanford, CA 94305, USA.
2Max-Planck-Institute for Quantum Optics, Hans-Kopfermann-Str. 1, 85748
Garching, Germany
$\dagger$Correspondence and requests for materials should be addressed to
<EMAIL_ADDRESS>

## Section I: Modelling Partially Reflecting Element in microring resonators

_Coupled mode equations for intracavity element_ : Here we develop coupled
mode theory for a general intracavity element in a whispering galley
resonator. As shown in the Fig S1, the counter directional propagating modes
(clockwise and counter clockwise propagating modes) $b(\omega)$ and
$a(\omega)$ can be formulated as:

$\displaystyle a_{in}(\omega)=a_{out}(\omega)\times e^{i\phi(\omega)}$ (1a)
$\displaystyle b_{in}(\omega)=b_{out}(\omega)\times e^{i\phi(\omega)}$ (1b)

The counter direction propagating modes can be coupled through the intracavity
element with its optical characteristics $|t|^{2}+|r|^{2}=1$ with following
coupled mode relation:

$\displaystyle a_{out}(\omega)=t\times a_{in}(\omega)+r\times b_{in}(\omega)$
(2a) $\displaystyle b_{out}(\omega)=-r^{*}\times a_{in}(\omega)+t\times
b_{in}(\omega)$ (2b)

Solving for the eigenvalue which is the coupling strength strength leads to
the coupling $\theta=arctan(\frac{|r|}{|t|})$.

The resonance conditions with and without the PRE are

$\displaystyle\phi_{\pm}(\omega)=2n\pi\pm\theta$ (3)

$\displaystyle\phi(\omega_{n})=2n\pi$ (4)

We use taylor expansion of $\phi_{\pm}(\omega)$ to obtain:

$\displaystyle\phi_{\pm}(\omega)=\phi_{\pm}(\omega_{n})+(\omega-\omega_{n})\phi^{\prime}(\omega_{n})=2n\pi\pm\theta$
(5)

$\displaystyle(\omega_{\pm}-\omega_{n})=\pm\frac{\theta}{\phi^{\prime}(\omega_{n})}$
(6)

where

$\displaystyle\phi^{\prime}(\omega_{n})=\frac{2\pi}{FSR}$ (7)

Hence the final resonance splitting ($\Gamma$) is

$\displaystyle\Gamma=\omega_{+}-\omega_{-}=\frac{arctan\left(\frac{|r|}{|t|}\right)}{\pi}\times
FSR$ (8)

Figure S1: Schematic of coupled mode theory for an intracavity element in
whispering gallery mode resonator

Figure S2: Measured integrated dispersion of three inverse designed FP
resonators with cavity length of 100 $\mu m$ and 480 $\mu m$

Figure S3: Schematic of side coupled inverse designed FP resonator with
selective reflection bands. This FP is coupled through a directional coupler
due to the flexibility in tuning of coupling strength by tuning the gap size
instead of targeting a particular through coupling strength. The gap size is
chosen to be 100 $nm$.

## Section II: Modelling TPA in optical parametric oscillations simulations

_Coupled mode equations for general cavities_ : Here we develop coupled mode
theory for a general resonator in the presence of TPA. We assume that the
resonator is described by modes $\textbf{E}_{\mu}(\textbf{x})$ and resonant
frequencies $\omega_{\mu}$ which satisfy:

$\displaystyle\nabla\times\nabla\times\textbf{E}_{\mu}(\textbf{x})=\mu_{0}\varepsilon_{0}\omega_{\mu}^{2}n^{2}(\textbf{x})\textbf{E}_{\mu}(\textbf{x})$
(9)

Furthermore, in the limit of a lossless cavity, the modes satisfy the
following orthonormality relationship:

$\displaystyle\int_{\Omega}n^{2}(\textbf{x})\textbf{E}_{\mu}^{*}(\textbf{x})\cdot\textbf{E}_{\nu}(\textbf{x})d^{3}\textbf{x}=\delta_{\mu,\nu}$
(10)

where $\Omega$ denotes the cavity region. The starting point for developing
the coupled mode theory is the Maxwell’s equations for the electric field:

$\displaystyle\nabla\times\nabla\times\textbf{E}(\textbf{x},t)+\mu_{0}\varepsilon_{0}n^{2}(\textbf{x})\frac{\partial^{2}}{\partial
t^{2}}\textbf{E}(\textbf{x},t)=-\mu_{0}\frac{\partial}{\partial
t}\textbf{J}_{\text{TPA}}(\textbf{x},t)-\mu_{0}\frac{\partial^{2}}{\partial
t^{2}}\textbf{P}_{\text{NL}}(\textbf{x},t)$ (11)

As a model for TPA, we setup a current that varies nonlinearly with
$\textbf{E}(\textbf{x},t)$:

$\displaystyle\textbf{J}_{\text{TPA}}(\textbf{x},t)=\begin{cases}\sigma_{\text{TPA}}\big{(}\textbf{E}(\textbf{x},t)\cdot\textbf{E}(\textbf{x},t)\big{)}\textbf{E}(\textbf{x},t)&\text{if
}\textbf{x}\in\Omega_{\text{NL}}\\\ 0&\text{otherwise}\end{cases}$ (12)

where $\Omega_{\text{NL}}$ is the region inside the cavity volume where the
nonlinear material is. Similarly, $\textbf{P}_{\text{NL}}(\textbf{x},t)$ has
the third order nonlinear form — we assume an isotropic polarization:

$\displaystyle\textbf{P}_{\text{NL}}(\textbf{x},t)=\begin{cases}\varepsilon_{0}\chi\big{(}\textbf{E}(\textbf{x},t)\cdot\textbf{E}(\textbf{x},t)\big{)}\textbf{E}(\textbf{x},t)&\text{if
}\textbf{x}\in\Omega_{\text{NL}}\\\ 0&\text{otherwise}\end{cases}$ (13)

Next, we write the electric field as a superposition of cavity modes:

$\displaystyle\textbf{E}(\textbf{x},t)=\underbrace{\bigg{(}\frac{\hbar\omega_{0}}{2\varepsilon_{0}}\bigg{)}^{1/2}\sum_{\mu}A_{\mu}(t)e^{-i\omega_{\mu}t}\textbf{E}_{\mu}(\textbf{x})}_{\tilde{\textbf{E}}(\textbf{x},t)}+\text{c.c.}$
(14)

where $\tilde{\textbf{E}}(\textbf{x},t)$ has only positive frequencies. We
note that $A_{\mu}(t)$ are dimensionless in this expression and
$|A_{\mu}(t)|^{2}$ is the number of photons in the $\mu^{\text{th}}$ cavity
mode at time $t$. Furthermore, we also assume that all the frequencies in this
term are very close to a central frequency $\omega_{0}$. With this form of the
field, the contribution of the following nonlinear expression in the electric
field to frequencies close to $\omega_{0}$ can be easily calculated:

$\displaystyle\big{(}\textbf{E}(\textbf{x},t)\cdot\textbf{E}(\textbf{x},t)\big{)}\textbf{E}(\textbf{x},t)\approx
2|\tilde{\textbf{E}}(\textbf{x},t)|^{2}\tilde{\textbf{E}}(\textbf{x},t)$ (15)

Here we have only retained the intensity dependent refractive index
contributions, and the parametric couplings in the nonlinear term, and ignored
harmonic generation. Defining

$\displaystyle\tilde{\textbf{J}}_{\text{TPA}}(\textbf{x},t)=\begin{cases}2\sigma_{\text{TPA}}|\tilde{\textbf{E}}(\textbf{x},t)|^{2}\tilde{\textbf{E}}(\textbf{x},t)&\text{if
}\textbf{x}\in\Omega_{\text{NL}}\\\ 0&\text{otherwise}\end{cases}\ \text{and}\
\tilde{\textbf{P}}_{\text{NL}}=\begin{cases}2\varepsilon_{0}\chi|\tilde{\textbf{E}}(\textbf{x},t)|^{2}\tilde{\textbf{E}}(\textbf{x},t)&\text{if
}\textbf{x}\in\Omega_{\text{NL}}\\\ 0&\text{otherwise}\end{cases}$ (16)

we obtain:

$\displaystyle\nabla\times\nabla\times\tilde{\textbf{E}}(\textbf{x},t)+\mu_{0}\varepsilon_{0}n^{2}(\textbf{x})\frac{\partial^{2}}{\partial
t^{2}}\tilde{\textbf{E}}(\textbf{x},t)=i\omega_{0}\mu_{0}\tilde{\textbf{J}}_{\text{TPA}}(\textbf{x},t)+\mu_{0}\omega_{0}^{2}\tilde{\textbf{P}}_{\text{NL}}(\textbf{x},t)$
(17)

Before proceeding further, we point out relationships between the model
parameters $\sigma_{\text{TPA}}$ and $\chi$ to the experimentally measured
$\beta_{\text{TPA}}$ (change in absorption coefficient per unit intensity) and
$n_{2}$ (change in refractive index per unit intensity). These parameters can
only be defined rigorously for plane-waves — supposes we consider a problem
where $n(\textbf{x})$ is uniform in space with value $n_{0}$ and consider
solutions of the form
$\tilde{E}(\textbf{x},t)=\hat{z}\mathcal{E}_{0}(x)e^{-i\omega t}$ to the above
equation then Eq. 17 reduces to:

$\displaystyle\frac{d^{2}\mathcal{E}(x)}{dx^{2}}+\frac{\omega_{0}^{2}}{c^{2}}\bigg{[}n_{0}^{2}+\bigg{(}2\chi+\frac{2i\sigma_{\text{TPA}}}{\omega_{0}\varepsilon_{0}}\bigg{)}|\mathcal{E}(x)|^{2}\bigg{]}\mathcal{E}(x)=0$
(18)

Furthermore, assuming $\mathcal{E}(x)$ to be very close to a plane-wave
solution $e^{ik_{0}x}$ where $k_{0}=\omega n_{0}/c$, we may calculate the
intensity in the electromagnetic field $I$ via:

$\displaystyle I(x)=2\epsilon_{0}cn_{0}|\mathcal{E}(x)|^{2}$ (19)

Now, from Eq. 18 we may conclude that the local refractive index $n(x)$
depends on the intensity $I(x)$ via:

$\displaystyle n(x)\approx
n_{0}+\frac{\chi|\mathcal{E}(x)|^{2}}{n_{0}}=n_{0}+\frac{\chi}{2\varepsilon_{0}cn_{0}^{2}}I(x)$
(20)

which yields $\chi=2n_{0}^{2}c\varepsilon_{0}n_{2}$. Similarly, from the
imaginary part of the local refractive index, we can calculate the intensity
dependent absorption coefficient (dimensions of 1 / length):

$\displaystyle\alpha(x)=\frac{2\sigma_{\text{TPA}}}{c\varepsilon_{0}n_{0}}|\mathcal{E}(x)|^{2}=\frac{\sigma_{\text{TPA}}}{c^{2}\varepsilon_{0}^{2}n_{0}^{2}}I(x)$
(21)

which yields
$\sigma_{\text{TPA}}=c^{2}\varepsilon_{0}^{2}n_{0}^{2}\beta_{\text{TPA}}$.
Note that:

$\displaystyle\nabla\times\nabla\times\tilde{\textbf{E}}(\textbf{x},t)+\mu_{0}\varepsilon_{0}n^{2}(\textbf{x})\frac{\partial^{2}}{\partial
t^{2}}\textbf{E}(\textbf{x},t)\approx-2i\mu_{0}\varepsilon_{0}\omega_{0}n^{2}(\textbf{x})\bigg{(}\frac{\hbar\omega_{0}}{2\varepsilon_{0}}\bigg{)}^{1/2}\sum_{\mu}\frac{dA_{\mu}(t)}{dt}e^{-i\omega_{\mu}t}\textbf{E}_{\mu}(\textbf{x})$
(22)

and thus

$\displaystyle\frac{dA_{\mu}(t)}{dt}=\bigg{(}\frac{2\varepsilon_{0}}{\hbar\omega_{0}}\bigg{)}^{1/2}\frac{\omega_{0}(i\varepsilon_{0}\omega_{0}\chi-\sigma_{\text{TPA}})e^{i\omega_{\mu}t}}{\omega_{\mu}\varepsilon_{0}}\int\big{|}\tilde{\textbf{E}}(\textbf{x},t)\big{|}^{2}\textbf{E}_{\mu}^{*}(\textbf{x})\cdot\tilde{\textbf{E}}(\textbf{x},t)d^{3}\textbf{x}$
(23)

We now consider the evaluation of the term in the integral: By a
straightforward substitution of the expression for
$\tilde{\textbf{E}}(\textbf{x},t)$, we obtain:

$\displaystyle\int\big{|}\tilde{\textbf{E}}(\textbf{x},t)\big{|}^{2}\textbf{E}_{\mu}^{*}(\textbf{x})\cdot\tilde{\textbf{E}}(\textbf{x})d^{3}\textbf{x}=\bigg{(}\frac{\hbar\omega_{0}}{2\varepsilon_{0}}\bigg{)}^{3/2}\frac{1}{n_{0}^{4}V_{0}}\sum_{\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\alpha}-\omega_{\nu}-\omega_{\beta})t}$
(24)

where

$\displaystyle\Lambda_{\mu.\nu,\alpha,\beta}=\frac{\int_{\Omega_{\text{NL}}}\big{[}\textbf{E}_{\alpha}^{*}(\textbf{x})\cdot\textbf{E}_{\beta}(\textbf{x})\big{]}\big{[}\textbf{E}_{\mu}^{*}(\textbf{x})\cdot\textbf{E}_{\nu}(\textbf{x})\big{]}d^{3}\textbf{x}}{\int_{\Omega_{\text{NL}}}|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}$
(25)

and

$\displaystyle
V_{0}=\frac{1}{n_{0}^{4}}\bigg{[}{\int_{\Omega_{\text{NL}}}|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}\bigg{]}^{-1}=\frac{1}{n_{0}^{4}}\bigg{[}\frac{\big{(}\int_{\Omega}n^{2}(\textbf{x})\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})d^{3}\textbf{x}\big{)}^{2}}{\int_{\Omega_{\text{NL}}}|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}\bigg{]}$
(26)

Here $n_{0}$ is the refractive inside the cavity. We note that for tightly
confined modes, $n^{2}(\textbf{x})\approx n_{0}^{2}$ inside the integrals in
the above expressions, and we obtain:

$\displaystyle\Lambda_{\mu,\nu,\alpha,\beta}\approx\frac{\int\big{[}\textbf{E}_{\alpha}^{*}(\textbf{x})\cdot\textbf{E}_{\beta}(\textbf{x})\big{]}\big{[}\textbf{E}_{\mu}^{*}(\textbf{x})\cdot\textbf{E}_{\nu}(\textbf{x})\big{]}d^{3}\textbf{x}}{\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}\
\text{and}\
V_{0}\approx{\frac{\big{(}\int\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}(\textbf{x})d^{3}\textbf{x}\big{)}^{2}}{\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}}$
(27)

Consequently, we obtain a nonlinear ODE for $A_{\mu}(t)$:

$\displaystyle\frac{dA_{\mu}(t)}{dt}=\frac{\hbar\omega_{0}(i\varepsilon_{0}\omega_{0}\chi-\sigma_{\text{TPA}})}{2n_{0}^{4}\varepsilon_{0}^{2}V_{0}}\sum_{\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\alpha}+\omega_{\mu}-\omega_{\nu}-\omega_{\beta})t}$
(28)

Next, we introduce the parametric coupling strength $g_{0}$ (with dimensions
of 1 / time) as:

$\displaystyle
g_{0}=\frac{\hbar\omega_{0}^{2}\chi}{2n_{0}^{4}\varepsilon_{0}V_{0}}=\frac{\hbar\omega_{0}^{2}n_{2}c}{n_{0}^{2}V_{0}}$
(29)

as well as the dimensionless parametric loss parameter $\xi$:

$\displaystyle\xi=\frac{\sigma_{\text{TPA}}}{\varepsilon_{0}\omega_{0}\chi}=\frac{\beta_{\text{TPA}}c}{2\omega_{0}n_{2}}$
(30)

with which:

$\displaystyle\frac{dA_{\mu}(t)}{dt}=ig_{0}(1+i\xi)\sum_{\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\alpha}+\omega_{\mu}-\omega_{\nu}-\omega_{\beta})t}$
(31)

Note that uptil now, we have ignored any losses in the cavity mode except for
the two-photon absorption loss, and have treated the linear cavity modes to be
lossless. An important consideration while studying losses in silicon based
devices is that the material absorption depends on the number of carriers in
the structure, which can become large at high field intensities due to carrier
generation based on two-photon absorption. Thus the decay rate of each cavity
mode is written as $\kappa_{\mu}+\gamma_{\mu}N(t)$, where $N(t)$ is the number
of carriers per unit volume inside the material:

$\displaystyle\boxed{\frac{dA_{\mu}(t)}{dt}=-\frac{\kappa_{\mu}+\gamma_{\mu}N(t)}{2}A_{\mu}(t)+ig_{0}(1+i\xi)\sum_{\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\alpha}+\omega_{\mu}-\omega_{\nu}-\omega_{\beta})t}}$
(32)

A simple estimate of $\gamma_{\mu}$ can be made from Maxwell’s equations by
(i) ignoring the spatial distribution of the generated carriers and (ii)
Assuming that the relaxation time-scale of the generated carriers is much
faster than the dynamics of the cavity mode. Under these assumptions, the
carriers simply introduce a conductivity
$\sigma=Ne^{2}\tau/m_{\text{eff}}(1+\omega_{0}^{2}\tau^{2})$. This
conductivity adds a current $\sigma\textbf{E}(\textbf{x},t)$ to Maxwell’s
equations above, which introduces the following decay term in the RHS of
equation for $dA_{\mu}(t)/dt$ (Eq. 23) from which we obtain the following
expression for $\gamma_{\mu}$:

$\displaystyle\gamma_{\mu}=\frac{\omega_{0}}{\omega_{\mu}}\bigg{[}\frac{e^{2}\tau}{m_{\text{eff}}(1+\omega_{0}^{2}\tau^{2})}\bigg{]}\int_{\text{Si}}|\textbf{E}_{\mu}(\textbf{x})|^{2}d^{3}\textbf{x}$
(33)

Next, we consider the dynamics of carrier generation due to two-photon
absorption. The number of carriers generated per unit time, $F(t)$ can be
estimated by calculating the total two-photon absorbed power and dividing it
by $2\hbar\omega_{0}$ (energy of two photons being absorbed). Using the
effective current in Eq. 12:

$\displaystyle
F(t)=\frac{\sigma_{\text{TPA}}}{2\hbar\omega_{0}}\bigg{\langle}\int(\textbf{E}(\textbf{x},t)\cdot\textbf{E}(\textbf{x},t))^{2}d^{3}\textbf{x}\bigg{\rangle}_{\text{opt.
cycle}}$ (34)

Here, we are assuming that the dynamics of the carrier generation is
comparable to the dynamics of the cavity mode, but still much slower than an
optical cycle. Thus, we perform an averaging over the optical cycle.
Performing this averaging, we obtain:

$\displaystyle
F(t)=\frac{4\sigma_{\text{TPA}}}{2\hbar\omega_{0}}\int(\tilde{\textbf{E}}^{*}(\textbf{x},t)\cdot\tilde{\textbf{E}}(\textbf{x},t))^{2}d^{3}\textbf{x}$
(35)

Using the expansion for $\tilde{\textbf{E}}(\textbf{x},t)$ in terms of the
cavity mode amplitude and assuming the cavity mode to be tightly confined, we
obtain:

$\displaystyle
F(t)=\frac{\sigma_{\text{TPA}}\hbar\omega_{0}}{2\varepsilon_{0}^{2}n_{0}^{4}V_{0}}\sum_{\mu,\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\mu}^{*}(t)A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\mu}+\omega_{\alpha}-\omega_{\mu}-\omega_{n}u)t}$
(36)

We can then introduce the following differential equation for the carrier
which accounts for an exponential decay of the carriers of the order of the
dwell time $\tau_{d}$ of the carrier inside the material:

$\displaystyle\boxed{\frac{dN(t)}{dt}=-\frac{N(t)}{\tau_{d}}+g_{0}\xi\sum_{\mu,\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\mu}^{*}(t)A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\mu}+\omega_{\alpha}-\omega_{\nu}-\omega_{\beta})t}}$
(37)

_Fabry Perot resonator_ : Several simplifications to the above coupled mode
equations happens for Fabry perot resonators by accounting for phase-matching
conditions between different resonant modes. These simplifications are
important since they make the dynamics easier to simulate. Consider. Fabry-
perot mode formed by mirrors at $x=0$ and $x=L$ with mirros that impart a
frequency dependent phase-shift $\phi(\omega)$. The resonator mode can be
expressed as:

$\displaystyle\textbf{E}_{\mu}(\textbf{x})=\textbf{E}^{(+)}_{\mu}(\rho)e^{i\beta_{\mu}x}+\textbf{E}^{(-)}_{\mu}(\rho)e^{-i\beta_{\mu}x}e^{-i\phi_{\mu}}$
(38)

where $\textbf{E}_{\mu}^{(\pm)}(\rho),\beta_{\mu},\phi_{\mu}$ are the
waveguide modes, propagation constant and mirror phase evaluated at
$\omega=\omega_{\mu}$ and $\rho\equiv(x,y)$ is the transverse coordinate.
Furthermore,
$\textbf{E}_{\mu}^{(-)}=\big{[}\textbf{E}_{\mu}^{(+)}\big{]}^{*}$. We can use
this form of the modes to evaluate $\Lambda_{\mu,\nu,\alpha,\beta}$ defined
above. For this calculation, we will only keep terms in the evaluation
$\Lambda_{\mu,\nu,\alpha,\beta}$ which are phase matched. Mathematically,
phase-matching is a consequence of approximating the following integral in the
limit of long length $L$:

$\displaystyle\int_{0}^{L}e^{i\Delta\beta x}dx=\frac{2e^{i\Delta\beta
L/2}\sin\Delta\beta L/2}{\Delta\beta}\approx\begin{cases}L&\text{if
}\Delta\beta=0\\\ 0&\text{otherwise}\end{cases}$ (39)

Furthermore, we also assume that $\beta_{\mu}$ and $\phi_{\mu}$ can be
approximated well by a linear function of $\mu$ as far as phase-matching is
concerned. This essentially means that the length of the resonator $L$ is much
larger than the inverse of the linear contribution to $\beta_{\mu}$, while
being smaller than the non-linear contributions. This immediately allows us to
write the following approximations to various integrals arising in the
calculation of $\Lambda_{\mu,\nu,\alpha,\beta}$:

$\displaystyle\int_{0}^{L}e^{\pm
i(\beta_{\mu}+\beta_{\nu}+\beta_{\alpha}+\beta_{\beta})x}dx\approx 0$ (40a)
$\displaystyle\int_{0}^{L}e^{\pm
i(\beta_{\mu}+\beta_{\nu}+\beta_{\alpha}-\beta_{\beta})x}dx\approx 0$ (40b)
$\displaystyle\int_{0}^{L}e^{i(\beta_{\mu}+\beta_{\nu}-\beta_{\alpha}-\beta_{\beta})x}dx\approx
L\delta_{\mu+\nu,\alpha+\beta}$ (40c)

A simple calculation then yields:

$\displaystyle\Lambda_{\mu,\nu,\alpha,\beta}=\text{M}_{\mu,\nu,\alpha,\beta}\delta_{\mu+\alpha,\nu+\beta}+\text{N}_{\mu,\nu,\alpha,\beta}\delta_{\mu+\beta,\alpha+\nu}+\text{P}_{\mu,\nu,\alpha,\beta}\delta_{\mu+\nu,\alpha+\beta}$
(41)

where

$\displaystyle\text{M}_{\mu,\nu,\alpha,\beta}=L\frac{\int\big{[}(\textbf{E}_{\mu}^{(+)*}(\rho)\cdot\textbf{E}_{\nu}^{(+)}(\rho))(\textbf{E}_{\alpha}^{(+)*}(\rho)\cdot\textbf{E}^{(+)}_{\beta}(\rho))+\text{c.c.}\big{]}d^{2}\rho}{\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}$
$\displaystyle\text{N}_{\mu,\nu,\alpha,\beta}=L\frac{\int\big{[}e^{i(\phi_{\alpha}-\phi_{\beta})}(\textbf{E}_{\mu}^{(+)*}(\rho)\cdot\textbf{E}_{\nu}^{(+)}(\rho))(\textbf{E}_{\alpha}^{(-)*}(\rho)\cdot\textbf{E}^{(-)}_{\beta}(\rho))+\text{c.c.}\big{]}d^{2}\rho}{\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}$
$\displaystyle\text{P}_{\mu,\nu,\alpha,\beta}=L\frac{\int\big{[}e^{i(\phi_{\alpha}-\phi_{\nu})}(\textbf{E}_{\mu}^{(+)*}(\rho)\cdot\textbf{E}_{\nu}^{(-)}(\rho))(\textbf{E}_{\alpha}^{(-)*}(\rho)\cdot\textbf{E}_{\beta}^{(+)}(\rho))+\text{c.c.}\big{]}d^{2}\rho}{\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}$
(42)

Another standard approximation that is needed to obtain an LLE like form is to
ignore dependence of
$\text{M}_{\mu,\nu,\alpha,\beta},\text{N}_{\mu,\nu,\alpha,\beta}$ and
$\text{P}_{\mu,\nu,\alpha,\beta}$ on the mode numbers and approximate them by
their value at the central mode. Doing so, we obtain that:

$\displaystyle\text{M}_{\mu,\nu,\alpha,\beta}\approx\text{N}_{\mu,\nu,\alpha,\beta}\approx\text{P}_{\mu,\nu,\alpha,\beta}\approx\Lambda_{0}=2L\frac{\int(\textbf{E}^{(+)*}_{0}(\rho)\cdot\textbf{E}_{0}^{(+)}(\rho))^{2}d^{2}\rho}{\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}}$
(43)

Using the form of the Fabry-perot modes, it can be seen that:

$\displaystyle\int|\textbf{E}_{0}^{*}(\textbf{x})\cdot\textbf{E}_{0}(\textbf{x})|^{2}d^{3}\textbf{x}=2L\bigg{(}1-\frac{\sin\phi_{0}}{\beta_{0}L}\bigg{)}\int\textbf{E}_{0}^{(+)*}(\rho)\cdot\textbf{E}^{(+)}_{0}(\rho)d^{2}\rho$
(44)

From this it follows that for long resonators ($\beta_{0}L\gg 1$),
$\Lambda_{0}\approx 1$. Consider now the nonlinear term in Eq. 32 — we now
make an assumption on the resonance frequencies $\omega_{\mu}$: we assume that
it is approximately linear with $\mu$:

$\displaystyle\omega_{\mu}=\omega_{0}+D_{1}\mu+\delta_{\mu}$ (45)

with $\delta_{\mu}\ll D_{1}$. We assume that $D_{1}$ is much larger than the
rate at which the system dynamics is happening, while the nonlinear
contribution is not. Then, we can perform a rotating wave approximation on the
time-variation in this equation which introduces a delta function in the
summation:

$\displaystyle\sum_{\mu,\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\alpha}+\omega_{\mu}-\omega_{\nu}-\omega_{\beta})t}\approx
e^{i\delta_{\mu}t}\sum_{\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}\delta_{\mu+\alpha,\nu+\beta}a_{\alpha}^{*}(t)a_{\nu}(t)a_{\beta}(t)$
(46)

where we have introduced $a_{\mu}(t)=A_{\mu}(t)e^{-i\delta_{\mu}t}$.
Substituting for $\Lambda_{\mu,\nu,\alpha,\beta}$, we obtain:

$\displaystyle\sum_{\nu,\alpha,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\alpha}^{*}(t)A_{\nu}(t)A_{\beta}(t)e^{i(\omega_{\alpha}+\omega_{\mu}-\omega_{\nu}-\omega_{\beta})t}\approx
e^{i\delta_{\mu}t}\bigg{[}\sum_{\nu,\alpha,\beta}a_{\nu+\beta-\mu}^{*}a_{\nu}a_{\beta}+2a_{\mu}\sum_{\alpha}|a_{\alpha}|^{2}\bigg{]}$
(47)

We can thus rewrite Eq. 32 as:

$\displaystyle\boxed{\frac{da_{\mu}(t)}{dt}=-\bigg{[}i\delta_{\mu}+\frac{\kappa_{\mu}+\gamma_{\mu}N(t)}{2}\bigg{]}a_{\mu}(t)+ig_{0}(1+i\xi)\bigg{[}\sum_{\nu,\beta}a_{\nu+\beta-\mu}^{*}a_{\nu}a_{\beta}+2a_{\mu}\sum_{\alpha}|a_{\alpha}|^{2}\bigg{]}+\sqrt{\kappa_{0}\eta}s_{\text{in}}\delta_{\mu,0}}$
(48)

where
$g_{0}=-\Lambda_{0}(i\sigma_{\text{TPA}}+\varepsilon_{0}\omega_{0}\chi)/2\varepsilon_{0}V_{0}$.
Similarly simplifying the nonlinear term in Eq. 37:

$\displaystyle\sum_{\mu,\alpha,\nu,\beta}\Lambda_{\mu,\nu,\alpha,\beta}A_{\mu}^{*}A_{\alpha}^{*}A_{\nu}A_{\beta}e^{i(\omega_{\mu}+\omega_{\alpha}-\omega_{\nu}-\omega_{\beta})t}\approx\Lambda_{0}\bigg{[}\sum_{\mu,\nu,\beta}a_{\mu}^{*}a_{\nu+\beta-\mu}^{*}a_{\nu}a_{\beta}+2\bigg{(}\sum_{\mu}|a_{\mu}|^{2}\bigg{)}^{2}\bigg{]}$
(49)

yielding

$\displaystyle\boxed{\frac{dN(t)}{dt}=-\frac{N(t)}{\tau_{d}}+g_{0}\xi_{0}\bigg{[}\sum_{\mu,\nu,\beta}a_{\mu}^{*}a_{\nu+\beta-\mu}^{*}a_{\nu}a_{\beta}+2\bigg{(}\sum_{\mu}|a_{\mu}|^{2}\bigg{)}^{2}\bigg{]}}$
(50)

Next, we put these equations in a dimensionless form for ease of analysis. We
introduce the pump threshold power (in photons per second) $p_{\text{th}}$
ignoring TPA:

$\displaystyle p_{\text{th}}=\frac{\kappa^{2}}{8g_{0}\eta}$ (51)

where $\kappa$ is the average decay rate of the resonators. Introducing
normalized mode amplitudes $\alpha_{\mu}=a_{\mu}\sqrt{2g_{0}/\kappa}$ and
normalized time $\tau=\kappa t/2$, we obtain:

$\displaystyle\frac{d\alpha_{\mu}(\tau)}{d\tau}=-(\lambda_{\mu}+\theta_{\mu}N(\tau))\alpha_{\mu}(\tau)+(i-\xi)\bigg{[}\sum_{\nu,\beta}\alpha_{\nu+\beta-\mu}^{*}(\tau)\alpha_{\nu}(\tau)\alpha_{\beta}(\tau)+2\alpha_{\mu}(\tau)\bigg{(}\sum_{\nu}|\alpha_{\nu}(\tau)|^{2}\bigg{)}\bigg{]}+\delta_{\mu,0}f_{\text{in}}(\tau)$
(52)

where $\lambda_{\mu}=(2i\delta_{\mu}+\kappa_{\mu})/\kappa$ and
$\theta_{\mu}=\gamma_{\mu}/\kappa$. Similarly:

$\displaystyle\frac{dN(\tau)}{d\tau}=-\lambda_{d}N(\tau)+\xi_{N}\bigg{[}\sum_{\mu,\nu,\beta}\alpha_{\mu}^{*}\alpha_{\nu+\beta-\mu}^{*}\alpha_{\nu}\alpha_{\beta}+2\bigg{(}\sum_{\mu}|\alpha_{\mu}|^{2}\bigg{)}^{2}\bigg{]}$
(53)

where $\lambda_{d}=2/\kappa\tau_{d}$ and $\xi_{N}=\kappa\xi/2g_{0}$.


# The relative importance of structure and dynamics on node influence in
reversible spreading processes

Junyi Qu  School of Physics and Electronic Science, East China Normal
University, Shanghai 200241, China  Ming Tang<EMAIL_ADDRESS>School of
Physics and Electronic Science, East China Normal University, Shanghai 200241,
China  Shanghai Key Laboratory of Multidimensional Information Processing,
East China Normal University, Shanghai 200241, China  Ying Liu
<EMAIL_ADDRESS>School of Computer Science, Southwest Petroleum
University, Chengdu 610500, China  Shuguang Guan<EMAIL_ADDRESS>School of Physics and Electronic Science, East China Normal University,
Shanghai 200241, China

###### Abstract

The reversible spreading processes with repeated infection widely exist in
nature and human society, such as gonorrhea propagation and meme spreading.
Identifying influential spreaders is an important issue in the reversible
spreading dynamics on complex networks, which has been given much attention.
Except for structural centrality, the nodes’ dynamical states play a
significant role in their spreading influence in the reversible spreading
processes. By integrating the number of outgoing edges and infection risks of
node’s neighbors into structural centrality, a new measure for identifying
influential spreaders is articulated which considers the relative importance
of structure and dynamics on node influence. The number of outgoing edges and
infection risks of neighbors represent the positive effect of the local
structural characteristic and the negative effect of the dynamical states of
nodes in identifying influential spreaders, respectively. We find that an
appropriate combination of these two characteristics can greatly improve the
accuracy of the proposed measure in identifying the most influential
spreaders. Notably, compared with the positive effect of the local structural
characteristic, slightly weakening the negative effect of dynamical states of
nodes can make the proposed measure play the best performance. Quantitatively
understanding the relative importance of structure and dynamics on node
influence provides a significant insight into identifying influential nodes in
the reversible spreading processes.

###### pacs:

87.19.X-,89.75.Hc,87.23.Ge

##  I Introduction

Identifying influential spreaders is an important topic in the spreading
dynamics on complex networks Koschützki et al. (2005); Lü et al. (2016); Pei
et al. (2020). Identifying these nodes can help us to formulate precise
marketing strategies Leskovec et al. (2007), control the spread of public
opinions Bovet and Makse (2019); Lin et al. (2018), prevent the catastrophic
disruptions in power girds Motter and Lai (2002); Albert et al. (2004) and
control the outbreak of epidemics Pastor-Satorras and Vespignani (2002);
Scarpino and Petri (2019); Zhou and Liu (2008).

In recent years, there have been a lot of researches devoted to identifying
the most influential spreaders Morone and Makse (2015); Pei et al. (2018); Hu
et al. (2018); Lokhov and Saad (2017); Zheng et al. (2021). As the network
structure has a significant impact on the spreading processes, some structural
centralities are used to measure the spreading influence of a node or node
set, such as degree, closeness centrality, betweenness centrality and
$k$-shell index Lü et al. (2016); Pei et al. (2020).

However, the complex interplays between network structure and spreading
dynamics sometimes make the structural centralities be unable to effectively
identify the influential nodes accurately in complex networks. Poux-Médard et
al. (2020); Erkol et al. (2020); Aral and Dhillon (2018). Integrating the
dynamical states of nodes in the spreading process can promote the accuracy of
measures in identifying the most influential spreaders Klemm et al. (2012).
Most of previous researches focus on irreversible spreading dynamics which has
a final state. But, some spreading processes are reversible with repeated
infections, such as meme spreading on online social networks Gleeson et al.
(2014), gonorrhea spreading propagation on sexual intercourse networks Pastor-
Satorras and Vespignani (2001a) and the rise and fall of stock prices on
financial stock market networks Stavroglou et al. (2019). These reversible
processes can be described by the SIS spreading dynamics Barzel and Barabási
(2013); Pastor-Satorras and Vespignani (2001b, c), where a node transfers
between infection state and recovery state. Qu $et$ $al$. proposed a single-
node control model to evaluate the influence of nodes and a structure-dynamics
combined centrality to identify influential nodes in the reversible spreading
system Qu et al. (2020). They found that taking the neighbors’ centrality and
dynamical states into account can identify influential nodes more accurately
in real-world networks than the benchmark centralities.

Although the interplay between network structure and dynamics determines the
influence of nodes, the relative importance of the two factors in identifying
influential nodes is still unknown. By considering the positive effects of the
neighbors’ local structures and the negative effects of the neighbors’
dynamical states, a new structure-dynamics combined centrality is proposed to
identify influential spreaders and evaluate the relative importance of
structure and dynamics on the determining node influence. In this new measure,
the number of outgoing edges and the infection risks of neighbors are used to
represent the local structures and dynamical states of the neighbors
respectively. Simulation results on real-world networks show that considering
neighbors’ local structures and dynamical states can help improving the
accuracy of the index in identifying influential spreaders. Notably, compared
with the number of outgoing edges of neighbors, slightly weakening the impact
of neighbors’ infection risks can greatly promote the accuracy of the index to
identify influential spreaders. There is an optimal range, within which the
local structures and dynamical states of neighbors enhance the accuracy of the
influence ranking index the most.

The rest of the paper is organized as follows. Section II introduces a single-
node control model to quantify the influence of a node in SIS spreading
processes. In Section III, based on the neighbors’ local structures and
dynamical states , we propose a node centrality index in the SIS spreading
dynamics. Section IV shows the numerical simulation results and analyses on
real-world networks. Section V is the conclusion and discussion.

##  II Single-node control model in the SIS spreading processes

In our previous research, we have proposed a single-node control method to
evaluate the influence of nodes in the reversible spreading system Qu et al.
(2020). In SIS spreading processes, a node has only two possible states:
susceptible (S) state and infected (I) state. An infected node transmits a
disease to each of its direct susceptible neighbors with rate $\beta$ and
recovers to susceptible state with rate $\mu$. After a long-time transient
process, the proportion of infected nodes in the network will be stable, which
is recorded as $\varphi_{0}$. Then we artificially set one node to be in the
infected state constantly. This will cause a change (improvement) in the
proportion of infected nodes in the entire network. After certain time steps,
the proportion of the I-state nodes reaches a new steady state and the
proportion of infected nodes is recorded as $\varphi$. We define
$\Delta\varphi=\varphi-\varphi_{0}$ as the spreading influence, also called
the spreading efficiency of the controlled node. The larger the
$\Delta\varphi$ is, the higher influence of the controlled node has.

In simulations, we use the synchronous updating method Shu et al. (2016). The
recovery rate is $\mu=0.1$ and the effective transmission rate
$\lambda=\beta/\mu$ should not be too large because when it is too large, the
propagation can easily spread to the entire network regardless of the node
under control Liu et al. (2015). In that case the difference between the
controlled nodes is diminished, and the influence of each node is
indistinguishable. Meanwhile, the effective transmission rate $\lambda$ should
not be too small. If it is too small, the propagation will be limited to the
local neighbourhood of the controlled node and cannot affect the entire
network. We obtain the SIS epidemic threshold $\lambda_{c}$ using the
susceptibility Ferreira et al. (2012); Shu et al. (2015); Xu et al. (2019)

$\chi(\lambda)=N\frac{<\rho^{2}>-<\rho>^{2}}{<\rho>},$ (1)

and choose an appropriate transmission rate $\lambda>\lambda_{c}$. Under this
effective transmission rate, the influence of nodes can be distinguished.

##  III The proposed centrality NSRC

It has been pointed out that compared with degree and $k$-shell index,
considering the contribution of neighbors to the initial spreader can
significantly improve the accuracy of ranking measures in the SIR-like
spreading processes Liu et al. (2016). If a neighbor $j$ of node $i$ only has
neighbors that are also the direct neighbor of node $i$, when we control node
$i$, the contribution of neighbor $j$ is less because it cannot assist node
$i$ to spread the message to the entire network.

Assume that the network $G=(V,E)$ is unweighted and undirected composed of
$|V|=N$ nodes and $|E|=M$ edges. The adjacency matrix of the network is
denoted as $\textbf{A}=(a_{ij})$, where $a_{ij}=1$ means there is a direct
connection between node $i$ and node $j$, otherwise $a_{ij}=0$. We further
consider the two-step neighbor $l$ of node $i$ connected to the node $i$
through the nearest neighbour node $j$. If $a_{ij}=1$, $a_{jl}=1$ and
$a_{il}\neq 1$, then we define $b_{jl}=1$ to represent an outgoing edge of
node $j$. The edge between the nearest neighbour node $j$ and node $l$ is the
outgoing edge of node $j$ relevant to node $i$. For each of the neighbor node
$j$, there is $k_{j}^{out}=\sum_{l\in\Gamma_{j}}{b_{jl}}$. Take the network in
Fig. 1 for example. Node 1 has four neighbors $\\{2,3,4,5\\}$ and the outgoing
edges of its nearest neighbour node 2 are the edges that connect node 2 and
the node set $\\{6,7,8,9\\}$. The outgoing edges of nodes 3, 4, 5 connect to
the node sets $\\{10,11\\}$, $\\{12,13,14\\}$ and $\\{15,16\\}$ respectively.
Considering the characteristics of SIS spreading processes, we introduce the
number of outgoing edges of neighbors into the node influential index with a
weighting parameter $a$. The idea is the more outgoing edges the neighbor
nodes have, the more they help the considered node to spread information or
message to the entire network Liu et al. (2016).

In the SIS spreading dynamics, when the system reaches a steady state, a node
can still be infected repeatedly. In order to maximize the spreading to the
entire network, we consider the risk that a node is infected. In Fig. 1, the
infection risk of node 3 is 0.5. When node 1 is controlled and spreads out
through node 3, because node 3 has a risk of infection, the influence of node
1 will be reduced. As the infection risk has a negative impact on the
influence of origin node, the term with $\rho$ has a minus sign. Parameter $b$
is an adjustable weight indicating the relative importance of local dynamical
states (i.e., infection risks).

Figure 1: Illustration of calculating the NSRC. The decimal number labeled to
the node number is the probability of being in the infected state. The NSRC of
node 1 is calculated as
$C(1)=k_{1}+a(k_{2}^{out}+k_{3}^{out}+k_{4}^{out}+k_{5}^{out})-b(k_{2}^{out}\rho_{2}+k_{3}^{out}\rho_{3}+k_{4}^{out}\rho_{4}+k_{5}^{out}\rho_{5})=4+a(4+2+3+2)-b(1.6+1+1.5+0.6)=4+11a-4.7b$

As considering neighbors’ local structures and dynamical states can
effectively improve the accuracy of ranking measure for node influence on them
Qu et al. (2020), we here investigate how to quantitatively evaluate the
importance of these two factors. We propose an index to rank node influence in
SIS reversible spreading processes called neighboring local structures and
infection risks coupled centrality (NSRC), which is defined as

$C(i)=k_{i}+a\sum_{j\in\Gamma_{i}}{k_{j}^{out}}-b\sum_{j\in\Gamma_{i}}{k_{j}^{out}\rho_{j}},$
(2)

where $k_{i}$ is the benchmark centrality (i. e., degree), $j\in\Gamma_{i}$ is
the direct neighbor of node $i$. The second item on the right side of Eq. (6)
represents the local structures of the neighbor node $j$, which is the number
of outgoing edges through all neighbors of node $i$. The third item represents
the contribution of dynamical state of neighbor node $j$, which is the
infection risk of node $j$. In the spreading process, when the system reaches
its steady state, a node can still be infected repeatedly. In order to
maximize the spread of infection to the entire network, we consider
controlling a node $i$ in I state. Suppose the neighboring node $j$ is in the
I state with infection probability $\rho_{j}$ in steady state. In simulations,
if the neighboring node $j$ is in I state at this time step, node $i$ will
have no effect on node $j$ and the influence of node $i$ will be reduced. As
the $\rho_{j}$ has a negative impact on the influence of node $i$, the term
with $\rho_{j}$ has a minus sign. For simplicity, only the one-step neighbors
of node $i$ and their outgoing edges are considered. $a$ and $b$ are two
weight parameters, which represent contribution of the the local structures
and the dynamical states of the neighbors respectively. The larger $a$
indicates the higher importance of the local structures and the larger $b$
indicates that the neighbors’ infections are more important. The combination
of $a$ and $b$ can be coordinated by their weights which will be discussed in
detail later.

In order to measure the accuracy of NSRC in identifying influential spreaders,
the following two evaluation methods are used. The imprecision function is
defined as

$\varepsilon(p)=1-\frac{\phi(p)}{\phi_{eff}(p)},$ (3)

where $p$ is the fraction of network size $N$. $\phi(p)$ is the average
spreading efficiency of $pN$ nodes with the highest NSRC, and $\phi_{eff}(p)$
is the average spreading efficiency of the $pN$ nodes with the highest
spreading efficiency. When $p$ is fixed, a smaller $\varepsilon(p)$ value
indicates a more accurate centrality in identifying influential spreaders.

The imprecision function reflects the accuracy of the NSRC in identifying the
most influential nodes. In order to measure the ability of NSRC in ranking all
nodes in the network, the Kendall’s $\tau$ correlation coefficient is used,
which is given by Kendall (1938)

$\tau=\frac{\sum_{i<j}\operatorname{sgn}[(x_{i}-x_{j})(y_{i}-y_{j})]}{\frac{1}{2}N(N-1)},$
(4)

where $\operatorname{sgn}(x)$ is a sign function. $\operatorname{sgn}(x)$
returns 1 if $x>0$, -1 if $x<0$, and 0 if $x=0$, respectively. $N$ is the
number of nodes in the list. $x_{i}$ and $x_{j}$ are the ranking orders of
node $i$ and node $j$ in sorting sequence 1. $y_{i}$ and $y_{j}$ are the
ranking orders of node $i$ and node $j$ in sorting sequence 2. If node $i$ and
node $j$ have a concordant order in sorting sequence 1 and 2,
$(x_{i}-x_{j})(y_{i}-y_{j})>0$. If node $i$ and node $j$ have a discordant
order in sorting sequence 1 and 2, $(x_{i}-x_{j})(y_{i}-y_{j})<0$. If node $i$
and node $j$ have a same order in sorting sequence 1 or 2,
$(x_{i}-x_{j})(y_{i}-y_{j})=0$. We rank the nodes by the proposed NSRC as the
sorting sequence 1 and rank the nodes by the simulated spreading efficiency as
sorting sequence 2. We compute the correlation coefficient of sorting sequence
1 and 2. The larger correlation coefficient is, the more concordant the
centrality and the spreading efficiency are. For an explicit comparison, we
calculate the improved $\tau$ ratio of the NSRC over the reference centrality,
which is

$\eta=\left\\{\begin{aligned}
\frac{\tau_{c}-\tau_{0}}{\tau_{0}}&&{\tau_{0}>0};\\\
\frac{\tau_{c}-\tau_{0}}{-\tau_{0}}&&{\tau_{0}<0};\\\
0&&{\tau_{0}=0},\end{aligned}\right.$ (5)

where $\tau_{c}$ is the correlation coefficient between the NSRC and spreading
efficiency, $\tau_{0}$ is the correlation coefficient between the neighboring
dynamical information centrality (NDIC) and spreading efficiency.

##  IV Numerical simulations

We apply the NSRC to eight real-world networks as listed in Table 1. The real-
world networks are: Netsci (collaboration network of network scientists),
Hamster (friendships and family network of a website), Router (the router
network collected by the Rocketfuel Project), Blog (communication network
between users on the MSN website), PGP (an encrypted communication network),
CA_hep (collaboration network of arxiv in high-energy physics theory), Astro
(collaboration network of astrophysics scientists) and Email (email
communication network of University College London) Newman (2006); Spring et
al. (2002); Boguna et al. (2004); Newman (2001); Kitsak et al. (2010).

Network | $N$ | $E$ | $<k>$ | $k_{max}$ | $H_{k}$ | $r$ | $C$ | $\lambda_{c}$ | $\lambda$
---|---|---|---|---|---|---|---|---|---
Net379 | 379 | 914 | 4.8 | 34 | 1.663 | -0.082 | 0.741 | 0.20 | 0.28
Hamster | 2000 | 16097 | 16.1 | 273 | 2.719 | 0.023 | 0.540 | 0.024 | 0.04
Router | 5022 | 6258 | 2.5 | 16 | 5.503 | -0.138 | 0.012 | 0.10 | 0.20
Blog | 3982 | 6803 | 3.4 | 189 | 4.038 | -0.133 | 0.284 | 0.10 | 0.20
PGP | 10680 | 24340 | 4.6 | 206 | 4.153 | 0.240 | 0.226 | 0.06 | 0.12
CA_Hep | 8638 | 24806 | 5.7 | 65 | 2.261 | 0.239 | 0.482 | 0.08 | 014
Astro | 14845 | 119625 | 16.1 | 360 | 2.820 | 0.228 | 0.670 | 0.02 | 0.10
Email | 12625 | 20362 | 3.2 | 576 | 34.249 | -0.387 | 0.109 | 0.01 | 0.06
Table 1: Characteristics of the real-world networks studied in this paper.
Theses characteristics include the number of nodes ($N$), number of edges
($E$), average degree ($<k>$), maximum degree ($k_{max}$), degree
heterogeneity ($H_{k}$), degree assortativity ($r$), clustering
coefficient($C$), epidemic threshold ($\lambda_{c}$) and disease-transmission
rate $\lambda$ used in the SIS spreading processes.

###  IV.1 The relative importance of neighbors’ local structures and dynamical
states

Next we focus on the most influential spreaders which are very critical in the
SIS spreading processes. We investigate the relative importance of the
neighbors’ local structures and dynamical states based on the top 5% nodes and
top 20% nodes ranked by the NSRC. Figures 2 and 3 show the the imprecisions of
the NSRC to identify the top 5% and 20% nodes under different values of the
weight parameters $a$ and $b$ respectively. In the Fig. 2, where the top 5%
influential spreaders are considered, and the blue areas in which the
imprecision of NSRC is relatively small correspond to the parameter range
$a\geq b$ in the lower right part of the diagram. In the Netsci, Router, Blog,
Email and PGP networks respectively, there is an optimal area at the lower
part of the diagonal $a=b$. In the Hamster, CA_Hep and Astro networks
respectively, the imprecision is relatively small in the area of $a<b$. There
is a saturation effect in this area where the imprecision does not change much
with the change of $a$ and $b$. When we consider the top ranked 20% nodes by
the NSRC, the results shown in Fig. 3 are similar to that of the top ranked 5%
nodes shown in Fig. 2. In the Netsci and Email networks, there are optimal
areas near the diagonal while in the other six networks, the optimal area is
the lower right part of $a>b$. This implies that for the proposed index NSRC
based on neighbors’ local structures and dynamical states to be more accurate,
the weight of the infection risks $b$ should not exceed the weight of the
neighbors’ local structures $a$.

Figure 2: The imprecision of NSRC as a function of $a$ and $b$ in eight real-
world networks for $p=5\%$. The real-world networks are Nstsci (a), Hamster
(b), Router (c), Blog (d), CA_Hep (e), Email (f), PGP (g), Astro (h).
Figure 3: The imprecision of NSRC as a function of $a$ and $b$ in eight real-
world networks for $p=20\%$. The real-world networks are Nstsci (a), Hamster
(b), Router (c), Blog (d), CA_Hep (e), Email (f), PGP (g), Astro (h).

To further explore the impact of the local structures and dynamical states of
neighbor nodes on the accuracy of the proposed centrality, we fix the value of
the parameter $a$ and investigate the imprecisions of NSRC as a function of
the parameter $b$. Figure 4 shows that when the weight parameter $a$ of the
neighbors’ local structure is fixed, the imprecision function has an optimal
value or saturation effect with the change of $b$. When $a=0.2$, in the Blog,
Email and PGP networks, the imprecisions are lowest at around
$b=0.1,b=0.2,b=0.2$ respectively. In the other five networks, the imprecisions
increase with the increase of $b$, and there is a saturation effect when
$b\leq a=0.2$. The results for $a=0.4$ and $a=0.8$ support the same
conclusion, as shown in Fig. 4. These results indicate that the same
conclusion also applies when $a=0.4$ and $a=0.8$ in Fig. 4. These results
indicate that the parameters $a$ and $b$ have a synergistic effect in
determining the accuracy of the centrality based on the neighbors’ local
structures and dynamical states. When the accuracy is the highest, the optimal
parameter $b$ changes with the parameter $a$.

Figure 4: The imprecision of NSRC as a function of $b$ in eight real-world
networks for $p=5\%$. The real-world networks are Nstsci (a), Hamster (b),
Router (c), Blog (d), CA_Hep (e), Email (f), PGP (g), Astro (h).

Through the experimental results in eight real-world networks, it is found
that regardless of the value of the parameter $a$, the imprecision function
will have a minimum value if $b$ is not greater than $a$. To understand the
above phenomenon, we can rewrite the centrality index as

$C(i)=k_{i}+a\sum_{j\in\Gamma_{i}}{k_{j}^{out}(1-\frac{b}{a}\rho_{j})}.$ (6)

In the form of Eq. (6), we can discuss the relationship between the items more
conveniently. The optimal area of $b<a$ indicates that although both the
positive effect of the outgoing edges and the negative effect of the infection
risks of neighbors’ have a significant contribution to the accuracy of the
ranking measure, the relative importance of them determines the ranking
accuracy. Specifically speaking, compared to the positive effect of the number
of neighbors’ outgoing edges, slightly weakening the negative effect of
neighbors’ infection risks can get better identification results.

The above results can be explained as follows. For a considered node, if two
of its neighbors have the same number of outgoing edges, a higher infection
risk of its neighbor $\rho_{j}$ means a limited impact towards the spreading
processes. However, a larger $\rho_{j}$ often implies that node $j$ tends to
be closely connected to the considered node and its neighbor set (i.e., node i
and its neighbor node $j$ having more common neighbors) Boguná et al. (2013);
Castellano and Pastor-Satorras (2012). For example, as there are more common
neighbors between node 1 and its neighbor node 3, the node 3 has a higher
infection risk (i. e., infection probability) than neighbor node 5. When the
considered node 1 is controlled as an I-state node, the neighbor node with
more common neighbors can enhance the influence of this neighbor node’s
outgoing propagation to a certain extent due to the echo chamber effect of
local close structures Zhang et al. (2014); Chen et al. (2018); Wang et al.
(2015). Combining the above two effects, an appropriate parameter combination
of $b$ being slightly less than $a$ can effectively and quantitatively
distinguish the impact of neighbors with different infection risks. For
example, for node 3 and node 5 in Fig. 1, under the echo chamber effect, the
infection risk increases from 0.5 to 0.8, and 0.3 to 0.65, respectively. When
$b=0$, the infection risks of these two nodes is not considered. Since the
number of outgoing edges is the same, the two nodes have the same contribution
to the spreading capacity of the controlled node 1. When $b=1$, although the
infection risk of node 5 is significantly greater than that of node 3, the
negative effect of high infection risk towards the spreading range of the node
1 is overvalued due to the enhancement of the echo chamber effect of local
close structures. But in fact we need to consider the impact of the echo
chamber effect and appropriately decrease the weight of the neighbor’s
infection risk, thus the value of $b$ should be less than $a$.

###  IV.2 Compared with NDIC

In our previous study, we proposed a similar ranking index based on
neighboring dynamical information (NDIC) Qu et al. (2020). Compared with the
NDIC, the newly proposed NSRC classifies the neighbors’ edges and especially
addresses the impact of outgoing edges. For a better illustration, we compared
the two indicators by the Kendall’s $\tau$ correlation and calculated the
improved $\tau$ ratio $\eta$ of the new index over the existing NDIC. It has
been discussed in Ref. Qu et al. (2020) that when $a=0.2$, the accuracy of
NDIC is the highest or a saturation effect appears when $b=0.5$.

Figure 5: The improved ratio $\eta$ of NSRC over NDIC as a function of $p$.
The weight parameters of $a$ and $b$ are both $a=0.2,b=0.5$. The real-world
networks are Nstsci (a), Hamster (b), Router (c), Blog (d), CA_Hep (e), Email
(f), PGP (g), Astro (h).
Figure 6: The improved ratio $\eta$ of NSRC over NDIC as a function of $p$.
The weight parameters of $a$ and $b$ for NSRC and NDIC are $a=0.2,b=0.1$ and
$a=0.2,b=0.5$ respectively. The real-world networks are Nstsci (a), Hamster
(b), Router (c), Blog (d), CA_Hep (e), Email (f), PGP (g), Astro (h).

In this combination of $a$ and $b$, the increased correlation ratio $\eta$ of
NSRC over NDIC is greater than 0 in most cases as shown in Fig. 5. We also see
that the NSRC displays a greater increase in ranking accuracy than the NDIC in
the all networks, when an appropriate combination of parameters is fixed, as
shown in Fig. 6. If we choose an optimal parameter combination of $a=0.4$ and
$b=0.2$, the improved $\tau$ ratio of NSRC over NDIC is much more significant
as shown in Fig. 7. In the Router network, the improved $\tau$ ratio $\eta$ is
as high as 90%. In the Blog network, $\eta$ is close to 50%. In the other six
networks, $\eta$ is up to 10% to 30%.

Figure 7: The improved ratio $\eta$ of NSRC over NDIC as a function of $p$.
The weight parameters of $a$ and $b$ for NSRC and NDIC are $a=0.4,b=0.2$ and
$a=0.2,b=0.5$ respectively. The real-world networks are Nstsci (a), Hamster
(b), Router (c), Blog (d), CA_Hep (e), Email (f), PGP (g), Astro (h).

Next, we investigate the dependence of $\eta$ on the transmission rate
$\lambda$. The transmission rate should not be too large when distinguishing
the node influence by using the single-node control model. Because when
$\lambda$ is too large, a node with low centrality can also spread the disease
to the entire network. In this case the influence of nodes cannot be
distinguished. We thus use the 1 to 3 times of the epidemic threshold as the
transmission rate Liu et al. (2015). The improved ratio $\eta$ of NSRC over
NDIC as a function of transmission rate $\lambda$ in eight real-world networks
are shown in Fig. 8. The weight parameters are set to be $a=0.4,b=0.2$ and
$a=0.2,b=0.5$ for NSRC and NDIC respectively. In the Netsci, Router, and Blog
networks, the improved $\tau$ ratio is about 80%. In the rest five networks,
$\tau$ is around 10% to 30%. These results indicate that the ranking accuracy
of NSRC is better than that of NDIC at various transmission rates.

Figure 8: The improved ratio $\eta$ of NSRC over NDIC as a function of
$\lambda$ in eight real-world networks. The imprecisions for $p=5\%$(red
column) and $20\%$(blue column) are compared in each network. The real-world
networks are Nstsci (a), Hamster (b), Router (c), Blog (d), CA_Hep (e), Email
(f), PGP (g), Astro (h).

##  V Conclusions

In this paper, we quantitatively analyzed the relative importance of the
neighbors’ local structures and dynamical states in ranking node influence in
the SIS-like spreading processes. We proposed a centrality index NSRC that
combines the positive effect of the neighbors’ outgoing edges and the negative
effect of the neighbors’ infection risks with two weighting factors $a$ and
$b$. A large number of neighbors’ outgoing edges have a positive effect on the
spreading influence of a node, while the high infection risks of the neighbor
nodes reduce the influence of the considered node. A larger value of $a$ or
$b$ means a stronger positive local structure effect or negative local
dynamical effect. Through a large number of simulations on eight real-world
networks, we found that when the positive and negative effects are taken
appropriately, under appropriate values of $a/b$, the proposed centrality
index NSRC can improve the ranking accuracy significantly.

Specifically, to identify the most influential nodes accurately, the weight of
the neighbors’ dynamical states $b$ should not be greater than the weight of
the neighbors’ local structures $a$. Although the high infection risk of a
neighbor node can inhibit the spreading capacity of the controlled node, the
close local structure with more common neighbors between the considered node
and its neighbor node can induce the echo chamber effect, and thus slightly
promote the spreading processes. Therefore, an optimal combination of
parameters with $b<a$ can make the NSRC evaluate the spreading influence of a
node more accurately. This is helpful to distinguish the impact of the
neighboring nodes having the same local structures but different infection
probabilities. By analyzing the simulation results of real-world networks, it
is found that when the outgoing edges of neighbor nodes are taken into
account, the accuracy of NSRC is greatly improved over the centrality measure
NDIC. After further optimizing the parameters, the improved ranking
correlation of NSRC over NDIC can be as high as 10% to 100%. Compared with the
degree centrality of the neighbors’, the neighbors’ outgoing edges can be more
accurate in reflecting the importance of neighbors’ local structures in the
spreading processes. Further more, when we change the transmission rate
$\lambda$, the NSRC has a higher accuracy than the existing index in ranking
influential spreaders in the SIS-like spreading processes.

In the reversible processes, we found that considering the positive effects of
neighbors’ outgoing edges and the negative effects of neighbors’ infection
risks can effectively improve the accuracy of indexes in identifying
influential spreaders. This deepens our understanding of identifying
influential spreaders in the reversible processes. The benchmark centrality
used in NSRC in our paper is degree, which can be substituted by any other
benchmark centralities, such as $k$-shell index centrality Kitsak et al.
(2010) and the collective influence index Morone and Makse (2015). The
proposed ranking index NSRC can be applied in other reversible dynamics and
dynamical systems with steady state such as synchronization processes Barzel
and Barabási (2013) and cascading failures Lin et al. (2020).

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 11875132, 11975099, 82161148012, 11835003, 61802321), the Natural
Science Foundation of Shanghai (Grant No. 18ZR1411800), and the Science and
Technology Commission of Shanghai Municipality (Grant No. 14DZ2260800).

## References

  * Koschützki et al. (2005) D. Koschützki, K. A. Lehmann, L. Peeters, S. Richter, D. Tenfelde-Podehl, and O. Zlotowski, in _Network Analysis_ (Springer, 2005), pp. 16–61.
  * Lü et al. (2016) L. Lü, D. Chen, X. Ren, Q. Zhang, Y. Zhang, and T. Zhou, Physics Reports 650, 1 (2016).
  * Pei et al. (2020) S. Pei, J. Wang, F. Morone, and H. A. Makse, Journal of Complex Networks 8, cnz029 (2020).
  * Leskovec et al. (2007) J. Leskovec, L. A. Adamic, and B. A. Huberman, ACM Transactions on the Web (TWEB) 1, 5 (2007).
  * Bovet and Makse (2019) A. Bovet and H. A. Makse, Nature Communications 10, 1 (2019).
  * Lin et al. (2018) Y.-T. Lin, X.-P. Han, B.-K. Chen, J. Zhou, and B.-H. Wang, Frontiers of Physics 13, 1 (2018).
  * Motter and Lai (2002) A. E. Motter and Y. Lai, Physical Review E 66, 065102 (2002).
  * Albert et al. (2004) R. Albert, I. Albert, and G. L. Nakarado, Physical Review E 69, 025103 (2004).
  * Pastor-Satorras and Vespignani (2002) R. Pastor-Satorras and A. Vespignani, Physical Review E 65, 036104 (2002).
  * Scarpino and Petri (2019) S. V. Scarpino and G. Petri, Nature Communications 10, 1 (2019).
  * Zhou and Liu (2008) J. Zhou and Z.-H. Liu, Frontiers of Physics 3, 331 (2008).
  * Morone and Makse (2015) F. Morone and H. A. Makse, Nature 524, 65 (2015).
  * Pei et al. (2018) S. Pei, F. Morone, and H. A. Makse, in _Complex spreading phenomena in social systems_ (Springer, 2018), pp. 125–148.
  * Hu et al. (2018) Y. Hu, S. Ji, Y. Jin, L. Feng, H. E. Stanley, and S. Havlin, Proceedings of the National Academy of Sciences 115, 7468 (2018).
  * Lokhov and Saad (2017) A. Y. Lokhov and D. Saad, Proceedings of the National Academy of Sciences 114, E8138 (2017).
  * Zheng et al. (2021) K. Zheng, Y. Liu, Y. Wang, and W. Wang, arXiv preprint arXiv:2101.02335 (2021).
  * Poux-Médard et al. (2020) G. Poux-Médard, R. Pastor-Satorras, and C. Castellano, Physical Review Research 2, 023332 (2020).
  * Erkol et al. (2020) Ş. Erkol, D. Mazzilli, and F. Radicchi, Physical Review E 102, 042307 (2020).
  * Aral and Dhillon (2018) S. Aral and P. S. Dhillon, Nature Human Behaviour 2, 375 (2018).
  * Klemm et al. (2012) K. Klemm, M. Á. Serrano, V. M. Eguíluz, and M. San Miguel, Scientific Reports 2, 1 (2012).
  * Gleeson et al. (2014) J. P. Gleeson, J. A. Ward, K. P. O’sullivan, and W. T. Lee, Physical Review Letters 112, 048701 (2014).
  * Pastor-Satorras and Vespignani (2001a) R. Pastor-Satorras and A. Vespignani, Physical Review E 63, 066117 (2001a).
  * Stavroglou et al. (2019) S. K. Stavroglou, A. A. Pantelous, H. E. Stanley, and K. M. Zuev, Proceedings of the National Academy of Sciences 116, 10646 (2019).
  * Barzel and Barabási (2013) B. Barzel and A. Barabási, Nature Physics 9, 673 (2013).
  * Pastor-Satorras and Vespignani (2001b) R. Pastor-Satorras and A. Vespignani, Physical Review Letters 86, 3200 (2001b).
  * Pastor-Satorras and Vespignani (2001c) R. Pastor-Satorras and A. Vespignani, Physical Review E 63, 066117 (2001c).
  * Qu et al. (2020) J. Qu, M. Tang, Y. Liu, and S. Guan, Chaos, Solitons & Fractals 140, 110197 (2020).
  * Shu et al. (2016) P. Shu, W. Wang, M. Tang, P. Zhao, and Y. Zhang, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 063108 (2016).
  * Liu et al. (2015) Y. Liu, M. Tang, T. Zhou, and Y. Do, Scientific Reports 5, 9602 (2015).
  * Ferreira et al. (2012) S. C. Ferreira, C. Castellano, and R. Pastor-Satorras, Physical Review E 86, 041125 (2012).
  * Shu et al. (2015) P. Shu, W. Wang, M. Tang, and Y. Do, Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 063104 (2015).
  * Xu et al. (2019) Y. Xu, M. Tang, Y. Liu, Y. Zou, and Z. Liu, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 103141 (2019).
  * Liu et al. (2016) Y. Liu, M. Tang, T. Zhou, and Y. Do, Physica A: Statistical Mechanics and its Applications 452, 289 (2016).
  * Kendall (1938) M. G. Kendall, Biometrika 30, 81 (1938).
  * Newman (2006) M. E. Newman, Physical Review E 74, 036104 (2006).
  * Spring et al. (2002) N. Spring, R. Mahajan, and D. Wetherall, ACM SIGCOMM Computer Communication Review 32, 133 (2002).
  * Boguna et al. (2004) M. Boguna, R. Pastorsatorras, A. Diazguilera, and A. Arenas, Physical Review E 70, 056122 (2004).
  * Newman (2001) M. E. Newman, Proceedings of the National Academy of Sciences 98, 404 (2001).
  * Kitsak et al. (2010) M. Kitsak, L. K. Gallos, S. Havlin, F. Liljeros, L. Muchnik, H. E. Stanley, and H. A. Makse, Nature Physics 6, 888 (2010).
  * Boguná et al. (2013) M. Boguná, C. Castellano, and R. Pastor-Satorras, Physical Review Letters 111, 068701 (2013).
  * Castellano and Pastor-Satorras (2012) C. Castellano and R. Pastor-Satorras, Scientific Reports 2, 1 (2012).
  * Zhang et al. (2014) H. Zhang, J. Xie, M. Tang, and Y. Lai, Chaos: An Interdisciplinary Journal of Nonlinear Science 24, 043106 (2014).
  * Chen et al. (2018) X. Chen, R. Wang, M. Tang, S. Cai, H. E. Stanley, and L. A. Braunstein, New Journal of Physics 20, 013007 (2018).
  * Wang et al. (2015) W. Wang, M. Tang, H. Zhang, and Y. Lai, Physical Review E 92, 012820 (2015).
  * Lin et al. (2020) Z. Lin, M. Feng, M. Tang, Z. Liu, C. Xu, P. M. Hui, and Y. Lai, Nature Communications 11, 1 (2020).