Nonlinear dynamics of wave packets in tunnel-coupled harmonic-oscillator traps

We consider a two-component linearly-coupled system with the intrinsic cubic nonlinearity and the harmonic-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 Josephson oscillations (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 JO state. Under still stronger nonlinearity, the regular JO initiated by the GS input carry over into a chaotic dynamical state. For the DM input, the chaotization happens at smaller powers than for the GS, which is followed by SSB at a slightly stronger nonlinearity. In the system with the defocusing nonlinearity, SSB does not take place, and dynamical chaos occurs in a small area of the parameter space. 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 spectrum is found analytically in the limits of weak and strong inter-component coupling, and numerically in the general case. Under the action of the coupling, the existence region of the confined modes shrinks for GSs and expands for DMs. In the full nonlinear system, the existence region for confined modes is identified in the numerical form. They are constructed too by means of the Thomas-Fermi approximation, in the case of the defocusing nonlinearity. Lastly, particular (non-generic) exact analytical solutions for confined modes, including vortices, in one- and two-dimensional asymmetric linearized systems are found. They represent bound states in the continuum.

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 of the field across a barrier separating the cores. The dual-core schemes are equally relevant to optics and BEC [44]. In particular, 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 new static and dynamical states in it. In particular, the symmetric system combining the attractive cubic terms of the SPM type and (optionally) HO potential acting in each component, with linear coupling between them, gives rise to spontaneous symmetry breaking (SSB) of two-component states [18,45]. In the case of repulsive SPM acting in each component and nonlinear repulsion between them (cross-phase modulation, XPM), it was found [15] that the linear mixing shifts the miscibility-immiscibility transition [46] in the trapped condensate. Further, effects of nonintegrability may be stronger in the linearly coupled system, because the linear coupling makes the system of one-dimensional NLSEs nonintegrable, even in the absence of the HO potential [47], although the integrability is kept by the system with the linear coupling added to the SPM and XPM terms with equal strengths (the Manakov's nonlinearity) [48].
The objective of the present work is two-fold. First, in Section II we aim to analyze the onset of chaos, as well as SSB, in the symmetric linearly-coupled system, with both attractive and repulsive signs of the SPM terms, starting from an input in the form of the ground state (GS), or the first excited state (the dipole mode (DM), represented by a spatially odd wave function) of the HO, which is initially launched in one core (component), while the other one is empty. This type of the input is, in particular, experimentally relevant in optics [45]. As a nonchaotic dynamical regime in this case, one may expect Josephson oscillations (JO) of the optical field [56]- [59], [47], [45], or of the BEC wave function [49]- [55], between the cores. As a measure of the transition to dynamical chaos in the system, we use a relative spread of the power spectrum of oscillations produced by simulations of the coupled NLSEs. Naturally, the chaotization sets in above a certain threshold value of the nonlinearity strength, and the chaos is much weaker in the case of the repulsive nonlinearity. In the dynamical state initiated by the GS, the SSB takes place prior to the onset of the dynamical chaos, while the DM input undergoes the chaotization occurs first, followed by the SSB at a slightly stronger nonlinearity.
The second objective, which is presented below in Section III, is to construct stationary GS and DM in the asymmetric (half-trapped ) linearly-coupled system, with the HO potential applied to one component only. The latter system can be realized in the experiment, applying, for instance, the trapping potential only to one core of the double waveguide 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 the narrow core emulating the component carrying a tightly confining potential, cf. Ref. [60]. A remarkable peculiarity of 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 one is free to escape. In addition, the asymmetry between the coupled cores makes it necessary to take into regard a difference in the chemical potentials or propagation constants (in terms of the BEC and optics, respectively) between them, which is represented by parameters ω in Eq. (10), see below. Some results for the half-trapped system are obtained in an analytical form, in the weak-and strong-coupling limits, as well as by means of the Thomas-Fermi approximation (TFA), and full results are produced numerically. In addition, particular solutions for localized states of the linear half-trapped system (including vortex states in its 2D version) are found in an exact form. The exact solutions belong to the class of bound states in the continuum [61][62][63], alias embedded [64] ones, which are spatially confined modes existing, as exceptional states, with the carrier frequency falling in the continuous spectrum. The paper is concluded by Section IV.

II. THE SYMMETRIC SYSTEM
A. 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-domain propagation in an optical waveguide with propagation distance z and transverse coordinate x: Here, coefficients σ > 0 or σ < 0 represents the strength of the focusing or defocusing SPM in each core, while λ and σg represent the linear mixing and XPM interaction, respectively. By means of rescaling, the strength of the HO trapping potential is set to be Ω = 1 (unless it is zero in one core). In other words, x is measured in units of the respective HO length. This implies that the unit of the transverse coordinate in the optical waveguides takes typical values in the range of 10 − 30 µm, hence the respective unit of the propagation distance (the Rayleigh/diffraction distance corresponding to the OH length) is estimated to be between 1 mm and 1 cm, for the carrier wavelength ∼ 1 µm. In the matter-wave realization of the system, typical units of x and time (replacing z in Eq. (1) are ∼ 10 µm and 10 ms, respectively. The system conserves two dynamical invariants., viz., the total norm (or power, in terms of optics), and the Hamiltonian (in which Ω = 1 is set), 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 work comply 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-well potential, in the case when the wave functions in two wells are linearly coupled by tunneling across the potential barrier, see, e.g., Refs. [65][66][67]. Nevertheless, the exact form of the system and its solutions are different.
In the limit of a small amplitude A 0 of the input, linearized equations (1) with Ω = 1 admit exact solutions for inter-core JO of the ground and dipole states (exact solutions for higher-order states can be readily found too): 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 HO lengths is sufficient to display all details of the solutions. Standard numerical methods were used, viz., the split-step fast-Fourier-transform scheme for simulations of the evolution governed by Eqs. (1) and (10), and the relaxation algorithm 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 0 of the input leads to nonlinear deformation of the oscillations, and eventually to the onset of dynamical chaos. A typical example of an essentially nonlinear but still regular JO dynamical regime, produced by numerical simulations of Eq. (1) with g = 0 (no XPM interaction), in interval 0 < z < Z, with the DM input, taken as per Eq. (5) at z = 0, is presented in Fig. 1. In particular, the left bottom panel of the figure displays oscillations of peak intensities of the fields, and, as a characteristic of the dynamics, in the right bottom panel we plot power spectra, |P (κ)| 2 , |Q(κ)| 2 , produced by the Fourier transform of the peak intensities: where κ is a real propagation constant. Very slow decay of the peak intensities, observed in the former panel, is a manifestation of the system's nonintegrability (in this connection, we again stress that the total norm is conserved in the course of the simulations). The regularity of the dynamical regime displayed in Fig. 1 is clearly demonstrated by its spectral structure, which exhibits a single narrow peak at κ peak ≈ 2. The peak's width, ∆κ 0.1, which corresponds to the relative width, ∆κ/κ peak ≈ 0.05, is comparable to the spread of the Fourier transform, corresponding to Z = 100 in Fig. 1. It can be estimated as δκ = 2π/Z ≈ 0.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 0 , give rise to chaotic ("turbulent") dynamical states with a broad dynamical spectrum, see a typical example in Fig. 2. Note that both the regular and chaotic dynamical pictures displayed in Figs. 1 and 2 extend over the distance estimated to be˜10 Rayleigh (diffraction) lengths corresponding to the width of the DM input. This estimate is sufficient to make conclusions about the character of the dynamics.
Systematic simulations with the GS input, provided by Eq. (4) at z = 0, produce similar results (not shown here in detail). In particular, as well as in the case of the DM input, amplitudes A 0 = 1 and 4 initiate, severally, quasi-regular and chaotic evolution. The results for the transition from regular JO to chaotic dynamics, initiated by the GS and DM inputs, taken as per Eqs. (4) and (5) at z = 0, are summarized in charts plotted in the plane of λ, A 2 0 in Fig. 3. They display heatmaps of values of the parameter which quantifies the sharpness of the central peak in the spectrum of the dynamical state: where the integration in the numerator is performed over the section of the central spectral peak selected according to the standard definition of the full width at half-maximum: |P (κ FWHM )| 2 = (1/2) (P (κ)) max . Values of Sharpness close to 1 imply the domination of a single sharp peak, such as one in Fig. 1, which corresponds to a regular dynamical regime, while decrease of this parameter indicates a transition to a broad spectrum, which is a telltale of the onset of chaotic dynamics -see, e.g., Fig. 2. Figure 3 clearly demonstrates decay of the central-peak's sharpness, i.e., transition to dynamical chaos, with the increase of the input's intensity, A 2 0 , in the case of the attractive SPM, σ = +1. Such a trend is not straightforward in the opposite case of the self-repulsion, σ = −1. In particular, the chaotization is not observed at all in the latter case 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 that work, demonstrate that no transition to dynamical chaos takes place at all values of parameters. The fact that an area of weak chaotization is, nevertheless, observed in the right panels of Fig. 3 is explained by the above-mentioned circumstance, 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 2 0 (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 a quasi-linear behavior.
Lastly, the comparison of the left and right panels in Fig. 3 suggests that the chaotization sets in faster in dynamical regimes initiated by the DM input, in comparison to their counterparts originating from the GS, for the same values of the input's intensity, A 2 0 . The difference between the GS and DM dynamical regimes is salient for relatively small values of λ. It may be explained by the fact that attractive SPM naturally tends to form a stable bright soliton from the GS input, which then maintains regular motion in the HO potential [19]. On the other hand, spatially odd bright solitons 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 the components, 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 without the trapping potential) suggests that the full system will be closer to integrability and farther from the onset of chaos. 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 λ, A 2 0 plane. In particular, the GS input generates "turbulent" behavior only at A 2 0 200, being limited to λ 0.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) and antisymmetry (oddness), respectively). This is a manifestation of the general effect which is well known, in diverse forms, in linearly coupled dual-core systems with intrinsic attractive SPM [44]. In particular, SSB of stationary states in systems with the HO trapping potential acting in both cores was addressed in Refs. [18] and [55], while  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 The SSB occurs as a transition from Θ = 0 to Θ = 0 with the increase of A 2 0 at some critical point, which is a generic property of stationary states in dual-core systems with intrinsic attractive nonlinearity [44,47], while here we consider it in the dynamical setting. On the other hand, in the case of the repulsive nonlinearity, σ = −1 in Eq. (1), SSB of stationary states takes place in this system only at g > 1 [55], while we here focus on the most relevant case of g = 0. Accordingly, the present system with σ = −1 does not feature SSB.
For σ = +1, the SSB boundaries in the parameter planes of the GS and DM solutions are shown by bold black lines in the top and bottom left panels of Figs. 3, respectively. The SSB bifurcation of the dynamical states under the consideration is of the supercritical, alias forward, type [68], in terms of dependence Θ A 2 0 , 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 2 0 increases with λ, as the symmetry is maintained by the linear coupling, hence stronger coupling needs stronger nonlinearity to break the symmetry. Note also that the transition from Θ = 0 to Θ ≈ 1 (a strongly asymmetric state) is steeper at larger λ.
It is worthy to note that, as clearly shown by the SSB boundary (the black line) in the top left panel of Fig. 3, the SSB of the GS mode happens prior to the transition to the dynamical chaos. This conclusion agrees with known results showing that SSB of stationary modes does not, normally, lead to chaotization of the system's dynamics [44,47]. On the other hand, the bottom left panel in Fig. 3 demonstrates a different situation for the DM states, which exhibit the chaotization prior to the SSB, although the separation between these transitions is small. This conclusion agrees with the fact that, as clearly seen in Fig. 4, the SSB in DM states occurs at values of the input's amplitude essentially higher 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 written as cf. Eq. (1). Here, as said above, Ω = 1 is set in the first equation, and the propagation-constant mismatch, ω (in terms of BEC, it represents a difference in the chemical potentials between the two wave functions) is a common feature of asymmetric systems. Stationary solutions to Eq. (10) are looked for as 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 + 1 2 µV + 1 2 Most results are produced below disregarding the XPM coupling between the components (g = 0). Nevertheless, the XPM terms are included when collecting numerical results for the threshold of the existence of bound states, displayed in Fig. 14.
A. The linearized system: analytical and numerical results

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 where the pre-exponential constants are determined by the normalization condition, which we adopt in this section. The eigenvalues corresponding to eigenmodes (14) and (15) are Proceeding to dynamical states, in the lowest approximation with respect to small λ the evolution of the v field is driven by the respective linearized equation in system (10), Obviously, Eq. (18) gives rise to emission of propagating waves ("radiation"), in the form of where the phase velocity is (in terms of the spatial-domain propagation in the optical waveguide, it is actually the beam's slope). The expansion of the area in the (x, z) plane occupied by the radiation field is bounded by the group velocity, |V gr | = |k| ≡ 2V 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 the radiation complies with the prediction based on the group velocity. The emission of radiation into the v core gives rise to a gradual decay of the amplitude in the u core, due to the conservation of the total norm, see Eq. (16). An example of the decay is displayed in Fig. 6, for a small initial amplitude of the GS input, A 0 = 0.1 (which corresponds to the quasi-linear dynamical regime), and a relatively large coupling constant, λ = 1, which makes the transfer of the norm (power) from u to v faster. On the other hand, the same input with large A 0 makes the u-v coupling a weak effect, in comparison with the dominant nonlinearity, hence the input mode in the u core seems quite stable, as shown in Fig. 7.  GS,DM < 0 the radiation is not generated by Eq. (19), hence two-component bound states may exist in this case. Our next objective is to find, for the coupled system with λ > 0, threshold values (ω GS,DM ) thr of the mismatch parameter ω in Eqs. (12) and (13), such that the bound states of the GS and DM types exist at ω > (ω GS,DM ) thr , respectively. Obviously, (ω GS ) thr = 1/2 and (ω DM ) thr = 3/2 in the limit of λ = 0, see Eq. (17). First, we aim to find lowest-order corrections to the GS and DM eigenvalues (17) for small λ. Then, a shift of the respective thresholds can be identified by setting µ GS,DM = 0. In the limit of λ = 0, the GS and DM wave functions are taken as per Eqs. (14) and (15), respectively. With small λ, the first-order solution for V (x), viz., GS,DM (x), has to be found from the inhomogeneous equation, that follows from Eq. (13), in which µ = 0 is set: Straightforward integration of Eq. (22), with expressions (14) and (15) substituted on the right-hand side, yields 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, as a feedback from component V . According to the standard rule of quantum mechanics, which deals with the linear Schrödinger equation [69], in the first approximation of the perturbation theory, when V (x) is replaced by expression (23) or (24), the result is with coefficients 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 The analytical predictions are compared to numerical results, obtained from a solution of the linearized variant of Eqs. (12) and (13), in Fig. 8. Note that the linear u-v coupling facilitates the formation of the DM bound state, by lowering (ω DM ) thr , but impedes to form the GS, by making the respective threshold, (ω DM ) thr , higher. It is seen that for the DM the analytical approximation is essentially more accurate than for the GS. An example of a bound-state solution 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 the DM at this point agrees with the right panel of Fig. 8.

The analysis for large values of λ
In the opposite limit of large coupling constant λ, an analytical approximation for the the discrete eigenvalues can be developed too. In this case, |µ| may also be large, ∼ λ. In the zero-order approximation, one may neglect the derivative term in Eq. (13), to obtain V ≈ − (λ/µ) U . Then, substituting this relation back into the originally neglected derivative term, one obtains a necessary correction to this relation: The subsequent substitution of this expression in the linearized equation (12) leads to an equation for U which is tantamount to the usual stationary linear Schrödinger equation with the HO potential: Then, the standard solution for the quantum-mechanical HO yields an equation which determines the spectrum of the eigenvalues: where n = 0, 1, 2, ... is the quantum number. Taking into regard that λ is now a large parameter, Eq. (30) produces a final result for the spectrum, The spectrum remains equidistant in the current approximation, while further corrections ∼ 1/λ give rise to terms ∼ (1/2 + n) 2 , which break this property. As concerns the existence threshold for the bound states, Eq. (31) predicts ω thr ≈ −2λ. The coefficient in this relation is not accurate, as the derivation is not valid for small |µ|, but the implication 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 compared to 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.

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 exceptional solutions, which, for an arbitrary value of the linear-coupling constant, λ, exist at a single, specially selected, value of the mismatch parameter, ω, and with a single value of the eigenvalue, µ. First, it is possible to find an exact mode which is a fundamental one (GS) in the V component, and a second-order mode in U : with amplitude U 0 defined by condition P = 1, see Eq. (16). We stress that, as seen in Eq. (32), this exact solution may only have µ > 0, i.e., it is BIC (a bound state in the continuous spectrum [61][62][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, with one 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 a DM in the V component and a third-order mode in U . In the normalized form, i.e., with P = 1 (as per Eq. (16)), the solution is which also exists at a single value of ω, and with a single eigenvalue, µ > 0. Note that both exact solutions, given by Eqs. (32) and (33), may exists at positive and negative values of ω, as well as at ω = 0 (in Eqs. (32) and (33), ω = 0 at λ 2 = 9/2 and λ 2 = 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, exact solutions were produced for a specially designed form of the trapping potential. The exact solutions of the linearized system may be tried as inputs in simulations of the full nonlinear system based on Eq. (10), with ω selected as per the solutions. The simulations, performed with moderate values of the initial amplitude, 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 the untrapped (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 of the figure. Still larger amplitudes of the input lead to irregular oscillations and conspicuous emission of radiation (not shown here in detail).
It may happen that the linearized system of Eqs. (12) and (13) produces isolated BIC/embedded modes in a numerical form at other values of parameters. The present work does not aim to carry our comprehensive search for such solutions. On the other hand, it is relevant to mention that a straightforward two-dimensional (2D) extension of the present system readily produces exceptional exact solutions for BIC/embedded modes of both fundamental (GS, alias zero-vorticity) and vortex types.
While the realization of this model in optics in not straightforward, it may be implemented for matter waves in a dual-core "pancake-shaped" holder of BEC [18,71]. In polar coordinates (r, θ), particular exact solutions of linear equations (34), with all integer values of the vorticity, S = 0, 1, 2, ..., are found as with arbitrary amplitude U 0 (the 2D solution of the GS type corresponds to S = 0 in Eq. (35)). As well as in 1D solutions (32) and (33), µ takes only positive values in the 2D solution, hence it also represents states of the BIC/embedded type.
B. The nonlinear half-trapped system

The Thomas-Fermi approximation (TFA)
In the presence of the self-defocusing nonlinearity, σ = −1 in Eqs. (12) and (13), it is relevant to apply TFA to finding the GS of the half-trapped system, omitting the second derivatives in both equations [2], and keeping condition µ < 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 and Eq. (12) amounts to an algebraic equation for the squared amplitude W ≡ −V 2 /µ > 0, viz., where and the applicability condition for TFA is easily shown to be m 1. The TFA solution exists under condition ξ 2 0 > 0 (see Eq. (39)), i.e., in a finite interval (bandgap) of values of the propagation constant, Outside of the bandgap, i.e., at µ < µ 0 < 0, the TFA solution does not exist. In the bandgap, Eq. (37) produces a usual GS profile, with a single value of W corresponding to each ξ 2 from region ξ 2 < ξ 2 0 (see Fig. 12, which displays a typical example of the TFA-predicted GS and its comparison with the numerically found counterpart). The solution vanishes at the border points, ξ = ±ξ 0 , and is equal to zero at ξ 2 > ξ 2 0 , 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 −µ 0 ≈ λ + ω/2, which is close to the largest value of −µ predicted by Eq. (31) for the GS (n = 0) in the linearized half-trapped system. Finally, the TFA solution may be cast in a simple explicit form close to the edge of the bandgap, i.e., at In this case, Eqs. (37) and (36) simplify to (2)), Note that relation (45) satisfies the anti-Vakhitov-Kolokolov criterion, dP/d (δµ) > 0, which is a necessary condition for stability of localized states supported by self-repulsive nonlinearity [73], as is the case in the present setting (the Vakhitov-Kolokolov criterion proper, dP/d (δµ) < 0, is a well-known necessary condition for stability of states in models with self-attraction [74][75][76]).

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 the coupled half-trapped system, shown in Fig. 8, suggests to identify existence boundaries of the same states in the full nonlinear system. Numerical data, necessary for the delineation of the existence region of the GSs and DMs in the nonlinear system, were collected by solving Eqs. (12) and (13) for spatially even and odd modes with the fixed total 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 results are summarized by the heatmap in Figs. 13 (for the GS) and 14 (for the DM), which show threshold values of the mismatch parameter, defined as per Eq. (21).
Nontrivial parametric areas for the GS and DM solutions are identified as domains with, respectively, (ω GS ) thr < 1/2 and (ω DM ) thr < 3/2. The former one is surrounded by red lines in Fig. 13, and its counterpart for the DMs is located between red lines in Fig. 14. In these areas, the coupled half-trapped system maintains stable localized GSs at ω < 1/2, and DMs at ω < 3/2, while in the absence of the coupling they may only exist at ω > 1/2 and ω > 3/2, respectively.
Note that the vertical cross-section of the heatmap in Fig. 13, drawn through σ = 0 (along which the system remains linear), that starts from λ = 0, belongs to the area with (ω GS ) thr > 1/2, in agreement with the result for the linear 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 regions surrounded by red lines. Naturally, the self-attraction (σ > 0) helps to create such states starting from arbitrarily  12) and (13). For given values of the nonlinearity and linear-coupling coefficients, σ and λ, the stable GS (ground state), subject to the normalization condition P = 1 (see Eq. (16)), exist above the threshold, i.e., at ω ≥ (ωGS) thr . Positive and negative values of σ correspond to the attractive and repulsive sign of the self-interaction, respectively. The nontrivial region is one confined by red lines, in which the result is (ωGS) thr < 1/2, i.e., the nonlinearity and linear coupling help to maintain self-trapped GSs in 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 = 0 and g = 1 (the left and right panels, respectively). In this case, the nontrivial region, located between the red boundaries, is one with (ωDM) thr < 3/2. small values of λ, while the self-repulsion (σ < 0) is able to do it in the region separated by a gap from λ = 0. Nevertheless, at λ > 0.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 nontrivial area, with (ω DM ) thr < 3/2, also in agreement with Eq. (27) and Fig. 8(b). Figure 14 demonstrates that the nonlinearity, generally, impedes the maintenance of the localized DMs. This conclusion is supported, in particular, by the comparison of panels (a) and (b) in the figure, which shows that the addition of the attractive XPM nonlinearity with 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 coupled fields, which are subject to the action of the HO (harmonic-oscillator) trapping potential and cubic self-attraction or repulsion. The system can be implemented in nonlinear optics and BEC. In the symmetric system, with identical HO potentials applied to both components, we focus on the consideration of JO (Josephson oscillations) in the system, by launching, 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 (spontaneous symmetry breaking) between the linearly coupled components in the dynamical JO state. At a higher power, the nonlinearity causes a transition from regular JO, initiated by the GS input, to chaotic dynamics. This transition is identified through consideration of spectral characteristics of the dynamical regime. The input in the form of the DM 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 the chaotization 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 identification of the system's linear spectrum, i.e., a parameter region in which the linearized system maintains trapped binary (twocomponent) 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 leads to the shrinkage of the spectral band in which the GS exists, and expansion of the existence band for the DM. The existence region for trapped states in the full nonlinear system is identified numerically, and such states are constructed analytically by means of the TFA (Thomas-Fermi approximation). In addition, exceptional solutions of the linearized system of the BIC (bound-state-in-continuum), alias embedded, type were found in the exact analytical form, in both the 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 relevant direction for the extension of the analysis is a systematic study of the 2D system.