Nonlinear dynamics of wave packets in tunnel-coupled harmonic-oscillator traps
NNonlinear dynamics of wave packets in tunnel-coupled harmonic-oscillator traps
Nir Hacker and Boris A. Malomed , Department of Physical Electronics, School of Electrical Engineering, Faculty of Engineering,and Center for Light-Matter interaction, Tel Aviv University, Tel Aviv 69978, Israel and Instituto de Alta Investigaci´on, Universidad de Tarapac´a, Casilla 7D, Arica, Chile
We consider a two-component linearly-coupled system with the intrinsic cubic nonlinearity and theharmonic-oscillator (HO) confining potential. The system models binary settings in BEC and optics.In the symmetric system, with the HO trap acting in both components, we consider Josephsonoscillations (JO) initiated by an input in the form of the HO’s ground state (GS) or dipole mode(DM), placed in one component. With the increase of the strength of the self-focusing nonlinearity,spontaneous symmetry breaking (SSB) between the components takes place in the dynamical JOstate. Under still stronger nonlinearity, the regular JO initiated by the GS input carry over into achaotic dynamical state. For the DM input, the chaotization happens at smaller powers than for theGS, which is followed by SSB at a slightly stronger nonlinearity. In the system with the defocusingnonlinearity, SSB does not take place, and dynamical chaos occurs in a small area of the parameterspace. In the asymmetric half-trapped system, with the HO potential applied to a single component,we first focus on the spectrum of confined binary modes in the linearized system. The spectrumis found analytically in the limits of weak and strong inter-component coupling, and numericallyin the general case. Under the action of the coupling, the existence region of the confined modesshrinks for GSs and expands for DMs. In the full nonlinear system, the existence region for confinedmodes is identified in the numerical form. They are constructed too by means of the Thomas-Fermiapproximation, in the case of the defocusing nonlinearity. Lastly, particular (non-generic) exactanalytical solutions for confined modes, including vortices, in one- and two-dimensional asymmetriclinearized systems are found. They represent bound states in the continuum.
I. INTRODUCTION
The combination of a harmonic-oscillator (HO) trapping potential and cubic nonlinearity is a ubiquitous settingwhich occurs in diverse microscopic, mesoscopic, and macroscopic physical settings. A well-known realization is offeredby Bose-Einstein condensates (BECs) with collisional nonlinearity [1–3], loaded in a magnetic or optical trap – see,e.g., Refs. [4]-[14]. A similar combination of the effective confinement, approximated by the parabolic profile of thelocal refractive index, and the Kerr term is relevant as a model of optical waveguides [20]-[23]. Models of the sametype appear in other physical systems too, such as networks of Josephson oscillators [24].The character of states created by the interplay of the intrinsic nonlinearity and externally applied trapping potentialstrongly depends on the sign of the nonlinearity. In the case of the self-attraction (or self-phase-modulation, SPM, interms of optics [25]), localized modes, similar to solitons, arise spontaneously. On the other hand, self-repulsion tendsto create flattened configurations, which, in turn, may support dark solitons in various static and dynamical states[6]-[14]. A specific situation occurs in a system combining repulsion and a weak parabolic potential with an additionalspatially periodic one (it represents an optical lattice in BEC [27], or a photonic crystal in optics and plasmonics[28–30]): the interplay of the periodic potential with the self-repulsion gives rise to bright gap solitons [27, 31], whoseeffective mass is negative. For this reason, gap solitons are expelled by the HO potential, but are trapped by theinverted one, that would expel modes with positive masses [32].Another noteworthy feature of the dynamics of one-dimensional (1D) nonlinear fields trapped in confining potentialsis the degree of its nonintegrability . The generic model for such settings is provided by the nonlinear Schr¨odingerequation (NLSE, alias the Gross-Pitaevskii equation, in terms of BEC [1–3]) with the cubic term, whose integrabilityin the 1D space [33] is broken by the presence of the HO potential. Nevertheless, systematic numerical simulations,performed for the repulsive sign of the cubic nonlinearity, have demonstrated that the long-time evolution governed byNLSE with the HO potential term does not leads to establishment of spatiotemporal chaos (“turbulence”), which wouldbe expected in the case of generic nonintegrability [34, 35]. Instead, the setup demonstrates quasi-periodic evolution,represented by a quasi-discrete power spectrum, in terms of a multi-mode truncation (Galerkin approximation) [14].This observation is specific for the harmonic (quadratic) confining potential, while anharmonic ones quickly lead tothe onset of clearly observed chaos [36, 37]. Thermalization of the model with the HO potential was recently explored,in the framework of a stochastically driven dissipative Gross-Pitaevskii equation, in Ref. [38].The interplay of the cubic nonlinearity and trapping potentials also occurs in two-component systems, whichrepresent, in particular, binary BEC [12, 15, 16, 18]. Here, we aim to consider the system with linear coupling betweenthe components. In optics, if two modes carrying orthogonal polarizations of light propagate in the same waveguide,linear mixing between them is induced by a twist of the guiding structure, see, e.g., Ref. [39]. In BEC, different atomic a r X i v : . [ n li n . PS ] F e b states which correspond to the interacting components can be linearly coupled by a resonant electromagnetic wave[40]-[43]. Another realization of linearly-mixed systems is offered by dual-core waveguides, coupled by tunneling ofthe field across a barrier separating the cores. The dual-core schemes are equally relevant to optics and BEC [44]. Inparticular, experiments with temporal solitons in dual-core optical fibers were reported in recent work [45].The coupling between the components enhances the complexity of the system and makes it possible to find newstatic and dynamical states in it. In particular, the symmetric system combining the attractive cubic terms of theSPM type and (optionally) HO potential acting in each component, with linear coupling between them, gives riseto spontaneous symmetry breaking (SSB) of two-component states [18, 45]. In the case of repulsive SPM actingin each component and nonlinear repulsion between them (cross-phase modulation, XPM), it was found [15] thatthe linear mixing shifts the miscibility-immiscibility transition [46] in the trapped condensate. Further, effects ofnonintegrability may be stronger in the linearly coupled system, because the linear coupling makes the system ofone-dimensional NLSEs nonintegrable, even in the absence of the HO potential [47], although the integrability iskept by the system with the linear coupling added to the SPM and XPM terms with equal strengths (the Manakov’snonlinearity) [48].The objective of the present work is two-fold. First, in Section II we aim to analyze the onset of chaos, as well asSSB, in the symmetric linearly-coupled system, with both attractive and repulsive signs of the SPM terms, startingfrom an input in the form of the ground state (GS), or the first excited state (the dipole mode (DM), represented bya spatially odd wave function) of the HO, which is initially launched in one core (component), while the other oneis empty. This type of the input is, in particular, experimentally relevant in optics [45]. As a nonchaotic dynamicalregime in this case, one may expect Josephson oscillations (JO) of the optical field [56]-[59], [47], [45], or of the BECwave function [49]-[55], between the cores. As a measure of the transition to dynamical chaos in the system, we use arelative spread of the power spectrum of oscillations produced by simulations of the coupled NLSEs. Naturally, thechaotization sets in above a certain threshold value of the nonlinearity strength, and the chaos is much weaker in thecase of the repulsive nonlinearity. In the dynamical state initiated by the GS, the SSB takes place prior to the onsetof the dynamical chaos, while the DM input undergoes the chaotization occurs first, followed by the SSB at a slightlystronger nonlinearity.The second objective, which is presented below in Section III, is to construct stationary GS and DM in theasymmetric ( half-trapped ) linearly-coupled system, with the HO potential applied to one component only. The lattersystem can be realized in the experiment, applying, for instance, the trapping potential only to one core of the doublewaveguide for matter waves in BEC (e.g., by focusing laser beams, which induced the trapping, on the single core).In optics, a similar setup may be built as a coupler with two widely different cores, narrow and broad ones, with thenarrow core emulating the component carrying a tightly confining potential, cf. Ref. [60]. A remarkable peculiarityof such a system is that the linear coupling mixes completely different types of the asymptotic behavior at | x | → ∞ :the trapped component is always confined, in the form of a Gaussian, by the HO potential, while the untrapped oneis free to escape. In addition, the asymmetry between the coupled cores makes it necessary to take into regard adifference in the chemical potentials or propagation constants (in terms of the BEC and optics, respectively) betweenthem, which is represented by parameters ω in Eq. (10), see below. Some results for the half-trapped system areobtained in an analytical form, in the weak- and strong-coupling limits, as well as by means of the Thomas-Fermiapproximation (TFA), and full results are produced numerically. In addition, particular solutions for localized statesof the linear half-trapped system (including vortex states in its 2D version) are found in an exact form. The exactsolutions belong to the class of bound states in the continuum [61–63], alias embedded [64] ones, which are spatiallyconfined modes existing, as exceptional states, with the carrier frequency falling in the continuous spectrum. Thepaper is concluded by Section IV. II. THE SYMMETRIC SYSTEMA. The coupled equations
As outlined above, we consider systems modeled by a pair of coupled NLSEs for complex amplitudes u ( x, z ) and v ( x, z ) of two interacting waves. In the normalized form, the equations are written in terms of the spatial-domainpropagation in an optical waveguide with propagation distance z and transverse coordinate x : i ∂u∂z + 12 ∂ u∂x + λv − Ω x u + σ (cid:0) | u | + g | v | (cid:1) u = 0 , (1) i ∂v∂z + 12 ∂ v∂x + λu − Ω x v + σ (cid:0) | v | + g | u | (cid:1) v = 0 . Here, coefficients σ > σ < λ and σg represent the linear mixing and XPM interaction, respectively. By means of rescaling, the strength of the HOtrapping potential is set to be Ω = 1 (unless it is zero in one core). In other words, x is measured in units of therespective HO length. This implies that the unit of the transverse coordinate in the optical waveguides takes typicalvalues in the range of 10 − µ m, hence the respective unit of the propagation distance (the Rayleigh/diffractiondistance corresponding to the OH length) is estimated to be between 1 mm and 1 cm, for the carrier wavelength ∼ µ m. In the matter-wave realization of the system, typical units of x and time (replacing z in Eq. (1) are ∼ µ mand 10 ms, respectively.The system conserves two dynamical invariants., viz ., the total norm (or power, in terms of optics), P = (cid:90) + ∞−∞ (cid:2) | u ( x ) | + | v ( x ) | (cid:3) dx, (2)and the Hamiltonian (in which Ω = 1 is set), H = (cid:90) + ∞−∞ (cid:34) (cid:32)(cid:12)(cid:12)(cid:12)(cid:12) ∂ u∂x (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) ∂ v∂x (cid:12)(cid:12)(cid:12)(cid:12) (cid:33) − σ (cid:0) | u | + | v | (cid:1) − σg | u | | v | + 12 x (cid:0) | u | + | v | (cid:1) − λ ( uv ∗ + u ∗ v ) (cid:21) dx, (3)where ∗ stands for the complex conjugate. The remaining scaling invariance of Eq. (1) makes it possible to either set | σ | = 1, or keep nonlinearity coefficient σ as a free parameter, but fix P ≡
1. All simulations performed in this workcomply with the conservation of P and H , up to the accuracy of the numerical codes.It is relevant to mention that the present two-component system resembles nonlinear models with a double-wellpotential, in the case when the wave functions in two wells are linearly coupled by tunneling across the potentialbarrier, see, e.g., Refs. [65–67]. Nevertheless, the exact form of the system and its solutions are different.In the limit of a small amplitude A of the input, linearized equations (1) with Ω = 1 admit exact solutions forinter-core JO of the ground and dipole states (exact solutions for higher-order states can be readily found too): u (GS)JO ( x, t ) = A exp (cid:18) − x − iz (cid:19) cos ( λz ) , (4) v (GS)JO ( x, t ) = iA exp (cid:18) − x − iz (cid:19) sin ( λz ) ,u (DM)JO ( x, t ) = A x exp (cid:18) − x − iz (cid:19) cos ( λz ) , (5) v (DM)JO ( x, t ) = iA x exp (cid:18) − x − iz (cid:19) sin ( λz ) . Expressions given by Eqs. (4) and (5) at z = 0 are used below as inputs in simulations of the full nonlinear equations(1). The simulations presented in this section were performed in domain | x | <
10. This size tantamount to 20 HOlengths is sufficient to display all details of the solutions. Standard numerical methods were used, viz ., the split-stepfast-Fourier-transform scheme for simulations of the evolution governed by Eqs. (1) and (10), and the relaxationalgorithm for finding solutions of stationary equations, such as Eqs. (12) and (13), see below.
B. The transition from regular to chaotic dynamics
Increase of amplitude A of the input leads to nonlinear deformation of the oscillations, and eventually to the onsetof dynamical chaos. A typical example of an essentially nonlinear but still regular JO dynamical regime, produced bynumerical simulations of Eq. (1) with g = 0 (no XPM interaction), in interval 0 < z < Z , with the DM input, takenas per Eq. (5) at z = 0, is presented in Fig. 1. In particular, the left bottom panel of the figure displays oscillationsof peak intensities of the fields, (cid:8) U ( z ) , V ( z ) (cid:9) ≡ max x (cid:110) | u ( x, z ) | , | v ( x, z ) | (cid:111) , (6) FIG. 1. A typical example of a regular Josephson dynamical regime, initiated by the DM (dipole-mode) input launched inthe u component (as given by Eq. (5) with z = 0 and A = 1). The solution is produced by simulations of Eq. (1) with λ = σ = Ω = 1, g = 0. Plots in the top row display the evolution of components u ( x, z ) and v ( x, z ). Left bottom: Theevolution of the peak intensities of both components, U ( z ) ≡ max x (cid:8) | u ( x, z ) | (cid:9) and V ( z ) ≡ max x (cid:8)(cid:2) | v ( x, z ) | (cid:3)(cid:9) . Rightbottom: The power spectrum of oscillations of the two components, defined as per Eq. (7). The spectra are virtually identicalfor both components. and, as a characteristic of the dynamics, in the right bottom panel we plot power spectra, | P ( κ ) | , | Q ( κ ) | , producedby the Fourier transform of the peak intensities: { P ( κ ) , Q ( κ ) } = (cid:90) Z (cid:8) U ( z ) , V ( z ) (cid:9) exp ( − iκz ) dz, (7)where κ is a real propagation constant. Very slow decay of the peak intensities, observed in the former panel, is amanifestation of the system’s nonintegrability (in this connection, we again stress that the total norm is conserved inthe course of the simulations).The regularity of the dynamical regime displayed in Fig. 1 is clearly demonstrated by its spectral structure, whichexhibits a single narrow peak at κ peak ≈
2. The peak’s width, ∆ κ (cid:39) .
1, which corresponds to the relative width,∆ κ/κ peak ≈ .
05, is comparable to the spread of the Fourier transform, corresponding to Z = 100 in Fig. 1. It canbe estimated as δκ = 2 π/Z ≈ .
06. Note the overall symmetry between the two components in Fig. 1 (in particular,their spectra are identical in the bottom panel).The simulations with the same input, but larger values of A , give rise to chaotic (“turbulent”) dynamical stateswith a broad dynamical spectrum, see a typical example in Fig. 2. Note that both the regular and chaotic dynamicalpictures displayed in Figs. 1 and 2 extend over the distance estimated to be ˜10 Rayleigh (diffraction) lengthscorresponding to the width of the DM input. This estimate is sufficient to make conclusions about the character ofthe dynamics.Systematic simulations with the GS input, provided by Eq. (4) at z = 0, produce similar results (not shown here indetail). In particular, as well as in the case of the DM input, amplitudes A = 1 and 4 initiate, severally, quasi-regularand chaotic evolution. FIG. 2. The same as in Fig. 1, but for a typical example of a chaotic dynamical regime, initiated by the DM input (5) with alarger amplitude, A = 4. Note different scales on vertical axes in two plots in the left panel. The results for the transition from regular JO to chaotic dynamics, initiated by the GS and DM inputs, taken as perEqs. (4) and (5) at z = 0, are summarized in charts plotted in the plane of (cid:0) λ, A (cid:1) in Fig. 3. They display heatmapsof values of the parameter which quantifies the sharpness of the central peak in the spectrum of the dynamical state:Sharpness ≡ (cid:82) FWHM | P ( κ ) | dκ (cid:82) ∞ | P ( κ ) | dκ , (8)where the integration in the numerator is performed over the section of the central spectral peak selected accordingto the standard definition of the full width at half-maximum: | P ( κ FWHM ) | = (1 /
2) ( P ( κ )) max . Values of Sharpnessclose to 1 imply the domination of a single sharp peak, such as one in Fig. 1, which corresponds to a regular dynamicalregime, while decrease of this parameter indicates a transition to a broad spectrum, which is a telltale of the onset ofchaotic dynamics – see, e.g., Fig. 2.Figure 3 clearly demonstrates decay of the central-peak’s sharpness, i.e., transition to dynamical chaos, with theincrease of the input’s intensity, A , in the case of the attractive SPM, σ = +1. Such a trend is not straightforwardin the opposite case of the self-repulsion, σ = −
1. In particular, the chaotization is not observed at all in the lattercase at small values of λ . This conclusion agrees with findings reported in work [14] for the single-component NLSE(which corresponds to λ = 0), with the HO potential and σ = −
1. Computations of the spectrum, reported in thatwork, demonstrate that no transition to dynamical chaos takes place at all values of parameters. The fact that anarea of weak chaotization is, nevertheless, observed in the right panels of Fig. 3 is explained by the above-mentionedcircumstance, that the addition of the linear coupling to a pair of NLSEs destroys their integrability (in free space).On the other hand, the increase of the sharpness with the increase of λ , observed at relatively small values of A (which is observed in an especially salient form in the left bottom panel of Fig. 3) also has a simple explanation:the increase of λ makes the linear terms in the system dominating over nonlinear ones, thus tending to maintain aquasi-linear behavior.Lastly, the comparison of the left and right panels in Fig. 3 suggests that the chaotization sets in faster in dynamicalregimes initiated by the DM input, in comparison to their counterparts originating from the GS, for the same values FIG. 3. Heatmaps of values of sharpness (8) of the central spectral peak quantifying proximity of the system’s dynamics to theregular regime. The maps are plotted in the plane of the linear-coupling strength, λ , and intensity of the input, A , which islaunched in one component. Top left: The GS input, given by Eq. (4) at z = 0, in the case of the self-attraction ( σ = +1). Topright: The same, but in the case of self-repulsion ( σ = − z = 0. Bottom right: The same as in the bottom left panel, but in the case of self-repulsion( σ = − g = 0 is set in Eq. (1) (no XPM interaction between the components). Black curves cutting the leftpanels in their lower areas designate the onset of SSB (spontaneous symmetry breaking), signalized by appearance of θ (cid:54) = 0,see Eq. (9). of the input’s intensity, A . The difference between the GS and DM dynamical regimes is salient for relatively smallvalues of λ . It may be explained by the fact that attractive SPM naturally tends to form a stable bright soliton fromthe GS input, which then maintains regular motion in the HO potential [19]. On the other hand, spatially odd brightsolitons do not exist in free space, which impedes transformation of the DM input into a regular dynamical state.As said above, the heatmaps are displayed in Fig. 3 for g = 0, i.e., in the absence of the XPM coupling between thecomponents, which is the case for dual-core couplers. On the other hand, in the case of the Manakov’s nonlinearity,i.e., g = 1, the above-mentioned integrability of such a system of NLSEs with the linear coupling [48] (but withoutthe trapping potential) suggests that the full system will be closer to integrability and farther from the onset ofchaos. Indeed, numerical results collected from simulations of Eq. (1) with σ = g = +1 (not shown here in detail)demonstrate a much smaller chaotic area in the (cid:0) λ, A (cid:1) plane. In particular, the GS input generates “turbulent”behavior only at A (cid:38) λ (cid:46) .
15, cf. Fig. 3(a).
C. Spontaneous symmetry breaking (SSB) between the coupled components
A noteworthy feature of the dynamical state presented in Fig. 2 is breaking of the symmetry between fields u and v (while the patterns initiated by the GS and DM inputs keep their parities, i.e., spatial symmetry (evenness) andantisymmetry (oddness), respectively). This is a manifestation of the general effect which is well known, in diverseforms, in linearly coupled dual-core systems with intrinsic attractive SPM [44]. In particular, SSB of stationary statesin systems with the HO trapping potential acting in both cores was addressed in Refs. [18] and [55], while Fig. 2 FIG. 4. Plots of the SSB (spontaneous symmetry breaking) in the dynamical states initiated by the GS (ground-state) andDM (dipole-mode) inputs: the asymmetry parameter, defined as per Eq. (9), is plotted versus the intensity of the input, A ,for three different fixed values of the linear-coupling constant, λ , as indicated in the figure. demonstrates the symmetry breaking in the dynamical JO state.The SSB effect may be quantified, as usual [47], by asymmetry of the dynamical states, which we define asΘ ≡ (cid:82) ∞ | P ( κ ) | dκ − (cid:82) ∞ | Q ( κ ) | dκ (cid:82) ∞ | P ( κ ) | dκ + (cid:82) ∞ | Q ( κ ) | dκ . (9)The SSB occurs as a transition from Θ = 0 to Θ (cid:54) = 0 with the increase of A at some critical point, which is a genericproperty of stationary states in dual-core systems with intrinsic attractive nonlinearity [44, 47], while here we considerit in the dynamical setting. On the other hand, in the case of the repulsive nonlinearity, σ = − g > g = 0.Accordingly, the present system with σ = − σ = +1, the SSB boundaries in the parameter planes of the GS and DM solutions are shown by bold blacklines in the top and bottom left panels of Figs. 3, respectively. The SSB bifurcation of the dynamical states underthe consideration is of the supercritical , alias forward , type [68], in terms of dependence Θ (cid:0) A (cid:1) , as shown in Fig.4 for the dynamical states initiated by the GS and DM inputs. It is observed that, naturally, the critical value of A increases with λ , as the symmetry is maintained by the linear coupling, hence stronger coupling needs strongernonlinearity to break the symmetry. Note also that the transition from Θ = 0 to Θ ≈ λ .It is worthy to note that, as clearly shown by the SSB boundary (the black line) in the top left panel of Fig. 3, theSSB of the GS mode happens prior to the transition to the dynamical chaos. This conclusion agrees with known resultsshowing that SSB of stationary modes does not, normally, lead to chaotization of the system’s dynamics [44, 47]. Onthe other hand, the bottom left panel in Fig. 3 demonstrates a different situation for the DM states, which exhibitthe chaotization prior to the SSB, although the separation between these transitions is small. This conclusion agreeswith the fact that, as clearly seen in Fig. 4, the SSB in DM states occurs at values of the input’s amplitude essentiallyhigher than those which determine the SSB threshold of the GS solutions. III. THE HALF-TRAPPED SYSTEM
The asymmetric system of linearly-coupled NLSEs, with the HO potential included in one equation only, is writtenas iu z + 12 u xx + λv − x u + σ (cid:0) | u | + g | v | (cid:1) u = − ωu, (10) iv z + 12 v xx + λu + σ (cid:0) | v | + g | u | (cid:1) v = 0 , cf. Eq. (1). Here, as said above, Ω = 1 is set in the first equation, and the propagation-constant mismatch, ω (interms of BEC, it represents a difference in the chemical potentials between the two wave functions) is a commonfeature of asymmetric systems. Stationary solutions to Eq. (10) are looked for as { u, v } = { U ( x ) , V ( x ) } exp ( − iµz ) , (11)with real propagation constant − µ (in BEC, with z replaced by t , µ is the chemical potential), and real functions U ( x ) and V ( x ) satisfying equations( µ + ω ) U + 12 d Udx + λV − x U + σ (cid:0) U + gV (cid:1) U = 0 , (12) µV + 12 d Vdx + λU + σ (cid:0) V + gU (cid:1) V = 0 . (13)Most results are produced below disregarding the XPM coupling between the components ( g = 0). Nevertheless, theXPM terms are included when collecting numerical results for the threshold of the existence of bound states, displayedin Fig. 14. A. The linearized system: analytical and numerical results
1. Emission of radiation in the untrapped component
In the linear limit, σ = 0, two decoupled equations in system (10) with λ = 0 produce completely different results:all excitations of component u stay confined in the HO trap, while the v component with any µ > − ω freely expands.In particular, in the limit of λ = 0, obvious bound-state GS and DM solutions of the linearized version of Eqs. (12)and (13), with zero v component, are U (0)GS ( x ) = 1 π / exp (cid:18) − x (cid:19) , V (0)GS = 0 , (14) U (0)DM ( x ) = √ π / x exp (cid:18) − x (cid:19) , V (0)DM = 0 , (15)where the pre-exponential constants are determined by the normalization condition, (cid:90) + ∞−∞ (cid:2) U ( x ) + V ( x ) (cid:3) dx = 1 , (16)which we adopt in this section. The eigenvalues corresponding to eigenmodes (14) and (15) are µ (0)GS = 1 / − ω ; µ (0)DM = 3 / − ω. (17)Proceeding to dynamical states, in the lowest approximation with respect to small λ the evolution of the v field isdriven by the respective linearized equation in system (10), iv z + 12 v xx = − λU (0)GS , DM ( x ) exp (cid:16) − iµ (0)GS , DM z (cid:17) . (18)Obviously, Eq. (18) gives rise to emission of propagating waves (“radiation”), in the form of v ∼ λ exp (cid:16) ikx − iµ (0)GS , DM z (cid:17) ,at resonant wavenumbers k = ± (cid:113) µ (0)GS , DM , provided that µ (0)GS , DM is positive, i.e., v rad ∼ λ exp (cid:18) ± i (cid:113) µ (0)GS , DM ( x − V ph z ) (cid:19) , (19)where the phase velocity is V ph = k ≡ ± (cid:113) µ (0)GS , DM (20)(in terms of the spatial-domain propagation in the optical waveguide, it is actually the beam’s slope). The expansionof the area in the ( x, z ) plane occupied by the radiation field is bounded by the group velocity, | V gr | = | k | ≡ V ph .An illustration of this dynamics is presented in Fig. 5. Straight red lines designate the wave-propagation directions,which exactly agree with the phase velocity predicted by Eq. (20), and the expansion of the area occupied by theradiation complies with the prediction based on the group velocity. X X
FIG. 5. Simulations of the evolution of the linearized half-trapped system (10), displayed in the untrapped component byplotting Re( v ( x, z )). Left: Emission of radiation generated by the GS (ground-state) populating the trapped component, u ( x, z ) (see Eq. (14)), with ω = 0 .
25 in Eq. (10). Right: The same, but for the radiation generated by the DM (dipole mode)in the trapped component (see Eq. (15)), with ω = 0. In both cases, the linear-coupling constant is λ = 0 . FIG. 6. Left: The evolution of fields u and v in the half-trapped system, as produced by simulations of Eq. (10) with ω = 0and coupling constant λ = 1, initiated by the DM input in the u core, taken as per Eq. (14) with A = 0 .
1. Here and inFig. 7, spurious left-right asymmetry of the radiation field in the v component is an illusion produced by plotting. Right: Theevolution of the peak intensities of both components, max x (cid:8) | u ( x, t ) | (cid:9) and max x (cid:8) | v ( x, t ) | (cid:9) . The emission of radiation into the v core gives rise to a gradual decay of the amplitude in the u core, due to theconservation of the total norm, see Eq. (16). An example of the decay is displayed in Fig. 6, for a small initialamplitude of the GS input, A = 0 . λ = 1, which makes the transfer of the norm (power) from u to v faster. On the other hand, thesame input with large A makes the u - v coupling a weak effect, in comparison with the dominant nonlinearity, hencethe input mode in the u core seems quite stable, as shown in Fig. 7.0 FIG. 7. The same as in Fig. 6, but for a large amplitude of the DM input in the u component, A = 8.
2. The shift of the GS and DM existence thresholds at small values of coupling constant λ At µ (0)GS , DM < λ >
0, threshold values ( ω GS , DM ) thr of the mismatchparameter ω in Eqs. (12) and (13), such that the bound states of the GS and DM types exist at ω > ( ω GS , DM ) thr , (21)respectively. Obviously, ( ω GS ) thr = 1 / ω DM ) thr = 3 / λ = 0, see Eq. (17).First, we aim to find lowest-order corrections to the GS and DM eigenvalues (17) for small λ . Then, a shift ofthe respective thresholds can be identified by setting µ GS , DM = 0. In the limit of λ = 0, the GS and DM wavefunctions are taken as per Eqs. (14) and (15), respectively. With small λ , the first-order solution for V ( x ), viz ., V ( x ) ≡ λV (1)GS , DM ( x ), has to be found from the inhomogeneous equation, that follows from Eq. (13), in which µ = 0is set: d dx V (1)GS , DM = − U (0)GS , DM ( x ) . (22)Straightforward integration of Eq. (22), with expressions (14) and (15) substituted on the right-hand side, yields V (1)GS ( x ) = −√ π / (cid:34) x √ (cid:18) x √ (cid:19) + (cid:114) π exp (cid:18) − x (cid:19)(cid:35) , (23) V (1)DM ( x ) = 2 π / erf (cid:18) x √ (cid:19) , (24)1 FIG. 8. The analytically predicted and numerically found threshold values of the mismatch parameter, ω , above which the GSand DM solutions (the left and right panels, respectively) are produced, for the half-trapped system, by the linearized versionof Eqs. (12) and (13), vs. coupling constant λ . The analytical results are produced by Eq. (27). For the GS, they are shown(in the inset) only for relatively small values of λ , as in this case the analytical approximation is inaccurate at larger λ . where erf( x ) is the standard error function, which is an odd function of x .Next, the small perturbation in Eq. (12), represented by term λV , produces a small shift δµ of the eigenvalue, asa feedback from component V . According to the standard rule of quantum mechanics, which deals with the linearSchr¨odinger equation [69], in the first approximation of the perturbation theory, when V ( x ) is replaced by expression(23) or (24), the result is δµ GS , DM = − λ (cid:90) + ∞−∞ V (1)GS , DM ( x ) U (0)GS , DM ( x ) dx ≡ − I GS , DM λ , (25)with coefficients I GS = − . , I DM = 4 , (26)where the former and latter ones are computed, respectively, in a numerical form and analytically (note opposite signs of these coefficients). Then, the accordingly shifted threshold values sought for are( ω GS ) thr = 1 / − I GS λ , ( ω DM ) thr = 3 / − I DM λ . (27)The analytical predictions are compared to numerical results, obtained from a solution of the linearized variant ofEqs. (12) and (13), in Fig. 8. Note that the linear u - v coupling facilitates the formation of the DM bound state, bylowering ( ω DM ) thr , but impedes to form the GS, by making the respective threshold, ( ω DM ) thr , higher. It is seen thatfor the DM the analytical approximation is essentially more accurate than for the GS. An example of a bound-statesolution of linearized equations (12) and (13) of the DM type, numerically found at λ = 0 .
225 and ω = 1 .
4, i.e., below the threshold value ( ω DM ) thr = 1 .
5, corresponding to the limit of λ →
0, is displayed in Fig. 9. The existence of theDM at this point agrees with the right panel of Fig. 8.
3. The analysis for large values of λ In the opposite limit of large coupling constant λ , an analytical approximation for the the discrete eigenvaluescan be developed too. In this case, | µ | may also be large, ∼ λ . In the zero-order approximation, one may neglectthe derivative term in Eq. (13), to obtain V ≈ − ( λ/µ ) U . Then, substituting this relation back into the originallyneglected derivative term, one obtains a necessary correction to this relation: V ≈ − λµ U + λ µ d Udx . (28)The subsequent substitution of this expression in the linearized equation (12) leads to an equation for U which istantamount to the usual stationary linear Schr¨odinger equation with the HO potential: (cid:18) µ − λ µ + ω (cid:19) U + 12 (cid:18) λ µ (cid:19) d Udx − x U = 0 . (29)2 FIG. 9. A bound state of the DM (dipole-mode) type in the half-trapped system, found as a numerical solution of Eqs. (12) and(13) with ω = 1 . λ = 0 . σ = 0 (the linearized version). The eigenvalue corresponding to this solution is µ ≈ − . Then, the standard solution for the quantum-mechanical HO yields an equation which determines the spectrum ofthe eigenvalues: µ − λ µ + ω = (cid:18)
12 + n (cid:19) (cid:115) λ µ , (30)where n = 0 , , , ... is the quantum number. Taking into regard that λ is now a large parameter, Eq. (30) producesa final result for the spectrum, µ ≈ − λ − ω √ (cid:18)
12 + n (cid:19) . (31)The spectrum remains equidistant in the current approximation, while further corrections ∼ /λ give rise to terms ∼ (1 / n ) , which break this property. As concerns the existence threshold for the bound states, Eq. (31) predicts ω thr ≈ − λ . The coefficient in this relation is not accurate, as the derivation is not valid for small | µ | , but theimplication is that, for large λ , ω thr drops to negative values with a large modulus, ∼ − λ .The prediction of the GS and DM eigenvalues, given by Eq. (31) with n = 0 and n = 1, respectively, is comparedto numerically found counterparts in Fig. 10, which shows proximity between the analytical and numerical results.The plots do not terminate in the displayed domain, i.e., they do not reach the existence boundary.
4. Exact solutions for one- and two-dimensional bound states in the continuum (BIC) in the linear system
A remarkable property of the coupled half-trapped system, represented by the linearized version of Eqs. (10) and(12), (13), is that it admits particular spatially-confined solutions in an exact analytical form. These are exceptionalsolutions, which, for an arbitrary value of the linear-coupling constant, λ , exist at a single, specially selected, value ofthe mismatch parameter, ω , and with a single value of the eigenvalue, µ . First, it is possible to find an exact modewhich is a fundamental one (GS) in the V component, and a second-order mode in U : U ( x ) = U (cid:20)(cid:18) λ − (cid:19) + x (cid:21) exp (cid:18) − x (cid:19) ,V ( x ) = − λU exp (cid:18) − x (cid:19) ,U = π − / (cid:0) λ + 4 λ + 1 / (cid:1) − ,ω = 94 − λ , µ = 12 (cid:18) λ + 12 (cid:19) , (32)3 FIG. 10. The top row: left and right panels display the analytically predicted eigenvalues given by Eq. (31) with n = 0 and n = 1, for the ground state (GS) and dipole mode (DM), respectively, of the half-trapped system, and their counterpartsproduced by the numerical solution of linearized coupled equations (12) and (13), as functions of the linear-coupling constant, λ , and mismatch parameter, ω . The bottom row: cross sections of the respective top panels along the diagonal connectingpoints ( λ, ω ) = (0 .
10) and (10 , λ . with amplitude U defined by condition P = 1, see Eq. (16). We stress that, as seen in Eq. (32), this exact solutionmay only have µ >
0, i.e., it is BIC (a bound state in the continuous spectrum [61–63]), alias an embedded mode [64].It is worthy to note that this BIC mode and additional ones, presented below, are found in the coupled system, withone component trapped in the HO potential.In addition to the above spatially-even solution, an odd one of the BIC type is available too. It is composed of aDM in the V component and a third-order mode in U . In the normalized form, i.e., with P = 1 (as per Eq. (16)), thesolution is U ( x ) = U x (cid:20)(cid:18) λ − (cid:19) + x (cid:21) exp (cid:18) − x (cid:19) ,V ( x ) = − λU x exp (cid:18) − x (cid:19) ,U = 2 π − / (cid:0) λ + 4 λ + 3 / (cid:1) − ,ω = 114 − λ , µ = 12 (cid:18) λ + 32 (cid:19) , (33)which also exists at a single value of ω , and with a single eigenvalue, µ >
0. Note that both exact solutions, given byEqs. (32) and (33), may exists at positive and negative values of ω , as well as at ω = 0 (in Eqs. (32) and (33), ω = 0at λ = 9 / λ = 11 /
2, respectively).These exact solutions for BIC states in the two-component systems are somewhat similar to those found in Ref. [70],which addressed a system of spin-orbit-coupled linear Gross-Pitaevskii equations for a binary BEC. In that work, exactsolutions were produced for a specially designed form of the trapping potential.4
FIG. 11. The evolution initiated by the exceptional (BIC) exact solution (32) of the system of linearized equations (12) and(13), as produced by simulations of the full nonlinear half-trapped system (10), with the attractive SPM, σ = 1 in the top row,and repulsive σ = − g = 0, λ = 7 / ω = 1 /
2, which are related as per Eq. (32).
The exact solutions of the linearized system may be tried as inputs in simulations of the full nonlinear system basedon Eq. (10), with ω selected as per the solutions. The simulations, performed with moderate values of the initialamplitude, produce a robust state with steady internal oscillations, as shown in Fig. 11 for the attractive ( σ = 1)and repulsive ( σ = −
1) SPM nonlinearity. In addition, the simulations demonstrate weak emission of radiation in theuntrapped ( v ) component. In the case of the self-attraction (the left panel in Fig. 7), the radiation is almost invisible.The self-repulsion in the v component, naturally, enhances the emission, which becomes visible in the right panel ofthe figure. Still larger amplitudes of the input lead to irregular oscillations and conspicuous emission of radiation (notshown here in detail).It may happen that the linearized system of Eqs. (12) and (13) produces isolated BIC/embedded modes in anumerical form at other values of parameters. The present work does not aim to carry our comprehensive search forsuch solutions. On the other hand, it is relevant to mention that a straightforward two-dimensional (2D) extension ofthe present system readily produces exceptional exact solutions for BIC/embedded modes of both fundamental (GS,alias zero-vorticity) and vortex types.The 2D extension of the linearized form of Eq. (10) is iu z + 12 ( u xx + u yy ) + λv − (cid:0) x + y (cid:1) u = − ωu, (34) iv z + 12 ( v xx + v yy ) + λu = 0 . While the realization of this model in optics in not straightforward, it may be implemented for matter waves in adual-core “pancake-shaped” holder of BEC [18, 71]. In polar coordinates ( r, θ ), particular exact solutions of linear5equations (34), with all integer values of the vorticity, S = 0 , , , ... , are found as u = U r S (cid:20) (cid:0) λ − − S (cid:1) + r (cid:21) exp (cid:18) − iµz + iSθ − r (cid:19) ,v = − λU r S exp (cid:18) − iµz + iSθ − r (cid:19) ,ω = 12 (cid:0) S − λ (cid:1) , µ = 12 (cid:0) λ + 1 + S (cid:1) , (35)with arbitrary amplitude U (the 2D solution of the GS type corresponds to S = 0 in Eq. (35)). As well as in1D solutions (32) and (33), µ takes only positive values in the 2D solution, hence it also represents states of theBIC/embedded type. B. The nonlinear half-trapped system
1. The Thomas-Fermi approximation (TFA)
In the presence of the self-defocusing nonlinearity, σ = − µ <
0, which is necessary for the existence of a generic (non-BIC) localized state in the V component. Then, Eq. (13)with g = 0 (the XPM coupling is omitted here) is solved as U = ( V /λ ) (cid:0) V − µ (cid:1) , (36)and Eq. (12) amounts to an algebraic equation for the squared amplitude W ≡ − V /µ > viz ., mW ( W + 1) + ξ W = ξ − ξ , (37)where m ≡ − µλ > , ξ ≡ − xµ , (38) ξ ≡ − µ (cid:2) λ − µ ( ω + µ ) (cid:3) , (39)and the applicability condition for TFA is easily shown to be m (cid:28)
1. The TFA solution exists under condition ξ > bandgap ) of values of the propagation constant,0 < − µ < (cid:112) λ + ω / ω/ ≡ − µ . (40)Outside of the bandgap, i.e., at µ < µ <
0, the TFA solution does not exist. In the bandgap, Eq. (37) produces ausual GS profile, with a single value of W corresponding to each ξ from region ξ < ξ (see Fig. 12, which displays atypical example of the TFA-predicted GS and its comparison with the numerically found counterpart). The solutionvanishes at the border points, ξ = ± ξ , and is equal to zero at ξ > ξ , so that the derivative of the TFA solution, dW/dξ , is discontinuous at the border points, which is a usual peculiarity of the TFA [2, 72].Note that, in the limit of large λ , the bandgap’s width, as given by Eq. (40), is − µ ≈ λ + ω/
2, which is close tothe largest value of − µ predicted by Eq. (31) for the GS ( n = 0) in the linearized half-trapped system. Finally, theTFA solution may be cast in a simple explicit form close to the edge of the bandgap, i.e., at0 < δµ ≡ µ − µ (cid:28) − µ , (41) ξ ≈ − (cid:0) /µ (cid:1) (cid:112) λ + ω / δµ (42)In this case, Eqs. (37) and (36) simplify to U ≈ (cid:26) (cid:0) µ / (cid:1) (cid:0) ξ − ξ (cid:1) , at ξ < ξ , ξ > ξ , (43) V ≈ (cid:26) (cid:0) λ / (cid:1) (cid:0) ξ − ξ (cid:1) , at ξ < ξ , ξ > ξ . (44)6 FIG. 12. A typical example of the GS solution predicted by TFA (Thomas-Fermi approximation) for the half-trapped system,as per Eqs. (36)-(39), for σ = − g = 0, λ = 8, µ = −
4, and its comparison to the numerically found counterpart. Therespective value of parameter m (see Eq. (38)), which should be small for the applicability of TFA, is m = 0 . Expressions (42)-(44) make it possible to calculate the total power of the GS (see Eq. (2)), P ≈ (cid:0) λ + µ (cid:1) (cid:18) λ + ω (cid:19) / ( − µ ) − / ( δµ ) / . (45)Note that relation (45) satisfies the anti-Vakhitov-Kolokolov criterion , dP/d ( δµ ) >
0, which is a necessary conditionfor stability of localized states supported by self-repulsive nonlinearity [73], as is the case in the present setting (theVakhitov-Kolokolov criterion proper, dP/d ( δµ ) <
0, is a well-known necessary condition for stability of states inmodels with self-attraction [74–76]).
2. Existence boundaries for nonlinear states
The shrinkage of the existence region for GS solutions, and its expansion for DM ones, in the linear version of thecoupled half-trapped system, shown in Fig. 8, suggests to identify existence boundaries of the same states in thefull nonlinear system. Numerical data, necessary for the delineation of the existence region of the GSs and DMs inthe nonlinear system, were collected by solving Eqs. (12) and (13) for spatially even and odd modes with the fixedtotal power, P = 1 (as per Eq. (16)), while the linear-coupling coefficient, λ , and the nonlinearity coefficient, σ ,were varied, the latter one taking both positive and negative values (for the self-attraction/repulsion). The resultsare summarized by the heatmap in Figs. 13 (for the GS) and 14 (for the DM), which show threshold values of themismatch parameter, defined as per Eq. (21).Nontrivial parametric areas for the GS and DM solutions are identified as domains with, respectively, ( ω GS ) thr < / ω DM ) thr < /
2. The former one is surrounded by red lines in Fig. 13, and its counterpart for the DMs is locatedbetween red lines in Fig. 14. In these areas, the coupled half-trapped system maintains stable localized GSs at ω < /
2, and DMs at ω < /
2, while in the absence of the coupling they may only exist at ω > / ω > / σ = 0 (along which the systemremains linear), that starts from λ = 0, belongs to the area with ( ω GS ) thr > /
2, in agreement with the result for thelinear system, which shows that the coupling impedes the existence of the GS, see Fig. 8(a) and Eq. (27). Further,Fig. 13 demonstrates that the SPM nonlinearity of either sign facilitates the creation of GS in parameter regionssurrounded by red lines. Naturally, the self-attraction ( σ >
0) helps to create such states starting from arbitrarily7
FIG. 13. The heatmap of threshold values of the mismatch parameter, ( ω GS ) thr , in the half-trapped system, based on Eqs. (12)and (13). For given values of the nonlinearity and linear-coupling coefficients, σ and λ , the stable GS (ground state), subject tothe normalization condition P = 1 (see Eq. (16)), exist above the threshold, i.e., at ω ≥ ( ω GS ) thr . Positive and negative valuesof σ correspond to the attractive and repulsive sign of the self-interaction, respectively. The nontrivial region is one confinedby red lines, in which the result is ( ω GS ) thr < /
2, i.e., the nonlinearity and linear coupling help to maintain self-trapped GSsin the parameter area where the decoupled system, with λ = 0, cannot create such states. FIG. 14. The same as in Fig. 13, but for DMs (dipole modes), obtained as numerical solutions to Eqs. (12) and (13) with g = 0and g = 1 (the left and right panels, respectively). In this case, the nontrivial region, located between the red boundaries, isone with ( ω DM ) thr < / small values of λ , while the self-repulsion ( σ <
0) is able to do it in the region separated by a gap from λ = 0.Nevertheless, at λ > .
24 the detrimental effect of the linear coupling cannot be outweighed by the SPM nonlinearity.In Fig. 14, the vertical cross section corresponding to the linear system ( σ = 0) entirely belongs to the nontrivialarea, with ( ω DM ) thr < /
2, also in agreement with Eq. (27) and Fig. 8(b). Figure 14 demonstrates that thenonlinearity, generally, impedes the maintenance of the localized DMs. This conclusion is supported, in particular, bythe comparison of panels (a) and (b) in the figure, which shows that the addition of the attractive XPM nonlinearitywith the to SPM conspicuously reduces the remaining nontrivial area.
IV. CONCLUSION
The objective of this work is to analyze new effects in the symmetric and asymmetric systems of linearly coupledfields, which are subject to the action of the HO (harmonic-oscillator) trapping potential and cubic self-attraction orrepulsion. The system can be implemented in nonlinear optics and BEC. In the symmetric system, with identical HOpotentials applied to both components, we focus on the consideration of JO (Josephson oscillations) in the system, bylaunching, in one component, an input in the form of the GS (ground state) or DM (dipole mode) of the HO potential.On the basis of systematically collected numerical data, we have identified two transitions in the system’s dynamics,which occur with the increase of the input’s power in the case of the self-attraction. The first is SSB (spontaneoussymmetry breaking) between the linearly coupled components in the dynamical JO state. At a higher power, thenonlinearity causes a transition from regular JO, initiated by the GS input, to chaotic dynamics. This transition isidentified through consideration of spectral characteristics of the dynamical regime. The input in the form of the8DM undergoes the chaotization at essentially smaller powers than the dynamical regime initiated by the GS input,which is followed by the SSB at slightly higher powers. In the case of self-repulsion, SSB does not occur, while thechaotization takes place in a weak form, in a small part of the parameter space.In the half-trapped system, with the HO potential acting on a single component, a nontrivial issue is identificationof the system’s linear spectrum, i.e., a parameter region in which the linearized system maintains trapped binary (two-component) modes. This problem is solved here analytically in the limit cases of weak and strong linear coupling,and in the numerical form in the general case. In particular, the linear coupling between the components leadsto the shrinkage of the spectral band in which the GS exists, and expansion of the existence band for the DM. Theexistence region for trapped states in the full nonlinear system is identified numerically, and such states are constructedanalytically by means of the TFA (Thomas-Fermi approximation). In addition, exceptional solutions of the linearizedsystem of the BIC (bound-state-in-continuum), alias embedded, type were found in the exact analytical form, in boththe 1D and 2D settings, the 2D solution being found with an arbitrary value of the vorticity.The work may be extended by considering inputs in the form of higher-order HO eigenstates. Another relevantdirection for the extension of the analysis is a systematic study of the 2D system.
ACKNOWLEDGMENT
This work was supported, in part, by the Israel Science Foundation through grant No. 1286/17. [1] C. J. Pethick and H. Smith,
Bose-Einstein Condensation in Dilute Gases (Cambridge University Press, Cambridge, 2002).[2] L. P. Pitaevskii and S. Stringari,
Bose-Einstein Condensation (Oxford University Press, Oxford 2003).[3] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonz´alez,
Emergent Nonlinear Phenomena in Bose-EinsteinCondensates: Theory and Experiment (Springer: Heidelberg, 2008).[4] B. I. Schneider and D. L. Feder, Numerical approach to the ground and excited states of a Bose-Einstein condensed gasconfined in a completely anisotropic trap, Phys. Rev. A , 2232-2242 (1999).[5] S. K. Adhikari, Numerical solution of the two-dimensional Gross-Pitaevskii equation for trapped interacting atoms, Phys.Lett. A , 91-96 (2000).[6] T. Busch and J. R. Anglin, Motion of dark solitons in trapped Bose-Einstein condensates, Phys. Rev. Lett. , 2298-2301(2000).[7] Y. S. Kivshar, T. J. Alexander, and S. K. Turitsyn, Nonlinear modes of a macroscopic quantum oscillator, Phys. Lett. , 225-230 (2001).[8] T. J. Alexander and L. Berg´e, Ground states and vortices of matter-wave condensates and optical guided waves, Phys.Rev. E , 026611 (2002).[9] G. Huang, J. Szeftel and S. Zhu, Dynamics of dark solitons in quasi-one-dimensional Bose-Einstein condensates, Phys.Rev. A , 053605 (2002).[10] N. G. Parker, N. P. Proukakis, M. Leadbeater, and C. S. Adams, Soliton-sound interactions in quasi-one-dimensionalBose-Einstein condensates, Phys. Rev. Lett. , 220401 (2003).[11] D. E. Pelinovsky, D. J. Frantzeskakis, and P. G. Kevrekidis, Oscillations of dark solitons in trapped Bose-Einstein conden-sates, Phys. Rev. E , 016615 (2005).[12] V. A. Brazhnyi and V. V. Konotop, Stable and unstable vector dark solitons of coupled nonlinear Schr¨odinger equations:Application to two-component Bose-Einstein condensates, Phys. Rev. E , 026616 (2005).[13] N. G. Parker, N. G. Proukakis, and C. S. Adams, Dark soliton decay due to trap anharmonicity in atomic Bose-Einsteincondensates, Phys. Rev. A , 033606 (2010).[14] T. Bland, N. G. Parker, N. P. Proukakis, and B. A. Malomed, Probing quasi-integrability of the Gross-Pitaevskii equationin a harmonic-oscillator potential, J. Phys. B: At. Mol. Opt. Phys. , 205303 (2018).[15] M. I. Merhasin, B. A. Malomed, and R. Driben, Transition to miscibility in a binary Bose-Einstein condensate induced bylinear coupling, J. Physics B , 877-892 (2005).[16] H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and R. Carretero-Gonzalez, Bright-dark solitoncomplexes in spinor Bose-Einstein condensates, Phys. Rev. A , 033612 (2008).[17] D.-S. Wang, X.-H. Hu, J. Hu, and W. M. Liu, Quantized quasi-two-dimensional Bose-Einstein condensates with spatiallymodulated nonlinearity, Phys. Rev. A , 025604 (2010).[18] Z. Chen, Y. Li, B. A. Malomed, and L. Salasnich, Spontaneous symmetry breaking of fundamental states, vortices, anddipoles in two and one-dimensional linearly coupled traps with cubic self-attraction, Phys. Rev. A , 033621 (2016).[19] Y. S. Kivshar and B. A, Malomed, Dynamics of solitons in nearly integrable systems, Rev. Mod. Phys. , 763-915 (1989).[20] S. Raghavan and G. P. Agrawal, Spatiotemporal solitons in inhomogeneous nonlinear media, Opt. Commun. , 377-382(2000). [21] D. A. Zezyulin, G. L. Alfimov, and V. V. Konotop, Nonlinear modes in a complex parabolic potential, Phys. Rev. A ,013606 (2010).[22] E. G. Charalampidis, P. G. Kevrekidis, D. J. Frantzeskakis, and B. A. Malomed, Dark-bright solitons in coupled NLSequations with unequal dispersion coefficients, Phys. Rev. E , 012924 (2015).[23] T. Mayteevarunyoo, B. A. Malomed, and D. V. Skryabin, One- and two-dimensional modes in the complex Ginzburg-Landau equation with a trapping potential, Opt. Exp. , 8849-8865 (2018).[24] M. Leib, F. Deppe, A. Marx, R. Gross, and M. J. Hartmann, Networks of nonlinear superconducting transmission lineresonators, New J. Phys. , 075024 (2012).[25] Y. S. Kivshar and G. P. Agrawal, Optical Solitons: From Fibers to Photonic Crystals (Academic Press: San Diego, 2003).[26] T. Dauxois, M. Peyrard,
Physics of Solitons (Cambridge University Press: Cambridge, 2006).[27] O. Morsch and M. Oberthaler, Dynamics of Bose-Einstein condensates in optical lattices, Rev. Modern Phys. , 179-215(2006).[28] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (PrincetonUniversity Press: Princeton, 2008).[29] M. Skorobogatiy, and J. Yang, Fundamentals of Photonic Crystal Guiding (Cambridge University Press: Cambridge, 2008).[30] E. A. Cerda-Mendez, D. Sarkar, D. N. Krizhanovskii, S. S. Gavrilov, K. Biermann, M. S. Skolnick, P. V. Santos, Exciton-polariton gap solitons in two-dimensional lattices, Phys. Rev. Lett. , 146401 (2013).[31] V. A. Brazhnyi and V. V. Konotop, Theory of nonlinear matter waves in optical lattices, Modern Phys. Lett. B , 627-651(2004).[32] H. Sakaguchi and B. A. Malomed, Dynamics of positive- and negative-mass solitons in optical lattices and inverted traps,J. Phys. B , 1443-1459 (2004).[33] V. E. Zakharov, S. V. Manakov, S. P. Novikov, and L. P. Pitaevskii, Theory of Solitons: Inverse Scattering Method (Nauka:Moscow, 1980; English translation: Consultants Bureau, New York, 1984).[34] V. Zakharov, F. Dias, and A. Pushkarev, One-dimensional wave turbulence, Phys. Rep. , 1-65 (2004).[35] I. E. Mazets and J. Schmiedmayer, Thermalization in a quasi-one-dimensional ultracold bosonic gas, New J. Phys. ,055023 (2010).[36] S. P. Cockburn, A. Negretti, N. P. Proukakis, and C. Henkel, Comparison between microscopic methods for finite-temperature Bose gases, Phys. Rev. A , 043619 (2011).[37] P. Grisins and I. E. Mazets, Thermalization in a one-dimensional integrable system, Phys. Rev. A , 053635 (2011).[38] K. F. Thomas, M. J. Davis, and K. V. Kheruntsyan, Thermalization of a quantum Newton’s cradle in a one-dimensionalquasicondensate, Phys. Rev. A , 023315 (2021).[39] M. Decker, M. Ruther, C. E, Kriegler, J. Zhou, C. M. Soukoulis, S. Linden, and M. Wegener, Strong optical activity fromtwisted-cross photonic metamaterials, Opt. Lett. , 2501-2503 (2009).[40] R. J. Ballagh, K. Burnett, and T. F. Scott, Theory of an output coupler for Bose-Einstein condensed atoms, Phys. Rev.Lett. , 1607-1611 (1997).[41] P. ¨Ohberg and S. Stenholm, Internal Josephson effect in trapped double condensates, Phys. Rev. A , 3890-3895 (1999).[42] D. T. Son and M. A. Stephanov, Domain walls of relative phase in two-component Bose-Einstein condensates, Phys. Rev.A , 063621 (2002).[43] S. D. Jenkins and T. A. B. Kennedy, Dynamic stability of dressed condensate mixtures, Phys. Rev. A , 053607 (2003).[44] B. A. Malomed, editor: Spontaneous Symmetry Breaking, Self-Trapping, and Josephson Oscillations , (Springer-Verlag:Berlin and Heidelberg, 2013).[45] N. V. Hung, L. X. T. Tai, I. Bugar, M. Longobucco, R. Buzcy´nski, B. A. Malomed, and M. Trippenbach, Reversibleultrafast soliton switching in dual-core highly nonlinear optical fibers, Opt. Lett. , 5221-5224 (2020).[46] V. P. Mineev, The theory of the solution of two near-ideal Bose gases, Zh. Eksp. Teor. Fiz. , 263-272 (1974) [Englishtranslation: Sov. Phys. - JETP , 132-136 (1974)].[47] B. A. Malomed, Solitons and nonlinear dynamics in dual-core optical fibers, In: Handbook of Optical Fibers (G.-D. Peng,Editor: Springer, 2018), pp. 421-474.[48] M. V. Tratnik and J. E. Sipe, Bound solitary waves in a birefringent optical fiber, Phys. Rev. A , 2011-2017 (1988).[49] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, Quantum dynamics of an atomic Bose-Einstein condensate in adouble-well potential, Phys. Rev. A , 4318 (1997).[50] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Quantum coherent atomic tunneling between two trapped Bose-Einstein condensates, Phys. Rev. Lett. , 4950 (1997).[51] M. Albiez, R. Gati, J. F¨olling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Direct observation of tunneling andnonlinear self-trapping in a single bosonic Josephson junction, Phys. Rev. Lett. , 010402 (2005).[52] Y. Shin, G.-B. Jo, M. Saba, T. A. Pasquini, W. Ketterle, and D. E. Pritchard, Optical weak link between two spatiallyseparated Bose-Einstein condensates, Phys. Rev. Lett. , 170402 (2005).[53] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, The a.c. and d.c. Josephson effects in a Bose–Einstein condensate,Nature , 579-583 (2007).[54] Z. Chen, Y. Li, and B. A. Malomed, Josephson oscillations of chirality and identity in two-dimensional solitons in spin-orbit-coupled condensates, Phys. Rev. Research , 033214 (2020).[55] H. Sakaguchi and B. A. Malomed, Symmetry breaking in a two-component system with repulsive interactions and linearcoupling, Comm. Nonlin. Sci. Num. Sim. , 105496 (2020).[56] C. Par´e and M. F(cid:32)lorja´nczyk, Approximate model of soliton dynamics in all-optical couplers, Phys. Rev. A , 6287 (1990). [57] A. I. Maimistov, Propagation of a light pulse in nonlinear tunnel-coupled optical waveguides, Kvant. Elektron. , 758(1991) [English translation: Sov. J. Quantum Electron. , 687 (1991)].[58] I. M. Uzunov, R. Muschall, M. Goelles, Y. S. Kivshar, B. A. Malomed, and F. Lederer, Pulse switching in nonlinear fiberdirectional couplers, Phys. Rev. E , 2527-2536 (1995).[59] M. Abbarchi, A. Amo, V. G. Sala, D. D. Solnyshkov, H. Flayac, L. Ferrier, I. Sagnes, E. Galopin, A. Lemaˆıtre, G. Malpuech,and J. Bloch, Macroscopic quantum self-trapping and Josephson oscillations of exciton polaritons, Nature Phys. , 275(2013).[60] N.-C. Panoiu, B. A. Malomed, and R. M. Osgood, Jr., Semidiscrete solitons in arrayed waveguide structures with Kerrnonlinearity, Phys. Rev. A , 013801 (2008).[61] F. H. Stillinger and D. R. Herrick, Bound states in continuum, Phys. Rev. A , 446-454 (1975).[62] A. Kodigala, T. Lepetit, Q. Gu, B. Bahari, Y. Fainman, and B. Kante, Lasing action from photonic bound states incontinuum, Nature
1, 196-199 (2017).[63] B. Midya and V. V. Konotop, Coherent-perfect-absorber and laser for bound states in a continuum, Opt. Lett. 43, 607-610(2018).[64] A. R. Champneys, B. A. Malomed, J. Yang, and D. J. Kaup, “Embedded solitons”: solitary waves in resonance with thelinear spectrum, Physica D , 340-354 (2001).[65] M. Albiez, R. Gati, J. Folling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Direct observation of tunneling andnonlinear self-trapping in a single bosonic Josephson junction, Phys. Rev. Leyy. 95, 010402 (2005).[66] D. Ananikian and T. Bergeman, Gross-Pitaevskii equation for Bose particles in a double-well potential: Two-mode modelsand beyond, Phys. Rev. A , 013604 (2006).[67] Y. V. Kartashov, V. V. Konotop, and V. A. Vysloukh, Dynamical suppression of tunneling and spin switching of aspin-orbit-coupled atom in a double-well trap, Phys. Rev. A , 063609 (2018).[68] G. Iooss and D. D. Joseph, Elementary Stability Bifurcation Theory (Springer, New York, 1980).[69] L. D. Landau and E. M. Lifshitz,