Emulating the local Kuramoto model with an injection-locked photonic crystal laser array

The Kuramoto model is a mathematical model for describing the collective synchronization phenomena of coupled oscillators. We theoretically demonstrate that an array of coupled photonic crystal lasers emulates the Kuramoto model with non-delayed nearest-neighbor coupling (the local Kuramoto model). Our novel strategy employs indirect coupling between lasers via additional cold cavities. By installing cold cavities between laser cavities, we avoid the strong coupling of lasers and realize ideal mutual injection-locking with effective non-delayed dissipative coupling. First, after discussing the limit cycle interpretation of laser oscillation, we demonstrate the synchronization of two indirectly coupled lasers by numerically simulating coupled-mode equations. Second, by performing a phase reduction analysis, we show that laser dynamics in the proposed device can be mapped to the local Kuramoto model. Finally, we briefly demonstrate that a chain of indirectly coupled photonic crystal lasers actually emulates the one-dimensional local Kuramoto chain. We also argue that our proposed structure, which consists of periodically aligned cold cavities and laser cavities, will best be realized by using state-of-the-art buried multiple quantum well photonic crystals.

www.nature.com/scientificreports/ coupled-mode equations, we numerically demonstrate the synchronization (mutual injection-locking) of two lasers. Furthermore, we confirm that strong-coupling between the two lasers is actually prohibited by the presence of the additional cold cavity. Second, in the same way as in our previous paper 43 , we perform a phase reduction analysis to calculate the phase equations of motion for two indirectly coupled lasers 2,44,45 . The obtained phase equations of motion indicate that the phase dynamics of lasers indirectly coupled via cold cavities is equivalent to the local Kuramoto model. Finally, we demonstrate that a one-dimensional chain of indirectly coupled PhC lasers can emulate the one-dimensional local Kuramoto chain 46 . We also argue that our proposed device can be realized best by using buried multiple quantum well (MQW) PhC cavities [47][48][49][50] , where MQWs are locally embedded in a PhC slab. With this state-of-the-art technology, laser and cold cavities can be periodically aligned on a PhC chip.

Lasers as limit cycle oscillators
Here, we review limit cycle interpretation for laser oscillation. In general, using complex field α and carrier number N, single-mode laser dynamics are, in the nonrotating frame, described by the following rate equations: [51][52][53] where P is the pumping rate for carriers, and ω c is the resonance frequency of the laser cavity. Decay rates γ c and γ are photon and carrier decay rates, respectively. Note that, in this paper, by employing the quantum optics convention, the electric field rotates as α(t) = α(0)e −iω c t , which is opposite to the rotation in conventional coupled-mode equations, [ α(t) = α(0)e iω c t ]. The coefficient β represents the fraction of photons spontaneously emitted into a lasing mode, and it is called the spontaneous emission coupling coefficient 51 . For simplicity, we neglect the linewidth enhancement factor in the rate equations (1) and (2) in the main text. In Section 5 in the supplemental material, we discuss the effect of the linewidth enhancement factor on synchronization, which may be negligible in quantum-dot lasers but generally has non-negligible effects in semiconductor lasers. Here, it is worth noting that, in Eqs. (1) and (2), the terms 1 2 βγ Nα and −βγ � N|α| 2 represent the stimulated emission, while there are no spontaneous emission terms. The effect of spontaneous emission will be included in the rate equations through a field noise term, if necessary. It is also important to note that Eqs. (1) and (2) hold only for a low β(≪ 1) , which is usually the case in most lasers. The rate equations (1) and (2) are known to exhibit Hopf bifurcation, which is equivalent to lasing, when the pump rate reaches a lasing threshold P = P th = γ c /β.
In this paper, for further simplification, we consider the case where the photon lifetime is much longer than the carrier lifetime ( γ c ≪ γ � ), which is called the class-A condition 54 . With this assumption, we adiabatically eliminate the carrier degree of freedom as Ṅ = 0 55,56 . The adiabatic elimination of the carrier dynamics reduces the rate equations (1) and (2) to Equation (3) is the well-known Stuart-Landau equation 2,57 , which is also called the Van der Pol equation 21,58 . Importantly, parameter ε in Eq. (3) is the pump parameter defined as which indicates that the Hopf bifurcation (lasing) again occurs when ε exceeds zero. Actually, when ε > 0 , the field amplitude |α| [see Fig. 1a] increases with an increase in the pump parameter as Therefore, in Eq. (3), the linear γ c εα/2 and nonlinear term βγ c |α| 2 α/2 can be interpreted as gain and gain saturation, respectively. Here, it is important to stress that the laser oscillation itself is interpreted as limit cycle oscillation, and thus the resonance frequency of the laser cavity ω c is the oscillation frequency of the limit cycle. As limit cycle oscillation emerges only in a nonlinear dissipative system with energy injection, lasing is achieved with the cavity decay, pumping, and gain saturation (nonlinearity).
Finally, we briefly comment on the effect of photon-carrier dynamics on synchronization properties, which will be important in real PhC cavity lasers. Since PhC cavity lasers are semiconductor lasers, their carrier lifetime is much longer than the photon lifetime (sometimes called class-B lasers 54 ), and the relaxation oscillation appears around lasing threshold 59,60 . Therefore, in real PhC cavity lasers, the adiabatic elimination approximation of the carrier degree of freedom cannot be justified, and we need to directly simulate the rate equations (1) and (2). Fortunately, we found that Eqs. (1) and (2) quantitatively provide the same results as the Stuart-Landau equation as long as phase dynamics are concerned, which can also be confirmed with the phase equation of motion for the class-B rate equations. See Section 4 in the supplemental material.

Synchronization of two lasers
Coupled-mode equations. Now, we consider the device shown in Fig. 1b, where the two lasers (L1 and L2) are indirectly coupled via the cold cavity (C1). The corresponding coupled-mode equations of motion representing field dynamics are given by where α 1,2 and E 1 represent fields in the laser cavity and cold cavity, respectively. Additionally, ω 1,2 and 1 respectively represent the resonance frequencies of the laser cavities (L1 and L2) and coldcavity (C1). Similarly, γ 1,2 and Ŵ 1 are the field decay rates of the laser-(L1,2) and cold cavity (C1), respectively. The parameter β 1,2 is the spontaneous emission coupling coefficient, while ε 1,2 is the pump parameter for laser L1 and L2. Finally, the two coupling strengths between the cavities are denoted by g 1 and g 2 . For simplicity, in the rest of this paper, we use β 1 = β 2 = 0.001 and ε 1 = ε 2 = 1.0 , which is above the lasing threshold. Furthermore, we use the same values for the normalized decay rates of the laser cavity and cold cavity: Ŵ 1 = γ 2 = γ 1 ≡ 1 , where γ 1 is interpreted as a dimensionless parameter for numerical simulations.

Time evolutions.
By showing field time evolutions described by the coupled-mode Eqs. (6)-(8), we demonstrate the synchronization of two lasers (mutual injection locking). Since the typical laser frequency, which is on the order of hundreds of terahertz, we perform the rotating-frame transformation for all fields, for example, as α 1 e −iω s t → α 1 . With this rotating frame transformation, we shift the resonance frequencies of the cavities as Importantly, there is an arbitrariness in the absolute frequencies, and only the relative frequencies are important. Thus, the frequency of the rotating frame, ω s , is arbitrary, and only the relative values between ω 1 , ω 2 , and 1 matter.
First, Fig. 2a shows the time evolutions of the real parts of the fields Re[α 1 (t)] (black) and Re[α 2 (t)] (blue) for the lasers L1 and L2, respectively, without coupling between cavities g 1 = g 2 = 0 . Without coupling between the cavities, there is no photon in cold-cavity C3, and thus E 1 (t) = 0 . As we expect, in Fig. 2a, the fields in laser L1 and L2 oscillate with their own frequencies: ω ′ 1 = 1γ 1 and ω ′ 2 = 1.01γ 1 . Second, we introduce coupling between the cavities as g 1 = g 2 = 0.1γ 1 in Fig. 2b, where the green curve is the time evolution of the real part of the field in cold-cavity C1 Re[E 1 ] . Figure 2b indicates that the two indirectly coupled laser oscillations exhibit synchronization (mutual injection locking), which is the main result of this paper. Furthermore, the synchronization phase is anti-phase, which is called anti-phase synchronization. Importantly, thanks to cold-cavity C3, normal-mode splitting associated with strong coupling between the two lasers is prohibited, which is confirmed from the fact that the frequency of the synchronized oscillations does not depend on the initial states of two lasers (not shown). In fact, when the system is in the strong-coupling regime, depending on the initial states of two lasers, for example, they form a "bonding" or "anti-bonding" mode, and their frequencies become lower or higher than the original oscillation frequencies 43 . Importantly, no matter how weak the coupling is, directly coupled Figure 1. (a) Laser oscillation is interpreted as limit cycle oscillation in the nonrotating frame, where the laser frequency ω corresponds to the frequency of the limit cycle. With the amplitude |α| = √ ε/β and phase φ = ωt of the laser, we define the limit cycle orbit as (x(φ), y((φ))) = √ ε/β(− cos φ, sin φ) , where x and y are the real Re[α] and imaginary parts Im[α] of the field, respectively. (b) Illustration of two PhC lasers (L1 and L2) indirectly coupled via a cold cavity (C1). The three cavities are evanescently coupled with coupling strengths g 1 and g 2 .  (8), anti-phase synchronization always occurs for any initial state, while if the signs of the two couplings are opposite such as g 2 = −g 1 , in-phase synchronization always occurs (not shown) The sign of a coupling constant depends on the overlap integral of cavity fields and may vary depending on the distance between cavities. In the device design in this paper, since all the distances between cavities are designed to be equal, all the signs of coupling constants can be assumed to be the same. In any case, the property that a synchronization phase does not depend on initial phases of lasers is of importance because the initial phase of PhC lasers cannot be controlled experimentally.

Scientific Reports
Synchronization tree. In Fig. 2c, we show the mean frequencies of the two laser oscillations ω ′ 1 and ω ′ 2 as a function of the coupling between cavities g 1,2 . Since, in general, limit cycle oscillations are quasi-periodic when coupling strength is lower than the critical strength of synchronization, we need to use their mean frequencies obtained with peak detection. Figure 2c clearly indicates that the mean frequencies symmetrically approach each other with an increase in the coupling strength g and that they merge as ω ′ 1 =ω ′ 2 = 1.005γ 1 at the critical strength g = 0.05γ 1 . In fact, the frequency ω ′ 1,2 = 1.005γ 1 is the mean frequency of ω ′ 1 = 1γ 1 and ω ′ 2 = 1.01γ 1 without coupling. Note that the synchronization tree shown in Fig. 2c is approximately symmetric for ω ′ 1 and ω ′ 2 , which is because the parameters are almost the same for L1 and L2.
Furthermore, in Fig. 2d, we plot the mean frequency ω ′ 1,2 as a function of the resonance frequency of the cold cavity 1 , where the coupling strengths are fixed as g 1 = g 2 = 0.07γ 1 while 1 is swept from ω 1 − 1γ 1 to ω 1 + 1γ 1 . Figure 2d indicates that the effective coupling strengths between the cavities can be tuned by changing the resonance frequency of the cold cavity 1 . Intuitively, as the cold-cavity's frequency deviates from the resonance frequencies of the two lasers, the effective coupling strengths decrease. In PhC cavities, the tuning of cavity coupling strength, which is determined by the distance between cavities, is almost impossible. Meanwhile the tuning of the cold cavity's resonance frequency is technically available with the carrier-injection 61,62 or thermooptic techniques 63,64 , and thus the synchronization tree shown in Fig. 2d could be measured.

Phase equations of motion
In this section, as we did in Ref. 43 , by performing the phase reduction analysis 2,65 for Eqs. (6)-(8), we attempt to obtain phase equations of motion. In our case, the phase of limit cycle oscillation is nothing else but the phase of a laser φ as illustrated in Fig. 1a, and thus the interpretation of corresponding phase equations of motion is also straightforward. Furthermore, we show that the determination of phase equations of motion is of importance in terms of mapping our model to the local Kuramoto model. The price to pay for obtaining phase equations of motion is the adiabatic elimination of the field in cold-cavity C1, which is required to transform the indirectly coupled system to a directly coupled model with dissipative coupling.
Adiabatic elimination approximation. The adiabatic elimination of the cold-cavity field degree of freedom Ė 1 = 0 requires that field E 1 rapidly decays compared with the laser field α 1,2 , and thus E 1 adiabatically follows α 1 and α 2 . The time-scale of a variable is generally characterized by its decay rate. Therefore, the conventional adiabatic elimination of field E 1 requires that the decay rate Ŵ 1 must be larger than the decay rates of α 1 and  www.nature.com/scientificreports/ α 2 , as shown in "Lasers as limit cycle oscillators" section, which is not the case, for example, when we consider Ŵ 1 = γ 1,2 as in Fig. 2. However, importantly, the time scale of the laser field α 1,2 is not characterized solely by γ 1,2 . Now, it is important to define the effective decay rates for α 1 , α 2 , and E 1 , including both oscillation frequencies and pump parameters, can be negative due to gain when the pump power is above the threshold ε 1,2 ≥ 0 . According to Ref. 56,66 , when Re[� 1 ] > 0 and Re[ 1,2 ] ≤ 0 , the field E 1 is a "stable" mode that rapidly decays, while the laser fields α 1 and α 2 are unstable modes that do not decay but govern the slow dynamics of the system, which allows putting Ė 1 = 0 (adiabatic elimination). In fact, the unstable mode α 1,2 "enslaves" the stable mode E 1 and plays a role as an "order parameter" (the slaving principle 56,66 ). Now, setting Ė 1 = 0 for Eq. (7), we eliminate the cold-cavity field degree of freedom as By substituting Eq. (9) into Eqs. (6) and (8), we obtain approximated equations of motion: To confirm the validity of this adiabatic elimination approximation, in Fig. 3a, we show synchronization dynamics calculated both with the original equations of motion (6)-(8) and approximated equations (10) and (11). In Fig. 3a, coupling with g 1,2 = 0.1γ 1 is switched on at t = 0 for uncoupled steady-state laser oscillations, and thus the time evolutions of fields represent synchronization dynamics from the unsynchronized to synchronized state. The upper panel in Fig. 3a shows only the synchronization dynamics calculated with the original equations of motion (6)- (8). Meanwhile, in the lower panel, synchronizations calculated with the original equations of motion (solid lines) overlap those calculated with the approximated equations of motion (dashed lines), which clearly indicates that two time evolutions are almost indistinguishable and that the adiabatic elimination approximation is surprisingly good. Note that, to clearly show the synchronization dynamics in Fig. 3a, we used shifted frequencies ω ′ 1 = 0.2γ 1 , ω ′ 2 = ω ′ 1 + �ω = 0.21γ 1 and � ′ 1 = ω ′ 1 = 0.2γ 1 , which are lower than those Fig. 1. As we commented in "Coupled-mode equations" section, these shifts of the resonance frequencies do not change the physics, because only the relative relationship between the resonance frequencies is important. Since the field in the cold cavity was adiabatically eliminated, Eqs. (10) and (11) represent directly coupled lasers. Furthermore, in Eqs. (10) and (11), the effective couplings represented by −(2g 1 g 2 / Ŵ 1 )α 2 and −(2g 1 g 2 / Ŵ 1 )α 1 are non-energyconserving dissipative couplings, which intuitively explains why normal-mode splitting does not appear in our model. Additionally, in Eqs. (10) and (11), the effective dissipative coupling does not have the time delay. www.nature.com/scientificreports/ Finally, we comment on synchronization with a large coupling strength. We found that Eqs. (10) and (11) fail to reproduce synchronization dynamics when g 1,2 ≥ γ 1,2 , Ŵ 1 , which is because the adiabatic elimination approximation cannot describe coherent intensity oscillation between cavities associated with this parameter region [see Section 2 in the Supplemental Material (SM)]. Therefore, the complete conditions required for the adiabatic elimination approximation are Here, it is also important to stress that, although the adiabatic elimination fails to describe synchronization dynamics, even when g 1,2 ≥ γ 1,2 , Ŵ 1 , stable synchronization itself can occur and the adiabatic elimination approximation well reproduces the steady-state synchronized oscillations (see Section 2 in the SM). Furthermore, even when the coupling is extremely strong, for example, g 1,2 = 10γ 1 , we can observe stable synchronization, where no normal-mode splitting is present (not shown). This insensitivity to coupling strength will be advantageous in terms of real device designs, because adjusting the value of weak coupling strength is technically difficult 43 . Furthermore, if coupling is sufficiently strong, we may prove synchronization from spectral shapes, which is discussed again in "Discussion" section.
Phase reduction analysis. Now, we perform the phase reduction analysis for equations of motion (10) and (11), which were obtained with the adiabatic elimination approximation. Here, we make use of the consequence of the phase reduction theory without going into the theoretical detail, which is briefly provided in Section 1 in the SM (further details can be found in our recent paper 43 and in Refs 2,44 ). The objective of the phase reduction analysis is to obtain the phase equations of motion for the phases of the laser L1 ( φ 1 ) and L2 ( φ 2 ) represented as where Ŵ 12 (φ) and Ŵ 21 (φ) are called the phase-coupling functions. For the approximated equations of motion (10) and (11), we found that Ŵ 12 (φ) and Ŵ 21 (φ) can be analytically calculated as Finally, the phase difference between the two lasers ψ ≡ φ 2 − φ 1 follows the following simple equation of motion: where �ω ≡ ω 2 − ω 1 is the frequency difference between the two lasers already defined in "Coupled-mode equations" section. Here, Ŵ a (ψ) ≡ Ŵ 21 (ψ) − Ŵ 12 (−ψ) is the anti-symmetric part of the phase coupling function Ŵ 21 (ψ) , which is shown in Fig. 3b. For a negligible laser frequency difference �ω ≃ 0 , since Ŵ a (π) = 0 and Ŵ ′ a (π) < 0 hold for Eq. (16), phase locking occurs at the phase ψ = φ 2 − φ 1 = π , which is anti-phase synchronization as expected from the simulations [see the arrows in Fig. 3b]. Meanwhile, since Ŵ a (0) = 0 and Ŵ ′ a (0) > 0 hold for φ = 0 , the phase φ = 0 is an unstable fixed point. Of course, for a non-negligible frequency difference �ω = 0 , the synchronization phase shifts from π . The phase equations of motion predict not only the synchronization phase but also the critical coupling strength of synchronization. For Eq. (16) to have a phase-locking solution, the condition −4g 1 g 2 / Ŵ 1 ≤ �ω ≤ 4g 1 g 2 / Ŵ 1 must be satisfied. For the oscillation frequency difference �ω = 0.01γ 1 and cold-cavity decay rate Ŵ 1 = 1γ 1 , which are assumed in Fig. 2c, synchronization occurs when the coupling strengths reach g 1 = g 2 = 0.05γ 1 [see Fig. 2c] because the above phase-locking condition is satisfied with these parameters as 4g 1 g 2 / Ŵ 1 = 0.01γ 1 = �ω.
Furthermore, the analytically calculated phase coupling functions in Eq. (15) are also of importance for mapping our model to the local Kuramoto model. In fact, the phase equations of motion are explicitly written as where g ij ≡ 2g i g j / Ŵ 1 ( g ij =g ji ) is the effective coupling strength. The phase equations of motion (17) and (18) are straightforwardly extended to a one-dimensional chain or two-dimensional array as where N j represents the nearest neighbour sites of the ith site. Importantly, the coupled phase oscillator described by Eq. (19) is equivalent to the local Kuramoto model 3,25 . Note that, in the original Kuramoto model, the sign of the coupling is minus as −g 21 , and thus in-phase synchronization occurs.
In conclusion, with the aide of the phase reduction theory, we proved that an array of lasers with cold-cavitymediated coupling can emulate the nearest-neighbor coupled Kuramoto model (the local Kuramoto model). Note that, of course, the strict mapping of given coupled-mode equations to the local Kuramoto model (19) requires an adiabatic elimination condition similar to Eq. (12).

Array configuration
Although the investigation of rich physics emerging from coupled phase oscillators is beyond the scope of this paper, we briefly simulate a one-dimensional chain of indirectly coupled PhC lasers and demonstrate that our device can actually reproduce collective dynamics predicted for the one-dimensional local Kuramoto chain 46 . The chain of indirectly coupled PhC lasers is schematically illustrated in Fig. 4a, where eleven laser cavities and ten cold cavities are alternately aligned. Of course, the configuration of cavities to realize the local Kuramoto chain is not limited to that shown in Fig. 4a, and various configurations can be imagined. For a one-dimensional chain, in principle, even the periodic boundary condition may be implemented with a ring-like configuration. For simplicity, for all the laser cavities, we assume β i = 0.001 , ε i = 1.0 , and γ i ≡ 1 . Similarly, all the cold cavities have the same resonance frequencies and photon decay rate: � i = 1γ 1 and Ŵ i = γ 1 ≡ 1 for all i. In Section 6 in the supplemental material, we demonstrate large-scale synchronization when the parameter values of all laser and cold cavities are slightly different. Furthermore, as in "Synchronization of two lasers" section, we assume that all the coupling constants have the same strengths: g i = g for all i. Meanwhile, the resonance frequencies of the eleven laser cavities are randomly distributed around a mean frequency ω i = 1γ 1 [for the actual values of ω i , please see the caption of Fig. 4]. Note that since all the laser and cold cavities have similar resonance frequencies and the coupling strengths are smaller than the cavity decay rates, an adiabatic elimination condition similar to Eq. (12) is satisfied, and thus corresponding simple phase equations of motion are expected to exist. By directly simulating the full coupled-mode equations corresponding to the configuration shown in Fig. 4a, we calculated the mean frequencies of the laser oscillations as a function of the coupling strength g [see the synchronization tree in Fig. 4b]. As Fig. 4b indicates, with an increase in coupling strength g, synchronized clusters are gradually formed, and finally all clusters merge into a single fully synchronized cluster at g ≃ 0.07γ 1 [see G on Fig. 4b]. Similarly to Ref. 46 , when two [at A, B, C in Fig. 4b] or three [at D, E in Fig. 4b] adjacent oscillators (or clusters) have close oscillation frequencies, they form a new synchronized cluster with an increase in coupling strength. When adjacent clusters have largely different frequencies, while non-adjacent clusters have similar frequencies, the non-adjacent clusters form a synchronized cluster. In fact, the synchronization denoted by F in Fig. 4b consists of the non-adjacent oscillators (clusters) L1-3 and L7-11. Furthermore, in Fig. 4c, we show a    Fig. 4b. The fact that both synchronization trees have almost the same structures indicates that our proposed device will actually emulate the local Kuramoto model. Finally, the time evolutions of the laser oscillations without ( g = 0 ) and with coupling ( g = 0.1γ 1 ) are shown in the left and right panels of Fig. 4d, respectively. When there is no coupling, as we expect, the laser oscillations are totally uncorrelated, while all the laser oscillations are fully synchronized with coupling g = 0.1γ 1 . Interestingly, in this fully synchronized state [see the right panel in Fig. 4d], the phases are opposite between the even and odd sites of the lasers oscillations. Therefore, even in the one-dimensional chain, a pair of adjacent laser oscillations exhibit anti-phase synchronization. Note that, in Fig. 4c,d, the "de-synchronization" discovered in 46 was not observed, which may be due to the small number of oscillators or, more interestingly, could be associated with anti-phase synchronization. We also comment on the offsets of the synchronization phases in the fully-synchronized oscillations shown in the right panel of Fig. 4d, where the synchronization phases slightly differ depending on the pair of the synchronized oscillations. We found that, with a further increase in coupling strength, these offsets of the synchronization phases disappear and that all the pairs of synchronized oscillations become indistinguishable.

Discussion
Here, we discuss several details that will be of importance in a real device design and experiments. In experiments, the easiest method to observe synchronization may be the spectral measurement of laser emissions. Since limit cycle oscillation frequencies are equivalent to laser oscillation frequencies, synchronization can be directly confirmed by the number of emission peaks in a measured spectrum. Namely, if a spectrum has a single emission peak, two lasers are synchronized, while if there are two emission peaks, they are not. Although coupling strengths between cavities are usually fixed in a device, it is still possible to actively tune the resonance frequency of a cold cavity [61][62][63][64]67 and effectively change coupling strengths as shown in Fig. 2d. In this context, the spectral shape of laser emissions will be of interest. Below the lasing threshold, since laser cavities will behave as "cold cavities", their emission spectrum is expected to exhibit normal-mode splitting. Meanwhile, above the lasing threshold the emission spectrum exhibits a single peak due to synchronization. Therefore, we may prove synchronization from the pump-power dependence of the change in spectral shape . Another promising experimental strategy to prove synchronization may be to pump two lasers independently and tune the respective laser frequencies by making use of the carrier-induced blue shift 49 . This strategy can be easily realized with spatially separated optical pumps or two electrodes for electric pumping.
The buried MQW PhC laser technique in the design of a real device is reported in Refs. 47,48 , where the PhC slab and buried PhC are composed of InP and InGaAsP/InGaAs, respectively. Furthermore, buried MQW PhC lasers can be pumped optically or electrically. If the photon lifetime of buried MQW PhC lasers is assumed to be 1/γ 1 = 1 ps ( ∼160 GHz), the frequency difference between two lasers corresponding to �ω = 0.01γ 1 , which is assumed in the simulations in Fig. 2, is �ω = 0.01γ 1 ∼ 1.6 GHz. This laser frequency difference may seem to be severe for experimental realization (even with the state-of-the-art fabrication technology, the frequency difference between cavities may be about 50 GHz 68 ), but we found that, qualitatively, the same synchronization can occur for a larger frequency difference. For example, synchronization with a laser frequency difference �ω = 0.1γ 1 is discussed in Section 3 in the SM. Furthermore, we found that even if all the parameters of the three cavities including β 1,2 are moderately different, synchronization can occur (not shown).

Conclusion and outlook
To conclude, we theoretically proposed a design of indirectly-coupled PhC cavity lasers that emulates the local Kuramoto model. In this study, we reinterpreted the injection-locking phenomenon of lasers as the synchronization of limit cycle oscillations. Furthermore, our design prevents laser oscillations from forming normal-modes (strong-coupling) with indirect coupling via additional cold cavities and realizes effective dissipative coupling without time-delay. Experimentally, this proposed structure will best be realized best by using buried MQW PhC cavities. First, after modelling laser oscillation with the Stuart-landau equation, we numerically demonstrated the synchronization of two indirectly-coupled PhC lasers using the coupled-mode equations of motion. Second, by applying the phase reduction theory to the two indirectly coupled lasers, we obtained corresponding phase equations of motion, which are equivalent to the local Kuramoto model. Finally, we briefly discussed synchronization dynamics for a one-dimensional chain of indirectly coupled PhC lasers and demonstrated that the proposed device can actually emulate the local Kuramoto chain.
For future perspectives, first of all, the one-dimensional local Kuramoto model briefly investigated in "Array configuration" section, already comprises rich physics that were actively investigated by detailed numerical simulations 46 and renormalization group analysis 25,69 . Furthermore, very recently, Ref. 70 demonstrated that even topological phenomena emerge in the one-dimensional chain of limit cycles. Thanks to the scalability of PhC cavities, the extension of the one-dimensional chain of PhC lasers to a two-dimensional array is straightforward, which is the realization of the celebrated two-dimensional local Kuramoto model 3,[24][25][26][27][28]71 . Compared with the in-phase synchronization case, large-scale anti-phase synchronization has not yet been drawing attention. For instance, as Ref. 72 indicates that a large anti-phase synchronization network is not possible, large-scale anti-phase synchronization itself may be of fundamental interest. Another important direction will be the inclusions of classical and quantum noise effects in indirectly coupled PhC lasers, which will provide spectral information. As we briefly discussed in "Discussion" section, we may prove synchronization in terms of the pump power dependence of spectral shape. In this direction, it is also easy to construct a quantum model corresponding to our coupled-mode equations. In fact, the quantum counterpart of the classical Stuart-Landau model is the Scully-lamb master equation 73,74 . Therefore, even the effect of quantum noises on synchronization 21 www.nature.com/scientificreports/ be tested with the proposed device. Moreover, since the synchronization problem in the local Kuramoto model is analogous to the energy minimization problem in the XY model, our device may be used for simulating the spin system in statistical physics 72,76 . Finally, from the standpoint of practical application, an injection-locked (synchronized) PhC laser array can be employed as a single-mode high-power PhC laser. Even though every PhC laser unavoidably has a different oscillation frequency, in the fully synchronized state, they behave as a laser with a single frequency. Furthermore, this type of a laser will also have high coherence because all laser phases are locked in the synchronized state.

Methods
All the time evolutions were obtained by integrating the coupled-mode equations of motion with the conventional Runge-Kutta method. The synchronization trees were calculated as the mean oscillation frequencies of time evolutions. The calculation of the mean frequencies is based on the peak detection technique. To precisely determine the mean frequencies, long time evolutions (typically 6000γ −1 1 ) were required.