**Analysis of a Trapped Bose–Einstein Condensate in Terms of Position, Momentum, and Angular-Momentum Variance**

#### **Ofir E. Alon 1,2**


Received: 30 September 2019; Accepted: 27 October 2019; Published: 1 November 2019

**Abstract:** We analyze, analytically and numerically, the position, momentum, and in particular the angular-momentum variance of a Bose–Einstein condensate (BEC) trapped in a two-dimensional anisotropic trap for static and dynamic scenarios. Explicitly, we study the ground state of the anisotropic harmonic-interaction model in two spatial dimensions analytically and the out-of-equilibrium dynamics of repulsive bosons in tilted two-dimensional annuli numerically accurately by using the multiconfigurational time-dependent Hartree for bosons method. The differences between the variances at the mean-field level, which are attributed to the shape of the BEC, and the variances at the many-body level, which incorporate depletion, are used to characterize position, momentum, and angular-momentum correlations in the BEC for finite systems and at the limit of an infinite number of particles where the bosons are 100% condensed. Finally, we also explore inter-connections between the variances.

**Keywords:** Bose–Einstein condensates; density; position variance; momentum variance; angular-momentum variance; harmonic-interaction model; MCTDHB

**PACS:** 03.75.Hh; 03.75.Kk; 67.85.Bc; 67.85.De; 03.65.-w

#### **1. Introduction**

Bose–Einstein condensates (BECs) made of ultra-cold atoms offer a wide platform to study many-body physics [1–5]. Here, there is a growing interest in the so-called particle limit [6–16], in which the interaction parameter (i.e., the product of the interaction strength times the number of particles) is kept fixed while the number of particles is increased to infinity. At the particle limit, the energy per particle, density per particle, and reduced density matrices [17] per particle computed at the many-body level of theory boil down to those obtained in mean-field theory [7–10,14,16], despite the fact that the respective many-boson wavefunctions are (much) different [13,15]. It turns out that variances of many-particle operators are a useful tool to characterize correlations (namely, differences between respective many-body and mean-field quantities) that exist even when the interacting bosons are 100% condensed [11,12].

The variance of a many-particle operator of a trapped BEC generally depends on the trap shape, strength and sign of the interaction and, in out-of-equilibrium problems, on time. Consequently, the difference between variances computed at the many-body and mean-field levels of theory also depends on these variables and, of course, on the observable under examination. The first examples [11,12] concentrated on one-dimensional problems and the position and momentum variances, and investigated conditions and mechanisms for the differences between the respective many-body and mean-field variances at the particle limit. In two spatial dimensions, further types of trap topologies come into play, and respective many-body and mean-field variances can exhibit

additional phenomena, such as opposite anisotropy [18] and distinct (effective) dimensionality [19]. The many-body variance of a trap BEC has been applied to extract excitations [20], analyze the range of inter-particle interaction [21], examine the effects of asymmetry of a double-well potential [22], and to assess numerical convergence [23,24].

So far, only the position and momentum variances were studied for BECs in rather general traps. In [25,26], the angular-momentum variance is studied for BECs in two-dimensional isotropic traps, and scenarios were the mean-field angular-momentum variance has less [25] or more [26] symmetry (in terms of its conservation) than the many-body angular-momentum variance are identified. Going beyond these works, in the present work we study, analytically and numerically, the angular-momentum variance of a trapped BEC in a two-dimensional anisotropic trap for static and dynamic scenarios, and analyze the difference between the many-body and mean-field variances for finite systems and at the limit of an infinite number of particles. Furthermore, we also study the respective position and momentum variances, and thereby offer a comprehensive characterization of the BEC in terms of its variances. This would allow us to put forward inter-connections between the variances.

Let us elaborate on the strategy of exposition chosen in the paper. We first study the ground state of a many-particle model which is exactly solvable, i.e., integrable, both at the many-body and mean-field levels of theory. A couple of symmetries are also used in the analysis. These would allow us to obtain exact and transparent results for any number of particles and particularly to analyze the variances at the particle limit. The merit of analytical closed-form results and, in the context of interacting bosons, their explicit evaluation at the limit of an infinite number of particles is obvious. Then, as is the usual case in many realistic systems, we continue to explore a set-up which is not integrable, and more so, examine its out-of-equilibrium dynamics which is rather complicated already at the mean-field level of theory, let alone at the many-body level of theory. The later necessitates the state-of-the art numerical tools for the accurate integration of the Schrödinger equation and a careful interpolation of properties to the particle limit. All in all, we show below that the combination of analytics and numerics, i.e., of completely opposite methodologies, provides substantial and complementary novel knowledge on the position, momentum, and angular-momentum variances of anisotropic trapped BECs in two spatial dimensions.

The structure of the paper is as follows. In Section 2 we study the position, momentum, and angular-momentum variances of the ground state within an exactly solvable model, the anisotropic harmonic-interaction model. In Section 3 we study numerically the time-dependent variances of an out-of-equilibrium BEC sloshing in a tilted annulus. Summary and outlook are given in Section 4. Finally, Appendix A discusses translations of variances and inter-connections of the latter.

#### **2. The Anisotropic Harmonic-Interaction Model**

Solvable models of particles interacting by harmonic forces, or, briefly, the harmonic-interaction model (and its variants), have drawn including in the BEC literature much attention [27–41]. Here we consider the anisotropic two-dimensional harmonic-interaction model

$$\hat{H}(\mathbf{r}\_{1},\ldots,\mathbf{r}\_{N}) = \sum\_{j=1}^{N} \left[ \left( -\frac{1}{2}\frac{\partial^{2}}{\partial x\_{j}^{2}} + \frac{1}{2}\omega\_{x}^{2}x\_{j}^{2} \right) + \left( -\frac{1}{2}\frac{\partial^{2}}{\partial y\_{j}^{2}} + \frac{1}{2}\omega\_{y}^{2}y\_{j}^{2} \right) \right]$$

$$+\lambda\_{0}\sum\_{1\le j$$

where *λ*<sup>0</sup> is the interaction strength; positive values imply attraction and negative repulsion. Without loss of generality we take *ω<sup>y</sup>* > *ωx*, namely, that the trap is tighter along the y-axis than along the x-axis (the trap anisotropy satisfies *<sup>ω</sup><sup>y</sup> ωx* > 1). Here and hereafter ¯*h* = *m* = 1.

Transforming from Cartesian to Jacobi coordinates,

$$\begin{aligned} Q\_{k,x} &= \frac{1}{\sqrt{k(k+1)}} \sum\_{j=1}^{k} \left( \mathbf{x}\_{k+1} - \mathbf{x}\_{j} \right), \quad Q\_{k,y} = \frac{1}{\sqrt{k(k+1)}} \sum\_{j=1}^{k} \left( y\_{k+1} - y\_{j} \right), \quad 1 \le k \le N-1, \\\ Q\_{N,x} &= \frac{1}{\sqrt{N}} \sum\_{j=1}^{N} \mathbf{x}\_{j}, \quad Q\_{N,y} = \frac{1}{\sqrt{N}} \sum\_{j=1}^{N} y\_{j} \end{aligned} \tag{2}$$

the many-body solution for the ground state is given by

$$\Psi(\mathbf{Q}\_1, \dots, \mathbf{Q}\_N) = \left(\frac{\omega\_X}{\pi}\right)^{\frac{1}{2}} \left(\frac{\omega\_Y}{\pi}\right)^{\frac{1}{4}} \left(\frac{\Omega\_x}{\pi}\right)^{\frac{N-1}{4}} \left(\frac{\Omega\_y}{\pi}\right)^{\frac{N-1}{4}}$$

$$\times e^{-\frac{1}{2}\left(\Omega\_x \sum\_{j=1}^{N-1} Q\_{jx}^2 + \omega\_X Q\_{N,x}^2\right)} \times e^{-\frac{1}{2}\left(\Omega\_y \sum\_{j=1}^{N-1} Q\_{jy}^2 + \omega\_Y Q\_{N,y}^2\right)}$$

$$= \Psi(\mathbf{r}\_1, \dots, \mathbf{r}\_N) = \left(\frac{\omega\_X}{\pi}\right)^{\frac{1}{2}} \left(\frac{\omega\_Y}{\pi}\right)^{\frac{1}{4}} \left(\frac{\Omega\_x}{\pi}\right)^{\frac{N-1}{4}} \left(\frac{\Omega\_y}{\pi}\right)^{\frac{N-1}{4}}$$

$$\times e^{-\frac{\frac{\omega\_Y}{\pi} \sum\_{j=1}^{N} x\_j^2 - \beta x\_j \sum\_{1 \le j < k}^{N} x\_j x\_k} \times e^{-\frac{\omega\_Y}{\pi} \sum\_{j=1}^{N} y\_j^2 - \beta y \sum\_{1 \le j < k}^{N} y\_j y\_k} \tag{3}$$

where

$$
\Omega\_{\mathbf{x}} = \sqrt{\omega\_{\mathbf{x}}^2 + 2N\lambda\_0} \qquad \Omega\_{\mathbf{y}} = \sqrt{\omega\_{\mathbf{y}}^2 + 2N\lambda\_0} \tag{4}
$$

are the interaction-dressed frequencies of the relative-motion degrees-of-freedom, and

$$\begin{aligned} \mathfrak{a}\_{\mathbf{x}} &= \Omega\_{\mathbf{x}} + \mathfrak{f}\_{\mathbf{x}\prime} & \mathfrak{f}\_{\mathbf{x}} &= \frac{1}{N} \left( \omega\_{\mathbf{x}} - \Omega\_{\mathbf{x}} \right) \\ \mathfrak{a}\_{\mathbf{y}} &= \Omega\_{\mathbf{y}} + \mathfrak{f}\_{\mathbf{y}\prime} & \mathfrak{f}\_{\mathbf{y}} &= \frac{1}{N} \left( \omega\_{\mathbf{y}} - \Omega\_{\mathbf{y}} \right) \end{aligned} \tag{5}$$

are parameters arising in the transformation from Jacoby coordinates back to Cartesian coordinates. Equation (4) prescribes the range of interactions for which the system is trapped, *λ*<sup>0</sup> > − *ω*2 *x* 2*N* , i.e., from moderate repulsion to any attraction. Clearly, the many-body solution (3) in two spatial dimensions factorizes to a product of respective one-dimensional many-body solutions.

All properties of the ground state can in principle be obtained from Ψ, such as the energy, densities, and reduced density matrices, see [29]. Here, as mentioned above, we concentrate on variances and their inter-connections. The many-particle position *X*ˆ = ∑ *N j*=1 *xj* , *Y*ˆ = ∑ *N j*=1 *y<sup>j</sup>* variance per particle is given by

$$\frac{1}{N}\Delta\_{\hat{X}}^{2} = \frac{1}{2\omega\_{\chi}}, \qquad \frac{1}{N}\Delta\_{\hat{Y}}^{2} = \frac{1}{2\omega\_{y}}.\tag{6}$$

Due to the symmetry of center-of-mass separation in the Hamiltonian (1), the many-particle position variance per particle is independent both of the interaction strength and the number of bosons in the system. Similarly, the many-particle momentum *P*ˆ*<sup>X</sup>* = ∑ *N j*=1 1 *i ∂ ∂x<sup>j</sup>* , *P*ˆ *<sup>Y</sup>* = ∑ *N j*=1 1 *i ∂ ∂y<sup>j</sup>* variance per particle is given by

$$
\frac{1}{N} \Delta\_{\mathbb{P}\_{\mathbb{X}}}^2 = \frac{\omega\_{\mathbb{X}}}{2}, \qquad \frac{1}{N} \Delta\_{\mathbb{P}\_{\mathbb{Y}}}^2 = \frac{\omega\_{\mathbb{Y}}}{2}. \tag{7}
$$

reflecting the minimal uncertainty product <sup>1</sup> *N* ∆ 2 *X*ˆ 1 *N* ∆ 2 *P*ˆ*X* = <sup>1</sup> *N* ∆ 2 *Y*ˆ 1 *N* ∆ 2 *P*ˆ *Y* = <sup>1</sup> 4 of the interacting system in the anisotropic harmonic trap.

The many-particle angular-momentum *L*ˆ *<sup>Z</sup>* = ∑ *N j*=1 1 *i xj ∂ ∂y<sup>j</sup>* − *y<sup>j</sup> ∂ ∂x<sup>j</sup>* variance per particle is, at least for bosons, a less familiar and more intricate quantity. After some lengthy but otherwise straightforward algebra it is given by

$$\frac{1}{N}\Delta\_{L}^{2} = \frac{1}{4} \frac{\left(\Omega\_{\text{y}} - \Omega\_{\text{x}}\right)^{2}}{\Omega\_{\text{y}}\Omega\_{\text{x}}} \left(\frac{N-1}{N}\right)^{2} \left[\left(1 + \frac{1}{N-1}\frac{\Omega\_{\text{y}}}{\omega\_{\text{y}}}\right) \left(1 + \frac{1}{N-1}\frac{\Omega\_{\text{x}}}{\omega\_{\text{x}}}\right) \right. \\ \tag{8}$$

$$+ \left(\frac{\Omega\_{\text{y}}}{\omega\_{\text{y}}} - 1\right) \left(\frac{\Omega\_{\text{x}}}{\omega\_{\text{x}}} - 1\right) \right] + \frac{1}{4N} \frac{\left[\left(\omega\_{\text{y}} - \Omega\_{\text{y}}\right) - \left(\omega\_{\text{x}} - \Omega\_{\text{x}}\right)\right] \left[\left(\omega\_{\text{y}} + \Omega\_{\text{y}}\right) - \left(\omega\_{\text{x}} + \Omega\_{\text{x}}\right)\right]}{\omega\_{\text{y}}\omega\_{\text{x}}}$$

where we have made use of the bosonic permutational symmetry, the structure of Ψ,

$$\hat{L}\_Z \Psi = -\frac{1}{i} \left\{ \left[ (\mathbf{a}\_y - \beta\_y) - (\mathbf{a}\_x - \beta\_x) \right] \left( \sum\_{j=1}^N \mathbf{x}\_j y\_j \right) + (\beta\_y - \beta\_x) \left( \sum\_{j=1}^N \mathbf{x}\_j \right) \left( \sum\_{k=1}^N y\_k \right) \right\} \Psi,\tag{9}$$

and the inverse coordinate transformations

$$\begin{aligned} x\_N &= \frac{1}{\sqrt{N}} Q\_{N,x} + \sqrt{\frac{N-1}{N}} Q\_{N-1,x\prime} & y\_N &= \frac{1}{\sqrt{N}} Q\_{N,y} + \sqrt{\frac{N-1}{N}} Q\_{N-1,y\prime} \\ x\_{N-1} &= \frac{1}{\sqrt{N}} Q\_{N,x} - \frac{1}{\sqrt{N(N-1)}} Q\_{N-1,x} + \sqrt{\frac{N-2}{N-1}} Q\_{N-2,x\prime} \\ y\_{N-1} &= \frac{1}{\sqrt{N}} Q\_{N,y} - \frac{1}{\sqrt{N(N-1)}} Q\_{N-1,y} + \sqrt{\frac{N-2}{N-1}} Q\_{N-2,y} \end{aligned} \tag{10}$$

to evaluate the various integral terms contributing to (8).

The angular-momentum variance per particle of the ground state (3) depends on the dressed frequencies, Ω*<sup>x</sup>* and Ω*y*, and the number of particles *N*. Namely, unlike the respective position and momentum variances it depends explicitly on the interaction strength and the number of particles. 1 *N* ∆ 2 *L*ˆ *Z* is, of course, non-zero only for anisotropic traps [for isotropic traps, from *ω<sup>y</sup>* = *ω<sup>x</sup>* we get Ω*<sup>y</sup>* = Ω*<sup>x</sup>* and expression (8) then vanishes]. For non-interacting bosons, Equation (8) boils down to 1 *N* ∆ 2 *L*ˆ *Z* = <sup>1</sup> 4 (*ωy*−*ωx*) 2 *ωyωx* = <sup>1</sup> 4 *ωy ωx* −1 2 *ωy ωx* , the value for a single particle in the anisotropic trap <sup>1</sup> 2*ω*<sup>2</sup> *x x* <sup>2</sup> + <sup>1</sup> 2*ω*<sup>2</sup> *yy* 2 , which only depends on the trap anisotropy. Opposite to the non-vanishing of the angular-momentum variance, we note that the expectation value of the angular-momentum operator, <sup>1</sup> *N* <sup>h</sup>Ψ|*L*<sup>ˆ</sup> *<sup>Z</sup>*|Ψi, vanishes for any anisotropy *<sup>ω</sup><sup>y</sup> ωy* , interaction strength *λ*0, and number of particles *N*. This is straightforward to see since Ψ is even under reflection of all coordinates *X* → −*X* and separately of *Y* → −*Y*, whereas *L*ˆ *<sup>Z</sup>* is odd under reflection.

The anisotropic harmonic-interaction model (1) can be solved analytically at the mean-field level of theory as well, like in [29], also see [41]. Starting from the ansatz where each and every boson resides in one and the same orbital, the mean-field solution is given by

$$\begin{split} &\Phi^{GP}(\mathbf{r}\_{1},\ldots,\mathbf{r}\_{N}) = \\ &= \left(\frac{\sqrt{\omega\_{\mathbf{x}}^{2}+2\Lambda}}{\pi}\right)^{\frac{N}{4}} \left(\frac{\sqrt{\omega\_{\mathbf{y}}^{2}+2\Lambda}}{\pi}\right)^{\frac{N}{4}} e^{-\frac{1}{2}\sqrt{\omega\_{\mathbf{x}}^{2}+2\Lambda}\sum\_{j=1}^{N}x\_{j}^{2}} \times e^{-\frac{1}{2}\sqrt{\omega\_{\mathbf{y}}^{2}+2\Lambda}\sum\_{j=1}^{N}y\_{j}^{2}} \\ &= \Phi^{GP}(\mathbf{Q}\_{1},\ldots,\mathbf{Q}\_{N}) \\ &= \left(\frac{\sqrt{\omega\_{\mathbf{x}}^{2}+2\Lambda}}{\pi}\right)^{\frac{N}{4}} \left(\frac{\sqrt{\omega\_{\mathbf{y}}^{2}+2\Lambda}}{\pi}\right)^{\frac{N}{4}} e^{-\frac{1}{2}\sqrt{\omega\_{\mathbf{x}}^{2}+2\Lambda}\sum\_{k=1}^{N}Q\_{k,\mathbf{x}}^{2}} \times e^{-\frac{1}{2}\sqrt{\omega\_{\mathbf{y}}^{2}+2\Lambda}}\Sigma\_{k,\mathbf{y}}^{N}\mathbf{Q}\_{k,\mathbf{y}}^{2} \end{split} \tag{11}$$

where Λ = (*N* − 1)*λ*<sup>0</sup> is the interaction parameter and Λ > − *ω*2 *x* 2 the condition for a trapped solution. Like the many-body solution (3), the mean-field solution (11) in two spatial dimensions factorizes to a product of respective one-dimensional mean-field solutions.

The many-particle position variance computed at the mean-field level is given by

$$\frac{1}{N}\Delta\_{\mathbb{X},GP}^2 = \frac{1}{2\sqrt{\omega\_{\mathbb{X}}^2 + 2\Lambda}}, \qquad \frac{1}{N}\Delta\_{\mathbb{Y},GP}^2 = \frac{1}{2\sqrt{\omega\_{\mathbb{Y}}^2 + 2\Lambda}}.\tag{12}$$

and seen to be dressed by the interaction. Similarly, the many-particle momentum variance computed at the mean-field level is dressed by the interaction and given by

$$\frac{1}{N}\Delta\_{\mathbb{P}\_X,GP}^2 = \frac{\sqrt{\omega\_x^2 + 2\Lambda}}{2}, \qquad \frac{1}{N}\Delta\_{\mathbb{P}\_Y,GP}^2 = \frac{\sqrt{\omega\_y^2 + 2\Lambda}}{2}. \tag{13}$$

Interestingly, because the mean-field solution (11) is made of Gaussian functions, it satisfies the minimal uncertainty product <sup>1</sup> *N* ∆ 2 *X*ˆ ,*GP* 1 *N* ∆ 2 *<sup>P</sup>*ˆ*X*,*GP* <sup>=</sup> <sup>1</sup> *N* ∆ 2 *Y*ˆ,*GP* 1 *N* ∆ 2 *P*ˆ *<sup>Y</sup>*,*GP* <sup>=</sup> <sup>1</sup> 4 as well.

The many-particle angular-momentum variance computed at the mean-field level is given by

$$\frac{1}{N} \Delta\_{L,GP}^2 = \frac{1}{4} \frac{\left(\sqrt{\omega\_y^2 + 2\Lambda} - \sqrt{\omega\_x^2 + 2\Lambda}\right)^2}{\sqrt{\omega\_y^2 + 2\Lambda}\sqrt{\omega\_x^2 + 2\Lambda}},\tag{14}$$

where we have made use of the structure and symmetries of Φ*GP* ,

$$\hat{L}\_Z \Phi^{GP} = -\frac{1}{i} \left( \sqrt{\omega\_y^2 + 2\Lambda} - \sqrt{\omega\_x^2 + 2\Lambda} \right) \sum\_{j=1}^N \mathbf{x}\_j \mathbf{y}\_j \Phi^{GP},\tag{15}$$

to arrive at the final expression.

The relation between the mean-field and many-body variances deserves a discussion. Their difference is used to define position, momentum, and angular-momentum correlations in the system. For the position and momentum variances, the following ratios hold,

$$\frac{\frac{1}{N}\Delta\_{\hat{X},GP}^{2}}{\frac{1}{N}\Delta\_{\hat{X}}^{2}}=\frac{1}{\sqrt{1+\frac{2\Lambda}{\omega\_{x}^{2}}}},\qquad\frac{\frac{1}{N}\Delta\_{\hat{Y},GP}^{2}}{\frac{1}{N}\Delta\_{\hat{Y}}^{2}}=\frac{1}{\sqrt{1+\frac{2\Lambda}{\omega\_{y}^{2}}}},$$

$$\frac{\frac{1}{N}\Delta\_{\hat{P}\_{X,GP}}^{2}}{\frac{1}{N}\Delta\_{\hat{P}\_{X}}^{2}}=\sqrt{1+\frac{2\Lambda}{\omega\_{x}^{2}}},\qquad\frac{\frac{1}{N}\Delta\_{\hat{P}\_{Y,GP}}^{2}}{\frac{1}{N}\Delta\_{\hat{P}\_{Y}}^{2}}=\sqrt{1+\frac{2\Lambda}{\omega\_{y}^{2}}},\tag{16}$$

obviously for any number of particles *N*. These ratios simply imply that, since repulsion (Λ < 0) broadens the position density, the many-body position variance is smaller than the corresponding mean-field one for repulsive interaction, and vise verse for attraction (Λ > 0). Inversely, since repulsion narrows the momentum density, the many-body momentum variance is larger than the corresponding mean-field one for repulsive interaction, and vise versa for attraction. Furthermore, both the position and momentum variances per particle exhibit the same anisotropies as the respective densities for any interaction parameter Λ, namely, if <sup>1</sup> *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> ,*GP* <sup>&</sup>gt; <sup>1</sup> *N* ∆ 2 *<sup>Y</sup>*ˆ,*GP* then <sup>1</sup> *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> <sup>&</sup>gt; <sup>1</sup> *N* ∆ 2 *Y*ˆ is satisfied and, analogously, if <sup>1</sup> *N* ∆ 2 *<sup>P</sup>*ˆ*X*,*GP* <sup>&</sup>lt; <sup>1</sup> *N* ∆ 2 *P*ˆ *<sup>Y</sup>*,*GP* then <sup>1</sup> *N* ∆ 2 *P*ˆ*X* < <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* is satisfied. We shall return to these relations and the anisotropy of the variance in the numerical example below.

We now extend the above discussion to the particle limit, in which the energy per particle, densities per particle, and reduced densities per particle at the mean-field and many-body levels of theory coincide, see in the context of the harmonic-interaction model [16]. Particularly, the system of

bosons becomes 100% condensed. The results (16) for the position and momentum variances hold at the particle limit as well, owing to the center-of-mass separability for any number of particles, namely, lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *X*ˆ ,*GP* lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *X*ˆ = <sup>r</sup> 1 1+ <sup>2</sup><sup>Λ</sup> *ω*2 *x* , lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *Y*ˆ,*GP* lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *Y*ˆ = <sup>r</sup> 1 1+ <sup>2</sup><sup>Λ</sup> *ω*2 *y* , lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *P*ˆ *X*,*GP* lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *P*ˆ *X* = q 1 + <sup>2</sup><sup>Λ</sup> *ω*2 *x* , and lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *P*ˆ *Y* ,*GP* lim*N*→<sup>∞</sup> 1 *N* ∆ 2 *P*ˆ *Y* = r

1 + <sup>2</sup><sup>Λ</sup> *ω*2 *y* . For the angular-momentum variance the limit has to be taken explicitly for each of the terms in (8). First are the frequencies (4), for which we have at the limit of an infinite number of bosons when Λ is held fixed

$$\lim\_{N \to \infty} \Omega\_{\mathbf{x}} = \sqrt{\omega\_{\mathbf{x}}^2 + 2\Lambda}, \qquad \lim\_{N \to \infty} \Omega\_{\mathbf{y}} = \sqrt{\omega\_{\mathbf{y}}^2 + 2\Lambda}. \tag{17}$$

Then, the angular-momentum variance takes on the appealing form

$$\lim\_{N \to \infty} \frac{1}{N} \Delta\_{\mathbb{Z}}^2 = \frac{1}{4} \frac{\left(\sqrt{\omega\_y^2 + 2\Lambda} - \sqrt{\omega\_x^2 + 2\Lambda}\right)^2}{\sqrt{\omega\_y^2 + 2\Lambda}\sqrt{\omega\_x^2 + 2\Lambda}} \left[1 + \left(\sqrt{1 + \frac{2\Lambda}{\omega\_y^2}} - 1\right)\left(\sqrt{1 + \frac{2\Lambda}{\omega\_x^2}} - 1\right)\right]. \tag{18}$$

Comparing (18) to the mean-field expression (14), it is instrumental to prescribe their ratio at the limit of an infinite number of particles (where, as mentioned above, the density per particle and other properties coincide),

$$\frac{\lim\_{N\to\infty}\frac{1}{N}\Delta\_{L\_{Z,GP}}^2}{\lim\_{N\to\infty}\frac{1}{N}\Delta\_{L\_Z}^2} = \frac{1}{1 + \left(\sqrt{1 + \frac{2\Lambda}{\omega\_y^2}} - 1\right)\left(\sqrt{1 + \frac{2\Lambda}{\omega\_x^2}} - 1\right)}\tag{19}$$

which is always smaller than 1 for interacting bosons in the anisotropic trap. Furthermore, we see that for attractive interaction the many-body variance can become much larger than the mean-field quantity in the anisotropic trap, signifying the growing necessity of the many-body treatment, even when the system is 100% condensed. This concludes our investigation of a solvable anisotropic many-boson model in which the variances of the momentum, position, and angular-momentum many-particle operators can be computed and investigated analytically and their values at the many-body and mean-field levels of theory compared and contrasted.

#### **3. Bosons in an Annulus Subject to a Tilt**

In most scenarios of interest, the position, momentum, and angular-momentum variance cannot be computed analytically. This in many cases is the situation when symmetries are lifted. Moreover, even when the variances can be computed for the ground state, like in the previous Section 2, their values for an out-of-equilibrium scenario are rarely within analytical reach. This would be the situation of the present investigation.

Bosons in rings, annuli, and shells have attracted considerable attention [42–66]. Here we consider weakly interacting bosons initially prepared in the ground state of a two-dimensional annulus. The annulus is then suddenly slightly tilted, leading to an out-of-equilibrium dynamics in an anisotropic setup. We build on and extend the study of bosons' dynamics in an annulus within an isotropic setup [19] (for which, e.g., the angular-momentum variance is 0). We analyze the BEC dynamics in terms of its time-dependent variances and other quantities of relevance, see Figures 1–7 below.

We consider the out-of-equilibrium dynamics governed by the time-dependent many-particle Schrödinger equation in two spatial dimensions, *H*ˆ (**r**1, . . . ,**r***N*)Ψ(**r**1, . . . ,**r***N*; *t*) = *i ∂*Ψ(**r**<sup>1</sup> ,...,**r***N*;*t*) *∂t* . The bosons are initially prepared in the ground state of the annulus, see Figure 1 in [19]. The trap potential is given by *V*ˆ(**r**) = 0.05**r** <sup>4</sup> + *V*0*e* − **r** 2 <sup>2</sup> , with a barrier of heights *V*<sup>0</sup> = 5 and 10 throughout this work. The interaction between the bosons is repulsive and taken to be *λ*0*W*(**r** − **r** 0 ) = *λ*0*e* − (**r**−**r** 0 ) 2 <sup>2</sup> ,

where the interaction strengths are *λ*<sup>0</sup> = 0.02 and 0.04 throughout this work. The form and extant of the interaction potential do not have a qualitative influence on the physics to be described below. At time *t* = 0 a linear term is added such that *V*(**r**) = 0.05**r** <sup>4</sup> + *V*0*e* − **r** 2 <sup>2</sup> + 0.01*x*. The physical meaning of the added potential is that a constant force pointing to the left is suddenly acting on the interacting bosons. Geometrically, the annulus can be considered to be slightly tilted to the left. Symmetry-wise, the isotropy of the potential is lifted and anisotropy sets in. All in all, the interacting bosons are not in their ground state any more and out-of-equilibrium dynamics emerges.

To compute the time-dependent many-boson wavefunction we use the multiconfigurational time-dependent Hartree for bosons (MCTDHB) method [67–69]. MCTDHB represents the wavefunction as a variationally optimal ansatz which is a linear-combination of all time-dependent permanents generated by distributing the *N* bosons over *M* time-adaptive orbitals. The quality of the wavefunction increases with *M* and convergence of quantities of interest is attained. The theory, applications, benchmarks, and extensions of MCTDHB are extensively discussed in the literature, see, e.g., Refs. [70–95]. Here we employ the numerical implementation in [96,97] both for preparing the ground state [98] (using imaginary-time propagation) and real-time dynamics. Finally, we mention that MCTDHB is the bosonic version of the nearly three-decades-established distinguishable-particle multiconfigurational time-dependent Hartree method frequently used (alongside its extensions) in molecular physics [99–105].

From the time-dependent wavefunction Ψ(**r**1, . . . ,**r***N*; *t*), here normalized to 1, we compute properties of interest. The reduced one-particle density matrix is defined as *ρ*(**r**,**r** 0 ; *t*) = *N* R *d***r**<sup>2</sup> · · · *d***r***N*Ψ<sup>∗</sup> (**r** 0 ,**r**2, . . . ,**r***N*; *t*)Ψ(**r**,**r**2, . . . ,**r***N*; *t*) = ∑*<sup>j</sup> nj*(*t*)*φ* ∗ *j* (**r** 0 ; *t*)*φj*(**r**; *t*), where {*φj*(**r**; *t*)} are the natural orbitals and {*nj*(*t*)} the natural occupations. The number of particles residing outside the condensed mode *φ*1(**r**; *t*), i.e., the total number of depleted particles, is given by ∑*j*><sup>1</sup> *nj*(*t*) = *N* − *n*1(*t*). Analogously, the reduced two-particle density matrix is given by *ρ*(**r**1,**r**2,**r** 0 1 ,**r** 0 2 ; *t*) = *N*(*N* − 1) R *d***r**<sup>3</sup> · · · *d***r***N*Ψ<sup>∗</sup> (**r** 0 1 ,**r** 0 2 ,**r**3, . . . ,**r***N*; *t*)× Ψ(**r**1,**r**2,**r**3, . . . ,**r***N*; *t*) = ∑*jpkq ρjpkq*(*t*)*φ* ∗ *j* (**r** 0 1 ; *t*)*φ* ∗ *p* (**r** 0 2 ; *t*)*φ<sup>k</sup>* (**r**1; *t*)*φq*(**r**2; *t*), from which the variance of a many-particle operator *<sup>A</sup>*<sup>ˆ</sup> = <sup>∑</sup>*<sup>j</sup> a*ˆ(**r**) is computed,

$$\begin{split} \frac{1}{N} \Delta\_{\hat{A}}^{2}(t) &= \frac{1}{N} \left( \langle \Psi(t) | \hat{A}^{2} | \Psi(t) \rangle - \langle \Psi(t) | \hat{A} | \Psi(t) \rangle^{2} \right) \\ &= \frac{1}{N} \left\{ \sum\_{j} n\_{j}(t) \int d\mathbf{r} \phi\_{j}^{\*}(\mathbf{r};t) \hat{\mathfrak{a}}^{2}(\mathbf{r}) \phi\_{j}(\mathbf{r};t) - \left[ \sum\_{j} n\_{j}(t) \int d\mathbf{r} \phi\_{j}^{\*}(\mathbf{r};t) \hat{\mathfrak{a}}(\mathbf{r}) \phi\_{j}(\mathbf{r};t) \right]^{2} \right. \\ &\left. + \sum\_{j \neq k \neq} \rho\_{j \text{p} \text{k} \text{q}}(t) \left[ \int d\mathbf{r} \phi\_{j}^{\*}(\mathbf{r};t) \hat{\mathfrak{a}}(\mathbf{r}) \phi\_{k}(\mathbf{r};t) \right] \left[ \int d\mathbf{r} \phi\_{p}^{\*}(\mathbf{r};t) \hat{\mathfrak{a}}(\mathbf{r}) \phi\_{q}(\mathbf{r};t) \right] \right\}. \end{split} \tag{20}$$

To compute the various terms for the position, momentum, and angular-momentum variance numerically we work in coordinate representation and operate on orbitals first with coordinate derivatives and then with coordinate multiplications. Thus, for the position operator *a*ˆ(**r**) = *x*ˆ and *a*ˆ 2 (**r**) = *x*ˆ <sup>2</sup> and likewise for *a*ˆ(**r**) = *y*ˆ, for the momentum operator *a*ˆ(**r**) = <sup>1</sup> *i ∂ ∂x* and *a*ˆ 2 (**r**) = <sup>−</sup> *<sup>∂</sup>* 2 *∂x* 2 and likewise for *a*ˆ(**r**) = <sup>1</sup> *i ∂ ∂y* , and for the angular-momentum operator *a*ˆ(**r**) = <sup>1</sup> *i x ∂ <sup>∂</sup><sup>y</sup>* − *y ∂ ∂x* and *a*ˆ 2 (**r**) = −*x* 2 *∂* 2 *∂y* <sup>2</sup> − *y* 2 *∂* 2 *∂x* <sup>2</sup> + 2*yx <sup>∂</sup> ∂y ∂ <sup>∂</sup><sup>x</sup>* + *x ∂ <sup>∂</sup><sup>x</sup>* + *y ∂ ∂y* . For the numerical solution we use a grid of 64<sup>2</sup> points in a box of size [−8, 8) × [−8, 8) with periodic boundary conditions. Convergence of the results with respect to the number of grid points has been checked using a grid of 128<sup>2</sup> points.

We begin with the dynamics of *N* = 10 bosons in the annulus. Following the sudden tilt of the potential, the bosons start to flow to the left. To quantify their sloshing dynamics, Figure 1 shows the time-dependent center-of-mass, <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*), for the two barrier heights, *<sup>V</sup>*<sup>0</sup> <sup>=</sup> <sup>5</sup> and *<sup>V</sup>*<sup>0</sup> <sup>=</sup> 10, and the two interaction strengths, *λ*<sup>0</sup> = 0.02 and *λ*<sup>0</sup> = 0.04 [we mention that <sup>1</sup> *N* <sup>h</sup>Ψ|*Y*ˆ|Ψi(*t*) = <sup>0</sup> due to the *<sup>Y</sup>* → −*<sup>Y</sup>* reflection symmetry]. The dynamics of <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*) appears to be almost periodic

and rather simple. We examine the amplitude and frequency of oscillations. It is useful to compare the amplitude of the center-of-mass motion with the radius of the (un-tilted) annulus. The radius of the density at its maximal value, *R*, is determined numerically using a computation with a resolution of 256<sup>2</sup> grid points as *R* = 1.75(0) for *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02, and *R* = 2.06(2) for *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02 [19]. From Figure 1 we see that the amplitude is about 13–25% of the radius, implying a mild sloshing of the density along the tilted annulus. The amplitude increases with the radius of the annulus and decreases with the interaction strength, where the latter implies that it is more difficult to compress the BEC for a stronger interaction. The decrease of the frequency of oscillations with *R* (*V*0) and increase with *λ*<sup>0</sup> are compatible with angular excitations, also see [19]. Last but not least, convergence with *M* is clearly seen. In fact, here already *M* = 1 orbitals accurately describe the center-of-mass dynamics for short and intermediate times, and *M* = 3 orbitals for all times.

**Figure 1.** Center-of-mass dynamics following a potential quench. The mean-field (*M* = 1 time-adaptive orbitals) and many-body (using *M* = 3, 5, 7, 10, 12, 14, and 15, 16 time-adaptive orbitals) time-dependent expectation value of the center-of-mass, <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*), of *<sup>N</sup>* <sup>=</sup> <sup>10</sup> bosons in the annuli with barrier heights and interaction strengths: (**a**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02; (**b**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.04; (**c**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02; and (**d**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04 following a sudden potential tilt by 0.01*x*. The corresponding depletions are plotted in Figure 2 and the respective position, momentum, and angular-momentum variances in Figures 3–5. See the text for more details. The quantities shown are dimensionless.

Figure 2 depicts the total number of depleted particles, *N* − *n*1(*t*), out of *N* = 10 bosons in the tilted annulus. During the dynamics, the depletion is rather small, ranging from less than 0.012 of a particle out of *N* = 10 particles (0.12%) for *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02 to less than 0.065 of a particle out of *N* = 10 particles (0.65%) for *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04. Generally, the depletion increases with the annulus radius and interacting strength, implying angular excitations, see [19]. Finally, convergence with *M* is clearly seen. Now, *M* = 3 orbitals nicely follow and *M* = 5 orbitals accurately describe the depletion dynamics, see Figure 2. The small amount of time-dependent depletion is in line with the accurate

description of the center-of-mass dynamics by *M* = 1 time adaptive orbitals, see Figure 1. Let us continue to the variances.

**Figure 2.** Depletion dynamics following a potential quench. The time-dependent total number of depleted particles, *N* − *n*1(*t*), of *N* = 10 bosons following a sudden potential tilt by 0.01*x* for annuli with barrier heights and interaction strengths (**a**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02; (**b**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.04; (**c**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02; and (**d**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04. *M* = 3, 5, 7, 10, 12, 14, and 15, 16 time-adaptive orbitals are used. The respective position, momentum, and angular-momentum variances are plotted in Figures 3–5. See the text for more details. The quantities shown are dimensionless.

Figure 3 plots the time-dependent many-particle position variance per particle, <sup>1</sup> *N* ∆ 2 *X*ˆ (*t*) and 1 *N* ∆ 2 *Y*ˆ (*t*), for the two barrier heights and two interaction strengths. There are several features that immediately are seen. First, since rotational symmetry is lifted, the dynamics of respective quantities along the x-axis and y-axis are different [note that at *t* = 0 the variances <sup>1</sup> *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> <sup>=</sup> <sup>1</sup> *N* ∆ 2 *Y*ˆ because the initial condition is the ground state of the un-tilted, isotropic annulus]. The mean-field (*M* = 1) and many-body (*M* ≥ 3) values are clearly separated from each other, and the former lie about 10–25% above the latter depending on the repulsion strength and barrier height, also see [11,19]. This is despite the small amount of depletion, see Figure 2. Furthermore, the many-body and mean-field variances do not cross each other, see Figure 3, indicating that the dynamics is mild and sufficiently close to the ground state and low-lying manifold of excited states (compare to [18] with interaction-quench dynamics in a single trap).

The mean-field position variance accounts for the geometry of the annulus and shape of the density and weakly depends on the interaction strength. The many-body position variance incorporates the (small amount of) depletion and hence strongly depends on the interaction strength. Both the mean-field and many-body variances oscillate with a relatively small amplitude, albeit with a different frequencies' content, see in this respect [20]. This amplitude slightly decreases with the repulsion strength, which correlates with the dependence of the center-of-mass dynamics on the interaction strength, see Figure 1. Moreover, the amplitude of oscillations of the y-axis variances is smaller than that of the x-axis variances, since the sloshing dynamics is primarily along the *x* direction. Last but not least is the so-called opposite anisotropy of the (position) variance [18]. During the dynamics, there can occur instances where <sup>1</sup> *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> <sup>&</sup>gt; <sup>1</sup> *N* ∆ 2 *Y*ˆ at the many-body level (*<sup>M</sup>* <sup>≥</sup> 3) whereas <sup>1</sup> *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> <sup>&</sup>lt; <sup>1</sup> *N* ∆ 2 *Y*ˆ at the mean-field level (*M* = 1) [or, in principle, vice versa, i.e., <sup>1</sup> *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> <sup>&</sup>lt; <sup>1</sup> *N* ∆ 2 *Y*ˆ at the many-body level whereas 1 *N* ∆ 2 *<sup>X</sup>*<sup>ˆ</sup> <sup>&</sup>gt; <sup>1</sup> *N* ∆ 2 *Y*ˆ at the mean-field level]. Examples for the former can be readily found for *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02, see Figure 3c,g around *t* = 70, and for *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04, see Figure 3d,h around *t* = 100, signifying among others that correlations 'win' over shape. Finally, we see that already *M* = 3 orbitals accurately describe the dynamics of the position variance.

We move to the momentum variance and also make contact with the results of the position variance. Figure 4 displays the many-particle momentum variance per particle, <sup>1</sup> *N* ∆ 2 *P*ˆ*X* (*t*) and <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* (*t*), for *V*<sup>0</sup> = 5, *V*<sup>0</sup> = 10 and *λ*<sup>0</sup> = 0.02, *λ*<sup>0</sup> = 0.04. Just like the results of the position variance, since rotational symmetry is lifted the dynamics of respective quantities along the x-axis and y-axis are different [the initial conditions imply <sup>1</sup> *N* ∆ 2 *P*ˆ*X* = <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* at *t* = 0]. The mean-field (*M* = 1) and many-body (*M* ≥ 3) values are, again, separated from each other, but now the former lie below the later, and there is only about 1–4% of a difference depending on the repulsion strength and barrier height, also see [12,19]. Thus, the momentum variance rather weakly depends on the (small amount of) depletion. This is because the matrix elements in (20) are typically smaller with the momentum operator than with the position operator. Yet, despite their small difference, the many-body and mean-field momentum variances do not cross each other, see Figure 4 (contrast with the interaction-quench dynamics in a single trap in [18]).

It is instructive to analyze the momentum-variance dynamics at short times. Whereas ∆ 2 *P*ˆ*X* (*t*) primarily increases, ∆ 2 *P*ˆ *Y* (*t*) mainly decreases. This matches the geometry of the sloshing dynamics in the tilted annulus, in which bosons from the 'north' and 'south' poles (on the y-axis) start to move to the left and accumulate in the 'west' pole (on the x-axis), and that the cross section of the rim of an annulus is enlarged when moving away from the center of the annulus. In other words, the dynamics of the momentum variances at short times when moving to the left reflects the relative localization of the bosons in the *x* direction and the effective broadening of the wavepacket along the *y* direction. Both the mean-field and many-body variances oscillate with a very small amplitude, note the scale on the y-axis in Figure 4. The high-frequency oscillations mark high-energy radial excitations across the (tight) annulus rim [19]. Like for the position variance, the amplitude of oscillations of the y-axis momentum variances is smaller than that of the x-axis momentum variances. Finally, we see that already *M* = 3 orbitals accurately describe the dynamics of the momentum variance; the difference to the *M* > 3 results is lower than 1%.

We now move to the angular-momentum variance and an interesting inter-connection with the momentum variance. Figure 5 presents the many-particle angular-momentum variance per particle, 1 *N* ∆ 2 *L*ˆ *Z* (*t*), for the two barrier heights, *V*<sup>0</sup> = 5 and *V*<sup>0</sup> = 10, and the two interaction strengths, *λ*<sup>0</sup> = 0.02 and *λ*<sup>0</sup> = 0.04. There are several features seen in the dynamics. Since rotational symmetry is lifted, 1 *N* ∆ 2 *L*ˆ *Z* 6= 0 expect for the initial conditions at *t* = 0 (the values of the minima for *t* > 0, see below, are close to but not 0). The dynamics of <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*) appears to be almost periodic and rather regular, more than that for the respective position and momentum variances, compare to Figures 3 and 4. On the other end, focusing on the dynamics of the center-of-mass in Figure 1, one can clearly observe correlation between the two quantities; Whenever <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*) has a minimum, i.e., the bosons are maximally localized to the left, <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*) has a maximum, and whenever <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*) has a maximum (which value is about 0), i.e., the bosons are momentarily, approximately equally distributed along the annulus, <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*) has a minimum (which value, as mentioned above, is close to 0). Furthermore, the frequencies of the two quantities as well as their relative amplitudes as a function of the barrier height and interaction strength are alike. These observations call for a dedicated analysis.

**Figure 3.** Position variance dynamics following a potential quench. The mean-field (*M* = 1 timeadaptive orbitals) and many-body (using *M* = 3, 5, 7, 10, 12, 14, and 15, 16 time-adaptive orbitals) time-dependent position variances per particle, <sup>1</sup> *N* ∆ 2 *X*ˆ (*t*) [left column, panels (**a**–**d**)] and <sup>1</sup> *N* ∆ 2 *Y*ˆ (*t*) [right column, panels (**e**–**h**)], of *N* = 10 bosons in the annuli with barrier heights and interaction strengths (**a**,**e**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02; (**b**,**f**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.04; (**c**,**g**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02; and (**d**,**h**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04 following a sudden potential tilt by 0.01*x*. The respective depletions are plotted in Figure 2. See the text for more details. The quantities shown are dimensionless.

**Figure 4.** Momentum variance dynamics following a potential quench. The mean-field (*M* = 1 time-adaptive orbitals) and many-body (using *M* = 3, 5, 7, 10, 12, 14, and 15, 16 time-adaptive orbitals) time-dependent momentum variances per particle, <sup>1</sup> *N* ∆ 2 *P*ˆ*X* (*t*) [left column, panels (**a**–**d**)] and <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* (*t*) [right column, panels (**e**–**h**)], of *N* = 10 bosons in the annuli with barrier heights and interaction strengths (**a**,**e**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02; (**b**,**f**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.04; (**c**,**g**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02; and (**d**,**h**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04 following a sudden potential tilt by 0.01*x*. The respective depletions are plotted in Figure 2. See the text for more details. The quantities shown are dimensionless.

**Figure 5.** Angular-momentum variance dynamics following a potential quench. The mean-field (*M* = 1 time-adaptive orbitals) and many-body (using *M* = 3, 5, 7, 10, 12, 14, and 15, 16 time-adaptive orbitals) time-dependent angular-momentum variance per particle, <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*), of *N* = 10 bosons in the annuli with barrier heights and interaction strengths (**a**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.02; (**b**) *V*<sup>0</sup> = 5, *λ*<sup>0</sup> = 0.04; (**c**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.02; and (**d**) *V*<sup>0</sup> = 10, *λ*<sup>0</sup> = 0.04 following a sudden potential tilt by 0.01*x*. The respective depletions are plotted in Figure 2. See the text for more details. The quantities shown are dimensionless.

To shed light on the above dynamics of the angular-momentum variance, see Figure 5, we analyze the translational properties of variances in Appendix A. Whereas the position variances and, trivially, the momentum variances, are translationally invariant, this invariance does not hold for the angular-momentum variance. If a wavepacket prepared in the origin has angular-momentum variance <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* , then several terms are added when the wavepacket is translated to the point (*a*, *b*) in plane, and angular-momentum variance is thereafter computed, see Equation (A3). Now, if this wavepacket is rotationally symmetric, i.e., <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* = 0, then several of the terms in (A3) vanish due to spatial symmetry and we are left with the appealing relation, <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* Ψ(*a*,*b*) = *a* 2 1 *N* ∆ 2 *P*ˆ *Y* Ψ + *b* 2 1 *N* ∆ 2 *P*ˆ*X* Ψ [Equation (A4)], connecting the angular-momentum variance of Ψ(*a*, *b*) localized at (*a*, *b*) and of Ψ at the origin. The meaning of this relation is that the momentum variances, <sup>1</sup> *N* ∆ 2 *P*ˆ*X* and <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* , together with the spatial translations along the y-axis and x-axis, respectively, determine the angular-momentum variance of a translated wavepacket (rotationally symmetric at the origin).

Returning to and combining Figure 5 for the angular-momentum variance, Figure 1 for the center-of-mass dynamics, and Figure 4e–h for <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* , we can now discuss and explain their inter-connection. Explicitly, the center-of-mass dynamics is analogous to translating the wavepacket along the x-axis (back and forth to the left), hence, according to Equation (A4), <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* is needed. This is why the dependencies of the frequency and amplitude of oscillations of <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*) on the barrier height and interaction strength nicely follow, respectively, those of <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*), compare

Figures 1 and 5. What is the role of <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* then? The momentum variance helps us understand the deviations between the many-body and mean-field results in Figure 5. We see that the maxima of the many-body <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*) (*M* > 3) are larger than the maxima of the mean-field <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*) (*M* = 1). The difference is about 7–25% (compare to the low depletion, Figure 2), depending on *V*<sup>0</sup> and *λ*0, and follows the respective trend of the many-body and mean-field results for <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* , see Figure 4e–h. We note that, although the wavepacket describing the bosons dynamics in the tilted annulus is not a translated, rotationally invariant wavepacket, and the values of deviations (in percents) between the many-body and mean-field results are actually larger for ∆ 2 *L*ˆ *Z* (at the maxima) than for ∆ 2 *P*ˆ *Y* , we find the above analytically based analysis to well explain the numerical findings and trends. Last but not least, a close inspection of the many-body and mean-field curves of the angular-momentum variance in Figure 5 shows that there are instances when they cross each other, i.e., one is smaller or larger than the other. This is in contrast with the non-crossing of the many-body and mean-field position and momentum variances, see Figures 3 and 4, respectively. Finally, we find that already *M* = 3 time-adaptive orbitals accurately describe the dynamics of the angular-momentum variance.

Our investigations are nearing their end, what is left to explore is the behavior of the position, momentum, and angular-momentum variances at the particle limit. Which of the above-described detailed findings, plotted in Figures 1–5 for a rather small (*N* = 10 bosons) yet weakly depleted BEC, survive this limit? To answer the question, we concentrate on the system with the higher barrier, *V*<sup>0</sup> = 10, and stronger interaction (for *N* = 10 bosons), *λ*<sup>0</sup> = 0.04. We hence fix the interaction parameter Λ = *λ*0(*N* − 1) = 0.36, and compute and compare the dynamics for *N* = 10, *N* = 100, and *N* = 1000 bosons using *M* = 3 time-adaptive orbitals. We have seen for *N* = 10 bosons that *M* = 3 time-adaptive orbitals accurately describe the variances. This implies that, keeping the interaction parameter Λ fixed while increasing the number of particles *N*, using *M* = 3 time-adaptive orbitals for calculating the variances will be (at least) as accurate as for *N* = 10 particles, see in this respect [24]. Before we proceed, a methodological remark. Examining the convergence of properties with the number of particles for *N* = 10, *N* = 100, and *N* = 1000 bosons is (still) far away from infinity, see in this respect [15]. We hence use, interchangeably, the term *en route* to the particle limit. We shall see below that, in effect, the particle limit is practically well achieved for the variances already for *N* = 1000 bosons.

Figure 6 prints the total number of depleted particles, *N* − *n*1(*t*), for *N* = 10, *N* = 100, and *N* = 1000 bosons for Λ = 0.36 and *V*<sup>0</sup> = 10 using *M* = 3 time-adaptive orbitals. Convergence of the number of depleted particles with *N* is nicely seen. Since *N* − *n*1(*t*) converges to a finite (and small) value with *N*, the bosons are becoming 100% condensed in the limit of an infinite number of particles, i.e., *n*1 (*t*) *<sup>N</sup>* → 1 as *N* → ∞, at least up to the maximal time of the computation, *t* = 100.

Figure 7 exhibits the position variances per particle, <sup>1</sup> *N* ∆ 2 *X*ˆ (*t*) and <sup>1</sup> *N* ∆ 2 *Y*ˆ (*t*), momentum variances per particle, <sup>1</sup> *N* ∆ 2 *P*ˆ*X* (*t*) and <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* (*t*), angular-momentum variance per particle, <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* (*t*), and the expectation value of the center-of-mass, <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*), for *<sup>N</sup>* <sup>=</sup> 10, *<sup>N</sup>* <sup>=</sup> 100, and *<sup>N</sup>* <sup>=</sup> <sup>1000</sup> bosons and for Λ = 0.36 and *V*<sup>0</sup> = 10 using *M* = 3 time-adaptive orbitals. Once again, convergence of each of the quantities with *N* is clearly seen. Yet, whereas the center-of-mass dynamics converges to the mean-field dynamics when the number of particles is increased, the variances exhibit many-body dynamics which converges nicely with *N*, but not to the respective mean-field dynamics. Beyond that, all the above results, for the frequencies, amplitudes, anisotropies, inter-connections, and particularly the differences between the many-body and mean-field position, momentum, and angular-momentum variances persist at the limit of infinite number of particles, despite the bosons becoming 100% condensed. This brings the present analysis to an end.

**Figure 6.** Depletion dynamics following a potential quench *en route* to the particle limit. The time-dependent total number of depleted particles, *N* − *n*1(*t*), of *N* = 10, *N* = 100, and *N* = 1000 bosons with interaction parameter Λ = *λ*0(*N* − 1) = 0.36 for an annulus with barrier height *V*<sup>0</sup> = 10 following a sudden potential tilt by 0.01*x*. The number of time-adaptive orbitals is *M* = 3. The respective position, momentum, and angular-momentum variances along with the expectation value of the center-of-mass are plotted in Figure 7. See the text for more details. The quantities shown are dimensionless.

**Figure 7.** *Cont*.

**Figure 7.** Position, momentum, and angular-momentum variance dynamics following a potential quench *en route* to the particle limit. The mean-field (*M* = 1 time-adaptive orbitals) and many-body (using *M* = 3 time-adaptive orbitals) time-dependent position variances per particle, (**a**) 1 *N* ∆ 2 *X*ˆ (*t*) and (**b**) 1 *N* ∆ 2 *Y*ˆ (*t*), momentum variances per particle, (**c**) 1 *N* ∆ 2 *P*ˆ*X* (*t*) and (**d**) 1 *N* ∆ 2 *P*ˆ *Y* (*t*), and angular-momentum variance per particle, (**f**) 1 *N* ∆ 2 *L*ˆ *Z* (*t*), of *N* = 10, *N* = 100, and *N* = 1000 bosons with interaction parameter Λ = *λ*0(*N* − 1) = 0.36 for an annulus with barrier height *V*<sup>0</sup> = 10 following a sudden potential tilt by 0.01*x*. (**e**) The time-dependent expectation value of the center-of-mass, <sup>1</sup> *N* <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψi(*t*). The respective depletions are plotted in Figure 6. See the text for more details. The quantities shown are dimensionless.

#### **4. Summary and Outlook**

In the present work we studied, analytically and numerically, the position, momentum, and especially the angular-momentum variance of interacting bosons trapped in a two-dimensional anisotropic trap for static and dynamic scenarios. Explicitly, we investigated the ground state of the anisotropic harmonic-interaction model in two spatial dimensions analytically and researched the out-of-equilibrium dynamics of repulsive bosons in tilted two-dimensional annuli numerically accurately by using the MCTDHB method. The differences between the variances at the mean-field level, which are attributed to the shape of the density per particle, and the respective variances at the many-body level, which incorporate a small amount of depletion outside the condensed mode, were used to characterize sometimes large position, momentum, and angular-momentum correlations in the BEC for finite systems and at the limit of an infinite number of particles where the bosons are 100% condensed. Finally, we also explored and utilized inter-connections between the variances, particularly between the angular-momentum and momentum variances, through the analysis of their translational properties.

There are many intriguing directions to follow out of which we list three below. First, variances of BECs in the rotating frame of reference in which high-lying excitations become low-energy excitations and even the ground state. Second, angular-momentum variance of a BEC flowing past an obstacle in which the mean angular-momentum variance vanishes. And third, variances in three-dimensional geometries lacking lower-dimensional analogs, such as a Möbius strip. In all these cases, whether considering a few interacting bosons or a BEC in the particle limit, interesting and exciting results are expected.

**Funding:** This research was funded by Israel Science Foundation (Grant No. 600/15).

**Acknowledgments:** This research was supported by the Israel Science Foundation (Grant No. 600/15). We thank Kaspar Sakmann, Anal Bhowmik, Sudip Haldar, and Raphael Beinke for discussions. Computation time on the BwForCluster, the High Performance Computing system Hive of the Faculty of Natural Sciences at University of Haifa, and the Cray XC40 system Hazelhen at the High Performance Computing Center Stuttgart (HLRS) is gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Variances and Translations**

Consider the many-particle translation operator in two spatial dimensions *e* <sup>−</sup>*i*(*P*ˆ*<sup>X</sup> <sup>a</sup>*+*P*<sup>ˆ</sup> *<sup>Y</sup>b*) , where *P*ˆ*<sup>X</sup>* = ∑ *N j*=1 *p*ˆ*x*,*<sup>j</sup>* and *P*ˆ *<sup>Y</sup>* = ∑ *N j*=1 *p*ˆ*y*,*<sup>j</sup>* . Its operation on a multi-particle wavefunction Ψ is given by *e* <sup>−</sup>*i*(*P*ˆ*<sup>X</sup> <sup>a</sup>*+*P*<sup>ˆ</sup> *<sup>Y</sup>b*)Ψ(*x*1, *<sup>y</sup>*1, . . . , *<sup>x</sup>N*, *<sup>y</sup>N*) = <sup>Ψ</sup>(*x*<sup>1</sup> <sup>−</sup> *<sup>a</sup>*, *<sup>y</sup>*<sup>1</sup> <sup>−</sup> *<sup>b</sup>*, . . . , *<sup>x</sup><sup>N</sup>* <sup>−</sup> *<sup>a</sup>*, *<sup>y</sup><sup>N</sup>* <sup>−</sup> *<sup>b</sup>*) <sup>≡</sup> <sup>Ψ</sup>(*a*, *<sup>b</sup>*). What are the implications on the variances when computed with respect to the translated wavefunction Ψ(*a*, *b*)?

For the position operator *X*ˆ = ∑ *N j*=1 *x*ˆ*j* (and equivalently for *Y*ˆ = ∑ *N j*=1 *y*ˆ*j* ) we have <sup>h</sup>Ψ(*a*, *<sup>b</sup>*)|*X*<sup>ˆ</sup> <sup>|</sup>Ψ(*a*, *<sup>b</sup>*)<sup>i</sup> <sup>=</sup> <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>|</sup>Ψ<sup>i</sup> <sup>+</sup> *Na* and <sup>h</sup>Ψ(*a*, *<sup>b</sup>*)|*X*<sup>ˆ</sup> <sup>2</sup> <sup>|</sup>Ψ(*a*, *<sup>b</sup>*)<sup>i</sup> <sup>=</sup> <sup>h</sup>Ψ|*X*<sup>ˆ</sup> <sup>2</sup> <sup>|</sup>Ψ<sup>i</sup> <sup>+</sup> <sup>2</sup>*Na*hΨ|*X*<sup>ˆ</sup> <sup>|</sup>Ψ<sup>i</sup> + (*Na*) 2 , implying that

$$\frac{1}{N} \Delta\_{\hat{\mathbb{X}}}^2 \Big|\_{\Psi(a,b)} = \frac{1}{N} \Delta\_{\hat{\mathbb{X}}}^2 \Big|\_{\Psi}. \tag{A1}$$

Trivially for the momentum operator *P*ˆ*<sup>X</sup>* (and equivalently for *P*ˆ *<sup>Y</sup>* = ∑ *N j*=1 *p*ˆ*y*,*<sup>j</sup>* ) one has

$$\frac{1}{N} \Delta\_{\hat{P}\_{\mathbb{X}}}^{2} \Big|\_{\mathbf{Y}(a,b)} = \frac{1}{N} \Delta\_{\hat{P}\_{\mathbb{X}}}^{2} \Big|\_{\mathbf{Y}'} \tag{A2}$$

i.e., both the position variance and momentum variance are translationally invariant.

For the angular-momentum operator *L*ˆ *<sup>Z</sup>* = ∑ *N j*=1 *x*ˆ*<sup>j</sup> p*ˆ*y*,*<sup>j</sup>* − *y*ˆ*<sup>j</sup> p*ˆ*x*,*<sup>j</sup>* the situation is more interesting. From <sup>h</sup>Ψ(*a*, *<sup>b</sup>*)|*L*<sup>ˆ</sup> *<sup>Z</sup>*|Ψ(*a*, *<sup>b</sup>*)<sup>i</sup> <sup>=</sup> <sup>h</sup>Ψ|*L*<sup>ˆ</sup> *<sup>Z</sup>*|Ψ<sup>i</sup> <sup>+</sup> *<sup>a</sup>*hΨ|*P*<sup>ˆ</sup> *<sup>Y</sup>*|Ψi − *<sup>b</sup>*hΨ|*P*ˆ*X*|Ψ<sup>i</sup> and <sup>h</sup>Ψ(*a*)|*L*<sup>ˆ</sup> <sup>2</sup> *Z* <sup>|</sup>Ψ(*a*)<sup>i</sup> <sup>=</sup> <sup>h</sup>Ψ|*L*<sup>ˆ</sup> <sup>2</sup> *Z* |Ψi + *a* 2 <sup>h</sup>Ψ|*P*ˆ<sup>2</sup> *Y* |Ψi + *b* 2 <sup>h</sup>Ψ|*P*ˆ<sup>2</sup> *X* <sup>|</sup>Ψ<sup>i</sup> <sup>+</sup> *<sup>a</sup>*hΨ|*L*<sup>ˆ</sup> *ZP*ˆ *<sup>Y</sup>* + *P*ˆ *YL*ˆ *<sup>Z</sup>*|Ψi − *<sup>b</sup>*hΨ|*L*<sup>ˆ</sup> *<sup>Z</sup>P*ˆ*<sup>X</sup>* + *P*ˆ*<sup>X</sup> L*ˆ *<sup>Z</sup>*|Ψi − <sup>2</sup>*ab*hΨ|*P*<sup>ˆ</sup> *<sup>Y</sup>P*ˆ*X*|Ψ<sup>i</sup> we have

$$\begin{split} \frac{1}{N} \Delta\_{L\_Z}^2 \Big|\_{\mathbf{Y}(a,b)} &= \frac{1}{N} \Delta\_{L\_Z}^2 \Big|\_{\mathbf{Y}} + a^2 \frac{1}{N} \Delta\_{\hat{\mathbf{P}}\_Y}^2 \Big|\_{\mathbf{Y}} + b^2 \frac{1}{N} \Delta\_{\hat{\mathbf{P}}\_X}^2 \Big|\_{\mathbf{Y}} \\ &+ a \left( \langle \Psi | \hat{L}\_Z \hat{\mathbf{P}}\_Y + \hat{\mathbf{P}}\_Y \hat{L}\_Z | \Psi \rangle - 2 \langle \Psi | \hat{L}\_Z | \Psi \rangle \langle \Psi | \hat{\mathbf{P}}\_Y | \Psi \rangle \right) \\ &- b \left( \langle \Psi | \hat{L}\_Z \hat{\mathbf{P}}\_X + \hat{\mathbf{P}}\_X \hat{L}\_Z | \Psi \rangle - 2 \langle \Psi | \hat{L}\_Z | \Psi \rangle \langle \Psi | \hat{\mathbf{P}}\_X | \Psi \rangle \right) \\ &- 2ab \left( \langle \Psi | \hat{\mathbf{P}}\_Y \hat{\mathbf{P}}\_X | \Psi \rangle - \langle \Psi | \hat{\mathbf{P}}\_Y | \Psi \rangle \langle \Psi | \hat{\mathbf{P}}\_X | \Psi \rangle \right). \tag{A3} \end{split} \tag{A3}$$

Equation (A3) deserves a discussion. In turn, even for the ground state of an interacting many-boson system in a rotationally symmetric [for which <sup>1</sup> *N* ∆ 2 *L*ˆ *Z* = 0 holds] but otherwise translated trap, the angular-momentum variance

$$\frac{1}{N} \Delta\_{\mathbb{Z}\_Z}^2 \Big|\_{\mathbb{Y}(a,b)} = a^2 \frac{1}{N} \Delta\_{\mathbb{P}\_Y}^2 \Big|\_{\mathbb{Y}} + b^2 \frac{1}{N} \Delta\_{\mathbb{P}\_X}^2 \Big|\_{\mathbb{Y}} \tag{A4}$$

differs at the many-body level and mean-field level of theory, i.e., when *a*, *b* 6= 0 and *λ*<sup>0</sup> 6= 0. This is, as can be seen in (A4), because of the respective many-body and mean-field momentum variances, 1 *N* ∆ 2 *P*ˆ*X* and <sup>1</sup> *N* ∆ 2 *P*ˆ *Y* . The analytical result (A4) is employed to analyze the numerical findings for the time-dependent angular-momentum variance in the main text. Generally in the absence of spatial symmetries, see Equation (A3), more terms contribute to the translated angular-momentum variance.

#### **References**


© 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Nonlinear Dynamics of Wave Packets in Tunnel-Coupled Harmonic-Oscillator Traps**

**Nir Hacker <sup>1</sup> and Boris A. Malomed 1,2,\***


**Abstract:** 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 carries 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.

**Keywords:** Bose-Einstein condensates; Josephson oscillations; spontaneous symmetry breaking; Thomas-Fermi approximation; dynamical chaos; ground states; perturbation theory

#### **1. Introduction**

The combination of a harmonic-oscillator (HO) trapping potential and cubic nonlinearity is a ubiquitous setting which occurs in diverse microscopic, mesoscopic, and macroscopic physical settings. A well-known realization is offered by 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 the local refractive index, and the Kerr term is relevant as a model of optical waveguides [15–18]. Models of the same type appear in other physical systems too, such as networks of Josephson oscillators [19].

The character of states created by the interplay of the intrinsic nonlinearity and externally applied trapping potential strongly depends on the sign of the nonlinearity. In the case of the self-attraction (or self-phase-modulation, SPM, in terms of optics [20]), localized modes, similar to solitons, arise spontaneously. On the other hand, self-repulsion tends to 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

**Citation:** Hacker, N.; Malomed, B.A. Nonlinear Dynamics of Wave Packets in Tunnel-Coupled Harmonic-Oscillator Traps. *Symmetry* **2021**, *13*, 372. https://doi.org/ 10.3390/sym13030372

Academic Editor: V. I. Yukalov and Rashid G. Nazmitdinov

Received: 06 February 2021 Accepted: 21 February 2021 Published: 25 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and a weak parabolic potential with an additional spatially periodic one (it represents an optical lattice in BEC [21], or a photonic crystal in optics and plasmonics [22–24]): the interplay of the periodic potential with the self-repulsion gives rise to bright gap solitons [21,25], whose effective mass is negative. For this reason, gap solitons are expelled by the HO potential, but are trapped by the inverted one that would expel modes with positive masses [26].

Another noteworthy feature of the dynamics of one-dimensional (1D) nonlinear fields trapped in confining potentials is the degree of its *nonintegrability*. The generic model for such settings is provided by the nonlinear Schrödinger equation (NLSE, alias the Gross–Pitaevskii equation, in terms of BEC [1–3]) with the cubic term, whose integrability in the 1D space [27] 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 by NLSE with the HO potential term does not lead to establishment of spatiotemporal chaos ("turbulence"), which would be expected in the case of generic nonintegrability [28,29]. 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 to the onset of clearly observed chaos [30,31]. Thermalization of the model with the HO potential was recently explored, in the framework of a stochastically driven dissipative Gross–Pitaevskii equation, in Ref. [32].

The interplay of the cubic nonlinearity and trapping potentials also occurs in twocomponent systems, which represent, in particular, binary BEC [12,33–35]. Here, we aim to consider the system with linear coupling between the 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. [36]. In BEC, different atomic states which correspond to the interacting components can be linearly coupled by a resonant electromagnetic wave [37–40]. 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 [41]. In particular, experiments with temporal solitons in dual-core optical fibers were reported in recent work [42].

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 [35,42]. In the case of repulsive SPM acting in each component and nonlinear repulsion between them (cross-phase modulation, XPM), it was found [33] that the linear mixing shifts the miscibility-immiscibility transition [43] in the trapped condensate. Furthermore, 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 [44], 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) [45].

The objective of the present work is two-fold. First, in Section 2 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 [42]. As a nonchaotic dynamical regime in this case, one may expect Josephson oscillations (JO) of the optical field [42,44,46–49], or of the BEC wave function [50–56], 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 3, 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 induce 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. [57]. 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 parameter *ω* in Equation (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* [58–60], alias *embedded* [61] 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 4.

#### **2. The Symmetric System**

#### *2.1. 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*:

$$\begin{aligned} i\frac{\partial u}{\partial z} + \frac{1}{2}\frac{\partial^2 u}{\partial x^2} + \lambda v - \frac{\Omega^2}{2} x^2 u + \sigma \left( |u|^2 + g|v|^2 \right) u &= 0, \\\\ i\frac{\partial v}{\partial z} + \frac{1}{2}\frac{\partial^2 v}{\partial x^2} + \lambda u - \frac{\Omega^2}{2} x^2 v + \sigma \left( |v|^2 + g|u|^2 \right) v &= 0. \end{aligned} \tag{1}$$

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 Equation (1) are ∼10 µm and 10 ms, respectively.

The system conserves two dynamical invariants., *viz*., the total norm (or power, in terms of optics),

$$P = \int\_{-\infty}^{+\infty} \left[ |u(\mathbf{x})|^2 + |v(\mathbf{x})|^2 \right] d\mathbf{x} \tag{2}$$

and the Hamiltonian (in which Ω = 1 is set),

$$H = \int\_{-\infty}^{+\infty} \left[ \frac{1}{2} \left( \left| \frac{\partial^2 u}{\partial x^2} \right|^2 + \left| \frac{\partial^2 v}{\partial x^2} \right|^2 \right) - \frac{\sigma}{2} \left( |u|^4 + |v|^4 \right) \right. $$

$$ -\sigma g |u|^2 |v|^2 + \frac{1}{2} \mathbf{x}^2 \left( |u|^2 + |v|^2 \right) - \lambda (uv^\* + u^\*v) \right] d\mathbf{x}, \tag{3} $$

where ∗ stands for the complex conjugate. The remaining scaling invariance of Equation (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. [62,63]. Nevertheless, the exact form of the system and its solutions are different.

In the limit of a small amplitude *A*<sup>0</sup> of the input, linearized Equation (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):

$$\begin{aligned} u\_{\rm{IO}}^{\rm{(GS)}}(x,t) &=& A\_0 \exp\left(-\frac{1}{2}x^2 - \frac{1}{2}iz\right) \cos(\lambda z), \\\\ v\_{\rm{IO}}^{\rm{(GS)}}(x,t) &=& iA\_0 \exp\left(-\frac{1}{2}x^2 - \frac{1}{2}iz\right) \sin(\lambda z), \\\\ u\_{\rm{IO}}^{\rm{(DM)}}(x,t) &=& A\_0 \mathbf{x} \exp\left(-\frac{1}{2}x^2 - \frac{3}{2}iz\right) \cos(\lambda z), \\\\ v\_{\rm{IO}}^{\rm{(DM)}}(x,t) &=& iA\_0 \mathbf{x} \exp\left(-\frac{1}{2}x^2 - \frac{3}{2}iz\right) \sin(\lambda z). \end{aligned} \tag{5}$$

Expressions given by Equations (4) and (5) at *z* = 0 are used below as inputs in simulations of the full nonlinear Equation (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 Equations (1) and (10), and the relaxation algorithm for finding solutions to stationary equations, such as Equations (12) and (13), see below.

#### *2.2. The Transition from Regular to Chaotic Dynamics*

Increase of amplitude *A*<sup>0</sup> 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 Equation (1) with *g* = 0 (no XPM interaction), in interval 0 < *z* < *Z*, with the DM input, taken as per Equation (5) at *z* = 0, is presented in Figure 1. In particular, the left bottom panel of the figure displays oscillations of peak intensities of the fields,

$$\left\{\mathsf{U}\_{\max}^{2}(z), V\_{\max}^{2}(z)\right\} \equiv \max\_{\mathsf{x}} \left\{|u(\mathsf{x}, z)|^{2}, |v(\mathsf{x}, z)|^{2}\right\},\tag{6}$$

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:

$$\{P(\kappa), Q(\kappa)\} = \int\_0^Z \left\{\mathcal{U}\_{\text{max}}^2(z), V\_{\text{max}}^2(z)\right\} \exp(-i\kappa z) dz,\tag{7}$$

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 during the simulations).

The regularity of the dynamical regime displayed in Figure 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 Figure 1. It can be estimated as *δκ* = 2*π*/*Z* ≈ 0.06. Note the overall symmetry between the two components in Figure 1 (in particular, their spectra are identical in the bottom panel).

**Figure 1.** A typical example of a regular Josephson dynamical regime, initiated by the DM (dipole mode) input launched in the *u* component (as given by Equation (5) with *z* = 0 and *A*<sup>0</sup> = 1). The solution is produced by simulations of Equation (1) with *λ* = *σ* = Ω = 1, *g* = 0. Plots in the top row display the evolution of components *u*(*x*, *z*) and *v*(*x*, *z*). Left bottom: The evolution of the peak intensities of both components, *U*<sup>2</sup> max(*z*) ≡ max *x* n |*u*(*x*, *z*)| 2 o and *V* 2 max(*z*) ≡ max *x* nh|*v*(*x*, *<sup>z</sup>*)<sup>|</sup> 2 io. Right bottom: The power spectrum of oscillations of the two components, defined as per Equation (7). The spectra are virtually identical for both components.

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 Figure 2. Note that both the regular and chaotic dynamical pictures displayed in Figures 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 Equation (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*<sup>0</sup> = 1 and 4 initiate, severally, quasi-regular and chaotic evolution.

**Figure 2.** The same as in Figure 1, but for a typical example of a chaotic dynamical regime, initiated by the DM input (5) with a larger amplitude, *A*<sup>0</sup> = 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 per Equations (4) and (5) at *z* = 0, are summarized in charts plotted in the plane of *λ*, *A* 2 0 in Figure 3. They display heatmaps of values of the parameter which quantifies the sharpness of the central peak in the spectrum of the dynamical state:

$$\text{Sharpness} \equiv \frac{\int\_{\text{FWHM}} |P(\kappa)|^2 d\kappa}{\int\_0^\infty |P(\kappa)|^2 d\kappa},\tag{8}$$

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)| <sup>2</sup> = (1/2)(*P*(*κ*))max. Values of Sharpness close to <sup>1</sup> imply the domination of a single sharp peak, such as one in Figure 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., Figure 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 Figure 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 Figure 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.

**Figure 3.** Heatmaps of values of sharpness (8) of the central spectral peak quantifying proximity of the system's dynamics to the regular regime. The maps are plotted in the plane of the linear-coupling strength, *λ*, and intensity of the input, *A* 2 0 , which is launched in one component. (**Top left**): The GS input, given by Equation (4) at *z* = 0, in the case of the self-attraction (*σ* = +1). (**Top right**): The same, but in the case of self-repulsion (*σ* = −1). (**Bottom left**): The same as in the top left panel, but produced by the DM input, given by Equation (5) at *z* = 0. (**Bottom right**): The same as in the bottom left panel, but in the case of self-repulsion (*σ* = −1). In all cases, *g* = 0 is set in Equation (1) (no XPM interaction between the components). Black curves cutting the left panels in their lower areas designate the onset of SSB (spontaneous symmetry breaking), signalized by appearance of *θ* 6= 0, see Equation (9).

Lastly, the comparison of the left and right panels in Figure 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 [64]. 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 Figure 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 abovementioned integrability of such a system of NLSEs with the linear coupling [45] (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 Equation (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 <sup>0</sup> & 200, being limited to *λ* . 0.15, cf. the top left panel in Figure 3.

#### *2.3. Spontaneous Symmetry Breaking (SSB) between the Coupled Components*

A noteworthy feature of the dynamical state presented in Figure 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 [41]. In particular, SSB of stationary states in systems with the HO trapping potential acting in both cores was addressed in Refs. [35,56], while Figure 2 demonstrates the symmetry breaking in the *dynamical* JO state.

The SSB effect may be quantified, as usual [44], by asymmetry of the dynamical states, which we define as

$$\Theta \equiv \frac{\int\_0^\infty |P(\kappa)|^2 d\kappa - \int\_0^\infty |Q(\kappa)|^2 d\kappa}{\int\_0^\infty |P(\kappa)|^2 d\kappa + \int\_0^\infty |Q(\kappa)|^2 d\kappa}. \tag{9}$$

The SSB occurs as a transition from Θ = 0 to Θ 6= 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 [41,44], while here we consider it in the dynamical setting. On the other hand, in the case of the repulsive nonlinearity, *σ* = −1 in Equation (1), SSB of stationary states takes place in this system only at *g* > 1 [56], 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 Figure 3, respectively. The SSB bifurcation of the dynamical states under consideration is of the *supercritical*, alias *forward*, type [65], in terms of dependence Θ *A* 2 0 , as shown in Figure 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 *λ*.

**Figure 4.** Plots of the SSB (spontaneous symmetry breaking) in the dynamical states initiated by the GS (ground state) and DM (dipole mode) inputs: the asymmetry parameter, defined as per Equation (9), is plotted versus the intensity of the input, *A* 2 0 , for three different fixed values of the linear-coupling constant, *λ*, as indicated in the figure.

It is worth noting that, as clearly shown by the SSB boundary (the black line) in the top left panel of Figure 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 [41,44]. On the other hand, the bottom left panel in Figure 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 Figure 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.

#### **3. The Half-Trapped System**

The asymmetric system of linearly coupled NLSEs, with the HO potential included in one equation only, is written as

$$\begin{aligned} i u\_z + \frac{1}{2} u\_{xx} + \lambda v - \frac{1}{2} x^2 u + \sigma \Big( |u|^2 + \mathbf{g} |v|^2 \Big) u &= -\omega u, \\\\ i v\_z + \frac{1}{2} v\_{xx} + \lambda u + \sigma \Big( |v|^2 + \mathbf{g} |u|^2 \Big) v &= 0, \end{aligned} \tag{10}$$

cf. Equation (1). Here, as said above, Ω = 1 is set in the first equation, and the propagationconstant 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 Equation (10) are looked for as

$$\{\mu, \upsilon\} = \{\mathcal{U}(\mathbf{x}), V(\mathbf{x})\} \exp(-i\mu z),\tag{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

$$(\mu + \omega)\mathcal{U} + \frac{1}{2}\frac{d^2\mathcal{U}}{dx^2} + \lambda\mathcal{V} - \frac{1}{2}\mathfrak{x}^2\mathcal{U} + \sigma\Big(\mathcal{U}^2 + \mathcal{g}\mathcal{V}^2\Big)\mathcal{U} = 0,\tag{12}$$

$$
\mu V + \frac{1}{2}\frac{d^2V}{dx^2} + \lambda U + \sigma \left(V^2 + gL^2\right)V = 0.\tag{13}
$$

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.

#### *3.1. The Linearized System: Analytical and Numerical Results*

#### 3.1.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 to the linearized version of Equations (12) and (13), with zero *v* component, are

$$\mathcal{U}\_{\rm GS}^{(0)}(\mathbf{x}) \;=\; \frac{1}{\pi^{1/4}} \exp\left(-\frac{\mathbf{x}^2}{2}\right) , V\_{\rm GS}^{(0)} = 0 ,\tag{14}$$

$$M\_{\rm DM}^{(0)}(\mathbf{x}) \quad = \ \frac{\sqrt{2}}{\pi^{1/4}} \mathbf{x} \exp\left(-\frac{\mathbf{x}^2}{2}\right), \\ V\_{\rm DM}^{(0)} = \mathbf{0}, \tag{15}$$

where the pre-exponential constants are determined by the normalization condition,

$$\int\_{-\infty}^{+\infty} \left[ \mathcal{U}^2(\mathbf{x}) + V^2(\mathbf{x}) \right] d\mathbf{x} = 1,\tag{16}$$

which we adopt in this section. The eigenvalues corresponding to eigenmodes (14) and (15) are

$$
\mu\_{\rm GS}^{(0)} = 1/2 - \omega; \; \mu\_{\rm DM}^{(0)} = 3/2 - \omega. \tag{17}
$$

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),

$$i v\_z + \frac{1}{2} v\_{\rm xx} = -\lambda \mathcal{U}\_{\rm GS,DM}^{(0)}(\boldsymbol{x}) \exp\left(-i \mu\_{\rm GS,DM}^{(0)} \boldsymbol{z}\right). \tag{18}$$

Obviously, Equation (18) gives rise to emission of propagating waves ("radiation"), in the form of *<sup>v</sup>* <sup>∼</sup> *<sup>λ</sup>* exp *ikx* − *iµ* (0) GS,DM*z* , at resonant wavenumbers *k* = ± q 2*µ* (0) GS,DM, provided that *µ* (0) GS,DM is positive, i.e.,

$$v\_{\rm rad} \sim \lambda \exp\left(\pm i \sqrt{2\mu\_{\rm GS,DM}^{(0)}} \left(x - V\_{\rm ph} z\right)\right),\tag{19}$$

where the phase velocity is

$$V\_{\rm ph} = \frac{k}{2} \equiv \pm \frac{1}{2} \sqrt{2\mu\_{\rm GS,DM}^{(0)}}\tag{20}$$

(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*| ≡ 2*V*ph.

An illustration of this dynamics is presented in Figure 5. Straight red lines designate the wave-propagation directions, which exactly agree with the phase velocity predicted by Equation (20), and the expansion of the area occupied by the radiation complies with the prediction based on the group velocity.

**Figure 5.** Simulations of the evolution of the linearized half-trapped system (10), displayed in the untrapped component by plotting Re(*v*(*x*, *z*)). (**Left**): Emission of radiation generated by the GS (ground state) populating the trapped component, *u*(*x*, *z*) (see Equation (14)), with *ω* = 0.25 in Equation (10). (**Right**): The same, but for the radiation generated by the DM (dipole mode) in the trapped component (see Equation (15)), with *ω* = 0. In both cases, the linear-coupling constant is *λ* = 0.05.

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 Equation (16). An example of the decay is displayed in Figure 6, for a small initial amplitude of the GS input, *A*<sup>0</sup> = 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*<sup>0</sup> 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 Figure 7.

3.1.2. The Shift of the GS and DM Existence Thresholds at Small Values of Coupling Constant *λ*

At *µ* (0) GS,DM < 0 the radiation is not generated by Equation (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 Equations (12) and (13), such that the bound states of the GS and DM types exist at

$$
\omega > (\omega\_{\text{GS,DM}})\_{\text{thr}} \tag{21}
$$

respectively. Obviously, (*ω*GS)thr = 1/2 and (*ω*DM)thr = 3/2 in the limit of *λ* = 0, see Equation (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 Equations (14) and (15), respectively. With small *λ*, the first-order solution for *V*(*x*), *viz*., *V*(*x*) ≡ *λV* (1) GS,DM(*x*), must be found from the inhomogeneous equation, which follows from Equation (13), in which *µ* = 0 is set:

$$\frac{d^2}{d\mathfrak{x}^2} V\_{\rm GS,DM}^{(1)} = -2L I\_{\rm GS,DM}^{(0)}(\mathfrak{x}).\tag{22}$$

Straightforward integration of Equation (22), with expressions (14) and (15) substituted on the right-hand side, yields

$$V\_{\rm GS}^{(1)}(\mathbf{x}) = -\sqrt{2}\pi^{1/4} \left[ \frac{\mathbf{x}}{\sqrt{2}} \text{erf}\left(\frac{\mathbf{x}}{\sqrt{2}}\right) + \sqrt{\frac{2}{\pi}} \exp\left(-\frac{\mathbf{x}^2}{2}\right) \right]. \tag{23}$$

$$V\_{\rm DM}^{(1)}(\mathbf{x}) = 2\pi^{1/4} \text{erf}\left(\frac{\mathbf{x}}{\sqrt{2}}\right),\tag{24}$$

where erf(*x*) is the standard error function, which is an odd function of *x*.

**Figure 6.** (**Left**): The evolution of fields *u* and *v* in the half-trapped system, as produced by simulations of Equation (10) with *ω* = 0 and coupling constant *λ* = 1, initiated by the DM input in the *u* core, taken as per Equation (14) with *A*<sup>0</sup> = 0.1. Here and in Figure 7, spurious left-right asymmetry of the radiation field in the *v* component is an illusion produced by plotting. (**Right**): The evolution of the peak intensities of both components, max *x* n |*u*(*x*, *t*)| 2 o and max *x* n |*v*(*x*, *t*)| 2 o .

**Figure 7.** The same as in Figure 6, but for a large amplitude of the DM input in the *u* component, *A*<sup>0</sup> = 8.

Next, the small perturbation in Equation (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 [66], in the first approximation of the perturbation theory, when *V*(*x*) is replaced by expression (23) or (24), the result is

$$
\delta\mu\_{\rm GS,DM} = -\lambda^2 \int\_{-\infty}^{+\infty} V\_{\rm GS,DM}^{(1)}(\mathbf{x}) \mathcal{U}\_{\rm GS,DM}^{(0)}(\mathbf{x}) d\mathbf{x} \equiv -I\_{\rm GS,DM} \lambda^2 \tag{25}
$$

with coefficients

$$I\_{\rm GS} = -3.414 \,\text{J}\_{\rm DM} = 4 \,\text{} \tag{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

$$(\omega\_{\rm GS})\_{\rm thr} = 1/2 - I\_{\rm GS} \lambda^2,\ (\omega\_{\rm DM})\_{\rm thr} = 3/2 - I\_{\rm DM} \lambda^2. \tag{27}$$

The analytical predictions are compared to numerical results, obtained from a solution of the linearized variant of Equations (12) and (13), in Figure 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.

**Figure 8.** The analytically predicted and numerically found threshold values of the mismatch parameter, *ω*, above which the GS and DM solutions (the left and right panels, respectively) are produced, for the half-trapped system, by the linearized version of Equations (12) and (13), vs. coupling constant *λ*. The analytical results are produced by Equation (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 *λ*.

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 Figure 9. The existence of the DM at this point agrees with the right panel of Figure 8.

**Figure 9.** A bound state of the DM (dipole mode) type in the half-trapped system, found as a numerical solution of Equations (12) and (13) with *ω* = 1.4, *λ* = 0.225, and *σ* = 0 (the linearized version). The eigenvalue corresponding to this solution is *µ* ≈ −0.036.

3.1.3. The Analysis for Large Values of *λ*

In the opposite limit of large coupling constant *λ*, an analytical approximation for 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 Equation (13), to obtain *V* ≈ −(*λ*/*µ*)*U*. Then, substituting this relation back into the originally neglected derivative term, one obtains a necessary correction to this relation:

$$V \approx -\frac{\lambda}{\mu} \mathcal{U} + \frac{\lambda}{2\mu^2} \frac{d^2 \mathcal{U}}{d\mathbf{x}^2}. \tag{28}$$

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:

$$\left(\mu - \frac{\lambda^2}{\mu} + \omega\right) \mathcal{U} + \frac{1}{2} \left(1 + \frac{\lambda^2}{\mu^2}\right) \frac{d^2 \mathcal{U}}{d\lambda^2} - \frac{1}{2} \mu^2 \mathcal{U} = 0. \tag{29}$$

Then, the standard solution for the quantum-mechanical HO yields an equation which determines the spectrum of the eigenvalues:

$$
\mu - \frac{\lambda^2}{\mu} + \omega = \left(\frac{1}{2} + n\right)\sqrt{1 + \frac{\lambda^2}{\mu^2}}\tag{30}
$$

where *n* = 0, 1, 2, . . . is the quantum number. Taking into regard that *λ* is now a large parameter, Equation (30) produces a final result for the spectrum,

$$
\mu \approx -\lambda - \frac{\omega}{2} + \frac{1}{\sqrt{2}} \left( \frac{1}{2} + n \right). \tag{31}
$$

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, Equation (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 Equation (31) with *n* = 0 and *n* = 1, respectively, is compared to numerically found counterparts in Figure 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.

**Figure 10.** ( **The top row**): left and right panels display the analytically predicted eigenvalues given by Equation (31) with *n* = 0 and *n* = 1, for the ground state (GS) and dipole mode (DM), respectively, of the half-trapped system, and their counterparts produced 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 connecting points (*λ*, *ω*) = (0.10) and (10, 0). The results shown in these plots are relevant for relatively large values of *λ*.

3.1.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 Equations (10), (12) and (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*:

$$\mathcal{U}(\mathbf{x}) = \mathcal{U}\_0 \left[ \left( \lambda^2 - \frac{1}{2} \right) + \mathbf{x}^2 \right] \exp \left( -\frac{\mathbf{x}^2}{2} \right),$$

$$V(\mathbf{x}) = -2\lambda \mathcal{U}\_0 \exp \left( -\frac{\mathbf{x}^2}{2} \right),$$

$$\mathcal{U}\_0^2 = \pi^{-1/2} \left( \lambda^4 + 4\lambda^2 + 1/2 \right)^{-1},$$

$$\omega = \frac{9}{4} - \frac{\lambda^2}{2}, \; \mu = \frac{1}{2} \left( \lambda^2 + \frac{1}{2} \right), \tag{32}$$

with amplitude *U*<sup>0</sup> defined by condition *P* = 1, see Equation (16). We stress that, as seen in Equation (32), this exact solution may only have *µ* > 0, i.e., it is BIC (a *bound state in the continuous spectrum* [58–60]), alias an *embedded mode* [61]. 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 Equation (16)), the solution is

$$\mathcal{U}(\mathbf{x}) = \mathcal{U}\_0 \mathbf{x} \left[ \left( \lambda^2 - \frac{3}{2} \right) + \mathbf{x}^2 \right] \exp \left( -\frac{\mathbf{x}^2}{2} \right),$$

$$V(\mathbf{x}) = -2\lambda \mathcal{U}\_0 \mathbf{x} \exp \left( -\frac{\mathbf{x}^2}{2} \right),$$

$$\mathcal{U}\_0^2 = 2\pi^{-1/2} \left( \lambda^4 + 4\lambda^2 + 3/2 \right)^{-1},$$

$$\omega = \frac{11}{4} - \frac{\lambda^2}{2}, \; \mu = \frac{1}{2} \left( \lambda^2 + \frac{3}{2} \right), \tag{33}$$

which also exists at a single value of *ω*, and with a single eigenvalue, *µ* > 0. Note that both exact solutions, given by Equations (32) and (33), may exists at positive and negative values of *ω*, as well as at *ω* = 0 (in Equations (32) and (33), *ω* = 0 at *λ* <sup>2</sup> = 9/2 and *λ* <sup>2</sup> = 11/2, respectively).

These exact solutions for BIC states in the two-component systems are somewhat similar to those found in Ref. [67], 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 Equation (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 Figure 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 Figure 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).

**Figure 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 *σ* = −1 in the bottom one. Other parameters are *g* = 0, *λ* = 7/2 and *ω* = 1/2, which are related as per Equation (32).

It may happen that the linearized system of Equations (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.

The 2D extension of the linearized form of Equation (10) is

$$\begin{aligned} i u\_z + \frac{1}{2} (u\_{xx} + u\_{yy}) + \lambda v - \frac{1}{2} \left(\mathbf{x}^2 + y^2\right) u &= -\omega u\_r \\\\ i v\_z + \frac{1}{2} (v\_{xx} + v\_{yy}) + \lambda u &= 0 \end{aligned} \tag{34}$$

Although 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 [35,68]. In polar coordinates (*r*, *θ*), particular exact solutions of linear Equation (34), with all integer values of the vorticity, *S* = 0, 1, 2, . . . , are found as

$$u = \mathcal{U}\_0 r^S \left[ \left( \lambda^2 - 1 - S \right) + r^2 \right] \exp \left( -i\mu z + iS\theta - \frac{r^2}{2} \right),$$

$$v = -2\lambda \mathcal{U}\_0 r^S \exp \left( -i\mu z + iS\theta - \frac{r^2}{2} \right),$$

$$\omega = \frac{1}{2} \left( 5 + S - \lambda^2 \right), \; \mu = \frac{1}{2} \left( \lambda^2 + 1 + S \right), \tag{35}$$

with arbitrary amplitude *U*<sup>0</sup> (the 2D solution of the GS type corresponds to *S* = 0 in Equation (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.

#### *3.2. The Nonlinear Half-Trapped System*

#### 3.2.1. The Thomas–Fermi Approximation (TFA)

In the presence of the self-defocusing nonlinearity, *σ* = −1 in Equations (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, Equation (13) with *g* = 0 (the XPM coupling is omitted here) is solved as

$$
\mathcal{U} = (V/\lambda) \Big(V^2 - \mu\Big),
\tag{36}
$$

and Equation (12) amounts to an algebraic equation for the squared amplitude *W* ≡ −*V* 2/ *µ* > 0, *viz*.,

$$m\mathcal{W}(\mathcal{W}+1)^3 + \mathfrak{F}^2\mathcal{W} = \mathfrak{f}\_0^2 - \mathfrak{f}^2,\tag{37}$$

where

$$m \quad \equiv \quad -\frac{2\mu}{\lambda^2} > 0,\ \sharp \equiv -\frac{\mathfrak{x}}{\mu'} \tag{38}$$

$$\mathfrak{E}\_0^2 \equiv -\frac{2}{\mu^3} \left[ \lambda^2 - \mu(\omega + \mu) \right],\tag{39}$$

and the applicability condition for TFA is easily shown to be *m* 1. The TFA solution exists under condition *ξ* 2 <sup>0</sup> > 0 (see Equation (39)), i.e., in a finite interval (*bandgap*) of values of the propagation constant,

$$0 < -\mu < \sqrt{\lambda^2 + \omega^2/4} + \omega/2 \equiv -\mu\_0. \tag{40}$$

Outside of the bandgap, i.e., at *µ* < *µ*<sup>0</sup> < 0, the TFA solution does not exist. In the bandgap, Equation (37) produces a usual GS profile, with a single value of *W* corresponding to each *ξ* 2 from region *ξ* <sup>2</sup> < *ξ* 2 0 (see Figure 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 *ξ* <sup>2</sup> > *ξ* 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,69].

**Figure 12.** A typical example of the GS solution predicted by TFA (Thomas–Fermi approximation) for the half-trapped system, as per Equations (36)–(39), for *σ* = −1, *g* = 0, *λ* = 8, *µ* = −4, and its comparison to the numerically found counterpart. The respective value of parameter *m* (see Equation (38)), which should be small for the applicability of TFA, is *m* = 0.125.

Note that in the limit of large *λ*, the bandgap's width, as given by Equation (40), is −*µ*<sup>0</sup> ≈ *λ* + *ω*/2, which is close to the largest value of −*µ* predicted by Equation (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

$$0 \quad < \quad \delta\mu \equiv \mu - \mu\_0 \ll -\mu\_0 \tag{41}$$

$$
\xi\_0^2 \approx -\left(4/\mu\_0^3\right)\sqrt{\lambda^2 + \omega^2/4}\delta\mu\tag{42}
$$

In this case, Equations (37) and (36) simplify to

$$\mathcal{U}^2 \approx \begin{cases} \ \left(\mu\_0^2/2\right)\left(\xi\_0^2 - \xi^2\right), \text{ at } \xi^2 < \xi\_0^2\\ \ 0 \text{ at } \xi^2 > \xi\_{0'}^2 \end{cases} \tag{43}$$

$$V^2 \approx \begin{cases} \ (\lambda^2/2) \left(\xi\_0^2 - \xi^2\right), \text{ at } \xi^2 < \xi\_{0'}^2\\ \ 0 \text{ at } \xi^2 > \xi\_0^2. \end{cases} \tag{44}$$

Expressions (42)–(44) make it possible to calculate the total power of the GS (see Equation (2)),

$$P \approx \frac{16}{3} \left(\lambda^2 + \mu\_0^2\right) \left(\lambda^2 + \frac{\omega^2}{4}\right)^{3/4} \left(-\mu\_0\right)^{-7/2} \left(\delta\mu\right)^{3/2}.\tag{45}$$

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 [70], 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 [71–73]).

#### 3.2.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 the coupled half-trapped system, shown in Figure 8, suggests identifying 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 Equations (12) and (13) for spatially even and odd modes with the fixed total power, *P* = 1 (as per Equation (16)), while the linearcoupling 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 Figure 13 (for the GS) and Figure 14 (for the DM), which show threshold values of the mismatch parameter, defined as per Equation (21).

**Figure 13.** The heatmap of threshold values of the mismatch parameter, (*ω*GS)thr, in the half-trapped system, based on Equations (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 Equation (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.

**Figure 14.** The same as in Figure 13, but for DMs (dipole modes), obtained as numerical solutions to Equations (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.

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 Figure 13, and its counterpart for the DMs is located between red lines in Figure 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 Figure 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 the left panel of Figure 8 and Equation (27). Furthermore, Figure 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 small values of *λ*, while the self-repulsion (*σ* < 0) can 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 Figure 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 Equation (27) and the right panel of Figure 8. 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 SPM conspicuously reduces the remaining nontrivial area.

#### **4. Conclusions**

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 (harmonicoscillator) 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 (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 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.

**Author Contributions:** Formulation of the model and analytical considerations, B.A.M.; numerical computations, N.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Israel Science Foundation, grant number 1286/17.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


*Article*
