Modulation Instability-Induced Multimode Squeezing in Quadratic Frequency Combs

TL;DR

This paper demonstrates how modulation instability in lithium niobate microring resonators can induce multimode squeezing in quadratic frequency combs, advancing integrated quantum photonics for quantum information and measurement.

Contribution

It introduces a universal framework for analyzing multimode squeezing and reveals modulation instability as a key mechanism for multimode entanglement in quadratic frequency combs.

Findings

Modulation instability enables formation of frequency combs and multimode squeezing.

A bipartite entanglement criterion for quantum frequency combs is proposed.

On-chip multimode squeezing facilitates collective noise suppression.

Abstract

Lithium niobate (LN) microring resonators, characterized by an exceptionally high second-order nonlinear coefficient and superior electro-optic tunability, serve as an outstanding platform for the precise control of integrated quantum frequency combs (QFCs). In this study, we introduce a bipartite entanglement criterion to investigate the pairwise entanglement characteristics of QFCs generated via the spontaneous parametric down-conversion (SPDC) process in lithium niobate microring resonators operating below threshold. Furthermore, we propose a universal framework for analyzing multimode squeezing in quadratic frequency combs, enabling the realization of ultrabroadband and high-degree multimode squeezing. We further reveal the underlying physical mechanism: modulation instability (MI), regulated by temporal walk-off control, not only enables the formation of frequency combs but also induces multimode squeezing in the corresponding resonant modes. This study uncovers the previously unexplored role of on-chip multimode squeezing in quadratic frequency combs while facilitating collective noise suppression across multiple modes, thus holding substantial potential for advancing quantum precision measurement and quantum information processing.

Equations

$$A_{\mathrm{eff}}=\frac{{\left({\iint }_{-\infty }^{+\infty }|F(y,z){|}^{2}dydz\right)}^2}{{\iint }_{-\infty }^{+\infty }|F(y,z){|}^{4}dydz}.$$

$${\omega }_l={\omega }_0+D_{1}l+\frac{D_2}{2}{l}^2+\cdots ={\omega }_0+\sum _{k=1}^{\infty }{D}_k\frac{l^k}{k!},$$

$$D_{\text{int}}={\omega }_l-{\omega }_0-D_{1}l,$$

$${\Delta }_p={\omega }_p-{\Omega }_p.$$

$$\Delta =(2{\omega }_0-{\omega }_p+{\Delta }_p)/2.$$

$$\gamma =|t_2{|}^2{\nu }_f=(1-|t_1{|}^2){\nu }_f,\mu ={\alpha }_{c}L{\nu }_f,$$

$${\alpha }_c\approx \frac{Q_{0}R{D}_1}{{\omega }_0},$$

$$Q={\omega }_0{\tau }_p=\frac{{\omega }_0}{\Gamma },$$

$$\mu =\frac{{\omega }_0}{Q_0},\gamma =\Gamma -\mu =\frac{{\omega }_0}{Q_{\mathrm{ex}}},$$

$$\frac{1}{Q}=\frac{1}{Q_0}+\frac{1}{Q_{\mathrm{ex}}},$$

$$\Delta k=\frac{{\omega }_{p}n({\omega }_p)}{c}-\frac{{\omega }_{s}n({\omega }_s)}{c}-\frac{{\omega }_{i}n({\omega }_i)}{c}=0.$$

$${\hat{x}}_j=\frac{1}{\sqrt{2}}({\hat{a}}_j^{†}{e}^{\mathrm{i}{\theta }_j}+{\hat{a}}_j{e}^{-\mathrm{i}{\theta }_j}),{\hat{y}}_j=\frac{1}{\sqrt{2}}(\mathrm{i}{\hat{a}}_j^{†}{e}^{\mathrm{i}{\theta }_j}-\mathrm{i}{\hat{a}}_j{e}^{-\mathrm{i}{\theta }_j}).$$

$${\left(\delta {\hat{x}}_s,\delta {\hat{x}}_i,\delta {\hat{y}}_s,\delta {\hat{y}}_i\right)}^{\mathrm{T}}=P{\left(\delta {\hat{a}}_s,\delta {\hat{a}}_s^{†},\delta {\hat{a}}_i,\delta {\hat{a}}_i^{†}\right)}^{\mathrm{T}}.$$

$${\hat{x}}_{\pm }=\frac{1}{\sqrt{2}}({\hat{x}}_s\pm {\hat{x}}_i),{\hat{y}}_{\pm }=\frac{1}{\sqrt{2}}({\hat{y}}_s\pm {\hat{y}}_i),$$

$$\delta {\hat{X}}_{\pm }={\left(\delta {\hat{y}}_{+},\delta {\hat{x}}_{+},\delta {\hat{y}}_{-},\delta {\hat{x}}_{-}\right)}^{\mathrm{T}}=G{\left(\delta {\hat{x}}_s,\delta {\hat{x}}_i,\delta {\hat{y}}_s,\delta {\hat{y}}_i\right)}^{\mathrm{T}}.$$

$$S_{{\hat{X}}_{\pm }}(\omega )$$

$$=\frac{T_r\cdot S_a(\omega )\cdot T_r^{\mathrm{T}}+(T_r\cdot S_a(\omega )\cdot T_r^{\mathrm{T}}{)}^{\mathrm{T}}}{2},$$

$$C_s=((\Delta {\hat{y}}_{+}{)}^2+(\Delta {\hat{x}}_{-}{)}^2-|O|\ge 0,$$

$$\hat{H}$$

$$-\frac{1}{2}\sum _{k,n,j}\hslash g_0\left({\hat{a}}_k{\hat{a}}_n{\hat{b}}_j^{†}+{\hat{a}}_k^{†}{\hat{a}}_n^{†}{\hat{b}}_j\right)\delta [{\Omega }_j-{\Omega }_k-{\Omega }_n],$$

$$\{\begin{array}{rl}\frac{\mathrm{d}{\alpha }_u}{\mathrm{d}t}& =\left[-\mathrm{i}({\omega }_u-{\Omega }_u)-\Gamma \right]{\alpha }_u+\mathrm{i}{g}_0\sum _{k,j}{\beta }_j{\alpha }_k^{\ast }\delta [{\Omega }_j-{\Omega }_k-{\Omega }_u],\\ \frac{\mathrm{d}{\beta }_v}{\mathrm{d}t}& =\left[-\mathrm{i}({\omega }_v-{\Omega }_v)-\Gamma \right]{\beta }_v+\mathrm{i}{g}_0\sum _{k,n}{\alpha }_n{\alpha }_k\delta [{\Omega }_v-{\Omega }_k-{\Omega }_n]\\ & +\sqrt{2\gamma }{\beta }^{in}\delta [{\Omega }_v-{\Omega }_p].\end{array}$$

$$\{\begin{array}{rl}\frac{\mathrm{d}a}{\mathrm{d}t}& =\left(-\mathrm{i}\Delta -\Gamma \right)a-\mathrm{i}\frac{k_1^{\prime \prime }L{\nu }_f}{2}\frac{{\partial }^{2}a}{\partial {\tau }^2}+\mathrm{i}{g}_{0}b\star a^{\ast },\\ \frac{\mathrm{d}b}{\mathrm{d}t}& =\left(-\mathrm{i}{\Delta }_p-\Gamma \right)b-\Delta k^{\prime }L{\nu }_f\frac{\partial b}{\partial \tau }-\mathrm{i}\frac{k_2^{\prime \prime }L{\nu }_f}{2}\frac{{\partial }^{2}b}{\partial {\tau }^2}\\ & +\mathrm{i}{g}_{0}a\star a+\sqrt{2\gamma }{B}_{\text{in}},\end{array}$$

$$\{\begin{array}{rl}\frac{\mathrm{d}\delta {\hat{a}}_u}{\mathrm{d}t}& =\left[-\mathrm{i}({\omega }_u-{\Omega }_u)-\Gamma \right]\delta {\hat{a}}_u+\sqrt{2\gamma }\delta {\hat{a}}_u^{\text{in}}+\sqrt{2\mu }\delta {\hat{a}}_u^{\text{loss}}\\ & +\mathrm{i}{g}_0\sum _{k,j}({\alpha }_k^{\ast }\delta {\hat{b}}_j+{\beta }_j\delta {\hat{a}}_k^{†})\delta [{\Omega }_j-{\Omega }_k-{\Omega }_u],\\ \frac{\mathrm{d}\delta {\hat{b}}_v}{\mathrm{d}t}& =\left[-\mathrm{i}({\omega }_v-{\Omega }_v)-\Gamma \right]\delta {\hat{b}}_v+\mathrm{i}{g}_0\sum _{k,n}{\alpha }_k\delta {\hat{a}}_n\delta [{\Omega }_v-{\Omega }_k-{\Omega }_n]\\ & +\sqrt{2\gamma }\delta {\hat{b}}_v^{\text{in}}+\sqrt{2\mu }\delta {\hat{b}}_v^{\text{loss}}.\end{array}$$

$$\{\begin{array}{rl}\delta \hat{\mathbf{A}}& =\left(\begin{array}{ccc}\delta {\hat{a}}_{-N/2+1},& \cdots & ,\delta {\hat{a}}_{N/2}\end{array}\right),\\ \delta \hat{\mathbf{B}}& =\left(\begin{array}{ccc}\delta {\hat{b}}_{-N/2+960},& \cdots & ,\delta {\hat{b}}_{N/2+959}\end{array}\right),\\ \hat{\delta }& =(\delta \hat{\mathbf{A}},\delta \hat{\mathbf{B}},\delta {\hat{\mathbf{A}}}^{†},\delta {\hat{\mathbf{B}}}^{†}{)}^{\text{T}}.\end{array}$$

$$\frac{\mathrm{d}\hat{\delta }}{\mathrm{d}t}=M_a\cdot \hat{\delta }+U_a^{\mathrm{in}}\cdot {\hat{\delta }}^{\mathrm{in}}+U_a^{\mathrm{loss}}\cdot {\hat{\delta }}^{\mathrm{loss}},$$

$$M_a=\left(\begin{array}{cccc}{\chi }_A& \mathrm{i}{g}_{0}N{F}_N^{-1}\text{diag}((F_N\alpha {)}^{\ast })F_N& \mathrm{i}{g}_{0}N{F}_N^{-1}\text{diag}(F_N\beta )F_N^{\ast }& 0E_N\\ \mathrm{i}{g}_{0}N{F}_N^{-1}\text{diag}(F_N\alpha )F_N& {\chi }_B& 0E_N& 0E_N\\ -\mathrm{i}{g}_{0}N{F}_N^{-1}\text{diag}(F_N{\beta }^{\ast })F_N^{\ast }& 0E_N& {\chi }_A^{\ast }& -\mathrm{i}{g}_{0}N{F}_N^{-1}\text{diag}(F_N^{\ast }\alpha )F_N\\ 0E_N& 0E_N& -\mathrm{i}{g}_{0}N{F}_N^{-1}\text{diag}(F_N{\alpha }^{\ast })F_N& {\chi }_B^{\ast }\end{array}\right).$$

$$\hat{L}(t)=\sum _{l=-N/2+1}^{N/2}{\eta }_l{\hat{a}}_l(t)+\sum _{l=-N/2+960}^{N/2+959}{\eta }_l{\hat{b}}_l(t).$$

$$S_x(\omega )=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}{E}_{2N}& E_{2N}\\ -\mathrm{i}{E}_{2N}& \mathrm{i}{E}_{2N}\end{array}\right)\cdot M^{\text{in}}(\omega )\cdot \frac{1}{\sqrt{2}}\left(\begin{array}{cc}{E}_{2N}& -\mathrm{i}{E}_{2N}\\ E_{2N}& \mathrm{i}{E}_{2N}\end{array}\right).$$

$$\{\begin{array}{rl}{V}^{+}(\omega )& ={\int }_{-\infty }^{\infty }⟨{\hat{L}}_{\text{out}}^{(+)}(t){\hat{L}}_{\text{out}}^{(+)}(t+\tau )⟩e^{-\mathrm{i}\omega \tau }d\tau ,\\ V^{-}(\omega )& ={\int }_{-\infty }^{\infty }⟨{\hat{L}}_{\text{out}}^{(-)}(t){\hat{L}}_{\text{out}}^{(-)}(t+\tau )⟩e^{-\mathrm{i}\omega \tau }d\tau .\end{array}$$

$$\{\begin{array}{rl}|{\eta }_{l,a}|& =\sqrt{|U(\omega {)}_{l,k}{|}^2+|U(\omega {)}_{l+2N,k}{|}^2},\\ {\phi }_{l,a}& =\text{atan2}\left(U(\omega {)}_{l+2N,k},U(\omega {)}_{l,k}\right),\\ |{\eta }_{l,b}|& =\sqrt{|U(\omega {)}_{l+N,k}{|}^2+|U(\omega {)}_{l+3N,k}{|}^2},\\ {\phi }_{l,b}& =\text{atan2}\left(U(\omega {)}_{l+3N,k},U(\omega {)}_{l+N,k}\right).\end{array}$$

$${\hat{H}}_0=\hslash \left({\omega }_s{\hat{a}}_s^{†}{\hat{a}}_s+{\omega }_i{\hat{a}}_i^{†}{\hat{a}}_i+{\omega }_p{\hat{b}}^{†}\hat{b}\right).$$

Full text

††journal: Chaos Solitons $\&$ Fractals

\affiliation

[add1]organization=State Key Laboratory of Photonics and Communications, School of Information Science and Electronic Engineering $\&$ School of Integrated Circuits, addressline=Shanghai Jiao Tong University, city=Shanghai, postcode=200240, country=China

\affiliation

[add2]organization=State Key Laboratory of Precision Spectroscopy, addressline=East China Normal University, city=Shanghai, postcode=200062, country=China

\affiliation

[add3]organization=Microwave Photonics Technology Key Laboratory of Sichuan Province, city=Chengdu, state=Sichuan Province, country=China

\affiliation

[add4]organization=School of Physics and Optoelectronic Engineering, Institute of Laser Engineering, addressline=Beijing University of Technology, city=Beijing, postcode=100124, country=China

\affiliation

[add5]organization=College of Electronic Science, National key laboratory of satellite navigation technology, addressline=National University of Defense Technology, city=Changsha, postcode=410073, country=China

\affiliation

[add6]organization=National Key Laboratory of Metrology and Calibration, addressline=Beijing Changcheng Institute of Metrology $\&$ Measurement, city=Beijing, postcode=100095, country=China

ttfootnotetext: These authors contributed equally to this work.

{graphicalabstract}

{highlights}

Lithium niobate microring quantum frequency combs achieve controllable entanglement

Unified theory enables full-regime multimode squeezing optimization in OPO

Modulation instability drives simultaneous comb-squeezing in quadratic systems

Unveils quantum functionality for scalable information and precision metrology

1 Introduction

Squeezed optical states, as indispensable nonclassical resources in quantum information science, enable noise suppression beyond the standard quantum limit [1, 2], thereby facilitating advances in secure quantum communication and computation. The rapid progress in integrated photonics is catalyzing a shift toward chip-scale quantum optical systems [3], where microring resonators have emerged as crucial platforms due to their strong optical confinement and enhanced nonlinear interactions [4]. Recent breakthroughs across diverse material platforms have significantly advanced microring-based quantum light sources. In silicon nitride (SiN) microrings, below-threshold spontaneous four-wave mixing has been employed to generate quadrature squeezed vacuum states and photon-number-difference squeezed states [5], as well as to realize heralded single-photon sources [6] and scalable multi-user quantum networks [7]. Silicon carbide (SiC) microrings utilize soliton microcomb dynamics to achieve multimode entanglement [8], while high-purity quantum sources have also been demonstrated in gallium nitride [9], aluminum gallium arsenide [10], and integrated nanophotonic platforms [11]. Complementing these advances, quantum frequency comb control techniques have evolved—from non-equilibrium driving of strongly coupled photonic dimer systems [12] to parity-time (PT) symmetric mode-selective pumping strategies [13]—further underscoring the central role of quantum frequency combs in emerging applications such as quantum computing, communication, and precision metrology.

Lithium niobate (LN) offers distinct advantages over materials primarily governed by third-order ($\chi^{(3)}$) nonlinearities, such as SiN, owing to its strong second-order ($\chi^{(2)}$) nonlinear response [14]. This intrinsic property enables direct exploitation of $\chi^{(2)}$ processes, particularly SPDC, allowing for the efficient generation of entangled photon pairs with a significantly lower pump threshold than that required in SiN-based platforms [15, 16]. In addition, LN facilitates quantum frequency conversion of single photons between the telecom band and the visible or mid-infrared range [17, 18], enabling spectral interfacing across disparate quantum systems and meeting the stringent compatibility requirements of multimode quantum architectures. The heterogeneous integration of LN thin films with silicon-based photonic circuits [19] further paves the way for multifunctional quantum photonic chips. Moreover, lithium niobate exhibits a pronounced Pockels effect, enabling high-speed and low-loss electro-optic modulation, which plays a vital role in optical communication systems. Driven by these advantages, this work focuses on the theoretical investigation of multimode squeezing in quadratic frequency combs generated by LN microring resonators.

In this study, we develop an integrated LN microring resonator platform that enables the generation of QFCs through precise dispersion and coupling engineering [20]. To elucidate the multimode squeezing mechanisms of quadratic frequency combs, we categorize our analysis into below-threshold and above-threshold regimes. Since below-threshold squeezing inherently occurs in a pairwise fashion between modes with quantum correlations in quadratures, it gives rise to bipartite entanglement. We therefore employ a bipartite entanglement criterion to rigorously quantify the degree of entanglement between each mode pair. Our results show that under below-threshold conditions, increasing the pump power enables a controllable and continuous redshift of the squeezed frequency components. We further investigate the frequency-dependent squeezing characteristics and reveal the interaction between the detection frequency and the optimal readout phase [21].

Above the optical parametric oscillator (OPO) threshold, the coupled-mode equations for the signal-idler pair break down due to cascaded $\chi^{(2)}$ nonlinear process. We therefore extend the model to a multimode formalism and derive quadratic coupled mean-field equations to describe the classical evolution of the optical frequency comb (OFC). Pronounced temporal walk-off in OPOs [22, 23] induces a novel MI band [24], wherein parametric amplification of the subharmonic field leads to enhanced pulse compression and synchronized temporal drift [25]. As the classical mean-field framework fails to capture the quantum correlations underlying threshold behavior [26], we employ a supermode basis—quadrature-weighted combinations of frequency modes [27]—to characterize multimode squeezing and intermodal quantum correlations [28]. In contrast to synchronously pumped OPOs (SPOPOs), where fluctuations are externally driven [29], quantum noise in our quadratic frequency comb arises intrinsically from $\chi^{(2)}$ nonlinearity, linking the below-threshold quantum state to the mean-field solution.

During the detuning scan that culminates in soliton crystal formation, we observe pronounced quadrature squeezing, with ultrabroadband multimode correlations emerging between the pump and subharmonic modes. This phenomenon originates from the intricate interplay between optical nonlinearity and dispersion, wherein minor perturbations are exponentially amplified, yielding substantial MI gain at selective frequencies. While MI underpins the formation of optical frequency combs, it concurrently drives extensive multimode squeezing across resonant modes. In contrast to the multipartite entanglement analysis in Ref. [30], our work provides a unified characterization of quadrature squeezing and spectral correlations across 800×2 modes that potentially contribute to the underlying multimode dynamics. We emphasize that a complete description must account for both pump and subharmonic field modes, as neglecting inter-band coupling can obscure critical squeezing phenomena. The general framework we develop for analyzing multimode squeezing in quadratic frequency combs is broadly applicable and holds strong potential for emerging applications such as high-throughput optical interconnects [31] and photonic machine learning accelerators [32].

2 Microring Resonator Simulation Model

Lithium niobate, renowned for its superior $\chi^{(2)}$ nonlinearity, enables a high-frequency pump photon ($\Omega_{p}$) to spontaneously split into two lower-frequency photons—a process known as SPDC. These resultant photons, designated as the signal ($\Omega_{s}$) and idler ($\Omega_{i}$) photons, form an entangled pair. Notably, SPDC is exclusively initiated by vacuum fluctuations and inherently adheres to both energy conservation and phase matching conditions.

In this study, we introduce an advanced on-chip microcavity architecture grounded in OPO theory, which ingeniously incorporates an enhanced coupling module to facilitate highly efficient nonlinear interactions. As depicted in Fig. 1(a), the micro-ring resonator comprises a ring waveguide coupled to a bus waveguide, with a z-cut LN ridge waveguide serving as the nonlinear medium. The waveguide is engineered with a silicon dioxide (SiO2) lower cladding to suppress substrate leakage and enhance vertical optical confinement, while an air upper cladding is employed to further tailor the mode field distribution. The optical field circulates along the azimuthal direction of the micro-ring resonator, and the effective cross-sectional area ($A_{\mathrm{eff}}$) of the ring waveguide [33], quantified by Equation (3), characterizes the extent of field localization within the resonator; a reduced effective mode area signifies superior confinement of the optical field.

[TABLE]

Here, $F(y,z)$ delineates the modal distribution within the LN and SiO2 material system. In our theoretical framework, this distribution is presumed invariant along the longitudinal propagation direction within the resonator, remaining independent of the evolution of the optical field along the cavity. This assumption is instrumental in facilitating the application of the separation of variables technique for solving the coupled mode equations.

The coupling region and the corresponding input–output schematic are illustrated in Fig. 1(b). The geometry of this region directly determines the coupling rate between the ring resonator and the bus waveguide, enabling precise control and forming a solid foundation for rigorous analysis of the resonator’s input–output characteristics. Figure 1(c) shows the cross-sectional profile of the ridge waveguide, with all relevant structural parameters clearly annotated. In this work, we focus on the electromagnetic properties of the fundamental transverse electric (TE) mode. The resonance profile of the cavity follows a Lorentzian lineshape. Due to the frequency dependence of the refractive index $n(\omega)$, resonance peaks exhibit spectral shifts as a function of optical frequency. Accurate control over key cavity parameters—including detuning, dispersion, coupling strength, and loss rate—within the designed platform requires a comprehensive and physically consistent theoretical framework.

2.1 Detuning and Dispersion

This section delves into the fundamental mechanisms of detuning and dispersion. We introduce the relative mode number $l$ ($l\in\mathbb{Z}$) to index the modes adjacent to the half-harmonic mode $\omega_{0}$ (corresponding to $l=0$). By performing a Taylor expansion of the resonance mode frequency around $\omega_{0}$, we obtain:

[TABLE]

where $D_{1}=2\pi\nu_{f}$, with $\nu_{f}$ representing the free spectral range (FSR). The parameter $D_{2}$ governs the group velocity dispersion  [34], where $D_{2}>0$ corresponds to anomalous dispersion, while $D_{2}<0$ characterizes normal dispersion. Higher-order dispersion terms ($D_{k}$ for $k\geq 3$) are neglected in this model. The integrated dispersion is defined as:

[TABLE]

which can be accurately approximated by a quadratic polynomial in the vicinity of $\omega_{0}$ [35].

The cold-cavity pump detuning is given by:

[TABLE]

Meanwhile, we set the detuning at the subharmonic mode to $\Delta=\omega_{0}-\Omega_{0}$. Here, $\Omega_{l}$ denote the comb frequency corresponding to the mode number $l$. Owing to $\Omega_{p}=2\Omega_{0}$, we obtain:

[TABLE]

2.2 Coupling, Loss and Gap

We analyze the fundamental relationship between coupling rate ($\gamma$), intrinsic loss rate ($\mu$), and resonator quality factor ($Q$). Following the virtual channel framework in Refs. [36, 37], intrinsic losses are modeled as a phantom channel with beam-splitter-like behavior, under the constraint $\gamma\ll\nu_{f}$, $\mu\ll\nu_{f}$, where $\nu_{f}$ denotes the free spectral range.

As illustrated in Fig. 1(b), the coupling and loss rates are given by:

[TABLE]

where $|t_{1}|^{2}$ and $|t_{2}|^{2}$ are the energy reflection and transmission coefficients, respectively, and $L=2\pi R$ is the resonator perimeter with $R=(D_{\text{in}}+D_{\text{out}})/4$. Here, $D_{\text{in}}$ and $D_{\text{out}}$refer to the microring’s inner and outer diameters, respectively. The absorption coefficient $\alpha_{c}$ (in m*-1*) is approximated as

[TABLE]

with $Q_{0}$ the intrinsic quality factor and $\omega_{0}$ the resonance frequency.

To evaluate the effect of external coupling, we fix $b_{2}=0$ and vary the ring–bus gap. The ratio $r=\gamma/\mu$ defines the coupling regime: under-coupled ($r<1$), critically coupled ($r=1$), and over-coupled ($r>1$). We focus on the over-coupled regime, which enables enhanced power extraction despite reduced intracavity field strength.

In LN resonators, the loaded quality factor governs field buildup and energy confinement:

[TABLE]

where $\tau_{p}$ is the photon lifetime and $\Gamma$ the total linewidth. The intrinsic and coupling rates are

[TABLE]

with $Q_{\mathrm{ex}}$ the external quality factor [38]. The total $Q$ satisfies

[TABLE]

linking the coupling gap to $Q_{\mathrm{ex}}$ and thereby enabling tunable external coupling.

2.3 Actual Simulation

Our objective is to engineer a dispersion-flattened waveguide by optimizing the structural parameters of the LN-integrated microcavity, such as its geometry and dimensions, to achieve phase matching in the SPDC process. This condition is fulfilled when the following equation holds:

[TABLE]

Here, $\Delta k$ quantifies the phase mismatch. For our designed micro-ring resonator, at 1560 nm the first-order dispersion is given by $k^{\prime}_{1}=\left.\frac{dk}{d\omega}\right|_{\omega_{0}}=7879.9ps/m,$ and the group velocity dispersion is $k^{\prime\prime}_{1}=\left.\frac{d^{2}k}{d\omega^{2}}\right|_{\omega_{0}}=-0.0219ps^{2}/m.$ Therefore, the FSR we have set can be expressed as: $\nu_{f}=1/(k^{\prime}_{1}L)=201.976$GHz. At 780 nm, the corresponding values are $k^{\prime}_{2}=\left.\frac{dk}{d\omega}\right|_{\omega_{p}}=8085.6ps/m,$ and $k^{\prime\prime}_{2}=\left.\frac{d^{2}k}{d\omega^{2}}\right|_{\omega_{p}}=0.3624ps^{2}/m.$ We employ COMSOL Multiphysics to simulate the complete waveguide geometry and its integrated dispersion profile. A fifth-order polynomial fit is then applied to closely match the integrated dispersion $D_{\text{int}}$ as depicted in Figure 2. Additionally, Lumerical FDTD simulations are conducted to elucidate the relationship between the gap and the coupling coefficient. Assuming an intrinsic quality factor of $Q_{0}=3.7\times 10^{6},$ the intrinsic loss rate is computed as $\mu=3.27\times 10^{8}\ \text{rad/s}.$ With a coupling ratio $r=\frac{\gamma}{\mu}=1.222,$ indicative of an over-coupled configuration, the external quality factor is determined to be $Q_{ex}=Q_{0}/r=3.03\times 10^{6}.$ Based on the FDTD simulation outcomes, the gap is calculated to be 490 nm.

Figure 1(c) illustrates the structure of an LN microring resonator, featuring a radius of $R=100\,\mu\text{m}$, with both the waveguide and bus widths set at $W_{R}=W_{B}=2000\,\text{nm}$. The outer height is $H=600\,\text{nm}$, while the inner height is $h=410\,\text{nm}$, and the design incorporates an angle of $\theta=75^{\circ}$. Simulation results for this configuration indicate a pump mode frequency ($\omega_{p}$) of 2.41406×$10^{15}$ rad/s, a fundamental mode (mode 0) frequency ($\omega_{0}$) of 1.20688×$10^{15}$ rad/s, and an effective area ($A_{\text{eff}}$) of $0.997\,\mu\text{m}^{2}$.

3 Multimode Squeezing

In this section, we present a comprehensive analysis of the multimode squeezing properties of quadratic frequency combs based on LN resonators. Below the OPO threshold, two-mode vacuum squeezing arises in correlated pairs, allowing the use of the Duan bipartite entanglement criterion to quantitatively evaluate the entanglement and squeezing of multiple two-mode squeezed states. Above the threshold, the evolution of the quadratic frequency comb is described by the quadratic-coupled mean-field equation. We further apply supermode theory in conjunction with the Analytical Bloch-Messiah Decomposition (ABMD) to characterize its multimode squeezing structure. Moreover, we identify that the MI gain driven by group-velocity mismatch plays a critical role in mediating the generation of multimode squeezing.

3.1 Bipartite Entanglement Criterion

The dynamics of the bipartite entangled system are detailed in Appendix A. To rigorously quantify the entanglement between the signal and idler modes, we employ the bipartite entanglement criterion from Ref. [39] to evaluate the entanglement measure $C_{s}$. The position ($\hat{x}_{j}$) and momentum ($\hat{y}_{j}$) quadrature operators for each mode are elegantly formulated in terms of the annihilation and creation operators $\hat{a}_{j}$ and $\hat{a}_{j}^{\dagger}$:

[TABLE]

By adjusting the detection angles of the signal and idler beams ($\theta_{s},\theta_{i}$), we derive:

[TABLE]

Introduce the sum and subtraction basis:

[TABLE]

and the fluctuation vector is formulated as:

[TABLE]

The spectral noise density matrix $S_{\hat{X}_{\pm}}(\omega)$ is determined as follows:

[TABLE]

where $T_{r}=G\cdot P=\frac{1}{2}\begin{pmatrix}-\mathrm{i}e^{-\mathrm{i}\theta_{s}}&\mathrm{i}e^{\mathrm{i}\theta_{s}}&-\mathrm{i}e^{-\mathrm{i}\theta_{i}}&\mathrm{i}e^{\mathrm{i}\theta_{i}}\\ e^{-\mathrm{i}\theta_{s}}&e^{\mathrm{i}\theta_{s}}&e^{-\mathrm{i}\theta_{i}}&e^{\mathrm{i}\theta_{i}}\\ -\mathrm{i}e^{-\mathrm{i}\theta_{s}}&\mathrm{i}e^{\mathrm{i}\theta_{s}}&\mathrm{i}e^{-\mathrm{i}\theta_{i}}&-\mathrm{i}e^{\mathrm{i}\theta_{i}}\\ e^{-\mathrm{i}\theta_{s}}&e^{\mathrm{i}\theta_{s}}&-e^{-\mathrm{i}\theta_{i}}&-e^{\mathrm{i}\theta_{i}}\end{pmatrix}.$

The Duan criterion is expressed as follows [40]:

[TABLE]

where $(\Delta\hat{x}_{-})^{2}=S_{\hat{X}_{\pm}}(\omega)(4,4)$, $(\Delta\hat{y}_{+})^{2}=S_{\hat{X}_{\pm}}(\omega)(1,1)$, and $O=\cos(\theta_{s}-\theta_{i})$. If the Duan criterion is violated, i.e., $C_{s}<0,$ the bipartite modes are confirmed to be entangled. A lower value of $C_{s}$ indicates a stronger degree of quantum entanglement.

3.2 Frequency-Dependent Squeezing

Our platform inherently accommodates the entanglement analysis of multiple bipartite mode pairs. However, given the similarity of results across different modes, we focus our analysis on the first mode ($l=1$) for conceptual clarity. Simulations are conducted at 20 °C, with thermal fluctuations disregarded. We compute the steady-state intracavity amplitudes of the pump and signal modes as functions of the injected pump amplitude using Mathematica [Figs. 3(a), 3(b)]. The pump detuning is fixed at $\Delta_{p}=3\times 10^{9}~\mathrm{rad/s}$. Below a critical threshold, the signal mode remains unexcited; above it, a nonzero amplitude emerges, defining the threshold power $P_{\mathrm{th}}$.

Above threshold, the intracavity pump power saturates, while the signal power increases quasi-linearly. A markedly distinct behavior arises at large pump detuning ($\Delta_{p}=3.03\times 10^{11}~\mathrm{rad/s}$), where the subharmonic field acquires positive detuning ($\Delta=4.78\times 10^{8}~\mathrm{rad/s}$). Under these conditions, the signal power curve exhibits pronounced bending above threshold [Fig. 3(d)], while the intracavity pump power falls below the signal power—opposite to the behavior observed at low pump detuning.

Figure 4(a) plots $P_{\mathrm{th}}$ as a function of mode number $l$ under high detuning. At lower $l$, the threshold increases with $l$, indicating that modes further detuned from $\omega_{0}$ are more readily excited—an inversion of the trend seen under weak detuning. This stems from modified phase-matching conditions, each dictating a distinct excitation threshold. Benefiting from lithium niobate’s high nonlinearity and favorable dispersion profile, the threshold power remains in the milliwatt regime—markedly lower than in conventional platforms.

We further implement a pump-tunable entanglement scheme. As shown in Fig. 4(b), increasing the pump amplitude induces a continuous redshift in the frequency maximizing entanglement. At threshold, this frequency approaches zero, accompanied by a substantial bandwidth broadening. The optimal frequency lies in the few-GHz domain, closely mirroring the multimode squeezing scenario discussed below. This spectral tunability highlights the potential for real-time optimization via pump power modulation.

Simulations with fixed $\Delta_{p}=3\times 10^{9}~\mathrm{rad/s}$ and $B_{\mathrm{in}}=1.5\times 10^{8}~\mathrm{V/m}$, and varying coupling strength $r$, reveal analogous behavior [Fig. 4(c)]: progressive redshift of the optimal entanglement frequency, bandwidth broadening, and enhanced squeezing. This arises from overcoupling, which augments the SPDC gain and broadens the phase-matching bandwidth.

To capture spectral dependence, we define the readout angle $\phi=\theta_{s}-\theta_{i}$. As shown in Fig. 4(d), the optimal readout angle varies monotonically with frequency, enabling its continuous adaptation. This frequency-dependent squeezing is highly pertinent for precision metrology, such as gravitational wave detection, offering a broadband, low-loss alternative to narrowband filtering cavities, and thus facilitating more efficient quantum-enhanced measurements.

3.3 Quadratic Coupled Mean-Field Equation

To investigate the frequency comb dynamics dominated by $\chi^{(2)}$ nonlinearity, we now derive the quadratic coupled mean-field equations. Starting from the coupled-mode equations for the signal and idler modes, we extend them to a multimode framework while neglecting quantum fluctuations. We first derive the total Hamiltonian of the multimode system as follows:

[TABLE]

where $g_{0}$ is the $\chi^{(2)}$ nonlinear coefficient, assumed uniform for all modes. This leads to the following coupled-mode equation for an arbitrary mode $u$ of the subharmonic field and an arbitrary mode $v$ of the pump field:

[TABLE]

The range of $u$ is defined as $\left[-\frac{N}{2}+1,\frac{N}{2}\right]$, while the range of $v$ is $\left[-\frac{N}{2}+960,\frac{N}{2}+959\right]$. This choice ensures that the central mode number of the pump field is 959. Here, $N$ represents the total number of modes considered for both the subharmonic and pump fields. The function $\delta[\cdot]$ denotes the Dirac delta function, which takes a value of 1 when its argument is zero and 0 otherwise.

We define $\boldsymbol{\alpha}=\left\{\alpha_{-N/2+1},\dots,\alpha_{N/2}\right\}^{\text{T}}$ and $\boldsymbol{\beta}=\left\{\beta_{-N/2+960},\dots,\beta_{N/2+959}\right\}^{\text{T}}$, allowing the system of $2N$ equations to be rewritten as two equations. Due to the computational complexity of the summation terms, we transform the frequency-domain convolution into a time-domain summation by applying the inverse discrete Fourier transform (IDFT) to the equations. We presume $a=\text{IDFT}[\boldsymbol{\alpha}],\quad b=\text{IDFT}[\boldsymbol{\beta}]$. According to Ref. [24], the nonlinear effect compensates for the spectral phase induced by walk-off, allowing a LN resonator with significant walk-off to still be modeled using the mean-field approximation. Under this approximation, the evolution of the subharmonic and pump fields can be described by two coupled mean-field equations:

[TABLE]

where all parameters are consistent with those of the LN micro-ring resonator employed for QFC generation. Here, $t$ denotes the cavity’s evolution time, while the fast-time variable $\tau$ spans the interval $\left[-\frac{\tau_{s}}{2},\frac{\tau_{s}}{2}\right]$, with $\tau_{s}$ representing the round-trip time. The appearance of the second-order partial derivative with respect to $\tau$ results from applying the IDFT to $D_{\text{int}}$, while the first-order derivative term in $\tau$ arises from the walk-off effect, which is induced by group velocity mismatch. We introduce the entrywise product operator $\star$, which denotes element-wise multiplication between corresponding rows of column vectors [41]. Clearly, Equation (20) admits a trivial all-zero solution for $a$, whereas $b$ does not exhibit such a trivial solution, corresponding to the sub-threshold region in Figure 3. Both trivial and nonzero solutions can exhibit MI gain, which amplifies random fluctuations, leading to the exponential growth of sidebands and ultimately forming a frequency comb [42]. By solving the mean-field equation using the split-step Fourier method, we can effectively simulate the dynamical evolution of the comb teeth.

3.4 Modulation Instability-Induced Frequency Comb

In this subsection, we simulate the evolution of a quadratic OFC. To assist readers in comprehending the underlying mechanism, a schematic diagram illustrating the formation of the quadratic OFC is presented in Figure 5. The pump mode initiates the generation of multiple subharmonic modes via OPO processes. These subharmonic modes subsequently regenerate pump components through second-harmonic generation (SHG) and sum-frequency generation (SFG). Facilitated by cavity-enhanced cascaded nonlinear interactions, this feedback loop drives the exponential proliferation of frequency components, ultimately resulting in the formation of a densely spaced, uniformly distributed quadratic OFC [43]. We set $N=800$, corresponding to the consideration of subharmonic field modes ranging from $-399$ to $400$ and pump field modes spanning from $560$ to $1359$, yielding a total of 1600 modes. Figure 6 presents the simulation outcomes at an evolution time of 700 ns, using a time step $dt=0.1$ ps to ensure the robustness of the split-step Fourier method. To initiate the optical parametric oscillation process, random Gaussian white noise is seeded in each optical mode, with an injected pump amplitude of $B_{\text{in}}=1.9\times 10^{8}$ V/m.

Figure 6(a) delineates the temporal evolution of the aggregate power in both the subharmonic and pump fields. Initially, the subharmonic field exhibits negligible power due to the limited excitation of modes. As the evolution time reaches the marked interval in the figure, optical parametric amplification becomes prominent, leading to the excitation of a multitude of comb modes and a precipitous increase in total power. Simultaneously, the temporal gain-clipping mechanism is activated: due to group velocity mismatch, the signal pulse with the greatest temporal overlap with the pump pulse experiences the maximum gain, resulting in pulse compression. As the pump becomes depleted, gain saturation sets in, and the pulse centroid undergoes a shift driven by nonlinear acceleration [25]. Once spectral broadening and compression, as well as gain and loss, reach dynamic equilibrium—and temporal synchronization between the pump and subharmonic fields is established—the system ultimately settles into steady state.

Figures 6(c) and 6(d) respectively exhibit the time-domain pulse and the OFC for the subharmonic field, while Figures 6(e) and 6(f) illustrate the corresponding time-domain pulse and OFC for the pump field. Under the stipulated simulation conditions, the time-domain profile reveals a bi-periodic pulse arrangement, corresponding to a spectral spacing of 2-FSR. It is noteworthy that, in contrast to soliton crystals, modulation instability combs are acutely sensitive to parameter variations, and the inter-pulse intervals are typically non-uniform. Extended simulation runs have corroborated the robust stability and coherence of the generated frequency comb.

Figure 6(b) depicts the evolution of the subharmonic field intensity $|a(t,\tau)|^{2}$ as a function of the slow evolution time $t$ and the fast time variable $\tau$. The image reveals a series of evenly spaced, uniformly sloped “bright stripes”. To dissect this phenomenon more profoundly, we extracted three temporal cross-sections, as shown in Figure 6(g). Notably, as time progresses, the subharmonic and pump field pulses exhibit a synchronized drift at a constant velocity, a behavior consistent with the observations reported in Ref. [24]. The drift is significantly slower than the walk-off parameter $\Delta k^{\prime}$, indicating that the linear spectral phase induced by the walk-off effect is nonlinearly compensated. Fundamentally, this compensation enables group velocity locking between the subharmonic and pump waves [44, 45]. This drift phenomenon manifests exclusively in the presence of walk-off and can be negated via a suitable transformation to a moving reference frame [46].

Furthermore, we observed that increasing the injected pump power leads to an augmentation in the number of pulses per period. For instance, when $B_{\text{in}}=3\times 10^{8}$ V/m, four pulses per period are observed; further increasing $B_{\text{in}}$ to $4\times 10^{8}$ V/m results in the emergence of eight pulses, as illustrated in Figure 7. Our interpretation is twofold: first, a higher pump power accelerates gain saturation, thereby triggering competition among multiple pulses; second, intense nonlinear phase shifts accumulate within each pulse, inducing self-steepening in the time domain and ultimately causing the pulse to fragment into multiple sub-pulses [47]. Moreover, although a single-pulse state can be accessed by tuning the detuning parameter, it merely corresponds to a saddle point—an unstable critical state—whereas multi-soliton states exhibit greater stability due to their more uniform energy distribution. By sweeping the detuning from $-1\times 10^{10}~\mathrm{rad/s}$ to $2.6\times 10^{10}~\mathrm{rad/s}$, we obtain soliton crystal states, as illustrated in Figures 7(d) and 7(e). These states demonstrate excellent coherence and robustness against perturbations. The corresponding frequency combs are consistent with the nonlinear processes depicted in Figure 5, confirming the underlying formation mechanism.

3.5 Supermode Theory

To unveil the structural characteristics of the collective quantum fluctuations in the quadratic OFC, we extend the quantum fluctuation equations to a multimode framework, yielding:

[TABLE]

To transform the system of $N$ coupled equations into a matrix form, we define the vectors:

[TABLE]

These vector representations enable us to express the quantum fluctuation dynamics in a compact matrix form as follows:

[TABLE]

where $U_{a}^{\mathrm{in}}=\sqrt{2\gamma}E_{{4N}}$, $U_{a}^{\mathrm{loss}}=\sqrt{2\mu}E_{{4N}}$. To simplify the summation terms in the elements of the $M_{a}$ matrix, we introduce the Fourier transform matrix $F_{N}$, defined as $F_{N}=\frac{1}{N}\begin{pmatrix}W^{0}&W^{0}&\cdots&W^{0}\\ W^{0}&W^{1}&\cdots&W^{N-1}\\ \vdots&\vdots&\ddots&\vdots\\ W^{0}&W^{N-1}&\cdots&W^{(N-1)(N-1)}\end{pmatrix}$, where $W$ is given by $W=e^{-\mathrm{i}\frac{2\pi}{N}}.$ By left-multiplying Equation (23) with $F_{N}$ and then subsequently left-multiplying by $F_{N}^{-1}$, we can express $M_{a}$ as Equation (24), where $\chi_{A}=-\mathrm{i}\text{diag}(\Delta_{-N/2+1},\dots,\Delta_{N/2})-\Gamma E_{N},\chi_{B}=-\mathrm{i}\text{diag}(\Delta_{-N/2+960},\dots,\Delta_{N/2+959})-\Gamma E_{N}$. Owing to $\hat{\delta}^{\text{out}}=-\hat{\delta}^{\text{in}}+U_{a}^{\text{in}}\hat{\delta}$, we can obtain $\hat{\delta}^{\text{out}}=M^{\text{in}}(\omega)\hat{\delta}^{\text{in}}+M^{\text{loss}}(\omega)\hat{\delta}^{\text{loss}}$, where $M^{\text{in}}(\omega)=U_{a}^{\text{in}}(\mathrm{i}\omega E_{4N}-M_{a})^{-1}U_{a}^{\text{in}}-E_{4N}$.

To calculate the noise variance spectrum of any local oscillator, we define a supermode decay operator [26]:

[TABLE]

Its Hermitian conjugate operator creates a superposition of photons between the longitudinal modes of the frequency comb, where $\eta_{l}=|\eta_{l}|e^{\mathrm{i}\varphi_{l}}$ defines the amplitude and phase components of the supermode. After performing the transformation as in Equation (12), we can obtain the spectral noise density matrix:

[TABLE]

It is worth noting that, for computational convenience, the inherent cavity loss term has been excluded in Equation (26). This simplification is justified, as cavity losses only affect the extraction of the deformed supermodes, without impacting their composition. By applying ABMD, we can identify the maximum squeezing based on the given spectral noise density matrix. Combining the methods outlined in references [48] and [49], we can implement the ABMD as follows: $S_{x}(\omega)=U(\omega)D(\omega)V^{\dagger}(\omega)$. Further details can be found in Appendix C.

In order to express the squeezing characteristics of cavity-deformed supermodes, we define the corresponding Hermitian quadrature operators as follows: $\hat{L}^{(+)}(t)=\frac{1}{\sqrt{2}}\left(\hat{L}^{\dagger}(t)+\hat{L}(t)\right),\hat{L}^{(-)}(t)=\frac{\mathrm{i}}{\sqrt{2}}\left(\hat{L}^{\dagger}(t)-\hat{L}(t)\right).$ Thus, the spectral noise of the quantum fluctuations can be expressed using the Wiener-Khinchin theorem [48]:

[TABLE]

Here, $V^{+}(\omega)$ and $V^{-}(\omega)$ represent the squeezing or anti-squeezing of the cavity-deformed supermodes in the absence of intrinsic cavity loss. The composition of the maximum squeezing cavity-deformed supermode corresponding to $U(\omega)$ can be derived from the columns of $U(\omega)$. When $U(\omega)$ is a real matrix, the components of the $k$-th maximum squeezing cavity-deformed supermode satisfy:

[TABLE]

We can now elucidate the collective quantum fluctuation behavior in the context of cavity-deformed supermodes and their squeezing dynamics.

3.6 Analysis of Multimode Squeezing

By combining the quadratic optical frequency comb generated in Section 3.4 with the supermode formalism developed in Section 3.5, we employ the ABMD to analyze the system’s multimode squeezing characteristics. To elucidate the underlying physical mechanism responsible for the emergence of multimode squeezing during the evolution of the quadratic frequency comb, we begin with the theoretical framework of MI gain. Specifically, MI yields both a trivial gain associated with the zero solution—corresponding to the unexcited parametric field—which is independent of group-velocity mismatch and determined solely by dispersion and pump strength, and a nontrivial MI gain driven by walk-off. The latter plays a critical role in the formation of the frequency comb. Details on MI theory and the gain of the trivial solution can be found in Appendix B. Figure 8(a) shows the evolution of nontrivial MI gain as a function of cavity detuning, which directly corresponds to the distribution of the morphing supermode coefficients shown in Fig. 8(b), satisfying the normalization condition $\sum_{l}|\eta_{l}|^{2}=1.$ In this process, periodic modes far from the critical frequency become unstable, leading to their disappearance and the emergence of new modes closer to the MI gain maximum. As the peak MI gain frequency evolves with detuning $\Delta$, previously stable sidebands eventually become sufficiently detuned from the critical frequency to become unstable, giving rise to new modes that align more closely with the evolving MI gain peak. This indicates that MI not only drives the formation of the quadratic frequency comb but also initiates the excitation of its corresponding squeezed modes.

Figure 8(c) displays the spectral profiles and supermode decomposition coefficients for both the Turing pattern and soliton crystal states. In both cases, strong multimode squeezing is observed, with squeezing levels reaching 21.01 dB and 37.06 dB, respectively. The supermode coefficients are predominantly concentrated around the subharmonic field modes, corresponding to their higher power occupancy. After the formation of the soliton crystal, the maximum of the supermode decomposition shifts to the central subharmonic mode, reflecting enhanced localization of squeezing in that mode.

To underscore the collective nature of multimode squeezing—where all modes within a spectral band must be incorporated to avoid the loss of quantum correlations due to inter-band coupling—we conduct a separate analysis of the subharmonic and pump fields. As illustrated in Fig. 8(d), the morphing supermode coefficients of the pump field are entirely concentrated in a single mode ($l=971$), owing to minimal detuning and favorable phase-matching conditions. Interestingly, the subharmonic field alone sustains significant squeezing, in some cases even exceeding the global squeezing level (see Fig. C.1). In stark contrast, the pump field exhibits negligible squeezing when considered in isolation, rendering the multimode squeezing framework ineffective. This divergence stems from the fact that, even without subsequent SHG or SFG, the initial OPO process driven by the pump field can establish multimode squeezing in the subharmonic field. However, neglecting the modes involved in SHG and SFG disrupts the quantum nonlinear couplings essential for sustaining multimode entanglement. As a predominantly classical drive, the pump field itself lacks the capacity to generate quantum squeezing, and thus cannot form a coherent multimode squeezed system on its own.

Figure 9(a) shows the evolution of multimode squeezing at zero Fourier frequency throughout the full detuning scan. Since the linearized fluctuation model is only valid on a stable mean-field background, our analysis focuses on the Turing pattern and soliton crystal regimes. In the Turing region, the degree of squeezing exhibits persistent oscillations, consistent with the periodic intracavity power variations during detuning. In contrast, in the soliton crystal regime, the squeezing level tracks the shift of the most strongly squeezed frequency with detuning. For large detuning, multimode squeezing vanishes entirely, indicating the breakdown of the soliton crystal structure. Using the ABMD, we extract the spectral distribution of multimode squeezing, as shown in Fig. 9(b). The spectrum typically features two dominant peaks whose positions are highly sensitive to detuning. As detuning increases, only a single isolated squeezing peak remains near zero Fourier frequency. This behavior stands in stark contrast to pump-power tuning, where the squeezing peak positions remain nearly unchanged.

Multimode squeezed states, enabled by the ultrabroadband quantum noise suppression inherent to LN, significantly enhance the scalability and practicality of quantum technologies. These states allow for the generation of multimode squeezed light fields, enabling parallel access to quantum resources and supporting the advancement of distributed quantum computing. They also facilitate efficient implementation of multi-node quantum key distribution and quantum repeater protocols. In multiparameter quantum metrology, multimode squeezing offers substantial improvements beyond classical and standard quantum noise limits, with applications ranging from optical clock synchronization to biological imaging. Prospectively, as photonic chip integration continues to mature, this technology is poised to drive large-scale deployment in quantum networks and deep-space communication.

4 Conclusion

We demonstrate a lithium niobate micro-ring platform that leverages $\chi^{(2)}$ nonlinearity to generate entangled QFCs with milliwatt-level thresholds. Through dispersion and coupling engineering, we achieve tunable frequency dependent squeezing and quantify entanglement across modes. Leveraging the quadratic coupled mean-field equation, we conduct an in-depth exploration of the evolutionary dynamics of multimode frequency combs. A supermode analysis reveals ultra-broadband quadrature squeezing, essential for scalable quantum networks and parallel quantum processing. We further propose and elucidate a physical mechanism wherein MI, beyond triggering the emergence of quadratic frequency combs, also drives multimode quadrature squeezing across the relevant modes.This work establishes LN as a superior material for integrated quantum photonics, bridging nonlinear dynamics and quantum optics for chip-scale applications.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 62075129), Microwave Photonics Technology Key Laboratory of Sichuan Province (Grant No.2023-04) and the Science and Technology on Metrology and Calibration Laboratory (Grant No. JLKG2024001B002).

Appendix A Bipartite Entanglement Dynamics

This section delves into the quantum dynamics of the resonator and presents a detailed analysis of two-mode quadrature squeezing.

A.1 System Hamiltonian

The system Hamiltonian rigorously characterizes the optical nonlinear dynamics in the quantum regime, with the resonant modes described by the annihilation operators $\hat{a}_{s}$, $\hat{a}_{i}$, and $\hat{b}$. Here, $\hat{a}_{s}$ and $\hat{a}_{i}$ correspond to the signal and idler modes in the 1560 nm band, while $\hat{b}$ denotes the pump mode in the 780 nm band.

The free Hamiltonian $\hat{H}_{0}$ is given by:

[TABLE]

For SPDC, the nonlinear Hamiltonian is formulated as:

[TABLE]

Here, $g_{0}$ represents the second-order nonlinear coefficient of LN and can be expressed as [50]:

[TABLE]

where $\chi^{(2)}$ denotes the second-order nonlinear susceptibility tensor of LN, $\epsilon_{0}$ and $\epsilon_{1}$ represent the vacuum permittivity and relative permittivity, respectively, $A_{\text{eff}}$ is the effective mode area, and $R$ denotes the radius of the microcavity. For simplicity, higher-order nonlinear effects beyond second-order nonlinearity are neglected.

The total Hamiltonian can be expressed as:

[TABLE]

A.2 Heisenberg-Langevin Equations

The Heisenberg-Langevin equation [51] elegantly unifies the Heisenberg equation of motion with Langevin noise terms, providing a comprehensive description of the dynamical evolution of an open quantum system influenced by quantum fluctuations. This framework can be succinctly formulated as follows, with the signal mode serving as an illustrative example:

[TABLE]

Here, the considered modes are presumed to exhibit analogous field distributions and are governed by a unified total loss rate $\Gamma$, which encapsulates the intrinsic dissipation rate $\mu$ and the coupling rate $\gamma$, satisfying the relation $\Gamma=\gamma+\mu$. The annihilation operators $\hat{a}^{\text{in}}$ and $\hat{a}^{\text{loss}}$ respectively represent the input and loss modes of the resonator. Within this theoretical construct, the loss mode is postulated to reside in the vacuum state, while the incident signal and idler modes are likewise assumed to be in vacuum states. The expectation value of the input pump mode is articulated as:

[TABLE]

where $P_{\mathrm{in}}$ (Watt) is the pump laser power.

By invoking the rotating wave approximation, wherein $\hat{a}_{j}e^{-i\Omega_{j}t}$ replaces $\hat{a}_{j}$ $(j=s,i)$ and $\hat{b}e^{-i\Omega_{p}t}$ replaces $\hat{b}$, the Heisenberg-Langevin equations governing the dynamical evolution of the pump, signal, and idler modes are elegantly expressed as:

[TABLE]

where $\Delta_{s}$ and $\Delta_{i}$ represent the normalized cold cavity detunings corresponding to the signal and idler modes, respectively. Since $\Omega_{p}=\Omega_{s}+\Omega_{i}$, the phase terms in the rotating wave approximation mutually cancel, simplifying the equations.

A.3 Steady-State Equations

We implement a linearization procedure by decomposing each field operator $\hat{a}_{j}$ into its steady-state mean value $\alpha_{j}$ and a fluctuation term $\delta\hat{a}_{j}$, expressed as $\hat{a}_{j}=\alpha_{j}+\delta\hat{a}_{j}$ (The pump mode is set as $\hat{b}=\beta+\delta\hat{b}$).

Under steady-state conditions, $\alpha_{j}$ remains constant, allowing us to impose $\delta\hat{a}_{j}=0$ and $\frac{d\alpha_{j}}{dt}=0,$ thereby deriving the steady-state Heisenberg-Langevin equations. Within this framework, the input fields of the signal and idler modes, along with the loss modes of the pump, signal, and idler, are all assumed to be in the vacuum state, satisfying:

[TABLE]

Without loss of generality, we take the phase of the external pump as the reference, yielding the following definitions: $\alpha_{j}=A_{j}\mathrm{e}^{\mathrm{i}\theta_{j}}$, $\beta=B\mathrm{e}^{\mathrm{i}\theta_{\mathrm{p}}}$, $\beta^{\mathrm{in}}=B_{\mathrm{in}}\mathrm{e}^{\mathrm{i}\theta_{\mathrm{in}}}$. For simplicity, we set $A_{s}=A_{i}=A$ and $\Delta_{s}=\Delta_{i}=\Delta$. Utilizing these parameters, we obtain the following expressions:

[TABLE]

By solving Equation (39), we derive the steady-state relationship among the input pump $B_{\text{in}}$, $B$, and $A$.

A.4 Quantum Fluctuations

To rigorously examine the quantum properties of the signal and idler modes, we derive the quantum fluctuation equations from Eq. (35). Given that the steady-state solutions have been systematically determined, these fluctuation equations are obtained by subtracting the steady-state components from Eq. (35). In this framework, the pump field is treated as a classical parameter, thereby neglecting fluctuations in the pump mode ($\delta\hat{b}=0$), while higher-order fluctuation terms are disregarded for analytical tractability.

We formulate the fluctuation vector for the signal and idler modes as follows:

[TABLE]

where $\theta_{j}$ represents the phase of the steady-state mean value, expressed as $\alpha_{j}=A_{j}e^{i\theta_{j}}.$ The time evolution of the fluctuations $\delta\hat{a}_{j}$ is dictated by the following set of linearized equations:

[TABLE]

where $U_{a}^{\mathrm{in}}=\mathrm{diag}\left(\sqrt{2\gamma},\sqrt{2\gamma},\sqrt{2\gamma},\sqrt{2\gamma}\right)$, $U_{a}^{\mathrm{loss}}=\mathrm{diag}\left(\sqrt{2\mu},\sqrt{2\mu},\sqrt{2\mu},\sqrt{2\mu}\right)$. The matrix $M_{a}$ arises from the linearization process, with its elements dictated by the mean field amplitudes and detuning parameters. Below the threshold (i.e., $A=0$), $M_{a}$ can be expressed as:

[TABLE]

The frequency-domain dynamics of these fluctuations are obtained via a Fourier transform. Incorporating this transformation with the resonator’s input-output relations, $\hat{a}^{\text{out}}=-\hat{a}^{\text{in}}+\sqrt{2\gamma}\hat{a},$ provides a comprehensive characterization of the system’s response:

[TABLE]

where $E$ is the identity matrix.

Thus, the output spectral noise density matrix is formulated as:

[TABLE]

where the matrix M_{c}=\left(\begin{array}[]{cccc}0&1&0&0\\ 0&0&0&0\\ 0&0&0&1\\ 0&0&0&0\end{array}\right).

A.5 Two-Mode Quadrature Squeezing

The Duan criterion [39, 40] enables a precise quantitative characterization of pairwise entanglement and quadrature squeezing in the sub-threshold regime. Figure A.1(a) illustrates the process of forming multiple pairs of entangled states via spontaneous parametric down-conversion (SPDC) below the oscillation threshold in a LN microring resonator OPO. As shown in Figure A.1(b), when measuring the idle frequency mode at a specific angle $\theta_{i}$, the signal mode exhibits a squeezing effect in the $-\theta_{i}$ direction, with statistical properties superior to the vacuum state, thereby surpassing the standard quantum limit. In contrast, the signal mode demonstrates anti-squeezing along the orthogonal direction, $\pi/2-\theta_{i}$. This phenomenon is further validated in Figure A.2, where we observe that for a given $\theta_{i}$, the maximum value of positive entanglement (corresponding to anti-squeezing) occurs at $\theta_{s}=\pi/2-\theta_{i}$ or $\theta_{s}=3\pi/2-\theta_{i}$.

Appendix B Modulation Instability Gain

The dynamics of the OPO can be described by an infinite-dimensional map of the field amplitudes, which, as shown in Ref. [42, 44], can be reduced to a single mean-field equation governing the subharmonic field $a$:

[TABLE]

Here, we assume the phase-matching condition is satisfied and $\otimes$ denotes the convolution operation. The temporal response function is given by $I(\tau)=\mathrm{IDFT}[\hat{I}(\Omega)]$, where

[TABLE]

The frequency $\Omega$ denotes the offset angular frequencyrelative to $\omega_{0}$.

It is straightforward to verify that the mean-field equation admits a time-independent constant solution of the form $a_{0}(t,\tau)=|a_{0}|e^{\mathrm{i}\phi}$. When $|a_{0}|=0$, this corresponds to the trivial steady-state solution of the below-threshold OPO. Substituting this ansatz into the mean-field equation, we obtain the following condition:

[TABLE]

Solving the above equation yields the steady-state intensity of the subharmonic field:

[TABLE]

To analyze the stability of the constant solution against the growth of new frequency components, we perform a modulation instability (MI) analysis. Specifically, we consider the following ansatz for the subharmonic field:

[TABLE]

where $A_{1}$ and $A_{2}$ represent small perturbations around the steady-state amplitude $a_{0}$. Substituting this ansatz into the mean-field equation and retaining only the first-order terms in $A_{1}$ and $A_{2}$, we project onto the corresponding frequency components. This yields the following set of linearized differential equations governing the evolution of the three field amplitudes:

[TABLE]

The coupled evolution of the perturbations can be recast into a compact matrix form as:

[TABLE]

For the nontrivial steady-state solution, the eigenvalues of the matrix $\mathbf{M}$ are given by:

[TABLE]

where the field amplitude $|a_{0}|$ is related to the input field $B_{\mathrm{in}}$ via Eq. (48), and $\iota_{\pm}(\Omega)=\hat{I}(\Omega)\pm\hat{I}^{*}(-\Omega).$

The nontrivial constant solution loses stability when the MI gain satisfies $\mathrm{Re}[\lambda_{+}]>0$.

Similarly, the trivial zero solution demonstrates MI when the eigenvalues have a positive real part.

[TABLE]

The instability of the trivial solution is evidently not governed by the detuning parameter but is instead induced by group velocity dispersion. Figure B.1 illustrates the MI gain associated with the trivial solution.

Appendix C Analytical Bloch-Messiah Decomposition

By integrating the methodologies proposed in Ref. [48] and Ref. [49], the Analytical Bloch-Messiah Decomposition (ABMD) of $S_{x}(\omega)$ (the propagation matrix after similarity transformation) can be carried out through a systematic analytical procedure centered at $\omega=0$, and extended across the frequency domain. The process begins with the polar decomposition of the zero-frequency transformation

[TABLE]

where $P$ is a real, symmetric, positive-definite symplectic matrix, and $Y$ is a real orthogonal matrix. This decomposition serves as the foundation for identifying the canonical modes.

The matrix $P$ is then partitioned into block form,

[TABLE]

which allows the construction of a complex symmetric matrix

[TABLE]

The Takagi (or Autonne) decomposition is then performed on $M$, yielding

[TABLE]

where $\Lambda$ is diagonal with non-negative entries and $W$ is a unitary matrix. This decomposition can be explicitly computed via the singular value decomposition

[TABLE]

followed by constructing the unitary matrix

[TABLE]

With $W$ obtained, we build the orthogonal symplectic matrix $U$ as

[TABLE]

and defines the diagonal matrix $\Xi=\Lambda+\sqrt{\Lambda\Lambda+E}$, from which the full diagonal scaling matrix $D=\begin{pmatrix}\Xi&0\\ 0&\Xi^{-1}\end{pmatrix}$is constructed. The final component of the canonical form is given by

[TABLE]

completing the decomposition $S_{x}(0)=U(0)D(0)V^{\dagger}(0)$ at zero frequency.

To extend this decomposition perturbatively around $\omega=0$, the derivative of the symplectic transformation is computed as

[TABLE]

and used to define the auxiliary matrix

[TABLE]

This matrix is then partitioned into four submatrices:

[TABLE]

From these submatrices, the Hermitian and anti-Hermitian generator matrices $H$ and $K$ are constructed. The elements of $H_{1}$ and $K_{1}$ satisfies

[TABLE]

To compute $H_{2}$ and $K_{2}$, we solve the coupled equations:

[TABLE]

which yield the remaining components of the generator matrices. Combining these, we construct

[TABLE]

The variation of the diagonal matrix $D$ with frequency is determined from the diagonal elements of $Q_{1}$, using the expression

[TABLE]

which enables extrapolation of $D(\omega)$ in the vicinity of $\omega=0$.

To update the transformation matrices for nearby frequencies, a perturbative Magnus expansion is applied. The matrices $U(\omega+\delta\omega)$ and $V(\omega+\delta\omega)$ are approximated as

[TABLE]

which allows the decomposition $S_{x}(\omega)=U(\omega)D(\omega)V^{\dagger}(\omega)$ to be analytically propagated across the frequency domain through iterative application of the differential steps.

According to the ABMD, optimal squeezed states can be achieved by appropriately tailoring either the output or input cavity modes. Figure C.1 presents a spectral comparison between the multimode squeezing levels obtained under different conditions and those derived by considering only the subharmonic field. Remarkably, the multimode squeezing achieved by considering solely the subharmonic modes can be substantial—occasionally even surpassing the squeezing obtained when all modes are included. This counterintuitive result arises because incorporating additional, non-squeezed modes may introduce greater squeezing degradation.