# Weighted Graph Nodes Clustering via Gumbel Softmax

Deepak Bhaskar Acharya  Computer Science5000 Technology Dr
NWHuntsvilleAlabamaUsa35805<EMAIL_ADDRESS>and  Huaming Zhang  Computer
Science5000 Technology Dr NWHuntsvilleAlabamaUsa35805<EMAIL_ADDRESS>
(2020)

###### Abstract.

Graph is a ubiquitous data structure in data science that is widely applied in
social networks, knowledge representation graphs, recommendation systems, etc.
When given a graph dataset consisting of one graph or more graphs, where the
graphs are weighted in general, the first step is often to find clusters in
the graphs. In this paper, we present some ongoing research results on graph
clustering algorithms for clustering weighted graph datasets, which we name as
Weighted Graph Node Clustering via Gumbel Softmax (WGCGS for short). We apply
WGCGS on the Karate club weighted network dataset. Our experiments demonstrate
that WGCGS can efficiently and effectively find clusters in the Karate club
weighted network dataset. Our algorithm’s effectiveness is demonstrated by (1)
comparing the clustering result obtained from our algorithm and the given
labels of the dataset; and (2) comparing various metrics between our
clustering algorithm and other state-of-the-art graph clustering algorithms.

Graph Neural Networks, Gumbel-Softmax, Weighted Graph Clustering.

††journalyear: 2020††copyright: acmcopyright††conference: 2020 ACM Southeast
Conference; April 2–4, 2020; Tampa, FL, USA††booktitle: 2020 ACM Southeast
Conference (ACMSE 2020), April 2–4, 2020, Tampa, FL, USA††price: 15.00††doi:
10.1145/3374135.††isbn: 978-1-4503-7105-6/20/03††ccs: Computing methodologies
Cluster analysis

##  1\. Introduction

Graph clustering concerns the discovery of densely related node groups in a
graph. Here an edge between two nodes normally implies any underlying
similarity or affinity between nodes, while the absence of an edge suggests
dissimilarity and distance. Given the (often noisy) similarity/dissimilarity
observations embedded in the network, graph clustering thus attempts to infer
groups of closely connected nodes. Graph analysis using the technique of
machine learning has been recognized as the strength of graphs is enormous
(Keyulu, ), i.e., graphs can be used to denote a vast variety of structures in
different fields, including (Hamilton, ) social networks, (Hamaguchi, )
knowledge graphs, (Khalil, ) natural science (physical systems
(BattagliaPascanu, ) and protein-protein interaction networks (Shariat, )) and
many other research areas. In several applications, such as community
detection, user profiling, graph clustering, is an essential subroutine. In
machine learning, statistics, and pattern recognition, clustering is a subject
of active research. Clustering, which is used almost everywhere in computing,
is a simple problem. It requires splitting the data into groups of similar
objects, and most machine learning applications are the fundamental issue
behind it. Data mining leads to datasets with many different attributes being
clustered. Therefore the relevant clustering algorithms need particular
computational requirements. Several algorithms have recently emerged that meet
these requirements and have been used successfully for real-life data mining
challenges.

Clustering is the unsupervised mechanism by which natural clusters are found
such that objects from the same cluster are identical, and objects from
separate clusters are distinct. If similarity relationships between objects
are interpreted in clustering as a simple, weighted graph where objects are
vertices and similarities between objects are edge weights, clustering reduces
the graph clustering problem. However, in many applications, the edge/non-edge
between each pair of nodes can reflect very different confidence levels of
node similarity. In certain circumstances, edge observation (the absence of an
edge) can be created by accurate measurements. Therefore, it is a clear
indication that the same clusters should be assigned to the two nodes.In other
cases, the findings will be very uncertain, and thus, no knowledge about the
cluster configuration will give an edge or the lack of it. The findings
between certain node pairs can hold no information as an exceptional example,
so they are practically unnoticed.

By extending Acharya and Zhang (acharya2020community, ) Community Detection
Clustering technique via Gumbel Softmax (WGCGS), in this research paper, we
propose the new method for Weighted Graph Nodes Clustering via Gumbel Softmax.
The article presented by using (acharya2020community, ) focuses on the
unweighted graph datasets. In this research, we recommend finding the cluster
of the nodes when the graph is weighted.

The remaining of the paper is organized as follows. In Section 2, we introduce
the background and related works in more detail. In Section 3, we introduce
our main method. In Section 4, we present the experimental results. In section
5, we conclude our work with the future enhancements.

##  2\. Background and Related work

###  2.1. Community Detection Clustering via Gumbel Softmax

Different algorithms have been proposed, and there is extensive research on
community detection. Detecting communities is a method for defining related
groups that can be difficult, depending on the graph network’s size and scale.
Acharya and Zhang (acharya2020community, ) provide a general understanding of
their paper on working with unsupervised graph datasets and identifying graph
nodes’ clusters.

In general, let the adjacency matrix be $A_{n\times n}$, where ’n’ represents
the total number of nodes in the graphs dataset. The $W_{C}$ matrix size is
$n\times k$, where k indicates the number of clusters. After obtaining a
matrix of the size $k\times k$ then perform the $W_{C}^{t}AW_{C}$ operation
and call this matrix as $R_{k\times k}$.

The resulting $R_{k\times k}$ matrix displays the cluster strength, where the
primary diagonal shows the strength of the data within the cluster groups, and
the other $R$ matrix components show the strength of the data between the
different clusters. Then as a discrete distribution of probabilities, apply
the softmax function to the $R$ matrix obtained to express inputs.
Mathematically be defined as follows:

(1) $Softmax(x_{i})=\frac{exp(x_{i})}{\sum_{j=0}^{m}{exp(x_{j})}}\ for\
i=1,2,3,......,m.$

###  2.2. Gumbel Softmax Approach on Graph Features Selection and Extraction

In the suggested approach for the citation datasets, the article features
selection and extraction for Graph Neural Networks, Acharya and Zhang
(bhaskar2019feature, ) selected and extracted Graph Neural Network (GNN). By
applying the feature selection and extraction technique to GNNs using gumbel
softmax, they conduct tests using different reference datasets: Cora, Pubmed,
and Citeseer. Acharya and Zhang demo an example of the Citeseer dataset, where
375 of 3703 features are selected and ranked according to prominent features.
The proposed deep learning method works well with reduced features, decreasing
the number of features that the dataset initially had by about 80-85 percent.
The experiment results show that the accuracy declines steadily when using the
selected classification features falling within the ranges 1-75, 76-150,
151-225, 226-300, and 301-375.

In general, let the input feature matrix be $X_{n\times f}$ where ’n’ and ’f’
represent the total number of nodes and features in the graph’s dataset. The
resulting matrix is $M_{n\times k}$, where ’k’ represents the features we
select from the ’f’ features.

Acharya and Zhang select the features in a graph citation dataset using the
gumbel softmax function and the proposed approach. Gumbel softmax’s
distribution is (Shixiang, ; maddison2016concrete, )’ a continuous
distribution over the simplex that can approximate categorical distribution
samples.’ A categorical distribution is a one-hot vector that defines the
maximum probability of one and all the other probability as zero.

The two layer Graph Convolution Network (GCN) used in the experiment is
defined as

(2) $GCN(X,A)=Softmax(A(ReLu(AXW_{G}W_{1}))W_{2})$

To verify the selected features and calculate the accuracy for classification
they use the following two layer Graph Convolution Network as defined below

(3) $GCN(X,A)=Softmax(A(ReLu(AXW_{G}^{\prime}))W_{2})$

$A$: Adjacency matrix of the undirected graph G.
$X$: Input feature matrix.
$W_{G}$: gumbel softmax feature selection / feature extraction matrix.
$W_{G}^{\prime}$: feature selection / feature extraction matrix obtain ed from
the result of Equation 2.
$W1,W2$: Layer-specific trainable weight matrix.
$ReLu$ : Activation function ReLu(.) = max(0,.).

##  3\. Proposed Method

The proposed approach in the article Community Detection Clustering via Gumbel
Softmax (CDCGS) clusters of graph dataset nodes and the article Feature
Selection and Extraction for Graph Neural Networks (FSEGNN) explains how
prominent features for graph citation datasets can be selected and extracted.
The dataset used is a graph in each of the papers. There are nodes and edges
in the graphs, and they have relations between the nodes that can have some
knowledge about the similarities or dependencies between the nodes. The two
papers above do not provide insight into the network of weighted graphs. We
suggest a new method of Weighted Graph Nodes Clustering via Gumbel Softmax
(WGCCS) to cluster the graph nodes as an extension to the CDCGS and FSEGNN.

We consider the network of ’n’ nodes and the adjacency matrix ’A’ in our WGCCS
method. An adjacency matrix is used to represent a finite graph. The elements
of the matrix show whether the pairs of vertices in the graph are adjacent or
not. Then we cluster the graph into k clusters, applying our method.

Gumbel softmax distributions can be used to estimate the sampling process of
discrete data if we have a stochastic neural network with discrete variables.
Using backpropagation, the network can then be trained, where the efficiency
of the network will depend on the temperature parameter choice.

The method used in the experiment is defined as below:

(4) $Weighted-Graph-Cluster(Adj)=Softmax(W_{C}^{t}(Adj)W_{C})$

In the Equation 4, ’$Adj$’ indicates the Adjacency matrix of the undirected
weighted graph G, where weight of graph edges indicates the similarity between
the nodes of the graph and ’$W_{C}$’ indicates the gumbel cluster weight
matrix.

We have the gumbel cluster weight matrix $W_{C}$ of size $N\times k$ once
model training is completed, where $N$ is the number of graph nodes, and $K$
is the number of clusters. Then, we look at the sum of the row elements, where
each row sums up to 1, and the row number’s maximum index is the cluster to
which the graph node belongs. For example, let’s consider k= 2, i.e., we’re
trying to cluster the dataset into two cluster groups. Let us say that the
data is [0.65 0.35] for the first row, where 0.65 is $0^{th}$, and 0.35 is
$1^{st}$. By looking at the maximum value in the row, we get the corresponding
maximum value index, and one can easily say that the data belongs to cluster
0. If [0.39 0.61] is the data for the second row, 0.39 will be $0^{th}$, and
0.61 will be $1^{st}$ index. By looking at the maximum value in the row, we
get the corresponding maximum value index, and one can easily say that the
data belongs to cluster 1.

The resulting ${k\times k}$ size matrix derived from the 4 equation is
compared to the identity matrix ($I_{k\times k}$) to find network losses. Here
the cluster groups represent diagonal components. We may then form the cluster
groups of the data points by using our proposed approach. A detailed analysis
is given in the result section, taking into account an example of a dataset
weighted karate club network.

##  4\. Experiment Results

###  4.1. Datasets

Zachary’s Weighted Karate Club Network is a well-known dataset describing the
relationships in a university karate club used by Wayne W. Zachary in his
paper "An Information Flow Model for Conflict and Fission in Small Groups."
This dataset is known for its simple description of the community structure,
which occurs because it is possible to cluster network nodes into strongly
linked sets. Focused on Mr. Hi, the karate instructor, and John A, the club
president, Zachary’s Karate Club network can be split into two groups. The
network accurately predicts how, following a disagreement over the pay, the
karate club splits into two new clubs and creates a rift within Hr. Hi and
John A. The network shows which members of the club will join the new club by
analyzing community members’ meetings outside the club’s background in 33 of
34 instances.

The initial dataset of Zachary’s Karate Club is weighted by multiple
friendship measures. Several visualizations of Zachary’s Karate Club have been
created since the first article was published in 1977. Michelle Girvan and
Mark Newman (girvan2002community, ) used Zachary’s Karate Club again in 2002
to demonstrate community structure in their article "Community structure in
social and biological networks." Figure 1 is a network created by a computer
that contains each member of the club.

Figure 1. Weighted Karate Club Network

###  4.2. Results

The weighted network of Zachary karate club is the network we can think of to
evaluate our methodology. We also know that the two teacher nodes are dominant
graph nodes 0 and 33, as the weighted Zachary network is labeled. Any other
network nodes can either join the teacher community of node 0 or the teacher
community of node 33. Until clustering, the original network is seen in figure
2. Second, we run our strategy to cluster the network into 2-clusters like the
Zachary karate club’s ground-truth community structure has 2-clusters. We
compare our proposed method results with the results of algorithms such as
greedy optimization of leading eigenvector of the community matrix(ECM)
(newman2006finding, ), edge betweenness(EB) (PhysRev, ), short random
walks(RW) (pons2005computing, ), Infomap community finding(ICF) (Rosvall_2009,
), multi-level optimization of modularity(MOM) (Blondel_2008, ), propagating
labels(PL) (Raghavan_2007, ). Metrics such as adjusted randomized index (ARI)
(sklearn_api, ), homogeneity (HOMO), normalized mutual information (NMI),
completeness (COMP), and v-measure V-MES) (rosenberg2007v, ) has also been
computed to complete the analysis. ARI = 1, HOMO = 1, NMI = 1, COMP = 1, V-MES
= 1 and the modularity calculation = 0.371 are the corresponding metric
metrics if we cluster the weighted Zachary network into 2-clusters (WGCGS(2C))
using WGCGS. The value of ARI, HOMO, NMI, COMP, and V-MES indicates that the
proposed algorithm clusters the nodes to their respective groups as per the
target label. Figure 2 represents Zachary’s network of 2-cluster karate clubs
resulting from our approach proposed.

Whenever we consider the graph dataset, the higher the modularity value
indicates a stronger assignment to the cluster. To obtain higher modularity
for the stronger assignment, we clustered the weighted Zachary karate club
network into four clusters. 4-clustered Zachary network (WGCGS(4C)) gives us a
0.4197 modularity metric (Figure 3). Considering that students joining node 0
as cluster group 1 and node 33 as cluster group 2, node 0 and node 33 are
already connected to several nodes, and nodes 23, 24, 25, 27, 28, 31 (cluster
3) and 4, 5, 6, 10, 16 (cluster 4) mean that students are not joining cluster
group 1 or 2. To test the classification of metadata groups by various
algorithms, we used the best scores of the Zachary karate club network along
with different algorithms. Our WGCGS approach outperforms the metric results
of the ECM, EB, RW, MOM, and PL algorithms, and we presented them in Table
4.2.

The result obtained from the 4 equation is compared to the loss function. The
loss function used for our experiment is the identity matrix ($I_{k\times
k}$), where the cluster group is represented by each diagonal member.

Figure 2. Weighted Karate Club Network with 2-Clusters

Figure 3. Weighted Karate Club Network with 4-Clusters Table 1. Zachary karate
club network metrics for the different algorithms and our method (WGCGS)

o XXXXXXXXX Metric | ECM | EB | RW | ICF | MOM | PL | WGCGS(2C) | WGCGS(4C)
---|---|---|---|---|---|---|---|---
ARI | 0.5120 | 0.5125 | 0.2620 | 0.7021 | 0.4619 | 0.4714 | 0.5414 | 1
NMI | 0.6770 | 0.6097 | 0.3729 | 0.6994 | 0.5866 | 0.5282 | 0.6873 | 1
HOMO | 1 | 0.8850 | 0.5588 | 0.8535 | 0.8535 | 0.6905 | 1 | 1
COMP | 0.5118 | 0.4651 | 0.2798 | 0.5925 | 0.44.68 | 0.4277 | 0.5235 | 1
V-MES | 0.6770 | 0.6097 | 0.3729 | 0.6994 | 0.5866 | 0.5282 | 0.6872 | 1

##  5\. Conclusion and Enhancements

We introduced the Weighted Graph Nodes Clustering via Gumbel Softmax strategy
for the weighted karate club network. The experimental results illustrate the
efficiency by selecting appropriate parameter values and checking the
resulting clustering accuracy. The approach currently available is just as
diverse as the implementations for weighted undirected graph clustering. We
are currently experimenting with applying the principle of clustering for
Graphs on few more weighted graph datasets and a directed graph dataset.

## References

  * [1] D. B. Acharya and H. Zhang.  Community detection clustering via gumbel softmax.  SN Computer Science, 1(5):1–11, 2020.
  * [2] D. B. Acharya and H. Zhang.  Feature selection and extraction for graph neural networks.  In Proceedings of the 2020 ACM Southeast Conference, ACM SE ’20, page 252–255, New York, NY, USA, 2020. Association for Computing Machinery.
  * [3] P. Battaglia, R. Pascanu, M. Lai, D. J. Rezende, et al.  Interaction Networks for Learning about Objects, Relations and Physics.  In Advances in neural information processing systems, pages 4502–4510, Barcelona, Spain, 2016.
  * [4] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre.  Fast unfolding of communities in large networks.  Journal of Statistical Mechanics: Theory and Experiment, 2008(10):P10008, Oct 2008.
  * [5] L. Buitinck, G. Louppe, M. Blondel, F. Pedregosa, A. Mueller, O. Grisel, V. Niculae, P. Prettenhofer, A. Gramfort, J. Grobler, R. Layton, J. VanderPlas, A. Joly, B. Holt, and G. Varoquaux.  API design for machine learning software: experiences from the scikit-learn project.  In ECML PKDD Workshop: Languages for Data Mining and Machine Learning, pages 108–122, 2013.
  * [6] A. Fout, J. Byrd, B. Shariat, and A. Ben-Hur.  Protein Interface Prediction using Graph Convolutional Networks.  In Advances in Neural Information Processing Systems, pages 6530–6539, Long Beach, USA, 2017.
  * [7] M. Girvan and M. E. Newman.  Community structure in social and biological networks.  Proceedings of the national academy of sciences, 99(12):7821–7826, 2002.
  * [8] T. Hamaguchi, H. Oiwa, M. Shimbo, and Y. Matsumoto.  Knowledge Transfer for Out-of-knowledge-base Entities: A Graph Neural Network Approach.  In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, Melbourne, Australia, 2017.
  * [9] W. Hamilton, Z. Ying, and J. Leskovec.  Inductive Representation Learning on Large Graphs.  In Advances in Neural Information Processing Systems, pages 1024–1034, Long Beach, USA, 2017.
  * [10] E. Jang, S. Gu, and B. Poole.  Categorical Reparameterization with Gumbel-softmax.  In ICLR, Toulon, France, 2017.
  * [11] E. Khalil, H. Dai, Y. Zhang, B. Dilkina, and L. Song.  Learning Combinatorial Optimization Algorithms over Graphs.  In Advances in Neural Information Processing Systems, pages 6348–6358, Long Beach, USA, 2017.
  * [12] C. J. Maddison, A. Mnih, and Y. W. Teh.  The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables.  In ICLR, Toulon, France, 2017.
  * [13] M. E. Newman.  Finding community structure in networks using the eigenvectors of matrices.  Physical review E, 74(3):036104, 2006.
  * [14] M. E. J. Newman and M. Girvan.  Finding and evaluating community structure in networks.  Phys. Rev. E, 69:026113, Feb 2004.
  * [15] P. Pons and M. Latapy.  Computing communities in large networks using random walks (long version), 2005.
  * [16] U. N. Raghavan, R. Albert, and S. Kumara.  Near linear time algorithm to detect community structures in large-scale networks.  Physical Review E, 76(3), Sep 2007.
  * [17] A. Rosenberg and J. Hirschberg.  V-measure: A conditional entropy-based external cluster evaluation measure.  pages 410–420, 2007.
  * [18] M. Rosvall, D. Axelsson, and C. T. Bergstrom.  The map equation.  The European Physical Journal Special Topics, 178(1):13–23, Nov 2009.
  * [19] K. Xu, W. Hu, J. Leskovec, and S. Jegelka.  How Powerful Are Graph Neural Networks?  In ICLR, New Orleans, USA, 2019.


# The Sun’s Shallow Meridional Circulation

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

The Sun’s global meridional circulation is evident as a slow poleward flow at
its surface. This flow is observed to carry magnetic elements poleward -
producing the Sun’s polar magnetic fields as a key part of the 11-year sunspot
cycle. Current theories for the sunspot cycle assume that this surface flow is
part of a circulation which sinks inward at the poles and turns equatorward at
depths below 100 Mm. Here we use the advection of the Sun’s convection cells
by the meridional flow to map the flow velocity in latitude and depth. Our
measurements show the largest cells clearly moving equatorward at depths below
35 Mm - the base of the Sun’s surface shear layer. This surprisingly shallow
return flow indicates the need for substantial revisions to solar/stellar
dynamo theory.

Sun: dynamo, Sun: rotation, Sun: surface magnetism

##  1 INTRODUCTION

The Sun’s meridional circulation has been a part of solar magnetic dynamo
theory for half a century. A poleward meridional flow from the Sun’s mid
latitudes was invoked in the earliest models (long before the flow was
actually measured) to transport magnetic elements from decaying sunspot
regions to the poles where they would erode the opposite polarity magnetic
field from the old sunspot cycle and build up the polar fields of the new
sunspot cycle (Babcock, 1961). This surface meridional flow, along with its
latitudinal structure and variation in time, is now well observed (Topka et
al., 1982; Komm et al., 1993; Hathaway et al., 1996; Hathaway & Rightmire,
2010, 2011) and its role in the surface magnetic flux transport is well
established (DeVore & Sheeley, 1987; van Ballegooijen et al., 1998; Schrijver
& Title, 2001). Dynamo theories over the last decade and a half have assumed
that the mass traveling poleward in the surface layers sinks inward at the
poles and returns to the equator along the base of the Sun’s convection zone
at a depth of 200 Mm. In these theories this slow, dense, equatorward flow is
responsible for the equatorward drift of sunspot activity (Dikpati &
Choudhuri, 1994; Dikpati & Charbonneau, 1999; Nandy & Choudhuri, 2001). Dynamo
models based on this deep meridional circulation have recently been used to
predict the currently emerging sunspot cycle, albeit with disparate results
(Dikpati et al., 2006; Choudhuri et al., 2007).

The meridional flow is difficult to measure. Its amplitude ($\sim 10-20\rm{\
m\ s}^{-1}$) is more than an order of magnitude weaker than that of the other
major flows observed on the surface of the Sun. The axisymmetric longitudinal
flow, differential rotation, has a dynamic range of $\sim 200\rm{\ m\ s}^{-1}$
and the non-axisymmetric cellular convection flows have typical velocities of
several hundred $\rm{\ m\ s}^{-1}$. The meridional flow is observed to be
poleward from the equator with peak flow speeds in the mid latitudes. The flow
amplitude measured from Doppler shifts of spectral lines formed at the surface
is $\sim 20\rm{\ m\ s}^{-1}$ (Hathaway et al., 1996; Ulrich, 2010) while the
amplitude found by measuring the motions of the small magnetic features is
$\sim 12\rm{\ m\ s}^{-1}$ (Komm et al., 1993; Hathaway & Rightmire, 2010,
2011). The meridional flow can also be measured using the methods of local
helioseismology which yield a peak velocity of $\sim 20\rm{\ m\ s}^{-1}$ that
appears to be constant with depth down to $\sim 26$ Mm (Giles et al., 1997;
Schou & Bogart, 1998). Recently, however, two new measurement methods have
indicated a decrease in amplitude with depth. A method using global
helioseismology (Mitra-Kraev & Thompson, 2007) found a meridional flow that
decreased with depth and became equatorward at a depth of only 40 Mm - but
with a large range of error. A method using the movement of the larger solar
convection cells, supergranules, also found a meridional flow that decreased
with depth but without precise depth information and without detection of a
return flow (Hathaway et al., 2010). While some helioseismic studies indicate
a poleward meridional flow at depths well below 26 Mm, Duvall & Hanasoge
(2009) found that those methods are prone to systematic errors and Gough &
Hindman (2010) conclude that the flow below 30 Mm remains unknown. Furthermore
Beckers (2007) has suggested that projection effects may have also compromised
some of the local helioseismology results and concludes that the meridional
flow velocity may decrease with depth.

Here we measure the meridional flow by tracking the motions of supergranules,
but extend the analysis to include larger cells with deeper roots. Numerical
models of compressible convection with radiative transfer in the near surface
layers (Stein & Nordlund, 2000) clearly show that small cells dominate at the
surface and larger structures are found at increasing depth. Supergranules
cover the surface of the Sun and have a broad range of sizes that sample a
corresponding range of depths. We measure the motion of the pattern of
supergranules by analyzing data acquired by the Michelson Doppler Imager (MDI)
(Scherrer et al., 1995) on the ESA/NASA Solar and Heliospheric Observatory
(SOHO) satellite in 1996 and 1997.

##  2 DATA PREPARATION

The data consist of 1024x1024 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. We
average them over 31 min with a Gaussian weighting function which filters out
any velocity components that vary on time scales less than about 16 min. We
then map these temporally filtered images onto a 1024x1024 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 for the most recent determinations of the
orientation of the Sun’s rotation axis (Beck & Giles, 2005; Hathaway &
Rightmire, 2010). We analyze data obtained during two 60+ day periods of
continuous coverage - one in 1996 from May 24 to July 24 and another in 1997
from April 14 to June 18.

We also generate and analyze simulated data to assist in our determination of
the representative depths. 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
(Hathaway et al., 2010). The cells are advected in longitude by differential
rotation and in latitude by meridional flow, both of which vary with cell
size.

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 supergranule convection cells. The latitudinal movement
of these supergranules yields a measurement of the Sun’s meridional
circulation.

##  3 DATA ANALYSIS AND RESULTS

We determine the motions of the cellular patterns in longitude and latitude
for strips of data by finding the displacement of the maximum in the cross-
correlation with similar strips from images acquired at time lags of 2, 4, 8,
16, 24, and 32 hr. Each strip is 11 pixels or $\sim 2\arcdeg$ high in latitude
and 600 pixels or $\sim 105\arcdeg$ long in longitude. We repeat this
procedure for 796 latitude positions between $\pm 70\arcdeg$ latitude and for
each hour over the 60 days of each MDI dataset and 60 days of simulated data.
This cross-correlation method was first used to determine the equatorial
rotation rate by Duvall (1980) who concluded that larger cells live longer and
rotated faster. Beck & Schou (2000) used a 2D Fourier transform method and
found a rotation rate that increased with the wavelength of the features.

We calculate the average differential rotation and meridional flow profiles
and fit them with 4th order polynomials in $\sin\theta$ where $\theta$ is the
heliographic latitude (Figure 2). The rotational velocity increases with
increasing time lags to a maximum at 24 hr but then decreases at 32 hr. The
meridional flow velocity decreases with time lag and, at the 32 hr time lag,
reverses sign. The cellular flows that live long enough to be positively
correlated 32 hrs later are moving equatorward. This is a clear detection of
the meridional return flow. The individual data points have a standard error
of $\sim 1\rm{\ m\ s}^{-1}$ but the vast majority of points indicate an
equatorward flow significantly bigger than this. The curves fit through the
data points indicate an equatorward return flow of $1.8\rm{\ m\ s}^{-1}$ with
a standard error of $<0.2\rm{\ m\ s}^{-1}$.

Figure 2:  Flow profiles as functions of latitude from the 1996 SOHO/MDI data
(top row), the 1997 SOHO/MDI data (middle row), and the simulated data (bottom
row). The flow velocities measured at each latitude are shown with colored
dots for each time lag as indicated in the figure. The solid lines with the
same color coding represent the 4th order polynomial fits to each profile. The
meridional flow decreases in amplitude with increasing time lag and reverses
direction for 32 hr lags. The rotation rate increases as the time lag
increases up to 24 hr then drops at 32 hr.

We determine the characteristic convection cell wavelengths for the different
time lags using the wavelength dependence of the differential rotation and
meridional flow profiles used in the simulation. The differential rotation
measurements are largely reproduced with a relatively simple rotation
velocity, $u$, relative to an inertial (sidereal) frame of reference given by

$\displaystyle u(\theta,\lambda)$ $\displaystyle=$
$\displaystyle[(1980-246\sin^{2}\theta-365\sin^{4}\theta)\cos\theta]$ (1)
$\displaystyle[1.0+0.046\tanh(\lambda/35)]\rm{\ m\ s}^{-1}$

where the wavelength, $\lambda$, is given in Mm. (Note that the prograde
velocities plotted in Figure 2 are relative to a frame of reference rotating
at the Carrington rotation rate with $u_{C}(\theta)=1991\cos\theta$.) The
meridional flow measurements are largely reproduced with a northward velocity,
$v$, given by

$\displaystyle v(\theta,\lambda)$ $\displaystyle=$
$\displaystyle[(65\sin\theta-78\sin^{3}\theta)\cos\theta]$ (2)
$\displaystyle[\tanh((35-\lambda)/20)]\rm{\ m\ s}^{-1}$

where the reproduction of the return flow is particularly sensitive to the
zero crossing occurring at $\sim 35$ Mm.

Figure 3 shows the equatorial differential rotation velocity relative to the
surface and the meridional flow velocity at $30\arcdeg$ latitude along with
the data points from the MDI (averaged for the two years) and simulation data
analyses. The data points from the simulation virtually coincide with those
from the MDI data except at the 32 hr time lag and for the northward velocity
at the 16 hr time lag. No doubt, better fits could be obtained with more
complicated flow profiles. It is apparent, however, that the differential
rotation velocity must decrease for wavelengths greater than $\sim 35$ Mm and
that the reversal in the meridional flow direction must be more abrupt.

Figure 3 also shows that the points do not fall right on the curves for the
input flow profiles. We attribute this to two different processes. The two
latitudes chosen for this figure represent the latitudes at which each flow
reaches its maximum. Since the convection cells span a finite range of
latitudes the measured values should be less than these maxima. However, the
differential rotation signal is subject to an additional line-of-sight
projection effect (Hathaway et al., 2006) which makes the Doppler features
appear to rotate faster. This effect raises the measured values to
increasingly higher values with increasing wavelength for the prograde
velocity.

Figure 3:  The simulation meridional flow speed at $30\arcdeg$ latitude (top)
and equatorial differential rotation relative to the surface (bottom) as
functions of wavelength are shown by the solid lines. The observed values from
the cross-correlation analysis with the MDI data and the simulated data
(filled circles and crosses respectively) are shown at their characteristic
wavelengths. Error bars centered on each symbol represent $2\sigma$ errors.
The dashed line shows the theoretical limit to the rate of increase in
rotation rate at the equator and the open circles show the rotation velocity
determined from global helioseismology by Schou et al. (1998) both assuming
that the wavelength equals the anchoring depth of the cells.

##  4 CONCLUSIONS

The relationship between the wavelength of a convection cell and the depth at
which it is anchored or steered is well constrained by the stability of the
surface shear layer and observations of the rotation rate with depth from
global helioiseismology. An 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; Gilman & Foukal, 1979; Hathaway, 1982).
However, a rotation rate which increases inward faster than that given by the
conservation of angular momentum is dynamically unstable (Chandrasekhar,
1961). 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. This gradient (the dashed line in Figure 3)
and the helioseismicly determined rotation rates (open circles in Figure 3)
are matched by the simulation input rotation profile if we assume that the
convection cells are anchored at depths equal to their widths. Shallower cells
would give unstable gradients. The mass density increases nearly quadratically
with depth so it is reasonable to expect the cells to be advected by the flows
near their deepest extent.

This association between cell wavelength and depth indicates that the poleward
meridional flow seen at the surface reverses at a depth of 35 Mm - the base of
the surface shear layer where the rotation rate reaches its maximum. Although
this shallow return flow violates the assumptions of flux transport dynamos
(Dikpati & Choudhuri, 1994; Dikpati & Charbonneau, 1999; Nandy & Choudhuri,
2001; Dikpati et al., 2006; Choudhuri et al., 2007), it was predicted by
numerical simulations of the effects of rotation on supergranules (Hathaway,
1982), it is in agreement with global helioseismology, and it helps to
reconcile other observations.

The surface has the slowest rotation and the fastest meridional flow. Small
magnetic elements rotate faster than the surface and have poleward meridional
flow which is slower than the surface (Komm et al., 1993; Hathaway &
Rightmire, 2010, 2011). Both velocity components for the small magnetic
elements are matched at a depth of about 15 Mm. While supergranules do have a
broad range of sizes, their spectrum exhibits an excess of power at
wavelengths of 30-35 Mm (Hathaway et al., 2000). Our results indicate that
cells this size have depths roughly equal to the depth of the surface shear
layer. This hardly seems coincidental but rather suggests an intimate
connection between the characteristic size of supergranules and the depth of
the surface shear layer.

A meridional circulation confined to the surface shear layer would also
explain why numerical simulations of the solar convection zone below this
surface shear layer have had difficulty producing the observed flows (Miesch
et al., 2000). In particular, the meridional circulations in these simulations
are highly structured in latitude and highly variable in time. The source of
this structure and variability can be attributed to the small number ($\sim
100$) of convection cells that populate the simulated volume. A meridional
circulation driven by ($\sim 10000$) supergranules - the convection cells that
populate the surface shear layer - is far more likely to be less structured
and variable.

This detection of a shallow equatorward return flow for the Sun’s meridional
circulation indicates the need for a reassessment of solar dynamo theory. The
flux transport dynamo models, all of which assume and require a deep
meridional flow, apparently cannot be correct. While other dynamo models
exist, the majority of these place the dynamo action at the base of the
convection zone with a possible secondary and less organized dynamo in the
surface shear layer. Upon reassessment we may find that the surface shear
layer plays a far more important role in the global dynamo.

The author would like to thank NASA for its support of this research through a
grant from the Heliophysics Causes and Consequences of the Minimum of Solar
Cycle 23/24 Program to NASA Marshall Space Flight Center. He is also indebted
to Ron Moore and Lisa Rightmire who read and commented on the manuscript.
SOHO, is a project of international cooperation between ESA and NASA.

## References

  * Babcock (1961) Babcock, H. W. 1961, ApJ 133, 572
  * 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
  * Beckers (2007) Beckers, J. M. 2007, Sol. Phys. 240, 3
  * Chandrasekhar (1961) Chandrasekhar 1961, Hydrodynamic and Hydromagnetic Stability Oxford University Press, London
  * Choudhuri et al. (2007) Choudhuri, A. R., Chatterjee, P., & Jiang, J. 2007, Phys. Rev. Lett. 98, 131103-1
  * Corbard & Thompson (2002) Corbard, T., & Thompson, M. J. 2002, Sol. Phys. 205, 211
  * DeVore & Sheeley (1987) DeVore, C. R., & Sheeley, N. R., Jr. 1987, Sol. Phys. 108, 47
  * Dikpati & Choudhuri (1994) Dikpati, M., & Choudhuri, A. R. 1994, A&A 291, 975
  * Dikpati & Charbonneau (1999) Dikpati, M., & Charbonneau, P. 1999, ApJ 518, 508
  * Dikpati et al. (2006) Dikpati, M., de Toma, G., & Gilman, P. A. 2006, Geophys. Res. Lett. 33, L05102
  * Duvall (1980) Duvall, T. L., Jr. 1980, Sol. Phys. 66, 213
  * Duvall & Hanasoge (2009) Duvall, T. L., Jr. & Hanasoge, S. M. 2009, arXiv:0905.3132
  * Foukal & Jokipii (1975) Foukal, P., & Jokipii, R. 1975, ApJ 199, L71
  * Giles et al. (1997) Giles, P. M., Duvall, T. L., Jr., Scherrer, P. H., & Bogart, R. S. 1997, Nature 390, 52
  * Gilman & Foukal (1979) Gilman, P. A., & Foukal, P. 1979, ApJ 229, 1179
  * Gough & Hindman (2010) Gough, D. O., & Hindman, B. W. 2010, ApJ 714, 960
  * Hathaway (1982) Hathaway, D. H. 1982, Sol. Phys. 77, 341
  * Hathaway et al. (1996) Hathaway, D. H. et al. 1996, Science 272, 1306
  * 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 & Rightmire (2010) Hathaway, D. H. & Rightmire, L. 2010, Science 327, 1350
  * Hathaway & Rightmire (2011) Hathaway, D. H. & Rightmire, L. 2011, ApJ 729, 80
  * 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
  * Komm et al. (1993) Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys. 147, 207
  * Miesch et al. (2000) Miesch, M. S., Elliott, J. R., Toomre, J., Clune, T. L., Glatzmaier, G. A., Gilman, P. A. 2000, ApJ 532, 593
  * Mitra-Kraev & Thompson (2007) Mitra-Kraev, U. & Thompson, M. J. 2007, Astron. Nachr. 328, 1009
  * Nandy & Choudhuri (2001) Nandy, D. & Choudhuri, A. R. 2012, Science 296, 1671
  * Scherrer et al. (1995) Scherrer, P. H., et al. 1995, Sol. Phys. 162, 129
  * Schou & Bogart (1998) Schou, J. & Bogart, R. S. 1998, ApJ 504, L131
  * Schou et al. (1998) Schou, J. et al. 1998, ApJ 505, 390
  * Schrijver & Title (2001) Schrijver, C. J. & Title A. M. 2001, ApJ 551, 1099
  * Stein & Nordlund (2000) Stein, R. F. & Nordlund, Å., Sol. Phys. 192, 91
  * Topka et al. (1982) Topka, K., Moore, R., LaBonte, B. J., & Howard, R. 1982, Sol. Phys. 79, 231
  * Ulrich (2010) Ulrich, R. K. 2010, ApJ 725, 658
  * van Ballegooijen et al. (1998) van Ballegooijen, A. A., Cartledge, N. P., & Priest, E. R. 1998, ApJ 501, 866