You are an expert at summarizing long articles. Proceed to summarize the following text: 
the pristine setting of atomic bose - einstein condensates ( becs ) has offered a significant playground for the examination of numerous physical concepts  @xcite . one of the flourishing directions has been at the interface of the theory of nonlinear waves and such atomic ( as well as optical ) systems , concerning , in particular , the study of matter - wave solitons  @xcite . such coherent nonlinear structures have been observed in experiments and studied extensively in theory . prototypical examples of pertinent studies include but are not limited to bright  @xcite , dark  @xcite and gap  @xcite matter - wave solitons , as well as vortices  @xcite , solitonic vortices and vortex rings  @xcite .    among these diverse excitations , dark solitons in repulsive becs have enjoyed a considerable amount of attention in effectively one - dimensional settings due to numerous experiments leading ( especially , more recently ) to their well - controlled creation  @xcite . moreover , numerous works considered dark solitons in higher - dimensional settings in order to experimentally explore their instability leading to the formation of vortex rings and vortex lines , as illustrated , e.g. , in refs .  @xcite . more recently , such excitations have also been observed in fermionic superfluids  @xcite . on the other hand , multi - dimensional ( i.e. , non - planar ) variants of dark solitons , were studied as well cf . e.g. , chapters 7 and 8 of ref . examples of such structures are radially symmetric solitons , such as the ring dark solitons ( rds ) , initially proposed in atomic becs in ref . @xcite ( see also refs . @xcite for subsequent studies ) . nevertheless , such structures were studied to a far lesser extend , arguably due to their generic identified instability , which was found to give rise to states such as alternating charge vortex polygons . notice that , as was recently shown in ref . @xcite , stabilization of rds is possible upon employing a proper potential barrier . the three - dimensional non - planar generalization of dark solitons , namely spherical shell solitons , has received very limited attention . the analysis of ref . @xcite in the case of box potentials ( rather than the typical parabolic trap setting of becs ) revealed their potential existence , but did not pursue their stability analysis . moreover , the adiabatic dynamics of these structures was studied analytically via a variational approach in ref . @xcite without providing , however , relevant numerical results . nevertheless , it is important to note that in an experimental realization of vortex ring states , the transient observation of spherical shell solitons was reported  @xcite . in addition , they were also suggested as a potential outcome of collisions of vortex rings in the work of ref . @xcite . it is the purpose of the present manuscript to revisit these states and explore not only their existence in isotropic traps , but also perform their systematic bogolyubov - de gennes ( bdg ) stability analysis . this is done by projecting perturbations onto eigenfunctions of the linear problem in the angular directions . using this basis , the relevant computations are greatly simplified , and the full 3d spectrum computation is reduced to an effectively small set of one - dimensional problems ( for the different modes ) . it should be mentioned here that the use of spherical harmonics as a basis in the angular variables for states bearing radial symmetry is rather common in terms of numerical schemes for linear and nonlinear schrdinger ( nls ) equations see , e.g. , ref . @xcite for a a recent example . earlier efforts along these lines in the context of linear schrdinger equations can be found , e.g. , in the works of refs . @xcite . in 2d becs , analogous decompositions of azimuthal modes for radial states including rdss and vortices  @xcite have appeared , and in 3d similar possibilities are starting to emerge  @xcite .    in this work , we find that indeed spherical shell solitons can be identified as stable 3d solutions of the gross - pitaevskii equation , sufficiently close to the linear ( small - amplitude ) limit of small chemical potential from which they bifurcate , as is shown below . moreover , we provide a theoretical analysis of their dynamics as `` effective particles '' . this approach allows to characterize the opposite limit of large chemical potential . our numerical computations of existence and stability , complemented by direct numerical simulations of the spherical shell soliton dynamics , not only corroborate the above two limits , but also provide a systematic way to interpolate between the two analytically tractable regimes . our presentation is organized as follows . first , in sec . [ setup ] , we introduce the model and describe the implementation of the linear stability analysis . our analytical approaches for the spherical shell dark solitons in the above mentioned analytically tractable limits are then presented in sec . [ sec : sub : part ] . next , in sec . [ results ] , we present our systematic numerical results . finally , our conclusions and a number of open problems for future consideration are given in sec . [ conclusion ] . in the appendix , a similar numerical decomposition method but using cylindrical ( rather than spherical ) symmetry is also discussed whose results are also given for comparison in sec .  [ results ] . in the framework of lowest - order mean - field theory , and for sufficiently low - temperatures , the dynamics of a 3d repulsive bec , confined in a time - independent trap @xmath0 , is described by the following dimensionless gross - pitaevskii equation ( gpe )  @xcite : @xmath1 where @xmath2 is the macroscopic wavefunction of the bec and @xmath3 is the chemical potential ( subscripts denote partial derivatives ) . here , we consider a harmonic trap of the form : @xmath4 where @xmath5 , @xmath6 and @xmath7 are the trapping frequencies along the @xmath8 plane and the vertical direction @xmath9 , respectively . note that the potential has rotational symmetry with respect to the @xmath9-axis . in our numerical simulations , we focus on the fully symmetric ( isotropic , spherically symmetric ) case with @xmath10 , in which we can benchmark our numerical methods in effective 1d and 2d against those of the fully 3d setting . we examine the dark spherical shell solitons , hereafter referred to as dss . this single radial node state exists in this isotropic case from the linear limit onwards as a stationary state of the form @xmath11 where @xmath12 is the spherical radial variable . at the linear limit , the waveform is an eigenmode of the quantum harmonic oscillator with chemical potential @xmath13 , and spatial profile : @xmath14 where the basis @xmath15 denotes the 3d harmonic oscillator quantum states in cartesian coordinates .    the dss is generally expected to be unstable for high chemical potentials due to transverse modulational ( snaking - type ) instability , in a way similar to its planar and rds counterparts , as summarized , e.g. , in refs . @xcite . however , for smaller chemical potentials , and especially near the linear limit , it is relevant to analyze stability of the dss and identify corresponding instabilities and where they may arise . in what follows ( cf . [ results ] below ) , we will turn to numerical computations , based on a fixed point iteration and our spherical harmonic decomposition of the spectral stability bdg problem , in order to efficiently identify the intervals of existence and stability of the relevant modes .      in this section , we discuss our analysis of the linear stability of the dss states utilizing their effective 1d nature , namely their spherical symmetry . in spherical coordinates @xmath16 , the laplacian can be decomposed into radial and angular parts , denoted by @xmath17 and @xmath18 respectively , namely @xmath19 where @xmath20 the @xmath18 operator has eigenstates given by the spherical harmonics @xmath21 with eigenvalues @xmath22 , i.e. , @xmath23 for states with spherical symmetry ( @xmath24 ) , the @xmath18 part of the laplacian is not relevant for identifying the stationary state . thus , the stationary dss state @xmath25 of eq . ( [ gpe ] ) , defined in the domain @xmath26 , satisfies the following radial equation : @xmath27 now , let us consider the linear stability ( bdg ) spectrum of such a stationary state . we consider small - amplitude perturbations expanded using the complete basis of @xmath21 ( the spherical harmonic eigen - basis ) as follows : @xmath28 , \label{perturbation}\ ] ] where the space- and time - dependent coefficients @xmath29 and @xmath30 determine the radial and temporal evolution of the perturbation . then , substituting eq . ( [ perturbation ] ) into eq . ( [ gpe ] ) and retaining up to linear terms in the expansion , one notes that all the @xmath31 modes are mutually independent due to the spherical symmetry of the steady state . for mode @xmath31 , the coefficients @xmath32 and @xmath33 of the expansion obey the following evolution equations : @xmath34 -i b_t&=&-\frac{1}{2 } \delta_r b+ \frac{\ell(\ell+1)}{2r^2}b + vb + 2 |\psi_0|^2b + \psi_0^{*2}a-\mu b , \end{array}\ ] ] where we have dropped for simplicity the subscripts @xmath35 . as we are interested in the stability of the steady - state , we now expand the coefficients using @xmath36 where @xmath37 is an eigenvalue of the following matrix @xmath38 with @xmath39 m_{12}&=&-i\psi_0 ^ 2,\\[1.0ex ] m_{21}&=&i\psi_0^{*2},\\[1.0ex ] m_{22}&=&i\left(-\frac{1}{2 } \delta_r + \frac{\ell(\ell+1)}{2r^2 }   + v   + 2|\psi_0|^2-\mu\right ) . \end{array}\ ] ] one can , therefore , compute the full spectrum by computing that of each mode @xmath31 independently , and then putting them all together to `` reconstruct '' the full spectrum . in what follows we use @xmath40 for values of the chemical potential up to @xmath41 . the reduction of dimensionality with the decomposition on a complete basis set is also available in 2d for states with rotational symmetry up to a topological charge along the @xmath9 axis . for clarity , this is shown in the appendix . it is worth mentioning that a large class of bec states belongs to this class . finally , we discuss the application of the degenerate perturbation method ( dpm )  @xcite to the dss state . this method is employed for the determination of the dss spectrum near the linear limit , and also for the identification of the nature of instabilities . the method can be most readily understood from the 3d version of the matrix @xmath42 , where one can separate the free ( linear ) part with known eigenstates and eigenvalues and the rest of the terms can be treated as perturbations near the linear limit . note that there are 2@xmath432 block matrices within @xmath42 and the complete basis consists of both `` positive '' and `` negative '' eigenstates . this terminology is associated with quantum harmonic oscillator eigenstates leading to positive or negative eigenvalues for the linear analogs of the operators @xmath44 in eq . ( [ m ] ) i.e. , with @xmath45 set to @xmath46 . in dpm , the functional dependence of the eigenvalues on the chemical potential @xmath3 is dominated by the interplay of equi - energetic , degenerate states . the result of the nonlinear perturbations ( @xmath47 ) in @xmath48 within eq . ( [ m ] ) is to cause the degenerate eigenvalues to depart from their respective linear limit , as the norm of the solution ( mass ) increases . the free part provides the basis that is formed from `` up '' states and `` down '' states , of the form @xmath49 and @xmath50 with eigenvalues @xmath51 and @xmath52 respectively , where @xmath53 is the eigenenergy of the 3d quantum harmonic oscillator for eigenstate @xmath54 . the relevant states for eigenvalues near the spectrum of im(@xmath37)=0 , 1 and 2 are listed in table . then , the effect of the nonlinearity on these eigenmodes and the corresponding eigenvalue corrections are calculated , as discussed in ref . @xcite ( and originally spearheaded in this context in ref . below , this is compared to the numerical spectrum in sec . [ results ] .    . relevant states for eigenvalues near the spectrum of im(@xmath37)=0,1 and 2 for the degenerate perturbation theory . here @xmath55 stands for eigenstates with all possible distinct permutations with the quantum numbers @xmath56 and @xmath57 . [ para ] [ cols="^,^,^",options="header " , ] we now turn to the limit of large chemical potential @xmath3 . in this so - called thomas - fermi ( tf ) limit of large @xmath3 , the ground state of the gpe is approximated as @xmath58 @xcite . a natural way to obtain a reduced dynamical description of the dss analogously to what is done in lower - dimensional settings is to develop an effective particle picture for the dss in the tf limit . here , following ref . @xcite , we focus on the equilibrium radius @xmath59 of the dss , as a prototypical diagnostic that we can compare to our numerical results . this equilibrium radius is determined as the critical value of the dss radius for which the restoring force due to the harmonic trap is counterbalanced by the effective force exerted due to the curvature of the dss . ( left panels ) and argument @xmath60 ( right panels ) of a dark soliton shell state at @xmath61 ; the top row illustrates the @xmath62 plane , while the middle row the @xmath63 plane . the bottom panel shows an isocontour plot of the same state . , title="fig:",width=302 ]   ( left panels ) and argument @xmath60 ( right panels ) of a dark soliton shell state at @xmath61 ; the top row illustrates the @xmath62 plane , while the middle row the @xmath63 plane . the bottom panel shows an isocontour plot of the same state . , title="fig:",width=226 ]    . bottom panel : number of particles as a function of chemical potential @xmath3 , showing the continuation of states from the linear limit to the nonlinear regime . , width=302 ]    the first approach ( motivated also by the work of ref . @xcite ) focuses on the stationary state using the ansatz @xmath64 , where one obtains an equation for the dss profile @xmath65 ( on top of the tf ground - state @xmath66 ) given by : @xmath67 where @xmath68 and primes denote derivatives with respect to @xmath69 . inspired by the profile of the 1d and 2d equivalents of the dss , namely the dark soliton and the rds , we seek a stationary dss solution in the form of @xmath70 $ ] . multiplying both sides of eq . ( [ neweq ] ) by @xmath71 and integrating over @xmath69 from @xmath72 to @xmath73 ( bearing in mind that the contribution of the integral from @xmath74 to @xmath75 is exponentially small ) , we find that the equilibrium radius is given by @xmath76 where @xmath77 . this result , as in the case of the dark soliton ring in 2d  @xcite , slightly overestimates the actual value for @xmath78 that , as we will see below , is accurately predicted by the second particle picture approach , which retrieves the precise value for @xmath78 to be @xmath79 ; this is confirmed by the numerical results shown below , in sec . [ results ] . the second approach relies on energy conservation and it is based on the analysis of ref . @xcite . in this approach , it is argued that the equation of motion can be derived by a local conservation law ( i.e. , an adiabatic invariant ) in the form of the energy of a dark soliton under the effect of curvature and of the density variation associated with it . more specifically , knowing that the energy of the 1d dark soliton centered at @xmath80 is given by @xmath81  @xcite , the generalization of the relevant adiabatic invariant quantity in a 3d domain bearing density modulations reads : @xmath82 \\[1.0ex ]          & = & 4 \pi r_0 ^ 2 \left [ \frac{4}{3 } ( \mu - v(r_0)-\dot{r}_0 ^ 2)^{3/2 } \right ] , \end{array } \label{ppnew1}\ ] ] where @xmath83 and @xmath84 are the ( radial ) location and velocity of the dss with initial conditions @xmath85 and @xmath86 . taking a time derivative on both sides of eq . ( [ ppnew1 ] ) and assuming that the dss has no initial speed ( @xmath87 ) , we obtain a newtonian particle equation of motion for the dss of the following form : @xmath88 . \label{ppextra}\ ] ] from this equation , the equilibrium position for the dss yields @xmath89 , a result consistent with the numerical findings of the next section . we now proceed to test these predictions , the stability analysis spectrum , and the equilibrium positions for the dss .    . ( a)-(b ) spectra computed using our effectively 1d method ( in the radial variable , using a spherical harmonic decomposition in the transverse directions ) . identical results are obtained when using the 2d variant of the method in cylindrical coordinates ( and decomposing only the azimuthal direction in suitable fourier modes ) . both 1d and 2d numerics use a spacing of @xmath90 . ( c ) spectra computed using the full 3d numerics with spacing @xmath91 . orange ( grey ) lines depict the real part of @xmath37 ( i.e. , unstable parts of the eigenvalues ) , while blue ( dark ) lines the imaginary part of @xmath37 . , title="fig:",width=264 ] . ( a)-(b ) spectra computed using our effectively 1d method ( in the radial variable , using a spherical harmonic decomposition in the transverse directions ) . identical results are obtained when using the 2d variant of the method in cylindrical coordinates ( and decomposing only the azimuthal direction in suitable fourier modes ) . both 1d and 2d numerics use a spacing of @xmath90 . ( c ) spectra computed using the full 3d numerics with spacing @xmath91 . orange ( grey ) lines depict the real part of @xmath37 ( i.e. , unstable parts of the eigenvalues ) , while blue ( dark ) lines the imaginary part of @xmath37 . , title="fig:",width=264 ] . ( a)-(b ) spectra computed using our effectively 1d method ( in the radial variable , using a spherical harmonic decomposition in the transverse directions ) . identical results are obtained when using the 2d variant of the method in cylindrical coordinates ( and decomposing only the azimuthal direction in suitable fourier modes ) . both 1d and 2d numerics use a spacing of @xmath90 . ( c ) spectra computed using the full 3d numerics with spacing @xmath91 . orange ( grey ) lines depict the real part of @xmath37 ( i.e. , unstable parts of the eigenvalues ) , while blue ( dark ) lines the imaginary part of @xmath37 . , title="fig:",width=264 ] let us start by providing some of the basic features of the dss . a typical dss state at @xmath61 is shown in fig . [ state ] , where its spherical symmetry is apparent . a radial plot of the wavefunction in the tf limit is depicted in the top panel of fig . [ lc ] . as the chemical potential is increased the number of atoms ( `` mass '' ) @xmath92 increases as depicted in the bottom panel of fig . [ lc ] . we now examine the stability of the dss . the corresponding spectra illustrating the imaginary part of the eigenvalues [ dark ( blue ) lines ) and the real part [ light ( orange ) lines ] are shown in fig . [ f : spectrum ] . panel ( a ) depicts the spectrum obtained via our numerical method over a wide range of chemical potentials . identical results were obtained when using the method in 2d , based on a cylindrical coordinate decomposition and a representation of the azimuthal variable dependence in the corresponding fourier modes . for reasons of completeness , a summary of this variant of the method is presented in the appendix . panels ( b ) and ( c ) depict a zoomed in region to contrast the results between these lower dimensional methods and the full 3d numerics . panel ( b ) corresponds to the spectrum as computed by means of the 1d spherical coordinate decomposition method using spherical harmonics with spatial spacing of @xmath90 , while panel ( c ) depicts the results from direct full 3d numerics using an spatial spacing of @xmath91 . a close comparison between the panels suggests small discrepancies of the imaginary parts at im@xmath93 and im@xmath94 , and a spurious unstable mode with re@xmath95 . these discrepancies ( and their amendment as the mesh resolution @xmath96 decreases see below ) are likely due to the large spatial spacing @xmath96 that needs to be chosen such that the 3d calculations are still `` manageable'' in which case the dss radial profile is under - resolved . it is crucial to note that the 3d calculations are far more expensive than their lower dimensional ( 1d or 2d ) counterparts considered above .    . the larger frequency corresponds to the oscillatory mode of the dss while the smaller one corresponds to the breathing mode of the background density . , width=264 ]    to confirm that the discrepancies are indeed caused by the mesh resolution , when we further decrease the lattice spacing to @xmath97 ( results not shown here ) and recalculate the spectrum in 3d , the branches at im@xmath93 indeed split into two non - degenerate curves and the spurious modes at im@xmath94 and re@xmath98 begin to get suppressed . this suggests that the discrepancies are indeed attributable to resolution effects and that the spectrum computed with the effectively 1d numerical method is indeed accurate .    ) from the full 3d numerical solution for @xmath99 . the initial condition corresponds to the stationary dss at its equilibrium radius . second set of panels : evolution of a dss shifted from its equilibrium position undergoing oscillatory motion over one oscillation period @xmath100 . the white dashed circle denotes the starting location of the dss at the beginning of the period . third panel : radial profiles for the oscillating dss over half a period . bottom panel : oscillations of the dss radius vs.  time . blue dots represent the radius extracted from full 3d numerics while the solid red line is the best ( least - square ) fit to a linear combination of oscillations with frequencies @xmath101 and @xmath102 ( see red dots in fig . [ dpm ] ) corresponding , respectively , to the frequencies of the dss oscillations and the breathing mode of the background cloud . for a movie depicting the stable oscillations of the dss please see the supplemental material . , title="fig:",width=257 ] ) from the full 3d numerical solution for @xmath99 . the initial condition corresponds to the stationary dss at its equilibrium radius . second set of panels : evolution of a dss shifted from its equilibrium position undergoing oscillatory motion over one oscillation period @xmath100 . the white dashed circle denotes the starting location of the dss at the beginning of the period . third panel : radial profiles for the oscillating dss over half a period . bottom panel : oscillations of the dss radius vs.  time . blue dots represent the radius extracted from full 3d numerics while the solid red line is the best ( least - square ) fit to a linear combination of oscillations with frequencies @xmath101 and @xmath102 ( see red dots in fig . [ dpm ] ) corresponding , respectively , to the frequencies of the dss oscillations and the breathing mode of the background cloud . for a movie depicting the stable oscillations of the dss please see the supplemental material . , title="fig:",width=257 ] ) from the full 3d numerical solution for @xmath99 . the initial condition corresponds to the stationary dss at its equilibrium radius . second set of panels : evolution of a dss shifted from its equilibrium position undergoing oscillatory motion over one oscillation period @xmath100 . the white dashed circle denotes the starting location of the dss at the beginning of the period . third panel : radial profiles for the oscillating dss over half a period . bottom panel : oscillations of the dss radius vs.  time . blue dots represent the radius extracted from full 3d numerics while the solid red line is the best ( least - square ) fit to a linear combination of oscillations with frequencies @xmath101 and @xmath102 ( see red dots in fig . [ dpm ] ) corresponding , respectively , to the frequencies of the dss oscillations and the breathing mode of the background cloud . for a movie depicting the stable oscillations of the dss please see the supplemental material . , title="fig:",width=257 ] ) from the full 3d numerical solution for @xmath99 . the initial condition corresponds to the stationary dss at its equilibrium radius . second set of panels : evolution of a dss shifted from its equilibrium position undergoing oscillatory motion over one oscillation period @xmath100 . the white dashed circle denotes the starting location of the dss at the beginning of the period . third panel : radial profiles for the oscillating dss over half a period . bottom panel : oscillations of the dss radius vs.  time . blue dots represent the radius extracted from full 3d numerics while the solid red line is the best ( least - square ) fit to a linear combination of oscillations with frequencies @xmath101 and @xmath102 ( see red dots in fig . [ dpm ] ) corresponding , respectively , to the frequencies of the dss oscillations and the breathing mode of the background cloud . for a movie depicting the stable oscillations of the dss please see the supplemental material . , title="fig:",width=257 ]    the comparison between the stability spectrum predicted by the degenerate perturbation method ( dpm ) and the one obtained numerically near the linear limit is shown in fig . as the figure shows , the dpm captures the behavior of the spectrum close to the linear limit . let us now confirm the stability of the dss in this limit . we have performed a direct numerical integration of the full 3d dynamics of the dss for @xmath99 , where it is predicted to be stable cf . [ f : spectrum ] and [ dpm ] . this evolution is depicted in the top set of panels of fig . also , in this figure , we depict the stable oscillations of a perturbed dss by displacing its initial location away from the equilibrium radius . the evolution for such a dss performing oscillations is displayed in the second set of panels in fig . [ rk4s ] . a movie for this evolution is provided in the supplemental material . to more clearly visualize these oscillations , we depict on the third panel of fig . [ rk4s ] the transverse density cuts at four different stages over half a period of the oscillation . these density cuts also reveal the excitation of the background cloud which performs breathing oscillations . from the density cuts we extracted the ( radial ) location of the dss along the evolution which displays a beating - type quasiperiodic oscillation as depicted by the blue dots in the bottom panel of fig . in fact , by ( least - square ) fitting a linear combination of two harmonic oscillations ( see red line ) we extract the frequencies @xmath101 and @xmath102 . upon closer inspection , these frequencies correspond , respectively , to the frequency of the dss oscillations and the frequency of the breathing mode of the background cloud . these frequencies are also depicted in fig . [ dpm ] clearly showing that they indeed belong to the stability spectrum of the stationary dss state . it is important to stress that the observed dss oscillations are supported by the stability of the dss at its equilibrium position . therefore , despite strong perturbations of the background that undergoes breathing oscillations , the dss robustly persists and performs stable oscillations about its equilibrium position . as the spectrum indicates , see figs . [ f : spectrum ] and [ dpm ] , for larger values of the chemical potential the dss becomes unstable . to identify this instability , we have looked at the eigenvectors of the dpm . there are high degeneracies in the eigenvalues in this case , rendering harder the identification of the nature of the instability . nevertheless , the dpm is still helpful in that regard . for instance , we know from the eigenvectors that the mode with the largest slope decreasing at im@xmath93 , which causes instabilities , is a linear combination of `` up '' states , and the `` down '' states are not involved . therefore , the instability is due to a bifurcation rather than a collision with `` negative '' energy modes ; the latter may lead to oscillatory instabilities , as discussed , e.g. , in ref . the former may be associated with symmetry breaking features and exhibits the instability via real eigenvalue pairs . direct numerical integration of eq . ( [ gpe ] ) for chemical potentials @xmath103 and 10.6 , shows that the first instability is a bifurcation of a six vortex lines ( vl6 ) cage , while the rest of the instabilities are dominated by a cubic six vortex rings ( vr6 ) cage . the first instability suggests a bifurcation of vl6 from the interactions of the dss and a dark soliton state with three plane nodes perpendicular to the @xmath8-plane , separated by angles of @xmath104 ; the dark soliton state ( ds3 ) can be written as @xmath105 and a linear eigenenergy @xmath106 . a two - mode analysis  @xcite using this state and the dss state yields a prediction of the bifurcation as occurring at @xmath107 , which is in good agreement with the numerical result @xmath108 of fig . [ f : spectrum ] . a similar treatment of the vr6 cage turns out to be more challenging due to the difficulty in designing the corresponding dark soliton state . instead , we have looked at the ds4 state , which can form an eight vortex line ( vl8 ) cage . the two - mode analysis then predicts a bifurcation at @xmath109 , which is in fair agreement with the second instability at @xmath110 from fig . [ f : spectrum ] . the numerical results , however , show that the dynamically resulting state is a vr6 cage , not a vl8 cage . presumably the vl8 cage can also cause instabilities of the dss but is less robust than the vr6 cage , possibly due to the intersections of the vortex lines , which are absent in vr6 . further analysis with higher values of @xmath111 ( not shown here ) also suggests the progressively lesser relevance of the dsm state in causing instabilities . on the other hand , the vr6 state seems to be very robust and dominates the instability in a wide range of chemical potentials from the onset of this instability around @xmath112 . the vl6 and vr6 cages are shown in fig . [ vl6vr6 ] .    ) dynamics at @xmath113 and @xmath114 , respectively . in particular , the vr6 cage seems to be robust and appears in the dynamics for a wide range of chemical potentials from the onset of the instability around @xmath112 . , title="fig:",width=158 ] ) dynamics at @xmath113 and @xmath114 , respectively . in particular , the vr6 cage seems to be robust and appears in the dynamics for a wide range of chemical potentials from the onset of the instability around @xmath112 . , title="fig:",width=158 ]    finally , we present the numerical results for the equilibrium location @xmath59 of the dss as a function of the chemical potential @xmath3 and compare the results with those of the two particle pictures in the large density limit ( see sec .  [ sec : sub : part ] ) . a plot of @xmath115 as a function of  @xmath3 is depicted in fig . [ rc ] . note that the first particle picture predicts a slightly larger equilibrium @xmath59 , while the particle picture based on energy conservation of the dss as an adiabatic invariant agrees very well with the numerical results . it is interesting to note that the asymptotic behavior only sets in at about @xmath116 , deep inside the tf regime . although this would be a rather computationally demanding parametric range when using 3d ( or even 2d ) methods , our effective 1d approach based on a spherical harmonic decomposition of the angular dependencies , allows us to reach such large values of @xmath3 . , as a function of @xmath3 ( solid red line ) . note that the numerical values reach the asymptotic value of @xmath117 ( dashed black line ) , as per the second particle picture , when @xmath3 is large . the first particle picture slightly overestimates @xmath59 ( dash - dotted blue line ) . , width=302 ] in this work , we have revisited the theme of spherical shell dark solitons . we have employed analytical methods to study these structures both in the vicinity of the linear limit , using a bifurcation analysis , and in the thomas - fermi limit , where they can be treated as effective particles . for the numerical simulations , we have used continuation methods , as well as a quasi - one - dimensional radial computation method , decomposing the angular dependence of the perturbations in spherical harmonics , to determine their stability . we have found that spherical shell dark solitons can , in fact , be dynamically robust , spectrally stable solutions of the 3d gross - pitaevskii equation within a certain interval of chemical potentials ( and atom numbers ) ; this occurs sufficiently close to the linear limit of relatively small chemical potentials . their instabilities , emerging when the chemical potential increases , have been elucidated and their dynamical outcome , namely the breakup into vortex ring and vortex line cages , has been showcased . we believe that these results , may not only be used to partially explain the transient observation of these structures in earlier experiments  @xcite and numerical simulations  @xcite , but may also pave the way for their consideration in future experiments . our work suggests a number of interesting future research directions , concerning , in particular , related coherent structures , their physical relevance and properties . for instance , this study suggests large scale stability computations for a wide range of pertinent steady states with radial symmetry , including ones with different numbers of radial nodes as well as a number of different external potentials . furthermore , it would be relevant to generalize and apply the numerical ( and analytical ) techniques used herein to multicomponent bec settings . for example , one can study `` symbiotic states '' composed of spherical shell dark - bright solitons , in the spirit of the dark - bright 2d rings of ref . @xcite . on the other hand , states emerging robustly from the instabilities of these spherically ( or cylindrically ) symmetric ones , such as the vortex ring ( and line ) cages are worthwhile to investigate further in their own right . efforts along these directions are currently in progress and will be presented in future publications . w.w .  acknowledges support from nsf ( grant nos . dmr-1151387 and dmr-1208046 ) . p.g.k .  gratefully acknowledges the support of nsf - dms-1312856 , as well as from the us - afosr under grant fa950 - 12 - 1 - 0332 , and the erc under fp7 , marie curie actions , people , international research staff exchange scheme ( irses-605096 ) . r.c.g .  gratefully acknowledges the support of nsf - dms-1309035 . we thank the texas a&m university for access to their ada cluster . the greatest advantage of the 1d method involving the decomposition into spherical harmonics is that one can compute the spectrum for large chemical potentials without excessive computational costs . the reason that typical computations for large chemical potentials become computationally prohibitive , is that in this limit the density is large and , thus , the width of the relevant localized structures ( proportional to @xmath118 ) becomes much smaller than the domain size ( proportional to @xmath119 ) ; hence , a large number of mesh points is necessary to resolve properly these configurations . nonetheless , it is worth mentioning that , for the projection method in 1d , the form of states that one can study is also restricted by the spherical symmetry ; therefore , it is also relevant to employ an analog of this method in 2d , so as to encompass a wider class of states . an example concerns states with a topological charge @xmath120 along the @xmath9 axis , i.e. , stationary states of the form @xmath121 . note that in this framework , one can study states including but not limited to the ground - state , planar , ring , or spherical shell dark soliton states , vortex lines , as well as coplanar or parallel vortex rings and hopfions  @xcite , in both isotropic and anisotropic traps . in 2d , and for cylindrical coordinates @xmath122 , the laplacian is decomposed as follows : @xmath123 where its components read : @xmath124 the @xmath125 operator has eigenstates @xmath126 with eigenvalues @xmath127 , i.e. , @xmath128 a stationary state @xmath129 , defined in the domain @xmath130 , and bearing a topological charge @xmath120 , satisfies the following equation : @xmath131 as in the 1d case , we construct the linear stability problem as follows . let @xmath132 a rotationally symmetric stationary state , up to a topological charge @xmath120 , along the @xmath9 axis perturbed using the complete basis of @xmath126 : @xmath133 \right ] . % \label{perturbation2 } \nonumber\ ] ] substituting this expansion into the gpe of eq . ( [ gpe ] ) , and upon linearizing and matching the basis expansion on both sides , one can see that the different @xmath111 modes are mutually independent . an equivalent derivation as in 1d shows that @xmath37 are the eigenvalues of the matrix @xmath134 where @xmath135 m_{12}&=&-i\psi_0 ^ 2,\\[1.0ex ] m_{21}&=&i\psi_0^{*2},\\[1.0ex ] m_{22}&=&i\left(-\frac{1}{2 } \delta_h + \frac{(m - s)^2}{2\rho^2 }   + v   + 2|\psi_0|^2-\mu\right ) . \end{array}\ ] ] as before , one can therefore compute the full spectrum by computing each mode @xmath111 independently and then putting them together . in the results based on this method and presented in sec . [ results ] , we use @xmath136 . it is worth commenting here that the matrix for the 2d eigenvalue problem with a charge @xmath120 is less symmetric than matrices corresponding to the 1d projection and the 2d case with no charge . therefore , the full 2d problem with charge appears to demand more computational work . however , it should be noted that , in the 2d case , the eigenvalues of the full matrix @xmath42 corresponding to the projections over the @xmath137 and @xmath138 modes are related via an orthogonal transformation and , thus , the two set of eigenvalues are complex conjugates of each other . therefore , it is only necessary to compute the non - negative @xmath111 set and , hence , the charged and the uncharged states actually have similar computational complexity . p. g. kevrekidis , d. j. frantzeskakis , and r. carretero - gonzlez ( eds . ) , _ emergent nonlinear phenomena in bose - einstein condensates . theory and experiment _ ( springer - verlag , berlin , 2008 ) ; r. carretero - gonzlez , d. j. frantzeskakis , and p. g. kevrekidis , nonlinearity * 21 * , r139 ( 2008 ) . j. denschlag , j.e . simsarian , d.l . feder , c.w . clark , l.a . collins , j. cubizolles , l. deng , e.w . hagley , k. helmerson , w.p . reinhardt , s.l . rolston , b.i . schneider , and w.d . phillips , science * 287 * ( 2000 ) 97101 . b. p. anderson , p. c. haljan , c. a. regal , d. l. feder , l. a. collins , c. w. clark , and e. a. cornell , phys .  rev . lett .  * 86 * , 2926 ( 2001 ) . p. engels and c. atherton , phys . lett .  * 99 * , 160405 ( 2007 ) . i. shomroni , e. lahoud , s. levy and j. steinhauer , nature phys . * 5 * , 193 ( 2009 ) . s. middelkamp , p. g. kevrekidis , d. j. frantzeskakis , r. carretero - gonzlez , and p. schmelcher , physica d * 240 * , 1449 ( 2011 ) . w. wang , p. g. kevrekidis , r. carretero - gonzlez , d. j. frantzeskakis , t. j. kaper , and m. ma , phys . a * 92 * , 033611 ( 2015 ) .
in this work we study spherical shell dark soliton states in three - dimensional atomic bose - einstein condensates . 
 their symmetry is exploited in order to analyze their existence , as well as that of topologically charged variants of the structures , and , importantly , to identify their linear stability bogolyubov - de gennes spectrum . 
 we compare our effective 1d spherical and 2d cylindrical computations with the full 3d numerics . 
 an important conclusion is that such spherical shell solitons can be stable sufficiently close to the linear limit of the isotropic condensates considered herein . 
 we have also identified their instabilities leading to the emergence of vortex line and vortex ring cages . 
 in addition , we generalize effective particle pictures of lower dimensional dark solitons and ring dark solitons to the spherical shell solitons concerning their equilibrium radius and effective dynamics around it . in this case too , we favorably compare the resulting predictions such as the shell equilibrium radius , qualitatively and quantitatively , with full numerical solutions in 3d .