# Stability of Multi-hump Optical Solitons

## Abstract

We demonstrate that, in contrast with what was previously believed, multi-hump solitary waves can be stable. By means of linear stability analysis and numerical simulations, we investigate the stability of two- and three-hump solitary waves governed by incoherent beam interaction in a saturable medium, providing a theoretical background for the experimental results reported by M. Mitchell, M. Segev, and D. Christodoulides \[Phys. Rev. Lett. 80, 4657 (1998)\].

Self-guided optical beams, or spatial optical solitons, are the building blocks of all-optical switching devices where light itself guides and steers light without fabricated waveguides . In the simplest case, a spatial soliton is created by one beam of a certain polarization and frequency, and it can be viewed as a self-trapped mode of an effective waveguide it induces in a medium . When a spatial soliton is composed of two (or more) modes of the induced waveguide , its structure becomes rather complicated, and the soliton intensity profile may display several peaks. Such solitary waves are usually referred to as multi-hump solitons; they have been found for various nonlinear models of coupled fields .

In realistic (nonintegrable) physical models, solitary waves can become unstable demonstrating self-focusing, decay, or a nonlinearity-driven transition to a stable state, if the latter exists . All these scenarios of soliton evolution are initiated by exponentially growing perturbations and they are attributed to linear instability. It is usually believed that all types of multi-hump solitary waves are linearly unstable, except for the special case of neutrally stable solitons in the integrable Manakov model . On the contrary, recent experimental results  indicate the possibility of observing stationary structures resembling multi-hump solitary waves. This naturally poses a question: Were those observations only possible because of short propagation distance and a small instability growth rate? Definitely, the experimental results challenge the conventional view on multi-hump solitary waves in different models of nonlinear physics.

The purpose of this Letter is twofold. First, we study the origin of multi-hump solitons supported by incoherent interaction of two optical beams in a photorefractive medium. We find that multi-hump solitons appear via bifurcations of one-component solitons and due to the process of hump multiplication, when the intensity profile of a composite soliton changes from single- to multi-humped with increasing power. Second, we perform numerical stability analysis of two- and three-hump solitary waves and also find analytically the instability threshold for two-hump solitons. We reveal that two-hump solitary waves are linearly stable in a wide region of their existence, whereas all three-hump solitons are linearly unstable, and that even linearly stable multi-hump solitons may not survive collisions.

In the experiments , spatial multi-hump solitary waves were generated by incoherent interaction of two optical beams in a biased photorefractive crystal. The corresponding model has been derived by Christodoulides et al. , and it is described by a system of two coupled nonlinear equations for the normalized beam envelopes, $`u(x,z)`$ and $`w(x,z)`$, which for the purpose of our current analysis can be written in the following form :

$$\begin{array}{c}i\frac{u}{z}+\frac{1}{2}\frac{^2u}{x^2}+\frac{u(|u|^2+|w|^2)}{1+s(|u|^2+|w|^2)}u=0,\hfill \\ i\frac{w}{z}+\frac{1}{2}\frac{^2w}{x^2}+\frac{w(|u|^2+|w|^2)}{1+s(|u|^2+|w|^2)}\lambda w=0,\hfill \end{array}$$
(1)   
where the transverse, $`x`$, and propagation, $`z`$, coordinates are measured in the units of $`(L_d/k)^{1/2}`$ and $`L_d`$, respectively, $`L_d`$ is a diffraction length, and $`k`$ is the wavevector in the medium. The parameter $`\lambda `$ is a ratio of the nonlinear propagation constants, and $`s`$ is an effective saturation parameter. For $`s0`$, the system (1) reduces to the integrable Manakov equations .

We look for stationary, $`z`$-independent, solutions of Eqs. (1) with both components $`u(x)`$ and $`w(x)`$ real and vanishing as $`|x|\mathrm{}`$. Different types of such two-component localized solutions, existing for $`0<\{\lambda ,s\}<1`$, can be characterized by the total power, $`P(\lambda ,s)=P_u+P_w`$, where the partial powers, $`P_u=_{\mathrm{}}^{\mathrm{}}|u|^2𝑑x`$ and $`P_w=_{\mathrm{}}^{\mathrm{}}|w|^2𝑑x`$, are integrals of motion. If one of the components is small, i.e. $`w/u\epsilon `$, Eqs. (1) become decoupled and, in the leading order, the equation for the $`u`$-component has a solution $`u_0(x)`$ in the form of a fundamental, $`sech`$-like, soliton with no nodes. The second equation can then be considered as an eigenvalue problem for the “modes” $`w_n(x)`$ of a waveguide created by the soliton $`u_0(x)`$ with the effective refractive index profile $`u_0^2(x)/[1+su_0^2(x)]`$. Parameter $`s`$ determines the total number of guided modes and the cut-off value for each mode, $`\lambda _n(s)`$. Therefore, a two-component vector soliton $`(u_0,w_n)`$ consists of a fundamental soliton and an $`n`$th-order mode of the waveguide it induces in the medium. Henceforward we denote such a composite solitary wave by its “state vector”: $`|0,n`$.

On the $`P(\lambda )`$ diagram (for fixed $`s`$), continuous branches representing $`|0,n`$ solitons emerge at the points of bifurcations $`\lambda _n(s)`$ of one-component solitons (see Fig. 1). It is noteworthy that the first-order mode is in fact the lowest possible mode of the waveguide induced by the fundamental soliton $`u_0(x)`$. This is because the state $`|0,0`$, node-less in both components, can exist only in the degenerate case $`\lambda =1`$, when Eqs. (1) have a family of equal-width solutions $`u_0=A(x)\mathrm{sin}\theta `$ and $`w_0=A(x)\mathrm{cos}\theta `$, with arbitrary $`\theta `$, and amplitude $`A`$ satisfying the scalar equation, $`dA/dx=\pm s^1[\mathrm{log}(1+sA^2)s(1s)A^2]^{1/2}`$.

Additionally, indefinitely many families of vector solitons $`|m,n`$, where $`mn0`$, can be formed as bound states of phase-locked $`|0,n`$ solitons . Although such states do contribute to the rich variety of the multi-hump solitons existing in our model, we exclude them from our present consideration.

Families of vector solitons can be found by numerical relaxation technique. Some results of our calculations are presented in Fig. 1, for $`|0,1`$ and $`|0,2`$ solitons found at $`s=0.8`$. Observing the modification of soliton profiles with changing $`\lambda `$ (see inset in Fig. 1), one can see that the modal description of two-component solitons is valid only near bifurcation points. For $`\lambda \lambda _n`$, the amplitude of an initially small $`w`$-component grows and the soliton-induced waveguide deforms. It is this purely nonlinear effect that gives rise to the existence of multi-hump solitons. In particular, two- and three-hump solitons are members of the soliton families $`|0,1`$ (branch A-B-C) and $`|0,2`$ (branch D-E-F) originating at different bifurcation points. At $`\lambda \lambda _n(s)`$, while the $`w`$-component remains small, all $`|0,n`$ solitons are single-humped, as shown in Figs. 1(a,d). As the amplitude of $`w`$ grows with increasing $`\lambda `$, the total intensity profile, $`I(x)=u_0^2(x)+w_n^2(x)`$, develops $`(n+1)`$ humps \[see Figs. 1(b,e)\], and at sufficiently large $`\lambda `$ the $`u`$-component itself becomes multi-humped \[Figs. 1(c,f)\]. The separation distance between the soliton humps tends to infinity as $`\lambda 1`$.

To analyze the linear stability of multi-hump solitons, we seek solutions of Eqs. (1) in the form of weakly perturbed solitary waves: $`u(x,z)=u_0(x)+\epsilon [F_u(x,z)+iG_u(x,z)]`$ and $`w(x,z)=w_n(x)+\epsilon [F_w(x,z)+iG_w(x,z)]`$, where $`\epsilon 1`$. Setting $`F_{u,w}f_{u,w}(x)e^{\beta z}`$, $`G_{u,w}g_{u,w}(x)e^{\beta z}`$, one can obtain the following eigenvalue problem (EVP)

$`\widehat{}_1\widehat{}_0\stackrel{}{g}=\mathrm{\Lambda }\stackrel{}{g},`$ (2)   
$`\widehat{}_0\widehat{}_1\stackrel{}{f}=\mathrm{\Lambda }\stackrel{}{f}.`$ (3)   
Here $`\stackrel{}{g}(g_u,g_w)^T`$, $`\stackrel{}{f}(f_u,f_w)^T`$, $`\mathrm{\Lambda }=\beta ^2`$, and

$`\widehat{}_{0,1}=\left(\begin{array}{cc}\frac{1}{2}\frac{d^2}{dx^2}+1a_{0,1}& b_{0,1}\\ b_{0,1}& \frac{1}{2}\frac{d^2}{dx^2}+\lambda c_{0,1}\end{array}\right),`$

where $`a_0=c_0=I/(1+sI)`$, $`b_0=0`$, $`a_1=a_0+2u_0^2/(1+sI)^2`$, $`c_1=c_0+2w_n^2/(1+sI)^2`$, and $`b_1=2u_0w_n/(1+sI)^2`$.

Because $`\widehat{}_1\widehat{}_0`$ and $`\widehat{}_0\widehat{}_1`$ are adjoint operators with identical spectra, we can consider the spectrum of only one of these operators, e.g. $`\widehat{}_1\widehat{}_0`$. Considering the complex $`\mathrm{\Lambda }`$-plane, it is straightforward to show that $`\mathrm{\Lambda }(\mathrm{},\lambda ^2)`$ is a continuum part of the spectrum with unbounded eigenfunctions. Stable bounded eigenmodes of the discrete spectrum (the so-called soliton internal modes ) can have eigenvalues only inside the gap, $`\lambda ^2<\mathrm{\Lambda }<0`$. The presence of either positive or complex $`\mathrm{\Lambda }`$ implies soliton instability, because in this case there always exists at least one eigenvalue of the soliton spectrum with $`\mathrm{Re}\beta >0`$.

Numerical solution of the EVP (2) shows that both $`|0,1`$ and $`|0,2`$ types of solitary wave solutions can be stable in a certain region of their existence domain, see Fig. 2. In the case of $`|0,1`$ solitons, the appearance of the instability is related to the fact that close to the curve where the total intensity $`I`$ becomes two-humped \[dashed line in Fig. 2\], a pair of internal modes split from the continuum into the gap. As $`\lambda `$ grows, the corresponding, purely imaginary, eigenvalues $`\beta =\pm i\sqrt{|\mathrm{\Lambda }(\lambda )|}`$ tend to zero, and at a certain critical value $`\lambda =\lambda _{\mathrm{cr}}(s)`$, they coincide at $`\beta =0`$. At this point, an eigenmode with positive eigenvalue $`\mathrm{\Lambda }`$ emerges, thus generating linear instability (see Fig. 3) with the instability growth rate $`\beta =\sqrt{\mathrm{\Lambda }(\lambda )}`$. For $`|0,2`$ solutions, the dynamics of internal modes can not be related in any obvious way with a change in the spatial solitary profiles, nevertheless the scenario of the instability development is similar to that for two-hump solitons. The dependence of $`\beta `$ on $`\lambda `$, for $`|0,1`$ and $`|0,2`$ soliton families giving rise to two- and three-hump solitary waves, is shown in Fig. 3 for $`s=0.3`$ and $`s=0.8`$, respectively. A decline in the instability growth rate as $`\lambda 1`$ (see Fig. 3) is caused by the fact that, in this limit, all multi-hump solitons decompose into a number of the neutrally stable $`|0,0`$ solitons separated by infinitely growing distance. Numerical analysis in the close vicinity of this limit is unfeasible due to lack of computational accuracy.

Note that, within the gap of the continuous spectrum, there exist several soliton internal modes not participating in the development of the linear instability. Analysis of their origin and influence on the soliton dynamics is beyond the scope of the present Letter.

With the aid of analytical asymptotic technique , it is possible to show that a perturbation mode with small but positive eigenvalue, and therefore the linear instability of a general localized solution $`(u,w)`$, appears if the functional $`J(u,w),`$ defined as

$$J=\frac{P_u}{2s}\frac{P_w}{\lambda }\frac{P_w}{2s}\frac{P_u}{\lambda }+\frac{P_u}{s}\frac{P_w}{\lambda }\frac{P_w}{s}\frac{P_u}{\lambda },$$
(4)   
changes its sign. The threshold condition $`J=0`$ is, in fact, the Vakhitov-Kolokolov stability criterion , generalized for the case of two-parameter vector solitons. In this case, it does not necessarily give a threshold of leading instability . Therefore, the presence of other instabilities (which are not associated with the condition $`J=0`$ and can have stronger growth rates) is still possible, as in some other cases .

For two-hump solitons, we have been able to locate the critical curve in $`(\lambda ,s)`$-plane corresponding to the condition $`J=0`$. Superimposing this curve onto the numerically calculated values $`\lambda _{\mathrm{cr}}(s)`$, we have found a remarkable agreement between the numerical and analytical instability thresholds, as shown in Fig. 2. This gives us the first example of the generalized Vakhitov-Kolokolov criterion for the instability threshold of vector multi-hump solitary waves. For the whole family of $`|0,2`$ solutions, including three-hump solitons, it appears that $`J0`$ throughout the entire existence region. Thus, appearance of instability of three-hump solutions is not associated with the change of the sign of the functional $`J`$.

To analyze long-term evolution of multi-hump solitary waves, we perform numerical simulations of the beam propagation for $`|0,n`$ solitons within the existence domain $`\lambda _n<\lambda <1`$, at fixed $`s`$. First, we use no perturbation so that the soliton instability can only develop from numerical noise. As long as the soliton maintains its single-humped shape \[see corresponding profiles in Fig. 1(a,d)\], it remains almost insensitive to numerical noise. Moreover, while the $`|0,1`$ solitons do become two-humped at $`\lambda <\lambda _{\mathrm{cr}}`$, they still remain stable in a wide domain of their parameters until the linear instability threshold is reached. On the contrary, $`|0,2`$ solitons remain single-humped up to the instability threshold value $`\lambda =\lambda _{\mathrm{cr}}`$, so that all three-hump solitons are indeed unstable. Above the instability threshold (i.e. for $`\lambda _{\mathrm{cr}}<\lambda <1`$), a two-hump soliton splits into two independent single-humped beams as a result of the instability developed from noise \[see Fig. 4(a)\], whereas a three-hump soliton exhibits a more complex symmetry-breaking instability, as shown in Fig. 4(b).

Next, we propagate two-hump (at $`s=0.3`$) and three-hump (at $`s=0.8`$) solitons perturbed by an eigenmode with the largest instability growth rates, i.e. $`\beta _{\mathrm{max}}0.055`$ and $`\beta _{\mathrm{max}}0.153`$, respectively. We find that in the presence of $`6\%`$ amplitude perturbation, the diffraction-induced decay of a soliton can be stabilized by the nonlinearity, whereas its splitting is significantly speeded up by the perturbation, compared with splitting due to a numerical noise.

To make a link between our stability analysis and experiment, we note that for the experiment  the diffraction length is defined as $`L_d=2/sb`$ and nonlinearity of the medium (SBN:60 crystal) is characterised by the parameter $`b=kr_{\mathrm{eff}}n_b^2E_0`$, where $`r_{\mathrm{eff}}`$ is the effective electro-optic coefficient ($`=280`$ pm/V), $`n_b`$ is the background refractive index $`(=2.3)`$, and $`E_0`$ is the applied electric field ($`2`$x$`10^5`$ V/m). For strong saturation we have $`s1`$ and $`L_d0.2`$mm. Now, the characteristic instability length $`z_{\mathrm{cr}}`$ can be defined through the maximum growth rate $`\beta _{\mathrm{max}}`$ and, as a result, for two-hump solitons at $`s=0.3`$ we obtain $`z_{\mathrm{cr}}12.18`$ mm. These estimates indicate that the instability, if it exists, could be detected for two-hump solitons within the experimental setup of Ref.  and therefore stable two-hump solitons have been indeed observed.

Importantly, three-hump solitons so far generated in the experiment belong to a different class of vector solitons which, in our notation, can be identified as $`|1,2`$ states. The extensive numerical analysis of soliton states $`|1,2`$  shows that all such solitons are linearly unstable. However, the observation of this instability is beyond the experimental parameters of Ref. .

The complex structure of multi-hump solitons and nonintegrability of the model (1) result in a variety of collision scenarios, which are quite dissimilar to the collisions of multi-hump solitons of the exactly integrable Manakov system . For instance, even linearly stable vector solitons do not necessarily survive soliton collisions. In Figs. 5(a,b) we show two examples of non-elastic interaction of linearly stable $`|0,1`$ and $`|0,2`$ solitons.

In conclusion, we have analyzed, analytically and numerically, the stability of multi-hump optical solitons in a saturable nonlinear medium. We have found that multi-hump solitons are members of an extended class of vector solitons which can be linearly stable in a wide region of their existence, although they may be destroyed in collisions. We believe that this is an important physical result that calls for a revision of our understanding of the structure and stability of many types of multi-hump solitary waves in nonintegrable multi-component models, usually omitted in the analysis because of their a priori assumed instability.

We thank M. Segev, D. Christodoulides, M. Mitchell, and A. Buryak for useful discussions. Yu.K. and E.O. are members of the Australian Photonics Cooperative Research Centre. D.S. acknowledges support from the Royal Society of Edinburgh and British Petroleum.