\section{Introduction} \label{sec:intro}

Jets are ubiquitous in astrophysical accreting systems. Two types of jets have been observed, namely continuous jets and episodic ones. In addition to the temporal difference,  the radiation spectrum, polarization, and velocities of these two types of jet are also different, as summarized in \citet{Fender04}. Episodic jets and radiation flares are often found to be physically associated with each other. One such example is Sgr A*, the supermassive black hole in the Galactic center. Flares are observed in multi-wavebands, ranging from radio, submm, infrared, and X-ray \citep{Baganoff01,Genzel2003,Eckart06,Trippe2007,Yusef-Zadeh2006,Marrone2008}. The flares are strongest in IR and X-ray \citep{Baganoff01,Genzel2003} and they occur simultaneously \citep{Genzel2003}, followed by flares at submm and radio wavebands \citep{Eckart06,Marrone2008}. The peak of the light curve at 43 GHz is found to lead that at 22 GHz by 20-40 minutes \citep{Yusef-Zadeh2006}. This time lag is regarded as the evidence for the ejection and adiabatic expansion of plasmoid from the accretion flow \citep{Yusef-Zadeh2006,Marrone2008}. This interpretation is confirmed by the 7 mm VLBA observations which found possible image of the ejected blob 4.5 hour later after the IR flare, and the speed of the blob is estimated to be 0.4 speed of light \citep{Rauch2016}.

The phenomena of multi-waveband flares with frequency-dependent time lag and associated ejection are common in black hole sources. For example, large radio flares were detected in three AGNs, namely M81  \citep{King2016}, blazar PKS 1510–089 \citep{Park2019}, and M87 \citep{Hada2014}, preceded by an X-ray flare, $\gamma$-ray flare, and $\gamma$-ray flare in these three sources, respectively. Discrete blobs associated with flares were clearly detected by VLBA or VLBI, with the speed of the blob measured to be 0.5 and 0.2 speed of light in the case of M81 \citep{King2016} and M87 \citep{Hada2014}. Similar phenomenon occurs in a black hole X-ray binary, Cygnus X-1 \citep{Wilms2007}.

Flares and episodic ejection (as jet knots) have also been observed in protostellar accretion systems \citep{Lee2020,Wolk2005,Feigelson1999}. Although consensus has not been reached on the mechanisms of the ejection and flares, both of them are believed to be powered by magnetic activities in the accretion flow \citep{Feigelson1999}, so it is natural to conjecture that they are physically associated with each other. More interestingly, even in the case of the Sun, large solar flares and coronal mass ejections (CMEs) also occur together \citep{2000JGR...105.2375L,Zhang2005}.

Many works have been published to explain the origin of black hole flares, usually by invoking synchrotron radiation or Compton scattering by relativistic electrons. These electrons are either accelerated by shocks in jets \citep{Markoff2001} or by magnetic reconnection in accretion flows \citep{Yuan2004,Yuan09,Dodds2009,Dodds2010,Nathanail20,Nathanail2021,Ripperda20,Ripperda2021,Dexter20,Dexter2021,2020MNRAS.499.1561Z}.
These works usually focus only on the flares, with little discussions on episodic ejections and their associations with flares.

An exception is \citet{Yuan09}. In this work, by analogy with the coronal mass ejections in the Sun, they propose that the turbulent motion and differential rotation of the black hole accretion flow will deform the magnetic field lines; consequently, magnetic reconnection occurs and flux ropes are formed in the coronal region of the accretion flow, which will then be ejected by the magnetic compression force.

These processes must be highly dynamic but assumptions and simplifications have to be adopted in the analytical treatment of \citet{Yuan09}.
In the present work, we perform three dimensional  magnetohydrodynamical numerical simulations to examine the  scenario proposed in \citet{Yuan09}, investigating the formation of episodic jets and their association with flares. The paper is structured as follows. In {\S\ref{sec:model}
, we briefly describe our model of general relativity magnetohydrodynamical simulations of hot accretion flows. In {\S\ref{sec:floats}
, we present the main results obtained from our simulations, including the formation of flux rope, its ejection, and the periodicity of the whole processes. In {\S\ref{sec:observation}
, we discuss the applications of our model to observations, including the periodic flares detected in Sgr A* and other sources, the hot spot in Sgr A* detected by the GRAVITY observations, the velocity and the possible periodicity of the ejected blobs in Sgr A*, M81, and M87, and finally the periodic ejection and flares in protostellar accretion systems.
The comparison between our present work and previous ones are presented in {\S\ref{sec:comp}}. Finally, we summarize our work in {\S\ref{summary}}

\begin{figure*}[htbp!]

   \plottwo{figures/figs1a.jpeg}
   {figures/figs1b.jpeg}
    \caption{Snapshots of properties of SANE and MAD in the $x-z$ plane. The color indicates {plasma} $\beta$, black lines denote magnetic field lines .}
    \label{fig:figs1}
\end{figure*}

\section{Model} \label{sec:model}
All black hole sources mentioned in the introduction section are powered by hot accretion flows \citep{Ho2008,YuanNar14}. We therefore perform a set of three-dimensional general relativity magnetohydrodynamical numerical simulations of a hot accretion flow around a black hole using the ATHENA++ code. Spherical coordinates are adopted. We have simulated two kinds of accretion flows, namely SANE (standard and normal evolution) and MAD (magnetically arrested disk), with the magnetization of the accretion flow in MAD much stronger than SANE \citep{Igum03,Nar03,Tchek11}. Two values of black hole spin are considered, $a=0$ and $a=0.98$. The details of the simulations are presented in Appendix A. We find that the value of $a$ is not important for our problem, so in the following we simply refer the models as SANE and MAD. Snapshots of the distributions of magnetic field lines and plasma $\beta$ (defined as the gas-to-magnetic pressure ratio) of the models are shown in Figure~\ref{fig:figs1}.

\section{Results} \label{sec:floats}
\subsection{Formation  of flux ropes}
 Solar flares and CMEs have been intensively studied, and the physical mechanism is known to be magnetic reconnection \citep{Zhang2005,Gou2019}. After reconnection, some plasma will be enclosed within helical magnetic field lines and form a “flux rope’’. The flux rope will then be ejected by strong magnetic pressure force and form the CMEs. It was proposed in \citet{Yuan09} that similar processes should also occur in the black hole accretion flows and result in the formation of episodic jets.

\begin{figure*}
\includegraphics[width=.305\linewidth]{figures/fig2-1.jpeg}\quad\includegraphics[width=.305\linewidth]{figures/fig2-2.jpeg}\quad\includegraphics[width=.39\linewidth]{figures/fig2-3.jpeg}
\includegraphics[width=.305\linewidth]{figures/fig2-4.jpeg}\quad\includegraphics[width=.305\linewidth]{figures/fig2-5.jpeg}\quad\includegraphics[width=.39\linewidth]{figures/fig2-6.jpeg}

\includegraphics[width=.305\linewidth]{figures/fig2-7.jpeg}\quad\includegraphics[width=.305\linewidth]{figures/fig2-8.jpeg}\quad\includegraphics[width=.39\linewidth]{figures/fig2-9.jpeg}

\caption{The poloidal magnetic field structure of a flux rope identified in MAD and at simulation time of T=20000{$r_g/c$} at various azimuthal angles $\varphi$ (denoted by the grid cell number $i$ at the top of each plot) in the $x-z$ plane. The magnetic field is denoted by black solid lines with arrows indicating the direction of the field while the color indicates the logarithmic value of the ratio of poloidal and toroidal magnetic field strength. The flux rope is located at $(x,z)=(5-10,5-10)$, depending on the value of $\varphi$.  In this example, the flux rope extends in the $\varphi$-direction for roughly $(63-41)\times5.6=123^{\circ}$.}
\label{fig:fluxrope}
\end{figure*}

\begin{figure*}
\includegraphics[width=0.65\columnwidth]{Figsmiki/densstrSane00T17400}
\includegraphics[width=0.65\columnwidth]{Figsmiki/densstrSane00T18400}
\includegraphics[width=0.65\columnwidth]{Figsmiki/densstrSane00T19400}
\caption{Six reconnection layers from our SANE simulation (marked with blue circles). They will evolve and result in the occurrence of magnetic reconnection and formation of six magnetic islands in $100r_g/c$, shown by the first column of Figure \ref{fig:isla}.  The density is shown by colour and magnetic field by solid lines, with arrows indicating the direction of the field.}
\label{islarec}
\end{figure*}

\begin{figure*}
\includegraphics[width=0.65\columnwidth]{Figsmiki/densstrSane00T19400b}
\includegraphics[width=0.65\columnwidth]{Figsmiki/densstrSane00T19430b}
\includegraphics[width=0.65\columnwidth]{Figsmiki/densstrSane00T19480b}
\caption{Evolution of the magnetic field structure around a reconnection layer shown in the rightmost panel of Figure \ref{islarec}, i.e, at $t=19400r_g/c$. Formation of magnetic island and its further motion can be found in the bottom panels of Fig. \ref{fig:isla}.  }
\label{islarec2}
\end{figure*}

 Stimulated by the above scenario, we have analyzed the simulation data of accretion flows, looking for flux ropes.  We do find that many such helical structure magnetic field lines (i.e., the flux ropes) exist both within and at the coronal region of the accretion flow, from radius as small as 4$r_g$ up to $\sim 30~r_g$, where $r_g$ stands for gravitational radius. The flux ropes are found to keep forming throughout the evolution of the accretion flows in all four models once they reach the steady state.

 In the present paper, we have tried to find as many as possible flux ropes. We then measure the properties of individual flux ropes and present the typical results.

 One example of flux rope identified in the case of MAD at time T=20000~$r_g/c$ is shown in Figure~\ref{fig:fluxrope}. The poloidal magnetic field structure  is presented in the $x-z$ plane at nine various azimuthal angles\footnote{ The time-averaged magnetic field configuration at a larger scale will show an ``hourglass'' shape, as shown by Figure 5 in \citet{Yang21}.}. As we can see from the figure, in the two-dimensional plane, these flux ropes will be represented by circular magnetic field line patterns --- “magnetic islands”. Note that if a constant-$\varphi$ slice through the flux rope were not perpendicular to the axis of the rope, we would see a whirl in poloidal field lines instead of a circular magnetic island.

\begin{figure*}
\includegraphics[width=2\columnwidth]{figures/phizratBpBphiSANET18500shifted}
  \includegraphics[width=2\columnwidth]{figures/phizplanedensSANET18500shifted
    \caption{The spatial distribution of the magnetic field strength and density when a flux rope is present in the SANE model. Top panel: Distribution of the ratio of poloidal and toroidal magnetic field strength in the $\varphi-z$ plane with a constant cylindrical radius $r= 30.29$ above the equatorial plane at $T = 18500~r_g/c$. The poloidal field reaches a local minimum at the center of the flux rope, which is shown as a slice through the azimuthal extensions of the magnetic islands. Bottom panel: Distribution of density in the same $\varphi-z$ plane. It is clear that the density of gas in the flux rope ({\bf middle of which is located around $(\varphi,z)=(25^\circ ,20)$ in this case}) is higher than the surrounding medium. Note that only parts of the arcs of the flux ropes lie in this $\varphi-z$ plane, because the cylindrical radius of the flux rope is not constant.  }
    \label{fig:banddensity}
\end{figure*}

To understand the physical mechanism of the formation of flux ropes, after we find a flux rope, we have traced back to the time and location of the flux rope formation. The motion of the flux ropes is traced by examining their movement in different $r-\theta$ planes with various angles $\varphi$. The center of the flux rope is identified by the minimum of the poloidal magnetic field and the circular pattern of the magnetic field lines, as we will explain at the end of section 3.1. Figures \ref{islarec} \& \ref{islarec2} show the tracing results. Three panels in Figure \ref{islarec} show the magnetic field structure in the $r-\theta$ plane at three snapshots. In each panel we can see two reconnection layers, located above and below the equatorial plane of the accretion flow respectively\footnote{We note that the appearance of reconnection layer and subsequent occurence of magnetic reconnection are found to occur randomly, not necessarily simultaneously above and below the equatorial plane.}. These reconnection layers will continue to evolve, result in the occurrence of magnetic reconnection and formation of magnetic islands after a time of $\sim 100r_g/c$. As an example, the evolution process of the reconnection layer located above the equatorial plane in the rightmost (i.e, $t=19400$) snapshot of Figure \ref{islarec} is illustrated in Figure \ref{islarec2}. The six magnetic islands finally formed are shown by the first column of Figure \ref{fig:isla}. We emphasize that similar steps of evolution from the reconnection layer to the magnetic island can be traced in all the cases of island formation captured in our simulations. It is interesting to note that similar to our case, magnetic reconnection is also invoked to explain the formation of solar prominence \citep{Chenpf2020}.

We speculate that the physical mechanism underlying the above processes is likely as follows.  In the inner region of the accretion flow, the field lines are tangled, due to the MHD turbulence driven by magnetorotational instability \citep{Balbus1991}. Thus the field lines with opposite polarity can come close enough, resulting in the occurrence of reconnection. In addition, the Parker instability and wind launched from the accretion flow will make the magnetic field lines emerge out of the accretion flow into the corona. The differential rotation and turbulent motion of the accretion flow where the field lines are rooted twist the field lines in the coronal region and make the reconnection occur there as well. This scenario of the formation of flux ropes is fully consistent with that proposed by \citet{Yuan09}. It has also been confirmed by other recent MHD numerical simulations \citep{Nathanail20,Nathanail2021,Ripperda20,Ripperda2021}, and it is confirmed again by our present analysis. The detailed comparison with these works will be presented in Section \ref{sec:comp}.

We find that, the formation of flux rope occurs not only in the coronal region, but also within the accretion flow, as also pointed out in \citet{Nathanail20}.  The reconnection rate is proportional to the Alfven speed. Since $\beta \gg 1$ in the accretion flow, the reconnection rate is much smaller there compared to the coronal region. Combined with the fact that the energy density of the magnetic field in the corona is also stronger than that in the accretion flow \citep[e.g.,][]{Yang21}, we expect that the energy release
during reconnection occurs slowly and weakly thus may not correspond
to observed strong flares. Because of these reasons, in the present paper, we only focus on the
flux ropes formed in the coronal region of the accretion flow.

Because of the reconnection, the poloidal magnetic field at the center of the flux rope should become weaker compared to its surrounding medium. This is confirmed by the top panel of Figure~\ref{fig:banddensity}. In fact, the movement of flux ropes is followed by tracing the trajectories of those minima and the circular magnetic field line pattern. The released magnetic energy should be converted into the thermal and kinetic energy of plasma, and be used to accelerate electrons whose radiation will be responsible for flares. In addition, the density within the flux rope is larger than its surrounding medium, as shown in the bottom panel of Figure~\ref{fig:banddensity}. We speculate the reason as follows. The weakening of the magnetic field results in the decrease of magnetic pressure within the flux rope; consequently the pressure from the surrounding medium compresses the flux rope and increases its density and pressure until a new pressure balance within and outside of the flux rope is established.  Combining
the top and bottom panels of the figure, we can find that the plasma $\beta$ within the flux rope should be higher while the magnetization parameter $\sigma\equiv B^2/(\rho c^2)$ should be lower than its surrounding medium. These results are consistent with \citep{Ripperda20}.

To present quantitative results, the change of various components of the magnetic field and density as a function of $\varphi$ at a given cylindrical radius are given in Figure \ref{fig:densityB}. Note
that since the cylindrical radius of the flux rope is different for different $\varphi$, in this figure, the flux rope is present only at two places, i.e., at $\varphi\sim 100^{\circ}$
and $\varphi\sim 350^{\circ}$. At these two
places, the poloidal field reaches its minima while the density reaches its maxima.

\begin{figure}
  \includegraphics[width=1\columnwidth]{figures/bandrhoalongphiT18500.png}
    \caption{The $\varphi$-profiles of various components of magnetic field and density of the flux rope shown in Fig. \ref{fig:banddensity} at about $z =27r_g$
in SANE, with the star mark denoting the same
magnetic island shown in Fig. \ref{islarec}. Note that since the cylindrical radius of the flux rope at different $\varphi$ is not a constant, in this figure the flux rope is present only at two places.}
    \label{fig:densityB}
\end{figure}

\begin{figure}
  \includegraphics[width=1\columnwidth]{figures/spiralkaa.png}
    \caption{The movement of a flux rope (i.e., a magnetic island) in the SANE simulation, denoted by slices crossing the center of the flux rope at different times and azimuthal angle $\varphi$.  We can see that with the increase of the radial distance of the island from the black hole, its corresponding $\theta$-angle is decreasing, indicating a  ``collimation'' effect.}
    \label{fig:spiral}
\end{figure}

\begin{figure}

  \includegraphics[width=1\columnwidth]{figures/fig3a.jpeg}
    \caption{The time-dependent location (represented by spherical radius) of a plasmoid in the case of MAD.}
    \label{fig:movement}
\end{figure}

\begin{figure*}

   \plottwo
      {figures/fig3b.jpeg}
      {figures/fig3c.jpeg}
    \caption{ (Left:) Various velocity components in units of speed of light of the plasmoid shown in Fig. \ref{fig:movement} as a function of time. (Right): Various velocity components in units of speed of light of a plasmoid in the case of SANE as a function of time.}
    \label{fig:blobvel}
\end{figure*}

\subsection{Ejection of flux ropes}

Although flux ropes can be formed at both small and large radii, we find that there is an important difference between them.  At small radii, within roughly $r\sim 10-15\, r_g$, although the flux ropes can be formed, few of them are ejected out. Many flux ropes just stay within the accretion flow and finally fall onto the black hole. This is likely because of the strong gravitational force of the black hole, similar to the existence of the ``stagnation radius’’ within which the motion of the matter is always toward the black hole. This is indicated by the direction of the velocity vectors close to the black hole shown in \citet{Yang21} (refer to Fig. 3 therein). Such a result is also consistent with the absence of wind in the innermost region of the accretion flow \citep{Yang21}. For other flux ropes, although they are found to be able to propagate outward for some distance, they quickly disappear because of the strong differential rotation of the accretion flow at small radii which tears the flux ropes.

Beyond 10-15 $r_g$, however, we find that the flux ropes  can be ejected out. This is because both the gravity and the differential rotation become weaker. The flux rope is found to become entangled soon after its ejection, because its two end-points are anchored in two different radii of the accretion disk which have different angular velocities. Such an entanglement would lead to the ejection of the flux rope due to the kink instability within an orbital timescale, similar to the solar and space physics cases \citep{Wang2016,2021JPlPh..87f9016S}. The short lifetime is consistent with previous analytical argument \citep{Broderick2006} and perhaps also related to the short lifetime of three flares detected by high resolution GRAVITY observations \citep{Gravity2018}.

We trace the motion of the broken flux rope (i.e., plasmoid) by examining its movement. The result is shown in Fig. \ref{fig:spiral}. We can see from the figure that the plasmoid is spiraling outward, reminding us of the helical structure of magnetic field in a continuous jet. It is interesting to note that the corresponding  $\theta$-angle of the plasmoid decreases with the increasing distance to the black hole, indicating a collimation effect, which is again similar to the behavior of continuous jets.

To give a quantitative result, we show by Figure \ref{fig:movement}  the radial location of a plasmoid in MAD as a function of time. The left panel of Fig. \ref{fig:blobvel} shows various components of the velocity as a function of time for this plasmoid.  It reaches as high as 0.3 speed of light at T=23000~$r_g/c$ from 0.05 speed of light at T=20000~$r_g/c$. As a comparison, we show by the right panel of Figure \ref{fig:blobvel} the evolution of various component of the velocity in the case of SANE. Compared to MAD, the typical values of both the magnitude of velocity and acceleration are much smaller.
\begin{figure}

  \includegraphics[width=1\columnwidth]{figures/fig4b}
    \caption{Forces along the radial direction acting on the magnetic island shown in the left plot of Figure~\ref{fig:blobvel}. The values of $\theta$ and $\varphi$ are indicated atop the upper right corner of the figure. The vertical solid line denotes the location of the magnetic island. We can see from the figure that the Lorentz force is the dominant force that accelerates the motion of the flux rope. }
    \label{fig:force}
\end{figure}

To understand the physical reason for the acceleration of the plasmoid, we have analyzed the forces acting on the flux rope shown in Figure~\ref{fig:movement} (and the left panel of Figure~\ref{fig:blobvel}). They include gravitational force, centrifugal force, Lorentz force, and the gradient of gas pressure. The result is shown in Figure~\ref{fig:force}. We can see from the figure that the Lorentz force is the dominant force that accelerates the plasmoid. The dominant component of the Lorentz force is the gradient of the magnetic pressure, which is enhanced by the occurrence of reconnection below the plasmoid. From this result, we can understand why the acceleration in the SANE case is much weaker than MAD; it is because the plasma $\beta$ in SANE is more than ten times larger than in MAD (refer to Figure~\ref{fig:figs1} and Figure 4 in \citet{Yang21}). Such a scenario of the magnetic acceleration of the ejected plasmoid confirms the prediction in \citet{Yuan09}, and is very similar to the case of solar CMEs \citep{Zhang2005}.

\begin{figure*}[htbp]
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T17500}
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T17700}
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T17900}\\
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T18500}
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T18800}
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T19100}\\
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T19500}
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T19800}
\includegraphics[width=0.55\columnwidth]{Figsmiki/densstrSane00T20100}

\caption{{Three successive cases of magnetic islands formation at T=17500, 18500 and 19500$~r_{\mathrm g}/c$, just after the reconnection events shown in Fig.~\ref{islarec}, and their subsequent evolution. The data is based on SANE simulation. Colour denotes density while solid lines magnetic field. The dashed black and red lines are spherical coordinate lines passing through the centers of two magnetic islands located at the crossings of the circular and straight dashed lines of the same color. The value of $\varphi$ is indicated atop the upper right corner of each panel. }}
\label{fig:isla}
\end{figure*}

\begin{figure}[t]

     \includegraphics[width=.48\textwidth]{figures/fig5a.jpeg}
     \includegraphics[width=.48\textwidth]{figures/Fig5b.jpeg}

     \caption{Properties of magnetic field at two different evolution times. The time (“T”) and the number of grid in the $\varphi$-direction (“i”) are shown at the top of each plot. The magnetic field is denoted by black solid lines with arrows indicating the direction of the field. The color indicates the logarithmic value of the ratio of poloidal and toroidal magnetic field strength. In the top panel, two flux ropes are located at $(x,z)=(8,7)$ and (20,28); in the bottom panel two flux ropes are located at (21,12) and (35,20). The magnetic field in the region where the two flux ropes are located is poloidally and toroidally dominated in the top and bottom panels, respectively. At the very center of the two flux ropes in the bottom panel, the poloidal field is weaker than the toroidal field because of the occurrence of reconnection (see text for details). The figure is plotted based on the simulation data of MAD.}

     \label{fig:fig5}
\end{figure}

\subsection{Periodicity of the formation and ejection of flux ropes}

We find from the simulation data that, in both SANE and MAD cases, every time when magnetic islands
are formed, their corresponding azimuthal angles $\varphi$ are usually different; but their corresponding radii are almost
always within the range of $20\la r\la 30 r_g$. More interestingly, the formation of flux rope seems to occur quasi-periodically, once the simulations have reached the steady state. Such a phenomenon does not depend on model parameters or numerical resolution. As an example, we show in Figure \ref{fig:isla}  three successive cases of magnetic island formation and their subsequent evolution. The formation time of the three islands is $T=17500, 18500$ and $19500${$~r_{\mathrm g}/c$}, respectively. More cases before $T=17500$ and after $T=19500${$~r_{\mathrm g}/c$} exist but are not shown in the figure due to the spatial limitation. This implies a period of $\sim 1000 r_g/c.$  Such a period matches the orbiting timescale of both
SANE and MAD at $20-30 r_g$ well, which are slightly larger than the Keplerian timescale since the angular velocities of both SANE and MAD are sub-Keplerian \citep{YuanNar14}.  Since the formation of flux rope at relatively large radii is accompanied by plasmoid ejection, such a periodicity implies that, in addition to flares, the ejection of plasmoids also occurs periodically. At smaller radii, the reconnection (thus flare) is also found to occur periodically, although not accompanied by ejections.

The presence of periodicity combined with its orbital timescale strongly implies that, compared to turbulence, the differential rotation of the accretion flow seems to play a more important role in twisting the magnetic field lines and resulting in the occurrence of magnetic reconnection. This scenario is enhanced by the fact that there seem to be no difference in terms of the periodicity of the
formation of flux ropes in SANE and MAD cases, although MAD is much less turbulent compared to SANE. In contrast, the occurrence of CMEs in the Sun mainly depends on solar activity. Although various periodicities exist for CMEs and flares \citep[e.g.,][]{Lou2003}, the periods are not related to the rotation period of the Sun. This is because, unlike the differentially rotating accretion flow, the geometry of Sun is spherical and the motion of the foot points of the field loops is driven by solar dynamo through local convection rather than by the rotation of the Sun.

In addition to the periodicity, another interesting finding is that in some cases, when a plasmoid moves outward, it may catch up and merge with a previously produced plasmoid having a slower speed. The merger of two magnetized islands must be accompanied by magnetic reconnection. This will results in  strong radiation so strong flares
are expected to occur. \citet{Yuan2012} have proposed this scenario as the mechanism of Gamma-Ray bursts.

\section{Applications to observations} \label{sec:observation}

\subsection{Periodic flares and the hot spot in Sgr A* and other black holes}

The formation of the flux rope occurs because of magnetic reconnection, which must be accompanied by electron acceleration. The synchrotron radiation of the accelerated electrons will be responsible for the observed flares in various wavebands, ranging from IR to X-ray and $\gamma$-ray flares in black hole sources mentioned above.

Near-infrared flares in Sgr A* were recently resolved by high resolution GRAVITY observations, with a compact ``hot spot’’ detected orbiting around the black hole with orbiting radius of 6-10 $r_g$ \citep{Gravity2018} . This hot spot is produced by  magnetic reconnection found in our simulations at small radii of the accretion flow. The flares are found to be polarized, with the polarization signature being consistent with orbital motion of a hot spot in a strong poloidal magnetic field \citep{Gravity2018}. In our simulation, we find that at the place of reconnection (i.e., the place of ``hot spot’’), which is also roughly where the flux rope is located, the magnetic field could be dominated by either poloidal or toroidal components, depending on time. This is shown by the two panels of Figure~\ref{fig:fig5}. The top panel of this figure shows a case that the field is dominated by the poloidal component, which is consistent with the observation. But at another time shown by the bottom panel of this figure, the field is dominated by toroidal component. This prediction could be verified by future observations. The detailed phenomenological modeling of the IR and X-ray light curves of Sgr A* indicates that to explain the simultaneity and symmetry of the X-ray and near-infrared light curves, the magnetic field strength is required to significantly decrease during the flare \citep{Dodds2010}. This is consistent with the weakening of the magnetic field during reconnection which we have described above and explicitly shown in Figure~\ref{fig:banddensity}.

From the over 2-hour long infrared light curve of Sgr A*, several periods of variability are identified by observations, with period ranging from 17 to 40 minute \citep{Genzel2003,Trippe2007,Genzel2010,Eckart2006}. The observed periods are consistent with orbital period at 3-5$r_g$ for the $3.6\times10^6$ solar mass black hole in Sgr A*, assuming the accretion flow in Sgr A* is described by an MAD. Given that the typical lifetime of a “hot spot” is shorter than one orbit timescale as we have explained before, the detected periodicity is best explained by the periodic occurrence of reconnection events found in our simulation at small radii.  The observed period is not a fixed value but has a range,  consistent with our finding that reconnection occurs in a range of radii.

\subsection{Ejected blobs in Sgr A*, M81, and M87: origin, velocity, and the possible periodicity}

In our simulations, the flux ropes (i.e., plasmoid) formed beyond 10-15 $r_g$ will be ejected outward from the surface of the accretion flow. This explains the origin of the radio-emitting blobs observed in Sgr A*, M81, M87, and other black hole sources. The observationally measured ejection velocity of 0.4, 0.5, and 0.2 of the speed of light for the ejected plasmoids in Sgr A*, M81 and M87 is also quantitatively similar to what we have obtained in the case of MAD, as shown in Figure~\ref{fig:blobvel}, which is about 0.3 of the speed of light. The ``MAD’’ nature of the accretion flow in M81 and M87 is consistent with most recent studies \citep{Shi2021,EHT2021,Yuan22}.

Examining the periodicity of the ejection of plasmoids is difficult since it requires very high spatial resolution and sensitivity of the telescope and the suitable observational targets. The high-resolution Global Millimeter-VLBI-Array observation to M87 jet reveals some discrete blobs at distance of about 100 $r_g$ from the central black hole \citep{Algaba2021}. By measuring the spacing between two adjacent blobs from the image and adopting a projected moving speed of 0.2 of the speed of light \citep{Hada2014}, we can estimate that the time interval of two adjacent ejections is about one year. This is in good consistency with our predicted period of 1000 $r_g/c \sim$ 1 yr for the $6\times10^9$ solar mass supermassive black hole in M87 \citep{Gebhardt2011}.

Most recently, two quasi-periodic oscillations (QPOs) are detected in the TeV blazar PKS 1510-089 by the Fermi Large Area Telescope \citep{Roy2021}. One is a 3.6-day QPO that lasted for five cycles with a moderate significance of $\sim 3.5 \sigma$, another is a 92-day QPO that lasted for seven cycles, with a significance of $\sim 7.0\sigma$. As discussed in \citet{Roy2021}, models like binary black holes and a procession jet are unlikely; two attractive models are a rotating hot spot in the inner region of the accretion flow and ``equispaced magnetic islands'' inside the jet. The remaining problem for the first model is that since the accretion flow in a blazar should have an almost face-on orientation, it may be hard for the rotating hot spot to produce enough variability.

Obviously, these two scenarios are very similar to what we have proposed in the present work. The only difference is that, for the first ``hot spot'' model, we think the variability is not caused by a single hot spot orbiting around the black hole, but several hot spots periodically appeared in the inner region of the accretion flow since the typical lifetime of a hot spot is shorter than one orbital time thus difficult to explain the observed several cycles.  Moreover, our ``modified'' hot spot model has the advantage of overcoming the ``not-enough variability'' difficulty.

\subsection{Possible application to protostellar accretion systems}

Although our simulations are done for hot accretion flows around black holes, we speculate that the basic physics relevant to flares and ejections should also work for accretion disks in protostars. This is why very similar phenomena have been observed in these systems. Specifically, in this case, the period should correspond to the orbital period at the inner edge of the accretion flow where it is truncated by the magnetosphere of the protostar, which is a few days. Periodic ejections are observed from the spacings of the jet knots \citep{Lee2020} in those systems, and the period can be estimated from the spatial intervals and the jet velocities. Unfortunately, the current instrumental limit of spatial resolution $\sim 0.1^{"}$ sets a limit on the period detection no shorter than about 1 yr \citep{Lee2020}. On the other hand, the Chandra light curve does show obvious periodicity with the period of $2\times10^5$ second or a few days, according to \citet{Wolk2005} and \citet{Flaccomio2012} (Ref. to Fig. 3 in \citet{Wolk2005} and Fig. 4 in \citet{Flaccomio2012}), consistent with our prediction.

\section{Comparison with previous works}
\label{sec:comp}

The most relevant previous work is \citet{Yuan09}. In that work, by analogy with the coronal mass ejection in the Sun, an analytical model was proposed to explain the formation of episodic jets and their association with flares. The basic scenario proposed in that work, i.e., flux ropes are formed due to magnetic reconnection and they are then ejected out by the magnetic forces, are fully confirmed by the present numerical simulations. Due to the limitation of analytical calculations, although the whole physical processes involved are highly dynamic, many assumptions and simplifications have to be adopted in \citep{Yuan09}. For example, the formation of flux rope is purely an assumption there. The quantitative calculation to the ejected velocity of the flux rope is also based on assumptions on the distributions of magnetic field and density of the coronal gas. Compared to \citet{Yuan09}, our current work is much more realistic and comprehensive.  Moreover, important new findings are obtained, e.g., the periodic formation of the flux rope and its ejection, and the quantitative prediction of the velocity of the ejected blobs.

Recently several MHD numerical simulation works focusing on the analysis of reconnection in accretion flows have been published. It is thus useful to compare our simulations with these works. In the two-dimensional general relativity ideal magnetohydrodynamic numerical simulations by \citet{Nathanail20}, the formation and evolution of current sheets are investigated. Special attention is paid to the effect of different initial magnetic field configurations. Similar to our work, they find that plasmoids are formed both within and at the surface of the accretion flow due to magnetic reconnection, which is driven by the turbulent motion of the accretion flow. Moreover, they find that plasmoids are formed during the reconnection, which is again consistent with our work.

In the works by another group, the two and three dimensional MHD simulations including weak Ohmic resistivity are performed, focusing on the modeling of magnetic reconnection and the formation of palsmoids \citep{Ripperda20,Ripperda2021}. Since an explicit resistivity is included which can approximately mimic the kinetic effects, their simulations are more suitable than ours to describe magnetic reconnection. Both SANE and MAD are studied as in our work. They find that within about $20r_g$ of the accretion flow current sheets and plasmoids are ubiquitous features that form regardless of the initial magnetic field, the magnetization of accretion flows, and the spin of the black hole.

Different from our work and \citet{Nathanail20}, another mechanism for the formation of plasmoid is identified in \citet{Ripperda20}, namely the tearing instability of the current sheet. When the reconnection layer becomes thin enough, the reconnection layer can break up and produce chains of plasmoids. These small plasmoids can merge and grow to macroscopic scales of the order of a few Schwarzschild radius, and be advected along the jet’s sheath or into the disk. We do not find such an instability in our simulations.  For this instability to occur, the Lundquist number $S=v_{\rm A} L/\eta$ needs to be larger than $10^4$. Here $v_{\rm A}$ being the Alfven speed, L being the typical length of the current sheet, and $\eta$ being the resitivity \citep{McKinney2012}. In our simulations, $S\lll10^4$   because of the high numerical resistivity caused by the relatively low resolution of our three-dimensional simulations. In reality both turbulence of the accretion flow and tearing instability in current sheet should be able to drive reconnection and the formation of plasmoids, as seems to be confirmed by a recent work \citep{Rosenberg2021}. In this work, resistive MHD equations are solved, reconnection and plasmoids are found to be formed both via merging current sheets due to turbulent motions of the accretion flow and through reconnection due to current sheet instabilities.

\section{Summary}\label{summary}

Episodic ejections  of blobs (episodic jets) have been observed in both black hole and protostellar accretion systems and they are often associated with radiation flares. Notable sources include Sgr A*, M81, and M87, among others. \citet{Yuan09} has proposed an analytical MHD model for episodic jet and its association with flares. In the present work, we develop that work by performing three dimensional GRMHD  numerical simulation of hot accretion flows. Both SANE and MAD are considered. We have analyzed the simulation data and obtained the following results.

\begin{itemize}
    \item[(1)] We find once the simulations have reached the steady state, flux ropes are keep forming both within and at the coronal region of the accretion flow, from radius as small as 4$r_g$ up to $\sim 30r_g$.  One example of a flux rope in the case of MAD is shown in Figure \ref{fig:fluxrope}. It extends in the $\varphi$ direction for about 120$^{\circ}$.

    \item[(2)] We find that the flux rope is formed due to the reconnection of magnetic field lines, driven by the differential rotation and turbulent motion of the accretion flow. This scenario is consistent with that proposed in \citet{Yuan09} and recent MHD numerical simulations. The detailed comparison between the present work and those previous ones are presented in Section \ref{sec:comp}.

    \item[(3)] We find that flux ropes formed at small and large radii have different fate. At small radii, inside of $\sim$ 10-15$r_g$, few flux ropes can be ejected out, likely because of the strong gravitational force of the black hole, similar to the absence of wind at small radii. Beyond that radius, flux ropes can be ejected out. The ejected flux ropes are found to break within one orbital timescale due to the kink instability.

    \item[(4)] The broken flux ropes (i.e., plasmoids)  spiral outward, with radial velocities shown in Figure \ref{fig:blobvel}. A clear acceleration in both cases is found and the  velocity can reach $0.3c$ and $0.08c$ for MAD and SANE, respectively. We have analyzed the forces and found that the acceleration is mainly due to the gradient of the magnetic pressure, as shown by Figure \ref{fig:force}, again consistent with \citet{Yuan09}.

    \item[(5)] The whole processes mentioned above, i.e., the formation of the flux ropes and their subsequent ejection, occur periodically, with the period being the orbiting time at the radius where the flux ropes are formed. This implies that differential rotation of the accretion flow plays a more important  role than turbulence in twisting the magnetic field lines and result in the formation of reconnection.

    \item[(6)] Magnetic reconnection should accelerate some electrons, their synchrotron emission explains observed flares and the origin of the ``hot spot'' in Sgr A* detected by GRAVITY.

    The periodicity detected from the IR light curves of Sgr A* and $\gamma-$ray light curves of blazar PKS 1510-089 are explained by the periodic formation of flux ropes.

    \item[(7)] The observed ejection of blobs is explained by the ejection of flux ropes formed at relatively large radii. In Sgr A$^*$, M81, and M87, the measured velocity of the ejected blobs agrees with the velocity predicted in our simulations of MAD. The predicted periodicity of the  ejection seems to be supported by some observations of M87 and PKS 1510-089.

    \item[(8)] Finally, we speculate that our results may also be able to interpreting flares and episodic ejection observed in protostellar systems. The measured period in the X-ray light curve is roughly consistent with the orbital period at the inner edge of the accretion flow where it is truncated by the magnetosphere of the protostar.

\end{itemize}

\section*{Acknowledgments}
We thank the useful discussions with Drs. P.F. Chen, J. Lin, and Y.M. Wang on the solar flares and coronal mass ejections, and C. White for his help of using ATHENA++ code. The {anonymous} referee is acknowledged for
constructive suggestions and comments. This work has made use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory. We also thank ASIAA, Taipei for the use of their visualization servers. This work was supported by CAS President’s International Fellowship for Visiting Scientists grant 2020VMC0002 (MC), Polish NCN grant 2019/33/B/ST9/01564 (MC), Natural Science Foundation of China grants 12133008, 12192220, and 12192223 (FY, HY), and Ministry of Science and Technology of Taiwan grants 108-2112-M-001-009 and 109-2112-M-001-028 (HS)


\section{INTRODUCTION}

In the recent decade, there has been a growing demand for mobile robots for a person following in the human-robot collaboration works such as assisting robots in industry, healthcare, surveillance, elderly care, and monitoring \cite{islam2019person}. Conventional robot tracking systems used lasers, which provide a cluster of distance measurements, to perceive the environment as well as to detect and track the target person \cite{islam2019person}. However, laser data are usually noisy and difficult to use for re-identification of the target. With the advent of deep learning, there has been a rapid growth of visual tracking algorithms \cite{VOT_TPAMI, kristan2018sixth, kristan2019seventh, kristan2020eighth, Kristan2021ninth}, which makes robot person following using vision-based trackers one of the popular and intelligent approaches \cite{chen2017integrating, chen2017person}.
The primary challenge is to accurately localize the target person of interest and take necessary control action to follow. Inherently, the task is quite challenging because it has to continuously recognize the target person in the presence of occlusion, scale change, rotation, etc. Also, it has to carefully adapt to the changes in appearances of the target persons due to the lighting variations.

Currently, the prominent and intelligent trend in the vision-based person following is using robust visual object trackers, which are integrated into the robots to provide them with perception ability. So, the performance of a person following robot highly depends on the robustness and speed of the object tracker. Several object trackers are available in the literature that has made tremendous progress by leveraging the expressive power of the deep learning algorithms \cite{ma2015hierarchical, mayer2021learning, chen2020siamese, ma2020rpt, zhao2021trtr, yan2021learning, chen2021transformer}. Furthermore, their progress has been fueled by several popular benchmark datasets that are well established in covering diverse real-world tracking challenges like rotation, occlusion, low-resolution, small-scale objects, etc \cite{muller2018trackingnet, fan2019lasot, huang2019got, wu2013online}.

The currently available benchmark datasets rarely cover scenarios of person tracking in uniform crowd scenarios,i.e. tracking the specific person while people nearby him are wearing similar clothes. This is commonly seen in the places like industries, sports stadiums, and schools as shown in Fig. \ref{fig:sample}. These scenarios are also vastly found in the Gulf region and many Asiatic countries. Tracking a person in these environments based on the RGB image is a challenging task for the state-of-the-art (SOTA) visual trackers as they perceive the surroundings based on appearance features. Hence, tracking a person wearing the same garment in a uniform crowd becomes a difficult task compromising the performance of the person-tracking robot. To address this challenge,  we present a new tracking system that uses an attention-based module to learn deeper representative features from color and depth information.

\begin{figure}[t]
    \centering
    \includegraphics[width=.95\linewidth]{images/sample.pdf}
    \caption{Sample images showing uniform crowd. $1^{st}$ row shows common crowd appearance in real-life and $2^{nd}$ row shows person following experimentation in uniform crowd environment.}
    \vspace{-7mm}
    \label{fig:sample}
\end{figure}

We argue that in person tracking in the uniform crowd, the color information is not enough because the target and the non-targets all look the same. In addition, the problem becomes more challenging when the target is occluded by a similar-appearing distractor. Thus, to address these common problems and to come up with a robust person-following robot system, we integrate depth information into the robust tracking architecture\cite{yan2021learning}. The addition of extra depth information in the tracking pipeline allows the best  segregation of the background and the foreground irrespective of high similarity \cite{lukezic2019cdtb}. In summary, this work has three contributions:
\begin{itemize}
    \item We propose a fully autonomous robot person following system with our new RGB-D tracker. It can track the person under the challenges like uniform crowd and occlusion with similar-appearing people.
    \item We conducted intensive experiments to evaluate different SOTA trackers’ performance for a person following in a uniform crowd environment.
    \item We introduce a new person tracking dataset that includes the challenges of a uniform crowd.
\end{itemize}

\section{METHODOLOGY}
\label{sec:methodology}
\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{images/block_diagram.pdf}
    \caption{System Block Diagram}
    \vspace{-7mm}
    \label{fig:block_diagram}
\end{figure}

This section describes the pipeline and components of our designed robot person following system as shown in Figure. \ref{fig:block_diagram}.Our proposed robot person following system is a fully autonomous system containing two main components: the perception module (which includes face identification, person detection, and tracking) and motion control of the robot. Section \ref{subsec:perception} explains the process of how the robot continuously perceives the target person of interest, and Section \ref{subsec:control} illustrates the motion control of the robot.

\subsection{Robot Perception}
\label{subsec:perception}
The perception module enables a robot to perceive the environment and identify its objective continuously. The main aim of the module is to allow visual servoing to work autonomously with minimal wrong behavior. Several individual components work together to allow robotic perception to work seamlessly. They are an identification module that initially identifies the specific target person's body based on face recognition and the tracking module, which tracks the target person in the crowd.

\subsubsection{Person Identification}
\label{subsec:identification}
The faces are visible initially thus to recognize the target person's body in the scene, the overall algorithm involves three parts, face identification, person detection, and face-person matching. To achieve fast and robust person identification, we employed a pre-trained yoloV5Face \cite{qi2021yolo5face} to detect the faces and a mobile facenet \cite{chen2018mobilefacenets} for feature extraction. The mobile facenet embeds the face region into the 128-dimensional feature vector. Face verification is achieved by finding the smallest euclidean distance between the feature vectors of the known target face and detected faces with a threshold of 0.9.

Lastly, for locating the whole body of the target person in the scene, we have used a pre-trained yoloV5 detector \cite{glenn_jocher_yolov5} as a person detector to detect the people and then computed the Intersection Over Union between the regions of the detected people and the identified target face. The body with the highest IOU is considered the target person's body and is used for tracking.

\subsubsection{Tracking}
\label{subsec:tracking}
\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{images/Picture1.png}
    \caption{Proposed RGBD tracker fusing RGB and depth data}
    \vspace{-7mm}
    \label{fig:architecture}
\end{figure}

For effective person tracking in the uniform crowd environment, we utilize depth information along with RGB. We argue that RGB only is not sufficient to discriminate the visual features in a uniform crowd environment. We propose a new tracker model, which is named Deep Tracking with RGB-D (\textbf{DTRD}),  that fuses the RGB and depth images for tracking. The tracker finds discriminative features from RGB and depth together to enhance its ability to distinguish the foreground target from the highly similar background distractors. Also, the tracker architecture adds a transformer module to enhance the feature learning and give attention to a specific target of interest as shown in Fig. \ref{fig:architecture}. The model contains three main components: a convolutional backbone for feature extraction, a transformer module for attention-based feature learning, and a bounding box regression head to predict the location of the target.

\textbf{Backbone:} The backbone adopts widely used vanilla ResNet \cite{he2016deep} convolutional layers for feature extraction for its expressivity and generality. The input of the trackers contains two pairs of images: template image $z$ and search image $s$ both containing an RGB and depth image. The pair of RGB and depth images are merged into a 4 channel image, $z \in \mathcal{R}^{H_z \times W_z \times 4}, x \in \mathcal{R}^{H_x \times W_x \times 4}$ and then passed through the backbone.  After the backbone, the input images are projected to two feature maps $f_z \in \mathcal{R}^{H_z/s \times W_z/s \times C}$ and $f_x \in \mathcal{R}^{H_x/s \times W_x/s \times C}$ where $s, C , H, W$ are the stride, number of feature channel, height and width respectively.

\textbf{Transformer:} The transformer contains six encoder-decoder blocks and a positional embedding layer. Before being fed into the transformer module, the feature maps are flattened and concatenated into a feature vector with a length of $H_x/s \times W_x/s + H_z/s \times W_z/s $ and a dimension of C. The transformer embeds the feature vector into a one-dimensional vector with a length of $d$.

\textbf{Regression Head:} The regression head is three fully connected layers to predict the four coordinates of the bounding box. The regression head directly predicts the top-left corner and the bottom-right corner of the target bounding box.

\subsection{Motion Control}
\label{subsec:control}
The designed robot-person following system is a visual servoing system in which the object tracker plays the most important role in locating the target. Thus, to keep the robot following the target person in an adaptive manner we have utilized a simple PI controller as the robot's motion controller. Also, we argue that a good object tracker should be able to compensate for the controller to improve the robot person following. The PI controller is defined in equation \ref{eqs:picontrol}.
\vspace{-4mm}
\begin{equation}
    u(t) = K_p e(t) + K_i \int e(t)dt
    \vspace{-2mm}
    \label{eqs:picontrol}
\end{equation}

Where the $K_p, K_i$ are the control gains, $e(t)$ is the error between the current measure and required measure, and $u(t)$ is the velocity output. For both the linear velocity and angular velocity, two PI controllers have been used.

\section{EXPERIMENT}
\label{sec:experiment}
In this section, we detail the developed dataset, tracker training, and intensive experiments performed to compare and evaluate the robot person following systems in a uniform crowd environment in real-time. Section \ref{subsec:dataset} presents the developed dataset and our tracker training detail. Section \ref{subsec:realtime} explains the real-time robot person tracking experimentation setup. Section \ref{subsec:result} and  \ref{subsec:discussion} wraps up the results of the experiment with an insightful observation.

\subsection{Dataset and Training Detail}
\label{subsec:dataset}
The primary problem for our proposed tracker was the lack of a uniform crowd tracking dataset(a similar-appearing distractor keeps on distracting in the view) that also contains depth data. Therefore, 45 uniform crowd videos (avg. 1891 frames/video) with synchronized RGB and depth data were collected and properly annotated initially with CVAT\footnote{\href{https://cvat.org/}{\textit{https://cvat.org/}}}(shown in $2^{nd}$ row of Fig.\ref{fig:sample}). The dataset records several scenarios such as occlusion, illumination variations, motion blur, etc.

We trained our tracker with $70\%$ of the dataset with initial weights from the STARK \cite{yan2021learning} spatial model. To properly embed depth, the first convolution layer of the backbone was modified to take four-channel image input and was trained from scratch. The training was done for 10 epochs with IOU and L1 loss using an ADAMW optimizer with a learning rate of 0.0001 for the model and 0.00001 for the backbone. The trained tracker was evaluated on the remaining $30\%$ of the dataset to confirm tracking ability.

\subsection{Real-time Tracking}
\label{subsec:realtime}
To compare the system following the target, we selected 2D SOTA trackers(DiMP\cite{bhat2019learning}, ATOM\cite{danelljan2019atom}, SiamBAN\cite{chen2020siamese}, STARK\cite{yan2021learning}, KeepTrack\cite{mayer2021learning}) that only use color image and our proposed 3D tracker that uses color and depth image.

\subsubsection{Experiment Setup}
\begin{figure*}[t]
    \centering
    \begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/targetmovement.png}
    \caption{Target moving path}
    \label{fig:target_moving_path}
    \end{subfigure}
    \begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/scenario1.png}
    \caption{Single distractor}
    \label{fig:single_distractor_path}
    \end{subfigure}
    \begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/scenario2.png}
    \caption{Two distractors (cross walking)}
    \label{fig:two_distractor_cross_path}
    \end{subfigure}
    \begin{subfigure}{.245\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/scenario3.png}
    \caption{Two distractors (parallel walking)}
    \label{fig:two_distractor_parallel_path}
    \end{subfigure}
    \caption{Pre-defined path for experimentation (Top view). The red solid line is where the target is occluded by the distractor/s.}
    \vspace{-3mm}
    \label{fig:path}

\end{figure*}

We used the pioneer3 AT robot as the mobile platform with ZEDm RGB-D camera placed one meter above the robot and a laptop with GPU as the control unit. We used an OptiTrack system\footnote{\href{https://optitrack.com/}{\textit{https://optitrack.com/}}} to measure the actual positions of the robot and people in the scene.

To make an effective and fair comparison between different trackers, we set a predefined rectangular path of the target person's movement as shown in Fig. \ref{fig:target_moving_path} and a trial is complete once the person reaches its initial position. The initial positions of the robot and target person are fixed at the left bottom corner (of the Fig. \ref{fig:target_moving_path}) with the robot placed at the corner and the person is 2 meters in front of it.

For studying the impact of different levels of the uniform crowd on a person following system,  we conducted experiments by adding a different number of distractors in the scene to examine the system performance with different trackers. Designed experiments include no distractor scene which provides a baseline evaluation of the system; One distractor scene where distractor walk cross to occlude the target as shown in Fig. \ref{fig:single_distractor_path}; Two distractor scenes where distractors walk cross or parallel to intervene the robot as shown in Fig. \ref{fig:two_distractor_cross_path} and \ref{fig:two_distractor_parallel_path}.

To avoid discrepancies during metric evaluation the distractor stops once the robot starts following him while the target person continues its path. A complete experiment protocol has been conducted as in Algorithm \ref{alg:realtime_experiment} to robustly and fairly examine the system's ability to follow the target.
\vspace{-2mm}
\begin{algorithm}
\caption{: Real-time experiment protocol}
\begin{algorithmic}[1]
\For{subject = A, B}
    \For{tracker = DiMP, ATOM, SiamBAN, STARK, KeepTrack, DTRD}
        \For{distractor = 0, 1, 2}
            \For{trial = 1, 2, 3}
                \State Perform robot person tracking experiment
            \EndFor
        \EndFor
    \EndFor
\EndFor
\end{algorithmic}
\label{alg:realtime_experiment}
\end{algorithm}
\vspace{-3mm}

\subsubsection{Experiment Evaluation Metric}
\label{subsec:experiment_evaluation_metric}
In previous research, we found the majority of work used distances error and trajectory matching to quantitatively and qualitatively evaluate person following. These metrics are insufficient to evaluate the systems' performance, especially in the case of the intervention from the distractors. Hence, we also add a new metric, following success, in the experiment for robust and quantitative evaluation as listed in Table \ref{tab:metrics}.

\begin{table*}[!t]
\centering
\caption{Evaluation metrics for robot person following}
\label{tab:metrics}
\begin{tabular}{p{0.25\textwidth}|p{0.65\textwidth}}
\hline
\multicolumn{1}{c|}{Metrics} & \multicolumn{1}{c}{Descriptions}                                                                                \\ \hline
Distance Error (DE)          & The difference in measurement from the tracker and the ground truth from the OptiTrack. Higher is the worst.                                  \\ \hline
Following Success (FS)       & Duration of the robot is successfully following the target person.  FS = {[}0, 1{)}. Higher is the best.                                                 \\ \hline
Frame rate Per Second (FPS)  & Frame processing time for the tracker which indicates the system's response time.
\\ \hline
\end{tabular}
\vspace{-3mm}
\end{table*}

\subsection{Results}
\label{subsec:result}
We report results on four evaluation metrics based on intensive experimentation.

\begin{figure*}[t]
    \centering
    \begin{subfigure}{.3\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/distanceerror.png}
    \caption{Distance error}
    \label{fig:distanceerror_curve}
    \end{subfigure}
    \begin{subfigure}{.3\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/FS_Curve.png}
    \caption{Person Following Success}
    \label{fig:fs_curve}
    \end{subfigure}
    \begin{subfigure}{.3\textwidth}
    \centering
    \includegraphics[width=\linewidth]{images/fps.png}
    \caption{Frame Processing Time}
    \label{fig:fps_curve}
    \end{subfigure}
    \caption{Experiment results}
    \vspace{-5mm}
\end{figure*}

\noindent\textbf{\textit{DE} Result:}
$DE$ quantifies how well the tracker predicts the target's state and signifies the error in distance measured from the tracker to the actual distance measured from the OptiTrack. From Fig. \ref{fig:distanceerror_curve} we can see that the overall $DE$ is around $0.5 m$ when there is no distractor in the field. But for all the SOTA trackers the errors have increased significantly with distractor presented scenarios making an average distance error greater than $1.2 m$. However, in our proposed DTRD tracker, it is around $0.6 m$ only for the two-distractor(cross) scenario. Whereas for two-distractor(parallel) it is slightly higher than other scenarios but still has the lowest error compared to all other trackers.

\noindent\textbf{\textit{FS} Result:}
$FS$ evaluates how well the tracker can drive the robot following the target person. This metric is calculated by comparing the distances traveled by the robot and the target person. The metric is derived based on the fact that, if the robot losses the target early, the distance traveled by the robot is short resulting in a lower $FS$ value. However, a higher $FS$ value signifies that the tracker is providing a robust perception to the robot despite any distractions. In Fig. \ref{fig:fs_curve}, when there are no distractors the $FS$ value for all the trackers is very high. However, when there are distractors all the 2D trackers have an $FS$ value of around 0.2 indicating they lost the target very early. But our proposed tracker has a higher $FS$ value than other trackers, showcasing robustness to the distractor.

\noindent\textbf{\textit{FPS} Result:}
This metric reveals the frame processing time of trackers which has a great impact on the response time and stability of the robot in real-time ($FPS$ should be above $20$). From Fig. \ref{fig:fps_curve}, we can see that all of the trackers fulfilled the constrained, except for KeepTrack (avg $FPS$ $13$).

\subsection{Discussion}
\label{subsec:discussion}
From the experimental results, it is evident that the visual-based robot person following system is able to keep track of the target person, referring to the reported higher performance of all trackers in the no-distractor scenario. However, the system performance significantly degrades when the distractors intervene. But our proposed tracker provides better perception ability than other 2D trackers. The reasons are two folds. Firstly, all the SOTA 2D trackers are only relying on appearance thus failing in uniform crowd scenarios and our DTRD tracker improved performance by exploiting the depth as well as using an attention-based module to learn better features from the images. Secondly, we observed that our DTRD tracker performs well in the uniform crowd when there is some distance between the target and distractors, but can get distracted otherwise. This is evidenced by the $FS$ value (around 0.6) of the DTRD in two-distractor scenarios showing the DTRD resisted the distractor's first intervention but failed at the second one because the gap was small.

Nevertheless, a more advanced design of trackers is required in the future to overcome the challenges of a uniform crowd environment. Also, trajectory prediction using visual cues is needed to make the robot handle the target out-of-view situation. It is to note that there is a great potential for depth to solve the challenges of tracking.

\section{Conclusion}
\label{sec:conclusion}
In this work, we present a new person tracking method exploiting the RGB and depth data that can overcome the deceptive challenge introduced by similar-appearing distractors. The proposed system is robust to different challenges of occlusion and distraction by the uniform crowd and thus provides resilient perception ability to the robot. Our tracking method also outperforms the SOTA 2D trackers in three different robust quantitative evaluation metrics.

\bibliographystyle{IEEEtran}


\section{Introduction}

Relation between brain function and cortical folding has long been intriguing due to its importance in offering profound insight into brain functional operation mechanisms, brain development and evolution, as well as being of clinical applications \cite{troiani2020use}. Gyri and sulci, the convex and concave cortical folding patterns, are two basic anatomical units of cortical morphology. It has been found that gyri/sulci align well with cortical functional boundaries \cite{fischl2008cortical}, at least in primary functional cortex. Specific gyro-sulcal patterns, such as the $plis$ $de$ $passage$ gyri \cite{regis2005sulcal}, topology of local gyri and sulci (such as H-shape sulci in frontal lobe \cite{troiani2020use}, H-shape gyri in fusiform gyri \cite{cachia2018interindividual} and the ventromedial prefrontal cortex (vmPFC) in DMN region \cite{lopez2019human}), were suggested to correspond to subtypes of a specific function (such as cognitive and perception ability associated with frontal lobe and reading ability associated with fusiform gyri) or relate to the placement of a specific functional region (such as functional activation location in vmPFC). In \cite{liu2019cerebral}, functional MRI signals were suggested to have stronger low frequency components on gyri, whereas sulci have signals of stronger high frequency components. Neuronal difference in their distributions \cite{hilgetag2005developmental}, types \cite{rash2019gliogenesis}, wiring patterns \cite{xu2010axons} and even morphology \cite{hilgetag2005developmental} were reported to be present between gyri and sulci, which was further suggested to lead to differentiated cell-level functioning patterns. In addition, expression of a huge number of genes was reported to be different between gyri and sulci \cite{de2015discrete}. Among them, interestingly, the same group of genes can lead to the formation of gyri/sulci on gyrencephalic species while mediating brain function distribution patterns on lissencephaly species \cite{de2015discrete}, suggesting a possible deeply rooted bond between cortical folding pattern and brain functions. Importantly, the functional differentiation between gyri and sulci could also provide clues to artificial neural network design. For example, gyri were suggested to serve as functional hubs and sulci as local processors, whereas sulci were isolated by gyri on the cortical sheet
\cite{deng2014functional}. Recalling that graph architecture of an artificial neural network in the sweet spot resembles the real neural networks \cite{you2020graph}, such a gyro-sulcal (hub-nonhub) spatial distribution might suggest an alternative approach for neural network design.

In spite of the importance in anatomy-function relation, no explicit explanation or generic theory of such a relation can be satisfactorily applied across brain regions and species \cite{fischl2008cortical}\cite{frost2012measuring}. Also, very few computational models are found \cite{liu2019cerebral} to provide a precise prediction from brain function to cortical folding. The reasons could be attributed to many folds: 1) the mapping between brain function and folding pattern is far from linear; 2) inter-individual variabilities are huge in the layout of both gyro-sulcal patterns and functional regions \cite{fischl2008cortical}\cite{frost2012measuring}; 3) brain function was usually delineated by atlases where “abrupt” boundaries were used to segregate functional “units” or “networks”, the functional attributes within which were considered to be of the same. However, the transition of a brain function from one to another was suggested to be “smooth” and “gradual”, which also gains supports from the related genetic expression \cite{huntenburg2018large}; 4) the spatial patterning of brain function distribution was not fully used to predict cortical folding.

To answer the aforementioned questions, as a preliminary effort, we proposed to used artificial neuronal networks to investigate the predictive ability of brain functional gradients to cortical folding patterns at macro scale. By means of resting-state fMRI (rsfMRI) and structural MRI (sMRI) in Human Connectome Project (HCP) dataset \cite{van2013wu}, we estimated the global functional gradients, embedding axes that encode the “gradual” differences in mesh vertices' connectivity patterns and can be present on a cortical sheet (top-left in Figure 1), where the prediction of gyro-sulcal segmentation (top-right in Figure 1) was performed. Technically, we proposed to use a mesh convolution model\cite{hu2022subdivision} based on the classic U-net architecture, which has been widely and successfully applied to medical image segmentation and is a close analogy with our
task. In addition, to enhance the interpretability of the mapping from functional gradients to cortical folding maps, we added a channel attention block \cite{hu2018squeeze} to the head of the U-net mesh convolution model, where functional gradients, as multi-channels, were re-weighted by the attention block before being fed to mesh convolution. The weights can encode the contribution of functional gradients to cortical folding prediction. Overall, the main contributions of this work are:
\begin{itemize}
    \item We proposed a new concept of adopting functional gradients to encode functional connectivities and present them on the cortex, such that mesh convolution can be adopted to predict gyri/sulci from spatial patterns of functional connectivity embedding.
    \item Channel attention was designed to enhance the interpretability of the nonlinear mapping between brain function and cortical folding.
    \item We explored and defined the activation maps of different layers and found novel clues to function-anatomy relation, which is valuable for both real and artificial neural networks.
\end{itemize}

\section{Related Works}
\label{gen_inst}

\subsection{Mesh Convolution}
Meshes are composed of three distinct types of geometric primitives: vertices, edges, and faces. The classification on meshes is a basic research area which has made great progress with the development of deep learning. DiffusionNet \cite{sharp2022diffusionnet} introduces a general-purpose approach to deep learning on 3D surfaces. The networks can be discretized on various geometric representations such as triangle meshes or point clouds, and can even be trained on one representation and then applied to another. HodgeNet \cite{smirnov2021hodgenet} uses many features that make it an attractive alternative for learning from meshes. During inference, its structure resembles that of most spectral geometry processing algorithms: construct a useful operator, and compute features from its spectrum. However, the above two approaches, as well as many other methods, including Pointnet \cite{qi2017pointnet}, MeshCNN  \cite{hanocka2019meshcnn}, PD-MeshNet \cite{milano2020primal}, MeshWalker \cite{lahav2020meshwalker}, PFCNN \cite{yang2020pfcnn}, need geometric features of meshes, such as coordinates of vertices, which is not applicable to our task that learns and predicts "texture" features on meshes. BrainSurfCNN \cite{ngo2020connectomic} and SubdivNet \cite{hu2022subdivision} provides a solution to the problem. BrainSurfCNN is a method which uses spherical convolutional kernel \cite{jiang2019spherical} to predict the task fMRI contrasts from resting-state functional connectivity ones under a regression framework. The work in \cite{jiang2019spherical} uses only one-ring neighborhood for convolution, which limits the extraction of long-range features. SubdivNet \cite{hu2022subdivision} uses a regular and uniform downsampling scheme to establish a fine-to-coarse mesh hierarchy. The convolutions efficiently support stride and large dilation, allowing the model to better capture long-range features. At the same time, the SubdivNet is a general approach, in which one can easily use 3D CNN on the mesh as similar as using the 2D CNN on the image without geometric features. In our method, we choose the convolution operation of SubdivNet to learn the "texture" features, rather than geometric features, for the purpose of predicting cortical folding patterns from brain cortical functional gradients.

\subsection{Attention Mechanisms}
Attention mechanisms were introduced into computer vision with the aim of effectively finding the salient regions by imitating the human visual system. Recently, a large number of attention mechanisms and approaches for deep learning are proposed to improve the performance of classification. Channel attention adaptively recalibrates the weight of each channel, and determines what to pay attention to. Hu and colleagues \cite{hu2018squeeze} proposed the concept of channel attention and presented SENet for this purpose. The core of SENet is a squeeze-and-excitation (SE) block which is used to collect global information, capture channel-wise relationships and improve representation ability. An SE block can be added after each residual unit due to their low computational resource requirements, which is suitable for our model. Wang and colleagues \cite{2019ECA} proposed the efficient channel attention (ECA) block which uses a 1D convolution to determine the interaction between channels, instead of dimensionality reduction, which can directly model correspondence between weight vectors and inputs. However, it is noticed that the ECA block is not very suitable for models with a small number of channels
in our experiments (as demonstrated in the Table 1 and Figure 2). \cite{ferrari2018computer} proposed the convolutional block attention module (CBAM) which stacks channel attention and spatial attention in series. However, the method of spatial attention which only fits 2D image can not be directly used in mesh. As a result, the SEnet was chose in our work.

\subsection{Predicting Gyro-sulcus From fMRI}
So far, there are very few works in predicting cortical folding pattern from fMRI via a deep learning model. Liu and colleagues \cite{liu2019cerebral} designed a convolutional neural network (CNN) based classifier, which can differentiate gyral and sulcal fMRI signals with a reasonable accuracy (0.67 on resting state fMRI). However, in this model, vertices on surfaces are samples (with a gyrus/sulcus label) independent from each other, whereas functional regions are not spatially isolated and independent. Either functional network studies or functional gradient ones demonstrated that remote or neighboring regions cooperate with each other to elicit and maintain a brain function \cite{huntenburg2018large}\cite{damoiseaux2006consistent}. Therefore, functional connectivity induced gradients that encode the spatial embedding of brain neuronal interactivity could be more suitable for prediction of spatial layout of gyri/sulci. Also, the relatively lower classification accuracy in \cite{liu2019cerebral} could be ascribed to the direct use of noisy fMRI signals as predictive features, whereas connectivities could mitigate such side effects and direct the attention to brain functions at a global scale rather than a vertex-wise one. Zhao and colleagues \cite{zhao2021exploring} proposed a novel Hierarchical Interpretable Autoencoder (HIAE) to explore the hierarchical functional difference between gyri and sulci. The hierarchical features learned by autoencoder can be embedded into a one-dimensional vector which interprets the features as spatial-temporal patterns. However, its major aim is embedding hierarchical and spatial-temporal functional features into a one-dimensional vector, whereas the functional difference between gyri and sulci was only a \emph{post hoc} analysis.

\begin{figure}[htp]
  \centering
  \includegraphics[width=14cm]{figure1.png}
  \caption{Attention mesh convolution network architecture}
\end{figure}
\section{Method}
\label{headings}

In general, we designed a novel attention mesh convolution method based on U-net architecture to predict the gyro-sulcal segmentation map from brain functional gradient maps. The model’s inputs are multi-channel functional gradients on icosahedral mesh \cite{baumgardner1985icosahedral}, while its output is a two-channel binary map on the same icosahedral mesh (Figure 1).

\subsection{Functional Gradients and Cortical Folding Patterns}

In brief, for each subject, a functional connectivity matrix was calculated using the correlation coefficient for all pairs between vertices (32,492 per hemisphere) on HCP grayordinate mesh system, introduced in section 4.1. For each row of this “dense” functional connectome matrix, the values of the top 10\% of connections were retained, whereas all others were zeroed. The negative values in the remaining connections were zeroed as well. The similarity of the connectivity matrix $\mathbf{A}$ was ensured by calculating the cosine distance between all pairs of rows, yielding a symmetric affinity matrix $\mathbf{L}$. On this matrix, diffusion embedding \cite{margulies2016situating} that allows local and long distance connections to be more effectively projected into a common space, was employed to yield unitless components, also known as functional “gradients”. The positions of nodes on these embedding axes encode the differences in nodes’ connectivity patterns. In this work, the number of gradients was set to 5. Results of other numbers of gradients were compared. It is noted that each functional gradient is associated with an eigen value that suggests its dominance. It has been demonstrated that gradients decomposed from different subjects are consistent, whereas their dominance could vary across subjects \cite{vos2020brainspace}. Hence, the algorithm in \cite{vos2020brainspace} used the group-mean matrix to yield the group-wise gradients, that were further used as templates to identify the cross-subject correspondence of gradients from individuals. The top-left panel of Figure 1 illustrates the aligned 5 functional gradient maps from one subject. Surface mean curvature \emph{k} was used to segment surface to gyri (\emph{k}>0) and sulci (\emph{k}<0) (the right-top panel of Figure 1).

\subsection{Channel Attention Block}

In order to boost the representation power of a network, the channel attention block SEnet \cite{hu2018squeeze} was applied in front of mesh convolution. In this attention block, a self-gating mechanism based on channel dependence yields the excitation of each channel (functional gradient in this work), such that channels of input were then re-weighted and fed to the following layers. This weight vector $\mathrm{w}$ could encodes the respective contribution of functional gradients to folding pattern prediction.

As shown in Figure 1 and Equation(1\&2), the $\mathbf{F}$ $(\cdot)$ is global average pooling function which compresses the characteristics of each channel into a real number, where $\mathrm{c}$ denotes the channel number, $N_{f}$ denotes the number of mesh faces, $G_{c}$ denotes the channel's gradient contrasts. Then, a $\mathrm{1*c}$ vector was sent to $\mathbf{F}_{\mathrm{ex}}(\cdot, \mathrm{w})$ which was implemented by a fully connected layer. In our method, we encoded 5 channels to 2 channels and then decoded them back to 5, in which $\mathrm{w}$ is a learnable parameter. After these operations, $\mathrm{w}$ gained the weight for every input channel, and then $\mathbf{F}_{\mathrm{scale}}(^\cdot, ^\cdot)$ scaled the $\mathrm{1*c}$ vector and make the gradient contrasts on every mesh face multiplied by the corresponding channel weight. Finally, the re-weighted channels were fed to the mesh U-Net.
\begin{equation}\label{eq1}
z = \mathbf{F}\left(\mathbf{G}_{c}\right)=\frac{1}{N_f} \sum_{i=1}^{N_f}  G_{c}(i)
\end{equation}
\begin{equation}\label{eq2}
\mathbf{s}=\mathbf{F}_{e x}(\mathbf{z}, \mathbf{W})=\sigma(g(\mathbf{z}, \mathbf{W}))=\sigma\left(\mathbf{W}_{2} \delta\left(\mathbf{W}_{1} \mathbf{z}\right)\right)
\end{equation}

\subsection{Mesh CNN with U-net Architecture and Loss Function}

The model is based on the U-Net architecture \cite{ronneberger2015u} using the mesh convolutional kernel proposed in \cite{hu2022subdivision}. As shown in Figure 1, after the attention block, a 'meshlinear' operation was used to generate 64 channels. This operation is a linear operation on meshes which resembles '1D' convolution on images. Then, the results above were sent to the classic U-Net framework. The kernel size was set to 3 on all mesh convolution operations, which was designed in \cite{hu2022subdivision}, and the 'mean pooling' was used for downsampling. A difference between the classic U-Net framework and ours is that a 'meshlinear' operation was used after the third 'pooling' operation (yellow arrow in Figure 1), and we found this change improves the accuracy of prediction. Finally, the deep network yielded two output channels. The loss was calculated by cross-entropy cost function, and the final result was obtained by an 'argmax' function, in which each mesh face was classified as sulcus or gyrus (Equation(3)).

\begin{equation}\label{eq3}
L=-[y \log \widehat{y}+(1-y) \log (1-\widehat{y})]
\end{equation}

The segmentation accuracy was defined in Equation(4), in which the ${N_f}$ denotes the number of mesh faces, $\mathbf{prd}_{i}$ and $\mathbf{lab}_{i}$ denotes the $i^{th}$ face's prediction and label, respectively. $sign$ function outputs 0 and 1, according to whether $\mathbf{prd}_{i}$ is equal to $\mathbf{lab}_{i}$ or not.

\begin{equation}\label{eq4}
Acc =\frac{1}{N_f} \sum_{i=1}^{N_f}  sign(\mathbf{prd}_{i},\mathbf{lab}_{i})
\end{equation}
\section{Experiments}
\subsection{Dataset}

We used the 3-Tesla T1-weighted MRI and resting-state fMRI (rsfMRI) data from Human Connectome Project (HCP) dataset \cite{van2013wu}. These data have been minimally pre-processed upon release \cite{glasser2013minimal}. Important imaging parameters are: T1-weighted MRI: TR=2400 $ms$, TE=2.14 $ms$, flip angle=8 $deg$, image matrix=260×311×260 and resolution=0.7×0.7×0.7 $mm^3$. Diffusion MRI: TR=5520 $ms$, TE=89.5 $ms$, flip angle=78 $deg$, FOV=210×180 $mm^2$, matrix=168×144, resolution=1.25×1.25×1.25 $mm^3$, echo spacing=0.78 $ms$; rsfMRI: TR=0.72 $s$, TE=33.1 $ms$, flip angle=52 $deg$, in-plane FOV=208×180 $mm^2$, matrix=90×104, 220 $mm$ FOV, 72 slices, 1200 time points, 2.0 $mm$ isotropic voxels, BW = 2290 $Hz/Px$.
T1-weighted MRI was used to reconstruct a cortical surface, since it provides more precise anatomical information due to its higher spatial resolution. Preprocessing of T1-weighted MRI includes skull stripping, tissue-segmentation and white matter surface reconstruction. The surfaces were transferred to the grayordinate system, a common surface template, \emph{via} cross-subject registration to standard volume and surface spaces. It is noted that the aligned surface of each subject was re-sampled such that vertices have cross-subject correspondence. Preprocessing of fMRI signals includes: ICA-denoise, spatial artifact/distortion removal. Within‐subject cross‐modal registration was applied such that each grayordinate vertex was associated with a preprocessed fMRI signal. More details of data acquisition and preprocessing are referred to \cite{glasser2013minimal}. Our experiments included 348 samples for training/testing (0.8/0.2).

\begin{figure}[!htp]
  \centering
  \includegraphics[width=10cm ]{figure2.png}
  \caption{Comparison experiments between ground truth and predicated results \emph{via} different methods. A 'U' shaped gyral area is enlarged for inspection in detail, where yellow circles highlight the unsatisfied predictions. Abbreviations: G for gradients}
\end{figure}

\subsection{Implementation Details}

Firstly, to adapt to the requirement in Subdivnet, we re-meshed the original surface and transferred the gradient features as well as class labels (0 for gyri and 1 for sulci) from mesh vertices to the mesh faces. The batch size and learning rate was set to 6 and 0.02, respectively. 'Bilinear' upsample operation was used. The accuracy on testing samples converges (75.85\%) after about 20 epochs (from 15 to 21). Then, we chose the checkpoint of the best accuracy on test data set and continued fine-tuning for several epochs, with the learning rate set to 0.008.

\subsection{Comparison with Mesh Convolution Baselines and Other Models}

Since our model was modified from Subdivnet \cite{hu2022subdivision}, it was used as a mesh convolution baseline. The BrainSurfCNN \cite{jiang2019spherical} predicts the task fMRI contrasts from resting-state functional connectivity ones by using a regression scheme. For a fair comparison, we modified its' network architecture by replacing the 100 input channels with 5 channels, the 47 output channels with 2 channels in our model, as well as the loss function with cross-entropy loss, to implement the same segmentation task. The features on the mesh vertices were transferred to faces as well. A classic linear classification method, support-vector machine (SVM) \cite{hearst1998support} was also selected for comparison. Figure 2 shows the segmentation results from one random testing subject and Table 1 reports the comparison of results \emph{via} different methods by means of a variety of metrics. It is found that
nonlinear models overwhelmingly outperform linear SVM in all metrics. Besides the intrinsic nonlinear relation between function and folding that could be significantly better encoded and decoded by nonlinear models, the definition of neighboring mask, a key factor for convolution operations on a mesh, could play an important role in the improvement. The importance of the convolution mask was also supported by a better performance of Subdivnet than BrainSurfCNN, since
BrainSurfCNN, unlike Subdivnet, only used a small neighborhood as convolution mask which limits the extraction of long-range features.
It is also found that sulci were often oversegmented by Subdivnet and BrainSurfCNN (yellow circles in Figure 2) in highly bent gyral regions, where the cortical morphology is far more complex than others.

\subsection{Ablation Study}
We considered the effects of different types of attention models and the number of functional gradient channels. It is found in Figure 2 and Table 2 that adding channel attention block or switching to another one does not significantly improve the performance in MIoU, Precision, F1 and Recall, whereas the accuracy could be affected by the selection of a channel attention block. The choice or absence of a channel attention could yield wrong segmentation in "straight" gyral/sulcal regions (the left yellow circles in third row of Figure 2). Neither fewer (4G) nor more (6G) functional gradients improve the segmentation performance. Mis-classification with 4G and 6G was usually found on the border between gyri and sulci (yellow circles in the fourth row of Figure 2).

\begin{table}[!h]
  \caption{Comparison with other methods. The MIoU, Precision, F1-score, Recall and Accuracy are reported.}
  \label{sample-table}
  \centering
  \begin{tabular}{llllll}
    \toprule

    Method     & MIoU     & Precision  & F1-score & Recall & Accuracy \\
    \midrule
    Our Method & 0.603  & 0.752  &0.751 &0.750 &75.85\%   \\
    Subdivnet     & 0.595 & 0.744   &0.745 &0.746	&75.11\%  \\
    BrainSurfCNN     & 0.464	&0.657	&0.622	&0.626	&65.03\%  \\
    Linear SVM     &0.295	&0.488	&0.371	&0.499	&59.09\% \\
    \bottomrule
  \end{tabular}
\end{table}

\begin{table}[!h]
  \caption{Ablation study on attention block and gradients. The MIoU, Precision, F1-score, Recall and Accuracy are reported. \textcolor{red}{Red} and \textcolor{blue}{blue} denote the best and the second-best results, respectively. Abbreviations: G for gradients.}
  \label{sample-table}
  \centering
  \begin{tabular}{lllllll}
    \toprule

    Type     & MIoU     & Precision  & F1-score & Recall & Accuracy \\
    \midrule
    SEnet+5Gs     &\textcolor{red}{0.603}	&\textcolor{red}{0.752}	&\textcolor{red}{0.751}	&\textcolor{red}{0.750}	&\textcolor{red}{75.85\%} 	 \\
    ECAnet+5Gs     &\textcolor{blue}{0.599} 	&\textcolor{blue}{0.750}	&\textcolor{blue}{0.747}	&0.745	&\textcolor{blue}{75.61\%}  \\
    No Attention+5Gs &0.598 &0.748 &0.747 &\textcolor{blue}{0.746} &75.45\% & \\
    SEnet+4Gs     &0.587  &0.744    &0.737 &0.734	&74.89\% \\
    SEnet+6Gs     &0.580 	&0.732	&0.732	&0.732	&73.86\%  \\
    \bottomrule
  \end{tabular}
\end{table}
\section{Discussion}
\begin{figure}[!htp]
  \centering
  \includegraphics[width=14cm]{figure3.png}
  \caption{(a) Activation maps in different deep layers that resemble the fronto-parietal functional network, not detected by functional gradients. (b) Similarity between functional gradients (G for short) and curvature map, and the weights of Gs learnt from attention block. (c) The segmentation activation maps (red for more activated regions) in the last two-channel layer and cortical landmarks (cyan bubbles), including gyral peaks and sulcal pits, that are independent from curvature.}
\end{figure}

The attention block re-weights the five gradients. The latter gradients (3G and 5G) gain relatively higher weights (Figure 3(b)), and their similarity to curvature map (measured by Pearson correlation) is significantly lower than the dominant gradients (1G and 2G). These results suggest that the implicit relation between function and folding might not be encoded in dominant gradients. It is also interesting to find that some functional networks, which are not explicitly represented by input gradients, can be identified in activation maps in deep layers. As an example, the fronto-parietal network, one of the resting-state networks \cite{smith2013functional} not present in the 5 gradients, were found in both encoding and decoding layers, as shown in Figure 3(a). Finally, the activation maps in the last two-channel layer, were related to folding landmarks, gyral peaks (local convexity maxima) and sulcal pits (local concavity maxima) \cite{li2022cortical}, which were obtained independent of curvature used to yield the output gyro-sulcal maps. These landmarks (blue bubbles in Figure 3(c), group-wise mean locations across subjects) are not located on the highly activated regions (red regions). On average, activation values on these landmarks are 0.46±0.62 and 0.60±0.55 (\emph{p}<0.05, \emph{via} unpaired \emph{t}-test) on other regions (expect the gray background), suggesting that these landmarks might mark the borders of folding patterns that are highly correlated with brain functions.

\section{Conclusion}
In this work, we developed an attention mesh convolution model based on the U-net architecture to predict cortical gyro-sulcal segmentation maps from brain functional gradients derived from resting-state fMRI. The model architecture, such as input channel, output channel, feature types and attention block, can be easily extended or customized according to specific needs. Unlike the previous study, where functional signals from gyral/sulcal vertices were independently used as samples to a CNN model and their frequency band was found to be a discriminative feature \cite{liu2019cerebral}, the mesh convolution in our model considers the spatial organization of functional gradients and folding patterns on a cortical sheet. The prediction performance \emph{via} our model outperforms the state-of-the-art. In addition, we found that the dominant functional gradients contribute less to folding prediction. On the activation maps of the last layer, we found that some well-studied cortical landmarks are on the borders of, rather than within, the highly activated regions. These results and findings suggest that a specifically designed artificial neural network can improve the precision of the mapping between brain functions and cortical folding patterns, and can provide valuable clues of brain anatomy-function relation for neuroscience and artificial neural network design.

{
\small
\bibliographystyle{unsrt}