Bifurcation analysis of stationary solutions of two-dimensional coupled Gross-Pitaevskii equations using deflated continuation
E. G Charalampidis, N. Boullé, P. E. Farrell, P. G. Kevrekidis
aa r X i v : . [ n li n . PS ] N ov Bifurcation analysis of stationary solutions of two-dimensional coupledGross-Pitaevskii equations using deflated continuation
E. G. Charalampidis, ∗ N. Boull´e, † P. E. Farrell, ‡ and P. G. Kevrekidis
3, 2, § Mathematics Department, California Polytechnic State University, San Luis Obispo, CA 93407-0403, USA Mathematical Institute, University of Oxford, Oxford, UK Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, MA 01003-4515, USA (Dated: December 3, 2019)Recently, a novel bifurcation technique known as the deflated continuation method (DCM) wasapplied to the single-component nonlinear Schr¨odinger (NLS) equation with a parabolic trap in twospatial dimensions. The bifurcation analysis carried out by a subset of the present authors shedlight on the configuration space of solutions of this fundamental problem in the physics of ultracoldatoms. In the present work, we take this a step further by applying the DCM to two coupledNLS equations in order to elucidate the considerably more complex landscape of solutions of thissystem. Upon identifying branches of solutions, we construct the relevant bifurcation diagrams andperform spectral stability analysis to identify parametric regimes of stability and instability andto understand the mechanisms by which these branches emerge. The method reveals a remarkablewealth of solutions: these do not only include some of the well-known ones including, e.g., from theCartesian or polar small amplitude limits of the underlying linear problem but also a significantnumber of branches that arise through (typically pitchfork) bifurcations. In addition to presentinga “cartography” of the landscape of solutions, we comment on the challenging task of identifying all solutions of such a high-dimensional, nonlinear problem.
PACS numbers:
I. INTRODUCTION
The study of Bose-Einstein condensates (BECs) has offered a versatile playground for the examination of a diversehost of mesoscopic quantum phenomena for more than two decades [1, 2]. Among them, nonlinear wave dynamicalfeatures are of particular interest in the form of bright [3] and dark [4] solitons, vortices [5–7] and vortex lines andvortex rings [8]. As discussed, e.g., in the compendium of Kevrekidis et al. [9], many of these patterns spontaneouslyemerge in the nonlinear dynamics of BEC systems and subsequently play a critical role in their dynamical evolutionincluding their density and phase profiles, as well as their regular or chaotic/turbulent phenomena.Although the majority of studies have considered the one component (single species) case, multi-component BECsare also of considerable interest; see, e.g. [10] for a recent review. These may consist of mixtures of, for instance,different spin states of a particular atom (pseudo-spinor systems) [11, 12] or different Zeeman sub-levels of the samehyperfine level (so-called spinor condensates) [13–15]. Within these multi-component generalizations, various coherentstructures have been realized experimentally [16–21], with one of the most notable arguably being the dark-brightsolitons (and their dark-dark cousins). Experimental realizations have also been extended to spinor BECs [14, 15, 22–24], where recently also solitonic states have been observed [25].In parallel to experimental and theoretical studies, the development of numerical methods can also enhance ourunderstanding of single- and multi-component BECs. One such example that a subgroup of the present authorsrecently adapted to single-component atomic BECs [26] is the method of deflation and, more specifically, the deflatedcontinuation method (DCM) [27, 28]. Given one numerically exact solution of the system (computed, in our work,by means of Newton’s method), the aim of deflation is to construct a new problem where applying Newton’s methodwill no longer converge to the already obtained solution. Hence, if the solver is applied again, and upon successfulconvergence (within user-prescribed tolerances), it will discover an additional solution that was previously unknown.To set the stage, let F : U → V be a nonlinear map between Banach spaces whose roots are sought. Suppose that ∗ Email: [email protected] † Email: [email protected] ‡ Email: [email protected] § Email: [email protected] φ is an isolated root of F . Then, one can construct a new operator G : U → V such that: G ( φ ) . = (cid:18) k φ − φ k U + 1 (cid:19) F ( φ ) , (1)where k · k U is the norm on U . The essential idea is that k φ − φ k − U approaches infinity as φ → φ faster than F ( φ )approaches 0, hence avoiding the convergence to φ of the fixed-point iteration when applied to G . The addition ofunity ensures that the deflated problem G ( φ ) = 0 behaves like the original problem F far away from φ .In the present work, we apply the method of deflation and DCM to examine the possible solutions of the multi-component atomic BEC system. This will enable us to obtain significant insights into the pattern formation of thissystem. The present analysis obtains solutions that, to the best of our knowledge, were previously unknown, offeringa significant tool for understanding the landscape of nonlinear waveforms featured by the system. In addition, bycomputing the spectral stability of the resulting solutions one can identify not only which of these solutions arepotentially stable (and where this is parametrically so) but also bifurcations and further solutions arising from theDCM-identified solutions.Our presentation is structured as follows. In section II, we provide an overview of the existence and stabilityproblems for the two-component system. We also discuss an important twist on the deflation method to properlyaccount for symmetries of the problem (to avoid discovering multiple copies of the same solution that are e.g. relatedby rotation). In section III, we discuss our numerical results for the different branches of the system. Finally, insection IV, we summarize our findings and present some conclusions and future directions. The Appendix presentssome technical details regarding the eigenvalue computations. While the emphasis of our analysis will be on theexistence of the branches and their bifurcation diagrams, it is worthwhile to highlight that to perform the spectralanalysis, we will utilize the state-of-the-art capabilities of FEAST, an eigenvalue solver combining highly desirablefeatures of accuracy, efficiency and robustness for problems such as the one considered herein. II. THE MODEL AND SETUP
In this work, the model of interest is a two-component nonlinear Schr¨odinger (NLS) system in (2+1)-dimensions(two spatial and one temporal variable) given by i ∂ Φ − ∂t = − D − ∇ Φ − + (cid:0) g | Φ − | + g | Φ + | (cid:1) Φ − + V ( r )Φ − , (2a) i ∂ Φ + ∂t = − D + ∇ Φ + + (cid:0) g | Φ − | + g | Φ + | (cid:1) Φ + + V ( r )Φ + , (2b)where ∇ = ∂ /∂x + ∂ /∂y is the Laplacian operator, and D ± and g ij , i, j = { , } (with g = g ) are thedispersion and interaction coefficients, respectively. The function V ( r ) describes the external harmonic confinementand takes the form of V ( r ) = 12 Ω | r | , (3)with the parameter Ω capturing its strength and | r | = x + y . In the realm of atomic BECs, the functionsΦ ± ( r , t ) : D × R + ∪ { } 7→ C in Eqs. (2a) and (2b) represent the macroscopic wave functions with D ⊆ R being the(two-dimensional) spatial domain. For the reduction of the original three-dimensional (in space) BEC problem to thelower-dimensional setting of Eqs. (2a)-(2b), see, e.g., the discussion of [9].Stationary solutions to Eqs. (2a) and (2b) with chemical potentials µ ± are found by assuming the standing waveansatz Φ ± ( r , t ) = φ ± ( r ) e − iµ ± t , φ ± ( r ) : D C . (4)Upon inserting Eq. (4) into Eqs. (2a) and (2b) we arrive at a boundary value problem consisting of two coupled(elliptic) nonlinear partial differential equations F − (( φ − , φ + ) , µ − ) . = − D − ∇ φ − + (cid:0) g | φ − | + g | φ + | (cid:1) φ − + V ( r ) φ − − µ − φ − = 0 , (5a) F + (( φ − , φ + ) , µ + ) . = − D + ∇ φ + + (cid:0) g | φ − | + g | φ + | (cid:1) φ + + V ( r ) φ + − µ + φ + = 0 , (5b)together with zero Dirichlet boundary conditions, i.e., φ ± ( r ) | ∂D = 0. The latter conditions are rather inconsequential(meaning that zero Neumann or other boundary conditions would function equally well) because we will consider wideenough domains that the confining potential of Eq. (3) forces the density to tend to vanish well before we reach theedge of the computational domain. Furthermore F (( φ − , φ + ) , ( µ − , µ + )) := { F − (( φ − , φ + ) , µ − ) , F + (( φ − , φ + ) , µ + ) } (6)represents the set of equations considered in the deflated continuation method (DCM).We apply the DCM to Eqs. (5a) and (5b) in order to find steady-state solutions φ ± ( r ) for various values of thebifurcation parameter µ + . This algorithm for performing bifurcation analysis is based on Newton’s method and theappropriate choice of a deflation operator G to compute multiple solutions to a nonlinear partial differential equation.It should be noted that the work of [26] for the single-component NLS equation constructed deflated problems via G ( φ ) = (cid:18) k| φ | − | φ | k U + 1 (cid:19) F ( φ ) (7)in order to deflate the group orbit { e iθ φ | θ ∈ [0 , π ) } . That is, steady-state solutions to the (single-component) NLSequation are not isolated: if φ is solution then so is e iθ φ for any θ ∈ [0 , π ). The deflation operator given by Eq. (1)is not appropriate, as the solutions are not isolated, due to the Lie group of symmetries.In this work, we overcome this problem by further extending the deflation operator to eliminate rotations of solutionsto the NLS equation. That is, we wish to deflate the group orbit { φ θ : ( x, y ) φ ( x cos θ − y sin θ, x sin θ + y cos θ ) | θ ∈ [0 , π ) } . (8)To do so, we define the rotationally invariant transformation of a function u denoted by ˆ u asˆ u ( x, y ) . = 12 π Z π u ( x cos θ − y sin θ, x sin θ + y cos θ ) dθ. (9)Let ( φ − , φ ) be a solution to Eqs. (5a) and (5b), computed by Newton’s method. We construct the deflated problemfor finding the steady-state solutions to the two-component NLS system as follows: G (( φ − , φ + ) , ( µ − , µ + )) = (cid:13)(cid:13)(cid:13)(cid:16) [ | φ − | − \ | φ − | (cid:17) (cid:16) [ | φ + | − \ | φ | (cid:17)(cid:13)(cid:13)(cid:13) U + 1 F (( φ − , φ + ) , ( µ − , µ + )) , (10)where the norm for U is the L ( D ; R ) norm. Since the amplitude is invariant under phase shift and the transformationdefined by Eq. (9) is unchanged by rotations, the modified problem defined by Eq. (10) ensures nonconvergence to asolution in the group orbit of ( φ − , φ + ).Having identified the steady-state solutions φ ± ( r ), we perform a stability analysis by assuming the perturbationansatz e Φ − ( r , t ) = e − iµ − t h φ − + ε (cid:16) a ( r ) e iωt + b ∗ ( r ) e − iω ∗ t (cid:17)i , (11a) e Φ + ( r , t ) = e − iµ + t h φ + ε (cid:16) c ( r ) e iωt + d ∗ ( r ) e − iω ∗ t (cid:17)i , (11b)where ω ∈ C is the eigenfrequency, ε ≪ ∗ indicates complex conjugation. Uponinserting Eqs. (11) into Eqs. (2a) and (2a), we obtain at order O ( ε ) an eigenvalue problem written in matrix form as ρ abcd = A A A A − A ∗ − A − A ∗ − A ∗ A ∗ A A A − A ∗ − A − A ∗ − A abcd , (12)with eigenvalue ρ = − ω , eigenvector V = [ a b c d ] T , and matrix elements given by A = − D − ∇ + (cid:0) g | φ − | + g | φ | (cid:1) + V − µ − , (13a) A = g (cid:0) φ − (cid:1) , (13b) A = g φ − (cid:0) φ (cid:1) ∗ , (13c) A = g φ − φ , (13d) A = − D + ∇ + (cid:0) g | φ − | + 2 g | φ | (cid:1) + V − µ + , (13e) A = g (cid:0) φ (cid:1) . (13f)If we decompose the eigenfrequencies ω into their real and imaginary parts according to ω = ω r + i ω i , then thefollowing cases regarding the classification of steady-states φ ± in terms of their stability can be distinguished. If ω i = 0 holds for all ω , then the steady-state φ ± is classified as spectrally stable. For practical considerations, in ournumerical results presented in the next section, we assume that φ ± is stable if ω i < − holds for all eigenfrequencies;this case scenario is depicted by solid blue lines in the bifurcation diagrams presented therein. On the contrary, if ω i = 0, then two types of instabilities can be identified: • If ω r = 0 holds, then φ ± is classified as exponentially unstable and characterized by a pair of imaginary eigen-frequencies. • If ω r = 0, then φ ± is classified as oscillatorily unstable and characterized by a complex eigenfrequency quartet.These two scenarios are depicted by dashed-dotted red and green lines respectively in the bifurcation diagrams thatfollow to highlight the nature of the dominant unstable mode. If a change in the dominant instability type happens,e.g., from an exponentially unstable to an oscillatory unstable steady-state solution, that change will be highlighted bya change between the respective colors. The eigenvalue problem of Eq. (12) is very challenging to solve efficiently andaccurately; technical details of how to solve these eigenproblems, using, e.g., FEAST, are described in Appendix A.For our numerical computations presented below, we consider the spatial domain D = ( − , . Hereafter, wefix D − = D + ≡ . µ − = 1, g = 1 . g = 0 .
97 and g = 1, motivated by relevant values ofcoefficients previously used in the case of Rb [9]. The case of Rb is one where the two components are normallyimmiscible [1, 2, 9], that is to say it is energetically favored for them not to occupy the same location in space, a keydriving force in the pattern formation that we will observe below. Slight deviations of the interaction coefficients fromthe above values will not change the essence of our results, provided that one stays within this weakly immiscibleregime.Fixing the above parameters, we perform a natural parameter continuation in the value of µ + , up to µ + = 1 .
355 (or µ + = 1 . | φ + | = 1 g (cid:2) µ + − V ( r ) − g | φ − | (cid:3) , (14)at those points on the plane where µ + > V ( r ) − g | φ − | holds whereas the solution is zero elsewhere. However,as µ + becomes large enough during the continuation process, and due to the nonlinear coupling between the twocomponents, the “ − ” component gradually becomes smaller in its amplitude and eventually vanishes giving (atthis limit) a Thomas-Fermi solution exactly the same as the single-component case with a Thomas-Fermi radius of R TF ≈ .
23 (or ≈ . µ + wherea bound pair forms but with | φ + | ≪
1, the associated Thomas-Fermi solution and radius of the “ − ” component againcould be reduced to the single-component case with a radius in particular of R TF ≈ .
07 (note that µ − = 1 and is keptfixed). To put it otherwise, in the weakly immiscible regime considered herein, there exist two distinct Thomas-Fermiconfigurations, one dominated by the “+” component with the “ − ” component being absent and vice-versa. In eitherof these limits nevertheless, the coherent structures sit comfortably inside the chosen computational domain.Two spatial discretizations were used in this work. First, the DCM was applied to a finite element discretizationof Eqs. (5a) and (5b) using FEniCS [29] on a relatively coarse mesh (as this requires the solution of many nonlinearsystems with Newton’s method). The resulting solutions were then used as initial guesses for a second-order accuratecentered finite difference scheme on a uniform two-dimensional grid of equidistant points with lattice spacing ∆ x ≡ ∆ y = 0 .
08. The underlying linear system arising at each Newton step was solved by using the induced dimensionreduction IDR(s) algorithm [30] (see also [31] for its applicability on a similar computation). The solutions obtainedby the latter numerical scheme were subsequently used for calculating their spectra.In the bifurcation diagrams presented in the following section, we use the total number of atoms (or the sum of the(squared) L -norms of the respective fields) as our diagnostic defined via N t = N + + N − , N ± = Z D | φ ± | dxdy, (15)as well as the absolute total-number-of-atoms difference∆ N t = (cid:12)(cid:12)(cid:12) N (a) t − N (b) t (cid:12)(cid:12)(cid:12) , (16)between the total number of atoms of branches (a) and (b). III. NUMERICAL RESULTS
We begin the presentation of our results with Fig. 1. Specifically, each panel shows the densities (top row) of φ − (top left) and φ + (top middle) together with the associated phases (bottom left and middle) as well as the real(top right) and imaginary parts (bottom right) of the eigenfrequencies computed using the highly accurate FEASTeigenvalue solver (see Appendix A). The state in Fig. 1(a) corresponds to a bound mode consisting of a ground statein the “ − ” component and a soliton necklace in the “ + ” component (hereafter, we will refer to the “ − ” and “ + ”components as first and second components, respectively). The latter state was identified in the single-componentNLS equation (see [26] and references therein) and is generally unstable. It is relevant to note that this hexagonalstate bears alternating phases between its constituent “blobs” (as is customary generally in such dipolar, quadrupolaretc. states) and this hexagonal symmetry bears an imprint also on the spatial profile of the first component due tothe weak immiscibility between the components for the case of Rb, as explained above.According to the bifurcation diagram shown in the bottom left panel of Fig. 1, this bound state is stable from µ + ≈ .
223 (where the second component is formed, i.e., it is nontrivial) to µ + ≈ . µ + . In addition,and as per the eigenvalue continuations in Fig. 1(a) a real eigenvalue –or equivalently an imaginary eigenfrequency–(see the top right panel therein) passes through the origin at µ + ≈ . µ + ≈ . µ + ≈ [1 . , . µ + ≈ .
13 and µ + ≈ .
263 giving birth to a vortex-bright (VB) soliton dipole branchof Fig. 2(b) which has been studied in [33] (see also [34] for further analysis of this transverse instability) as well asthe tripole branch (of alternating vortices) of Fig. 2(c) in the second component (see also the bottom right panel ofthe figure). This progression of first the dipole (of VB solitons), then the tripole, then an aligned quadrupole etc.is strongly reminiscent of the corresponding process of breakup of a stripe into multi-vortex states due to transverseinstability in single-component BECs, as discussed, e.g., in [35]. The dipole branch itself appears dynamically stableexcept for a narrow interval of oscillatory instability, once again analogous to the one-component case [35]. On theother hand, the tripole branch is classified as exponentially unstable, inheriting the exponential instability of its parentDB stripe branch. While we are not aware of an experimental manifestation of states bearing multiple vortex-brightsolitary waves, the ability to tune the experimentally realized DB solitons [16–21] should, in principle, enable therealization of such observations, including, e.g., through the transverse instability of a dark-bright solitonic stripe. Itis relevant to also mention that a single VB solitary wave has been previously realized experimentally in the workof [36].A similar bifurcation pattern, i.e., a cascade of pitchfork bifurcations appears in Fig. 3. In particular, Fig. 3(a)corresponds to the crossed DB soliton [37], i.e., a fundamental state in the first component and the | , i (c) ∝ -10010 00.20.40.60.8 00.10.20.3 00.10.2-10 0 10-12-28 -0 -10 0 10 -0 1.23 1.27 1.31 1.3500.020.040.06-10010 00.20.40.60.8 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.31 1.33 1.3500.010.020.03 FIG. 1: A bifurcation diagram of a bright-soliton-necklace state. In particular, the clustered panels (a) present the | φ − | (topleft) and | φ + | (top middle), together with the respective phases (bottom left and middle panels therein). In addition, theassociated spectra are presented where the real (top right) and imaginary (bottom right) parts of the relevant eigenfrequencies ω are depicted (for the real part of these eigenfrequencies only the lowest ones are shown). Similarly, the clustered panels (b)correspond to the bifurcated state of a stripe and two vortex dipoles in the second component emerging at µ + ≈ . N t (left) and atom number difference ∆ N t (right) as functionsof µ + (see text for its definition). The vanishing of the latter is used to signal the bifurcation point. Recall that solid blue linesdenote stability, while red and green dash-dotted ones exponential and oscillatory instability, respectively, in the bifurcationdiagram here and in what follows. xye − Ω( x + y ) / quadrupolar state in the second component as per its Cartesian representation at the linear limit. See,e.g., [26] for this representation; here, we remind the reader for completeness that Cartesian eigenstates of the linearlimit can be denoted as | m, n i (c) . = φ m,n ∼ H m ( √ Ω x ) H n ( √ Ω y ) e − Ω r / , (17)where H m,n stands for the Hermite polynomials with m, n >
0, the quantum numbers of the harmonic oscillator. Theassociated energy of such a state (i.e., its eigenvalue) is E m,n . = ( m + n + 1)Ω [cf. [26]]. In the polar representation, -10010 00.20.40.6 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.08 1.16 1.24 1.3200.050.10.15-10010 00.050.1 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.16 1.24 1.3200.010.02-10010 00.20.4 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.28 1.3 1.3200.050.1 FIG. 2: Same as Fig. 1 but for the dark-bright (DB) soliton stripe branch. All densities and respective phases presented inpanels (a)-(c) are shown for values of µ + of µ + = 1 .
32. Note that bifurcations happen at µ + ≈ .
13 (b) and µ + ≈ .
263 (c),respectively. The emerging branches from the bifurcations of the DB stripe of the top panel are the VB dipole of the middlepanel (b) and the VB soliton tripole of the lower panel (c). -10010 00.20.40.60.8 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.16 1.24 1.3200.020.040.060.08-10010 00.10.2 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.2 1.26 1.3200.020.040.06 -10010 00.20.40.60.8 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.28 1.3 1.3200.020.040.06
FIG. 3: Same as Fig. 1 but for the crossed dark-bright (DB) soliton branch. Densities and phases shown in panels (a) and(b) correspond to a value of µ + of µ + = 1 .
32 whereas the ones presented in (c) correspond to µ + = 1 . µ + ≈ .
162 and branch (c) at µ + ≈ . -10010 00.20.4 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.16 1.24 1.3200.050.1-10010 00.20.4 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.23 1.26 1.29 1.3200.050.1-10010 00.20.4 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.26 1.28 1.3 1.3200.020.04-10010 00.10.20.3 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.26 1.28 1.3 1.3200.020.040.06 FIG. 4: Same as Fig. 1 but for the dark-bright (DB) ring soliton branch. Densities and phases shown in panels (a), (b) and(c) correspond to a value of µ + of µ + = 1 .
32 whereas the one presented in (d) corresponds to µ + = 1 . µ + ≈ . µ + ≈ . µ + ≈ . µ + ) are shown in Fig. 5. FIG. 5: Continuation of Fig. 4. The left and right panels present the total number of atoms and total number of atomsdifference both as functions of µ + . Notice the “spike” shown in the right panel as well as the fact that the ∆ N t of branches (b)and (c) from (a) become equal, i.e., equidistant from (a), at µ + ≈ .
273 corresponding to the intersection of their respectivecurves (see text for details). we also have | k, l i (p) . = φ k,l = q k,l ( r ) e ilθ (18)with eigenvalues E k,l . = (1 + | l | + 2 k )Ω. Here, l and k stand for the eigenvalue of the ( z -component of the) angu-lar momentum operator and the number of radial zeros of the corresponding eigenfunction q k,l , respectively. Theeigenfunction’s radial part can be denoted by q k,l ∼ r l L lk (Ω r ) e − Ω r / , (19)where L lk are the associated Laguerre polynomials. Interestingly once again, in the pattern of Fig. 3(a) note theastroid pattern of the first component induced by its phase immiscibility with the second component. The relevantstate bifurcates from the linear limit (i.e., in the absence of the second component) and is dynamically stable for µ + ≈ [1 . , . µ + ≈ [1 . , . µ + ≈ . µ + ≈ .
162 shown in Fig. 3(b). The latter waveform is oscillatorily unstable overthe interval of µ + we consider herein. A subsequent pitchfork bifurcation happens later at µ + ≈ .
279 (see also thebottom right panel of the figure) giving birth to the state of Fig. 3(c). This state is exponentially unstable (due tothe instability of its parent branch of Fig. 3(a)) and bears one vortex of charge 2 in the middle surrounded by anotherfour vortices in the second component.Subsequently, we focus on the DB ring soliton state of Fig. 4(a) and its bifurcations. This state (and in particular,with the first and second components reversed) has been identified and studied in [37] together with its stability.The state of Fig. 4(a) is generically unstable except for parametric intervals of stability µ + ≈ [1 . , . µ + ≈ [1 . , . µ + ≈ . µ + ≈ . µ + ≈ . µ + ≈ .
222 (a double imaginary eigenfrequency emanates from this collision) results in the emergenceof the two DB soliton stripes shown in Fig. 4(b). This daughter branch inherits the stability of the parent branchand maintains its stability for a very short parametric interval (see the bottom right inset of the left panel of Fig. 5showcasing the exchange of stability). Then, the two DB soliton stripes branch becomes oscillatorily unstable past µ + ≈ . µ + ≈ . µ + ≈ [1 . , . µ + of µ + ≈ . µ + ≈ .
247 (see also the upper left inset of the left panel of Fig. 5). Finally, the parent branch of Fig. 4(a) undergoesone more bifurcation at µ + ≈ .
251 giving birth to the VB hexagon state [37] (again, the emerging unstable modecorresponds to a double imaginary eigenfrequency). This branch is shown in Fig. 4(d) and classified as exponentiallyunstable although past µ + = 1 . µ + ≈ . -10010 00.10.20.3 00.51 00.10.2-10 0 10-12-28 -0 -10 0 10 -0 0.8 1 1.2 1.400.10.2-10010 00.10.2 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.82 1 1.18 1.3600.10.2-10010 00.050.1 00.51 00.10.2-10 0 10-12-28 -0 -10 0 10 -0 1.16 1.2 1.2400.050.1 -10010 00.010.02 00.51 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.34 1.36 1.3800.050.1 FIG. 6: Same as Fig. 1 but for the quadrupolar ring-dark soliton (RDS) branch. The densities and respective profiles areshown for µ + = 1 . µ + = 1 .
16 (b), µ + = 1 . µ + = 1 .
377 (d), respectively. The branches of (b), (c), and (d)emerge at values of µ + of µ + ≈ . µ + ≈ .
135 and µ + ≈ . FIG. 7: Continuation of Fig. 6. The left and right panels present the total number of atoms and total number of atomsdifference both as functions of µ + . Note that the first component of the branch (b) vanishes past µ + ≈ .
373 as well as (c)does past µ + ≈ .
244 and (d) past µ + = 1 . happens is offered in the following. Based on the definition of ∆ N t [cf. Eq. (16)], we measure the distance of N t [cf.Eq. (15)] of a bifurcating branch from a reference branch. Indeed, although the N t for branch (b) in Fig. 4 is smallerthan the one for the reference branch (a) over the entire interval in µ + considered therein, this difference appears positive in the right panel of Fig. 5, based on the definition of the relevant diagnostic. Furthermore, the “spike” ofthe curve therein corresponding to the difference of N t between branches (c) and (a) happens because the N t of theformer branch (as this emerges from branch (b)) is smaller than the N t of branch (a) (and larger than the N t ofbranch (b)) until its curve (as a function of µ + ), intersects the corresponding one of branch (a) at µ + ≈ . N t of branch (c) is larger than the one of branch (a).Hereafter, the configurations become more complex. The quadrupolar type solution in the first component describedas | , i (c) in its Cartesian classification trapping a ring-dark soliton (RDS) in the second component is shown inFig. 6(a). This branch emerges, i.e., bearing a non-trivial solution in its second component at µ + ≈ .
791 and isclassified as exponentially unstable over the parametric interval of µ + ≈ [0 . , . µ + ≈ .
815 giving birth to the daughter branch of Fig. 6(b). This branchbears a two-dark soliton stripes waveform (or | , i (c) as per its Cartesian classification) in the second componentwhose total number of atoms is less than its parent branch (see the left panel of Fig. 7). In addition, this state isclassified as exponentially unstable throughout its interval of existence where the first component vanishes eventuallyat µ + ≈ . µ + ≈ .
135 where thebranch of Fig. 6(c) bearing a hexapolar mode in the second component [26] is born. This branch, i.e., the quadrupolar-vortex-hexagon branch is classified as exponentially unstable over µ + ≈ [1 . , . µ + ≈ . µ + increases, the parent branch of Fig. 6(a) undergoes onemore pitchfork bifurcation giving birth to the quadrupolar vortex-octagons state of Fig. 6(d). Such a pattern, i.e.,the subsequent bifurcations of the RDS state was briefly mentioned in the one-component case [26]. In our case, thebranch of Fig. 6(d) has a narrow interval of existence of µ + ≈ [1 . , . µ + considered herein (see the left panel ofFig. 9). However, there is again a cascade of pitchfork bifurcations, the first of which occurs at µ + ≈ .
16 as is alsohighlighted in the right panel of Fig. 9. In particular, the first bifurcating branch is shown in Fig. 8(b) where thesolution in the second component can be thought of as a symmetry-broken (between the axes) variant of panel (a)that can be represented as a | , i (c) + | , i (c) Cartesian state at the linear limit. In the present two-component case,the branch of Fig. 8(b) is exponentially unstable as well (see the left panel of Fig. 9). The next pitchfork bifurcationoccurs at µ + ≈ .
242 where Fig. 8(a) gives birth to the bound mode shown in Fig. 8(c) with the second componentfeaturing an intriguing combination of a ring and a vortex necklace. Although this solution has been identified asoscillatorily unstable in [26] in the single-component NLS, the daughter branch shown in Fig. 8(c) is classified tobe exponentially unstable (see the left panel of Fig. 9). It is interesting to highlight here that although in somecases (such as in Fig. 1(a) as we saw before) the presence of a second component plays a stabilizing role, in otherssuch as the one of Fig. 8(c), it adds further unstable eigendirections to a particular waveform. Finally, the parent3 -10010 00.20.4 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.14 1.2 1.26 1.3200.050.1 -10010 00.20.4 00.20.40.6 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.16 1.24 1.3200.030.060.09 -10010 00.10.20.3 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.26 1.29 1.3200.040.08-10010 00.20.4 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.28 1.3 1.3200.030.060.09
FIG. 8: Same as Fig. 1 but for the quadrupolar–two dark rings branch. The branches (a)-(d) are presented for a value of µ + of µ + = 1 .
32. The branch of panel (a) undergoes a cascade of pitchfork bifurcations occurring at µ + ≈ . µ + ≈ .
242 and µ + ≈ . FIG. 9: Continuation of Fig. 8. The left and right panels present the total number of atoms and total number of atomsdifference both as functions of µ + . According to the left panel, the N t of branch (b) is less than the one of its parent branch(a). Similarly to Fig. 5 the ∆ N t of branches (b) and (c) from (a) become equal, i.e., equidistant from (a), at µ + ≈ . branch of Fig. 8(a) undergoes one further pitchfork bifurcation at µ + ≈ .
277 giving birth to the daughter branch ofFig. 8(d). This branch forms four vortices in its inner ring and is generally exponentially unstable (over the intervalwe considered herein) except for a narrow interval of µ + ≈ [1 . , . N t of the bifurcating branchof Fig. 8(b) is less than the one of its parent branch [cf. Fig. 8(a)], the respective total-number-of-atoms differencebetween those two branches is positive as is shown in the right panel of Fig. 9 per the definition of Eq. (16).The state of Fig. 10(a) involves a symmetry-broken state where the first component features an elliptical type of aRDS, while the second component represents a dipolar ( | , i (c) in the notation of the Cartesian Hermite eigenstates)structure bearing a π phase shift between the two matter wave blobs. This branch is highly unstable and in fact,classified as exponentially unstable over µ + ≈ [0 . , . µ + ≈ . µ + where the first componentvanishes, this branch undergoes a similar symmetry-breaking bifurcation as in the previous cases in this work givingbirth to the branch of Fig. 10(b) (see also the bottom right panel of Fig. 10). This daughter branch involves a seriesof density blobs in the first component centered along the nodal line of the second component. The latter features avortex quadrupole (this is hard to discern at the level of the density but more transparent at the level of the phaseof this component) which is exponentially unstable as is also shown in the inset in the bottom left panel of Fig. 10;its interval of existence over µ + is rather narrow (the first component is non-vanishing for µ + ≈ [1 . , . µ + ≈ . , . µ + ≈ .
99 where the first two will correspondto reverse pitchfork bifurcations (see, Figs. 12(c) and 13(g) next for an example of such a bifurcation).On the other hand, and as per the branch of Fig. 11(a), the first component features similarly a RDS whereas thesecond component involves a hexapolar double soliton necklace. The latter state (which also deforms the RDS into apattern with hexagonal symmetry due to the nonlinear interactions) was identified in the previous work of [26] in thesingle-component case and classified as exponentially unstable. In the present two-component case, the bound modeof Fig. 11(a) is exponentially unstable over the reported interval of µ + ≈ [1 . , .
4] (see the bottom left panel ofFig. 11). The double solitonic necklace emerges at µ + ≈ . µ + ≈ .
39 giving birth (via a pitchfork bifurcation) to the branch of Fig. 11(b) (see the bottom right panelof Fig. 11). This branch, in turn, features a vortex necklace (notice the modification in the relevant phase profile fromthe purely real parent branch of Fig. 11(a)) and is classified as exponentially unstable over the reported interval of µ + ≈ [1 . , .
4] as is shown in the bottom left panel of Fig. 11.We turn now our focus on the branches presented in Figs. 12, 13 and 14. Specifically, the branch of Fig. 12(a)involves a RDS-necklace state whose first component is a RDS of octagonal type and bears small “blobs” in its density.These blobs are complementary (due to inter-component repulsion) to the ones of the solitonic necklace of the secondcomponent. This bound mode emerges from the linear limit at µ + ≈ .
177 and is classified (according to our spectralstability analysis) as exponentially unstable over µ + = [1 . , .
4] (see also the top panel of Fig. 14). In addition, thisbranch gives birth to the structures of Fig. 12(b) and 12(c), respectively through pitchfork bifurcations. In particular,the state of Fig. 12(b) emerges at µ + ≈ .
19 and features a vortex necklace in the second component (see also thebottom panel of Fig. 14). It can be discerned from the second component of Fig. 12(b) that the blobs of the necklaceof its parent branch [cf. Fig 12(a)] re-arrange themselves due to the emergence of vorticity therein. Based on ourspectral stability analysis, this state is classified as exponentially unstable over µ + ≈ [1 . , . µ + ≈ .
344 (see the top panel of Fig. 14). Furthermore, the soliton necklace in the second component5 -10010 00.20.40.60.8 00.10.2 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.61 0.84 1.07 1.300.10.2
FIG. 10: Same as Fig. 1 but for the ring-dark-soliton-dipolar (RDS-dipolar) branch. The densities and associated phases areshown for values of µ + of µ + = 0 . µ + = 1 .
29 (b), respectively in the first and second rows. Note that the branch ofpanel (b) emerges at µ + ≈ . µ + . of the branch of Fig. 12(a) undergoes a density re-arrangement in its blobs at µ + ≈ .
192 resulting in the bound stateof Fig. 12(c) whose first component is morphed into an elliptical RDS state and is classified as exponentially unstableover µ + ≈ [1 . , . µ + ≈ . µ + ≈ [1 . , . µ + ≈ .
322 are not shown) as is depicted in the top panel of Fig. 14. Subsequently, the branch of Fig. 12(a)undergoes two further pitchfork bifurcations at µ + ≈ .
27 and µ + ≈ .
37, respectively. The first bifurcating state isshown in Fig. 13(e) which involves a vortex necklace of octagonal shape in the second component. The latter statewas identified in the single-component NLS equation in [26] (and references therein) as a combination of a doublering configuration and a necklace. According to our stability analysis results shown in the top panel of Fig. 14, thisbranch is exponentially unstable for µ + ≈ [1 . , .
32] and oscillatorily unstable for µ + ≈ [1 . , . µ + ≈ . µ + ≈ [1 . , .
4] (results6 -10010 00.20.40.60.8 00.10.2 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.34 1.36 1.38 1.400.050.1-10010 00.20.40.60.8 00.10.2 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.39 1.395 1.400.050.1
FIG. 11: Same as Fig. 1 but for the hexagonal type RDS-necklace branch. Similarly, the densities and associated phases areshown for a value of µ + of µ + = 1 . µ + ≈ .
39 via a pitchfork bifurcation. Again,the bottom left and right panels present the total number of atoms and total-number-of-atoms difference both as functions of µ + . are not shown past µ + = 1 . µ + ≈ .
362 merges with the branch ofFig. 13(g) (see also the black open circles in the bifurcation diagram in the top panel of Fig. 14). This mergingof 12(c) with 13(g) is a canonical example of a reverse pitchfork bifurcation where the branch of Fig. 13(g) is theparent branch. This branch is classified as exponentially unstable. However, as µ + increases and at µ + ≈ . µ + ≈ [1 . , .
4] over which it appears (dominantlyto be) oscillatorily unstable as is shown in the top panel of Fig. 14 (this branch is not shown for values of µ + past µ + ≈ . -10010 00.20.40.60.8 00.10.2 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.2 1.3 1.400.050.1-10010 00.20.40.6 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.2 1.3 1.400.050.1-10010 00.20.40.60.8 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.2 1.25 1.3 1.3500.050.1-10010 00.20.4 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.3 1.31 1.3200.030.060.09 FIG. 12: Same as Fig. 1 but for the RDS-necklace branch. All densities and phases are shown for a value of µ + = 1 . µ + ≈ .
32. The branch of panel (a) emerges at µ + ≈ . µ + ≈ .
19 and µ + ≈ . µ + ≈ .
362 (see textfor details). Finally, the branch of panel (d) emerges at µ + ≈ .
291 from (c). -10010 00.10.20.3 00.20.40.6 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.29 1.32 1.35 1.3800.050.1-10010 00.20.40.6 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.37 1.38 1.39 1.400.050.1-10010 00.20.40.60.8 00.20.40.6 00.10.2-10 0 10-12-28 -0 -10 0 10 -0 1.2 1.3 1.400.050.1-10010 00.10.2 00.20.40.6 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.3 1.35 1.400.050.1 FIG. 13: Continuation of Fig. 12. Densities and phases are presented for values of µ + of µ + = 1 .
32 (e), µ + = 1 . µ + = 1 . µ + = 1 . µ + ≈ .
27 and vanishes at µ + ≈ .
381 whereas thebranch of (f) emerges at µ + ≈ .
37 and is shown up to µ + = 1 .
4. Finally, the branch of panel (g) emerges at µ + ≈ . µ + ≈ . FIG. 14: Continuation of Fig. 12 and also 13. The total number of atoms and total number of atoms difference are shown inthe top and bottom panels, respectively as functions of µ + . The black open circles in the top panel are explained in the text. in reverse i.e., with the second component playing the role of the first one and vice versa. Recall that in the case ofthe equal g ij (the so-called Manakov model integrable in one-dimension [38]), these states would be interchangeable,however the weak asymmetry (and corresponding immiscibility) distinguishes between these states and the earlierones. Indeed, Fig. 15(a) provides a reversed analogue of Fig. 1(a) whereas the branch of Fig. 15(b) can be similarlyconnected to Fig. 3(a). These branches involve multipoles (a hexapole and a quadrupole) coupled to a fundamentalstate. The branch of Fig. 15(c) presents a DB solitonic stripe, analogous to that of Fig. 2(a). Also, the branch ofFig. 15(d) is tantamount to Fig. 4(a), namely the ring dark-bright soliton, whereas the branch of Fig. 16(a) is stronglyreminiscent of the branch of Fig. 6(a) involving a RDS and the | , i (c) (in its Cartesian classification) state. In viewof these similarities and the slight asymmetries of the g ij , the precise instability details of the states of Fig. 15(a)-(d)and Fig. 16(a) are somewhat different, yet qualitatively similar to the examples that we have already studied above.It should be noted that we performed the continuation of all the above branches over µ + and stopped when one ofthe components was found to be below numerical precision, i.e., when effectively the states become single componentwaveforms with the other component being trivial. Indicatively, the intervals of existence (and stability) of the statesshown in Fig. 15 are µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ≈ [0 . , . µ + ).As we advance to panels (b)-(d) of Fig. 16, we encounter higher excited state combinations. For instance, inFig. 16(b), a double RDS in the second component is coupled to a single RDS in the first component (this boundstate emerges at µ + ≈ . l = 0 azimuthal index) of the associated Laguerre eigenmodes of the two-dimensional system. Inparticular, these are modes with k = 2 in the second component and k = 1 in the first component. Subsequently, inFig. 16(c), a fundamental state in the first component is coupled with the multipole | , i (c) + | , i (c) in the notationof the Cartesian Hermite eigenstates resulting in a bound state emerging at µ + ≈ .
283 as a stable branch until µ + ≈ . µ + ≈ .
286 and features a higher order example of a Cartesianexcited state. Here, a | , i (c) state in the second component is coupled to a | , i (c) one in the first component.This bound state features a multiplicity of exponentially unstable modes, in addition to a number of oscillatoryinstabilities in suitable parametric ranges. Finally, it should be noted that all the above states with the exceptionof that of Fig. 16(c) are very highly unstable, as can be seen from the right panels of Fig. 16. Yet, it is relevant tocomment that some of these highly excited states like the one depicted in Fig. 16(c) may feature intervals of spectralstability and the highly accurate numerical techniques used herein can identify these states and their correspondingstability intervals.0 -10010 00.10.2 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.3 0.4 0.5 0.600.020.040.06-10010 00.20.4 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.4 0.6 0.800.050.1-10010 00.20.40.6 00.10.20.3 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.55 0.67 0.79 0.9100.050.10.15-10010 00.20.40.6 00.10.2 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.54 0.61 0.68 0.7500.050.1 FIG. 15: Branches of solutions that are reversed versions of ones discussed previously in this work. Densities and phases as wellas associated spectra are shown for µ + = 0 .
52 (a) (reversed version of Fig. 1(a)), µ + = 0 . µ + = 0 .
68 (c) (reversed version of Fig. 2(a)), and µ + = 0 .
65 (d) (reversed version of Fig. 4(a)), respectively. -10010 00.20.40.60.8 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.78 1.08 1.3800.10.2 -10010 00.20.4 00.20.40.6 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.27 1.31 1.3500.050.1 -10010 00.20.40.60.8 00.10.20.3 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.29 1.31 1.33 1.3500.0050.010.015 -10010 00.20.40.6 00.050.1 00.10.2-10 0 10-10010 -0 -10 0 10 -0 1.29 1.31 1.33 1.3500.050.1 FIG. 16: Same as Fig. 15. Note that panel (a) depicted for µ + = 0 .
95 is a reversed version of Fig. 6(a). The branches of panels(b)-(d) are all presented for a value of µ + of µ + = 1 .
32 at which the first component does not become the trivial solution. -10010 00.20.4 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.78 0.96 1.14 1.3200.050.10.15-10010 00.050.10.15 00.20.40.60.8 00.10.2-10 0 10-10010 -0 -10 0 10 -0 0.86 0.95 1.04 1.1300.030.060.09 FIG. 17: Same as Fig. 1 but for the rotated dipolar branch. Densities and phases are shown in panels (a) and (b) for values of µ + of µ + = 1 and µ + = 1 .
1, respectively. The bound state involving the vortex quadrupole in the second component of panel(b) emerges at µ + ≈ . A. Further linear modes and DCM
The DCM was initialized at µ + = 0 . µ + and µ − to the nonlinear regime, i.e., for non-vanishing values of N + and N − .3To this end, we briefly discuss a few examples (among several others) shown in Figs. 17 and 18 that could be usedin addition to the results of the DCM method. We note in passing that we have found numerous additional branchesto the ones reported here (stemming from the linear limit). However, we refrain from discussing all of them here toavoid cluttering the manuscript with additional figures.We begin our discussion with the quadrupolar complex of Fig. 17(a). In particular, one can envision the | , i (c) state in the one component co-existing with its rotated version in the second component (due to the mutual repulsionbetween components and their immiscibility) as is shown in the figure. This bound state emerges at µ + ≈ .
76 andis classified as exponentially unstable over µ + ≈ [0 . , .
32] (we did not perform the continuation past µ + = 1 . µ + ≈ . µ + ≈ [0 . , .
89] and then oscillatorilyunstable until µ + ≈ .
172 where the first component vanishes, thus reducing to the single-component case. The singlecomponent vortex quadrupole features an interval of oscillatory instability but is otherwise dynamically stable; seee.g., [26] and references therein.To offer an anthology of possible states that we have been able to identify on the basis of continuation from thelinear limit, we refer the reader to Fig. 18. These (together with the branches shown in Fig. 17) constitute examplesof nonlinear eigenstates of the problem that we have continued from their respective linear limit up to µ + = 1 . not obtained via the DCM.Although typical examples of the relevant solutions are shown, we do not present their associated spectra: wesimply present them as representative examples of the richness of additional available nonlinear patterns that aphysical understanding of relevant limits can enable us to access. Nevertheless, we briefly comment on intervals ofstability (when applicable) and instability of these waveforms. Similarly as in Fig. 17(a), one can construct a dipolestate (in the form of | , i (c) in the first component) featuring its rotated version in the second component as is shownin Fig. 18(a). This state emerges at µ + ≈ .
681 and is generally exponentially unstable except for the interval of µ + ≈ [0 . , .
05] in which it appears to be oscillatorilly unstable. The bound state mode of Fig. 18(b) featuresa single charge vortex in the first component and the | , i (c) in the second component. This complex is classifiedas stable for µ + ≈ [1 . , . µ + ≈ [1 . , . µ + ≈ .
191 and is classified as stable for µ + ≈ [1 . , . µ + it becomes oscillatorily unstableup to µ + ≈ .
212 where it starts featuring a dominant exponentially unstable mode. Next, all solutions of panels(d)-(h) of Fig. 18 involve the | , i (c) state in the first component and the states in the second component bifurcateat µ + ≈ . µ + ≈ . µ + ≈ . µ + ≈ .
22, and µ + ≈ . | , i (c) in the second component whereas the | , i (c) in Fig. 18(e) traps a soliton necklace asthis is the case with the second component of Fig. 18(g) and (h). In Fig. 18(f), the solution that is trapped in thesecond component can be approximated in their linear limit by a linear combination of Cartesian eigenstates [26] (seealso the second component of Fig. 13(g)). The solutions of Figs. (d)-(h) are all classified as exponentially unstable. IV. CONCLUDING REMARKS AND FUTURE CHALLENGES
In the present work, we have shed some light on the wealth and complexity of the pattern formation that arises inthe context of two-component atomic Bose-Einstein condensates, revealed by a state-of-the-art numerical technique.Naturally, some of the states that we explored have been previously considered in earlier studies such as [33, 34, 37](among others) or represent two-component extensions of single-component ones. However, several of the statesfound have not previously been discussed in the literature, to the best of our knowledge. In addition to identifyingthese states, we explored their bifurcation structure and gave a systematic mapping of the parametric regions (as afunction of the chemical potential of the second component) where the states were potentially stable, exponentiallyor oscillatorily unstable.Naturally, there are numerous directions of further work to pursue. A clear starting point is that of seeking away for obtaining a complete “cartography” of the possible nonlinear states of the model. We saw that the DCMoffers a rich repository of such states. Similarly, we argued that the linear limit and its theoretical understandingprovides plenty more. Additionally bifurcation events from the existing states lead to even more. Yet, no toolboxat the moment appears to allow a systematic classification of the full span of the highly nonlinear states in thissystem. A summary of the current techniques and their advantages and disadvantages will be a useful tool forfurther efforts. Moreover, to date we have focused on two-dimensional states, but it is particularly relevant to studythree-dimensional configurations in both single but also multi-component settings. Another direction where recently4 -10010 00.20.40.6 00.20.40.6-10 0 10-10010 -0 -10 0 10 -0 -10010 00.20.4 00.20.40.6-10 0 10-12-28 -0 -10 0 10 -0 -10010 00.20.40.6 00.20.40.6-10 0 10-12-28 -0 -10 0 10 -0 -10010 00.20.4 00.20.4-10 0 10-10010 -0 -10 0 10 -0 -10010 00.20.4 00.20.40.6-10 0 10-10010 -0 -10 0 10 -0 -10010 00.20.4 00.20.4-10 0 10-10010 -0 -10 0 10 -0 -10010 00.20.4 00.10.20.3-10 0 10-10010 -0 -10 0 10 -0 -10010 00.20.4 00.10.2-10 0 10-10010 -0 -10 0 10 -0
FIG. 18: Other branches of solutions identified via continuation from the linear limit (see text for details). Densities andphases are shown for µ + = 1 (a) and µ + = 1 .
21 (b), µ + = 1 .
26 (c), µ + = 1 . µ + = 1 .
232 (e), µ + = 1 .
27 (f), µ + = 1 . µ + = 1 . µ + would give it more branches to track, and therefore more initial guesses touse at each continuation step. Second, branch switching algorithms (as implemented in, e.g., AUTO [39]) couldbe applied to branches identified via deflation, combining their advantages. Third, deflation can be combined withnested iteration to greatly enhance its robustness [40]. Fourth, the use of more robust nonlinear solvers (improvedline searches or trust region techniques) would make the solution of deflated problems more successful on average.Together, these extensions could substantially enhance the ability of deflation to reveal previously unknown solutionsto nonlinear partial differential equations and, thus, represent a direction worth pursuing in its own right. Acknowledgments
EGC is indebted to Hans Johnston (UMass) for his endless support and guidance throughout this work as wellas providing computing resources. He extends his deepest gratitude to Eric Polizzi (UMass) and Pavel Holoborodko(Advanpix) for substantial assistance regarding the eigenvalue computations using FEAST and the MultiprecisionComputing Toolbox for MATLAB, respectively. This work is supported by EPSRC grant EP/R029423/1 (PEF), theEPSRC Centre For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) (NB). Thismaterial is based upon work supported by the US National Science Foundation under Grants No. PHY-1602994 andDMS-1809074 (PGK). PGK also acknowledges support from the Leverhulme Trust via a Visiting Fellowship and theMathematical Institute of the University of Oxford for its hospitality during part of this work.
Appendix A: State-of-the-art eigenvalue computations using FEAST
In this appendix we briefly give more details of the eigenvalue computations associated with Eq. (12). Uponconvergence of Newton’s method (with 10 − relative and absolute tolerances, respectively), the stability matrix A ofEq. (12) is computed. For the finite difference discretization employed, the matrix A is a 357 , × ,
604 sparsematrix with 2 , ,
048 non-zero elements. To compute its eigenvalues, we initially used MATLAB’s eigs commandwhich employs a Krylov-Schur method [41] to compute a subset of the eigenvalues. Although MATLAB did not raiseany warnings during the computation (i.e., flag=0 ), in some of the cases the obtained spectrum contained spuriouseigenvalues. This finding was further investigated by calculatingmax ≤ i ≤ n max (cid:13)(cid:13)(cid:13) A W ( i )R − ρ i W ( i )R (cid:13)(cid:13)(cid:13) ℓ k A k ℓ , (A1)where n max is the number of eigenpairs computed, and W ( i )R is the i th right eigenvector associated with the eigenvalue ρ i . For instance, Eq. (A1) evaluated for the solution of the branch shown in Fig. 2(c) and for µ + = 1 .
32 gives ≈ . eigs computesfirst the LU decomposition (with full pivoting and scaling) of the matrix A , which is performed very accurately.However, the matrix A itself is ill-conditioned (as we will see subsequently) and any slight change in subspace vectorsin the Krylov-Schur algorithm will lead to a large change in resulting eigenvectors [43]. As a consequence, all theaccuracy is lost by using eigs without any warning raised.To further investigate this issue, we used the Multiprecision Computing Toolbox “Advanpix” [44] which implementsan extended precision version of eigs function called mpeigs . We performed the eigenvalue computation of the branchshown in Fig. 2(c) with 34 digits (over 121 distinct values of µ + ) which took approximately 3 months of computingtime on an Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60GHz processor with 64GB of RAM. The respective results areshown with red stars in the left and right panels of Fig. 19 corresponding to the real and imaginary parts of ω ,respectively. We further checked Eq. (A1) for this computation and for, e.g., µ + = 1 .
32 we found that the residual ofEq. (A1) for 100 eigenpairs ( ρ, W R ) was ≈ . × − . This computation clearly suggests that extended precisionis capable of diminishing any small perturbations in subspace vectors in the Krylov-Schur algorithm.6 FIG. 19: The eigenfrequencies of Eq. (12) obtained by using the FEAST eigenvalue solver and mpeigs are shown with blueopen circles and red stars for comparison.
However, a new algorithm (motivated by contour integration and density-matrix representation in quantum me-chanics) for solving eigenvalue problems known as FEAST was introduced in [45] (see also [46] for details). FEASTitself combines accuracy, efficiency and robustness while exhibiting natural parallelism at multiple levels. Recently,the algorithm was generalized and applied to non-Hermitian matrices [47] (and references therein). In our presentwork, we used FEAST extensively to calculate the spectra of all branches shown with 10 − relative tolerance (on theresiduals) as the stopping criterion (e.g., for the branch of Fig. 2(c) and for µ + = 1 .
32, Eq. (A1) gives ≈ . × − for 100 eigenpairs). It should be noted that FEAST converges quite quickly (3 to 5 iterations were required in mostof the cases we studied), thus providing us with a robust tool for solving large eigenvalue problems. Its executiontime varies depending on the number of eigenvalues inside a prescribed contour. Indicatively, the computation of theeigenvalues of the branch of Fig. 2(c) took approximately 5 hours.To further demonstrate the robustness of the FEAST algorithm, we calculate the condition number of a simpleeigenvalue ρ [43] defined by κ ( ρ ) = k W R k ℓ k W L k ℓ | ( W R , W L ) | , (A2)where W L is the left eigenvector associated with the eigenvalue ρ . Direct computation of Eq. (A2) using the left andright eigenvectors computed in FEAST (inside an elliptical contour) for µ + = 1 .
32, we get a value of max ≤ i ≤ n max κ ( ρ i ) ≈ . × , again for 100 eigenpairs. The spectra of the entire branch were computed using FEAST and are shown inFig. 19 with blue open circles. The agreement of the spectra obtained using mpeigs and FEAST is clearly evident.FEAST has clearly demonstrated its robustness and accuracy (see also the discussion in [47]) and is a powerful toolfor studying the spectra of very large ill-conditioned matrices. [1] C. J. Pethick and H. Smith, Bose-Einstein condensation in Dilute Gases , Cambridge University Press (Cambridge, 2002).[2] L. P. Pitaevskii and S. Stringari,
Bose-Einstein Condensation and Superfluidity , Oxford University Press (Oxford, 2016).[3] F. Kh. Abdullaev, A. Gammal, A. M. Kamchatnov, and L. Tomio, Int. J. Mod. Phys. B , 3415 (2005).[4] D. J. Frantzeskakis, J. Phys. A: Math. Theor. , 213001 (2010).[5] A. L. Fetter and A. A. Svidzinsky, J. Phys.: Cond. Mat. , R135 (2001).[6] A. L. Fetter, Rev. Mod. Phys. , 647 (2009).[7] P. G. Kevrekidis, R. Carretero-Gonz´alez, D. J. Frantzeskakis, and I. G. Kevrekidis, Mod. Phys. Lett. B , 1481 (2004).[8] S. Komineas, Eur. Phys. J.- Spec. Topics , 133 (2007).[9] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-Gonz´alez, The Defocusing Nonlinear Schr¨odinger Equation: fromDark Solitons to Vortices and Vortex Rings (SIAM, Philadelphia, 2015).[10] P. G. Kevrekidis, D. J. Frantzeskakis, Rev. Phys. , 140 (2016).[11] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. , 1539–1542 (1998).[12] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.-J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev.Lett. , 2027–2030 (1998).[13] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H.-J. Miesner, A. P. Chikkatur, and W. Ketterle, Nature , 345–348 (1998).[14] Y. Kawaguchi and M. Ueda, Phys. Rep. , 253–381 (2012).[15] D. M. Stamper-Kurn and M. Ueda, Rev. Mod. Phys. , 1191–1244 (2013). [16] C. Becker, S. Stellmer, P. Soltan-Panahi, S. D¨orscher, M. Baumert, E.-M. Richter, J. Kronj¨ager, K. Bongs, and K.Sengstock, Nature Phys. , 496–501 (2008).[17] S. Middelkamp, J. J. Chang, C. Hamner, R. Carretero-Gonz´alez, P. G. Kevrekidis, V. Achilleos, D. J. Frantzeskakis,P. Schmelcher, and P. Engels, Phys. Lett. A , 642–646 (2011).[18] C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer, Phys. Rev. Lett. , 065302 (2011).[19] D. Yan, J. J. Chang, C. Hamner, P. G. Kevrekidis, P. Engels, V. Achilleos, D. J. Frantzeskakis, R. Carretero-Gonz´alez,and P. Schmelcher, Phys. Rev. A , 053630(2011).[20] M. A. Hoefer, J. J. Chang, C. Hamner, and P. Engels, Phys. Rev. A , 041605(R) (2011).[21] D. Yan, J. J. Chang, C. Hamner, M. Hoefer, P. G. Kevrekidis, P. Engels, V. Achilleos, D. J. Frantzeskakis, and J. Cuevas,J. Phys. B: At. Mol. Opt. Phys. , 115301 (2012).[22] D. M. Stamper-Kurn and W. Ketterle, in Coherent Atomic Matter Waves , R. Kaiser, C. Westbroook and F. David (Eds.),Springer-Verlag (Berlin, 2001), p. 139.[23] M.-S. Chang, C. D. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, Phys.Rev. Lett. , 140403 (2004).[24] M.-S. Chang, Q. Qin, W. Zhang, L. You, and M. S. Chapman, Nat. Phys. , 111 (2005).[25] T. M. Bersano, V. Gokhroo, M. A. Khamehchi, J. D’Ambroise, D. J. Frantzeskakis, P. Engels, and P. G. Kevrekidis, Phys.Rev. Lett. , 063202 (2018).[26] E. G. Charalampidis, P. G. Kevrekidis, and P. E. Farrell, Commun. Nonlinear Sci. Numer. Simulat. , 482 (2018).[27] P. E. Farrell, A. Birkisson, and S. W. Funke, SIAM J. Sci. Comp. , 2026 (2016).[28] P. E. Farrell, C. H. L. Beentjes and ´A. Birkisson, arXiv:1603.00809 (2016).[29] A. Logg, K. A. Mardal, and G. N. Wells (eds.), Automated Solution of Differential Equations by the Finite Element Method (Springer-Verlag, Berlin, 2012).[30] P. Sonneveld and M. B. van Gijzen, SIAM J. Sci. Comput. , 1035 (2008).[31] K. J. H. Law, P. G. Kevrekidis, and L. S. Tuckerman, Phys. Rev. Lett. , 160405 (2010).[32] J. Yang, Stud. Appl. Math , 133 (2012).[33] M. Pola, J. Stockhofe, P. Schmelcher, and P. G. Kevrekidis, Phys. Rev. A , 053601 (2012).[34] P. G. Kevrekidis, W. Wang, R. Carretero-Gonz´alez, and D. J. Frantzeskakis, Phys. Rev. A , 063604 (2018).[35] S. Middelkamp, P. G. Kevrekidis, D. J. Frantzeskakis, R. Carretero-Gonz´alez, and P. Schmelcher, Phys. Rev. A ,013646 (2010).[36] B. P. Anderson, P. C. Haljan, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. , 2857 (2000).[37] J. Stockhofe, P. G. Kevrekidis, D. J. Frantzeskakis, and P. Schmelcher, J. Phys. B , 191003 (2011).[38] M. J. Ablowitz, B. Prinari, and A. D. Trubatch, Discrete and Continuous Nonlinear Schr¨odinger Systems , CambridgeUniversity Press (Cambridge, 2004).[39] E. Doedel and J. P. Kern´evez,
AUTO: Software for continuation and bifurcation problems in ordinary differential equations ,Tech. Rep., California Institute of Technology (1986).[40] J. H. Adler, D. B. Emerson, P. E. Farrell, and S. P. MacLachlan, SIAM J. Sci. Comput. Numerical Solution of Large Eigenvalue Problems , 115112 (2009).[46] P. T. P. Tang and E. Polizzi, SIAM J. Matrix Anal. Appl. , 354 (2014).[47] J. Kestyn, E. Polizzi, and P. T. P. Tang, SIAM J. Sci. Comput.38