#  Spatial Mapping of Electrostatics and Dynamics across 2D Heterostructures

Akshay A. Murthy1,2, Stephanie M. Ribet1, Teodor K. Stanev3, Pufan Liu1,
Kenji Watanabe4, Takashi Taniguchi5, Nathaniel P. Stern3, Roberto dos
Reis1,6,∗,
Vinayak P. Dravid1,2,6,∗

( 1Department of Materials Science and Engineering, Northwestern University,
Evanston, IL, USA
2International Institute of Nanotechnology, Northwestern University, Evanston,
IL, USA
3Department of Physics and Astronomy, Northwestern University, Evanston, IL,
USA
4Research Center for Functional Materials, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan
5International Center for Materials Nanoarchitectonics, National Institute for
Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan
6The NUANCE Center, Northwestern University, Evanston, IL, USA

∗Corresponding Authors:
Roberto dos Reis<EMAIL_ADDRESS>
Vinayak P. Dravid<EMAIL_ADDRESS>)

In situ electron microscopy is a key tool for understanding the mechanisms
driving novel phenomena in 2D structures. Unfortunately, due to various
practical challenges, technologically relevant 2D heterostructures prove
challenging to address with electron microscopy. Here, we use the differential
phase contrast imaging technique to build a methodology for probing local
electrostatic fields during electrical operation with nanoscale precision in
such materials. We find that by combining a traditional DPC setup with a high
pass filter, we can largely eliminate electric fluctuations emanating from
short-range atomic potentials. With this method, a priori electric field
expectations can be directly compared with experimentally derived values to
readily identify inhomogeneities and potentially problematic regions. We use
this platform to analyze the electric field and charge density distribution
across layers of hBN and MoS2.

KEYWORDS: In situ electron microscopy, heterostructure, transition metal
dichalcogenides, tunneling contacts, MoS2, differential phase contrast

## I INTRODUCTION

The diverse class of two-dimensional (2D) materials including graphene,
hexagonal boron nitride (hBN), and transition metal dichalcogenides (TMDs),
have produced a wide variety of captivating and technologically relevant
features for next generation optoelectronic architectures and quantum
information science.[1, 2, 3] These include unique memory devices,[4, 5]
Josephson junctions for solid state qubit architectures,[6, 7] photonic
structures,[8, 9, 10] and electronic devices.[11, 12, 13] Although charge
transport dynamics present at various 2D or contact interfaces contribute
significantly to device performance in such structures, the atomistic features
giving rise to macroscopic properties often remains unclear. This makes it
difficult to identify the optimal methods for improving performance and
reliability. In situ electrical biasing scanning/transmission electron
microscopy (S/TEM) is a particularly attractive method for probing property-
performance relationships at these atomic length scales. This approach
involves stimulating the sample of interest with an external electrical field
while simultaneously probing it with a high energy electron beam. It is thus
possible to directly ascertain the evolution of the atomic structure with
nanoscale spatial resolution, such as electromigration of atoms in MoS2.[14]

Beyond providing structural information, the multimodal nature of STEM makes
it possible to detect emergent phenomena. For instance, electric and/or
magnetic fields within the sample partially deflect the incident electron
probe leading to shifts in the intensity center of mass (ICOM) in the
diffraction patterns. The result of this beam-specimen interaction can be
measured at the sample plane through methods including differential phase
contrast (DPC).[15, 16, 17] Traditional DPC makes use of the idea that the
difference in measured diffraction intensities between detector regions
diametrically opposite one another is directly proportional to the deflection
of the electron probe along the common detector axis.[18, 19, 20, 21]
Similarly, first moment (FM)–STEM configuration calculates ICOM directly by
summing up the measured intensity at detector pixel. High-speed direct
detection cameras [22, 23, 24] have opened the possibility to capture the
scattered electron probe at each pixel in a real space image with roughly 10-3
second dwell times [25] and enabled both DPC methods.[26]

Here we investigate 2D layers that can be simultaneously biased electronically
while being probed with an electron beam. Using DPC to calculate variations in
the center of mass of the electron probe, we are able to map the resultant
electrostatic fields and charge densities that arise in 2D layers when an
external bias is applied. By combining a traditional DPC setup with a high
pass filter, we are able to selectively isolate the low spatial frequency
information associated with these long-range potentials. This method paves the
way for exploring a wide variety of 2D memory, quantum, and optoelectronic
heterostructure devices.

Figure 1: Schematic of STEM electrical field measurement setup (a) In this
approach, an electrical potential is applied across a 2D material system
leading to a deflection in the center of mass of the electron diffraction
pattern. (b) This momentum transfer experienced by the incident electrons in
the detector plane is used to calculate the electric field vector at each
probe position.

## II SAMPLE PREPARATION

Although fabricating TEM devices of 2D materials and heterostructures is
typically a challenging task, it is possible to construct such structures
using a process described previously.[14] In this setup, 75 nm thick gold
electrodes are patterned using electron beam lithography on a SiO2/Si (300 nm
thickness) substrate. This substrate is then spin-coated with a polycarbonate
solution (5% polycarbonate to 95% chloroform w/v) at 2000 RPM for 60 seconds.
After this coating is baked at 120∘C for 1 minute, the substrate is placed in
DI water until the polymer containing the gold electrodes naturally released
from the substrate. Once the polymer film had dried, it is carefully placed on
two strips of PDMS gel attached to a glass slide and suspended roughly 1 mm
above the glass slide. Using a micromanipulator stage, the gold electrodes are
aligned over a TEM window (TEMwindows.com) and brought into contact with the
TEM grid. This assembly is then gradually heated up to 150∘C for a few minutes
before the sample is baked at 170∘C for 15 minutes. The structure is then
placed in a chloroform bath for 6 hours to dissolve the polycarbonate film and
moved to an IPA bath afterwards for a few minutes to remove chloroform
residue.

Once a thin layer of mechanically exfoliated hBN (Figure S1) or a mechanically
transferred MoS2/hBN heterostructure is identified on a SiO2/Si substrate, a
similar process is then repeated to move the layers from the substrate to the
TEM grid. The resultant sample is left overnight in chloroform and moved to an
IPA bath afterwards. The sample is then annealed in a 95% Ar/5% H2 mixture at
200∘C for 2 hours to remove residue prior to STEM imaging. Finally, gold wires
are used to connect the electrode pads to the contact points on the
Nanofactory® TEM holder using indium as a bonding agent. Following fabrication
of the sample, Raman and photoluminescence (PL) spectra are taken to confirm
the presence of suspended layers of MoS2 and hBN (Figure S2). The current-
voltage characteristics that arise when a voltage is applied across this
sample is provided in (Figure S2) demonstrating the ability to apply an
external field across the sample.

## III EXPERIMENTAL

S/TEM images of 2D layers are acquired using a JEOL ARM 300F Grand ARM S/TEM
operated at 300 kV. The camera length is set to 8 cm and the condenser
aperture is chosen to result in a convergence semi-angle of 30 mrad with beam
current 20pA to minimize beam-induced damage. Application of a high convergent
angle (30 mrad) is intentionally chosen to optimize the measurement of a
constant field in the sample that causes a uniform shift of the unscattered
probe. The asymmetry in the measured pattern creates a difference signal that
is proportional to the field strength. 4D STEM data set containing over 40,000
convergent beam electron diffraction (CBED) patterns are acquired in a 680 x
80 mesh (pixel size of 5.96Å) across the sample using Gatan® K3-IS in counting
mode. A frame rate of 285 fps (0.0035s per diffraction) is chosen to
prioritize acquisition of high resolution CBED (1024 x 1024 pixels) while
maintaining sufficient sample stability under external bias.

## IV THEORY: Deducing In-Plane Electric Fields from Intensity Variations in
CBED Patterns

Intensity variations in the CBED pattern can be related to in-plane electric
fields through the relations discussed in this section. When an electric field
is present in the sample plane, the incident beam of electrons experience a
Coulombic force leading to a deflection in the CBED pattern parallel to the
detector plane. As discussed, DPC and FM-STEM offer methods to measure these
changes in the ICOM of the CBED pattern as a function of position. By
rearranging DeBroglie’s wave equation and normalizing the overall electron
probe intensity, the ICOM values calculated at each probe position can be
related to the expectation value of the momentum transfer,
$\left<p_{\bot}\right>$ (Eq. 1). In this equation, h represents Planck’s
constant and $\lambda$ represents electron wavelength.

$\left<p_{\bot}\right>=h\frac{\sin(I_{COM})}{\lambda}$ (1)

Further, [27] showed that based on Ehrenfest’s theorem, the measured in-plane
electric field is given by Eq. 2.[27] In this equation, $v$ represents the
velocity of electrons carrying a charge $e$ that interact with the sample for
a time equal to $v$/$z$, where $z$ represents the thickness of the sample.[27,
28]

$E_{\bot}=\left<p_{\bot}\right>\frac{v}{ez}$ (2)

Because this equation assumes that the intensity of the electron probe does
not vary with $z$, this relation requires specimen thicknesses to be
restricted to values less than 5 nm to avoid beam broadening effects.[27, 29]
In the case of thicker samples, sources of diffraction contrast can also
contribute to probe deflection and introduce imaging artifacts.[30]

Figure 2: DPC simulations for interrogating long-range electric fields (a-b)
Simulated DPC images taken from hBN structure subjected to an varying
electrical potential in the vertical direction detailed in Figure S4. The
detector configurations are provided in the insets. (c) An intensity line
profile taken from (a). (d-e) Simulated annular DPC images taken from same hBN
structure with detector configurations are provided in the insets. (f) An
intensity line profile taken from (c). From these line profiles, the simulated
data suggests that intensity deviations from a uniform contrast profile
resulting from atomic electric fields are significantly minimized in the
annular geometry.

Furthermore, as discussed by [31], it is apparent that low spatial frequency
information associated with long-range fields leads to a uniform displacement
of the CBED pattern.[31] Meanwhile, high spatial frequency information
associated with short-range fields leads to a redistribution in intensity
within the CBED pattern.[31] As such, one method to isolate information
associated with long-range fields is through the use of an annular DPC
detector.[32] Namely, by filtering out low frequency signal with an inner
collection angle approaching the convergence semi-angle of the CBED pattern
(Figure S5), the redistribution in signal associated with short-range
potentials can be largely diminished. Simulated differential phase contrast
images of hBN subjected to an external electric field demonstrate this aspect
(Figure 2). In the case of a traditional DPC image (Figure 2a-b), local atomic
electrostatic potentials and long-range external potentials can both
contribute to the measured signal. In annular DPC images, however, (parameters
listed in Methods) the contribution of local potentials can be selectively
diminished (line profiles provided in Figure 2c & f), which makes this method
potentially very useful for mapping long-range electric fields.

## V RESULTS AND DISCUSSION

Figure 3: Experimental scanning diffraction data acquired from the sample seen
in center dark field TEM image. (a) CBED patterns captured at point 1 marked
on the hBN sample. Colorbar intensity is in units of relative detection
intensity. CBED patterns taken as a function of indicated positions when (b)
0V and (c) 5V are applied to the hBN sample. The purple dots refer to the
center of mass of the CBED pattern taken from regions devoid of sample, while
the yellow dots represent the center of mass of the CBED pattern taken from a
particular position. The border color represents the degree of probe
deflection in the x-direction.(d) CBED patterns captured at point 1 marked on
the MoS2/hBN sample. Colorbar intensity is in units of relative detection
intensity. CBED patterns taken as a function of indicated positions when (e)
0V and (f) 5V are applied to the MoS2/hBN sample. Scale bars represent 2 mrad.
A clear shift in center of mass is observed when an external field is applied
for both samples.

Experimental results displaying ICOM values measured in a layer of hBN prior
to applying and after applying external bias are presented in Figure 3. Prior
to applying an external bias, the CBED pattern taken from spot 1 on the TEM
image is provided in Figure 3a. In order to better observe the beam
deflections present in these diffraction patterns, the center portion of
individual CBED patterns is provided in Figure 3b. In these electron
diffraction patterns, the green dot represents the ICOM of the vacuum probe
that has not interacted with the sample. The yellow dots represent the ICOM of
the CBED pattern taken from the specified positions of the hBN layer. The ICOM
values are calculated using FM-STEM. When this measurement is repeated in the
case where an external voltage of 5V is applied and an electric field is
generated along the direction of the indicated arrow (Figure 3c), larger
deflections in the ICOM values compared to the vacuum probe are detected. This
suggests that this external stimulation leads to the incident electron probe
experiencing additional momentum transfer.

To create a more electronically-relevant structure, a monolayer of MoS2 is
placed on top of the hBN layers. This semiconducting layer provides an
alternative current pathway and this architecture represents a tunneling
contact geometry previously explored for reducing contact resistance.[33, 34]
Once more, it is apparent that applying an external voltage leads to greater
beam deflection in this geometry (Figure 3d-f).

Eq. 2 is applied to relate shifts in the ICOM of the CBED pattern to the
electric field distribution. In the absence of an external electric field,
these fields are seen in Figure 4. The quasi-random pattern observed (Figure
4a) likely represents electric fields arising from atomic orbitals that are
captured by the uncorrected nanoscale electron probe.[27]

In order to calculate the electric field distribution when an external voltage
of 5V is applied, electrons detected at each probe position are filtered using
an annular, quadrant DPC virtual detector (Methods). Once again Eq. 2 is used
to calculate the electric field vector at each probe position. This map of the
electric field distribution when an external voltage of 5V is applied is
provided in Figure 4b. This calculated field using the annular detector
architecture agrees with our expectations for the uniform field distribution
that would arise when an electrical potential is applied across hBN.

For the MoS2/hBN sample, similar methods are used to characterize the
resultant electrical fields in the absence of an external voltage and the
long-range fields that arise when an external potential is applied. These,
too, agree with prior expectations and are indicative of a uniform long-range
field existing across the sample (Figure 4c & d).

Figure 4: Experimental electric field maps. In the absence of an applied bias,
electric field maps of the (a) hBN and (c) MoS2/hBN samples are provided.
Similar field maps in the case where an external bias of 5V is applied are
provided in (b) and (d) for the hBN and MoS2/hBN samples, respectively. In the
case of an applied field, the measured electric field vectors align with the
external potential gradient.

Although based on the geometry of the system (5V applied across 1$\mu$m
region), we expect the external field to have a magnitude on the order of 1
V/$\mu$m, the calculated field magnitude is found to be on the order of 1
V/nm. This discrepancy can potentially be explained by the fact that this
measured electric field includes contributions from the electrostatic field in
the sample plane as well as the leakage electric fields in the vacuum regions
directly above and beneath the sample. We model the relative intensity of this
leakage field by solving Laplace’s equation for the sample geometry as seen in
Figure S8. These calculations suggest that while the accumulated leakage field
scales with sample field, they can persist for several hundred nanometers both
above and beneath the sample. As such, rigorous analysis of the size of the
electric field would be needed to quantitatively extract the precise electric
field magnitudes within the sample, which is outside the scope of this study.

Finally, since Maxwell’s theory states that the divergence of the electric
field is proportional to the charge density (Methods), a local charge density
distribution of this MoS2/hBN region under applied electric field is provided
in Figure S9. From this figure, it is clear that a carrier density roughly
equivalent to 1-2 x $10^{13}cm^{-2}$ is observed when an external voltage is
applied. These values are reasonable for monolayer MoS2 and at such carrier
densities, the material is expected to exist in a conductive regime with a 4
probe resistivity $<$ 100 k$\Omega$/sq.[35]

In summary, we demonstrate the possibility for measuring the electrostatic
field distribution that arises in 2D materials as a function of external bias.
This technique involves the use of an annular DPC mask to relate shifts in the
center of mass of the diffraction pattern as a function of position to local
fields. Through simulated and experimental investigations, we determine that
such masks provide a method for selectively probing the low spatial frequency
information associated with long-range electric fields. This work demonstrates
that these tools offer both a route to identifying non-uniformities in
electric field as well as potential sources of noise or performance
degradation during device operation in 2D architectures moving forward.

## VI ACKNOWLEDGEMENTS

This material is based upon work supported by the National Science Foundation
under Grant No. DMR-1929356. This work made use of the EPIC, Keck-II, and SPID
facilities of Northwestern University’s NUANCE Center, which has received
support from the SHyNE Resource (NSF ECCS-2025633), the IIN, and
Northwestern’s MRSEC program (NSF DMR-1720139). A.A.M. gratefully acknowledges
support from the Ryan Fellowship and the IIN at Northwestern University.
T.K.S. was supported by the Office of Naval Research (N00014-16-1-3055). P.L.
was supported by Argonne National Laboratory. K.W. and T.T. acknowledge
support from the Elemental Strategy Initiative conducted by the MEXT, Japan,
Grant Number JPMXP0112101001, JSPS KAKENHI Grant Number JP20H00354 and the
CREST (JPMJCR15F3), JST. The authors thank Erik Lenferink for assistance with
sample characterization. Research reported in this publication was supported
in part by instrumentation provided by the Office of The Director, National
Institutes of Health of the National Institutes of Health under Award Number
S10OD026871. The content is solely the responsibility of the authors and does
not necessarily represent the official views of the National Institutes of
Health. This research was supported in part through the computational
resources and staff contributions provided for the Quest high performance
computing facility at Northwestern University which is jointly supported by
the Office of the Provost, the Office for Research, and Northwestern
University Information Technology.

## VII AUTHOR CONTRIBUTIONS

The manuscript was written through contributions of all authors. All authors
have given approval to the final version of the manuscript.

## VIII ADDITIONAL INFORMATION

Supplementary Information Supplementary Figs. 1–11.
Competing Financial Interests: The authors declare no competing financial
interests.

## IX METHODS

### A Material Synthesis

Single crystal flakes of monolayer MoS2 were grown through an atmospheric
pressure chemical vapor deposition process (CVD). This process involved the
use of 5 mg of MoO3 (Sigma-Aldrich) powder, which is spread evenly throughout
an alumina boat and capped by two face-down, 1 cm by 1 cm SiO2/Si substrates
(300 nm oxide). This boat is placed within a 1-inch diameter quartz tube at
the center of a tube furnace. Similarly, 300 mg of sulfur pieces (Alfa-Aesar)
were placed in another alumina boat located upstream and outside of the
furnace. Following a 15 min purge of the sealed quartz tube with argon gas,
the center of the furnace is heated to 700∘C over a period of 40 minutes and
held at that temperature for 3 minutes. During this process, the flow rate of
argon gas is maintained at 13 sccm. When the temperature of the furnace
reached 575∘C during this ramp step, the alumina boat containing the sulfur is
moved to a region inside the furnace where the temperature is approximately
150∘C. This is done with the help of a magnet. Following the growth process,
the furnace is allowed to cool naturally.

#### A.1 DPC STEM Simulations

Multislice simulations were performed using custom MATLAB® codes following the
methods laid out by Kirkland.[36] These simulations used similar microscope
parameters as those used in acquiring the experimental datasets along with
eight frozen phonon configurations. An external potential is applied by adding
a mesh grid with linearly varying potentials from 0 to 10kV to the first layer
of the multislice calculation. Traditional DPC images were acquired by
processing the data using the virtual quadrant detectors seen in Figure 2.
Electron signal in quadrants 1 and 3 were subtracted from one another to
acquire DPCx images and electron signal in quadrants 2 and 4 were subtracted
from one another to acquire DPCy images. For annular DPC images, the data is
processed using the virtual annular quadrant detectors seen in Figure 2 with
an inner collection semi-angle of 28 mrad and an outer collection semi-angle
of 32 mrad. DPCx and DPCy images were acquired through similar means.

#### A.2 Scanning Diffraction Image Processing

4D STEM datasets were processed using custom Python scripts. For analysis
associated with samples under no applied bias, the x and y components of the
intensity center of mass at each probe position (ICOM) were calculated
directly. For analysis associated with samples under an applied bias,
transmitted electron signal both prior to applying and after applying applied
bias is filtered using an annular quadrant detector. This virtual DPC detector
is constructed to have an inner collection semi-angle of 28 mrad and an outer
collection semi-angle of 32 mrad. In this case, DPC signal obtained under no
applied bias is used as the background image and subsequently subtracted from
the DPC signal obtained under applied bias to further eliminate the impact of
short-range atomic potentials. The net signal in the x and y annular quadrants
is used to approximate the x and y components of the ICOM.

Following the protocol detailed by Fang et al. and Müller et al,[27, 37]
electric field maps were calculated using the x and y components of the ICOM
values. These values were related to the expectation value of the momentum
transfer, $\left<p_{\bot}\right>$, through Eq. 1. $\left<p_{\bot}\right>$ is
calculated in units of eV s Å-1. Eq. 2 is used to relate the expectation value
of the momentum transfer to electric field. $E_{\bot}$ is represented in the
text using a quiver plot displaying both the $E_{x}$ and $E_{y}$ components
and represented in units of V Å-1.

The total charge density ($\rho$), which is represented in units of elementary
charge per Å2, is calculated using Eq. 3. In this equation,$\epsilon_{0}$
represents the vacuum permittivity.

$div(E_{\bot}z)=\frac{\rho}{\epsilon_{0}}$ (3)

#### A.3 Confocal Raman spectroscopy

Raman spectra were obtained using a Horiba LabRAM HR Evolution Confocal Raman
system under laser illumination with a wavelength of 532 nm and a power of
0.25mW.

## References

  * [1] Byungjin Cho et al.  “Chemical Sensing of 2D Graphene/MoS2 Heterostructure Device”  In _ACS App. Mater. Interfaces_ 7.30, 2015, pp. 16775–16780
  * [2] Yang Yang et al.  “Edge-Oriented MoS2 Nanoporous Films as Flexible Electrodes for Hydrogen Evolution Reactions and Supercapacitor Devices”  In _Adv. Mater._ 26.48, 2014, pp. 8163–8168
  * [3] Jie Jiang et al.  “2D MoS2 Neuromorphic Devices for Brain-Like Computational Systems”  In _Small_ 13.29, 2017, pp. 1700933
  * [4] Yunqiu Kelly Luo et al.  “Opto-Valleytronic Spin Injection in Monolayer MoS2/Few-Layer Graphene Hybrid Spin Valves”  In _Nano Lett._ 17.6, 2017, pp. 3877–3883
  * [5] Ruijing Ge et al.  “Atomristor: Nonvolatile Resistance Switching in Atomic Sheets of Transition Metal Dichalcogenides”  In _Nano Lett._ 18.1, 2018, pp. 434–441
  * [6] V.. Calado et al.  “Ballistic Josephson Junctions in Edge-Contacted Graphene”  In _Nat. Nanotechnol._ 10.9, 2015, pp. 761–4
  * [7] M. Kim et al.  “Strong Proximity Josephson Coupling in Vertically Stacked NbSe2-Graphene-NbSe2 Van Der Waals Junctions”  In _Nano Lett._ 17.10, 2017, pp. 6125–6130
  * [8] Feng Wang et al.  “Tunable GaTe-MoS2 Van Der Waals P–N Junctions With Novel Optoelectronic Performance”  In _Nano Lett._ 15.11, 2015, pp. 7558–7566
  * [9] Wenzhuo Wu et al.  “Piezophototronic Effect in Single-Atomic-Layer MoS2 for Strain-Gated Flexible Optoelectronics”  In _Adv. Mater._ 28.38, 2016, pp. 8463–8468
  * [10] Min Sup Choi et al.  “Lateral MoS2 P–N Junction Formed by Chemical Doping for Use in High-Performance Optoelectronics”  In _ACS Nano_ 8.9, 2014, pp. 9332–9340
  * [11] Dattatray J Late et al.  “Hysteresis in Single-Layer MoS2 Field Effect Transistors”  In _ACS Nano_ 6.6, 2012, pp. 5635–5641
  * [12] Jingli Wang et al.  “High Mobility MoS2 Transistor With Low Schottky Barrier Contact by Using Atomic Thick H‐BN as a Tunneling Layer”  In _Adv. Mater._ 28.37, 2016, pp. 8302–8308
  * [13] Xiaochi Liu et al.  “P-Type Polar Transition of Chemically Doped Multilayer MoS2 Transistor”  In _Adv. Mater._ 28.12, 2016, pp. 2345–2351
  * [14] A.. Murthy et al.  “Direct Visualization of Electric-Field-Induced Structural Dynamics in Monolayer Transition Metal Dichalcogenides”  In _ACS Nano_ 14.2, 2020, pp. 1569–1576
  * [15] H. Rose  “Nonstandard Imaging Methods in Electron Microscopy”  In _Ultramicroscopy_ 2, 1976, pp. 251–267
  * [16] H. Rose  “Image Formation by Inelastically Scattered Electrons in Electron Microscopy”  In _Optik_ 45.2, 1976, pp. 139–158
  * [17] N.. Dekkers and H.. Lang  “Differential Phase-Contrast in a STEM”  In _Optik_ 41.4, 1974, pp. 452–456
  * [18] N. Shibata et al.  “Differential Phase-Contrast Microscopy at Atomic Resolution”  In _Nat. Phys._ 8.8, 2012, pp. 611–615
  * [19] J.. Chapman, P.. Batson, E.. Waddell and R.. Ferrier  “The Direct Determination of Magnetic Domain Wall Profiles by Differential Phase Contrast Electron Microscopy”  In _Ultramicroscopy_ 3.2, 1978, pp. 203–14
  * [20] E.. Waddell and J.. Chapman  “Linear Imaging of Strong Phase Objects Using Asymmetrical Detectors in STEM”  In _Optik_ 54.2, 1979, pp. 452–456
  * [21] R. Close, Z. Chen, N. Shibata and S.. Findlay  “Towards Quantitative, Atomic-Resolution Reconstruction of the Electrostatic Potential via Differential Phase Contrast Using Electrons”  In _Ultramicroscopy_ 159, 2015, pp. 124–137
  * [22] Marco Battaglia, Devis Contarato, Peter Denes and Piero Giubilato  “Cluster Imaging With a Direct Detection CMOS Pixel Sensor in Transmission Electron Microscopy”  In _Nucl. Instrum. Methods Phys. Res., Sect. A_ 608.2, 2009, pp. 363–365
  * [23] K. Müller-Caspary et al.  “Comparison of First Moment STEM With Conventional Differential Phase Contrast and the Dependence on Electron Dose”  In _Ultramicroscopy_ 203, 2019, pp. 95–104
  * [24] Ivan Lazić, Eric G.. Bosch and Sorin Lazar  “Phase Contrast STEM for Thin Samples: Integrated Differential Phase Contrast”  In _Ultramicroscopy_ 160, 2016, pp. 265–280
  * [25] Colin Ophus  “Four-Dimensional Scanning Transmission Electron Microscopy (4d-Stem): From Scanning Nanodiffraction to Ptychography and Beyond”  In _Microsc. Microanal._ 25.3, 2019, pp. 563–582
  * [26] Ch Androulidakis et al.  “Strained Hexagonal Boron Nitride: Phonon Shift and Grüneisen Parameter”  In _Phys. Rev. B_ 97.24, 2018, pp. 241414
  * [27] Knut Müller et al.  “Atomic Electric Fields Revealed by a Quantum Mechanical Approach to Electron Picodiffraction”  In _Nat. Commun._ 5.1, 2014, pp. 5653
  * [28] Takehito Seki et al.  “Quantitative electric field mapping in thin specimens using a segmented detector: Revisiting the transfer function for differential phase contrast”  In _Ultramicroscopy_ 182, 2017, pp. 258–263  DOI: https://doi.org/10.1016/j.ultramic.2017.07.013
  * [29] Knut Müller-Caspary et al.  “Measurement of atomic electric fields and charge densities from average momentum transfers using scanning transmission electron microscopy” FEMMS 2015  In _Ultramicroscopy_ 178, 2017, pp. 62–80  DOI: https://doi.org/10.1016/j.ultramic.2016.05.004
  * [30] Ian MacLaren et al.  “On the origin of differential phase contrast at a locally charged and globally charge-compensated domain boundary in a polar-ordered material”  In _Ultramicroscopy_ 154, 2015, pp. 57–63  DOI: https://doi.org/10.1016/j.ultramic.2015.03.016
  * [31] Michael C Cao et al.  “Theory and practice of electron diffraction from single atoms and extended objects using an EMPAD”  In _Microscopy_ 67.suppl_1, 2017, pp. i150–i161  DOI: 10.1093/jmicro/dfx123
  * [32] J.. Chapman, I.. McFadyen and S. McVitie  “Modified differential phase contrast Lorentz microscopy for improved imaging of magnetic structures”  In _IEEE Transactions on Magnetics_ 26.5, 1990, pp. 1506–1511  DOI: 10.1109/20.104427
  * [33] A.. Murthy et al.  “Intrinsic Transport in 2D Heterostructures Mediated Through H-Bn Tunneling Contacts”  In _Nano Lett._ 18.5, 2018, pp. 2990–2998
  * [34] U. Chandni, K. Watanabe, T. Taniguchi and J.. Eisenstein  “Evidence for Defect-Mediated Tunneling in Hexagonal Boron Nitride-Based Junctions”  In _Nano Lett_ 15.11, 2015, pp. 7329–33
  * [35] X. Cui et al.  “Low-Temperature Ohmic Contact to Monolayer MoS2 by Van Der Waals Bonded Co/H-Bn Electrodes”  In _Nano Lett._ 17.8, 2017, pp. 4781–4786
  * [36] Earl J. Kirkland  “Advanced Computing in Electron Microscopy”  Springer, 2020
  * [37] Shiang Fang et al.  “Atomic Electrostatic Maps of 1D Channels in 2D Semiconductors Using 4D Scanning Transmission Electron Microscopy”  In _Nat. Commun._ 10.1, 2019, pp. 1127

Supplementary Information
Spatial Mapping of Electrostatics and Dynamics across 2D Heterostructures

Akshay A. Murthy1,2, Stephanie M. Ribet1, Teodor K. Stanev3, Pufan Liu1,
Kenji Watanabe4, Takashi Taniguchi5, Nathaniel P. Stern3, Roberto dos
Reis1,6,∗,
Vinayak P. Dravid1,2,6,∗

1Department of Materials Science and Engineering, Northwestern University,
Evanston, IL, USA
2International Institute of Nanotechnology, Northwestern University, Evanston,
IL, USA
3Department of Physics and Astronomy, Northwestern University, Evanston, IL,
USA
4Research Center for Functional Materials, National Institute for Materials
Science, 1-1 Namiki, Tsukuba 305-0044, Japan
5International Center for Materials Nanoarchitectonics, National Institute for
Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan
6The NUANCE Center, Northwestern University, Evanston, IL, USA

∗Corresponding Authors:
Roberto dos Reis<EMAIL_ADDRESS>
Vinayak P. Dravid<EMAIL_ADDRESS>

KEYWORDS: In situ electron microscopy, heterostructure, transition metal
dichalcogenides, tunneling contacts, MoS2, differential phase contrast

Table of contents:
S1: Height profile of thin hBN support layers
S2: Spectroscopy and I-V results from prepared TEM sample
S3: Collapse of suspended TEM sample fabricated through alternative means
S4: hBN structure used for DPC simulations
S5: Comparison of inner collection semi-angles for annular DPC detector
S6: TEM diffraction pattern taken from MoS2/hBN sample
S7: Variation in CBED with and without applied field
S8: Electric field simulation of TEM sample to understand impact of leakage
field
S9: Charge Density across MoS2/hBN

S1: Height profile of thin hBN support layers

Figure S1: Atomic force microscopy image taken from a flake of mechanically
exfoliated hBN used as the support substrate. The flake thickness was found to
be roughly 5 nm.

S2: Spectroscopy and I-V results from prepared TEM sample

Figure S2: (a) Colorized SEM image of suspended 2D heterostructure geometry.
The gold region is indicative of gold electrodes. The purple region indicates
the presence of h-BN, and the green region represents monolayer MoS2, (b)
Current (Ids) - Voltage (Vds) characteristics taken from sample indicate that
an electric field can be applied across the sample, (c) Photoluminescence (PL)
spectra taken from MoS2/h-BN regions that are supported on SiO2, supported on
gold, and suspended. Raman spectrum taken from the same (d) MoS2 and (e) h-BN
regions are also provided. Due to the inherent differences in charge doping
and strain, the position and magnitude of the PL peaks and the position of
Raman modes vary with substrate.

S3: Collapse of suspended TEM sample fabricated through alternative means

Figure S3: Various challenges arise when trying to construct an in situ TEM
sample of a 2D heterostructure when using other methods. In this case, the
electrodes were first transferred onto the TEM windows and 2D layers were
sequentially transferred on top. Unfortunately, the repeated thermal
processing involved in this setup prevents successful material suspension and
creates cracking and fracturing of the metal electrodes.

S4: hBN structure used for DPC simulations

Figure S4: (a) Projected potential of structure used for multislice
simulation. (b) External potential applied to simulated structure. Simulated
dark field images taken from structure with (c) and without (d) external
potential. Average CBED taken from structure with (e) and without (f) external
potential. The scale bar is consistent for the real space images in a-d as
well as the diffraction space images in e-f

S5: Comparison of inner collection semi-angles for annular DPC detector

Figure S5: (a) Normalized intensity line profiles taken from same region in
Figure 2 for inner collection semi-angles of 0 mrad, 15 mrad, 25 mrad, 28
mrad, and 30 mrad. (b) Normalized intensity lines profiles taken from same
region in Figure 2 replotted for inner collection semi-angles of 25 mrad, 28
mrad, and 30 mrad. Minimal intensity variation from atomic electrostatic
fields is observed for both 28 mrad and 30 mrad inner collection angles. 28
mrad was selected to provide sufficient latitude in terms of beam and detector
alignment in practical situations.

S6: TEM diffraction pattern taken from sample

Figure S6: TEM selected area diffraction pattern taken from MoS2/hBN sample
indicates that the two layers are present and offset by 5∘.

S7: Variation in CBED with and without applied field

Figure S7: Difference in experimental CBED patterns taken from MoS2/hBN region
at position 1 on Figure 3 between Vds=0V and Vds=5V conditions. This colormap
indicates that when the probe is smaller than the feature, the ICOM signal is
measured from the shifts of the unscattered beam as opposed to asymmetries
within the disc.

S8: Electric field simulation of TEM sample to understand impact of leakage
field

Figure S8: Simulated relative electric field intensities in the sample as well
as a 1$\mu$m vacuum region directly above and below the sample plane. For this
simulation, boundary condition values of 0V at the ground electrode and a
positive voltage at the source electrode (Vds) were used. The Laplace equation
was solved to determine the electric potential at each finite location. By
taking the in-plane spatial gradient of these values, the in-plane electric
field distribution was calculated.

S9: Charge Density across MoS2/hBN

Figure S9: Charge density map across MoS2/hBN region calculated from electric
field map in Figure 4d.


# Tunable Degenerate Optical Parametric Oscillation with Coupled
Microresonators

Nathalia B. Tomazio*,1,2,3  Luca O. Trinchão*,1,2  Eduardo S. Gonçalves*,1,2
Laís Fujii dos Santos *,1,4  Paulo F. Jarschel1,2  Felipe G. S. Santos1,2
Thiago P. Mayer Alegre1,2  Gustavo S. Wiederhecker1,2,${\dagger}$ * These
authors contributed equally to this work.
1 Photonics Research Center, Universidade Estadual de Campinas (UNICAMP),
Campinas, SP, Brazil
2 Gleb Wataghin Physics Institute, Universidade Estadual de Campinas
(UNICAMP), Campinas, SP, Brazil 3 Instituto de Física, Universidade de São
Paulo, São Paulo, Brazil
4 School of Electrical Engineering and Computer Science, University of Ottawa,
K1N 6N5, ON, Canada
<EMAIL_ADDRESS>
††journal: opticajournal††articletype: Research Article

Microresonator-based degenerate optical parametric oscillation (DOPO) has
recently been explored as a compelling platform for all-optical computing and
quantum information applications, such as truly random number generation and
the production of squeezed states of light. Emerging research has highlighted
the potential of coupled microresonators, or photonic molecules, as a novel
avenue for spectral engineering, unlocking an extra degree of freedom for the
optimization of four-wave mixing interactions. Here, we demonstrate DOPO
within the coupled modes of a silicon nitride triple-state photonic molecule.
Our design introduces a distinctive mechanism for spectral engineering, using
microheaters to individually tune the resonance spectral positions, thus
enabling dynamic local dispersion control within the coupled modes. We
successfully generate a DOPO signal with active efficiency control and explore
the optical mode spacing in the tens of gigahertz range to use native phase-
locked optical pumps driven by a radio-frequency source.

##  1 Introduction

Degenerate optical parametric oscillation (DOPO) has been investigated as a
novel strategy for various applications, including all-optical quantum random
number generation (RNG) [1, 2, 3], generation of squeezed states of light [4,
5] and coherent optical computing [6, 7, 8]. DOPO occurs when the signal and
idler fields generated by parametric amplification in an optical cavity are
degenerate in frequency. In $\chi^{(2)}$ and $\chi^{(3)}$-media, DOPO is
generated, respectively, via parametric down-conversion and four-wave mixing
(FWM) with a dual-pump scheme. In both cases, the DOPO phase properties arise
from the underlying phase-sensitive parametric amplification [9, 1]. Below the
oscillation threshold, the phase state of the signal and idler pair is
described by a squeezed vacuum state [4, 5]. Upon reaching the oscillation
threshold – where parametric gain exceeds losses – DOPO undergoes a
nonequilibrium phase transition to a coherent binary-phase state, which is
either in-phase or $\pi$ out-of-phase relative to the pump fields [1, 10]. In
this case, both phase states are equally likely, since oscillation is
initiated from vacuum fluctuations. This random binary phase state provides a
means to truly RNG [1] and is used to emulate spin states on photonic Ising
machines [8]. Recently, CMOS-compatible integrated microresonators based on
the $\chi^{(3)}$ nonlinearity have enabled the development of on-chip DOPO-
based applications [11, 2, 3, 7, 4, 5], offering scalability, compactness,
enhanced dispersion control, and high Q-factors, enabling efficient nonlinear
effects even at low-power CW pumping.

Near-field coupled microresonators, also known as photonic molecules, offer
extra degrees of freedom for spectral engineering[12, 13, 14, 15, 16, 17, 18,
19, 20]. These systems exploit evanescent coupling to induce hybridization of
the bare resonator modes, resulting in frequency-split coupled modes
(supermodes) that locally modify the dispersion. By coupling two resonators,
several authors have harnessed controllable avoided mode-crossings to engineer
phase-matching conditions for FWM, allowing generation of optical frequency
combs in the normal dispersion regime [21, 22, 23], purity enhancement and
preservation of squeezed states [5], and high-efficiency parametric conversion
[24]. Despite significant progress in this direction, to the best of our
knowledge, no nanophotonic device has yet demonstrated simultaneous tunability
for dispersion control and isolation of the target FWM interaction from
competing processes. Additionally, for dual-pumped FWM interactions, such as
DOPO, the participating modes of single or two coupled resonator platforms are
spaced by free spectral range (FSR) units, ranging from tenths to units of
THz. This wide frequency spacing inevitably demands two lasers for the dual-
pump scheme, which can be experimentally challenging, especially for those
applications requiring phase-locked pumps.

Figure 1: Photonic molecule design: (a) Photonic molecule design comprising an
array of three identical Si3N4 ring resonators with $50~{}\upmu$m of radii
(FSR $\sim 450$ GHz), a cross-section of $1~{}\upmu$m $\times\,0.8~{}\upmu$m
(width $\times$ height) and an inter-ring gap of 550 nm, cladded in a SiO2
substrate. Light is launched at the coupled microrings through the bus
waveguide (waveguide-to-ring gap of 600 nm) and can be collected via either
the through or the drop port. Inverted tapers at both the input and output
ports optimize the TE coupling. Integrated microheaters atop each ring allow
for spectral tunability. (b,c) Third-order parametric nonlinear processes that
occur in a dual-pumped (a) single-ring microcavity and (b) a three-coupled
ring photonic molecule. When phase-matching is achieved, DOPO is generated at
the central resonance between the two pumps. In (b), the normal dispersion
regime hinders the undesired ND-FWM interactions by phase-mismatching
supermodes at different mode numbers. (d) Measured linear transmission at the
through-port of the photonic molecule, exhibiting evenly spaced triplets
occurring every free-spectral range of the bare rings. The colored Lorentzian
dips indicate the Antisymmetric (AS, blue), Central (C, black), and Symmetric
(S, red) supermodes. The loaded quality factor of the C supermode was measured
to be roughly $Q_{L}=2.2\times 10^{5}$. (e) 2D simulation of the electric
field profile of each supermode using the finite-element method. This
simulation uses a representative ring radius for better visualization. (f)
Frequency splitting of supermodes as a function of the inter-ring gap for
identical rings, calculated via CMT [25]. The line serves as a guide for the
eye.

In this work, we propose an array of three identical coupled SiN
microresonators as a platform for DOPO, illustrated in Fig. 1(a). The mode
hybridization induced by the coupling yields a triplet of coupled modes that
repeats every FSR of the bare rings, so DOPO can be fully encompassed within
one triplet of coupled modes without relying on modes of other azimuthal
orders. For comparison, in single-resonator DOPOs, achieving phase-matching
for the three interacting modes often inadvertently phase-matches other
azimuthal modes as well, leading to parasitic and cascading non-degenerate FWM
(ND-FWM) processes, as illustrated in Fig. 1(b). In contrast, our photonic
molecule design can potentially isolate the DOPO from these competing effects.
This is due to the wide spectral separation between triplets, allowing local
phase-matching within a triplet to occur independently from phase-matching
across other triplets (Fig. 1(c)). By tuning the individual phases of the
coupled modes with integrated microheaters, we unlock locally tunable
dispersion control, enabling us to dynamically correct for phase mismatch and
optimize third-order nonlinear optical interactions within a triplet of
coupled modes. Moreover, by setting up inter-resonator gaps of a few hundred
nanometers, frequency splittings of tens of GHz can be achieved among the
coupled modes, without increasing cavity length [26, 27]. Such radio-frequency
(RF) range splittings allow for phase-coherent dual-tone pumping from a single
laser source using electro-optic modulation (EOM). This strategy confers
native phase-locking between the pump fields, considerably reducing the
complexity of the experimental setup for DOPO.

##  2 Design considerations

Our photonic molecule comprises an array of three identical (degenerate) SiN
microring resonators, fabricated by Ligentec SA, illustrated in Fig. 1(a).
Near-field evanescent coupling between the rings induces hybridization of
their individual modes, resulting in coupled modes that undergo frequency
splitting from their original uncoupled frequencies. This results in evenly
spaced triplets occurring in every FSR of the bare rings, shown in the
through-port transmission in Fig. 1(d). Moreover, integrated microheaters
placed atop each microring enable dynamic tuning of the individual resonant
frequencies via the thermo-optic effect. This capability allows for precise
control over the coupling dynamics, enabling us to spectrally engineer the
supermode splittings and locally tailor the group velocity dispersion (GVD)
within a set of supermodes.

Figure 2: GVD control and FWM: (a-c) Residual dispersion
($D_{int}=\omega_{j,\mu}-\omega_{j,0}-D_{1,j}\mu$) for the S (a), C (b), and
AS (c) supermodes, revealing their normal GVD regimes. Circles represent
measured dispersion; line curves are cubic fittings. (d) Theoretical
calculation of the S (red), C (black), and AS (blue) supermodes relative
frequencies as a function of a frequency deviation of ring (1). The diagonal
dashed curve represents the relative frequency deviation
$(\delta\omega-\omega_{\mathrm{C}})/2\pi$. (e) sGVD ($J_{2}$) calculated from
(d), indicating the transition from the normal to the anomalous local
supermode dispersion regime. The green hatched region in both (d, e) delimits
the parameter region explored in our experiments depicted in (g, h). (f)
Optical spectrum of FWM. The involved supermodes are colored following the
color convention (S, red; C, black; AS, blue), while the insert indicates the
energy conservation of the nonlinear process. A signal-to-idler conversion
efficiency of -20.4 dB is displayed. (g) Relative spectral position of
supermodes collected at the through port as a function of the current applied
to the microheater atop ring (1). The white dashed lines represent our
theoretical model fitted to the experimental data. (h) sGVD ($J_{2}$)
calculated from (g), indicating the transition from the normal to the
anomalous local supermode dispersion regime. The black solid curve indicates
the fitted model of our calculation in (e) and the color-coded scatter points
indicate the idler field power achieved in the stimulated FWM experiment.

The supermode frequency splitting can be determined by solving the coupled
amplitude basis $\\{a_{j}\\}$ and its associated eigenfrequencies $\omega_{j}$
via coupled mode theory in time (CMT) [25]:

$\displaystyle a_{\mathrm{AS}}=\frac{1}{2}(1~{}-\sqrt{2}~{}1)^{T},$
$\displaystyle\omega_{\mathrm{AS}}=\omega_{0}+\sqrt{2}J,$ (1a) $\displaystyle
a_{\mathrm{C}}=\frac{1}{\sqrt{2}}(-1~{}0~{}1)^{T},$
$\displaystyle\omega_{\mathrm{C}}=\omega_{0},$ (1b) $\displaystyle
a_{\mathrm{S}}=\frac{1}{2}(1~{}\sqrt{2}~{}1)^{T},$
$\displaystyle\omega_{\mathrm{S}}=\omega_{0}-\sqrt{2}J.$ (1c)

The supermode amplitudes are normalized such that $|a_{j}|^{2}$ corresponds to
the $j$-mode energy, for which $j=$ AS, C and S refer to the antisymmetric,
central, and symmetric supermodes, respectively. The electric field
distribution of the supermodes, simulated by the finite element method (COMSOL
Multiphysics), is depicted in Fig. 1(e). In the $a_{\mathrm{AS}}$ and
$a_{\mathrm{A}}$ supermodes, the electric field displays opposite phases
within the central resonator, which leads to their designation as
antisymmetric (AS) and symmetric (S) modes: in the AS supermode, the electric
field in the central microring is out-of-phase with respect to the outer
counterparts, while in the S supermode, the electric field in the three
microrings is in phase. The central supermode $a_{\mathrm{C}}$ mostly occupies
the outer rings and undergoes minimal frequency shifts due to coupling [14].
Such modal distributions have been shown to affect the FWM dynamics in coupled
resonator systems [28] (see Supplementary Material Section S1 for details).
Moreover, these supermodes deviate from the original uncoupled degenerate
resonance frequency $\omega_{0}$ by the factor $\sqrt{2}J$, where $J/2\pi$
represents the ring-to-ring coupling rate, adjustable via the inter-ring gap.
The frequency splitting between the supermodes is presented in Fig. 1(f) as a
function of the inter-ring gap distance. Our devices exhibit an inter-ring gap
of 550 nm, which yields a radio-frequency (RF) range splitting ($\sqrt{2}J\sim
10~{}\text{GHz}$). The calculation details of the inter-ring coupling
coefficient can be found in the Supplementary Material Section S2.

The supermodes’ rich dispersive dynamics can be described by the resonator
optical frequency expansion [29],

$\omega_{j,\mu}=\omega_{j,0}+D_{1,j}\mu+\frac{1}{2}D_{2,j}\mu^{2}+\dots$ (2)

where $\mu=m-m_{0}$ is the azimuthal mode number relative to a reference mode
($m_{0}$) with frequency $\omega_{j,0}$, and the coefficients $D_{1,j}$,
$D_{2,j}$ represent the FSR and the change in the FSR with respect to $\mu$
(GVD parameter), respectively. The index $j$ is included to account for the
supermode numbers ($j$ = AS, C, or S). Our device exhibits normal GVD regimes
($D_{2,j}<0$) as shown in Fig. 2(a-c), with $D_{2}/2\pi=-8.3$ (S), $-7.6$ (C)
and $-6.7$ (AS) MHz, at $\lambda_{0}\approx 1560$ nm. The GVD varies slightly
for each supermode and higher-order dispersion terms have been neglected. For
a given triplet (fixed mode number $\mu=\mu_{0}$), the supermode dispersion is
described by defining the supermode free spectral range (sFSR, $J_{1}$) as the
average frequency splitting and the supermode group velocity dispersion (sGVD,
$J_{2}$) as its variation within the triplet, similar to standard FSR and GVD
definitions across mode numbers [17] (for details on this treatment, see
Supplementary Information Section S4). In triple-state photonic molecules, up
to second-order supermode dispersion is defined: the $J_{2}$ sign indicates
whether the supermodes’ frequency spacing increases ($J_{2}>0$) or decreases
($J_{2}<0$) with frequency, describing an anomalous and normal local
dispersion regime, respectively. It can be calculated straightforwardly as
$J_{2}\equiv(\omega_{\mathrm{AS}}-\omega_{\mathrm{C}})-(\omega_{\mathrm{C}}-\omega_{\mathrm{S}})$.

While Eq. 1 holds for the case of three identical bare rings
($\omega^{(1)}=\omega^{(2)}=\omega^{(3)}=\omega_{0}$), individual frequency
deviations from the degenerate case can produce large spectral asymmetry among
the triplet of supermodes, modifying the coupled basis solution. In Fig. 2(d),
we show the supermode frequencies as a function of a deviation $\delta\omega$
in the eigenfrequency of the ring coupled to the bus waveguide, hereafter
referred to as ring (1) ($\omega^{(1)}=\omega_{0}+\delta\omega$,
$\omega^{(2)}=\omega^{(3)}=\omega_{0}$). For $\delta\omega=0$, the supermodes
are evenly spaced, as shown in the representative spectrum depicted in Fig.
1(d)). As $\delta\omega$ increases, the eigenfrequencies shift, resulting in
the avoided mode crossing curves observed in Fig. 2(d). These shifts in the
supermode frequencies can be exploited as a local dispersion control mechanism
within a triplet of supermodes [30]. We show in Fig. 2(e) that CMT predicts
that the supermode frequency spacing $J_{2}$ is controllable and can reach a
normal, zero, and anomalous local dispersion (sGVD). We experimentally
demonstrate this local dispersion control by employing integrated microheaters
atop each microresonator, allowing for controllable phase-shifting of the
bare-ring resonances through the thermo-optic effect. In Fig. 2(g), we
illustrate the tuning of the relative frequency spacing of the supermodes
achieved by varying the current supplied to the microheater on top of the ring
(1). This enables active supermode dispersion engineering in our devices,
allowing us to set up normal ($J_{2}<0$), zero ($J_{2}=0$), and anomalous
($J_{2}>0$) local dispersion regimes, as depicted in Fig. 2(h). Notably, this
local dispersion control (sGVD) was found to have a negligible impact on the
overall GVD of the supermodes.

To illustrate the control on sGVD achieved in our devices, we carried out
parametric wavelength conversion via stimulated FWM within a triplet of
coupled modes. In this process, the idler field is produced in the S supermode
by tuning the pump and signal fields to the C and AS supermodes, respectively,
as shown in the optimized optical spectrum in Fig. 2(f). The signal and pump
fields are generated using a high-extinction EOM (40 dB) bias to suppress the
optical carrier. The input (off-chip) peak power was set to 22 dBm ($3.5$
dB/facet insertion loss). For more experimental details, see Supplementary
Information Section S5. By varying the current applied to the microheater atop
ring (1), the sGVD is pretuned close to zero local dispersion regime, ensuring
the phase-matching condition for stimulated FWM. The color code inserted into
the experimental data points in Fig. 2(h) indicates the maximum power of the
idler field for each sGVD condition, indicated by the value of $J_{2}/2\pi$.
Parametric wavelength conversion is optimized close to the zero dispersion
condition $(J_{2}=0)$, at which the idler power was found to be $20.4$ dB
below the signal power level. Moreover, variations of $J_{2}/2\pi$ of around
300 MHz can cause an efficiency loss of nearly 3 dB, and beyond that range,
phase mismatch rapidly diminishes the process efficiency, and the idler
becomes obscured by FWM occurring in the fibers. Therefore, the sGVD control
offered by our photonic molecule design enables optimizing or suppressing the
phase-matching condition for the occurrence of FWM processes within a triplet
of supermodes.

##  3 DOPO generation

DOPO is generated in the photonic molecule by pumping the outer supermodes of
a triplet with a phase-coherent dual-pump scheme obtained by electro-optic
modulation of a single laser. The experimental setup used to study DOPO is
shown in Fig. 3(a). The high-extinction EOM 1(40 dB), is biased to ensure an
equal distribution of optical power between the carrier and the sidebands, and
the frequency of the RF drive source is adjusted to match the frequency
splitting between the AS and S supermodes ($\Omega/2\pi\approx 20$ GHz). The
optical carrier ($\omega_{\mathrm{b}}$) and one of the sidebands
($\omega_{\mathrm{r}}$), $\omega_{\mathrm{b}}>\omega_{\mathrm{r}}$, are used
to pump the supermodes $\omega_{\mathrm{AS}}$ and $\omega_{\mathrm{S}}$, as
illustrated in Fig. 3(b). The remaining sideband is not resonant with any
resonator mode and can be neglected in our experiment. To achieve high pump
powers, necessary to surpass the DOPO threshold, the pump fields are modulated
by another electro-optic intensity modulator (EOM 2), producing 6 ns pulses at
a repetition rate of 1.7 MHz. The pulse duration, approximately 17 times
larger than the cavity lifetime ($\tau=2Q_{L}/\omega\approx$ 0.37 ns), is long
enough to ensure a quasi-CW regime while significantly reducing thermo-optic-
induced phase shifts. The pulsed optical tones are amplified using an erbium-
doped fiber amplifier (EDFA) and coupled to the chip with the aid of lensed
fibers. From a measured input (off-chip) peak power of 8.2 W, we estimate in-
chip peak power of 1.45 W per pump field. Polarization controllers are used to
ensure the excitation of the fundamental transverse electric (TE) mode. The
output light is directed to an optical spectrum analyzer (OSA), and the DOPO
signal is detected by a photodetector after filtering out the pump wavelengths
(details of the filtering process can be seen at Section S6).

Figure 3: DOPO measurements: (a) Main components of the experimental setup for
the DOPO measurements. Tunable laser, EOM1: electro-optic modulator driven by
a sinusoidal RF-signal, EOM2: electro-optic modulator driven by a pulsed RF-
signal, EDFA: erbium-doped fiber amplifier, PC: polarization controllers, DUT:
device under test, OSA: optical spectrum analyzer, tunable narrow filter, and
PD: photodetector. (b) Schematics of the detuning $\Delta\omega_{\mathrm{AS}}$
($\Delta\omega_{\mathrm{S}}$) between the supermode frequency
$\omega_{\mathrm{AS}}$ ($\omega_{\mathrm{S}}$) and the pump frequency
$\omega_{\mathrm{b}}$ ($\omega_{\mathrm{r}}$). (c) Optical spectra collected
at the through port as a function of the normalized hot detuning
($\Delta\omega/\kappa_{\mathrm{C}}$) between the pump fields and the photonic
molecule resonances. The blue and red dashed lines indicate the pumped S and
AS supermodes. The colorbar was saturated to enable the visualization of the C
supermode oscillation, highlighted in white. (d) Optical spectrum at zero hot
detuning of an operating DOPO for zero (blue) and normal (pink) supermode
local dispersion regimes. (e,f) Transmission spectrum of the photonic molecule
around 1560 nm for the zero ($J_{2}=0$; blue) and normal ($J_{2}<0$; pink)
supermode local dispersion regimes.

DOPO occurs in the central supermode $\omega_{\mathrm{C}}$ as the dual-tone
pumps are simultaneously red-detuned across the AS and S supermodes
($\omega_{\mathrm{AS}}$ and $\omega_{\mathrm{S}}$). The optical spectra
collected at the through port as a function of hot pump detuning are shown in
Fig. 3(c). As the intracavity power increases, the supermode resonances are
red-shifted, mainly due to self- and cross-phase modulation (SPM and XPM,
respectively) contributions. The maximum DOPO power is observed when the hot
detunings of the pump fields approach zero, i.e.,
$\Delta\omega=\Delta\omega_{\mathrm{AS}}=\Delta\omega_{\mathrm{S}}=0$, with
$\Delta\omega_{\mathrm{AS}}\equiv\omega_{\mathrm{b}}-\omega_{\mathrm{AS}}$ and
$\Delta\omega_{\mathrm{S}}\equiv\omega_{\mathrm{r}}-\omega_{\mathrm{S}}$. Once
the pump fields are detuned beyond the bistable resonance edge, the intra-
cavity power abruptly drops, and the DOPO ceases to oscillate. The unloaded
supermodes’ spectral positions can be observed as the darkest regions in Fig.
3(c) due to the broadband spontaneous emission produced by the EDFA.

In Fig. 3(d), the DOPO spectral signature is displayed at zero hot detuning
$(\Delta\omega/\kappa_{\mathrm{C}}=0)$ for distinct sGVD regimes, fine-tuned
with the aid of the microheater on top of one of the outer rings. In
particular, in the normal dispersion regime ($J_{2}<0$, Fig. 3(f)), DOPO
exhibits a significant improvement compared to that in the zero dispersion
regime ($J_{2}=0$, Fig. 3(e)). In the normal sGVD regime, DOPO achieves power
levels 17 dB below the residual pump power, compatible with which it has been
achieved in single microrings [7, 11, 4]. The interplay between the SPM and
XPM effects in the supermodes plays a determinant role in the optimization of
DOPO. As the intracavity power reaches the nonlinear regime, the pumped S and
AS supermodes experience SPM and XPM, while the C supermode experiences XPM
contributions of the two pumps and negligible SPM. Since the effect of XPM is
stronger than that of SPM, the C supermode undergoes a larger redshift
compared to the S and AS supermodes, leading to an asymmetry in the frequency
spacing of the triplet of coupled modes. This is compensated for by presetting
the sGVD to the normal regime, as shown in Fig. 3(f). In this case, the local
normal dispersion imposed by the microheaters balances the nonlinear phase-
shifts, thus providing an optimum phase-matching condition for DOPO.

Figure 4: Time instability and OFC: Top (a,b) and bottom (c-f) figures present
the characterization of DOPO in time and frequency, respectively. (a)
Transmission of the filtered DOPO signal displaying amplitude fluctuation
across pulses. (b) Zoom-in on the highlighted DOPO pulse on (a) in blue,
overlaid with the input pump pulse (127 ns, orange) and the photonic lifetime
of the C supermode (0.37 ns, green). (c,e) Narrowband DOPO spectra centered at
the triplet of interest ($\lambda_{0}\sim$1560 nm) at normalized hot detunings
($\Delta\omega/\kappa_{\mathrm{C}}$) of 0.5 and 0.0, respectively. (d,f)
Broadband DOPO spectra at the same conditions of (c,e). Comb lines are
indicated in black.

We measured the DOPO’s temporal stability by filtering out the dual-tone pumps
from the pulsed output light. To remove the remaining from the DOPO signal, we
used a free-space Mach–Zehnder interferometer with the same FSR as the pump
tones frequency spacing, cascaded by a tunable narrowband filter to get a 23
dB rejection pump-to-DOPO ratio. For these measurements, we set the pulse
duration of the dual-tone pumps to 127 ns, which also provides a quasi-CW
regime within the microrings. As the dual-tone pump is tuned towards the
optical AS and S supermode resonances, oscilloscope traces reveal significant
variability in the amplitude reached for each pump pulse, as shown in Fig.
4(a). Inspecting a single pump pulse, as presented in Fig. 4(b), it is
noticeable that the DOPO pulse is shorter than the input pump pulse, with a
pulse duration of approximately 20 ns. This instability condition was verified
over the whole range of pump detuning and for a range of power levels above
the threshold. Our attempt to further explore the power-detuning parameter
space was hindered by the excitation of spurious comb lines located far away
from the pump frequencies, as illustrated in Fig. 4(d,f). When the pumped
power exceeds the threshold and DOPO initiates, two distant fundamental
supermodes around 1517 and 1605 nm ($\mu=\pm 12$) begin to oscillate at
similar power levels to DOPO itself (Fig. 4(c,d)). As the dual-tone pumps
further approach the AS and S resonances, more comb lines are progressively
populated. Such cascading FWM occurs, despite the normal GVD regime imposed on
the supermodes, because of the presence of avoided crossings between the
fundamental and higher-order TE modes of the bare rings. Experimental evidence
for the coupling to higher-order TE modes can be found in the fluctuations
seen in the dispersion curves of Fig. 2(a-c), and in the excited comb lines
beyond the triplet signature of our devices at 1517 and 1605 nm showed in
Section S7. The suppression of higher-order modes can be achieved by
exploiting design strategies that prevent their excitation by using nonlinear
adiabatic couplers [31] or that eliminate them by using selective loss
mechanisms [32] or coupling engineering [33]. In future work, our triple-state
photonic molecule design can be further improved by incorporating these
strategies.

The instability of Kerr-based DOPO in single microresonators has been
thoroughly investigated based on Lugiato-Lefever bifurcation dynamics[34, 4,
35]. It is known that in single microrings, the stability of the DOPO depends
on the suppression of non-degenerate FWM tones in the modes near the pump
wave[34, 35], which are reduced by spectral engineering of the photonic
molecules. Although a detailed experimental analysis of stability is needed to
identify the sources of the dynamics seen in Fig. 4(a,b), the temporal
duration of the DOPO signal indicates that thermal effects may also be
involved in our observations [36, 37]. We also do not exclude the possibility
of intrinsic instability due to bifurcations inherent in the Lugiato-Lefever
equations (LLEs) in coupled resonators [38]. However, our experimental
findings create an opportunity to investigate the intricate behavior of dual-
pumped coupled microresonators beyond the dimer regime [5].

##  4 Conclusions

We have realized DOPO in a coupled triple-mode photonic device. In our design,
the DOPO process happens entirely within the three coupled modes, without
depending on modes of different azimuthal mode numbers. This helps to isolate
DOPO from other nonlinear optical effects, such as non-degenerate FWM, which
takes place even in the normal GVD regime. Using integrated microheaters, we
could dynamically adjust the phases of the coupled modes, and thus optimize
phase-matching conditions for stimulated FWM and DOPO. We have demonstrated
that this active local dispersion control can be achieved independently of the
overall GVD regime of the microcavities, creating more opportunities for
spectral engineering. Furthermore, our design produces mode splittings with
frequency difference in the RF range, enabling dual-pump schemes by electro-
optic modulation of a single laser. This makes the experimental setup much
easier by eliminating the need for extra lasers and amplifiers and providing
native phase-locked pumps. The distant comb lines that compete with the DOPO
are due to mode interactions and can be reduced with future design
improvements that incorporate, for example, nonlinear adiabatic couplers [33].
Our triple-mode photonic device has great potential for DOPO applications,
such as genuine random number generation, coherent optical computing, and
generation and characterization of squeezed states of light.

Funding This work was supported by São Paulo Research Foundation (FAPESP)
through grants 18/15577-5, 18/15580-6, 18/21311-8, 18/25339-4, 20/04686-8,
21/10334-0, 22/06267-8, 23/09412-1, and Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior - Brasil (CAPES) (Finance Code 001)

Acknowledgments The authors would like to acknowledge Prof. Felippe Barbosa
for insightful discussions on the DOPO behavior.

Disclosures The authors declare no conflicts of interest.

Data Availability Statement Data underlying the results presented in this
paper are not publicly available at this time but may be obtained from the
authors upon reasonable request.

Supplemental document The Supplementary Material provides additional data and
characterization of the DOPO presented in the main text.

## References

  * [1] A. Marandi, N. C. Leindecker, K. L. Vodopyanov, and R. L. Byer, “All-optical quantum random bit generation from intrinsically binary phase of parametric oscillators,” Optics express 20, 19322–19330 (2012).
  * [2] Y. Okawachi, M. Yu, K. Luke, _et al._ , “Quantum random number generator using a microresonator-based kerr oscillator,” Optics letters 41, 4194–4197 (2016).
  * [3] Y. Okawachi, B. Y. Kim, Y. Zhao, _et al._ , “Dynamic control of photon lifetime for quantum random number generation,” Optica 8, 1458–1461 (2021).
  * [4] Y. Zhao, Y. Okawachi, J. K. Jang, _et al._ , “Near-degenerate quadrature-squeezed vacuum generation on a silicon-nitride chip,” Physical Review Letters 124, 193601 (2020).
  * [5] Y. Zhang, M. Menotti, K. Tan, _et al._ , “Squeezed light from a nanophotonic molecule,” Nature communications 12, 1–6 (2021).
  * [6] A. Marandi, Z. Wang, K. Takata, _et al._ , “Network of time-multiplexed optical parametric oscillators as a coherent ising machine,” Nature Photonics 8, 937–942 (2014).
  * [7] Y. Okawachi, M. Yu, J. K. Jang, _et al._ , “Demonstration of chip-based coupled degenerate optical parametric oscillators for realizing a nanophotonic spin-glass,” Nature communications 11, 1–7 (2020).
  * [8] T. Inagaki, K. Inaba, R. Hamerly, _et al._ , “Large-scale ising spin network based on degenerate optical parametric oscillators,” Nature Photonics 10, 415–419 (2016).
  * [9] C. J. McKinstrie and S. Radic, “Phase-sensitive amplification in a fiber,” Optics Express 12, 4973–4979 (2004).
  * [10] C. Nabors, S. Yang, T. Day, and R. Byer, “Coherence properties of a doubly resonant monolithic optical parametric oscillator,” JOSA B 7, 815–820 (1990).
  * [11] Y. Okawachi, M. Yu, K. Luke, _et al._ , “Dual-pumped degenerate kerr oscillator in a silicon nitride microresonator,” Optics letters 40, 5267–5270 (2015).
  * [12] Y. A. Espinel, F. G. Santos, G. O. Luiz, _et al._ , “Brillouin optomechanics in coupled silicon microcavities,” Scientific reports 7, 43423 (2017).
  * [13] B. E. Little, S. T. Chu, H. A. Haus, _et al._ , “Microring resonator channel dropping filters,” Journal of lightwave technology 15, 998–1005 (1997).
  * [14] M. A. Popović, C. Manolatou, and M. R. Watts, “Coupling-induced resonance frequency shifts in coupled dielectric multi-cavity filters,” Optics express 14, 1208–1222 (2006).
  * [15] X. Zeng and M. A. Popović, “Design of triply-resonant microphotonic parametric oscillators based on Kerr nonlinearity,” Optics Express 22, 15837 (2014).
  * [16] M. C. Souza, G. F. Rezende, L. A. Barea, _et al._ , “Spectral engineering with coupled microcavities: active control of resonant mode-splitting,” Optics letters 40, 3332–3335 (2015).
  * [17] A. Tusnin, A. Tikan, K. Komagata, and T. J. Kippenberg, “Nonlinear dynamics and kerr frequency comb formation in lattices of coupled microresonators,” Communications Physics 6, 317 (2023).
  * [18] A. Pal, A. Ghosh, S. Zhang, _et al._ , “Linear and nonlinear coupling of twin-resonators with kerr nonlinearity,” arXiv preprint arXiv:2404.05646 (2024).
  * [19] M. Singh, X. Zhang, B. Zhang, _et al._ , “Efficient wideband tunable radio frequency–optical conversion via triply resonant photonic molecules,” Optics Letters 49, 3930–3933 (2024).
  * [20] J. Wang, K. Liu, A. Isichenko, _et al._ , “Integrated programmable strongly coupled three-ring resonator photonic molecule with ultralow-power piezoelectric control,” Optics Letters 48, 2373–2376 (2023).
  * [21] X. Xue, Y. Xuan, P.-H. Wang, _et al._ , “Normal-dispersion microcombs enabled by controllable mode interactions,” Laser & Photonics Reviews 9, L23–L28 (2015).
  * [22] S. Fujii, Y. Okabe, R. Suzuki, _et al._ , “Analysis of Mode Coupling Assisted Kerr Comb Generation in Normal Dispersion System,” IEEE Photonics Journal 10, 1–11 (2018).
  * [23] B. Y. Kim, Y. Okawachi, J. K. Jang, _et al._ , “Turn-key, high-efficiency Kerr comb source,” Optics Letters 44, 4475–4478 (2019). Publisher: Optica Publishing Group.
  * [24] C. M. Gentry, X. Zeng, and M. A. Popović, “Tunable coupled-mode dispersion compensation and its application to on-chip resonant four-wave mixing,” Optics letters 39, 5689–5692 (2014).
  * [25] H. A. Haus and W. Huang, “Coupled-mode theory,” Proceedings of the IEEE 79, 1505–1518 (1991).
  * [26] L. A. Barea, F. Vallini, P. F. Jarschel, and N. C. Frateschi, “Silicon technology compatible photonic molecules for compact optical signal processing,” Applied Physics Letters 103, 201102 (2013).
  * [27] A. Tikan, J. Riemensberger, K. Komagata, _et al._ , “Emergent nonlinear phenomena in a driven dissipative photonic dimer,” Nature Physics 17, 604–610 (2021).
  * [28] X. Zeng and M. A. Popović, “Design of triply-resonant microphotonic parametric oscillators based on kerr nonlinearity,” Optics Express 22, 15837–15867 (2014).
  * [29] J. Liu, V. Brasch, M. H. P. Pfeiffer, _et al._ , “Frequency-comb-assisted broadband precision spectroscopy with cascaded diode lasers,” Optics Letters 41, 3134 (2016).
  * [30] T. Herr, V. Brasch, J. Jost, _et al._ , “Mode spectrum and temporal soliton formation in optical microresonators,” Physical review letters 113, 123901 (2014).
  * [31] A. Kordts, M. H. Pfeiffer, H. Guo, _et al._ , “Higher order mode suppression in high-q anomalous dispersion sin microresonators for temporal dissipative kerr soliton formation,” Optics letters 41, 452–455 (2016).
  * [32] X. Ji, J. K. Jang, U. D. Dave, _et al._ , “Exploiting ultralow loss multimode waveguides for broadband frequency combs,” Laser & Photonics Reviews 15, 2000353 (2021).
  * [33] F. A. Sabattoli, H. El Dirani, L. Youssef, _et al._ , “Suppression of parasitic nonlinear processes in spontaneous four-wave mixing with linearly uncoupled resonators,” Phys. Rev. Lett. 127, 033901 (2021).
  * [34] A. Roy, R. Haldar, and S. K. Varshney, “Analytical model of dual-pumped kerr micro-resonators for variable free spectral range and composite frequency comb generation,” Journal of Lightwave Technology 36, 2422–2429 (2018).
  * [35] T. Hansson and S. Wabnitz, “Bichromatically pumped microresonator frequency combs,” Phys. Rev. A 90, 013811 (2014).
  * [36] M. Gao, Q.-F. Yang, Q.-X. Ji, _et al._ , “Probing material absorption and optical nonlinearity of integrated photonic materials,” Nature communications 13, 3323 (2022).
  * [37] M. He and K. Jamshidi, “Conditions for dual-pumped optical parametric oscillation in kerr microresonators,” Phys. Rev. Appl. 20, 054036 (2023).
  * [38] A. Ghosh, A. Pal, L. Hill, _et al._ , “Controlled light distribution with coupled microresonator chains via kerr symmetry breaking,” arXiv preprint arXiv:2402.10673 (2024).
  * [39] T. Herr, K. Hartinger, J. Riemensberger, _et al._ , “Universal formation dynamics and noise of kerr-frequency combs in microresonators,” Nature photonics 6, 480–487 (2012).
  * [40] Y. K. Chembo and N. Yu, “Modal expansion approach to optical-frequency-comb generation with monolithic whispering-gallery-mode resonators,” Physical Review A 82, 033801 (2010).
  * [41] X. Xue, Y. Xuan, P.-H. Wang, _et al._ , “Normal-dispersion microcombs enabled by controllable mode interactions,” Laser & Photonics Reviews 9, L23–L28 (2015).
  * [42] L. Fujii, M. Inga, J. H. Soares, _et al._ , “Dispersion tailoring in wedge microcavities for kerr comb generation,” Optics letters 45, 3232–3235 (2020).

## Supplementary Information

##  S1 DOPO in photonic molecules

The slow time evolution of the supermodes’ amplitudes involved in DOPO can be
described by the coupled Lugiato-Lefever equations (LLE) [39, 40, 28], written
in the supermode basis:

$\displaystyle\frac{\partial a_{1}}{\partial\tau}$
$\displaystyle=-\left(\frac{\kappa_{1}}{\kappa_{2}}+i\frac{\Delta_{1}}{\kappa_{2}}\right)a_{1}+i\left(f^{111}_{1}|a_{1}|^{2}+2f^{221}_{1}|a_{2}|^{2}+2f^{331}_{1}|a_{3}|^{2}\right)a_{1}+if^{232}_{1}a_{2}^{2}a_{3}^{\ast}+s_{\text{b}},$
(S1a) $\displaystyle\frac{\partial a_{2}}{\partial\tau}$
$\displaystyle=-\left(1+i\frac{\Delta_{2}}{\kappa_{2}}\right)a_{2}+i\left(2f^{112}_{2}|a_{1}|^{2}+f^{222}_{2}|a_{2}|^{2}+2f^{332}_{2}|a_{3}|^{2}\right)a_{2}+i2f^{123}_{2}a_{1}a_{2}^{\ast}a_{3},$
(S1b) $\displaystyle\frac{\partial a_{3}}{\partial\tau}$
$\displaystyle=-\left(\frac{\kappa_{3}}{\kappa_{2}}+i\frac{\Delta_{3}}{\kappa_{2}}\right)a_{3}+i\left(2f^{113}_{3}|a_{1}|^{2}+2f^{223}_{3}|a_{2}|^{2}+f^{333}_{3}|a_{3}|^{2}\right)a_{3}+if^{212}_{3}a_{2}^{2}a_{1}^{\ast}+s_{\text{r}}.$
(S1c)

The equations are normalized by $\tau=\kappa_{2}t/2$, where $\kappa_{j}$ is
the total loss of the $j^{\text{th}}$supermode. Here, the index $j=1,2,$ and
$3$ is used to denote the AS, C, or S supermode, respectively. The normalized
amplitude of the $j^{\text{th}}$ supermode is given by
$a_{j}=\sqrt{2g_{0}/\kappa_{2}}A_{j}e^{-i\Delta_{j}t}$, where $g_{0}$ is the
Kerr coefficient, $A_{j}$ is the unnormalized amplitude of the $j^{\text{th}}$
supermode, $\omega_{j}$ is its frequency, and $\Delta_{j}$ is its detuning.
Furthermore, $s_{b}$ and $s_{r}$ represent normalized blue and red pumps,
respectively. The third, fourth and fifth terms in the RHS of Eq. S1 denote
the effects of SPM and XPM, while the sixth term accounts for the FWM
interactions. Nonlinear interactions are weighted by the spatial overlap
between the involved modes, calculated via the overlap integrals [40]:

$f^{\alpha\beta\gamma}_{\mu}=\frac{\omega_{\mu}}{\omega_{0}}\frac{\int\epsilon^{2}\,(\vec{\mathcal{E}}^{\ast}_{\mu}\cdot\vec{\mathcal{E}}_{\alpha})\,(\vec{\mathcal{E}}^{\ast}_{\beta}\cdot\vec{\mathcal{E}}_{\gamma})\,dV}{\int\epsilon^{2}\,||\vec{\mathcal{E}}_{2}||^{4}\,dV}.$
(S2)

Here $\vec{\mathcal{E}}_{j}$ denotes the spatial mode profile of the
$j^{\text{th}}$ supermode. Notice that, due to the distinct mode distributions
of the supermodes shown in Fig. 1(e), the overlap integrals
$f_{\mu}^{\alpha\beta\gamma}$ can not be approximated by 1. To estimate the
nonlinear overlaps, we perform a change of basis to write the supermodes
$\vec{\mathcal{E}}_{j}(\vec{r})$ as a linear combination of the bare modes
$\vec{e}_{j}(\vec{r})$. As an example, the integrand in the numerator of the
FWM coefficient $f^{123}_{2}$ is calculated as follows:

$\begin{split}\left(\vec{\mathcal{E}}_{2}^{\ast}\cdot\vec{\mathcal{E}}_{1}\right)\left(\vec{\mathcal{E}}_{2}^{\ast}\cdot\vec{\mathcal{E}}_{3}\right)&=\frac{1}{8}\left(\vec{e}_{3}-\vec{e}_{1}\right)^{\ast}\cdot\left(\vec{e}_{1}-\sqrt{2}\vec{e}_{2}+\vec{e}_{3}\right)\left(\vec{e}_{3}-\vec{e}_{1}\right)^{\ast}\cdot\left(\vec{e}_{1}+\sqrt{2}\vec{e}_{2}+\vec{e}_{3}\right)\\\
&\approx\frac{|\vec{e}_{1}|^{4}+|\vec{e}_{3}|^{4}-2\left(\vec{e}_{1}^{\ast}\cdot\vec{e}_{2}\right)^{2}-2\left(\vec{e}_{3}^{\ast}\cdot\vec{e}_{2}\right)^{2}}{8},\end{split}$
(S3)

where we have discarded the field products involving rings 1 and 3 because
there is no spatial overlap between the outer rings. We chose to normalize all
$f_{\mu}^{\alpha\beta\gamma}$ by the self-overlap of the central mode
($f_{2}^{222}$):

$\begin{split}|\vec{\mathcal{E}}_{\text{2}}|^{4}&=\frac{1}{4}\left[\left(\vec{e}_{3}-\vec{e}_{1}\right)^{\ast}\cdot\left(\vec{e}_{3}-\vec{e}_{1}\right)\right]^{2}\\\
&\approx\frac{1}{4}\left(|\vec{e}_{1}|^{2}+|\vec{e}_{3}|^{2}\right)^{2}\\\
&\approx\frac{|\vec{e}_{1}|^{4}+|\vec{e}_{3}|^{4}}{4}.\end{split}$ (S4)

Since

$\int\epsilon^{2}\left(|\vec{e}_{1}|^{4}+|\vec{e}_{3}|^{4}\right)dV\gg
2\int\epsilon^{2}\left[2\left(\vec{e}_{1}^{\ast}\cdot\vec{e}_{2}\right)^{2}+2\left(\vec{e}_{3}^{\ast}\cdot\vec{e}_{2}\right)^{2}\right]dV,$
(S5)

we find $f^{123}_{2}\approx 1/2$. Following the same reasoning, we obtain the
remaining nonlinear coupling coefficients:

$\begin{pmatrix}f^{111}_{1}&f^{221}_{1}&f^{331}_{1}&f^{232}_{1}\\\
f^{112}_{2}&f^{222}_{2}&f^{332}_{2}&f^{123}_{2}\\\
f^{113}_{3}&f^{223}_{3}&f^{333}_{3}&f^{212}_{3}\end{pmatrix}=\begin{pmatrix}3/4&1/2&3/4&1/2\\\
1/2&1&1/2&1/2\\\ 3/4&1/2&3/4&1/2\end{pmatrix}.$ (S6)

Notice that some of the SPM and XPM terms become more prominent than the FWM
terms, and XPM can no longer be described as twice as strong as SPM. At the
same time, according to Eq. S4, the effective mode volume of the reference
mode of the molecule is twice that of the single ring, which means that the
threshold power for DOPO will be higher in the molecule.

##  S2 Inter-ring coupling coefficient of a photonic molecule

The calculation presented in this section follows a treatment of coupled mode
theory (CMT) described in ref [25] for weakly coupled and lossless
microcavities. In this treatment, the electric field in the microcavity is
described as a linear superposition of its mode fields:

$\vec{E}=\sum_{i}a_{i}\,\vec{e}_{i},$ (S7)

in which $a_{i}$ and $\vec{e}_{i}$ are the amplitude and field profile of the
$i$-th mode. Using this trial solution in Maxwell’s equations, the coupled
mode equations for two-coupled single-mode microcavities become:

$\frac{d\vec{a}}{dt}=i\,\mathbb{H}\,\vec{a},$ (S8)

where $\vec{a}=(a_{1}\>a_{2})^{T}$ is a vector containing the mode amplitudes
of each microcavity. The matrix $\mathbb{H}$ accounts for the eigenfrequencies
of the individual modes (diagonal elements) and the coupling between them
(off-diagonal elements). Its elements are calculated as follows:

$H_{ij}=-\frac{\omega_{0}}{2}\int\vec{e_{i}}^{*}\cdot\delta\epsilon_{i}\>\vec{e}_{j}\,dV,$
(S9)

where $\vec{e}_{i,j}$ and $\omega_{0}$ are, respectively, the mode field
profile and angular frequency of each microcavity and
$\delta\epsilon_{i}=\epsilon-\epsilon_{j}$ is the difference between the
spatially-varying dielectric constant and that experienced by the $j$-th
microcavity in the absence of the $i$-th microcavity. Fig. S1(a) illustrates
how the quantity $\delta\epsilon_{i}$ is calculated. The mode field profiles
($\vec{e}_{i,j}$) are normalized so that the square of their amplitudes
($|a_{i,j}|^{2}$) is equal to the energy carried by the modes:

$U_{i,j}=|{\vec{a}}_{i,j}|^{2}\,\frac{1}{2}\,\int\epsilon|{\vec{e}}_{i,j}|^{2}dV\equiv|{\vec{a}}_{i,j}|^{2},$
(S10)

which implies:

$\vec{e}_{i,j}=\frac{\tilde{\vec{e}}_{i,j}}{\sqrt{\int\frac{1}{2}\epsilon|\tilde{\vec{e}}_{i,j}|^{2}dV}}.$
(S11)

with $\tilde{\vec{e}}_{i,j}$ being the E-field in V/m units.

The Eqs. S8, S10, S11 and S9 for a two-coupled microring system were
implemented in Comsol Multiphysics. The mode field profiles ($\vec{e}_{i,j}$)
of each microring were obtained through a mode analysis study performed in a
2D axisymmetric component, where the microring cross section is defined. The
field profile of the TE fundamental mode of a silicon nitride (SiN) microring
is shown in the inset of Fig. S1(b). Once the desired mode solutions were
found, they were projected into a 3D component using Comsol extrusion
operators. As shown in Fig. S1(b), the 3D component encloses only the coupling
domain, where the field overlap between the modes of each microring is non-
zero.

Figure S1: Inter-ring coupling calculation: (a) Dielectric constant
distributions used in the calculation of the coupling matrix $\mathbb{H}$,
represented in the transverse plane of the smallest gap between two
microrings. $\epsilon$ accounts for the overall dielectric constant,
$\epsilon_{1}$ accounts for the overall dielectric constant in the absence of
the other microring, and $\delta\epsilon_{2}$ is the difference between
$\epsilon$ and $\epsilon_{1}$. $\varepsilon_{0}$, $\varepsilon_{1}$ and
$\varepsilon_{2}$ are the dielectric constants of the surrounding medium,
microring 1 and microring 2, respectively. (b) 3D component enclosing the
coupling domain between two microrings. The fundamental TE mode of each
microring, found through a mode analysis study in Comsol, is shown in the
inset. They were projected into the 3D component using Comsol extrusion
operators to calculate the field overlap integrals of the matrix $\mathbb{H}$.
The microrings used in this simulation are SiN microrings buried in a SiO2
substrate with 50 $\mu$m of radius and 1 $\mu$m x 0.8 $\mu$m of width x
height.

##  S3 Device characterization

Our photonic molecule design comprises three coupled microring resonators,
illustrated in the micrograph in Fig. S2(a). The central ring is coupled to
both outer rings, each coupled solely to the central ring and a single (input
or end) waveguide. Light is launched into the chip using lensed fibers, while
inverse tapers optimize the coupling for TE-polarized light. The loaded
Q-factor of the central supermode was found to be
$Q_{L}=\omega_{0}/\Delta\omega=2.2\times 10^{5}$ by fitting the narrow
through-port transmission around the central supermode, shown in Fig. S2(b),
using the equation:

$T(\omega)=\frac{4(\omega-\omega_{0})^{2}+(\kappa_{i}-\kappa_{e})^{2}}{4(\omega-\omega_{0})^{2}+(\kappa_{i}+\kappa_{e})^{2}}.$
(S12)

After determining the intrinsic and extrinsic losses, $\kappa_{i}$ and
$\kappa_{e}$, respectively, the resonance linewidth $\Delta\omega$ was
calculated through:

$\Delta\omega=\kappa_{i}+\kappa_{e}.$ (S13)

Figure S2: Photonic molecule characterization: (a) Micrograph of the device,
showing the photonic molecule, the bus waveguide, and the integrated
microheaters. (b) Transmission spectrum around the triplet of interest, used
to calculate the quality factor. Individual supermodes are fitted with
Lorentzian curves. The fitting of the central supermode yields a loaded
quality factor of $Q_{L}=2.2\times 10^{5}$. (c) Effective index
$n_{\text{eff}}$ as a function of the wavelength for the TE (orange) and TM
(green) optical modes supported in the microring resonators employed in our
photonic molecules. The fundamental (two upper curves) and first high-order
modes (two lower curves) are displayed, alongside the ring waveguide (SiN) and
cladding (SiO2) refractive indexes. The insert shows the electric field
profile of each mode, simulated using the finite element method. (d)
Transmission at the device’s through port as the input light is polarized to
excite mainly the TE (orange) and TM (green) modes.

Moreover, a 2D axisymmetric simulation is employed to calculate the optical
modes supported by the individual microring resonators. Fig. S2(c) shows the
effective index and mode profile of the eigenmode solutions, revealing the
support of the fundamental and the first higher order TE and TM optical modes.
In the photonic molecule, these modes can couple across the rings to produce
higher-order or TM-polarized supermodes whose presence can potentially disturb
the dispersion locally by avoided mode crossing [41]. This process can occur
even if the coupling between the waveguide and the photonic molecule is weaker
for these parasitic supermodes. In Fig. S2(d), we present the through-port
transmission for both TE and TM polarized light, where higher-order modes are
not visible. However, their presence is evident from the fluctuations in the
dispersion curves of the supermodes shown in Fig. 2(a-c).

##  S4 General description of supermode dispersion

To study the dispersion of an optical resonator, the standard procedure
involves expanding the mode frequencies $\omega$ around a reference mode
($m_{0}$) with frequency $\omega_{0}$, using the relative azimuthal mode
$\mu=m-m_{0}$ [42]. For photonic molecules, the dispersion can be individually
calculated for each supermode by fixing the supermode index $j$:

$\omega_{j,\mu}=\omega_{j,0}+D_{1,j}\mu+\frac{1}{2}D_{2,j}\mu^{2}+\dots$ (S14)

Here, the coefficients $D_{1,j}=\frac{\partial\omega_{j,\mu}}{\partial\mu}$
and $D_{2,j}=\frac{\partial^{2}\omega_{j,\mu}}{\partial\mu^{2}}$ denote the
average spectral spacing (FSR) and its variation in the $\mu$-space (GVD
parameter). Alternatively, it might be interesting to define the supermode
dispersion within a given azimuthal mode number ($\mu$ fixed) [17]. In a
photonic molecule composed of an array of $N$ identical (degenerate)
resonators, the frequency splitting induced by the linear evanescent coupling
yields N-tuples that repeat for each azimuthal number $m$. Expanding the
supermode frequency in relation to the supermode index results in a new
dispersion relation:

$\omega_{j,\mu}=\omega_{j_{0},\mu}+J_{1,j}(j-j_{0})+\frac{1}{2}J_{2,j}(j-j_{0})^{2}+\dots$
(S15)

In this case, $J_{1,j}=\frac{\partial\omega_{j,\mu}}{\partial j}$ represents
the average frequency splitting,
$J_{2,j}=\frac{\partial^{2}\omega_{j,\mu}}{\partial j^{2}}$ denotes its
variation along the $N$-tuple supermode, and so on. These define the supermode
FSR (sFSR, $J_{1,j}$) and supermode GVD (sGVD, $J_{2,j}$), similar to the
standard dispersion definitions in the $\mu$-space. For $N=3$ ($j=\\{\text{S,
C, and AS}\\}$, the design explored in this work), up to the second-order
supermode dispersion is defined, which is sufficient to characterize normal
($J_{2,j}<0$) and anomalous ($J_{2,j}>0$) supermode regimes, locally within a
triplet indexed by the mode number $\mu$.

##  S5 Four-wave mixing experimental details

Stimulated FWM is performed using the same setup illustrated in Fig. 3(a). We
demonstrate parametric frequency conversion by adjusting the pump and signal
fields across the C and AS supermodes, producing an idler field at the S
supermode, as shown in Fig. 2(f). Dual-pump tones are generated using a high-
extinction EOM ($\Omega\sim 5$ GHz, half the supermode frequency splittings)
biased to suppress the carrier, while an EDFA delivers 22 dBm off-chip peak
power. Moreover, a second EOM, driven by a pulsed source, is integrated to the
setup to pulse the optical pump fields, mitigating the thermo-optic effect.
Fig. S3(c) displays a typical spectrum before the chip, showing FWM occurring
in the fibers. In this scenario, the higher-wavelength field generated in the
fiber coincides with the idler field of the photonic molecule FWM, although
this does not occur in the DOPO measurements.

In Fig. S3(a, b), the optical spectra are presented for CW (a) and pulsed
(quasi-CW, b) regimes as the dual-tone pump is simultaneously finely detuned
using the laser’s built-in piezo control. Since the pulse width of 50 ns (135
times longer than the photonic lifetime) is much faster than the thermo-optic
response time [36], red phase-shift and bistability are drastically reduced in
the pulsed regime. In Fig. S3(d), we characterize the conversion efficiency
relative to the input power. Initially, the conversion efficiency increases
linearly with the pump power, but it saturates around 23 dBm due to induced
phase-mismatch. Since the supermodes present different spatial distributions,
they experience nonuniform phase-shift contributions from thermo-optic and
Kerr effects. Additionally, the idler phase shift undergoes significant XPM
contributions from both pumps, whereas the pumps experience SPM and XPM. Since
the effect of XPM is stronger than SPM, the idler field exerts a larger red-
shift. Consequently, the supermodes achieve the bistability edge at distinct
detunings, terminating the FWM process early, as shown in Fig. S3(e). This
effect becomes more pronounced as the input power increases, affecting the
efficiency trend of the process.

Figure S3: FWM: (a,b) FWM optical spectra collected at the drop port as a
function of the hot detuning for cw (a) and pulsed (b) pumps. The dashed lines
indicate the supermode frequencies, red-shifted by thermo-optic and Kerr phase
shifts. (c) Typical optical spectrum before the chip, showing the surpassed
carrier and FWM occurring in the fibers. (d) Signal-to-idler conversion
efficiency as a function of the CW input power. The efficiency follows a
linear trend until it saturates. (e) Power of each supermode as a function of
the detuning in the CW regimes. The dashed lines highlight the detuning at
which the power is maximum for each curve.

##  S6 DOPO signal filtering

The spectral proximity of the supermodes involved in the DOPO process
(separation of $\sim 10~{}\text{GHz}$), while advantageous for the pumping
scheme, introduces experimental challenges. This, associated with the high
contrast between the two pumps and the DOPO signal, with a best-case scenario
of 17 dB, complicates the filtering of the desired signal. This challenge is
particularly crucial when analyzing the temporal stability of DOPO, as shown
in Fig. 4(a,b). Even minor contributions from the pumps to the detected signal
can overshadow the signal of interest when using a photodiode. Therefore, to
enhance the relative intensity of the DOPO signal compared to the pumps, we
employed a series of optical filters: first, we used a KYLIA WT-MINT delay
line interferometer with variable FSR as an optical filter by setting the
transmission minima to coincide with the pump positions. This approach was
used because the interferometer is capable of withstand the high power from
the chip, reducing the pump intensity to levels comparable to the DOPO signal.
Subsequently, a Yenista XTM-50 tunable optical filter with adjustable
bandwidth was used to further attenuate the pump fields and any remaining ASE
from the EDFA. As a result, the DOPO signal was enhanced to 23 dB above the
filtered pumps (Fig. S4), making it suitable for measurement with the
photodiode.

Figure S4: DOPO filtering: Narrow (a) and broadband (b) DOPO spectra
displaying the narrow filtering achieved to isolate the DOPO signal resonance
frequency. Not filtered (light blue) and filtered spectra (dark blue) are
shown.

##  S7 Optical frequency comb

Cascading FWM occurs despite the normal GVD regime imposed on the supermodes
of our photonic molecule. Here, we show that this phenomenon is attributed to
the presence of avoided crossings between the fundamental and higher-order TE
modes of the bare rings.

As the dual-tone pumps approach the antisymmetric (AS) and symmetric (S)
resonances, comb lines are progressively populated, accompanied by the DOPO
signal. Fig. S5(a) shows the excitation of the first comb lines around 1517
and 1605 nm ($\mu=\pm 12$) as soon as DOPO begins to oscillate. The first sign
of DOPO is evidenced in Fig. S5(b) by an asymmetry in the cavity resonance,
seen as a result of the broadband spontaneous emission of the erbium-doped
fiber amplifier. The intensity of these initial comb lines is orders of
magnitude smaller than the intensity of the pumps, indicating that the
oscillation at the $\omega_{c}$ frequency is due to pump interaction rather
than secondary cascaded nonlinear processes. As the detuning between the pumps
and AS and S resonances is further reduced, additional comb lines are
populated and an increase in the DOPO amplitude is observed (Fig. S5(c, d)).
This dynamic process is possibly a source of the instability measured in the
DOPO amplitude.

Fig. S5(e, f) provide a zoomed view of the spectrum around the first comb
lines. In this scenario, modes beyond the triplet signature of our devices are
excited, providing experimental evidence of undesired coupling to higher-order
TE modes. This coupling leads to avoided crossings that enable phase-matching
for the initial FWM effect, causing energy to be transferred to other triplets
rather than the one being pumped.

Figure S5: OFC (a) Broadband and (b) narrowband spectra collected as the pumps
approach the resonances. The position of first lines spectrally distant from
the pumps are highlighted in (a). The first sign of DOPO is seen in (b) as an
asymmetry in the cavity resonance. (c) Broadband and (d) narrowband spectra
collected as the detuning is reduced. (e, f): Zoomed view of the spectrum
around the first comb lines around 1517 nm (e) and 1605 nm (f), showing
excited modes beyond the triplet signature expected for the fundamental TE
supermodes.