TEMPORAL CAVITY SOLITONS IN A DELAYED MODEL OF A DISPERSIVE CAVITY RING LASER

. Nonlinear localised structures appear as solitary states in systems with multistability and hysteresis. In particular, localised structures of light known as temporal cavity solitons were observed recently experimentally in driven Kerr-cavities operating in the anomalous dispersion regime when one of the two bistable spatially homogeneous steady states exhibits a modulational instability. We use a distributed delay system to study theoretically the formation of temporal cavity solitons in an optically injected ring semiconductor-based ﬁber laser, and propose an approach to derive reduced delay-diﬀerential equation models taking into account the dispersion of the intracavity ﬁber delay line. Using these equations we perform the stability and bifurcation analysis of injection-locked continuous wave states and temporal cavity solitons


Introduction
Temporal localised structures of light propagating along the axial direction in nonlinear cavities attracted significant theoretical and experimental attention in the last decade due to their potential applications for optical data storage and transmission [7,12,13,15,16].Typically they appear in the vicinity of parameter regions where a branch of spatially homogeneous steady states exhibits S-shaped hysteresis loop [1,20,22].Similarly to the solitons of nonlinear Schrödinger equation [40], dissipative optical localised structures known also as temporal cavity solitons (TCSs) are localised in time and in longitudinal direction.They can be studied by direct numerical simulations of complex Ginzburg-Landau-type equations [8,11] or alternatively as stationary solutions of properly constructed ordinary differential equations in the co-moving reference frame [30,34,35].Although this approach allows for a detailed bifurcation analysis of TCSs, complex Ginzburg-Landau models are hardly applicable to account accurately for some important physical effects in realistic laser devices, such as those containing intracavity semiconductor medium [19,33].This is why travelling wave-type models [3,37] are commonly used to model the dynamics of semiconductor devices.However, since the traveling wave models are rather complicated and their analysis is usually limited to direct numerical simulations, an alternative and more simple approach to the analysis of multimode semiconductor lasers was proposed in [32,33,36] based on the use of delay differential equations (DDEs).DDE laser models can be derived from the travelling wave equations under certain non-restrictive simplifying physical assumptions and proved to be a viable alternative to the standard models based on partial differential equations.In addition to asymptotic stability analysis [25,33] the DDE approach allows for numerical study of continuous wave (CW) and periodic intensity regimes using well-developed Floquet theory and software packages such as DDE-BIFTOOL [10,17,23,25,32,33,36].
In this paper, using DDE models we investigate the properties of TCSs in an optically injected laser containing a semiconductor optical amplifier and a long dispersive fiber delay line [19,25].Physical mechanisms of the formation of these solitons are similar to those reported earlier in the Lugiato-Lefever equation (LLE) with anomalous dispersion [14] and in driven nonlinear cavities with dissipative nonlinearity and diffraction [6,21,24,26].The LLE is an externally driven, damped nonlinear Schrödinger equation for the electric field in a passive optical cavity subject to weak optical injection slightly detuned from the resonant frequency of the cavity.Although the LLE allows one to perform numerical bifurcation analysis of stationary homogeneous and TCS solutions [13], the simplifying assumptions behind the LLE limit its applicability.For example, the LLE fails to describe the bistability between two TCS branches corresponding to neighbouring longitudinal cavity modes [9], which requires large detunings.On the other hand, this bistability is well captured in the framework of the travelling wave equation approach [11].Similarly to the travelling wave equations DDE models are free from the small detuning approximation.Furthermore, they allow for an adequate account of the characteristic features of intracavity semiconductor optical amplifier and inclusion of the dispersion of the fiber delay line by incorporating the distributed delay term in the model equations [19].Therefore, DDEs provide a more relevant and rigorous tool to study TCS in the physical system considered here than the Ginzburg-Landau type models.In the framework of time delay models, we propose another way to study moving CSs in optical resonators with the help of differential equations involving nonlinear delay terms that complements the existing studies of the effect of feedback loop described by linear delay term in form of Pyragas control [27], which can be applied to the models of broad-area lasers [18,28,31], as well as to nonlinear optical cavities such as fiber resonators or disk microresonators subjected to delayed optical feedback, where a LLE with time delay can be used to study the drift of TCSs [29].
In this study, using the distributed DDE laser model of reference [19] and following the approach described in [39], we perform linear stability analysis of CW states in the limit of large delay.We report the presence of modulational instability (MI) of the upper part of CW solutions branch in the strong anomalous dispersion regime.Above the MI threshold, we demonstrate numerically the formation of stable TCSs.Thereafter, we derive a reduced DDE model with a single delay that preserves the effect of the chromatic dispersion on the dynamics of the ring laser.With the help of the reduced model we perform numerical continuation and stability analysis of the periodic TCS solutions in the reduced model using DDE-BIFTOOL [5].Among other things, we find a narrow region of multistability between TCSs of different width.

Distributed DDE model
Schematic presentation of an optically injected long cavity laser under consideration is given in Figure 1 (left).The laser containing a short semiconductor optical amplifier (SOA) section, linear frequency selective spectral filter, and a long dispersive fiber delay line, is subject to external optical injection from a single-mode laser.A similar device without external injection operating in the so-called Fourier domain mode-locked regime was studied experimentally in [25].A distributed DDE model to describe this device was derived in reference [19].
Here we write the distributed DDE model with an additional term describing single-mode coherent external injection: Here A(t) and B(t) are the electric field envelopes at the entrance and at the output of the fiber delay line, respectively; G(t) is the cumulative saturable gain in the SOA section, and T is the cavity round-trip time.The parameter γ is the spectral filter bandwidth, ω describes the detuning between the central frequency of the filter and the injection frequency, κ < 1 is the round trip attenuation factor due to the linear nonresonant losses in the SOA and output of radiation from the cavity, α is the linewidth enhancement factor, ϕ describes the detuning between injection frequency and the considered cavity mode, γ g is the the carrier relaxation rate in SOA, g 0 is the pump parameter, and σ is proportional to the length of the fiber delay line.The parameter η describes the injection rate.Since it is assumed that the reference frequency coincides with the frequency of the injection, the injection term does not depend on time.For the complex resonance frequency Ω c , see equation (2.4).

Modelling dispersive fiber delay line
The distributed DDE model was derived in reference [19] following the procedure described in [32,33,36] with the help the so-called lumped element approach.In this approach the travelling wave equations are solved for each of the laser sections and the resulting solutions are combined to obtain a closed set of model DDEs.In particular, in [19] it was assumed that the dispersion of the fiber delay line is created by a Lorentzian absorption line with central frequency Ω and full width at half maximum Γ, which is strongly detuned from the laser transition.Neglecting all nonlinearities, propagation of light in the fiber delay line was described by the set of two differential equations for the electric field envelope E(t, z) and polarisation associated with the absorption line P L (t, z): where z is the coordinate along the fiber line, t = t − z/v gr relates the retarded time t to the laboratory t , and v gr is the group velocity at the carrier frequency.The parameter ζ is proportional to the oscillator strength of the Lorentzian absorption line, the latter is associated with the complex frequency Ω c .The distributed delay term P (t − T ), which appears in the system (2.1)-(2.2) through equation (2.3), was derived in [19] by integration of equations (2.4).An important drawback of the distributed DDE model (2.1)-(2.2) is, however, that its numerical simulations are time consuming for small (positive) Γ due to the slow decay of the Bessel function.Below we re-examine the derivation of the distributed delay term and describe an approximate approach, which allows us to replace this term by a single Lorentz-type ODE for P (t).This will be done by invoking the Laplace transformation to solve equation (2.4) and by using the Padé approximants when performing the inverse Laplace transformation.

Derivation of the distributed delay term
The Laplace transformed system (2.4), assuming for simplicity that the polarisation relaxes between the pulses and therefore taking P (z, t)| t≤0 = 0, reads Here the Laplace image f (ω) of a causal function f (t) vanishing for t < 0 is defined by  where ω ∈ C is the complex frequency having "sufficiently positive" imaginary part.The parameter C ∈ R is large enough so that all singularities of f (ω) are below the integration contour R + iC in the complex plane.Since typical f (ω), being analytically continued, has poles in the lower half-plane (on the real axes at most), it is sufficient to set C = +0.The integration contour can then be pushed down as in Figure 2a.Integrating equation (2.5) along the longitudinal coordinate z and performing inverse Laplace transform we obtain where z 1 is the longitudinal coordinate at the entrance of the fiber delay line.Taking z = z 2 , where z 2 is the delay line output coordinate, we write the relation between input and output fields, E(z 1 , t) and E(z 2 , t): where σ = ζ(z 2 − z 1 ).Introducing the notations E(z 1 , t) ≡ A(t) and E(z 2 , t) ≡ B(t) we can rewrite equation (2.8) in the form: with Green's function G(t) vanishes for t < 0 and is defined by for t > 0. To calculate the Bromwich integral in equation (2.10) we use the contour shown in Figure 2(a) and expand the first exponential term.Equations (2.9) and (2.10) yield equation (2.3) [19].In the next sections we propose two approximations of Green's function G(t) leading to DDE laser models with a single delay T .

Approximation of Green's function I
When the dispersion of the fiber delay line is sufficiently weak, the main contribution to the integral in equation (2.8) comes from the frequency domain satisfying the inequality |ω − Ω c | σ.However, a naively truncated Taylor expansion of the first exponential in equation (2.10) yields unphysical behaviour, e.g., polynomial divergence for Γ = 0 and t → ∞.A more relevant approach is to use Padé approximants [2].The simplest with a shifted pole at ω = Ω − i(Γ + 1 2 σ).Since σ, is assumed to be positive, the pole still is in the lower half-plane of the complex plane.The resulting Green's function G(t) is then causal and vanishes for t → ∞.
Substituting the [1/1] Padé approximant (2.11) into equation (2.10) and integrating along the contour shown in Figure 2a we obtain which gives the following ODE for the quantity P (t) defined by (2.9): (2.12) Equation (2.12) can be used instead of the second equation in (2.3) to close the system (2.1)-(2.2).The above scheme can be easily generalised by using higher order Padé approximants to construct causal approximations of Green's function G(t) as illustrated in Figure 2b.
Since the condition |ω − Ω c | σ is quite restrictive for adequate description of the laser dynamics as discussed in Sections 3.2 and 4.2, in the next subsection we propose an alternative and more accurate model to approximate equation (2.3).This model also contains the distributed delay term and also leads to a single ODE for the quantity P .

Approximation of Green's function II
Here we assume that the laser field has a narrow spectrum localised near the injection frequency ω 0 = 0 and apply the Padé approximation around this frequency.Rewriting equation (2.7) in the form replacing the first exponent under the integral with its [1/1] Padé approximant, and taking z = z 2 we obtain: where with new Green's function ) where Note that the factor e 2iΘ describing an additional phase shift and losses induced by the fiber delay line can be easily eliminated by including it in the parameter ϕ and κ in equation (2.16).Furthermore, similarly to equation (2.12), equation (2.18) suggests that the variable P (t) exhibits faster decay than e −Γt .Such decay existing also for Γ = 0 is somewhat similar to the Landau damping.

Stability analysis in the limit of large delay
Steady state solutions of equations (2.1)-(2.3)correspond to the injection-locked CW laser regimes with the rotation frequency equal to that of the injection field.Since the reference frequency is chosen to coincide with the injection frequency they can be written in the form A = A 0 e iϕ0 and G = G 0 with time independent A 0 , G 0 , and ϕ 0 .Using these relations we get from equation (2.3) P = P 0 = e −σ/[Γ+iΩ] − 1 A 0 .Then, substituting the expressions for A, G, and P into equations (2.1)-( 2.3) we come to a system of transcendental equations for A 0 , ϕ 0 , and G 0 , which can be solved numerically.For example, an S-shaped hysteresis loop formed by three injection-locked CW states calculated in the absence of chromatic dispersion, σ = 0, is shown on Figure 1 (right).
Since our aim is to study formation of TCSs in the long laser with round-trip time T 100, we can perform linear stability analysis of steady state solutions of equations (2.1)-(2.3) in the limit of large delay using the method described in [39].Although originally this method was developed for DDEs with a single discrete delay, it was demonstrated in [19,38] that it is easily extended for the analysis of the distributed delay laser model.We linearise the system (2.1)-(2.3)near the steady states A = A 0 + δAe λτ , G = G 0 + δGe λτ , and P = P 0 + δP e λτ with the relation δP = δA e −σ/[Γ+λ+iΩ] − 1 obtained from (2.3).Substituting these expressions into (2.1) and (2.2) we obtain the following characteristic equation for the eigenvalues λ describing the stability of the steady state solution: where Y = e −λT is the exponential term that comes from the delayed variables A(t − T ) and P (t − T ) and with .
In the limit of large delay time T → ∞ the eigenvalues belonging to the pseudo-continuous spectrum can be represented in the form λ = iµ + Λ T + O(1/T 2 ) with real µ [39].Keeping only the single leading term iµ in a 0 (λ), a 1 (λ), a 2 (λ) and two leading terms iµ + Λ T in Y (λ), we obtain from (3.1) two branches of pseudo-continuous spectrum given by Real parts of these eigenvalue branches are shown in Figure 3 for η = 0.006 and σ = 0. Left, central, and right panels in this figure belong, respectively, to lower, middle, and upper parts of S-shaped CW branch shown in Figure 1 (right).We see that for the chosen parameter values the CW solutions belonging to the upper and lower parts of CW branch are stable.The middle part of the branch is always unstable.Finally, we note that the necessary condition of the MI of CW solutions of a distributed DDE model (2.1)-(2.3)without optical injection, η = 0, was obtained in [19].In the next subsection we compare this condition with those obtained with the simplified models (2.1)-(2.2),(2.12) and (2.16)-(2.18).In Section 4.1, we use these conditions to locate MI of the injection-locked CW state and find a stable temporal cavity soliton.

Modulational instability
In this section we compare analytical conditions for the appearance of long wavelength MI of CW solutions of the distributed DDE model with those obtained using the reduced DDE models.For simplicity we consider the case when the external injection is absent, η = 0.In this case the system (2.1)-(2.3)has a phase shift symmetry so that for any given solution A(t), P (t), and G(t) the phase shifted solution A(t)e iφ , P (t)e iφ , and G(t) with arbitrary constant φ also a satisfies the system.Due to this symmetry CW solutions of this system take the form A(t) = A 0 e iνt , P (t) = P 0 e iνt , and G(t) = G 0 and the characteristic equation (3.1) remains invariant under the transformations w → w − ν and Ω → Ω + ν.In addition, for ν = 0 the amplitudes of the CW solutions |A 0 |, |P 0 |, and G 0 can be obtained explicitly [19].
Owing to the phase shift symmetry of equations (2.1)-(2.3)with η = 0 one of the two eigenvalue branches Λ ± (µ) satisfies the condition Λ(0) = 0. Long wavelength MI takes place when the second derivative of this branch d 2 Re Λ(µ)/dµ 2 changes its sign at the point µ = 0 from negative to positive, see Figure 4. Using equation (3.2) and analytical expressions for the coefficients a 0 (λ), a 1 (λ), and a 2 (λ) entering the characteristic equation (3.1), a necessary condition for appearance of MI of CW solution with the frequency ν in a laser with anomalously dispersive fiber delay line and η = 0 can be obtained [19]: Here the second-order dispersion coefficient is given by Let us now for simplicity consider the particular limit when Γ = 0.In this limit the necessary and sufficient MI condition for CW solutions of the distributed DDE model (2.1)-(2.3)takes a relatively simple form [19]: where D 2 is obtained by substitution Γ = 0 into (3.4) and F represents auxiliary terms responsible for appearance of MI of the side modes which exists even in the normal dispersion regime: .
Note that unlike (3.3) the condition (3.5) accounts for the MI of the side modes that can take place even in the normal dispersion regime.Since in this study we are mainly interested in the formation of TCSs in the anomalous dispersion regime, we consider the maximum gain mode, which cannot undergo MI in the normal dispersion regime.For this mode we can set w = ν = 0. Then for Γ = 0 we have D 2 = 2σ/Ω 3 and the MI condition (3.5) takes the simplified form Next, we consider the MI of the CW solutions of the reduced DDE models (2.1)-(2.2),(2.12) and (2.16)-(2.18).Since in this case the coefficients of the characteristic equation (3.1) become very cumbersome, we do not provide an explicit formula for Y (µ).However, in the particular case Γ = 0, when only Landau damping remains in equations (2.12) and (2.18), we can obtain rather simple analytical conditions of the MI.
For the reduced DDE model (2.1)-(2.2),(2.12) we get the following MI condition For the maximum gain mode with w = ν = 0 instead of (3.7) we obtain the condition which differs from (3.6) by an additional term σ 2 4Ω 2 .Due to the presence of this term the condition (3.8) holds only for moderate values of σ, while for large σ the absolute value of the left-hand side of this condition is always small and tends to zero as σ → ∞.In particular, from (3.8) we can obtain a necessary condition Ω ≥ −γ 3 √ 3α/2 for the MI in the anomalous dispersion regime (Ω < 0).On the contrary, the condition (3.3) is always satisfied in the anomalous dispersion regime when σ is sufficiently large.
Finally, for the second reduced DDE model (2.16)-(2.18)with Γ = 0 we get the MI condition, which coincides exactly with (3.5) obtained with the distributed DDE model.It is worth noting that for Γ > 0 the exact agreement with the distributed DDE model is lost, however, the reduced DDE model is regular enough for Γ = 0 due to Landau damping in all our simulations.Therefore, we can conclude that the second reduced DDE model (2.16)-(2.18)reproduces better the MI of the original distributed DDE model.

Temporal cavity solitons
In this section, we perform numerical analysis of the original distributed DDE model (2.1)-(2.3)and the reduced DDE models (2.1)-(2.2),(2.12) and (2.16)-(2.18).For each of these three models we locate MI of the upper part of S-shaped injection-locked CW branch and find numerically periodic TCSs in the range of parameter values that are in the vicinity of this instability.Furthermore for reduced DDE models we perform numerical bifurcation analysis of TCSs using the software package DDE-BIFTOOL [5].Note, that in our numerical examples given below we choose w = 0 and ϕ = ϕ c + δϕ in equations (2.1)-(2.3),where is a phase shift that places the maximum gain mode of the laser at the reference (zero) frequency, which coincides with the frequency of the injected field.In this case synchronisation of the laser output to the external injection takes place at zero frequency and the bistability can be observed near the boundaries of the locking range.

Distributed DDE model
In reference [13] a formation of TCSs was reported in the anomalous dispersion regime, when the upper part of the S-shaped CW branch exhibits a MI.Since unlike the experimental setup of this paper, semiconductor based fiber laser can hardly be described by simple mean-field models, like the LLE, in this section we use the distributed DDE model to find the parameter range where TCSs can exist.As it was demonstrated analytically in reference [19] (see also Sect.3.2) in the absence of external injection, η = 0, CW solutions of (2.1)-(2.3)become modulationally unstable when anomalously dispersive (Ω < 0) fiber delay line is sufficiently long, i.e.,  1.Other panels: Real parts of two branches of the pseudo-continuous spectrum ReΛ ± (µ) for the bottom fold point at η ≈ 0.0051 (top-center), upper fold point at η ≈ 0.006 (top-right), and for three steady states observed at η ≈ 0.0058, where short wavelength MI appears on the upper part of the curve (bottom-right), the steady state from the middle part is unstable (bottom-center), whereas the steady state from the bottom part is stable (bottom-left). the parameter σ is sufficiently large.When the injection rate η becomes nonzero the MI remains in the system but requires even larger values of σ.The pseudo-continuous spectrum of the upper branch steady state is shown in Figure 5, where top-right and bottom-right panels illustrate long-and short-wavelength MI, respectively.Note that for the parameter values of this figure threshold in σ is approximately 3 times higher than given by (3.6) for η = 0.
Stable TCS solution of the distributed DDE laser model (2.1)-( 2.3) with the repetition period close to the cavity round trip time is shown in Figure 6.This solution was obtained numerically for small injection amplitude, η = 0.0058, by seeding a sufficiently large localised perturbation of the CW solution, which belongs to the lower part of the S-shaped bistability curve shown in Figure 5.A large spike of the field amplitude |A(t)| in Figure 6 is accompanied by a weak drop of the gain variable G(t).Unlike Figures 5 and 6, the next two figures, 7 and 8 correspond to the case of strong optical injection η ≈ 0.4, when the locking frequency is quite far from the center of the locking cone ϕ = ϕ c − 1.65.In this case the numerically calculated TCSs correspond to much narrower and higher spikes of the electric field than in the case of weak injection, and the gain depletion is much stronger as well, see Figure 8.It follows from our numerical simulations that the larger is the injection rate the higher is the dispersion necessary to achieve MI and the formation of TCSs.Therefore, experimental observation of the TCSs in a strongly injected laser would require the use of sufficiently long fiber delay line.
In the next two subsections we perform bifurcation analysis of the TCS solutions with the help of the reduced DDE models (2.1)-(2.2),(2.12) and (2.16)-(2.18).3) for σ = 6000, ϕ = ϕ c − 1.65, α = 4.5, g 0 = 1.38, and other parameters are as in Figure 5.Other panels: Real parts of two branches of the pseudo-continuous spectrum ReΛ ± (µ) for the steady states observed at η ≈ 0.43, where MI appears on the upper part of the curve (right), and the steady state from the bottom part is stable (center).

Reduced DDE model I
In this subsection, we briefly describe the results obtained using the simplest reduced DDE model (2.1)-(2.2),(2.12).In particular, we show that this model has certain deficiencies that could be improved by using the second DDE model (2.16)-(2.18),which is studied in the next subsection.Since some of the analysis of this subsection does not appear to be in full qualitative agreement with that of the distributed DDE model, we leave a detailed discussion of the results to the next subsection, where similar results are demonstrated to be more trustworthy.
According to the condition (3.8) in order to achieve MI we have to choose the parameter Ω sufficiently close to zero.For that, we have used Ω = −0.3 that was first tested in the end of the previous subsection for the distributed DDE model (2.1)-(2.3)as shown on Figure 9, and successfully located a MI of the injection-locked CW state for σ = 0.061.Using DDE-BIFTOOL [5], we can observe the same S-shaped branch of injection-locked CW states on Figure 10, where at η ≈ 0.009 the steady state from the bottom part of the curve is stable, the steady state from the upper part is modulationally unstable, and the curves of the pseudo-continuous spectrum are very similar to the ones observed in the full model.Nevertheless, one can immediately spot the difference  5. Other panels: Real parts of two branches of the pseudo-continuous spectrum ReΛ ± (µ) for the steady states observed at η ≈ 0.0058, where MI appears on the upper part of the curve (right), and the steady state from the bottom part is stable (center).We see that comparing to Figure 5, where injection strength is similar, the losses induced by Lorentzian absorption line deform the spectrum around µ = ±0.3 at the point of MI. in the stability properties of the steady states in this DDE model: the upper part of the CW branch is always modulationally unstable, whereas there is a Hopf bifurcation on the bottom part, which is not observed in the full model.Similarly to LLE, the branch of TCS solutions is born in the fold bifurcation point of the bottom part of the CW branch, and we observe a little bit different bifurcation scenario for the stable TCSs than in LLE, which we will discuss in the next subsection in more detail.In particular, we can observe on Figure 11 multistable TCSs of different widths for η = 0.009, however we note that the shape of TCS is strongly distorted in comparison to the shape of TCSs observed in the distributed DDE model, see Figure 6.

Reduced DDE model II
In this subsection we perform numerical investigation of the reduced DDE model (2.16)-(2.18)that, according to Section 3.2, should be valid for larger σ and |Ω| than in the case of equations (2.1)-(2.2),(2.12),discussed in the previous subsection.Another advantage of this model is that we could naturally eliminate linear losses and phase shift induced by the Lorentzian absorption line and simplify our analysis, because the increase of 2), (2.12) for σ = 0.061, Γ = 0.0195 obtained using DDE-BIFTOOL [5].Other parameters are as in Figure 9.
Here, the simple S-shaped curve is formed by injection-locked CW states, whereas the branch of the TCS solutions starts at the rightmost fold bifurcation point of the S-shaped curve, and a spiralling behaviour is observed, which leads to the appearance of multiple stable branches of TCSs of various types, which are depicted on Figure 11.We discuss this behaviour in more details in the analysis demonstrated on Figure 12 of a more adequate reduced DDE model.Other panels: Spectra of the steady states for η ≈ 0.009 for the bottom part (center) and the upper part (left) of the S-shaped curve.Next, we can obtain periodic TCSs as the direct numerical solution of (2.16)-(2.18)exactly as in Section 4.1 by choosing as an initial condition the stable background steady state from the bottom part of the CW branch perturbed with a narrow Gaussian pulse e −(t+T /2) 2 /100 for t ∈ [−T, 0].Alternatively, one can perturb optical injection η with similar pulse, which is more realistic from the experimental point of view.On the other hand, optical injection is usually performed using a single-mode CW laser, which does not normally emit short pulses, hence we also explore a possibility to excite the TCSs by controlling the parameters of the system similarly to  1. Right: Bifurcation diagram obtained using DDE-BIFTOOL [5].Here, the simple S-shaped curve is formed by injection-locked CW states, whereas the branch of the TCS solutions starts at the rightmost fold bifurcation point of the S-shaped curve, and a spiralling behaviour is observed, where several stable parts of the TCSs branch shown on Figure 13 are connected via fold bifurcations and unstable parts, and the width of the TCS increases after each fold bifurcation.[15].In that paper, several branches of mode-locked periodic solutions with different number of pulses on a period were observed to be multistable with a zero |A 0 | = 0 (laser off) background state below the lasing threshold in a two-section passively mode-locked ring semiconductor laser, where each branch appeared and disappeared for different values of injection current, and by increasing and decreasing injection current in the gain section g 0 above and bellow the lasing threshold one can easily switch between the branches and obtain different number of localised pulses on a period.Since in our system for η = 0.006 we observe a bistability between a bottom steady state and a non-localised periodic solution that is born from the Hopf bifurcation of the upper steady state under MI (see Fig. 12), trajectory of (2.16)-(2.18)usually converges to one of these solutions.We propose the following technique to excite stable TCSs for the injection strength parameter η = η b inside the region of the S-shaped bistability.First, we choose a constant initial function A(t) = A in , G(t) = G in , P (t) = P in for t ∈ [−T, 0] and solve equations (2.16)-(2.18)for t ∈ [0, t 1 ] with sufficiently strong injection η = η 1 > η b , for which there is only one stable steady state that belongs to the upper part of the CW branch.We choose t 1 so that |A(t 1 )| is sufficiently large, then reduce injection strength to η = η 1 < η b , for which there is also only one stable steady state that belongs to the bottom part of the CW branch, and solve the system for t ∈ [t 1 , t 2 ].Here, t 2 is chosen so that |A(t 2 )| is sufficiently small, and finally we set η = η b and solve the system for t ∈ [t 2 , t 3 ], and at the final moment of time the trajectory is close enough to a desired attractor.This technique turns out to be quite sensitive to the choice of the parameters, however we could successfully excite single TCS per round trip by taking η b = 0.006, η 1 = 0.007, η 2 = 0.001, t 1 = 40000, t 2 = 47600 for the parameters used in Figure 12.
Some results of bifurcation analysis of (2.16)-(2.18)obtained using DDE-BIFTOOL are presented in Figure 12 (right).Unlike the bifurcation diagram of Figure 10 obtained for the reduced DDE model (2.1)-(2.2),(2.12), here we observe a better qualitative agreement with the diagram of the distributed DDE model shown in Figure 5 (top-left).In particular, the bottom part of the S-shaped curve formed by injection-locked CW state is always stable in the considered parameter range, whereas the upper part of the curve exhibits a MI at η ≈ 0.006.Nevertheless, the the periodic TCS solutions looks similar to that shown in Figure 10.This branch is born at the saddle-node bifurcation of the bottom part of the branch of the injection-locked CW states at η ≈ 0.0061, stabilised via another fold bifurcation at η ≈ 0.006, and looses stability via a fold bifurcation at η ≈ 0.00607.Further evolution of TCSs could be seen in a magnified portion of the bifurcation diagram on Figure 13 (left).Here, another TCS gains stability at η ≈ 0.00605, which is a wider TCS, displayed on Figure 13 (right), then it looses stability at η ≈ 0.006056, and we observe a spiralling behaviour with multiple TCSs of increasing widths at η ≈ 0.006056.Spiralling behaviour is not observed in the LLE model [13], however it bares similarities to the bifurcation scenario for the spatial cavity soliton branches demonstrated earlier in the models of broad-area semiconductor devices [4].We note that the shape of the main TCS pulse, shown on Figure 13 (right), looks regular and much closer to the shape of the pulses obtained in the full model (2.1)-( 2.3) and demonstrated on Figures 6, 8 than to the pulses observed in the simplified model (2.1)-(2.2),(2.12) that are discussed in the previous subsection and depicted on Figure 11.

Conclusion
In this paper, we have proposed for the first time an efficient methodological approach for numerical bifurcation analysis of periodic TCSs in dispersive time-delay systems that was previously available only for envelope PDEs.We have studied a model of an optically injected ring cavity laser containing a semiconductor gain medium and long fiber delay line.The most basic and general prerequisite for the existence of TCSs in this model is the presence of anomalous dispersion leading to a MI of the upper branch of the hysteresis loop formed by CW solutions locked to external injection frequency.It was demonstrated recently [19] that the effect of chromatic dispersion on the laser dynamics can be accurately accounted using a distributed DDE model (2.1)-(2.3).Here we have revisited the derivation of the this model and developed with the help of Padé approximants two reduced DDE models, which preserve the causality property of the original model.We have shown that although both the reduced models can exhibit a MI in the anomalous dispersion regime, only the second reduced model (2.16)-(2.18)reproduces well the MI threshold of the original distributed DDE model when the dispersion is sufficiently strong.We have demonstrated numerically that all the three models, the distributed DDE one and the two reduced models, can support TCSs in a certain parameter range and used the software package DDE-BIFTOOL to continue TCS branches in the parameter space.In particular, we have spotted a spiralling behaviour of the TCS branch similar to that reported earlier in spatially-distributed models with transverse diffraction and dissipative nonlinearity [4,6,21].Moreover, we have discussed possible techniques for excitation and control of TCSs in realistic experimental setups and the role of hysteresis.We have observed that the slow increase and decrease of the control parameter value leads to the attraction of the laser output trajectory to one of the two bistable non-localised periodic patterns, and have demonstrated how to switch the operation of the laser to the periodic localised TCS regime by a sufficiently fast and precise change of the control parameter that takes into account the properties of the hysteresis loop.Finally, the results of our analysis indicate that the dynamics of the distributed DDE model proposed in [19] can be well captured by the reduced DDE model (2.16)- (2.18).Therefore, we expect that being much simpler than the distributed DDE laser model these equations can become useful tool for studying the effect of chromatic dispersion on the dynamics of different types of optoelectronic devices, and, in particular, mode-locked semiconductor lasers.

Figure 2 .
Figure 2. (a) A transformed integration contour in the complex ω plane for the Bromwich integral (2.10).(b) Comparison of the Bessel function in the kernel of equation (2.3) to the expressions derived from equation (2.10) using diagonal Padé approximants.

Figure 3 .
Figure 3. Real parts of two branches of pseudo-continuous spectrum ReΛ ± (µ) for the injectionlocked CW states of (2.1)-(2.2) lying on the bottom (left), middle (center), and the upper (right) part of the S-shaped curve shown in the right panel of Figure 1.η = 0.006, other parameters are as in Figure 1.

Figure 5 .
Figure 5. Top-left: S-shaped curve, formed by injection-locked CW states of (2.1)-(2.3)for Ω = −13, Γ = 0.001, σ = 2000, ϕ = ϕ c − 0.2, and other parameters are as in Figure1.Other panels: Real parts of two branches of the pseudo-continuous spectrum ReΛ ± (µ) for the bottom fold point at η ≈ 0.0051 (top-center), upper fold point at η ≈ 0.006 (top-right), and for three steady states observed at η ≈ 0.0058, where short wavelength MI appears on the upper part of the curve (bottom-right), the steady state from the middle part is unstable (bottom-center), whereas the steady state from the bottom part is stable (bottom-left).

Figure 6 .
Figure 6.Stabilised transient of the TCS solution for the field amplitude |A(t)| (left) and gain G(t) (center), obtained using numerical solution of (2.1)-(2.3)for T = 400, η = 0.0058 and other parameters as in Figure 5. Long-time transient of the field amplitude |A(t)| in a spacetime representation (right), where the time trace is divided into intervals of length T , and for each subsequent interval a horizontal line representing the color code of the field amplitude is plotted at the vertical position corresponding to the number of round-trips passed before this interval.

Figure 7 .
Figure 7. Left: S-shaped curve (left) formed by injection-locked CW states of (2.1)-(2.3)for σ = 6000, ϕ = ϕ c − 1.65, α = 4.5, g 0 = 1.38, and other parameters are as in Figure5.Other panels: Real parts of two branches of the pseudo-continuous spectrum ReΛ ± (µ) for the steady states observed at η ≈ 0.43, where MI appears on the upper part of the curve (right), and the steady state from the bottom part is stable (center).

Figure 8 .Figure 9 .
Figure 8. Stabilised transient of the TCS solution (left) for the field amplitude |A(t)| (solid ) and gain G(t) (dashed ), obtained using numerical solution of (2.1)-(2.3)for T = 400, η = 0.4 and other parameters as in Figure 7. Long-time transient of the field amplitude |A|(t) in a space-time representation (right), see Figure 6 for details.

Figure 10 .
Figure 10.Left: Bifurcation diagram of the reduced DDE model (2.1)-(2.2),(2.12) for σ = 0.061, Γ = 0.0195 obtained using DDE-BIFTOOL[5].Other parameters are as in Figure9.Here, the simple S-shaped curve is formed by injection-locked CW states, whereas the branch of the TCS solutions starts at the rightmost fold bifurcation point of the S-shaped curve, and a spiralling behaviour is observed, which leads to the appearance of multiple stable branches of TCSs of various types, which are depicted on Figure11.We discuss this behaviour in more details in the analysis demonstrated on Figure12of a more adequate reduced DDE model.Other panels: Spectra of the steady states for η ≈ 0.009 for the bottom part (center) and the upper part (left) of the S-shaped curve.

Figure 12 .
Figure 12.Left: Real parts of two branches of the pseudo-continuous spectrum ReΛ ± (µ) for the injection-locked CW state exhibiting MI in the reduced DDE model (2.16)-(2.18).Here, Ω = −2, σ = 9, g 0 = 1.33, κ = 0.25, η = 0.006, Γ = 0, ϕ = ϕ c − 0.252 and other parameters are as in Figure1.Right: Bifurcation diagram obtained using DDE-BIFTOOL[5].Here, the simple S-shaped curve is formed by injection-locked CW states, whereas the branch of the TCS solutions starts at the rightmost fold bifurcation point of the S-shaped curve, and a spiralling behaviour is observed, where several stable parts of the TCSs branch shown on Figure13are connected via fold bifurcations and unstable parts, and the width of the TCS increases after each fold bifurcation.

Figure 13 .
Figure 13.Left: Magnified portion of the bifurcation diagram from Figure 12 for η ∈ [0.006045, 0.00607].Right: Profiles of various types of TCSs for the field amplitude |A(t)| observed in the bifurcation diagram for η ≈ 0.006056, where the TCS gets wider while we continue the trajectory along the spiral.