text:
# Collective modes, chaotic behavior and self-trapping in the dynamics of three coupled Bose-Einstein condensates

## I Introduction

A remarkable progress in the experimental design has been done since the first direct observation of Bose-Einstein condensation in diluite atomic gas PEXP . One of most promising development concerns the construction of experimental devices in which condensates, achieved within complex geometries, interact with each other giving rise to quantum effects that are observable at the macroscopic level COLL  -INGU . In this respect, one should recall, for example, the (superfluid) boson Josephson-junction arrays obtained by means of optical lattices that trap weak interacting BECs in periodic arrays of potential wells JJC . In parallel with the experimental work, an increasing attention has been devoted theoretically to study the dynamical behavior of low-energy states in arrays of Bose-Einstein condensates (BEC) where the number of lattice sites (namely the potential wells occupied by the condensate) is very large LTR ; TROMB .

Opposite situations, corresponding to ‘lattices’ formed by two or three interacting wells, have been investigated as well in various recent papers (see Refs. MILB  -FP2 ). In particular, the two-well system (dimer) has been analysed thorougly from both the semiclassical (mean-field) and the purely quantum viewpoint in Refs. SM1 ; SM2  and FPZ1 ; MACZ ; FP2 , respectively. Such investigations have revealed how the nontrivial structure of dimer phase space causes many significant phenomena such as the symmetry-breaking effect (issuing oscillation modes that are isoenergetic but inequivalent), the onset of $`\pi `$-phase oscillations, the self-trapping of boson populations, and, quantum-mechanically, the occurrence of (parameter-dependent, nondegenerate) doublets in the energy spectrum entailing periodic self-trapping.

The connection between the quantum and the semiclassical picture of many interacting bosonic wells has been illustrated in Ref. FPZ1  and, in view of the closed link between an array of interacting BECs and the Bose-Hubbard (BH) model FPZ1 ; CIR , in Ref. AP  within the BH model theory. We wish to observe that the semiclassical approach (corresponding to describing condensates within the Bogoliubov approximation) is appropriate for interacting wells with macroscopic boson populations. However, the systems recently obtained, where condensate is distributed among many wells are able to provide mesoscopic numbers of bosons per wells INGU , might be best modeled by using the space-mode approximation MILB ; PW  stemmed from the second-quantized boson field theory.

Compared with the dimer nonlinear behavior, indeed the system of three interacting bosonic wells (trimer model) exhibits a new phenomenology and provides a novel research topic, both theoretically and experimentally. In fact, despite the simplicity of the model, the trimer dynamics is affected by a strong inner instability leading to the chaos onset. This feature originates from the combination of the model nonlinearity with the nonintegrable nature that distinguishes trimer dynamics within extended regions of its phase space. In respect of this, the solutions of the trimer dynamical equations \[displaying nonlinear periodic oscillations (dimerlike orbits) with possible self-trapping effect\] recognized in Refs. NEM ; FP3  must be viewed as special subregimes in which dynamics is integrable. Concerning the instability of trimer dynamics, we wish to recall that the presence of homoclinic chaos has been revealed for the asymmetric trimer in Ref. hg , while, at the quantum level, the survival of breather configurations has been investigated in Ref. ff .

In the present paper we perform a systematic analysis of the trimer dynamics directed to ascertain that the nonintegrability generates chaos and that this dominates the interactions of three BECs modeled in a standard semiclassical picture. A rich scenario of dynamical behaviors emerges from our analysis, which confirms the extremely structured character of trimer dynamics and represents the natural prosecution of the work of Ref. FP3  in which we focused our attention on the special dimerlike (integrable) regime of trimer and on the self-trapping effect.

Indeed the analysis we perform appears to be topical in relation to the study of the dynamics of solitons TROMB  and of vortices CAPE  on one-dimesional chains of BECs as well as of experimental architectures obtained recently INGU . The observation of such chain excitations, in fact, must take into account the possible destructive action of inner instabilities whose influence is clearly manifested in the trimer chain. We emphasize the fact that the trimer chain is the simplest possible situation in which interacting BECs turns out to be governed by nonintegrable equations. The dimer dynamics, in fact, is completely integrable.

The paper layout is the following. In Sec. II we review the derivation of the space-mode Hamiltonian for three coupled wells from the quantum field theory of bosonic fluids in the dilute-gas approximation, and present the semiclassical picture that describes coupled boson wells with macroscopic populations. In Sec. III we identify the set of fixed points of the trimer Hamiltonian equations and show how such points are associated to periodic solutions (representing collective modes) owing to the dynamically conserved total boson number $`N`$. Sec. IV is devoted to study the second variation (with $`N=const`$) of the energy function for such extremal configurations in order to make explicit their nature. In Sec. V we perform a dimensional reduction of trimer dynamics by defining a new set of canonical variables whereby the constant of motion $`N`$ is incorporated in the dynamical equations. This paves the way to the implementation of the Poincaré sections’ method. The chaos onset in the trimer dynamics is investigated in Sec. VI. First, in a qualitative way, by constructing the Poincaré sections for trajectories whose initial conditions are chosen in proximity of the fixed-point configurations identified in Sec. III. Then, quantitatively, by measuring the maximum Lyapunov exponent of such trajectories. Sec. VII contains concluding remarks and comments on future work.

## II Trimer dynamics

The model for a chain (or more complex structures) of $`M`$ interacting bosonic wells can be derived from the quantum field theory for boson fluids (with a nonlocal $`\psi ^4`$ term) by implementing the M-(space)mode approximation CIR . If the boson fluid is a diluite gas of $`N`$ interacting bosons trapped in an external potential $`V_e`$ then its dynamics is generated by the local boson-field Hamiltonian PW ; LEG ; DGPS

$`\widehat{H}={\displaystyle 𝑑𝐫\widehat{\psi }^{^{}}(𝐫)\left[V_e\frac{_{\mathrm{}^2^2}}{^{2m}}+\frac{_{U_0}}{^2}\widehat{\psi }^{^{}}(𝐫)\widehat{\psi }(𝐫)\right]\widehat{\psi }(𝐫)}`$
where $`m`$ is the boson mass, $`U_0=4\pi \mathrm{}^2a/m`$ takes into account the two-body interaction strength, and $`a`$ is the $`s`$wave scattering length. The field $`\psi (𝐫)`$ \[$`\widehat{\psi }^+(𝐫)`$\] is the Heisenberg field operator that destructs \[creates\] bosons at position $`𝐫`$. The nonlinear term has been written in the usual normal ordered form. In order to work out the trimer Hamiltonian we state some assumptions: First, $`V_e`$ is a three-well symmetric-shaped potential. Second, the lowest energy level of each well (within the approximation of $`V_e`$ in terms of a single-well parabolic potential) must be well separated from the higher energy levels PW . Third, the binary particle interactions is not strong enough to significatively change the latter assumptions.

Some furher approximations are necessary to make explicit the ground-state structure in $`\widehat{H}`$. Let $`𝐫_i`$ ($`i=1,2,3`$) be the locations of the minima of $`V_e`$ and let $`V_j=V(𝐫𝐫_j)`$ be the parabolic approximation to the potential in the $`j`$th minimum, so that $`V_e(𝐫)V(𝐫𝐫_j)`$ when $`𝐫𝐫_j`$. Also, let us introduce the eigenstates $`u_j`$ that represent the normalized single-boson ground-states with energy $`E_0`$ of the local parabolic potential $`V_j`$. These states are only approximatively orthogonal because of $`d^3𝐫\overline{u}_ju_k=\delta _{jk}+R_{jk}`$ but the residue $`R_{jk}`$ is a quantity exponentially suppressed depending on the overlap between $`u_j`$ and $`u_k`$. The analysis is restricted to those potentials for which $`R_{jk}<<1`$ thus making negligible such contributions.

The picture of the system thus resulting suggests that the boson field operator can be expanded as $`\widehat{\psi }(𝐫,t)=\mathrm{\Sigma }_i\overline{u}_i(𝐫)\widehat{a}_i(t)`$, where $`\widehat{a}_i(t)`$ \[$`\widehat{a}_i^+(t)`$\] is the annihilation \[creation\] boson operator (associated to the space-mode state $`u_i`$) that satisfy the commutation relation $`[\widehat{a}_i(t),\widehat{a}_k^+(t)]=\delta _{ik}`$. By substituting this expression in the many-body hamiltonian one can obtain, to the lowest order in the overlap between the single-well modes, the quantum trimer Hamiltonian FPZ1 ; NEM

$$H_3=\mathrm{\Sigma }_i[U(n_i1)v]n_i\frac{_T}{^2}\mathrm{\Sigma }_{i,j}\left(a_i^{^{}}a_j+a_j^{^{}}a_i\right),$$

where the site index $`i,j=1,2,3`$, $`v=E_0`$, and the operators $`n_ia_i^{}a_i`$ count the number of bosons at site $`i`$. Also, in $`H_3`$

$$T=2d^3𝐫\overline{u}_j[V_jV_e]u_{j\pm 1},U=\frac{_{U_0}}{^2}d^3𝐫|u_j|^4,$$

represent the (interwell) hopping amplitude and the strength of the Coulomb on-site repulsion, respectively MILB ; PW .

The quantum dynamics involved by Hamiltonian $`H_3`$ can be cast in a classical form by representing the system quantum state through a trial state $`|Z`$ written in terms of Glauber’s states $`|z_i`$ (defined by $`a_i|z_i=z_i|z_i`$). By implementing the time-dependent variational principle (TDVP) and the procedure discussed in Refs. AP , ZFG  on $`|Z=\mathrm{\Pi }_i|z_i`$, one obtains the effective Hamiltonian FPZ1

$$_3(Z,Z^{}):=Z|H_3|Z$$

$$\mathrm{\Sigma }_{j=1}^3[U|z_j|^4v|z_j|^2\frac{_T}{^2}(z_j^{}z_{j+1}+c.c.)],$$
(1)
($`j=14`$ on the trimer chain) with the equations

$`i\mathrm{}\dot{z}_1`$ $`=`$ $`(2U|z_1|^2v)z_1{\displaystyle \frac{_T}{^2}}(z_2+z_3),`$
$`i\mathrm{}\dot{z}_2`$ $`=`$ $`(2U|z_2|^2v)z_2{\displaystyle \frac{_T}{^2}}(z_3+z_1),`$ (2)
$`i\mathrm{}\dot{z}_3`$ $`=`$ $`(2U|z_3|^2v)z_3{\displaystyle \frac{_T}{^2}}(z_1+z_2).`$
Such equations for $`z_j`$ (notice that $`z_jZ|a_j|Z`$ and $`z_j^{}=Z|a_j^+|Z`$) can be calculated from $`_3`$ via the Poisson brackets $`\{z_k^{},z_j\}=i\delta _{k,j}/\mathrm{}`$ furnished by the TDVP method. Those for $`z_j^{}`$ are easily obtained by complex conjugation. Various aspects concerning the special dimeric subregime of Eqs. (2), in which trimer dynamics is integrable, have been studied in Refs. NEM , FP3 .

## III Fixed points and periodic orbits

The distinctive features that characterize the dynamics of a given Hamiltonian system can be recognized by exploring the structure of its phase space. The first step to do this is to locate the fixed points. The latter are derived in the present section based on Eqs. (2), whereas the nature of such points is studied in the next section.

The fixed-point equations for the trimer are obtained by setting $`\dot{z}_j0`$ in Eqs. (2). Since the trimer dynamics is constrained by the constant of motion $`N=|z_1|^2+|z_2|^2+|z_3|^2`$, one must incorporate explicitly the restriction to the phase-space submanifold defined by $`N=const`$ through a Lagrange multiplier $`\chi `$. This just requires that one considers the variations of $`_3\chi N`$ in place of that of $`_3`$. The resulting equations are

$`0=\left(2U|z_j|^2\mu +{\displaystyle \frac{_T}{^2}}\right)z_j{\displaystyle \frac{_T}{^2}}Z,`$ (3)
where $`j=1,2,3`$, and $`Z:=z_1+z_2+z_3`$, $`\mu :=\chi +v`$ have been introduced . Therefore, any fixed-point, namely any vector $`(z_i)`$ $`(\eta _1,\eta _2,\eta _3)`$ that satisfies Eqs. (3), provides, at the same time, a dynamically active solution of Eqs. (2) represented by $`z_j(t)=\eta _j\mathrm{exp}[i\chi t/\mathrm{}]`$. Despite its time dependence, such a periodic orbit represents a (one-dimensional) extremal configuration of $`_3`$ on the hypersurface $`N=const`$. The solutions of Eqs. (3) can be grouped in various classes the first of which is represented by the case

$$Z=\mathrm{\Sigma }_jz_j=0.$$

The remaining cases are obtained from Eqs. (38)-(39), where $`z_i`$ can be replaced with the real quantities $`x_i`$ (see appendix A). Such cases are given by

$$\begin{array}{c}x_1=x_2=x_30,\\ \\ x_1=x_2x_3=\mathrm{\hspace{0.17em}0},\\ \\ x_1=x_2x_3\mathrm{\hspace{0.17em}0}.\end{array}$$

Such configurations are discussed below COM1 .

### III.1 Ground-state configurations

When $`x_1=x_2=x_3`$, Eqs. (38) and (39) are satisfied. Based on Eq. (37) combined with the conserved quantity $`N`$, one finds that $`x_j=\pm \sqrt{N/3}`$ which provides

$$z_i=\sqrt{N/3}\mathrm{exp}(i\mathrm{\Phi }),(i=1,2,3),$$
(4)
and $`3\mu =2UN3T`$. The energy of such configurations (they are shown to represent the ground-state in Sec. III) is given by

$$E_{gs}=\frac{_1}{^3}UN^2TN.$$
(5)
The ground-state phase $`\mathrm{\Phi }`$ is arbitrary since it represents a symmetry of the model. Also, the fact that $`z_j=|z_j|\mathrm{exp}(i\varphi _j)`$ have the same phase $`\varphi _j\mathrm{\Phi }`$ reproduces the symmetry breaking phenomenon that distinguishes the minimum energy state and, particularly, the vanishing of the phase difference (between closed points) in superfluid media. Such two features, which are known to characterize the superfluid ground-state of BH lattice model, naturally extend to the trimer model of condensates in our semiclassical picture. As regards Eqs. (2), they are solved by $`z_j(t)=\sqrt{N/3}\mathrm{exp}[i(\mathrm{\Phi }+\chi t/\mathrm{})]`$, with $`\chi =\mu v`$, that describes the ground-state collective mode.

### III.2 Vortexlike configurations

The situation in which $`Z=0`$ leads to a special configuration. In this case Eqs. (3) reduce to $`0=\left(2U|z_j|^2\mu +T/2\right)z_j`$ implying, in turn, that

$$|z_j|^2=\frac{_{2\mu T}}{^{4U}},j,\mu =\frac{_{4NU+3T}}{^6}$$

where the value of $`\mu `$ is derived from $`\mathrm{\Sigma }_i|z_i|^2=N`$. As a consequence of the independence of $`|z_j|`$ on the site index, the condition $`Z=0`$ is realized only if the phases of $`z_j=|z_j|\mathrm{exp}(i\theta _j)`$ are such that $`\theta _j(k)=2\pi jk/3+\mathrm{\Phi }_0`$ where $`k=1,2`$ and $`\mathrm{\Phi }_0`$ is an arbitrary phase. This configuration represents a particular case of the vortex state discussed for the Bose-Hubbard model on a $`M`$-well chain lattice ($`M>2`$) in Refs. CAPE  and FPZ1 . The energy associated to the vortex states

$$z_j(k)=\sqrt{N/3}\mathrm{exp}[i\theta _j(k)]$$
(6)
is given by

$$E_v=\frac{_{UN^2}}{^3}TN\mathrm{cos}\left(\frac{_{2\pi }}{^3}k\right)=\frac{_{UN^2}}{^3}+\frac{_{TN}}{^2},$$
(7)
while dynamics issued from Eq. (2) is described by the solution $`z_j(t)=z_j(k)\mathrm{exp}[it(\mu v)/\mathrm{}]`$.

### III.3 Configurations with a single depleted well

These configurations are characterized by the presence of a single depleted well (SDW). Without losing generality, one can choose the second well, so that one has $`x_1,x_30`$, and $`x_2=0`$. This case can be faced based on Eqs. (36): to satisfy $`E_2(𝐱)=0`$ one must impose $`x_3=x_1`$, which entails, in turn, that the equation $`E_1(𝐱)=0`$ is equivalent to the equation $`E_3(𝐱)=0`$. The latter, together with the constraint $`N=x_1^2+x_3^2`$, implies that $`x_1=\pm \sqrt{N/2}=x_3`$,$`x_2=0`$, which provides the fixed points

$$z_1=\sqrt{N/2}e^{i\mathrm{\Phi }},z_2=0,z_3=\sqrt{N/2}e^{i(\mathrm{\Phi }+\pi )},$$
(8)
with $`\mu 2Ux_1^2+T/2`$. Their energy is

$$E_{dw}=\frac{_1}{^2}UN^2+\frac{_1}{^2}TN.$$
(9)
Permutations of site indices $`j`$ of $`z_j`$ furnish other five fixed points of the same type. In view of state (8), it is worth noting that SDW configurations have the same structure of the $`\pi `$-states occurring in the dimer dynamics SM1 , where the phases of each well keep a constant phase difference $`\pi `$ in the course of time evolution. Solutions $`z_i(t)=z_i(0)\mathrm{exp}[it(\mu v)/\mathrm{}]`$ –with $`z_i(0)`$ given by Eqs. (8)– exhibit this property.

### III.4 Dimerlike configurations

For such states the variables $`x_j`$ (recall that $`x_j^2`$ is the boson number of the $`j`$-th well) satisfy the condition $`x_i=x_jx_k`$. Three cases are thus obtained through index permutation. In general, since three fixed points turn out to be associated with each dimerlike configuration (namely to each choice $`x_i=x_jx_k`$), nine fixed points are finally found in the dimeric class.

To deal with an explicit case, we shall consider the fixed points of the case $`x_1=x_2x_3`$. Owing to its complexity, their derivation is described in appendix B. Here we discuss the results obtained. Fixed points depend on the parameter $`\tau :=T/NU`$ and are representable as points $`(x_1,x_2,x_3)𝐑^3`$ on a sphere due to the constraint $`N=\mathrm{\Sigma }_ix_i^2`$. Their expressions read

$`\begin{array}{c}A_1:=(a_1,a_1,(a_1/|a_1|)\sqrt{N2a_1^2}),\\ \\ A_2:=(a_2,a_2,(a_2/|a_2|)\sqrt{N2a_2^2}),\\ \\ A_3:=(a_3,a_3,(a_3/|a_3|)\sqrt{N2a_3^2}),\end{array}`$ (15)
where

$$a_3=\pm \left[\frac{Np^2}{1+2p^2}\right]^{\frac{1}{2}},a_\nu =\pm \left[\frac{N}{2+q_\nu ^2}\right]^{\frac{1}{2}},$$
(16)
and $`\nu =1,2`$. Parameters $`q_\nu `$ and $`p`$ are defined as implict functions of parameter $`\tau `$ through systems (42) and (43), respectively. In particular, the cubic equation

$$\tau (2+q^2)(2+q)+4q(1+q)=0,$$
(17)
derived from Eq. (42), with $`q[1,0]`$ furnishes the $`\tau `$ functions $`q_1(\tau )`$, $`q_2(\tau )`$ corresponding to the real roots of Eq. (17). Notice that $`1q_1(\tau )<q_2(\tau )0`$ for $`\tau <\tau ^{}`$, whereas, for $`\tau >\tau ^{}`$, there are no solutions. The value $`\tau ^{}`$, where $`q_1(\tau ^{})=q_2(\tau ^{})`$ (and $`A_1A_2`$), is calculated in appendix B. Instead, system (43), with $`p[1/2,0]`$, always exhibits a single solution $`p(\tau )`$ that is carried out from the equation

$$\tau (1+2p^2)(1+2p)+4p(1+p)=0.$$
(18)
The dependence on physical parameters $`T`$, $`U`$, $`N`$ of $`A_1(q_1)`$, $`A_2(q_2)`$, $`A_3(p)`$ by means of parameter $`\tau `$ is thus established.

Points $`A_i`$ found in this way generate, by varying $`\tau `$, three curves on the sphere with $`N=const`$ \[in view of the sign $`\pm `$ in Eqs. (16) they actually are six\]. These become eighteen when considering the fixed points generated by index permutations. This process is described in appendix B, where the actual number of dimerlike fixed points is shown to reduce to twelve. Such curves \[parametrized by $`\tau `$ via $`p(\tau )`$ and $`q_\nu (\tau )`$\] can be proven to never intersect with one another except for the special case $`\tau =\tau ^{}`$ where $`A_1A_2`$. This coalescence effect is discussed below.

If the values of $`q_\nu `$ and $`p`$ for some given $`\tau \tau ^{}`$ are carried out explicitly, the energy for $`A_i`$

$$E_d=U[N^2+6a_i^44Na_i^2]T[a_i^22a_i\sqrt{N2a_i^2}],$$
(19)
($`i=1,2,3`$) is obtained via formulas (16).

We conclude illustrating the physical situations that correspond to configurations $`A_1`$, $`A_2`$ and $`A_3`$ when $`\tau `$ changes. Also, we compare them to the pure-dimer scenario SM1 ; FPZ1 . Let us start with $`\tau 0`$. One has $`q_1=1`$ and $`q_2=p=0`$ that entail

$`\begin{array}{c}A_1(1):=\pm (\sqrt{N/3},\sqrt{N/3},\sqrt{N/3}),\\ \\ A_2(0):=(\sqrt{N/2},\sqrt{N/2},\mathrm{\hspace{0.17em}0}),\\ \\ A_3(0):=(0,\mathrm{\hspace{0.17em}0},\sqrt{N}),\end{array}`$ (25)
respectively. By increasing $`\tau `$, $`A_1(q_1)`$ and $`A_2(q_2)`$ get closer and closer ($`A_1=A_2`$ for $`\tau \tau _{}`$). When $`\tau >\tau _{}`$ only the fixed point $`A_3(p)`$ survives. In particular, $`\tau \mathrm{}`$ implies $`p=1/2`$, so that \[from the third equation of Eqs. (15)\] one obtains

$$A_3(\frac{_1}{^2})=\pm (\sqrt{N/6},\sqrt{N/6},\sqrt{2N/3}).$$
(26)
Since $`A_3(p)`$ \[as well as $`B_3(p)`$, $`C_3(p)`$, obtained via index permutation (see appandix B)\] is shown to be a maximum in Sec. IV, indeed $`A_3(p)`$ appears to be comparable with the $`\tau `$-dependent maximum of the (pure) dimer model FPZ1 , where a unique well ends up to be filled when $`\tau 0`$. Instead, when $`\tau `$ is increased, no merging of $`A_3(p)`$ with other maxima \[e.g., $`B_3(p)`$, $`C_3(p)`$\] happens as that observed in the (pure) dimer model. In this model, in fact, a (macroscopic) coalescence effect takes place \[see Ref. SM1 ; FPZ1 \] since two maxima merge in a unique one when $`\tau >1`$ \[the opposite effect (bifurcation) occurs for $`\tau <1`$\]. As shown in Appendix A, such effects involving maxima pairs do not distinguish trimer dynamics.

A different macroscopic effect happens, however,in the trimer phase space. This is caused by the merging of distinct saddle points (e. g., $`A_1`$, $`A_2`$) discussed above, that disapper for $`\tau >\tau ^{}`$. Since chaotic behavior develops around saddle points, their coalescence should entail an evident local regularization of dynamical behavior.

As a final remark, we wish to observe that state $`A_1(1)`$ exhibits the same per-well boson distribution of the ground-state. Such two states differ uniquely owing to the phase of the third well. Similarly, state $`A_1(0)`$ and the $`\pi `$-like state $`(\sqrt{N/2},\sqrt{N/2},0)`$ just differ owing to the opposite phase of the second well. Despite their identical boson distribution, such situations will exhibit very different behavior thus confirming the profound influence of the interwell phase differences in distinguishing dynamical states.

## IV Character of fixed points

In this section we consider the stability character of fixed points just identified. Such a character is recognized by studying the second variation of the energy function on the hyper-surface defined by $`N=const`$. Explicitly, one should analyze the signature of the quadratic form associated to the Hessian of $`_3`$ (with $`N=const`$) in each fixed point. Being this process rather technical, we mainly develop it in Appendix D. Below, after showing how the separation of $`_3`$ in two independent (local) subhamiltonians simplifies remarkably the stability analysis, we summarize the results obtained.

To simplify $`_3`$, it is advantageous to introduce the local variables $`\xi _j=z_jv_j`$, where $`v_j`$ are the coordinates of some given fixed point. Then, neglecting third and fourth order terms, $`_3`$ takes the form

$$_3=_3(v)+\mathrm{\Sigma }_{i=1}^3\left(2U|v_j|^2\mu T\right)|\xi _i|^2$$

$$+\mathrm{\Sigma }_{j=1}^3U(v_j^{}\xi _j+v_j\xi _j^{})^2+\frac{_T}{^4}\mathrm{\Sigma }_{ik}|\xi _i\xi _k|^2,$$
(27)
which undergoes the further semplification $`v_i=`$ $`v_i^{}x_i`$ when one recalls that the phase factor of $`v_i𝐂`$ is a constant that can be absorbed by $`\xi _j`$. This fact allows us to separate $`_3`$ in a $`q`$-dependent part and a $`p`$-dependent part, with $`\xi _j:=q_j+ip_j`$. By making explicit the latter definition in $`_3`$, we find

$$_3=_0+h(q;6U,T,\mu )+h(p;2U,T,\mu ),$$
(28)
where $`_0:=_3(v)`$, and

$$h(q;6U,T,\mu ):=\underset{ij}{}(M_q)_{ij}q_iq_j,$$
(29)
$$h(p;2U,T,\mu ):=\underset{ij}{}(M_p)_{ij}p_ip_j.$$
(30)
Dynamical matrices $`M_q`$ and $`M_p`$ are defined as

$$M_q=\frac{T}{2}\left[\begin{array}{ccc}\mathrm{\Delta }_1& 1& 1\\ 1& \mathrm{\Delta }_2& 1\\ 1& 1& \mathrm{\Delta }_3\end{array}\right],$$
(31)
$$M_p=\frac{T}{2}\left[\begin{array}{ccc}\delta _1& 1& 1\\ 1& \delta _2& 1\\ 1& 1& \delta _3\end{array}\right],$$
(32)
with $`\mathrm{\Delta }_j:=2(\mu 6Ux_j^2)/T`$ and $`\delta _j:=2(\mu 2Ux_j^2)/T`$. Hence the diagonalization of the (local) quadratic form associated to $`_3`$ can be performed in a separate way on (29) and (30). Further simplifications come from the fact that matrix $`M_p`$ is proven to always have a vanishing eigenvalue (see Appendix C), while, due to the conserved quantity $`N=\mathrm{\Sigma }_i|z_i|^2=const`$, the induced local constraint $`q_1x_1+q_2x_2+q_3x_3=0`$ makes $`h(q;2U,T,\mu )`$ dependent only on two variables $`q_i`$ (see Appendix D).

The stability character for ground-states, vortex states, SDW states, and dimerlike states is studied explicitly in Appendix D. The calculation of the $`_3`$ second variation and, when this is necessary, of Hessian eigenvalues provide the following scenario.

\- States with $`x_1=x_2=x_3`$ are energy minima;

\- Vortex configurations are saddle points;

\- SDW configurations are saddle points;

\- Dimeric states exhibit two saddle points and one maximum point for $`0<\tau <\tau ^{}`$.

In the last case the merging of the two saddle points to form a regular point is enacted for $`\tau \tau ^{}`$. Thus, for $`\tau ^{}<\tau `$, a single maximum survives.

## V Mapping of dynamics on the reduced phase space

We develop a both qualitative and quantitative analysis of the chaotic behavior of the trimer based on the Poincaré section method. After performing a qualitative study of the periodic-orbit instability, we effect a quantitative analysis by measuring the very chaos’ indicator, namely the Lyapunov exponent of every single orbit choosen near a period one. To this end, we introduce a more suitable coordinate system involving a symplectic reduction of dynamics.

Hamiltonian (1) depends on three complex degrees of freedom and commutes with the boson number $`N`$. This fact as well as the structure of coupling term in $`_3`$ permit us to reduce from six to four the number of (real) canonical coordinates. The numerical integration of Eqs. (2) thus furnishes a system picture consisting of a trajectory in a four-dimensional ($`4d`$) reduced phase-space $`𝒫𝐂^3`$. In $`𝒫`$, a Poincaré section (PS) is the figure made by the points where a trajectory cut a $`2d`$ reference plane. The new set of canonical coordinates used to construct the PS embodies explicitly the conserved quantity $`\mathrm{\Sigma }_{i=1}^3|z_i|^2=N`$. Complex coordinates $`z_i=\sqrt{n_i}\mathrm{exp}(i\theta _i)`$ are replaced by

$$\begin{array}{cc}\varphi _1=\theta _2\theta _1,\hfill & \xi _1=(n_2+n_3n_1)/N\hfill \\ & \\ \varphi _2=\theta _3\theta _2,\hfill & \xi _2=(n_3n_1n_2)/N\hfill \\ & \\ \psi =(\theta _3+\theta _1)/2,\hfill & N=n_1+n_2+n_3,\hfill \end{array}$$
(33)
which obey the canonical Poisson brackets

$$\begin{array}{cc}\{\varphi _i,\xi _j\}=2\delta _{ij}/N,\hfill & \{\varphi _i,N\}=0,\hfill \\ & \\ \{\psi ,\xi _j\}=0,\hfill & \{\psi ,N\}=1.\hfill \end{array}$$
(34)
With such new variables COM2  $`_3`$ becomes

$$=\frac{_2}{^{UN^2}}(\xi _1,\xi _2,\varphi _1,\varphi _2)=$$

$$\xi _1^2+\xi _2^2\xi _1\xi _2\xi _1+\xi _2+_0$$

$$\tau \sqrt{(1\xi _1)(\xi _1\xi _2)(1+\xi _2)}\times $$

$$\left[\frac{\mathrm{cos}\varphi _1}{\sqrt{(1+\xi _2)}}+\frac{\mathrm{cos}\varphi _2}{\sqrt{(1\xi _1)}}+\frac{\mathrm{cos}\varphi _{12}}{\sqrt{(\xi _1\xi _2)}}\right]$$
(35)
in which $`\varphi _{12}:=\varphi _1+\varphi _2`$, $`_0=12\mu /(UN)`$ (the associated Hamiltonian equations are contained in Appendix E). In terms of coordinates $`(\xi _1,\xi _2,\varphi _1,\varphi _2)`$ the ground-state configuration (4) and vortexlike fixed points (6) correspond to

$$(\frac{_1}{^3},\frac{_1}{^3},\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0}),(\frac{_1}{^3},\frac{_1}{^3},\frac{_2}{^3}\pi k,\frac{_2}{^3}\pi k),$$

($`k=1,2`$), respectively. Dimeric fixed points \[consider, e. g., $`A_1`$ in Eq. (15)\] are given by

$$(\xi _1,2\xi _11,0,\pi ),$$

where $`\xi _1=12a_1^2(\tau )/N`$ and $`a_1(\tau )`$ is defined by Eqs. (16), whereas SDW states can be expressed as $`(0,0,\phi ,\pi \phi )`$ for $`n_2=0`$.

Operationally, the motion equations (2) are numerically integrated by using a first-order bilateral sympletic scheme; the algorithm precision is checked by monitoring the conserved quantities, that is, the system energy and total number of bosons. Trajectories can be traced in the phase space $`𝒫`$ in terms of $`\xi _1`$, $`\xi _2`$, $`\varphi _1`$, and $`\varphi _2`$. For any given value of the reduced energy, Hamiltonian (35) defines a 3d hypersurface in $`𝒫`$. The 2d surface used to construct a PS then is obtained by firming the value of $`\xi _2`$ to a constant. Hence, PS is made by coordinates $`(\xi _1,\varphi _1)`$ of the points’ set in which trajectories cut the selected 2d surface.

### V.1 Discussion of numerical results

We present the results of the numerical analysis aimed at investigating the phenomenology of the trimer dynamics in proximity of fixed points -these identify with the periodic orbits that stationarize $`_3`$ with $`N=const`$\- calculated in Sec. III. The Hamiltonian parameters and the total boson number chosen for the numerical simulations are

$$U=T=\mathrm{\hspace{0.17em}1},v=0,N=\mathrm{\hspace{0.17em}10},$$

respectively. Simulations have been carried out by using an integration time step of order $`1\times 10^4`$, while the total number of time steps employed in constructing each orbit is of order $`2^{28}`$. For each periodic orbit a PS has been selected by setting $`\xi _2=const`$ and $`d\xi _2/dt>0`$; on this section we have considered samples of about 100 initial conditions (IC). Also, the maximum Lyapunov exponent (MLE) has been measured for each trajectorie of the samples associated to the extremal periodic orbits.

Ground-state. The PS of this case is fixed by the condition $`\xi _2=1/3`$. About 100 initial conditions have been choosen in the $`\xi _1\varphi _1`$ plane, placed close to the point $`(\xi _1,\varphi _1)=(1/3,0)`$. The energy of the corresponding trajectories is $`E0.47(UN^2/2)`$. Fig. 1 shows that all the trajectories are regular in the phase space $`𝒫`$. Their regular character appears to be consistent with the periodic character shown by the oscillations of populations $`n_i=|z_i|^2`$ (see Fig. 2) of each condensate.

Vortexlike initial conditions. The choice $`\xi _2=0.3`$, $`E0.77(UN^2/2)`$ distinguishes the PS of the present case. The initial conditions for the trajectories have been chosen close to the fixed points $`(\xi _1,\varphi _1)=`$ $`(1/3,2k\pi /3)`$, $`k=1,2`$ (vortex state). Since no difference distinguishes the PS with $`k=1`$ and that with $`k=2`$, we restrict our attention to $`k=1`$. Fig. 3 shows the presence of both regular and chaotic trajectories \[Fig. 4 supplies two examples, one for each orbit type\]. It shows as well that the PS points related to chaotic trajectories are distributed in a region well separated from that occupied by points generating regular orbits. In particular, the PS displayed in Fig. 3 \[together with other PS involving slightly different $`\xi _2`$ ($`1/3`$)\] suggests that the vortex-state fixed point is basically surrounded by chaotic orbits. In Fig. 3, to reach the nearest regular orbits starting from $`\xi _1=1/3`$, $`\varphi _1=2\pi /3`$, a finite variation of both $`\xi _1`$ and $`\varphi _1`$ is necessary.

When this change is carried out the trimer-population oscillations change in a significant way. The time evolution of condensate populations $`n_i`$ ($`i=1,2,3`$) related to the nonchaotic orbit of Fig. 4 is illustrated in Fig. 5 and confirms its regular character. By considering IC closer and closer to the vortex state position such a character is progressively lost. This is shown in Fig. 6 that plots populations $`n_i`$, as a function of time, for the chaotic trajectory of Fig. 4.

The regular orbit of Fig. 4 involves an evident self-trapping effect provided their IC are enough far from $`(\xi _1,\varphi _1)=(1/3,2\pi /3)`$. This is clearly manifested in Fig. 5 where $`n_i(t)`$’s oscillate in such a way that $`n_2(t)<n_1(s)`$, $`n_3(s)`$, $`t`$, $`s`$: a stable gap, in fact, separates $`n_2`$ from $`n_1`$, $`n_3`$. On the contrary, no stable gap is involved, in general, by chaotic orbits (see Fig. 6) that develop large oscillations on the whole range of $`n_i`$. The fact that, in average, $`n_i1/3`$ is the only feature inherited by the vortex state.

SDW-like intial conditions. In this case the choice $`\xi _2=0.002`$ and $`E1.09(UN^2/2)`$ fixes the PS that is represented in Fig. 7. The zoom of the section reveals regular trajectories in $`𝒫`$ placed near SDW fixed points the latter being characterized by $`\xi _2=0=\xi _1`$, $`\varphi _1=\phi `$, $`\varphi _2=\pi \phi `$, where we have set ($`\phi `$ can be choosen arbitrarily) $`\phi =0.57\pi `$. Fig. 8 describes the PS of a chaotic orbit chosen among those of Fig. 7: The PS points are distributed in two, well separated, basins in a quite evident way. The interpretation of such an effect is the following: after recalling that setting $`\xi _20`$ implies that $`n_31/2`$, one deduces that, concerning the points of the PS, the values allowed for $`n_1`$ are either $`n_11/2`$ or $`n_10`$, which involves either $`n_20`$ or $`n_21/2`$, respectively.

One therefore recognizes the presence of an inversion population phenomenon beetwen $`n_1`$ and $`n_2`$. Interestingly, no intermediate values seems to be permitted. The corresponding scenario is given, on a shorter time interval, in Fig. 9 where the nonperiodic oscillations of populations $`n_1(t)`$, $`n_2(t)`$, and $`n_3(t)`$ are compared and a population inversion effect involving $`n_1`$, $`n_2`$ gets going. The population oscillations referred to a regular orbit of those contained in the zoom of Fig. 7 are shown in Fig. 10.

Dimerlike initial conditions. As proven in Secs. III and IV, the fixed points of this case consist of two saddle points and a maximum. The conditions $`\xi _2=0.295`$, $`E0.73(UN^2/2)`$ and $`\xi _2=0.005`$, $`E0.91(UN^2/2)`$ firm the PS associated to the first saddle (see Fig. 11) with coordinates $`\xi _1=\xi _2=0.295`$, $`\varphi _1=\varphi _2=\pi `$, and to the second saddle (see Fig. 15) with coordinates $`\xi _1=\xi _2=0.005`$, $`\varphi _1=\varphi _2=\pi `$, respectively. In both the cases, the PS’s exhibit both regular and chaotic trajectories. Concerning Fig. 11 (first saddle point), the coexistence of such regimes is confirmed in Fig. 12, where the PS of a regular trajectory is compared with the PS of a chaotic one. Figs. 16 shows analogous quantities referred to the second saddle point. It is worth noting that in the latter case the neighbourhood of the saddle point is characterized IC issuing regular orbits, whereas the first saddle is surrounded by IC generating chaotic motions.

In Fig. 15, regular orbits reside either on the bottom or on the top of the figure, but dynamics never connects the top orbits with the bottom ones. This feature is confirmed by the time behavior (see Fig. 17) of $`n_i(t)`$ related to the regular trajectory of Fig. 16. In Fig. 17, the presence of the gap beetwen $`n_2(0.5)`$ and $`n_1`$, $`n_3`$ $`(5)`$ indicates a macroscopic self-trapping effect. The latter differs from the dimeric self-trapping reviewed in Sec. II-A, where $`z_1=z_3n_1=n_3`$, in that $`n_1`$, $`n_3`$ develop independent oscillations. The scenario just described no longer holds for the chaotic orbits: Fig. 18 shows an intermittent effect of population inversion beetwen $`n_1`$ and $`n_2`$. This reflects the fact the points of the chaotic orbit of Fig. 16 are distributed intermittently both in the higher and in the lower part of the PS.

The trajectories near the maximum (we consider the case $`\xi _1=0.999`$, $`\xi _2=0.998`$, $`\varphi _1=\varphi _2=\pi `$,) appear to be regular, as one can deduce from Fig. 19 that shows PS close to a maxium fixed point with $`\xi _2=0.8`$ and $`E1.7(UN^2/2)`$. In this case, populations $`n_1(t)`$, $`n_2(t)`$, and $`n_3(t)`$ display in Fig. 20 a periodic effect of population-inversion involving $`n_1`$, $`n_2`$, whereas the fact that $`n_1`$, $`n_2`$ $`<<n_3`$ entails an evident self-trapping phenomenon.

The averages of the MLEs calculated for the chaotic orbits are 0.74 (Fig. 3), 0.42 (Fig. 7), 0.29 (Fig. 11), and 0.96 (Fig. 15). The evaluation of the MLE for the regular trajectories described in the previous examples (as well as for those of Fig. 1 and Fig. 19) exhibit the expected decreasing behavior thus suggesting that a weak stochasticity occurrs on such trajectories.

## VI Conclusions

In this paper we have focused our attention on the structure of the phase space of ’classic’ trimer, that is the mean-field form of the model describing three interacting BECs. Our analysis puts in light, on the one hand, the remarkable complexity that characterizes the trimer dynamics (by comparison with the integrable dynamics of the dimer system), on the other hand, the phenomena that are expected to characterize the $`M`$-well chain of interacting condensates. In view of the recent experimental results, the phenomenology of this system seems more and more viable to experimental observations. Trimer dynamics has been investigated within a semiclassical Hamiltonian picture, reviewed in Sec. II and developed in previous papers, based on a coherent-state representation of the trimer quantum state.

The identification of the set of fixed points of trimer Hamiltonian equations and the fact that such points are associated to periodic solutions (collective modes) of several types represents the initial, central result of our paper. The presence of the constraint $`N=const`$ entails that the states that stationarize the Hamiltonian are not isolated points but periodic orbits (one-dimensional manifolds). The solutions thus found enlarge the set of exact solutions NEM ; FP3  pertaining the dimerlike integrable subregimes of trimer and exhibiting a parameter-dependent self-trapping effect.

Based on the second variation (with $`N=const`$) of the energy function around its fixed points, the character of the latter has been recognized in Sec. IV and appendix D thereby revealing the presence of several saddle points and maxima, in addition to the expected ground-state. Numerical simulations and the PS method have furnished a wide scenario of trimer dynamical behaviors whose possible chaotic character has been tested by measuring the maximum Lyapunov exponent. We summarize the results of our dynamical simulations.

i) The orbits that have ICs close to ground-state exhibit a regular behavior with periodic oscillations of populations $`n_i`$.

ii) On the contrary, orbits with vortexlike IC (namely based at points close to vortex fixed points) are, in general, chaotic. Regular orbits, however, are found at sufficient distance from vortex fixed points. For such orbits a stable gap separates the oscillations of $`n_i`$ from those of $`n_j`$, $`n_k`$ ($`j,ki`$) thus generating self-trapping; $`n_j`$, $`n_k`$ show independent oscillations (in Sec. VI configurations with $`i=2`$, $`j=1`$, $`k=3`$ has been considered). The gap disappears for chaotic orbits.

iii) Orbits with SDW-like ICs also display both regular and chaotic behaviors but their ICs are not separated spatially. The regular orbits we have considered keep memory the IC since one of the three wells remains almost empty (pure dimer) while the other two undergo regular oscillations. Such states identify essentially with $`\pi `$-like states and manifest a stable character. An example of chaos emergence, which starts with a macroscopic population inversion entailing the filling of the (initially) depleted well, has been detected by assuming various SDW-like ICs involving a chaotic orbits.

iv) Regular orbits generated by dimerlike ICs (related to the second saddle point) exhibit periodic oscillations of $`n_i`$’s with an evident self-trapping. Chaotic orbits, instead, give rise to oscillations displaying self-trapping on short time intervals and intermittent population-inversion effects. ICs close to the maxima further generate regular orbits with self-trapping. The corresponding states display the presence of a unique almost filled well ff .

The scenario just depicted supplies a rich account of properties, behaviors and possible observable effects issuing from trimer dynamics and suggests promising future developments. An effect that might have a macroscopic character is related to the expected chaos suppression (see Sec. IIID) caused by the $`A_1`$-$`A_2`$ coalescence for $`\tau =\tau _{}`$. The basic configurations (ground-state, vortex states, SDW states, and dimerlike states involving both saddle and maximum points) recognized in the present work, and the complexity of dynamical regimes, both chaotic and regular, that develop in their neighbourhoods deserve further investigations in two directions at least. First, classically, one should carry out a systematic study (requiring huge computational resources) of long-time behavior of dynamical states of interest to disclose further macroscopic effects. For the same reason, a larger number of ICs (together with the trajectories thus issued) should be considered near fixed-point configurations.

We point out the fact that predictions on the dynamics of phases $`\varphi _i`$ should be important in relation to phase-interference experiments LATT . This aspect, which has not been deepened in the present work, requires a separate analysis and further numerical study directed to detect phenomena exhibiting phase coherence and their stability in proximity of states endowed with ordered phase configurations such as vortex states ($`\varphi _i=2\pi k/3`$, $`k=1,2`$), SDW states ($`\varphi _i\varphi _j=\pi `$, $`ij`$, $`|z_k|=0`$ with $`i,jk`$), and dimer configurations ($`\varphi _i=\varphi _j`$, $`z_kz_i=z_j`$ with $`i,jk`$).

Second, in view of the possibility of realizing systems with small per-well populations, the pure quantum approach to trimer dynamics (along the same lines of previous work directed to study the spectral properties of dimer) seems to be quite natural. The study of quantum trimer might put in evidence unexpected effects caused by the competition of chaotic (classical) behavior and integrable (quantum) behavior on the borderline of appropriate mesoscopic regimes where the transition from quantum to classical dynamics takes place.

###### Acknowledgements.

We are indebted to M. Pettini for useful observations and to P. Buonsante for his comments. The work of V.P. has been supported by I.N.F.M. (Italy) and M.U.R.S.T. (within the Project ”Theory of coherent fluids”). R.F. gratefully acknowledges the financial support of M.U.R.S.T. (within the Project ”Complex problems in Statiscal Mechanics and Field Theory”) and I.N.F.N. (Italy), as well as Prof. Arimondo and his group for useful discussions.

## Appendix A Derivation of fixed points

Eqs. (3) can be simplified by noting their invariance under the global symmetry transformation $`z_{\mathrm{}}z_{\mathrm{}}\mathrm{exp}(i\mathrm{\Phi })`$ and the fact that $`z_{\mathrm{}}/Z𝐑`$ (whenever $`Z0`$) with $`\mathrm{}=1,2,3`$. Then one can set $`z_{\mathrm{}}x_{\mathrm{}}\mathrm{exp}(i\mathrm{\Phi })`$, where the $`x_{\mathrm{}}`$ are real numbers and $`\mathrm{\Phi }`$ is an arbitrary phase, thus reducing Eqs. (3) to a system of the three real equations

$$E_j(𝐱):=\left(2Ux_j^2\mu +\frac{_T}{^2}\right)x_j\frac{_T}{^2}X0,$$
(36)
with $`j=1,2,3`$ and $`X:=x_1+x_2+x_3`$. When the condition $`x_i0`$ for $`i=1,2,3`$ is imposed, the latter equations can be recast in terms of an equivalent system of three equations one of which fixes the Lagrange multiplier $`\mu `$, while the other two, now formulated in a $`\mu `$-independent form, determine $`x_1`$, $`x_2`$, $`x_3`$ thanks to the further condition $`N=const`$. In fact, the sum of the quantities $`E_j(𝐱)/x_j`$ can be set equal to zero provided $`x_j0`$ thus giving the equation

$$3\mu =2UN\frac{_T}{^2}\mathrm{\Sigma }_i\frac{Xx_i}{x_i}.$$
(37)
Moreover, from $`E_1(𝐱)/x_1E_2(𝐱)/x_2=0`$ and $`E_3(𝐱)/x_3E_2(𝐱)/x_2=0`$ one obtains

$$0=(x_2x_1)\left[2U(x_2+x_1)+\frac{TX}{2x_1x_2}\right],$$
(38)
$$0=(x_2x_3)\left[2U(x_2+x_3)+\frac{TX}{2x_3x_2}\right],$$
(39)
completed by the condition $`N=x_1^2+x_2^2+x_3^2`$.

## Appendix B Derivation of fixed points of dimer type

Because of the identification $`x_1=x_2`$ (characterizing dimerlike fixed points), Eqs. (38) and (39) become a unique equation that can be written as

$$2U(x_1^2x_3^2)=\frac{T}{2}\left[\mathrm{\hspace{0.17em}1}+\frac{x_3}{x_1}2\frac{x_1}{x_3}\right].$$
(40)
Such an expression suggests two possible ways to parametrize $`x_1`$, $`x_3`$

$$\begin{array}{cc}x_1=\sqrt{N}R\mathrm{ch}\alpha ,\hfill & x_3=\sqrt{N}R\mathrm{sh}\alpha ,\hfill \\ & \\ x_1=\sqrt{N}R\mathrm{sh}\beta ,\hfill & x_3=\sqrt{N}R\mathrm{ch}\beta ,\hfill \end{array}$$
(41)
both allowing for the elimination of $`N`$ in the constraint $`2x_1^2+x_3^2N`$. They also provide two independent class of solutions that make explicit the three roots involved by the cubic character of Eq. (40). In the first case one finds

$$R^2=\frac{1q^2}{2+q^2},R^2=\frac{\tau }{4}\left[1+q\frac{2}{q}\right],$$
(42)
\[the first formula comes from the constraint on the total number of bosons, the second one comes from Eq. (40)\], whereas the second choice gives

$$R^2=\frac{1p^2}{1+2p^2},R^2=\frac{\tau }{4}\left[\mathrm{\hspace{0.17em}2}p1\frac{1}{p}\right],$$
(43)
where $`\tau :=T/UN`$, $`q=\mathrm{th}\alpha `$, $`p=\mathrm{th}\beta `$, and $`\alpha ,\beta 𝐑`$. System (42) reduces to the cubic equation

$$\tau (2+q^2)(2+q)+4q(1+q)=0,$$
(44)
which, provided $`q[1,+1]`$ in order to ensure the condition $`R^20`$, supplies either two or none solutions, depending on the fact that $`\tau <\tau _{}`$, $`\tau >\tau _{}`$. By solving the system one finds that the two roots $`q_\nu (\tau )`$, $`\nu =1,2`$ range in $`[1,0]`$ and fulfil the conditions

$$1q_1q_20,\mathrm{for}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}\tau \tau ^{},$$

with $`q_1(\tau )=q_2(\tau )`$ for $`\tau =\tau ^{}`$. The parameter $`\tau _{}`$ is identified by imposing, in addition to Eq. (44), the requirement that the two $`q`$-dependent functions of Eqs. (42) (that is, the two right-hand sides) are tangent at some point

$$\frac{6q}{(2+q^2)^2}\frac{\tau }{4}\left[1+\frac{2}{q^2}\right].$$
(45)
Equations (44) and (45), solved numerically, supply the value $`\tau ^{}0.29718`$.

Such a structure does not characterize system (43), which always exhibits a single solution for some appropriate value in the sector $`p[1/2,0]`$ obtained from the equation

$$\tau (1+2p^2)(1+2p)+4p(1+p)=0.$$
(46)
In view of the restriction $`q,p<0`$, from definitions (41) one deduces that the fixed-point coordinates are such that $`x_1x_3<0`$. The three solutions $`q_1(\tau )`$, $`q_2(\tau )`$, $`p(\tau )`$ just obtained correspond, within the space of coordinates $`\{(x_1,x_2,x_3)\}𝐑^3`$, to three vectors expressed as

$`\begin{array}{c}A_1:=(a_1,a_1,(a_1/|a_1|)\sqrt{N2a_1^2}),\\ \\ A_2:=(a_2,a_2,(a_2/|a_2|)\sqrt{N2a_2^2}),\\ \\ A_3:=(a_3,a_3,(a_3/|a_3|)\sqrt{N2a_3^2}),\end{array}`$ (52)
where

$$a_3=\pm \left[\frac{Np^2}{1+2p^2}\right]^{\frac{1}{2}},a_\nu =\pm \left[\frac{N}{2+q_\nu ^2}\right]^{\frac{1}{2}},$$
(53)
with $`q_\nu `$, $`p`$ solving Eqs. (44), (46). One can easily check that the two $`\tau `$-dependent curves $`A_1`$ and $`A_2`$ can be seen as two branches of a unique curve based at the common point $`A_1(q_1)=A_2(q_2)`$ for $`\tau =\tau ^{}`$, where they join smoothly.

Dimerlike fixed points become nine when considering the further set of points generated by index permutations

$$\begin{array}{c}B_1:=(a_1,(a_1/|a_1|)\sqrt{N2a_1^2},a_1),\\ \\ B_2:=(a_2,(a_2/|a_2|)\sqrt{N2a_2^2},a_2),\\ \\ B_3:=(a_3,(a_3/|a_3|)\sqrt{N2a_3^2},a_3),\end{array}$$
(54)
$$\begin{array}{c}C_1:=((a_1/|a_1|)\sqrt{N2a_1^2},a_1,a_1),\\ \\ C_2:=((a_2/|a_2|)\sqrt{N2a_2^2},a_2,a_2),\\ \\ C_3:=((a_3/|a_3|\sqrt{N2a_3^2},a_3,a_3),\end{array}$$
(55)
related to the subcases $`x_1=`$ $`x_3`$ $`x_2`$ and $`x_3=`$ $`x_2`$ $`x_1`$. Such nine curves become actually six. Our previous observation on considering $`A_1`$ and $`A_2`$ as a unique curve, in fact, readily extends the other curves $`B_1`$, $`B_2`$, and $`C_1`$, $`C_2`$. Owing to the double choice $`\pm `$ in formula (16) such curves are twelve.

In order to visualize the dimerlike fixed points, one can plot their position vectors $`A_{\mathrm{}},B_{\mathrm{}},C_{\mathrm{}}`$ $`\mathrm{}=1,2,3`$ within the three-dimensional space by varying their parameters $`a_{\mathrm{}}`$ ($`\mathrm{}=1,2,3`$) in the appropriate range. As noted above, in the case $`\tau <\tau ^{}`$, for $`q`$ ranging in the intervall $`[1,0]`$, Eq. (44) exhibits two solutions, whereas, for any $`\tau `$, Eq. (46) admits one solution with $`p[1/2,0]`$. The corresponding ranges of variation for $`a_\nu `$ $`(\nu =1,2)`$ and $`a_3`$ are

$$\sqrt{\frac{N}{3}}|a_\nu |\sqrt{\frac{N}{2}},0|a_3|\sqrt{\frac{N}{6}},$$

respectively. The representation of the position of dimeric fixed points on the sphere $`\mathrm{\Sigma }_ix_i^2=N`$ provides arcs that never intersect the one with the other when $`\tau `$ ranges in $`[0,\mathrm{}]`$.

The double sign $`\pm `$ in parametrizations (53) entails that each curve $`P_j`$ ($`P=A,B,C`$) is formed by two disjoint curves. For $`P_3`$ ($`P=A,B,C`$) one has

$$a_3[\sqrt{N/6},0[,a_3]0,\sqrt{N/6}],$$

while for $`P_1`$ and $`P_2`$ ($`P=A,B,C`$) the two disjoint branches are originated by the mappings $`a_\nu P_\nu `$ with

$$a_\nu [\sqrt{N/2},\sqrt{N/3}],a_\nu [\sqrt{N/3},\sqrt{N/2}].$$

## Appendix C

The diagonalization process of $`M_p`$ allows one to prove that one of the eigenvalues is always zero. ¿From the standard condition $`\mathrm{det}(M_p\lambda )=0`$ one obtains the eigenvalue equation (upon introducing $`\mathrm{\Lambda }:=2\lambda /T`$)

$$\mathrm{\Lambda }^3+(\mathrm{\Sigma }_j\delta _j)\mathrm{\Lambda }^2+(\delta _1\delta _2+\delta _2\delta _3+\delta _3\delta _13)\mathrm{\Lambda }+$$

$$+\delta _1\delta _2\delta _3+2\mathrm{\Sigma }_i\delta _i=0$$
(56)
where $`\delta _1\delta _2\delta _3+2\mathrm{\Sigma }_j\delta _j`$, upon setting $`\delta _j(x_{\mathrm{}}+x_k)/x_j`$ owing to Eqs. (36), can be shown to vanish in virtue of the identity

$$\mathrm{\Pi }_j^{}\frac{x_{\mathrm{}}+x_k}{x_j}=2+\mathrm{\Sigma }_j^{}\frac{x_{\mathrm{}}+x_k}{x_j}.$$

The superscript symbol $``$ recalls that the indices $`\mathrm{}`$, $`k`$, and $`j`$ must differ the one from the other. Hence, as a general result, the diagonalization of $`M_p`$ entails the presence of a zero eigenvalue consistent with the analysis of dynamics in the reduced phase space we develop in the sequel. In view of the invariance of the matrix trace, one also finds

$$\mathrm{\Sigma }_{j=1}^3\lambda _j=\frac{_T}{^2}\mathrm{\Sigma }_{j=1}^3\delta _j$$

$$\lambda _1+\lambda _2\frac{_T}{^2}\mathrm{\Sigma }_j\frac{x_{\mathrm{}}+x_k}{x_j}(3\mu +2UN)$$

while the two roots

$$\lambda _\nu =\frac{_T}{^4}\left[\mathrm{\Delta }\pm \sqrt{(\mathrm{\Delta })^24(\mathrm{\Delta }_{123}3)}\right]$$
(57)
where, $`\mathrm{\Delta }=\mathrm{\Sigma }_i\delta _i`$ and $`\mathrm{\Delta }_{123}:=\mathrm{\Sigma }_i\delta _1\delta _2\delta _3/\delta _i`$, and $`\nu =1,2`$, \[ $`\lambda _1`$ ($`\lambda _2`$) is associated with $``$ ($`+`$) in the formula\], are derived from the quadratic equation that emerges from Eq. (56) when removing a factor $`\mathrm{\Lambda }`$.

## Appendix D

This appendix is devoted to recognize, case by case, the minima, the maxima and the hyperbolic points within the four (class of) states identified as fixed points.

Concerning the eigenvalues of matrix $`M_q`$ one must take into account the restriction on the displacements $`\xi _i`$ from $`v_i`$ induced by the constraint $`N=\mathrm{\Sigma }_i|z_i|^2=const`$. After recalling that the phase of $`v_i`$ can be absorbed, for each $`i`$, by $`\xi _i`$ due to its arbitrariness, the substitution $`v_ix_i`$ implies that

$$N=\mathrm{\Sigma }_i|\xi _i+x_i|^2=\mathrm{\Sigma }_i\left[|\xi _i|^2+2q_ix_i+x_i^2\right],$$

which, in turn, entails $`\mathrm{\Sigma }_i|\xi _i|^2+2q_ix_i=0`$, namely –to first order– the plane equation

$$q_1x_1+q_2x_2+q_3x_3=0.$$

It represents the restriction on the displacements that variables $`\xi _i`$’s are allowed to effect. Substituting $`q_i`$ with $`q_i=(x_rq_r+x_sq_s)/x_i`$, where $`r,si`$ (and the choice of $`i`$ depends in general on the condition $`x_i0`$) finally gives

$$h(q;6U,T,\mu )$$

$$\left[2\frac{x_rx_s}{x_i^2}(6Ux_i^2\mu )+\frac{T(x_r+x_sx_i)}{x_i}\right]q_rq_s+$$

$$+\mathrm{\Sigma }_{ji}\left[12Ux_j^2\frac{\mu }{x_i^2}(x_i^2+x_j^2)+T\frac{x_j}{x_i}\right]q_j^2,$$
(58)
### D.1  Ground-state case

These fixed points are characterized by the fact that $`x_j=\pm \sqrt{N/3}`$, for $`j=1,2,3`$, and $`\mu =2UN/3T`$. By inserting this solutions in Eq. (28) one obtains the Hamiltonian written as

$`_3=E_{gs}`$ $`+`$ $`{\displaystyle \underset{j=1}{\overset{3}{}}}\left[({\displaystyle \frac{4}{3}}UN+T)q_j^2+Tp_j^2\right]+`$
$``$ $`{\displaystyle \frac{T}{2}}{\displaystyle \underset{ij=1}{\overset{3}{}}}\left(p_ip_j+q_iq_j\right),`$
where $`E_{gs}`$ is the ground state energy defined previously. Then, by taking the constraint $`N=\mathrm{\Sigma }_{j=1}^3|z_j|^2`$ into account, one obtains $`\mathrm{\Sigma }_{j=1}^3[q_j^2+p_j^2\pm 2\sqrt{N/3}q_j]=0`$. For little displacements from $`q_j=0=p_j`$ the latter equation reduces to $`\mathrm{\Sigma }_{j=1}^3q_j=0`$, which implies that

$`_3E_{gs}`$ $`+`$ $`({\displaystyle \frac{8}{3}}UN+3T)(q_1^2+q_2^2+q_1q_2)+`$ (59)
$`T{\displaystyle \underset{i=1}{\overset{3}{}}}p_i^2{\displaystyle \frac{T}{2}}{\displaystyle \underset{ij=1}{\overset{3}{}}}p_ip_j.`$
The eigenvalues of the Hessian corresponding to the $`q`$-dependent and the $`p`$-dependent part of $`_3`$ are $`\{(3\tau +8/3)UN/2,(9\tau +8)UN/2\}`$ and $`\{0,3T/2,3T/2\}`$, respectively. They are positive, coherent with the fact this is a minimum.

### D.2  Vortex case

The conditions $`|x_i|^2=N/3`$ and $`\mu =(4NU+3T)/6`$ characterize vortex configurations. One thus finds

$$h(p;2U,T,\mu )=T(p_1+p_2+p_3)^2/2<0,$$

whereas from Eq. (58) one gets

$$h(q;6U,T,\mu )=8UN(q_1^2+q_2^2+q_1q_2)/3>0$$

whose eigenvalues are always postive. Vortex configurations are therefore saddle points.

### D.3 SDW case

The previous analysis shows that the (fixed point) configurations in which one of the three well is depleted (e. g., well $`i=2`$) is such that $`x_2=0`$, $`x_1=x_3=\pm \sqrt{N/2}`$ ($`\pi `$-state structure), and $`\mu =NU+T/2`$. Site index permutations allows one to obtain two further, similar cases. In these points, Hamiltionan (28) can be written as

$`_3=E_{dw}`$ $`+`$ $`UN\left[2(q_1^2+q_3^2)q_2^2\right]+`$
$``$ $`{\displaystyle \frac{T}{2}}\left[({\displaystyle \underset{i=1}{\overset{3}{}}}q_i)^2+({\displaystyle \underset{i=1}{\overset{3}{}}}p_i)^2\right],`$
In this case, the constraint on the total number of particles supplies the constraint $`_{j=1}^3(p_j^2+q_j^2)\pm \sqrt{2N}(q_1q_3)=0`$, which reduces to $`q_1q_30`$ when $`q_j,p_j0`$. Then $`_3`$ takes the form

$`_3E_{dw}+UN(4q_1^2q_2^2)UNp_2^2+`$
$`{\displaystyle \frac{T}{2}}[(2q_1+q_2)^2+(p_1+p_2+p_3)^2],`$ (60)
whose Hessian matrix is endowed with the eigenvalues $`\{(65\tau \pm \sqrt{5(2012\tau +5\tau ^2)})UN/4\}`$ and $`\{0,(2+3\tau \pm \sqrt{44\tau +9\tau ^2})UN/4\}`$ for the $`q`$-dependent part and the $`p`$-dependent part, respectively. The analysis of the signature of such eigenvalues leads to identify the fixed points of the empty-well case with saddle points.

### D.4  Dimeric case

In the dimeric case the conditions on the coordinates are $`x:=x_1=x_2x_3=:y`$ (none vanishing). Furthermore, one has to impose the condition $`(2Ux_j^2\mu )x_j=T/2_{kj}x_k`$ on the Lagrange multiplier $`\mu `$ that becomes

$$\mu =(2UNT(1+y/x+x/y))/3.$$
(61)
Let us start by analyzing the Hessian eigenvalues of Hamiltonian $`p`$ sector namely of Eq. (30). These eigenvalues are given by Eq. (57) and one has that whenever the term $`\sigma :=(\delta _1\delta _2+\delta _2\delta _3+\delta _3\delta _13)`$ is greater then zero, these eingevalues are negative. To verify the condition $`\sigma >0`$ one can proceed in the following way. First, one substitutes in $`\sigma `$ the $`\delta _j`$ written in term of thier definitions; namely in term of $`\mu `$, $`T`$, $`U`$ and the fixed points coordinates $`x,y`$, Second, one writes $`\mu `$ in terms of $`x`$, $`y`$ and etc. \[Eq. (61)\]. Finally, one eliminates the dependence on $`x`$, $`y`$ and $`T`$ in $`\sigma `$ by choosing either the parametrization

$$y/x=q,x=\pm \sqrt{N/(2+q^2)}$$

combined with Eq. (44) to express the $`q`$-dependence of $`\tau `$, or the parametrization

$$x/y=p,x=\pm \sqrt{Np^2/(1+2p^2)}$$

combined with equation (46) to express the $`p`$-dependence of $`\tau `$. The expression for $`\sigma `$ achieved in such a way depends on $`q`$ or $`p`$, respectively. One can show that both the expressions are always positive in the range of definition of $`q`$ ($`[1,0]`$) and $`p`$ ($`[1/2,0]`$), which means that eigenvalues (57) are all negative.

As usual, for working out the Hessian eigenvalues of Hamiltonian (29) related to the $`q`$ part of the original one, it is necessary to take into account the constraint $`\mathrm{\Sigma }_i|z_i|^2=N=const`$. The latter, in the present case, becomes $`q_3(q_1+q_2)x/y`$. By means of this condition, one can reduce the dimension of the eigenvalues problem related to Hamiltonian (29) from 9 to 4. The Hamiltonian $`h_r(q_1,q_2,U,\mu ,T)`$ thereby obtained can be further simplified through the substitutions $`\mu \mu =[2UNT(1+y/x+2x/y)]/3`$ and $`y/x=q`$ or $`x/y=p`$, depending on the parametrization one adopts. With the first choice ($`y/x=q`$), and relying on Eq. (44), one finds two Hessian eigevalues one of which is always positive, whereas the other has an illdefined sign in the domain $`q[1,0]`$. By using the second parametrization $`x/y=p`$ and Eq. (46), both the eigenvalues thus obtained can be proven to be negative for $`p[1/2,0]`$. In summary, in the dimeric case, for $`0<\tau <\tau ^{}`$, one has two saddle points and one maximum point; for $`\tau ^{}<\tau `$, instead, fixed points reduce to a single maximum point.

## Appendix E

The fixed-point configurations corresponding to the change of Hamiltonian variables $`z_i`$, $`z_i^{}`$ $``$ $`\xi _a`$ $`\varphi _a`$, $`N`$, $`\psi `$ are obtained from the equations of motion rewritten in terms of the new variables (see below). Coordinates transformation (33) can exhibit (isolated) singular points in which they are not invertible. Dimer configurations in the presence of an empty well provides an explicit example where transformations (33) are illdefined. In fact, fixed points

$$(z_i)=(\sqrt{N/2}\mathrm{exp}i\varphi ,0,\sqrt{N/2}\mathrm{exp}i(\varphi +\pi ))$$

correspond to the set $`\{(0,0,\chi ,\pi \chi )|\chi \}`$, in the new description. As regards dynamical applications, fortunately, this problem is bypassed because the trajectories chosen in proximity of the periodic orbits associated to this kind of fixed points (as well as the PS used to study the dynamics near the same fixed points) do not contain, by construction, such pathological points. Upon setting $`\varphi _{12}:=\varphi _1+\varphi _2`$, $`s:=UNt`$, the equations of motion on the reduced phase space are given by

$$\frac{d\varphi _1}{ds}=12\xi _1+\xi _2+\frac{\tau }{2}[\sqrt{\frac{1+\xi _2}{\xi _1\xi _2}}\mathrm{cos}\varphi _2+$$

$$\frac{(12\xi _1+\xi _2)\mathrm{cos}\varphi _1}{\sqrt{(1\xi _1)(\xi _1\xi _2)}}\sqrt{\frac{1+\xi _2}{1\xi _1}}\mathrm{cos}\varphi _{12}]$$
(62)
$$\frac{d\varphi _2}{ds}=1\xi _1+2\xi _2+\frac{\tau }{2}[\sqrt{\frac{1\xi _1}{\xi _1\xi _2}}\mathrm{cos}\varphi _1+$$

$$\frac{(1\xi _1+2\xi _2)\mathrm{cos}\varphi _2}{\sqrt{(\xi _1\xi _2)(1+\xi _2)}}\sqrt{\frac{1\xi _1}{1+\xi _2}}\mathrm{cos}\varphi _{12}]$$
(63)
$$\frac{d\xi _1}{ds}=\tau \sqrt{(1\xi _1)(\xi _1\xi _2)}\mathrm{sin}\varphi _1+$$

$$\tau \sqrt{(1\xi _1)(1+\xi _2)}\mathrm{sin}\varphi _{12}$$
(64)
$$\frac{d\xi _2}{ds}=\tau \sqrt{(\xi _1\xi _2)(1+\xi _2)}\mathrm{sin}\varphi _2+$$

$$\tau \sqrt{(1\xi _1)(1+\xi _2)}\mathrm{sin}\varphi _{12}.$$
(65)