Multimode squeezing in soliton crystal microcombs

Soliton microcombs are self-organized pulses of light sustained in driven Kerr microresonators, intensively studied for applications in integrated photonic technologies and for their rich nonlinear dynamics. In this work, we theoretically study the collective dynamics of the quantum ﬂuctuations of soliton microcombs. We ﬁnd that the mean ﬁeld of a dissipative Kerr soliton crystal is accompanied by pulses of squeezed multimode vacuum and derives its operational stability from the strong detuning of the below-threshold parametric process. We present a photonic architecture that enables independent control of the above-and below-threshold states and achieves a high degree of squeezing ( > 15 dB) in the output waveguide with realistic losses. Our work elucidates the quantum dynamics of formation and annihilation in dissipative Kerr soliton systems, and establishes a pathway for the realization of a practical integrated source of multimode squeezedlight.

The optical mean field of the Kerr frequency comb (i.e., "classical" comb) has been modeled with great success by the Lugiato-Lefever equation (LLE) [10], a nonlinear Schrödinger equation that includes dissipation, drive, and detuning. The most commonly studied configuration is that where a single coherent pump laser supplies parametric gain to populate the comb lines through stimulated four-wave mixing [11]. Spontaneous light generation is absent from the model; to seed threshold processes such as optical parametric oscillation (OPO), random noise must be added into every optical mode. Thus, the LLE can reveal neither the properties of the quantum state of the comb nor the coherent dynamics of the threshold processes that drive the formation of the comb itself.
The quantum state of soliton microcombs has received little attention, and studies to date have approached it via an extension of pairwise mode analysis [12,13] which is used to describe signal-idler quantum frequency combs [14][15][16]. However, the signal-idler basis cannot be expected to reveal the structure of the collective quantum fluctuations in a soliton microcomb due to the extended modal coupling (Fig. 1), suggesting the need for a multimode analysis. An example of a well-studied multimode system is the synchronously pumped optical parametric oscillator (SPOPO), where an external mode-locked classical comb source drives a quantum comb via a second-order (χ (2) ) nonlinear process: the SPOPO is naturally described in a basis of supermodes (i.e., superpositions of frequency modes) [17,18], and this basis reveals multimode quadrature squeezing [19,20]. Three key features distinguish the soliton microcomb from commonly studied squeezed light sources. First, Kerr cavity soliton systems typically feature (and rely on) significant modal dispersion, which contributes strongly detuned parametric processes. Second, the third-order (χ (3) ) nonlinearity in soliton microcombs introduces four-wave mixing nonlinear terms not present in χ (2) systems, most notably Bragg scattering (frequency translation of a photon). In degenerate and signal-idler squeezing schemes, these system properties are often considered as parasitic processes that degrade squeezing [21,22]. Finally, the quantum fluctuations in a soliton microcomb are driven not by an external source [23], but by a coherent comb that is itself generated in situ via the same Kerr nonlinearity, intimately linking the mean-field and the below-threshold states.
In this work, we apply a multimode quadrature squeezing analysis to the soliton microcomb and its formation. We show that multimode analysis is necessary to understand Kerr microcomb threshold processes beyond the single-pump regime, and that such analysis can predict the properties of the post-threshold mean field, such as spatiotemporal oscillations. We find that the quantum state of the soliton crystal microcomb is highly squeezed across the entire range of its existence, and that the passage from a soliton crystal to a single-soliton state is a coherent threshold process associated with asymptotic growth in squeezing. Finally, we describe how the soliton crystal can be engineered as a practical source of multimode quadrature-squeezed light. Top: a continuous-wave pump laser is evanescently coupled into the microring through a bus waveguide and used to generate the soliton crystal. In the temporal domain, the below-threshold state (red) co-propagates with the coherent soliton pulses (blue). Bottom: in the frequency mode basis, the above-and below-threshold modes form two subsets, allowing the quantum fluctuations to be studied in isolation.

SUPERMODE DECOMPOSITION
The quantum fluctuations of the soliton microcomb are described by a many-body quartic Hamiltonian, where all resonator modes are coupled by a four-photon interaction through the χ (3) nonlinearity of the medium. In a soliton crystal microcomb [8,24], the resonator modes can be partitioned into two sets: (i) abovethreshold modes populated via stimulated four-wave mixing and comprising the classical dissipative Kerr soliton state, and (ii) the below-threshold modes populated only via spontaneous four-wave mixing [13]. Due to this robust partition guaranteed by energy and momentum conservation, the soliton crystal forms a unique system for studying in isolation the quantum fluctuations of dissipative Kerr solitons. To model the quantum state, the simplest approach is to assume that the above-threshold modes are in a classical coherent state [25] (complex amplitudes A m ) and that this coherent state drives parametric processes in the below-threshold modesâ j . This assumption constitutes the linearization of the quartic Hamiltonian into a computationally tractable quadratic Hamiltonian: (1) where δ FWM = δ ( j +k−m−n) is the four-wave mixing mode number matching condition, g 0 is the nonlinear coupling coefficient, and h.c. is the Hermitian conjugate. This approximation has been used to predict the second-order photon correlations that exist across the below-threshold modes of a soliton crystal state, measured pairwise in the resonator basis [13].
Starting with the linearized Hamiltonian, we obtain the coupled-mode equation where κ µ is the total decay rate of mode µ [13], δ µ is the detuning from the rotating frame (set by the group velocity of the soliton), andb in,µ are the bath operators. In Eq. (2), the first term accounts for the modal detunings (dispersion); the second term represents pair generation; the third term describes cross-phase modulation (XPM) and Bragg scattering; and the last term is the coupling to the bath. From Eq. (2), it is evident that the presence of multiple pump modes A j generates multimode coupling, resulting in collective comb dynamics that cannot be understood through pairwise mode analysis. An example of such a collective effect is the temporal envelope of the below-threshold comb, obtained from the steady-state solution of Eq. (2). The temporal shape of the quantum fluctuations does not mimic that of the mean field, but rather has a split shape (Fig. 1). The origin of this peculiar feature will, in a later section, be understood through the supermode decomposition.
To calculate the maximally squeezed supermodes of the system, we rewrite the Heisenberg equations in the basis of the quadrature operators of each mode: r (t) = (x 1 (t), ..., x n (t)| y 1 (t), ..., y n (t)) T , where x n = 1 √ 2 (a † n + a n ), and y n = i √ 2 (a † n − a n ). Input-output relations can be written for the quadrature operators as r out (t) = r in (t) + √ r (t), where is a diagonal matrix of the cavity decay rates κ µ . In the Fourier basis, the input and output fields are related by a transfer matrix, S(ω): which can be diagonalized by Euler decomposition [26,27]: The matrix D(ω) is diagonal with corresponding antisqueezing and squeezing levels associated with maximally squeezed orthogonal supermodes encoded in U (ω): the columns define the linear combination of quadratures, which can be mapped to a local oscillator for homodyne detection of the squeezing for ω, where U (ω) is real. This is always the case for ω = 0. For systems with terms corresponding to detuning or Bragg scattering, the ideal local oscillator configuration will depend on the Fourier frequency [26,28], e.g., the ideal local oscillator configuration for measuring maximum squeezing at zero Fourier frequency may measure sub-optimal squeezing across the rest of the spectrum.

MULTIMODE ANALYSIS OF COMB FORMATION
We begin with the supermode analysis of the stages of the microcomb that precede the soliton. The mean-field Hamiltonian coupling terms A m are obtained via an LLE simulation. We use system parameters consistent with silicon-carbide-on-insulator microrings [29,30] used in a recent experimental study of the quantum correlations in soliton crystals in the telecommunications C-band [13]: nonlinear coupling g 0 /2π = 3.4 Hz, a free spectral range (FSR) of D 1 /2π = 350 GHz, integrated quadratic dispersion D 2 /2π = 30 MHz, and a loaded quality factor of Q = 1.5 · 10 6 with critical coupling to the bus waveguide. The loss rate is assumed equal for all modes, denoted as κ (i.e., κ µ = κ). We note that these parameters are similar to soliton devices in many material platforms such as silicon nitride [8,11,24], lithium niobate [31,32], and tantala [33]. Figure 2(a) shows the evolution of intra-cavity mode amplitudes A m (t) under adiabatic red-tuning of the pump laser through the formation of primary and secondary combs. The formation of the primary comb produces a spatial rolls pattern, and the subsequent formation of the secondary comb generates spatiotemporal oscillations [ Fig. 2 A powerful feature of supermode analysis is the ability to reveal the multimode nature of threshold processes. Figure 2(c) illustrates the supermode analysis near the primary comb threshold: here, maximally squeezed supermodes consist of signal-idler pairs [14,15] described by supermodes 1 √ 2 (â −µ +â µ ) and i √ 2 (â −µ −â µ ), reflecting the amplitudes and phases of the local oscillators that could be used for homodyne detection. The well-known result of this single-pump Hamiltonian is that phase matching dictates which pair reaches threshold first, and predicts the spacing of the subsequent primary comb [11]. We now turn to the formation of the secondary comb: with the multiple nonzero-amplitude modes of the primary comb, the connectivity in the Hamiltonian increases beyond pairwise mode interactions. The squeezing spectrum just before the secondary-comb threshold is shown in Fig. 2(d). The maximally squeezed supermode approaches the threshold at a nonzero Fourier frequency ω = , indicating the presence of detuning in the multimode parametric gain process. To understand the origin of this detuning in the squeezing spectrum, we can examine the spectral composition â † µ (ω)â µ (ω) of the below-threshold modes, which, near threshold, is dominated by emission into the first supermode. The spectral composition [ Fig. 2(d)] reveals that the supermode consists of two quantum subcombs of equal and opposite detuning from the primary comb spacing. This detuning in the squeezing spectrum predicts the RF beatnote 2 that accompanies the formation of the secondary comb [11], giving rise to the spatiotemporal oscillation [34] seen in the numeric solution of the LLE [ Fig. 2(b)]. Thus, the supermode analysis sheds light on the formation process (contrasting the pair-wise mode analysis of secondary-comb formation [11]), and furthermore predicts the dynamics of the ensuing secondary comb.

SOLITON CRYSTAL AND ITS ANNIHILATION
After the formation of the secondary comb, the system enters the chaotic-modulation instability (MI)-regime. We will not discuss the MI state: the Hamiltonian is time dependent, and a subspace of below-threshold modes cannot be clearly delineated, thereby complicating the quantum analysis. To induce soliton-pair crystallization from the MI state, we introduce a −100 MHz perturbation at µ = +2 (as described in Ref. [8]). In the LLE simulation [ Fig. 3(a)], the MI state can be seen to end in a low-noise 2-FSR soliton crystal state, followed by an abrupt transition to the single soliton.
We analyze the quadrature squeezing for the below-threshold state as it evolves with detuning across the soliton crystal step. The state features significant spatiotemporal precession with respect to the resonator FSR, and care must be taken to perform the squeezing analysis in the stationary frame of the mean-field solution for each detuning. This procedure is described in Supplement 1. The inset of Fig. 3(a) shows the evolution of maximum squeezing, which exceeds 20 dB for all detunings. The maximally squeezed supermodes extend across the entire comb [ Fig. 3(b)] and resemble the Hermite-Gauss (HG) modes, the eigenmodes of SPOPO [18]. This reflects the all-to-all coupling in the Hamiltonian.
The evolution of the complete supermode eigenspectrum across the soliton crystal step is shown in Fig. 3(c). For all detunings, two supermodes show strong and comparable levels of squeezing while the rest have levels below 0.5 dB. This two-mode dominance of the eigenspectrum is a unique feature of the soliton crystal, unobserved in other multimode systems studied to date [18,23,27,35], although HG-like squeezed supermodes have been observed in soliton propagation through a χ (3) fiber [27]. The supermode structure of the squeezed vacuum explains the puzzling contrast between its temporal profile and that of the mean field, shown in the inset of Fig. 3(c): the strong contribution of the quasi-HG 1 supermode (whose Fourier transform to the time domain is also bi-modal) induces the temporal splitting of the squeezed vacuum pulse. Furthermore, an anti-crossing in the squeezing values of the two dominant supermodes is observed [ Fig. 3(c)], in contrast to the signal-idler pairs preceding primary comb formation, where no modal interaction is revealed in the squeezing spectra [ Fig. 2(c)]. The apparent interaction of the supermodes at the anti-crossing is further evidenced by the hybridization of the supermode shapes. Even at zero Fourier frequency, the phase profile of the quasi-HG 0 mode [ Fig. 3(b)] reveals a contribution of the odd quasi-HG 1 mode. This anti-crossing is observed universally for higher-order soliton crystals as well (Supplement 1 Fig. 2). Its physical significance will be investigated in a future work.
The end of the soliton crystal step is accompanied by asymptotic growth in the squeezing of one supermode, indicating that the dissociation of the soliton crystal state is a coherent threshold process. This can be anticipated from the fact that the passage from the crystallized two-soliton state to the single soliton results in the breaking of C 2 symmetry and correspondingly the onset of OPO in the modes µ = 1 (mod 2). In this regard, soliton crystal annihilation is unique among other state transitions, such as one-by-one disappearance of pulses in multi-soliton states [7] and transitions between soliton crystals with defects [24], where symmetry breaking does not occur. The annihilation of the soliton crystal is also the only known example (to our knowledge) of a threshold process that results in a reduction of the mean-field intensity: according to the LLE simulation [ Fig. 3(d)], one soliton disappears without energy transfer to the other.  (2) (τ ) for a single resonator mode (µ = 1) reveals the twin-comb signature through a temporal oscillation with period π/ .
As shown in Fig. 4(a), the strong squeezing across the step culminating in asymptotic growth is universally present across the soliton crystal existence condition. For all powers, the threshold is reached at a non-zero Fourier frequency [inset of Fig. 4(a)], indicating a strongly detuned threshold process. This is confirmed through calculating the spectral composition of the squeezed vacuum near threshold, which shows that every resonator mode displays a strongly split spectrum [ Fig. 4(b)]. The near-threshold state is thus composed of twin quantum frequency combs offset by from the soliton crystal rotating frame. Equal intensities of the twin combs is guaranteed by energy and momentum conservation. Twin combs would be directly observable in the second-order autocorrelation g (2) (τ ), manifesting as a temporal oscillation with a period of π/ [ Fig. 4(c)]. The autocorrelation peaks at g (2) (0) = 3 and exhibits significant coherence broadening, as expected for near-threshold OPO [13,36].
A degenerate parametric oscillator reaches parametric oscillation in a continuous process in which damping precludes the formation of a coherent cat state and instead produces a classical mixture of coherent states [25], manifesting as random phase selection of the above-threshold OPO. If the soliton crystal did not exhibit quantum twin comb behavior ( = 0) at threshold, the supermode analysis would lend itself to an analogous picture: the post-threshold state is a classical mixture of either one of the soliton pulses disappearing, corresponding to the modes µ = 1 (mod 2) possessing phase of zero or π with respect to modes µ = 0 (mod 2). However, since = 0, such a simple interpretation is not possible: the passage through threshold must be accompanied by the spectral collapse of the quantum twin comb into a rotating frame. This cannot be explained within the framework of the linearized model, which we show here to predict unbounded growth of twin combs at non-zero . The annihilation of the soliton crystal thus represents a clear opportunity for experimental and theoretical exploration of the breakdown of the linearization assumption in nonlinear Kerr resonators.

SOLITON CRYSTAL SQUEEZED LIGHT SOURCE
Until now, we have considered squeezing in the absence of parasitic loss channels-all of the quantum light generated inside the resonator is collected with unity efficiency. While in principle one may realize resonators with an arbitrarily high waveguidecoupling rate relative to the intrinsic loss, peak-efficiency Kerr soliton devices operate near the critical coupling point, which limits the outcoupled squeezing to 3 dB. Soliton operation in the over-coupled regime is associated with a severe increase in power requirements, as the OPO threshold scales quadratically with resonator losses. For instance, the outcoupling-limited squeezing of 10 dB (15 dB) requires an escape efficiency of 90% (96.8%), and corresponds to an increase of the OPO threshold by 25 (250) times, putting into question the practicality of this approach. In this section, we present a device architecture that overcomes this limitation.
The photonic architecture based on a 2-FSR soliton crystal is shown in Fig. 5(a). The proposed device consists of the squeezing resonator, engineered for anomalous dispersion to support soliton formation, and critically coupled to the bus waveguide (κ i = κ c = κ/2) for efficient in-coupling of the pump. The auxiliary resonator is designed to have an FSR twice larger than that of the squeezing resonator, and its modes overlap in frequency with the odd [µ = 1 (mod 2)] azimuthal modes of the squeezing resonator. The auxiliary resonator is strongly over-coupled to the drop waveguide (coupling rate κ aux κ i ), and is coupled to the squeezing resonator with strength J such that κ aux J > κ. This corresponds to the regime known in cavity quantum electrodynamics as Purcell enhancement: the odd modes of the squeezing resonator are coupled to the auxiliary resonator and experience a decay rate (κ out ) into the drop waveguide. If the magnitude of Purcell enhancement (κ out /κ i ) is well below the finesse F of the squeezing resonator (typical F is 10 3 −10 4 ), the even modes in the squeezing resonator are unaffected by the auxiliary resonator. Thus, rather than disturb the formation of the soliton crystal, the auxiliary ring stabilizes its formation through the suppression of undesirable OPO processes. The choice of J and κ aux thus provides control over the outcoupling rate of the below-threshold modes without negatively impacting the above-threshold state. The magnitude of squeezing in the drop waveguide as a function of κ out /κ is shown in Fig. 5(b). The squeezing calculation in presence of loss is presented in Supplement 1. In the regime of small κ out , the squeezing is limited by the outcoupling efficiency into the drop waveguide, and increases with κ out . The intrinsic squeezing, however, drops with increasing κ out , since the classical soliton crystal state remains unchanged and thus results in weaker effective drive and increased distance to threshold. In the regime where κ out dominates other losses, the outcoupled squeezing is limited by the intrinsic squeezing of the system.
The same conditions responsible for the formation of the twin quantum frequency combs near the soliton crystal annihilation threshold also provide resilience of the squeezing strength against the addition of the outcoupling rate κ aux , rendering the soliton crystal an attractive source of squeezed light. This is illustrated by contrast with the squeezing of a typical, non-detuned squeezed source [ Fig. 5(b)]. Ramping up κ out without altering the pump power, the maximum outcoupled squeezing is ≈ 3 dB at κ out = 2κ, beyond which it rapidly decays. In contrast, the peak outcoupled squeezing of the below-threshold soliton crystal state reaches its maximum of 10 dB at κ out = 20κ. The resilience of squeezing to added losses is the consequence of the detuning of the squeezing process: the broadening of the below-threshold modes associated with growing κ out reduces the effective drive strength but simultaneously reduces the loss-normalized detuning of the belowthreshold modes with respect to the pump modes. The detuning thus acts as a "squeezing strength" reservoir, and the auxiliary resonator enables the extraction of 10 dB of outcoupled squeezing without increasing the pump power.
We have so far described the effect of the auxiliary resonator while holding constant the power of the pump laser. Since, as noted above, the auxiliary resonator stabilizes the soliton crystal state, the state can be synthesized in the squeezing resonator at much higher pump power to achieve squeezing levels limited by the escape efficiency. We illustrate this in Fig. 5(c) for escape efficiency η = 0.98 (κ out = 50 κ), which corresponds to the outcouplinglimited squeezing of ≈ 17 dB in the waveguide. In this condition, the soliton crystal may be captured with a 100 mW pump, over a greatly extended detuning range. The outcoupled squeezing grows steadily along the soliton crystal step (consistent with the growing above-threshold comb power), asymptotically approaching the outcoupling limit.
A few additional characteristics of the photonic molecule architecture, as they pertain to the realization of a practical highly squeezed multimode source, are worth noting.
Construction of the local oscillator. The photonic molecule conveniently separates the squeezed vacuum and the above-threshold soliton crystal into separate waveguides. The latter may then be manipulated independently of the squeezed light to construct the local oscillator. Since the temporal bandwidth of the generated squeezed vacuum matches closely that of the above-threshold modes, the local oscillator may be readily generated by first electrooptically interleaving [37] the above-threshold soliton crystal, followed by pulse-shaping the resulting state to retain coherent light with the phases and intensities corresponding to the desired supermode. This process is illustrated in Fig. 5(a). The entirety of the local oscillator preparation may be implemented on-chip, with the recent advances in on-chip electro-optic frequency comb generation [38], low-loss phase-shifting [39], and, if necessary, amplification [40].
Squeezing bandwidth. While the above-threshold modes have linewidths of the order of 100 MHz (enabling low-power OPO and formation of the soliton crystal), the bandwidth of the squeezed light is dictated by the total linewidth of the belowthreshold modes, and thus is magnified by a factor of κ out /κ. The photonic molecule thus may be used to generate squeezed light of bandwidth potentially much greater than that of the squeezing resonator [ Fig. 5(d)]. The local oscillator compositions corresponding to the maximum squeezing at ω = 0 in Fig. 5(d) are shown in Fig. 5(e).
Fourier frequency of maximal squeezing. A local oscillator with optimal squeezing may be prepared in the corresponding supermode as long as U (ω) is real. However, this is guaranteed only for ω = 0. In the case where U (ω) is complex, a local oscillator stationary in the rotating frame cannot be used to measure the maximum squeezing [23]. However, for κ out κ, maximum squeezing shifts to ω = 0, because the loss-normalized system detunings are reduced. This fortuitously renders the squeezing in the photonic molecule configuration amenable to the straightforward local oscillator measurement.

DISCUSSION
We have applied a multimode squeezing analysis to the 2-FSR soliton crystal microcomb, which captures the essential multimode features of the below-threshold state. We have considered the simplest realistic configuration for soliton crystal formation. Additional work must be done to analyze the effects of Raman scattering [27,41] and higher-order dispersion on the squeezing structure. Application of the squeezed supermode decomposition to higher-order soliton crystal states (N circulating pulses) results in N − 1 orthogonal supermodes for each quasi-HG order. For instance, a 7-FSR crystal features six quasi-HG 0 and six quasi-HG 1 supermodes of nearly identical squeezing levels (see Supplement 1). Thus, the squeezing spectrum of a soliton crystal source can be tailored beyond two prominent supermodes. However, the pattern of symmetry breaking in the annihilation of higher-FSR soliton crystals is expected to be qualitatively different, as multiple pulses (and thus multiple decay pathways) are present.
To experimentally test the findings of this work, several techniques may be applied. Measuring the evolution of second-order photon correlations [13] across the step would validate the existence of the threshold process as well as the formation of twin quantum combs. Ultrafast imaging [7] could be used to observe the soliton crystal decay and potentially identify whether the macroscopic state evolution reveals the distinction between the annihilation of soliton crystals and non-crystallized multisoliton states. Finally, the photonic molecule device proposed here for direct measurement of out-coupled squeezing can be readily fabricated using current experimental capabilities [8,13,24,[31][32][33].
This study may also serve as a starting point for exploring the process by which the linearized model breaks down. The unbounded growth of twin quantum combs in a soliton crystal transitioning to the single-soliton state is an unreconciled discontinuity within the linearized model. The experimental observation of this discontinuity may give insights into the quantum parametric processes in Kerr microcombs beyond the linearized model [42].
Funding. Defense Advanced Research Projects Agency (LUMOS and QuICC Programs).