2.1.

Two Self-Sustained Microresonators and Interaction Model

As shown in Fig. 1(a), the system is composed of two optical microcavities with different resonant frequencies, ${\omega }_{10}$ and ${\omega }_{20}$, coupled at the strength $g$. The $j$’th ($j=\mathrm{1,2}$) cavity is self-sustained by the internal gain described by the factor $G_{\mathrm{p},j}$ and the intrinsic dissipation at the rate ${\kappa }_j$. In the presence of nonlinear gain, a self-Kerr type modulation ${\delta }_j{(a_j^{†}{a}_j)}^2/2$ is present, with $a_j$ being the annihilation operator and the factor ${\delta }_j$ describing the self-Kerr effect.³⁶ With the gain saturation³⁷ or the multiphoton absorption,^38–40 the effective dissipation of the $j$’th mode is modeled as $K_j={\kappa }_j/2+R_j⟨a_j^{†}{a}_j⟩$, where $R_j$ is the nonlinearity factor.

Fig. 1

Schematic diagram of the system. (a) Two detuned and self-sustained optical microcavities with different resonant frequencies, ${\omega }_{10}$ and ${\omega }_{20}$, which are directly coupled at strength $g$. (b)–(d) Frequency spectra of the coupled cavities, showing three different long-term states: unsynchronized, limit cycle (LC), and synchronized (Sync.). Light blue represents the noise backgrounds from which the first- and second-order synchronizations are distinguished.

The dissipative evolution of the system is described by the Lindblad density-matrix equation ($\hslash =1$ hereafter),

2.2.

Synchrony Solution in Static Case

We focus on the phase difference and the transient frequencies of two modes in the coherent-state representation.⁶ In this representation, the complex amplitude ${\alpha }_j=⟨a_j⟩$ is parameterized as $r_j\sqrt{G/R}{e}^{-i{\varphi }_j}$, where $r_j$ and ${\varphi }_j$ are the amplitude and the phase, respectively. Let $\varphi ={\varphi }_1-{\varphi }_2$ be the phase difference, which is the preserved degree of freedom, and let ${\omega }_j={\dot{\varphi }}_j$ be the transient frequency for the $j$’th mode.

Following the standard Wigner function formalism,⁴³ the mode equation is described by ${(\dot{\mathbf{\Lambda }},\overline{\dot{\mathbf{\Lambda }}})}^{⊺}=\mathbf{f}(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})$, where $\mathbf{\Lambda }=({\alpha }_1,{\alpha }_2)$ and $\mathbf{f}(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})$ denote the quasiprobability drift flow of the two modes (see Supplementary Material for details). The synchrony solution is achieved when $\mathbf{f}(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})=0$, and a fixed point ${\mathbf{\Lambda }}_{\mathrm{s}}(\tilde{\mathrm{\Delta }},\tilde{\delta },\tilde{g})$ emerges in the parametric space (see Supplementary Material for details). In Fig. 2, we plot the phase differences and the transient frequencies for different $\tilde{g}$. Three categories of long-term behaviors are discovered. When the coupling strength is low, $\tilde{g}=0.3$, for example, the phase difference $\varphi$ accumulates to infinity quickly [see Fig. 2(a1)], and the transient frequencies ${\omega }_1$ and ${\omega }_2$ are effectively separated [Fig. 2(a2)], showing two separated modes in the frequency spectrum [Fig. 1(b)]. When the coupling strength is turned higher, $\tilde{g}=0.398$, for example, the phase difference $\varphi$ vibrates around the stationary point but does not accumulate [Fig. 2(b1)], and the frequencies ${\omega }_1$ and ${\omega }_2$ breathe slowly around the stationary frequency [Fig. 2(b2)], generating a limit cycle state. A stationary mode is localized around the original two cavity modes, and a pair of weak limit cycle modes can be found symmetrically detuned from the stationary modes [Fig. 1(c)]. Finally, with a high enough coupling strength, such as $\tilde{g}=0.4$, the phase difference $\varphi$ stabilizes [Fig. 2(c1)], and the frequencies ${\omega }_1$ and ${\omega }_2$ also converge to a single value [Fig. 2(c2)], reaching the synchronized state with a single mode in the frequency spectrum [Fig. 1(d)].

Fig. 2

Long-term evolutions of the two cavity modes under different coupling strengths. Three different categories are shown: (a) the unsynchronized ($\tilde{g}=0.3$), (b) limit cycle ($\tilde{g}=0.398$), and (c) synchronized states ($\tilde{g}=0.4$). (a1)–(c1) Phase difference; (a2)–(c2) transient frequencies; (a3)–(c3) trajectory encircling types (black cross as the axis); and (a4)–(c4) dynamical potential near the synchrony point. In all figures, the given detuning $\tilde{\mathrm{\Delta }}=0.3$ and Kerr factor $\tilde{\mathit{\delta }}=0.1$.

It is noted that the temporal translational symmetry (TTS) is preserved in the synchronized state because both the amplitudes and phase difference remain invariant while the symmetry is broken in the unsynchronized and limit cycle states. The discrete topological character number as the average encircling number is further defined as

2.3.

Analysis of Different Transition Types

When we further study the maximum of the frequency differences, $\max |{\omega }_1-{\omega }_2|$, two types of synchronization transitions are found. As shown in Fig. 3(a), when the coupling strength $\tilde{g}$ is low, the maximal frequency difference varies slowly. At a critical strength ${\tilde{g}}_{\mathrm{c}}$, it suddenly falls to zero, which shows the characteristics of the first-order transition. As shown in Fig. 3(b), the maximal frequency difference continuously decreases to zero but has a discontinuity in its derivative at ${\tilde{g}}_{\mathrm{c}}$, showing the feature of the second-order transition. In addition, the noise spectrum is also calculated in long-term motions. For the synchronized spectrum in Fig. 1(d), the background noise has coinciding peaks with synchronized frequencies in the first-order transition, whereas the noise has shifted-away peaks in the second-order transition (see Supplementary Material for details).

Fig. 3

Parameter dependence of the synchronization. (a), (b) Maximum of the frequency differences, $\max |{\omega }_1-{\omega }_2|$, versus the coupling strength $\tilde{g}$, with ($\tilde{\mathrm{\Delta }}=0.2$, $\tilde{\mathit{\delta }}=0.1$) in (a) and ($\tilde{\mathrm{\Delta }}=0.3$, $\tilde{\mathit{\delta }}=0.1$) in (b); inset shows the derivative. (c) Phase diagram in the $(\tilde{\mathrm{\Delta }},\tilde{g})$ plane with the Kerr factor $\tilde{\mathit{\delta }}=0.1$. The inaccessible (gray), limit cycle (dark blue), and synchronized (light blue) regimes are marked. The red cross stands for the triple phase point $({\tilde{\mathrm{\Delta }}}_{\mathrm{T}},{\tilde{g}}_{\mathrm{T}})$. (d) The triple phase point $({\tilde{\mathrm{\Delta }}}_{\mathrm{T}},{\tilde{g}}_{\mathrm{T}})$ depending on the Kerr factor $\tilde{\mathit{\delta }}$.

In order to study the critical coupling strength ${\tilde{g}}_{\mathrm{c}}$ and the transition behaviors in its vicinity, a real-valued dynamical potential $V(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})$ is defined (see Supplementary Material for details). Only if the dynamical potential has a local minimum, the fixed point ${\mathbf{\Lambda }}_{\mathrm{s}}(\tilde{\mathrm{\Delta }},\tilde{\delta },\tilde{g})$ emerges and remains stable, and thus indicates the existence of a synchronized state (see Supplementary Material for details). In the vicinity of the fixed point, the dynamical potential can be expanded as $V(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})=V({\mathbf{\Lambda }}_{\mathrm{s}},{\overline{\mathbf{\Lambda }}}_{\mathrm{s}})-\frac{1}{2}[(\mathrm{\Delta }\overline{\mathbf{\Lambda }},\mathrm{\Delta }\mathbf{\Lambda })·\mathbf{J}·{(\mathrm{\Delta }\overline{\mathbf{\Lambda }},\mathrm{\Delta }\mathbf{\Lambda })}^{†}+\mathrm{H.c.}]$, where $\mathbf{J}(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})=\partial \mathbf{f}(\mathbf{\Lambda },\overline{\mathbf{\Lambda }})/\partial (\mathbf{\Lambda },\overline{\mathbf{\Lambda }})$ is the Jacobian matrix, $\mathrm{\Delta }\mathbf{\Lambda }$ is the arbitrarily small displacement from the fixed point, and H.c. is the Hermitian conjugate. The displacement $\mathrm{\Delta }\mathbf{\Lambda }$ signifies the breaking of the TTS. At $\mathrm{\Delta }\mathbf{\Lambda }=0$, the TTS is preserved. The stability near the fixed point is thus governed by the eigenvalues of the Jacobian. When the largest real part of the $\mathbf{J}$ eigenvalues [known as the largest Lyapunov exponent $\mathcal{L}({\mathbf{\Lambda }}_{\mathrm{s}})$] is positive, the fixed point is unstable and vice versa (see Supplementary Material for details).⁴⁵ The critical coupling strength ${\tilde{g}}_{\mathrm{c}}$ is then taken at $\mathcal{L}({\mathbf{\Lambda }}_{\mathrm{s}})=0$.

In the three-dimensional space $(r_1,r_2,\varphi )$, the Jacobian $\mathbf{J}$ has purely real $3\times 3$ components, and thus the complex eigenvalues must come in pairs. If the largest Lyapunov exponent $\mathcal{L}$ equals one of the eigenvalues, the dynamical potential is simplified as

In the light of the static analysis above, the phase diagram in the $\tilde{\delta }$-cross section is plotted in Fig. 3(c), where three regions of different long-term behaviors are marked. The synchronized and limit cycle regimes are specified according to the existence of a single and double local minima of the dynamical potentials, respectively. The inaccessible (unsynchronized) regime corresponds to the saddle nodes in the dynamical potentials, and thus neither the synchronized state nor the limit cycle state can survive in this regime. The transition from the unsynchronized state to the synchronized state is of first-order and has a variant topological character number. The transition from the limit cycle state to the synchronized state is of second-order and has an invariant topological character number (see Supplementary Material for details). It is also found that a triple phase point emerges at ${\tilde{g}}_{\mathrm{c}}={\tilde{g}}_{\mathrm{T}}$ and $\tilde{\mathrm{\Delta }}={\tilde{\mathrm{\Delta }}}_{\mathrm{T}}$, where ${\tilde{\mathrm{\Delta }}}_{\mathrm{T}}$ is the minimal detuning required for the second-order transition. This point corresponds to the solution where two eigenvalues of the Jacobi matrix $\mathbf{J}$ are zeros. For the detuning $\tilde{\mathrm{\Delta }}<{\tilde{\mathrm{\Delta }}}_{\mathrm{T}}$ ($\tilde{\mathrm{\Delta }}>{\tilde{\mathrm{\Delta }}}_{\mathrm{T}}$), the first-order (second-order) transition happens around ${\tilde{g}}_{\mathrm{c}}$ (black solid line). The triple phase point relies crucially on the strength of the Kerr effect. In Fig. 3(d), we plot ${\tilde{\mathrm{\Delta }}}_{\mathrm{T}}$ and ${\tilde{g}}_{\mathrm{T}}$ with respect to the Kerr factor $\tilde{\delta }$, showing monotone increasing and decreasing dependence, respectively. When the factor $\tilde{\delta }$ increases, the self-tuning ability of the Kerr effect is strengthened, and the second-order synchronization under a larger detuning and a weaker coupling becomes possible.

2.4.

Hysteresis Behavior

The two types of synchronization transitions present distinct hysteresis behaviors near the critical coupling strength ${\tilde{g}}_{\mathrm{c}}$. In Fig. 4, the frequency differences $|{\omega }_1-{\omega }_2|$ and their maxima, $\max |{\omega }_1-{\omega }_2|$, are plotted, to indicate when the real-time coupling strength $\tilde{g}(\tau )$ slowly increases (forward) and then decreases (backward). In the first-order transition regime, whatever direction $\tilde{g}(\tau )$ moves, the synchronization emerges or disappears at the same ${\tilde{g}}_{\mathrm{c}}$ as derived in the static analysis [see Fig. 4(a)]. The maximal frequency differences in the forward and backward evolutions are identical to those depicted in Fig. 4(c), coinciding with Fig. 3(a). In the second-order transition regime, although the synchronization also emerges at ${\tilde{g}}_{\mathrm{c}}$ in the forward trip, it does not disappear at the same critical point in the backward trip [see Fig. 4(b)]. Actually, the synchronization survives far below ${\tilde{g}}_{\mathrm{c}}$, even into the statically inaccessible region depicted in Fig. 3(c). Further calculation of the maximal frequency differences reveals a hysteresis loop in this case [see Fig. 4(d)]. The forward half of the loop remains the same as the curve in Fig. 3(b), whereas the backward half is beneath it. The critical coupling strength under the static model does not apply under a dynamical model with the second-order transition. As explained in the Supplementary Material,⁴⁷ the $\tilde{g}(\tau )$ passes ${\tilde{g}}_{\mathrm{c}}$ with the emergence of new non-zero eigenvalues proportional to $-\dot{\tilde{g}}{\mathbf{J}}^{-1}$. This mechanism is thus attributed to the singularity of ${\mathbf{J}}^{-1}$ and the altering direction of real-time $\tilde{g}(\tau )$. The hysteresis property breaks the minimal coupling required for the synchronization, and the consequent temporal nonreciprocity enables the reading out of coupling history as has been done in the ferromagnetic materials.⁴⁸

Fig. 4

Hysteresis behavior in frequency difference. (a), (b) Frequency differences $|{\omega }_1-{\omega }_2|$ versus the evolution time $\tau$ in the first- and second-order transition regimes. Insets: the real-time evolution of the coupling strength $\tilde{g}(\tau )$. (c), (d) Maxima of the frequency differences, $\max |{\omega }_1-{\omega }_2|$ versus $\tilde{g}(\tau )$. For each plot, the Kerr factor $\tilde{\delta }=0.1$; the detuning $\tilde{\mathrm{\Delta }}=0.2$ in (a) and (c), and $\tilde{\mathrm{\Delta }}=0.3$ in (b) and (d).