11institutetext: Rongjiao Ji, Alessandra Micheletti 22institutetext:
Universitá degli Studi di Milano, 22email<EMAIL_ADDRESS><EMAIL_ADDRESS>33institutetext: Natasa Krklec Jerinkic
44institutetext: University of Novi Sad, 44email<EMAIL_ADDRESS>55institutetext: Zoranka Desnica 66institutetext: 3Lateral DOO, 66email:
<EMAIL_ADDRESS>
# Group pattern detection of longitudinal data using functional statistics

Rongjiao Ji  Alessandra Micheletti  Natasa Krklec Jerinkic  Zoranka Desnica

###### Abstract

Estimations and evaluations of the main patterns of time series data in groups
benefit large amounts of applications in various fields. Different from the
classical auto-correlation time series analysis and the modern neural networks
techniques, in this paper we propose a combination of functional analysis of
variance (FANOVA) and permutation tests in a more intuitive manner for a
limited sample size. First, FANOVA is applied in order to separate the common
information and to dig out the additional categorical influence through paired
group comparison, the results of which are secondly analyzed through
permutation tests to identify the time zones where the means of the different
groups differ significantly. Normalized kernel functions of different groups
are able to reflect remarkable mean characteristics in grouped unities, also
meaningful for deeper interpretation and group-wise classification. In order
to learn whether and when the proposed method of FANOVA and permutation F-test
works precisely and efficiently, we compare the estimated kernel results with
the ground truth on simulated data. After the confirmation of the model’s
efficiency from simulation, we apply it also to the RAVDESS facial dataset to
extract the emotional behaviors of humans based on facial muscles contractions
(so-called action units (AU) technically in computer graphics), by comparing
the neutral performances with emotional ones.

###### Keywords:

functional analysis of variance (FANOVA), permutation F-test, emotion,
expression evolution, action units

##  1 Introduction

Many natural, economical, or technological phenomena can be represented as
real functions $f:\mathbb{R}\to\mathbb{R}$ where the domain can be time or
another one-dimensional variable, and data of many copies of these functions
can be collected, usually in a discrete format.

We here consider the case where we collect data of a set of random functions
describing a real phenomenon, which is characterized by the presence of
groups, where each group is characterized by a different mean underlying
pattern that we want to detect.

Examples of such situations are: the collection of temperature evolution in
different regions of the world, in the same time window, where we want to
study the underlying temperature evolution functions characterizing each
region and their main differences; the collection of stock price evolution,
with stocks related with goods of different economical sectors, or different
segments of the same sector, and we want to extract and compare the underlying
evolution patterns of the sector or segments; etc.

In this paper we will face the problem of underlying pattern detection and
comparison of groups of random functions, using a specific driving example:
the detection of human emotions from functions describing the movements of
facial muscles during a video.

The study of human facial expressions and emotions never stops in our daily
life while we communicate with others. Following the increased interest in
automatic facial behavior analysis and understanding, the need of a semantic
interpretation of the evolution of facial expressions and of human emotions
has become of interest in recent years Fridlun2014 . In this paper, based on a
work cooperated with the Serbian company 3Lateral, which has special expertise
on building visual styles and designs in animation movies, we want to explore
functional statistical instruments to identify the emotions while analyzing
the expressions through recorded videos of human faces. The final aim of this
research is to use this information to better and more realistically establish
virtual digital characters, able to interact autonomously with real humans.

The data that we consider are multivariate longitudinal data, showing the
evolution in time of different face muscles contraction. Functional Data
Analysis (FDA) offers the possibility to analyze the entire expression
evolution process over time and to gain detailed and in-depth insight into the
analysis of emotion patterns. The basic idea in functional data analysis is
that the measured data are noisy observations coming from a smooth function.
Ramsay and Silverman ramsay1997functional  describe the main features of FDA,
that can be used to perform exploratory, confirmatory or predictive data
analysis. Ullah and Finch UllahShahid2013  published a systematic review on
the applications of functional data analysis, where they included all areas
where FDA was applied.

Here we will apply models of functional regression and functional analysis of
variance (FANOVA) to identify the mean functional patterns related with
different emotions. Additionally we setup a technique, based on permutation
tests, to identify the time zones where the identified patterns are
significantly different. The identified patterns can then be used to perform
an automatic detection and classification of the emotions expressed by a human
face.

The methods here proposed are also tested on simulated data to highlight their
effectiveness in identifying the right underlying patterns and the correct
time zones in which the patterns differ, also in relationship with the level
of noise present in the collected data.

In spite of the specific motivating case study here considered, the methods
here proposed and the related software can be applied to many other similar
situations in different fields of application.

##  2 Methods

###  2.1 Functional Data Analysis

Functional Data Setting Let $\Omega$ be a population space of statistical
units ($\omega\in\Omega$) and
$s_{n}=\\{\omega_{1},\omega_{2},\dots,\omega_{n}\\}$ be a random sample of
size $n$ drawn from $\Omega$. Let $\mathbf{Y}$ be a $L_{2}$-continuous
multivariate stochastic process defined on $\Omega$ associating to a set of
$D$ curves, $D\geq 2$, where each curve in one dimension defined on a given
finite interval $\mathcal{T}=[0,T],0<T<\infty$, i.e.,
$\mathbf{Y}(\omega)=[\mathbf{y}_{\omega,1}(t),\dots,\mathbf{y}_{\omega,D}(t)]^{T}$
(denoted as $\mathbf{Y}=(\mathbf{y}_{1},\dots,\mathbf{y}_{D})^{T}$ for
simplicity), $t\in\mathcal{T}$, such that

$\lim\limits_{h\to
0}\mathbb{E}[||\mathbf{Y}(t+h)-\mathbf{Y}(t)||^{2}]=\lim\limits_{h\to
0}\int_{0}^{T}\sum\limits_{d=1}^{D}\mathbb{E}[(\mathbf{y}_{d}(t+h)-\mathbf{y}_{d}(t))^{2}]dt=0.$

The observations of $\mathbf{Y}$ on the sample $s_{n}$ are called multivariate
functional data, providing independent realizations as samples
$\\{\mathbf{Y}(\omega_{1}),\dots,\mathbf{Y}(\omega_{n})\\}$ with multivariate
curves jacques2012clustering . The mean functions of $\mathbf{y}_{d}$ and
$\mathbf{Y}$ are

$\mathbb{E}[\mathbf{y}_{d}(t)]=\mu_{d}(t),t\in[0,T],\text{ and
}\mathbb{E[\mathbf{Y}]}=\mathbf{\mu}(t)=(\mu_{1}(t),\dots,\mu_{d}(t),\dots,\mu_{D}(t))^{T}.$

The covariance operator of $\mathbf{Y}$ is then defined as an integral
operator $\mathcal{C}$ with kernel

$\mathbf{C}(t,s)=\mathbb{E}[(\mathbf{Y}(t)-\mathbf{\mu}(t))\otimes(\mathbf{Y}(s)-\mathbf{\mu}(s))],\forall
t,s\in[0,T],$

where $\otimes$ is the tensor product on $\mathbb{R}^{D}$. Thus,
$\mathbf{C}(t,s)$ is a $D\times D$ matrix with elements

$\mathbf{C}(t,s)[i,j]=Cov(\mathbf{Y}_{i}(t),\mathbf{Y}_{j}(s)),i,j=1,\dots,D.$

Functional Principal Component Analysis (FPCA) Under the hypothesis of
$L_{2}$-continuity, $\mathcal{C}$ is an Hilbert-Schmidt operator, i.e., it is
compact and self-adjoint. Meanwhile, the spectral analysis of $C$ provides a
countable set of positive eigenvalues $\\{\eta_{j,j\geq 1}\\}$ associated to
an orthonormal basis of multivariate eigen-functions $\\{\mathbf{f}_{j,j\geq
1}\\}$, $\mathbf{f}_{j}=[\mathbf{f}_{j,1},\dots,\mathbf{f}_{j,D}]^{T}$, such
that

$C\mathbf{f}_{j}=\eta_{j}\mathbf{f}_{j},$

with $\eta_{1}\geq\eta_{2}\geq\dots$ and
$\langle\mathbf{f}_{i},\mathbf{f}_{j}\rangle_{\\{L_{2}([0,T])\\}^{D}}=\delta_{i,j}$
with $\delta_{i,j}=1$ if $i=j$, otherwise $\delta_{i,j}=0$
berrendero2011principal . Assume that the principal components $\eta_{j}$ of
$\mathbf{Y}$ are zero-mean random variables, as the projections of
$\mathbf{Y}$ on the $j$th eigenfunction of $\mathbf{C}$:

$\eta_{j}=\int_{0}^{T}\langle\mathbf{Y}(t)-\mathbf{\mu}(t),\mathbf{f}_{j}(t)\rangle
dt=\int_{0}^{T}\sum\limits_{d=1}^{D}(Y_{d}(t)-\mu_{d}(t))\mathbf{f}_{j,d}(t)dt,$
(1)

which means $\mathbf{Y}(t)$ varies around the mean function with random
amplitude variations in the directions of $\\{\mathbf{f}_{j}\\}_{j\geq 1}$.
Therefore, following Karhunen-Loeve expansion, it holds that

$\mathbf{Y}(t)=\mathbf{\mu}(t)+\sum\limits_{j\geq
1}\eta_{j}\mathbf{f}_{j}(t)\text{ and
}\mathbf{y}_{d}(t)=\mu_{d}(t)+\sum\limits_{j\geq
1}\eta_{j}\mathbf{f}_{j,d}(t)$ (2)

are the best approximations of this form under the mean square criterion.

Functional Registration In practice, it is essential to align observed
functions into a common registered internal timeline to isolate the phase
variability of the functions, but keeping, at the same time, the amplitude-
phase unchanged to maintain the information under interests. The phase
variation is normally represented by a random change of time scale, which is
mostly a non-linear transformation. We use the warping functions
$T_{i}:[0,T]\rightarrow[0,T],i=1,\dots,n$, assuming that they are increasing
functions independent of amplitude variation. They map unregistered
chronological time $t_{i}^{*}$ to registered internal time $t$ so that
$T_{i}^{-1}(t_{i}^{*})=t$, with $E[T_{i}(t)]=t$. The observed time-warped
curves, represented through a Karhunen-Loeve expansion based on a functional
basis $\mathbf{f}_{j,d}$, are

$\tilde{y}_{i,d}(t)=y_{i,d}(T_{i}^{-1}(t_{i}^{*}))=\mu_{d}(T_{i}^{-1}(t_{i}^{*}))+\sum\limits_{j\geq
1}\eta_{j}\mathbf{f}_{j,d}(T_{i}^{-1}(t_{i}^{*})).$

The literature is rich of proper methods to do this, and a more in-depth
review would be beyond the scope of this report. Briefly speaking, we followed
the principal components based registration method wrobel2019register  with a
generative process, by alternating the process between using the generalized
FPCA to calculate template functions and estimating smooth time warping
functions that map observed curves to templates 111Codes are available in the
R package ”registr” wrobel2018register .

Approximation into finite basis of functions Often in practice, data are
observed at discrete time points and with some noise. In order to get the
functional feature of data, researchers normally use some smoothing and
interpolation methods to rebuild the consecutive function from the discrete
data, while assuming that the latent function belongs to a finite dimensional
functional space spanned by some functional basis. In this way, functional
data objects are constructed from time series data by specifying a set of
basis functions and a set of coefficients defining a linear combination of
these basis functions. For non-periodic time series, B-spline basis functions
which contain polynomial segments jointed end-to-end at argument values are
frequently considered, benefiting the most from fast computational ability and
great flexibility ramsay1997functional .

###  2.2 Multiple Multivariate Function-on-Scalar Regression

Let $y_{g,d,k}(t)$ denote the evolution of one specific variate
$d\in\\{1,\dots,D\\}$ in sample $k\in\\{1,\dots,K\\}$ for one group
$g\in\\{0,\dots,G\\}$. We can assume that

$y_{g,d,k}(t)=\mu_{d,0}(t)+\alpha_{d,g}(t)+\epsilon_{g,d,k}(t),$ (3)

where $\mu_{d,0}(t)$ is the grand mean function independent from groups. The
term $\alpha_{d,g}(t)$ is the additional effect function of group $g$ on the
considered variate $d$, while $\epsilon_{g,d,k}(t)$ represents the unexplained
zero-mean variation, specific of the $k$-th sample within group $g$. In order
to identify the parameter functions uniquely, we require the constraint
$\sum\limits_{g=0}^{G}\alpha_{d,g}(t)=0,\forall t.$

By grouping the samples from the same group, we can define a
$(GK+1)\times(G+1)$ design matrix $\mathbf{Z}_{d}$ for univariate model, with
suitable 0 and 1 entries, as described in (ramsay1997functional, , Section
9.2), and rewrite Equation (3) in matrix form

$\mathbf{y}_{d}(t)=\mathbf{Z}_{d}\mathbf{\beta}_{d}(t)+\mathbf{\epsilon}_{d}(t),\forall
t,$

where

$\mathbf{y}_{d}(t)=\begin{bmatrix}y_{0,d,1}(t)\\\ \vdots\\\ y_{0,d,K}(t)\\\
y_{1,d,1}(t)\\\ \vdots\\\ y_{1,d,K}(t)\\\ \vdots\\\ y_{G,d,1}(t)\\\ \vdots\\\
y_{G,d,K}(t)\\\
0\end{bmatrix},\mathbf{Z}_{d}=\begin{bmatrix}1&1&0&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\cdots&\vdots&\vdots\\\ 1&1&0&0&\cdots&0&0\\\
1&0&1&0&\cdots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\cdots&\vdots&\vdots\\\
1&0&1&0&\cdots&0&0\\\ &&&&\vdots&&\\\ 1&0&0&0&\cdots&0&1\\\
\vdots&\vdots&\vdots&\vdots&\cdots&\vdots&\vdots\\\ 1&0&0&0&\cdots&0&1\\\
0&1&1&1&\cdots&1&1\end{bmatrix},\mathbf{\beta}_{d}(t)=\begin{bmatrix}\mu_{d,0}(t)\\\
\alpha_{d,0}(t)\\\ \alpha_{d,1}(t)\\\ \vdots\\\
\alpha_{d,G}(t)\end{bmatrix}=\begin{bmatrix}\beta_{d,0}(t)\\\
\beta_{d,1}(t)\\\ \beta_{d,2}(t)\\\ \vdots\\\
\beta_{d,G+1}(t)\end{bmatrix},\mathbf{\epsilon}_{d}(t)=\begin{bmatrix}\epsilon_{0,d,1}(t)\\\
\vdots\\\ \epsilon_{0,d,K}(t)\\\ \epsilon_{1,d,1}(t)\\\ \vdots\\\
\epsilon_{1,d,K}(t)\\\ \vdots\\\ \epsilon_{G,d,1}(t)\\\ \vdots\\\
\epsilon_{G,d,K}(t)\\\ 0\end{bmatrix}.$

Here $\mathbf{\beta}_{d}$ is defined as the corresponding set of $G+1$
regression functions for mean and $G$ additional groups, and
$\mathbf{\epsilon}_{d}$ is the matrix of noise with the same dimension of
$\mathbf{y}_{d}$.

Following the notations above, we can write the multivariate matrix form for
the multiple regression functions, as

$\mathbf{Y}(t)=\mathbf{Z}\mathbf{\beta}(t)+\mathbf{\epsilon}(t),$

where

$\mathbf{Y}(t)=\begin{bmatrix}\mathbf{y}_{1}(t)\\\ \vdots\\\
\mathbf{y}_{d}(t)\\\ \vdots\\\
\mathbf{y}_{D}(t)\end{bmatrix},\mathbf{Z}=\begin{bmatrix}\mathbf{Z}_{1}\\\
\vdots\\\ \mathbf{Z}_{d}\\\ \vdots\\\
\mathbf{Z}_{D}\end{bmatrix}^{T},\mathbf{\beta}(t)=\begin{bmatrix}\mathbf{\beta}_{1}(t)\\\
\vdots\\\ \mathbf{\beta}_{d}(t)\\\ \vdots\\\
\mathbf{\beta}_{D}(t)\end{bmatrix},\mathbf{\epsilon}(t)=\begin{bmatrix}\mathbf{\epsilon}_{1}(t)\\\
\vdots\\\ \mathbf{\epsilon}_{d}(t)\\\ \vdots\\\
\mathbf{\epsilon}_{D}(t)\end{bmatrix}.$

Under the functional least squares fitting criterion on $R^{GDK\times T}$
space, we can estimate the parameters

$\hat{\mathbf{\beta}}(t)=\hat{\mathbf{B}}\Theta(t),$ (4)

where $\Theta(t)=(\theta_{1}(t),\dots,\theta_{q}(t),\dots,\theta_{Q}(t))^{T}$
denotes a linear combination of basis functions with $J_{\theta\theta}$ as the
inner product of each pair of functional basis, and $\hat{\mathbf{B}}$ is the
estimated coefficient matrix with

$vec(\mathbf{B}^{T})=[(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})^{T}(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})]^{-1}(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})^{T}vec(J_{\theta\theta}^{\frac{1}{2}}\mathbf{A}^{T}).$

Here, $\mathbf{A}$ is the coefficient matrix of samples $\mathbf{Y}(t)$ with a
fixed B-spline functional basis and more computation details in matrix format
are shown in Appendix 3.1.

###  2.3 Functional ANOVA and F-type tests

After establishing the model for functional evolution processes and estimating
the additional parameter functions affected by group categories, next goal is
to investigate whether and when group labels are significantly influencing the
change of variables. ANOVA-related methods are suitable tools to determine the
effects of independent categorical variables on continuous dependent
variables, by comparing inter-group differences between centroids represented
as the vectors of the mean values of the dependent variables. Moreover, in
order to construct feasibly more general hypothesis tests for group means,
contrasts, defined as a linear combination of regression coefficients where
the values of the linear combination sum up to 0, are applied by specifying
the values on different pairs of groups gradually.

For each subgroup of variable $d\in\\{1,\dots,D\\}$ and group
$g\in\\{0,\dots,G\\}$, assume that the $k$-th sample curve in this subgroup
$y_{d,g,k}(t)\sim\mathbb{N}(\mu_{d,g}(t),\sigma^{2}),k\in\\{1,\dots,K\\}$ with
group mean $\mu_{d,g}(t)$ and constant variance $\sigma^{2}$, and therefore
the subgroup sample mean
$\overline{Y}_{d,g,\cdot}(t)=\frac{1}{K}\sum\limits_{k=1}^{K}y_{d,g,k}(t)\sim\mathbb{N}(\mu_{d,g}(t),\frac{\sigma^{2}}{K})$.
Here denote $\mathbf{\mu}(t)$ as the sequence of subgroup means $\mu_{d,g}(t)$
for all $d\in\\{1,\dots,D\\}$ and $g\in\\{0,\dots,G\\}$. We construct one test
for each pair of groups in comparison (one control group $g=0$ and another
group under study $g=\tilde{g}$ for one variable $d$) and in $i$-th test
($i\in\\{1,\dots,DG\\}$) we design a vector
$\mathbf{a}=[0,\dots,0,1,0,\dots,0,-1,0,\dots,0]$ in length $D(G+2)$ (same as
the length of $\mathbf{\beta}(t)$ in last section) with the entries pointed to
the pair of groups in comparison being 1 and -1 respectively, i.e.,
$a_{d(G+2)+2}=1$, $a_{d(G+2)+2+\tilde{g}}=-1$, and 0 in entries remained. We
define the contrast in $i$-th test as
$c=\mathbf{a}^{T}\mathbf{\mu}(t)=\alpha_{d,0}(t)-\alpha_{d,\tilde{g}}(t)=(\alpha_{d,0}(t)+\mu_{d,0})-(\alpha_{d,\tilde{g}}(t)+\mu_{d,0})$
for $DG$ different group comparisons, where the sum of weigths in $\mathbf{a}$
is zero, and then the difference of subgroup sample means
$\overline{Y}_{d,0,\cdot}(t)-\overline{Y}_{d,\tilde{g},\cdot}(t)$ is an
unbiased estimator of the contrast $c$. Following this setting, the null
hypothesis at each time point $t$ for the control group and the group under
study $\tilde{g}$ in variable $d$ is

$\mathcal{H}_{0}:\mathbf{a}^{T}\beta(t)=0\Leftrightarrow\alpha_{d,0}(t)-\alpha_{d,\tilde{g}}(t)=0,$

and the related F statistics under null hypothesis is

$F_{d,\tilde{g}}(t)=\frac{(\overline{Y}_{d,0,\cdot}(t)-\overline{Y}_{d,\tilde{g},\cdot}(t))^{2}}{S_{p}^{2}\frac{2}{K}},$

where the sum of sample variances is
$S_{p}^{2}=\frac{1}{N-DG}\sum\limits_{g=1}^{G}\sum\limits_{d=1}^{D}\sum\limits_{k=1}^{K}(y_{d,k,g}(t)-\overline{y}_{d,g,\cdot}(t))^{2}$
with sample size $N$. When assuming the functional data $y_{d,k,g}(t)$ as
realizations of Gaussian processes, for any given $t\in\mathcal{T}$,
$F_{d,\tilde{g}}(t)\sim F_{1,N-DG}$, where $F_{1,N-DG}$ denotes
$F$-distribution with $1$ and $N-DG$ degrees of freedom.

The statistic $F_{d,\tilde{g}}(t)$ should be compared with a critical value of
$F_{1,N-DG}$ for any predetermined significant level $\alpha$. The critical
value of $F_{1,N-DG}(\alpha)$ represents the $\alpha$-fraction of rejection
area, where $F_{1,N-DG}(\alpha)$ denotes the upper $100\alpha$ percentile of
$F_{1,N-DG}$. We write the $1-\alpha$ quantile of the distribution of F
statistic by $c_{1-\alpha}$, and given $c_{1-\alpha}$, we reject the null
hypothesis when $F_{d,\tilde{g}}(t)>c_{1-\alpha}$, which indicates the
dissimilarity of group means under comparison. First, we can test the
statistic at all points of $\mathcal{T}$ using the same critical value from
standard F-Distribution table with correlated degrees of freedom. Since
complex assumptions and convoluted formula are not often satisfied in practice
for most classic statistical tests, we consider also a nonparametric test,
permutation test, which is light on assumptions, widely applicable and
intuitive friendly. The permutation distribution is obtained by shuffling the
labels within the paired groups in comparison and calculate the F statistic
$F_{d,\tilde{g}}(t)$ many times to approximate to the statistic’s
distribution. More inference of permutation test can be seen in Appendix 4.

###  2.4 Kernel-based Classification

Figure 1: Post hoc analysis by classifying from a given video

After the explanation of functional linear regressions and test methods, we
are able to obtain the ground mean functions and the additional group-wise
functions revealing the latent group effects on variable curves. Most
importantly, we can detect the time zones where the additional effects impact
significantly in statistical sense. To further examine the efficiency of the
information, a natural attempt would be to apply a post hoc analysis: do
classifications using the knowledge obtained and check the correctness. Before
the classification, four steps are needed in our framework to project the
input functions into normalized scores in lower dimensions:

  * •

First, centralize the input evolution by subtracting the estimated mean
function for each sample $k$ on variable $d$.

$\tilde{y}_{d,k}(t)=y_{d,k}(t)-\mu_{d,0}(t)$

  * •

Secondly, use the centralized evolution to make inner product with $G$
estimated addtional group-wise functions (also called kernels for simplicity)
in turn, and thus the evolution function is projected to a vector space with
$G$ dimensions. Let $S_{g,d,k}$ represent the score of group $g$ on variable
$d$ in sample $k$:

$S_{g,d,k}=(\tilde{y}_{d,k}(t),\hat{\alpha}_{g,d}(t))=\int_{0}^{T}w_{g,d}\tilde{y}_{d,k}(t)\hat{\alpha}_{g,d}(t)dt=\int_{I}w_{g,d}\tilde{y}_{d,k}(t)\hat{\alpha}_{g,d}(t)dt,$

where $I$ represents the union of relevant intervals, and $w_{g,d}$ is the
weights in the integral process, constant in time.

  * •

Furthermore, normalize the scores across all groups and variables from one
sample $k$, so that $\tilde{S}_{g,d,k}\in[0,1]$, in order to keep the relative
rate/percentage across groups.

    * –

If $S_{g,d,k}<0$ for some $g$, move them to positive numbers by subtracting
the minimum ($min\\{0,min_{g}S_{g,d,k}\\}$).

    * –

Then, rescale them to be smaller than $1$ across groupss.

$\tilde{S}_{g,d,k}=\frac{S_{g,d,k}-min\\{0,min_{g}S_{g,d,k}\\}}{\sum\limits_{g=1}^{G}[S_{g,d,k}-min\\{0,min_{g}S_{g,d,k}\\}]}.$

In this way $\tilde{S}_{g,d,k}$ represents the probability that variable $d$
in sample $k$ behaves as in (has been generated by) group $g$. If the sample
$k$ that we are considering is labeled with $g$, more importance is given to
the variables that assign a highest probability $\tilde{S}_{g,d,k}$ to the
group $g$.

  * •

Finally, combine the scores across variables for sample $k$ (labeled with $g$)
by a linear convex combination with the trained weights $r_{l,g}$, such that,
if the label of sample $k$ is $g$, then $S_{g,k}$ is the maximum score.

$S_{g,k}^{comb}=r_{1,g}\tilde{S}_{g,1,k}+r_{2,g}\tilde{S}_{g,2,k}+\dots+r_{L,g}\tilde{S}_{g,L,k},$

with $r_{l,g}\in[0,1]$ and $\sum\limits_{l=1}^{L}r_{l,g}=1$. For sample $k$,
$[S_{1,k},S_{2,k},\dots,S_{8,k}]$ contains the scores of all groups under
study, and all the elements of this vector fall in region $[0,1]$.

For classification based on scores, we apply SVM method on scores
$[S_{1,k},S_{2,k},\dots,S_{G,k}]$ to classify the samples and examine the
kernels efficiency. A further study on the classification methods is important
but is skippped here for space saving because the classification is not the
focus of this paper, although we will explore further in next study.

##  3 Details of Methods

###  3.1 Approximation into a finite basis of functions

Assume that each curve $y_{g,d,k}(d\in\\{1,\dots,D\\}$ for $k$-th video can be
expressed as a linear combination of basis functions
$\Theta(t)=(\theta_{1}(t),\dots,\theta_{q}(t),\dots,\theta_{Q}(t))^{T}$. We
define $J_{\theta\theta}$ to represent the inner product of each pair of
functional basis.

$y_{g,d,k}(t)=\sum\limits_{q=1}^{Q}a_{g,d,k,q}\theta_{q}(t)=\mathbf{a}_{g,d,k}\Theta(t),$

$\mathbf{y}_{g,d}(t)=\mathbf{A}_{g,d}\Theta(t),\ \
\mathbf{A}_{g,d}=[\mathbf{a}_{g,d,1},\dots,\mathbf{a}_{g,d,K}]^{T}$

$\mathbf{Y}_{g}(t)=\mathbf{A}_{g}\Theta(t),\ \
\mathbf{A}_{g}=[\mathbf{A}_{g,1},\dots,\mathbf{A}_{g,d},\dots,\mathbf{A}_{g,D}]^{T}$

$\mathbf{Y}=\mathbf{A}\Theta(t),\ \
\mathbf{A}=[\mathbf{A}_{1},\dots,\mathbf{A}_{g},\dots,\mathbf{A}_{G}]^{T}$

$\text{with
}\mathbf{A}_{g}=[\mathbf{A}_{g,1},\dots,\mathbf{A}_{g,d},\dots,\mathbf{A}_{g,D}]^{T},\mathbf{A}_{g,d}=\begin{bmatrix}a_{g,d,1,1}&\dots&a_{g,d,1,Q},\\\
\vdots&\vdots&\vdots\\\ a_{g,d,K,1}&\dots&a_{g,d,K,Q}\end{bmatrix}.$

Similarly, the coefficient function
$\mathbf{\beta}=[\beta_{1}(t),\dots,\beta_{D}(t)]^{T}$ can be decomposed with
the same functional basis in the following way

$\mathbf{\beta}=\mathbf{B}\Theta(t),$ $\text{with
}\mathbf{B}=[\mathbf{B}_{1},\dots,\mathbf{B}_{d},\dots,\mathbf{B}_{D}]^{T},\mathbf{B}_{d}=\begin{bmatrix}b_{1,d,1}&\dots&b_{1,d,Q},\\\
\vdots&\vdots&\vdots\\\ b_{G,d,1}&\dots&b_{G,d,Q}\end{bmatrix}.$

Consider that in multivariate linear regression setting

$\mathbf{Y}=\mathbf{Z}\mathbf{\beta}+\mathbf{E},$

by substituting the functional curves with coefficient matrix and functional
basis, we have

$\mathbf{A}\Theta(t)=\mathbf{Z}\mathbf{B}\Theta(t)+\mathbf{E}$

Notice that $\mathbf{B}$ is the coefficient that we want to estimate. Here by
applying least square method on $R^{GDK\times T}$ space, we have

$\displaystyle J(B)$
$\displaystyle=\|E\|_{L_{2}}^{2}=\sum\limits_{\begin{subarray}{c}g\in\\{1,\dots,G\\},\\\
d\in\\{1,\dots,D\\},\\\
k\in\\{1,\dots,K\\}\end{subarray}}[y_{g,d,k}(t)-\hat{y}_{g,d,k}(t)]_{L_{2}}^{2}$
(5)
$\displaystyle=\int\|\mathbf{A}\Theta(t)-\mathbf{Z}\mathbf{B}\Theta(t)\|_{F}^{2}dt$
(6)
$\displaystyle=\int\|(\mathbf{A}-\mathbf{Z}\mathbf{B})\Theta(t)\|_{F}^{2}dt$
(7) $\displaystyle=\int
trace[((\mathbf{A}-\mathbf{Z}\mathbf{B})\Theta(t))^{T}((\mathbf{A}-\mathbf{Z}\mathbf{B})\Theta(t))]dt$
(8) $\displaystyle=\int
trace[\Theta(t)^{T}(\mathbf{A}-\mathbf{Z}\mathbf{B})^{T}(\mathbf{A}-\mathbf{Z}\mathbf{B})\Theta(t)]dt$
(9)
$\displaystyle=trace[\int\Theta(t)^{T}(\mathbf{A}-\mathbf{Z}\mathbf{B})^{T}(\mathbf{A}-\mathbf{Z}\mathbf{B})\Theta(t)dt]$
(10)
$\displaystyle=trace[(\mathbf{A}-\mathbf{Z}\mathbf{B})^{T}(\mathbf{A}-\mathbf{Z}\mathbf{B})\int\Theta(t)^{T}\Theta(t)dt]$
(11)
$\displaystyle=tr[((\mathbf{A}-\mathbf{Z}\mathbf{B})J_{\theta\theta}^{\frac{1}{2}})^{T}(\mathbf{A}-\mathbf{Z}\mathbf{B})J_{\theta\theta}^{\frac{1}{2}}$
(12)
$\displaystyle=vec[(\mathbf{A}-\mathbf{Z}\mathbf{B})J_{\theta\theta}^{\frac{1}{2}}]^{T}vec[(\mathbf{A}-\mathbf{Z}\mathbf{B})J_{\theta\theta}^{\frac{1}{2}}]$
(13)
$\displaystyle=\|vec[(\mathbf{A}-\mathbf{Z}\mathbf{B})J_{\theta\theta}^{\frac{1}{2}}]^{T}\|_{2}^{2}$
(14)
$\displaystyle=\|vec(J_{\theta\theta}^{\frac{1}{2}}\mathbf{A}^{T})-vec[J_{\theta\theta}^{\frac{1}{2}}\mathbf{B}^{T}\mathbf{Z}^{T})\|^{2}$
(15)
$\displaystyle=\|vec(J_{\theta\theta}^{\frac{1}{2}}\mathbf{A}^{T})-(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})vec(\mathbf{B}^{T})\|^{2},$ (16)

Let $J(\mathbf{B})=\mathbf{0}$, we can derive that

$vec(\mathbf{B}^{T})=[(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})^{T}(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})]^{-1}(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})^{T}vec(J_{\theta\theta}^{\frac{1}{2}}\mathbf{A}^{T}),$

together with

$\displaystyle\mathbf{C}_{mat}$ $\displaystyle=(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})^{T}(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})$
$\displaystyle=(\mathbf{Z}^{T}\mathbf{Z})\otimes({J_{\theta\theta}^{\frac{1}{2}}}^{T}J_{\theta\theta}^{\frac{1}{2}})$
$\displaystyle\mathbf{D}_{mat}$ $\displaystyle=(\mathbf{Z}\otimes
J_{\theta\theta}^{\frac{1}{2}})^{T}vec(J_{\theta\theta}^{\frac{1}{2}}\mathbf{A}^{T})$

##  4 Permutation Test

Let’s denote $\Phi_{n}$ as the set of all possible permutations of
${1,\dots,n}$, $M_{n}$ as the cardinality of $\Phi_{n}$. A permutation
$\phi\overset{\Delta}{=}(\phi_{1},\dots,\phi_{n})\in\Phi$.

$X_{n}^{\phi}$: the permuted version of $X_{n}$,
$X_{n}^{\phi}\overset{\Delta}{=}(X_{1}^{\phi},\dots,X_{n}^{\phi})$.
$Fr_{n}^{\phi}=Fr_{n}(X_{n}^{\phi})$ to denote the F-test statistic computed
based on the permutations $X_{n}^{\phi}$.

Let $F_{Fr_{n}^{\phi}}(s)$ be the permutation distribution function of
$Fr_{n}^{\phi}$ , defined as

$F_{Fr_{n}^{\phi}}(s)=M_{n}^{-1}\sum\limits_{\phi\in\Phi_{n}}1_{Fr_{n}(X_{n}^{\phi})\leq
s}.$

We write the $1-\alpha$ quantile of $F_{Fr_{n}^{\phi}}$ by $c_{1-\alpha,n}$,
defined as $c_{1-\alpha,n}\overset{\Delta}{=}\inf\\{t:F_{Fr_{n}^{\phi}}(s)\geq
1-\alpha\\}$.

Given the quantile $c_{1-\alpha,n}$, the permutation test rejects the null
hypothesis when $Fr_{n}>c_{1-\alpha,n}$.

##  5 Acknowledgements

This work was funded by European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska Curie grant agreement No 812912 for the
project BIGMATH. Computational resources have been provided by the INDACO
Platform, which is a project of High Performance Computing at the University
of Milan.

## References

  * [1] José Ramón Berrendero, Ana Justel, and Marcela Svarc.  Principal components for multivariate functional data.  Computational Statistics & Data Analysis, 55(9):2619–2634, 2011\.
  * [2] A.J. Fridlund.  Human facial expression: An evolutionary view.  Academic Press, 2014.
  * [3] Julien Jacques and Cristian Preda.  Clustering multivariate functional data.  In COMPSTAT 2012, pages 353–366, 2012.
  * [4] J. Ramsay and B.W. Silverman.  Functional data analysis.  Springer, 1997.
  * [5] S. Ullah and F. F. Caroline.  Applications of functional data analysis: A systematic review.  BMC medical research methodology, 2013.
  * [6] J. Wrobel.  Register: Registration for exponential family functional data.  Journal of Open Source Software, 3(22):557, 2018.
  * [7] J. Wrobel, V. Zipunnikov, J. Schrack, and J. Goldsmith.  Registration for exponential family functional data.  Biometrics, 75(1):48–57, 2019.


# Switchback Deflections Beyond the Early Parker Solar Probe Encounters

R. Laker1, T. S. Horbury1, L. Matteini1, S. D. Bale2,3, J. E. Stawarz1, L. D.
Woodham1, T. Woolley1
1Imperial College London, Blackett Laboratory, South Kensington, SW7 2AZ
2Physics Department, University of California, Berkeley, CA 94720-7300, USA
3Space Sciences Laboratory, University of California, Berkeley, CA 94720-7450,
USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Switchbacks are Aflvénic fluctuations in the solar wind, which exhibit large
rotations in the magnetic field direction. Observations from Parker Solar
Probe’s (PSP’s) first two solar encounters have formed the basis for many of
the described switchback properties and generation mechanisms. However, this
early data may not be representative of the typical near-Sun solar wind,
biasing our current understanding of these phenomena. One defining switchback
property is the direction in which the magnetic fields deflects as the
spacecraft samples them. During the first solar encounter, this was primarily
in the tangential direction for the longest switchbacks, which has since been
discussed as evidence, and a testable prediction, of several switchback
generation methods. In this study, we re-examine the deflection direction of
switchbacks during the first eight PSP encounters to confirm the existence of
a systematic deflection direction. We first identify switchbacks exceeding a
threshold deflection in the magnetic field and confirm a previous finding that
they are arc-polarized. In agreement with earlier results from PSP’s first
encounter, we find that groups of longer switchbacks tend to deflect in the
same direction for several hours. However, in contrast to earlier studies, we
find that there is no unique direction for these deflections, although several
solar encounters showed a non-uniform distribution in deflection direction
with a slight preference for the tangential direction. This result suggests a
systematic magnetic configuration for switchback generation, which is
consistent with interchange reconnection as a source mechanism, although this
new evidence does not rule out other mechanisms, such as the expansion of wave
modes.

###### keywords:

Sun: magnetic fields – Sun: heliosphere – solar wind

††pubyear: 2021††pagerange: Switchback Deflections Beyond the Early Parker
Solar Probe Encounters–References

##  1 Introduction

Switchbacks are Alfvénic fluctuations that represent folds in the magnetic
field, rather than changes in magnetic polarity (Balogh et al., 1999;
Neugebauer & Goldstein, 2013; McManus et al., 2020), and have been observed
sporadically in the fast solar wind with previous spacecraft (Balogh et al.,
1999; Horbury et al., 2018). Following observations from Parker Solar Probe
(PSP; Kasper et al., 2019; Bale et al., 2019), switchbacks have become an
active area of research, with many studies describing their properties during
the first two solar encounters (Dudok de Wit et al., 2020; Horbury et al.,
2020; Farrell et al., 2020; Larosa et al., 2021).

One interesting characteristic identified in this early data was the
deflection direction of the magnetic field inside switchbacks. Although Dudok
de Wit et al. (2020) demonstrated that the magnetic field deflections were
isotropic around the Parker spiral direction (Parker, 1958), they did note
that the longest switchbacks in the times series data seemed to have a
preferential direction. For an interval during PSP’s first solar encounter,
Horbury et al. (2020) found that the longest switchbacks deflected in the
tangential direction in the same sense as the Sun’s rotation. Kasper et al.
(2019) also reported a tangential global flow deflection of the solar wind in
the same direction, with the concurrence of these observations leading to the
suggestion that they may be related (Fisk & Kasper, 2020).

Interest in switchbacks from the first encounters prompted many theories
regarding their creation on the Sun, for example, interchange reconnection
(Yamauchi et al., 2004; Fisk & Kasper, 2020); coronal jets (Sterling & Moore,
2020); or flux ropes (Drake et al., 2021). There have also been several
mechanisms proposed to create switchbacks in-situ via interactions between
solar wind streams (Ruffolo et al., 2020; Schwadron & McComas, 2021) or
through the growth of various wave modes because of solar wind expansion
(Squire et al., 2020; Zank et al., 2020; Liang et al., 2021; Mallet et al.,
2021; Shoda et al., 2021). Schwadron & McComas (2021) and Fisk & Kasper (2020)
both argue that their super-Parker spiral and interchange reconnection ideas
naturally explain both the observed tangential deflections from Horbury et al.
(2020) and a global solar wind flow (Kasper et al., 2019). These systematic
deflections are not an explicit prediction of in-situ generation methods,
which involve expanding Alfvénic waves and turbulence (Squire et al., 2020;
Mallet et al., 2021; Shoda et al., 2021), as suggested by Fisk & Kasper
(2020).

Therefore, evidence, or lack thereof, for tangential deflections in
switchbacks can differentiate between theories of their origin. Crucially,
much of the evidence for this property is based on observations from the first
solar encounter, when PSP was connected to a small equatorial coronal hole
(Bale et al., 2019; Badman et al., 2020). These conditions are not
representative of all observations taken to date. Therefore, in this paper we
revisit the idea of tangential deflections in subsequent solar encounters and
aim to remove any ambiguity around this switchback property.

##  2 Methodology

###  2.1 Definition of a Switchback

Following the work of several authors (Dudok de Wit et al., 2020; Horbury et
al., 2020; Laker et al., 2021; Fargette et al., 2021), we consider switchbacks
as deflections of the magnetic field away from the Parker spiral direction .
This represents a stable reference direction that is a based on the physical
speed of the solar wind, and is insensitive to large-scale magnetic field
deflections associated with the switchbacks themselves. We estimate the Parker
spiral direction using the spacecraft position and the solar wind velocity
filtered using a Butterworth filter (of order 2) with a cut-off at $3$ hours.
In particular, we use proton core fits to the SPAN-i measurements (Woodham et
al., 2021) for the solar wind velocity estimates, and fill any data gaps using
moments of the proton distributions measured by the Solar Probe Cup (SPC,
Kasper et al., 2016; Case et al., 2020).

The smooth velocity profile is not influenced by switchbacks that have much
smaller durations than the 3 hour cut-off (Dudok de Wit et al., 2020). A
systematic under/over estimation of the Parker spiral direction may introduce
unwanted results when investigating any preferential deflection direction of
switchbacks. To avoid this, we compare the Parker spiral direction to periods
of manually identified solar wind with minimal switchback activity, finding no
systematic offset in either the $\vec{R}$-$\vec{T}$ or $\vec{T}$-$\vec{N}$
planes in the RTN coordinate system111$\vec{R}$ is radial direction from the
Sun to the spacecraft, $\vec{N}$ is the component of the solar north direction
perpendicular to $\vec{R}$, and $\vec{T}$ completes the right-handed set.

In this paper, we define a switchback as a magnetic deflection with an angle
greater than 45∘ from the Parker spiral direction, as shown in Fig. 1. We
ignore fluctuations with a duration less than 10 seconds and merge switchbacks
that are separated by less than 20 seconds. In order to include complete
observations of the switchbacks, rather than just the measurements at
deflections above 45∘, we define the switchback boundaries as the point at
which the deflection angle crosses a threshold of $30^{\circ}$. This method
may exclude smaller fluctuations from our analysis, but it allows for the
identification of individual switchbacks. In contrast, studying each
individual magnetic vector above a certain threshold will provide a
distribution heavily weighted towards the longest duration switchbacks.

###  2.2 Clock Angle

Figure 1: An example switchback identified using the methods outlined in
Section 2.1. (a) and (b) show the magnetic field components, along with a
dashed line representing their mean average across the switchback. We show the
deflection angle from the Parker spiral in (c), along with our thresholds to
define a switchback at $30^{\circ}$ and $45^{\circ}$. (d) shows how the clock
angle varies inside a switchback. We then compare different methods for
obtaining a single value for clock angle per switchback (defined in the text).

When discussing the deflection direction of switchbacks, we restrict ourselves
to magnetic field data from the FIELDS suite (Bale et al., 2016), since rapid
deflections in the magnetic field direction can move the proton velocity
distribution out of the SWEAP plasma instruments’ field of view (e.g., see
Woolley et al., 2020; Woodham et al., 2021).

Following the work of Horbury et al. (2020) and Laker et al. (2021), we use
the term ‘clock angle’ ($\Phi_{CA}$) to describe the direction of switchback
deflection. The clock angle is the angle made by the magnetic field vector in
a plane normal to the local Parker spiral direction. We define $0^{\circ}$,
$90^{\circ}$, $180^{\circ}$, $270^{\circ}$ as pointing towards +$\vec{N}$,
+$\vec{T}$, -$\vec{N}$, -$\vec{T}$, respectively. To simplify this procedure,
we force the Parker spiral direction to point away from the Sun, which ensures
that a deflection with a +$\vec{T}$ component gives
$0^{\circ}<\Phi_{CA}<180^{\circ}$, regardless of the magnetic field polarity.
We can calculate the clock angle for every magnetic field vector in the
switchback, meaning there are several methods to measure a single value for
$\Phi_{CA}$ per switchback. We have opted for a method that averages
$\Phi_{CA}$ when the deflection angle is within $90^{\circ}\pm 10^{\circ}$, or
a maximum if this threshold is not reached. This method samples the magnetic
field when it is perpendicular to the local Parker spiral direction, which
gave the most consistent results. Fig. 1 demonstrates that there is little
difference between our definition of clock angle ($\Phi_{CA}$) and the clock
angle of the average magnetic components ($\Phi_{avg}$) or the average of the
clock angle for each magnetic vector in the switchbacks
($\left\langle\Phi_{CA}\right\rangle$).

##  3 Results

###  3.1 Arc-Polarized Fluctuations

Figure 2: Two dimensional histograms for the magnetic vectors inside all
switchbacks over the eight encounters. Top panels show the difference between
the clock angle of each vector and our definition of clock angle for an
individual switchback ($\Phi_{CA}$, from Section 2.2) versus normalized time
(a) and deflection angle (b). Panels (b) and (d) are normalized by column,
which shows that there is little variation in $\Phi_{CA}$. Bottom panels show
the difference between $|\vec{B}|$ and the mean, which shows that $|\vec{B}|$
is constant inside switchbacks.

Several authors have already stated that switchbacks are ‘arc-polarized’
(Bruno et al., 2004; Matteini et al., 2014; Horbury et al., 2020), meaning the
magnetic field vectors inside a switchback lying on the arc of a sphere of
radius $|\vec{B}|$. We first test whether this property is consistent across
all available encounters, which is a crucial step, since our results are based
on the assumption that each switchback has a unique deflection direction.
After identifying switchbacks using the methods in Section 2, we plot
superposed epoch analysis for $|\vec{B}|$ in the bottom panels of Fig. 2.
Switchbacks show little variation in $|\vec{B}|$ in both normalized time
(panel (c)) and deflection angle (panel (d)), which agrees with previous
results (Krasnoselskikh et al., 2020), confirming the Alfvénic nature of
switchbacks across all encounters.

Fig. 2 also shows $\Phi_{CA}$ inside all switchbacks, which again shows little
variation compared to our definition. Therefore, we can conclude that
switchbacks are arc-polarized, i.e the fluctuations lie on an arc of the
$|\vec{B}|$ sphere and return along the path they deflected from, as shown by
the symmetry about $0.5$ normalized time (panel (a)). Alfvénic fluctuations do
not have to exhibit this property, and may in principle keep following the arc
past $180^{\circ}$ deflection angle, which would case a $180^{\circ}$ change
in $\Phi_{CA}$. Such a result may reflect an inherent symmetry of the physical
switchback structure, which would be the case for a typical ‘S-shaped’ curve
(Kasper et al., 2019).

###  3.2 Switchback deflections

Figure 3: Contextual solar wind plasma data (a-d) for the third PSP encounter,
with the $\Phi_{CA}$ of identified switchbacks plotted in (e), where scatter
size is proportional to switchback duration (similar to Horbury et al.
(2020)). The histogram in (f) shows the distribution of all switchback
deflections in the encounter, which indicates a tendency for tangential
deflections.

Since switchbacks are arc-polarized, with the magnetic field returning along
the same path, a single value can describe their deflection direction. In Fig.
3, we plot the deflection direction in terms of $\Phi_{CA}$ for each
switchback in the encounter. It is clear that an individual switchback can
have any value of $\Phi_{CA}$, which is a key feature that holds over all the
encounters. However, when considering the distribution of deflection
directions across all switchbacks within an encounter, Fig 3(f), we find that
$\Phi_{CA}$ is more likely to be around $100^{\circ}$ and $270^{\circ}$. This
anisotropy is also present for the longest switchbacks, which seem to follow
the overall distribution.

We summarize the switchback deflections for encounters 1 to 8 by wrapping the
histogram in Fig 3(f) onto polar axes, as shown in Fig 4. The Alfvénic nature
of switchbacks allows us to distinguish between positive (red) and negative
(blue) polarity background field as this will influence the direction of
velocity deflection - i.e. for positive polarity the direction of magnetic
deflection is opposite to the velocity deflection, which would rotate the
$\Phi_{CA}$ by $180^{\circ}$. A null hypothesis of a uniform distribution is
plotted for each histogram, which was rejected using Chi-Square test with a 1%
confidence level for all of the observed $\Phi_{CA}$ distributions across the
encounters. Therefore, we conclude that switchback deflections are not
uniformly distributed, with a preference towards tangential, rather than
normal, deflections being clear in Fig 4.

Interestingly, this preference is not purely tangential (i.e. $90^{\circ}$ or
$270^{\circ}$), with negative polarities showing a tendency for an axis of
symmetry along $\sim 135^{\circ}$ or $\sim 315^{\circ}$ (e.g. Fig. 4 encounter
2) and positive having symmetry along $\sim 225^{\circ}$ or $\sim 45^{\circ}$
(e.g. Fig. 4 encounter 7). This statement holds true even when using other
definitions of clock angle (Sec. 2.2). While the origins of this effect are
not clear, it could be related to a systematic magnetic configuration either
at the generation site (e.g. a tilted coronal hole boundary), or as the
switchback propagates through the solar wind. It is important to note that not
all encounters show this offset symmetry axis (e.g. encounter 1 and 4),
suggesting it may be related to the large-scale solar wind structure.

While Fig. 4 shows the distribution for all switchbacks, Fig. 3(e) allows us
to see the distribution of the longest switchbacks. Similar to previous
reports (Horbury et al., 2020; Laker et al., 2021), these longer switchbacks
exhibit ‘clustering’ - i.e. successive longer switchbacks tend to deflect in
the same direction. Such behavior occurs throughout all the encounters,
although there is no specific $\Phi_{CA}$ at which this happens. During the
first encounter, Horbury et al. (2020) showed that longer switchbacks
deflected in $+\vec{T}$, which was then subsequently linked to the transverse
flow in Kasper et al. (2019). However, when looking at the $\Phi_{CA}$ of
switchbacks in encounter 1 (Fig. 4) we see that there is no preference for
$+\vec{T}$ over $-\vec{T}$, meaning that the earlier result from Horbury et
al. (2020) was most likely during a ‘cluster’ of longer switchbacks.

Figure 4: Distribution of all switchback clock angles for each PSP encounter,
split into positive (red) and negative (blue) polarity groups. We show an
expected uniform distribution for both positive (dashed circle) and negative
(solid circle) polarities. This demonstrates that switchback deflections are
not isotropic, with a preference for deflections in the tangential direction
($90^{\circ}$ and $270^{\circ}$).

##  4 Discussion

As Fisk & Kasper (2020) highlighted, a preference for tangential magnetic
deflections is a natural consequence of interchange reconnection (Fisk &
Schwadron, 2001; Fisk, 2005). In this generation mechanism, open magnetic
field lines are dragged over closed loops at the Sun, which then reconnect and
launch switchbacks into the solar wind. Therefore, this process is
intrinsically linked to the differential rotation of the corona and the solar
surface, which would create switchbacks with orientations in the $\pm\vec{T}$
directions (Fisk & Kasper, 2020, Fig. 1). Such a mechanism can also explain
the clustering behavior in Fig. 3, with successive switchbacks being launched
from the same reconnection site on the Sun. A stable magnetic configuration
surrounding the reconnection site would then lead to successive switchbacks
having similar deflection directions (Zank et al., 2020; Fisk & Kasper, 2020),
with the distributions in Fig. 4 reflecting the magnetic field orientation at
the Sun.

The super-Parker spiral concept from Schwadron & McComas (2021) also predicts
a preferential tangential deflection for switchbacks. Here, the magnetic field
foot-point passes from slow to fast solar wind, creating a kink in the
magnetic field with a velocity deflection in the same direction as the Sun’s
rotation. Assuming an Alfvénic correlation, this would create magnetic field
deflections with a $+\vec{T}$ component for negative polarity and $-\vec{T}$
for positive polarity. Although this was observed during a period of the first
solar encounter (Horbury et al., 2020), Fig. 4 demonstrates that the
preference for tangential deflections is not restricted to a single direction
(e.g. not only $+\vec{T}$ component for negative polarity). It is important to
note that this argument does not consider the propagation of switchbacks,
which could account for a change in orientation.

The observed clustering and preference for tangential deflections are harder
to reconcile for those mechanisms that are based on the expansion of randomly
orientated in-situ fluctuations (Squire et al., 2020; Mallet et al., 2021;
Shoda et al., 2021). There would have to be some systematic effect, either in
the seeding fluctuations or on the switchbacks as they propagate. For example,
if the Parker spiral could influence the evolution of switchbacks, then this
would provide a systematic effect on switchbacks, regardless of their
generation mechanism.

Recently, Bale et al. (2021) and Fargette et al. (2021) demonstrated that
patches of switchbacks are organized into structures that are consistent with
funnels closer to the Sun, most likely related to super-granules, which can
exhibit large-scale magnetic field deflections. Therefore, it will be
interesting to see if simulations can study the influence of large-scale
deflections on switchbacks.

However, there are several ways to explain the clustering behavior because all
PSP observations are inherently single point measurements. The most obvious is
that there are successive switchbacks with longer duration that have similar
magnetic deflection directions. As discussed in Section 4, this is naturally
explained by intermittent events being launched from the Sun, but is harder to
explain with expanding waves. While there are some suggestions that plasma is
hotter inside switchbacks (Woodham et al., 2021), there is no agreed plasma
signature that can be used to identify switchbacks. Therefore, with single
point measurements, we may misindentify a single coherent structure as several
smaller switchbacks in the PSP data. These ‘switchbacks’ would actually be
part of the same structure, which could explain the similarity between their
deflections directions. Such an issue could be explored by identifying
switchback edges through the presence of wave activity (Agapitov et al., 2020;
Mozer et al., 2020; Larosa et al., 2021) and magnetic field dropouts (Farrell
et al., 2020), or by finding a unique plasma signature for the inside of a
switchback.

The arc-polarized behavior demonstrated in Section 3.1 implies a symmetry
within the switchback structure, since the magnetic field deflects back along
the same path. Such an observed property is a potential constraint for wave
models of switchback generation. This can be visualized with an S-shaped
magnetic field line in 2D, as shown in Kasper et al. (2019) and is reminiscent
of a ‘phase-steepened Alfvén wave’ (Tsurutani & Ho, 1999; Tsurutani et al.,
2002; Tsurutani et al., 2018). Since this arc-polarized behavior exists for
all clock angles, this symmetry in the magnetic field must be present in three
dimensions.

The spacecraft can pass through the 3D structure of the switchback in
different ways, as discussed in Zank et al. (2020). Horbury et al. (2020) and
(Laker et al., 2021) suggested that there is a relationship between the
deflection direction of the switchback and the direction in which PSP cuts
through it. Therefore, if switchbacks have an elongated structure, then the
deflection direction will influence the duration of the switchback in the
spacecraft time series (Macneil et al., 2020; Laker et al., 2021). With
increasing spacecraft tangential speed, it was expected that PSP would travel
through the side of switchbacks (assuming they are long and thin (Laker et
al., 2021)), leading to durations being less affected by deflection direction.
We do not observe this effect in the later encounters, suggesting that the
link between deflection direction and duration is not as simple as first
thought. For example, considering the most simple ‘S-shaped’ curve in the
magnetic field (e.g. Kasper et al., 2019) and applying a threshold of
$45^{\circ}$ for the magnetic deflection, as shown in Section. 2.1. This will
produce a long, thin shape that is aligned away from the background magnetic
field, similar to Laker et al. (2021) that can have a different length and
alignment based on the ‘tightness’ of the S shape. This simple example
demonstrates the non-trivial relationship between the true switchback
structure in the solar wind and what we observe in the PSP time series, which
will continue to be a major challenge in the future.

##  5 Conclusions

Many of the current switchback generation theories were based on observations
from the first few solar encounters from PSP, with the deflection direction of
switchbacks being used as evidence. However, the solar wind conditions in the
first encounter were not representative of later encounters, which are now
also closer to the Sun. Therefore, in this study, we have investigated the
deflection directions of switchbacks across the first eight encounters. We
have identified switchbacks as deflections in the magnetic field, and
confirmed previous results that they have almost constant $|\vec{B}|$ and
deflection direction, meaning they can be described as arc-polarized. We
showed that the longer switchbacks still exhibit clustering behavior, where
successive switchbacks have similar orientations over periods of hours.
Although switchbacks could deflect in any direction, we demonstrated that they
have a tendency to deflect in the tangential, rather than normal, direction.
Such a systematic observation, and the clustering behavior, would be natural
consequences of some switchback generation methods, such as interchange
reconnection. More research is needed to rule out other generation theories,
such as the evolution of switchbacks and the effect of large-scale magnetic
deflections closer to the Sun. We found no obvious trend in switchback
direction or clustering with radial distance from the Sun, although this may
be revealed as PSP travels closer in future encounters.

## Acknowledgements

RL was supported by an Imperial College President’s Scholarship, TSH by STFC
ST/S000364/1, TW by ST/N504336/1, LDW by ST/S000364/1. JES is supported by the
Royal Society University Research Fellowship URF\R1\201286. The SWEAP and
FIELDS teams acknowledge support from NASA contract NNN06AA01C. This work has
made use of the open source and free community-developed space physics
packages HelioPy (Stansby et al., 2021) and SpiceyPy (Annex et al., 2021).

## Data Availability

The data used in this study are available at the NASA Space Physics Data
Facility (SPDF): https://spdf.gsfc.nasa.gov.

## References

  * Agapitov et al. (2020) Agapitov O. V., et al., 2020, The Astrophysical Journal, 891
  * Annex et al. (2021) Annex A., et al., 2021, AndrewAnnex/SpiceyPy: SpiceyPy 4.0.1, Zenodo, doi:10.5281/zenodo.4883901
  * Badman et al. (2020) Badman S. T., et al., 2020, The Astrophysical Journal Supplement Series, 246, 23
  * Bale et al. (2016) Bale S. D., et al., 2016, Space Science Reviews, 204, 49
  * Bale et al. (2019) Bale S. D., et al., 2019, Nature, 576, 237
  * Bale et al. (2021) Bale S. D., et al., 2021, The Astrophysical Journal, 923, 174
  * Balogh et al. (1999) Balogh A., Forsyth R. J., Lucek E. A., Horbury T. S., Smith E. J., 1999, Geophysical Research Letters, 26, 631
  * Bruno et al. (2004) Bruno R., Carbone V., Primavera L., Malara F., Sorriso-Valvo L., Bavassano B., Veltri P., 2004, Annales Geophysicae, 22, 3751
  * Case et al. (2020) Case A. W., et al., 2020, The Astrophysical Journal Supplement Series, 246, 43
  * Drake et al. (2021) Drake J. F., et al., 2021, Astronomy & Astrophysics
  * Dudok de Wit et al. (2020) Dudok de Wit T., et al., 2020, The Astrophysical Journal Supplement Series, 246, 39
  * Fargette et al. (2021) Fargette N., et al., 2021, The Astrophysical Journal, 919, 96
  * Farrell et al. (2020) Farrell W. M., MacDowall R. J., Gruesbeck J. R., Bale S. D., Kasper J. C., 2020, The Astrophysical Journal Supplement Series, 249, 28
  * Fisk (2005) Fisk L. A., 2005, The Astrophysical Journal, 626, 563
  * Fisk & Kasper (2020) Fisk L. A., Kasper J. C., 2020, The Astrophysical Journal, 894, L4
  * Fisk & Schwadron (2001) Fisk L. A., Schwadron N. A., 2001, The Astrophysical Journal, 560, 425
  * Horbury et al. (2018) Horbury T. S., Matteini L., Stansby D., 2018, Monthly Notices of the Royal Astronomical Society, 478, 1980
  * Horbury et al. (2020) Horbury T. S., et al., 2020, The Astrophysical Journal Supplement Series, 246, 45
  * Kasper et al. (2016) Kasper J. C., et al., 2016, Space Science Reviews, 204, 131
  * Kasper et al. (2019) Kasper J. C., et al., 2019, Nature, 576, 228
  * Krasnoselskikh et al. (2020) Krasnoselskikh V., et al., 2020, The Astrophysical Journal, 893, 93
  * Laker et al. (2021) Laker R., et al., 2021, Astronomy & Astrophysics
  * Larosa et al. (2021) Larosa A., et al., 2021, Astronomy & Astrophysics
  * Liang et al. (2021) Liang H., Zank G. P., Nakanotani M., Zhao L.-L., 2021, The Astrophysical Journal, 917, 110
  * Macneil et al. (2020) Macneil A. R., Owens M. J., Wicks R. T., Lockwood M., Bentley S. N., Lang M., 2020, Monthly Notices of the Royal Astronomical Society, 494, 3642
  * Mallet et al. (2021) Mallet A., Squire J., Chandran B. D. G., Bowen T., Bale S. D., 2021, Evolution of Large-Amplitude Alfvén Waves and Generation of Switchbacks in the Expanding Solar Wind, doi:10.3847/1538-4357/ac0c12
  * Matteini et al. (2014) Matteini L., Horbury T. S., Neugebauer M., Goldstein B. E., 2014, Geophysical Research Letters, 41, 259
  * McManus et al. (2020) McManus M. D., et al., 2020, The Astrophysical Journal Supplement Series, 246, 67
  * Mozer et al. (2020) Mozer F. S., et al., 2020, The Astrophysical Journal Supplement Series, 246, 68
  * Neugebauer & Goldstein (2013) Neugebauer M., Goldstein B. E., 2013, in AIP Conference Proceedings. American Institute of PhysicsAIP, pp 46–49, doi:10.1063/1.4810986
  * Parker (1958) Parker E. N., 1958, The Astrophysical Journal, 128, 664
  * Ruffolo et al. (2020) Ruffolo D., et al., 2020, The Astrophysical Journal, 902, 94
  * Schwadron & McComas (2021) Schwadron N. A., McComas D. J., 2021, The Astrophysical Journal, 909, 95
  * Shoda et al. (2021) Shoda M., Chandran B. D. G., Cranmer S. R., 2021, The Astrophysical Journal, 915, 52
  * Squire et al. (2020) Squire J., Chandran B. D. G., Meyrand R., 2020, The Astrophysical Journal, 891, L2
  * Stansby et al. (2021) Stansby D., et al., 2021, Heliopython/Heliopy: HelioPy 0.15.4, Zenodo, doi:10.5281/zenodo.5090511
  * Sterling & Moore (2020) Sterling A. C., Moore R. L., 2020, The Astrophysical Journal Letters, 896
  * Tsurutani & Ho (1999) Tsurutani B. T., Ho C. M., 1999, Reviews of Geophysics, 37, 517
  * Tsurutani et al. (2002) Tsurutani B. T., et al., 2002, Geophysical Research Letters, 29, 86
  * Tsurutani et al. (2018) Tsurutani B. T., Lakhina G. S., Sen A., Hellinger P., Glassmeier K. H., Mannucci A. J., 2018, Journal of Geophysical Research: Space Physics, 123, 2458
  * Woodham et al. (2021) Woodham L. D., et al., 2021, A&A, 650
  * Woolley et al. (2020) Woolley T., et al., 2020, Monthly Notices of the Royal Astronomical Society, 498, 5524
  * Yamauchi et al. (2004) Yamauchi Y., Suess S. T., Steinberg J. T., Sakurai T., 2004, Journal of Geophysical Research: Space Physics, 109
  * Zank et al. (2020) Zank G. P., Nakanotani M., Zhao L.-L., Adhikari L., Kasper J., 2020, The Astrophysical Journal, 903, 1