Designs toward synchronization of optical limit cycles with coupled silicon photonic crystal microcavities

A driven high-Q Si microcavity is known to exhibit limit cycle oscillation originating from carrierinduced and thermo-optic nonlinearities. We propose a novel nanophotonic device to realize synchronized optical limit cycle oscillations with coupled silicon (Si) photonic crystal (PhC) microcavities. Here, coupled limit cycle oscillators are realized by using coherently coupled Si PhC microcavities. By simulating coupled-mode equations, we theoretically demonstrate mutual synchronization (entrainment) of two limit cycles induced by coherent coupling. Furthermore, we interpret the numerically simulated synchronization in the framework of phase description. Since our proposed design is perfectly compatible with current silicon photonics fabrication processes, the synchronization of optical limit cycle oscillations will be implemented in future silicon photonic circuits.


INTRODUCTION
Synchronization is a universally observed phenomenon in nature [1]. In fact, the observation of synchronization has a long history, which may go back to the 17th century with Huygens's discovery of synchronization of two pendulum clocks. In the 19th century, Lord Rayleigh reported the unison of sounds in acoustical systems. The first modern experimental studies of synchronization were performed by Appleton and van der Pol in the early 20th century using electrical and radio engineering techniques [2,3]. On the other hand, for a clear understanding of synchronization, we had to wait until the late 20th century, when phase description of limit cycles was developed by Winfree and Kuramoto [4,5]. Limit cycle oscillation emerges from a nonlinear dissipative system and well models various rhythm and self-pulsing phenomena. Since limit cycles have stable orbits, they are different from harmonic oscillations in conservative systems. The main idea of phase description is to describe limit cycle dynamics solely with a (generalized) phase degree of freedom. The phase description was found to be a powerful tool for understanding not only single limit cycle dynamics but also synchronization phenomena. In fact, for an intuitive understanding of mutual synchronization (entrainment) of coupled limit cycles, phase description provides a powerful tool called the phase coupling function. Furthermore, phase description is not limited to two oscillators, and it can also be used to analyze an ensemble of coupled oscillators, which is called the Kuramoto model. Nowadays, the phase analysis of synchronization is an indispensable tool to understand various synchronization phenomena in physics, chemistry, biology, and physiology. In biology, the numerous examples of synchronization phenomena range from the circadian rhythm to firefly synchronization [1]. In physics, synchronization phenomena in several systems has only recently been discussed. The most famous example may be the Josephson junction array, which is known to be described by the Kuramoto model [6][7][8]. In photonic systems, synchronization has been demonstrated with coupled lasers, microcavity polaritons, optomechanical oscillators, and trapped ions [9][10][11][12][13][14][15][16][17][18]. Furthermore, very recently, a frequency comb was interpreted in terms of synchronization [19].
In this paper, we propose a novel nanophotonic system with standard silicon (Si) photonic crystal (PhC) technologies that realizes synchronization of optical limit cycles. In our previous paper [20], we experimentally investigated the detailed properties of stochastic limit cycle oscillation (self-pulsing) in a single driven high-Q Si PhC microcavity. Here, we extended the previous study to coupled driven Si PhC microcavities. First, by numerically simulating coupled-mode equations, we demonstrate that introducing coherent field coupling between two cavities gives rise to synchronization (entrainment) of two limit cycle oscillations. Interestingly, we found that the synchronization phase (for example, in-and anti-phase synchronizations) can be controlled by the phase difference between two laser inputs. Second, we qualitatively interpreted the numerically demonstrated synchronization in the framework of the phase description (phase reduction theory). For this purpose, we calculated the phase coupling function, which plays a central role in phase description [5,21]. The obtained phase coupling function intuitively explains the origin of the synchronization and the synchronization phase. Finally, we demonstrated synchronization in a realistic coupled cavity device, which has moderately different cavity resonance frequencies.
PhC cavity structures largely enhances carrier-induced and thermo-optic optical nonlinearities with their very high-Q value and nanoscale mode-volume [22][23][24]. Employing the enhanced carrier-induced and thermo-optic nonlinearities in high-Q PhC cavities, optical bistability [25][26][27][28][29], limit cycle oscillation [30][31][32], and excitability [31,32] were demonstrated. Furthermore, recently, coupled PhC cavities has been actively investigated to realize, for instance, slow-light [33], Fano resonance [34,35], unconventional photon blockade [36], and self-pulsing coupled nanolasers [37][38][39]. Here, Si PhC cavities are advantageous also for studying synchronization of optical limit cycles from the standpoint of measurements and their controllability. In particular, for measurements, the real-time dynamics of light outputs are easily obtained with conventional optical setups. Meanwhile, the input power and frequency of a driving laser are easily controlled. Furthermore, since the proposed coupled Si PhC cavity device does not require any active material, and is based on the standard Si fabrication technique, its integration with other Si photonic devices is easy. Thus, it will be easy to implement the demonstrated limit cycle synchronization for future silicon photonic information processing and optical communications [40]. Ultimately, an array of Si PhC cavities will work as a one-dimensional nearest-neighbor coupled Kuramoto model.

LIMIT CYCLE IN A SINGLE HIGH-Q SI PHC CAVITY
First, we review limit cycle oscillation emerging from a single high-Q Si PhC cavity, which we investigated in our previous paper [20]. We consider a single Si L3-type cavity with two waveguides as schematically shown in Fig. 1(a), which is the same as in Ref [20]. The PhC slab is a two-dimensional hexagonal lattice, and the cavity is introduced by removing three air-holes. Note that, in the sample used in Ref. [20], several air-holes around the cavity region were carefully modulated to achieve larger Q value than that of the conventional L3 cavity [41]. The cavity, which has resonance frequency ω c , is driven by a laser input with frequency ω L and power P through the input waveguide. When input power exceeds a critical value, the output light exhibits limit cycle oscillation (self-pulsing) originating from nonlinear field, Detuning δ=ω Detuning δ=ω FIG. 1. (a) Schematic of a high-Q Si PhC microcavity with two waveguides. (b) Self-pulsing (SP) and bistable (BS) regions as functions of laser input power P and detuning δ(= ωL − ωc). (c) Input power P and detuning δ dependence of limit cycle oscillation frequency Ω. In (b) and (c), blue and red filled circles represent parameters used for the cavity C1 and C2 in Fig. 2, respectively. (d), (e) Time evolutions of the light output intensity |α(t)| 2 (left), carrier n(t) (right), and thermal component θ(t) (right) for P = 0.6κ 2 (d) and 1.0κ 2 (e). In (d) and (e), we used δ = −2κ.
carrier, and thermal dynamics. Now, we write up the coupled-mode equations describing field, carrier, and thermal dynamics in the nonlinear Si PhC cavity, which were proposed in Ref. [42,43] and also used in our previous paper [20]. Electric field α, normalized carrier density n, and normalized thermal effect θ follow the coupled-mode equationṡ where the detuning δ is defined as δ = ω L − ω c . The thermal effect θ is proportional to a temperature difference between the internal and external regions of the cavity. It is important to note that the variables n and θ are normalized so that the nonlinear coefficients before n and θ in Eq. (1) are unity. The nonlinear coefficients f , ξ, β, and η represent free-carrier absorption (FCA), two-photon absorption (TPA), heating with linear photon absorption, and FCA-induced heating, respectively. The small Kerr nonlinearity is neglected in the coupled-mode equations. In the rest of this paper, we use f = 0.0244, ξ = 8.2, β = 0.0296, and η = 0.0036, which are the same as in Ref [43]. Although a precise determination of the values of the nonlinear coefficients is difficult, exact values are not necessary, and the qualitative reproduction of the observed limit cycle oscillation is sufficient. For the lifetimes of the three variables, we set 1/2κ = 300 ps (Q ∼ 3.5 × 10 5 ), 1/γ = 200 ps, and 1/Γ = 100 ns. As discussed in Ref. [25,44], in the L3-type PhC cavity, due to the small cavity region, fast carrier diffusion makes the carrier lifetime comparable to the field lifetime. The details of our model are described in the Supplemental Material in Ref. [20].
Here, we briefly discuss the steady-state properties of coupled-mode equations (1)-(3). Here, α ss , n ss , and θ ss represent the steady state values of the field, carrier, and thermal effect, respectively. By settingα = 0,ṅ = 0, andθ = 0 in Eqs (1)-(3), an algebraic equation for I ss = |α ss | 2 is obtained as Depending on input power P and detuning δ, the algebraic equation (4) has one or two solutions for I. We numerically solve Eq. (4) and obtain I ss . Using I ss , we respectively calculate n ss and θ ss as n ss = κ γ ξI 2 ss and θ ss = Using n ss and θ ss , we can write the complex electric field α ss as Second, to check the stabilities of the steady states, we perform a linear stability analysis for coupled-mode equations (1)- (3). For this purpose, decomposing the complex field α as α = x + iy, we rewrite Eqs. (1)-(3) aṡ where the vector x is defined as x = (x, y, n, θ). Now, a 4×4 Jacobian matrix corresponding to the dynamical system in Eq. (7) is given by Now, a small fluctuation δx followsδx J δx, where δx ≡ x − x ss with the steady state values x ss = (x ss , y ss , n ss , θ ss ). We calculate the eigenvalues of J (x) at the steady states x = x ss for various input power P and detuning δ. When the pair of the eigenvalues of the Jacobian J (x) have positive real values, the steady state x ss becomes unstable, which leads to limit cycle oscillation (the Hopf bifurcation) [5,45]. We show nontrivial regions as functions of input power P and detuning δ in Fig. 1(b), where the bistable and limit cycle (self-pulsing) region are indicated by "BS" and "SP", respectively. In the SP+BS region, one steady state is stable, while the other is unstable. In this paper, since we are interested in limit cycle oscillation, we focus solely on the SP region. We also note that the Jacobian matrix Eq. (8) is used again in Section 4. Additionally, in Fig. 1(c), we show the input power P and detuning δ dependence of the limit cycle's frequency Ω, which were obtained from numerical time evolutions. Figure 1(c) indicates that the limit cycle's frequency decreases with increasing pump power P . Now, we directly simulate coupled-mode equations (1)-(3). The real-time evolutions of light output I(t) = |α(t)| 2 (left), carrier n(t) (right), and θ(t) (right) are shown in Fig. 1(d), where the detuning is δ = −2κ, and laser input powers are P = 0.6κ 2 (d) and 1.0κ 2 (e). In Fig. 1(d), which is for P = 0.6κ 2 , all the variables reach steady states when t 1000 ns, and there is no self-pulsing. Meanwhile, for P = 1.0κ 2 [see Fig. 1(e)], all the variables clearly exhibit temporal periodic oscillations (limit cycle oscillation) with a frequency of Ω/2π = 11 MHz. In fact, in Fig. 1(b), the values δ = −2κ and P = 1.0κ 2 are represented as a blue filled circle in the SP region. Meanwhile, the values δ = −2κ and P = 0.6κ 2 are outside the SP region . In the rest of this paper, we show only light output I(t) = |α(t)| 2 because it is the only measurable valuable in experiments.
Finally, we comment on the origin of limit cycle oscillation in Si PhC microcavities. In a minimum model that exhibits limit cycle oscillation, we set η = 0 and f = 0 in Eqs. (1) and (3), which have only quantitative effects. Furthermore, the exponent of the term κξ|α| 4 in Eq. (2) is not essential, because limit cycle oscillation appears even if this term is replaced with κξ|α| 2 . In fact, limit cycle oscillation requires only that the signs of the nonlinear energy shifts be opposite for carrier and thermal components; that the carrier lifetime be comparable to or even shorter than the photon lifetime, with the thermal lifetime much longer than the photon lifetime; and that β be much smaller than ξ, approximately β/ξ γ θ /γ n , to make carrier-and thermal-induced energy-shifts comparable. Due to the large time-scale difference and the opposite sign of the nonlinear energy shifts, a delayed positive feedback instantaneously occurs when the effective cavity frequency returns to the frequency of the laser input, which leads to self-pulsing.

COUPLED LIMIT CYCLE DYNAMICS
. The proposed device with two coupled Si PhC cavities is sketched in Fig. 2, where the two cavities are labelled as C1 and C2. Since the cavities are evanescently coupled, coupling strength g depends on the distance between the two cavities. To drive the two cavities, we separate a single laser source into two inputs using on-chip Si wire waveguides instead of two laser sources. This process is very important for temporally fixing the relative phase difference φ L between the two laser inputs. Actually, we show that the relative phase difference φ L has a crucial impact on synchronization. The design shown in Fig. 2 has two output waveguides, which are used to measure light outputs from C1 and C2. The coherent coupling g is introduced through tunneling of evanescent fields, which is linear coupling. The light outputs from the two cavities are extracted from the two output waveguides, while the two cavities are driven though the two input waveguides. The two laser inputs from a single laser source are separated with on-chip Si wire waveguides and a beam splitter. The phase difference of the laser inputs φL is adjusted by the optical path lengths of the on-chip Si wire waveguides.
coupled cavities aṡ Equations (9)-(12) and Eqs. (12)-(14) represent dynamics for C1 and C2, respectively. The coherent field coupling between C1 and C2 is represented by the coupling strength g. In Eq. (12), the term e iφ L represents a phase factor originating from the phase difference between the two laser inputs. It is worth noting that φ L is the phase associated with the field, and thus is not directly related to a limit cycle's phase, which is introduced in Section 4. The values of the nonlinear coefficients f , ξ, β, and η are the same as those in Fig. 1. The cavity detuning is defined as δ 1,2 ≡ ω L − ω 1,2 , where ω 1,2 is the resonance frequency of the cavity. To observe synchronization, there must to be a small frequency difference in two limit cycles. However, in our proposal, the two cavities are designed to be identical because natural disorders or unavoidable fabrication errors will introduce an intrinsic parameter and resonance frequency difference between the two cavities. In this section, for the demonstration of synchronization, we consider a rather ideal device. Namely, only the cavity resonance frequencies are slightly different: The other parameters are the same for the cavity C1 and C2: 1/2κ 1 = 1/2κ 2 = 300 ps, 1/γ 1 = 1/γ 2 = 200 ps, and 1/Γ 1 = 1/Γ 2 = 100 ns. Additionally, we derive the two cavities with the same input powers, P 1 = P 2 = κ 1 . In Section 5, we consider a more realistic device, where the resonance frequencies of the two cavities are moderately different. Figure 3(a) shows time evolutions of the light output |α 1,2 (t)| 2 without g = 0 (top) and with coupling g = 0.02κ 1 (middle and bottom), which are the central results of this paper. In Fig. 3 With coupling (g = 0.02κ 1 and ΦL = 0) With coupling (g = 0.02κ 1 and ΦL = �) Simulated time evolutions of the light output intensity |α1,2(t)| 2 without (top) and with coupling g = 0.02κ1 (middle and bottom). The middle and bottom time evolutions are for φL = 0 and π, respectively. (b) The average frequenciesΩ1,2 of the two limit cycle oscillations |α1,2(t)| 2 for φL = 0 (upper) and π (lower) as a function of the coupling strength g. The critical coupling strengths of synchronization are gc = 0.0115 and 0.011 for φL = 0 and π, respectively. In these simulations, we used 1/2κ1 = 1/2κ2 = 300 ps, 1/γ1 = 1/γ2 = 200 ps, 1/Γ1 = 1/Γ2 = 100 ns, φ L = π, respectively. Note that this coupling strength (g = 0.02κ 1 ) is much smaller than the cavity decay rates (g κ 1,2 ), so the two coupled cavities are in the weak-coupling regime. Without coherent coupling (g = 0) [ Fig. 3(a), top], the two limit cycle oscillations are completely decoupled and thus have their own frequencies: Ω 1 /2π = 11 and Ω 2 /2π = 11.75 MHz for C1 and C2, respectively. On the other hand, when small coupling g = 0.02κ 1 is introduced, the time evolution dramatically changes as shown in the middle and bottom time evolutions in 3(a), where the two limit cycle oscillations are perfectly synchronized (entrained) with each other. Furthermore, we notice that synchronization is in-phase for φ L = 0 (middle), while "anti-phase" for φ L = π (bottom). Now, we briefly discuss a synchronization time. By turning on coupling for "steady-state" uncoupled limit cycle oscillations (not shown), we found that the synchronization time for g = 0.02κ 1 is about 500 ns, which corresponds to approximately the five periods of the limit cycle oscillations. In fact, while the oscillation frequency of limit cycles is typically about ∼10 MHz in Fig. 3(a), the strength of coherent coupling g = 0.02κ 1 (1/g = 1/0.02κ 1 = 30 ns) corresponds to 5.3 MHz. We also found that when the coupling strength is increased to g = 0.05κ 1 , the synchronization time actually becomes comparable to the period of the limit cycles (not shown).
In addition to the time evolutions, we show in Fig. 3(b) the mean frequenciesΩ 1,2 of the two limit cycles as a function of coupling strength g, where the upper and lower graphs are for φ L = 0 and π, respectively. We used the mean frequencies because the oscillations are not perfectly periodic when the coupling strength is smaller than a critical value for synchronization. Figure 3(b) clearly shows that as the coupling strength g increases, the mean frequenciesΩ 1 andΩ 2 approach each other and merge when g reaches the critical value g c : g c = 0.0115κ 1 and 0.011κ 1 for φ L = 0 and π, respectively. Furthermore, for both φ L = 0 and π, the frequencies of the synchronized limit cycles are the same: Ω 1 = Ω 2 , which is called 1:1 synchronization. Finally, we comment on the fact that the value of g c is not the same for φ L = 0 and π. We found that when the frequency difference between uncoupled limit cycles becomes smaller, synchronization occurs with the same critical values of g c for both φ L = 0 and π, respectively.
In Appendix A, we show a simulation of an intermediate phase difference φ L = 0.5π, which does not exhibit synchronization with coupling strength g = 0.02κ 1 . Furthermore, we discuss synchronization in a near-strong coupling region (g = κ 1,2 ) in Appendix B. Compared to the very small coupling g = 0.02κ 1 considered in this section, the near-strong coupling region may be technically easy to realize.

PHASE DESCRIPTION
In Section 3, we demonstrated synchronization of limit cycle oscillations in two cavities by directly simulating time evolutions. For a qualitative understanding of the synchronization, phase reduction theory provides a powerful tool called the phase coupling function [5,21,46]. In particular, the phase coupling function can explain why the in-or anti-phase synchronization occurs depending on the phase difference between the two laser inputs. In this section, after a brief review of phase reduction theory, we numerically derive the phase equation of motion and phase coupling function for coupled-mode equations (9)-(14).

General phase description for a single limit cycle
The key idea in phase reduction theory is to describe limit cycle dynamics solely with a generalized phase degree of freedom. First, we consider the phase description for general single limit cycle dynamics and introduce a scalar "phase field" φ(x). Let us consider a general dynamical system that exhibits limit cycle oscillation: where f (x) is a general function. In the phase description, the phase field φ(x) is defined in such a way thatφ where Ω is the frequency of the limit cycle oscillation. If there is no perturbation, dynamics converge on the orbit of the limit cycle and follow the very simple equation of motionφ = Ω, where φ without any argument represents the phase on the limit cycle's orbit. For simplicity, we denote its orbit as χ(φ). When the dynamical system [Eq. (15)] is perturbed by a force If the perturbation p(x) is sufficiently weak, x is approximated as a point on the limit cycle's orbit, x χ(φ). With this approximation, Eq. (17) is further simplified aṡ where Z(φ) ≡ ∇ x=χ(φ) φ(x) is called "sensitivity" [4]. Here the capital P (φ, t) is defined as P (φ, t) ≡ p(χ(φ), t). Equation (18) is called the phase equation of motion, and it plays a central role in phase reduction theory. Actually, with Eq. (18), the perturbed limit cycle dynamics are described solely by the phase degree of freedom φ. Therefore, our next step is to numerically determine the sensitivity Z(φ) for our dynamical system described by coupled-mode equations (1)- (3). Fortunately, to numerically obtain Z(φ), we can use the adjoint method [46,47], which employs the fact that Z(φ) satisfies the following equation of motion: where J (χ(Ωt)) is the transpose of the Jacobian of a dynamical system. In our case, the Jacobian matrix J is already given in Eq. (8). Since Eq. (19) is unstable for forward time integration due to the minus sign before J , we need to perform backward time integration as dZ(−Ωt )/dt = J (χ(−Ωt ))Z(−Ωt ) with t = −t. Additionally, the numerically obtained Z(φ) was normalized as Z(φ) · f (χ(φ)) = Ω, which is equivalent to Eq. (16). Figure 4(a) shows numerically obtained Z i (φ), where the index i represents x, y, n, and θ. Importantly, parameters used for calculating Z(φ) are the same as those in Fig. 1(e). From Fig. 4(a), we notice that there is a scale difference between the four components and the shape of Z(φ) is very complicated compared with, for example, the sensitivity of the simple Stuart-Landau model [46].

Phase coupling function
Here, we extend phase equation of motion (18) to two coupled limit cycles. Let us consider two weakly coupled dynamical systems, both of which exhibit limit cycle oscillations: where δf 1,2 (x 1,2 ) is a deviation from the "standard" oscillator f (x) [Eq. (15)], while g 12 (x 1 , x 2 ) and g 21 (x 2 , x 1 ) represent coupling between the two systems. Rewriting with the phase coordinate φ 1,2 of the standard oscillator and taking the terms δf i (x i ), g 12 (x 1 , x 2 ), and g 21 (x 2 , x 1 ) as perturbations, the phase equations of motion corresponding to Eqs (20) and (21) are given byφ where the upper-case symbols represent the functions of the standard oscillator's phase φ 1,2 , which is given byẋ 1,2 = f (x 1,2 ). For further simplification of Eqs. (22) and (23), we transform φ 1,2 into the rotating frame of the standard oscillator as ψ 1,2 ≡ φ 1,2 − Ωt, where Ω is the standard oscillator' oscillation frequency. Additionally, we perform an approximation for the coupled phase equations of motion by averaging over one period of the standard oscillator. With these procedures, Eqs. (22) and (23) becomė Here, the frequency shift δΩ 1,2 and the phase coupling function Γ ij (ψ) are given by and respectively. Finally, the phase difference between the two oscillators, ψ = ψ 2 − ψ 1 , follows the following simple equation:ψ where ∆Ω ≡ δΩ 2 − δΩ 1 and Γ a (ψ) ≡ Γ 21 (ψ) − Γ 12 (−ψ). In fact, Γ a (ψ) is the anti-symmetric part of the phase coupling function. A synchronization phase ψ sync is required to satisfy Γ a (ψ sync ) = 0 and Γ a (ψ sync ) < 0, where the prime represents the derivative. For example, if Γ a (0) = 0 and Γ a (0) < 0, the phase difference ψ is locked to ψ = 0 by negative feedback, which is in-phase synchronization. Therefore, the shapes of the phase coupling function allow an intuitive interpretation of a synchronization phase.
The purpose of this phase rotation is to define a common standard oscillator and its common phase φ for the two limit cycles. In fact, except for the coupling terms, Eqs (29)-(31) and Eqs (32)- (34) are the same equations of motion represented as f (x) [Eq. (7)]. Meanwhile, the coupling function g 12 (x 1 , x 2 ) and g 21 (x 2 , x 1 ) are given by and respectively. Additionally, for simplicity, we use the limit cycle in the cavity C1 as a standard oscillator, and thus we put δΩ 1 = 0. Since the parameters for the standard oscillator are the same as those in Fig. 1(e), we can use the sensitivity Z shown in Fig. 4(a). Representing the coupling function g ij (x j ) with the standard oscillator's phase coordinate as G ij (φ j ), we numerically integrate Eq. (27). Figure 4(b) and (c) show the anti-symmetric parts of the phase coupling function Γ a (ψ) ≡ Γ 21 (ψ) − Γ 12 (−ψ) for φ L = 0 and π, respectively. Here, the power of the phase description is that the complex limit cycle dynamics represented by coupled-mode equations (9)-(14) are reduced to a simple phase equation of motion (28). In fact, the origin of synchronization is understood only in this phase coordinate. Figure 4(b) and (c) clearly indicate that when φ L = 0 (b), Γ a (0) = 0 and Γ a (0) < 0 hold, and thus in-phase locking occurs. Meanwhile when φ L = π (c), Γ a (π) = 0 and Γ a (π) < 0 hold, and thus anti-phase locking occurs. Here, Fig. 4(c) is a mirror image of Fig. 4(b) about the x-axis, which is intuitive because the signs of Eqs (35) and (36) are opposite for φ L = 0 and π. In our case, since the phase coupling function for φ L = 0 [see Fig. 4(b)] resembles the sine function, in-and anti-phase synchronizations will occur for φ L = 0 and π, respectively. The surprise is that although the two cavities are in the weak-coupling regime (g κ 1,2 ), the phase φ L in Eqs (35) and (36) strongly modifies synchronization behavior. In fact, since coherent coupling between fields has a (relative) phase degree of freedom, in a coupled-cavity system, it is always important to take the phase into account.
The parameters used for calculating phase coupling functions in Fig 4(b) and (c) are again the same as those in Fig. 1(e). Note that the shape of the phase coupling function depends on parameters used for phase reduction. For example, if we perform phase reduction with the parameters for the limit cycle in C2 shown in Fig. 3(a), although the detailed shape of the phase coupling function changes from that in Fig. 4(b) (not shown), both have qualitatively the same quasi-sinusoidal shapes.
Thus, anti-phase synchronization for φ L = π is not a general result, which depends on models and parameters. Meanwhile, we found that in-phase synchronization for φ L = 0 seems to be general.
Finally, we discuss why the linear coherent coupling g ij (x j ) gives rise to the nonlinear phase coupling function Γ 1,j (φ j ) shown in Fig. 4. The mathematical answer is the transformation of the coordinate from the Cartesian coordinates x into the phase coordinate of the limit cycle φ. Namely, on the phase coordinate, the linear coupling g ij (x j ) appears as a nonlinear ) as a function of the phase ψ ≡ ψ2 − ψ1 for φL = 0 (left) and π (right). The arrows indicates phase locking points. To calculate Z(φ) and Γa(ψ), we used the same parameter values as those in Fig. 1 (e).
function G ij (φ j ). Since limit cycle oscillation itself originates in a nonlinear dissipative system, the transformation from the Cartesian to the phase coordinate is also nonlinear. We can also interpret our synchronization phenomenon as analogous to injection locking [48] or mutual injection locking [49] in laser physics. In injection locking, coupling between slave and master lasers is usually provided by partially transmitting mirrors, which is definitely linear coupling. Therefore, although the coupling itself is linear, synchronization occurs with the modulation of the slave laser's field by the master laser. Similarly to injection locking, in our system, the coherent coupling g allows the oscillating light in the cavity C1 to modulate the light in the cavity C2. Thus, synchronization is interpreted as a response of the limit cycle in the cavity C2 (C1) to the modulation from C1 (C2).

SYNCHRONIZATION OF TWO MODERATELY DIFFERENT LIMIT CYCLES
Until now, we have considered synchronization in rather ideal systems, where the two limit cycles are almost identical and only their cavity resonance frequencies are slightly different. Thus, it is still questionable whether or not realistic Si PhC cavity devices are able to exhibit synchronization of limit cycle oscillations. Even with state-of-the-art fabrication technology, fabrication errors or natural disorders cause, for example, unavoidable resonance frequency differences in cavities. Therefore, in this section, we consider a more realistic device, where two cavities have a moderate resonance frequency difference.
First, we set photon lifetimes for the two cavities as 1/2κ 1 = 1/2κ 2 = 100 ps, which correspond to Q ∼ 1.0 × 10 5 . Importantly, compared with the simulations in previous sections, we slightly decreased the photon lifetime from 300 to 100 ps. This is because it is technically easier to reduce the difference in cavity resonance frequencies for a shorter photon lifetime (a lower Q value). We use the same values as in previous sections for the nonlinear coefficients: f = 0.0244, ξ = 8.2, β = 0.0296, and η = 0.0036. For these parameters, the SP and BS regions are represented by the diagram shown in Fig. 5(a). Additionally, we show the detuning and input power dependence of the limit cycles' frequency Ω in Fig. 5(b), which is more complicated than Fig. 1(c). In fact, the oscillation frequency does not monotonically decrease with increasing input power, because there is an increase in the oscillation frequency at P 2.6κ 2 , and this jump might be related to the onset of fast photon-carrier oscillation [30]. Second, we introduce a moderate difference to the cavity resonance frequencies as ω 2 − ω 1 = 7κ 1 . Finally, we also set the value of the coupling strength as g = 0.4κ 1 , which is (d) Without coupling (g = 0) With coupling (g = 0.4κ and ΦL = 0) Input laser power P [κ 1 2 ] Laser frequency Detuning δ=ω Input laser power P [κ 2 ] Detuning δ=ω FIG. 5. Synchronization in realistic cavities with photon lifetime 1/2κ1 = 1/2κ2 = 100 ps and moderate resonance frequency difference ω2 − ω1 = 7κ1. (a) Self-pulsing (SP) and bistable (BS) regions as functions of laser input power P and detuning δ = ωL−ωc. (b) Input power P and detuning δ dependence of the limit cycle's frequency Ω. In (a) and (b), the blue and red filled circles represent parameters used for the cavity C1 and C2, respectively. (c) Transmission spectra of the coupled cavities obtained as steady state light output intensities |α1,2(ωL)| 2 as a function of laser input frequency ωL. The spectra were obtained with the laser input power fixed as P1 = P2 = 0.001κ 2 1 . (d) Simulated time evolutions of the light output intensity |α1,2(t)| 2 without g = 0 (upper) and with coupling g = 0.02κ1 (lower). To simulate the time evolutions, we used P1 = P2 = 9κ 2 1 , δ1 ≡ ωL − ω1 = −1.0κ1, and δ2 ≡ ωL − ω2 = −8.0κ1. much stronger than in Section 3.
We show the spectra of the two cavities in Fig. 5(c), which was obtained by sweeping the laser frequency ω L from ω 1 − 20κ 1 to ω 1 + 20κ 1 and plotting the steady state outputs |α 1 | 2 and |α 2 | 2 with very low input power P 1 = P 2 = 0.001κ 2 1 so as not to induce any nonlinearity. In Fig. 5(c), the dashed curves are the spectra without coupling, g = 0; the solid blue and red curves are the spectra with coupling, g = 0.4κ 1 . Comparing the spectra with and without coupling, we notice that the moderately large coupling strength (g = 0.4κ) induces the signatures of coupling as peaks [see the two arrows in Fig. 5(c)] in the spectral tails, but does give rise to normal-mode splitting. Therefore, the system is still in the weakcoupling regime, and we are able to consider coupling as perturbation. Here, though the coupling is moderately strong, the system is in weak-coupling because of the large frequency difference between the two cavities ∆ω ≡ ω 2 − ω 1 = 7κ 1 . Recall that strong-coupling requires not only g > κ 1 but also g > |∆ω|. Note also that this value of the resonance frequency difference ∆ω = 7κ 1 is experimentally available with state-of-the-art fabrication technology [50,51]. Additionally, in Appendix C, we briefly discuss the configuration of two PhC cavities to realize the coupling strength g = 0.4κ 1 with finite-difference time-domain (FDTD) simulations.
To drive the cavities, we set the detuning values between the cavity resonance and laser frequency as δ 1 = ω 1 −ω L = −1.0κ 1 , which leads to δ 2 = ω 2 −ω L = −8.0κ 1 . Both cavities are driven by inputs with the same power P 1 = P 2 = 9.0κ 2 1 . These parameters are represented by the blue (C1) and red (C2) filled circles in the diagram in Fig. 5(a) and (b), which indicate that both cavities exhibit self-pulsing (limit cycle oscillation). Now, in the same way as in Fig. 3(a), we show the time evolution of light output intensity |α 1,2 (t)| 2 in Fig. 5(d) with (upper) and without coupling (lower). First, we discuss time evolution without coupling, g = 0. Without coupling, both cavities exhibit limit cycle oscillations with their own frequencies: Ω 1 /2π = 5.8 and Ω 2 /2π = 9.2 MHz for the cavity C1 and C2, respectively. On the other hand, with coupling, g = 0.4κ 1 , the time evolution clearly shows synchronization of the two oscillations. However, the profile of synchronized oscillations is very different from that in Fig. 3(a). For instance, the profile of |α 2 (t)| 2 is strongly modified by introducing the large coupling.
In the synchronized state [see the lower time evolutions in Fig. 5(d)], the oscillation periods of the limit cycle oscillations for |α 1 (t)| 2 and |α 2 (t)| 2 are identified as T 1 = T 2 = 210 ns, which corresponds to Ω 1 /2π = Ω 2 /2π = 4.77 MHz. However, from |α 2 (t)| 2 on the lower panel in Fig. 5(d), we notice that the limit cycle orbit for C2 is strongly modified by coherent coupling compared with the uncoupled orbit. In fact, in terms of the Poincaré section, the period of the limit cycle orbit |α 2 (t)| 2 with coupling will be T 2 105 ns and the corresponding frequency is Ω 2 /2π 9.54 MHz, which is the double of 4.77 MHz. This co-existence of the two frequencies and temporal profile for |α 2 (t)| 2 [see Fig. 5(d)] may be signatures of period doubling bifurcation [45], which is an interesting theme for future investigation both from the theoretical and experimental standpoints.
Finally, we comment on the phase difference of the two laser inputs, φ L . In this section, we have shown the simulation only for φ L = 0 because we found that synchronization occurs only for near-zero phase φ L 0. Actually, when φ L is not close to zero, synchronization does not occur even with g > κ. This result is related to the large frequency difference between the two uncoupled limit cycles (Ω 1 = 5.8 and Ω 2 /2π = 9.2 MHz). In fact, if the parameters for the cavity C1 and C2 are similar , the frequency difference between two uncoupled limit cycles is smaller, and synchronization occurs both for φ L = 0 and π. Therefore, to realize synchronization in realistic coupled cavities with a moderate frequency difference, it is important to adjust the phase difference of laser inputs to near-zero (φ L = 0), which will be achieved by adjusting optical path lengths with, for example, on-chip Si wire waveguides.

DISCUSSION AND FUTURE PERSPECTIVE
First, we argue that the proposed scheme of synchronization is not limited to Si PhC cavities, but applicable to a wide range of limit cycle oscillations in nanophotonic systems, such as nanolasers [37,39], semiconductor microcavities [31,32], and microring resonators [42,43,[52][53][54]. Actually, the coherent field coupling is easily implemented in these nanophotonic devices, which will lead to synchronization of optical limit cycles. In particular, since coupled-mode equations (1)-(3) were originally proposed for modelling optical limit cycles in Si microring resonators [42,43], our synchronization scheme is easily applicable to them. In terms of the tunability of various physical parameters such as resonance frequencies, Si microring resonators may be advantageous over PhC structures. In particular, a different type synchronization dynamics was investigated with coupled microring resonators in Ref [55]. Second, we discuss a future perspective of limit cycle synchronization in Si PhC cavities. One can naturally imagine the extension of the two coupled Si PhC cavities to an array of coupled cavities as illustrated in Fig. 6. In principle, the coupled PhC cavity array illustrated in Fig. 6 could behave as a one-dimensional (1D) nearest-neighbor coupling (local) Kuramoto oscillator. The 1D local Kuramoto model has been theoretically investigated by numerical simulation [56] and renormalization group analysis [57], which have predicted various nontrivial collective phenomena, including a synchronization state, a phase slip at the onset of de-synchronization, and coupling-induced chaos. From the standpoint of device application, the predicted chaotic state in the 1D local Kuramoto model could be used for photonic reservoir computing [58].

CONCLUSION
In conclusion, we have theoretically demonstrated synchronization of optical limit cycles with driven coupled silicon (Si) photonic crystal (PhC) cavities, where limit cycle oscillation emerges from carrier-and thermal-induced nonlinearities. Introducing coherent field coupling between two cavities synchronizes (entrains) two limit cycle oscillations. First, we quantitatively demonstrated synchronization by directly simulating the time evolutions of coupled-mode equations. We found that synchronization phase depends on the phase difference of two laser inputs. Second, the numerically simulated synchronization was qualitatively interpreted in the framework of phase description. In particular, we calculated phase coupling functions, which intuitively explain why the synchronization phase depends on the phase difference between the two laser inputs. Finally, we discussed synchronization in a realistic coupled cavity device, where the resonance frequencies of the two cavities are moderately different. Since our proposed design is perfectly compatible with conventional Si fabrication processes, synchronization of optical limit cycles will be easy to implement in future silicon photonic devices and can be extended to an array of coupled cavities.  Fig. 3(a).

DISCLOSURES
The authors declare no conflicts of interest.
APPENDIX A: SIMULATIONS FOR φL = 0.5π In this appendix, we show simulations when the phase difference in laser inputs is φ L = 0.5π in Fig. 3. The simulation described in Section 3 in the main text was performed only for φ L = 0 and π, which exhibited in-and anti-phase synchronization, respectively. Thus, it is of natural interest to discuss the intermediate case φ L = 0.5π. Figure 7(a) shows the time evolution of light output |α 1,2 (t)| 2 for φ L = 0.5π with coherent coupling g = 0.02κ 1 . In fact, in Fig. 7(a), all the parameters except φ L = are the same as in Fig. 3(a). Surprisingly, even though the coupling strength is the same as in Fig. 3(a), no synchronization is observed in 7(a). We found that this result can be explained in terms of a phase coupling function. Similarly to 4(b) and (c) , we show the anti-symmetric part of the phase coupling function in Fig. 7(b). Interestingly, Γ a (ψ) for φ L = 0.5π never crosses the zero axis, and thus there is no phase locking point. This explains why phase synchronization does not occur for φ L = 0.5π with the small coupling strength (g = 0.02κ 1 ).
Even for φ L = 0.5π, if the coupling strength is further increased, for example, to g 0.2κ 1 , synchronization occurs (not shown). However, this synchronization with a large coupling strength may not be interpreted as 1:1 synchronization, because, there is no smooth transition of the limit cycles' average frequencies from the independent to synchronized state for φ L = 0.5π. In summary, for φ L = 0.5π, 1:1 synchronization does not occur, but m:n synchronization can occur with a large value of coupling.

APPENDIX B: SYNCHRONIZATION IN STRONG COUPLING REGIONS
We demonstrate synchronization in the strong coupling region of two coupled cavities. For this purpose, we investigate a near-strong coupling region where g = κ 1 = κ 2 and g > |∆ω| with ∆ω ≡ ω 2 − ω 1 = 0.5κ 1 . Although synchronization with weak coherent coupling is theoretically interesting, the very weak coupling strength g = 0.02κ 1 for the high-Q cavities 1/2κ 1 = 1/2κ 2 = 300 ps used in Section 3 may not be easy in real PhC devices, it is important to discuss coupled limit cycles in the strong coupling region.
In this Appendix, for the detuning, we set δ 1 = ω L −ω 1 = 0 and δ 2 = ω L −ω 2 = −0.5κ 1 for the cavity C1 and C2, respectively. The other parameters except for the coupling strength and detuning are the same as those in Fig. 3. Thus, for laser input powers, we used P 1 = P 2 = κ 2 1 . Importantly, when there is no coupling (g = 0), no limit cycle oscillations appear with these parameters as indicated by the time evolutions shown in Fig. 8(b). In fact, the parameters for C1 and C2 are, respectively, indicated by the blue and red filled circles on the trivial region in Fig. 8(a).
On the other hand, when a near-strong coupling (g = κ 1 ) is introduced, there are still no limit cycle oscillations for φ L = 0, while synchronized limit cycle oscillations, surprisingly, appear for φ L = π. These phase φ L dependent results can be understood by considering the normal-modes of the fields α 1 and α 2 . When φ L = 0, the electric fields α 1 and α 2 form a "bonding" normal-mode. Since the frequency (energy) of the bonding normal-mode ω b is lower than the original cavity frequencies ω 1 and ω 2 (ω b < ω 1 , ω 2 ), the corresponding detuning δ a (≡ ω L −ω a ) is still outside the SP region and no limit cycle oscillations appear [see the green filled square in Fig. 8(a)]. Meanwhile, when φ L = π, the electric fields α 1 and α 2 form an "anti-bonding" normal-mode whose frequency (energy) ω a is higher than the original cavity frequencies: ω a > ω 1 , ω 2 . Therefore, the corresponding detuning δ a (≡ ω L − ω a ) is lower than δ 1,2 and enters the SP (self-pulsing) region [see the black filled square in Fig. 8(a)]. Note that to plot the green and black filled square in Fig. 8(a), we used the frequencies of the bonding ω b and anti-bonding normal-modes ω a given by and ω a = 1 2 respectively. We also note that, the temporal profile of the synchronized oscillations in Fig.  8(d) exhibits anti-phase synchronization, but is quantitatively different from those in Fig.  3(a) for φ L = π due to the large coupling strength g = κ 1 .
In conclusion, even in the strong-coupling region, synchronization can be realized, but we have to take the effect of normal-mode splitting into account. However, note that, the weakcoupling region is more interesting than the strong-coupling region from the standpoint of synchronization physics. This is because synchronization in the strong-coupling region can be interpreted simply as a limit cycle oscillation of a normal-mode that appears in both cavities.

APPENDIX C: ESTIMATING COHERENT COUPLING STRENGTH USING FDTD SIMULATIONS
Here, we briefly discuss the estimation of a coherent coupling strength g in a realistic coupled PhC cavity structure with three dimensional (3D) finite-difference time-domain (FDTD) simulations. The general idea of the FDTD simulation is as follows. We excite the cavity C1 at t = 0 and probe the time evolution of the energy of the electromagnetic field in the cavity C2. The time evolution of the electromagnetic field energy in the cavity C2 exhibits damped oscillation, where the oscillation originates from the coherent coupling (tunneling), Without coupling (g = 0) With coupling (g = 1κ 1 and ΦL = 0) With coupling (g = 1κ 1 and ΦL = �) Input laser power P [κ 2 ] t (ns) . 8. (a) Simulations for a near strong-coupling region (g = κ1 = κ2). For the cavity C1 and C2, we use δ1 = ωL − ω1 = 0 and δ2 = ωL − ω2 = −0.5κ1, respectively. The other parameters are the same as those in Fig. 3. (a) Self-pulsing (SP) and bistable (BS) regions as functions of laser input power P and detuning δ(= ωL −ωc). The blue and red filled circles represent parameters used for the cavity C1 and C2, respectively. Meanwhile the green and black filled squares represent the expected parameters for "bonding" (for φL = 0) and "anti-bonding" normal-modes (for φL = π) formed by the cavity fields, respectively. (b,c,d) Simulated time evolutions of the light output intensity |α1,2(t)| 2 without (b) and with coupling g = κ1 = κ2 (c,d). The phases are φL = 0 and π for (c) and (d), respectively.
while the damping is associated with the finite lifetimes of the cavities. Therefore, the oscillation frequency of the probed damped oscillation of the electromagnetic field energy in C2 corresponds to the coupling strength g. Note that as the distance between two cavities increases, the oscillation period becomes longer, which results in a very long computation time. In the 3D FDTD simulations in this Appendix, the maximum time range of time evolution was 1 ns, which already took a few days for computation. Figure 9(a) illustrates a coupled PhC cavity structure and simulated electric field distribution at t = 9 ps after pulse excitation to C1, while (b) represents the time evolutions of electric field energy probed in C1 and C2. Figure 9(b) indicates that the oscillation phases of the electric field energies in C1 and C2 are opposite, which is the evidence of the coherent energy transfer between the two cavities. Furthermore, from the oscillation period of the time evolution in Fig. 9(b), we can extract the coupling strength as 1/g = 6/2π [ps].
We investigated how the coupling strength changes depending on the configuration of two cavities denoted by ∆x and ∆y, which represent the x-and y-direction spacing between the two cavities measured by the number of air-holes, respectively [see the the bidirectional arrows in Fig. 9(a)]. For example, the configuration of the coupled PhC cavities shown in Fig. 9(a) is denoted as (∆x, ∆y) = (3,4). In Table I, we summarized the coupling strength for the five different configurations of two coupled cavities. For all 3D FDTD simulations, the total calculation area was 18.4 µm×18.4 µm×4 µm, while the lattice constant was 435 nm. The cavity C2 was on-resonantly excited with a pulse electric field whose pulse width is 3 ps. Depending on the configuration of the two cavities, their cavity lifetime vary from 100 to 200 ps, which must be at least longer than the period of the coherent oscillation (2π/g). Due to the limited computation time, the weakest coupling strength in Table I was 1/g = 380/2π 60 [ps] for the configuration (∆x, ∆y) = (7,7). Compared with the field decay rate 1/κ 1 = 200 ps assumed in Section 5, this coupling strength 1/g 60 is still   Fig.  9(a)]. Actually, 2π/g is the period of coherent oscillation of electromagnetic field energies.
stronger than the field decay rate: g > κ 1 . Although our FDTD simulations failed to calculate the weak-coupling region for 1/κ 1 = 200, Table I provides a hint of cavity configuration to realize the weak-coupling region. In fact, Table I indicates that the coupling strength seems to exponentially decreases with an increase in the distance between the cavities. Therefore, it is natural to expect that the coupling strength g = 0.4κ 1 with 1/κ 1 = 200, which is assumed in Section 5, will be soon achieved by slightly increasing the distance between cavities, for example, as (∆x, ∆y) = (8,8). Of course, weak coupling is also realized by decreasing the cavity photon lifetime. However, we found that, when the photon lifetime is further decreased, the numerical integration of coupled-mode equations (1)-(3) become unstable due to the too large time scale difference between the field, carrier, and thermal components [20]. Finally, we comment on an alternative strategy proposed in [51], which is worth considering for the design of coupled PhC cavities. In fact, Ref. [51] demonstrated a barrier engineering technique for robustly tailoring a coupling strength between cavities, which employs the modulation of the radius of air-holes between two cavities.