\section{Introduction}
Dark solitonic structures emerge ubiquitously in systems that combine dispersion with a self-defocusing nonlinearity.
Numerous physical examples encompassing these properties prominently include atomic Bose-Einstein condensates (BECs)~\cite{kevrekidis2015defocusing,Frantzeskakis_2010}, and nonlinear Kerr (and generalized non-Kerr) media~\cite{KIVSHAR199881}. 
More broadly, solitonic structures also arise in mechanical lattices of coupled pendula~\cite{PhysRevLett.68.1730}, electrical transmission lines~\cite{PhysRevE.49.828}, microwave magnetic thin films~\cite{PhysRevLett.70.1707,PhysRevLett.84.4697}, acoustic waveguides of Helmholtz resonators~\cite{PhysRevE.91.023204}, surface water waves~\cite{PhysRevLett.110.124101}, and nematic liquid crystals~\cite{Piccardi:11}.
Outside of these closed systems, solitonic waves are even present in media balancing gain and loss, such as polariton fluids in semiconductor microcavities~\cite{carusotto_2013}. The presence of solitary waves in such disparate contexts attests to the breadth of their relevance and the generality of the central mechanisms leading to their creation. In particular, atomic BECs hold a central role amongst these systems owing to their highly tunable interaction Hamiltonians, the control of dimensionality and potential landscapes, as well as the diverse available mechanisms for creating excitations~\cite{pethick,stringari,kevrekidis2015defocusing}.

There are different ways BECs are used to create solitonic states. For example, dragging a laser beam through a BEC~\cite{PhysRevLett.99.160405}, phase-imprinting~\cite{Becker2008}, and matter-wave interference~\cite{PhysRevLett.101.130401,PhysRevA.81.063604}. More recent experiments showed control of both the position and velocity of soliton~\cite{Fritsch2020creating}, with contemporary theoretical studies exploring 3-dimensional (3d) solitary waves in detail~\cite{mateo2015stability}. The intriguing dynamics of the resulting structures, including their breakup into vortical patterns~\cite{PhysRevLett.86.2926,Shomroni2009} (even in superfluid Fermi gases~\cite{PhysRevLett.113.065301}) is a subject of wide investigation and has a wide variety of potential applications ranging from atomic matter-wave interferometers~\cite{Negretti_2004} to producing two-level qubit systems~\cite{PhysRevA.95.053618}. Much of this experimental and theoretical effort is summarized in Refs.~\cite{Frantzeskakis_2010,kevrekidis2015defocusing}.

In the present work, we revisit dark solitary waves 3d repulsively interacting BECs highly elongated along the longitudinal axes $\ez$. These waves are associated with a reduction in the local atomic density accompanied with a modification of the condensate phase. In elongated systems, vortex rings have a similar density distribution but with a ring-shaped phase singularity~\cite{komineas,kevrekidis2015defocusing}. The simplest example is a planar (or kink) soliton where the density is reduced to be two dimensional on a plane normal to the soliton's direction of propagation. By contrast, a solitonic vortex is a highly anisotropic vortex with a linear phase-singularity~\cite{kevrekidis2015defocusing}.  
The potential conversion of kink solitons into these excitations is of wide interest~\cite{PhysRevA.68.043617,PhysRevA.65.043612,Becker_2013,tylutki2015solitonic,PhysRevA.94.023627,PhysRevLett.113.065302}.

We begin by comparing experimentally observed solitary waves with 3d numerical simulations computed with no adjustable parameters. Motivated
by this comparison, we theoretically explore the related 3d case of a longitudinally traveling solitonic wave in a system that is transversely confined but longitudinally infinite. Then, we compute the relevant Bogolyubov-de Gennes (BdG) modes to identify the onset of transverse instability and to quantify the growth rate of the unstable mode. This process confirms the prediction of Ref.~\cite{PhysRevLett.89.110401} that the interval of stability of kink solitons increases as their speed approaches the local speed of sound c. Finally in the unstable regime, we numerically follow the exponential growth of the unstable mode to its ultimate conversion to vortex structures.

\begin{figure*}[t]
\includegraphics{Fig_PHE2Din0}
	\caption{Integrated GPE wavefunction density (top) and phase cross-section (bottom) evaluated 

	immediately after a $t_p = 140\ \mu\mathrm{s}$ phase imprint.
	(a) and (b) plot  $|\psi(x,z)|^2$ and $|\psi(y,z)|^2$ respectively.}
    \label{fig:PHE2D0}
\end{figure*}

\section{Experimental Motivation: Phase Engineering and In-Trap Motion}

In the present work, we model the experimental procedure in reference~\cite{Fritsch2020creating} which describes a new wavefunction engineering approach for launching solitons in atomic BECs and show that dark planar solitons are first created and then relax into a solitonic vortex.

The experiments started with $N=2.4\times10^5$ $^{87}\text{Rb}$ BECs confined in an optical dipole trap with frequencies $(\omega_x,\omega_y,\omega_z)=2\pi \times (94.73,153.23,9.1)\ \text{Hz}$ and chemical potential $\mu=h\times1.1$~kHz. 
We created solitonic excitations by combining
conventional phase imprinting with density engineering, i.e., wavefunction engineering, yielding solitonic excitations with a wide range of initial velocities~\cite{Fritsch2020creating}.

We consider the fully $3d$ Gross-Pitaveskii model (GPE), complementing these experiments~\cite{pethick,stringari}
\begin{equation}
  i\hbar\frac{\partial\psi}{\partial t}= \left[-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{r})+g|\psi|^2 \right] \psi.
  \label{eq:3DGP}
\end{equation}
Here, $\psi$ is normalized to the total atom number $N_{\mathrm{3D}}$,
interaction constant and potential are defined as $g=4 \pi \hbar^2 a_s/m$ and $V(\mathbf{r})=m(\omega_x^2 x^2+\omega_y^2 y^2+\omega_z^2 z^2)/2+V_\mathrm{step}(t) \mathcal{H}(z)$, respectively, where $\mathcal{H}(z)$ represents the Heaviside function.
We fix the magnitude of the step potential $V_\mathrm{step}(t)$ to be $0$ except during the phase-imprinting pulse when it is $h\times5.5\ \mathrm{kHz}$.

We simulate the experiments using the parameters described above and create solitons via the phase imprinting method: applying a potential generated with an infinitely sharp edge (smaller than the BEC healing length) to half of the BEC.
We note that the density engineering method is implemented in the experiment to circumvent technical limitations, given that the BEC healing length is typically much smaller than the optical resolution of the imaging systems used for phase imprinting, which is not the case in numerical simulations.
We found that the idealization of a step potential centered at $z=0$ and amplitude $V_{\mathrm{step}}$, is sufficient to reproduce the experimental results.
This indirectly demonstrates the efficacy of the experimental approach.

\begin{figure*}
     \includegraphics[width=\textwidth,angle =0]{Fig_PHE2Dm}
       \caption {
Integrated GPE wavefunction density (top) along with phase (bottom) cross-sections evaluated $400\ \mathrm{ms}$ after a $t_p = 140\ \mu\mathrm{s}$ phase imprint.
(a) plots the $|\psi|^2$ projection and associated cross-sections in the $\ex-\ez$ plane and (b) plots $\ey-\ez$ quantities.
 }
    \label{fig:PHE2D}
\end{figure*}

\begin{figure}
     \includegraphics[width=\columnwidth,angle =0]{Fig_PHE1Dm}
     \caption{ Oscillation of the solitonic vortices (a) and oscillation frequencies (b).
    The black symbols on the grey colored integrated GPE wavefunction density in $\ez - t$ plane of (a) represent experimental data.
    Every set of experimental parameters was repeated three times, and a symbol is present for each clearly identified solitonic excitation. In cases in which more than one solitonic excitation was found in a single image we chose the position of the deepest depletion. The simulations use $t_p=140\ \mu\mathrm{s}$. 
     The symbols, red squares and black circles, in (b) represent GPE and experimental calculations, respectively, where the error bar represents the single-$\sigma$ uncertainty  obtained from the sinusoidal fits.
     }
    \label{fig:PHE1D}
\end{figure}


Figure~\ref{fig:PHE2D0} shows the 3d density projected onto the $\ex-\ez$ and $\ey-\ez$ planes, i.e., $|\psi(x,z)|^2$ and $|\psi(y,z)|^2$, and the associated planar phase cuts, immediately following the application of $V_\text{step}$ for a phase imprinting time $t_p = 140\ \mu\text{s}$. 
The phase-imprinted planar dark soliton subsequently evolves into a vortex ring as shown in ~\cite{supple}.
This ring eventually decays~\cite{PhysRevA.98.033609} into a pair of solitonic vortices aligned along the direction of stronger transverse confinement and exhibits a large amplitude oscillation ~\cite{supple}. 
The weak interaction of the resulting
structures with the background of excited phonons leads to effective (apparent) dissipative longitudinal dynamics and to the selection of two separate trajectories, one for each solitonic vortex as shown in Figs.~\ref{fig:PHE2D} and~\ref{fig:PHE1D}. 
{In the experimental results, over long times, the amplitude
of oscillation appears to be growing, possibly reflecting the (weak)
anti-damping effect due to the thermal fraction~\cite{prouka}.}
Both the numerical and experimental results (black crosses) find an oscillation frequency $\omega_s$

far from the ``canonical'' result of $\omega_s =  \omega_z / \sqrt{2}$ for planar solitons~\cite{PhysRevLett.84.2298,Frantzeskakis_2010}. 

Instead, the solitonic vortex oscillation frequency is 
significantly reduced~\cite{mateo2015stability,Fritsch2020creating}, e.g.,
 by a factor of about 3 in comparison to $\omega_z$, as shown in the bottom panel of Fig.~\ref{fig:PHE1D}. More generally, in line with~\cite{mateo2015stability}, we find that the relevant reduction is a function of the chemical potential, and the general agreement between the experimental data and the numerical simulations supports 
 solitonic vortex pair formation from the method of phase imprinting.
 
This experiment-theory correspondence prompts us to further explore the stability of planar solitons and their transformation into other structures.

\section{Stability of planar dark solitons}

In this section we study the stability of dark solitary waves using a dimensionless version of Eq.~\eqref{eq:3DGP}
\begin{equation}
  i\frac{\partial\psi'}{\partial t'}= \left[-\frac{1}{2}(\nabla')^2+V'(\mathbf{r}')+g'|\psi'|^2 \right] \psi',
  \label{eq:3DGPd}
\end{equation}
in terms of the dimensionless quantities $t'=\omega_{\perp} t$, $\mathbf{r}'=\mathbf{r}/a_{\perp}$, $\psi'=a_{\perp}^{3/2} \psi$, $V(\mathbf{r}')=[l_x^2 (x')^2+l_y^2 (y')^2+l_z^2 (z')^2]/2$, $g'=4 \pi a_s/a_{\perp}$, and $\mu'=\mu/(\hbar \omega_{\perp})$.
Here, we use the transverse degrees of freedom to define the natural units of frequency $\omega_{\perp}=\sqrt{\omega_x \omega_y}$ and length $a_{\perp}=\sqrt{\hbar/(m \omega_{\perp})}$. 
In addition, we introduce anisotropy parameters $l_x=\omega_x/\omega_{\perp}$, $l_y=\omega_y/\omega_{\perp}$ and $l_z=\omega_z/\omega_{\perp}$.

In line with the last section, we use a chemical potential $\mu= h\times 1.1\ \mathrm{kHz}$, longitudinal frequency $\omega_z/(2\pi)=10\ \mathrm{Hz}$, and transverse frequency ranging from $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$ to $1\ \mathrm{kHz}$. In what follows, we take $\omega_x=\omega_y$ for simplicity.

The 1D counterpart of Eq.~\eqref{eq:3DGP} 
with $V=0$ supports the dark soliton solution~\cite{kevrekidis2015defocusing}
\begin{equation}
 \psi_k'(z',t')=\sqrt{n_0} \left[B \tanh(\zeta)+i A\right]e^{i(kz'-\omega t')}.
  \label{eq:1D_dark}
\end{equation}
Here, $n_0=\mu'$, $\zeta=\sqrt{n_0}B[z'-z_0(t')]$, $z_0=vt'+Z_0$, $v=A\sqrt{n_0}+k$, $\omega=k^2/2+n_0$ and $A^2+B^2=1$, and restrict ourselves to $k=0$ solutions (of a stationary background). The remaining properties have the standard meaning,
i.e., $v$ is the soliton speed, $n_0$ the background density,
$z_0$ marks the center of mass (and $Z_0$ is its initial value), 
while $B$ is associated with the
inverse solitonic width~\cite{kevrekidis2015defocusing}.

\subsection{Stationary solutions}

In the presence of a trap, the wavefunction for a 3d stationary solution is approximately  
\begin{equation}
  \psi'(\mathbf{r'},t')=\psi_{ds}(z')\Phi(\mathbf{r'})e^{-i \mu' t'},
  \label{eq:1D_darks}
\end{equation}
  where $\Phi(\mathbf{r'})$ represents the ground state Thomas-Fermi (TF) wavefunction~\cite{pethick,stringari} and 
  $\psi_{ds}$ represents the stationary dark soliton from
  Eq.~(\ref{eq:1D_dark}), associated with $A=v=0$. 
\begin{figure*}
     \includegraphics{Fig_DSS}
       \caption{Condensate GPE density of the Eq.~\eqref{eq:ste} (top) and phase cross-section (bottom) show the stationary 3D dark soliton solution for $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$, where $\omega_x=\omega_y$ and $v=0$. The lower panel shows the soliton phase jump from  $\pi$ (dark) to 0 (light). 
       (a) and (b) plot $|\psi'(x',0,z')|^2$ and $|\psi'(0,y',z')|^2$ respectively.
       }
    \label{fig:DSS}
\end{figure*}
Using Eq.~\eqref{eq:1D_darks} as an initial guess for a root finding algorithm, we obtain numerically exact solutions by solving
\begin{equation}
   \mathcal{F}(\psi') = \left[-\frac{1}{2}\nabla^2+V(\mathbf{r'})+g'|\psi'|^2 -\mu' \right]\psi'=0.
      \label{eq:ste}
  \end{equation}
Fig.~\ref{fig:DSS} shows the solution of Eq.~\eqref{eq:ste} for $\omega_{\perp}/(2\pi)=100$~Hz. We now 
explore both the stable and unstable modes of $\psi'$ by identifying the eigenvalues associated with the BdG equations~\cite{pethick,stringari,kevrekidis2015defocusing} 
for the parameter range of $\omega_{\perp}$.
For relevant details, see the Appendix \ref{app:ES}. An eigenvalue $\lambda$ has a  real part denoted by $\lambda_r$ and an imaginary part denoted by $\lambda_i$. It is important to recall that a non-vanishing $\lambda_r$, given the  Hamiltonian nature of our GPE (\ref{eq:3DGPd}), is tantamount to the presence of a dynamical instability.
 %
  %
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_DSEigim}
       \caption{Imaginary eigenvalues $\lambda_i$ obtained from the BdG analysis  plotted as a function of $\omega_{\perp}$ for $\omega_x=\omega_y$. 
       The vertical dashed black line represents the mode responsible for instability of the planar
dark soliton, while blue lines represent stable modes. In the limit of large $\omega_{\perp}$, the 1d analytically available results, indicated by the red short lines, are retrieved (see text for further details).
   
     
      }
    \label{fig:DSEigim}
\end{figure}
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_DSEigr}
       \caption{The largest real eigenvalue $\lambda_r$, 
       associated with the
       instability of stationary, planar dark soliton, as a function of $\omega_{\perp}$ for $\omega_x=\omega_y$. The red circles are obtained from GPE dynamics and blue line represents BdG analysis. 
      
     
     
       }
    \label{fig:DSEigr}
\end{figure}

There exist two limiting cases for the stability problem of the dark solitary waves in a trapped condensate, more
specifically, the TF limit at large chemical
potential $\mu'$ and the linear limit at small $\mu'$. 
As is well-known~\cite{kevrekidis2015defocusing} 
in the 1d limit of the GPE model, the anomalous mode ($\omega_0$), corresponding to negative energy is associated with the oscillation of the single dark soliton with a frequency of $\approx l_z /\sqrt{2}$~\cite{PhysRevLett.84.2298}.
The other modes $\omega_n$ (associated with the
vibration modes of the background cloud) are expected to approach the frequencies of $l_z\sqrt{n(n+1)/2}$~\cite{stringari}. On the other hand, in the linear limit of vanishing density $|\psi'|^2 \rightarrow 0$, the eigenvalue problem of Eq.~\eqref{eq:3DGP} reduces to that of a {\it linear} quantum harmonic oscillator with corresponding energies of eigenmodes $|k,l,m\rangle$ 
\begin{align}
E_{k,l,m} &= \left(k+ \frac{1}{2}\right) l_x+ \left(l+ \frac{1}{2}\right) l_y+
\left(m + \frac{1}{2}\right) l_z.
\end{align}
A single dark soliton at $z'=0$ (pertaining to the $|0,0,1\rangle$ in the linear limit) has energy $E_{0,0,1} = 
(l_x + l_y+ 3 l_z)/2$.
To explore the validity of these predictions, we rescale the imaginary eigenvalues by $l_z$ and show the corresponding results in Fig.~\ref{fig:DSEigim}. 
Indeed, our numerically calculated imaginary eigenvalues confirm the expectations on the basis of the 1d
predictions, as our numerical continuation results
approach progressively the 1d TF limit as $\omega_{\perp}$
increases. We note that the chemical potential controls the approach to the TF limit for the present setting. The mode responsible for the instability of the planar
dark soliton appears as a nearly vertical dashed black line in 
Fig.~\ref{fig:DSEigim}.


   \begin{figure*}
 \includegraphics[width=\textwidth,angle =0]{Fig_DSEigv415}
       \caption{(left)  The eigenvector $|\varphi|^2$ corresponds to the largest real eigenvalue of Eq.~\eqref{eq:gpeeig} and its phase profile for the parameter $\omega_{\perp}/(2\pi)=415$~Hz and $\omega_x=\omega_y$. (right panel) The density $|\psi'+\epsilon \varphi|^2$, where $\psi'$ is the steady state wavefunction and $\epsilon=20$. Lower panels in (a) and (b) show the phase cross-section.
     
      }
    \label{fig:DSEigv415}
\end{figure*}
The analysis of real eigenvalues displayed in Fig.~\ref{fig:DSEigr} shows that, as $\omega_{\perp}$  decreases (going away from the above mentioned quasi-1d limit), there exists an instability window for smaller values of the relevant parameter. The critical $\omega_{\perp}$, namely $\omega_{\perp}^c$, at which $\lambda_r$ crosses zero is at $2\pi\times 418\ \mathrm{Hz}$. Figure~\ref{fig:DSEigv415} shows the eigenvector corresponding to the largest real eigenvalue of Eq.~\eqref{eq:gpeeig} for the parameter $\omega_{\perp}/(2\pi)=415\ \mathrm{Hz}$ and $\omega_x=\omega_y$. Indeed, on the basis of the relevant unstable eigenmode, we expect the planar dark solitary wave structure to break its axial symmetry and reshape itself in the form of a vortex ring as a result of the relevant growing mode, recalling that the latter only breaks up towards a pair of solitonic vortices considerably later in the dynamical evolution. It is interesting to highlight that the above instability manifestation
appears to
be a natural 3d generalization of the destabilization
of a 2d planar dark soliton towards the formation of
a vortex dipole, explored from a bifurcation 
theory perspective, e.g., in~\cite{PhysRevA.77.053610,PhysRevA.82.013646}.

The most unstable eigenvalue (the one bearing the largest real part) can be computed from the GPE dynamics starting from the perturbed steady state solutions of Eq.~\eqref{eq:ste}~\cite{PhysRevA.105.053303}.
We evolve such waveforms starting with a small perturbation $\delta \psi$ 
[$\approx 10^{-10}$ relative to the $\mathrm{max}(|\psi'|)$].
As expected from the linear
stability analysis, we find that this perturbation grows exponentially as $\psi'=\psi_0'\exp(\lambda_r t')$ with an instability growth time $\tau=-\log(|\delta \psi|)/\lambda_r$.
The computed growth rate $\lambda_r$ via this ``direct integration'' approach,
indicated by circles on Fig.~\ref{fig:DSEigr}, matches well with the largest eigenvalues obtained from the spectral analysis of the BdG equations, thus confirming the robustness of our numerical findings.
The GPE dynamics for different $\omega_{\perp}$ shows that the strength of transverse confinement determines the resulting wave structures induced by the transverse instability of dark solitary waves.
We generally find that tighter confinement favors the formation of a solitonic vortex, while looser traps result in a long-lived vortex ring.
Figs.~\ref{fig:DSE} and \ref{fig:DSVR} show the vortex ring and solitonic vortex generated during the GPE dynamics with $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$ and $\omega_{\perp}/(2\pi)=250\ \mathrm{Hz}$.
Their existence can be further confirmed from the 3D condensate densities shown in Fig.~\ref{fig:DS3Dview}.
The single solitonic vortex formation from the imprinted dark soliton for higher $\omega_{\perp}$ is due to the asymmetric flow (see Fig.~\ref{fig:DSEigv415}) arising from the transverse instability that leads to the bending of density on the $\ey-\ez$ plane \cite{PhysRevLett.113.065302}.
   %
      %
   \begin{figure*}
 \includegraphics[width=\textwidth,angle =0]{Fig_DSE}
       \caption{Condensate densities (top) and phase cross-section (bottom) show the existence of a vortex ring at $t'=15$ after evolving the steady state dark soliton solution. Here, $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$ and $\omega_x=\omega_y$. (a) and (b) plot $|\psi(x',0,z')|^2$ and $|\psi(0,y',z')|^2$ respectively.} 
    \label{fig:DSE}
\end{figure*}
   \begin{figure*}
 \includegraphics[width=\textwidth,angle =0]{Fig_DSVR}
       \caption{Condensate densities (top) and phase cross-section (bottom) show a solitonic vortex at $t'=50$ after evolving the steady state dark soliton solution. Here, $\omega_{\perp}/(2\pi)=250\ \mathrm{Hz}$ and $\omega_x=\omega_y$. (a) and (b) plot $|\psi'(x',0,z')|^2$ and $|\psi'(0,y',z')|^2$ respectively.
      
       }
    \label{fig:DSVR}
\end{figure*}
   %
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_DS3Dview}
       \caption{The three dimensional view of the snapshots shown in Fig.~\ref{fig:DSE} (left) and Fig.~\ref{fig:DSVR} (right). The red and blue colors represent cuts of constant density and vorticity.
     }
    \label{fig:DS3Dview}
\end{figure}

\subsection{Traveling solutions}

   \begin{figure*}
     \includegraphics[width=\textwidth,angle =0]{Fig_vDSS}
       \caption{Condensate densities (top) and phase cross-section (bottom) of numerically exact traveling planar dark soliton solution for $\omega_{\perp}/(2\pi)=400$~Hz, where $\omega_x=\omega_y$ and $v=0.1c_s$. (a) and (b) plot $|\psi(x',0,z')|^2$ and $|\psi(0,y',z')|^2$ respectively.}
    \label{fig:vDSS}
\end{figure*}
  %
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_eigvelo.pdf}
       \caption{The largest real eigenvalues leading to the unstable dynamics of moving solitons as a function of $\omega_{\perp}$.The lines of the main figure (top to bottom) represent the cases, $v=(0,0.1,0.2,0.25,0.3)\times c_s$ obtained from the eigenvalue problem, while circles represent the growth rate obtained from the GPE dynamics. The inset shows the critical chemical potential $\mu_c$ as a function of $v/c_s$.
       }
    \label{fig:eigvelo}
\end{figure}

We now consider the case of traveling solitary waves. In this case, we find the steady solution of the 3D-GPE equation  
in the so-called co-traveling frame variant of
Eq.~\eqref{eq:ste} for which
\begin{align*}
 \mathcal{F}(\psi') &= \left[-\frac{1}{2}\nabla^2+iv\frac{\partial}{\partial z'}+V(\mathbf{r'})+g'|\psi'|^2 -\mu' \right] \psi' =0
\end{align*}
and $v$ represents the soliton velocity 
i.e., by seeking a numerically exact stationary solution
in the co-traveling frame along $\ez$, we obtain a
traveling solution of the original frame.
In this case, we fix the longitudinal length $L_z=120 a_{\perp}$ and remove the trap along the z-direction, so that
genuine traveling can be feasible along this direction. Fig.~\ref{fig:vDSS} shows the steady state solution for $v=0.1c_s$, where $c_s=\sqrt{\mu'}$ is the sound velocity, and $\omega_{\perp}/(2\pi)=400\ \mathrm{Hz}$.
This shows a phase jump which corresponds to the gray soliton~\cite{Frantzeskakis_2010} at the $\mathbf{r'}=0$ plane. 
The additional phase jump 
at the boundary is due to the implemented periodic boundary condition. 


Importantly, the computation of the traveling 
planar dark solitary wave solution is accompanied by
the solution of the BdG equations to 
identify the spectral stability of the relevant
state. 
The largest real eigenvalues leading to the unstable dynamics of moving solitons are shown in Fig.~\ref{fig:eigvelo} as a function of $\omega_{\perp}$. As the velocity increases, the $\lambda_r/l_z$ decreases and, in this way, leads to the  increase in the critical transverse frequency $\omega_c$, which representing the transition point of stable to unstable dynamics. The corresponding chemical potential $\mu_c = \mu/(\hbar \omega_{\perp}^c)$ is shown in the inset of Fig.~\ref{fig:eigvelo} as a function of $v/c_s$. This corroborates, through detailed 3D computations, 
the earlier dynamical results of~\cite{PhysRevLett.89.110401}. On the other hand, we find from the GPE dynamics that the transverse instability leads to the formation of solitonic vortices as in the case of $v=0$. Fig.~\ref{fig:DS3Dviewv02} shows the solitonic vortex at $t'=80$ formed during the dynamical evolution of the condensate density for $v=0.2c_s$.

   \begin{figure}
 \includegraphics[width=0.8\columnwidth,angle =0]{Fig_DS3Dviewv02}
       \caption{The three dimensional view of condensate density $|\psi'|^2$ at $t'=80$ for $v=0.2c_s$. The red and blue colors represent cuts of constant density and vorticity.
     
      }
    \label{fig:DS3Dviewv02}
\end{figure}

\section{Conclusions and outlook}

Motivated by experimental observations in 3D BECs indicative of the dynamical instability of planar dark solitons, we analyzed the stability of both stationary and traveling planar dark solitons.
We first illustrated the comparison of our 3D time-dependent numerical computations with the experimentally observed time-evolution. 
This indicated that, for the experimental parameters, the phase imprinting of a planar dark soliton gave rise to a vortex ring
that subsequently decayed into a solitonic vortex. 
This finding motivated us to explore the stability of planar dark solitons as a function of the transverse confinement, thereby quantifying the transition from the 3D to the 1D regime. 
Indeed, for sufficiently strong transverse confinement, 
our BdG analysis identified the quasi-1d limiting values of the corresponding excitation frequencies. This analysis further captured the destabilization of planar dark solitons in the 3D regime.
We identified the most rapidly growing unstable mode in two distinct ways{; one by the GPE direct simulation and another by the BdG analysis,} and monitored the resulting dynamics leading to the formation of solitonic vortices or vortex rings (depending on the transverse confinement).
 
Although this problem has been studied from a range of perspectives, there are still numerous open questions. 
Our primary focus here
is on the implications of the transverse confinement.
However, these structures are not exact solutions for non-zero longitudinal confinement, therefore their stability and properties as a function of longitudinal confinement 
merits further study. 
Moreover, a careful inspection of the movie in~\cite{supple} raises
intriguing questions on the scattering interactions between 
solitonic vortices in confined geometries.
In the same vein, the behavior of solitonic structures in box potentials~\cite{PhysRevLett.110.200406} is of considerable
current interest.
In this case, one can study how structures that are near-exact solutions in the bulk abruptly impinge on the system's edge.
Further possibilities could involve localized barriers~\cite{PhysRevA.94.063612} or even random potentials~\cite{doi:10.1073/pnas.1615004114} where the soliton's fate and trajectory are unknown.

\begin{acknowledgments}
We gratefully acknowledge G.~H.~Reid, and K.~Rodriguez for  meticulously reading the manuscript, R.~Carretero-Gonz\'{a}lez for fruitful discussions; and Mason Porter for long ago initiating our exchanges on this topic.
This material is based upon work supported by the US
National Science Foundation under Grants No. PHY-2110030 and DMS-2204702 (P.G.K.).
I.B.S. and A.F. were partially supported by the National Institute of Standards and Technology, as well as the NSF through: the Physics Frontier Center at the Joint Quantum Institute (PHY-1430094) and the Quantum Leap Challenge Institute for Robust Quantum Simulation (OMA-2120757).
\end{acknowledgments}


\section{Introduction}
Dark solitonic structures emerge ubiquitously in systems that combine dispersion with a self-defocusing nonlinearity.
Numerous physical examples encompassing these properties prominently include atomic Bose-Einstein condensates (BECs)~\cite{kevrekidis2015defocusing,Frantzeskakis_2010}, and nonlinear Kerr (and generalized non-Kerr) media~\cite{KIVSHAR199881}. 
More broadly, solitonic structures also arise in mechanical lattices of coupled pendula~\cite{PhysRevLett.68.1730}, electrical transmission lines~\cite{PhysRevE.49.828}, microwave magnetic thin films~\cite{PhysRevLett.70.1707,PhysRevLett.84.4697}, acoustic waveguides of Helmholtz resonators~\cite{PhysRevE.91.023204}, surface water waves~\cite{PhysRevLett.110.124101}, and nematic liquid crystals~\cite{Piccardi:11}.
Outside of these closed systems, solitonic waves are even present in media balancing gain and loss, such as polariton fluids in semiconductor microcavities~\cite{carusotto_2013}. The presence of solitary waves in such disparate contexts attests to the breadth of their relevance and the generality of the central mechanisms leading to their creation. In particular, atomic BECs hold a central role amongst these systems owing to their highly tunable interaction Hamiltonians, the control of dimensionality and potential landscapes, as well as the diverse available mechanisms for creating excitations~\cite{pethick,stringari,kevrekidis2015defocusing}.

There are different ways BECs are used to create solitonic states. For example, dragging a laser beam through a BEC~\cite{PhysRevLett.99.160405}, phase-imprinting~\cite{Becker2008}, and matter-wave interference~\cite{PhysRevLett.101.130401,PhysRevA.81.063604}. More recent experiments showed control of both the position and velocity of soliton~\cite{Fritsch2020creating}, with contemporary theoretical studies exploring 3-dimensional (3d) solitary waves in detail~\cite{mateo2015stability}. The intriguing dynamics of the resulting structures, including their breakup into vortical patterns~\cite{PhysRevLett.86.2926,Shomroni2009} (even in superfluid Fermi gases~\cite{PhysRevLett.113.065301}) is a subject of wide investigation and has a wide variety of potential applications ranging from atomic matter-wave interferometers~\cite{Negretti_2004} to producing two-level qubit systems~\cite{PhysRevA.95.053618}. Much of this experimental and theoretical effort is summarized in Refs.~\cite{Frantzeskakis_2010,kevrekidis2015defocusing}.

In the present work, we revisit dark solitary waves 3d repulsively interacting BECs highly elongated along the longitudinal axes $\ez$. These waves are associated with a reduction in the local atomic density accompanied with a modification of the condensate phase. In elongated systems, vortex rings have a similar density distribution but with a ring-shaped phase singularity~\cite{komineas,kevrekidis2015defocusing}. The simplest example is a planar (or kink) soliton where the density is reduced to be two dimensional on a plane normal to the soliton's direction of propagation. By contrast, a solitonic vortex is a highly anisotropic vortex with a linear phase-singularity~\cite{kevrekidis2015defocusing}.  
The potential conversion of kink solitons into these excitations is of wide interest~\cite{PhysRevA.68.043617,PhysRevA.65.043612,Becker_2013,tylutki2015solitonic,PhysRevA.94.023627,PhysRevLett.113.065302}.

We begin by comparing experimentally observed solitary waves with 3d numerical simulations computed with no adjustable parameters. Motivated
by this comparison, we theoretically explore the related 3d case of a longitudinally traveling solitonic wave in a system that is transversely confined but longitudinally infinite. Then, we compute the relevant Bogolyubov-de Gennes (BdG) modes to identify the onset of transverse instability and to quantify the growth rate of the unstable mode. This process confirms the prediction of Ref.~\cite{PhysRevLett.89.110401} that the interval of stability of kink solitons increases as their speed approaches the local speed of sound c. Finally in the unstable regime, we numerically follow the exponential growth of the unstable mode to its ultimate conversion to vortex structures.

\begin{figure*}[t]
\includegraphics{Fig_PHE2Din0}
	\caption{Integrated GPE wavefunction density (top) and phase cross-section (bottom) evaluated 

	immediately after a $t_p = 140\ \mu\mathrm{s}$ phase imprint.
	(a) and (b) plot  $|\psi(x,z)|^2$ and $|\psi(y,z)|^2$ respectively.}
    \label{fig:PHE2D0}
\end{figure*}

\section{Experimental Motivation: Phase Engineering and In-Trap Motion}

In the present work, we model the experimental procedure in reference~\cite{Fritsch2020creating} which describes a new wavefunction engineering approach for launching solitons in atomic BECs and show that dark planar solitons are first created and then relax into a solitonic vortex.

The experiments started with $N=2.4\times10^5$ $^{87}\text{Rb}$ BECs confined in an optical dipole trap with frequencies $(\omega_x,\omega_y,\omega_z)=2\pi \times (94.73,153.23,9.1)\ \text{Hz}$ and chemical potential $\mu=h\times1.1$~kHz. 
We created solitonic excitations by combining
conventional phase imprinting with density engineering, i.e., wavefunction engineering, yielding solitonic excitations with a wide range of initial velocities~\cite{Fritsch2020creating}.

We consider the fully $3d$ Gross-Pitaveskii model (GPE), complementing these experiments~\cite{pethick,stringari}
\begin{equation}
  i\hbar\frac{\partial\psi}{\partial t}= \left[-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{r})+g|\psi|^2 \right] \psi.
  \label{eq:3DGP}
\end{equation}
Here, $\psi$ is normalized to the total atom number $N_{\mathrm{3D}}$,
interaction constant and potential are defined as $g=4 \pi \hbar^2 a_s/m$ and $V(\mathbf{r})=m(\omega_x^2 x^2+\omega_y^2 y^2+\omega_z^2 z^2)/2+V_\mathrm{step}(t) \mathcal{H}(z)$, respectively, where $\mathcal{H}(z)$ represents the Heaviside function.
We fix the magnitude of the step potential $V_\mathrm{step}(t)$ to be $0$ except during the phase-imprinting pulse when it is $h\times5.5\ \mathrm{kHz}$.

We simulate the experiments using the parameters described above and create solitons via the phase imprinting method: applying a potential generated with an infinitely sharp edge (smaller than the BEC healing length) to half of the BEC.
We note that the density engineering method is implemented in the experiment to circumvent technical limitations, given that the BEC healing length is typically much smaller than the optical resolution of the imaging systems used for phase imprinting, which is not the case in numerical simulations.
We found that the idealization of a step potential centered at $z=0$ and amplitude $V_{\mathrm{step}}$, is sufficient to reproduce the experimental results.
This indirectly demonstrates the efficacy of the experimental approach.

\begin{figure*}
     \includegraphics[width=\textwidth,angle =0]{Fig_PHE2Dm}
       \caption {
Integrated GPE wavefunction density (top) along with phase (bottom) cross-sections evaluated $400\ \mathrm{ms}$ after a $t_p = 140\ \mu\mathrm{s}$ phase imprint.
(a) plots the $|\psi|^2$ projection and associated cross-sections in the $\ex-\ez$ plane and (b) plots $\ey-\ez$ quantities.
 }
    \label{fig:PHE2D}
\end{figure*}

\begin{figure}
     \includegraphics[width=\columnwidth,angle =0]{Fig_PHE1Dm}
     \caption{ Oscillation of the solitonic vortices (a) and oscillation frequencies (b).
    The black symbols on the grey colored integrated GPE wavefunction density in $\ez - t$ plane of (a) represent experimental data.
    Every set of experimental parameters was repeated three times, and a symbol is present for each clearly identified solitonic excitation. In cases in which more than one solitonic excitation was found in a single image we chose the position of the deepest depletion. The simulations use $t_p=140\ \mu\mathrm{s}$. 
     The symbols, red squares and black circles, in (b) represent GPE and experimental calculations, respectively, where the error bar represents the single-$\sigma$ uncertainty  obtained from the sinusoidal fits.
     }
    \label{fig:PHE1D}
\end{figure}


Figure~\ref{fig:PHE2D0} shows the 3d density projected onto the $\ex-\ez$ and $\ey-\ez$ planes, i.e., $|\psi(x,z)|^2$ and $|\psi(y,z)|^2$, and the associated planar phase cuts, immediately following the application of $V_\text{step}$ for a phase imprinting time $t_p = 140\ \mu\text{s}$. 
The phase-imprinted planar dark soliton subsequently evolves into a vortex ring as shown in ~\cite{supple}.
This ring eventually decays~\cite{PhysRevA.98.033609} into a pair of solitonic vortices aligned along the direction of stronger transverse confinement and exhibits a large amplitude oscillation ~\cite{supple}. 
The weak interaction of the resulting
structures with the background of excited phonons leads to effective (apparent) dissipative longitudinal dynamics and to the selection of two separate trajectories, one for each solitonic vortex as shown in Figs.~\ref{fig:PHE2D} and~\ref{fig:PHE1D}. 
{In the experimental results, over long times, the amplitude
of oscillation appears to be growing, possibly reflecting the (weak)
anti-damping effect due to the thermal fraction~\cite{prouka}.}
Both the numerical and experimental results (black crosses) find an oscillation frequency $\omega_s$

far from the ``canonical'' result of $\omega_s =  \omega_z / \sqrt{2}$ for planar solitons~\cite{PhysRevLett.84.2298,Frantzeskakis_2010}. 

Instead, the solitonic vortex oscillation frequency is 
significantly reduced~\cite{mateo2015stability,Fritsch2020creating}, e.g.,
 by a factor of about 3 in comparison to $\omega_z$, as shown in the bottom panel of Fig.~\ref{fig:PHE1D}. More generally, in line with~\cite{mateo2015stability}, we find that the relevant reduction is a function of the chemical potential, and the general agreement between the experimental data and the numerical simulations supports 
 solitonic vortex pair formation from the method of phase imprinting.
 
This experiment-theory correspondence prompts us to further explore the stability of planar solitons and their transformation into other structures.

\section{Stability of planar dark solitons}

In this section we study the stability of dark solitary waves using a dimensionless version of Eq.~\eqref{eq:3DGP}
\begin{equation}
  i\frac{\partial\psi'}{\partial t'}= \left[-\frac{1}{2}(\nabla')^2+V'(\mathbf{r}')+g'|\psi'|^2 \right] \psi',
  \label{eq:3DGPd}
\end{equation}
in terms of the dimensionless quantities $t'=\omega_{\perp} t$, $\mathbf{r}'=\mathbf{r}/a_{\perp}$, $\psi'=a_{\perp}^{3/2} \psi$, $V(\mathbf{r}')=[l_x^2 (x')^2+l_y^2 (y')^2+l_z^2 (z')^2]/2$, $g'=4 \pi a_s/a_{\perp}$, and $\mu'=\mu/(\hbar \omega_{\perp})$.
Here, we use the transverse degrees of freedom to define the natural units of frequency $\omega_{\perp}=\sqrt{\omega_x \omega_y}$ and length $a_{\perp}=\sqrt{\hbar/(m \omega_{\perp})}$. 
In addition, we introduce anisotropy parameters $l_x=\omega_x/\omega_{\perp}$, $l_y=\omega_y/\omega_{\perp}$ and $l_z=\omega_z/\omega_{\perp}$.

In line with the last section, we use a chemical potential $\mu= h\times 1.1\ \mathrm{kHz}$, longitudinal frequency $\omega_z/(2\pi)=10\ \mathrm{Hz}$, and transverse frequency ranging from $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$ to $1\ \mathrm{kHz}$. In what follows, we take $\omega_x=\omega_y$ for simplicity.

The 1D counterpart of Eq.~\eqref{eq:3DGP} 
with $V=0$ supports the dark soliton solution~\cite{kevrekidis2015defocusing}
\begin{equation}
 \psi_k'(z',t')=\sqrt{n_0} \left[B \tanh(\zeta)+i A\right]e^{i(kz'-\omega t')}.
  \label{eq:1D_dark}
\end{equation}
Here, $n_0=\mu'$, $\zeta=\sqrt{n_0}B[z'-z_0(t')]$, $z_0=vt'+Z_0$, $v=A\sqrt{n_0}+k$, $\omega=k^2/2+n_0$ and $A^2+B^2=1$, and restrict ourselves to $k=0$ solutions (of a stationary background). The remaining properties have the standard meaning,
i.e., $v$ is the soliton speed, $n_0$ the background density,
$z_0$ marks the center of mass (and $Z_0$ is its initial value), 
while $B$ is associated with the
inverse solitonic width~\cite{kevrekidis2015defocusing}.

\subsection{Stationary solutions}

In the presence of a trap, the wavefunction for a 3d stationary solution is approximately  
\begin{equation}
  \psi'(\mathbf{r'},t')=\psi_{ds}(z')\Phi(\mathbf{r'})e^{-i \mu' t'},
  \label{eq:1D_darks}
\end{equation}
  where $\Phi(\mathbf{r'})$ represents the ground state Thomas-Fermi (TF) wavefunction~\cite{pethick,stringari} and 
  $\psi_{ds}$ represents the stationary dark soliton from
  Eq.~(\ref{eq:1D_dark}), associated with $A=v=0$. 
\begin{figure*}
     \includegraphics{Fig_DSS}
       \caption{Condensate GPE density of the Eq.~\eqref{eq:ste} (top) and phase cross-section (bottom) show the stationary 3D dark soliton solution for $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$, where $\omega_x=\omega_y$ and $v=0$. The lower panel shows the soliton phase jump from  $\pi$ (dark) to 0 (light). 
       (a) and (b) plot $|\psi'(x',0,z')|^2$ and $|\psi'(0,y',z')|^2$ respectively.
       }
    \label{fig:DSS}
\end{figure*}
Using Eq.~\eqref{eq:1D_darks} as an initial guess for a root finding algorithm, we obtain numerically exact solutions by solving
\begin{equation}
   \mathcal{F}(\psi') = \left[-\frac{1}{2}\nabla^2+V(\mathbf{r'})+g'|\psi'|^2 -\mu' \right]\psi'=0.
      \label{eq:ste}
  \end{equation}
Fig.~\ref{fig:DSS} shows the solution of Eq.~\eqref{eq:ste} for $\omega_{\perp}/(2\pi)=100$~Hz. We now 
explore both the stable and unstable modes of $\psi'$ by identifying the eigenvalues associated with the BdG equations~\cite{pethick,stringari,kevrekidis2015defocusing} 
for the parameter range of $\omega_{\perp}$.
For relevant details, see the Appendix \ref{app:ES}. An eigenvalue $\lambda$ has a  real part denoted by $\lambda_r$ and an imaginary part denoted by $\lambda_i$. It is important to recall that a non-vanishing $\lambda_r$, given the  Hamiltonian nature of our GPE (\ref{eq:3DGPd}), is tantamount to the presence of a dynamical instability.
 %
  %
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_DSEigim}
       \caption{Imaginary eigenvalues $\lambda_i$ obtained from the BdG analysis  plotted as a function of $\omega_{\perp}$ for $\omega_x=\omega_y$. 
       The vertical dashed black line represents the mode responsible for instability of the planar
dark soliton, while blue lines represent stable modes. In the limit of large $\omega_{\perp}$, the 1d analytically available results, indicated by the red short lines, are retrieved (see text for further details).
   
     
      }
    \label{fig:DSEigim}
\end{figure}
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_DSEigr}
       \caption{The largest real eigenvalue $\lambda_r$, 
       associated with the
       instability of stationary, planar dark soliton, as a function of $\omega_{\perp}$ for $\omega_x=\omega_y$. The red circles are obtained from GPE dynamics and blue line represents BdG analysis. 
      
     
     
       }
    \label{fig:DSEigr}
\end{figure}

There exist two limiting cases for the stability problem of the dark solitary waves in a trapped condensate, more
specifically, the TF limit at large chemical
potential $\mu'$ and the linear limit at small $\mu'$. 
As is well-known~\cite{kevrekidis2015defocusing} 
in the 1d limit of the GPE model, the anomalous mode ($\omega_0$), corresponding to negative energy is associated with the oscillation of the single dark soliton with a frequency of $\approx l_z /\sqrt{2}$~\cite{PhysRevLett.84.2298}.
The other modes $\omega_n$ (associated with the
vibration modes of the background cloud) are expected to approach the frequencies of $l_z\sqrt{n(n+1)/2}$~\cite{stringari}. On the other hand, in the linear limit of vanishing density $|\psi'|^2 \rightarrow 0$, the eigenvalue problem of Eq.~\eqref{eq:3DGP} reduces to that of a {\it linear} quantum harmonic oscillator with corresponding energies of eigenmodes $|k,l,m\rangle$ 
\begin{align}
E_{k,l,m} &= \left(k+ \frac{1}{2}\right) l_x+ \left(l+ \frac{1}{2}\right) l_y+
\left(m + \frac{1}{2}\right) l_z.
\end{align}
A single dark soliton at $z'=0$ (pertaining to the $|0,0,1\rangle$ in the linear limit) has energy $E_{0,0,1} = 
(l_x + l_y+ 3 l_z)/2$.
To explore the validity of these predictions, we rescale the imaginary eigenvalues by $l_z$ and show the corresponding results in Fig.~\ref{fig:DSEigim}. 
Indeed, our numerically calculated imaginary eigenvalues confirm the expectations on the basis of the 1d
predictions, as our numerical continuation results
approach progressively the 1d TF limit as $\omega_{\perp}$
increases. We note that the chemical potential controls the approach to the TF limit for the present setting. The mode responsible for the instability of the planar
dark soliton appears as a nearly vertical dashed black line in 
Fig.~\ref{fig:DSEigim}.


   \begin{figure*}
 \includegraphics[width=\textwidth,angle =0]{Fig_DSEigv415}
       \caption{(left)  The eigenvector $|\varphi|^2$ corresponds to the largest real eigenvalue of Eq.~\eqref{eq:gpeeig} and its phase profile for the parameter $\omega_{\perp}/(2\pi)=415$~Hz and $\omega_x=\omega_y$. (right panel) The density $|\psi'+\epsilon \varphi|^2$, where $\psi'$ is the steady state wavefunction and $\epsilon=20$. Lower panels in (a) and (b) show the phase cross-section.
     
      }
    \label{fig:DSEigv415}
\end{figure*}
The analysis of real eigenvalues displayed in Fig.~\ref{fig:DSEigr} shows that, as $\omega_{\perp}$  decreases (going away from the above mentioned quasi-1d limit), there exists an instability window for smaller values of the relevant parameter. The critical $\omega_{\perp}$, namely $\omega_{\perp}^c$, at which $\lambda_r$ crosses zero is at $2\pi\times 418\ \mathrm{Hz}$. Figure~\ref{fig:DSEigv415} shows the eigenvector corresponding to the largest real eigenvalue of Eq.~\eqref{eq:gpeeig} for the parameter $\omega_{\perp}/(2\pi)=415\ \mathrm{Hz}$ and $\omega_x=\omega_y$. Indeed, on the basis of the relevant unstable eigenmode, we expect the planar dark solitary wave structure to break its axial symmetry and reshape itself in the form of a vortex ring as a result of the relevant growing mode, recalling that the latter only breaks up towards a pair of solitonic vortices considerably later in the dynamical evolution. It is interesting to highlight that the above instability manifestation
appears to
be a natural 3d generalization of the destabilization
of a 2d planar dark soliton towards the formation of
a vortex dipole, explored from a bifurcation 
theory perspective, e.g., in~\cite{PhysRevA.77.053610,PhysRevA.82.013646}.

The most unstable eigenvalue (the one bearing the largest real part) can be computed from the GPE dynamics starting from the perturbed steady state solutions of Eq.~\eqref{eq:ste}~\cite{PhysRevA.105.053303}.
We evolve such waveforms starting with a small perturbation $\delta \psi$ 
[$\approx 10^{-10}$ relative to the $\mathrm{max}(|\psi'|)$].
As expected from the linear
stability analysis, we find that this perturbation grows exponentially as $\psi'=\psi_0'\exp(\lambda_r t')$ with an instability growth time $\tau=-\log(|\delta \psi|)/\lambda_r$.
The computed growth rate $\lambda_r$ via this ``direct integration'' approach,
indicated by circles on Fig.~\ref{fig:DSEigr}, matches well with the largest eigenvalues obtained from the spectral analysis of the BdG equations, thus confirming the robustness of our numerical findings.
The GPE dynamics for different $\omega_{\perp}$ shows that the strength of transverse confinement determines the resulting wave structures induced by the transverse instability of dark solitary waves.
We generally find that tighter confinement favors the formation of a solitonic vortex, while looser traps result in a long-lived vortex ring.
Figs.~\ref{fig:DSE} and \ref{fig:DSVR} show the vortex ring and solitonic vortex generated during the GPE dynamics with $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$ and $\omega_{\perp}/(2\pi)=250\ \mathrm{Hz}$.
Their existence can be further confirmed from the 3D condensate densities shown in Fig.~\ref{fig:DS3Dview}.
The single solitonic vortex formation from the imprinted dark soliton for higher $\omega_{\perp}$ is due to the asymmetric flow (see Fig.~\ref{fig:DSEigv415}) arising from the transverse instability that leads to the bending of density on the $\ey-\ez$ plane \cite{PhysRevLett.113.065302}.
   %
      %
   \begin{figure*}
 \includegraphics[width=\textwidth,angle =0]{Fig_DSE}
       \caption{Condensate densities (top) and phase cross-section (bottom) show the existence of a vortex ring at $t'=15$ after evolving the steady state dark soliton solution. Here, $\omega_{\perp}/(2\pi)=100\ \mathrm{Hz}$ and $\omega_x=\omega_y$. (a) and (b) plot $|\psi(x',0,z')|^2$ and $|\psi(0,y',z')|^2$ respectively.} 
    \label{fig:DSE}
\end{figure*}
   \begin{figure*}
 \includegraphics[width=\textwidth,angle =0]{Fig_DSVR}
       \caption{Condensate densities (top) and phase cross-section (bottom) show a solitonic vortex at $t'=50$ after evolving the steady state dark soliton solution. Here, $\omega_{\perp}/(2\pi)=250\ \mathrm{Hz}$ and $\omega_x=\omega_y$. (a) and (b) plot $|\psi'(x',0,z')|^2$ and $|\psi'(0,y',z')|^2$ respectively.
      
       }
    \label{fig:DSVR}
\end{figure*}
   %
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_DS3Dview}
       \caption{The three dimensional view of the snapshots shown in Fig.~\ref{fig:DSE} (left) and Fig.~\ref{fig:DSVR} (right). The red and blue colors represent cuts of constant density and vorticity.
     }
    \label{fig:DS3Dview}
\end{figure}

\subsection{Traveling solutions}

   \begin{figure*}
     \includegraphics[width=\textwidth,angle =0]{Fig_vDSS}
       \caption{Condensate densities (top) and phase cross-section (bottom) of numerically exact traveling planar dark soliton solution for $\omega_{\perp}/(2\pi)=400$~Hz, where $\omega_x=\omega_y$ and $v=0.1c_s$. (a) and (b) plot $|\psi(x',0,z')|^2$ and $|\psi(0,y',z')|^2$ respectively.}
    \label{fig:vDSS}
\end{figure*}
  %
   \begin{figure}
 \includegraphics[width=\columnwidth,angle =0]{Fig_eigvelo.pdf}
       \caption{The largest real eigenvalues leading to the unstable dynamics of moving solitons as a function of $\omega_{\perp}$.The lines of the main figure (top to bottom) represent the cases, $v=(0,0.1,0.2,0.25,0.3)\times c_s$ obtained from the eigenvalue problem, while circles represent the growth rate obtained from the GPE dynamics. The inset shows the critical chemical potential $\mu_c$ as a function of $v/c_s$.
       }
    \label{fig:eigvelo}
\end{figure}

We now consider the case of traveling solitary waves. In this case, we find the steady solution of the 3D-GPE equation  
in the so-called co-traveling frame variant of
Eq.~\eqref{eq:ste} for which
\begin{align*}
 \mathcal{F}(\psi') &= \left[-\frac{1}{2}\nabla^2+iv\frac{\partial}{\partial z'}+V(\mathbf{r'})+g'|\psi'|^2 -\mu' \right] \psi' =0
\end{align*}
and $v$ represents the soliton velocity 
i.e., by seeking a numerically exact stationary solution
in the co-traveling frame along $\ez$, we obtain a
traveling solution of the original frame.
In this case, we fix the longitudinal length $L_z=120 a_{\perp}$ and remove the trap along the z-direction, so that
genuine traveling can be feasible along this direction. Fig.~\ref{fig:vDSS} shows the steady state solution for $v=0.1c_s$, where $c_s=\sqrt{\mu'}$ is the sound velocity, and $\omega_{\perp}/(2\pi)=400\ \mathrm{Hz}$.
This shows a phase jump which corresponds to the gray soliton~\cite{Frantzeskakis_2010} at the $\mathbf{r'}=0$ plane. 
The additional phase jump 
at the boundary is due to the implemented periodic boundary condition. 


Importantly, the computation of the traveling 
planar dark solitary wave solution is accompanied by
the solution of the BdG equations to 
identify the spectral stability of the relevant
state. 
The largest real eigenvalues leading to the unstable dynamics of moving solitons are shown in Fig.~\ref{fig:eigvelo} as a function of $\omega_{\perp}$. As the velocity increases, the $\lambda_r/l_z$ decreases and, in this way, leads to the  increase in the critical transverse frequency $\omega_c$, which representing the transition point of stable to unstable dynamics. The corresponding chemical potential $\mu_c = \mu/(\hbar \omega_{\perp}^c)$ is shown in the inset of Fig.~\ref{fig:eigvelo} as a function of $v/c_s$. This corroborates, through detailed 3D computations, 
the earlier dynamical results of~\cite{PhysRevLett.89.110401}. On the other hand, we find from the GPE dynamics that the transverse instability leads to the formation of solitonic vortices as in the case of $v=0$. Fig.~\ref{fig:DS3Dviewv02} shows the solitonic vortex at $t'=80$ formed during the dynamical evolution of the condensate density for $v=0.2c_s$.

   \begin{figure}
 \includegraphics[width=0.8\columnwidth,angle =0]{Fig_DS3Dviewv02}
       \caption{The three dimensional view of condensate density $|\psi'|^2$ at $t'=80$ for $v=0.2c_s$. The red and blue colors represent cuts of constant density and vorticity.
     
      }
    \label{fig:DS3Dviewv02}
\end{figure}

\section{Conclusions and outlook}

Motivated by experimental observations in 3D BECs indicative of the dynamical instability of planar dark solitons, we analyzed the stability of both stationary and traveling planar dark solitons.
We first illustrated the comparison of our 3D time-dependent numerical computations with the experimentally observed time-evolution. 
This indicated that, for the experimental parameters, the phase imprinting of a planar dark soliton gave rise to a vortex ring
that subsequently decayed into a solitonic vortex. 
This finding motivated us to explore the stability of planar dark solitons as a function of the transverse confinement, thereby quantifying the transition from the 3D to the 1D regime. 
Indeed, for sufficiently strong transverse confinement, 
our BdG analysis identified the quasi-1d limiting values of the corresponding excitation frequencies. This analysis further captured the destabilization of planar dark solitons in the 3D regime.
We identified the most rapidly growing unstable mode in two distinct ways{; one by the GPE direct simulation and another by the BdG analysis,} and monitored the resulting dynamics leading to the formation of solitonic vortices or vortex rings (depending on the transverse confinement).
 
Although this problem has been studied from a range of perspectives, there are still numerous open questions. 
Our primary focus here
is on the implications of the transverse confinement.
However, these structures are not exact solutions for non-zero longitudinal confinement, therefore their stability and properties as a function of longitudinal confinement 
merits further study. 
Moreover, a careful inspection of the movie in~\cite{supple} raises
intriguing questions on the scattering interactions between 
solitonic vortices in confined geometries.
In the same vein, the behavior of solitonic structures in box potentials~\cite{PhysRevLett.110.200406} is of considerable
current interest.
In this case, one can study how structures that are near-exact solutions in the bulk abruptly impinge on the system's edge.
Further possibilities could involve localized barriers~\cite{PhysRevA.94.063612} or even random potentials~\cite{doi:10.1073/pnas.1615004114} where the soliton's fate and trajectory are unknown.

\begin{acknowledgments}
We gratefully acknowledge G.~H.~Reid, and K.~Rodriguez for  meticulously reading the manuscript, R.~Carretero-Gonz\'{a}lez for fruitful discussions; and Mason Porter for long ago initiating our exchanges on this topic.
This material is based upon work supported by the US
National Science Foundation under Grants No. PHY-2110030 and DMS-2204702 (P.G.K.).
I.B.S. and A.F. were partially supported by the National Institute of Standards and Technology, as well as the NSF through: the Physics Frontier Center at the Joint Quantum Institute (PHY-1430094) and the Quantum Leap Challenge Institute for Robust Quantum Simulation (OMA-2120757).
\end{acknowledgments}

