---
abstract: 'We investigate competition between two phase transitions of the second kind induced by the self-attractive nonlinearity, *viz.*, self-trapping of the leaky modes, and spontaneous symmetry breaking (SSB) of both fully trapped and leaky states. We use a one-dimensional mean-field model, which combines the cubic nonlinearity and a double-well-potential (DWP) structure with an elevated floor, which supports leaky modes (quasi-bound states) in the linear limit. The setting can be implemented in nonlinear optics and BEC. The order in which the SSB and self-trapping transitions take place with the growth of the nonlinearity strength depends on the height of the central barrier of the DWP: the SSB happens first if the barrier is relatively high, while self-trapping comes first if the barrier is lower. The SSB of the leaky modes is characterized by specific asymmetry of their radiation tails, which, in addition, feature a resonant dependence on the relation between the total size of the system and radiation wavelength. As a result of the SSB, the instability of symmetric modes initiates spontaneous Josephson oscillations. Collisions of freely moving solitons with the DWP structure admit trapping of an incident soliton into a state of persistent shuttle motion, due to emission of radiation. The study is carried out numerically, and basic results are explained by means of analytical considerations.'
author:
- 'Krzysztof B. Zegadlo$^{1}$, Nir Dror$^{2}$, Marek Trippenbach$^{3}$, Miroslaw A. Karpierz$^{1}$ and Boris A. Malomed$^{2}$'
title: 'Spontaneous symmetry breaking of self-trapped and leaky modes in quasi-double-well potentials'
---

Introduction
============

Usually, the ground state (GS) of quantum-mechanical systems exactly follows the symmetry of the underlying Hamiltonian [@LL], while excited states may realize different representations of the same symmetry (a different situation is exemplified by the Jahn-Teller effect in molecules, which makes the GS of the electron subsystem spatially asymmetric, thus breaking the symmetry of the respective Hamiltonian [@RE]). In particular, for the double-well potential (DWP), which is dealt with in the present work, the GS wave function is spatially even, while the first excited state is odd. This is not necessarily true in many-body settings. In that context, the mean-field description of atomic Bose-Einstein condensates (BECs) is provided by the Gross-Pitaevskii equation (GPE) [@BEC], which includes the cubic term accounting for attractive forces between colliding atoms. Essentially the same is the nonlinear Schrödinger equation (NLSE) modeling the propagation of optical signals in Kerr-nonlinear media [NLS]{}. If the self-focusing nonlinearity is strong enough, it gives rise to the phase transition in the form of *spontaneous symmetry breaking* (SSB) of the GS [@book]. In its simplest manifestation, which is provided by the DWP, the SSB implies that one well traps a larger atomic density or field power than the other. This effect also implies the breakup of the basic principle of quantum mechanics, according to which the GS cannot be degenerate, as the SSB gives rise to a pair of two mutually symmetric GSs in the DWP, with the maximum of the wave function found in either potential well (as mentioned above, the Jahn-Teller effect gives rise to a qualitatively similar situation). The same DWP setting admits a symmetric state coexisting with the asymmetric ones, but, above the SSB point, the symmetric wave function no longer represents the GS, being unstable against symmetry-breaking perturbations. In the course of the spontaneous transition from the unstable symmetric state to a stable asymmetric one, the choice between the two mutually degenerate asymmetric states is determined by random perturbations, which push the system to build the maximum of the wave function in the left or right potential well. The SSB is a ubiquitous phenomenon, with well-known manifestations in nonlinear optics, BEC, superfluidity, superconductivity, ferromagnetism, etc. [book]{}.

The concept of the SSB in nonlinear systems of the NLS type was, plausibly, introduced for the first time in 1979 by E. B. Davies [@Davies], who addressed a nonlinear extension of the Schrödinger equation for a pair of quantum particles with an isotropic interaction potential. In this context, the SSB was predicted as the breaking of the rotational symmetry in the GS. Another early work, which predicted the SSB in a relatively simple form, addressed the *self-trapping model*, based on a system of linearly coupled ordinary differential equations including the self-attractive cubic terms [@Scott].

In the effectively one-dimensional geometry, the SSB can be studied in the framework of the scaled NLSE/GPE with potential $H(x)$ of the DWP type, for the amplitude of the electromagnetic wave, or the single-particle wave function, $\psi \left( x,z\right) $: $$i\frac{\partial \psi }{\partial z}=-\frac{1}{2}\frac{\partial ^{2}\psi }{%
\partial x^{2}}-\left\vert \psi \right\vert ^{2}\psi +H(x)\psi ,
\label{NLSE}$$where $z$ is the propagation distance in optics, or time in the GPE. This equation can be reduced to a system of coupled ordinary differential equations for two amplitudes, $u_{1,2}(z)$, by means of the *tight-binding approximation* [@tight], which replaces $\psi (x,z)$ by a linear superposition of two stationary wave functions, $\phi $, corresponding to the states trapped separately in either potential well, with their centers located at $x=\pm a$ [@Ananikian]:$$\psi \left( x,z\right) =u_{1}(z)\phi \left( x-a\right) +u_{2}(z)\phi \left(
x+a\right) .  \label{12}$$

The analysis of the SSB in BEC and similar models based on Eq. (\[NLSE\]) was initiated in Refs. [@Milburn] and [@Smerzi]. In this case, in the framework of the mean-field approximation, the symmetry breaking is the phase transition of the second kind (alias the supercritical bifurcation, which does not admit hysteresis [@Joseph]). Further, GPE (\[NLSE\]) was extended by adding an extra (free) spatial coordinate, which transforms the DWP into a two-dimensional dual-core structure [@Warsaw]. In such a setting, the self-attractive nonlinearity gives rise to matter-wave solitons, which self-trap in the free direction [@soliton]. The SSB destabilizes symmetric solitons and replaces them by asymmetric ones, provided that the norm of the wave function (which determines the effective strength of the intrinsic nonlinearity) exceeds a critical value [Warsaw]{}. In the latter case, the mean-field symmetry breaking is a phase transition of the first kind (alias a subcritical bifurcation [@Joseph]), which includes hysteresis. The subcritical transition is typical to solitons in dual-core waveguides with the Kerr self-focusing [soliton-dual-core]{}. The same type of the transition may be featured by CW (continuous-wave) states in dual-core systems with non-Kerr nonlinearities [@Snyder].

In addition to the analysis of static symmetric and asymmetric modes, dynamical regimes, most typically in the form of oscillations of the mean-field wave function between two wells of the DWP structure, were analyzed too. Following the analogy with Josephson oscillations of the wave function of Cooper-paired electrons in superconducting tunnel junctions [Ustinov]{}, the possibility of oscillations in *bosonic Josephson junctions* was predicted [@junction]. The simplest dynamical model of the Josephson oscillations in bosonic systems was derived by means of the tight-binding approximation ([@superconductor]).

Experimental manifestations of the SSB have been observed in both BEC and photonics. Self-trapping of a macroscopically asymmetric state of the atomic condensate of $^{87}$Rb atoms, loaded into the DWP, as well as Josephson oscillations in that setting, were reported in Ref. [@Markus] (in that case, the effective nonlinearity is self-repulsive, therefore the respective SSB occurs not in the symmetric GS, but rather in the antisymmetric first excited state). The SSB of laser beams coupled into an effective transverse DWP created in the self-focusing photorefractive medium was demonstrated in Ref. [@photo]. Other experimentally observed SSB effects in optics are spontaneously established asymmetric regimes of operation of coupled lasers [@lasers; @NatPhot; @NaturePhot2; @Chili], and breaking of chiral symmetry in metamaterials [@Kivshar].

In addition to usual bound states, one may work with quasi-localized modes in potentials which do not admit complete trapping in linear quantum mechanics, but give rise to leaky bound states, alias quasi-bound ones. The combination of such a potential and self-attractive nonlinearity makes it possible to transform the leaky states into truly bound ones [Carr1,Carr2]{}. This possibility, in turn, suggests another setting, which is the subject of the present work: DWP structures embedded into a potential barrier. In the linear limit, this structure support solely symmetric leaky modes, that may be transformed into self-trapped ones with the help of the cubic self-attraction. The main feature of the system which, to the best of our knowledge, was not explored before, is competition between two different mean-field phase transitions of the second kind, driven by the nonlinearity: the SSB and transition to the self-trapping. Realization of the competition in stationary states of the DWP system is the main subject of the present work. We demonstrate that, depending on parameters of the DWP structure and nonlinearity strength, either transition may happen first, with the growth of the nonlinearity. Another essential problem addressed in the paper is a dynamical one, namely, Josephson oscillations in the DWP structure, initiated by the instability of the symmetric mode, and collisions of free solitons with the structure.

It is relevant to mention that, in terms of the BEC realizations, the present setting represents macroscopic quantum states, with the phase transitions between them being quasi-classical ones, considered in the framework of the mean-field approximation. The validity of this approximation is usually justified by the large number of atoms in the condensate [@BEC]. The consideration of a few-body state in the DWP can give rise to quantum phase transitions, such as those in the Lipkin-Meshkov-Glick model, which applies in this case [@LMG]. As suggested by a recent analysis of the three-dimensional many-body quantum gas with repulsive binary interactions, which is pulled to the center by a potential $\sim -1/r^{2}$ [@Grisha], the quantum phase transition may produce results similar to but different from their mean-field counterparts [@HS]. In particular, the GS predicted by the mean-field may be replaced by a metastable state in the quantum many-body theory [@Grisha]. In any case, the consideration of truly quantum phase transitions in the DWP structure is a subject for a separate work.

The subsequent presentation is structured as follows. The model is elaborated in Section II. Results of the analysis of symmetric and spontaneously emerging asymmetric trapped and leaky modes in it are summarized in Section III. Detailed results are obtained in a numerical form, and their basic features are explained by means of an analytical approach. Both the trapped and leaky modes undergo the SSB transition with the increase of the norm, the symmetric modes getting unstable above the transition point. The nonlinear evolution of the unstable modes, which features Josephson oscillations, is studied by means of systematic simulations in Section IV. A related possibility is capture of incident solitons by the DWP structure into shuttle states. This possibility is studied in a systematic form in Section V. The paper is concluded by Section VI.

The model
=========

The underlying dynamical model, based on Eq. (\[NLSE\]), gives rise to stationary modes with propagation constant $k$ (in BEC, $-k$ is the chemical potential),  $$\psi \left( x,z\right) =e^{ikz}u\left( x\right) ,  \label{psi}$$where real modal functions $u(x)$ satisfy the equation

$$-ku+\frac{1}{2}u^{\prime \prime }+u^{3}=H\left( x\right) u.
\label{stationaryNLSE}$$

Solutions with $k>0$ represent self-trapped localized states, while $k<0$ corresponds to delocalized leaky modes, which do not vanish at $x\rightarrow
\pm \infty $. The states of these two types are characterized, respectively, by convergent and divergent norms, $N=\int_{-\infty }^{+\infty
}u^{2}(x)dx$ (proportional to the total power of the light beam in optics, or the total number of atoms in BEC), and Hamiltonian (energy),$$E=\int_{-\infty }^{+\infty }\left[ \frac{1}{2}\left( u^{\prime }\right) ^{2}-%
\frac{1}{2}u^{4}+2H(x)u^{2}\right] dx.  \label{H}$$

The DWP can be readily implemented in the experiment. In optics, waveguides with this structure are fabricated with the help of the implanting technique [@Opt-DWP], while in BEC the DWP setting can be created by means of electromagnetic fields [@BEC-DWP]. In the present work, calculations are reported for the rectangular DWP profile with the elevated floor:

$$H\left( x\right) =\left\{
\begin{array}{ll}
A, & ~\mathrm{at}~~|x|<0.5, \\
2, & ~\mathrm{at}~~3<|x|<7, \\
0, & ~\mathrm{elsewhere},%
\end{array}%
\right.   \label{potential}$$

where $A>0$ is the height of the inner potential barrier, see Fig. [potential-fig]{}, while height $H_{\max }=2$ of the outer barriers is fixed by scaling. Values of lengths adopted in Eq. (\[potential\]) adequately represent the generic situation, as demonstrated by additional numerical results (not shown here in detail). Indeed, it is shown below that the symmetry-breaking and self-trapping transitions, and the competition between them, crucially depend on the tunneling transparency of the central barrier and the nonlinearity strength, i.e., the barrier height, $A$, and total norm, $N$. These are two control parameters which are subject to the variation in the subsequent analysis. For the same reason, the main findings are not sensitive to a particular shape of the DWP. In particular, the rectangular form of the DWP, adopted in Eq. (\[potential\]), which may be essential for some dynamical effects, such as temporal scaling in the relaxation of perturbations [@Azbel], produces results which are essentially the same as generated by smooth DWP profiles (for the self-trapping of the leaky modes in a single potential well, this property was known before [@Carr1]). As concerns the necessity of having the elevated potential floor, it may make the experimental creation of the structure easier, as the floor" is naturally built by overlap of fringes of two potential barriers which determine the DWP (in previously reported experimental realizations of the DWP in BEC [@Markus], the bottom level of the potential had to be depressed, because those DWPs trapped condensates with the repulsive interactions, on the contrary to the case of the self-attraction, considered herein).

![The symmetric double-well potential for leaky modes (quasi-bound states), defined as per Eq. (\[potential\]).[]{data-label="potential-fig"}](potential.pdf)

Obviously, in the absence of the self-attractive cubic term, potential ([potential]{}) cannot support any bound state in the respective linear model, while weakly delocalized quasi-bound states are possible. Indeed, a straightforward estimate of the tunneling coefficient for the tall barriers separating the inner and outer parts of structure (\[potential\]) yields $$T\simeq \exp \left( -\sqrt{2H_{\max }-q^{2}}W\right) \approx \allowbreak
4.4\times 10^{-4},  \label{T}$$where $H_{\max }=2$ and $W=4$ are the height and width of the potential barriers, as per Eq. (\[potential\]), and $q\equiv \pi /l=\pi /6$ is the wavenumber of the lowest quasi-bound state in the potential box of width $%
l=6 $.

For the study of collisions of moving solitons with the DWP structure, the height of the inner barrier is fixed to be $A=0.5$, while the height, $H_{0}$, and width, $W_{0}$, of the outer barriers will be varied, to allow clearer observations of different collision scenarios: $$H\left( x\right) =\left\{
\begin{array}{ll}
0.5, & ~\mathrm{at}~~|x|<0.5, \\
H_{0}, & ~\mathrm{at}~~3<|x|<3+W_{0}\equiv \Lambda /2, \\
0, & ~\mathrm{elsewhere}.%
\end{array}%
\right.  \label{P2}$$

Detailed consideration of the SSB in the leaky modes will require an explicit calculation of small-amplitude radiation tails attached to those modes outside of the barriers, i.e., at $|x|>7$, see Eq. (\[potential\]). For this purpose, the DWP is embedded into a broad free-space domain, $%
|x|<L/2$, with zero boundary conditions (b.c.): $$u(|x|=L/2)=0.  \label{bc}$$

Spontaneous symmetry breaking (SSB) of leaky and trapped modes
==============================================================

The structure of symmetric and asymmetric modes
-----------------------------------------------

Numerical solutions of Eq. (\[stationaryNLSE\]) were obtained by means of the shooting and Newton-matrix methods. While the former one makes it possible to find all relevant solutions independently of an input trial function, the latter method can be applied to obtain solution with high accuracy, provided that the initial guess is taken not too far from the final result. For the setting addressed in this paper, the combination of both algorithms is the most efficient way of obtaining stationary solutions.

As is typical for the SSB in systems with self-focusing nonlinearity, it was found that the GS is spatially symmetric below the bifurcation point ($k<k_{%
\mathrm{bif}}$) and asymmetric above it, at $k>k_{\mathrm{bif}}$. The symmetric state exists at $k>k_{\mathrm{bif}}$ too, but in that case it is not a GS, and is no longer stable. As mentioned above, all mean-field phase transitions exhibited by the present system are of the second kind, featuring no hysteresis or bistability between symmetric and asymmetric modes.

Generic examples of unstable symmetric and stable asymmetric states of both trapped ($k>0$) the leaky ($k_{\mathrm{bif}}<k<0$) types are displayed in Figs. \[symmetric and asymmetric modes\](a) and (b), respectively. In the latter case, the delocalized tails of the leaky mode are extremely small, with amplitude $$u_{\mathrm{rad}}^{(\max )}\approx 1.2\times 10^{-4}.  \label{rad}$$Taking into regard the value of the amplitude of the delocalized mode at its center, $U_{0}\approx 0.40$, tunneling coefficient (\[T\]) predicts amplitude $u_{\mathrm{rad}}^{(\max )}\sim TU_{0}\simeq 1.8\times 10^{-4}$, in reasonable agreement with its numerically found counterpart (\[rad\]). For given $k$, the spatial period of the tail in the free space is expected to be $\lambda =\pi \sqrt{2/|k|}\approx \allowbreak 15$ for $k=-0.085$ in Fig. \[symmetric and asymmetric modes\](b), while the numerical solutions demonstrates a close value, $\lambda \approx 13.5$ (a small deviation from the predicted value may be explained by the proximity of the tail to the outer barrier).

 

The formally diverging contribution of the tails to the total norm of the leaky mode is negligible, $N_{\mathrm{rad}}\simeq \left( L/2\right) \left[
u_{\mathrm{rad}}^{(\max )}\right] ^{2}$, where $L\approx 20$ is the total size of the free-space part of the configuration displayed in Fig. [symmetric and asymmetric modes]{}(b). Indeed, using estimate (\[rad\]) for the tail’s amplitude yields $N_{\mathrm{rad}}\sim 10^{-7}$, therefore the leaky modes have definite values of the norm, as indicated in the caption to Fig. \[symmetric and asymmetric modes\](b).

For the sake of direct comparison between leaky and trapped modes, in Fig. \[symmetric and asymmetric modes\](c) we display stationary states with the same propagation constant as in Fig. \[symmetric and asymmetric modes\](b), but in the case when they are confined by impenetrable (infinitely tall) outer potential barriers in Eq. (\[potential\]). It is seen that, similar to the situation shown in Fig. \[symmetric and asymmetric modes\](b), the smaller and larger values of $N$ correspond to the broken asymmetric and unbroken symmetric modes, respectively. However, in the infinitely deep potential box the profiles of the wave functions are, naturally, narrower and taller.

The symmetric modes displayed in Fig. \[symmetric and asymmetric modes\] feature split peaks, due to the fact that the inner potential barrier in Eq. (\[potential\]) is relatively high. On the other hand, the same potential structure with an essentially smaller barrier’s height $A$ supports single-peak symmetric modes, see Fig. \[unstable evolution\] below.

SSB of radiation tails in leaky modes
-------------------------------------

As mentioned above, the SSB of trapped modes is a known effect, which was previously studied in other forms [@book; @NaturePhot2]. A new phenomenon reported here is the SSB of leaky modes, which include nonvanishing tails extending into the free space outside of the DWP structure. Even if the tails have small amplitudes, it is interesting to analyze their structure in asymmetric modes, as this issue was not considered previously. To this end, we here focus on the symmetric and asymmetric states in the setting based on the DWP (\[potential\]) with $A=0.5$ and $\Lambda =14$, embedded into a broad domain of size $L=128$, see Eqs. (\[potential\]), (\[P2\]), and (\[bc\]). In this case, the asymmetric modes are found at $k\geq -0.100$.

A characteristic example of *asymmetric tails*, i.e., left and right ones with unequal amplitudes, is shown in Fig. \[tails\](a) for $k=-0.075$. Further, separately calculated total norms of the right and left tails, in regions $\Lambda /2<x<L/2$ and $-\Lambda /2<x<-L/2$, along with the norm of the tails in the coexisting unstable symmetric leaky mode, are displayed, as functions of $k$, in Fig. \[tails\](b). This dependence exhibits two notable features. First, the asymmetry between the right and left tails emerges at $k=-0.100$ and gradually increases with the increase of $k$ (i.e., decrease of $|k|$), even if each tail’s norm vanishes in the limit of $k$$\rightarrow 0$ (when the transition to the self-trapped mode takes place, and the tails vanish). This feature is illustrated by Fig. \[tails\](c), which displays the asymmetry measure vs. $k$:$$\theta (k)\equiv \frac{\int_{\Lambda /2}^{L/2}u^{2}(x;k)dx-\int_{-\Lambda
/2}^{-L/2}u^{2}(x;k)dx}{\int_{\Lambda /2}^{L/2}u^{2}(x;k)dx+\int_{-\Lambda
/2}^{-L/2}u^{2}(x;k)dx}~.  \label{theta_tail}$$Second, the dependence of the tails’ norms on $k$ shows strong oscillations, which is explained by the commensurability-incommensurability transitions between the wavelength of the radiation tail and the total size of the free-space domains, $L/2-\Lambda /2$. Indeed, the radiation wavenumber given by the free-space dispersion relation for linearized equation ([stationaryNLSE]{}), $q=\sqrt{-2k}$, determines the the radiation half-wavelength, $\pi /q$, which, in the case of the commensurability, satisfies relation $(\pi /q)n=L/2-\Lambda /2$, with $n=1,2,3,...$ . Thus, maxima of the radiation amplitude are expected at discrete values of the propagation constant,  $$k_{n}=-2\left[ \pi n/\left( L-\Lambda \right) \right] ^{2}.  \label{k}$$As shown in Fig. \[tails\](b), Eq. (\[k\]) quite accurately predicts positions of the tail-norm peaks for $n=2,3,4,5,6,7,$ and $8$, for $%
L-\Lambda =114$, which corresponds to the present case \[$n=1$ yields $%
k_{1}\approx -1.5\times 10^{-3}$, for which the tail’s amplitudes are too small to discern the corresponding maximum, while Eq. (\[k\]) with $n=9$ predicts $k_{9}\approx -0.123$, for which asymmetric modes do not exist).

SSB diagrams
------------

Getting back to the consideration of the SSB for the entire system, systematic results are presented by means of plots $N(k)$ and $E(N)$ \[the Hamiltonian is defined by Eq. (\[H\])\] for symmetric and asymmetric modes, which are displayed in Figs. \[bifurcation diagrams\](a-c) and (d-f) for three different values of height $A$ of the inner barrier of potential structure (\[potential\]). Note that the $N(k)$ curves obey the Vakhitov-Kolokolov criterion [@VK; @VK2], which is necessary but not sufficient for the stability of modes supported by the self-attractive nonlinearity (it does not catch the instability of the symmetric modes coexisting with asymmetric ones). It is relevant to compare these plots with their counterparts,$$N_{\mathrm{sol}}=2\sqrt{2k},~E_{\mathrm{sol}}=-(1/3)\left( 2k\right) ^{3/2}%
\text{,}  \label{sol}$$for the NLS solitons in the free space, given by Eq. (\[psi\]) with$$u_{\mathrm{sol}}(x)=\sqrt{2k}\,\mathrm{sech}\left( \sqrt{2k}x\right) ,
\label{soliton}$$which are displayed by dashed curves in Figs. \[bifurcation diagrams\](a-f).

The SSB in the families of trapped and leaky states is characterized by the asymmetry ratio,

$$\Theta \equiv N^{-1}\left[ \int_{0}^{+\infty }u^{2}(x)dx-\int_{-\infty
}^{0}u^{2}(x)dx\right]  \label{theta}$$

\[cf. a similar definition for the tails, given by Eq. (\[theta\_tail\])\], which is shown as a function of $N$ in Figs. \[bifurcation diagrams\](g-i). These plots clearly identify the SSB-onset points, at which the symmetric mode gets destabilized, and simultaneously a stable asymmetric one emerges. In accordance with what is said above, the bifurcation is of the supercritical type [@Joseph], i.e., the emerging branches of the asymmetric states immediately go forward". Conclusions about the stability and instability of the solution branches displayed in Fig. \[bifurcation diagrams\] were produced by means of the well-known method [@VK2] based on numerical computation of (in)stability eigenvalues (imaginary parts of eigenfrequencies) for small perturbations added to the stationary solutions, using the respective linearized (Bogoliubov - de Gennes) equations. In particular, the instability of those symmetric states which coexist with asymmetric ones is always represented by a single pair of purely imaginary eigenfrequencies.

The fact that the asymmetric modes, when they exists, have smaller values of $E$ for given $N$, hence they realize the system’s GS \[see Figs. [bifurcation diagrams]{}(d-f)\], can be easily understood: having their center shifted from the layer occupied by the positive potential ($|x|<0.5$) to a region where $H=0$ \[$0.5<|x|<3$, see Figs. \[potential-fig\] and [symmetric and asymmetric modes]{}\], they obviously reduce the integral value of $E$. The same argument explains why, for given $k$, the asymmetric modes feature smaller $N$: having smaller overlap with the region of $H>0$, they need a smaller norm to compensate the shift of $k$ towards negative values induced by the positive potential, see Eq. (\[stationaryNLSE\]).

Note that the families of states displayed in Figs. \[bifurcation diagrams\](a-c) comprise both $k<0$ and $k>0$, i.e., the leaky and trapped modes alike. In particular, the SSB bifurcation occurs at $k>0$ in panels \[bifurcation diagrams\](a,b), and at $k<0$ in (c) (in the latter case, the SSB sets in at $k=-0.100$, as shown for the same system in Fig. [tails]{}). A noteworthy fact is that the system features the competition of the two different mean-field phase transitions driven by the increase of $N$, i.e., the strength of the self-attraction: the transition from the quasi-bound (leaky) state to the self-trapped one, which was previously found in single-well elevated potentials [@Carr1; @Carr2], and the SSB in the DWP structure. Thus, in the cases shown in Figs. \[bifurcation diagrams\](a,b) the self-trapping transition happens first (at smaller $N$), while in Fig. \[bifurcation diagrams\](c) the SSB takes place prior to the onset of the self-trapping.

The values of the propagation constant and norm at the SSB bifurcation point are shown, as a function of the inner-barrier’s height $A$ \[see Eq. ([potential]{})\], in Fig. \[bifurcation point\]. In panel (\[bifurcation\_k\]), the boundary between the SSB occurring with the delocalized and trapped modes ($k_{\mathrm{bif}}<0$ and $k_{\mathrm{bif}}>0$, respectively) is located at $A\approx 0.30$, the same value corresponding to the boundary designated by the square symbol in panel \[bifurcation\_N\]). That is, the SSB happens first (at smaller $N$) at $A>0.30$, while the transition to the self-trapping precedes the SSB at $A<0.30$.

The fact that the SSB happens with the trapped and leaky modes, respectively, at small and large $A$, as seen in Fig. \[bifurcation point\], is easy to explain: small $A$ implies strong linear coupling between the wave functions in the two barely separated wells, hence very large $N$ is required to induce the SSB, being far greater than the value of $N$ needed for the onset of the self-trapping, which is determined by the fixed parameters of the outer barrier in the DWP structure (\[potential\]). On the contrary, large $A$ implies weak linear coupling between the strongly separated wells, hence the respective strength of the nonlinearity (measured by $N$), required for the SSB, is much smaller than the value necessary for the commencement of the self-trapping. These arguments clearly suggest that the same sequence of the SSB and self-trapping phase transitions should take place in generic DWP structures.

These arguments can be cast in a more definite form, if the central barrier in Eq. (\[potential\]) is approximated by $H_{\mathrm{central}}(x)=A\delta
(x)$, and the outer barriers are made impenetrable, similar to what is shown in Fig. \[symmetric and asymmetric modes\](c). These conjectures replace the present model by the one for an infinitely deep potential box split by delta-functional barrier, which is the simplest model of the SSB [NaturePhot2,Chili]{}. In the limit of large $A$, the latter model predicts the following value of the norm at the SSB bifurcation point,$$N_{\mathrm{bif}}\approx 8\pi ^{2}/\left( 3l^{2}A\right)   \label{analyt}$$(see caption to Fig. 2 in Ref. [@NaturePhot2]), where $l/2$ is the coordinate of the point at which the wave function vanishes (half-width of the infinitely deep box). In particular, one can identify $l\simeq 9$ for $%
A=0.5$ in Fig. \[symmetric and asymmetric modes\](b), hence Eq. ([analyt]{}) yields $N_{\mathrm{bif}}~\simeq 0.65$, while the respective numerical  value in Fig. \[bifurcation point\](b) is $N_{\mathrm{bif}%
}~\approx 0.62$, which implies a reasonable agreement for the present (not really large) value of $A$.

Evolution of unstable symmetric states
======================================

Conclusions concerning the stability and instability of the symmetric and asymmetric modes, presented in Fig. \[bifurcation diagrams\], were verified, in addition to the computation of eigenvalues for small perturbations, by direct simulations, performed by dint of the finite-difference algorithm. The instability development of unstable symmetric states was catalyzed by adding small initial symmetry-breaking perturbations to them. This was done for the unstable symmetric states with both $k>0$ and for $k<0$, i.e., self-trapped and leaky ones. Typical examples are displayed in Fig. \[unstable evolution\], in which the unstable symmetric modes feature single- and double (split)-peak shapes at small and large values of the inner-barrier’s height, $A=0.05$ and $A=0.5$, respectively (in the top and bottom rows of the figure).

As mentioned above, the instability of symmetric states is accounted for by pure imaginary eigenfrequencies, hence the originally developing instability is not oscillatory. However, the nonlinearity makes the unstable dynamics oscillatory, as seen in Fig. \[unstable evolution\]. In other words, the unstable symmetric modes spontaneously develop bosonic Josephson oscillations. Close to the instability onset, the effective oscillation period is very large (as it diverges precisely at the onset point), gradually decreasing deeper into the instability region. The dynamical symmetry breaking induced by the weak and moderate instability is incomplete, leading to periodic oscillations between the original symmetric state and a new asymmetric one, as observed in Fig. \[unstable evolution\](a-e). Stronger instability causes complete symmetry breaking, replacing the symmetric state by an irregularly vibrating asymmetric mode, as seen in Fig. \[unstable evolution\](f).

It is relevant to note the difference between the oscillatory regimes generated by the instability of self-trapped and leaky modes. Indeed, while Fig. \[unstable evolution\](c) demonstrates that the shape of the oscillating mode is sharp in the former case, the shape is fuzzy in Fig. \[unstable evolution\](e) because it involves an essential radiation component in the case when the underlying unstable mode is a leaky one.

Collisions of free solitons with the quasi-double-well potential structure {#sec:CollisionResults}
==========================================================================

In addition to the analysis of the stationary states performed above, it is also relevant to consider collisions of free solitons with the elevated DWP structure. For this purpose, the initial soliton was created as the tilted (moving) version of the static one given by Eqs. (\[psi\]) and ([soliton]{}),$$\psi \left( x,z\right) =\sqrt{2k}\exp \left( i\left( k-c^{2}/2\right)
z+icx\right) \text{\textrm{sech}}\left( \sqrt{2k}\left( \left(
x-x_{0}\right) -cz\right) \right) ,  \label{tilt}$$where $c$ is a real tilt (velocity), and $x_{0}$ is the initial position of the soliton, chosen far enough from the localized potential structure (we here take $x_{0}=-15$). Generic findings are produced here for incident solitons with $k=3,$ the corresponding norm being $$N=2\sqrt{2k}\approx 4.90,  \label{N}$$according to Eq. (\[sol\]), other values of $N$ giving similar results.

Figure \[ThetaVsHeightMapping\] presents a parameter chart for three different outcomes of the collisions, produced by varying the height of the outer barriers, $H_{0}$ in Eq. (\[P2\]), and tilt $c$ in Eq. (\[tilt\]). In the region designated as (1), i.e., for $c$ small and/or $H_{0}$ large enough, the incident soliton bounces back the left outer barrier, as shown in Fig. \[FundCollision1\] for $H_{0}=0.6$ and $c=0.8$. At larger $c$, in a relatively narrow region (2) of Fig. \[ThetaVsHeightMapping\], the soliton gets captured inside the potential structure, which resembles the previously explored possibility of capturing an incident Bragg-grating soliton by a cavity formed by two locally repulsive defects [@cavity]. Similar to that setting, the trapped soliton performs shuttle motion between the outer and inner potential barriers, as shown in Fig. \[FundCollision2\], for $H_{0}=0.6$ and $c=1.15$. The shuttle dynamical regime is an essential addition to the stationary states revealed by the analysis presented above. At still larger $c$, but yet staying within the boundaries of the narrow region (2), the shuttle motion of the trapped soliton becomes irregular, see Fig. \[FundCollision5\] for $H_{0}=1.3$ and $c=1.8$.

![The chart of outcomes of collisions of free solitons, launched with tilt $c$ \[see Eq. (\[tilt\])\], with potential structure (\[P2\]), in the plane of $\left( H_{0},c\right) $. The norm of the incident soliton is fixed as per Eq. (\[N\]). In this figure, the width of the outer potential barriers is $W_{0}=4$. In region (1), the soliton bounces back from the left outer barrier. Region (2) refers to trapping the soliton inside the potential structure, where it performs shuttle motion. In region (3), the incident soliton passes the structure. The dashed line is the analytical prediction produced by Eq. ([c\_cr]{}).[]{data-label="ThetaVsHeightMapping"}](ThetaVsHeightMapping.pdf){width="3.2in"}

In region (3), the initial tilt, $c$, is large enough for the soliton to pass the potential structure, see an example in Fig. \[FundCollision3\] for $H_{0}=0.6$ and $c=1.5$. Another option admitted by this scenario is shown in Fig. \[FundCollision4\], where the soliton passes the left outer barrier, bounces back from the inner one, and escapes in the reverse direction, passing the left outer barrier again. Naturally, this collision pattern is common for lower $H_{0}$, when the outer barriers are *lower* than the inner one, which plays the role of a strong “bouncer". Furthermore, at $c>1.85$, the incident soliton splits into two fragments, one escaping and the other one staying in a chaotically evolving trapped state, see an example in Fig. \[FundCollision6\], for $H_{0}=1.2$ and $%
c=1.92$.

Collision scenarios were also explored by varying width $W_{0}$ of the outer barriers, while keeping their height constant, $H_{0}=1$, as well as characteristics of the inner barrier and the distance between the barriers, see Eq. (\[P2\]). The respective results, for the same incident soliton as used above ($k=3$, $N=4.90$), are summarized in Fig. [ThetaVsWidthMapping]{}, where regions (1)-(3) have the same meaning as their counterparts in Fig. \[ThetaVsHeightMapping\].

In the latter case, the results may be classified into three outcomes of the collision, depending on $W_{0}$. The first outcome occurs at $0<W_{0}<0.8$. It is characterized by a rapidly growing region of the shuttle motion of the trapped soliton \[region (2)\], and *sharp transitions* between the three evolution scenarios, $(1)\leftrightarrow (2)$, $(2)\leftrightarrow (3)$. The second outcome, which was observed in the range of $0.8<W_{0}<1.8$ (in this case, $W_{0}$ is, roughly, close to the width of the incident soliton), is distinguished by *gradual transitions* between the scenarios. That means that, for certain values of $c$, the soliton does not fully bounce from the barrier, nor penetrates it, but rather splits into two segments, one of which escapes, while the other remains trapped. In this case, both the upper and lower boundaries of region (2) represent tilts at which the soliton is split into equal fragments. An example of such an outcome is shown in Fig. \[SplitSolitonCollision\], for $W_{0}=1.2$ and $c=1.66$. The third outcome, which is observed at $W_{0}>1.8$, is distinguished by the fact that the variation of $W_{0}$ almost does not affect the soliton’s motion. In contrast to the second outcome, and similar to the first one, the respective transitions between the regions are sharp.

![The same as in Fig. \[ThetaVsHeightMapping\], but varying the width of the outer potential barriers, $W_{0}$, while their height is fixed, $H_{0}=1$. The dashed curve shows the analytical prediction for the boundary of area (1) given by Eq. (\[parabola\]). []{data-label="ThetaVsWidthMapping"}](ThetaVsWidthMapping.pdf){width="3.2in"}

![Outcomes of the collisions, for $W_{0}=1.2$ and $c=1.66$. In this borderline example, the soliton hits the right (second) outer potential barrier and splits into two fragments, one escaping and the other one staying trapped in the potential structure.[]{data-label="SplitSolitonCollision"}](fundcollisionmum3wo1p2theta1p66.pdf){width="3.2in"}

It is easy to explain the parabolic boundary of region (1) in Fig. [ThetaVsHeightMapping]{}, as well as the boundary of the same region in Fig. \[ThetaVsWidthMapping\], using the perturbation theory for NLS solitons, which treats them as quasi-particles with mass $N$ [@RMP]. Indeed, the kinetic energy of the soliton is $E_{\mathrm{kin}}=(N/2)c^{2}$, while the height of the outer potential barrier for the quasi-particle, $E_{0}$, can be obtained from the third term in expression (\[H\]), assuming that the soliton’s center is located at the midpoint of the potential barrier:$$E_{0}=H_{0}\int_{-W_{0}/2}^{+W_{0}/2}u_{\mathrm{sol}}^{2}\left( x^{\prime
}\right) dx^{\prime }=H_{0}N\tanh \left( \frac{1}{4}W_{0}N\right) ,
\label{E0}$$where $u_{\mathrm{sol}}$ is the solitons’s profile (\[soliton\]), $%
x^{\prime }\equiv x-\left( 3+W_{0}/2\right) $ \[see Eq. (\[P2\])\], and the result is expressed in terms of the soliton’s norm, as per Eq. (\[sol\]). Next, equating $E_{\mathrm{kin}}$ to $E_{0}$ predicts that the boundary between the rebound and passage corresponds to$$c_{\mathrm{cr}}=\sqrt{2H_{0}\tanh \left( \frac{1}{4}W_{0}N\right) }.
\label{parabola}$$In particular, the respective prediction for dependence $c_{\mathrm{cr}%
}(H_{0})$ corresponding to the case displayed in Fig. [ThetaVsHeightMapping]{}, with $W_{0}=4$ and $N$ fixed as per Eq. (\[N\]), simplifies to$$c_{\mathrm{cr}}=\sqrt{2H_{0}}.  \label{c_cr}$$Figure \[ThetaVsHeightMapping\] demonstrates that Eq. (\[c\_cr\]) predicts the parabolic boundary of region (1) quite accurately, a discrepancy at large $H_{0}$ being explained by the fact that the collision with the tall barrier causes a deformation of the soliton’s shape. Further, the full analytical expressions (\[parabola\]) predicts the boundary of the same region in Fig. \[ThetaVsWidthMapping\] well enough too.

As concerns the trapping regime in area (2) of Figs. [ThetaVsHeightMapping]{} and \[ThetaVsWidthMapping\], it is explained by the fact that, while slowly passing the left barrier, and then passing the inner one, the soliton with the initial tilt slightly exceeding value ([parabola]{}) suffers radiation losses due to its deceleration and acceleration. The losses cause a drop in the kinetic energy below the level necessary for clearing the right barrier. A feature which relates the trapping and splitting to the leaky modes considered above is tunneling of the radiation across the outer potential barriers, which can be seen, in particular, in Figs. \[FundCollisions\](b,e,f).

Conclusions
===========

We have extended the known possibility of the stabilization of leaky modes in quasi-trapping potentials by means of the self-attractive nonlinearity. Unlike the previously studied single-well potential, we have introduced the spatially symmetric DWP (double-well potential) with the elevated floor, embedded into the potential barrier. The setting can be implemented in nonlinear optical waveguides and BEC. The new possibility offered by this system is the competition of two phase transitions of the second kind, described in the mean-field approximation: the onset of the self-trapping of the leaky modes, and the SSB (spontaneous symmetry breaking) of both true bound states and leaky modes, under the action of the self-attractive nonlinearity. With the increase of the norm of the wave field (which determines the nonlinearity strength), the former or latter transition happens first if the central barrier of the DWP structure is, respectively, low or tall. These conclusions are generic, as they do not depend on details of the DWP structure. New effects are revealed by the consideration of the SSB of the leaky modes: asymmetry of radiation tails, which are parts of these modes, and the commensurability-incommensurability interplay between the radiation wavelength and the total size of the system, into which the DWP is embedded. Systematic results have been produced in the numerical form, and their basic features were explained with the help of analytical considerations. The simulations demonstrate that unstable symmetric modes initiate Josephson oscillations. Collisions of freely moving solitons with the DWP structure were studied in a systematic form too, revealing various generic outcomes of the collisions, boundaries between which were explained in the analytical form. In particular, in addition to the stationary states with the unbroken and broken symmetry, the collisions reveal the dynamical mode, in the form of a soliton which performs persistent shuttle motion in the DWP structure.

Relevant possibilities for the extension of the analysis reported in this work may be offered by two-component systems, as well as by a two-dimensional generalization of the present setting. On the other hand, the results produced by the competition of the mean-field phase transitions suggest that it may be interesting too to consider quantum phase transitions in a many-body bosonic state with attractive inter-particle interactions, loaded into the DWP.

Acknowledgment
==============

M.T. and B.A.M. acknowledge partial support from the National Science Center of Poland in the framework of the HARMONIA Program, No. 2012/06/M/ST2/00479. K.B.Z. acknowledges support from the National Science Center of Poland through Project ETIUDA No. 2013/08/T/ST2/00627.

[99]{} D. Landau and E. M. Lifshitz, *Quantum Mechanics* (Nauka Publishers: Mosocw, 1974).

R. Englman, *The Jahn-Teller Effect in Molecules and Crystals* (Wiley-Interscience: London, 1972).

S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. **80**, 1215 (2008); H. T. C. Stoof, K. B. Gubbels, and D. B. M. Dickrsheid, *Ultracold Quantum Fields* (Springer: Dordrecht, 2009).

Y. S. Kivshar and G. P. Agrawal, *Optical Solitons: From Fibers to Photonic Crystals* (Academic Press: San Diego, 2003).

*Spontaneous Symmetry Breaking, Self-Trapping, and Josephson Oscillations*, Editor: B. A. Malomed (Springer-Verlag: Berlin and Heidelberg, 2013).

E. B. Davies, “Symmetry breaking in a non-linear Schrödinger equation", Commun. Math. Phys. **64**, 191-210 (1979).

J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott, “The discrete self-trapping equation", Physica D **16**, 318-338 (1985).

G. L. Alfimov, P. G. Kevrekidis, V. V. Konotop, and M. Salerno, “Wannier functions analysis of the nonlinear Schrödinger equation with a periodic potential", Phys. Rev. E **66**, 046608 (2002).

D. Ananikian and T. Bergeman, “Gross-Pitaevskii equation for Bose particles in a double-well potential: Two-mode models and beyond", Phys. Rev. A **73**, 013604 (2006).

G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, “Quantum dynamics of an atomic Bose-Einstein condensate in a double-well potential", Phys. Rev. A **55**, 4318-4324 (1997).

A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, “Quantum coherent atomic tunneling between two trapped Bose-Einstein condensates", Phys. Rev. Lett. **79**, 4950-4953 (1997).

G. Iooss and D. D. Joseph, *Elementary Stability and Bifurcation Theory* (Springer: New York, 1980).

M. Matuszewski, B. A. Malomed, and M. Trippenbach, “Spontaneous symmetry breaking of solitons trapped in a double-channel potential", Phys. Rev. A **75**, 063621 (2007).

V. M. Pérez-García, H. Michinel, and H. Herrero, “Bose-Einstein solitons in highly asymmetric traps", Phys. Rev. A **57**, 3837-3842 (1998).

C. Paré and M. Florjańczyk, “Approximate model of soliton dynamics in all-optical couplers", Phys. Rev. A **41**, 6287-6295 (1990); A. I. Maimistov, “Propagation of a light pulse in nonlinear tunnel-coupled optical waveguides", Sov. J. Quantum Electron. **21**, 687-690 (1991); P. L. Chu, B. A. Malomed, and G. D. Peng, “Soliton switching and propagation in nonlinear fiber couplers: analytical results", J. Opt. Soc. Am. B **10**, 1379-1385 (1993); N. Akhmediev and A. Ankiewicz, “Novel soliton states and bifurcation phenomena in nonlinear fiber couplers", Phys. Rev. Lett. **70**, 2395-2398 (1993); B. A. Malomed, I. M. Skinner, P. L. Chu, and G. D. Peng, “Symmetric and asymmetric solitons in twin-core nonlinear optical fibers", Phys. Rev. E **53**, 4084-4091 (1996).

A. W. Snyder, D. J. Mitchell, L. Poladian, D. R. Rowland, and Y. Chen, “Physics of nonlinear fiber couplers", J. Opt. Soc. Am. B **8**, 2102-2118 (1991).

G. Schön and A. D. Zaikin, “Quantum coherent effects, phase transitions, and the dissipative dynamics of ultra small tunnel junctions", Phys. Rep. **198**, 237-412 (1990).

A. V. Ustinov, “Solitons in Josephson junctions", Physica D **123**, 315-329 (1998).

S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, “Coherent oscillations between two weakly coupled Bose-Einstein condensates: Josephson effects, $\pi $ oscillations, and macroscopic quantum self-trapping", Phys. Rev. A **59**, 620-633 (1999); A. Smerzi, and S. Raghavan, Phys. Rev. A **61**, 063601 (2000); R. Burioni, D. Cassi, I. Meccoli, M. Rasetti, S. Regina, P. Sodano, and A. Vezzani, “Bose-Einstein condensation in inhomogeneous Josephson arrays", Europhys. Lett. **52**, 251-256 (2000); Gati and M. K. Oberthaler, “A bosonic Josephson junction", J. Phys. B At. Mol. Opt. Phys. **40**, R61-R89 (2007).

M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, “Direct observation of tunneling and nonlinear self-trapping in a single bosonic Josephson junction", Phys. Rev. Lett. **95**, 010402 (2005).

P. G. Kevrekidis, Z. Chen, B. A. Malomed, D. J. Frantzeskakis, and M. I. Weinstein, “Spontaneous symmetry breaking in photonic lattices: Theory and experiment", Phys. Lett. A **340**, 275-280 (2005).

T. Heil, I. Fischer, W. Elsässer, J. Mulet, and C. R. Mirasso, “Chaos synchronization and spontaneous symmetry-breaking in symmetrically delay-coupled semiconductor lasers". Phys. Rev. Lett. **86**, 795-798 (2001).

P. Hamel, S. Haddadi, F. Raineri, P. Monnier, G. Beaudoin, I. Sagnes, A. Levenson, and A. M. Yacomotti, “Spontaneous mirror-symmetry breaking in a photonic molecule", Nature Phot. **9**, 311-315 (2015).

B. A. Malomed, Symmetry breaking in laser cavities", Nature Phot. **9**, 287-289 (2015).

B. A. Malomed, Spontaneous symmetry breaking in nonlinear systems: An overview and a simple model", in: *Nonlinear Dynamics: Materials, Theory and Experiments* (Springer Proceedings in Physics, vol. 173, 2016), ed. by M. Tlidi and M. Clerc, pp. 97-112.

M. Liu, D. A. Powell1, I. V. Shadrivov, M. Lapine, and Y. S. Kivshar, “Spontaneous chiral symmetry breaking in metamaterials", Nature Comm. **5**, 4441 (2014).

N. Moiseyev, L. D. Carr, B. A. Malomed, and Y. B. Band, “Transition from resonances to bound-states in nonlinear systems: application to Bose-Einstein condensates", J. Phys. B **37**, L193 (2004).

L. D. Carr, M. J. Holland, and B. A. Malomed, Macroscopic quantum tunneling of Bose-Einstein condensates in a finite potential well", J. Phys. B **38**, 3217-3231 (2005).

H. J. Lipkin, N. Meshkov, and A. J. Glick, Validity of many-body approximation methods for a solvable model. I. Exact solutions and perturbation theory“, NucI. Phys. **62**, 188-198 (1965); V. V. Ulyanov and O. B. Zaslavskii, New methods in the theory of quantum spin systems”, Phys. Rep. **216**, 179-251 (1992).

G. E. Astrakharchik and B. A. Malomed, Quantum versus mean-field collapse in a many-body system", Phys. Rev. A **92**, 043632 (2015).

H. Sakaguchi and B. A. Malomed, Suppression of the quantum-mechanical collapse by repulsive interactions in a quantum gas", Phys. Rev. A **83**, 013607 (2011).

D. Kip, “Photorefractive waveguides in oxide crystals: fabrication, properties, and applications", Appl. Phys. B - Lasers and Optics **67**, 131-150 (1998); F. Chen, X.-L. Wang, and K.-M. Wang, “Development of ion-implanted optical waveguides in optical materials: A review", Opt. Materials **29**, 1523-1542 (2007).

T. Schumm, S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer, and P. Kruger, “Matter-wave interferometry in a double well on an atom chip", Nature Phys. **1**, 57-62 (2005); S. Hofferberth, I. Lesanovsky, B. Fischer, J. Vedu, and J. Schmiedmayer, “Radiofrequency-dressed-state potentials for neutral atoms", *ibid*. **2**, 710-716 (2006).

M. Ya. Azbel and B. A. Malomed, “Possibility for time-scale-invariant relaxation in tunneling", Phys. Rev. Lett. **71**, 1617-1620 (1993).

M. Vakhitov and A. Kolokolov, “Stationary solutions of the wave equation in a medium with nonlinearity saturation", Radiophys. Quantum Electron. **16**, 783 (1973).

L. Bergé, “Wave collapse in physics: principles and applications to light and plasma waves", Phys. Rep. **303**, 259-370 (1998); G. Fibich and G. Papanicolaou, “Self-focusing in the perturbed and unperturbed nonlinear Schrödinger equation in critical dimension",SIAM J. Appl. Math. **60**, 183-240 (1999); E. A. Kuznetsov and F. Dias, “Bifurcations of solitons and their stability",Phys. Rep. **507**, 43-105 (2011).

P. Y. P. Chen, B. A. Malomed, and P. L. Chu, “Trapping Bragg solitons by a pair of defects", Phys. Rev. E **71**, 066601 (2005).

Y. S. Kivshar and B. A. Malomed, “Dynamics of solitons in nearly integrable systems", Rev. Mod. Phys. **89**, 763-915 (1989).
