\section{Introduction}

The concept of "Rubber Duck Debugging" originated from a story in the book \textit{The Pragmatic Programmer} \cite{rubberduckbook}. It encompasses the idea that a rubber duck can assist the process of "debugging" a program by serving as an embodied listener when a person explains the code aloud step-by-step. This concept is related to the technique of Think-Aloud Problem Solving (TAPS) \cite{2010Ku}, which consists of the verbalization of one's thoughts while solving tasks. Prior literature found that verbalization is an effective strategy to improve learning and task performance \cite{1985stinessen,1994michelene}, but this activity takes place as an individual process. Rubber duck debugging reframes the act of problem-solving as a thinking-aloud interaction, where the second interactant is an unresponsive listener.

Social robots can assist people in therapy for rehabilitation \cite{butchart2021child, casas2021social} and mental health \cite{rabbitt2015integrating, kabacinska2021socially}, or support children with special needs \cite{scassellati2012robots, dickstein2018socially, martinez2020socially}. More specifically, the presence of a robot in a TAPS session was seen to improve learning gains in young students \cite{2018ramachandran}. Moreover, work by \citet{Gratch2006} showed that a responsive listener agent, exhibiting non-verbal behavior, leads subjects to spend more time talking and use more words, when compared to an unresponsive listener agent.

In this work, we explored a "rubber duck problem-solving session" and the impact of replacing the inanimate object with an active listener robot.
Prior work has established different methods to generate listener behaviors, formulating rule-based heuristics extracted from human conversation data \cite{2010truong,2021oertel} or using deep learning (DL) models \cite{2019huang,2019wijnen}. We aim to examine if listener behaviors generated through either a heuristic or a DL model lead to different outcomes in a TAPS session. This allows us to explore the sensitivity of the user to robot behavior when the focus of the session is on solving the task rather than the social interaction with the robot.

\begin{figure}
    \centering

    \includegraphics[width=0.37\textwidth]{figures/setup_LIGHTSHORT.png}
        \caption{Overview of the experimental set-up. Participants were invited to think aloud while completing tasks on the tablet. They directed their words to A) a Furhat robot which displays either a heuristic (\textit{NaïveL}) or learned (\textit{DataL}) listening behavior or B) an inanimate object (giant rubber duck).
    }
    \label{fig:setup}
    \vspace{-1em}
\end{figure}

In a between-subject study with three conditions, we evaluated how the two different robot listening behaviors and an inanimate object - a rubber duck - affect the outcome and perception of two distinct think-aloud problem-solving tasks, a deductive logic quiz and an open-ended question. To the best of our knowledge, this is the first study characterizing how the presence of an active listener robot affects the perception of a TAPS session. Considering the reported positive effects of verbalization during problem-solving and the presence of robots in learning scenarios, we explore the potential of robots as "social assistants" in everyday problem-solving tasks. In sum, our main contributions are: (1) investigating the impact of an active listener robot in a thinking-aloud task; (2) developing and evaluating a rule-based listening behavior and a deep-learning based listening behavior system from publicly available human-human dyadic conversational data; (3) exploring how methods and datasets used for developing listener behaviors in one social setting - a conversation - perform in the context of attentive listening to assist a problem-solving session.

\section{Related Work}
\label{sec:relwork}

This research is based on previous work in the areas of social robots in think-aloud tasks and generating active listening behavior.
\begin{description}[align=left, leftmargin=0em, labelsep=0.2em, font=\textbf, itemsep=0em,parsep=0.3em]

\item[Social Robots in Thinking-aloud tasks:] Prior literature has explored how a robot can play an auxiliary role in problem-solving. \citet{2019wijnen} studied how the presence of a social robot can affect explanatory behavior in children, which is associated with deeper learning. The authors concluded that the presence of a robot to which children could direct their words was effective in triggering explanatory behavior (longer duration of utterances, more relevant content). Closer to our work, \citet{2018ramachandran} explored how the presence of a social robot while engaging in a TAPS task could improve learning gains in young students, versus a no-robot condition. The robot would intervene with tutoring notes and reminders to keep thinking aloud. They concluded that the robot platform could support the thinking-aloud process and enhance immediate and persistent learning gains. The authors suggest that the robot represented an embodied social entity towards which students could direct their words, which is in line with the findings from \citet{2019wijnen}. To the best of our knowledge, there is no similar work on adult subjects.

 We expand on prior work by comparing the effects of the presence of an active listener robot to those of an inanimate object in a "rubber duck debugging"/TAPS session with adults. In order to assess if different robot listening behaviors lead to different outcomes, we implemented two types of active listener behaviors, both based on previous work.

\item[Generating Listening Behavior:]

Listening behavior is characterized by backchannels (BCs) \cite{1970yngve}, which signal attentiveness or emotion to a speaker. These include body language, such as head nods or smiles, and short vocalizations (\textit{'uh-huh', 'hmm'}). Backchannel placement is an important component of conversations, as it can skew the perception of the interactants~\cite{2010poppe}.

Significant amount of work has studied the generation of listening behavior in robots.
One approach is to develop rule-based methods that are informed by prosodic features, such as voice pitch and pauses \cite{2010truong,2000ward}, or by behavior observed in human-human conversational data \cite{2021oertel}. As an alternative to designing hand-crafted heuristics, the use of automated methods has been widely adopted\cite{2013dekok, 2016kawahara}. \citet{1996okato} designed a Hidden Markov Model (HMM) based on prosodic patterns to detect when to emit a BC. \citet{2010morency} also used an HMM with multimodal (audiovisual) input features to predict BC timing.

 More recently, deep learning techniques have been deployed to inform the design of listening behaviors. Examples include Recurrent Neural Network (RNN) models, which capture the temporal dependencies of continuous signals (i.e., retain "memory" of previous inputs). \citet{2019ruede} used Long-Short Term Memory (LSTM) layers to predict BCs, with acoustic features and word history as input features. \citet{2019huang} also tested LSTMs and Gated Recurrent Unit (GRU) layers for the development of idle behavior in virtual agents, with multimodal input features. Similarly, \citet{2021murray} made use of LSTMs and multimodal input, and suggested a method for data augmentation that positively impacts BC prediction.
 While the use of RNNs has been widely adopted in sequential decision-making tasks, such as the generation of listener behavior, we highlight other interesting approaches, such as semi-supervised learning \cite{jain2021exploring}, or reinforcement learning \cite{2019hussain,2019rudovic}.

In the present work, we developed two different robot listener behavior generation policies - one rule-based, one DL-based - and evaluated how these impacted various aspects of think-aloud problem-solving tasks: engagement with thinking-aloud, task performance, and user experience.

\end{description}

\section{Implementation of Listener Behavior} \label{sec:listener}

We set out to study how different listening behaviors can support a "rubber duck" problem-solving session. The listening behavior consisted of the robot's \textbf{gaze} and \textbf{backchanneling} (vocal or non-vocal). For \textit{vocal} utterances, the specific sound was selected randomly from a set of pre-defined utterances (\textit{e.g., 'hmm', 'ahh'}). The \textit{non-vocal} backchannel was realized through a head nod.

For developing backchanneling behavior models, we included two different approaches, as we intended to shed light the sensitivity of the user to listening behaviors when the interaction with the robot is only \textit{auxiliary} (not \textit{collaborative}) to task completion.
We tested each of the two models as a different robot listening condition in our user study.

As suggested by prior work~\cite{jain2021exploring}, the decision-making process for generating a backchannel was divided into two stages: first, the \textit{timing} of the backchannel was determined; in a second step, the \textit{type} of backchannel was selected. This separation allows for handling the natural imbalance in conversational data (few examples of each backchannel type and more datapoints for "no backchannel").

Backchannel timing and type were either heuristically determined (naïve listener model, \textit{NaiveL}) or learned through the application of a deep learning model (data-driven model, \textit{DataL}). We detail the implementation of these behaviors below.

\subsection{Naïve Listener}
\label{sec:naivel}

To develop the robot's behavior for the naïve listener condition, we decided to use a frequently deployed heuristic, originally developed by \citet{2000ward}. This work used a corpus of human conversational data in English to implement a heuristic that finds the most appropriate \textit{timing} for a backchannel in a dyadic interaction. The heuristic is based on prosodic cues, i.e. pitch. The authors define a set of conditions for producing a backchannel, which we reproduce \textit{verbatim}:
\textit{\begin{itemize}
    \item[\textbf{(P1)}] a region of pitch less than the 26th-percentile pitch level and
    \item [\textbf{(P2)}] continuing for at least 110 milliseconds,
    \item [\textbf{(P3)}] coming after at least 700 milliseconds of speech,
    \item [\textbf{(P4)}] providing you have not output backchannel feedback within the preceding 800 milliseconds,
    \item [\textbf{(P5)}] after 700 milliseconds wait.
\end{itemize}
}

To run this model in our user study, we computed the pitch distribution of the user's voice over the last $50$ s of audio and calculated the pitch region percentile over the last $10$ ms.
The backchannel \textit{type} was randomly selected to be either vocal or non-vocal when the five conditions shown above were met.

\subsection{Data-Driven Listener}
\label{sec:datal}
We explored the use of machine learning to determine backchannel timing and type. While we maintain the input modality in the form of audio features, as in the naïve condition, we wanted to explore if the learned behavior, by having higher potential to adapt to different speakers, could lead to more positively evaluated interactions \cite{2010morency}.

\begin{table*}[]
\centering
\tabcolsep=0.11cm
\caption{Performance metrics on model used (averaged over all validation folds). Both models used augmented data for training. Comparison with performance of the rule-based model (\textit{NaïveL} condition) on the same dataset.}
\label{tab:modelperf}
\begin{tabular}{c|c|c|c}
\hline
\textbf{Module}                & \textbf{Model}                & \textbf{Hyperparameters}                                                                                & \textbf{Performance}                                         \\ \hline
\multirow{3}{*}{\textbf{Timing}} & \multirow{3}{*}{\textit{GRU}} & \multirow{3}{*}{\textit{\begin{tabular}[c]{@{}c@{}}Lookback: 5,activation: sigmoid, \\ batch size: 16,  dropout: 0.0, \\ loss function: focal,optimizer: Adam\end{tabular}}} & \textbf{Macro Accuracy}: 0.95, \textbf{Precision}: 0.52, \textbf{Recall}: 0.51, \textbf{F1}:0.50  \\
   &    &   & \textbf{Margin Accuracy}: 0.95, \textbf{Precision}: 0.59, \textbf{Recall}: 0.76,  \textbf{F1}:0.65 \\
                               &                               &                                                                                & \textbf{BC Prediction Deviation}: 0.83                                       \\ \hline
\multirow{3}{*}{\textbf{Type}} & \multirow{3}{*}{\textit{GRU}} & \multirow{3}{*}{\textit{\begin{tabular}[c]{@{}c@{}}Lookback: 10, activation: sigmoid, \\ batch size: 32, dropout: 0.2,\\  loss function: MSE, optimizer: SGD\end{tabular}}}   & \textbf{Macro Accuracy}: 0.64, \textbf{Precision}: 0.37, \\    &    &         &       \textbf{Recall}: 0.39, \textbf{F1}:0.35                                                         \\
                               &                               &                                                                                &      \\  \hline \hline
\multirow{2}{*}{\textbf{Timing}} & \multirow{2}{*}{\textit{NaïveL}} & \multirow{2}{*}{\citet{2000ward}}   & \textbf{Macro Accuracy}: 0.95, \textbf{Precision}: 0.48, \\    &    &         &       \textbf{Recall}: 0.50, \textbf{F1}:0.49
\end{tabular}
\end{table*}

\subsubsection{Problem Formulation}
We formulated the problem of generating backchanneling in a TAPS session as a sequential decision-making problem. At any time-step $t$, the environment (user's voice features) is captured as state variable $s_t \in S$. The robot takes actions $a_t \in A$. It first chooses between a binary output: \textit{performing BC} or \textit{doing nothing}. If the robot chooses \textit{performing BC}, $s_t$ is again used to inform \textit{which type} of BC to perform. The potential actions $a_t$ that follow are \textit{vocal BC} (vocal utterance), \textit{non-vocal BC} (nodding), or \textit{both}.

\begin{description}[align=left, leftmargin=0em, labelsep=0.2em, font=\textbf, itemsep=0em,parsep=0.3em]
\item[Dataset:]
The collection of conversational data is frequently the first step to developing listening behaviors. In this work, we wanted to test if we could instead use data collected in a different setting from that in which we deployed our model. Some interesting corpora of dyadic interactions in English include CCDb \cite{ccdb}, IEMOCAP \cite{iemocap} and D64 \cite{D64}, among others \cite{2017corpora}. We were looking for unscripted, non-topic-bounded interactions that were publicly available. As per these criteria, the dataset used for training was Cardiff's Conversation Database (CCDb)\cite{ccdb}.

The CCDb consists of 30 dyadic conversations, each lasting around five minutes. The 30 conversations are between 16 different speakers (12 M, 4 F), with ages ranging from 25-56 years old. We used the subset of conversations that were annotated for facial expressions and utterances to train the listening behavior policy. We extracted data from the perspective of each participant, totaling 80 minutes of conversational data.

We split the individual feature streams into moments where the participant was the \textit{speaker} or the \textit{listener}. As we are focusing on listener behavior, we extracted features only in the \textit{listener} portions of the conversation. As training data, we used the audio features from the \textit{speaker} (other participant), while extracting the annotated backchannels performed by the \textit{listener} participant as ground truth for our model. The conversations from the CCDb dataset include extensive annotations. We filtered the annotations to provide positive training instances in moments where vocal backchannels or head motion (nodding) were present.

\item[State space:]

As input features for our model, we extracted speech features from the speaker, including 13-dimensional mel-frequency cepstrum coefficients (MFCC) and 4-dimensional prosody features as used in prior works \cite{2019ruede,2021murray,jain2021exploring}. The MFCC features are computed every 30 ms with a sliding hamming window of 400 ms. Prosody features include pitch (fundamental frequency) and yin-energy, as well as the first derivative of these variables. The final 34-item feature vector is composed of the mean and standard deviation of each of these features. The dataset frequency is 2 Hz (one sample every 0.5 s) and was normalized before training.

\item[Action Space:]
The robot could perform the same actions as those performed by the \textit{NaïveL} robot (Sec. \ref{sec:naivel}) - vocal utterance, nod of varying amplitude, or a simultaneous combination of both.
\end{description}

\subsubsection{Training policies for listening behavior}
To train sufficiently good behavior policies, we explored different model architectures and techniques such as data augmentation\cite{2021murray}. Below, we explain the models tested and the selection process, as well as model performance.

\begin{description}[align=left, leftmargin=0em, labelsep=0.2em, font=\textbf, itemsep=0em,parsep=0.3em]

\item[Data Augmentation:]

\citet{2021murray} suggest a method to improve the robustness of robot listening behaviors by using audio data augmentation. The authors show that the model trained on this data outperforms a rule-based and a random model as reported by users. We emulated the proposed method by making use of \textit{masking} techniques in the time and frequency domains. In our work, training instances (audio features) were partially masked in one or both domains, chosen at random. Both the original and the deformed instances were used for training.

\item[Models:]
We explored the use of deep learning techniques to develop listening behavior. We separated our model into two components according to the two-staged decision-making process, the \textbf{Timing} and the \textbf{Type} modules. The \textbf{Timing} module learns the appropriate timing to perform backchanneling, and the \textbf{Type} module decides which type of feedback - vocal, non-vocal (nodding) or both - the robot should emit. The first module is a binary classification task (\textit{do BC} or \textit{do nothing}). In the case of a positive (\textit{'do BC'}) output, the second module is activated (a multiclass classification model).

The use of neural networks has proven to be a versatile tool to learn relevant features automatically \cite{2015mueller}. Multiple authors have approached this problem by making use of models that consider previous internal states. LSTMs are a common approach \cite{2019ruede,2021murray}, although these architectures are prone to overfitting. Alternatively, Gated Recurrent Units (GRUs) can be used \cite{2019huang}. These are less complex and usually preferred over LSTMs for small training datasets.

\item[Model Hyperparameter Tuning:]
In order to determine the best model, we tested different combinations of models and parameters. For both the \textit{Timing} and \textit{Type} modules, we trained single-layer LSTM and GRU models, which were followed by a dropout layer and a final dense layer. Each model was either trained on the original or the augmented dataset, considering varying previous timesteps, i.e. lookback sizes (5, 10 and 15, respectively 2.5, 5 or 7.5 seconds of data), activation functions (sigmoid, ReLU and softmax), batch sizes (8, 16, 32), dropout rates (0, 0.2 and 0.4), L2 regularization (0.1, 0.01, 0.001, 0.0001), loss functions (focal loss, MSE, binary cross-entropy or hinge loss) and optimizers (SGD and Adam).

\item[Evaluation Metrics:]

Performance on both modules was evaluated based on macro \textbf{accuracy} (number of correct predictions over all predictions), \textbf{precision} (probability that emitted BCs match the ground-truth data), \textbf{recall} (probability that a ground-truth BC is predicted by the model) and \textbf{F1-score} (harmonic mean of precision and recall).

Additionally, for the \textbf{Timing} module and in line with similar work \cite{2012dekok,2019ruede}, we calculated the above metrics considering a tolerance margin of $[-500, 500]$ ms; finally, we controlled for robot overreaction with a \textbf{Backchannel frequency deviation} metric. For each participant $i$, we considered the relative difference between the predicted number of backchannels and its original value:
\begin{equation}
\Delta{BC}^{i}_{dev} = |Y^i_{true} - Y^i_{pred}| / Y^i_{true}
\end{equation}

\item[Models' Selection and Performance:]
For training, data was split into 8 folds, with non-overlapping participant data in a 6:1:1 train-validation-test split. Each model was trained for a maximum of 100 epochs, but early stopping strategies were deployed for both loss and validation loss, in order to prevent overfitting. The use of a dropout layer and L2 regularization entailed the same purpose.

Each module (\textit{Timing} and \textit{Type}) was trained with an analogous methodology. We fine-tuned hyperparameters based on macro accuracy and F1-score. The final candidates for each model type (LSTM or GRU, augmented and not augmented training data) were then trained using k-fold cross-validation.

For each module, we deployed the model that was most frequently in the top three best performing models for each of the metrics mentioned above. Both models are single-layer GRUs trained on the augmented dataset (see performance on Table \ref{tab:modelperf}, along with performance for the \textit{NaïveL} heuristic on the same dataset). The \textbf{Timing} model outperforms other similar model architectures for BC prediction \cite{2019ruede,2021murray}. Interestingly, the modules have different lookback periods, with the \textbf{Type} model input considering a higher number of past internal states.
\end{description}

\subsection{Gaze Behavior}

Gaze behavior was kept constant across the two listening conditions, as direct eye contact was predicted to occur infrequently during the problem-solving sessions. We implemented the robot's gaze behavior based on the work by \citet{2001garau}. They found that adapting an avatar's gaze and head pose to the turn-taking pattern of a conversation outperformed a random gaze agent. The agent alternates between two actions, looking \textit{at} the user or \textit{away} from them (gaze aversion). The duration of each action is informed by previous work on dyadic interactions \cite{1976argyle,1967kendon} and varies depending on whether the agent is listening or speaking. For this work, we only used the listening timing indications.

\section{Study Design}\label{sec:problem}

We designed a between-subject study with three conditions in order to evaluate the effects of different previously proposed robot listener behaviors in the outcome of a TAPS activity. Two conditions concerned robot listener behaviors: the first, a rule-based heuristic where non-verbal utterances and head nods are crafted according to previous literature (as described in Section \ref{sec:naivel}); the second, based on deep learning models which learn listener behaviors through data input from human-human dyadic interactions (as described in Section \ref{sec:datal}). The two listening behaviors were chosen to gauge if different complexity in the model used to generate the listener behavior affects the outcomes of TAPS activities. We compare these effects with the baseline condition, the presence of an inanimate object - a rubber duck - to whom words are directed. The three conditions are detailed in Section \ref{subsec:cond}.

\subsection{Hypotheses}

\label{subsec:hyp}

We were interested in evaluating how a robot can affect the think-aloud process in a "rubber duck debugging" session. Thus, we defined our outcomes based on how the user interacted with the task. We measure user \textit{engagement} in thinking-aloud and user \textit{perception} of the problem-solving session (or user experience). \textbf{Engagement} with the think-aloud process measures how the user talks when thinking aloud, e.g. number and duration of verbal utterances.
\textbf{User experience} addresses how the participant perceives solving the task, including self-assessed cognitive load (or "mental effort") and enjoyment (or "user experience index"). These measures are detailed in Section \ref{subsec:metrics}.

\citet{2019wijnen} found that the presence of a robot leads to better explaining behavior while learning in children, and \citet{2018ramachandran} found that children engage more in thinking aloud utterances when a robot is present. Consequently, we hypothesized:
\begin{itemize}
	\item \textbf{H1a:} \textit{Engagement in thinking aloud is higher when thinking aloud to an actively listening robot than to an inanimate object.}
\end{itemize}
Moreover, \citet{2010morency} report that a multimodal probabilistic model is better at predicting backchannels than hand-crafted rules. As such, we expected the data-driven condition to generate listener behavior that leads to a more organic interaction:
\begin{itemize}
	\item \textbf{H1b:} \textit{Listener behaviors learned from data lead to higher degrees of engagement in thinking-aloud, when compared to rule-based heuristics.}
\end{itemize}

Following \textbf{H1a, H1b}, and because the verbalization of one's thoughts during problem-solving has been reported to have affect how users solve a task \cite{1985stinessen,1994michelene}, we also hypothesize about user behavior during the task:
\begin{itemize}

    \item \textbf{H2a:} \textit{User behavior during task completion will be different in the conditions when a robot is present. Users will take longer to solve tasks, but also perform better.}
    \item \textbf{H2b:} \textit{Following \textbf{H1b}, users with the DL-based listener robot condition will take a longer time completing the tasks than users in the rule-based listener robot condition.}

\end{itemize}

Finally, we wanted to evaluate the subjective experience of the user during the problem-solving sessions. Literature is sparse regarding if thinking aloud affects how users perceive the task, and if different listening behaviors skew that perception. As such, we took an exploratory approach to answer the following research questions:
\begin{itemize}
	\item \textbf{RQ1a:} \textit{Is user experience affected by an active listener robot when compared to an inanimate object?}
	\item \textbf{RQ1b:} \textit{Is user experience affected by the type of listening behavior the robot exhibits?}
\end{itemize}

\subsection{Conditions}
\label{subsec:cond}

We sought to evaluate how the presence and behavior of a social robot impacted the outcome of a "rubber duck" problem-solving session. We designed a between-subject study with three conditions:
\begin{description}[align=left, leftmargin=0em, labelsep=0.2em, font=\textbf, itemsep=0em,parsep=0.3em]
\item[Rubber Duck (RDuck):] To reproduce a "rubber duck debugging" scenario, participants in the \textit{RDuck} condition were instructed to think aloud directing their words to a rubber duck.
\item[Naïve Listener Robot (NaiveL):] The social robot displayed listening behavior generated by hand-crafted heuristics. Information on the detailed implementation is given in Section \ref{sec:naivel}.
\item[Data-driven Listener Robot (DataL):] The social robot displayed listening behavior that was learned through machine learning methods using a human-human conversational corpus. The implementation details are provided in Section \ref{sec:datal}.

By including two different robot listening behaviors, we explored the transferability of two implementation methods of different complexity into a new social context, as well as the sensitivity of the user to differences between these behaviors.

\subsection{Problem-solving Tasks}
\label{subsec:task}

We designed two different tasks in order to evaluate distinct contexts for thinking-aloud problem-solving and how the presence of a robot can impact the perception of those tasks and the outcomes.

\subsubsection{\textbf{Logic Deduction Quiz (LogicQ)}} Tasks using deductive logic are common when developing thinking-aloud protocols \cite{2009leighton}. The questionnaire used here was adapted from \textit{IQ Test Labs}\footnote{https://www.intelligencetest.com/questions/logic/1.html}. For each question, 5 possible answers were shown as illustrated in Figure \ref{fig:logicd}. Participants were able to skip questions but could not see if their answers were correct. There was no time limit per question.

\subsubsection{\textbf{NASA Moon Survival Task (NASA)}}  In the NASA task \cite{1970hall}, the participant is placed in an hypothetical life-threatening scenario (spaceship crashed on the moon) and asked to rank a list of 15 available items according to their contribution towards survival.

The LogicQ task allowed for an objective evaluation of performance by counting correct answers to the quiz, whereas the NASA task was framed as an open-ended question involving strategic thinking. We analyzed the effects of each task separately, as they elicit different cognitive processes, and tasks similar to these give rise to different protocols in thinking-aloud studies \cite{1994vansomeren}.

\begin{figure}
    \centering

    \includegraphics[width=0.25\textwidth]{figures/q1.png}
        \caption{Layout example for one of the questions from the deductive logic quiz, from \textit{IQ Test Labs}.
    }
    \label{fig:logicd}

\end{figure}
\subsection{Measures}
\label{subsec:metrics}
We established a set of metrics that explore the outcomes of the "rubber duck" problem-solving session. \textbf{Engagement} in thinking-aloud and \textbf{Task-related} metrics quantify the objective effects of the different conditions on how the user acted, whereas \textbf{User experience} attempts to quantify the subjective experience of the user while solving the tasks.

\subsubsection{Engagement metrics} To assess \textbf{\textit{H1}}, we measured the user's \textbf{Number} and \textbf{Duration of verbal utterances}, as well as \textbf{Speech-to-silence ratio} (time spent talking divided by time spent silent, during a task) while solving each task. This provided an objective measure of the user's engagement with the think-aloud process. This data was collected from Voice Activation Detection (VAD) through lapel microphones.

\subsubsection{Task-related metrics} For evaluation of \textbf{\textit{H2}}, we defined metrics that relate to how the user behaved while solving the task. We measured performance in the \textit{LogicQ} task by monitoring the number of \textbf{Questions answered (QAns)}, and number of \textbf{Correct answers}. For the \textit{NASA} task, we considered \textbf{time spent} before submitting the final answer.

\subsubsection{User experience metrics} Finally, to investigate \textbf{\textit{RQ1}}, we measured self-reported \textbf{Cognitive load} and \textbf{User experience index (UEI)} of the problem solving session. This was assessed in a questionnaire taken \textit{after each task}. Cognitive load was measured using the NASA-TLX scale \cite{1988nasatlx}. The perception of the problem-solving session was evaluated with three dimensions from the short version of the User Engagement Scale (focused attention, perceived usability and reward) \cite{2018obrien}.

\subsubsection{Other measures}
We also collected data on participants' demographics (age, gender identity, country of origin), as well as the personality trait of \textbf{Extraversion} from the Big Five Inventory\cite{john1991big,2007rammstedt} in the pre-experiment questionnaire. In the post-experiment questionnaire, we assessed the \textbf{Impact of thinking aloud}, which asked users to rate how they thought thinking aloud helped them solve the tasks. This measure was obtained through five questions with a 5-point Likert scale. A Cronbach's Alpha of $\alpha=0.84$ indicated good reliability. The exact questions are provided as supplemental material. We further measured \textbf{Extraversion} and \textbf{Impact of thinking aloud} since \textit{personality} and \textit{natural inclination for speaking aloud} can impact user behavior during the tasks. Finally, we considered the \textbf{Order} of task presentation, since thinking aloud can suffer from a "cold start effect" \cite{1997gibson} (as described in Section \ref{subsec:exp}).

Finally, to explore if differences in user behavior could be due to the perception of the robot listening behaviors, the participants with conditions where the robot was present were asked about the robot's \textbf{Listening behavior} and \textbf{Closeness} \cite{2021murray}, \textbf{Social attributes} (from RoSAs \cite{rosas}), as well as previous experience with robots.

\subsection{Experimental procedure}
\label{subsec:exp}

Participants were randomly assigned to one of the three conditions. After giving informed consent, participants read the experiment instructions. In all conditions, they were instructed to speak aloud to "Alex", their "rubber duck", in a study that aims to understand cognitive processes behind problem-solving tasks. The instructions also clarified that Alex did not have any insight into the solutions of any of the two tasks. Afterwards, participants were asked to fill the pre-experiment questionnaire (demographic and personality information). Before moving to the experiment room, participants were given 2 minutes to respond to questions from the LogicQ task, in order to get familiar with the types of questions presented in the quiz.

Participants then moved to the experimental room. They sat at the desk, facing a tablet and "Alex". "Alex" was either the rubber duck or a Furhat robot, depending on the assigned condition. Figure \ref{fig:setup} illustrates the set-up. Audio features and voice activity were collected through a lapel microphone. Before starting the first task, participants were introduced to "Alex" and were asked to introduce themselves to it. This was prompted by the robot (or by the researcher in the \textit{RDuck} condition) and aimed to a) eliminate the "cold start effect' of talking aloud \cite{1997gibson} and, for the \textit{NaïveL} and \textit{DataL} conditions, b) eliminate the novelty effect \cite{2016kennedy} when first speaking to the robot. Participants were given $7$ minutes to complete each task (for the NASA task, they were free to submit their answer before that time). Between the two tasks, they were asked to fill out a questionnaire about the task they had just solved (collection of user experience measures). Task order was counterbalanced between participants, within each condition. After the second task, participants filled out the second task-based questionnaire and a final questionnaire (collection of perception of the robot and impact of thinking aloud measures).
After a final debriefing participants received a voucher valued at around 10 USD.

\paragraph{\textbf{System implementation}} As our social robot, we used a Furhat robot\footnote{https://furhatrobotics.com/} with the William voice from CereProc\footnote{https://cereproc.com/}, including its set of backchannels. For the nod movement, the "amplitude" - range of the up and down movement - was randomly sampled from a uniform distribution. During the nod, the robot paused at the lowest point for $0.5$ s.  Furhat was controlled through a computer using a NVIDIA GeForce RTX 2080 SUPER and an Intel® Core™ i9.

\paragraph{\textbf{Participants}}

Participants were recruited through posters, flyers, social media platforms and word of mouth. A total of 101 participants were recruited with ages ranging from 14-76 years (M = 26.4, SD = 7.6). 53 participants identified as male and 48 as female, with a total of 29 different nationalities.
\end{description}

\section{Results} \label{sec:results}

Data from eleven participants had to be excluded due to substantial software or hardware problems. The remaining 90 participants were distributed across conditions with the following demographics: $RDuck - N = 30$ $(14F,16M)$, ages $ 27.8 \pm 10.2$; $NaiveL - N = 29$ $(9F,20M)$, ages $26.3 \pm 7.4$; $DataL - N = 31$ $(18F,13M)$, ages $25.3 \pm 5.3$. Out of the 60 participants which interacted with the robot, 20 indicated they had interacted with a robot before. In this section, we analyze the listening behaviors developed and report our findings for each measure presented in Section \ref{sec:problem}.

\subsection{User Engagement in Thinking Aloud}

We evaluated how much participants engaged in thinking aloud utterances in two separate tasks with an ANCOVA to examine the effects of \textit{condition} on the \textbf{Number} and \textbf{Duration of utterances}, as well as \textbf{Speech-to-silence ratio}, after controlling for \textbf{Extraversion} and \textbf{Impact of Thinking Aloud}. We also control for the \textbf{Order} of presentation of tasks, due to the "cold start effect" mentioned in Section \ref{sec:problem} \cite{1997gibson}. Table \ref{tab:vad} summarizes our findings.

\subsubsection{NASA task}

The effect of the \textit{condition} on the \textbf{Number of utterances} was not significant ($F(2,83)=0.33,p=0.71$). The covariate \textbf{Order} was a significant predictor for the \textbf{Duration of utterances} and \textbf{Speech-to-silence ratio} ( $F(1,83)=4.16,p=0.04$ and $F(1,83)=5.84, p=0.02$, resp.). Doing the NASA task first decreased mean \textbf{duration of utterances} ($3.30\pm1.31 s$ vs. $3.72\pm1.85 s$ if doing the LogicQ task first) and \textbf{speech-to-silence ratio} ($1.46\pm1.01$ vs. $1.88\pm1.74$).

\subsubsection{LogicQ task}
The ANCOVA test showed no significant effects of the \textit{condition} on the engagement of the user with the task (\textbf{Number of utterances}, $F(2,83)=0.21,p=0.81$, \textbf{Duration of utterances}, $F(2,83)=1.15,p=0.32$, \textbf{Speech-to-silence ratio}, $F(2,83)=1.18,p=0.31$).

\begin{table
\centering
\tabcolsep=0.11cm
\caption{User Engagement metrics, for both tasks ($M\pm SD$).}
\label{tab:vad}
\begin{tabular}{c|c|c|c|c}
\hline
\textbf{Task}                   & \textbf{Condition} & \textit{\textbf{Number}} & \textit{\textbf{Duration}} & \textit{\textbf{Speech-to-silence}} \\ \hline
\multirow{3}{*}{\textbf{NASA}}  & \textit{RDuck}     & $0.17\pm0.04$            & $3.50\pm1.34$              & $1.73\pm1.12$                       \\
                                & \textit{NaiveL}    & $0.16\pm0.04$            & $3.70\pm1.46$              & $1.73\pm1.19$                       \\
                                & \textit{DataL}     & $0.16\pm0.04$            & $4.37\pm2.66$              & $2.55\pm2.44$                       \\ \hline
\multirow{3}{*}{\textbf{LogicQ}} & \textit{RDuck}     & $0.16\pm0.05$            & $2.90\pm0.87$              & $1.21\pm0.78$                       \\
                                & \textit{NaiveL}    & $0.16\pm0.04$            & $3.28\pm1.32$              & $1.27\pm0.85$                       \\
                                & \textit{DataL}     & $0.17\pm0.03$            & $3.31\pm1.10$              & $1.52\pm1.11$                       \\
\end{tabular}
\end{table}

\subsection{Task-related metrics}

Table \ref{tab:performance} shows task-related measures across participants within each condition. For the NASA task, we consider \textbf{Time taken} until submission of final object ranking (in seconds, with a maximum of $420 s$). For the LogicQ task, we looked at the \textbf{Number of Questions Answered (QAns)} and \textbf{Number of Correct Answers} as a way to measure how the different conditions affected how the user completes the LogicQ task. We controlled for the \textbf{Impact of Thinking Aloud} in both tasks.

\subsubsection{NASA task}

An ANCOVA revealed no significant effect of \textit{condition} on \textbf{Time Taken} solving the task ($F(2,84) = 1.67, p=0.19$).

\subsubsection{LogicQ task}
A one-way ANOVA revealed no significant effect of \textit{condition} on \textbf{Correct answers} ($F(2,84) = 0.21, p=0.81$), or on \textbf{Number of questions answered} ($F(2,84) = 0.14, p=0.86$).

\begin{table
\centering
\tabcolsep=0.11cm
\caption{Task Performance metrics, for both tasks ($M\pm SD$).}
\label{tab:performance}
\begin{tabular}{c|c|c|c|c}
\hline
\textbf{Task}                   & \textbf{Condition} & \textit{\textbf{QAns}} & \textit{\textbf{Correct}} & \textit{\textbf{Time Taken (s)}} \\ \hline
\multirow{3}{*}{\textbf{NASA}}  & \textit{RDuck}     & -                      & -                             & $360\pm75$                       \\
                                & \textit{NaiveL}    & -                      & -                             & $349\pm80$                       \\
                                & \textit{DataL}     & -                      & -                             & $320\pm99$                       \\ \hline
\multirow{3}{*}{\textbf{LogicQ}} & \textit{RDuck}     & $8.41\pm2.93$         & $0.011\pm0.005$                & -                                \\
                                & \textit{NaiveL}    & $8.93\pm3.50$         & $0.010\pm0.36$                & -                                \\
                                & \textit{DataL}     & $8.58\pm3.50$         & $0.009\pm0.005$                & -
\end{tabular}
\end{table}

\subsection{User Experience metrics}
\begin{table
\centering
\tabcolsep=0.11cm
\caption{User Experience metrics, for both tasks ($M\pm SD$).}
\label{tab:ues}
\begin{tabular}{c|c|c|c}
\hline
\textbf{Task}                   & \textbf{Condition} & \textit{\textbf{Cog. Load}} & \multicolumn{1}{c}{\textit{\textbf{UEI}}} \\ \hline
\multirow{3}{*}{\textbf{NASA}}  & \textit{RDuck}     & $3.47\pm0.87$               & $3.42\pm0.31$                               \\
                                & \textit{NaiveL}    & $3.53\pm0.73$               & $3.39\pm0.31$                              \\
                                & \textit{DataL}     & $3.11\pm0.67$               & $3.40\pm0.33$                              \\ \hline
\multirow{3}{*}{\textbf{LogicQ}} & \textit{RDuck}     & $4.18\pm0.88$               & $3.40\pm0.29$                               \\
                                & \textit{NaiveL}    & $4.08\pm0.81$               & $3.26\pm0.30$                              \\
                                & \textit{DataL}     & $4.11\pm0.84$               & $3.27\pm0.34$
\end{tabular}
\end{table}
Participants' self-reported experience solving the tasks was evaluated through measures of \textbf{Cognitive load} and \textbf{User Experience Index} ($UEI$). We tested the effect of condition on these metrics with an ANCOVA, controlling for the \textbf{Order} of tasks and \textbf{Impact of Thinking Aloud}. See Table \ref{tab:ues} for more detail.

\subsubsection{NASA task}
The effect of the \textit{condition} on self-reported \textbf{Cognitive load} is not significant, $F(2,84)=2.43, p=0.09$.
The covariate \textbf{Order} of task presentation significantly effected (decreased) \textbf{UEI}, $F(1,84)=8.11,p=0.005$ (NASA first, \textbf{UEI} $=3.32\pm0.31$; LogicQ first, \textbf{UEI} $=3.39\pm0.32$).

\subsubsection{LogicQ task}

A one-way ANCOVA showed no significant effects of the \textit{condition} on User Experience metrics (\textbf{Cognitive Load}, $F(2,86)=0.05,p=0.95$, \textbf{UEI}, $F(2,86)=2.05,p=0.14$).

\input{figures/testfigure}

\subsection{Robot Behavior}

In addition to the main measures, we also looked into how the robot acted and how it was perceived in the \textit{NaïveL} and \textit{DataL} conditions. This helps to illustrate the two listener behavior models implemented, but can also correlate with user behavior during the TAPS sessions under different conditions.

Backchannel frequency per minute (see Fig. \ref{fig:robotbehavior}) is given by the count of backchannels executed by the robot every 60 seconds. It is calculated as a sliding window with hop length of $15 s$.
The \textit{NaïveL} robot displayed similar values for vocal and non-vocal backchanneling, but the frequency decreased as task time increases. The \textit{DataL} robot, on the other hand, favors non-vocal backchannels (nodding), which is displayed in much higher frequency, as well as a combination of both nodding and vocal utterances.

\subsubsection{User Reporting of Robot Behavior}

\input{figures/test2figure}

Fig. \ref{fig:robotbehavior_listening} shows participants' evaluation of the robot's listening behavior. A one-way ANOVA showed no significant differences between conditions for both the \textbf{Closeness} and \textbf{Listening behavior} dimensions ($F(1,58)=2.07,p=0.16$ and $F(1,58)=0.56, p=0.46$, respectively). The robot's social attributes \cite{rosas} were also not significantly different between the \textit{NaiveL} and \textit{DataL} conditions (one-way ANOVA for the \textbf{Competence} dimension, $F(1,58)=0.46,p=0.50$, and Wilcoxon rank sum test for the \textbf{Warmth}, $W=543.5,p=0.10$ and \textbf{Discomfort}, $W=410,p=0.55$, dimensions, after a Shapiro–Wilk test of normality revealed these variables did not have normal distributions, $p<0.05$).

\section{Discussion} \label{sec:discussion}

In order to capture the impact of a "robot duck," we analyzed various aspects of a thinking-aloud problem-solving session. Overall, we cannot say that the presence of the robot impacted the "rubber duck" TAPS sessions, nor that the two listening behavior models caused significant differences in task performance or user experience. Our results reject \textbf{H1a,b} - the presence of the robot, in either condition (\textit{NaïveL} or \textit{DataL}), did not appear to affect user engagement in thinking aloud. This result is in contrast to the findings in \citet{2018ramachandran}, where children thinking aloud to a robot spoke more than those in a no-robot condition. We also reject \textbf{H2a,b} - task performance was not affected, and the time spent to complete the tasks did not significantly change regardless of the condition.
Finally, \textbf{RQ1} was aimed at the impact of the presence of the robot on subjective user experience. \citet{2013jung} found that users in human-robot teams with robots that backchannel show lower cognitive load when solving tasks. However, we found no significant effects of condition on how users self-reported the sessions. Below, we advance potential explanations for the findings in this study.

The data-driven model appeared to greatly favor non-verbal feedback, with rare instances of vocal-only backchannels.
However, notably, the users rank both behaviors similarly in the dimensions of \textbf{Listening behavior}, \textbf{Closeness} and social attributes (\textbf{Competence, Warmth} and \textbf{Discomfort}). These findings differ from those reported in \citet{2021murray}, which found that a data-driven listener behavior was perceived as a better listener than a rule-based model. This might indicate that the user becomes fully involved with the task, paying little attention to the robot. This phenomenon has been observed in a setting where children playing a game stop listening to the robot's suggestions after some time \cite{2022ali}.

High user involvement in task completion could also explain why we do not find significant differences between the \textit{RDuck} condition and the two robot conditions. Relevant prior work on human dishonesty in the presence of a robot \cite{2022petisca} found that there were no differences in how users cheated between being alone or in the presence of an unaware robot. We hypothesize that there was a mismatch between the expectations of the users and robot behavior, leading users to believe that the robot did not truly understand what they were saying. Furhat's anthropomorphic appearance elicits the development of expectations about the system's performance \cite{2000nass}, causing cognitive anthropomorphization \cite{2022sacino}. Previous studies also showed different reactions to robot failure depending on the robot's appearance \cite{2020kontogiorgos} and functionality framing \cite{2020washburn}.
Inadequate backchannel timing thus may have had a detrimental effect on the expectations of the user, leading them to believe Furhat was unaware of the meaning of their utterances - and diluting the differences between the three conditions tested. This is consistent with comments left by the users in the post-experiment questionnaire, e.g., \textit{" (...) the robot is responsive however, can't replace a human interaction as I didn't get the emotional feedback.”}(P72).

The order of presentation of tasks is a significant predictor of measures such as \textbf{Speech-to-silence ratio} and self-reported \textbf{User Experience}. As mentioned in Section \ref{sec:problem}, an initial inhibition is to be expected when thinking aloud \cite{1997gibson}. Further, some participants reported that the two tasks differed in how organic it felt to say their thoughts aloud, with many feeling that it was harder to do it during \textit{LogicQ} quiz. Deductive logic skills also vary among humans, and this may explain why the \textbf{Impact of thinking aloud} is a predictor of involvement measures for the \textit{LogicQ} task and not the \textit{NASA} task.

On a final note, we observe that much work in the development of robot listener behavior requires the collection of specific conversational data in chosen human contexts that are then replicated in human-robot interactions \cite{2010morency, 2019huang,2021murray}. Given the wide range of social settings and contextual characteristics within which HRI can be applied, one might wonder: is it sustainable to continue leveraging social HRI research by collecting task-specific datasets? In this work, we developed two fully-autonomous listening behaviors models that were informed by human-human conversational data collected in a contrasting social setting to that in which they were deployed. In spite of this, both behaviors lead to successful interactions of the user with the robot, with no significant detriment in any of the aspects of the sessions when compared to an inactive listener. This opens the door to continue exploring transferability of datasets to new social settings.

\subsection{Limitations and Future Work}
\label{sec:limit}

This study carries some limitations. Functionality misframing, or mismatch of user expectations towards the robot, could be reduced by providing participants with more interaction time with the robot before task completion, in order to adjust expectations. Since participants spent most of the time looking at the tablet, some reactions from the robot may have been missed and impacted the perception of robot behavior. Further, some participants described the robot as "intimidating", which could have impacted their comfort in thinking aloud.

Future work may help elucidate some questions raised by the results reported here. Age and cultural background may impact user behavior, and these aspects were not taken into account.
There is also potential to explore if different-looking social robots affect the sessions, or the impact of a human listener. It would be interesting to further investigate if the outcomes of the TAPS session change with a different set of tasks (e.g. programming or debugging tasks).

Regarding listener behavior implementation, we did not explore multimodal data input (such as facial expressions, gaze tracking or head pose). Future work may explore these aspects, in line with what was developed in \cite{2019huang,2021murray}. The robot could also display more complex listening behavior, such as asking general questions to the user (e.g.\citet{1966weizenbaum}'s ELIZA). Finally, and while we explored transferability of datasets when implementing robot listener behavior, we note that behavior expectations may differ in a human-human versus a human-robot interaction \cite{moshkina2014social,gambino2020building}. This is yet another aspect that calls for further exploration.

\section{Conclusion} \label{sec:conclusion}

This study explored human-robot interaction in a new setting and broadened knowledge on the role and usefulness (or lack thereof) of social robots as adjacent tools for human task completion. Further, the development and/or implementation of two fully-autonomous robot listener behaviors provided an interesting contribution on the transferability of heuristics and conversational datasets into new interactive contexts. Our findings indicate that the presence of a social robot displaying elementary listener behavior is not sufficient to elicit substantial differences in human behavior and perception of a "rubber duck" think-aloud problem-solving session, when compared to an inanimate object. Additional studies are required to inform the design of socially-assistive robots in problem-solving - as the future brings us closer to social robots, an optimized "robot duck" could play an important role in the problem-solving tasks which populate our days.


\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}