Two-dimensional phase-space picture of the photonic crystal Fano laser

TL;DR

This paper simplifies the complex dynamics of a photonic crystal Fano laser to a 2D phase-space model, revealing new insights into its stability and operation.

Contribution

It reduces a 5D dynamic model to a 2D model, providing a clearer understanding of the laser's stability and the origin of instabilities.

Findings

Laser operation is confined to two degrees-of-freedom after initial transition.

The 1D model reduces to a 2D model after the exceptional point.

A new unstable periodic orbit explains the origin of instability.

Abstract

The recently realized photonic crystal Fano laser constitutes the first demonstration of passive pulse generation in nanolasers [Nat. Photonics $\boldsymbol{11}$, 81-84 (2017)]. We show that the laser operation is confined to only two degrees-of-freedom after the initial transition stage. We show that the original 5D dynamic model can be reduced to a 1D model in a narrow region of the parameter space and it evolves into a 2D model after the exceptional point, where the eigenvalues transition from being purely to a complex conjugate pair. The 2D reduced model allows us to establish an effective band structure for the eigenvalue problem of the stability matrix to explain the laser dynamics. The reduced model is used to associate a previously unknown origin of instability with a new unstable periodic orbit separating the stable steady-state from the stable periodic orbit.

Figures

Tables

Table 1. Table 1: Laser parameters used in all numerical simulations.

Parameter name	Symbol	Value
Transparency carrier density	$N_0$	$1\cdot {10}^{24}{m}^{-3}$
Parity of the cavity mode	$P$	$1$
Linewidth enhancement factor	$\alpha$	$1$
Internal loss factor	${\alpha }_i$	$1000m^{-1}$
Lasing cavity length	$L$	$5\mu m$
Carrier lifetimes	${\tau }_s$, ${\tau }_c$	$0.5$ns
Laser cavity volume	$V_{LC}$	$1.05\mu m^3$
Nanocavity volume	$V_{NC}$	$0.243\mu m^3$
Nanocavity resonance	${\lambda }_r$	$1554nm$
Reference refractive index	$n_{ref}$	$3.5$
Group refractive index	$n_{grp}$	$3.5$
Differential gain	$g_N$	$5\cdot {10}^{-20}{m}^2$
Waveguide confinement factor	$\Gamma$	$0.5$
Nanocavity confinement factor	${\Gamma }_C$	$0.3$
Left mirror reflectivity	$R_1$	$1$
Nanocavity-waveguide coupling	${\gamma }_C$	$1.14{\mathrm{ps}}^{-1}$
Nanocavity total passive decay rate	${\gamma }_T$	$1.21{\mathrm{ps}}^{-1}$

Equations

$$r_1({\omega }_s)r_R({\omega }_s,{\omega }_c)e^{i2k({\omega }_s,N_s)L}=1,$$

$$\frac{d|a_{+}|(\tau)}{d\tau}=-\frac{\gamma_{L}|a_{+}|(\tau)}{G_{NC}N_{0}}+\frac{\Gamma|a_{+}|(\tau)(n_{w}(\tau)-n_{s})}{2\Gamma_{C}}+\frac{\gamma_{L}}{G_{NC}N_{0}}|a_{+}|(\tau)\mathrm{Re}\bigg{(}\frac{A_{-}(\tau)}{r_{R}A_{+}(\tau)}\bigg{)}$$

$$\frac{d|a_{-}|(\tau)}{d\tau}=-\frac{P\gamma_{C}}{G_{NC}N_{0}}|a_{-}|(\tau)\mathrm{Re}\bigg{(}\frac{A_{+}(\tau)}{A_{-}(\tau)}\bigg{)}-\frac{\gamma_{T}|a_{-}|(\tau)}{G_{NC}N_{0}}+\frac{|a_{-}|(\tau)(n_{c}(\tau)-1)}{2}$$

$$\frac{d\Delta\phi(\tau)}{d\tau}=-\frac{\alpha}{2}(n_{c}(\tau)-1)+\frac{\Gamma\alpha(n_{w}(\tau)-n_{s})}{2\Gamma_{C}}-\frac{\Delta\omega}{G_{NC}N_{0}}-\mathrm{Im}\bigg{(}\frac{r_{R}P\gamma_{C}{A}^{2}_{+}(\tau)+\gamma_{L}{A}^{2}_{-}(\tau)}{r_{R}G_{NC}N_{0}A_{-}(\tau)A_{+}(\tau)}\bigg{)}$$

$$\frac{d{n}_w(\tau )}{d\tau }=\frac{-|a_{+}(\tau ){|}^2(n_w(\tau )-1)-n_w(\tau )+j_c}{G_{NC}{N}_0{\tau }_s}$$

$$\frac{d{n}_c(\tau )}{d\tau }=\frac{-|a_{-}(\tau ){|}^2(n_c(\tau )-1)-n_c(\tau )}{G_{NC}{N}_0{\tau }_c}$$

$$\vec{\psi }(\tau )=\{|a_{+}(\tau )|,|a_{-}(\tau )|,\Delta \varphi (\tau ),n_w(\tau ),n_c(\tau )\}$$

$$\frac{d\vec{\psi }(\tau )}{d\tau }=V(\vec{\psi }(\tau ))$$

$$\frac{d^2\vec{\psi }(\tau )}{d{\tau }^2}={\nabla }_{\vec{\psi }}V(\vec{\psi }(\tau ))\frac{d\vec{\psi }(\tau )}{d\tau }=\mathbf{A}(\vec{\psi }(\tau ))\frac{d\vec{\psi }(\tau )}{d\tau }$$

$$\vec{\psi }=$$

$$\frac{d\vec{\psi }}{d\tau }=\{$$

$$V_4(|a_{-}|,n_w),V_5(|a_{-}|,n_w)\}.$$

$$\bigg{(}\frac{d^{2}\vec{\psi}(\tau)}{d\tau^{2}}\bigg{)}_{i}=\bigg{\{}\frac{\partial V_{i}}{\partial|a_{-}|}\frac{d|a_{-}|}{d\tau}+\frac{\partial V_{i}}{\partial n_{w}}\frac{dn_{w}}{d\tau}\bigg{\}},$$

$$\frac{d|a_{-}(\tau )|}{d\tau }=V_2(|a_{-}(\tau )|,n_w(\tau )),$$

$$\frac{d{n}_w(\tau )}{d\tau }=V_4(|a_{-}(\tau )|,n_w(\tau )),$$

$$P_{+}(t)=2{ϵ}_0{n}_{ref}{c}_0|A_{+}(t)+P{A}_{-}(t){|}^2,$$

$$\frac{d\vec{\psi }}{d\tau }=c_1(\tau )\overrightarrow{v_1}(\tau )+c_2(\tau )\overrightarrow{v_2}(\tau )$$

$$\vec{r}^{\prime}(\tau)=\bigg{\{}\frac{d|a_{-}|}{d\tau},\frac{dn_{w}}{d\tau},\frac{d\vec{\psi}}{d\tau}\bigg{\}},$$

$$\frac{d\vec{\psi }}{d\tau }=\frac{\partial \vec{\psi }}{\partial |a_{-}|}\frac{d|a_{-}|}{d\tau }+\frac{\partial \vec{\psi }}{\partial n_w}\frac{d{n}_w}{d\tau }.$$

$$\vec{r}^{\prime}_{1}=\bigg{\{}1,0,\frac{\partial\vec{\psi}}{\partial|a_{-}|}\bigg{\}},\hskip 18.49411pt\vec{r}^{\prime}_{2}=\bigg{\{}0,1,\frac{\partial\vec{\psi}}{\partial n_{w}}\bigg{\}},$$

$$\left[\begin{array}{c}{\vec{v}}_1^{\prime }\\ {\vec{v}}_2^{\prime }\end{array}\right]={\left[\begin{array}{cc}{v}_{12}& v_{14}\\ v_{22}& v_{24}\end{array}\right]}^{-1}\left[\begin{array}{ccccccc}{v}_{12}& v_{14}& v_{11}& v_{12}& v_{13}& v_{14}& v_{15}\\ v_{22}& v_{24}& v_{21}& v_{22}& v_{23}& v_{24}& v_{25}\end{array}\right]$$

Full text

††preprint: APS/123-QED

I Introduction

Integrated photonic circuits require energy efficient, fast and compact light sources miller_device_2009 . Particularly promising candidates to realize them are photonic crystal (PhC) lasers due to their flexibility in design and precise control of the cavity properties akahane_high-q_2003-1 ; tran_directive_2009 . PhC lasers can be electrically driven and allow for modulation in the GHz range matsuo_ultralow_2013 ; jang_sub-microwatt_2015 . Moreover, they have been shown to exhibit very rich dynamics, e.g. spontaneous symmetry breaking hamel_spontaneous_2015 . Recently, a new type of PhC laser has been proposed mork_photonic_2014 where one of the mirrors arises due to a Fano resonance fano_effects_1961 ; limonov_fano_2017 . Furthermore, this laser has been demonstrated to be able to generate a self-sustained train of pulses at GHz frequencies, a property that has been observed only in macroscopic lasers thus far yu_demonstration_2017 . Generation of pulses by an ultracompact laser is of interest for applications in future on-chip optical signal processing.

The configuration of the Fano laser is shown in Fig. 1. The active material may be composed of several layers of InAs quantum dots or quantum wells and is incorporated inside the InP PhC membrane. The laser cavity is composed of a PhC line-defect waveguide blocked with a PhC mirror on the left side forming a broad-band mirror, whereas the right mirror is due to the Fano interference between the nanocavity and the waveguide. The Fano resonance arises due to the interference of a discrete mode of the nanocavity with the continuum of PhC waveguide modes. The spectral width of the resonance is determined by the quality factor of the nanocavity enabling the realization of a narrow-band mirror. The dynamic operation of the laser is modeled using a combination of coupled-mode theory and conventional laser rate equations rasmussen_theory_2017 . The model has been used to demonstrate that there are two regimes of operation, the continuous-wave regime and a self-pulsing regime rasmussen_theory_2017 . Particularly, it has been shown that as the real part of any of the eigenvalues of the underlying stability matrix, evaluated at the steady-state, becomes positive, the relaxation oscillation becomes undamped resulting in the laser becoming unstable and self-pulsing behaviour setting in rasmussen_theory_2017 ; rasmussen_modes_2018 . However, it does not fully explain the origin of instability in the whole parameter space of the laser as there exists a region in which the laser can become unstable even when all the steady-state eigenvalues are negative rasmussen_theory_2017 .

The purpose of this work is to analyze not only the steady-state eigenvalues of the stability matrix of the dynamic model, but also the instantaneous eigenvalues during the laser operation. Moreover, we determine the ’minimal’ model for the laser that is required to explain the dynamics in different regimes, thereby obtaining an alternative perspective on the dynamics of the Fano laser. We thus demonstrate that the laser operation can be effectively modeled by a one-dimensional (1D) system of differential equations in a limited region of the parameter space when the steady-state eigenvalues are purely real and that it evolves into an effective two-dimensional (2D) system beyond the steady-state exceptional point, when the eigenvalues form complex conjugate pairs. These findings are used to determine the origin of the instability that is observed when the steady-state eigenvalues are negative. We notice that the analysis of instabilities and chaos in injection-locked lasers, e.g., using bifurcation analysis has been very successful mork_chaos_1992 ; krauskopf_bifurcation_2000 ; wieczorek_dynamical_2005 ; erzgraber_bifurcation_2007 . Here, we use it to analyse origin of laser instability when the steady-state is stable and set an important goal to identify reduced systems for getting further physical insight.

The manuscript is organized as follows. In Sec. II we introduce the model used to describe the laser dynamics. In Sec. III we show that the laser operation can be understood by means of a 2D phase-space picture and we analyze the steady-state and instantaneous eigenvalues of the stability matrix. In Sec. IV we exploit the simplified 2D model to associate a self-pulsing operation, when a steady-state is stable to a generalized Hopf (Bautin) bifurcation which is characterized by a bifurcation of two periodic orbits and an equilibrium point (steady-state) govaerts_numerical_2000 ; kuznetsov_elements_2004 .

II Dynamic model of the Fano laser

We next briefly describe the procedure required to establish the dynamic model of the Fano laser; for more details refer to rasmussen_theory_2017 . The complex field is decomposed into the fields propagating to the right and left from the reference plane, see Fig. 1. By combining the boundary conditions for both fields, we can arrive at the oscillation condition mork_photonic_2014 :

[TABLE]

where $r_{1}$ and $r_{R}$ are the broad-band (left) and the narrow-band (right) reflection coefficients, respectively. $r_{R}$ is determined using the coupled-mode theory fan_temporal_2003 ; wonjoo_suh_temporal_2004 ; kristensen_theory_2017 , while $r_{1}$ is the reflection coefficient due to the PhC band gap and has to be transformed towards the common reference plane using standard transmission line theory tromborg_transmission_1987 . $k$ is the complex wavenumber of the waveguide, $L$ is the length of the lasing cavity, and $\omega_{c}$ is the resonance frequency of the nanocavity. The condition in Eq. (1), is solved for $(\omega_{s},N_{s})$, which are the steady-state lasing frequency and carrier density, respectively. They serve as expansion points of the dynamic model. There are multiple solutions of Eq. (1) mork_photonic_2014 ; rasmussen_modes_2018 among which the one with the lowest modal threshold gain is chosen. The wavenumber, $k$, accounts for dispersion of the refractive index of the PhC membrane and the gain of the active material.

Subsequently, the boundary condition is solved for the left propagating field and then the term $1/r_{1}(\omega)e^{i2k(\omega_{s},N_{s})L}$ is Taylor expanded around the steady-state operation point $(\omega_{s},N_{s})$ and a first-order differential equation for the right-propagating complex field envelope evaluated at the reference plane $A_{+}(t)$ is derived using the Fourier transform. In the special case of an open waveguide considered here, the coupled-mode equation for the field in the nanocavity can be directly reformulated as an equation for the left-propagating complex field envelope evaluated at the reference plane $A_{-}(t)$. The equations for $A_{+}(t)$ and $A_{-}(t)$ are complemented with the traditional rate equations for carrier densities in the waveguide and the nanocavity.

Since the variables introduced above differ by orders of magnitude, we introduce dimensionless near-unity variables in order to improve numerical stability. Moreover, detuning from the expansion point frequency $\omega_{s}$ results in time-varying real and imaginary parts of $A_{+}(t)$ and $A_{-}(t)$ at the steady-state. Because of that, the differential equations for $A_{+}(t)$ and $A_{-}(t)$ are separated into equations for amplitudes and phase evolutions by the following substitution: $A_{+}(t)=A_{0+}|a_{+}|(\tau/(G_{NC}N_{0}))e^{i\phi_{+}(\tau/(G_{NC}N_{0})}$, $A_{-}(t)=A_{0-}|a_{-}|(\tau/(G_{NC}N_{0}))e^{i\phi_{-}(\tau/(G_{NC}N_{0})}$, where $A_{0+}$ and $A_{0-}$ are the normalization constants, $a_{-}(\tau)$ and $a_{+}(\tau)$ are the normalized complex field envelopes, $\tau=tG_{NC}N_{0}$ is the normalized time, $G_{NC}=\Gamma_{C}v_{g}g_{N}$, and $v_{g}$ is the group velocity. The system depends solely on the phase difference $\Delta\phi(\tau)=\phi_{-}(\tau)-\phi_{+}(\tau)$; thus by subtracting the equations for phase evolutions $\phi_{+}(\tau)$, $\phi_{-}(\tau)$ and exploiting linearity of differentiation, these equations can be combined into one. This leads us to the following system of five differential equations describing the dynamics of the laser:

[TABLE]

Here, $n_{w}(\tau)$ and $n_{c}(\tau)$ are the carrier densities in the waveguide and the nanocavity, respectively, normalized with respect to the transparency carrier density, $N_{0}$. $\gamma_{L}=v_{g}/(2L)$ is the inverse of the cavity roundtrip time, $n_{s}$ is the steady-state carrier density obtained from the oscillation condition normalized with respect to $N_{0}$, $\Delta\omega=\omega_{c}-\omega_{s}$ is the detuning of the steady-state lasing frequency, $\omega_{s}$, from the cavity resonance frequency, $\omega_{c}$, and $j_{c}$ is the normalized effective pumping current, which includes the injection efficiency.

Subsequently, we linearize the problem by calculating the total derivative of Eqs. (2) with respect to $\tau$. The system of equations describing the laser dynamics in Eqs. (2) can be expressed in the short form as a function $V(\cdot)$ of the state vector $\vec{\psi}(\tau)$:

[TABLE]

By taking the total derivative of $V(\vec{\psi}(\tau))$, we obtain a directional derivative along the curve parameterized by $\tau$:

[TABLE]

Consequently, $d\vec{\psi}(\tau)/d\tau$ in Eq. (3b) is interpreted as the velocity of the state vector and is expressed as a function of the current position of the state vector in Eq. (2). $d^{2}\vec{\psi}(\tau)/d\tau^{2}$ is interpreted as the acceleration of the state vector, see Eq. (4). Matrix $\bm{A}$ is the so-called Jacobian matrix; its eigenvalues $\lambda$ are used to determine the stability of the laser when evaluated at the steady-state. The system is stable if all eigenvalues have negative real parts. On the other hand, if any eigenvalue has a positive real part, the system is unstable. The matrix $\bm{A}$ is purely real, but not symmetric as we separated the complex field envelopes into the magnitudes $|a_{+}(\tau)|$, $|a_{-}(\tau)|$ and the phase difference $\Delta\phi$. Therefore, the matrix is non-Hermitian and we have to distinguish between right $\vec{v}$ and left $\vec{w}$ eigenvectors which are normalized so that $\bm{W}^{T}\bm{V}=\bm{I}$ is satisfied morse_methods_1953 ; ibanez_adiabaticity_2014 . The columns of $\bm{W}$ and $\bm{V}$ are the left and the right eigenvectors, and $\bm{I}$ is the identity matrix. Furthermore, eigenvalues of the matrix $\bm{A}$ can be purely real or form complex conjugate pairs arfken_mathematical_2005 . In the following sections, we use Eqs. (2) and (5) to investigate the origin of instability in case of a stable steady-state and to show that the original laser model can be simplified to a system of two differential equations.

III Two dimensional phase-space picture

III.1 Steady-state eigenvalues

Above the threshold, the laser can exhibit two types of operation, the continuous wave and the self-pulsing operation rasmussen_theory_2017 ; yu_demonstration_2017 . Figure 2 shows the real and imaginary parts of the two steady-state eigenvalues of the Jacobian matrix, $\bm{A}$ with the largest real parts plotted versus $\Delta\omega_{c}$, which is the detuning of the cavity resonance frequency, $\omega_{c}$ from the resonance frequency of the isolated cavity, normalized with respect to $\gamma_{T}$. It is noted that $\Delta\omega_{c}$ defines $\Delta\omega$ in Eq. (2) through the oscillation condition, Eq. (1), and is controlled externally. As our case study, we choose $\Delta\omega_{c}$ marked with the blue line in Fig. 2.

Interestingly, it has been observed in rasmussen_theory_2017 that in the vicinity of the $\Delta\omega_{c}$ marked by the blue dashed line in Fig. 2 the laser can exhibit the continuous wave or the self-pulsing operation depending on the initial condition despite its steady-state eigenvalues having negative real parts and thus suggesting stable operation of the laser. However, origin of this instability has not been explained and is examined in Sec. IV. On the other hand, when $\Delta\omega_{c}$ is increased beyond $1.26\gamma_{T}$, the real parts become positive, the relaxation oscillation becomes undamped, the laser becomes unstable, and the state approaches a stable periodic orbit for any initial condition rasmussen_theory_2017 . All the following figures are obtained for the parameters listed in Table 1, while the pumping current is set to $J=1.2J_{thr}$, where $J_{thr}$ is the minimum threshold current.

III.2 Exceptional points

It is interesting to observe in Fig. 2 that for $\Delta\omega_{c}$ lower than $-1.72\gamma_{T}$, the real part of the two eigenvalues split and the eigenvalues become purely real, see Fig. 2b). At $\Delta\omega_{c}=-1.72\gamma_{T}$, the two eigenvalues coalesce and not only the eigenvalues are identical at this point, but so are the eigenvectors as well dembowski_experimental_2001 ; heiss_exceptional_2004 ; berry_physics_2004 ; heiss_physics_2012 ; liertzer_pump-induced_2012 . This constitutes an exceptional point which is also known as a symmetry breaking point for a non-Hermitian system bender_generalized_2002 ; ruter_observation_2010 ; feng_non-hermitian_2017 ; el-ganainy_non-hermitian_2018 . However, exceptional points are a general phenomenon observed in optical waveguides klaiman_visualization_2008 , unstable laser resonators berry_mode_nodate , coupled PhC nanolasers kim_direct_2016 , quantum systems lefebvre_resonance_2009 , electronic circuits stehmann_observation_2004 and mechanical resonators xu_topological_2016 . They only require non-Hermiticity of the system for their existence heiss_repulsion_2000 ; berry_physics_2004 ; heiss_exceptional_2004 . We emphasize that exceptional points may arise upon coalescence of eigenvectors/eigenvalues of any matrix, e.g. a Hamiltonian matrix rotter_non-hermitian_2009 , a S-parameter matrix chong_symmetry_2012 and an impedance/admittance matrices hanson_exceptional_2004 , to name a few. Exceptional points have also been linked to a self-pulsing mechanism in distributed feedback lasers bandelow_theory_1993 ; wenzel_mechanisms_1996 in which case the self-pulsing mechanism was attributed to dispersive quality factor self-switching similarly as in the case of the Fano laser yu_demonstration_2017 . In the present case, exceptional points arise due to dissimilar decay rates, $\gamma_{C}$, $\gamma_{L}$ and phenomenologically introduced gain terms $|a(\tau)|(n(\tau)-1)$, Eqs. (2a) and (2b). They play an analogous role to the loss and gain usually introduced as an imaginary part of the refractive index in parity-time symmetric systems feng_single-mode_2014 ; hodaei_parity-time-symmetric_2014 .

III.3 Two-dimensional phase-space

In Fig. 3 we plot three trajectories of $n_{c}$, marked in red, green and blue versus $n_{w}$, $|a_{-}|$ and obtained for the three different initial conditions. The trajectories are parameterized by $\tau$. It is found that there are actually three different stages of the laser operation: the initial transition stage, the later transition stage and the self-pulsing stage. This is in contrast to the previously reported picture of two stages: the transition stage and the self-pulsing stage. The red, green and blue dots mark the initial conditions in Fig. 3. It is found that at first they lie above a yellow surface, this is the initial transition stage which lasts only a few picoseconds. After a very short initial transition stage, the state reaches the surface at the time instant marked with the orange dashed line in Fig. 3b). The state stays on the surface within the later transition stage and the self-pulsing/continuous wave stage. The later transition stage is when the state approaches the stable steady-state or the stable periodic orbit. Eventually, the state reaches the stable periodic orbit at the time instant marked with the green dashed line in Fig. 3b). The state stays at the orbit unless perturbed; this stage is called the self-pulsing stage and takes place at the edge of the yellow surface.

Thus, it is found that once the state reaches the yellow surface, the state is confined to the surface. The phenomenon of data collapse to a surface happens also for the other two variables, $|a_{+}|$, $\Delta\phi$. Since the state always lies on the surface after a very short initial time, we conclude that two degrees of freedom are sufficient to specify the state after the initial transition stage and the propagation of the state is locally restricted to two directions. We note that this phenomenon is a general feature of a dynamical system close to a Hopf bifurcation and is called a reduction to the center manifold kuznetsov_elements_2004 ; seydel_practical_2010 ; guckenheimer_nonlinear_1983 . The dimension of the center manifold is strictly related to the number of steady-state eigenvalues, the real parts of which cross zero kuznetsov_elements_2004 ; seydel_practical_2010 ; guckenheimer_nonlinear_1983 . In Fig 2a), we have seen that in the present case there are two eigenvalues with real parts crossing zero, while all the remaining eigenvalues have negative real parts giving rise to a stable manifold. Thus, the center manifold is two-dimensional as it is confirmed by the yellow curved surface in Fig. 3a). The dynamics in the remaining three directions quickly approach the surface during the initial transition stage.

In Fig. 3, $n_{w}$ and $|a_{-}|$ are the two degrees of freedom, while all the remaining degrees of freedom {$|a_{+}|$, $\Delta\phi$, $n_{c}$} of the state vector $\vec{\psi}$ are expressed as functions of the variables $n_{w}$ and $|a_{-}|$ after the initial transition stage:

[TABLE]

Similarly, equations for each component of the velocity vector $d\vec{\psi}/d\tau$ in Eqs. (2) and (3b) can be expressed as functions of $n_{w}$ and $|a_{-}|$:

[TABLE]

By taking the total derivative of $d\vec{\psi}(\tau)/d\tau$, we obtain:

[TABLE]

which, when compared with Eq. (4), indicates that the laser dynamics can be locally approximated by a $2\times 2$ Jacobian matrix $\bm{A}$.

The functions of $n_{w}$ and $|a_{-}|$ in Eq. (5a) and the yellow surface in Fig. 3 are approximated by polynomials. In order to do that, we solve the system of differential equations in Eqs. (2) for varying initial conditions. Each solution then corresponds to a different trajectory plotted versus $n_{w}$ and $|a_{-}|$. All of these trajectories are seen to lie on the surface after the initial transition stage similarly to Fig. 3. Next, we fit a polynomial with all the trajectories excluding the initial transition stage. Then the polynomial describes the surface $n_{c}$ in terms of $n_{w}$ and $\mathopen{|}a_{-}\mathclose{|}$. Similarly, we can approximate the surfaces for $\mathopen{|}a_{+}\mathclose{|}(\mathopen{|}a_{-}\mathclose{|},n_{w})$ and $\Delta\phi(\mathopen{|}a_{-}\mathclose{|},n_{w})$. We emphasize that in order to keep the original coordinate system of the variables, we exclude the part of the trajectory in the initial transition stage and fit a polynomial with the remaining parts of all the trajectories. Thus, we fit the polynomials once the state has reached the center manifold. Having obtained these surfaces, we can determine any state in the phase space within the periodic orbit once $n_{w}$ and $|a_{-}|$ are known without any need of solving the five-dimensional system of equations in Eq. (2).

Moreover, in order to describe the dynamics on these surfaces, we need to solve a system of two differential equations describing the two degrees of freedom, $n_{w}$ and $|a_{-}|$. These equations are the components of the velocity vector in Eq. (III.3):

[TABLE]

Once Eq. (7) is solved, the remaining degrees of freedom {$|a_{+}|$, $\Delta\phi$, $n_{c}$} can be determined using the polynomials.

III.4 Instantaneous eigenvalues

We compute the instantaneous eigenvalues of the matrix $\bm{A}$ in Eq. (4) for the state vectors $\vec{\psi}$, Eq. (5a) over the whole surfaces. In Fig. 3, we have seen that we can define three surfaces for $|a_{+}|$, $\Delta\phi$, $n_{c}$ plotted versus $n_{w}$ and $|a_{-}|$. After the initial transition stage, the state always lies on these surfaces. Thus, all the dynamics are confined to these surfaces. Then, each point of these surfaces can be substituted into the matrix $\bm{A}$ and the instantaneous eigenvalues of the state at this position are obtained. Figure 4 shows the real and imaginary parts of the three instantaneous eigenvalues with the largest real parts over the whole surface as well as along the trajectory marked with the green dashed line. The two remaining eigenvalues have significantly smaller real parts, and thus are not included in Fig. 4 as the contribution from the corresponding eigenvectors decays rapidly.

Figure 4a) shows that the pair of eigenvalues marked with orange and yellow has considerably larger real parts than the third eigenvalue (purple) over the major part of the surface. The third eigenvalue is only comparable to the other two eigenvalues along the line $|a_{-}|=0$, but it is still smaller and never becomes positive. The negative real parts of the purely real third eigenvalue (purple) and the remaining complex conjugate pair of eigenvalues (not shown) signifies that the contribution of the corresponding eigenvectors in a reconstruction of the solution decays very quickly. This is what is observed in Fig. 3 in the initial transition stage. Afterwards, once the state is on the surface, the contribution from the three corresponding eigenvectors is negligible and the state description is dominated by the eigenvectors corresponding to the two eigenvalues with the largest real part (orange and yellow).

The real parts of the pair of eigenvalues marked with orange and yellow are seen to dominate for large values of $|a_{-}|$; this is where the pulse is released. In Fig. 4c),d,e)) we show the instantaneous eigenvalues in the vicinity of the pulse along the green trajectory in a), b) when the state has already reached a limit cycle. On the right axis we plot the pulse power in the straight-port defined as:

[TABLE]

where $c_{0}$ is the speed of light and $\epsilon_{0}$ is vacuum permittivity.

Figure 4 shows that when the state moves along the $n_{w}$ axis (just after the previous pulse has been released and before a new one) the three eigenvalues with the largest real parts are purely real. As the limiting value of $n_{w}$ is reached, Fig. 4a), one of the eigenvalues (orange) starts to rapidly increase, Fig. 4c), while the third eigenvalue (purple) drops rapidly. Just as the pulse is released, the second eigenvalue (yellow) rapidly increases and collapses with the first one (yellow) at the exceptional point. Therefore, it is found that, as the pulse grows, the pair of eigenvalues transitions from being purely real to being complex conjugate when crossing the exceptional point. As the pulse power decreases, the complex conjugate pair of eigenvalues coalesce at the second exceptional point and transitions back to the pair of purely real eigenvalues. Thus, most of the pulse is observed to be bounded by the two instantaneous exceptional points with positive/negative real part of the eigenvalue at the beginning/end of the pulse, respectively. Interestingly, two more exceptional points are observed as the pulse is decaying, see Fig. 4e).

Within one period the state traverses a loop in the phase space of the model. We have seen that four exceptional points are crossed within a single loop when the laser state is a periodic orbit. When the exceptional point is approached, the eigenvectors exhibit a characteristic phase jump and are phase shifted relative to each other by $\pm i$ gunther_projective_2007 ; heiss_repulsion_2000 . Therefore, during an evolution along any trajectory in the diminishing vicinity of an exceptional point, eigenvectors will acquire a phase shift of $\pm i$ keck_unfolding_2003 ; rotter_non-hermitian_2009 ; muller_exceptional_2008 . In menke_state_2016 , it has been shown that this effect is preserved as long as the exceptional point is inside the loop or crossed by it. Therefore, it is only a four-fold loop around an exceptional point or a single loop around four exceptional points that will restore an original scenario for the eigenvectors concerned heiss_collectivity_1998 ; heiss_phases_1999 ; dembowski_experimental_2001 ; heiss_exceptional_2004 ; heiss_physics_2012 . Since the laser is operating in the periodic orbit in our case, in order to remain periodic it has to cross four exceptional points within one period in phase space.

III.5 Reconstruction of the solution

At most two out of five instantaneous eigenvalues have positive real parts. Thus, after the initial transition stage, the eigenvectors which correspond to the two dominating eigenvalues can be used to reconstruct the solution of Eq. (4) as follows:

[TABLE]

where $\vec{v}_{1,2}$ are the instantaneous right eigenvectors, and $c_{1,2}$ are the amplitudes of the corresponding eigenvectors. These amplitudes can be reconstructed from a solution $d\vec{\psi}(\tau)/d\tau$ using the left eigenvectors as $c_{1,2}=\vec{w}_{1,2}^{T}(\tau)d\vec{\psi}(\tau)/d\tau$. In the following we show that the two eigenvectors can be used to approximate the two tangential vectors to the surface pointing along the $|a_{-}|$, $n_{w}$ coordinate lines. This confirms that the solution can be approximately expanded in the two eigenvectors.

The tangential vector to the surface $z=f(x,y)$ along the parametric curve $\vec{r}(t)=\{x(t),y(t),z(t)\}$ on this surface is expressed as $\vec{r}^{\prime}(t)=\{x^{\prime}(t),y^{\prime}(t),z^{\prime}(t)\}$, where $z^{\prime}(t)=\nabla f\cdot\vec{u}$, $\vec{u}=\{x^{\prime}(t),y^{\prime}(t)\}$. In our case, the tangential vectors to the surfaces, which approximate the components of the state vector $\vec{\psi}$, Eq. (5a), are expressed as:

[TABLE]

where

[TABLE]

The tangential vectors $\vec{r}^{\prime}(\tau)$ along the parameterized trajectory can be decomposed into a linear combination of the tangential vectors to the surface along its coordinates $|a_{-}|$, $n_{w}$:

[TABLE]

It is observed that the tangential vectors to the surface are composed of the components of the velocity vector, Eq. (10), and the velocity vector can be expanded into the two eigenvectors, see Eq. (9). Since the five-dimensional (5D) matrix $\bm{A}$ is real, the top two eigenvalues ($\lambda_{1}$ and $\lambda_{2}$) and eigenvectors ($\vec{v_{1}}$ and $\vec{v_{2}}$) are either real or form a complex conjugate pair. Then, we change these eigenvectors to point along the original coordinate lines, $|a_{-}|$, $n_{w}$, as follows:

[TABLE]

Then, the two vectors $\vec{v}_{1}^{\prime}$ and $\vec{v}_{2}^{\prime}$ are determined at positions of the state vector approximated by the polynomials, see Eq. (5a), and separated by equidistant steps. The vectors are purely real and are plotted over the whole surfaces $|a_{+}|$, $\Delta\phi$, $n_{c}$, see Fig. 5. Subsequently, these vectors are scaled by the distance between the steps in the state vector along each direction in order to avoid an overlap and create a square grid pattern. If these vectors create an ideal square grid then they can perfectly reconstruct the tangential vectors in Eq. (12). A small discrepancy is only found in Fig. 5(c) for small values of $|a_{-}|$, which can be explained by the third eigenvalue becoming comparable to the dominating pair of eigenvalues at these points, see Fig. 4. However, the two vectors $\vec{v}_{1}^{\prime}$ and $\vec{v}_{2}^{\prime}$ are found to approximate the tangential vectors over the whole surface as observed in Figure 5. Thus, the two degree-of-freedom picture is justified over the whole surface and shown to precisely reconstruct $d\vec{\psi}(\tau)/d\tau$. Therefore, the system of five nonlinear differential equations can be reduced to only two differential equations after the initial transition stage. The other three dimensions are functions of $n_{w}$ and $|a_{-}|$ and are presently approximated by polynomials. We note that the instantaneous eigenvalues/eigenvectors are not needed to reduce dimensionality of the system, but they provide an additional insight into the solution. Furthermore, the fact that the two instantaneous eigenvectors approximate the tangential vectors to the surfaces proves that the system dimensionality can be reduced to two.

Moreover, we note that although a 2D model can be used to describe the laser dynamics after the initial transition stage, there exists a parameter region in which even a 1D model is sufficient to replace the original 5D model after the initial transition stage. One may observe in Fig. 2 that for a large negative detuning $\Delta\omega_{c}$, the steady-state eigenvalues undergo transition from a complex conjugate pair of eigenvalues to two purely real eigenvalues. Then, one of the eigenvalues decreases rapidly, and the other one approaches zero. Therefore, for detunings $-2.05\gamma_{T}<\Delta\omega_{c}<-1.72\gamma_{T}$, there is a single steady-state eigenvalue that dominates and thus, the velocity vector can be described by a single eigenvector, see Eq. (9). In this case, the laser dynamics can be described by a single differential equation after the transition stage in which the contribution from the other four eigenvalues rapidly decays. For detunings $\Delta\omega_{c}<-2.05\gamma_{T}$, the lasing mode ceases to exist mork_chaos_1992 ; mork_photonic_2014 ; rasmussen_theory_2017 . Thus, as the detuning $\Delta\omega_{c}$ increases, the steady-state eigenvalues transition from being purely real to a complex conjugate pair and the system evolves from a 1D to a 2D system.

IV Origin of the laser instability

IV.1 Detection of periodic orbits

In what follows, we use the simplified 2D model, Eq. 7, to explain the origin of the laser instability that may be observed even when all real parts of the steady-state eigenvalues are negative.

At first, the phase space of the Fano laser is scanned in search for periodic orbits. We choose our initial conditions as follows: $1)$ $|a_{-}|^{\mathrm{initial}}$ is set to the steady-state value and $2)$ $n_{w}^{\mathrm{initial}}$ is varied over the whole phase space along the purple line as shown in Fig. 6b). For each initial condition, we then compute the trajectory by solving Eq. (7) up to the point when $|a_{-}|(T)=|a_{-}|^{\mathrm{initial}}$, where $T$ is the time corresponding to one cycle. Some of these trajectories are shown in Fig. 6(b) in different colours. Subsequently, we evaluate the shift $\Delta x=n_{w}^{\mathrm{initial}}-n_{w}^{\mathrm{cycle}}$ in the state vector after the time $T$. If the shift between the initial state and the state after one cycle is zero then we are at a periodic orbit or steady-state. On the other hand if it is non-zero it means that the state is approaching or departing from the steady-state/periodic orbit.

Figure 6a) shows the numerically evaluated shift $\Delta x$ in $n_{w}$ after one cycle of the trajectory. It is seen that there are five crossings with zero. These crossings are marked with blue, red and black dots. The single black dot indicates the steady-state, while the pairs of blue and red dots indicate periodic orbits. The outer periodic orbit, marked with a pair of blue dots, has been observed before and is known to be stable for a pair of complex conjugate steady-state eigenvalues with a positive real part rasmussen_theory_2017 . Here, it is seen that for a strong enough perturbation of the initial conditions from the steady-state value, the state can still reach the outer periodic orbit despite all the steady-state eigenvalues having a negative real part and thus the steady-state being stable and attracting the state. Furthermore, we find an additional periodic orbit marked with a pair of red dots in Fig. 6a). The newly found periodic orbit separates the steady-state and the outer periodic orbit.

IV.2 Stability of the orbits

We now prove the stability of the newly found orbit using the simplified 2D model. This is done by calculating the Floquet multipliers, $\lambda_{f}$, which tell us how the solution behaves in the vicinity of the periodic orbit, i.e., whether it diverges/converges from/towards the orbit iooss_elementary_1997 ; glendinning_stability_1994 . In order to compute the Floquet multipliers, we first obtain the fundamental solution matrix, $\bm{\Phi}(\tau)$, which can be determined using $d\bm{\Phi}(\tau)/d\tau=\bm{A}(\tau)\bm{\Phi}(\tau)$ and satisfies d\vec{\psi}/d\tau=\bm{\Phi}(\tau)d\vec{\psi}/d\tau\bigl{|}_{\tau=0} with $\bm{\Phi}(0)=\bm{I}$. The Floquet multipliers are the eigenvalues of $\bm{\Phi}(\tau)$ evaluated at $\tau=T$, where $T$ is the period of the orbit.

If the Floquet multipliers are within the unit circle in the complex plane, the orbit is stable, otherwise it is unstable. The Floquet multipliers for the outer periodic orbit are $\lambda_{f1}=0.04$ and $\lambda_{f2}=1$ confirming its stability. On the other hand, the Floquet multipliers of the newly found orbit are $\lambda_{f1}=2.31$ and $\lambda_{f2}=1$, proving that this orbit is unstable. We note that for a periodic orbit there is always one of the Floquet multipliers for which $\lambda_{f}=1$ and the corresponding eigenvector is tangential to the periodic orbit. This neutral stability accounts for the possibility of drift along the periodic orbit glendinning_stability_1994 .

Furthermore, we study the stability of the newly found orbit with variation of the detuning, $\Delta\omega_{c}$. It is observed that the Floquet multiplier crosses the unit circle along the real axis in the complex plane. This indicates an exchange of instability seydel_practical_2010 . Indeed, as $\Delta\omega_{c}$ decreases from $\Delta\omega_{c}=0.52\gamma_{T}$, see Fig. 2, the newly found unstable periodic orbit increases in size. Eventually, it collapses with the stable periodic orbit. Both orbits disappear due to a fold bifurcation kuznetsov_elements_2004 and only the stable (time-independent) steady-state remains present in the phase space. On the other hand, when $\Delta\omega_{c}$ increases, the new unstable orbit decreases in size. Eventually it collapses with the stable steady-state resulting in the steady-state becoming unstable. This happens in the vicinity of $\Delta\omega_{c}=1.52\gamma_{T}$, which is the critical bifurcation point and as $\Delta\omega_{c}$ is increased further, the real part of the steady-state eigenvalues becomes positive. Since the cycle is present before the bifurcation point, i.e. before the real parts of the steady-state eigenvalues become positive, the bifurcation at this point is called a subcritical Hopf bifurcation.

Both bifurcations are marked in the phase diagram in Fig. 7. It shows that as $\Delta\omega_{c}$ is decreased from large values along the dark green arrow, at first the system undergoes a supercritical Hopf bifurcation at the dashed yellow line. There, the real part of the steady-state eigenvalues crosses zero and becomes positive resulting in the steady-state point becoming unstable and a stable limit cycle being present after the bifurcation point. As we decrease $\Delta\omega_{c}$ further, the system undergoes a subcritical Hopf bifurcation at the dotted red line. Here, the steady-state point becomes stable again, while the stable periodic orbit coexists with an unstable periodic orbit. Eventually, as $\Delta\omega_{c}$ is further decreased, the unstable periodic orbit collides with the stable one, and both orbits disappear through a fold bifurcation leaving the stable steady-state point as the only solution kuznetsov_elements_2004 ; seydel_practical_2010 . Analogous behaviour is observed for lower pump currents including $J=1.2J_{\mathrm{thr}}$, however in this case the laser is below threshold for large detuning $\Delta\omega_{c}$.

Occurrence of both Hopf bifurcations, supercritical and subcritical, is a signature of a Bautin bifurcation, also known as generalized Hopf bifurcation govaerts_numerical_2000 ; kuznetsov_elements_2004 . A Bautin bifurcation is characterized by the presence of two orbits and an equilibrium point (steady-state) in phase space. We note that a Bautin bifurcation cannot be detected by merely monitoring the eigenvalues govaerts_numerical_2000 ; kuznetsov_elements_2004 . Upon the external parameter variation, $\Delta\omega_{c}$, an inner orbit may collide with an outer orbit and annihilate or exchange stability with an equilibrium point, as it has been observed. We note that since each $\Delta\omega_{c}$ results in different solutions $(\omega_{s},N_{s})$ of the oscillation condition in Eq. (1), we adjust the polynomial approximation of the surface in each case.

The stability of the orbits can also be assessed based on Fig. 6a). It is seen that if the model is initialized outside the orbit marked with the red dots, $n_{w}^{\mathrm{cycle}}$ will increase after one cycle compared to the initial value. Thus, the state is repelled away from the orbit. If the model is initialized inside the orbit, $n_{w}^{\mathrm{cycle}}$ will keep increasing with each cycle confirming that the newly found orbit is unstable.

V Conclusion

We demonstrate that after a fast initial transient the dynamics of the recently realized Fano laser yu_demonstration_2017 are confined to a 2D center manifold. The dimension of the center manifold follows the number of steady-state eigenvalues, the real parts of which cross zero. We show that there are two steady-state eigenvalues with real parts crossing zero, while the remaining three eigenvalues have negative real parts forming a stable manifold. The dynamics is attracting along the corresponding three directions and quickly tends to the curved surface, i.e. the center manifold, during the initial transition stage. Afterwards, the state vector is confined to the curved surface and can be solely described by two degrees of freedom. The surface geometry of the phase space can be approximated by the two eigenvectors of the linear stability matrix corresponding to the eigenvalues with the largest real parts. As the pulse develops, the instantaneous eigenvalues transition from a pair of purely real eigenvalues to a complex conjugate pair at the first exceptional point. The main part of the repeating pulse is bounded by two exceptional points with positive/negative real part of the eigenvalue at the beginning/end of the pulse, respectively. Moreover, the trajectory encounters four exceptional points during one period, ensuring that both, the eigenvalues and the eigenvectors, are periodic in $\tau$. Furthermore, we show that the 5D model used to describe the laser dynamics, after the initial transition stage, can be reduced to only 1D in part of the parameter space and evolves into a 2D model beyond the exceptional point of steady-state eigenvalues as the detuning $\Delta\omega_{c}$ increases. Moreover, we have used the simplified 2D model to associate the unknown source of laser instability with the newly found unstable periodic orbit, which arises due to a generalized Hopf (Bautin) bifurcation. These findings allow to better understand the laser dynamics and may lead to the design of new functionalities in nanolasers used for on-chip communications and sampling.

Acknowledgements.

The authors would like to thank T. S. Rasmussen for helpful discussions on the implementation of the dynamic model. This work was supported by Villum Fonden via the Centre of Excellence NATEC (grant 8692) and Research Grants Council of Hong Kong through project C6013-18GF.