¹

²

1 Introduction↩︎

Photon correlations are central to quantum optics and underpin a wide range of quantum technologies. They are commonly characterized by the second-order correlation function \(g^{(2)}(\tau)\) introduced by Glauber [1]: \(g^{(2)}(0)=1\) for coherent light, \(g^{(2)}(0)=2\) for thermal or chaotic light (which exhibits greater bunching than coherent light) [2], and \(g^{(2)}(0)<1\) for antibunched (e.g., photon-blockaded) light. Classical setups, such as cascaded rotating ground glasses, can convert laser light into pseudothermal light with altered intensity fluctuation distributions achieving \(g^{(2)}(0) \sim 3\) [3]. The superbunching regime, \(g^{(2)}(0)\gg 2\), is of special interest: photons with such extreme correlations could be useful in applications such as ghost imaging, nonlinear optics, and electron sources [4]–[6].

Superbunching can occur through various quantum phenomena, including spontaneous parametric down-conversion [7]–[11], collective effects like superradiance [12]–[15], and intra-emitter quantum effects [16], [17]. These have been shown to achieve \(g^{(2)}(0) \sim 10^2\), though recently a specially-prepared nonlinear photonic crystal fiber has been claimed to reach \(g^{(2)}(0) \sim 10^4\) [18]. Broadly speaking, known methods rely on either strong optical nonlinearities to generate a bright squeezed vacuum [8], or the presence of quantum matter (e.g., collectively coherent emitters or quantum dots with cascaded emission pathways). This poses a challenge for using integrated quantum photonic devices to generate superbunched light, because such devices have small form factors and, preferably, utilize conventional optical media featuring only weak intrinsic nonlinearities.

However, superbunching by vacuum squeezing may not strictly require strong nonlinearity. An analogy may be drawn to the phenomenon of unconventional photon blockade (UPB) [19]–[27], whereby photon blockade (\(g^{(2)}(0) \rightarrow 0\)) is achieved in a system of coupled optical modes by carefully designing the mode excitation pathways and how they interfere with one another. UPB has been demonstrated experimentally [28], [29], and there have been recent theoretical proposals aimed at drastically reducing the nonlinearity strengths required [26], as well as increasing the antibunching lifetime (i.e., the range of \(\tau\) over which \(g^{(2)}(\tau)\) is suppressed) [27]. Could this interference-engineering approach be applied to the opposite purpose—superbunching as opposed to antibunching? If so, what is the largest \(g^{(2)}(0)\) that can be achieved, and what sort of photonic system is required?

In this paper, we predict that a system of coupled photonic modes with realistically weak Kerr nonlinearities can exhibit giant superbunching with \(g^{(2)}(0) > 100\), comparable to previous superbunching results based on strong nonlinearities, and possibly even higher. This is accomplished by engineering destructive interference among coherent driving paths to cancel the classical amplitude \(\langle \hat{a}_i \rangle\) in a target mode \(i\), leaving the photon emission governed entirely by the quantum fluctuations \(\hat{d}_i \equiv \hat{a}_i-\langle \hat{a}_i \rangle=\hat{a}_i\). Consequently, the pair correlation \(\langle \hat{d}_i\, \hat{d}_i \rangle\) dominates over the normal occupation \(\langle \hat{d}_i^\dagger \hat{d}_i \rangle\), giving rise to \[g_{ii}^{(2)}(0) - 2 \;\propto\; \frac{1}{\langle \hat{a}_i^\dagger \hat{a}_i\rangle} \gg 1. \label{eq:scaling95result}\tag{1}\] We develop an analytical theory for this mechanism, and verify it with numerical master equation simulations on a minimal four-mode model. This superbunching scheme is highly compatible with scalable integrated photonic platforms such as silicon carbide microring resonators [30], [31], as it does not rely on giant optical nonlinearities or collective emission, nor exotic (e.g., non-passive) couplings between the participating optical modes. Moreover, the required interferometric conditions can be achieved via standard active control parameters in photonic systems, such as laser frequency detuning and inter-cavity coupling strength. Such interference-aided superbunching may thus supply a workable design principle for compact quantum light sources.

2 Coupled-mode model↩︎

We begin by presenting our theoretical framework for studying superbunching. Consider a system of \(N\) coupled photonic modes, each of which may be subject to on-site Kerr nonlinearity, loss, and coherent driving. In the frame co-rotating with the driving frequency, the Hamiltonian is (setting \(\hbar = 1\)) \[\begin{align} \mathcal{H} = \sum_{j} &\left[-\Delta_j \hat{a}_j^\dagger \hat{a}_j + U_j \hat{a}_j^\dagger \hat{a}_j^\dagger \hat{a}_j \hat{a}_j + F_j \hat{a}_j^\dagger + F_j^* \hat{a}_j\right] \nonumber \\ &+ \sum_{j \neq k} J_{jk} \hat{a}_j^\dagger \hat{a}_k, \label{eq:H} \end{align}\tag{2}\] where, for each mode \(j\), \(\Delta_j\) is the frequency detuning (relative to the drive), \(U_j\) is the Kerr coefficient, \(F_j\) is the drive amplitude, \(J_{jk} = J_{kj} \in \mathbb{R}\) is the inter-mode coupling (in the absence of magneto-optic effects), and \(\hat{a}_j, \hat{a}_j^\dagger\) are the photon annihilation and creation operators. We also assign each mode a loss rate \(\gamma_j\), which enters via the Lindblad master equation \[\begin{align} \frac{d\hat{\rho}}{dt} &= -i[\hat{\mathcal{H}}, \hat{\rho}] + \sum_j \gamma_j \mathcal{D}[\hat{a}_j]\hat{\rho}, \tag{3} \\ \mathcal{D}[\hat{o}]\hat{\rho} &\equiv \hat{o}\hat{\rho} \hat{o}^\dagger - \{\hat{o}^\dagger \hat{o}, \hat{\rho}\}/2. \tag{4} \end{align}\]

The steady-state second-order correlation between sites \(i\) and \(j\), with time delay \(\tau\ge0\), is \[g^{(2)}_{ij}(\tau)=\lim\limits_{t\rightarrow + \infty} \frac{\langle\hat{a}_j^\dagger(t)\hat{a}_i^\dagger(t+\tau)\hat{a}_i(t+\tau)\hat{a}_j(t)\rangle}{\langle\hat{a}_j^\dagger(t)\hat{a}_j(t)\rangle\langle\hat{a}_i^\dagger(t+\tau)\hat{a}_i(t+\tau)\rangle}, \label{g2def}\tag{5}\] with the expectation values taken in the steady state. We will mainly be interested in the on-site correlation \(g^{(2)}_{ii}(0)\) and the mean photon number \(\langle n_i \rangle = \langle \hat{a}_i^\dagger \hat{a}_i \rangle\).

To separate out the classical response, let \[\hat{a}_j = \alpha_j + \hat{d}_j, \label{eq:displace}\tag{6}\] where \(\alpha_j \in \mathbb{C}\) is the classical steady-state amplitude and \(\hat{d}_j\) is a fluctuation operator satisfying \(\langle \hat{d}_j \rangle = 0\). If we take the Heisenberg equation \(i\dot{\hat{a}}_j = [\hat{a}_j, \mathcal{H}] - i(\gamma_j/2)\hat{a}_j+\hat{\xi}_j\), where \(\hat{\xi}_j\) is the quantum noise operator, applying the steady-state condition \(\dot{\alpha}_j = 0\) yields \[\left(-\Delta_j - i\frac{\gamma_j}{2}\right) \alpha_j + 2U_j |\alpha_j|^2 \alpha_j + \sum_k J_{jk} \alpha_k + F_j = 0. \label{eq:classical}\tag{7}\] This is a set of \(N\) coupled cubic equations whose solutions determine the classical mean field.

As shown in Appendix 6, after extracting this classical contribution, the remaining terms form an effective Hamiltonian \[\begin{align} \hat{\mathcal{H}}_\mathrm{eff} = \sum_j \bigg[&\tilde{\omega}_j \hat{d}_j^\dagger \hat{d}_j + U_j\!\left(\alpha_j^2 (\hat{d}_j^\dagger)^2 + (\alpha_j^*)^2 \hat{d}_j^2\right) \nonumber \\ &+ 2U_j\!\left(\alpha_j^* \hat{d}_j^\dagger \hat{d}_j^2 + \alpha_j (\hat{d}_j^\dagger)^2 \hat{d}_j\right) \nonumber \\ &+ U_j\, \hat{d}_j^\dagger \hat{d}_j^\dagger \hat{d}_j \hat{d}_j\bigg] + \sum_{j \neq k} J_{jk}\, \hat{d}_j^\dagger \hat{d}_k, \label{eq:Heff} \end{align}\tag{8}\] where \[\tilde{\omega}_j = -\Delta_j + 4U_j|\alpha_j|^2 \label{eq:omegatilde}\tag{9}\] is a renormalized on-site frequency detuning. This effective Hamiltonian enters via the master equation \[\frac{d\hat{\rho}}{dt} = -i[\hat{\mathcal{H}}_\mathrm{eff}, \hat{\rho}] + \sum_j \gamma_j \mathcal{D}[\hat{d}_j]\hat{\rho}. \label{eq:S95master}\tag{10}\] The Lindblad dissipator has the same form as in Eq. 4 , since the displacement 6 is a unitary operation on the bosonic Hilbert space.

To solve this quantum Langevin equation analytically, we can linearize \(\hat{\mathcal{H}}_\mathrm{eff}\) by retaining only terms quadratic in \(\hat{d}_j\) (in subsequent sections, we will discuss numerical simulations that include the cubic and quartic contributions). The problem then reduces to \[i\partial_t \boldsymbol{d} = \left(\boldsymbol{H} - i\frac{\boldsymbol{\Gamma}}{2}\right)\boldsymbol{d} + \boldsymbol{U}\,\boldsymbol{d}^{\dagger} - i\sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}, \label{eq:langevin}\tag{11}\] where \(\boldsymbol{d} = (\hat{d}_1, \ldots, \hat{d}_N)^T\) is the vector of fluctuation operators, \(\boldsymbol{\eta} = (\hat{\eta}_1, \ldots, \hat{\eta}_N)^T\) is the vacuum input noise satisfying \(\langle \hat{\eta}_j(t)\,\hat{\eta}_k^\dagger(t') \rangle = \delta_{jk}\,\delta(t - t')\), and \[\begin{align} \boldsymbol{H} &= \mathrm{diag}(\tilde{\omega}_1,\ldots,\tilde{\omega}_N) + [J_{jk}], \\ \boldsymbol{\Gamma} &= \mathrm{diag}(\gamma_1,\ldots,\gamma_N), \\ \boldsymbol{U} &= \mathrm{diag}(2U_1\alpha_1^2,\ldots,2U_N\alpha_N^2). \label{eq:matrices} \end{align}\tag{12}\] Here \(\boldsymbol{H}\) is the single-particle Hamiltonian matrix (encoding renormalized detunings and couplings), \(\boldsymbol{\Gamma}\) is the loss matrix, and \(\boldsymbol{U}\) is the anomalous (squeezing) coupling arising from the Kerr nonlinearity interacting with the classical field.

Equation 11 can be solved by the method of Green’s functions, as shown in Appendix 7. In terms of the retarded Green’s function \(\boldsymbol{G}_\omega = (\boldsymbol{H} - i\boldsymbol{\Gamma}/2 - \omega)^{-1}\), the steady-state two-time correlation functions are \[\begin{align} \langle \boldsymbol{d}(t)\,\boldsymbol{d}^{\dagger T}\!(0) \rangle &= \mathcal{F}_{t}^{-1}\!\left[ \boldsymbol{G}_\omega\, \boldsymbol{\Gamma}\, \boldsymbol{G}_\omega^{\dagger} \right], \tag{13} \\ \langle \boldsymbol{d}^{\dagger}(0)\,\boldsymbol{d}^{T}\!(t) \rangle &= \mathcal{F}_{t}^{-1}\!\left[ \boldsymbol{G}_\omega^{\dagger} \boldsymbol{U}^{\dagger} \boldsymbol{G}_{-\omega}\, \boldsymbol{\Gamma}\, \boldsymbol{G}_{-\omega}^{\dagger} \boldsymbol{U}\, \boldsymbol{G}_\omega \right], \tag{14} \\ \langle \boldsymbol{d}(t)\,\boldsymbol{d}^{T}\!(0) \rangle &= -\mathcal{F}_{t}^{-1}\!\left[ \boldsymbol{G}_\omega\, \boldsymbol{\Gamma}\, \boldsymbol{G}_\omega^{\dagger} \boldsymbol{U}\, \boldsymbol{G}_{-\omega} \right], \tag{15} \\ \langle \boldsymbol{d}^{\dagger}(0)\,\boldsymbol{d}^{\dagger T}\!(t) \rangle &= -\mathcal{F}_{t}^{-1}\!\left[ \boldsymbol{G}_\omega^{\dagger} \boldsymbol{U}^{\dagger} \boldsymbol{G}_{-\omega}\, \boldsymbol{\Gamma}\, \boldsymbol{G}_{-\omega}^{\dagger} \right], \tag{16} \end{align}\] where \(\mathcal{F}_{t}^{-1}[f(\omega)] = \int \frac{d\omega}{2\pi}\, e^{-i\omega t}\, f(\omega)\) is the usual inverse Fourier transform.

3 Superbunching point↩︎

We are interested in the \(g^{(2)}\) function given in Eq. 5 . To simplify the four-operator correlation function in the numerator, we use Wick’s theorem (which is valid for Gaussian states in the linearized theory). At \(\tau = 0\), \[\begin{align} \big\langle \hat{d}_i^\dagger \, \hat{d}_i^\dagger\, \hat{d}_i\, \hat{d}_i \big\rangle = \big|\big\langle \hat{d}_i\, \hat{d}_i \big\rangle\big|^2 + \big|\big\langle \hat{d}_i^\dagger\, \hat{d}_i \big\rangle\big|^2 + \big\langle \hat{d}_i^\dagger\, \hat{d}_i \big\rangle\, \big\langle \hat{d}_i^\dagger\, \hat{d}_i \big\rangle, \end{align}\] and so \[g_{ii}^{(2)}(0) = 2 + \left|\frac{\langle \hat{d}_i\, \hat{d}_i \rangle}{\langle \hat{d}_i^\dagger \hat{d}_i \rangle}\right|^2. \label{eq:S95g295opt}\tag{17}\] In the second term, the numerator is given by Eq. 15 : \[\langle \hat{d}_i\, \hat{d}_i \rangle = -\left[\mathcal{F}_0^{-1}\!\left(\boldsymbol{G}_\omega\,\boldsymbol{\Gamma}\,\boldsymbol{G}_\omega^\dagger\,\boldsymbol{U}\,\boldsymbol{G}_{-\omega}\right)\right]_{ii}, \label{eq:S95dd}\tag{18}\] which is of order \(U\alpha^2\).

As for the numerator, let us consider the special case where the system is tuned to a point where the classical amplitude vanishes at a site \(i\): \(\alpha_i = 0\). The mean photon number is then determined entirely by the zero-point fluctuations: \[\begin{align} \langle n_i \rangle &= \langle \hat{d}_i^\dagger \hat{d}_i \rangle \\ &= \left[\mathcal{F}_0^{-1}\!\left(\boldsymbol{G}_\omega^\dagger\,\boldsymbol{U}^\dagger\,\boldsymbol{G}_{-\omega}\,\boldsymbol{\Gamma}\,\boldsymbol{G}_{-\omega}^\dagger\,\boldsymbol{U}\,\boldsymbol{G}_\omega\right)\right]_{ii}, \label{eq:S95n95opt} \end{align}\tag{19}\] where on the second step we have used Eq. 14 . This is of order \(U^2\alpha^4\). Combining this with Eqs. 17 –15 gives \[g_{ii}^{(2)}(0) - 2 \;\propto\; \frac{1}{U^{2}\alpha^{4}} \;\propto\; \frac{1}{\langle n_i\rangle}. \label{eq:S95scaling}\tag{20}\]

We have thus established the phenomenon of superbunching in this weakly nonlinear system, along with an associated trade-off: it occurs around a point where the mean photon flux vanishes. Intuitively, it is precisely when the classical field vanishes that emission is dominated by correlated photons from parametric scattering. The practical significance of this phenomenon, therefore, depends on the photon detection limits and the maximum value of \(g^{(2)}\) within those limits.

Let us now look for a configuration of optical modes that can host such a superbunching point. Returning to Eq. 7 , we seek a solution with \(\alpha_i = 0\) (for some \(i\)) in the weak-nonlinearity limit \(U_j|\alpha_j|^2 \ll |\Delta_j|, \gamma_j\). By dropping the \(U\) terms, we bring the equation into the linear form \[\begin{align} \boldsymbol{L} &\begin{pmatrix}\alpha_1 \\ \vdots\\\alpha_N \end{pmatrix} + \begin{pmatrix}F_1\\\vdots\\F_N \end{pmatrix} = 0, \label{eq:zerocond} \\ L_{jk} &= z_j\,\delta_{jk} + J_{jk}, \\ z_j &\equiv -\Delta_j - i\gamma_j/2. \end{align}\tag{21}\] When only one mode \(d\) is driven (\(F_j = F\,\delta_{jd}\)), the zero-displacement condition on a signal site \(i \neq d\) reduces to \[L^{-1}_{id} = 0. \label{eq:Linvcond}\tag{22}\] By the adjugate formula \(L^{-1}_{id} = \mathrm{adj} ({\boldsymbol{L}})_{id}/\det\boldsymbol{L}\), this requires the \((i,d)\) cofactor of \(\boldsymbol{L}\) to vanish. A systematic analysis shows that for \(N = 2\) or \(N = 3\) modes, Eq. 22 cannot be satisfied unless at least one mode is totally lossless or has gain (see Supplemental Materials [32]). A simpler way to achieve such a solution is to use \(N = 4\) lossy modes in a ring arrangement, as depicted in Fig. 1 (a). We will focus on this configuration.

The behavior away from the optimal superbunching point can also be derived within this framework. Suppose \(\alpha_i\) is not merely nonzero but dominant: \[|\alpha_j|^2 \gg \langle \hat{d}_j^\dagger \hat{d}_j \rangle.\] The denominator of \(g^{(2)}\) then simplifies as \(\langle \hat{a}_j^\dagger \hat{a}_j \rangle \approx |\alpha_j|^2\). Expanding the numerator \(\langle \hat{a}_j^\dagger(0)\, \hat{a}_i^\dagger(\tau)\, \hat{a}_i(\tau)\, \hat{a}_j(0) \rangle\) and keeping only terms up to first order in \(U_j\alpha_j^2\), we obtain \[\begin{gather} \langle \hat{a}_j^\dagger(0)\, \hat{a}_i^\dagger(\tau)\, \hat{a}_i(\tau)\, \hat{a}_j(0) \rangle \\ \approx |\alpha_i\alpha_j|^2 - 2\operatorname{Re}\!\left\{\alpha_i^*\alpha_j^*\,\mathcal{F}_\tau^{-1}\!\left[\boldsymbol{G}_\omega\,\boldsymbol{\Gamma}\,\boldsymbol{G}_\omega^\dagger\,\boldsymbol{U}\,\boldsymbol{G}_{-\omega}\right]_{ij}\right\}. \end{gather}\] Hence, in this regime, \[g_{ij}^{(2)}(\tau) \approx 1 - \frac{2\operatorname{Re}\!\left\{\alpha_i^*\alpha_j^*\,\mathcal{F}_\tau^{-1}\!\left[\boldsymbol{G}_\omega\,\boldsymbol{\Gamma}\,\boldsymbol{G}_\omega^\dagger\,\boldsymbol{U}\,\boldsymbol{G}_{-\omega}\right]_{ij}\right\}}{|\alpha_i\alpha_j|^2}. \label{eq:S95g295away}\tag{23}\]

4 Numerical results↩︎

We now want to verify the superbunching predictions of the previous section. Instead of the linearized Eq. 11 , from which the analytic predictions were obtained, let us return to the full quantum master equation 10 , and solve it numerically with all nonlinear terms (cubic and quartic) included. For the model, we take the four-mode configuration of Fig. 1 (a) and choose \(d=1\) as the driven mode, and \(i = 2\) as the signal mode.

We target a set of parameters applicable to the realistic setup of Fig. 1 (b), which aims to realize the model with two coupled silicon carbide microring resonators, utilizing the clockwise and counterclockwise modes in each ring. The setup has three sets of inter-mode couplings: a tunable inter-ring coupling \(J\), and two fixed intra-ring couplings \(J' = 1\,\mu\text{eV}\) and \(J'' = 0.5\,\mu\text{eV}\). All modes have the same decay rate \(\gamma_j = 2\,\mu\text{eV}\) and Kerr coefficient \(U_i = 4.7\times10^{-6}\,\mu\text{eV}\). Notably, the Kerr coefficient matches the intrinsic optical nonlinearity of silicon carbide, without quantum dots or other inclusions. The coherent drive is fixed at \(F = 200\,\mu\text{eV}\), corresponding to an input power of \(P_\mathrm{in} = \omega_L F / 2\pi \approx 6\,\text{nW}\) at wavelength \(\lambda = 1550\,\text{nm}\). Our chosen parameters are based on the SiC waveguide simulations; for details, see the Supplemental Materials [32].

With the above parameters fixed, we can vary the inter-ring coupling \(J\) and detuning \(\Delta\) (the same for all modes), so as to satisfy the zero-displacement condition 22 . Note that these are standard active control parameters in integrated microring platforms; in particular, fine adjustments to the effective couplings can be achieved by thermo-optic or interferometric means [33]–[35] (see Supplemental Materials [32]).

To perform numerical calculations of \(g^{(2)}\), we first obtain the classical amplitudes \(\alpha_j\) by solving Eq. 7 using a modification of the Powell hybrid method as implemented in MINPACK [36]. The effective Hamiltonian \(\mathcal{H}_\mathrm{eff}\) and the Lindblad collapse operators are then constructed in a truncated Fock basis with per-mode photon-number cutoffs \(N_1=7\), \(N_2=4\), \(N_3=5\), and \(N_4=5\), which are later self-consistently validated to be significantly larger than the obtained steady-state per-mode photon numbers. Setting \(d\hat{\rho}/dt=0\) yields a sparse linear system on the vectorized density matrix (the Liouvillian superoperator), which is solved using the direct LGMRES iterative method implemented in QuTiP [37]. We then evaluate \(g^{(2)}(\tau)\) via the quantum regression theorem [37], [38]. As a final check, we verify that the results are insensitive to further increases in the Fock cutoffs \(N_j\).

It is worth noting that the aforementioned Fock space cutoffs refer to the displaced operators \(\hat{d}_j\), not the original \(\hat{a}_j\). Hence, the calculations are valid and tractable even away from the zero-displacement point, where the full photon occupation \(\langle \hat{n}_j \rangle\) may not be small.

Fig. 1 (c) shows the resulting values of \(g^{(2)}(0)\) as a function of the detuning \(\Delta\) and coupling perturbation \(\delta J\). A single sharp peak of superbunching is found at the point where \(\alpha_i = 0\). Its maximum reaches \(g^{(2)}(0) \approx 135\), which is well within the superbunched regime. At this point, the signal cavity has mean photon occupation \(\langle n_i \rangle \approx 6.7\times10^{-4}\), so for a realistic output-waveguide coupling rate of \(\gamma_\mathrm{out} \approx 0.5\,\mu\text{eV}\) (see Supplemental Materials [32]), the photon emission rate is \(R_\mathrm{em} = \langle n_i \rangle \gamma_\mathrm{out} / (2\pi\hbar) \approx 80\,\text{kHz}\), compatible with commercial single-photon detectors.

As we have seen, there is a trade-off between the degree of superbunching that can be achieved using this method and the total photon population, since the superbunching point coincides with a zero of the classical field. If one can tolerate lower photon emission rates (e.g., by having ultra-low detector noise levels), the trade-off may be further exploited to achieve extraordinarily large \(g^{(2)}(0)\).

In Fig. 2 (a), we plot both \(g^{(2)}(0)\) and \(\langle \hat{n}_i \rangle\) at the optimal superbunching point against the drive amplitude \(F\) (i.e., for each \(F\), the couplings and detuning are re-tuned to the optimal superbunching point). The results indicate that if one is able and willing to operate at very low photon count rates, it is possible (say) \(g^{(2)}(0) > 3\times10^4\) at \(F=50\,\mu\text{eV}\), a level of superbunching previously reported only in a macroscopic nonlinear photonic crystal fiber [18]. In principle, even larger values are accessible at weaker drives.

To verify that the numerically-observed trade-off follows the theory quantitatively, Fig. 2 (b) directly plots \(g^{(2)}(0)\) against \(\langle \hat{n}_i \rangle\), using the numerically-obtained values at each optimal superbunching point under varying \(F\), alongside the theoretical scaling relation \(g^{(2)}(0)-2 = C/\langle \hat{n}_i \rangle\) from Eq. 20 . Here, the coefficient \(C\) is chosen so that the value at \(F=250\,\mu\text{eV}\) agrees with the numerical result. Evidently, the scaling relation accurately predicts the numerical results, including in the giant superbunching regime \(g^{(2)}(0) \gg 10^2\).

These predictions may be directly compared to previous experimental superbunching results on other sytems. For optically-driven giant superbunching from a single perovskite quantum dot at cryogenic temperatures, \(g^{(2)}(0) \approx 9\) has been obtained with an emission rate of \(R_{\textrm{em}} \approx 533\,\text{kHz}\) (at \(100\,\text{nW}\) pump power), with a maximum of \(g^{(2)}(0) \approx 30\) at weaker pump powers [17]. For four-wave mixing in a nonlinear photonic crystal fiber, \(g^{(2)}(0) \approx 1200\) has been reported with a count rate of \(200\,\text{kHz}\) [39]. For cathodoluminescence from nanodiamond nitrogen-vacancy centers, \(g^{(2)}(0) \approx 49\) was achieved with \(R_{\textrm{em}} \sim 10\,\text{kHz}\) [40].

Therefore, our proposed approach supports two complementary operating regimes. First, at \(F=200\,\mu\text{eV}\), it can achieve larger \(g^{(2)}(0)\) than several previous compact sources while maintaining a comparable emission rate of \(R_{\textrm{em}} \approx 80\,\text{kHz}\). Second, if the drive is reduced to \(F=50\,\mu\text{eV}\), one can reach \(g^{(2)}(0)>3\times10^4\), well above established nonlinear photonic crystal fiber results [39] (and comparable only to a very recent claim based on another photonic crystal fiber [18])—albeit with a much lower count rate.

Finally, Fig. 2 (c) shows the time-delayed second-order correlation \(g^{(2)}(\tau)\) at the optimal superbunching point in Fig. 1 (c). The correlation exhibits a sharp peak at \(\tau = 0\) and decays on a timescale set by the cavity linewidth \(\gamma\), confirming that the superbunching is a genuine quantum correlation effect rather than arising from classical intensity fluctuations.

5 Conclusions↩︎

We have shown that photon superbunching can be achieved in a system of coupled optical modes subject to only weak Kerr nonlinearities. The essential physics is captured by Eq. 1 , which predicts a divergence in the second-order photon correlation as the nonlinearity strength \(U\) and the photon population \(\langle \hat{n}_i\rangle\) both approach zero. It is striking that a weaker nonlinearity induces a stronger quantum effect, which runs counter to the conventional wisdom that strong photon correlations demand strong interactions. The reason is that interference is utilized to suppress the classical amplitude, so as to allow quantum fluctuations to dominate.

We have presented an exemplary setup for realizing this phenomenon, based on coupled silicon carbide ring resonators with nanowatt-level driving, standard forms of active parameter tuning, and other design elements that are commonplace in integrated photonics [30], [31]. When the system is driven at \(F = 200\,\mu\text{eV}\), the photons in the signal cavity reach \(g^{(2)}(0) \approx 135\) with an emission rate \(80\,\text{kHz}\), which is measurable with commercially available single-photon detectors without employing any unusual ultralow-noise or ultrahigh-efficiency techniques. If the drive may be reduced to \(F = 50\,\mu\text{eV}\), we predict \(g^{(2)}(0) \gtrsim 3\times10^4\) [Fig. 2 (a)]; higher values of \(g^{(2)}(0)\) are accessible so long as one can tolerate a lower emission rate from such a light source.

So far as we know, this route to achieving arbitrary strong levels of superbunching has not previously been explored. Although our predictions were derived from a linearized theory, they are confirmed by quantum master equation simulations retaining all nonlinear terms.

Several directions merit further investigation. Engineering cavity couplings in larger lattices [41] may enable multiphoton superbunching through higher-order \(g^{(N)}\), which is relevant to quantum-enhanced metrology and ghost imaging beyond the two-photon level. Unlike conventional schemes that rely on intrinsic emitter properties, our mechanism requires only coupled nonlinear cavities and a coherent drive; however, combining both approaches, by embedding quantum emitters in coupled cavity networks, may allow collective emission dynamics to be further engineered by intercavity interference. This could be used to access correlation regimes inaccessible to either mechanism alone.

6 Effective fluctuation Hamiltonian↩︎

In this Appendix, we derive the effective Hamiltonian 8 , which serves as the basis of the superbunching prediction. We start from the system Hamiltonian 2 , which describes \(N\) coherently-driven coupled modes with Kerr nonlinearities, in a co-rotating frame; the quantum fluctuation operators \(\hat{d}_j\) are split out via Eq. 6 , with the mean fields given by Eq. 7 . To proceed, we expand the Hamiltonian, substituting Eq. 6 into each term of the Hamiltonian and collect contributions by order in \(\hat{d}_j\).

The first (detuning) term in the Hamiltonian yields \[\begin{align} -\Delta_j \hat{a}_j^\dagger \hat{a}_j &= -\Delta_j \left[|\alpha_j|^2 + \alpha_j^* \hat{d}_j + \alpha_j \hat{d}_j^\dagger + \hat{d}_j^\dagger \hat{d}_j\right]. \end{align}\] For the second term, which corresponds to the Kerr nonlinearity, we expand \((\alpha_j^* + \hat{d}_j^\dagger)^2(\alpha_j + \hat{d}_j)^2\) exactly: \[\begin{align} &U_j \hat{a}_j^\dagger \hat{a}_j^\dagger \hat{a}_j \hat{a}_j \nonumber \\ &= U_j \bigg[|\alpha_j|^4 + 2|\alpha_j|^2 \alpha_j^* \hat{d}_j + 2|\alpha_j|^2 \alpha_j \hat{d}_j^\dagger \nonumber \\ &\quad\qquad + 4|\alpha_j|^2 \hat{d}_j^\dagger \hat{d}_j + \alpha_j^2 (\hat{d}_j^\dagger)^2 + (\alpha_j^*)^2 \hat{d}_j^2 \nonumber \\ &\quad\qquad + 2\alpha_j^* \hat{d}_j^\dagger \hat{d}_j^2 + 2\alpha_j (\hat{d}_j^\dagger)^2 \hat{d}_j + (\hat{d}_j^\dagger)^2 \hat{d}_j^2\bigg]. \label{eq:S95kerr95expand} \end{align}\tag{24}\] For the coupling terms, we obtain \[\begin{align} J_{jk} \hat{a}_j^\dagger \hat{a}_k &= J_{jk} \left[\alpha_j^* \alpha_k + \alpha_j^* \hat{d}_k + \alpha_k \hat{d}_j^\dagger + \hat{d}_j^\dagger \hat{d}_k\right]. \end{align}\] Finally, the driving terms give \[\begin{align} F_j \hat{a}_j^\dagger + F_j^* \hat{a}_j &= F_j \alpha_j^* + F_j^* \alpha_j + F_j \hat{d}_j^\dagger + F_j^* \hat{d}_j. \end{align}\]

We now collect terms by order. The zeroth order terms correspond to a constant energy shift, which only affects the global phase and therefore may be ignored. As for the first order terms, we find that the coefficients exactly match the classical steady-state equation 7 , so as to form \[\begin{align} \mathcal{H}_1 = \sum_j \frac{i\gamma_j}{2}(\alpha_j \hat{d}_j^\dagger - \alpha_j^* \hat{d}_j) + \mathrm{h.c.}, \end{align}\] where “\(\textrm{h.c.}\)” stands for “Hermitian conjugate”. We should also note that the displacement transformation modifies the Lindblad dissipator: \[\gamma_j\mathcal{D}[\alpha_j + \hat{d}_j]\rho = \gamma_j\mathcal{D}[\hat{d}_j]\rho - i[\mathcal{H}_\mathrm{diss}^{(j)}, \rho],\] where \[\mathcal{H}_\mathrm{diss}^{(j)} = \frac{i\gamma_j}{2}\!\left(\alpha_j^* \hat{d}_j - \alpha_j \hat{d}_j^\dagger\right).\] Hence, the linear terms cancel exactly between the Hamiltonian and dissipative contributions: \[\mathcal{H}_1 + \sum_j \mathcal{H}_\mathrm{diss}^{(j)} = 0,\] which self-consistently satisfies \(\langle \hat{d}_j \rangle = 0\).

Thus, the fluctuation dynamics starts at quadratic order. Collecting these terms gives \[\mathcal{H}_2 = \sum_j \left[\tilde{\omega}_j\, \hat{d}_j^\dagger \hat{d}_j + U_j\!\left(\alpha_j^2 (\hat{d}_j^\dagger)^2 + \textrm{h.c.}\right)\right] + \sum_{j \neq k} J_{jk}\, \hat{d}_j^\dagger \hat{d}_k,\] where \(\tilde{\omega}_j\) is the renormalized frequency detuning given in Eq. 9 . The terms of order \(\alpha^2\) describe squeezing (parametric) processes induced by the classical field.

Next, we handle the the third- and fourth-order terms, which arise solely from the Kerr nonlinearity: \[\begin{align} \mathcal{H}_3 &= \sum_j 2U_j \left(\alpha_j^*\, \hat{d}_j^\dagger \hat{d}_j^2 + \alpha_j\, (\hat{d}_j^\dagger)^2 \hat{d}_j\right). \tag{25} \\ \mathcal{H}_4 &= \sum_j U_j\, \hat{d}_j^\dagger \hat{d}_j^\dagger \hat{d}_j \hat{d}_j. \tag{26} \end{align}\] Combining all non-constant orders, we obtain the effective Hamiltonian 8 . Note that this result is exact—we have not dropped any terms in the expansion.

7 The linearized Langevin equation↩︎

Here, we give the derivation leading to Eqs. 13 –16 . Starting from the effective Hamiltonian 8 and the master equation 10 , we retain only terms up to quadratic order in \(\hat{d}_j\). The quantum Langevin equations become \[\begin{align} i\partial_t \hat{d}_k &= \!\left(\tilde{\omega}_k - i\frac{\gamma_k}{2}\right)\! \hat{d}_k + 2U_k \alpha_k^2 \hat{d}_k^\dagger + \sum_l J_{kl} \hat{d}_l - i\gamma_k^{\frac{1}{2}}\hat{\eta}_k, \label{eq:S95langevin95single} \end{align}\tag{27}\] where \(\hat{\eta}_k\) is the vacuum input noise operator satisfying \(\langle \hat{\eta}_j(t)\, \hat{\eta}_k^\dagger(t') \rangle = \delta_{jk}\,\delta(t-t')\) and \(\langle \hat{\eta}_j(t)\, \hat{\eta}_k(t') \rangle = 0\). In Eq. 11 –12 , these equations were casted into matrix form by introducing \(\boldsymbol{d}\) (collecting the \(\hat{d}_j\) operators), as well as the matrices \(\boldsymbol{H}\), \(\boldsymbol{\Gamma}\), and \(\boldsymbol{U}\) (collecting the detunings/couplings, loss rates, and nonlinearity). Similarly, taking the Hermitian conjugate of Eq. 27 yields \[-i\partial_t \boldsymbol{d}^\dagger = (\boldsymbol{H} + i\boldsymbol{\Gamma}/2)\,\boldsymbol{d}^\dagger + \boldsymbol{U}^*\,\boldsymbol{d} + i\sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}^\dagger. \label{eq:S95langevin95adj}\tag{28}\] Eqs. 11 and 28 combine into a Bogoliubov equation \[\begin{gather} i\partial_t \begin{pmatrix} \boldsymbol{d} \\ \boldsymbol{d}^\dagger \end{pmatrix} = \begin{pmatrix} \boldsymbol{H} - i\boldsymbol{\Gamma}/2 & \boldsymbol{U} \\ -\boldsymbol{U}^* & -(\boldsymbol{H} + i\boldsymbol{\Gamma}/2) \end{pmatrix} \begin{pmatrix} \boldsymbol{d} \\ \boldsymbol{d}^\dagger \end{pmatrix} \\ - i \begin{pmatrix} \sqrt{\boldsymbol{\Gamma}}\boldsymbol{\eta} \\ \sqrt{\boldsymbol{\Gamma}}\boldsymbol{\eta}^\dagger \end{pmatrix}. \label{eq:S952N} \end{gather}\tag{29}\] To solve this, we define the Fourier transformed operators \[\hat{d}_k(\omega) = \int dt\, e^{i\omega t}\, \hat{d}_k(t).\] In the frequency domain, Eq. 29 becomes \[\begin{pmatrix} -i\boldsymbol{d}_\omega \\ i\boldsymbol{d}_{-\omega}^\dagger \end{pmatrix} = \boldsymbol{\mathcal{M}}_\omega^{-1} \begin{pmatrix} \sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}_\omega \\ \sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}_{-\omega}^\dagger \end{pmatrix}, \label{eq:S95fourier95sol}\tag{30}\] where \[\begin{align} \boldsymbol{\mathcal{M}}_\omega &= \begin{pmatrix} \boldsymbol{H} - i\boldsymbol{\Gamma}/2 - \omega & -\boldsymbol{U} \\ -\boldsymbol{U}^* & \boldsymbol{H} + i\boldsymbol{\Gamma}/2 + \omega \end{pmatrix}, \tag{31} \\ \langle \boldsymbol{\eta}_{\omega_1} \boldsymbol{\eta}_{-\omega_2}^{\dagger T} \rangle &= 2\pi\,\delta(\omega_1 + \omega_2)\,\mathbb{I}, \quad \langle \boldsymbol{\eta}_{\omega_1} \boldsymbol{\eta}_{\omega_2}^T \rangle = 0. \tag{32} \end{align}\] Next, let us define the retarded Green’s function as the inverse of the upper-left block in \(\boldsymbol{\mathcal{M}}_\omega\): \[\boldsymbol{G}_\omega = (\boldsymbol{H} - i\boldsymbol{\Gamma}/2 - \omega)^{-1}. \label{eq:S95Green}\tag{33}\] Note that for real \(\omega\), \[\begin{align} \boldsymbol{G}_{-\omega}^\dagger = (\boldsymbol{H} + i\boldsymbol{\Gamma}/2 + \omega)^{-1}. \end{align}\] We can then compute the inverse of \(\boldsymbol{\mathcal{M}}_\omega\) via the Schur complement. To leading order in the nonlinearity term \(\boldsymbol{U}\), Eq. 30 reduces to \[\begin{align} -i\boldsymbol{d}_\omega &\approx \boldsymbol{G}_\omega \sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}_\omega + \boldsymbol{G}_\omega\, \boldsymbol{U}\, \boldsymbol{G}_{-\omega}^\dagger \sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}_{-\omega}^\dagger, \tag{34}\\ i\boldsymbol{d}_{-\omega}^\dagger &\approx \boldsymbol{G}_{-\omega}^\dagger\, \boldsymbol{U}^\dagger\, \boldsymbol{G}_\omega\, \sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}_\omega + \boldsymbol{G}_{-\omega}^\dagger\sqrt{\boldsymbol{\Gamma}}\,\boldsymbol{\eta}_{-\omega}^\dagger. \tag{35} \end{align}\] Using Eq. 32 , we can calculate correlators like \[\begin{align} \langle \boldsymbol{d}_{\omega_1}\, \boldsymbol{d}_{-\omega_2}^T \rangle &= -2\pi\, \boldsymbol{G}_{\omega_1} \boldsymbol{\Gamma} \big[\boldsymbol{G}_{\omega_2}\, \boldsymbol{U}\, \boldsymbol{G}_{-\omega_2}^\dagger\big]^T \;\delta(\omega_1+\omega_2), \end{align}\] and, likewise, \(\langle \boldsymbol{d}_{-\omega_1}^\dagger\, \boldsymbol{d}_{\omega_2}^T \rangle\), \(\langle \boldsymbol{d}_{-\omega_1}^\dagger\, \boldsymbol{d}_{-\omega_2}^{\dagger T} \rangle\), and \(\langle \boldsymbol{d}_{\omega_1}\, \boldsymbol{d}_{-\omega_2}^{\dagger T} \rangle\). Finally, Fourier transforming back to the time domain yields the correlation functions 13 –16 .

Supplemental Materials for
Giant Photon Superbunching from Weak Nonlinearity
Y. Wang, X. Zheng, T. C. H. Liew, and Y. D. Chong

8 Minimum number of cavities for zero displacement↩︎

Here we show that a minimum of four cavities is required to satisfy the zero-displacement condition \(L^{-1}_{id} = 0\) [Eq. 22 of the main text] under the practical requirement that every cavity has nonzero loss (\(\gamma_j > 0\)) and the inter-cavity couplings \(J_{jk}\) are real and nonvanishing. We consider a general scenario in which all couplings and the complex detunings \(z_j = -\Delta_j - i\gamma_j/2\) are independently tunable.

8.1 Two cavities↩︎

For two cavities coupled by \(J\): \[\boldsymbol{L} = \begin{pmatrix} z_1 & J \\ J & z_2 \end{pmatrix}, \quad \boldsymbol{L}^{-1} = \frac{1}{z_1 z_2 - J^2} \begin{pmatrix} z_2 & -J \\ -J & z_1 \end{pmatrix}. \label{eq:S95L2}\tag{36}\] The off-diagonal element \(L^{-1}_{12} = -J/(z_1 z_2 - J^2)\) vanishes only if \(J = 0\) (trivially decoupled). The diagonal element \(L^{-1}_{11} = z_2/(z_1 z_2 - J^2)\) vanishes only if \(z_2 = 0\), i.e., \(\Delta_2 = 0\) and \(\gamma_2 = 0\), requiring a lossless on-resonance cavity. Neither option is practical.

8.2 Three cavities↩︎

For three cavities with all pairwise couplings \(J_1\), \(J_2\), \(J_3\): \[\boldsymbol{L} = \begin{pmatrix} z_1 & J_1 & J_3 \\ J_1 & z_2 & J_2 \\ J_3 & J_2 & z_3 \end{pmatrix}. \label{eq:S95L3}\tag{37}\] The cofactors entering \(\boldsymbol{L}^{-1}\) are \(L^{-1}_{ij} = \mathrm{adj} ({\boldsymbol{L}})_{ij}/\det\boldsymbol{L}\), where \[\begin{align} \mathrm{adj} ({\boldsymbol{L}})_{11} &= z_2 z_3 - J_2^2, \tag{38} \\ \mathrm{adj} ({\boldsymbol{L}})_{12} &= J_2 J_3 - z_3 J_1, \tag{39} \end{align}\] and the remaining cofactors follow by permutation. There are two cases:

(i) Diagonal zero. Setting \(\mathrm{adj} ({\boldsymbol{L}})_{11} = 0\) requires \(z_2 z_3 = J_2^2\). Taking modulus and argument separately: \(|z_2|\,|z_3| = J_2^2\) and \(\arg(z_2) + \arg(z_3) = 0\). Since \(\mathrm{Im}(z_j) = -\gamma_j/2 < 0\), the argument condition forces \(z_2 = z_3^*\), which demands that one cavity has gain precisely balanced against another’s loss, or that both cavities are lossless.

(ii) Off-diagonal zero. Setting \(\mathrm{adj} ({\boldsymbol{L}})_{12} = 0\) requires \(z_3 = J_2 J_3/J_1\), which is purely real. This is only possible when \(\gamma_3 = 0\)—the cavity must be entirely lossless.

Both requirements are impractical in typical photonic systems.

8.3 Four cavities in a ring↩︎

We now consider four cavities arranged in a ring with nearest-neighbor couplings \(J_{12}\), \(J_{23}\), \(J_{34}\), \(J_{41}\) and no cross-couplings: \[\boldsymbol{L} = \begin{pmatrix} z_1 & J_{12} & 0 & J_{41} \\ J_{12} & z_2 & J_{23} & 0 \\ 0 & J_{23} & z_3 & J_{34} \\ J_{41} & 0 & J_{34} & z_4 \end{pmatrix}. \label{eq:S95L4}\tag{40}\] Taking cavity 1 as the driven site and cavity 2 as the signal site, the zero-displacement condition is \(L^{-1}_{21} = 0\), i.e., \(\mathrm{adj} ({\boldsymbol{L}})_{21} = 0\). Evaluating the \((2,1)\) cofactor: \[\mathrm{adj} ({\boldsymbol{L}})_{21} = -\bigl[J_{12}(z_3 z_4 - J_{34}^2) + J_{41} J_{34} J_{23}\bigr]. \label{eq:S95adj21954cav}\tag{41}\] Setting this to zero gives \[z_3 z_4 = J_{34}^2 - \frac{J_{41}\, J_{34}\, J_{23}}{J_{12}}. \label{eq:S954cav95cond}\tag{42}\] The right-hand side is real, while the left-hand side \(z_3 z_4 = (-\Delta_3 - i\gamma_3/2)(-\Delta_4 - i\gamma_4/2)\) has imaginary part \(-(\Delta_3\gamma_4 + \Delta_4\gamma_3)/2\), which can be zeroed by choosing appropriate detunings. For identical cavities (\(z_j = z = -\Delta - i\gamma/2\)), the condition becomes \(z^2 = J_{34}^2 - J_{41} J_{34} J_{23}/J_{12}\), whose imaginary part \(\Delta\gamma = 0\) vanishes at \(\Delta = 0\), leaving the purely real condition \[\frac{\gamma^2}{4} = \frac{J_{41}\, J_{34}\, J_{23}}{J_{12}} - J_{34}^2. \label{eq:S95gamma95cond}\tag{43}\] This is readily satisfied when \(J_{41} J_{23}/J_{12} > J_{34}\), demonstrating that the four-cavity ring admits zero-displacement solutions with all cavities having nonzero loss. Four is therefore the minimum number of cavities that enables giant superbunching under practical conditions.

9 Physical realization using ring resonators↩︎

To achieve the four coupled optical modes of our model, we can employ either four coupled microring resonators, as shown in Fig. 3 (a), or two coupled microring resonators in which the clockwise (CW) and counterclockwise (CCW) modes of each ring serve as the four bosonic modes, as shown in Fig. 3 (b), or Fig. 1 (b) of the main text. We focus on the two-ring realization.

Here, the CW and CCW modes of each ring have the same natural frequency, and are coupled to each other with rates \(J'\) (second ring) and \(J''\) (first ring). In practice, such couplings can be achieved via deliberately-introduced subwavelength defects or surface roughness [42]–[45]. Assuming both rings have the same shape and size, the Langevin equations of motion for the four modes \(\hat{a}_1\) (CW, ring 1), \(\hat{a}_2\) (CCW, ring 1), \(\hat{a}_3\) (CW, ring 2), \(\hat{a}_4\) (CCW, ring 2) are \[i\frac{d}{dt}\begin{pmatrix}\hat{a}_1\\\hat{a}_2\\\hat{a}_3\\\hat{a}_4\end{pmatrix}= \begin{pmatrix} z_1 & J'' & 0 & J\\ J'' & z_1 & J & 0\\ 0 & J & z_2 & J'\\ J & 0 & J' & z_2 \end{pmatrix} \begin{pmatrix}\hat{a}_1\\\hat{a}_2\\\hat{a}_3\\\hat{a}_4\end{pmatrix} +\begin{pmatrix}i\sqrt{\gamma_{\mathrm{in}}}\,\alpha_{\mathrm{in}}\\0\\0\\0\end{pmatrix}, \label{eq:S95langevin95ring}\tag{44}\] where \(z_1 = -\Delta - i(\gamma+\gamma_{\mathrm{in}}+\gamma_{\mathrm{out}})/2\), \(z_2 = -\Delta - i\gamma'/2\), \(\Delta=\omega_c-\omega_L\) is the cavity–laser detuning, \(\gamma\) is the intrinsic loss rate, \(\gamma_{\mathrm{in}}\) is the input-waveguide coupling rate, \(\gamma_{\mathrm{out}}\) is the output-waveguide coupling rate, \(J\) is the inter-ring coupling, and \(\alpha_{\mathrm{in}}\) is the input field amplitude. The output field is \(\alpha_{\mathrm{out}}=\sqrt{\gamma_{\mathrm{out}}}\,\hat{a}_2\). The Kerr nonlinear Hamiltonian for the four modes is \(\mathcal{H}_{nl} = U\sum_{i=1}^{4} \hat{a}_i^{\dagger}\hat{a}_i^{\dagger}\hat{a}_i \hat{a}_i\), where \(U\) is the Kerr coefficient of each mode.

We perform Finite Element Method (FEM) simulations using material parameters drawn from the experimental literature. We assume each ring is composed of silicon carbide (SiC) with ring radius \(R=3\,\mu\)m, waveguide thickness \(t=0.35\,\mu\)m, and waveguide width \(w=0.8\,\mu\)m, operating at wavelength \(\lambda\approx1550\,\text{nm}\). The refractive indices are \(n_\mathrm{SiC}=2.45\) for the waveguides and \(n(\text{SiO}_2)=1.44\) for the cladding, and the nonlinear refractive index for SiC is \(n_2=4.8\times10^{-6}\,\mu\text{m}^2/\text{W}\). The Kerr coefficient is estimated by [46], [47] \[U=\frac{c\hbar^2\omega_c^2}{n_\mathrm{SiC}^2\,V}\,n_2,\] where the mode volume is approximated by the ring volume \(V=2\pi Rwt\). This gives \(U\approx4.7\times10^{-6}\,\mu\text{eV}\), consistent with the value used in the main text simulations.

Fig. 3 (c) shows the inter-ring coupling strength \(J\) as a function of gap distance \(d_\mathrm{gap}\), extracted from FEM simulations. In each simulation, the eigenvalue splitting between the symmetric and antisymmetric supermodes of the coupled system equals \(2J\). The FEM identifies an eigenmode near \(\lambda=1553\,\text{nm}\) with a simulated intrinsic quality factor \(Q_\mathrm{int}\approx1.5\times10^6\); accounting for additional loss channels including waveguide coupling (\(\gamma_\mathrm{in}, \gamma_\mathrm{out}\approx 0.5\,\mu\)eV), the loaded loss rate of ring-1 modes (\(\hat{a}_1\), \(\hat{a}_2\)) is \(\gamma+\gamma_\mathrm{in}+\gamma_\mathrm{out}\), while ring-2 modes (\(\hat{a}_3\), \(\hat{a}_4\)) have only the intrinsic rate \(\gamma'\) (loaded \(Q\approx4\times10^5\)). For simplicity, all four modes are assigned the same loss rate \(\gamma_j = 2\,\mu\)eV in the simulations, corresponding to the loaded quality factor of ring 1. The general theoretical model supports site-dependent \(\gamma_j\) and this simplification does not qualitatively affect the conclusions. The CWCCW couplings are set to \(J'=1\,\mu\text{eV}\) and \(J''=0.5\,\mu\text{eV}\), which are achievable via subwavelength surface defects [42]–[45].

In an actual experiment, the system parameters may deviate from nominal values, so it is desirable to have \(\Delta\) and \(J\) actively tunable, so as to achieve the superbunching condition:

• The detuning \(\Delta\) can be adjusted simply via the laser frequency \(\omega_L\).

• The coupling \(J\) is determined principally by the gap distance \(d_\mathrm{gap}\) [Fig. 3 (c)], but can be further fine-tuned by engineering the refractive index in and around the gap region. In integrated microcavities, such tuning can be implemented with an auxiliary control laser [48]; laser-induced heating and photoexcited carriers can modify the refractive index of semiconductors [49], [50], and hence the inter-ring coupling.

  In Fig. 4 (b), we illustrate the dependence of \(J\) on the regional refractive index change \(\Delta n\) within a spot of radius \(r=0.5\) \(\mu\)m centered on the inter-ring gap, assuming an initial gap size of \(d_{\textrm{gap}}=1.74\) \(\mu\)m. For realistic levels of \(\Delta n\), it is possible to achieve variations in \(J\) compatible with our proposal.

References↩︎