Dynamics of a class A nonlinear mirror mode-locked laser

TL;DR

This paper models the complex dynamics of a class-A nonlinear mirror mode-locked laser using delay differential equations, revealing stability conditions, instabilities, and the formation of asymmetric mode-locked pulses.

Contribution

It introduces a theoretical delay differential equation model to analyze the stability and pulse formation in a class-A nonlinear mirror mode-locked laser.

Findings

CW regimes can be destabilized by modulational, Turing, and square-wave instabilities.

Mode-locked pulses exhibit asymmetric exponential decay in their edges.

Harmonic mode-locking arises from asymmetric pulse interactions.

Abstract

Using a delay differential equation model we study theoretically the dynamics of a unidirectional class-A ring laser with a nonlinear amplifying loop mirror. We perform linear stability analysis of the CW regimes in the large delay limit and demonstrate that these regimes can be destabilized via modulational and Turing-type instabilities, as well as by an instability leading to the appearance of square-waves. We investigate the formation of square-waves and mode-locked pulses in the system. We show that mode-locked pulses are asymmetric with exponential decay of the trailing edge in positive time and faster-than-exponential (super-exponential) decay of the leading edge in negative time. We discuss asymmetric interaction of these pulses leading to a formation of harmonic mode-locked regimes.

Figures

Equations

$${\gamma }^{-1}\frac{dE}{dt}+E=\sqrt{{\kappa }_1}{e}^{g/2+i\theta }\tilde{f}\left[|E(t-T){|}^2,G\right]E(t-T),$$

$${\gamma }_g^{-1}\frac{dg}{dt}=g_0-g-\left(e^g-1\right)|E(t-T){|}^2.$$

$$\tilde{f}(|E{|}^2,G)=\sqrt{{\kappa }_2}{e}^{G/2}f(|E{|}^2),$$

$$f(|E{|}^2)=\frac{1}{2}\left(1-e^{ia|E{|}^2}\right),$$

$${\gamma }^{-1}\frac{dE}{dt}+E=\sqrt{\kappa }{e}^{g_0/[2(1+|E(t-T){|}^2)]+i\theta }f(|E(t-T){|}^2)E(t-T),$$

$$E=\sqrt{R}{e}^{i\Omega t},$$

$$\kappa e^{\frac{g_0}{1+R}}{\sin }^2\left(\frac{aR}{2}\right)=1+\frac{{\Omega }^2}{{\gamma }^2},$$

$$\tan \left(\frac{a{R}^2}{2}-T\Omega +\theta \right)+\frac{\gamma }{\Omega }=0.$$

$$g_0=\tilde{g}\equiv a(1+R{)}^2\cot \left(\frac{aR}{2}\right).$$

$$c_2{Y}^2+2c_1(\lambda )Y+c_0(\lambda )=0,$$

$$\lambda =i{\lambda }_0+\frac{{\lambda }_1}{T}+\mathcal{O}\left(\frac{1}{T^2}\right),{\lambda }_1={\lambda }_{11}+i{\lambda }_{12},$$

$${\lambda }_{11}^{\pm }=\mathrm{Re}{\lambda }_1^{\pm }=\mathrm{Re}\ln \left(Y_{\pm }^{-1}\right),$$

$$Y_{\pm }=\frac{-c_1(i{\lambda }_0)\pm \sqrt{c_1(i{\lambda }_0{)}^2-c_0(i{\lambda }_0)c_2}}{c_2}.$$

$${\left(\frac{d^2{\lambda }_{11}^{-}}{d{\lambda }_0^2}\right)}_{{\lambda }_0=0}=0.$$

$$Y_{+}^{-1}{|}_{{\lambda }_0=0}=-1,$$

$$2-\frac{g_{0}R}{(1+R{)}^2}+aR\cot \left(\frac{aR}{2}\right)=0.$$

$${\Omega }^2\pm \gamma \Omega aR={\gamma }^2.$$

$$Y_{-}{|}_{\Omega =0}$$

$$Y_{+}{|}_{\Omega =0}$$

$$\epsilon \frac{dE}{d\tau }+E=\sqrt{\kappa }{e}^{g_0/[2(1+|E(\tau -1){|}^2)]}f\left[|E(\tau -1){|}^2\right]E(\tau -1),$$

$$R_{n+1}=h(R_n),h(R_n)=\kappa e^{\frac{g_0}{1+R_n}}{\sin }^2\left(\frac{a{R}_n}{2}\right)R_n,$$

$$\kappa e^{\frac{g_0}{1+R^{\ast }}}{\sin }^2\left(\frac{a{R}^{\ast }}{2}\right)=1.$$

$$1-\frac{g_0{R}^{\ast }}{(1+R^{\ast }{)}^2}+a{R}^{\ast }\cot \left(\frac{a{R}^{\ast }}{2}\right)=-1,$$

$$E(t-T)=E(t-T+T_0)e^{i\Delta }\equiv E(t+\delta )e^{i\Delta }$$

$${\gamma }^{-1}\frac{dE}{dt}+E=\frac{\sqrt{\kappa }}{2}{e}^{g_0/[2(1+|E(t+\delta ){|}^2)]+i\Theta }f(|E(t+\delta ){|}^2)E(t+\delta ),$$

$${\gamma }^{-1}\frac{dE}{dt}+E=ϵE(t+\delta )e^{i\Theta },$$

$${\lambda }_k=-\gamma \left[1+\frac{W_k(-ϵ\gamma \delta e^{-\gamma \delta +i\Theta })}{\gamma \delta }\right],$$

$$c_0(\lambda )$$

$$c_1(\lambda )$$

$$c_2$$

$$\frac{{\left[a(1+R{)}^2\Omega +\gamma (\tilde{g}-g_0)\right]}^2}{{\left({\gamma }^2+{\Omega }^2\right)}^2}+\frac{2(1+R{)}^2{\Omega }^2\left[g_0^2+2a^2(1+R{)}^4-g_0^2\cos (aR)-2a{g}_0(1+R{)}^2\sin (aR)\right]}{R(\tilde{g}-g_0){\left({\gamma }^2+{\Omega }^2\right)}^2\left[1-\cos \left(aR\right)\right]}=0.$$

Full text

I Introduction

The possibility of generation of short light pulses by locking the longitudinal modes of a laser was discussed only a few years after the development of the laser in 1960. Mode-locking techniques can be classified into two major classes: (i) active mode-locking, based on an external modulation at a frequency close to the cavity free spectral range and (ii) passive mode-locking where an intracavity nonlinear component reduces losses for pulsed operation with respect to the those of continuous-wave (CW) regime. A standard theoretical approach to study the properties of mode-locked devices is based on direct integration of the so-called traveling wave equations describing space-time evolution of the electric field and carrier density in the laser sections Tromborg et al. (1994); Avrutin et al. (2000); Bandelow et al. (2001, 2006). Another, much simpler approach limited to small gain and loss approximation was developed by Haus in Haus (2000). To overcome this limitation the third modelling approach was suggested in Vladimirov and Turaev (2004); Vladimirov et al. (2004); Vladimirov and Turaev (2005) based on the lumped element method that allows to derive a delay differential equation (DDE) model for the temporal evolution of the optical field at some fixed position in the cavity. DDE models successfully describe the dynamics of multi-section mode-locked semiconductor lasers Viktorov et al. (2006); Rachinskii et al. (2006); Vladimirov et al. (2010); Arkhipov et al. (2013); Marconi et al. (2014); Viktorov et al. (2014), frequency swept light sources Slepneva et al. (2013, 2014); Pimenov et al. (2017), optically injected lasers Rebrova et al. (2011); Pimenov et al. (2014), semiconductor lasers with feedback Kelleher et al. (2017); Otto et al. (2012); Jaurigue et al. (2017), as well as some other multimode laser devices Viktorov et al. (2007).

In this work, we propose and study a DDE model for nonlinear amplifying loop mirror (NALM) mode-locked laser. Nonlinear mirror laser as a device for ultrafast light processing was proposed in Doran and Wood (1988). Mode-locked pulse formation mechanism in this laser is based on the asymmetric nonlinear propagation in a waveguide loop, where two counterpropagating waves aquire different intensity-dependent phase shifts caused by the Kerr nonlinearity. As a result of the interference of these waves the loop acts as a nonlinear mirror with the reflectivity dependent on the intensity of incident light. Such nonlinear mirror plays a role of a saturable absorber that leads to the appearance of a pulsed laser operation.

We develop a DDE laser model by assuming dispersion-free unidirectional operation inside the cavity and symmetrical beam splitting of the field into two counter-propagating fields at the entrance of the NALM. When the material variables are adiabatically eliminated, one obtains a single DDE for the complex electric field envelope, which, despite of its simplicity, gives a good insight on the laser dynamics. Using this equation we find different mode-locked regimes of laser operation including square waves, ultrashort optical pulses and their harmonics. The mode-locked pulses are always bistable with the non-lasing state, and, hence, can also be considered as temporal cavity solitons or nonlinear localized structures of light Leo et al. (2010); Marconi et al. (2014). The linear stability analysis of CW solutions corresponding to different longitudinal modes of the laser reveals modulational and Turing instabilities, as well as an instability leading to the emergence of square waves. We demonstrate analytically the mode-locked pulses are asymmetric with exponential decay of the pulse trailing tail in positive time and faster-than-exponential (super-exponential) decay of the leading tail in negative time. We find that the repulsive interaction of such asymmetric pulses leads to a harmonic mode-locking regime. In order to explain the experimentally observed square wave generation recently reported in Aadhi et al. (2018) in the figure-of-eight laser, we construct a one-dimensional map exhibiting a period doubling route to chaos.

II Model equation and CW solutions

A schematic of a NALM laser with a gain and a spectral filter in an unidirectional cavity coupled to a bidirectional loop with a second gain medium and a nonlinear element, is shown in Fig. 1. This scheme corresponds to experimentally implemented setups of mode-locked lasers with a high-Q microring resonator Kues et al. (2017) and an integrated waveguide Aadhi et al. (2018) acting as nonlinear elements. Using the approach of Vladimirov et al. (2004); Vladimirov and Turaev (2005) we write the following DDE model for the time evolution of the electric field amplitude $E$ and saturable gain $g$ in a laser shown in Fig. 1:

[TABLE]

[TABLE]

Here $\gamma$ is the spectral filter bandwidth, $\kappa_{1}$ is the linear attenuation factor describing nonresonant cavity losses in the unidirectional part of the laser cavity, $g_{0}$ is the pump parameter, $\gamma_{g}$ is the relaxation rate of the amplifying medium, $\theta$ describes the detuning between the central frequency of the filter and one of the cavity modes, and the delay parameter $T$ is equal to the cold cavity round trip time.

The function ${\tilde{f}}\left(|E|^{2},G\right)$ in Eq. (1) describes the reflectivity and the phase shift introduced by the NALM. It is given by ${\tilde{f}}\left(|E|^{2},G\right)=\sqrt{\kappa_{2}}e^{G/2}\left(e^{i\phi_{1}}-e^{i\phi_{2}}\right)/2$, where $G$ and $\kappa_{2}$ and describe amplification and linear losses inside the loop, while $\phi_{1}=\eta|E|^{2}/2$ and $\phi_{2}=\eta e^{G}|E|^{2}/2$ are the phase shifts of the clockwise and counter-clockwise propagating waves indicated by arrows in Fig. 1. Since the counter-clockwise propagating wave is amplified before passing through the nonlinear element its phase shift $\phi_{2}$ is $e^{G}$ times larger than that the clockwise propagating wave, which is amplified after passing through the nonlinear element. For the gain $G$ an equation similar to Eq. (2) can be written. Here, however, we assume for simplicity that the electric field intensity $|E|^{2}$ is small and gain medium inside the NALM operates far below the saturation regime. In this case $G$ can be considered to be a constant parameter. Furthermore, for $\phi_{1}=\eta|E|^{2}/2\ll 1$ we have $e^{\phi_{1}}\approx 1$ and the function ${\tilde{f}}$ can be rewritten in the form

[TABLE]

with

[TABLE]

where $a=\eta(e^{G}-1)/2$. Finally, replacing $e^{g}-1$ with $g$ in Eq. (2), eliminating the gain $g$ adiabatically, $g=g_{0}/\left[1+|E(t-T)|^{2}\right]$, and substituting this expression into Eq. (1) we get

[TABLE]

where $\kappa=\kappa_{1}\kappa_{2}e^{G/2}<1$, which means that the linear gain in the NALM is compensated by the linear cavity losses, and the nonlinear function $f(|E|^{2})$ is defined by Eq. (3). In the following we will show that despite being very simple, the model equation (4) is capable of reproducing such experimentally observed behaviours of nonlinear mirror lasers as mode-locking and square wave generation Aadhi et al. (2018).

The simplest solution of Eq. (4) is that corresponding to laser off state, $E=0$. Linear stability analysis of this trivial non-lasing solution indicates that it is always stable, which means that the laser is non-self-starting for all possible parameter values. Non-trivial CW solutions are defined by the relation

[TABLE]

where $R=|E|^{2}>0$ is the intensity and $\Omega$ is the frequency detuning of CW regime from the reference frequency coinciding with the central frequency of the spectral filter. $R$ and $\Omega$ are the solutions of a system of two transcendental equations

[TABLE]

[TABLE]

Multiple solutions of these equations corresponding to different longitudinal modes of the laser are illustrated in the left panel of Fig 2, where CW regimes correspond to the intersections of black closed curves obtained by solving Eq. (6) and thin gray lines calculated from Eq. (7). Black curves in right panel of Fig. 2 show the branches of non-trivial CW solutions as functions of the pump parameter $g_{0}$, while gray lines are defined by the condition

[TABLE]

The intersections of the gray lines with black CW branches indicate fold bifurcation points, each of which separates the corresponding branch into two parts with the lower part being always unstable.

III CW stability in the large delay limit

To study linear stability the CW solutions of Eq. (4) we apply the large delay limit approach described in Ref. Yanchuk and Wolfrum (2010). By linearizing Eq. (4) on a CW solution given by Eq. (5) we obtain the characteristic equation in the form:

[TABLE]

where $Y=e^{-\lambda T}$ and $c_{0,1,2}$ are given in the Appendix. In the large delay limit, $\gamma T\gg 1$, the eigenvalues belonging to the so-called pseudo-continuous spectrum can be represented in the form:

[TABLE]

with real $\lambda_{0}$, $\lambda_{11}$, and $\lambda_{12}$. Therefore, using the approximate relations $c_{0,1}(\lambda)\approx c_{0,1}(i\lambda_{0})$ and $Y\approx e^{-i\lambda_{0}T-\lambda_{1}}$ we can solve the characteristic equation to express $\lambda_{11}$ as a function of $\lambda_{0}$ Yanchuk and Wolfrum (2010)

[TABLE]

Two solutions $\lambda_{11}^{\pm}(\lambda_{0})$ given by Eq. (10) define two branches of pseudo-continuous spectrum shown in the left panel of Fig. 3. Due to the phase shift symmetry of the model equation (4) $E\to Ee^{i\varphi}$ with arbitrary constant $\varphi$ one of these branches is tangent to the $\lambda_{11}=0$ axis on the ($\lambda_{0}$,$\lambda_{11}$)-plane at the point $\lambda_{0}=0$, i.e., $Y_{-}^{-1}\rvert_{\lambda_{0}=0}=1$. The right panel of Fig. 3 shows the stability diagram of the CW solutions on the plane of two parameters: normalized frequency offset $\Omega/\gamma$ of a CW solution and pump rate $g_{0}$. CW solutions are stable (unstable) in the dark (light) gray domains. Black curve indicates the fold bifurcation, where two CW solutions merge and disappear. Below this curve calculated from Eqs. (6) and (8) there are no CW solutions while above it a pair of CW solutions is born with one them corresponding to smaller intensity $R$ being always unstable. Similarly to the case of Eckhaus instability Tuckerman and Barkley (1990); Wolfrum and Yanchuk (2006) the upper branch of CW solutions can be stable only within the so-called Busse balloon (region labelled “1” in Fig. 3) limited from below by modulational instability curve shown by blue line. The modulational instability curve is defined by the condition that one of the two branches of the pseudo-continuous spectrum, which satisfies the relation $\left(\lambda_{11}^{-}\right)_{\lambda_{0}=0}=\left(d\lambda_{11}^{-}/d\lambda_{0}\right)_{\lambda_{0}=0}=0$, changes the sign of its curvature at the point $\lambda_{0}=0$:

[TABLE]

This instability is illustrated sub-panels 2 and 6 of the left panel of Fig. 3. Explicit expression for the modulational instability condition (11) in terms of the model equation parameters is given in the Appendix.

It is seen from the right panel of Fig. 3 that the modulational instability curve is tangent to the fold bifurcation curve at $\Omega=0$ and that this curve becomes asymmetric with respect to the axis $\Omega=0$ sufficiently far away from the tangency point.

The upper boundary of the CW stability domain shown in the right panel of Fig. 3 consists of two parts separated by codimension-two point C+. The right part of this boundary lying between the points C+ and C is indicated by red line and corresponds to the so-called Turing-type (wave) instability Yanchuk and Wolfrum (2010), where one of the two branches of pseudo-continuous spectrum cross the $\lambda_{11}=0$ axis at two symmetric points with $\lambda_{0}\neq 0$ i.e., $\lambda_{11}^{+}=0$, at $\lambda_{0}=\pm\lambda_{0}^{*}$ with $\lambda_{0}^{*}>0$, see sub-panel 5 of the left panel of Fig. 3. The left part of the stability boundary lying between two symmetric codimension-two points C*±* corresponds to a flip instability leading to the appearance of square waves with the period close to $2T$. In the right panel of Fig. 3 the flip instability curve is shown by green line defined by the condition

[TABLE]

which can be rewritten in the form:

[TABLE]

The codimension-two points C*±* are defined by Eq. (13) together with additional conditions $\left(d^{2}\lambda_{11}^{\pm}/d\lambda_{0}^{2}\right)_{\lambda_{0}=0}=0$. Using the relation (13) the additional conditions can be rewritten as

[TABLE]

An implicit equation for the coordinates of C*±* on the ($\Omega$,$g_{0}$)-plane are obtained by solving Eq. (14) for $R$ and substituting the resulting solution into Eq. (13). It follows from Eq. (14) that the codimension-two points C*±* shown in the right panel of Fig (3) are symmetric with respect to $\Omega=0$ axis. It is seen from the right panel of Fig. 3 that the central longitudinal mode with $\Omega=0$ is the first one undergoing a flip instability with the increase of the pump parameter $g_{0}$. The larger is the frequency offset of the mode, the higher is the flip bifurcation threshold for this mode. When, however, positive (negative) frequency offset is sufficiently large the mode is already unstable with respect to Turing (modulational) instability at the flip instability point. In this case the flip bifurcation results in the appearance of unstable square waves. Note, that for any finite values of delay time, $\gamma T<\infty$, it is an Andronov-Hopf bifurcation of a CW solution that leads to the emergence of square waves. However, in the limit of infinite delay the period of square wave regime tends to infinity, which means that the imaginary parts of complex eigenvalues crossing the imaginary axis at the Andronov-Hopf bifurcation tend to zero. Hence, in the limit of infinite delay we refer to this bifurcation as a flip instability.

Let us consider the CW solution of Eq. (4) corresponding to the central longitudinal mode with zero detuning from the central frequency of the spectral filter, $\Omega=0$. For this solution the two quantities $Y_{\pm}$ in Eq. (10) take the form

[TABLE]

From Eq. (15) we get the relations $Y_{-}^{-1}\rvert_{\Omega=0,\lambda_{0}=0}=1$ and $\left|Y_{-}^{-1}\right|_{\Omega=0,\lambda_{0}\neq 0}<1$ meaning that the first branch of the pseudo-continuous spectrum of the central longitudinal mode is always stable and is tangent to the imaginary axis at $\lambda_{0}=0$. On the other hand, from Eq. (16) we see that the fold bifurcation of the central longitudinal mode defined by the condition $Y_{+}^{-1}\rvert_{\Omega=0,\lambda_{0}=0}=1$ coincides with Eq. (8). Similarly the flip instability responsible for the emergence of square-waves is defined by the condition $Y_{+}^{-1}\rvert_{\Omega=0,\lambda_{0}=0}=-1$ coinciding with Eq. (13). As it will be shown in the next section, the condition (13) defines also the period doubling bifurcation of a 1D map, which we construct in the next section to study the square wave formation in the DDE model Eq. (4).

IV 1D map and square waves

The existence of stable square wave solutions in Eq. (4) can be demonstrated by constructing a 1D map Chow and Mallet-Paret (1983) that exhibits a period doubling bifurcation corresponding to the emergence of square waves in the DDE model (4). To this end we rescale the time $\tau=t/T$ in Eq. (4) and obtain

[TABLE]

where in the large delay limit we have $\varepsilon\equiv 1/\gamma T\ll 1$. By discarding the time derivative term, which is proportional to the small parameter $\varepsilon$, we transform this equation into a complex map, which describes the transformation of the electric field envelope $E$ after a round trip in the cavity and resembles the well known Ikeda map Ikeda (1979). Then, taking modulus square of both sides we obtain a 1D map for the intensity $R$:

[TABLE]

where $R_{n}\equiv R(n)$ ($n=0,1,2\dots$) with fixed points $R_{n}=R^{*}$ satisfying the condition $R^{*}=h(R^{*})$:

[TABLE]

Graphical representation of the function $h$ is given in the left panel of Fig. 4. Note that since Eq. (19) is equivalent to Eq. (6) taken at $\Omega=0$ and $\Omega=0$ is a solution of Eq. (7), fixed points of the map (18) have the intensity $R$ coinciding with that of the central longitudinal mode, i.e. the CW solution of Eq. (4) with zero frequency offset $\Omega=0$ from the central frequency of the spectral filter. Furthermore, for sufficiently large $g_{0}$ a stable fixed point of the map (18) exhibits a period doubling bifurcation which is defined by

[TABLE]

together with (19). Since the relations (19) and (20) are equivalent to (6) and (12) evaluated at $\Omega=0$, the period doubling bifurcation point of the map (18) coincides in the limit of large delay with the flip instability to square waves of the central longitudinal mode having zero frequency offset, $\Omega=0$. It is seen from right panel of Fig. 4 that after the first period doubling bifurcation the 1D map (18) demonstrates a period doubling transition to chaos.

Period-doubling route to chaos obtained by numerical integration of the DDE model (4) is illustrated in left panel of Fig. 5, where local maxima of the electric field intensity time-trace are plotted versus increasing values of the pump parameter $g_{0}$. It is seen that the diagram in this figure is very similar to that obtained with the 1D map (18), cf. Fig. 4. Note, however, that the period doubling threshold is slightly higher in Fig. 5 than in Fig. 4. This can be explained by taking into consideration that in the DDE model (4) the threshold of the flip instability leading to the emergence of square waves increases with the absolute value of the frequency detuning $\Omega$. Therefore, we can conclude that in Fig. 5 the CW solution undergoing the period-doubling cascade must have a small non-zero frequency detuning, $\Omega\neq 0$. The first period doubling bifurcation in the left panel of Fig. 4 is responsible for the formation of square waves shown in the right panel of Fig. 4, while further period doublings give rise to more complicated square wave patterns with larger periods. Finally, we note that with the increase of the gain parameter $g_{0}$ new pairs of fixed points of the map (18) appear in saddle-node bifurcations. For example, the second pair of fixed points corresponds to the second (right) maximum of the function $h$ shown in the left panel of Fig. 4 and to additional branches of CW solutions visible in the upper right part of the right panel of Fig. 2. Although the linear stability analysis performed in the previous section is applicable to these additional high intensity CW solutions as well, here we restrict our consideration to moderate pumping levels.

V Mode-locking

Bifurcation diagram in the left panel of Fig. 6 was obtained by numerical integraton of Eq. (4). It is similar to that shown in Fig. 5, but spans a larger range of the pump parameter values. It follows from this diagram that with the increase of the pump parameter $g_{0}$ after a chaotic regime associated with the period doubling cascade the phase trajectory of the system jumps to a pulsed solution with time periodic laser intensity. This solution corresponding to a fundamental mode-locked regime with the repetition period close to the cavity round trip time, $T$, is illustrated in right panel of Fig. 6, where it is seen that the leading edge the pulses is much steeper than the trailing edge. In the following, we will show that unlike the trailing tail of the pulses, which decays exponentially, their leading tail demonstrates faster-than-exponential (super-exponential) decay in negative time.

Let us consider a mode-locked solution of Eq. (4) with the period $T_{0}$ close to the cavity round trip time $T$. For this solution satisfying the condition $E(t)=E(t+T_{0})e^{i\Delta}$ we can write

[TABLE]

where $\delta=T_{0}-T>0$ is the small difference between the solution period and the delay time. Substituting (21) into (4) we get a time advance equation

[TABLE]

where $\Theta=\theta+\Delta$. Note, that unlike the original DDE model (4), which has stable trivial solution $E=0$, the trivial solution of Eq. (22) is a saddle with one stable and an infinite number of unstable directions. The stable direction determines the decay rate of the trailing tail of mode-locked pulses, while unstable directions are responsible for the decay of the leading tail in negative time. To perform linear stability analysis of the trivial solution of Eq. (22) we write the following equation linearized at $E=0$:

[TABLE]

where linear time advance term in the right hand side proportional to small perturbation parameter $\epsilon$ describes an imperfection introduced by a slight asymmetry of the coupler between the laser cavity and the NALM. The spectrum of Eq. (23) is defined by

[TABLE]

where $W_{k}$ denotes the $k$th branch of multivalued Lambert function. In particular, in the limit $\epsilon\to 0$ the eigenvalue with the index $k=0$ is the only stable and negative one $\lambda_{0}\to-\gamma<0$. This eigenvalue determines the decay rate of the trailing tail of mode-locked pulses. The remaining eigenvalues with $k\neq 0$ have positive real parts diverging in the limit $\epsilon\to 0$, $\operatorname{Re}\lambda_{k}\to+\infty$. Among these unstable eigenvalues, depending on the value of $\Theta$, one of the two eigenvalues $\lambda_{\pm 1}$ with $k=\pm 1$ has the smallest real part. All other eigenvalues with $|k|>1$ have larger real parts. Hence, for small nonzero $\epsilon$ generically one of the two eigenvalues $\lambda_{\pm 1}$ determines the decay rate the pulse leading tail in negative time. The fact that this eigenvalue tends to infinity as $\epsilon\to 0$ suggests that this decay is super-exponential.

Finally let us discuss briefly the interaction of asymmetric mode-locked pulses shown in the right panel of Fig. 6. When integrating the model equation (4) numerrically it is possible to seed two or more non-equidistant pulses in the laser cavity as an initial condition. Then the pulses start to interact locally via their decaying tails. The asymmetric nature of these tails suggests that similarly to the case discussed in Vladimirov et al. (2018) local interaction of non-equidistant pulses will be very asymmetric as well. This can be is seen in Fig. 7 illustrating an interaction of two asymmetric pulses of Eqs. (4) on the time-round trip number plane. We see that two initially non-equidistant pulses repel each other and tend to be equidistantly spaced in the long time limit. Furthermore, when the two pulses are sufficiently close to one another, exponentially decaying trailing tail of the left pulse repels noticeably the right pulse, while super-exponentially decaying leading tail of the right pulse almost does not affect the position of the left pulse. When the pulses become equidistant the forces acting on a pulse from opposite directions balance each other and a stable harmonic mode-locking regime with two pulses per cavity round trip time is established.

VI Conclusion

We have considered a simple DDE model of unidirectional class-A ring NALM mode-locked laser. Linear stability analysis in the large delay limit revealed that similarly to the well-known Eckhaus instability only those CW solutions, which belong to the Busse balloon can be stable. This balloon is limited from below by the modulational instability boundary. We demonstrated that with the increase of the pump parameter CW regimes loose their stability either via a Turing-type instability, or through a flip instability leading to a formation of stable square waves. We have constructed a 1D map which describes the transition to square waves and their secondary bifurcations giving rise to a more complicated square wave patterns with larger and larger periods. We have shown that mode-locked pulses, which appear after a chaotic square-wave dynamics, have exponentially decaying trailing tail and a leading tail, which decays super-exponentially in negative time. When two or more pulses circulate in the laser cavity, the interaction of these pulses is repulsing and very asymmetric leading to a formation of harmonic mode-locked regimes. Noteworthy, that the mode-locked regime considered here is always non-self-starting, which means that the pulses are sitting on a stable laser-off solution. Hence, these pulses can be viewed as temporal cavity solitons having similar properties to spatial and temporal localized structures of light observed in bistable optical systems.

Authors thank D. Turaev for useful discussions. A.V.K. and E.A.V. acknowledge the support by Government of Russian Federation (Grant 08-08). A.G.V. acknowledges the support of the Fédération Doeblin CNRS and SFB 787 of the DFG, project B5.

VII Appendix

The coefficients $c_{0}(\lambda)$, $c_{1}(\lambda)$, and $c_{2}$ in the characteristic equation (9) are defined by

[TABLE]

The modulational instability condition (11) of the CW solutions of Eq. (4) can be rewritten in the following form:

[TABLE]