\subsection{Social Navigation}
In traditional robot motion planning frameworks, humans are usually treated as obstacles that the robot should avoid~\cite{fox1997}. However, in social navigation scenarios, humans are different from most other obstacles -- they have a purpose, follow social rules, and react to other agents in the environment, including the robot~\cite{kruse2013}. Researchers have explored various methods for social-aware motion planning.

Understanding and modeling reactive behaviors of humans is necessary for planning socially acceptable actions. A popular framework to model human navigation in response to other agents is the Social Force Model (SFM)~\cite{helbing1995, zanlungo2011, farina2017}, which frames pedestrians as affected by interactive forces that drive them towards their goals and away from each other. SFM has been applied to human-robot interaction and social navigation scenarios~\cite{ratsamee2013, Ferrer2014, mehta2016, Che2018}. Some work in the field of multi-agent planning also consider the reactive nature of moving agents. For example, optimal reciprocal collision avoidance (ORCA)~\cite{van2011} assumes that each agent takes some responsibility of resolving pair-wise conflicts, which result in global collision free paths. Guy et al. extended ORCA to model human navigation by introducing additional parameters such as reaction time and personal space~\cite{guy2010}. These reaction-based methods are popular for crowd simulation because of the high computational efficiency. However, they tend to be less accurate in predicting individual trajectories, and do not capture the stochasticity of human behavior.

Researchers have applied probabilistic methods to modeling human behavior and planning for social navigation. The Interactive Gaussian Process model was proposed in~\cite{trautman2010unfreezing}; it takes into account the variability of human behavior and mutual interactions. Trautman et al. used this model to navigate a robot through a crowded cafeteria~\cite{Trautman2013}. Inverse Reinforcement Learning (IRL)~\cite{ng2000, levine2012, abbeel2005, ziebart2008maximum} is another framework that allows a probabilistic representation of human behavior. IRL was applied in~\cite{kuderer2012, kim2016, kretzschmar2016, pfeiffer2016} to learn pedestrian models and plan for socially compliant navigation. In general, these approaches aim to generate ``human-like'' actions, with the underlying assumption that the robot should behave similarly to a human in order to be socially acceptable. We also use IRL as part of our model to predict humans' responses. Our work differs from previous approaches in that we model jointly the effects of both explicit communication and robot motions on human navigation behavior.

Recently, deep neural networks also gained attention in robot navigation applications. Chen et al. incorporated social norms and multi-agent reactions into a deep reinforcement learning framework to enable robot navigation in a pedestrian-rich environment~\cite{chen2017}. Tai et al. applied a generative adversarial imitation learning strategy to generate socially compliant navigation from raw sensor data~\cite{tail2018}. Deep learning frameworks have the advantage that they often impose less assumptions on human behavior. However, it is unclear whether they generate motions that are natural and transparent. In addition, the amount of data required for training these neural networks can be large and expensive to obtain. Our approach is more data efficient and explainable by leveraging modeling techniques.

\subsection{Expressive Motion in Robotics}
Besides social navigation, researchers have investigated the problem of planning interactive and communicative motions in the field of computational human-robot interaction~\cite{thomaz2016}. Dragan et al. formalized the idea of legibility in~\cite{dragan2013legibility} and proposed methods to generate legible motions for robot manipulators~\cite{Dragan2013}. Legible motions are motions that express the robot's intent to observers, and can be different from predictable motions~\cite{dragan2013legibility}. Sadigh et al. modeled a human-robot system jointly as a dynamical system, and proposed methods to plan motions that actively influence human behaviors~\cite{sadigh2016planning} and further actively gather information of internal state of humans~\cite{sadigh2016information,sadigh2018planning}. Similar ideas were explored in human-robot collaboration scenarios. Bestick et al. developed a method to enable a robot to purposefully influence the choices of a person in handover tasks to help the person avoid suboptimal grasps~\cite{bestick2016}. Nikolaidis et al. studied human-robot mutual adaptation, and proposed a decision framework to allow the robot to guide participants towards a better way of completing the task~\cite{nikolaidis2017}. In our work, we consider robot motions that can inform or affect people as a form of implicit communication, as people often have to infer about the robot's intent with context.

\subsection{Non-Verbal Explicit Communication Methods}
In this work, we define actions that communicate information deliberately and unambiguously as explicit communication. The most common form of explicit communication is verbal communication. Non-verbal communication can be made explicit, such as pointing gestures~\cite{huang2013, lohse2014}, head nod~\cite{hashimoto2007} and visual displays~\cite{baraka2018}.

Haptic feedback has been used as an explicit communication mechanism in human-robot interaction. Scheggi et al. designed a vibrotactile haptic bracelet to guide the user along trajectories that are feasible for human leader-robot follower formation tasks~\cite{scheggi2012, scheggi2014}. The bracelet can display three distinct signals that inform human to slow down and turn toward left or right. Sieber et al. used haptic feedback to assist a user teleoperating a team of robots manipulating an object~\cite{sieber2015} by communicating the status (stable vs. in transient) of the robot team. Here, we use haptic feedback to explicitly convey the robot's discrete intent to the user in collision avoidance scenarios during navigation.

Non-verbal communication, such as gaze, are considered as implicit if it indirectly conveys information~\cite{mutlu2009, admoni2014, moon2014}. In this work, we use the motion of the mobile robot, which does not impose any constraint on the appearance of the robot and can naturally be incorporated in social navigation.
\section{Introduction} \label{sec: intro}
\input{intro.tex}

\section{Background} \label{sec: background}
\input{background.tex}

\section{A Planning Framework \\for Socially-Aware Robot Navigation}\label{sec: planning framework}
\input{framework.tex}

\section{Hybrid Model of Human Navigation} \label{sec: human model}
\input{human.tex}

\section{Planning for Communication} \label{sec: planning}
\input{planning.tex}

\section{Simulation Results} \label{sec: simulation}
\input{simulation.tex}

\section{Experimental Results} \label{sec: experiment}
\input{experiments.tex}

\section{Conclusion and Future Work} \label{sec: discussion}
\input{discussion.tex}

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}

\subsection{Experimental Setup}
We again used a TurtleBot2 as the mobile robot platform. The robot was equipped with a Hokuyo URG-04LX-UG01 laser range finder, and an ASUS Zenbook UX303UB as the onboard computer. The planning algorithm is implemented on a desktop computer with Intel Core i7 processor and 16GB RAM. The two computers communicate with each other through a wireless network.

The onboard computer processes the laser range finder readings to localize the robot and estimate the position and velocity of nearby pedestrians. We localize the robot in a static map using laser-based Monte Carlo Localization~\cite{thrun2001}. The robot detects and tracks pedestrians using a leg detection algorithm provided by ROS~\cite{quigley2009}. The algorithm first attempts to segment leg objects from laser range finder readings with a pre-trained classifier, then pairs nearby legs and associates the detection with a tracker for each person. The localization and tracking results are sent to the desktop computer that runs our planning algorithm. The planning algorithm computes the plan and sends it back to the robot.

The explicit communication is displayed to the user via a wearable haptic interface that consists of a single vibration motor (Haptuator Mark II, shown in Fig.~\ref{fig: setup} (d)). The interface is capable of rendering distinct signals by modulating the vibration pattern. In this experiment, we display haptic cues with different vibration amplitudes and durations to indicate the robot's intent. Robot priority is represented by a single long vibration (duration 1.5s, max current 250 mA), and human priority is represented by 3 short pulses (0.2 sec vibration with 0.2 sec pause in between, max current 150 mA). To train the user interpreting the haptic cues, we first display cues in random order (3 trials for each vibration pattern) and explain the meaning. Then we test the user with randomized cues. We found that all users can recognize each cue correctly. The training process is done for each user before the actual experiment.

The experiment field is illustrated in Fig.~\ref{fig: setup}(a), which is a room of size 8 m$\times$10 m. Two tables are placed at the two ends of the field as the targets for the user. To distract the user from focusing solely on the robot, we place questionnaires on the tables, and ask the user to remember and answer arithmetic questions. The entire experiment is recorded with an overhead camera (GoPro Hero 4, recording at 60 Hz), and we post-process the video to extract the trajectories of the user and the robot. To facilitate tracking, we ask the user to wear a purple hat, and we attach an ArUco marker~\cite{Aruco2014} on the top of the robot.

\subsection{Experimental Design}

\begin{figure*}[t]
	\centering
    \includegraphics[width=\textwidth]{figures/stats}
    \caption{Comparison of different metrics for three experimental conditions. (a) Percentage of time that the user passed in front of the robot. (b) Average path length. A base length (6.2 m) is subtracted from all for visualization purpose. (c) The user's trust in the robot. Brackets indicate statistical significance (*$p < 0.05$, **$p < 0.01$, ***$p < 0.001$). }
    \label{fig: stats}
\end{figure*}

We conducted a $3 \times 1$ within-subject experiment, manipulating the communication mode (baseline, implicit only, explicit + implicit) and measuring users' path lengths, whether they pass the robot in the front, and their trust in the robot. Details of the experiment design and procedure are described below.

\prg{Manipulated Factors} We manipulate one factor: the \textbf{\textit{communication mode}}. The user's performance under each communication mode reflects the effectiveness of the communication method. In addition, the experimental trials consist of different \textit{task priorities} and \textit{robot starting positions}. We randomize the order of task priority so that the user does not know the robot's true intent beforehand. We also vary the starting position of the robot to create richer interaction scenarios and prevent the user from developing a fixed pattern (e.g. always pass in front of the robot regardless of communication). Overall, the experiment is divided into three sessions based on the communication mode:
\begin{itemize}
    \item \textbf{\textit{explicit + implicit}}: the robot communicates its intent to the user both explicitly (via haptic feedback) and implicitly (by changing speed and direction), using a model to predict the user's movement.
    \item \textbf{\textit{implicit only}}: the robot does not plan for explicit communication, but still changes speed and direction according to a model that predicts the user's movement.
    \item \textbf{\textit{baseline}}: the robot simply performs collision avoidance with the user without predicting the user's movement.
\end{itemize}

Each session consists of 20 trials. In each trial, the human user and the robot both move from a target position to a goal position, as shown in Fig.~\ref{fig: setup}(b). The robot is assigned one of the two task priorities in each trial: robot priority ($c_h : c_r = 1:5$) or human priority ($c_h : c_r = 2:1$). We also vary the starting position of the robot. The starting positions can be classified as \textbf{\textit{far}} (4 out of 20 trials), \textbf{\textit{close}} (4/20) and \textbf{\textit{normal}} (12/20), based on the distance to the user's starting position (Fig.~\ref{fig: setup}(b)). There are equal numbers of high-priority trials and low-priority trials for each type of starting position. The order of the trials for each communication mode is pseudo-randomized. We counter-balance the order of three communication modes among all users using a Latin square design. The first 10 trials are treated as training trials and not included for analysis. In the training trials, we will reveal the robot's true intent to the user after each trial. We never tell the robot's intent to the user in the experimental trials. The purpose of the training trials is to familiarize the user with the experimental procedure and the robot's behavior. After the experiment, we ask the user to rate her trust in the robot on a scale from 1 to 7, 1 being the least trust and 7 being the most trust. The exact question we use in the survey is ``Please rate your overall trust in the robot in the first/second/third condition''. There are 7 choices under each question and we explain to the user the meaning of the scale.

In addition to the orthogonal encounter scenario, we also tested the planner in other experimental setups. Fig.~\ref{fig: setup} (c) presents human-robot encounter at different angles ($+45^\circ$ and $-45^\circ$). We found that the planner worked similarly in these scenarios.

\prg{Dependent Measures} We measure the user's path length for each trial. Path length is commonly used to estimate the user's efficiency and effort in navigation tasks. We also measure whether the user passes the robot in the front, which is an indicator of whether the user understands the robot's intent correctly (the user is expected to pass the robot in the front in human priority trials and from the behind in robot priority trials, if they understand the robot's intent correctly. This measurement reflects the transparency of the robot's behavior. Finally, we measure the user's trust in the robot. Trust is measured using a post-experiment survey as described previously. In the session where explicit communication is allowed, we also record whether explicit communication happens and the time at which it happens.

\prg{Hypothesis} We hypothesize that:
\begin{enumerate}[label=\Roman*., widest=IV, align=left, leftmargin=*]
\item Using \textbf{\textit{explicit + implicit}} communication conveys the robot's intent more clearly than \textbf{\textit{implicit only}} and \textbf{\textit{baseline}}, such that users will elect to pass in front of or behind the robot as appropriate for a given priority.
\item The user's average path length is shorter when the robot plans for communication with the human model (\textbf{\textit{explicit + implicit}} and \textbf{\textit{implicit only}} modes).
\item The user's overall trust rating in the robot is higher when the robot plans for \textbf{\textit{explicit + implicit}} communication, compare to the other two conditions.
\end{enumerate}

\prg{Subject Allocation} A total of 12 people (7 males and 5 females, age from 23 to 33) participated in the experiment after giving informed consent, under a protocol that was approved by the Stanford University Institutional Review Board. We used a within-subjects design and counterbalanced the order of the three sessions.

\subsection{Analysis and Results}
\begin{table}[t]
\caption{Mean and standard deviation of measurements. For the first two measurements, the first row is Human Priority and the second row is Robot Priority.}
\label{table: exp measurements}
\centering
\begin{tabular}{|@{\hspace{3pt}}l@{\hspace{3pt}}|c|c|c|}
\hline
Measure & Imp+Exp & Imp & Base\\
\hline\hline
\multirow{2}{1.5cm}{\% Passing in Front} & $0.94\pm0.12$ & $0.72\pm0.36$ & $0.33\pm0.27$\\ \cline{2-4}
& $0.11\pm0.15$ & $0.50\pm0.37$ & $0.55\pm0.41$\\
\hline
\multirow{2}{1.5cm}{Path Length} & $6.42\pm0.20$ & $6.62\pm0.29$ & $6.79\pm0.22$\\ \cline{2-4}
& $6.64\pm0.07$ & $6.68\pm0.10$ & $6.91\pm0.18$\\
\hline
Trust  & $6.58\pm0.49$ & $4.58\pm0.64$ & $3.33\pm0.85$\\
\hline
\end{tabular}
\end{table}

Fig.~\ref{fig: stats} and Table~I summarizes major results. We describe the analysis and results in detail in the following paragraphs.

\prg{Understanding Robot Intent} To characterize how often the user correctly understood the robot's intent, we compute the percentage of trials that the user passed in front of the robot for each communication mode and task priority. A one-way repeated measures ANOVA revealed a significant effect on percentage passing in the front for the communication mode ($F(2, 33) = 14.61, p < .001$ for human priority trials, $F(2, 33) = 5.75, p = .007$ for robot priority trials). We performed the statistical test separately for robot and human priority trials. We performed a post-hoc analysis with Tukey HSD to determine pairwise differences. For human priority trials, results show that the baseline mode is significantly different from the implicit only mode ($p = .005$) and the implicit + explicit mode ($p < .001$). For robot priority trials, the implicit + explicit mode is significantly different from the implicit only mode ($p = 0.027$) and the baseline mode ($p = 0.01$). This supports Hypothesis I that the user can better understand the robot's intent with explicit communication.

\prg{Path Length} Path length is a commonly used metric to measure efficiency. We computed average path length (over all trials of each condition and robot intent) of each user. Similar to the metric of passing in the front, we performed a one-way repeated measure ANOVA analysis, and results show a significant effect ($F(2, 33) = 6.36, p = .005$ for human priority trials, $F(2, 33) = 14.18, p < .001$ for robot priority trials). Post-hoc comparisons reveal that, for human priority trials, the average path length of the implicit + explicit mode is significantly shorter than the baseline mode ($p = .003$). For robot priority trials, the average path length of the baseline mode is significantly longer than the implicit only mode ($p < .001$) and implicit + explicit mode ($p < .001$). The result suggests that the user can navigate more efficiently when the robot plans for communication (Hypothesis II).

\prg{Trust} Past research has shown that human would trust the robot more if it can express its intent via motion~\cite{dragan2015} or other forms of communication~\cite{baraka2018}, and we expect similar results in our experiment. Here users' trust in the robot was measured with a post-experiment survey, as explained in the previous subsection. A one-way repeated measure ANOVA analysis revealed a significant effect for communication mode ($F(2, 33) = 64.5, p < .001$). Post-hoc Tukey HSD analysis showed that all three pairs are significantly different from each other ($p < .001$ for all), which supports Hypothesis III.

\begin{figure}[t]
	\centering
    \includegraphics[width=0.4\textwidth]{figures/velocity_profile}
    \caption{Velocity profiles from sample trials. Top: human priority trials. Bottom: robot priority trials.}
    \label{fig: velocity profile}
\end{figure}

\subsection{Effect of Communication}
Fig.~\ref{fig: stats}(a) shows that when there was both implicit and explicit communication, users were able to understand the robot's intent most of the time and acted accordingly -- users passed in front of the robot 92\% of the time in human priority trials, and only 9\% of the time in robot priority trials. When there was only implicit communication, users were more confused, especially in robot priority trials (about 50\% of the time users passed in front of the robot despite that the robot's intent was not to yield). When the robot performed collision avoidance naively (baseline condition), users were the most confused, and acted more conservatively.

Fig.~\ref{fig: velocity profile} depicts the robot's velocity over time from sample trials of different communication modes and task priorities. In human priority trials (top row), when the robot planned for communication (implicit only and implicit + explicit modes), it started to slow down relatively early, compared to the baseline mode. By doing so, the robot implicitly expressed its intent. Indeed, significantly more users chose to pass in front of the robot with these communication modes (Fig.~\ref{fig: stats}(a)). The reason that the planner generated this behavior is that, with the human model, the planner was able to predict the effect of robot motions on the user's navigation behavior. More specifically, our human model would predict that, given the robot slows down early, the human is more likely to believe that the robot is yielding priority and will choose to pass first, which would generate higher expected rewards. This encouraged the robot to slow down earlier. On the other hand, the baseline planner did not perform any prediction and simply solved for a collision-free trajectory. As a result, the robot did not slow down until it got very close to the user, which caused the user to believe that it did not want to yield. The result suggests that legible motions emerge when reactions of the human are taken into consideration during planning.

Our results also suggest that users can navigate more efficiently when the robot planned for communication. Fig.~\ref{fig: stats}(b) shows that for human priority trials, users' path length is significantly shorter in the implicit + explicit mode than the baseline mode. For robot priority trials, both the implicit + explicit mode and the implicit only mode result in significantly shorter path lengths than the baseline mode. Knowing the robot's intent, users can better coordinate their movements with the robot to navigate more efficiently. In the implicit + explicit mode, users in general got feedback early in the trial. In the implicit only mode, users had to infer the robot's intent, and tended to be less efficient. This result does not show statistical significance; users' behaviors have high variance in the implicit only mode. In Fig.~\ref{fig: stats}(a), we observe the same high variance in the implicit only mode.

An important feature of the study is that we never explained to the user how the robot may move differently in three communication modes. One of the main purpose of the training trials is to let the user figure out the robot's movement pattern. We observe that users in general were able to adapt to the robot's communication mode in a few trials. Most users started conservatively, trying to keep relatively large distances from the robot. In the implicit + explicit mode, most users adapted to trust the explicit communication quickly. In the baseline mode, users made their decisions largely based on the initial distance from the robot, as they found it difficult to tell the robot's intent from its motion. Interestingly, in the implicit only mode, we observed that users' behaviors varied based on their strategy of interacting with the robot. Aggressive users tend to realize earlier that the robot would slow down when it wanted to yield priority; they were more likely to cut in front of the robot. This helped them to interpret the implicit communication. Timid users, on the other hand, tend to keep larger distances or let the robot pass first. Sometimes this behavior resulted in the robot not slowing down in human priority trials. Unlike legibility planning as in~\cite{dragan2013legibility}, our algorithm does not explicitly optimize legibility. In addition, the algorithm responds to human actions in real time. As a result, the robot would not wait if the user is too far, even in human priority trials, which could confuse the user. This explains in part the high variance observed in the implicit only mode.

Note that in the implicit + explicit mode, the robot may not always use explicit communication. Again, this is the result of optimizing the overall reward function. We assigned a negative reward for performing explicit communication to prevent over-communicating. In trials where the robot started too close to or too far  from the intersection point, it often chose to not use explicit communication because the chance of collision if extremely low.

Finally, we show that users gained more trust in the robot with the proposed planner in Fig.~\ref{fig: stats} (c). We found statistical significance between each pair of communication modes. When the robot communicated both implicitly and explicitly, users trusted the robot the most. This is also the condition where users understood the robot's intent the best (based on the previous analysis shown in Fig.~\ref{fig: stats}(a)). The correlation suggests that the easier users can understand the robot's intent, the more they trust the robot -- transparency is beneficial in social navigation. When users do not trust the robot, they tend to act conservatively and be less efficient (fewer trials where the user passed in front of the robot, and longer path lengths in the baseline mode). Note that the trust measurement in this work is exploratory, as we only asked users to rate their overall trust over the robot. In order to obtain more accurate results, a more specific and validated questionnaire should be used.

\subsection{Modeling Continuous Navigation Actions of a Human}

\prg{Inverse Reinforcement Learning}
We employ a data-driven approach to model humans' continuous navigation actions. Specifically, we use maximum entropy inverse reinforcement learning (MaxEnt IRL)~\cite{ziebart2008maximum, levine2012} to learn the distribution of human actions that matches a set of provided demonstrations in expectation. We assume the humans are approximately optimizing a reward function $R_{h}^\psi(\traj{h}, \traj{r}, E)$ for a given discrete class $\psi$. Under this model, the probability of the actions $\mathbf{u}_{h}$ is proportional to the exponential of the total reward of the trajectory ($Z$ is a normalization vector):
\begin{equation}
\begin{split}
p_\psi(\traj{h} | \traj{r}, E) &= \frac{1}{Z} \exp \left( R_{h}^\psi(\traj{h}, \traj{r}, E) \right) \\
&= \frac{1}{Z} \exp \left( \mathbf{w}_\psi \cdot \mathbf{f}(\traj{h}, \traj{r}, E) \right)
\end{split}
\label{eq: max ent}
\end{equation}

We parameterize the human reward function as a linear combination of features $\mathbf{f}$ that capture relevant properties of the navigation behavior. The features are weighted by $\mathbf{w}_\psi$, which can be learned by maximizing the overall likelihood of the provided demonstrations $\mathcal{D}$:
\begin{equation}
\argmax_{\mathbf{w}} \prod_{\mathbf{u}_{h} \in \mathcal{D}} p_\mathbf{w}(\traj{h} | \traj{r}, E)
\label{eq: irl}
\end{equation}

\noindent The weight vector $\mathbf{w}_\psi$ is a function of the behavior mode: $\psi$. Given a selected set of features such as distance to the other agents, heading, or velocity, we learn appropriate weights corresponding to the humans' reward functions for each discrete class $\psi$ from collected demonstrations. We will discuss the specific features used in this work in Sec.~\ref{sec: simulation}.

\subsection{Modeling Behavior Modes of a Human}
When modeling the distribution of behavior modes of a human, we consider the effect of both implicit communication (via robot movements) and explicit communication (via haptic feedback). We assume that the human will infer the robot's intent, and act cooperatively during the interaction. Mathematically, this suggests that the distribution of behavior modes of the human is related to her belief over the robot's intent, $\theta$. In our particular social navigation problem, we assume this belief is equal to probability of choosing a behavior mode:
\begin{equation}
p(\psi | \traj[-N:0]{r}, \traj[-N:0]{h}, \overline{u_c}, t_c) = p(\theta | \traj[-N:0]{r}, \traj[-N:0]{h}, \overline{u_c}, t_c)
\label{eq: belief model assumption}
\end{equation}
In other words, if the human believes that there is an 80\% chance that the robot will yield priority, she will behave in ``human priority'' mode with the same probability of 0.8. Non-linear relationships between $p(\psi)$ and $p(\theta)$ could also be used to model risk-averse (more willing to yield priority) and risk-seeking behaviors.

Using Bayes rule, we can transform Eq.~(\ref{eq: belief model assumption}) to:
\begin{equation}
\begin{split}
&p(\theta | \traj[-N:0]{r}, \traj[-N:0]{h}, \overline{u_c}, t_c)
\\
&\propto  p(\traj[-N:0]{r}, \overline{u_c}, t_c | \theta, \traj[-N:0]{h}) \cdot p(\theta | \traj[-N:0]{h}) \\
 & = p(\traj[-N:0]{r} | \theta, \traj[-N:0]{h})\cdot p(\overline{u_c}, t_c | \theta)\cdot p(\theta)
\end{split}
\label{eq: discrete model factorization}
\end{equation}

The second step in Eq.~(\ref{eq: discrete model factorization}) assumes conditional independence of the robot trajectory $\traj[-N:0]{r}$ and explicit communication $u_c$. With this factorization, we can separately model the effect of robot motion (implicit communication) and explicit communication on behavior modes of the human. The last term $p(\theta)$ is the prior on the robot's intent, which should be decided based on the application. In our implementation, there is an equal chance for the robot communicates different intents. So we choose a uniform prior.

The formulation casts the backward inference problem (from action to intent) into forward prediction problems (from intent to action). Next, we derive methods to compute each part of Eq.~\eqref{eq: discrete model factorization}.\\

\prg{Implicit Communication}
We assume that the human in general expects the robot to be rational and efficient according to some reward function. Applying the principle of maximum entropy, we model the human as expecting robot movements with probability:
\begin{equation}
    p(\traj[-N:0]{r} | \theta, \traj[-N:0]{h}) = \frac{\exp \left( R_r^\theta(\traj[-N:0]{r}, \traj[-N:0]{h}) \right)}{\int_{\mathbf{\xi}'_{r}} \exp \left( R_r^\theta(\mathbf{\xi}'_{r}, \traj[-N:0]{h}) \right) \mathrm{d} \mathbf{\xi}'_{r}}
    \label{eq: implicit model}
\end{equation}
where $R_r^\theta(\traj[-N:0]{r}, \traj[-N:0]{h})$ is a reward function that the human expects the robot to optimize, given its intent $\theta$. To compute the integration in the denominator of Eq.~\eqref{eq: implicit model}, we use the second-order Taylor series to approximate $R_r^\theta(\traj[-N:0]{r}, \traj[-N:0]{h})$.\\

\prg{Explicit Communication}
The belief of the human over the robot's intent is strongly affected by explicit communication, because the intent is directly conveyed. However, the effect of explicit communication should decay over time, as the robot's intent may change, and only short-term intents are communicated. Inspired by a model of human short-term verbal retention~\cite{peterson1959}, we propose:
\begin{equation}
    p(\overline{u_c}, t_c | \theta) = \frac{1}{Z} \left\{
    \begin{array}{ll}
         A \exp \left( -\frac{t - t_c}{M} \right) + 1 & \overline{u_c} = \theta \\
         1 & \overline{u_c} \neq \theta
    \end{array}
    \right.
    \label{eq: explicit model}
\end{equation}

\noindent Here $A$ and $M$ are parameters that determine the characteristic of the distribution, and $Z$ is a normalization factor.
The explicit communication initially reflects the true intent with very high probability. However, the inference strength decays over time, and eventually the communication becomes irrelevant to the robot's latest intent.

\begin{figure}[t]
	\centering
    \includegraphics[width=0.4\textwidth]{figures/explicit_model_demo}
    \caption{Inference on the robot's intent based on explicit communication. Assuming that $\overline{u_c} = \theta_1$ happened at $t_c = 0$}
    \label{fig: explicit model demo}
\end{figure}

To clarify this model, let us consider a simpler version of Eq.~\eqref{eq: belief model assumption}: $P(\theta | \overline{u_c}, t_c)$ -- to infer the robot's intent with only explicit communication. Assume that $\theta \in \{\theta_1, \theta_2\}$ is binary, $\overline{u_c}= \theta_1$, and $t_c = 0$, and use Bayes rule again:
\begin{equation}
\begin{split}
    p(\theta = \theta_1 | \overline{u_c}, t_c) &= \frac{p(\overline{u_c}, t_c | \theta_1) p(\theta_1)}{p(\overline{u_c}, t_c |\theta_1) p(\theta_1) + p(\overline{u_c}, t_c | \theta_2) p(\theta_2)}\\
    &= \frac{A \exp \left(- (t - t_c) / {M} \right) + 1}{A \exp \left(- (t - t_c) / {M} \right) + 2}
\end{split}
\end{equation}
The change of $p(\theta = \theta_1 | \overline{u_c}, t_c)$ over time is plotted in Fig.~\ref{fig: explicit model demo}. Initially, the robot's intent can be inferred with high probability given the explicit communication. As time passes, the inference strength decreases and eventually the belief over the robot's intent becomes uniform (0.5).

In the scenarios we consider, the human infers the robot's intent with both implicit and explicit communication. The combined effect can be computed with Eqs.~(\ref{eq: discrete model factorization})~--~(\ref{eq: explicit model}) given in this section.

\subsection{Robot Reward Function}
The overall reward $R_r$ that robot optimizes in Eq.~\eqref{eq: computing optimal u_r} consists of four parts that quantify \textit{robot efficiency}, \textit{human comfort}, \textit{safety}, and \textit{reward of explicit communication}:

\prg{Robot Efficiency} The robot should get as close as possible to its target with minimum effort:
    \begin{equation}
        R_{r,r} (\traj{r}) = \norm{\st[N]{x}{r} - x_{r}^{g}} +  \gamma_r \sum_{t=1}^{N} \norm{\st{u}{r}}^2,
    \end{equation}
    where $\gamma_r$ is a weight that balances effort and goal. This component determines how fast the robot approaches the goal. $\st[N]{x}{r}$ is the robot position at time step $N$, and $x_{r}^{g}$ is the goal position. $\norm{\cdot}$ denotes Euclidean norm

\prg{Human Comfort} We want the human to spend less effort to achieve her goal:
    \begin{equation}
        R_{r,h} (\traj{h}) = \norm{\st[N]{x}{h} - x_{h}^{g}} + \gamma_h \sum_{t=1}^{N} \norm{\st{u}{r}}^2.
    \end{equation}
    Similar to $\gamma_r$, $\gamma_h$ balances between human reaching her goal and her effort.

\prg{Safety} The robot should avoid collisions with the human:
    \begin{equation}
        R_{r,\text{safety}} (\traj{r}, \traj{h}) = \sum_{t=1}^{N} \exp \left( - \frac{\norm{\mathbf{x}_{h,t} - \mathbf{x}_{r,t}}^2}{D_r^2} \right).
    \end{equation}

\prg{Explicit Communication} Too much explicit communication will distract or annoy the human. Therefore, we set a constant reward for performing each explicit communication:
    \begin{equation}
        R_{r,EC} (u_c) = \gamma_{EC} \cdot \mathds{1}(u_c \neq \text{``none"})
    \end{equation}
where $\mathds{1}$ is the indicator function.

\noindent The components defined above are then combined linearly with weights $c_i$:
\begin{equation}
	R_r(\traj{r}, \traj{h}, \ensuremath{u_c}) = \sum_{i=0}^4 c_i \cdot R_{r,i}(\traj{r}, \traj{h}, \ensuremath{u_c})
	\label{eq: robot reward function}
\end{equation}

\subsection{Human-Aware Planning}
To solve for the optimal robot actions using Eq.~\eqref{eq: computing optimal u_r}, we need to compute the expectation of the robot reward over the distribution of $\mathbf{u}_h$. In the last section, we derived a hybrid model for human behavior. To speed up computation, here we only consider the most likely trajectory given each possible human intent. With this, the expectation in Eq.~\eqref{eq: computing optimal u_r} can be expressed as:
\begin{equation}
    \mathbb{E}_{\traj{h}} \left[ R_r \right] = \sum_{\psi} p(\psi) \cdot R_r(\traj{r}, \traj[\psi]{h},\ensuremath{u_c}),
    \label{eq: expected reward}
\end{equation}
where $\traj[\psi]{h}$ is the most likely human response given $\psi$, and can be computed as:
\begin{equation}
\traj[\psi]{h} = \argmax_{\traj{h}} \mathbf{w}_\psi \cdot \mathbf{f}(\traj{h}, \traj{r})
\label{eq: most likely actions}
\end{equation}
The distribution $p(\psi)$ is given by Eqs.~\eqref{eq: belief model assumption} and~\eqref{eq: discrete model factorization}.

We use gradient-based optimization to solve the optimal robot controls $\mathbf{u}_{r}^*$. This requires the gradient information of the expected reward in Eq.~\eqref{eq: expected reward}. Letting $\overline{R}_r = \mathbb{E}_{\traj{h}} \left[ R_r \right]$, we aim to find:
\begin{equation}
    \frac{\partial \overline{R}_r}{\partial \mathbf{u}_{r}} = \frac{\partial \overline{R}_r}{\partial \mathbf{u}_{r}} + \sum_{\psi} \frac{\partial \overline{R}_r}{\partial \mathbf{u}_{h}} \frac{\partial \mathbf{u}_{h}^{\psi}}{\partial \mathbf{u}_{r}},
\end{equation}
where $\mathbf{u}_{h}^{\psi}$ is the most likely human actions. As we have a symbolic representation of $R_r$, both $\frac{\partial \overline{R}_r}{\partial \mathbf{u}_{r}}$ and $\frac{\partial \overline{R}_r}{\partial \mathbf{u}_{h}}$ can be computed analytically. Following the implicit differentiation method discussed in~\cite{sadigh2016planning}, we compute the last unknown term using the following expression:
\begin{equation}
    \frac{\partial \mathbf{u}_{h}^{\psi}}{\partial \mathbf{u}_{r}} = \left[ \frac{\partial^2 R_h^{\psi}}{\partial \mathbf{u}_{h}^2} \right]^{-1} \left[ - \frac{\partial^2 R_h^{\psi}}{\partial \mathbf{u}_{h} \partial \mathbf{u}_{r}} \right]
\end{equation}
We have assumed a discrete set of communication actions $\mathcal{U}_c \cup \text{``none"}$.
To obtain the optimal explicit communication $u_c^*$, we enumerate all possible communication actions $u_c \in \mathcal{U}_c \cup \text{``none"}$ and solve the optimization for each $u_c$. This can be done in parallel because the optimizations do not depend on each other.

\begin{algorithm}[tb]
	\KwData{$\mathcal{U}_c$ - set of explicit communicative actions}
	\SetKwData{TimerI}{timer1}	\SetKwData{TimerII}{timer2}
	\SetKwFunction{SendHaptics}{SendHapticMessage}

    \SetKwFunction{GoalReached}{GoalReached}
    \SetKwFunction{GetTrackingInfo}{GetTrackingInfo}
    \SetKwFunction{UpdateModel}{UpdateModel}
    \SetKwFunction{PredictOneStep}{PredictOneStep}
    \SetKwFunction{GenerateInitGuess}{GenerateInitGuess}
    \SetKwFunction{ComputePlan}{ComputePlan}
	\SetKwFunction{Communicate}{Communicate}
    \SetKwFunction{ExecuteControl}{ExecuteControl}

    \SetKwProg{ForParallel}{for}{ do in parallel}{end}

    \Repeat {\GoalReached{$\st{x}{r}$, $x_{g, r}$, $\epsilon$}} {
    	\tcc{get current+predicted states}
    	$\st{s}{r}, \st{u}{r}, \st{s}{h} \gets$ \GetTrackingInfo{}\;
        $\st{\hat{s}}{r}, \st{\hat{s}}{h} \gets$ \PredictOneStep{$\st{s}{r}, \st{s}{h}$}\;

        \BlankLine
        \tcc{update the hybrid model}
        $p(\psi) \gets$ \UpdateModel{$p(\psi), \st{s}{r}, \st{u}{r}, \st{s}{h}, u'_c, t_c$}\;

        \BlankLine
        \tcc{find the optimal actions}
        $\mathbf{u}_{r}^{0}, \mathbf{u}_{h}^{0} \gets$ \GenerateInitGuess{$\st{\hat{s}}{r}, \st{\hat{s}}{h}, x_{g,r}, x_{g,h}$}\;

        \ForParallel {$u_c \in \mathcal{U}_c \cup \{null\}$ } {
        	$\mathbf{u}_r(u_c), \overline{R}_r(u_c) \gets$ \ComputePlan{$\mathbf{u}_{r}^{0}, \mathbf{u}_{h}^{0}, \st{\hat{s}}{r}, \st{\hat{s}}{h}, p(\psi), u_c$}\;
        }

        $u_c^* \gets \argmax_{u_c} \overline{R}_r(u_c)$ \;
        $\mathbf{u}_r^* = \mathbf{u}_r(u_c^*)$\;
        }

        \BlankLine
        \tcc{execute the actions}
        \If {$u_c^* \neq null$} {
        	\Communicate{$u_c^*$}\;
            $u'_c \gets u_c^*$\;
            $t_c \gets t$\;
        }

        \ExecuteControl{$\mathbf{u}_r^*$}\;

        \BlankLine
        $t \gets t + 1$\;

	\caption{Outline of the planning algorithm.}
	\label{alg: planning framework}
\end{algorithm}

The planning algorithm is outlined in Alg.~\ref{alg: planning framework}. At every time step, the algorithm first retrieves states of the robot and the human, and performs a one step prediction of robot and human states in order to compensate time spent for planning. Then it updates the belief over the behavior mode of the human $\psi$ with new observations using equations (\ref{eq: belief model assumption}) and (\ref{eq: discrete model factorization}). Before optimizing for robot actions, the algorithm needs to generate an initial guess for the initial state of the human. If there is a plan from the previous time step, the plan is used as the initial guess. When there is no previous plan (first time step, or first detection of human), we generate the initial guesses as follows: First, we compute the robot actions ignoring the human. We use a feedback control based policy to steer the robot towards its goal~\cite{astolfi1999}. We then compute the human actions using an attracting potential field at the goal position, and a repelling potential field centered at the robot. The algorithm then performs a set of optimizations, one for each explicit communication action. Taking advantage of modern multi-core CPUs, we can run these optimizations in parallel. Finally, the best explicit communication $u_c^*$ and robot movements (implicit communication) $\mathbf{u}_r^*$ are selected and executed.

We implemented the planning algorithm in C++, and used the software package NLopt~\cite{johnson2014} to perform numerical optimization. In our implementation, we chose a planning horizon of $N = 6$, and a time step of 0.5 seconds.

\subsection{Social Navigation Scenario}
We consider the scenario shown in Fig.~\ref{fig: intro} for our simulation and experimental evaluation: a human and a robot move in (approximately) orthogonal directions and encounter each other. In this scenario, they have to coordinate with each other and decide who should pass first to avoid collision. We chose this scenario over the face-to-face encounter scenario because in the face-to-face encounter, there is a relatively strong social norm that people pass each other on the right hand side. We prefer to start with a scenario with no such bias to allow us better focus on the effect of communication. The robot can explicitly communicate its intent: to yield priority to the human (human priority) or not (robot priority). To numerically evaluate the model and compute the optimal plan, we need to define features leading to a learned reward function, which can then be used to model humans' continuous actions Eq.~\eqref{eq: max ent} as well as behavior modes Eq.~\eqref{eq: implicit model} of the human.

We select features that are commonly used to characterize social navigation of pedestrians~\cite{hoogendoorn2003}:
\begin{description}[leftmargin=0cm]
\item [Velocity.] Pedestrians tend to maintain a constant desired velocity. We use a feature that sums up the squared velocity along the trajectory:
\begin{equation}
f_{velocity} = \sum_{t=1}^{N} \norm{\st{\dot{x}}{h}}^2
\end{equation}
Here, and in the following equations, $||\cdot||$ represents L2 norm.

\item [Acceleration.] To navigate efficiently, pedestrians usually avoid unnecessary accelerations:
\begin{equation}
f_{acceleration} = \sum_{t=1}^{N} \norm{\st{\ddot{x}}{h}}^2 = \sum_{t=1}^{N} \norm{\st{u}{h}}^2
\end{equation}

\item [Distance to goal.] Pedestrians typically try to get to the target position as close as possible given a time horizon:
\begin{equation}
f_{goal} = \norm{\st[N]{x}{h} - {x}_{\text{goal}}}
\end{equation}

\item [Avoiding static obstacles.] Pedestrians avoid static obstacles in the environment:
\begin{equation}
f_{obstacle} = \sum_{t=1}^{N} \exp \left( - \frac{\norm{\st{x}{h} - \st{o}{\text{closest}}}^2}{D_o^2} \right)
\end{equation}
where $o_{\text{closest},t}$ is the position of the closest obstacle to the human at time $t$, and $D_o$ is a scaling factor.

\item [Collision avoidance with the robot.] Pedestrians tend to avoid each other in navigation. We assume that they avoid the robot in a similar fashion:
\begin{equation}
f_{avoidance} = \sum_{t=1}^{N} \exp \left( - \frac{\norm{\st{x}{h} - \st{x}{r}}^2}{D_r^2} \right)
\end{equation}

\item [Avoiding the front side of the robot.] In addition to avoiding the robot, we observe that humans tend to not cut in front of the robot, especially when they think that the robot has priority. This behavior is captured by the feature:
\begin{equation}
f_{front} = \sum_{t=1}^{N} \exp \left( - (\st{x}{h} - \st{x}{f})^\intercal \mathbf{D}_f^{-1} (\st{x}{h} - \st{x}{f}) \right)
\end{equation}
where $\st{x}{f}$ is a position in front of the robot, and $\mathbf{D}_f$ scales the Gaussian and aligns it with the robot's orientation.
\end{description}

\begin{figure}[t]
	\centering
    \includegraphics[width=0.5\textwidth]{figures/features_vis}
    \caption{Visualization of (a subset of) the features. Warm color indicates high reward and cool color indicates low reward. Human trajectory is in black, and robot trajectory is in red. Arrows indicate the positions and moving directions at the specified time step.
    }
    \label{fig: features}
\end{figure}

The features (except for velocity and acceleration) are visualized in Fig.~\ref{fig: features}. Human and robot trajectories are from a demonstration we collected to train the IRL model. The figure shows that the human indeed avoided low reward regions (cool color) and navigated to the high reward region (warm color).

\begin{figure*}[t]
	\centering
    \includegraphics[width=0.95\textwidth]{figures/cross_val}
    \caption{Sample predictions using the learned reward functions and cross validation. (a)-(b) show an example where the prediction matches the actual measurement. (c)-(d) show an example where initially the prediction doesn't match the measurement. We use the model to re-predict human actions at each time step, and the prediction starts to match measurement after two seconds. (e) Cross validation of prediction accuracy. Type I corresponds to cases where the model initially predicts the same homotopy class as the demonstration ((a), (b)), and type II corresponds to cases where the model initially predicts a different homotopy class ((c), (d)). Blue indicates errors of the IRL model, and orange indicates errors of the social force (SF) model. }
    \label{fig: cross validation}
\end{figure*}

\subsubsection*{Rewards for the Behavior Mode Model}
In general, it is challenging to to define the expected reward function of the human $R_h^\psi(\traj[-N:0]{r}, \traj[-N:0]{h})$, as it is not measurable. Here we design two reward functions based on our observations. Data driven methods like IRL can potentially produce more accurate results. Our goal, however, is to derive an approximate model that can inform the planning algorithm. Therefore, we prefer to start with simpler designs:

    \prg{Human Priority ($\psi = \psi_h$)} When the robot's intent is to yield priority to the human, it should avoid moving directly towards the human and getting too close. We use a reward function:
    \begin{equation}
        R_h^{\psi_h} = -\sum_{t = -N}^{0} \frac{\hat{n}_t \cdot (\st{\dot{x}}{r} - \st{\dot{x}}{h})}{\max(\norm{\st{x}{r} - \st{x}{h}}^2, 1.0)}
        \label{eq: implicit reward human priority}
    \end{equation}
    where $\hat{n}_t$ is a unit vector that points from the robot to the human: $\hat{n}_t = (\st{x}{r} - \st{x}{h}) / {\norm{\st{x}{r} - \st{x}{h}}}$. Note that here $t \leq 0$ means that the summation is over the previous time steps.

    \prg{Robot Priority ($\psi = \psi_r$)}  When the robot has priority, the human would expect it to move with a desired velocity, regardless of the state of the human:
    \begin{equation}
        R_h^{\psi_r} = -\sum_{t = -N}^{0} \frac{|\st{v}{r} - v_d|}{\max(\norm{\st{x}{r} - \st{x}{h}}^2, 1.0)}
        \label{eq: implicit reward robot priority}
    \end{equation}
    Here $\st{v}{r}$ is the robot's speed and $v_d$ is the desired speed.

In Eqs.~(\ref{eq: implicit reward human priority}) and (\ref{eq: implicit reward robot priority}), the reward at each time step is divided by $\max(\norm{\st{x}{r} - \st{x}{h}}^2, 1.0)$ because when the robot is relatively far from the user, its actions carries less information about its intent. We use the $\max()$ operation to achieve numerically stable  results.

\subsection{Human Model Evaluation}
\prg{Data Collection for Learning Human Model}
We collected navigation demonstrations from four users. Each user was asked to walk back and forth between two target locations (setup is similar to that of Section~\ref{sec: experiment}). At the same time, a mobile robot (TurtleBot2, Open Source Robotics Foundation, Inc.) also moved between a set of target positions in the environment. The target positions were selected so that the robot would cross path with the user. The user was allowed to walk freely as in normal human environment, while avoiding collision with the robot.

We collected data for two scenarios: the human priority scenario and the robot priority scenario. When the user had priority, the robot would slow down and yield to the user upon encounter (we decreased the maximum allowable speed of the robot as it got close to the user). When the robot had priority, it would simply ignore the user. In each condition, the user was told the intent of the robot, so the user always knew the robot's true intent. To distract the user from focusing solely on the robot, we asked the user to remember two numbers at one target location, and answer arithmetic problems using the two numbers at the other target location.

For each user, we recorded 64 trials for each scenario (a trial is moving from one target location to the other). We used 80\% of the data to train the model, and the other 20\% as the testing set.

\prg{Prediction with Learned Reward Function}
Using MaxEnt IRL, we recover the human reward functions in human priority and robot priority scenarios. The reward functions can be used to compute the most likely continuous navigation actions using Eq.~\eqref{eq: most likely actions}. This computes the actions over a fixed time horizon. To predict the entire trajectory, we use an MPC approach: we compute the most likely actions starting from time step $k$, predict the human state at $k+1$, and recompute the most likely actions starting from $k+1$.

Fig.~\ref{fig: cross validation} presents two prediction examples. In the first example ((a)-(b)), the predicted trajectory matches the measured trajectory very well (same homotopy class, type I). In the second example ((c)-(d)), the model initially predicts a trajectory that passes the obstacle on a different side (different homotopy class, type II). Although the prediction does not match the measured trajectory, it is still reasonable. In this example, both the predicted and measured trajectories pass behind the robot. Passing on either side of the obstacle has little effect on the reward and is almost equally likely. When predicting at $t=2.0$ s, because the human already started to move toward one side of the obstacle, the model is able to
generate a more accurate prediction.

The cross validation result is shown in Fig.~\ref{fig: cross validation}(e). The prediction error is calculated as the average Euclidean distance between the predicted and measured trajectories:
\begin{equation}
    e = \frac{1}{N'} \sum_{t=1}^{N'} \norm{\st{x}{\text{pred}} - \st{x}{\text{meas}}}
\end{equation}
where $N'$ is the total number of time steps. We compute the prediction error separately for cases where the initial prediction is in the same homotopy class as the measurement (type I), and cases where the initial prediction and the measurement are in different homotopy classes (type II). It can be observed that the initial prediction error is much smaller for type I, but both types becomes more accurate as the human gets closer to the goal. Compared with a social force model we used previously~\cite{Che2018}, the model described here performed better for this scenario.

\prg{Evaluation of the Behavior Mode Model}
\begin{figure*}[t]
	\centering
    \includegraphics[width=1.0\textwidth]{figures/belief_model_examples}
    \caption{Demonstrations of the behavior mode model in different scenarios. (a) A scenario where the robot slowed down to let the user pass first. Top plots show the trajectories of the robot and the user, and a map of the environment at three different time steps. Bottom plot shows the user's belief over the robot's intent (human priority) over time, given different explicit communication at $t=0.5$ s, predicted by our model. (b) A scenario where the robot didn't slow down and passed first.}
    \label{fig: belief model examples}
\end{figure*}
In addition to continuous trajectories, our model can also predict behavior modes, or equivalently, the belief of the human over the robot's intent. As it is impossible to measure this belief in an experiment, we aim to show that the prediction is reasonable with a few test scenarios.

Fig.~\ref{fig: belief model examples}(a) presents a scenario where the robot slowed down to let the user pass first. The scenario is illustrated with the three plots in the top row, showing trajectories and the environment at different time steps. The data is from one demonstration we collected for training the IRL model. The plot in the bottom row shows the predicted belief $p(\theta=\theta_h)$ over time, given no explicit communication, and two different explicit communication actions at $t_c = 0.5$ s: communicating human priority ($u_c = \theta_h$) and communicating robot priority ($u_c = \theta_r$). In the case of no communication, the belief decreases initially, but rises up as the robot stops to yield to the user. This result matches the intuition that the user can infer the robot's intent (human priority) to some extent by observing its movement (slowing down). When the robot communicates human priority at the beginning, the belief jumps to a high value and maintains this value afterwards. An interesting test scenario is when the robot communicates robot priority, which is not its true intent. According to our model, initially the belief drops, indicating that the user believes the explicit communication. However, when the robot slows down to let the user pass first, the belief starts to increase, indicating that the user starts to believe otherwise as she observes the robot movements.

We also tested scenarios where the robot's intent is to give itself priority. An example is shown in Fig.~\ref{fig: belief model examples}(b). Similarly, we show the change in belief given different explicit communication in the bottom plot. These examples demonstrate that our model can capture the effect of both explicit and implicit communication, and the prediction matches our intuition.

\subsection{Case Study in Simulation}
\begin{figure*}[t]
	\centering
    \includegraphics[width=1.0\textwidth]{figures/simulation_result}
    \caption{(a) Simulated scenario using the proposed interactive planner. In this scenario, the robot slows down and explicitly communicates its intent (human priority) to the human. Top: trajectory snippets at five individual time steps. Black and red arrows represent the positions and moving directions of the human and the robot at the given time step. Thin, solid red lines represent the robot's planned trajectories, thin solid/dashed black lines represent predicted human trajectories. Bottom: Velocity profile of the robot and the human. (b) Various starting positions of the robot vs. whether explicit communication is used by the planner for the same scenario. Red region represent the starting positions where explicit communication is used.}
    \label{fig: simulation example}
\end{figure*}

Before conducting a user study with a physical mobile robot (described Section~\ref{sec: experiment}), we first tested our algorithm in simulation. The purpose of the simulation is two-fold: first to validate that the proposed algorithm can work in idea setups,
and second to select appropriate parameters that produce reasonable behaviors.

Fig.~\ref{fig: simulation example}(a) illustrates a simulated scenario where the human and the robot move along orthogonal paths. Here the simulated human user follows a predefined trajectory, and the trajectory is deterministic regardless of the robot's actions. While this is not realistic, the purpose is to test whether the planner can generate reasonable robot behaviors. We will further validate the effectiveness of the plan with real-world user studies. The subplots in the top row shows the scenario at 5 different time steps, and the subplot in the bottom row shows the speed of the user and the robot over time. We can observe that the robot starts to slow down at $t = 1$ s, and explicitly communicates its intent to the user at $t = 1.5$ s. The robot continues moving slowly to allow the user pass first, and then speeds up towards its goal.

\begin{figure*}[thpb]
	\centering
    \includegraphics[width=\textwidth]{figures/exp_setup}
    \caption{(a) Experimental field. (b-c) Different starting and goal positions for the robot. (d) The wearable haptic interface and the vibration motor.}
    \label{fig: setup}
\end{figure*}

To understand why the planner does this, we visualize the generated plan and the predictions of human movements, as shown in Fig.~\ref{fig: simulation example}(a). As the robot speeds up and approaches the intersection, the user becomes unsure of how to avoid the robot. This is reflected in the second subplot in the first row, where the user's possible future trajectories diverge. In this simulation, we set the reward function coefficients $c_h \gg c_r$, so that the planner cares more about the human user's efficiency and comfort. As a result, the planner explicitly communicates human priority to the user to minimize the chance that the user would slow down to yield to the robot.

We tested the planner with various robot starting positions and target positions. Fig.~\ref{fig: simulation example}(b) presents the relationship of robot starting position and whether explicit communication is generated by the planner, for one specific target position. We can see that explicit communication is only used if the robot starts within the banded region in red. When the robot starts too close to the intersection (upper-right grey region), the planner predicts that the user will let the robot pass first as it is closer to the intersection. When the robot starts too far (bottom-left grey region), the planner predicts that the user will pass first. Explicit communication is used only when it becomes ambiguous who should pass first.

We also studied the effect of the coefficients $c_h$ and $c_r$ in the reward function in Eq.~(\ref{eq: robot reward function}). Setting $c_h \gg c_r$ resulted in a submissive robot that would yield to the user when there was a potential collision. Conversely, setting $c_h \ll c_r$ resulted in aggressive robot behaviors. When $c_h \approx c_r$, the robot tends to be submissive. In our implementation, the robot is usually slower than the human, so it is often more efficient to let the user pass first.
Based on the simulation results, we decided to experimentally test two sets of parameters: $c_h : c_r = 1:5$ as \emph{robot-prioritized} scenario, and $c_h : c_r = 2:1$ as \emph{human-prioritized} scenario.


\section{Introduction}

It is well known that frustration in magnetic materials
enables exotic quantum states to be realized.
However, not many such quantum states are obtained
in a mathematically rigorous form.
The Shastry-Sutherland model\cite{Shastry_Sutherland_Physica1981}
is a member of the family of mathematical models
in which not all but only some of the eigenstates are exactly obtained
\cite{Majumdar_Ghosh_1969,Majumdar_1970,Caspers_1982,
AKLT_PRL1987,Lange_ZPhys1994,Long_Siak_1993,HN_VBS_PRB1996,
HN_JPSJ1997,HN_JPSJ1998,Tonooka_JPSJ2007}.
Among such systems, the Shastry-Sutherland model
became important after a good candidate material, SrCu$_{2}$(BO$_{3}$)$_{2}$,
was discovered\cite{Kageyama_PRL1999,SMiyahara_KUeda_PRL1999}.
The discovery was followed
by extensive theoretical and experimental studies.

Between the region with the exact dimer ground state
and
the weakly frustrated region
with the typical N$\acute{\rm e}$el-ordered ground state,
the existence of the plaquette-singlet phase
was pointed out in Ref.~\ref{AKoga_NKawakami_PRL2000}.
Various approaches\cite{YFukumoto_JPSJ2000,
Lauchli_PRB2002,Lou_arXiv1212_1999,
Corboz_Mila_PRB2013,Wang_Batista_PRL2018}
have theoretically attempted to clarify
the behavior of the system in the intermediate region.
In the numerical-diagonalization studies, unfortunately,
the maximum of the treated system sizes
\---
32 spin sites
in Ref.~\ref{Lauchli_PRB2002} to the best of our knowledge
\---
was not so large.
On the other hand,
the pressure dependence of the spin gap was observed
experimentally\cite{HOhta_JPhysChem2015,Zayed_NatPhys2017,
TSakurai_JPSJ87_ESR_HP}.
Among these studies, Ref.~\ref{TSakurai_JPSJ87_ESR_HP},
employing an electron spin resonance study under
high pressure and high field, recently reported
the behavior of the spin gap
around the phase transition at the edge of the dimer phase.
This experimental result becomes a significant motivation
to theoretically clarify the behavior of the spin gap
around the intermediate region
from the numerical-diagonalization
calculations for even larger systems.

Under the circumstances, the purpose of the present paper is
to report numerical-diagonalization results of even larger systems
and
to examine the behavior of the system
around the intermediate region
between the dimer and N$\acute{\rm e}$el-ordered regions.
The present study provides results
for 36-site and 40-site systems in addition to smaller systems
to deepen our understanding of this system.
In particular, numerical diagonalizations
of the 40-site system require large-scale parallel calculations
in an appropriate supercomputer.
We successfully detect the edge of the dimer phase
and the edge of the N$\acute{\rm e}$el-ordered phase
from the results of the two large sizes.
Between the two edges, we additionally detect
the third ratio of the boundary dividing the intermediate region
into two parts, each state of which has characteristics
of a correlation function that are different from each other.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg1.eps}
\end{center}
\caption{(Color)
Lattice structure of the orthogonal dimer system,
namely, the Shastry-Sutherland model.
Thick solid lines and thin solid lines denote
bonds for $J_1$ and $J_2$, respectively.
Panel (a) shows the finite-size clusters
of $N_{\rm s}=16$, 32, and 36 by blue, green, and red dotted lines,
respectively.
Panel (b) shows the finite-size clusters
of $N_{\rm s}=20$ and 40 by blue and red dotted lines,
respectively.
}
\label{fig1}
\end{figure}

The Hamiltonian studied here is given by
\begin{equation}
{\cal H}
=
\sum_{\langle i ,j\rangle : \ {\rm orthogonal~dimer}} J_{1}
\mbox{\boldmath $S$}_{i}\cdot\mbox{\boldmath $S$}_{j}
+
\sum_{\langle i ,j\rangle : \ {\rm square~lattice}} J_{2}
\mbox{\boldmath $S$}_{i}\cdot\mbox{\boldmath $S$}_{j}
.
\label{Hamiltonian}
\end{equation}
Here, $\mbox{\boldmath $S$}_{i}$
represents the $S=1/2$ spin operator at site $i$.
We consider the case of an isotropic interaction
in spin space in this study.
Site $i$ is assumed to characterize the vertex
of the square lattice.
The number of spin sites is denoted by $N_{\rm s}$.
The first term of Eq.~(\ref{Hamiltonian}) denotes
orthogonal dimer interactions represented
by thick solid bonds in Fig.~\ref{fig1}.
The second term of Eq.~(\ref{Hamiltonian}) represents
interactions forming the square lattice
represented by thin solid bonds in Fig.~\ref{fig1}.
We consider that the two interactions between the two spins
are antiferromagnetic,
namely, $J_{1} > 0$ and $J_{2} > 0$.
Energies are measured in units of $J_{1}$;
hereafter, we set $J_{1}=1$.
We denote the ratio $J_{2}/J_{1}$ by $r$.
Note here that when $r=0$, the system is an assembly of
isolated dimerized-spin models,
whereas the system is reduced
to the $S=1/2$ Heisenberg antiferromagnet
on the ordinary square lattice
in the limit $r\rightarrow\infty$.

We treat finite-size clusters with system size $N_{\rm s}$
under the periodic boundary condition.
In this study,
$N_{\rm s}=16$, 20, 32, 36, and 40 are treated;
finite-size clusters are shown in Fig.~\ref{fig1}.
Note here that $N_{\rm s}/4$ is an integer
and that all the clusters are regular squares,
although the squares for $N_{\rm s}=20$, 32, and 40 are tilted.
The regular-square clusters help us
capture well the two dimensionality of the present system.

We carry out our numerical diagonalizations
on the basis of the Lanczos algorithm
to obtain the lowest energy of ${\cal H}$
in the subspace belonging to $\sum _j S_j^z=M$.
Note here that the $z$-axis
is taken as the quantized axis of each spin.
It is widely believed that
numerical-diagonalization calculations are unbiased.
Thus, one can obtain reliable information about the system.
The energy is denoted by $E(N_{\rm s},M)$,
where $M$ is an integer;
in particular, we calculate the cases $M=0$ and $M=1$
because our attention is focused primarily
on the behavior of the spin gap given by
\begin{equation}
\Delta=E(N_{\rm s},1) - E(N_{\rm s},0) .
\label{gap_determining_finite_N}
\end{equation}
Some of the Lanczos diagonalizations were carried out
using MPI-parallelized code that was originally
developed in the study of Haldane gaps\cite{HNakano_HaldaneGap_JPSJ2009}.
The usefulness of our program was confirmed in large-scale
parallelized calculations\cite{HNakano_kgm_gap_JPSJ2011,
HNakano_s1tri_LRO_JPSJ2013,HN_TSakai_kgm_1_3_JPSJ2014,
HN_TSakai_kgm_S_JPSJ2015,HN_YHasegawa_TSakai_dist_s
huriken_JPSJ2015,HN_TSakai_dist_tri_JPSJ2017,
HN_TSakai_tri_NN_JPSJ2017,HN_TSakai_kgm45_JPSJ2018,
YHasegawa_HN_TSakai_dist_shuriken_PRB2018,
TSakai_HN_ICM018,HN_TSakai_S2HaldaneGap_JPSJ2018}.
Note here that the largest-scale
calculations in this study
have been carried out using either the K computer or Oakforest-PACS.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg2.eps}
\end{center}
\caption{(Color)
Finite-size energy differences $\Delta$
of the cases $N_{\rm s}=16$, 20, 32, 36, and 40
for some representative $r$ as a function of $1/N_{\rm s}$.
Black circles, red pluses, blue crosses, and green squares
denote results for $r=1.5$, 0.72, 0.69, and 0.66, respectively.
}
\label{fig2}
\end{figure}

Now, let us observe the $N_{\rm s}$-dependence of $\Delta/J_{\rm 1}$
for some representative cases of $r$;
the results are depicted in Fig.~\ref{fig2}.
One finds that for $r=1.5$, $\Delta/J_{\rm 1}$ significantly decreases
as $N_{\rm s}$ is increased.
The decreasing behavior of $\Delta/J_{\rm 1}$ is consistent with
that in the gapless N$\acute{\rm e}$el-ordered phase.
Our results for $r=1.5$ suggest
an almost linear dependence on $1/N_{\rm s}$.
On the other hand, for $r=0.72$ and 0.69,
$\Delta/J_{\rm 1}$ decreases with increasing $N_{\rm s}$ for small $N_{\rm s}$
but shows only a very weak $N_{\rm s}$-dependence for large $N_{\rm s}$.
For $r=0.66$,
$\Delta/J_{\rm 1}$ finally becomes almost constant for all ranges
of $N_{\rm s}$.
From these observations, it is considered
that clusters with $N_{\rm s}=36$ and 40 capture well
the behavior of large systems approaching the thermodynamic limit.
Therefore, focusing our attention on the results of $N_{\rm s}=36$ and 40,
we hereafter investigate the behavior of the present system.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg3.eps}
\end{center}
\caption{
$r$-dependence of the spin gap
for $N_{\rm s} = 36$ and 40.
Open squares and closed circles denote results
for $N_{\rm s} = 36$ and 40, respectively.
Inset is a zoom-in view of the region of $r$
around $r\sim 0.7$.
}
\label{fig3}
\end{figure}
Next, let us observe the $r$-dependence
of the spin gap $\Delta$ for finite-size clusters in detail;
the results of $N_{\rm s}=36$ and 40 are depicted in Fig.~\ref{fig3}.
First, one finds, in the region up to $r\sim 0.67$,
that $\Delta/J_{1}$ gradually decreases as $r$ is increased
and
that data for the two sizes agree well with each other.
The agreement strongly suggests that
the system-size dependence has already become weak
and
that finite-size results almost agree
with the corresponding values for the thermodynamic limit.
Next, in the region from $r\sim 0.68$ to $r\sim 0.70$,
on the other hand, $\Delta/J_{1}$ gradually increases
with increasing $r$.
The good agreement of data for the two sizes is still maintained,
although the $r$-dependence of whether it increases or decreases
has been changed.
In the region above $r\sim 0.70$,
$\Delta/J_{1}$ decreases once but increases again
as $r$ is increased.
The upturn of $\Delta/J_{1}$ is observed
for both $N_{\rm s}=36$ and 40.
The significant characteristic in this region
is that there appears a considerable system-size dependence:
that is, $\Delta/J_{1}$ for $N_{\rm s}=40$ is smaller than
that for $N_{\rm s}=36$ at a given $r$.
To find whether the nonzero spin gap exists or
is absent in the thermodynamic limit,
we will need to carry out further analysis.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg4.eps}
\end{center}
\caption{(Color)
Energy-level scheme and second derivatives of the energies
with respect to the ratio $r$.
Panels (a) and (b) show results of the energy-level scheme
for $N_{\rm s} = 36$ and 40, respectively.
Panels (c) and (d) show results of the second derivatives
for $N_{\rm s} = 36$ and 40, respectively.
Squares and diamonds denote results for $M=0$ and $M=1$, respectively.
}
\label{fig4}
\end{figure}

To deepen our understanding of the behavior of $\Delta$
showing a complex dependence of decreases and increases,
next, let us observe the $r$-dependence of $E(N_{\rm s},M)$;
the results are depicted in Fig.~\ref{fig4}(a) and (b).
First, one easily finds that
$E(N_{\rm s},0)$ in the region up to $r\sim 0.67$ is constant.
The constant value corresponds to the eigenenergy
of the rigorous dimer ground state.
This clearly indicates that
this region corresponds to the dimer phase with the exact ground state.
One finds, from the results of $E(N_{\rm s},0)$ above $r \sim 0.675$, that
another state different from the dimer state becomes the ground state.
No significant size dependence is observed
with respect to the boundary ratio $r \sim 0.675$.
This result agrees well with previously known estimates
of the phase boundary:
$r=0.677$ in Ref.~\ref{AKoga_NKawakami_PRL2000},
$r=0.678$ in Ref.~\ref{Lauchli_PRB2002},
$r=0.687$ in Ref.~\ref{Lou_arXiv1212_1999}, and
$r=0.675$ in Ref.~\ref{Corboz_Mila_PRB2013}.
In order to capture the behavior in the region $r \simge 0.675$,
we evaluate a numerical second derivative given by
$-\delta^{2}[E(N_{\rm s},M)]/\delta r^{2}=[2E(N_{\rm s},M)|_{r}
-E(N_{\rm s},M)|_{r+\delta r}-E(N_{\rm s},M)|_{r -\delta r}]/(\delta r)^{2}$;
the results are depicted in Fig.~\ref{fig4}(c) and (d).
It is known that the analysis based on second derivatives
is useful to detect the boundaries
of a target system\cite{You_second_deri_PRE2007,
Albuquerque_second_deri_PRB2010,HN_TSakai_tri_NN_JPSJ2017}.
In Fig.~\ref{fig4}(c) and (d),
the second derivatives can appropriately capture
the discontinuity at $r \sim 0.675$ for both $M=0$ and $M=1$.
In addition, the second derivatives
for $N_{\rm s}=36$
show another discontinuity around $r\sim 0.69$-0.70.
For $N_{\rm s}=40$, the second derivative of $M=1$ also shows
a discontinuity
around $r\sim 0.70$;
that of $M=0$ does not show a discontinuity but
it shows a peak at $r\sim 0.71$ instead.
The behavior around $r\sim 0.70$ is consistent with
the result $r=0.702$ for the edge of the plaquette-singlet phase
reported in the $N_{\rm s}=32$ diagonalization study
in Ref.~\ref{Lauchli_PRB2002}.
On the other hand, the present result of $r\sim 0.70$
differs from
$r=0.86$ in Ref.~\ref{AKoga_NKawakami_PRL2000},
$r=0.75$ in Ref.~\ref{Lou_arXiv1212_1999}, and
$r=0.765(15)$ in Ref.~\ref{Corboz_Mila_PRB2013}
as results for the edge of the plaquette-singlet phase.
The difference will be examined later.
From Fig.~\ref{fig4}, therefore, one can understand that
the observed changes in the dependence of $\Delta$
are due to the energy-level structure
of both $E(N_{\rm s},0)$ and $E(N_{\rm s},1)$.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg5.eps}
\end{center}
\caption{
$r$-dependence of the product of the system size and the spin gap
for $N_{\rm s} = 36$ and 40.
Open squares and closed circles denote results
for $N_{\rm s} = 36$ and 40, respectively.
}
\label{fig5}
\end{figure}
Next, let us examine the system-size dependence of the spin gap
in the region of large $r$.
When $r$ is infinitely large, the system is reduced to
the simple square-lattice antiferromagnet, showing that
the spin excitation is gapless
owing to the existence of the N$\acute{\rm e}$el order.
As a means of distinguishing
whether the system is gapped or gapless,
the method of observing the product of the system size and the spin gap
is known.
This method was successfully used in the study of
the plateau \--- the gap under the magnetic field \---
at the one-third height of the saturation
in the triangular-lattice Heisenberg antiferromagnet
with next-nearest-neighbor interactions\cite{HN_TSakai_tri_NN_JPSJ2017}.
The results of this analysis for the present system
with $N_{\rm s}=36$ and 40 are depicted in Fig.~\ref{fig5}.
One clearly finds that the results from the two sizes
agree with each other in the region down to $r\sim 0.75$.
When $r$ is further decreased, the results of $N_{\rm }=40$
clearly become larger than those of $N_{\rm }=36$.
The agreement in the behavior of $N_{\rm s} \Delta/J_{1}$
in the region of large $r$ suggests that  the finite-size spin gap
in this region exhibits $\Delta \propto 1/N_{\rm s}$,
which means that the system is gapless.
In the region below $r\sim 0.75$, on the other hand,
the finite-size spin gap does not have
the dependence $\Delta \propto 1/N_{\rm s}$.
Although the $N_{\rm s}$-dependence of $\Delta$
in the region between $r\sim 0.71$ and $r\sim 0.75$ is unclear
at the present stage, there are two possible scenarios.
One is that the system is gapped without any long-range orders.
The other is that the system is gapless, but
the $N_{\rm s}$-dependence of $\Delta$ is different
from $\Delta \propto 1/N_{\rm s}$
corresponding to the N$\acute{\rm e}$el-ordered phase.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg6.eps}
\end{center}
\caption{(Color)
Correlation functions $\langle S_{i}^{z} S_{j}^{z} \rangle$.
Panel (a) shows the positions of site $i$ and $j$.
For a given site $i$ denoted by the double circle,
$j$ is taken in three directions shown by squares,
triangles, and inversed triangles.
Panels (b), (c), and (d) show results
for $\langle S_{i}^{z} S_{j}^{z} \rangle$
as a function of the distance between $i$ and $j$
for $r=1.50$, 0.73, and 0.68, respectively.
Panels (e) and (f) show the $r$-dependence
of $\langle S_{i}^{z} S_{j}^{z} \rangle$
for the shortest-distance pair
along the direction represented by triangles.
Black and red symbols denote results
for $N_{\rm s} = 36$ and 40, respectively.
}
\label{fig6}
\end{figure}
To capture the change from the N$\acute{\rm e}$el-ordered phase
to the plaquette-singlet phase,
let us observe correlation functions in the ground state,
namely, $\langle S_{i}^{z} S_{j}^{z} \rangle$;
the results for both $N_{\rm s}=36$ and $N_{\rm s}=40$
are depicted in Fig.~\ref{fig6}.
The case $r=1.50$ in Fig.~\ref{fig6}(b) is a typical one
for the N$\acute{\rm e}$el-ordered phase.
All the results shown by triangles and inversed triangles
are positive; this feature is explained by the fact that
both $i$ and $j$ for a measured pair are in a common sublattice
among the two sublattices of the N$\acute{\rm e}$el-ordered state.
The results shown by the triangles and the inversed triangles
also indicate a gradual decay as the distance is increased.
The results shown by the squares indicate alternating signs.
This behavior suggests the staggered nature
of the N$\acute{\rm e}$el-ordered state.
Therefore, the characteristics of the N$\acute{\rm e}$el-ordered state
are well captured in Fig.~\ref{fig6}(b).
The results in Fig.~\ref{fig6}(d), on the other hand,
are completely different from those in Fig.~\ref{fig6}(b).
Among the results shown by the triangles and the inversed triangles,
only the shortest-distant datum by the inversed triangle is positive,
and the rest are negative.
Absolute values of $\langle S_{i}^{z} S_{j}^{z} \rangle$
for distances larger than two are very small.
These behaviors of the correlation functions
are different from those of the N$\acute{\rm e}$el-ordered state,
but are consistent with those of the plaquette-singlet state.
In this state, each plaquette singlet is located at a local square
involving the $J_{1}$ bond and is the one that has a component
of two-spin singlet in diagonal pairs of the square
among two possible singlet states of four spins.
In the results in Fig.~\ref{fig6}(c) for $r=0.73$,
the pattern of whether $\langle S_{i}^{z} S_{j}^{z} \rangle$ is
positive or negative is common with Fig.~\ref{fig6}(b)
and different from Fig.~\ref{fig6}(d).
A significant difference between Fig.~\ref{fig6}(c) and Fig.~\ref{fig6}(d)
is $\langle S_{i}^{z} S_{j}^{z} \rangle$
for the shortest-distant pair shown by the triangle;
its $r$-dependence is depicted in Fig.~\ref{fig6}(e) and (f).
In Fig.~\ref{fig6}(e),
the dependence reveals a discontinuity at $r\sim 0.675$
for $N_{\rm s}=36$ and 40;
in Fig.~\ref{fig6}(f), another discontinuity appears
for $N_{\rm s}=36$, which divides the region of $r$ into
negative $\langle S_{i}^{z} S_{j}^{z} \rangle$
and positive $\langle S_{i}^{z} S_{j}^{z} \rangle$ regions.
For $N_{\rm s}=40$ in Fig.~\ref{fig6}(f),
$\langle S_{i}^{z} S_{j}^{z} \rangle$ changes its sign
around  $r$ similar to that for the discontinuity of $N_{\rm s}=36$,
although $\langle S_{i}^{z} S_{j}^{z} \rangle$ for $N_{\rm s}=40$
is not discontinuous.
One possible scenario for the spin state in the larger-ratio region
is that,
if the system forms plaquette singlets located at a local square
involving the $J_{1}$ bond,
each plaquette singlet is
the other singlet state which does not include
a component of two-spin singlet in diagonal pairs of the square.
Therefore,
our calculations
suggest that the state for $r=0.73$ shows a behavior
that is different from that of the state for $r=0.68$
and that the behavior changes at common $r$ values
for $N_{\rm s}=36$ and 40.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg7.eps}
\end{center}
\caption{(Color)
$r$-dependence of correlation functions.
Squares and crosses denote results of $\langle S_{i}^{z} S_{j}^{z} \rangle$
for next-nearest-neighbor pair illustrated in left inset.
Diamonds and pluses denote results of $\langle S_{i}^{z} S_{j}^{z} \rangle$
for the longest-distance pair in the finite-size clusters.
The right inset shows the pair from the common centered site
by the double circle to the corner sites by green (red) closed triangles
of the dotted-line squares for $N_{\rm s}=36$ ($N_{\rm s}=40$).
Squares and diamonds are for $N_{\rm s}=36$;
crosses and pluses are for $N_{\rm s}=40$.
}
\label{fig7}
\end{figure}

To find out
whether or not the N$\acute{\rm e}$el-type long-range order survives,
next, let us observe the $r$-dependence
of correlation functions $\langle S_{i}^{z} S_{j}^{z} \rangle$ in detail;
the results are depicted in Fig.~\ref{fig7}.
One finds that as $r$ is decreased down to $r\sim 0.7$,
$\langle S_{i}^{z} S_{j}^{z} \rangle$ for the next-nearest-neighbor pair
gradually decreases, but its magnitude is not so small.
On the other hand, $\langle S_{i}^{z} S_{j}^{z} \rangle$
for the pair between the longest distance decreases more rapidly;
its magnitude becomes considerably small in the region below $r\sim 0.8$.
To capture the difference of $\langle S_{i}^{z} S_{j}^{z} \rangle$
between the next-nearest-neighbor pair and the longest-distance pair,
we examine $R_{\rm cf}$, defined as
the ratio of $\langle S_{i}^{z} S_{j}^{z} \rangle$ for $N_{\rm s}=40$
divided by the corresponding $\langle S_{i}^{z} S_{j}^{z} \rangle$
for $N_{\rm s}=36$ presented in Fig.~\ref{fig7};
the results are depicted in Fig.~\ref{fig8}.
One finds that $R_{\rm cf}$ for the longest-distance pair
significantly decreases below $r\sim 0.76$,
whereas the ratio for the next-nearest-neighbor pair
is maintained at $R_{\rm cf}\sim 1$.
This observation suggests that
the N$\acute{\rm e}$el-type long-range order survives
in the region above $r\sim 0.76$
and that the order may disappear in the region below $r\sim 0.76$,
where
the N$\acute{\rm e}$el-type short-range correlations still survive.
As previous estimates for the edge of the N$\acute{\rm e}$el-ordered phase,
recall $r=0.75$ in Ref.~\ref{Lou_arXiv1212_1999},
and $r=0.765(15)$ in Ref.~\ref{Corboz_Mila_PRB2013};
the present result, $r\sim 0.76$ for the edge of the region
where the N$\acute{\rm e}$el-type long-range order definitely exists,
agrees well with those previous estimates.

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8cm]{hnkn-fg8.eps}
\end{center}
\caption{
Ratios of $\langle S_{i}^{z} S_{j}^{z} \rangle$ for $N_{\rm s}=40$
divided by $\langle S_{i}^{z} S_{j}^{z} \rangle$ for $N_{\rm s}=36$.
Crosses and closed circles denote
results for the next-nearest-neighbor pair and
for the longest-distance pair presented in Fig.~\ref{fig7}.
}
\label{fig8}
\end{figure}

In summary,
we have studied
the Shastry-Sutherland model
by the Lanczos-diagonalization method.
The present study has presented
diagonalization results for 36-site and 40-site
clusters that have not been reported before.
Our numerical results have successfully
clarified the dependence of the spin gap
on the ratio of interactions.
Our calculations have successfully
captured the edge of the dimer phase
to be $J_{2}/J_{1}\sim 0.675$
and
the edge of the N$\acute{\rm e}$el-ordered phase
to be $J_{2}/J_{1}\sim 0.76$.
A noteworthy finding is a third specific ratio $J_{2}/J_{1}\sim 0.70$
which divides the intermediate region
between the ratios of the two edges
into two parts.
We have found
from observation of correlation functions
that the spin state in the smaller-ratio region
and the one in the
larger-ratio region
are different from each other.
The properties of the ground states in the two intermediate regions
should be further studied
from different viewpoints in future.
Such studies would greatly contribute
to our fundamental understanding of frustrated magnetism.

\begin{acknowledgment}

We wish to thank
Professor N.~Todoroki
for fruitful discussions.
This work was partly supported
by JSPS KAKENHI Grant Numbers
16K05418, 16K05419, 16H01080 (JPhysics), and 18H04330 (JPhysics).
Nonhybrid thread-parallel calculations
in numerical diagonalizations were based on TITPACK version 2
coded by H. Nishimori.
In this research, we used the computational resources of the K computer
provided by the RIKEN Advanced Institute for Computational Science
through the HPCI System Research projects
(Project ID: hp170018, hp170028, and hp170070).
We used the computational resources
of Fujitsu PRIMERGY CX600M1/CX1640M1(Oakforest-PACS)
provided by
Joint Center for Advanced High Performance Computing
through the HPCI System Research project
(Project ID: hp170207 and hp180053).
Some of the computations were
performed using the facilities of
the Department of Simulation Science,
National Institute for Fusion Science;
Institute for Solid State Physics, The University of Tokyo;
and Supercomputing Division,
Information Technology Center, The University of Tokyo.

\end{acknowledgment}