Computing stationary solutions of the two-dimensional Gross-Pitaevskii equation with Deflated continuation

In this work we employ a recently proposed bifurcation analysis technique, the deflated continuation algorithm, to compute steady-state solitary waveforms in a one-component, two dimensional nonlinear Schr\"odinger equation with a parabolic trap and repulsive interactions. Despite the fact that this system has been studied extensively, we discover a wide variety of previously unknown branches of solutions. We analyze the stability of the newly discovered branches and discuss the bifurcations that relate them to known solutions both in the near linear (Cartesian, as well as polar) and in the highly nonlinear regimes. While deflated continuation is not guaranteed to compute the full bifurcation diagram, this analysis is a potent demonstration that the algorithm can discover new nonlinear states and provide insights into the energy landscape of complex high-dimensional Hamiltonian dynamical systems.


Introduction
Over the past two decades, the pristine setting of atomic Bose-Einstein condensates (BECs) has enabled the exploration of numerous physical concepts [1,3] . One of the principal themes examined involves the interface between nonlinear wave dynamics and such atomic systems, concerning the study of so-called matter-wave solitons [4][5][6] . These coherent structures have not only been theoretically predicted but also in many cases experimentally observed. Some relevant examples include bright [8][9][10] , dark [11] and gap [12] matter-wave solitons. Higher dimensional analogues of these have also been studied; these permit a wide variety of interesting structures, such as vortices [13,14] , solitonic vortices and vortex rings [15] .
From a bifurcation theoretic perspective, the one-dimensional repulsive version of this problem is interesting but not particularly rich [16][17][18][19][20] in the customary experimental setting of a parabolic trapping potential. One can systematically trace the bifurcation of nonlinear branches from the corresponding linear ones. As the characteristic eigenvalue parameter, the so-called chemical potential, is detuned from the linear limit of the quantum harmonic oscillator, nonlinear states emerge via bifurcation from the trivial branch. However, no subsequent bifurcations occur on these branches even as the chemical potential approaches the other analytically tractable limit of large values, the Thomas-Fermi (TF) limit, where excited states feature the same profile with embedded single or multiple dark solitons.
However, the solutions of the two-dimensional repulsive case are far more intricate and intriguing. In this paper we revisit the two-dimensional repulsive case with a new numerical tool, the deflated continuation algorithm of Farrell et al. [21] . This specific problem is of significant interest for a number of reasons. It constitutes a prototypical example where numerous secondary bifurcations have been shown to occur [22,23] . These include symmetry-breaking bifurcations and provide the potential for genuinely complex order parameter states bearing vorticity. These bifurcations are also crucial in modifying the stability characteristics of the states at hand. Hence, as we will see below, this problem provides a very rich testbed for investigating the effectiveness of deflated continuation on an important and intensely studied problem [20,[22][23][24][25][26][27][28][29] . From a mathematical perspective, the analysis of the relevant bifurcations and the resulting vortex-bearing states are still a topic of recent investigation [30] .
At the heart of deflated continuation is a deflation technique for computing previously unknown solutions of differential equations [31] , and we discuss this first. Suppose we have a nonlinear problem F (φ) = 0 , and one solution φ 1 has been found via Newton's method from an initial guess φ . The central idea of deflation is to eliminate the known solution from consideration ; we construct a new nonlinear problem G (φ) = 0 that possesses all of the solutions of F (φ) = 0 , except for the known solution φ 1 . We say that the solution φ 1 has been deflated. More precisely, G (φ) = 0 has the property that Newton's method is guaranteed not to converge to φ 1 from any initial guess, even φ . Thus, if Newton's method is initialized from the same initial guess φ , and it converges, it will converge to a distinct solution φ 2 = φ 1 . The process can then be repeated by deflating each solution found until Newton's method no longer converges to any solution within a fixed number of iterations. In this way, many solutions to the same problem can be discovered from a single initial guess. However, there is no guarantee that all of the relevant solutions (for a given parameter value) will be obtained in this way, as Newton's method may fail to converge when initialized far from a solution.
The deflated problem is constructed via the application of a deflation operator to the residual. Suppose F : U → V is a nonlinear map between Banach spaces, and φ 1 is an isolated root of F , i.e. a root of multiplicity 1 with F ( φ 1 ) invertible 1 . In previous work [31] , G : U → V was constructed via where · U is the norm on U . The essential idea is that φ − φ 1 −2 U approaches infinity as φ → φ 1 faster than F ( φ) approaches 0, and hence lim inf ensuring that Newton's method will not converge to φ 1 when applied to G . In the present case, this idea must be modified slightly, as the solutions are no longer isolated: if φ 1 is a solution, then so is its phase shift ψ θ = e iθ φ 1 for any phase θ . Given knowledge of a solution φ 1 , we therefore wish to deflate the entire group orbit { ψ θ : θ ∈ [0, 2 π )}. This is achieved by constructing the deflated problem via where | φ| represents the amplitude of the complex-valued wavefunction φ. As the amplitude is invariant under phase shift, this modified deflation operator eliminates the entire group orbit, ensuring nonconvergence to any solution trivially related to a known solution. The norm chosen for U is the H 1 (D ; R ) norm, where D is the domain on which the PDE is posed.
This deflation idea is extended to an algorithm for bifurcation analysis, called deflated continuation, as follows. Suppose a set of solutions S ( μ) is known for the nonlinear problem F (φ; μ) = 0 at a given parameter value μ ∈ R , and we wish to compute the solutions at a modified parameter value S(μ + μ) . In the first phase of the algorithm, known branches are continued: each solution is used as initial guess for the nonlinear problem at μ + μ, and the resulting solution deflated. In the second phase, new solutions are sought: each φ ∈ S ( μ) is used again as initial guess for the nonlinear problem at μ + μ. In this way, we allow for the discovery of new branches that have bifurcated between μ and μ + μ. If a search is successful, the initial guess φ is used again. Once the algorithm has exhausted all initial guesses in S ( μ), it proceeds to compute S(μ + 2 μ) from S(μ + μ) , and so on. For a full description of the algorithm, see Appendix A .
There are two central advantages of deflated continuation. First, the algorithm is capable of discovering disconnected branches, those not continuously connected to known solutions. Second, the algorithm can scale to very fine discretizations of PDEs. The only expensive subroutine required by deflated continuation is the computation of the Newton step of the undeflated system; thus, if a good preconditioner is available for the Newton system, the bifurcation diagram of the system can be efficiently computed. This will be of central importance in future work on the corresponding three-dimensional problem (see e.g. [5,Chapter 4] ).
In the two-dimensional problem of interest here we have extensive knowledge of the solutions at the linear limit, as described in Section 2 . We augment deflated continuation with this knowledge by using these linear solutions as initial guesses near the bifurcation from the trivial branch. This augmentation identifies a handful of additional solutions that deflated continuation alone misses, as will be discussed in Section 4 .
We apply deflated continuation to this complex problem in the hope of gaining new insight into its families of solutions. We will complement the solutions identified via deflation with a stability analysis, aiming at a systematic map of the newly discovered bifurcations and associated stability changes. In Section 2 , we will provide a brief overview of the theoretical setup of both the existence and stability problem. This will explore the analytically tractable linear limit which will subsequently serve as a way of potentially seeding nonlinear solutions away from that limit. In Section 3 , we will offer a systematic classification of our numerical results. Based on a detailed comparison with earlier works including [23,28,32] , we will provide a wealth of novel families and associated bifurcations. Finally, in Section 4 , we summarize our main findings and discuss the directions of future study arising from this work.

The model and setup
We consider the nonlinear Schrödinger (NLS) equation in (2 + 1) dimensions (two spatial and one temporal) written in dimensionless form (see, e.g., [4,5] for details of the nondimensionalization) as Here ∇ 2 stands for the standard Laplace operator in two spatial dimensions and the external potential V ( r ) takes the standard harmonic form of V (r ) = 1 2 2 | r | 2 , with | r | 2 = x 2 + y 2 and the normalized trap strength . The latter represents the ratio of trappings along and transverse to the plane and should thus be 1; in what follows we will fix = 0 . 2 . In the context of BECs, the field (r , t ) : D × R + → C in Eq. (4) represents the macroscopic wave function with D ⊆ R 2 being the (two-dimensional) spatial domain. The sign of the nonlinear term is chosen to reflect the self-repulsive nature of the interatomic interactions considered herein. The construction of stationary solutions is based on the well-known ansatz with chemical potential μ ≥ . By substituting Eq. (5) into (4) , we obtain the stationary equation where F ( φ; μ) stands for the nonlinear set of equations utilized in our Newton solvers as well as in the deflated continuation method. The equation is discretized with piecewise linear finite elements for the real and imaginary components using FEniCS [2] . The problem is posed on the domain D = (−12 , 12) 2 and homogeneous Dirichlet conditions are imposed. This choice of domain is made to ensure that the influence of truncating the domain is negligible, as the support of the solutions remains far from the boundary. We use the following diagnostic to summarize the parametric dependence of each steady state branch on μ. The above integral represents the number of atoms in the BEC, considered as a function of the chemical potential μ. When N → 0, the nonlinearity of the problem becomes irrelevant and the states bifurcate from the respective linear limit. We will also use the difference in the number of atoms with respect to a given state N in order to highlight the emergence of bifurcations. This diagnostic, examined in [23] , transparently illustrates the origin of different branches. In the N → 0 limit, we can decompose the modes in Cartesian form [7] as being proportional to: where H m, n are the Hermite polynomials with m, n > 0 being the associated quantum numbers of the harmonic oscillator. These linear eigenfunctions have corresponding eigenvalues μ = E m,n := (m + n + 1) . On the other hand, the linear eigenfunctions mentioned above can also be expressed in polar coordinates via with eigenvalues μ = E k,l := (1 + | l| + 2 k ) . The parameters, or equivalently, quantum numbers l and k denote the eigenvalue of the ( z -component of the) angular momentum operator and the number of radial zeros of the corresponding eigenfunction q k, l , respectively. This radial part of the relevant eigenfunction is given by [33] q k,l ∼ r l L l where L l k are the associated Laguerre polynomials. The subscripts in | m, n ( c ) and | k, l ( p ) stand for the Cartesian (c) and polar (p) representations respectively. For instance, the ground state in the polar decomposition can be represented by | 0 , 0 (p) (i.e., in the absence of a zero crossing in the radial profile and without angular momentum), and the first excited state with angular momentum (the single charge vortex) by | 0 , ±1 (p) . Similarly, the second excited states can be represented by | 0 , ±2 (p) (a doubly charged vortex) and | 1 , 0 (p) (the ring dark soliton). In most of what follows, we will use the Cartesian notation, but we will occasionally resort to the polar decomposition where convenient.
Once the relevant solutions have been identified, we wish to investigate their stability. We assume the perturbation ansatz around a stationary solution φ 0 to be of the following form: where ω is the (complex) eigenfrequency, ε is a (formal) small amplitude of the perturbation, and the asterisk stands for complex conjugation. Inserting Eq. (11) into the NLS Eq. (4) , we obtain at order ε an eigenvalue problem written in the following matrix form: with eigenfrequencies ρ = −ω (the eigenvalues λ = iω), eigenvectors V = [ a b] T , and matrix elements given by The steady state φ 0 is classified as stable in this Hamiltonian system if no eigenfrequency ω = ω r + i ω i has a nonvanishing imaginary part ω i ; this is because, given the Hamiltonian nature, if ω is an eigenfrequency, so is −ω, ω and −ω . The scenario of stability will be depicted by a solid blue line in the bifurcation diagrams presented below. All stable stationary solutions to be discussed in Section 3 are conveniently summarized in Table 1 of the Appendix B . On the other hand, when the solution becomes unstable, two types of instabilities can be identified: i) exponential instabilities characterized by a pair of imaginary eigenfrequencies with zero real part, and ii) oscillatory instabilities characterized by a complex eigenfrequency quartet. These two scenarios are depicted by dashed-dotted red and green lines respectively in the bifurcation diagrams that follow to highlight the nature of the dominant unstable mode; a transition between these two colors will thus signal a change in the dominant instability type.
We are now ready to describe the different solutions produced by the deflation technique.

Numerical results
The simplest state of the system is its ground state |0, 0 (c) with eigenvalue μ = at the linear limit. As this state is generically stable [5] (i.e., for all values of μ from the linear limit to the Thomas-Fermi regime), no bifurcations occur from it. This solution is well-known and we do not examine it further.

Bifurcations from μ = 2
From the point of view of bifurcation analysis, the first interesting events occur at μ = 2 , with n + m = 1 . As is wellknown [22,23,34] , two branches bifurcate from this point ( Fig. 1 ). One is the 1-dark soliton stripe |1, 0 (c) , Fig. 1 (a). This is the two-dimensional generalization of the one-dimensional dark soliton analyzed, e.g., in [35] . The other is the single charge (i.e., unit vorticity) vortex of Fig. 1 (b). In Cartesian coordinates this is described as the linear combination | 1 , 0 ( c ) + i | 0 , 1 ( c ) ; in polar coordinates it is |0, 1 (p) . Among the two, the vortex is very robust, incurring no instabilities [5,22] . We note in passing that a corresponding analysis for the so-called bright vortex in the attractive interaction case can be found, e.g., in [36] ; there the vortex is only found to be stable for sufficiently small atom number.
As the chemical potential increases the stripe progressively approaches a rectilinear dark soliton that is well-known to be subject to a transverse (modulational) instability [37] . In fact, there is a whole cascade of such instabilities, arising in the form of pitchfork bifurcations from the stripe [5,22,23] . The first of these bifurcations gives rise to the emergence of the vortex dipole state, Fig. 1 (c). This state is well-known and has been studied experimentally [38,39] . It is dynamically stable except for a narrow interval of oscillatory instability (associated with a Hamiltonian Hopf bifurcation), as previously noted [22,23] . The next bifurcation gives rise to a vortex tripole, Fig. 1 (d). By this stage the stripe branch is unstable and the vortex tripole inherits this instability. This is a configuration with three vortices of alternating charge (+ , −, +) or (−, + , −) .
This has also been identified in experiments [40] and explored in simulations [41,42] . This pattern of bifurcations continues to higher excited states for larger values of μ; we do not pursue these bifurcations further, given their earlier analysis, e.g., in [22,23] .

Bifurcations from μ = 3
We now turn to the considerably more complicated case of bifurcations from μ = 3 with n + m = 2 . These bifurcations are also mostly well-known (although with some important twists to be discussed below), and are summarized in Figs. 2-  Somewhat surprisingly, the double tripole bifurcation (leading to a state with 6 vortices of alternating charges) corresponds to a bifurcating stable branch, at least when the latter states first emerge; for higher μ they become oscillatorily unstable, acquiring a complex eigenfrequency quartet. The subsequent bifurcations arise from the already unstable double soliton stripe branch and hence the bifurcating branches are also unstable. For a demonstration of the relevant stability properties, see the bottom left panel of Fig. 2 , while the bottom right uses the number of atom difference from the 2-dark soliton stripe branch as a diagnostic to display the occurrences of the different bifurcations (as the daughter branches depart from zero in this quantity).
The bifurcation of the 6 vortex state (numerically occurring at μ = 0 . 86 ; panel 2 (b)) precedes that of the 4 vortex one (numerically occurring around μ = 0 . 89 ; panel 2 (c)). The former bifurcation in the formulation of [23] comes from the combination of |2, 0 (c) with |0, 3 (c) , with a π /2 phase shift, while the latter emerges from the symmetry-breaking event involving |2, 0 (c) with |1, 2 (c) (again with a π /2 phase shift). It is also intriguing that the theoretical prediction of these bifurcations based on the two-mode theory of [23] occurs at 503 /113 ≈ 0.89 and 265 /61 ≈ 0.87, respectively, i.e., very close to the computationally obtained values, although in reverse order. Fig. 3 (a) considers |1, 1 (c) and its subsequent bifurcations. The first bifurcation from |1, 1 (c) gives rise to a branch which features a dark soliton stripe together with two same charge vortices, as shown in Fig. 3 (b). Nearly concurrent to this bifurcation is the emergence of a branch involving a "diagonal" set of 6 vortices; see Fig. 3 (c). As discussed in [23] , the branch of solutions in Fig. 2

(b) collides and disappears
(in a saddle-center bifurcation) with that of Fig. 3 (c), as μ is increased. A subsequent bifurcation gives rise to a branch with a vortex of charge l = −2 in the middle and four surrounding vortices of charge l = 1 , Fig. 3 (d). Explicit algebraic conditions for such states have been obtained via a generating function approach in the large-density Thomas-Fermi limit [43] , as have the rectilinear vortex states of Fig. 1 . It is also interesting to note, as illustrated in the bottom panel of Fig. 3 , the resulting branches out of these bifurcations also feature intervals of oscillatory instabilities. However, the principal branch |1, 1 (c) out of which these bifurcations arise is unstable, and hence all the bifurcating branches inherit this instability. This includes the 8 diagonal vortex state Fig. 3 (e). A more complex structure emerges in the context of solutions with radially symmetric density and their bifurcations that are analyzed in Fig. 4 . The first such solution to consider is the ring dark solution state, Fig. 4 (a), which is well-known and has been extensively studied (for a recent discussion, see [44] and the references therein). In Cartesian notation, this is | 2 , 0 ( c ) + | 0 , 2 ( c ) , while in radial notation this is the | 1 , 0 ( p ) state. A systematic study of the stability of this and related states from the linear limit onwards was conducted in [32] , showing that its degeneracy with the vortex quadrupole of Fig. 4 (b) leads to an immediate quadrupolar instability through a real eigenvalue pair for this mode. In contrast, the vortex quadrupole is generically stable, aside from a finite interval of oscillatory instability [23] . The ring dark soliton, progressively becomes more unstable to undulations with higher wavenumbers. The first to emerge is a hexapolar mode, giving rise (through a pitchfork bifurcation) to the vortex hexagons of Fig. 4 (c); these may for larger values of μ also possess oscillatory instabilities (as denoted in the bottom panel). This pattern continues with an octapolar mode leading to vortex octagons and so on.
Another mode bifurcating here from the linear limit is the doubly charged vortex | 0 , 2 ( p ) of Fig. 4 (d). This state is unstable from the linear limit onwards through a sequence of intervals of oscillatory instabilities originally examined in [45] and subsequently retraced in a variety of publications [29,46] . This state can also be represented via a combination of Cartesian eigenstates as | 2 , 0 ( c ) − | 0 , 2 ( c ) + 2 i | 1 , 1 ( c ) . Out of this branch bifurcates branch (e) which bears three vortices of the same charge in the periphery and one of opposite charge at the center (hence has the same total charge of 2). Intriguingly, due to the spherical symmetry of the solution, despite the pitchfork nature of the bifurcation, the branch (e) has a pair of eigenvalues of the linearization at the origin, and no genuinely imaginary eigenfrequencies. It is also interesting to mention here that the bifurcation of (e) essentially coincides (for our parametric resolution of steps of 0.01 in the chemical potential) with the stabilization against oscillatory instabilities of the charge 2 branch. Finally, the branch (f) emanating from the linear limit and being subject to oscillatory instabilities can be approximated by | 2 , 0 ( c ) + i | 1 , 1 ( c ) . This solution also seems to "harbor" 4 vortices, although 2 are more clearly observable close to the condensate center, while the two others are less discernible, merging with the background. It is worth noting here that the systematic classification of [32] identified the solutions stemming from the linear limit (including the ring of Fig. 4 (a), the multi-pole of Fig. 3 (a), the radially symmetric vortex of Fig. 4 (d) and the vortex necklace (in this case, a vortex quadrupole) of Fig. 4 (b)). Nevertheless, a remarkable feature of our analysis is that the branch of Fig. 4 (f) appears to have never been previously discussed, to the best of our knowledge. This branch is subject to oscillatory instabilities, as shown in the bottom panel of Fig. 4 .
The analysis presented heretofore considered bifurcations that have been analyzed in earlier works; we have successfully validated the method by comparing it to the known bifurcation diagram of this problem and in the process have unraveled novel branches of solutions, such as Figs. 2 (c), 3 (b) and 4 (f).

Bifurcations from μ = 4
We next examine the bifurcations emanating from μ = 4 with n + m = 3 , plotted in Figs. 5-7 . There are states that we can immediately recognize as emerging out of the Cartesian linear limit, the |3, 0 (c) state of Fig. 5 (a) and the |2, 1 (c) state of Fig. 5 (b). The |3, 0 (c) state undergoes subsequent bifurcations leading to real eigenvalues and symmetry-breaking destabilizations, e.g., the six vortex and dark soliton stripe state of Fig. 5 (o) (see, also, Fig. 6 of [23] ). This pitchfork bifurcation, due to the admixture of |3, 0 (c) with |1, 3 (c) with a π /2 shift (see the relevant theory of [23] ), leads to the destabilization of the former branch. Furthermore, the |2, 1 (c) mode gives birth to a state having two dark solitons and two vortices -through admixture with |4, 0 (c) again with a π /2 phase shift-. The bifurcating state appears to be exponentially unstable, inheriting the instability of its "parent branch", over the full parametric interval of μ and is depicted in Fig. 5 (c). In addition, the algorithm has discovered states that can be naturally expressed as linear combinations of the Cartesian eigenstates. The state depicted in Fig. 5 (d) can be expressed as | 3 , 0 ( c ) + i | 0 , 3 ( c ) in the linear limit and represents a lattice of 9 vortices of alternating charge and was previously characterized via algebraic conditions [43] . This solution appears to be oscillatorily unstable throughout our computations.
The states depicted in Fig. 5

(e) and (f) can be approximated by the linear combinations
To the best of our knowledge, these solutions have not been previously identified (cf.    [23] ), yet adhere to the formulation whereby more complex nonlinear solutions generalize linear combinations of simpler linear eigenstates. Furthermore, in a manner reminiscent of Fig. 4 (f), such solutions feature a pattern of vortices (along a line or on a cross, respectively) close to the center of the trap and another one further away. Both branches appear to be oscillatorily unstable, except for an interval where the former branch also develops exponential instabilities associated with real eigenvalues for μ ࣡ 1.21.  There are also states such as the ring-vortex state of Fig. 5 (g) (previously examined, e.g., in [28,29] ) and the well-known triple charge vortex of Fig. 5 (h) that are best described in polar notation, | 1 , 1 ( p ) and | 0 , 3 ( p ) respectively. Both branches are oscillatorily unstable. In the terminology of [32] , Fig. 5 (i) depicts a multi-pole, a real solution of the form q 0, 3 ( r )cos (3 θ ), whereas a complex combination of q 1, 1 ( r ) with q 0, 3 ( r )sin (3 θ ) will produce the vortex necklace of Fig. 5 (j). Furthermore, the state of Fig. 5 (i) undergoes a symmetry-breaking bifurcation -through an admixture with |0, 4 (c) with a π /2 phase shift -leading to the state of Fig. 5 (k). The latter possesses two isolated vortices (two more appear "hidden" in the region of vanishing density and can be discerned in the associated phase plot). It is exponentially unstable, inheriting the instability of its multi-pole parent branch. This is with the exception of a narrow window (in μ) where the oscillatory instability dominates. The branch of Fig. 5 (l) is associated with q 1, 1 ( r )cos ( θ ) i.e., another type of multi-pole (or | 2 , 1 ( c ) + | 0 , 3 ( c ) in the Cartesian format) and is referred to as a soliton in the terminology of [33] . However, this state becomes subject to instability similar to that leading to branch Fig. 4 (c), associated with a hexagonal mode. As a result, the daughter branch of Fig. 5 (m) emerges possessing four isolated vortices at the top and bottom and two additional charge 2 vortices along the x -axis.
The bifurcation diagram in Fig. 7 sheds light on the potential stability of the branches, as well as on the bifurcations that arise. More concretely, we observe that almost all the branches are now dynamically unstable, a feature that is not surprising given the highly excited nature of the states. Nevertheless, the branch |3, 0 (c) of Fig. 5 (a) possesses intervals of stability. The top left panel of Fig. 7 sheds light on relevant bifurcations including the fact that the 9-vortex state of Fig. 5 (d) emerges from the ring-vortex branch of Fig. 5 (g), while the vortex necklace of Fig. 5 (j) and the state of Fig. 5 (k) emerge as bifurcations from the soliton necklace of Fig. 5 (i). Finally, the 6-vortex plus soliton state of Fig. 5 (o) emerges as a bifurcation from the 3-dark soliton stripe state of Fig. 5 (a).
Arguably, of particular note are the complex patterns of Figs. 5 (n) and 6 (p), both emerging from the linear limit. The state of Fig. 5 (n) represents a vortex pattern involving 9 vortices (and total charge 3) where seven of the vortices are involved in an elaborate vortical H-shape, while the other two are distinct. This state (and the one shown in Fig. 6 (p)) can be classified as emerging from a complex (both literally and figuratively) combination at the linear limit. To unveil the relevant superposition, we project the mode -in the immediate vicinity of the linear limit -onto the fundamental modes | m, n (c) (again, with m + n = 3 ). This gives rise to the combination: α( | 3 , 0 (c) + i | 2 , 1 (c) ) + iα * ( | 0 , 3 (c) − i | 1 , 2 (c) ) with α being a suitable complex prefactor; in the case we considered, a typical value of α was found as: α ≈ −0 . 15 + 0 . 5 i . In a similar vein, the projection of the state of Fig. 6 (p) onto the fundamental modes suggests the combination: α( | 3 , 0 (c) − i | 1 , 2 (c) ) + α * ( | 0 , 3 (c) + i | 2 , 1 (c) ) with α ≈ 0 . 14 + 0 . 04 i in this case. To the best of our knowledge, such patterns have also not been previously identified and are direct by-products of the use of deflated continuation.
These bifurcations are best understood by plotting the diagnostic N , where the base solution is taken to be the 3-(planar) dark soliton branch |3, 0 (c) . The resulting bifurcation diagram is presented in the bottom panels of Fig. 7 .

Bifurcations from μ = 5
In the case of μ = 5 , the wealth of relevant states is even greater. We first discuss the states that are naturally expressed in Cartesian format. Fig. 8 (a)-(c) are associated with the | 2 , 2 ( c ) , |3, 1 (c) and | 0 , 4 ( c ) states respectively. All of these branches are potentially subject to exponential instabilities except for the branch (c) which bears oscillatory instabilities for μ ࣠ 1.235 and for μ ࣠ 1.06 the waveform is stabilized. Additional states can be produced from linear combinations of the Cartesian eigenstates. In particular, the states of Figs. 8 (d) and (e) can be approximated in their linear limit as | 4 , 0 ( c ) + | 2 , 2 ( c ) and | 4 , 0 ( c ) − | 0 , 4 ( c ) respectively. The former can be thought of as a double solution. Both solutions appear to be exponentially unstable in our computations. In addition, the branch of Fig. 8 (f) can be characterized by a linear combination of the . This branch appears to be oscillatorily unstable in our stability analysis.
Other states are more naturally classified in the polar representation. For instance, Fig. 8 (g) corresponds to the mode | 1 , 2 ( p ) , associated with an oscillatorily unstable double vortex at the origin, bearing also a ring dark soliton in the periphery. Similarly, Fig. 8 (h) can be represented as | 2 , 0 ( p ) and corresponds to a highly unstable double ring configuration. The vortex of charge four depicted in Fig. 8 (i) can be represented as | 0 , 4 ( p ) and is oscillatorily unstable; all of these states could be identified in the polar decomposition of [28,29] . Other states are generalizations of ones that we identified in Fig. 5 . Fig. 8 (j) depicts a multi-pole; this is described by q 0, 4 ( r )cos (4 θ ) in its linear limit and is subject to exponential instabilities.
Additional states can be produced from linear combinations of the polar eigenstates. Fig. 8 (k) can be represented by a complex superposition of the double ring configuration and the soliton necklace; it corresponds to the | 2 , 0 ( p ) + iq 0 , 4 (r) cos (4 θ ) state which is oscillatorily unstable. This is a canonical example of a vortex necklace in the terminology of [32] (and as can also be seen in the figure; cf. Fig. 2(d) of [32] ). Other states are more difficult to classify, although it is still plausible to classify them by using "exotic" combinations of Cartesian and polar eigenmodes. For instance, Fig. 8 (l) depicts a solution that can be described by q 0 , 4 (r) cos (4 θ ) − 2 | 2 , 2 ( c ) + | 0 , 4 ( c ) and is subject to exponential instabilities.
This real solution is a soliton necklace in the terminology of [32] . Fig. 8 (m) corresponds to the "curved" variant of Fig. 8 (e) (subject to an exponential instability, as well). The vortex necklace of Fig. 8 (n) appears to be subject only to oscillatory instabilities. This is an elaborate pattern, once again revealed by the technique of deflation, that may be thought of as consisting of 4 vortex triplets (in a Y shape) each of which adds a charge of 1 to the total charge of 4 within the structure. In terms of stability, we observe in Fig. 9 that the branch (c) is linearly stable near the linear limit.

Bifurcations from μ = 6
Finally, we examine an even more "exotic" set of solutions, the states bifurcating from the linear limit at μ = 6 . Here too we observe that most of the states are oscillatorily unstable for sufficiently large μ, with the exception of the branches of the solutions depicted in Fig. 10 (a), (g)-(i), (m) and (n) that seem to be dominated by an exponential instability.
Many of the bifurcating states can be naturally classified in Cartesian form. Fig. 10 (a) represents the |4, 1 (c) state, while Fig. 10 (b) represents the |5, 0 (c) state. The states depicted in Figs. 10 (c), (d), (e)-(g), (h) and (i) can be approximated by linear combinations of the Cartesian eigenstates; they can approximated as | 1 , solution with the two radii being concentric), | 3 , 2 ( c ) + | 1 , 4 ( c ) , | 5 , 0 ( c ) + | 3 , 2 ( c ) (triple solution), respectively. It should be noted that none of these solutions has previously appeared in a systematic fashion in the literature, to the best of our knowledge.
It is already clear at this stage that the full classification of the pertinent solutions becomes an extremely cumbersome task. This motivates the application of methods from bifurcation analysis to yield a comprehensive perspective that might not otherwise be available. In this case, deflated continuation has made it possible to unravel a wide variety of branches that had not been previously identified.  As regards the stability of the branches, Fig. 11 suggests that the branch in Fig. 10 (b), i.e., the |5, 0 (c) state is stable over a narrow parametric interval in μ.

Concluding remarks and future challenges
In this work we applied deflated continuation to the analysis of the solutions of a Bose-Einstein condensate in a twodimensional isotropic parabolic trap. This problem has a well-understood linear limit in which the eigenmodes can be identified in closed analytical form, offering guidance for what to expect within the nonlinear regime. There has been a wide range of publications on this system that have revealed a broad spectrum of nonlinear excitations and their stability characteristics. However, the number of solutions and their complexity increase significantly as the chemical potential increases and analytical calculations become rather tedious. Deflated continuation is therefore extremely useful in yielding insight into the problem. We have identified branches (from the linear limit) and configurations bearing complex vortex patterns as genuine (although often unstable) solutions of the system. We have also identified the bifurcations involving such configurations including, to give but one example, the elaborate multi-vortex patterns that arise from the symmetry breaking e.g, for l = 2 ( Fig. 3 (d) and (e)), l = 3 ( Fig. 4 (c)), and so on.
The use of deflated continuation enabled us to identify a wide range of branches that had not been previously numerically constructed or continued. Nevertheless, it should be highlighted that not all branches obtained herein were identified using this method alone; deflated continuation is not guaranteed to find all solutions to a given nonlinear problem. In particular, the solutions depicted in Figs. 2 (d), 3 (c), 5 (e) and (o) as well as Fig. 8 (h) and (k), and 10 (c)-(e) were obtained via alternative techniques. More specifically, the branches of 2 (d), 3 (c), 5 (o) were obtained via Newton's method with an initial guess consisting of the (unstable) parent branch solution itself and a small perturbation along the eigendirection of the associated unstable mode of the parent branch. For the rest of the above mentioned branches ( Figs. 5 (e), 8 (h) and (k), and 10 (c)-(e)) we made use of the underlying eigenfunctions at the linear limit. This serves as an important caveat: deflated continuation should not be thought of as a universal method that blindly reveals all solutions to a particular bifurcation problem, but as a useful tool that becomes even more powerful when combined with physical insight. In particular, a deep understanding of the underlying physics of the system (in this case knowledge of the linear limit in both Cartesian and polar coordinates) remains of paramount value in uncovering the complexity of its landscape of stationary solutions. It is conceivable that the robustness of deflated continuation could be improved by reducing the continuation step-size or by employing more robustly globalized nonlinear solvers, thus enabling the method to detect more branches. This issue merits further computational investigation. Nevertheless, the significant success of deflated continuation in the present setting suggests that it is well suited to identifying steady states in multi-dimensional PDEs with such energy landscapes, providing insights on the connections between different extrema (and saddle points) in them. In the context of BECs there exists a wide array of problems that are very much worth pursuing. A natural extension is to attempt to generalize the methodology to multi-component BECs [5,47] . Another important extension is to three-dimensional configurations in a single component, including dark solitons, vortex rings, vortex lines, Hopfions etc. [33,4 8,4 9] . A further generalization would be to three-dimensional multi-component settings, where structures such as skyrmions and merons arise [50][51][52] . Of large interest in recent years have also been directions associated with media that feature an inhomogeneous defocusing nonlinearity [53][54][55] . In these, different types of multi-poles and vortex clusters have been demonstrated [54,55] and seeking additional solutions via the deflated continuation method constitutes a natural approach. Such studies are currently in progress and will be reported in future publications. Fig. 10. Same as Fig. 8 but for states bifurcating from μ = 6 . First, third and fifth as well as second, fourth and sixth rows correspond to plots of the density profiles and phases, respectively, of the (a) | 4 ,   11. Same as Fig. 9 but for states bifurcating from μ = 6 . The number of atoms N as a function of μ.
In this appendix we give more details of the algorithm used for bifurcation analysis. We consider the equation where in our case F is given by (6) . The central task of bifurcation analysis is the calculation of the solutions φ ∈ U as a function of μ ∈ M , with a particular emphasis on identifying values of μ where the number of solutions varies. Deflated continuation is an algorithm for this task. We denote by S : M ⇒ U the map that associates the known solutions (a subset of U ) to a parameter value μ.
Suppose that a single solution φ 0 is known for an initial parameter value μ 0 . For concreteness, suppose that the initial chemical potential lies between the first and second bifurcation ( < μ < 2 ) and that the ground state φ 0 , the only solution to the equation, is known. That is, the initial data supplied to the algorithm is the fact that S(μ 0 ) = { φ 0 } . The central step in deflated continuation is the construction of an estimate for S(μ + μ) from knowledge of S ( μ). That is, using knowledge of the solutions at a previous (easier) value of the parameter, we attempt to identify a set of solutions for the next (harder) parameter value. This process repeats until the algorithm reaches the end of the interval of interest. The algorithm is described in Algorithm 1 , adapted from Farrell et al. [21] .
The central step of calculating S(μ + μ) from S ( μ) is divided into two phases. In the first phase (lines 5-9), all of the known branches are continued to the next parameter value. As each solution is found, it is recorded (line 8) and deflated (line 9). In the next phase (lines 10-18), new solutions are sought, i.e. those solutions on branches that have not yet been discovered. Each previous solution φ ∈ S ( μ) is used again as initial guess: as the continuation of each branch has already been deflated, if Newton's method converges, it will converge to a previously undiscovered branch (line 13). If a search step is successful, then the newly discovered solution is recorded and the initial guess that yielded the new branch is attempted once more. Once Newton's method has failed from each initial guess available, we post-process each φ ∈ S(μ + μ) to compute its stability as described in Section 2 (lines [19][20]. Alternatively, this could be done offline, as the bifurcation analysis algorithm does not use the results of the stability calculation. As a concrete example, we consider the application of deflated continuation to calculating the roots of the polynomial as a function of μ ∈ [ −10 , 10] . The other parameters are fixed to (α, β, γ ) = (5 , 1 / 2 , −5) . The bifurcation diagram for this problem is shown in Fig. 12 .

5:
for ˜ φ ∈ S(μ) do Continue known branches. 6: Apply Newton's method to G from initial guess ˜ φ. 7: if solution φ * found then 8: S(μ + μ) ← S(μ + μ) ∪ { φ * } Record solution. 9: Deflate solution φ * . 10: for ˜ φ ∈ S(μ) do Seek new branches. 11: success ← true 12: while success do 13: Apply Newton's method to G from initial guess ˜ φ. 14: if solution φ * found then New branch found. 15: Record solution. 16: Deflate solution φ * . 17: else 18: success ← false 19: for φ * ∈ S(μ + μ) do Post-processing. 20: Compute the stability of φ * . 21: μ ← μ + μ 22: return S in the second phase, the fold point of the upper branch is located (Newton converges to the solution x ≈ 1.512, from the initial guess x ≈ −0 . 8896 , once the solutions x ≈ −0 . 8917 and x ≈ −2 . 144 have been deflated). At the next iteration two solutions are identified, corresponding to the two upper branches. The algorithm then continues all four branches until μ = 0 . 6 ; in the continuation step from μ = 0 . 6 to μ = 0 . 61 , the lower two branches fail to converge, as they have been extinguished at a fold bifurcation. The algorithm then carries the upper two branches to μ = 10 without discovering any new branches, as the problem permits no more solutions. While this example is entirely straightforward, and the solutions of polynomials can be calculated far more efficiently with other algorithms, the algorithm carries over to calculating the bifurcation diagrams of partial differential equations. This example illustrates the essential properties of the algorithm: it is capable of computing disconnected bifurcation diagrams, and requires only the implementation of Newton's method for the problem under consideration. We therefore expect it will have a wide range of applicability to other physical problems.

Appendix B. Summary of stable states
In this appendix, all the stable solutions discussed in Section 3 are summarized in Table 1 . Note that we indicate the multiples of in the first column where the bifurcation happens in order to highlight the relevant bifurcation points and to clearly indicate the value of for which our table is given. Furthermore, the third and fourth columns of the table indicate the figure (such that these stable states can be located in Section 3 ) and the parametric interval of stability in μ, respectively. We point out that the upper bound in μ that we considered here is μ = 1 . 32 .