\section{Introduction}
The theory of wrapping one object around another, known as homotopy in topology~\cite{NashBook,Nakahara}, has played an important role in characterizing physical phenomena. Topological excitations~\cite{Nakahara,Chaikin03}, such as vortices in superconductors and superfluids or defects in liquid crystals, have been characterized according to the winding number by monitoring the phase change of the order parameter around a loop enclosing the excitation. More recently, homotopy has been used to classify topological matter by counting the winding in a target space as the Brillouin zone is traversed, leading to a periodic table of topological systems~\cite{Chiu2016,Stanescu_book}. When applied to spinor Bose-Einstein condensates (BECs), homotopy helps elucidate the topological mechanisms behind exotic excitations~\cite{Ueda14}.

Here we characterize a class of solitons in a spinor BEC on a spherical shell with a topological winding number coming from a mapping that wraps the order-parameter space around the real 2D space. Historically, similar solitons, call lumps, have been studied on a plane with specific boundary condition at infinity to make the space topologically equivalent to a sphere~\cite{Manton}. Placing the lump soliton on a spherical shell will naturally realize the setup and makes the topological meaning of the solitons more transparent.  Various solitons have been realized~\cite{Khaykovich02,Strecker02,Eiermann04,Marchant13,Lamporesi13,Medley14,Meier16,Aycock17,Bersano18,BH21,Hung23} or proposed~\cite{Wu04, Pedri05,Huang05,Li05,Damski10,Maucher11,Scott11,Pu15,Dutta17,Kengne21,Mostaan22} in cold atoms. As reviewed in Ref.~\cite{Lundblad22}, spherical-shell potentials for ultracold atoms can be realized by either coupling radio-frequency excitation to multi-component atoms in harmonic potentials to create a bubble trap or producing a phase-separation structure of multi-species atoms in a harmonic trap with the core occupied by one species while the spherical shell occupied by the desired species. We envision a spinor BEC on a spherical shell will be realized and analyze its topological solitons by finding the minimal-energy solutions with quantized winding numbers. While atomic spinor BECs have been realized and analyzed extensively~\cite{Kawaguchi12,SK13}, spherical bubble traps of BEC in microgravity have recently been achieved~\cite{Carollo22}.

In the limit of infinite contact interactions, analytic solutions of the lumps with quantized winding numbers can be constructed explicitly, which elucidate the wrapping between the real space and order-parameter space via the vector of the BEC wavefunction. When the interactions are finite, the lump solitons are solutions to a set of coupled nonlinear equations. In the strong interaction regime, we still find analytic solutions using the spherical harmonics. In the intermediate and weak interaction regimes, we found another set of lumps via numerical solutions. While the lumps carry quantized winding numbers, their energies are above the uniform ground state. Our calculations suggest higher-winding lumps are energetically stable against decomposing into lower-winding lumps, which is in contrast to quantum vortices in superconductors or superfluids~\cite{FetterBook,Pethick_book}. We mention there are other means for generating topological properties via manipulations in real space, such as inducing band topology by modulating real-space patterns~\cite{PhysRevLett.109.215302,PhysRevA.97.023618,Song20}.



\section{Theoretical framework}
According to Refs.~\cite{Kawaguchi12,SK13}, the mean-field energy functional of a three-dimensional $F=1$ spinor BEC is
\begin{eqnarray}\label{eq-F1}
\tilde{E}&=&\int d^3\tilde{x} \left\{\sum_{m=-1}^{1}\phi_m^*\left(-\frac{\hbar^2\nabla^2}{2M}-\mu_m-pm+qm^2\right)\phi_m \right. \nonumber \\
& &\left. +\frac{\tilde{c}_0}{2}n^2+\frac{c_1}{2}|\mathbf{F}|^2 \right\}.
\end{eqnarray}
Here $p$ and $q$ denote the linear and quadratic Zeeman splittings, $n=\sum_{m=-1}^{1}|\phi_m|^2$ is the total density, and $\mathbf{F}=(F_x, F_y, F_z)$ denotes the magnetization vector. We consider a simplified case with $c_1=0$, ignore the Zeeman splitting terms, and assume the populations of the three components are equal by introducing a uniform chemical potential $\mu$ for all components. Both $\mu$ and $c_0$ are positive for a stable BEC.


Next, we consider the $F=1$ spinor BEC is confined on a spherical surface $\tilde{S}^2$ with radius $R$ and replace $\tilde{x}^\nu\to Rx^\nu$ with $\nu=1, 2$. The derivatives are scaled accordingly. Moreover, $\tilde{\phi}_a\to\mathcal{M}\phi_a$ and $\tilde{E}\to E\tilde{E}_0$, where $\tilde{E}_0=\hbar^2\mathcal{M}^2/(2M)$ and $\mathcal{M}^2=\mu/c_0$. Following the derivation summarized in the Appendix, the dimensionless energy functional is
\begin{eqnarray}
E=\int_{S^2}d^2x \sqrt{g}\Big[g^{\mu\nu}\partial_\mu\phi_a \partial_\nu\phi_a+\frac{c_0}{2}(\phi_a^2-1)^2\Big].
\label{eq-en-1}
\end{eqnarray}
Here we define $c_0=2M \tilde{c}_0(R\mathcal{M})^2/\hbar^2$. The integral is over the unit sphere $S^2$ with the metric $g_{\mu\nu}$ and $g=\det(g_{\mu\nu})$. 





The energy functional is similar to the $O(3)$ non-linear sigma model \cite{Polyakov,Rajaraman}, which can support topological solitons called ``lumps'', as discussed in  Ref.~\cite{Manton}.
For the spinor BEC on a sphere, the vector field formed by $\phi^a(x)$ defines a map from the spatial two-sphere parameterized by $x=(x^1, x^2)$ to the sphere of $\phi^a\phi^a=1$ in the order-parameter space. The topology of this mapping is given by the homotopy group $\pi_2(S^2)=\mathbb{Z}$ \cite{NashBook,Nakahara}, which can be characterized by the integer winding number
\begin{eqnarray}\label{eq-WN}
N=\frac1{8\pi}\int_{S^2}d^2x \epsilon_{\mu\nu}\epsilon_{abc}\phi_a \partial_\mu\phi_b \partial_\nu\phi_c.
\end{eqnarray}
Here $\epsilon_{\mu\nu}$ and $\epsilon_{abc}$ are the Levi-Civita symbols of the real-space coordinates and order-parameter indices, respectively.


\section{Spherical lump soliton in the large $c_0$ limit}
To appreciate the topological properties of the lump soliton, we first consider the limit of $c_0\to\infty$, where exact solutions exist and can be visualized. 
The condition $\phi_a\phi_a=1$ must be satisfied to minimize the energy functional if $c_0\rightarrow\infty$. Thus, the above spinor BEC model becomes the genuine $O(3)$ non-linear sigma model on a spherical shell. 
In the literature~\cite{Manton}, the lumps are defined on a 2D plane $\mathbb{R}^2$, where the condition $\phi_a\to$ constant as $|x|\to\infty$ has been imposed to map the plane to an infinite sphere. A brief overview of the planar case is given in the Appendix. Here the spinor BEC is placed on a finite spherical shell.
We now look for configurations of $\phi_a(x)$ that minimize the total energy and also give rise to non-vanishing integers of the winding number $N$.

It is convenient to convert $\phi_a$ to a complex field according to the stereographic projection defined by
\begin{eqnarray}
R(x)=\frac{\phi_1+i\phi_2}{1+\phi_3},\quad
R^*(x)=\frac{\phi_1-i\phi_2}{1+\phi_3}.
\end{eqnarray}
The spherical coordinates with $x_1=\theta$ and $x_2=\varphi$ correspond to the metric on a unit sphere given by $g_{\mu\nu}=diag(h_1^2, h_2^2)$ with $h_1=1$ and $h_2=\sin\theta$.
We will focus on the case when the Bogomolny bound \cite{Belavin} is saturated with the following equations. (See the Appendix for a derivation.)
\begin{eqnarray}\label{eq-Req}
h_2\partial_1 R=-i h_1\partial_2 R,\quad
h_2\partial_1 R^*=i h_1\partial_2 R^*.
\label{eq-CR}
\end{eqnarray}
Since the equations are complex conjugates, we only focus on the first equation. 


From Eq.~\eqref{eq-CR} and $u=\tan\frac{\theta}{2}e^{i\varphi}$, it follows 
$\Big(\sin\theta\frac{\partial}{\partial\theta}+i\frac{\partial}{\partial\varphi}\Big)u=0$.
Thus, the generic lump solution on a unit sphere can be taken as
\begin{eqnarray}
R(u)=\frac{p(u)}{q(u)},
\end{eqnarray}
where $p(u)$ and $q(u)$ are polynomials of $u$ with no common factors. Clearly, $R(u)$ satisfies Eq.~(\ref{eq-CR}). In terms of $R(u)$, we find 
$\phi_1=\frac{2\text{Re}(R)}{1+|R|^2}$, $\phi_2=\frac{2\text{Im}(R)}{1+|R|^2}$, and $\phi_3=\frac{1-|R|^2}{1+|R|^2}$.
As an example, we consider a general solution describing a single lump centered at $(\theta_0,\varphi_0)$ with winding number $N=1$ given by
$R(u)=C(u-u_0)$.
Here $C$ is a complex constant and $u_0=\tan(\theta_0/2)e^{i\varphi_0}$. For simplicity, we can assume $C=1$ and $u_0=0$, so $R(u)=u$ that translates to the configuration of the $\phi_a$ fields as
\begin{eqnarray}\label{eq:N1Inf}
\phi_1=\sin\theta\cos\varphi,\quad \phi_2=\sin\theta\sin\varphi,\quad \phi_3=\cos\theta.
\end{eqnarray}
This is nothing but the identity map from the real-space $S^2$ to the order-parameter space $S^2$ defined by $\phi_a^2=1$. 
This simple soliton solution is visualized in the left panel of Figure \ref{fig-sphere}, where the vector field of $\phi_a$ on a unit sphere is represented by the arrows spreading around the whole surface. The configuration is similar to the electric field of a point charge or a hedgehog.

\begin{figure}
\centering
\includegraphics[width=0.45\columnwidth]{sphere.eps}
\includegraphics[width=0.45\columnwidth]{sphere-1.eps}
\caption{ Vector fields of $\phi_a$ for the $N=1$ lump soliton~\eqref{eq:N1Inf} (left) and the $N=2$ lump soliton~\eqref{eq:N2Inf} (right) on a unit sphere. }
\label{fig-sphere}
\end{figure}

A $N=2$ lump solution can be constructed by setting $R(u)=u^2$, which translates to the configuration of the $\phi_a$ fields as
\begin{eqnarray}\label{eq:N2Inf}
\phi_1&=&\frac{\sin^2\theta}{1+\cos^2\theta}\cos2\varphi,\quad
\phi_2=\frac{\sin^2\theta}{1+\cos^2\theta}\sin2\varphi, \nonumber \\ \phi_3&=&\frac{2\cos\theta}{1+\cos^2\theta}.
\end{eqnarray}
A direct calculation verifies that the winding number is indeed $N=2$ for this solution. The configuration of the $\phi_a$ fields of this lump solution is shown in the right panel of Figure \ref{fig-sphere}, where the order-parameter space wraps twice around the real space. The energies of the $N=1$ and $N=2$ lump solitons are $E_1=8\pi$ and $E_2=16\pi$, respectively. The linear energy dependence on $N$ is different from the quadratic dependence of the vortex energy on its winding number~\cite{FetterBook,Pethick_book} and implies energy neutrality by converting a $N=2$ lump to two $N=1$ lumps.




\begin{figure}
\centering
\includegraphics[width=\columnwidth]{fg.eps}
\caption{$f$ and $g$ of the $N=-1$ (Black line) and $N=-2$ (Red dashed line) solutions of Eq.~(\ref{eq-phi}) on a unit sphere as functions of $w=\cos\theta$. Here $c_0=1$.}
\label{fg}
\end{figure}

\section{Spherical lump solitons with finite $c_0$}
In realistic cases with finite $c_0$, the extremum condition of the energy functional leads to the following set of non-linear differential equations with $a=1,2,3$.
\begin{eqnarray}
-\frac{1}{\sqrt{g}}\partial_\mu\Big(\sqrt{g}g^{\mu\nu}\partial_\nu\phi_a\Big)+c_0(\phi_b\phi_b-1)\phi_a=0.
\label{eq-fB}
\end{eqnarray}
The extremum condition admits a trivial solution $\phi_a=v_a$, where $v_a$ is a constant unit vector in 3D. For the trivial solution, the energy is zero, which indicates it is the ground state without any excitation. 
We will present several solutions of the equations with non-zero winding numbers. However, finding solutions of non-linear differential equations is in general very difficult. For convenience, we rewrite the spherical Laplacian operator $\nabla_s^2\equiv-\frac{1}{\sqrt{g}}\partial_\mu\sqrt{g}g^{\mu\nu}\partial_\nu$ on a unit sphere as $-\Big(\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\sin\theta\frac{\partial}{\partial\theta}+\frac{1}{\sin^2\theta}\frac{\partial^2}{\partial\varphi}\Big)$.
The eigenfunctions of the spherical Laplacian operator are the spherical harmonics:
$\nabla_s^2Y_{lm}(\theta,\varphi)=l(l+1)Y_{lm}(\theta,\varphi)$.

In terms of $Y_{lm}(\theta,\varphi)$, we can find one solution of Eq.~(\ref{eq-fB}) inspired by the infinite-$c_0$ solution~\eqref{eq:N1Inf} as
\begin{eqnarray}
&&\phi_1=A\sin\theta\cos\varphi=\frac{A}{2}(Y_{1,1}+Y_{1,-1}),\nonumber\\
&&\phi_2=A\sin\theta\sin\varphi=\frac{A}{2i}(Y_{1,1}-Y_{1,-1}),\nonumber\\
&&\phi_3=A\cos\theta=AY_{1,0}.
\label{eq-s1}
\end{eqnarray}
Here $A$ is a real constant, and we ignore the normalization factors of $Y_{lm}$ at this moment. Then $\phi_a^2=A^2$, and Eq.~(\ref{eq-fB}) reduces to
$2\phi_a+c_0(A^2-1)\phi_a=0$.
Hence, we find $A=\sqrt{1-\frac2{c_0}}$, which requires that $c_0>2$. One can verify that the winding number for this solution is $N=1$, and its configuration is similar to Fig.~\ref{fig-sphere} (a).


Making use of the extremum condition \eqref{eq-fB}, the kinetic energy $E_{\text{kin}}=\int_{S^2}d^2x\phi_a\Big(-\partial_\mu\sqrt{g}g^{\mu\nu}\partial_\nu\Big)\phi_a$ takes the form $\int_{S^2}d^2x\sqrt{g}c_0(\phi_b^2-1)\phi_a^2$.
Therefore, the total energy with finite $c_0$ is
\begin{eqnarray}
E=\int_{S^2}\sqrt{g}\Big[\frac{c_0}{2}(1-\phi_a^2)(1+\phi_a^2)\Big].
\end{eqnarray}
For the $N=1$ lump soliton solution from Eq.~(\ref{eq-s1}), we have $\phi_b^2=A^2<1$, and the energy is 
$E=8\pi\Big(1-\frac{1}{c_0}\Big)>0$.
Therefore, the lump soliton is an excitation with positive energy, implying that the existence of a soliton solution requires some finite energy provided to the spinor BEC. Moreover, the excitation energy approaches the Bogomolny bound when $c_0\rightarrow\infty$. We remark that searching for a $N=2$ lump solution from an analogue of Eq.~\eqref{eq:N2Inf} with finite $c_0$ is challenging because there are five $Y_{2,m}$ functions but only three $\phi_a$ components. In the following, we will use a different ansatz to find more lump solutions. 




\begin{figure}
\centering
\includegraphics[width=\columnwidth]{vector.eps}
\caption{Projections of the $\phi_a$ fields from the north hemisphere onto the $xy$ plane for the $N=-1$ (left) and $N=-2$ (right) lump soliton solutions of Eq.~\eqref{eq-phi}. All projected vectors have been normalized to unit length.}
\label{fig-vec}
\end{figure}

Inspired by the $N=1$ lump solution from the spherical harmonics, we seek a more general solution of Eq.~(\ref{eq-fB}) with the following form
\begin{eqnarray}
\phi_1=f(\theta)\cos(n\varphi),\quad \phi_2=f(\theta)\sin(n\varphi),\quad \phi_3=g(\theta).
\label{eq-phi}
\end{eqnarray}
Here $f(\theta)$ and $g(\theta)$ are two unknown functions to be determined, and $n$ is some positive integer. By defining $w=\cos\theta$ and substituting the functional form into Eq.~(\ref{eq-fB}), we find two coupled equations with $h=f, g$:
\begin{eqnarray}
\frac{d}{dw}\Big[(1-w^2)\frac{dh}{dw}\Big]-\frac{n^2}{1-w^2}h-c_0(f^2+g^2-1)h=0. 
\label{eq-fg}
\end{eqnarray}
Since some coefficients of the equations are singular at $w=\pm1$, we impose the boundary condition allowed on a sphere that both $f(\pm1)$ and $g(\pm1)$ are finite.
The numerical solutions of $f(w)$ and $g(w)$  with $n=1$ and $c_0=1$ are shown in Figure \ref{fg}. 
Since the norm of the vector $\phi_a$ is no longer a constant, we introduce the unit vector $\hat{\phi}_a=\phi_a/\sqrt{\phi_a^2}$, which defines a map from the real-space $S^2$ to the unit sphere of $\hat{\phi}_a$. Then we can evaluate the winding number using Eq.~\eqref{eq-WN}, which becomes
$N
=\frac{n}{2}\int_{-1}^1 dw\,\sqrt{1-w^2}\frac{f(g f'-f g')}{(f^2+g^2)^{3/2}}$with $f'=df/dw$. For the $n=1$ solution of Eq.~\eqref{eq-phi}, we find $N=-0.97$, which is close to $N=-1$ and indicates the result as a $N=-1$ lump soliton on a sphere. 



\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{E-c0.eps}
\caption{The lump soliton energy $E=\tilde{E}/\tilde{E}_0$ on a unit sphere as a function of $c_0$ for the numerical solutions with $N=-1$ (black dots) and $N=-2$ (red squares) and the analytical solution with $N=1$ (green triangles). 
}
\label{fig-E-c0}
\end{figure}


For a better visualization of the $\phi_a$ fields, we make a projected vector-field plot on a plane. To this end, we first extract the $x$ and $y$ coordinates on the north hemisphere of the real-space $S^2$ and then project the corresponding three-component $\phi_a$ vectors to the $xy$ plane. To make all vectors equally visible, we normalize the projected $\phi_a$ vectors to unit length. The $N=-1$ result is shown in the left panel of Figure \ref{fig-vec}. Taking the south hemisphere produces the same projected vector-field plot. 
If one follows a loop in real space circling the north pole once, the vectors $\hat{\phi}_a$ will make an opposite $2\pi$ rotation, which is the signature of the winding number $N=-1$.

By plugging $n=2$ in Eq.~(\ref{eq-phi}), we find a higher-winding soliton solution. The numerical results of $f(w)$ and $g(w)$ are shown by the dashed lines in Figure \ref{fg}. From the plots, one can see that the $f$ and $g$ curves are steeper around $w=\pm1$.
The numerical value of the winding number of the $n=2$ solution is $N=-1.98$, which is close to $N=-2$ and confirms the solution as a $N=-2$ lump soliton on a sphere. 
The projected vector-field plot of the $N=-2$ lump soliton is shown in the right panel of Fig.~\ref{fig-vec}. One can see the double-wrapping of the vectors around a point, reflecting the $N=-2$ winding number. When comparing Fig.~\ref{fig-vec} with finite $c_0$ to Fig.~\ref{fig-sphere} with infinite $c_0$, one can see their resemblance, despite the opposite directions of the vectors due to the opposite signs of the winding numbers.

In principle, one can obtain lump solitons on a sphere with higher winding numbers by using larger values of $n$ in Eq.~\eqref{eq-phi}. However, the coefficients of Eq.~(\ref{eq-fg}) become more singular for larger values of $n$. Therefore, the higher-$n$ solutions will appear more singular around $w=\pm1$, which then lead to larger numerical errors for those high-winding solutions. 

The energies of the lump solitons are above the ground-state energy because they are topological excitations. Thus, some energy according to Eq.~\eqref{eq-en-1} is required to maintain a soliton configuration.
In Figure \ref{fig-E-c0}, we show the energies as functions of $c_0$ of the $N=-1$ and $N=-2$ lump solutions of Eq.~(\ref{eq-phi}) and also that of the $N=1$ lump solution of Eq.~(\ref{eq-s1}), which is valid only for $c_0>2$. All their energies increase with $c_0$. Both $N=-1$ and $N=-2$ solutions of Eq.~(\ref{eq-phi}) have lower energies than that of the $N=1$ solution of Eq.~(\ref{eq-s1}) when $c_0 >2$. Moreover, the $N=-2$ lump energy from Eq.~(\ref{eq-phi}) is higher than the $N=-1$ lump energy but less than twice of the $N=-1$ energy, which implies that it is not energetically favorable for a $N=-2$ lump to decay into two $N=-1$ lumps. This is again in contrast to quantum vortices with excitation energies proportional to the square of the winding number~\cite{FetterBook,Pethick_book}, which energetically favor vortices with lower winding numbers.


\section{Implications}
The condensate wavefunctions $\phi_a$ determines the corresponding densities, so the lump solitons may be viewed as particular density profiles of the three components. 
To induce a lump soliton, one may need to generate inhomogeneous and non-overlapping density profiles of the three components. Since the lumps are minimal-energy solutions with globally defined winding numbers, they will last longer while other non-soliton modes decay away. 
For example, one may use focused laser lights coupled to different components to induce local perturbations, similar to the methods for measuring critical velocity~\cite{Raman99,Desbuquois12,Weimer15} or exotic sounds (see Ref.~\cite{Hu22} for a review) in cold atoms. It may also be possible to use optical Feshbach resonance~\cite{Fatemi00,Theis04,Enomoto08,Bauer09,Blatt11,YanOFR13} or other means to locally tune the spin-dependent coupling $c_1$. 
Imaging the density profiles of different components of a spinor BEC may be performed by the method in Ref.~\cite{Bersano18}, where non-topological solitons of quasi-1D $F=1$ spinor BECs have been observed.



We remark on extensions of the framework. While the kinetic energy is typically negligible when compared to the interaction energy in bosonic systems, the Fermi energy of fermions may remain competitive in the presence of interactions~\cite{FetterBook}. A curvature-induced BCS-BEC crossover on a spherical shell due to the competition of the kinetic and interaction energies has been discussed~\cite{He22}. Whether the lumps can survive in fermionic systems is an intriguing question. On the other hand, the homotopy $\pi_3 (S^2)=\mathbb{Z}$~\cite{Nakahara} implies a $F=1$ spinor BEC may also exhibit topological behavior if the geometry is $S^3$, which may be possible by engineering the 4D space-time.

\section{Conclusion}
In conclusion, we have shown that spherical lump solitons are realizable a $F=1$ spinor BEC on a spherical shell because the vector field of the condensate wavefunctions wraps the order-parameter space around the real space. The homotopy is characterized by the winding number of the lump solitons. In addition to the exact solutions from spherical harmonics, we have demonstrated an ansatz for generating lump soliton solutions of arbitrary winding numbers in principle. The lump solitons on sphere not only offer concrete examples of topological excitations of interacting quantum systems in curved space but also open a way for investigating nonlinear physics via geometry and topology.

\begin{acknowledgments}
Y. H. was supported by the NNSF of China (No. 11874272) and Science Specialty Program of Sichuan University (No. 2020SCUNL210). C. C. C. was supported by the NSF (No. PHY-2011360).
\end{acknowledgments}


