Interaction of solitons and the formation of bound states in the generalized Lugiato-Lefever equation

TL;DR

This paper investigates how solitons in the generalized Lugiato-Lefever equation form bound states through their oscillatory tails, analyzing stability, effects of higher order dispersion, and noise, revealing dynamics of soliton interactions.

Contribution

It introduces a detailed analysis of soliton bound state formation and stability in the generalized Lugiato-Lefever equation, including effects of higher order dispersion and noise.

Findings

Bound states form via oscillatory tail interactions.

Higher order dispersion and noise influence soliton dynamics.

Reversible perturbations lead to soliton separation, while broken reversibility can stabilize or push solitons together.

Abstract

Bound states, also called soliton molecules, can form as a result of the interaction between individual solitons. This interaction is mediated through the tails of each soliton that overlap with one another. When such soliton tails have spatial oscillations, locking or pinning between two solitons can occur at fixed distances related with the wavelength of these oscillations, thus forming a bound state. In this work, we study the formation and stability of various types of bound states in the Lugiato-Lefever equation by computing their interaction potential and by analyzing the properties of the oscillatory tails. Moreover, we study the effect of higher order dispersion and noise in the pump intensity on the dynamics of bound states. In doing so, we reveal that perturbations to the Lugiato-Lefever equation that maintain reversibility, such as fourth order dispersion, lead to bound states that tend to separate from one another in time when noise is added. This separation force is determined by the shape of the envelope of the interaction potential, as well as an additional Brownian ratchet effect. In systems with broken reversibility, such as third order dispersion, this ratchet effect continues to push solitons within a bound state apart. However, the force generated by the envelope of the potential is now such that it pushes the solitons towards each other, leading to a null net drift of the solitons.

Tables

Table 1. Table 1: Several stable separation distances z n s subscript superscript 𝑧 𝑠 𝑛 z^{s}_{n} calculated using the full Hamiltonian ( U H subscript 𝑈 𝐻 U_{H} ), the interaction Hamiltonian ( U H I subscript 𝑈 subscript 𝐻 𝐼 U_{H_{I}} ), and numerical time-evolution simulations. Parameters as in Fig.  3 .

Pinning dist.	$U_H$	$U_{H_I}$	Exact
$z_1^s$	35.2979	35.6575	36.3181
$z_2^s$	44.3926	44.7891	45.4856
$z_3^s$	53.5116	53.931	54.6523
$z_5^s$	80.0	80.0	80.0
$z_{n+1}^s-z_n^s$	9.1167	9.1416	9.1668
$\left(z_{n+1}^s-z_n^s\right)/2\pi$	0.6892	0.6873	0.6854

Table 2. Table 2: Stable separation distances z n s subscript superscript 𝑧 𝑠 𝑛 z^{s}_{n} for θ = 2.2 𝜃 2.2 \theta=2.2 and ρ = 1.5 𝜌 1.5 \rho=1.5 calculated with the interaction potential are compared with the exact numerical solution, both for a situation with FOD ( d 4 = 0.15 subscript 𝑑 4 0.15 d_{4}=0.15 ) and a situation with TOD ( d 3 = 0.3 subscript 𝑑 3 0.3 d_{3}=0.3 ).

$d_4=0.15$	$z_4^s$	$z_9^s$	$z_{16}^s$	$z_{25}^s$
$U_{H_I}$	12.7976	24.2656	40.3503	61.0321
Exact	12.6436	24.2517	40.4122	61.119
$d_3=0.3$	$z_0^s$	$z_3^s$	$z_7^s$	$z_{13}^s$
$U_{H_I}$	14.788130	19.357768	25.453909	34.600589
Exact	14.966661	19.517315	25.601455	34.742282

Equations

$${\partial }_{t}A=-(1+i\theta )A+i{\partial }_x^{2}A+iA|A{|}^2+\rho$$

$$i\frac{d^{2}A}{d{x}^2}(x)-(1+i\theta )A(x)+iA(x)|A(x){|}^2+\rho =0,$$

$$I_0^3-2\theta I_0^2+(1+{\theta }^2)I_0={\rho }^2$$

$$I_{b,t}=\frac{2\theta }{3}\pm \frac{1}{3}\sqrt{{\theta }^2-3}.$$

$$\left[\begin{array}[]{c}{\rm Re}[A_{0}]\\ {\rm Im}[A_{0}]\end{array}\right]=\left[\begin{array}[]{c}\displaystyle\frac{\rho}{1+(I_{0}-\theta)^{2}}\\ \displaystyle\frac{(I_{0}-\theta)\rho}{1+(I_{0}-\theta)^{2}}\end{array}\right].$$

$$\mathcal{S}[A,\bar{A}]\equiv {\int }_{\mathbb{R}}{e}^{2t}L[A,\bar{A}]dt={\int }_{\mathbb{R}}{e}^{2t}{\int }_{\mathbb{R}}\mathcal{L}[A,\bar{A}]dxdt,$$

$$\begin{gathered} \mathcal{L}[A,\bar{A}]\equiv \frac{i}{2}\left(\bar{A}{\partial }_{t}A-A{\partial }_t\bar{A}\right) \\ -{\partial }_{x}A{\partial }_x\bar{A}-i\rho \left(\bar{A}-A\right)-\theta A\bar{A}+\frac{1}{2}{\bar{A}}^2{A}^2. \end{gathered}$$

$$\frac{\delta \mathcal{S}[A,\bar{A}]}{\delta \bar{A}}=0$$

$$\mathcal{H}[A,\bar{A}]={\partial }_{x}A{\partial }_x\bar{A}+i\rho \left(\bar{A}-A\right)+\theta A\bar{A}-\frac{1}{2}{\bar{A}}^2{A}^2,$$

$${\mathcal{H}}_{\mathrm{I}}[A,\bar{A}]=-\frac{1}{2}{A}^2{\bar{A}}^2.$$

$$H[A,\bar{A}]={\int }_{\mathbb{R}}\mathcal{H}[A,\bar{A}]dx.$$

$$A(x,z)=A_{-}(x)+A_{+}(x)-A_0,$$

$${\delta }_z\mathcal{S}(z)=\frac{d}{dz}{\int }_{\mathbb{R}}{e}^{2t}H(z)dt=0.$$

$$\frac{dH}{dz}=0.$$

$$U_H(z)=H[A(z)]={\int }_{\mathbb{R}}\mathcal{H}[A(z)]dx.$$

$$H_I=-\frac{1}{2}{\int }_{\mathbb{R}}{A}^2{\bar{A}}^{2}dx=-\frac{1}{2}{\int }_{\mathbb{R}}|A{|}^{4}dx,$$

$$U_{H_I}(z)=H_I[A(z)]=-\frac{1}{2}{\int }_{\mathbb{R}}|A(x,z){|}^{4}dx.$$

$$\begin{array}[]{l}d_{x}y_{1}=y_{3}\\ d_{x}y_{2}=y_{4}\\ d_{x}y_{3}=y_{2}+\theta y_{1}-y_{1}y_{2}^{2}-y_{1}^{3}\\ d_{x}y_{4}=-y_{1}+\theta y_{2}-y_{2}y_{1}^{2}-y_{2}^{3}+\rho.\end{array}$$

$${\lambda }^4+c_2{\lambda }^2+c_0=0.$$

$$\lambda =\pm \sqrt{(\theta -2I_0)\pm \sqrt{I_0^2-1}}.$$

$$q_0=\frac{1}{\sqrt{2}}\sqrt{\sqrt{c_0}+\theta -2I_0}$$

$$k_0=\frac{1}{\sqrt{2}}\sqrt{\sqrt{c_0}-(\theta -2I_0)}.$$

$$\mathrm{Re}[A]=\mathrm{Re}[A_0]+a_1{\mathrm{e}}^{q_{0}x}\cos (k_{0}x+{\phi }_1),$$

$${\lambda }_0^{tails}=\frac{2\sqrt{2}\pi }{\sqrt{\sqrt{c_0}-(\theta -2I_0)}}.$$

$${\partial }_{t}A=-(1+i\theta )A+i{\partial }_x^{2}A+d_3{\partial }_x^{3}A+i{d}_4{\partial }_x^{4}A+iA|A{|}^2+\rho$$

$$\begin{array}[]{l}d_{x}y_{i}=y_{i+2},\hskip 15.6491pti=1,...,6\\ d_{x}y_{7}=d_{4}^{-1}\left[y_{2}+\theta y_{1}-y_{5}-y_{1}(y_{1}^{2}+y_{2}^{2})\right]\\ d_{x}y_{8}=d_{4}^{-1}\left[-y_{1}+\theta y_{2}-y_{6}-y_{2}(y_{1}^{2}+y_{2}^{2})\right].\end{array}$$

$$d_4^2{\lambda }^8+2d_4{\lambda }^6+c_4{\lambda }^4+c_2{\lambda }^2+c_0=0,$$

$${\lambda }^{-}=\pm q^{-}\pm i{k}^{-}=\pm \sqrt{\frac{\sqrt{c_4}-1}{2d_4}}\pm i\sqrt{\frac{\sqrt{c_4}+1}{2d_4}}.$$

$$\Lambda =2\pi /k^{-}=2\pi \sqrt{\frac{2d_4}{\sqrt{c_4}+1}},$$

$$-v{d}_{x}A=-(1+i\theta )A+i{d}_x^{2}A+d_3{d}_x^{3}A+iA|A{|}^2+\rho ,$$

Full text

I Introduction

Localized dissipative structures (LSs) are a particular type of emergent structures appearing in spatially extended systems out of the thermodynamical equilibrium, i.e. systems that are internally dissipative and externally driven, commonly known as dissipative systems. They can be understood as part of one state embedded in a background of a different state and are formed due to a double balance between nonlinearity and dispersion (or diffraction) on the one hand, and driving and dissipation on the other Akhmediev_GTS . When they consist of a single localized element, such as a single peak or a localized patterned patch, LSs are also called dissipative solitons or simply solitons. They are unique once the system parameters are fixed and they can exhibit multistability with other LSs Akhmediev_GTS . These structures can display a variety of dynamical regimes such as periodic oscillations Umbanhowar ; Firth ; Vanag_Epstein ; leo_lendert , chaos leo_lendert ; Michaelis , or excitability gomila , and they are found in a wide variety of systems in nature, such as chemical reactions, plant ecology, biology and nonlinear optics Cross_Hohenberg ; Murray ; Hoyle . Moreover, several single LSs can also coexist and form bound states (BSs) that keep a stable separation distance between them Malomed91 ; Malomed93 ; Malomed94 ; Akhmediev ; Akhmediev_Soto_PRL ; Barashenkov ; Schapers .

In nonlinear optics, LSs arise in driven nonlinear optical cavities such as Fabry-Perot interferometers lugiato_spatial_1987 , fiber ring cavities Haelterman ; Leo_NaturePhot , and microresonators Coen_comb ; Chembo_comb . In this context, LSs are well described by a mean-field model, which in the case of Kerr nonlinearities reduces to the Lugiato-Lefever (LL) equation lugiato_spatial_1987 . In the last decade, this model has sparked new interest as it has been found to describe accurately the formation and dynamics of Kerr frequency combs in microresonators Coen_comb ; Chembo_comb . Optical frequency combs in microresonators show great promise thanks to their wide range of applications and their potential to be integrated on-chip kippen ; delhaye ; Okawachi ; hansch ; Scott ; Herr_kippenberg . Therefore, understanding the formation of single solitons, patterns, and bound states, is a crucial step in gaining new insights into optical frequency combs.

Several works by Malomed Malomed91 ; Malomed93 ; Malomed94 , Akhmediev Akhmediev ; Akhmediev_Soto_PRL , and Barashenkov Barashenkov have theoretically studied BSs in perturbed versions of the Nonlinear Schrödinger (NLS) equation using a variety of techniques, such as calculating effective interaction potentials. Even in non-integrable equations, which do not allow finding a closed analytical expression for the interaction potential, one can use numerical solutions of solitons to calculate such a potential. In Malomed91 , it was found that equilibrium distances of BSs correspond to the maxima and minima of that potential. The maxima correspond with the unstable equilibrium separations, while the minima indicate the stable ones. The interaction in the presence of so-called ”skew” terms (terms breaking spatial reversibility) has also been studied Malomed93 ; Malomed94 . In this particular case, the maxima of the potential, instead of the minima, determine the stable equilibrium separation distances for BSs Malomed93 . In even earlier works, it was shown that locking of solitons can be produced via radiative interaction or dispersive wave emission Mollenauer ; Kaup . Interaction of solitons in dissipative systems has also been studied using perturbation theory in the modified Swift-Hohenberg equation Tlidi_SHE ; Tlidi_SHE_interaction . We also note that BSs have been studied in two and three spatial dimensions Malomed98 , and that the solitons within a BS can undergo temporal oscillatory instabilities and move with respect to one another soto_crespo . Experimentally, BSs have been observed in spatial optical cavities Schapers , and more recently in microresonators Herr_kippenberg ; Wang_leo_jae_coen ; Wang_leo_jae_coen_Optica ; Brasch ; Vahala ; DelHaye .

In this paper, using similar techniques as in Malomed91 ; Malomed93 ; Malomed94 ; Akhmediev ; Barashenkov , we study how different LSs can interact and form BSs in the generalized LL equation due to their interaction via the oscillatory tails in the soliton profiles. To do so, we focus on the anomalous group velocity dispersion (GVD) regime. In Section II, we introduce the standard LL model and review how different types of LSs form. Next, in Section III, we show a variety of BSs and discuss the general mechanism by which the solitons can interlock at different separation distances. In the next Section, we then calculate an interaction potential using variational principles, and analyze how the shape of the tails of the LSs changes. Using both approaches, we discuss for which control parameters BSs of various separation distances can form. Afterwards, in Section V, we present different types of bifurcation structures for BSs. Section VI is devoted to the generalization of these results to the LL equation in the presence of high order dispersion effects, in particular third order dispersion (TOD) and fourth order dispersion (FOD), where we also perform simulations in the presence of noise analyzing how two solitons jump between neighboring stable distances (Section 7). Finally, in Section VIII, we end with some concluding remarks.

II Localized structures in the Lugiato-Lefever equation

The standard LL equation with one spatial dimension in the anomalous GVD regime reads

[TABLE]

where, $\rho$, $\theta$ are control parameters representing the normalized injection amplitude and the frequency detuning, respectively lugiato_spatial_1987 ; Haelterman .

Steady state solutions of (1) satisfy the following ordinary differential equation

[TABLE]

and they can be homogeneous steady states (HSS), or non-uniform solutions (with a spatial dependence) consisting of periodic patterns or LSs Coullet .

The HSS solutions $A_{0}$ are given by the classic cubic equation of dispersive optical bistability, namely

[TABLE]

where $I_{0}\equiv|A_{0}|^{2}$. For $\theta<\sqrt{3}$, (3) is monovaluate and hence the system is monostable. However, if $\theta>\sqrt{3}$, (3) is trivaluate and $A_{0}$ has three branches of solutions forming a S-shape bifurcation diagram. The different branches meet at saddle-nodes SNhom,1 and SNhom,2 located at

[TABLE]

In terms of the real and imaginary parts of the field A, the HSSs $A_{0}$ is given implicitly by

[TABLE]

For $\theta<\sqrt{3}$, the homogeneous solution, $A_{0}$, becomes unstable to perturbations of the form $\propto{\rm e}^{\Omega t-ikx}$ at $I_{0}=I_{c}=1$, with critical wavenumber $k=k_{c}=\sqrt{2-\theta}$. Above such threshold the real part of $\Omega$ is positive and perturbations grow. At $I_{c}$ a pattern solution is created either supercritically ($\theta<41/30$) or subcritically ($\theta>41/30$) lugiato_spatial_1987 . If the intracavity intensity $I_{0}$ is above the modulation instability (MI) threshold, $I_{c}=1$, the HSS is unstable to a range of wavevectors, with $k=k_{u}=\sqrt{2I_{0}-\theta}$ the one growing the fastest. For $\sqrt{3}<\theta<2$, while the homogeneous solution is trivalued, the lower branch ($A_{0}^{b}$) has the same instability threshold and critical wavenumber.

In the regime where patterns arise subcritically, the bistability between the low intensity solution ($A_{0}^{b}$, zero in Fig. 1) and the pattern solution ($P$, red branch in Fig. 1) allows the formation of LSs Woods_1999 ; gomila_Scroggi ; Burke_2006 . Together with the pattern, two branches of LSs bifurcate from the MI. One branch includes profiles with an odd number of peaks (blue curve), and the other one contains profiles with an even number of peaks (green curve). These two branches of solutions are intertwined in a so-called homoclinic snaking bifurcation structure: a sequence of saddle-node bifurcations cause the branches to oscillate back and forth across a parameter range called the snaking or pinning region Woods_1999 ; Burke_2006 . The LSs repeatedly add peaks on either side symmetrically at each back-and-forth oscillation, increasing the width of the LS by two wavelengths $2\pi/k_{c}$, until filling the whole size and connecting to the pattern solution $P$. This type of bifurcation structure in the LL equation is shown in Fig. 1 for a domain size $L=160$ and for a fixed value of the detuning $\theta=1.5$, and was discussed in more detail in Parra_Rivas_2 .

For higher values of the detuning ($\theta>2$), a region of operation that is more commonly used in frequency comb generation in microresonators, the situation is considerably different. The low intensity HSS, $A_{0}^{b}$, is now stable all the way up to SNhom,1, thus $I_{b}=\frac{2\theta}{3}+\frac{1}{3}\sqrt{\theta^{2}-3}$. When increasing the detuning the typical soliton becomes sharper and the oscillatory tails decrease in amplitude (see Fig. 2). Now $I_{c}=1$ corresponds to a Belyakov-Devaney transition Parra_Rivas_2 . For $I_{0}<I_{c}$, LSs still maintain oscillatory tails, although they are highly damped. However, as soon as $I_{0}>I_{c}$, the tails become monotonic.

III Origin of bound states

In Fig. 3, we show a single stable soliton, as well as two different stable BSs consisting of two solitons bound at different separation distances. The maxima of all solitons are indicated by red lines, showing that the maxima of each soliton in the BS closely correlate with the maxima of the oscillatory tails of the single soliton. Indeed, it was analyzed in a general way that different LSs interact through their tails and can form bound states in this way Elphick ; Aranson . Moreover, it was shown that such tail interaction is not only important between different LSs, but LSs can also interact with system boundaries through their tails Kozyreff1 ; Kozyreff2 .

This general mechanism of forming BSs through tail interaction has clear implications in the context of the LL equation. In Fig. 2, we illustrated that oscillatory tail amplitude damps increasing the frequency detuning $\theta$, weakening the interaction between solitons. As a result, BSs cannot be strongly locked together and are very sensitive to perturbations. The formation of BSs becomes even more cumbersome for values of the frequency detuning larger than 2 where tails become monotonic. The interaction via monotonic tails leads to a continuous, albeit very slow, movement of the solitons, so that they no longer form any stable BSs. In recent experiments, however, BSs have been found for $\theta>2$, showing that the standard LL equation does not capture all of the dynamics.

Any perturbation to the LL equation that introduces oscillatory tails into the soliton profile would be sufficient to explain the formation of BSs. In what follows, we will first discuss in detail BSs in the standard LL equation for $\theta<2$. Next, we will extend our study to BSs in the generalized LL equation in the presence of higher order dispersion effects, which is one way to introduce oscillatory tails GelensPRL ; Gelens_NL .

IV Interaction potential and oscillatory tails

IV.1 Variational methods

Here, using similar variational methods as in Malomed91 ; Malomed93 ; Malomed94 ; Akhmediev ; Barashenkov , we study how solitons interact and form BSs. The LL equation with one spatial dimension has the following generalized action functional:

[TABLE]

where $\bar{A}$ is the complex conjugate of the field $A$. The time-dependent exponential is linked to the dissipative nature of the system Barashenkov . In (6), the Lagrangian density is given by Firth_lagran :

[TABLE]

In this framework the LL equation corresponds to the Euler-Lagrange equation, derived from the least action principle defined by:

[TABLE]

where $\delta$ stands for functional derivative.

The Hamiltonian density is,

[TABLE]

and the interaction Hamiltonian density is given only by the last term

[TABLE]

The Hamiltonian function is then defined by

[TABLE]

We now use this formalism to derive effective interaction potentials of two solitons that are separated by a distance $z$. This potential is determined by the overlapping integral between the tail of one soliton and the core of the other soliton. We consider that the BS formed by two solitons that are widely separated can be described by the ansatz

[TABLE]

where $A_{-}(x)=A_{sol}(x-z/2)$ represents a soliton displaced by a distance $z/2$ to the left of the center of the domain at $x=0$, and $A_{+}=A_{sol}(x+z/2)$ a soliton displaced by a distance $z/2$ to the right, with $z$ being a time-dependent free parameter Barashenkov . In Ref. Akhmediev_Soto_PRL a different ansatz including the phase difference between solitons was considered. In the case considered here the phase is fixed by the external pump and there is no need of such free parameter. As there is not a known exact analytical solution for the soliton in the LL equation, we use the stationary solution profiles obtained numerically via a Newton method to calculate the potential.

IV.2 Full Hamiltonian

Using the ansatz (12), the action functional depends on $z$, i.e. $\mathcal{S}(z)\equiv\mathcal{S}[A(z),\bar{A}(z)]$. This time independent ansatz implies that the kinetic terms in the Lagrangian density (7) vanish and thus the action corresponds to the integral of the Hamiltonian density (9) Barashenkov . The extreme of the action correspond to stationary equilibrium distances or pinning distances:

[TABLE]

which is equivalent to

[TABLE]

This equation has the form of a constraint and its solutions $z=z_{n}$ correspond to the equilibrium separation distances for the BSs. $H$ as a function of $z$ defines an effective potential

[TABLE]

IV.3 Interaction Hamiltonian

Another constraint equation similar to (14) can be obtained using the interaction Hamiltonian (10)

[TABLE]

instead of the complete one given by (9). This Hamiltonian defines an effective interaction potential

[TABLE]

IV.4 Effective potentials

Fig. 4 shows the potential corresponding to the soliton shown in Fig. 3, calculated using the full Hamiltonian density (9) (top) and the interaction Hamiltonian density (10) (bottom), integrated over the range $[-L/2,L/2]$. In both cases the potentials oscillate in $z$ with a fixed period or wavelength $\Lambda_{U}$ and decays for increasing values of $z$. The insets in Fig. 4 show that even at large distances the oscillations are still present. The minima (maxima) of the potential correspond to stable (unstable) equilibrium separations $z_{n}^{s}$ ($z_{n}^{u}$), and therefore to stable or unstable BSs. These minima (maxima) are indicated by $\bullet$ ($\circ$) and by the red (black) dashed lines. The numerical values for some of these stable points using both calculated potentials are shown in Table 1. For comparison we have also added the exact separation distances calculated using numerical time-evolution simulations. Although the agreement is not perfect, the differences are relatively small and of the order of those found in Barashenkov , where a similar equation was used. They can be attributed to the approximation (12). The periodicity of this potential $\Lambda_{U}=z^{s}_{n+1}-z^{s}_{n}$, i.e. the difference between the position of two consecutive local minima (or maxima) can be calculated for every set of control parameters. Here, we find that $\Lambda_{U}\in(9.1,9.2)$.

IV.5 Wavelength of the oscillatory tails

Another approach to finding the separation $\Lambda_{U}$ between consecutive stable separation distances of two bound solitons, is to study the spatial dynamics of the system Champneys_homoclinic ; Gelens_NL . Defining the variables $y_{1}(x)={\rm Re}[A(x)]$, $y_{2}(x)={\rm Im}[A(x)]$, $y_{3}(x)=d_{x}{\rm Re}[A]$ and $y_{4}(x)=d_{x}{\rm Im}[A]$, (2) can be recast as the four dimensional dynamical system given by

[TABLE]

In this framework, LSs biasymptotic to the HSS $A_{0}^{b}$ correspond to homoclinic orbits to $A^{b}_{0}$ Holmes . By studying the linearization of the system (18) around $A_{0}^{b}$, and therefore the spectrum of eigenvalues of its Jacobian matrix, it is possible to understand how trajectories leave and approach $A_{0}^{b}$. The spatial eigenvalues of the Jacobian matrix are solutions of the characteristic polynomial

[TABLE]

with $c_{0}=\theta^{2}+3I_{0}^{2}-4\theta I_{0}+1$ and $c_{2}=4I_{0}-2\theta$. This equation is invariant under $\lambda\rightarrow-\lambda$ and, since the coefficients are real, it leads to eigenvalues whose configuration in the complex plain is symmetric with respect to both axes. The form of this equation is a consequence of spatial reversibility Devaney_A ; Homburg ; Knobloch15 . The eigenvalues satisfying (19) are

[TABLE]

For any value of $\theta$, $A^{b}_{0}$ is a saddle-focus for $I_{0}<I_{c}$, with eigenvalues $\lambda_{1,2,3,4}=\pm q_{0}\pm ik_{0}$, where

[TABLE]

and

[TABLE]

In this case, trajectories approach or leave $A_{0}^{b}$ in an oscillatory manner, and therefore, the tails of the soliton profile are oscillatory in space. Thus, in the linear regime, the oscillatory tails of the soliton profile are approximated by

[TABLE]

where $a_{1}$ and $\varphi_{1}$ are functions of the control parameters $(\theta,\rho)$ determined by the complete nonlinear dynamics, and they can be estimated by fitting the tails with function (23). For the control parameter values used here, $(\theta,\rho)=(1.5,1.11445)$, the frequency and the decay rate are respectively $k_{0}=0.6853$ and $q_{0}=0.1827$, and the fitted parameters are $(a_{1},\varphi_{1})=(1.45\cdot 10^{-7},-0.053)$. The imaginary part of the tails is described by a function analogous to (23).

From (22), one finds that the wavelength of the oscillatory decaying tails is given by

[TABLE]

This expression is plotted in Fig. 5 as a function of $\rho$ for $\theta=1.5$. The red diamonds correspond to the periodicity of the potential $\Lambda_{U}$, calculated using (15) in the pinning region of the homoclinic snaking. The fact that $\Lambda_{U}\approx\lambda_{0}^{tails}$ shows that both approaches (interaction potential vs. spatial eigenvalues) can be used to calculate stable separation distances from an initial stable separation distance $z^{s}_{n}$, by using the relation $z^{s}_{n+1}=z^{s}_{n}+\Lambda_{U}$.

V Bifurcation structure

In this Section, we explore the bifurcation structure of some of the BSs that can be formed combining LSs consisting of one or multiple peaks locked at different separation distances, such as those shown in Fig. 3. A detailed study of the bifurcation structure of BSs has been done in Multi_Burke in the context of the Swift-Hohenberg equation. In the previous Section, we showed that the equilibrium separations between two soliton states are given by the overlapping of the tail of one soliton with the core of the other one. The presence of stable BSs is intrinsically determined by the properties of the oscillatory tails. Therefore, for a fixed set of parameters $(\theta,\rho)$, one expects to have as many BSs as one can make combinations of existing LSs. In an infinite domain, the number of possible BSs is infinite, while in a finite size system the number of BSs is constrained by the domain size.

First, let us discuss the bifurcation structure of BSs consisting of two single solitons separated by $L/2$. As these pulses are equidistant to the nearest neighbor on either side, their behavior is identical to the behavior of single pulses on a periodic domain of size $L/2$. These type of structures are, therefore, organized in a snaking bifurcation diagram as shown in Fig. 6. Both pulses behave similarly as a single pulse previously shown in Fig. 1, and so each individual pulse of the BS adds extra peaks after each saddle-node bifurcation.

As shown in Fig. 3, BSs not only exist with a stable separation distance $L/2$, but solutions with various other stable separations $z^{s}_{n}$ (see Table 1) exist. Similarly, BSs consisting multiple peaks separated by $z^{s}_{n}$ also exist. Such solutions share the same bifurcation structure as the LSs separated by $L/2$, shown in Fig. 6, as long as both LSs within the BS remain far enough to only interact weakly and behave as independent entities.

The branches associated to BSs corresponding to LSs that are close to one another are, however, organized into large isolas, and do not connect directly with the pattern branches. Such a bifurcation structure of isolas is shown in Fig. 7 for BSs separated by $z^{s}_{1}$. This kind of behavior was reported in Multi_Burke to originate from insufficient accuracy in the continuation algorithm, which generates jumps between independent isolas, creating one large one as seen in Fig. 7. Despite improving the accuracy in the numerical continuation algorithm, we have not found such independent isolas in the LL equation, instead we always find several connected isolas in a structure similar to that of Fig. 7.

In addition to the symmetric BSs consisting of two LSs, there is a wide variety of BSs that are built up from any number of pulses randomly separated by the pinning distances $z^{s}_{n}$. Many of these states are not invariant under reflection symmetry, i.e. are asymmetric, and in the LL equation they move with constant velocity as a consequence of non-variational effects. If the separation between the peaks of these states is large, they form a snaking-type bifurcation diagram (see Fig. 6 in Multi_Burke ). However, when the separation between LSs is smaller, these BSs are also organized in isolas. Fig. 8 shows two examples of such isolas. Panel (a) shows the isola corresponding to a four-pulse BS whose peaks are separated by a distance $z^{s}_{2}$. Panel (b) shows a similar structure but for a 4-pulse BS where the inter-distance between peaks is also given by $z^{s}_{2}$.

VI Oscillatory tails induced by higher order dispersion

In this Section, we study the effect of higher order dispersion on the shape of the oscillatory tails of solitons, and as a result, on the stability and shape of BSs. Taking into account higher order dispersion effects is critical when considering cavities that operate close to zero GVD, a situation which has been targeted experimentally to achieve broader frequency combs Lamont ; Brasch . Theoretically, various works have addressed the influence of higher order dispersion on the formation and stability of solitons Coen_comb ; GelensPRA2007 ; Tlidi_FOD ; Milian1 ; Pedro_TOD ; Milian2 .

VI.1 Solitons in the generalized LL equation with higher order dispersion

The LL equation in the anomalous regime including dispersion terms of third and fourth order reads

[TABLE]

Fig. 9 shows several stable soliton solution profiles for $(\theta,\rho)=(2.2,1.5)$, including either fourth order dispersion (FOD, $d_{4}\neq 0$) or third order dispersion (TOD, $d_{3}\neq 0$), each time varying the strength of this higher order dispersion term. In both cases, adding these dispersion terms of higher order introduces oscillatory tails where previously (only second order dispersion, and $\theta>2$) they were highly damped. Initially (small FOD or TOD), these oscillations have a small amplitude and very short wavelength, but when increasing the higher order dispersion strength further these oscillatory tails become more pronounced and their wavelength increases. One critical difference between FOD and TOD is the symmetry of the resulting soliton profiles. While FOD retains reversibility in the system, TOD breaks such reversibility. As a result, solitons in systems with FOD remain symmetric, while in the presence of TOD they become asymmetric and start drifting (to the right for $d_{3}>0$).

VI.2 Interaction potential

As before, inserting the ansatz (12) into (17) the potential $U_{H_{I}}$ can be calculated as a function of the separation distance $z$. Also in the presence of higher order dispersion, this leads to an oscillatory interaction potential, where the extrema give an estimate of the equilibrium separation distances of solitons within BSs. Such approximation is in principle only valid when the separation between solitons in the BS is large enough. Therefore, in Table 2, we compare the predictions for the separation distances using the interaction potential with the exact numerical solution, finding good correspondence.

VI.3 Wavelength of the oscillatory tails: fourth order dispersion

Similarly as in Section II, using the variables $y_{1}={\rm Re}[A]$, $y_{2}={\rm Im}[A]$, $y_{3}=d_{x}{\rm Re}[A]$, $y_{4}=d_{x}{\rm Im}[A]$, $y_{5}=d_{x}^{2}{\rm Re}[A]$, $y_{6}=d_{x}^{2}{\rm Im}[A]$, $y_{7}=d_{x}^{3}{\rm Re}[A]$ and $y_{8}=d_{x}^{3}{\rm Im}[A]$, one can recast the stationary version of (25) to the following dynamical system:

[TABLE]

The eigenspectrum of the Jacobian of (26) evaluated at $A_{0}$ is given by the characteristic polynomial

[TABLE]

with $c_{4}=(1-2d_{4}\theta+4d_{4}I_{0})$. The solution of this polynomial consists of two sets of eigenvalues $\lambda_{1,2,3,4}=\pm q^{+}\pm ik^{+}$ and $\lambda_{5,6,7,8}=\pm q^{-}\pm ik^{-}$, as those shown in Fig. 10 for $d_{4}=0.02$ (a) and $d_{4}=0.15$(b). Due to the spatial reversibility the eigenspectrum is symmetric with respect to the axes Im$[\lambda]=0$ and Re$[\lambda]=0$. This eigenspectrum allows estimating the wavelength of the oscillatory tails by $2\pi/k^{-}$, where $k^{-}$ is the imaginary part of the spatial eigenvalue with smallest negative real part in absolute value ($q^{-}$). Fig. 11(a) shows that the wavelength of the oscillatory tails increases with $d_{4}$. Moreover, it shows that this wavelength (red line with crosses) corresponds well with an estimate of the periodicity of the interaction potential (blue diamonds), again confirming the validity of both approaches.

For low values of $d_{4}$, the modulus of the spatial eigenvalue $\lambda^{-}$ is very large (see Fig. 10), in such a way that the dominant terms in (27) are those with the highest order in $\lambda$. Due to this it is possible to obtain an analytical approximation for those eigenvalues by solving (27) only considering the highest order in $\lambda$. In this way, we obtain that $\lambda^{-}$ can be approximated by the expression

[TABLE]

This expression shows that the eigenvalue approaches infinity ($\lambda^{-}\rightarrow\infty$) when FOD becomes zero ($d_{4}\rightarrow 0$), a result that can be observed looking at the spectrum for different values of $d_{4}$. When $d_{4}\rightarrow 0$, these two eigenvalues $\lambda^{-}$, and its complex conjugate, tend to $\pm i\infty$. The approximation of the wavelength of the oscillatory tails using this expression (28) gives:

[TABLE]

and is shown in Fig. 11(a) by the black line with asterisks.

VI.4 Wavelength of the oscillatory tails: third order dispersion

Following the analysis for FOD, we now write down the LL equation with TOD and look for solutions that move rigidly at a velocity $v$ (to be determined), namely $A(x,t)=A(x-vt)$. From (25) with $d_{4}=0$ one has:

[TABLE]

which as in the previous section can be recast into:

[TABLE]

Steadily drifting soliton solutions of (30) can be found numerically, with arbitrary precision, using a Newton-Raphson method where the velocity $v$ is computed as part of the solution. Then, the eigenspectrum of the Jacobian of (31) around $A_{0}^{b}$ with the found velocity is given by the characteristic equation:

[TABLE]

The corresponding spatial eigenvalues are shown in Fig. 10 for $d_{3}=0.1$ (c) and $d_{3}=0.25$ (d). As reversibility is broken by TOD, the spatial eigenvalues no longer appear in quadruplets, namely complex eigenvalues appear in conjugate pairs while real eigenvalues appear non-symmetrically around the imaginary axis. Similar as in the FOD case, however, the wavelength of the oscillatory tails as determined by the dominant spatial eigenvalue $\lambda^{-}$ increases with TOD strength $d_{3}$, as shown in Fig. 11(b).

We can obtain an analytical approximation for the wavelength of the oscillatory tails $\Lambda$ using a similar analysis as in the previous FOD case, considering the highest order in $\lambda$:

[TABLE]

which is again plotted in Fig. 11(b) as a black line with asterisks. We note that the approximation (33) gives very accurate results over a broad range of values of the parameter $d_{3}$, being quite more accurate than the equivalent approximation in the case of FOD.

VII Effects of noise on bound states

In this section we study the dynamics of BSs in the presence of noise. The effect of a fluctuating driving term on the dynamics of BSs was studied in the NLS equation with a linear dissipative term, where the noise was uniform in space and varying in time Malomed95random . In contrast, in our case, noise is incorporated in the system through a term that describes a fluctuating pump intensity which is random both in time and space, and that has the form

[TABLE]

where $\xi(x,t)$ is a Gaussian white noise with zero mean $\langle\xi(x,t)\rangle$=0 and correlations

[TABLE]

where $\langle\cdots\rangle$ stands for the mean value.

In Ref. Malomed95random a Fokker-Planck equation was derived for the probability distribution of the amplitude of a single soliton. Using the stationary solution of such equation, a mean potential for the soliton-soliton interaction was calculated, showing that the bound distances were modified with respect to the deterministic case. In our case, the equilibrium separations are not substantially modified and we focus on the jumps between neighboring stable fixed points due to the noise.

Fig. 12 shows the time evolution of BSs consisting of two solitons starting from different initial separations, both in a reversible system (with FOD: $d_{4}=0.025$ (a)) and a system with broken reversibility (with TOD: $d_{3}=0.3$ (b)). Due to noise, the separation fluctuates around stable locking distances ($z^{s}_{n}$) indicated by the red lines, which is especially clear in the zooms. Sudden larger jumps increasing or decreasing the separation between peaks occur. They correspond to jumps between neighboring equilibrium positions due to the noise. The distance between the stable (red dashed lines) and unstable (black dashed lines) separations is given by $\Lambda_{U}/2=z^{s}_{n}-z^{u}_{n}$, which, for the parameters considered in the figure, corresponds to $\approx 0.440$ for the FOD case (a) and to $\approx 0.754$ for the TOD case (b). A general trend is that these jumps occur much more often as the separation distance increases, which is expected as the amplitude of the oscillatory tails, and thus the interaction strength, strongly decreases, making the “potential” barriers lower. There is, however, a significant difference between the reversible (FOD) system and the system with broken reversibility (TOD). In the reversible system two solitons always tend to increase their separation over long time scales, eventually approaching a BS where both solitons are equidistantly spaced over the domain of size $L$, which is illustrated by the two-soliton profiles taken at the trajectory and times indicated by the blue dots. In the presence of TOD, however, trajectories where two solitons of a BS approach are also commonly observed, and no net drift towards larger or smaller distances is clearly appreciated.

In the case of FOD, Fig. 12(a) shows two different time scales, a fast one corresponding to the jumps between neighboring stable positions, and a slow one describing how for long times the system evolves to larger separation distances. In an infinite system the separation between the two solitons will increase indefinitely. However, our system is periodic and the largest separation allowed is $z=L/2$. The fact that solitons are not point objects but have a spatial extension allow them to feel the envelope of other solitons. Since in this case solitons move down gradients, two solitons tend to move apart as a result of the interaction through their envelopes. In fact $z=L/2$ correspond to the minimum of such envelope and a minimum of the interaction potential.

In contrast, in the case of TOD solitons move up gradients and it has been shown that the maxima of the interaction potential (and no longer the minima) determine the stable separation distances Malomed93 . As a result, in principle, one could expect the system now to evolve to smaller separation distances corresponding to the maximum of the envelope of the potential. However, we do not clearly observe such dynamics (see Fig. 12(b)). Instead no net drift of the separation is observed in our stochastic simulations with TOD.

The reason behind the different behavior for FOD and TOD can be understood as follows: in addition to the effect of the envelope, another mechanism must be taken into account to explain the dynamics of two interacting solitons in the presence of noise. Note that the wells in the modulated potential shown in Fig. 4 are asymmetric: the barrier to larger distances is lower than to smaller distances. As a result, the noise is rectified and jumps to larger distances happen more often than to shorter distances, inducing a net movement separating the two solitons. This phenomenon is similar to the Brownian ratchet Magnasco . This mechanism reinforces the tendency for two solitons to separate in the reversible case (FOD). However, it opposes the envelope effect that pushes solitons towards each other in the presence of TOD (broken reversibility). We observe that the two opposite forces effectively cancel reducing much the net drift, and leading to solitons that seem to wander around randomly.

Considering a large number of different realizations of the noise, and starting from a given initial separation $z_{0}$, we have calculated the probability $P$ of finding the two solitons separated a distance $z$ after a time $t_{f}=30000$ (see Fig. 13). In the case of FOD (panels (a) and (b)) we observe a net drift of the mean position until half of the system size is reached. The variance of the fluctuations grows with the distance due to the potential wells becoming more shallow. In the case of TOD, the envelope and ratchet effects cancel each other, leading to a nearly diffusive motion of the solitons with zero mean displacement, and variance that grows as the square root of time. The fine structure of the probability distribution also reveals the stable and unstable positions. These are more clear in the case of TOD (Fig. 13d), where for the parameters of the figure the oscillatory tails are more pronounced (see Fig. 9).

VIII Conclusions

In this work, we have studied the interaction and formation of BSs of localized structures in the LL equation with and without additional higher order dispersion terms. Using variational techniques previously used in Malomed91 ; Malomed93 ; Malomed94 ; Akhmediev ; Barashenkov , we have derived an effective potential depending on the separation distance between solitons in the BS. The extrema of this potential determine the stable and unstable separation distances, and they are related with the overlapping of the oscillatory tails of one soliton and the other soliton’s core. Therefore, we have also analyzed BSs by studying the eigenvalues of the spatial dynamics for the LL equation, i.e. the way the system approaches or leaves the homogeneous steady state solutions. The periodicity of the potential is determined by the wavelength of a soliton’s oscillatory tails. We have found that the variational approach and the spatial eigenvalue analysis are consistent, and both give accurate estimates of the separation between two consecutive equilibrium distances of BSs. Moreover, the variational approach predicts the absolute allowed separation distances.

Next, we have calculated the bifurcation diagrams associated to different types of BSs. BSs that consist of two solitons, which are separated enough such that they interact very weakly through their oscillatory tails, behave as independent LSs, resulting in a snaking-type bifurcation diagram Multi_Burke . In contrast, when the separation is smaller such that tail interaction is more significant, BSs are organized in stacks of isolas that are no longer connected to the pattern solution branch. If we consider arrays of more than two solitons, the bifurcation structures obtained are usually isolas.

We have then extended our analysis to the LL equation with higher order dispersion effects, which cannot be ignored when operating close to zero second-order GVD. We considered two qualitatively different situations: (i) fourth order dispersion, which maintains reversibility, and (ii) third order dispersion, which breaks reversibility leading to asymmetric, moving solitons. Both types of higher order GVD can introduce oscillations into the spatial profile of the soliton’s tails, even for high values of the detuning, where in the standard LL equation no oscillatory tails were present ($\theta>2$). By calculating the interaction potentials and characterizing the spatial eigenvalues of the system, we have shown that a wide range of BSs locked at different separation distances come into existence.

Finally, in Section VII, we have studied the effects of adding white Gaussian noise to the pump intensity, which lead to random jumps between BSs of different separation distances. By pooling many time evolution simulations for different noise realizations, we have calculated the probability of finding two solitons at a given distance after a given time. We have found that there exists a critical difference between systems that are reversible (fourth order dispersion) and those that are not (third order dispersion). In the reversible case, noise tends to drive the solitons within a BS apart over time, eventually leading to a situation where a BS consisting of two solitons separated by half of the system width is the most probable one. In contrast, in the non-reversible case, the potential envelope and ratchet effect cancel out leading to a diffusive like behavior that leads the two solitons wandering around randomly.

Recently, BSs of solitons have been observed experimentally in semiconductor lasers Barland , and in microresonators in the context of frequency comb generation Herr_kippenberg ; Wang_leo_jae_coen ; Wang_leo_jae_coen_Optica ; Brasch ; Vahala ; DelHaye . Moreover, several binding mechanisms have been studied theoretically and experimentally in the context of passively driven nonlinear optical resonators, such as Gordon/Kelly sidebands, birefringence, and dispersive waves Wang_leo_jae_coen ; Wang_leo_jae_coen_Optica . We expect that our results will prove useful to help interpreting various of these experimental observations of BSs.