One-dimensional Townes solitons in dual-core systems with localized coupling

Shatrughna Kumar¹, Pengfei Li^2,3, and Boris A. Malomed^1,4 ¹Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering, and Center for Light-Matter Interaction, Tel Aviv University, P.O.B. 39040, Tel Aviv, Israel
²Department of Physics, Taiyuan Normal University, Jinzhong, 030619, China
³Institute of Computational and Applied Physics, Taiyuan Normal University, Jinzhong, 030619, China
⁴Instituto de Alta Investigación, Universidad de Tarapacá, Casilla 7D, Arica, Chile

Abstract

The recent creation of Townes solitons (TSs) in binary Bose-Einstein condensates [10, 11] and experimental demonstration of spontaneous symmetry breaking (SSB) in solitons propagating in dual-core optical fibers [26] draw renewed interest to the TS and SSB phenomenology in these and other settings. In particular, stabilization of TSs, which are always unstable in free space, is a relevant problem with various ramifications. We introduce a system which admits exact solutions for both TSs and SSB of solitons. It is based on a dual-core waveguide with quintic self-focusing and fused (localized) coupling between the cores. The respective system of coupled nonlinear Schrödinger equations gives rise to exact solutions for full families of symmetric solitons and asymmetric ones, which are produced by the supercritical SSB bifurcation (i.e., the symmetry-breaking phase transition of the second kind). Stability boundaries of asymmetric solitons are identified by dint of numerical methods. Unstable solitons spontaneously transform into robust moderately asymmetric breathers or strongly asymmetric states with small intrinsic oscillations. The setup can be used in the design of photonic devices operating in coupling and switching regimes.

I Introduction

Critical collapse is a fundamental phenomenon with well-known manifestations in optics, plasmas, Bose-Einstein condensates (BECs), etc. [1, 2, 3, 4, 5]. Inherently related to it is the notion of Townes solitons (TSs), which were first predicted (before the concept of solitons was introduced by Zabuski and Kruskal, and the term was coined by them [6]) as two-dimensional (2D) light beams propagating in an optical medium with the cubic self-focusing [7]. TSs play the role of boundaries (separatrices) between collapsing and decaying dynamical states.

These concepts are adequately represented by the ubiquitous nonlinear Schrödinger (NLS) equations for the wave amplitude $\psi(\mathbf{r},z)$ ,

i\psi_{z}=-\frac{1}{2}\nabla^{2}\psi-|\psi|^{n-1}\psi,

(1)

where, in terms of optics, $z$ and $\mathbf{r}$ are the propagation distance and set of transverse coordinates in the $D$ -dimensional space, and $n$ is the power of the self-focusing nonlinearity. This equation can be derived from the Hamiltonian,

H=\int\left(\frac{1}{2}\left|\nabla\psi\right|^{2}-\frac{2}{n+1}|\psi|^{n+1}\right)d^{D}\mathbf{r.}

(2)

To analyze the possibility of the onset of the collapse governed by Eq. (1), one can consider a localized state of size $L$ with amplitude $A$ . The conservation of the total power (norm), $P=\int\left|\psi(\mathbf{r})\right|^{2}d^{D}\mathbf{r}$ , imposes a relation on the scaling of $A$ and $L$ in the course the development of the collapse, which corresponds to the evolution towards $L\rightarrow 0$ : $A\sim\sqrt{P}L^{-D/2}$ . Taking this relation into account, a straightforward estimate shows that the gradient and self-attraction terms in the Hamiltonian scale, at $L\rightarrow 0$ , as

H_{\mathrm{grad}}\sim+PL^{-2},H_{\mathrm{self-attr}}\sim-P^{(n+1)/2}L^{-(D/2)(n-1)},

(3)

The comparison between these terms leads to the conclusion that the negative term dominates and leads to the supercritical collapse at $D(n-1)>4$ , starting from the input with arbitrarily small $P$ . In the opposite case, $D(n-1)<4$ , the positive term is the dominant one in Eq. (3), thus preventing the collapse and making it possible to construct stable solitons. At the border, i.e., at

D(n-1)=4,

(4)

the critical collapse, leading to emergence of a singularity in the solution, takes place after a finite propagation distance. In this case, the negative term dominates in Eq. (3), which happens if $P$ exceeds a critical value, $P_{\mathrm{cr}}$ , while at $P<P_{\mathrm{cr}}$ the positive term is dominating in Eq. (3), leading to decay of the input. In the case of $D=2$ and $n=3$ , the occurrence of the singularity at $P>P_{\mathrm{cr}}$ is a rigorous corollary of the virial theorem [8, 1, 4].

Stationary solutions to Eq. (1), $\psi(\mathbf{r},z)=\exp\left(ikz\right)u(\mathbf{r})$ , with propagation constant $k$ , are invariant with respect to the conformal transformation,

\mathbf{r}\rightarrow\rho\mathbf{r},\psi\rightarrow\rho^{-2/(n-1)}\psi,k\rightarrow\rho^{-2}k,P\rightarrow\rho^{D-4/(n-1)}P,

(5)

where $\rho$ is an arbitrary scaling factor. In the critical case (4), Eq. (5) demonstrates that the soliton family, which is one of the TS type, is degenerate, as its power takes a single value that does not depend on $k$ (e.g., $P_{\mathrm{TS}}\approx 5.85$ for $D=2$ and $n=3$ , its counterpart predicted by the variational approximation being $P_{\mathrm{TS}}=2\pi$ [9]).

As mentioned above, TSs play the role of separatrices between the collapsing and decaying solutions in the case when critical condition (4) holds. For this reason, TSs are unstable against small perturbations, although the initial growth of small perturbations is subexponential, accounted for by a pair of zero eigenvalues in the respective spectrum [1, 2, 3, 4]. The latter circumstance suggests that the instability grows slowly at the initial stage, making it recently possible to experimentally observe 2D TSs in a binary BEC [10, 11], which is modeled by the Gross-Pitaevskii equation in the form of Eq. (1) with $n=3$ . This first experimental demonstration of physical states of the TS type, achieved 57 years after they had been predicted [7], makes further studies of TSs in other physical settings a relevant objective. In particular, a remaining challenge is to elaborate possibilities for observing such modes in optics, where they were predicted for the first time.

In the 1D case, the critical collapse and TS family correspond to $n=5$ in Eq. (4), i.e., the quintic self-focusing [12, 13]:

i\psi_{z}=-\frac{1}{2}\psi_{xx}-|\psi|^{4}\psi,

(6)

where $x$ is the transverse coordinate in the planar waveguide. The quintic nonlinearity occurs in various optical materials, often in a combination with a cubic term. The effective coefficients before quintic and cubic terms can be accurately adjusted to target values, including the case of the quintic-only nonlinearity (nullifying the cubic coefficient), for the light propagation through colloidal suspensions of metallic nanoparticles [14, 15]. This possibility is provided by selecting an appropriate density of the nanoparticles and their size, as well as the wavelength of light. In terms of the light propagation in optical fibers, an equation similar to Eq. (6), with $x$ replaced by the reduced temporal variable $\tau\equiv t-z/V_{\mathrm{gr}}$ (here $t$ is time, and $V_{\mathrm{gr}}$ is the group velocity of the carrier wave), may also give rise to the critical collapse [16]. However, the collapse in optical fibers is easily suppressed by higher-order effects, such as the stimulated Raman scattering [17].

The family of TS solutions to Eq. (6) is commonly known,

\psi_{\mathrm{TS}}=\exp\left(ikz\right)\frac{\left(3k\right)^{1/4}}{\sqrt{\cosh\left(\sqrt{8k}x\right)}},

(7)

where real propagation constant $k$ takes values in the semi-infinite interval, $0<k<\infty$ . In accordance with what is said above, the family is degenerate, as the power of all the solutions takes the single value, which does not depend of $k$ :

P_{\mathrm{TS}}=\int_{-\infty}^{+\infty}\left|\psi\left(x\right)\right|^{2}dx=\frac{\sqrt{6}}{4}\pi\approx 1.92.

(8)

The well-known necessary stability condition for soliton families, if they are characterized by the dependence of the power on the propagation constant, is the Vakhitov-Kolokolov (VK) criterion [18, 1, 2, 3, 4], $dP/dk>0$ . Generally, it is only necessary but not sufficient for the stability of the solitons. For the degenerate TS family, with $dP/dk=0$ , the VK criterion yields a neutral result. Actually, as mentioned above, TSs are always unstable, while the formally neutral VK prediction corresponds to the above-mentioned fact that their instability against small perturbations is subexponential, corresponding to the pair of zero eigenvalues for small perturbations.

Dual-core waveguides, alias couplers, are systems which realize another fundamental effect, viz., spontaneous symmetry breaking (SSB) – in particular, in optics [19, 20, 21, 22, 23] and BEC [24, 25]. In the latter setting, SSB was experimentally realized about 20 years ago [25], while SSB in nonlinear optical couplers, predicted still earlier [20], was demonstrated experimentally for solitons in dual-core fibers 31 years later [26].

The well-elaborated theoretical predictions for the critical collapse and SSB, and the recent experimental works which have demonstrated these phenomena separately [10, 26], suggest a possibility to look for new physical effects produced by their interplay in dual-core waveguides. In particular, it is interesting to use such systems as a means for the stabilization of TSs. It is demonstrated below, by means of an exact solution, that this is possible indeed. The exact solution offers a reference point for the study of stabilized TS modes in a more general form.

In the simplest case, the system combining the critical collapse and SSB is based on the system of linearly coupled NLS equations with the quintic intra-core nonlinearity,

	$\displaystyle i\psi_{z}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\psi_{xx}-\|\psi\|^{4}\psi-\lambda\phi,$		(9)
	$\displaystyle i\phi_{z}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\phi_{xx}-\|\phi\|^{4}\phi-\lambda\psi,$		(10)

where $\lambda>0$ is a real coupling constant. It can be realized experimentally as a set of two parallel planar optical waveguides with the nonlinearity adjusted to the quintic form as outlined above.

Soliton solutions to Eqs. (9) and (10) with propagation constant $k$ are looked for as

\left\{\psi(x,z),\phi(x,z)\right\}=\exp\left(ikz\right)\left\{u(x),v(x)\right\},

(11)

where spatially even real functions $u(x)$ and $v(x)$ satisfy stationary equations

	$\displaystyle-ku$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}u}{dx^{2}}-u^{5}-\lambda v,$		(12)
	$\displaystyle-kv$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}v}{dx^{2}}-v^{5}-\lambda u,$		(13)

with boundary conditions $u(|x|\rightarrow\infty)=v(|x|\rightarrow\infty)=0$ . The solitons are characterized by powers of their two components,

P_{u}\equiv\int_{-\infty}^{+\infty}u^{2}(x)dx,P_{v}\equiv\int_{-\infty}^{+\infty}v^{2}(x)dx,

(14)

and the total power, $P\equiv P_{u}+P_{v}$ , the latter one being a dynamical invariant of the system.

It is well known that, in the case of the cubic self-focusing, with terms $|\psi|^{4}\psi$ and $|\phi|^{4}\phi$ in Eqs. (9) and (10) replaced by $|\psi|^{2}\psi$ and $|\phi|^{2}\phi$ , respectively, the competition between the intra-core self-focusing and linear inter-core coupling leads to the SSB bifurcation (alias symmetry-breaking phase transition), which destabilizes obvious symmetric soliton solutions, with $\psi=\phi$ , and replaces them by stable asymmetric ones [22, 23] when $P$ and $k$ exceed the corresponding critical values,

P_{\mathrm{SSB}}^{\mathrm{(cubic)}}=8\sqrt{\lambda}/3,k_{\mathrm{SSB}}^{\mathrm{(cubic)}}=5\lambda/3.

(15)

For the system of Eqs. (9) and (10), a similar analysis demonstrates that the SSB takes place at $k\geq k_{\mathrm{SSB}}^{\mathrm{(quintic)}}=5\lambda/4$ . On the other hand, all the symmetric solitons produced by Eqs. (9) and (10) assume a single value of the power, which is $2P_{\mathrm{TS}}=\pi\sqrt{3/2}$ , according to Eq. (8), therefore there is no special value of $P_{\mathrm{SSB}}^{\mathrm{(quintic)}}$ . Finally, the analysis of the Hamiltonian of the system of Eqs. (9) and (10), similar to that leading to Eq. (3), leads to the conclusion that the system with the quintic self-focusing and simplest linear coupling can only create unstable solitons of the TS type.

The objective of this work is to introduce a new system of linearly coupled NLS equations with the quintic self-focusing, which admits exact stable solutions for symmetric and asymmetric solitons (thus, the system provides an exact solution for the SSB phase transition). Unlike the system of Eqs. (9) and (10), the new one includes localized coupling (implying that a fused [32, 33] dual-core system is considered), and also a local attractive potential attached to the coupling spot, which helps to stabilize the solitons in the system. The model is introduced, with necessary details, in Section II. The exact solutions for solitons, which are symmetric and asymmetric with respect to the coupled cores, are reported in Section III. The analysis of the solitons’ stability makes it necessary to use numerical methods, with the results summarized in Section IV. In particular, it is found that asymmetric solitons, created by the SSB, are stable in a finite interval of values of the propagation constant. Evolution of those symmetric and asymmetric states which are unstable leads to establishment of either moderately asymmetric breathers (periodically oscillating localized modes) or strongly asymmetric states with small intrinsic vibrations. The work is concluded by Section V, where, in particular, we provide estimates of physical parameters for the experimental realization of the system, and solitons maintained by it, in optical couplers.

II The model

II.1 Local attractive potentials

II.1.1 The linear potential

A possibility for the stabilization of TSs in the 1D single-component NLS equation with the quintic self-focusing was elaborated in Ref. [27], where a delta-functional attractive potential (local defect) was added to Eq. (6):

i\psi_{z}=-\frac{1}{2}\psi_{xx}-|\psi|^{4}\psi-\varepsilon\delta(x)\psi,

(16)

with potential strength $\varepsilon>0$ . In optics, it may be realized as a narrow stripe with a higher refractive index embedded in the planar waveguide. A family of exact solutions for solitons pinned to the attractive defect is

\psi\left(x,z\right)=\exp\left(ikz\right)\frac{\left(3k\right)^{1/4}}{\sqrt{\cosh\left(2\sqrt{2k}(|x|+\xi\right)}},

(17)

cf. Eq. (7), with offset parameter $\xi>0$ determined by equation

\tanh\left(2\sqrt{2k}\xi\right)=\varepsilon/\sqrt{2k}.

(18)

Because $\tanh$ takes only values $<1$ , the solution for the pinned solitons, admitted by Eq. (18), exists in interval

\varepsilon^{2}/2<k<\infty.

(19)

For the soliton family given by Eqs. (17) and (18), the dependence of the power on the propagation constant is [27]

P(k)=\sqrt{\frac{3}{2}}\left[\pi-2\arctan\left(\sqrt{\frac{\sqrt{2k}+\varepsilon}{\sqrt{2k}-\varepsilon}}\right)\right].

(20)

This expression starts from $P=0$ at the left edge of the existence range (19), $k=\varepsilon^{2}/2$ , and attains the limit value (8) in the limit of $k\rightarrow\infty$ . Thus, the attractive delta-functional potential lifts the power degeneracy of the TS family (7).

Dependence $P(k)$ given by Eq. (20) satisfies the VK criterion. In agreement with this fact, it was found, by means of numerical methods, that the pinned-soliton family is completely stable [27]. The same result can be predicted by the estimate of the Hamiltonian, following the pattern of Eq. (3): the estimate of the extra term which accounts for the attractive potential in Eq. (16), $H_{\varepsilon}=-\varepsilon\left|\psi(x=0)\right|^{2}$ , yields $H_{\varepsilon}\sim-P/L$ , cf. Eq. (3). Then, at $P<P_{\mathrm{TS}}$ [note that Eq. (20) indeed yields values of $P$ which are smaller than $P_{\mathrm{TS}}$ given by Eq. (8)], the analysis of the estimate for $H$ , including term $H_{\varepsilon}$ , predicts a minimum of the Hamiltonian at a finite value of $L$ , which corresponds to stable solitons.

II.1.2 The nonlinear potential

For the sake of comparison with the model based on Eq. (16), it is relevant to mention one in which the attractive defect is not linear, but cubic, which was introduced in Ref. [28]:

i\phi_{z}=-\frac{1}{2}\phi_{xx}-\delta(x)|\phi|^{2}\phi.

(21)

In optics, it describes a planar waveguide with a narrow strip of strong local nonlinearity, which may be built by the respective distribution of the density of dopants that induce the local self-focusing [29]. In terms of BEC, Eq. (21), with $z$ replaced by $t$ , is the Gross-Pitaevskii equation with the local self-attractive nonlinearity, which may be imposed by a tightly focused laser beam through the optically-induced Feshbach resonance [30, 31]. Equation (21) gives rise to a family of pinned states,

\phi=\left(2k\right)^{1/4}\exp\left(ikz-\sqrt{2k}|x|\right),

(22)

whose power does not depend on $k$ , viz., $P=1$ , i.e., Eq. (21) is another 1D model which gives rise to the TS family [hence solutions (22) are unstable].

Further, it is possible to consider the 1D NLS equation which combines nonlinear terms from both equations (6) and (21) which give rise to the critical collapse, namely,

i\phi_{z}=-\frac{1}{2}\phi_{xx}\mp|\phi|^{4}\phi-\varepsilon\delta(x)|\phi|^{2}\phi.

(23)

Unlike Eq. (16), coefficient $\varepsilon>0$ in Eq. (23) cannot be eliminated by rescaling. Exact solutions to Eq. (23) for pinned solitons take the form [cf. Eq. (17)]

\left\{\begin{array}[]{c}\phi_{-}\left(x,z\right)\\ \phi_{+}\left(x,z\right)\end{array}\right\}=\left(3k\right)^{1/4}\exp\left(ikz\right)\left\{\begin{array}[]{c}1/\sqrt{\cosh\left(2\sqrt{2k}(|x|+\xi_{\_}\right)}\\ 1/\sqrt{\sinh\left(2\sqrt{2k}(|x|+\xi_{+}\right)}\end{array}\right\},

(24)

with offset $\xi_{\mp}$ defined by equation [cf. Eq. (18)]

\left\{\sinh\left(2\sqrt{2k}\xi_{-}\right),\cosh\left(2\sqrt{2k}\xi_{+}\right)\right\}=\sqrt{\frac{3}{2}}\varepsilon.

(25)

Here subscripts $-$ and $+$ correspond to the same signs in Eq. (23), where they represent the self-focusing/defocusing bulk nonlinearity, respectively. Note that, in the latter case, Eq. (25) produces a solution if the localized cubic self-focusing is stronger than the bulk defocusing, viz., $\varepsilon>\sqrt{2/3}$ . The calculation of the power for these soliton families yields

	$\displaystyle P_{-}(\varepsilon)=\sqrt{6}\left[\frac{\pi}{2}-\arctan\left(\sqrt{\frac{3}{2}}\varepsilon+\sqrt{\frac{3}{2}\varepsilon^{2}+1}\right)\right];$
	$\displaystyle P_{+}(\varepsilon)=\frac{1}{2}\sqrt{\frac{3}{2}}\ln\left(\frac{\varepsilon+\sqrt{2/3}}{\varepsilon-\sqrt{2/3}}\right).$		(26)

Both expressions (26) again do not depend on $k$ , i.e., the soliton families produced by Eqs. (24) and (25), with the combination of two nonlinear terms, also belong to the class of TS solutions, hence they are unstable too.

II.2 The system with fused coupling

As suggested by Refs. [34, 35], stable two-component solitons can be supported by a system with localized (fused [32, 33]) optical coupling, which is modeled by equations

	$\displaystyle i\psi_{z}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\psi_{xx}-\|\psi\|^{4}\psi-\delta(x)\phi-\varepsilon\delta(x)\psi,$		(27)
	$\displaystyle i\phi_{z}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\phi_{xx}-\|\phi\|^{4}\phi-\delta(x)\psi-\varepsilon\delta(x)\phi,$		(28)

cf. the model with the uniform coupling, provided by Eqs. (9) and (10). The fused coupling, approximated by the delta-functional terms with the coefficient scaled to be $1$ in Eqs. (27) and 28), can be realized as a bridge connecting two parallel waveguiding cores [34, 35]. The bridge also introduces a local increase of the thickness in each core, which may be modeled, in the simplest form, by the delta-functional potential terms with coefficient $\varepsilon>$ $0$ (the latter terms were not considered in Refs. [34, 35]). The local attractive potential draws the wave field to the spot at which the coupler is installed, thus helping to strengthen the coupling effect. In the fused-coupler setup, the attractive potential can be enhanced by a locally infused resonant dopant.

In addition to its potential use in photonic devices operating in the coupling regime [36, 37], the system based on Eqs. (27) and (28) offers a benefit for the theoretical study, as it provides exact analytical solutions for the symmetry-breaking phenomenology of TS families, as shown in the next section. Further, the analytical solutions provide clues for the understanding of the operation of a realistic setup, in which the ideal delta-function is replaced by its finite-width regularization, as is demonstrated below in the numerical part of the work.

III Analytical solutions

III.1 The solution ansatz

Stationary solutions to Eqs. (27) and (28), with real propagation constant $k$ , are looked for as per Eq. (11), with real functions $u(x)$ and $v(x)$ satisfying equations

	$\displaystyle-ku$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}u}{dx^{2}}-u^{5}-\delta(x)\left(v+\varepsilon u\right),$		(29)
	$\displaystyle-kv$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}v}{dx^{2}}-v^{5}-\delta(x)\left(u+\varepsilon v\right),$		(30)

cf. Eqs. (12) and (13). The solution to Eq. (16) represented by expression (17) suggests an ansatz for exact solutions to Eqs. (29) and (30):

	$\displaystyle u(x)$	$\displaystyle=$	$\displaystyle\frac{\left(3k\right)^{1/4}}{\sqrt{\cosh\left(\sqrt{8k}(\|x\|+\xi\right)}},$		(31)
	$\displaystyle v(x)$	$\displaystyle=$	$\displaystyle\frac{\left(3k\right)^{1/4}}{\sqrt{\cosh\left(\sqrt{8k}(\|x\|+\eta\right)}}.$		(32)

Here, possible asymmetry between the two components of the soliton in the symmetric system is represented by different positive offset parameters, $\xi\neq\eta$ .

The effect of the delta-functional terms in Eqs. (29) and (30) is represented by jumps (discontinuities) of the first derivatives of $u(x)$ and $v(x)$ at point $x=0$ , which are produced by the integration of Eqs. (29) and (30) in an infinitesimal vicinity of $x=0$ , while functions $u(x)$ and $v(x)$ are continuous at this point:

	$\displaystyle\frac{du}{dx}(x$	$\displaystyle=$	$\displaystyle+0)-\frac{du}{dx}(x=-0)=-2\left[v(x=0)+\varepsilon u(x=0)\right],$		(33)
	$\displaystyle\frac{dv}{dx}(x$	$\displaystyle=$	$\displaystyle+0)-\frac{dv}{dx}(x=-0)=-2\left[u(x=0)+\varepsilon v(x=0)\right].$		(34)

The substitution of expressions (31) and (32) in Eqs. (33) and (34) leads to the following system of equations for $\xi$ and $\eta$ :

	$\displaystyle 2k\left(c_{1}^{2}-1\right)c_{2}$	$\displaystyle=$	$\displaystyle\left(c_{1}+\varepsilon^{2}c_{2}+2\varepsilon\sqrt{c_{1}c_{2}}\right)c_{1}^{2},$		(35)
	$\displaystyle 2k\left(c_{2}^{2}-1\right)c_{1}$	$\displaystyle=$	$\displaystyle\left(c_{2}+\varepsilon^{2}c_{1}+2\varepsilon\sqrt{c_{1}c_{2}}\right)c_{2}^{2},$		(36)

where

c_{1}\equiv\cosh\left(\sqrt{8k}\xi\right),c_{2}\equiv\cosh\left(\sqrt{8k}\eta\right).

(37)

III.2 Symmetric states

First, for symmetric solutions, with $\xi=\eta$ , i.e., $c_{1}=c_{2}$ in Eq. (37), the matching conditions given by Eqs. (35) and (36) coalesce into a single equation,

\tanh\left(2\sqrt{2k}\xi\right)=\frac{1+\varepsilon}{\sqrt{2k}},

(38)

cf. Eq. (18). The constraint $\tanh\left(2\sqrt{2k}\xi\right)<1$ defines the existence range of the symmetric solitons,

\frac{1}{2}\left(1+\varepsilon\right)^{2}<k<\infty,

(39)

which is identical to range (19) in the case of $\varepsilon=0$ . The power of the symmetric soliton can be easily obtained from Eqs. (31), (32), and (38):

P_{\mathrm{symm}}(k)=\sqrt{6}\left[\pi-2\arctan\left(\sqrt{\frac{\sqrt{2k}+(1+\varepsilon)}{\sqrt{2k}-(1+\varepsilon)}}\right)\right].

(40)

This expression starts from $P_{\mathrm{symm}}=0$ at the left edge the existence range (39) and monotonously grows with the increase of $k$ , up to $P_{\mathrm{symm}}(k\rightarrow\infty)=\sqrt{6}\pi/2\equiv 2P_{\mathrm{TS}}$ , see Eq. (8). Thus, similar to Eq. (20, for any value of $\varepsilon$ Eq. (40) meets the VK criterion, $dP/dk>0$ . Furthermore, the symmetric solution given by Eqs. (31), (32), and (38) exists even for $-1<\varepsilon<0$ , when the local potential is moderately repulsive, as the effective attractive potential induced by the fused coupling is stronger.

III.3 Asymmetric states

III.3.1 The case of $\varepsilon=0$

It is easy to find a solution of Eqs. (35) and (36) with broken symmetry, i.e., $c_{1}\neq c_{2}$ , in the case of $\varepsilon=0$ , when the system does not include the delta-functional potential:

c_{1,2}^{2}=2k^{2}\pm 2k\sqrt{k^{2}-1}.

(41)

This solution exists, obviously, in the range of

1<k<\infty.

(42)

Note that condition $c_{1,2}>1$ , which follows from definition (37), holds for solution (41) in interval (42). Comparison with Eq. (39) demonstrates that, in the case of $\varepsilon=0$ , only the symmetric states exist in the narrow interval of $1/2<k<1$ , while the semi-infinite area (42) is populated by both the symmetric and asymmetric modes.

The analytical calculation of the power for the asymmetric solution, which is given, for $\varepsilon=0$ , by Eqs. (31), (32), and (41), yields a surprising exact result: it does not depend on $k$ [on the contrary to power (40) of the symmetric solitons in the same case of $\varepsilon=0$ ], keeping the value given by Eq. (8). It is easy to check that, naturally, this single value of the norm is identical to the norm of the symmetric solution at $k=1$ , which is given by the doubled expression (40) with $k=1$ and $\varepsilon=0$ . The fact that the family of the asymmetric solitons is degenerate, featuring the single value of the power, implies that it is another species of TS family, thus being certainly unstable.

The asymmetry of the stationary states is characterized by the power-imbalance parameter,

\Theta\equiv\pm\frac{P_{u}-P_{v}}{P_{u}+P_{v}},

(43)

where the opposite signs correspond to two branches of the asymmetric states, that are mirror images of each other, see Figs. 4, 5, and 8 below. The analytical solution (41) yields

	$\displaystyle\Theta(k,\varepsilon$	$\displaystyle=$	$\displaystyle 0)=\pm\frac{4}{\pi}\left[\arctan\left(\sqrt{2k^{2}+2k\sqrt{k^{2}-1}}+\sqrt{2k^{2}+2k\sqrt{k^{2}-1}-1}\right)\right.$		(44)
			$\displaystyle\left.-\arctan\left(\sqrt{2k^{2}-2k\sqrt{k^{2}-1}}+\sqrt{2k^{2}-2k\sqrt{k^{2}-1}-1}\right)\right].$		(44)

With the increase of $k$ from $1$ to $\infty$ , as per Eq. (42), $\Theta$ given by Eq. (44) monotonously varies from $0$ to $\pm 1$ , the latter limit implying that one core of the coupler gets asymptotically empty. Usually, the SSB transition is characterized by dependence $\Theta(P)$ , rather than $\Theta(k)$ , because the power and asymmetry are observable parameters of optical beams, unlike the “hidden” propagation constant. However, in the present case it is irrelevant to refer to $\Theta(P)$ , as the degenerate family of the asymmetric solitons of the TS type keeps the single value (8) of the power at $\varepsilon=0$ .

III.3.2 The general case, $\varepsilon\neq 0$

In the case of $\varepsilon\neq 0$ , equations (35) and (36) are too complex for the full analytical solution, Nevertheless, analysis of these equations makes it possible to identify the point at which SSB commences, i.e., a solution with an infinitesimal difference $c_{1}-c_{2}$ branching off from the symmetric one with $c_{1}=c_{2}\equiv c$ . The linearization of Eqs. (35) and (36) with respect to infinitesimal $\left(c_{1}-c_{2}\right)$ demonstrates that this bifurcation happens at

c=c_{\mathrm{SSB}}\equiv\sqrt{2+\varepsilon},~{}k\equiv k_{\mathrm{SSB}}=\frac{1}{2}\left(1+\varepsilon\right)\left(2+\varepsilon\right),

(45)

i.e., the asymmetric solitons exist in the range of

\frac{1}{2}\left(1+\varepsilon\right)\left(2+\varepsilon\right)<k<\infty.

(46)

In the case of $\varepsilon=0$ , this range is identical to one given above by Eq. (42).

It is also instructive to compare values of the power at $k=k_{\mathrm{SSB}}$ and $k\rightarrow\infty$ . It is easy to find

P\left(k=k_{\mathrm{SSB}}\right)=\sqrt{6}\left[\pi-2\arctan\left(\sqrt{\varepsilon+2}+\sqrt{\varepsilon+1}\right)\right]

(47)

[in the case of $\varepsilon=0$ , Eq. (47) gives the value identical to one given by Eq. (8)], while, in the limit of $k\rightarrow\infty,$ the power is the same as given above by Eq. (8). Then, an essential conclusion is that

	$\displaystyle P_{\mathrm{asymm}}\left(k=k_{\mathrm{SSB}}\right)$	$\displaystyle<$	$\displaystyle P_{\mathrm{asymm}}(k\rightarrow\infty)\text{,~{}at }\varepsilon>0,$
	$\displaystyle P_{\mathrm{asymm}}\left(k=k_{\mathrm{SSB}}\right)$	$\displaystyle\geq$	$\displaystyle P_{\mathrm{asymm}}(k\rightarrow\infty)\text{,~{}at }\varepsilon\leq 0.$		(48)

Actually, Eq. (48) implies that the VK criterion holds at $\varepsilon>0$ , and does not hold at $\varepsilon\leq 0$ . Therefore, the asymmetric solitons may be stable in the former case (with the attractive local potential), and are definitely unstable if the potential is repulsive. It is shown below that a part of the family of the asymmetric solitons is indeed stable for $\varepsilon>0$ .

IV Numerical results

IV.1 The setup for the numerical analysis

The exact solution for the asymmetric states, given by Eqs. (35) and (36), cannot be cast in a fully explicit form for $\varepsilon\neq 0$ , and the full stability analysis cannot be performed analytically beyond the verification of the VK criterion. These problems should be addressed by means of numerical methods. The numerical solutions and test of their stability are produced below in the framework of the system in which the ideal delta-function in Eqs. (52) and (28) is replaced by the usual Gaussian regularization,

\tilde{\delta}(x)=\left(\sqrt{\pi}\sigma\right)^{-1}\exp\left(-\frac{x^{2}}{\sigma^{2}}\right),

(49)

with small width $\sigma$ . We here fix $\sigma=0.02$ , which is sufficiently small in comparison with the transverse width of various trapped modes. Comparison between the exact analytical solutions produced above for the system with the ideal delta-function and numerical results obtained with the help of regularization (49) is demonstrated below in Figs. 6 and 7. Note also that the system with the ideal delta-function replaced by regularization (49) is appropriate for modeling real fused couplers, with a finite size of the coupling region.

Numerical solution of Eqs. (29) and (30) with $\delta(x)$ replaced by $\tilde{\delta}(x)$ was performed by means of the Newton-conjugate-gradient method, which readily provides robust convergence to stationary solutions of the system [38, 39]. Stability of the stationary solutions against small perturbations, taken in the usual form [39],

	$\displaystyle\psi_{1}(x,z)$	$\displaystyle=$	$\displaystyle\exp(ikz)\left[\exp(\gamma z)u_{1}(x)+u_{2}^{\ast}(x)\exp(\gamma^{\ast}z)\right],$
	$\displaystyle\phi_{1}(x,z)$	$\displaystyle=$	$\displaystyle\exp(ikz)\left[\exp(\gamma z)v_{1}(x)+v_{2}^{\ast}(x)\exp(\gamma^{\ast}z)\right],$		(50)

is determined by the linearized equations for perturbation amplitudes $u_{1,2}(x)$ and $v_{1,2}(x)$ and instability eigenvalue (growth rate) $\gamma$ :

$\displaystyle-(k-i\gamma)u_{1}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}u_{1}}{dx^{2}}-3u^{4}(x)u_{1}-2u^{4}(x)u_{2}-\delta(x)\left(v_{1}+\varepsilon u_{1}\right),$
$\displaystyle-(k+i\gamma)u_{2}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}u_{2}}{dx^{2}}-3u^{4}(x)u_{2}-2u^{4}(x)u_{1}-\delta(x)\left(v_{2}+\varepsilon u_{2}\right),$
$\displaystyle-(k-i\gamma)v_{1}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}v_{1}}{dx^{2}}-3v^{4}(x)v_{1}-2v^{4}(x)v_{2}-\delta(x)\left(u_{1}+\varepsilon v_{1}\right),$
$\displaystyle-(k+i\gamma)v_{2}$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\frac{d^{2}v_{2}}{dx^{2}}-3v^{4}(x)v_{2}-2v^{4}(x)v_{1}-\delta(x)\left(u_{2}+\varepsilon v_{2}\right).$	(51)

The ordinary stability condition is that all eigenvalues $\gamma$ produced by Eqs. (51) must have zero real parts. The spectrum of the eigenvalues was produced by means of the Fourier collocation method.

Predictions for the (in)stability of the stationary states, provided by the numerical computation of the eigenvalues for small perturbations, were verified by means of simulations of the perturbed evolution in the framework of Eqs. (27) and (28), again with $\delta(x)$ replaced by $\tilde{\delta}(x)$ . Characteristic examples of stable and unstable evolution are produced below.

IV.2 Stable and unstable families of symmetric and asymmetric bound states

First, results of the numerical solution of Eqs. (29) and (30), with $\delta(x)$ replaced by regularization (49), are collected, in the form of dependences $P(k)$ for the families of symmetric and asymmetric bound states, in Fig. 1. The plots also present conclusions for the stability, as obtained from the numerical solution of Eqs. (51). Typical examples of the symmetric and asymmetric solitons, both stable and unstable ones, are displayed in Fig. 2.

These results are reported for the strength of the attractive potential varying in the interval of

0.15\leq\varepsilon\leq 2.5.

(52)

At $\varepsilon<0.15$ , the findings are virtually identical to those for $\varepsilon=0.15$ . The case of $\varepsilon>2.5$ is not presented here, as in that case the strong attractive potential makes the bound state so narrow that its width is no longer much larger than the fixed width $\sigma=0.02$ of the regularized delta-function, adopted in Eq. (49); for this reason, the setup becomes essentially different from the one defined by Eqs. (27) and 28). This is shown, in particular, below in Fig. 7.

Refer to caption — Figure 1: Power $P$ of symmetric and asymmetric bound states (solitons) vs. the propagation constant $k$ for different values of strength $\varepsilon$ of the attractive potential in Eqs. (27) and (28), which are indicated in panels. Note essentially different scales of $k$ in different panels. Families of asymmetric states are represented by nearly horizontal branches originating at SSB (spontaneous-symmetry-breaking) bifurcation points. The results are produced by the numerical solution of Eqs. (29) and (30), where the ideal delta-function is replaced by the narrow Gaussian (49) with $\sigma=0.02$ . The families of symmetric states are stable and unstable before and after the bifurcation (to the left and right of the SSB points, respectively). The branches of the asymmetric states are stable in segments between the SSB points and stability boundaries designated by arrows. The stability segments are extremely narrow for $\varepsilon\leq 0.5$ , in which cases the arrows are located very close to the SSB points. The stability is identified according to eigenvalues produced by the numerical solution of Eqs. (51).

Shapes of the symmetric and asymmetric solitons, and their evolution following the increase of the propagation constant, $k$ , are displayed in Fig. 2. The figure clearly shows the change in the shape triggered by the onset of SSB at the point which is indicated by arrows.

Naturally, symmetric solitons are stable below the SSB point, and unstable above it. As shown above by means of the VK criterion, the asymmetric solitons are completely unstable at $\varepsilon=0$ . As an extension of that finding, Fig. 1 shows that the asymmetric states remain almost completely unstable at $\varepsilon\leq 0.5$ , while a visible stability interval, of width $\Delta k$ , appears above the SSB point at $\varepsilon>0.5$ .

The dependence of width $\Delta k$ of the stability interval for the asymmetric bound states on $\varepsilon$ is shown in Fig. 3(a). In addition, panel (b) of Fig. 3 shows values of the power for the asymmetric bound states at the midpoint of the respective stability intervals. These results, produced by the computation of stability eigenvalues, are confirmed by direct simulations of the perturbed evolution of the asymmetric bound states, as shown below in Figs. 9-12.

It is relevant to present the asymmetric-soliton families in terms of the asymmetry parameter $\Theta$ , defined as per (43). To this end, dependences $\Theta(k)$ and $\Theta(P)$ are plotted in Figs. 4 and 5, respectively. The virtually flat vertical shape of the $\Theta(P)$ lines for $\varepsilon<0.5\$ in Fig. 5 corresponds to the above analytical result obtained for $\varepsilon=0$ , according to which the power of the whole family of the asymmetric solitons takes the single value (8). It is relevant to stress that the $\Theta(P)$ curves which are not fully flat feature the forward-going (convex) shape, clearly indicating that the SSB bifurcation in the present system is of the supercritical type [40], i.e., it represents the phase transition of the second kind.

IV.3 Comparison of the numerical and analytical results

It is relevant to compare the analytical predictions, obtained above in the model based on Eqs. (27) and (28) with the ideal delta-function, and systematic numerical results produced for the system with the delta-function replaced by the regularized expression (49). First, in Fig. 6 we compare the analytical values of the power and propagation constant at the SSB point, as given by Eq. (15), and their numerically found counterparts, as functions of $\varepsilon$ in the interval (52) under the consideration. The top panel of Fig. 6 clearly demonstrates that the analytically predicted critical values of $P$ are virtually identical to the numerical ones. Further, the bottom panel shows that the analytical and numerical critical values of $k$ are almost identical at $\varepsilon<1$ , and a relatively small discrepancy appears at $\varepsilon>1$ . As mentioned above, it is explained by the fact that strong attraction to the central point by the local potential makes the trapped state so narrow that its width ceases to be much larger than the finite size $\sigma$ of the regularized delta-functional profile in Eq. (49).

Further comparison of the analytical and numerical results is provided by Fig. 7, in terms of the $P(k)$ curves for the families of symmetric solitons (the numerical curves are essentially the same as shown, also for the symmetric states, in Fig. 1). This picture again demonstrates close proximity of the analytical and numerical results at $\varepsilon<1$ , and gradual increase of the discrepancy at $\varepsilon>1$ .

It is also natural to compare the analytical dependence of the asymmetry parameter, $\Theta(k)$ , given for $\varepsilon=0$ by Eq. (44), and its numerically generated counterpart, which is shown in Fig. 8. It is seen that both dependences match perfectly, being also very close to the numerically generated $\Theta(k)$ dependence plotted for $\varepsilon=0.15$ in Fig. 4(a).

IV.4 Perturbed evolution of stable and unstable bound states

The predictions for the (in)stability, which are produced above in the analytical form by means of the VK criterion and by dint of the numerical solution of the eigenvalue problem based on Eq. (51), have been verified by direct simulations of Eqs. (27) and (28), with $\delta(x)$ again replaced by approximation (49) with $\sigma=0.02$ . First, the simulations of the symmetric solitons, which are reported in Figs. 9 and 10 for weaker and stronger attractive potentials, viz., $\varepsilon=1$ and $2.5$ , demonstrate, as expected, that these states remain stable at $k<k_{\mathrm{SSB}}$ , and develop spontaneous instability at $k>k_{\mathrm{SSB}}$ . A clear trend is that the instability makes amplitudes of the originally symmetric components different, and initiates Josephson oscillations between them. In accordance with the fact that the solution of Eq. (51) produces purely real instability eigenvalues $\gamma$ for the symmetric solitons at $k>k_{\mathrm{SSB}}$ , the oscillation frequency is small close to the SSB point [e.g., at $k=2.7$ in Fig. 9(b)]. As a result, the system establishes robust breathers with different amplitudes in their two components, as seen in Figs. 9(b,c) and 10(b,c).

Farther from the bifurcation (at larger values of $k$ ), the oscillations are much faster, and the system quickly develops a state with a large difference in amplitudes of the two components [e.g., at $k=3$ in Fig. 9(d), and at $k=7.5$ in Fig. 10(d)]. The latter states seem as robust strongly asymmetric ones with irregular small-amplitude intrinsic vibrations. They approximately resemble stationary bound states with strong asymmetry between the components; however, strictly stationary asymmetric states do not exist for the same values of parameters. For instance, in the two above-mentioned cases, corresponding to $k=3$ in Fig. 9(d) and $k=7.5$ in Fig. 10(d), the respective value of the power of the symmetric solitons is $P\approx 1.6$ in either case, as seen from Figs. 1(e) and (j), respectively; on the other hand, the same figures demonstrate that stationary asymmetric states do not exist for this $P$ – for instance, the largest available power of the asymmetric solitons is $\approx 1.5$ in the former case, and only $\approx 1.3$ in the latter one. Thus, the conclusion is that robust modes with strong asymmetry and large values of the power are admitted by the system, but in the weakly vibrating state.

Simulations of the evolution of asymmetric states corroborate that these solitons are also stable or unstable in the cases when the numerical solution of Eq. (51) produces, respectively, zero or nonzero values of the instability growth rate, Re $(\gamma)$ . Typical examples for the same values, $\varepsilon=1$ and $2.5$ , as used in Figs. 9 and 10, are displayed, severally, in Figs. 11 and 12. It is observed that the instability emerges in an oscillatory form, with a finite frequency, in accordance with the fact that the instability eigenvalues for the asymmetric states are found with nonzero imaginary parts, unlike the above-mentioned pure real unstable eigenvalues for the symmetric bound states.

In the case of relatively weak instability, the stationary asymmetric states are spontaneously transformed into apparently robust asymmetric breathers, which exhibit regular oscillations in Figs. 11(b,c) and 12(c). They are similar to the above-mentioned asymmetric breathers created by the instability of the symmetric states, as shown in Figs. 9(b,c) and 10(b,c).

On the other hand, strong instability, as seen in Figs. 11(d) and 12(d), creates states with a very strong asymmetry and small-amplitude irregular intrinsic vibrations. The robust trapped states of the latter type are similar to the above-mentioned strongly asymmetric modes with irregular internal oscillations which are observed in Figs. 9(d) and 10(d)

V Conclusion

This work aims to introduce a physically relevant model which admits the exact solution for two basic problems, viz., the stabilization of TSs (Townes solitons) and SSB (spontaneous symmetry breaking) in nonlinear couplers. The interest to these problems is drawn, in particular, by the recently reported first experimental demonstration of weakly unstable 2D TSs in BEC [10, 11] and SSB of optical solitons in dual-core fibers [26]. The present model represents the system of parallel 1D waveguides with the intrinsic quintic nonlinearity, which is critical in the 1D case, giving rise to the corresponding TSs. The waveguides form a fused coupler, with the linear connection introduced in a narrow region. In addition to that, the region of the local coupling carries the attractive potential acting in each waveguiding core. The system can be directly realized in optics, and it may be used in the design of photonic devices operating in coupling and switching regimes. The model admits full analytical solutions for symmetric and asymmetric solitons, revealing the explicit picture of the SSB in solitons and a straightforward scenario for the stabilization of the TSs, which are completely unstable in the uniform waveguides. It is found that the SSB bifurcation of the supercritical type, i.e., the phase transition of the second kind, destabilizes the symmetric bound states and gives rise to the asymmetric ones, which remain stable in a finite interval of values of the propagation constant. The evolution of those symmetric and asymmetric states which are unstable leads to robust moderately asymmetric breathers, which perform Josephson oscillations in the coupler, or (also robust) strongly asymmetric states with small-amplitude irregular internal vibrations.

It is relevant to estimate parameters of the fused optical coupler which can realize the proposed setup. An appropriate material is the above-mentioned colloidal suspension of silver nanoparticles, with the volume-filling factor $\simeq 1.5\times 10^{-5}$ and the corresponding quintic susceptibility $\chi^{(5)}\simeq$ $4\times 10^{-38}$ m ${}^{4}/$ V⁴ at a visible wavelength $\simeq 0.5$ $\mathrm{\mu}$ m, while the cubic nonlinearity is negligible [14, 15]. Appropriate thickness of each waveguide, as well as the width of the bridge connecting them, is $\simeq 3$ $\mathrm{\mu}$ m. Then, the predicted solitons, with the transverse size $\gtrsim 30$ $\mathrm{\mu}$ m, can be created by laser beams with power density about $10$ GW/cm² and total power $\sim 10$ kW. The propagation distance $\gtrsim 3$ cm is sufficient to create well-formed solitons.

As an extension of the analysis, it may be relevant to introduce a system with a double fused coupler, similar to the one considered in Ref. [35]. Such a system should make it possible to study double SSB effects – between the two cores of the coupler and, in the spatial direction, between two separated regions of the fused coupling.

Acknowledgment

This work was supported, in part, by the Israel Science Foundation through grant No. 1695/22.

References

[1] L. Bergé, Wave collapse in physics: principles and applications to light and plasma waves, Phys. Rep. 303, 259-370 (1998).
[2] C. Sulem and P.-L. Sulem, The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse (Springer, New York, 1999).
[3] V. E. Zakharov V. E. and E. A. Kuznetsov, Solitons and collapses: two evolution scenarios of nonlinear wave systems, Physics – Uspekhi 55, 535-556 (2012).
[4] G. Fibich, The Nonlinear Schrödinger Equation: Singular Solutions and Optical Collapse (Springer, Heidelberg, 2015).
[5] B. A. Malomed, Multidimensional Solitons (American Institute of Physics, Melville, NY, 2022).
[6] N. J. Zabusky and M. D. Kruskal, Interaction of “solitons” in a collisional plasma and the recurrence of initial states, Phys. Rev. Lett. 15, 240-243 (1965).
[7] Chiao R. Y., E. Garmire, and C. H. Townes, Self-trapping of optical beams, Phys. Rev. Lett. 13, 479-482 (1964).
[8] S. N. Vlasov, V. A. Petrishchev, and V. I. Talanov, Izv. Vyssh. Uchebn. Zaved. Radiofiz. 14, 1353-1363 (1971) [English translation: Radiophys. Quantum Electron. 14, 1062-1070 (1971)].
[9] M. Desaix, D. Anderson, and M. Lisak, Variational approach to collapse of optical pulses, J. Opt. Soc. Am. B 8, 2082-2086 (1991).
[10] B. Bakkali-Hassani B., C. Maury, Y.-Q. Zhou, E. Le Cerf, R. Saint-Jalm, P. C. M. Castilho, S. Nascimbene, J. Dalibard, and J. Beugnon, Realization of a Townes soliton in a two-component planar Bose gas, Phys. Rev. Lett. 127, 023603 (2021).
[11] B. Bakkali-Hassani, C. Maury, S. Stringari, S. Nascimbene, J. Dalibard, and J. Beugnon, The cross-over from Townes solitons to droplets in a 2D Bose mixture, New J. Phys. 25, 013007 (2023).
[12] F. Kh. Abdullaev and M. Salerno, Gap-Townes solitons and localized excitations in low-dimensional Bose-Einstein condensates in optical lattices, Phys. Rev. A 72, 033617 (2005).
[13] G. L. Alfimov, V. V. Konotop, and P. Pacciani, Stationary localized modes of the quintic nonlinear Schrödinger equation with a periodic potential 75, 023624 (2007).
[14] A. S. Reyna and C. B. de Araújo, Spatial phase modulation due to quintic and septic nonlinearities in metal colloids, Opt. Exp. 22, 22456-22469 (2014).
[15] A. S. Reyna and C. B. de Araújo, High-order optical nonlinearities in plasmonic nanocomposites – a review, Adv. Opt. Phot. 9, 720-774 (2017).
[16] L. Bergé, J. Juul Rasmussen, and J. Wyller, Dynamics of localized solutions to the Raman-extended derivative nonlinear Schrödinger equation, J. Phys. A: Math. Gen. 29, 3581-3595 (1996); J. S. Hesthaveny, J. Juul Rasmussen, L. Bergé, and J. Wyller, Numerical studies of localized wavefields governed by the Raman-extended derivative nonlinear Schrödinger equation, J. Phys. A: Math. Gen. 30, 8207-8224 (1997).
[17] Y. Kodama and A. Hasegawa, Nonlinear pulse propagation in a monomode dielectric guide, IEEE J. Quant. Elect. QE-23, 510-524 (1987).
[18] N. G. Vakhitov and A. A. Kolokolov, Stationary solutions of the wave equation in a medium with nonlinearity saturation, Radiophys. Quantum Electron. 16, 783-789 (1973); https://doi.org/10.1007/BF01031343.
[19] S. R. Friberg, Y. Silberberg, M. K. Oliver, M. J. Andrejko, M. A. Saifi, and P. W. Smith, Ultrafast all-optical switching in a dual-core fiber nonlinear coupler, Appl. Phys. Lett. 51, 1135-1137 (1987).
[20] E. M. Wright, G. I. Stegeman, and S. Wabnitz, Solitary-wave decay and symmetry-breaking instabilities in two-mode fibers. Phys. Rev. A 40, 4455-4466 (1989).
[21] 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, 2101-2118 (1991).
[22] M. Romagnoli, S. Trillo, and S. Wabnitz, Soliton switching in nonlinear couplers, Opt. Quantum Electron 24, S1237–S1267 (1992).
[23] B. A. Malomed, A variety of dynamical settings in dual-core nonlinear fibers, In: Handbook of Optical Fibers, Vol. 1, pp. 421-474 (G.-D. Peng, Editor: Springer, Singapore, 2019).
[24] 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).
[25] 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).
[26] V. H. Nguyen, L. X. T. Tai, I. Bugar, M. Longobucco, R. Buzcynski, B. A. Malomed, and M. Trippenbach, Reversible ultrafast soliton switching in dual-core highly nonlinear optical fibers, Opt. Lett. 45, 5221-5224 (2020).
[27] L. Wang, B. A. Malomed, and Z. Yan, Attraction centers and $\mathcal{PT}$ -symmetric delta-functional dipoles in critical and supercritical self-focusing media, Phys. Rev. E 99, 052206 (2019).
[28] B. A. Malomed and M. Ya. Azbel, Modulational instability of a wave scattered by a nonlinear center, Phys. Rev. B 47, 10402-10406 (1993).
[29] J. Hukriede, D. Runde, and D. Kip, Fabrication and application of holographic Bragg gratings in lithium niobate channel waveguides, J. Phys. D 36, R1 (2003).
[30] R. Yamazaki, S. Taie, S. Sugawa, and Y. Takahashi, Submicron spatial modulation of an interatomic interaction in a Bose-Einstein condensate, Phys. Rev. Lett. 105, 050405 (2010).
[31] L. W. Clark, L.-C. Ha, C.-Y. Xu, and C. Chin, Quantum dynamics with spatiotemporal control of interactions in a stable Bose-Einstein condensate, Phys. Rev. Lett. 115, 155301 (2015).
[32] A. Takagi, K. Jinguji, and M. Kawachi, Wavelength characteristics of ( $2\times 2$ ) optical channel-type directional couplers with symmetric or nonsymmetric coupling structures, J. Lightwave Tech. 10, 735-746 (1992).
[33] R. Ismaeel, T. Lee, B. Oduro, and Y. Jung, and G. Brambilla, All-fiber fused directional coupler for highly efficient spatial mode conversion, Opt. Exp. 22, 11610-11619 (2014).
[34] Y. Li, W. Pang, S. Fu, and B. A. Malomed, Two-component solitons under a spatially modulated linear coupling: Inverted photonic crystals and fused couplers, Phys. Rev. A 85, 053821 (2012).
[35] A. Harel and B. A. Malomed, Interactions of spatial solitons with fused couplers, Phys. Rev. A 89, 043809 (2014).
[36] C.-L. Chen, Foundations for Guided-Wave Optics (Wiley-Interscience, Hoboken, NJ, 2006).
[37] R. Halir, P. J. Bock, P. Cheben, A. Ortega-Monux, C. Alonso-Ramos, J. H. Schmid, J. Lapointe, D. X. Xu, J. G. Wanguemert-Perez, and I. Molina-Fernandez, Waveguide sub-wavelength structures: a review of principles and applications, Laser & Phot. Reviews 9, 25-49 (2015).
[38] J. Nocedal and S. J. Wright, Numerical Optimization (Springer, New York, 2006).
[39] J. Yang, Nonlinear Waves in Integrable and Nonintegrable Systems (SIAM, Philadelphia, 2010).
[40] G. Iooss and D. D. Joseph, Elementary Stability and Bifurcation Theory (Springer, Berlin, 1980).