# On the Collaborative Object Transportation Using Leader Follower Approach

Sumanta Ghosh1∗, Subhajit Nath1∗, Sarvesh Sortee1, Lokesh Kumar1, and Titas
Bera1 * Sumanta Ghosh and Subhajit Nath contributed equally to this
work.1Sumanta Ghosh, Subhajit Nath, Sarvesh Sortee, Lokesh Kumar, Titas Bera
with Robotics and Autonomous Systems Group, TCS Research, India
{g.sumanta1,nath.subhajit, sarvesh.sortee, kumar.lokesh7<EMAIL_ADDRESS>
###### Abstract

In this paper we address the multi-agent collaborative object transportation
problem in a partially known environment with obstacles under a specified goal
condition. We propose a leader follower approach for two mobile manipulators
collaboratively transporting an object along specified desired trajectories.
The proposed approach treats the mobile manipulation system as two independent
subsystems: a mobile platform and a manipulator arm and uses their kinematics
model for trajectory tracking. In this work we considered that the mobile
platform is subject to non-holonomic constraints, with a manipulator carrying
a rigid load. The desired trajectories of the end points of the load are
obtained from Probabilistic RoadMap-based planning approach. Our method
combines Proportional Navigation Guidance-based approach with a proposed Stop-
and-Sync algorithm to reach sufficiently close to the desired trajectory, the
deviation due to the non-holonomic constraints is compensated by the
manipulator arm. A leader follower approach for computing inverse kinematics
solution for the position of the end-effector of the manipulator arm is
proposed to maintain the load rigidity. Further, we compare the proposed
approach with other approaches to analyse the efficacy of our algorithm.

##  I INTRODUCTION

With the rapid pace of recent developments in robotics, its applications have
become extensive. Given the general advantages of cooperative mobile
manipulators over a single mobile manipulator; vast studies have been carried
out in recent years over the abilities of the same, mainly to accomplish
complex tasks like object transportation, assembling, packaging and welding[1,
2, 3]. Most of the tasks mentioned above are repetitive in nature and
performed in structured environment. Centralized control approaches are very
popular for solving this kind of problem though they suffered from high
computational complexities [4, 5, 6]. Decentralized and distributed control
strategies are proposed to mitigate this problem [7, 8, 9, 10]. A
decentralised approach to multi-robot manipulation in which a team of robots
encircles and traps an object and transports it by dragging and pushing to the
target configuration in an obstacle laden environment is presented in [11]. A
formalism of local motion planning for collaborative multi-robot manipulation
was developed in [12]. It exploits the deformability of soft objects during
manipulation. Formation based control approaches are very popular in
collaborative multi-agent load transportation.

Figure 1: Collaborative transportation of a rigid rod by two mobile
manipulators in a circular trajectory.

In [13], a sequential convex programming formulation is proposed to solve a
formation control problem for collaborative transportation of a rigid object
by two mobile manipulators in dynamic environment. Authors in [14] presented a
method for an optimal formation of a team of non-holonomic robots for a common
object manipulation. A distributed swarm algorithm with application in multi-
agent collaborative load transportation is developed in [15] where agents are
reaching an aggregate state maintaining the environment constraints. Authors
in [16] proposed a decentralized leader-follower hierarchical strategy for two
robots collaboratively transporting an object in a partially known environment
with obstacles.

In this work we consider the problem of object transportation using mobile
manipulators where two robots are transporting a rigid rod collaboratively
complying with object’s physical constraints in an environment with obstacles.
The desired trajectories of the end points of the object are generated using a
global sampling based path planner. We consider that the measure of success
while attempting to reach the desired goal location is not the the time taken
or the trajectory error, rather how close the distance between the end
effectors of the robots remains to the length of the load during the
execution, while staying as close to the desired trajectory as possible. We
propose a Stop-and-Sync algorithm, combined with Proportional Navigation
Guidance (PNG) based approach to steer the mobile manipulator close to desired
the trajectory, and a leader follower approach for solving inverse kinematics
to track the desired trajectory.

The main contribution of this paper is developing an efficient method for a
dual robot system to transport a rigid object from a specified initial
position to a final position, while avoiding collisions with obstacles
maintaining the synchronization between the agents.

The remainder of the paper is organized as follows, Section II introduces the
problem formulation and a overview of the method. In Section III we introduce
a Stop-and-Sync algorithm along with the PN Guidance based method. The
proposed leader follower approach is discussed at the later part of that
section. Section IV contains the simulation results. We conclude the paper
with discussion in Section V.

##  II Problem Formulation & Approach

We represent the position of a point $B$ measured from a point $A$, expressed
in frame $F$ as ${}^{A}\textbf{{p}}_{F}^{B}$. When the positions are expressed
in the world frame, for the ease of notation, we denote the position of point
$B$, measured from $A$, as ${}^{A}\textbf{{p}}^{B}$. Further, to represent the
position vector of a point $B$ measured from the origin of the world frame
expressed in the world frame we use $\textbf{{p}}^{B}$.

Figure 2:  Schematic Representation of load transport with two TurtleBot3
Waffle Pi mounted with the OpenMANIPULATOR-X arm. The specifications for the
same can be found here. Due to orientation constraints of the end-effector,
the rigid load is supported by two vertical fixtures connected by revolute
joints.

We consider a system of two mobile manipulators carrying a rigid cylindrical
load connected through their end-effector, as shown in Fig.(2). The aim is to
navigate the load from an initial position to a final while avoiding obstacles
and keeping the distance between the end-effectors same as the length of the
rigid load for all the time.

The robot which acts as the Leader is denoted by subscript $L$, while the
Follower is expressed by the subscript $F$. The desired end-effector
trajectory to be followed by each $a^{\text{th}}$ robot is denoted as:
$T_{a}$; whose projection on the $X-Y$ plane is denoted as $P_{a}$. We assume
that the desired end-effector trajectory for each robot has $n_{d}$ number of
way-points. These way-points, when projection onto the $X-Y$ plane, act as
target points for the corresponding mobile bases. Accordingly, these projected
desired target points are represented as:
$\textbf{d}^{a}_{P}=(x^{d,a}_{P},y^{d,a}_{P})$. The desired end-effector
positions at each time step are obtained by projection to find the nearest
point on $T_{a}$. We discuss more about it later in Section III. These
positions are represented as:
$\textbf{d}^{a}_{n}=(x^{T,a}_{n},y^{T,a}_{n},z^{T,a}_{n})$. The base locations
at each time step $k$ are given as:
$\textbf{b}^{a}_{k}=(x^{a}_{k},y^{a}_{k})$.

We assume that $T_{a}$ is generated from PRM-based path planner considering
the obstacles. We consider that the base of the mobile manipulators are non-
holonomic and the motion of the mobile platform is confined within XY plane
only. The mobile manipulator system consists of two main subsystems: the
platform or base of the mobile manipulator, and the manipulator arm. The goal
is to start the $a^{\text{th}}$ robot from a position $(x_{k}^{a},y_{k}^{a})$,
sufficiently close to $\textbf{d}^{a}_{P}$, compute the control inputs,
$(v^{a},\omega^{a})$ to the mobile base that lands the mobile manipulator to a
point $(x^{a}_{k+1},y^{a}_{k+1})$ within the reachability distance of the next
point $(x^{d}_{k+1},y^{d}_{k+1})$ on the desired trajectory. The reachability
distance of a point is defined as follows:

$\rho_{d}=\gamma*\rho_{l}$ (1)

Where, $\rho_{l}$ is the reachability of a manipulator defined as the volume
of the workspace that can be reached by the manipulator. $\gamma$ is a user
defined value that represents a tolerable fraction of the reachability.
Henceforth, the joint angles of the manipulator arm are calculated by solving
the inverse kinematics of the manipulator arm to the desired point
$\textbf{d}^{a}_{n}$ measured from the base of the mobile manipulator
expressed in the world frame depending on the role (leader or follower) of the
manipulator. The whole process is performed for every base position
$\textbf{b}^{a}_{k}$ to track the desired end-effector trajectory for a single
robot. Stop-and-Sync algorithm is introduced to maintain the collaboration in
between them.

###  II-A Path Planning

We consider that a 3D model of the environment with polyhedron obstacles is
present. The position, orientation, and the shape of all the static obstacles
are assumed to be known apriori.

(a)

(b)

Figure 3: Schematic of the nodes generated through PRM in (a), and the path
planning achieved in (b). The grayed areas represent the obstacles, while the
robots carrying the load are represented by the red squares.

For path planning, we use PRM on a 2D projection of the above stated obstacle
laden environment. To that end, given the initial and final configurations
$q_{init}$ and $q_{goal}$, a number of vertices are generated randomly over
$Q_{free}$ – which are then connected to its $c$ closest neighbours via edges
using a local planner. Henceforth, the shortest path between $q_{init}$ and
$q_{goal}$ through the obtained edges is extrapolated using Dijkstra’s
Algorithm. The execution of this algorithm in MATLAB is shown in Fig. 3.

###  II-B Forward Kinematics of Mobile Platform

For a mobile robot the forward kinematics is defined as the mapping between
the control variables (as a function of time, $t$) and the pose of the robot.
The approach to obtain a mapping between $\textbf{b}^{a}_{k}$ and the control
variables is the same as followed in []. In our case, the control variables
are the linear and the angular velocities of the mobile robot. Assuming the
trajectory covered by the mobile base from timesteps $k$ to $k+1$ to be
composed of small circular segments; for a differential drive robot starting
with a pose $[x(t),y(t),\theta(t)]^{T}$, moving with linear velocity $v(t)$
and angular velocity $\omega(t)$, the pose of the robot at $(t+\delta t)$ can
be determined from the forward kinematics,

$\begin{bmatrix}x(t+\delta t)\\\ y(t+\delta t)\\\ \theta(t+\delta
t)\end{bmatrix}=\textbf{{g}}^{B}(v(t),\omega(t))$ (2)

Figure 4: Schematic diagram of the Leader Follower Approach with Stop-and-Sync
Algorithm.

where the matrix $\textbf{g}^{B}$ was as derived in [17]. Appropriate limits
can be introduced to replace the circular trajectory with a straight line
path.

###  II-C Inverse Kinematics of Manipulator Arm

The forward kinematics of a manipulator provides the position and the
orientation of the end effector given the joint variables in our case the
joint angles, $\textbf{{r}}^{E}=\textbf{{f}}^{E}\left(\bm{\beta}\right)$. This
mapping from the joint space to the end effector position is highly nonlinear
in nature. The problem of solving the inverse mapping, from the end effector
position to the joint angles, is known as inverse kinematics of the
manipulator, $\bm{\beta}=\textbf{{f}}^{-1}\left(\textbf{{r}}^{E}\right)$.
Various analytic and numerical approaches exist for finding the solution of
the inverse kinematics problem. One efficient approach is to formulate the
problem as an optimization problem,

$\displaystyle\bm{\beta}^{\star}=\arg\min_{\bm{\beta}}\hskip
14.22636pt\left\|{}^{B}\textbf{{f}}^{E}(\bm{\beta})-^{B}\textbf{{r}}^{E}\right\|^{2}_{2}$
(3) $\displaystyle\hskip 22.76228pt\text{subject to
}\bm{\beta}_{min}\leq\bm{\beta}\leq\bm{\beta}_{max}$ (4) $\displaystyle\hskip
65.44142pt\dot{\bm{\beta}}_{min}\leq\dot{\bm{\beta}}\leq\dot{\bm{\beta}}_{max}$
(5)

The above optimization problem is non convex in nature. Sequential quadratic
programming [18] approach can be used to find the solution of this nonlinear
optimization problem. Finding a good solution of this optimization problem is
very hard as it may often get stuck at the local optima and therefore can fail
to find a feasible solution of the optimization problem if it exists. But this
approach gives more freedom to handle the joint limits, joint velocity limits
(defined in Eq.(3)-Eq.(4)) and other conditions by adding constraints in the
optimization problem. In our work we exploited this characteristic.

##  III Collaborative Object Transportation

###  III-A Trajectory Tracking

Trajectory Tracking for each of the robots is first carried out independently
as presented here. In this stage, it is considered that the robots ($R_{1}$
and $R_{2}$) do not have any information about each other’s position and the
collaborative object transportation problem is divided into two separate
planning problems for each robot.

Numerous approaches have been tried out over the years to attain trajectory
tracking of mobile manipulators. We focus our attention on one such approach,
which has proven to be simple, yet effective. The Proportional Navigation
Guidance Law[19]; widely used in homing missile applications, is used to
predict the possibility of collision of two objects based on the direct Line-
of-Sight between them. When the velocity vector of the pursuer changes at the
rate of rotation of the target, a collision is predicted.

Given that each of the mobile bases starts at their first respective projected
desired point, the subsequent desired point acts as a stationary target for
the mobile base at that particular time-step. Accordingly, the requisite Line
of Sight (LOS) is calculated at each time-step to update the angular velocity
$\omega$ of the base, while keeping the forward velocity, $v$ constant. As the
mobile base finally enters within the radius of the aforementioned
reachability distance $\rho_{d}$ of the current target point, the target
switches to the next projected desired point on the trajectory. This iterative
procedure is carried out until the robots reach their final destination.

The rate of change of LOS is taken as $\dot{\lambda}$. With respect to the
Proportional Guidance Law, we may then write:

$\alpha_{n}=N\dot{\lambda}v$ ; where $\alpha_{n}$ is the generated
acceleration perpendicular to the instantaneous velocity vector of the mobile
base, and N is a proportionality constant, taken as 6 in our case.

Assuming that the robots cover a circular arc trajectory at each time-step
(given $a_{n}\neq 0$), the generated acceleration acts with acentripetal
nature on the mobile base. Accordingly, the value of $\omega$ at that timestep
can be generated as:

$\omega=\alpha_{n}/v$

Having thus obtained the values: [v, $\omega$] for the current timestep, the
forward kinematics routine for the non-holonomic base, as explained before, is
applied to obtain the next position $\textbf{b}^{a}_{k}$ for each robot.

###  III-B Collaborative Load Transportation

Assuming that both the desired trajectories have same number of points
$n_{d}$, the mobile base has to traverse repeatedly at each time step, until
it enters the projected circle of reachability distance $\rho_{d}$ of a
particular target point, as defined previously.

This approach may cause each robot to generate different number of base
positions: $\textbf{b}^{a}$ between their respective subsequent projected
points within each desired base trajectory– especially if the distances
between the respective target points are significantly dissimilar in the
defined base trajectories. Assuming similar velocity limits in both robots,
this causes a loss of correspondence between the robots– thus violating the
load rigidity constraint in our assumptions, given that the bases move so far
apart as to hinder manipulation by the manipulator arms. To account for this
discrepancy, we adopt a _“Stop-and-Sync”_ Algorithm to morph the generated
velocity and base locations at each time step. This algorithm is outlined in
further detail below.

####  III-B1 Stop-and-Sync Algorithm

For each base position obtained by the robots at a particular time step, the
nearest point on the desired end-effector trajectory curve is obtained through
a projections of $\textbf{b}^{a}_{k}$ on $T_{a}$. Here, the desired end-
effector trajectory $T_{a}$ is approximated into a set of lines through linear
interpolation between the desired points. Henceforth, perpendicular
projections are drawn from the current base position to the aforementioned
interpolated curve. The computed intersection point between the projected
perpendicular, and the interpolated line then acts as the candidate point for
solving the inverse kinematics of the manipulator arm with respect to the
current base position; at the current time step.

When a robot base enters within the defined reachability radius, the
corresponding target point $p$ shifts to the next point $p+1$ on the desired
trajectory. At this time-step, the base location of the other robot is
checked. If the target point of the other robot has not changed in a similar
fashion yet, the former robot base is stopped; until the target point of the
other robot changes – at which point of time, the two robots start moving
together again towards their respective new target points. This approach
necessitates that the two robots start moving synchronously after target
switching occurs from the point $p$ to $p+1$ on both trajectories $P_{a}$. The
Stop-and-Sync algorithm is summarized in Algorithm - 1.

Algorithm 1 Stop-and-Sync Algorithm

1:Input:
$b^{1}_{0},b^{2}_{0}$: initial base positions of robots
$d^{1},d^{2}$: discretized points of desired base trajectories $B_{1}$ and
$B_{2}$

2:Output:
$V_{1},V_{2}$: generated velocities of the robots
$b^{1},b^{2}$: generated base positions of robots

3:$R_{tol}=\gamma*\rho_{l}$

4:$stop_{1}=0$ $stop_{2}=0$

5:$n_{d.p}=size(d_{1})$

6:while $p\neq n_{d.p}+1$ do

7: for each timestep $k$ do

8: for each agent $a$ do

9: if $stop_{a}\neq 1$ then

10: Calculate velocity: $V^{a}_{k}=[v^{*},\omega^{*}]$, and the base position:
$b^{a}_{k+1}$

11: else

12: $V^{a}_{k}=[0,0]$

13: $b^{a}_{k+1}$ = $b^{a}_{k}$

14: end if

15: $dist_{a}=\lVert{d_{p}^{a}}-b^{a}_{k+1}\rVert$

16: end for

17: if $dist_{1}\leq R_{tol}\And dist_{2}\leq R_{tol}$ then

18: $p=p+1$

19: $stop_{1}=0$

20: $stop_{2}=0$

21: else if $dist_{1}\leq R_{tol}\And dist_{2}>R_{tol}$ then

22: $stop_{1}=1$

23: else if $dist_{1}>R_{tol}\And dist_{2}\leq R_{tol}$ then

24: $stop_{2}=1$

25: end if

26: end for

27:end while

####  III-B2 Leader Follower Approach

In collaborative object transportation with rigid load, the load rigidity
constraint should be given more importance than the tracking error. We propose
a leader follower approach that uses the information of the end effector
position of one of the mobile manipulator (leader) to solve the inverse
kinematics of the other manipulator arm (follower) to satisfy the load length
condition, i.e. the Euclidean distance between the estimated end effector
position of the leader and the end effector of the follower with respect to
the world frame should be equal to the load length, $l$.

$\left\|\hat{\textbf{{r}}}^{E_{L}}-\textbf{{f}}\left(\bm{\beta}_{F}\right)\right\|_{2}=l$
(6)

When the two mobile manipulators reach inside the reachability distance, the
leader robot, $R_{L}$ computes the joint angles by solving the inverse
kinematics problem define in Eq. (3)-Eq. (5). Further, the end effector
position of the leader, $\hat{\textbf{{r}}}^{E_{L}}$ can also be estimated.
Then the joint angles of the follower can be determined by solving the
following constraint optimization problem,

$\displaystyle\bm{\beta}^{*}_{F}=\hskip
5.69046pt\arg\min_{\bm{\beta}_{F}}\hskip
11.38092pt\left\|{}^{B_{F}}\textbf{{f}}^{E_{F}}\left(\bm{\beta}_{F}\right)-^{B_{F}}\textbf{{r}}^{E_{F}}\right\|^{2}_{2}$
(7) $\displaystyle\hskip 25.6073pt\text{subject to}\hskip 8.5359pt\left|\hskip
2.84544pt\left\|\hat{\textbf{{r}}}^{E_{L}}-\textbf{{f}}\left(\bm{\beta}_{F}\right)\right\|_{2}-l\right|=0$
(8) $\displaystyle\hskip
73.97733pt\bm{\beta}_{min}\leq\bm{\beta}\leq\bm{\beta}_{max}$ (9)
$\displaystyle\hskip
73.97733pt\dot{\bm{\beta}}_{min}\leq\dot{\bm{\beta}}\leq\dot{\bm{\beta}}_{max}$
(10)

The leader follower algorithm is summarized in Algorithm -2. The schematic
diagram of the proposed approach is shown in fig. 4.

Algorithm 2 Leader Follower Approach

1:Input: $\bm{\beta}_{min}$, $\bm{\beta}_{max}$, $\dot{\bm{\beta}}_{min}$,
$\dot{\bm{\beta}}_{max}$, ${}^{B_{L}}\textbf{{r}}^{E_{L}}$,
${}^{B_{F}}\textbf{{r}}^{E_{F}}$, $l$.

2:Output: $\bm{\beta}^{\star}_{L}$,$\bm{\beta}^{\star}_{F}$

3:At $k+1$, $\forall k=0,1,\ldots,N-1$;

4:Leader: Determine $\bm{\beta}^{*}_{L}$ solving $(\ref{IK})-(\ref{IK2})$;

5:Leader: Estimate ${}^{B_{L}}\hat{\textbf{{r}}}^{E_{L}}$;

6:Follower: Determine $\bm{\beta}^{\star}_{F}$ from $(\ref{ik3})-(\ref{ik6})$.

Figure 5: The desired trajectories and their projection on the $X-Y$ plane of
the two robots.

##  IV Simulation and Results

###  IV-A Simulation Setup

The proposed algorithm is tested on a simulated environment in MATLAB to
capture their efficiency. For the mobile manipulator system, we use the
specifications for the TurtleBot3 Waffle Pi with the OpenMANIPULATOR-X
attached atop it. The pair of robots is assumed to carry a slender load of
length 0.65 m. The control variables $[v,\omega]$ are passed on to the robots
with a time gap, $\Delta t$ of $0.08$ seconds. For the manipulator arm we
incorporate the following constraints on joint angles $(\bm{\beta})$ in rad
and rate of change of the joint angles $\dot{\bm{\beta}}$ in rad/sec.:
$[-0.9\pi,-0.57\pi,-0.3\pi,-0.5\pi]^{\intercal}\leq\bm{\beta}^{\intercal}\leq[0.9\pi,0.5\pi,0.44\pi,0.65\pi]^{\intercal}$
and $\dot{\bm{\beta}}_{i}=4\hskip 8.5359pt\forall i\in\\{1,2,3,4\\}$. We
further limit the gripper orientation of the manipulator arm such that it
always remains parallel to the $X-Y$ plane while carrying the load for
hardware limitation.

Since most generated trajectories can be broken down into a straight line and
circular arcs, we have, for generality, considered the desired trajectories of
the two robots as shown in fig. 5 for simulation. The desired trajectory of
each robot has three parts, a semicircle in the middle and two straight lines
of length 1m at the both end of the semicircle. The radius of the semicircle
for inner and outer robot is taken 1 m and 1.65 m respectively. The height
along the z axis of each trajectory is taken as 0.2m. The initial positions of
the base of the mobile manipulators are $(1,-1)$ and $(1.65,-1)$. The goal
positions of the them are $(-1,-1)$ and $(-1.65,-1)$. The reachability
tolerance value $\gamma$ is taken as 0.4. The time delay in the system between
computing the solution of inverse kinematics of the leader and the follower is
considered negligible.

(a)

(b)

(c)

Figure 6: The obtained distances between the end-effectors of the two robots
for the three different methods. Fig. 6(a), fig. 6(b) and fig.6(c) represent
the results from PN Guidance, RRT generation and MPC based methods
respectively.

(a)

(b)

(c)

Figure 7: The trajectory tracking errors of two robots for the three methods.
Fig. 7(a), fig. 7(b) and fig. 7(c) represents the result from PN Guidance, RRT
generation and MPC based methods respectively

(a)

(b)

(c)

Figure 8: The linear velocities of the two robots for the three different
methods. Fig. 8(a), fig. 8(b) and fig.8(c) represent the results from PN
Guidance, RRT generation and MPC based methods respectively.

(a)

(b)

(c)

Figure 9: The angular velocities of the two robots for the three different
methods. Fig. 9(a), fig. 9(b) and fig.9(c) represent the results from PN
Guidance, RRT generation and MPC based methods respectively.

Figure 10: The joint angles of the two robots for the three methods. The top
plot portrays the joint angles of robot 1 while the bottom plot displays the
same for robot 2. Fig.9a, fig.9b and fig.9c represent the results for PNG
based, RRT generation based and MPC based approaches respectively.

###  IV-B Comparative Approaches

The presented approach is brought into comparison with two other approaches,
to compare efficacy of the results with respect to the ascertained parameters.

In the first approach, an attempt is made to generate a single Rapidly-
exploring random tree (RRT) from the given start node to the final goal node,
with the intermediate discretized points $d^{a}_{P}$ acting as the
transitional target nodes at each timestep. This process is carried out
iteratively until the tree spans the entire trajectory $T^{a}$.

Given a point $(x^{d}_{p},y^{d}_{p})$ on the desired trajectory, each robot
first generates 500 random sets of linear and angular velocities: $[v,\omega]$
within a specified range. Each velocity is registered upon the current base
location of each robot, and the corresponding final position is generated upon
the application of the velocities for time $\Delta t$, all the while
respecting the non-holonomic constraints of the robots. Among the generated
subsequent locations, the one nearest to the aforementioned desired point
$d^{a}_{P}$; and the equivalent velocity data: $[v^{*},\omega^{*}]$ is chosen
for the current timestep.

The computation of the control variables for the agent $a$, can thus, be
formulated as an optimization problem that minimizes the following objective
function:

$J=\left|d-\rho_{d}\right|$ (11)

where:

$d=\sqrt{(x^{d,a}_{p}-x^{a}_{k})^{2}+(y^{d,a}_{p}-y^{a}_{k})^{2}}$ (12)

The objective function indicates the control velocities needed to reach within
the reachability distance of the next desired point. The optimization problem
is defined below,

$v^{\star},\omega^{\star}=\arg\min_{v,\omega}J$ (13)

Once path generation is achieved in this manner, prior approaches of the Stop-
and-Sync Algorithm, and the Leader-Follower are incorporated to bring about
collaborative load transport by the robots.

The presented approach is also compared with a Sequential Linear Quadratic
Model Predictive [20] approach in simulation wherein each robot individually
tracks their own reference trajectory based on a desired object trajectory.
The model used is given by

$\dot{\boldsymbol{x}}=\begin{bmatrix}v\cos{\theta}\\\ \omega\sin{\theta}\\\
\dot{\boldsymbol{\beta}}\end{bmatrix}$ (14)

where the state
$\boldsymbol{x}\triangleq[x,y,\theta,\boldsymbol{\beta}^{T}]^{T}\in\mathbb{R}^{3+n}$
and
$\boldsymbol{u}\triangleq[v,\omega,\dot{\boldsymbol{\beta}^{T}}]^{T}\in\mathbb{R}^{2+n}$
are the linear and angular velocities for the base and joint velocities
respectively. For each robot, the cost function considered over a horizon of
length $N$ is

$J=\boldsymbol{ee}(\boldsymbol{x}_{N})^{T}\boldsymbol{Q}\boldsymbol{ee}(\boldsymbol{x}_{N})+\sum_{t=1}^{N-1}\boldsymbol{ee}(\boldsymbol{x}_{t})^{T}\boldsymbol{Q}\boldsymbol{ee}(\boldsymbol{x}_{t})+\boldsymbol{u}_{t}^{T}\boldsymbol{R}\boldsymbol{u}_{t}$
(15)

where the function
$\boldsymbol{ee}(\boldsymbol{x}_{t})=\boldsymbol{FK}(\boldsymbol{x}_{t})-\boldsymbol{r}(t)$
is the tracking error between the end effector position
$\boldsymbol{FK}(\boldsymbol{x}_{t})\in\mathbb{R}^{3}$ and reference
trajectory $\boldsymbol{r}(t)$ at time $t$. The matrices $\boldsymbol{Q}$ and
$\boldsymbol{R}$ are positive diagonal weighting matrices. Input and joint
angle constraints are incorporated in a soft manner by adding relaxed barrier
functions [21] to the cost equation. The parameters of all the stated
approaches are kept similar to the extent possible, for proper comparison.

###  IV-C Simulation Results

The simulation results are depicted in fig. 6 -10. The time taken to reach the
goal state from the initial position differs for the three approaches. Fig. 6
represents the distance between the end effectors of the mobile
manipulators$(l)$. This indicates the extension and compression the rigid load
may observe while the mobile manipulators are in action. The proposed leader
follower with stop and wait approach ensures minimum deviation of $l$ from the
desired load length $0.65$m which can also be seen in fig. 6(a) and fig. 6(b).
For SLQ MPC approach distortion of $l$ can be observed in fig. 6(c).

Trajectory tracking errors for the three methods are depicted in Fig. 7. In
our work, the Leader follows the trajectory $T_{L}$ independently, while the
Follower tries to remain close to the trajectory $T_{F}$, while giving greater
priority to maintaining a constant $l$. It is expected to observe tracking
error for the Follower as seen in fig. 7(a) though the tracking error for the
Leader is minimal. For RRT generation based approach it can be seen that the
tracking error of Leader is comparable to that of follower in fig. 7(b) though
both of them return better tracking error than the SLQ MPC approach.

The control input to mobile base $(v,\omega)$ is shown in fig. 13-fig.9. Due
to the Stop-and-Sync algorithm fluctuations can be noticed both in $v$ and
$\omega$ for PN Guidance and RRT generation based approach while for MPC based
approach it appears to be smooth. For the PN Guidance based approach, the
control variable $v$ remains constant for Robot-2 (which has to traverse a
longer distance for the outer semicircular path), while that of Robot-1 is
interrupted at regular intervals to accommodate synchronization with Robot-2.
Due to difference in orientations at the end of the circular arc trajectory,
both robots undergo seemingly unexpected intervals of variations in $v$ during
the second straight line phase. As Robot-1 has to navigate a tighter curve
with smaller radius, its values of $\omega$ are consistently higher. For RRT
generation based approach the variations in $v$ due to _Stop-and-Sync_ can be
spotted from the very beginning in both robots, due to randomizations in
velocity generation. The evaluation of the joint angles are shown in fig. 10
for the three methods.

###  IV-D Hardware Experimentations

The experiments run on the TurtleBot 3 Waffle Pi hardware can be accessed at
this playlist. Here, two robots are made to work collaboratively to transport
a rigid rod while following specified trajectories. Due to inherent hardware
constraints, a horizontal constraint is imposed upon the end-effectors to
prevent undue twisting of the load. As can be seen from the videos, the robots
were able to carry out load transportation successfully on an eclectic variety
of trajectories.

##  V CONCLUSIONS

In this paper, we have considered two robots collaboratively transporting a
rigid load, while subject to non-holonomic constraints. To that end, a Leader-
Follower Approach for the inverse kinematics solution; coupled with a Stop-
and-Sync Algorithm to keep the mobile bases in close proximity has been
proposed. The proposed approach, as can be ascertained from the results, shows
better results with respect to our salient metrics of load length, and
trajectory tracking error. Hardware experimentations also prove the efficacy
of the simulation results.

The proposed approach, due to its inherent nature, provides discontinuous
profiles of velocities. In future works, it shall be attempted to smoothen out
these profiles through a predictive localization technique. Interchangeability
of the agents acting out as Leader and Follower shall also be attempted to
obtain better results during trajectory tracking of tight curves.

## References

  * [1] Z. Feng, G. Hu, Y. Sun, and J. Soon, “An overview of collaborative robotic manipulation in multi-robot systems,” _Annual Reviews in Control_ , vol. 49, pp. 113–127, 2020.
  * [2] O. Khatib, K. Yokoi, K. Chang, D. Ruspini, R. Holmberg, and A. Casal, “Coordination and decentralized cooperation of multiple mobile manipulators,” _Journal of Robotic Systems_ , vol. 13, no. 11, pp. 755–764, 1996.
  * [3] P. Song and V. Kumar, “A potential field based approach to multi-robot manipulation,” in _Proceedings 2002 IEEE International Conference on Robotics and Automation (Cat. No. 02CH37292)_ , vol. 2. IEEE, 2002, pp. 1217–1222.
  * [4] C. Alford and S. Belyeu, “Coordinated control of two robot arms,” in _Proceedings. 1984 IEEE international conference on robotics and automation_ , vol. 1. IEEE, 1984, pp. 468–473.
  * [5] D. Rus, B. Donald, and J. Jennings, “Moving furniture with teams of autonomous robots,” in _Proceedings 1995 IEEE/RSJ International Conference on Intelligent Robots and Systems. Human Robot Interaction and Cooperative Robots_ , vol. 1. IEEE, 1995, pp. 235–242.
  * [6] K. Böhringer, R. Brown, B. Donald, J. Jennings, and D. Rus, “Distributed robotic manipulation: Experiments in minimalism,” in _Experimental Robotics IV: The 4th International Symposium, Stanford, California, June 30–July 2, 1995_. Springer, 1997, pp. 11–25.
  * [7] B. Donald, L. Gariepy, and D. Rus, “Distributed manipulation of multiple objects using ropes,” in _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No. 00CH37065)_ , vol. 1. IEEE, 2000, pp. 450–457.
  * [8] A. Marino, “Distributed adaptive control of networked cooperative mobile manipulators,” _IEEE Transactions on Control Systems Technology_ , vol. 26, no. 5, pp. 1646–1660, 2017.
  * [9] G. Habibi, Z. Kingston, W. Xie, M. Jellins, and J. McLurkin, “Distributed centroid estimation and motion controllers for collective transport by multi-robot systems,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2015, pp. 1282–1288.
  * [10] E. L. Zhu, Y. R. Stürz, U. Rosolia, and F. Borrelli, “Trajectory optimization for nonlinear multi-agent systems using decentralized learning model predictive control,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 6198–6203.
  * [11] J. Fink, M. A. Hsieh, and V. Kumar, “Multi-robot manipulation via caging in environments with obstacles,” in _2008 IEEE International Conference on Robotics and Automation_. IEEE, 2008, pp. 1471–1476.
  * [12] J. Alonso-Mora, R. Knepper, R. Siegwart, and D. Rus, “Local motion planning for collaborative multi-robot manipulation of deformable objects,” in _2015 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2015, pp. 5495–5502.
  * [13] J. Alonso-Mora, S. Baker, and D. Rus, “Multi-robot navigation in formation via sequential convex programming,” in _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 4634–4641.
  * [14] A. Abbaspour, K. Alipour, H. Z. Jafari, and S. A. A. Moosavian, “Optimal formation and control of cooperative wheeled mobile robots,” _Comptes Rendus Mécanique_ , vol. 343, no. 5-6, pp. 307–321, 2015.
  * [15] R. F. Carpio, L. Di Giulio, E. Garone, G. Ulivi, and A. Gasparri, “A distributed swarm aggregation algorithm for bar shaped multi-agent systems,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 4303–4308.
  * [16] M. Bujarbaruah, Y. R. Stürz, C. Holda, K. H. Johansson, and F. Borrelli, “Learning environment constraints in collaborative robotics: A decentralized leader-follower approach,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 1636–1641.
  * [17] K. S. Chong and L. Kleeman, “Accurate odometry and error modelling for a mobile robot,” in _Proceedings of International Conference on Robotics and Automation_ , vol. 4. IEEE, 1997, pp. 2783–2788.
  * [18] P. T. Boggs and J. W. Tolle, “Sequential quadratic programming,” _Acta numerica_ , vol. 4, pp. 1–51, 1995.
  * [19] S. A. Murtaugh and H. E. Criel, “Fundamentals of proportional navigation,” _IEEE spectrum_ , vol. 3, no. 12, pp. 75–85, 1966.
  * [20] M. Neunert, C. de Crousaz, F. Furrer, M. Kamel, F. Farshidian, R. Siegwart, and J. Buchli, “Fast nonlinear model predictive control for unified trajectory optimization and tracking,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 1398–1404.
  * [21] C. Feller and C. Ebenbauer, “Relaxed logarithmic barrier function based model predictive control of linear systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 3, pp. 1223–1238, Mar. 2017. [Online]. Available: https://doi.org/10.1109/tac.2016.2582040


# Supergranules as Probes of Solar Convection Zone Dynamics

David H. Hathaway  NASA Marshall Space Flight Center, Huntsville, AL 35812 USA
<EMAIL_ADDRESS>
###### Abstract

Supergranules are convection cells seen at the Sun’s surface as a space
filling pattern of horizontal flows. While typical supergranules have
diameters of about 35 Mm, they exhibit a broad spectrum of sizes from $\sim
10$ Mm to $\sim 100$ Mm. Here we show that supergranules of different sizes
can be used to probe the rotation rate in the Sun’s outer convection zone. We
find that the equatorial rotation rate as a function of depth as measured by
global helioseismology matches the equatorial rotation as a function of
wavelength for the supergranules. This suggests that supergranules are
advected by flows at depths equal to their wavelengths and thus can be used to
probe flows at those depths. The supergranule rotation profiles show that the
surface shear layer, through which the rotation rate increases inward, extends
to depths of $\sim 50$ Mm and to latitudes of at least $70\arcdeg$. Typical
supergranules are well observed at high latitudes and have a range of sizes
that extend to greater depths than those typically available for measuring
subsurface flows with local helioseismology. These characteristics indicate
that probing the solar convection zone dynamics with supergranules can
complement the results of helioseismology.

Sun: convection, Sun: rotation

##  1 INTRODUCTION

Supergranules were discovered in the 1950s by Hart (1954) but were best
characterized in the 1960s by Leighton et al. (1962) who gave them their name.
Leighton et al. (1962) showed that this cellular pattern of horizontal flows
covers the solar surface and that the boundaries of the cells coincide with
the chromospheric/magnetic network. Typical supergranules have diameters of
$\sim 35$ Mm and maximum flow speeds of $\sim 500$ m s-1. While the kinetic
energy spectrum has a distinct peak at wavelengths of $\sim 35$ Mm, the
spectrum includes cells at least three times larger and extends to much
smaller cells where the supergranule spectrum blends into the granulation
spectrum (Hathaway et al., 2000).

Larger cells live longer than smaller cells. Typical supergranules with
diameters of $\sim 30$ Mm live for $\sim 24$ hr (Simon & Leighton, 1964; Wang
& Zirin, 1989). Typical granules with diameters of $\sim 1$ Mm only live for
$\sim 5$ min (Title et al., 1989). Cells of intermediate size ($\sim 5-10$ Mm)
have intermediate lifetimes of $\sim 2$ hr (November et al., 1981).

The rotation rate of the supergranule pattern was first measured by Duvall
(1980) who cross-correlated the Doppler velocity pattern from equatorial
spectral scans obtained over several days. He found that the pattern rotates
about 3% faster than the photospheric plasma and faster rates are found for
the 24-hr time lags from day-to-day than for the 8-hr time lags from the
beginning to end of an observing day. He concluded that larger cells dominate
the longer time lags and that the observations are consistent with
supergranules embedded in a surface shear layer in which the rotation rate
increases with depth.

The presence of this shear layer was first suggested by Foukal & Jokipii
(1975) as a consequence of the conservation of angular momentum by convective
elements moving inward and outward in the near surface layers. Global
helioseismology inversions (Thompson et al., 1996; Schou et al., 1998)
indicate that the shear layer extends to depths of 35-50 Mm near the equator
but may disappear or reverse at latitudes above $\sim 55\arcdeg$. While local
helioseismology (Basu et al., 1999; Corbard & Thompson, 2002) does not probe
as deeply, it produces similar results which suggest that the rotation rate
may not continue to increase inward at higher latitudes.

Beck & Schou (2000) measured the rotation of the supergranule pattern using a
Fourier technique with space-based Doppler data from the ESA/NASA Solar and
Heliospheric Observatory (SOHO) Michelson Doppler Imager (MDI) (Scherrer et
al., 1995). They mapped the data onto heliographic coordinates, took the
Fourier transform in longitude of data from equatorial latitudes, and then the
Fourier transform in time of those spectral coefficients over 6 10-day
intervals in 1996 which had continuous coverage at a 15-min cadence. They
found that the larger cells do indeed rotate more rapidly than the smaller
cells, but with rotation rates that exceeded the peak internal rotation rate
at the base (50 Mm depth) of the surface shear layer as determined from global
helioseismology (Schou et al., 1998). This discovery led them to conclude that
supergranules must have wave-like properties in order to rotate faster than
the flows they are embedded in.

However, Hathaway et al. (2006) showed that line-of-sight projection effects
on a rigidly rotating fixed velocity pattern could reproduce the excess
rotation velocities found by Beck & Schou (2000). In projecting the vector
velocities onto the line-of-sight the function $\sin\phi$ (where $\phi$ is the
heliographic longitude relative to the central meridian) multiplies the
longitudinal flow velocities near the equator. Since the flows are largely
horizontal the $\sin\phi$ multiplier effectively pushes the peaks in the
Doppler pattern away from the central meridian and makes the pattern appear to
rotate more rapidly.

In fact, Schou (2003) largely removed the line-of-sight projection effects
from the equatorial Doppler data by dividing the data by a function that
approximated the function $\sin\phi$ and found that rotation velocities were
much more in line with those from global helioseismology (but noted that there
were still motions relative to this that suggested wave-like properties for
supergranules). Recently, Hathaway et al. (2010) found that the rotation
profiles as functions of latitude determined by the cross-correlation
technique used by Duvall (1980) for time lags from 2-hr to 16-hr could be
reproduced by cellular patterns that are advected by a differential rotation
with a peak velocity consistent with that found in the Sun’s surface shear
layer by global helioseismology.

Here we measure the rotation of the pattern of supergranules by analyzing the
same data used by Beck & Schou (2000) and by Schou (2003). We execute a series
of 2D Fourier transform analyses. We remove the line-of-sight projection
effects near the equator as was done by Schou (2003) and repeat the Fourier
analysis done by Beck & Schou (2000) to show that the equatorial rotation rate
of supergranules as a function of wavelength matches the equatorial rotation
rate as a function of depth determined from global helioseismology (Schou et
al., 1998). This “de-projection” can only be done at the equator and is only
approximate since the flows are not purely horizontal. We determine the
rotation rates at other latitudes by repeating the Fourier transform analysis
on the raw Doppler data (without the removal of projections effects) and using
a data simulation to support our conclusions.

##  2 DATA PREPARATION

The data consist of $1024^{2}$ pixel images of the line-of-sight velocity
determined from the Doppler shift of a spectral line due to the trace element
nickel in the solar atmosphere. The images are acquired at a 1 min cadence and
cover the full visible disk of the Sun. We average the data over 31 min with a
Gaussian weighting function which filters out variations on time scales less
than about 16 min, and sample that data at 15 min intervals. We then map these
temporally filtered images onto a $1024^{2}$ grid in heliographic latitude
from pole to pole and in longitude $\pm 90\arcdeg$ from the central meridian
(Figure 1). This mapping accounts for the position angle of the Sun’s rotation
axis relative to the imaging CCD and the tilt angle of the Sun’s rotation axis
toward or away from the spacecraft. Both of these angles include modification
in line with the most recent determinations of the orientation of the Sun’s
rotation axis (Beck & Giles, 2005; Hathaway & Rightmire, 2010). We analyze
data obtained during a 60 day period of continuous coverage in 1996 from May
24 to July 22.

Line-of-sight projection effects influence the results so we also generate and
analyze simulated data to assist in our determination of the actual rotation
as a function of latitude and depth in the Sun’s outer convection zone. We
construct the simulated data from an evolving spectrum of vector spherical
harmonics in such a manner as to reproduce the spatial, spectral, and temporal
behavior of the observed cellular flows. The amplitudes of the spectral
coefficients are constrained by matching the observed velocity spectrum
(Hathaway et al., 2000) with the radial flow component constrained by the disk
center to limb variation in the RMS Doppler signal(Hathaway et al., 2002).

Figure 1:  Heliographic map details of the line-of-sight (Doppler) velocity
from SOHO/MDI (top) and from the data simulation (bottom). Each map detail
extends $90\arcdeg$ in longitude from the central meridian on the left and
about $35\arcdeg$ in latitude from the equator (the thick horizontal line).
The mottled pattern is the Doppler signal (blue for blue shifts and red for
red shifts) due to the supergranules.

The cells are given finite lifetimes and made to rotate by adding changes to
the phases of the spectral coefficients (Hathaway et al., 2010). The rotation
rates are constrained by matching the observed rotation rates as functions of
latitude and wavelength.

The rotation rates of the cells in the simulated data are given by a fairly
simple function with the latitudinal, $\theta$, variation separated from the
wavelength, $\lambda$, variation such that

$\Omega(\theta,\lambda)/2\pi=f(\theta)[1+g(\lambda)]$ (1)

with

$f(\theta)=454-51\sin^{2}\theta-92\sin^{4}\theta\ \rm{nHz}$ (2)

and

$g(\lambda)=0.045\tanh{\lambda\over 31}\left[2.3-\tanh{(\lambda-65)\over
20}\right]/3.3$ (3)

where the wavelength, $\lambda$, is given in Mm.

##  3 EQUATORIAL ROTATION RATE

We determine the equatorial rotation rate of the supergranules by 2D Fourier
transforms with and without removing line-of-sight projection effects. The
line-of-sight projection effects can be minimized using the method described
by Schou (2003). The mapped Doppler velocities near the equator are divided by
a function, ${\rm sgn}(\phi)\sqrt{(}\sin^{2}\phi+0.01),$ which approximates
the geometric factor, $\sin\phi$, that multiplies the longitudinal velocity in
producing the Doppler signal. This “de-projected” signal and the raw Doppler
signal are both then apodized near the limb and then Fourier transformed over
longitude for the 50 latitude positions that straddle the equator. These
spectral coefficients are then Fourier transformed in time over six 10-day
intervals.

The rotation rate as a function of wavelength is determined by first finding
the temporal frequency of the centroid of the spectral power for that
wavelength. This temporal frequency is then divided by the longitudinal
wavenumber to give the synodic rotation rate which is then converted to a
sidereal rotation rate by adding a correction based on the rate of change of
the ecliptic longitude during the observations ($\sim 1\arcdeg\ {\rm
day}^{-1}$).

The results of these analyses are shown in Fig. 2 along with the rotation rate
with depth from a global helioseismology analysis by Schou et al. (1998). The
sidereal rotation rate of the de-projected supergranules as a function of
longitudinal wavelength matches the rotation rate as a function of depth
through the outer half of the solar convection zone and the simulation matches
the MDI observations. The raw Doppler data give faster rotation rates than the
de-projected data at all wavelengths (in both the MDI and the simulated data)
with larger increases for larger cells.

Figure 2:  The equatorial rotation rate for the Doppler pattern as a function
of longitudinal wavelength. Results for the raw data are shown in the upper
panel. Results for the de-projected data are shown in the lower panel. MDI
results are shown by the black dots (with $2\sigma$ error bars for wavelengths
longer than 20 Mm). Simulated data results are shown by open circles. The
equatorial rotation rate as a function of depth from global helioseismology
(Schou et al., 1998) is shown by the large red dots while the rotation profile
used in the simulation is shown by the solid black lines.

The rotation rates seen with both de-projected datasets fall slightly below
both the helioseismology results and the input profile for the simulated data.
Experiments with the simulated data suggest that this can be attributed to the
presence of a small radial flow component. This component is projected into
the line-of-sight along the equator by multiplying by a projection factor
$\cos\phi$. When the Doppler signal is de-projected by dividing by $\sin\phi$
this small signal can become large and influence the results in this manner.

We also see that the rotation rates for the raw simulated data fall
systematically below the MDI results for wavelengths $>50$ Mm. This too may be
a result of the radial flows.

These results do however lead to a key conclusion - that supergranules with
sizes from 10 Mm to 100 Mm are advected by flows within the convection zone at
depths equal to their widths.

##  4 ROTATION PROFILES

We determine the rotation rate of the supergranules as functions of latitude
and wavelength (depth) by repeating the 2D Fourier transforms on the raw
Doppler data for a series of latitude strips. Each strip is 11 pixels or $\sim
2\arcdeg$ high in latitude and offset from the previous strip by 5 pixels. We
repeat the procedure for 172 positions between $\pm 75\arcdeg$ latitude and
for each of six 10-day intervals for both MDI and simulated data. Latitudinal
rotation profiles are obtained for a series of cell wavelengths by averaging
the profiles for all waveumbers that produce wavelengths within 5 Mm of the
target wavelength. These profiles are then averaged between hemispheres and
smoothed with a 9-point binomial smoothing kernal which then limits the data
to $\pm 70\arcdeg$ latitude. The results for averages from the MDI datasets
and the from simulated data are shown in Fig. 3.

Figure 3:  The differential rotation profiles (average longitudinal velocity
relative to a frame of reference rotating at the Carrington sidereal rate of
456 nHz) for the cells with 30 Mm wavelengths are shown in the upper panel
with a solid line for the MDI data, a dashed line for the simulated data, and
$2\sigma$ error limits on the MDI data indicated by thin solid lines. The
differences between the differential rotation at the other wavelengths from
that found with the 30 Mm wavelength cells are shown in the lower panel using
colors for the different wavelengths - violet for 10 Mm, blue for 20 Mm, green
for 50 Mm, and red for 70 Mm.

The rotation velocity appears to increase at all latitudes with increasing
wavelength. However, the small increase at 70 Mm wavelength (in particular
near the equator) can be attributed to projection effects. The actual rotation
rate of the convection cells in the simulation is given by equations 1-3 and
the solid lines in Fig 2. which give a rotation rate that increases with
wavelength to a maximum at $\sim 50$ Mm at all latitudes.

The smallest cells, $10\pm 5$ Mm, are difficult to resolve at high latitudes
but clearly show slower rotation than larger cells over their observed
latitude range. Cells with wavelengths of $20\pm 5$ Mm can be resolved at all
latitudes and have significantly lower rotation rates than the larger cells.

While this method with the raw Doppler data is subject to systematic offsets
due to line-of-sight projection effects, the same offsets are present in the
simulated data and are much smaller for the smaller cells. We conclude that
the surface shear layer, in which the rotation rate increases inward, extends
to latitudes of at least $70\arcdeg$.

##  5 CONCLUSIONS

We conclude that supergranules are anchored or steered in the subsurface flows
at depths equal to their wavelengths. This is a simple explanation for the
match with the rotation rate with depth from global helioiseismology (Figure 2
lower panel) and is based on well known physics. This is consistent with
numerical simulations of convection in the outermost 16 Mm of the Sun by Stein
et al. (2011) who show that flow structures at different depths have diameters
about equal to the depth itself.

This conclusion is, however, somewhat surprising given the much smaller
estimates for the depth of typical supergranules previously determined from
the visibility of their internal flows using local helioseismology. Duvall
(1998) estimated a depth of 8 Mm for typical supergranules while Zhao &
Kosovichev (2003) estimated a depth of 15 Mm. This suggests that local
helioseismology is less sensitive to these deeper (and slower) flows and that
this new method of probing the convection zone with supergranules can probe
flows at greater depths.

We also conclude that the surface shear layer extends to a depth of $\sim 50$
Mm at all latitudes. The increase in rotation rate with depth has long been
suggested by observations and is attributed to the conservation of angular
momentum for fluid elements moving inward and outward in the near surface
layers (Foukal & Jokipii, 1975; Hathaway, 1982). Measurements of this rotation
rate increase from helioseismology (Schou et al., 1998; Corbard & Thompson,
2002) indicate that it follows this critical gradient to depths of 10-15 Mm
and reaches a maximum rotation rate at depths of 35-50 Mm. However, many
helioseismology results also suggest that the shear layer disappears at
latitudes above about $50\arcdeg$. The results reported here, using
supergranules, indicate that the shear layer extends to the highest latitudes
probed with this MDI data - $\sim 70\arcdeg$.

The rotation increase with depth given by Equation 3 follows the critical
gradient (with $\partial\ln\Omega/\partial\ln r=-1$) given by angular momentum
mixing near the surface but then drops below that value at greater depths as
the cell turn-over times become longer and the convective flows adjust to the
solar rotation.

Global helioseismology can measure the internal rotation rate to great depths
but it gives less reliable results at high latitudes. Local helioseismology
can measure the non-axisymmetric flows (as well as the axisymmetric meridional
flow) but only in the near surface layers. Using supergranules of different
sizes to probe the flows in the Sun’s convection zone extends these
measurements to greater depths and higher latitudes. This new method of
probing solar convection zone dynamics should provide information
complementary to that obtained with helioseismology.

The author would like to thank NASA for its support of this research through
grants from the Heliophysics Causes and Consequences of the Minimum of Solar
Cycle 23/24 Program and the Living With a Star Program to NASA Marshall Space
Flight Center. He is indebted to Ron Moore, Lisa (Rightmire) Upton, and an
anonymous referee whose comments greatly improved the manuscript. He would
also like to thank the American taxpayers who support scientific research in
general and this research in particular. SOHO, is a project of international
cooperation between ESA and NASA.

## References

  * Basu et al. (1999) Basu, S., Antia, H. M., & Tripathy, S. C. 1999, ApJ 512, 458
  * Beck & Giles (2005) Beck, J. G., & Giles, P. 2005, ApJ 621, L153
  * Beck & Schou (2000) Beck, J. G., & Schou, J. 2000, Sol. Phys. 193, 333
  * Corbard & Thompson (2002) Corbard, T., & Thompson, M. J. 2002, Sol. Phys. 205, 211
  * Duvall (1980) Duvall, T. L., Jr. 1980, Sol. Phys. 66, 213
  * Duvall (1998) Duvall, T. L., Jr. 1998, in Structure and Dynamics of the Interior of the Sun and Sun-like Stars ed. S. Korzennik (Noordwijk, Netherlands: ESA SP 418), 581
  * Foukal & Jokipii (1975) Foukal, P., & Jokipii, R. 1975, ApJ 199, L71
  * Hart (1954) Hart, A. B. 1954, MNRAS 114, 17
  * Hathaway (1982) Hathaway, D. H. 1982, Sol. Phys. 77, 341
  * Hathaway et al. (2000) Hathaway, D. H., Beck, J. G., Bogart, R. S., Bachmann, K. T., Khatri, G., Petitto, J. M., Han, S., & Raymond, J. 2000, Sol. Phys. 193, 299
  * Hathaway et al. (2002) Hathaway, D. H., Beck, J. G., Han, S., & Raymond, J. 2002, Sol. Phys. 205, 25
  * Hathaway & Rightmire (2010) Hathaway, D. H. & Rightmire, L. 2010, Science 327, 1350
  * Hathaway et al. (2006) Hathaway, D. H., Williams, & Cuntz, M. 2006, ApJ 644, 598
  * Hathaway et al. (2010) Hathaway, D. H., Williams, P. E., Dela Rosa, K., & Cuntz, M. 2010, ApJ 725, 1082
  * Leighton et al. (1962) Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ 135, 474
  * November et al. (1981) November, L. J., Toomre, J., Gebbie, K. B. & Simon, G. W. 1981, ApJ 245, L123
  * Scherrer et al. (1995) Scherrer, P. H., et al. 1995, Sol. Phys. 162, 129
  * Schou et al. (1998) Schou, J. et al. 1998, ApJ 505, 390
  * Schou (2003) Schou, J. 2003, ApJ 596, L259
  * Simon & Leighton (1964) Simon, G. W. & Leighton, R. B. 1964, ApJ 140, 1120
  * Stein et al. (2011) Stein, R. F., LagerFjärd, A., Nordlund, Å., & Georgobiani, D. 2011, Sol. Phys. 268, 271
  * Thompson et al. (1996) Thompson, M. J., et al. 1996, Science 272, 1300
  * Title et al. (1989) Title, A. M., Tarbell, T. D., Topka, K. P., Ferguson, S. H., Shine, R. A. & The SOUP Team 1989, ApJ 336, 475
  * Wang & Zirin (1989) Wang, H. & Zirin, H. 1989, Sol. Phys. 120, 1
  * Zhao & Kosovichev (2003) Zhao, J. & Kosovichev, A. G. 2003, in GONG+ 2002. Local and Global Helioseismology: the Present and Future ed. H. Sawaya-Lacoste (Noordwijk, Netherlands: ESA SP 517), 417