\section{Introduction}

In the hierarchical structure formation scenario, galaxy clusters are
formed by continuous merging of smaller infalling groups and finally
achieve dynamical equilibrium \citep[e.g.][]{ps74, mbb+09, cog13}.
As the largest gravitational bound systems in the universe, clusters
of galaxies consist of a large number of member galaxies embedded in
the intracluster medium (ICM) and the dark matter halo.
Often seen inside a galaxy cluster are the substructures in the
spatial distribution of galaxies \citep[e.g.,][]{ds88,wh13}.
Because the hot gas in the ICM observed in the X-ray band is in
kinematics equilibrium under the assumption of the energy density
equipartition between the ICM and galaxies, many substructures have
also been detected in the X-ray images  \citep[e.g.,][]{kbp+01}.

The dynamical state of galaxy clusters describes if clusters are
relaxed dynamically or are still in a merging process of subclusters.
Clusters of galaxies are therefore broadly divided into relaxed
clusters and dynamically disturbed clusters. Relaxed clusters roughly
reach the virial equilibrium, showing a symmetrical structure in the
distribution of galaxies and the intracluster gas around their
brightest cluster galaxy (BCG).
While dynamically disturbed clusters, such as the Bullet cluster
\citep{mgd+02}, obviously deviate from the virial equilibrium with
subclusters around two or more very bright galaxies
\citep[e.g.,][]{cd96,bbr+02,wh13} or clear substructures in X-ray
images \citep[e.g.,][]{mefg95,mjfv08,yda+16}.
The dynamical state of galaxy clusters is ideally expressed by the
three-dimensional (3D) velocity distribution of member galaxies
\citep[e.g.][]{cd96,ets+10,evn+12} or the ICM
\citep[e.g.][]{db06,lyt+15,lyt+16,yda+16}. In practice, the dynamical
state can be roughly indicated by the 1D radial velocity or redshift
distribution
\citep[e.g.,][]{yv77,ds88,wb90,ssg99,hmp+04,hph+09,rph18,yds+18}, or
the 2D positions of member galaxies in the sky plane projected from
their real spatial distribution
\citep[e.g.,][]{gb82,wb90,fk06,rbp+07,as10,evn+12,wh13,ltl+18} or the
projected hot gas distribution shown in X-ray or microwave-band images
of clusters \citep[e.g.,][]{mefg95,kbp+01,bfs+05,cbs12,cds+18}.

Substructures in such a projected two dimensional distribution of
member galaxies or hot gas have been parameterized quantitatively to
describe the dynamical state of galaxy clusters, such as the asymmetric
factors or the probability distribution based on the optical data
\citep[e.g.][]{wb90,rbp+07,wh13}, or the concentration index
\citep[e.g.,][]{mefg95,srt+08}, the centroid shift \citep[e.g.,][]{
  mefg95,pfb+06,mjfv08} and the power ratio \citep[e.g.,][]{bt95,
  bt96, bpa+10} based on the data of the X-ray images, see details in
Section~\ref{knownpara}.

\citet{wh13} derived relaxation parameters for 2092 rich clusters
based on photometric data by using the asymmetry, the ridge flatness
and the normalized deviation of the smoothed optical maps, which is
the largest sample of clusters to our knowledge with dynamical
parameters.
\citet{kgm+19} found a good agreement between the galaxy density maps
and the X-ray brightness maps, and therefore estimated the dynamical
parameters for 890 clusters by the method given by \citet{wh13}.
Recently, \citet{cds+18} also tried to parameterize the morphological
characterization of synthetic maps of the Sunyaev-Zel\'dovich (SZ)
effect for a sample of 258 simulated clusters.
\citet{ltl+18} estimated optical substructure parameters for 72
clusters which have {\it Chandra} X-ray images and/or {\it Planck}
Sunyaev-Zel\'dovich maps by using optical imaging and spectroscopic
data.
\citet{rlt+18} showed the close correlation between dynamical
parameters (difference between the optical and X-ray center, the
projected offset of the BCG from other cluster centroids) and the
offset from the scaling relations for clusters.
\citet{zhk+20} also calculated the offsets between BCGs and gas center
as traced by SZ effect ({\it SPT}) and/or X-ray ({\it Chandra} and
{\it XMM-Newton}) images as dynamical parameter for 288 massive
clusters.

The dynamical parameters have been involved to many studies, e.g. the
cosmological constraints based on galaxy clusters with different X-ray
morphologies \citep[e.g.][]{mefg95},
the evolution of galaxy clusters in the cosmological model of
structure formation \citep[e.g.][]{bt96},
the presence of radio halos and relics in merging clusters and the
formation of mini-halos in relaxed clusters
\citep[e.g.,][]{ceg+10,ceb+13,ccb+15, yhw15},
the mass estimation \citep[e.g.,][]{rlt11} and radial mass profile
\citep[e.g.,][]{bap+19} for clusters of galaxies,
the scaling relations for cluster mass estimations
\citep[e.g.,][]{omb+06,crb+07,zfb+08,zjc+13},
the dynamical-state dependence of galaxy luminosity functions
\citep[e.g.,][]{wh15},
and activity of super massive black holes in BCGs
\citep{kvc+15,yhw16}.

In fact, any definition of dynamical parameters may be biased because
we cannot get the real 3D distribution and velocities. Nevertheless,
the identification and analysis of substructures of any projected
galaxy distribution or X-ray images should give the lower limit of
dynamical parameters ideally derived from the unavailable real 3D
distribution.  Therefore, it is still valuable to get dynamical
parameters from the X-ray images of galaxy clusters.
The {\it Chandra} satellite has observed about 1000 clusters of
galaxies in the X-ray band, and the high resolution {\it Chandra}
image is ideally to reveal the substructures from the ICM.
In this paper, we derive a set of dynamical parameters for 964
clusters based on the archival data of the {\it Chandra X-ray
  Observatory}.

\begin{figure}
\centering
\includegraphics[angle=270,width=0.42\textwidth]{z_distribution.eps}
\caption{Redshift distribution for 911 of 964 clusters with
    available redshift.}
\label{redshift}
\end{figure}

In Section 2, we describe how we get and process the {\it Chandra}
data. In Section 3, we define and calculate various dynamical
parameters. Comparison of the derived parameters and discussions on
the impact of dynamical state on the scaling relations are presented
in Section 4. A summary is given in Section 5. Throughout this paper,
we assume a flat $\Lambda$CDM cosmology taking $H_0=70$ km~s$^{-1}$
Mpc$^{-1}$, $\Omega_m=0.3$ and $\Omega_{\Lambda}=0.7$.

\section{Data collection and processing}

\subsection{Chandra data for galaxy clusters}
\label{datacoll}
We collect the observation information for clusters of galaxies
archived by the {\it Chandra} in two approaches.

First, we get the information for all observations which categorized
as ``Clusters of galaxies'' in the {\it Chandra} archival
database\footnote{http://cda.harvard.edu/chaser/}. As a result, more
than 2000 observations are obtained, often several observations are
made for one cluster. For clusters observed more than one time, we
generally choose the observation with largest photon number. However,
we select the observation for nearby clusters (mostly with $z<0.05$)
not only with a sufficient photon number but also with a good coverage
of CCDs. After removing reduplicate observations, we obtained data for
915 clusters.

Second, we cross match several widely used cluster catalogs: the Abell
\citep{aco89}, MCXC \citep{pap+11}, redMaPPer \citep{rrb+14}, CAMIRA
\citep{o14}, WH15 \citep{whl12,wh15} and WHY18 \citep{why18}, with
archived {\it Chandra} observations to find serendipitous observed
clusters with criteria as following: (1) the search radius from the
nominal center of a {\it Chandra} observation is set as 10 arcminutes
so that a cluster in the field of view of a {\it Chandra} ACIS-I
observation with a side-length of 16.9 arcminutes can be covered, (2)
the exposure time is larger than 10 kiloseconds and the average count
rate is larger than 2 Hz, so that there are enough photons to work on,
(3) the data status should be archived, (4) observations in all
science categories of {\it Chandra} are selected except for ``Clusters
of galaxies'' (the category has been done), (5) the exposure
instrument is limited with ACIS and without any gratings, (6) the
observation types are constrained as GO (Guest Observer) and GTO
(Guaranteed Time Observation), and the exposure mode is set as TE
(Timed Exposure). With these steps, after merging the outputs for
different cluster catalogs and checking the actual coordinates in
X-ray images, we finally obtain another 49 serendipitously observed
clusters.

\label{imgproc}
\begin{figure*}
\centering
\includegraphics[angle=270,width=55mm]{fig2/01.eps}
\includegraphics[angle=270,width=55mm]{fig2/02.eps}
\includegraphics[angle=270,width=55mm]{fig2/03.eps}\\[2mm]
\includegraphics[angle=270,width=55mm]{fig2/04.eps}
\includegraphics[angle=270,width=55mm]{fig2/05.eps}
\includegraphics[angle=270,width=55mm]{fig2/06.eps}\\[2mm]
\caption{Example images for 6 galaxy clusters, and images for all
  clusters are available on the web page
  http://zmtt.bao.ac.cn/galaxy\_clusters/dyXimages/. Cluster name,
  redshift and contour levels are written on the top of each
  panel. Surface brightness is logarithmically indicated as the color
  bar on the right of each panel. We select a clean region around
  every clusters to calculate the mean brightness $\langle S_{bg}
  \rangle$, and its fluctuations $\sigma$. The white cross stands for
  the center $(x_0,y_0)$ of the best fitted model. The adaptive
  dynamical parameters are calculated within the white ellipse, which
  also defines the actual region for the morphology parameters of
  clusters. The red circle in the bottom-left corner of each panel
  indicates the smooth scale generally with a diameter of 30 kpc. The
  scale of 500 kpc is plotted on the bottom-right corner to indicate
  the physical scale of clusters, which is also the radius size used
  to calculate the concentration index, the centroid shift and the
  power ratio.}
\label{examples}
\end{figure*}

In total, we have a sample of 964 galaxy clusters with good {\it
  Chandra} X-ray image data (see Table 1), which mostly have a lower
redshift less than 0.5 (see Figure~\ref{redshift}) but some high
redshift clusters up to $z\sim1.5$ (e.g., SPT-CL J2040-4451) do come
out. The estimated masses of some clusters \citep[e.g.,][]{pap+11} are
in the range from $10^{13} M_{\odot}$ (e.g., MCXC J1242.8+0241) to
$2.0\times10^{15} M_{\odot}$ (e.g., MACS J0417.5-1154).

\input tab1.tex

\subsection{Image processing}

The {\it Chandra} data for these clusters of galaxies are processed in
the standard manner by using CIAO 4.9 \citep{fma+06} and calibration
database CALDB 4.7.2. Data files of all observations are downloaded
with CIAO command {\tt download\_chandra\_obsid}, and reprocessed by
{\tt chandra\_repro} to make sure the latest calibration products are
used. To be consistent with previous works \citep[e.g.,][]{srt+08},
the photons are filtered in 0.5-5 keV band. To remove flare events, we
use the tool {\tt lc\_clean} to cut out data in time intervals where
the photon count rate deviates from their mean value more than
20\%. For few clusters with giant flares, the mean count rate of
photons is not a good reference, thus time intervals for giant flares
are removed manually. Point sources are detected automatically with
the CIAO routine {\tt wavdetect} with a proper threshold which should
not be too high to miss some faint point sources, or not too low to
lead a false detection of some fluctuations as being point
sources. The size of point source model should not be too small to
miss some relatively large member or foreground galaxies, or too large
to confuse the core of relaxed clusters as point sources. We check the
list of point sources very carefully and then remove the real point
sources. We use the tool {\tt dmfilth} to fill ``holes'' of subtracted
point sources according to the brightness of ambient regions. Images
of clusters are also exposure corrected and background subtracted
accordingly. Since the pixel stands for different physical scales for
different redshifts, the X-ray images of clusters are smoothed to 30
kpc by a Gaussian function, except for 4 very close clusters NGC 1395,
Virgo, MCXC J1242.8+0241 and MCXC J1315.3-1623 with $z<0.01$ for which
the smooth scale are set to 10 kpc. For few clusters without available
redshifts, the images are smoothed to 5 arcseconds. The examples of
final X-ray images are shown in Figure~\ref{examples}.

\section{Dynamical parameters from X-ray images}
\label{dp}

In this section, we obtain three widely used dynamical parameters,
i.e., the concentration index $c$, the centroid shift $\omega$, the
power ratio $P_3/P_0$, and compare our results with those available in
literature. Then, we define two adaptive dynamical parameters, i.e.,
the profile parameter $\kappa$ and the asymmetry factor $\alpha$, in
the best fitted elliptical region. We combine them to make the
morphology index $\delta$, which is an excellent indicator for the
dynamical state of galaxy clusters.

\subsection{Three widely used dynamical parameters}
\label{knownpara}

\subsubsection{The concentration index, $c$}
\label{conindex}
Relaxed clusters usually host a very luminous cool core in their
center \citep[e.g.,][]{fnc84,f94,mbb+12}, while the core of disturbed
clusters generally have been destroyed by violent merger
events. \citet{srt+08} defined the concentration index as the ratio of
X-ray fluxes integrated in the centroid and the whole regions of
galaxy clusters to quantify their dynamical state. The concentration
index $c$ is calculated in two circular regions with the core radius
of 100 kpc and the outer radius of 500 kpc
\citep[e.g.,][]{ceg+10,ceb+13}, i.e.,
\begin{equation}
  c=\frac{S_{100~\rm kpc}}{S_{500~\rm
      kpc}}=\frac{\sum\limits_{R<100\rm~kpc}f_{\rm
      obs}(x_i,y_i)}{\sum\limits_{R<500\rm~kpc}f_{\rm obs}(x_i,y_i)},
\label{c}
\end{equation}
where $f_{\rm obs}(x_i,y_i)$ means the observed X-ray flux at pixel
$(x_i,y_i)$. With the definition in Equation~\ref{c}, we set the two
radii as being 100 kpc and 500 kpc, and calculate the concentration
index $c$ for clusters in our sample, as listed in Table 1. The
concentration index for some low-redshift clusters (e.g., the Vigro
cluster) cannot be calculated due to the size of 500 kpc is much out
of the coverage of the {\it Chandra} CCD, and they are denoted as
``--'' in the full Table~1.

\begin{figure}
  \begin{center}
\includegraphics[angle=0,width=0.23\textwidth]{fig3/ceg_c.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig3/der+16a_c.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig3/der+16b_c.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig3/der+16b_c1.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig3/ajf+17a_c.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig3/ajf+17b_c.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig3/zrs+17_c.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig3/lty+18_c.eps}
\caption{Comparison for the concentration indexes that we obtained
  with those in literatures. The solid line in panels a and b and
  dotted line in panel d indicate equivalent values in X and Y axes
  because the values in both axes are calculated with the same radii
  for the center and outer regions. The Spearman rank-order
  correlation coefficient $r_{\rm s}$ and the relevant significance
  $p_{\rm s}$ are labelled in the right-bottom corner of each
  panel. Data for the X-axis are our values and for the Y-axis are
  obtained from: C10=\citet{ceg+10}, C13=\citet{ceb+13},
  D16=\citet{der+16}, A17=\citet{ajf+17}, Z17=\citet{zrs+17} and
  L18=\citet{lty+18}.}
\label{ccorr}
\end{center}
\end{figure}

Although the concentration index has been derived by a lot of authors
from {\it Chandra} and/or {\it XMM-Newton} images, the radii of the
central and the outer regions are chosen differently. In
Figure~\ref{ccorr}, we compare our values with those calculated in
literatures. \citet{ceg+10} defined the concentration index as
$c=S_{100\rm~kpc}/S_{500\rm~kpc}$ and calculated the $c$ for 32
clusters based on the {\it Chandra} images and for 7 clusters in a
later paper \citep{ceb+13}. Since we set the same core and outer radii
as \citet{ceg+10,ceb+13}, we obtain almost the same values (see
Figure~\ref{ccorr}a). Very good correlations can be found between our
values of $c=S_{100\rm~kpc}/S_{500\rm~kpc}$ and those obtained by
\citet{der+16} in Figure~\ref{ccorr}b for 25 clusters who also defined
the core radius of 100 kpc and the outer radius of 500 kpc, and in
Figure~\ref{ccorr}c and Figure~\ref{ccorr}d for 19 clusters with the
core radius of $0.15~R_{500}$ and the outer radius of $0.5~R_{500}$
which we also calculated. Here $R_{500}$ is the radius of a galaxy
cluster within which the matter density of a cluster is 500 times of
the critical density of the universe. \citet{ajf+17} worked on a large
sample of 214 clusters with the {\it Chandra} data, and calculated the
concentration index with two definitions, i.e.,
$c=S_{40\rm~kpc}/S_{400\rm~kpc}$ and
$c=S_{0.15~R_{500}}/S_{R_{500}}$. Again, good consistences are found
between our values of $c=S_{100\rm~kpc}/S_{500\rm~kpc}$ and those in
\citet{ajf+17}, see Figure~\ref{ccorr}e and Figure~\ref{ccorr}f. The
outliers, i.e., A115, A2440, A3716 and RXC J1414.2+7115, in
Figure~\ref{ccorr}e and Figure~\ref{ccorr}f are merging clusters with
a bi-model brightness distribution (see the image of A115 in
Figure~\ref{examples} as an example). For these clusters we take the
best-fitted center (see Section 3.2.1) located between the two
subclusters similar to \citet{ceg+10}, while \citet{ajf+17} took the
X-ray peak of one subcluster as the center, which cause the value
difference. There are 53 clusters in common among our sample and
\citet{zrs+17}, inverse correlation appears with large scatter in
Figure~\ref{ccorr}g because they define reversely as $c_{\rm
  L}=S_{R_{500}}/S_{0.2~R_{500}}$. \citet{lty+18} calculated
$c=S_{40\rm~kpc}/S_{400\rm~kpc}$ for 41 clusters from the {\it
  Chandra} images. We also get a good agreement between our values and
results of \citet{lty+18} for 41 clusters as shown in
Figure~\ref{ccorr}h. The Spearman rank-order correlation coefficient
$r_{\rm s}$, and the relevant significance $p_{\rm s}$, are marked
inside each panel to show the reliability of correlations
\citep[see][p. 640]{ptvf92}. Here the small value of $p_{\rm s}$
indicates a significant correlation between the two values.

The concentration index can be calculated within the widely used,
mass-related radius $R_{500}$, but there are some limitations for a
large sample of clusters. First, the radius $R_{500}$ is the byproduct
during the mass estimation for clusters, which (1) requires redshifts
that for 53 clusters in our sample are missing, (2) needs good quality
of X-ray image with sufficient photons for spectrum analysis and (3)
is estimated with the assumption of spherically symmetry and virial
equilibrium that is invalid for many clusters as seen in this
paper. Second, the aperture with a radius of $R_{500}$ is also too
large for the {\it Chandra} CCD coverage for many clusters. The median
$R_{500}$ for 1,742 clusters in \citet{pap+11} is about 0.8 Mpc, the
corresponding aperture exceeds 10 arcminutes when $z<0.15$ (note that
the scale of the whole field of view for ACIS-I is 16.9'$\times$16.9'
and ACIS-S is 8.3'$\times$50.6').

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=0.23\textwidth]{fig4/ceg_w.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig4/wbsa13_w.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig4/wbc13_w.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig4/der+16a_w.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig4/der+16b_w.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig4/der+16b_w1.eps}
\caption{Comparison for the centroid shifts we calculated with those in
  literatures. The solid line in the panels a and d indicates
  equivalent values in X and Y axes if the radius is set in the same
  value of 500 kpc, the dotted line in panel f is also an equivalent
  line but for the radius of $0.5~R_{500}$. The correlation parameter
  $r_{\rm s}$ and $p_{\rm s}$ are labelled in the right-bottom corner
  of each panel. Data in the X-axis are our values and in the Y-axis
  are obtained from: C10= \citet{ceg+10}, C13=\citet{ceb+13},
  W13a=\citet{wbsa13}, W13b=\citet{wbc13} and D16=\citet{der+16}.}
\label{wcorr}
\end{center}
\end{figure}

\subsubsection{The centroid shift, $\omega$}
\label{censhift}
The observed X-ray image peak of merging clusters can be off from
their fitted center significantly \citep[e.g.,][]{mfg93,kbp+01}, while
the deviation is usually inconspicuous for relaxed
clusters. \citet{pfb+06} defined the centroid shift $\omega$ as the
standard deviation of the projected separation between the X-ray
brightness peak and the model fitted center, which is computed in a
series of circular aperture centered on the X-ray brightness peak from
$0.05~R_{\rm ap}$ to $R_{\rm ap}$ in steps of $0.05~R_{\rm ap}$
\citep[see also][]{omb+06}, thus
\begin{equation}
\omega=[\frac{1}{n-1}\sum\limits_{i}(\Delta_{i}-\langle \Delta
  \rangle)^2]^{\frac{1}{2}}\times \frac{1}{R_{\rm ap}}.
\label{w}
\end{equation}
Here $R_{\rm ap}=500~\rm{kpc}$ and $n=20$, $\Delta_{i}$ is the
distance between the X-ray brightness peak and the fitted model center
of the $i$th aperture, $\langle\Delta\rangle$ is the mean value of all
$\Delta_{i}$ \citep[e.g.,][]{ceg+10,ceb+13}. We calculate in this work
the centroid shift for clusters in central region with radius of 500
kpc, see Equation~\ref{w} and values are listed in Table 1.

Previously, the centroid shift of clusters has been calculated by
several authors with the similar formula, but for different
regions. In Figure~\ref{wcorr}, the values of centroid shift $\omega$
that we obtained are compared with those calculated in literatures.
\citet{ceg+10} derived the centroid shift for 32 clusters in the
aperture with a radius of 500 kpc, and 7 new clusters later in
\citet{ceb+13}. Good correlation is found between our values and
parameters in \citet{ceg+10,ceb+13}, as shown in Figure~\ref{wcorr}a.
\citet{wbsa13} calculated the centroid shift for 80 clusters by using
the {\it XMM-Newton} images in the annulus region of 0.1--1
$R_{500}$. \citet{wbc13} worked on a larger sample of clusters in
$R_{500}$ which contains 126 clusters based on {\it Chandra} and {\it
  XMM-Newton} data. We get good correlations, with large scatter,
between our values and parameters for 57 clusters in \citet{wbsa13} as
shown in Figure~\ref{wcorr}b, and for 101 clusters in \citet{wbc13} in
Figure~\ref{wcorr}c. \citet{der+16} calculated the centroid shift for
25 clusters in 500 kpc, and for 19 clusters in 0.5 $R_{500}$ with the
{\it Chandra} images. We find clear correlations between our values
and those in \citet{der+16}, as shown in Figure~\ref{wcorr}d and
Figure~\ref{wcorr}e, with large uncertainties. We also calculate the
centroid shift for 19 clusters in \citet{der+16} within the radius of
$0.5~R_{500}$, and compare our values to their results in
Figure~\ref{wcorr}f, and find the data are around the equivalent line
with a better correlation.

In Figure~\ref{wcorr}a and \ref{wcorr}d, the centroid shift are
calculated within the same radius of 500 kpc, but our results seem to
be systematically smaller than those obtained by \citet{ceg+10,ceb+13}
and \citet{der+16}, especially for relaxed clusters. This is mainly
caused by the smoothed images we use. For example, Z2701 and RXC
J1504.1-0248 (marked in Figure~\ref{wcorr}a) have the largest
deviations to the equivalent line. The logarithmic value of the
centroid shift for Z2701 is equal to -4.37 when the image is smoothed
to 30 kpc, but it reduces to -3.47 when the image is smoothed to 10
kpc, and to -2.50 for unsmoothed image which matches the value of
-2.44 obtained by \citet{ceg+10} very well. The logarithm of the
centroid shift for RXC J1504.1-0248 is equal to -4.00 when the image
is smoothed to 30 kpc, but it is reduced to -3.31 when the image is
smoothed to 10 kpc, and to -2.06 for unsmoothed image which is
slightly smaller than the -2.33 obtained by \citet{ceg+10}. The smooth
scale can therefore affect the value of centroid shift
significantly. Considering that the pixel size stands for different
physical scales for clusters with different redshifts, we suggest that
the X-ray images of clusters should be smoothed to a certain physical
size such as the 30 kpc we used.

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=0.23\textwidth]{fig5/ceg+10_p.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig5/wbsa13_pa.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig5/wbsa13_pb.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig5/wbc13_p.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig5/der+16a_p.eps}
\includegraphics[angle=0,width=0.23\textwidth]{fig5/der+16b_p.eps}\\[2mm]
\includegraphics[angle=0,width=0.23\textwidth]{fig5/der+16b_p1.eps}
\caption{Comparison for the power ratios that we calculated with those
  in literatures. The solid line in the panels a and e indicates
  equivalent values in X and Y axes if the radius is set in the same
  value of 500 kpc, the dotted line in panel g is also an equivalent
  line but for the radius of $0.5~R_{500}$. The correlation parameter
  $r_{\rm s}$ and $p_{\rm s}$ are labelled in the right-bottom corner
  of each panel. Data in the X-axis are our values and in the Y-axis
  are taken from: C10=\citet{ceg+10}, W13a=\citet{wbsa13},
  W13b=\citet{wbc13} and D16=\citet{der+16}.}
\label{pcorr}
\end{center}
\end{figure}

\subsubsection{Power ratio, $P_3/P_0$}
\label{powratio}
Because disturbed clusters generally have more remarkable fluctuations
of surface brightness than relaxed clusters, \citet{bt95} defined the
power ratio as dimensionless morphological parameters from the
two-dimensional multipole expansion of the projected gravitational
potential of clusters within $R_{\rm ap}=500\rm kpc$. The moments,
$P_{m}$, are defined as follows:
\begin{equation}
P_{0}=[a_{0}\ln (R_{\rm ap})]^2,
\label{p0}
\end{equation}
\begin{equation}
P_{m}=\frac{1}{2m^{2}R_{\rm ap}^{2m}}(a_{m}^2+b_{m}^2).
\label{pm}
\end{equation}
The moments $a_{m}$ and $b_{m}$ are calculated using
\begin{equation}
a_{m}=\int_{r\le R_{\rm ap}}f_{\rm obs}(x_i,y_i)(r)^{m}\cos(m\theta)dx_idy_i,
\label{am}
\end{equation}
and
\begin{equation}
b_{m}=\int_{r\le R_{\rm ap}}f_{\rm obs}(x_i,y_i)(r)^{m}\sin(m\theta)dx_idy_i,
\label{bm}
\end{equation}
where $f_{\rm obs}(x_i,y_i)$ and $\theta$ have the same meaning as
before.  $P_{3}/P_{0}$ is the power ratio, which was found to be
related to substructures
\citep[e.g.,][]{bfs+05,bpa+10,ceg+10,lfj+17,cds+18}.  Following the
custom, we take the radius at 500 kpc to calculate the power ratio,
see Equation~\ref{p0}-\ref{bm}, and list the results in Table 1.

The values of power ratio that we obtained are compared with those
calculated in literatures in Figure~\ref{pcorr}. Our values show good
consistence to those obtained by \citet{ceg+10} for 32 clusters in
central region with radius of 500 kpc, as shown in
Figure~\ref{pcorr}a. \citet{wbsa13} worked out power ratios for 80
clusters by using the {\it XMM-Newton} images in $R_{500}$, and also
defined $(P_3/P_0)_{\rm max}$ as the maximum value in different
annuluses along the radial direction \citep[see details
  in][]{wbsa13}. We get 57 clusters in common and find good
correlations between our results and those in \citet{wbsa13}, see
Figure~\ref{pcorr}b for $P_3/P_0$ and Figure~\ref{pcorr}c for
$(P_3/P_0)_{\rm max}$. By combining the {\it XMM-Newton} and {\it
  Chandra} data, \citet{wbc13} calculated power ratios for 126 galaxy
clusters in $R_{500}$. The Figure~\ref{pcorr}d shows a clear
correlation between values calculated by us and those in \citet{wbc13}
for 101 clusters in common. \citet{der+16} derived power ratios for 25
clusters in 500 kpc and 19 clusters in $0.5~R_{500}$ from the {\it
  Chandra} images. Good consistencies can be found between our results
and values for clusters from \citet{der+16} in Figure~\ref{pcorr}e and
Figure~\ref{pcorr}f. In Figure~\ref{pcorr}g, we show the correlation
between power ratios in $0.5~R_{500}$ calculated by us and those
obtained by \citet{der+16}.

In Figure~\ref{pcorr}a, \ref{pcorr}e and \ref{pcorr}g, the power
ratios in X-axis and Y-axis are calculated in the same radius, so that
the values obtained by us match well with those calculated by
\citet{ceg+10} and \citet{der+16} around the equivalent line. However,
we get slightly smaller power ratios for relaxed clusters than those
in \citet{ceg+10} and \citet{der+16} because we use the smoothed
images. For example, the logarithm of the power ratio of A267, which
has the largest deviation and labelled in Figure~\ref{pcorr}a, is
equal to -9.28$\pm$0.34 when the image is smoothed to 30 kpc, and it
is reduced to -7.24$\pm$0.48 when we use the unsmoothed image which
matches very well with the value of 7.75$\pm$0.83 obtained by
\citet{ceg+10}.

\subsection{A new morphological parameter for estimating the dynamical state}
In this subsection, we introduce and calculate a new morphological
parameter for estimating the dynamical state of galaxy cluster.

As seen in section \ref{knownpara}, the concentration index $c$, the
centroid shift $\omega$ and the power ratio $P_3/P_0$, which have been
widely used, are mostly calculated in the region of a fixed radius of
500 kpc or $R_{500}$ in recent papers. Obviously the fixed size to 500
kpc is not the best choice because clusters have various
sizes. Second, these three dynamical parameters are defined and
calculated in a circular region but the morphologies for most clusters
are elliptical or even irregular. Third, the redshift information is
needed but may not be available to define the 500 kpc or $R_{500}$,
especially for clusters identified in X-ray band or through SZ
effect. For clusters containing two or more subclusters, the central
region for calculating the dynamical parameters can only focus on the
main subcluster, which is not a proper indicator for the dynamical
state of the whole cluster (such as A115 in Figure~\ref{examples}).

The dynamical state of clusters should be derived in regions adaptive
to their size. The concerned region should be more reasonably taken as
being an elliptical region rather than a circular region. The image
edge of clusters in practice often is shown by the $\langle
S_{bg}\rangle+3\sigma$ contours, as the white ellipse in Figure
\ref{examples}. The dynamical state of clusters should be independent
to the exposure time of X-ray images. In the following, we use two
parameters to describe the global morphological properties of clusters
which are less influenced by the quality of X-ray images. We carry a
few steps toward the goal as following.

\subsubsection{Model fitting for the cluster shape}
\label{modelfit}
To describe the global morphology of a galaxy cluster, we fit the
X-ray surface brightness distribution with an elliptical 2-dimensional
$\beta$-model \citep{cf76,gad+14}:
\begin{equation}
f_{\rm mod}(x_i,y_i)=f_{\rm mod}(r)=A(1+(\frac{r}{r_0})^2)^{-\beta}+C,
\label{Ir}
\end{equation}
where
\begin{equation}
r(x_i,y_i)=\frac{\sqrt{x^2(1-\epsilon)^2+y^2}}{1-\epsilon},
\label{r}
\end{equation}
and
\begin{equation}
  \begin{split}
    x=(x_i-x_0)\cos\theta+(y_i-y_0)\sin\theta,\\
    y=(y_i-y_0)\cos\theta-(x_i-x_0)\sin\theta.
    \label{xy}
\end{split}
\end{equation}
Here $(x_0,y_0)$ are the coordinates of the center of the model, $A$
is the model amplitude, $r_0$ means the core radius, $\beta$ is the
power law index, and $C$ is a constant adjusting the count number for
the average background. $\epsilon$ is the ellipticity of cluster,
$\theta$ is the position angle of cluster, defined as the direction of
major axis from north to east. All these model parameters can be
determined by the least $\chi^2$ fitting.

We fitted the model to the X-ray image of 964 clusters of galaxies,
and obtained the model parameters which define the shape of galaxy
clusters. Among them, the most important is the index $\beta$, which
define steepness of the brightness profile. Others are not very
important for the dynamical state. After the model fitting, the
residual maps show more clearly the disturbed gas distribution. We
have tried some quantitative description of such disturbed state, and
found that the key factor is the asymmetry factor $\alpha$ (described
below), which in fact can be calculated free from the model fitting,
though the center has to be determined by the model-fitting.

\subsubsection{Two key morphology parameters derived from Chandra X-ray images}

By looking at hundreds of X-ray images of clusters and the residual
images, we find that relaxed clusters generally have steeper radial
distribution of brightness than disturbed ones. Further more, relaxed
clusters generally have roundish morphologies while merging clusters
could show more elongated shapes. Thus, we define a new profile
parameter $\kappa$ for the brightness distribution as being
\begin{equation}
\kappa=\frac{1+\epsilon}{\beta}.
\label{kappa}
\end{equation}
Here, $\beta$ and $\epsilon$ are the power law index and the
ellipticity for the fitted $\beta$-model.

Disturbed clusters usually are more significantly asymmetric than
relaxed ones \citep[e.g.,][]{wod88,ozf+10,zof+10,rme13,wh13}. Here we
use the asymmetry factor $\alpha$ to reveal the dynamical state of
galaxy clusters, which is defined as being
\begin{equation}
\alpha=\frac{\sum\limits_{x_i,y_i}{[f_{\rm
        obs}(x_i,y_i)-f_{\rm
        obs}(x_i',y_i')]^2}}{\sum\limits_{x_i,y_i}{f^2_{\rm
      obs}(x_i,y_i)}} \times100\;  \%,
\label{alpha}
\end{equation}
where $f_{\rm obs}(x_i',y_i')$ stands for observed flux at the
symmetry pixel of $(x_i,y_i)$ with respect to the cluster center
$(x_0,y_0)$.

\input tab2.tex

\begin{figure*}
  \begin{center}
\includegraphics[angle=0,width=0.75\textwidth]{ka.eps}\\
\caption{{\it Left panel:} Distribution for a test sample of 125
  clusters in the $\kappa-\alpha$ space. The relaxed clusters are
  presented as circles while disturbed clusters as crosses. The dash
  line stands for the best fitted line for all sources in the test
  sample, the green solid line is the adopted boundary that separates
  relaxed and disturbed clusters with highest success rate. {\it Right
    panel:} Similar to the left panel but for 964 clusters.}
\label{ka}
\end{center}
\end{figure*}

\label{combpara}
Because the dynamical state cannot be represented with a single
parameter and an easy threshold (see Figure~\ref{ka}), the combination
of the profile parameter $\kappa$ and the asymmetry factor $\alpha$
can indicate the dynamical state of galaxy clusters more
properly. Here we define the morphology index $\delta$ as being:
\begin{equation}
\delta=A\kappa+B\alpha+C.
\label{delta}
\end{equation}

To find the appropriate values of the coefficients $A$, $B$ and $C$,
we take a test sample of 125 clusters with known dynamical states
qualitatively classified as relaxed or disturbed by \citet{me12} who
built a statistically complete sample from their MAssive Cluster
Survey \citep[MACS:][]{eeh01,ebd+07,eem+10} and the Brightest Cluster
Sample \citep[BCS \& eBCS:][]{eeb+98,eea+00} based on data of the {\it
  ROSAT All-sky Survey} \citep[{\it RASS}:][]{t93}, with a flux larger
than $1\times10^{-12}\rm erg ~s^{-1}cm^{-2}$ in 0.1-2.4 keV and a
luminosity $L_{\rm RASS}>5\times10^{44}\rm erg~s^{-1}$. They collected
129 clusters and got X-ray and optical images for 108 clusters of
them, and dynamically classified these clusters into 4 groups. We
directly take the dynamical information of the 108 clusters from
\citet{me12} and define clusters in group 1 as relaxed clusters and
those in group 2-4 as disturbed clusters. In addition, we take 17
extra clusters from the BCS and eBCS sample \citep{eeb+98,eea+00} that
satisfy $f_{\rm RASS}>1\times10^{-12}\rm erg ~s^{-1}cm^{-2}$ and
$L_{\rm RASS}>5\times10^{44}\rm erg~s^{-1}$ and have been observed by
{\it Chandra} later, and dynamically classify them with the same
standard in \citet{me12}. In total, we have 125 clusters with known
dynamical state as listed in Table 2. This test sample is X-ray flux
complete and volume complete to $z\sim0.4$ \citep{me12}.

\begin{figure*}
\begin{center}
\includegraphics[angle=-90,width=0.45\textwidth]{2dcorr_sample.eps}
\hspace{.4in}
\includegraphics[angle=-90,width=0.45\textwidth]{2dcorr_all.eps}
\caption{{\it Left panel:} One-dimensional histogram distributions and
  two-dimensional correlations of four kinds of dynamical proxies for
  125 clusters with known dynamical states from literatures. Crosses
  and red histograms indicate disturbed clusters, while circles and
  blue histograms stand for relaxed clusters. The Spearman rank-order
  correlation coefficient $r_s$ and the relevant significance level of
  the correlation $p_s$ are labelled in each panel. A small value of
  $p_s$ indicates a significant correlation ($r_s$>0) or
  anti-correlation ($r_s$<0). {\it Right panel:} Similar to the left
  panel but for the sample of 964 clusters.
\label{corrsct}}
\end{center}
\end{figure*}

The $\kappa-\alpha$ distribution for these 125 clusters are plotted in
the left panel of Figure~\ref{ka}. Relaxed clusters are separated from
disturbed ones. First, in the $\kappa-\alpha$ space we find the
best-fitted line as
\begin{equation}
\rm{log_{10}}(\alpha)=0.94\kappa-2.43.
\label{kafit}
\end{equation}
The line of demarcation has the form of
\begin{equation}
{\rm log_{10}}\alpha=A_1\kappa+B_1,
\label{demarc}
\end{equation}
which is perpendicular to the best-fitted line with
$A_1=-1/0.94=-1.07$ and $B_1=-0.31$, and has the highest success rate
to discriminate the relaxed and disturbed clusters. For the test
sample the success rate is 115/125=88\%. The morphology index $\delta$
of galaxy clusters is therefore defined as the distance to the best
demarcation line in $\kappa-\alpha$ space, that is
\begin{equation}
\begin{aligned}
\delta &=0.68\rm{log_{10}}(\alpha)+0.73\kappa+0.21,
\label{delta}
\end{aligned}
\end{equation}
which can quantitatively indicate the dynamical state of galaxy
clusters just from the X-ray image without the redshift
information. We calculate the two parameters $\kappa$ and $\alpha$ and
hence the morphology index $\delta$ for all 964 clusters, as listed in
Table 1.

\section{Comparison and applications of the dynamical parameters}

\begin{figure}
\begin{center}
\includegraphics[angle=0,width=0.38\textwidth]{wh13_delta.eps}
\caption{Correlations between the X-ray morphology index, $\delta$,
  which we obtained in this paper and optical dynamical parameter,
  $\Gamma$, obtained by \citet{wh13} for 190 clusters. The Spearman
  rank-order correlation coefficient $r_s$ and the correlation
  significance $p_s$ are marked in the bottom left corner.}
\label{gcorr}
\end{center}
\end{figure}

\begin{figure*}
\begin{center}
\includegraphics[angle=0,width=0.65\textwidth]{radio.eps}
\caption{Applications of the morphology index $\delta$ to the scaling
  relations for 35 clusters with giant-halos (open circles) and 12
  clusters with mini-halos (solid points).  The correlation is
  enhanced and the scatter is reduced when the dynamical parameter
  $\delta$ is involved ({\it right panel}) compared to the case
  without involvement ({\it left panel}). The dashed line stands for
  the best fitted line. The intrinsic scatter $\sigma$, the Spearman
  rank-order correlation coefficient $r_s$ and the correlation
  significance $p_s$ are marked in the bottom right corner.}
\label{halo}
\end{center}
\end{figure*}

\subsection{Comparison for the four dynamical parameters}
\label{paracorr}
Here we discuss the correlations among the four X-ray dynamical
parameters, the concentration index, the centroid shift, the power
ratio and the morphology index. To reflect the goodness of these
parameters more clearly, the test sample and the whole sample of
clusters are plotted and examined separately in Figure~\ref{corrsct}.

The left panels of Figure~\ref{corrsct} show that relaxed (blue) and
disturbed (red) clusters are separated in four parameters, indicating
that these four parameters are sensitive diagnostics for dynamical
state. However, as shown in the histograms they all are continuously
distributed without a clear boundary for the two states, which is
reasonable in practice. If one has to set a criterion to separate them
and get the fraction of disturbed clusters, that is 48.8\% if the
criterion is taken as $\delta>0$, close to the value of $\sim$50\%
obtained by previous works based on X-ray data
\citep[][]{bfs+05,crb+07,hmr+10}. In 2D parameter diagrams, clear
correlations are shown between each pair of two dynamical proxies,
though some of them have large scatter. The small $p_{s}$ in each
panel means the correlations are significant. In the right panels,
these dynamical parameters for the whole sample of 964 clusters are
found to be well correlated, and also distributed continuously from
very disturbed state to the relaxed state with a peak in
between. Strongest correlations appear between ${\rm log}_{10}(c)$,
$\delta$ and ${\rm log}_{10}(\omega)$, which implies that the three
parameters are similarly sensitive to the dynamical state.

\subsection{Correlation between the morphology index derived from X-ray image and the relaxation factor derived from the optical data}

Based on the projected distribution of member galaxies, \citet{wh13}
derived dynamical parameters for 2092 rich clusters. The optical image
of clusters is smoothed to 20 kpc with a Gaussian kernel and weighted
by luminosities of member galaxies. They got the relaxation factor
$\Gamma$ by combining the asymmetry factor, the ridge flatness and the
normalized deviation based on a test cluster sample with known
dynamical states from literatures. We get 190 clusters by cross
matching the optical catalogue in \citet{wh13} with our X-ray
sample. The Figure \ref{gcorr} shows the distribution of the optical
relaxation factor $\Gamma$ and the morphology index $\delta$. It is
clear that relaxed clusters are both well recognized in the view of
optical and X-ray dynamical parameters as the data are concentrated in
the upper-left corner. However, parameters are very scattered for
disturbed clusters, which means that both hot gas and member galaxies
can be very disordered but are not necessary to move or
distribute coincidentally.

\begin{figure*}
\begin{center}
\includegraphics[angle=-90,width=0.95\textwidth]{mass_offset.eps}\\[3mm]
\includegraphics[angle=-90,width=0.98\textwidth]{mass_scaling.eps}\\[3mm]
\includegraphics[angle=-90,width=0.98\textwidth]{mass_scaling_lx.eps}\\
\caption{{\it Top panels:} Relations between dynamical parameters and
  ratio of mass estimated from X-ray and SZ-effect data for 316
  clusters. The solid line is for the equality of the two estimated
  masses. {\it Middle panels:} Influence of dynamical parameters on
  the $M_{\rm X,500}-M_{\rm SZ,500}$ scaling relations for 316
  clusters. The dashed line is the best fitted line, and the intrinsic
  scatter $\sigma$, the Spearman rank-order correlation coefficient
  $r_s$ and the correlation significance $p_s$ are marked in the
  bottom right corner. {\it Bottom panels:} Similar to the middle
  panels but for the $L_{\rm X,500}-M_{\rm SZ,500}$ relation.  The
  25\% most relaxed clusters ($\delta<-0.15$) are denoted as black
  circles.}
\label{Moff}
\end{center}
\end{figure*}

\subsection{The scaling relations and the fundamental plane for radio giant-halos and mini-halos in clusters}

\begin{figure*}
\begin{center}
\includegraphics[angle=-90,width=0.95\textwidth]{mass_offset1.eps}\\[3mm]
\includegraphics[angle=-90,width=0.98\textwidth]{mass_scaling1.eps}\\
\caption{Similar to Figure~\ref{Moff} but for masses estimated through
  weak lensing analyses.}
\label{Moff-wl}
\end{center}
\end{figure*}

Large scale diffuse radio sources in clusters are associated with ICM
rather than member galaxies, which can be classified into giant-halos,
relics and mini-halos \citep[see][as an observational
  review]{fgg+12}. It has been widely believed that the formation of
diffuse radio sources are related to merging process of their host
clusters. Giant-halos and relics are usually discovered in disturbed
clusters, while mini-halos are detected around the cool core of
relaxed clusters. Giant-halos and mini-halos may related to the
(re-)acceleration process by turbulence in the ICM, while relics are
related to the shock resulted by major merger \citep[see][for
  detail]{bj14}. The scaling relation has been established between
radio power of giant-halos and mini-halos and the mass or mass proxies
of the host clusters \citep[e.g.,][]{lhb+00,bcds09,ceb+13,yhw15}.

We collected the radio powers $P_{1.4~\rm GHz}$ in
$10^{24}\rm~W~Hz^{-1}$, the X-ray luminosities $L_{500}$ within
$R_{500}$ in $10^{44}\rm~erg~s^{-1}$ from \citet[][see the references
  therein]{yhw15}. Now we check if the morphology index $\delta$ can
reduce the scatter of data and enhance the correlation for the scaling
relation. In the left panel of Figure~\ref{halo}, the radio power of
giant-halos (red open circles) for 35 clusters and mini-halos (blue
solid points) for 12 clusters are plotted against X-ray luminosity of
host clusters. The mini-halos have statistically less radio power than
giant-halos but follow the same scaling relation with the X-ray
luminosity. We fit all these data with weight of data errors, see the
Appendix in \citet{yhw15}. The scaling relation has a Spearman
rank-order correlation coefficient $r_s=0.52$. In the right panel, the
morphology index is involved and the scaling relation show a larger
correlation index ($r_{\rm s}=0.68$) and a reduced scatter
($\sigma=0.25$), which means that giant-halos and mini-halos follow
the similar fundamental plane between the radio power, cluster mass
and the morphology index.

\subsection{Influence of dynamical state on mass estimation for clusters}

The gravitational mass of clusters is an essential parameter for a lot
of researches, which is however affected by the dynamical state if
they are estimated from X-ray images. \citet{pap+11} published a
meta-catalogue contains 1,743 clusters based on previous works from
data of {\it ROSAT} All Sky Survey. The masses $M_{\rm X,500}$ and
characteristic radius $R_{500}$ of clusters are estimated
homogeneously by using X-ray data. \citet{paa+16} detected 1,653
clusters and estimated masses $M_{\rm SZ, 500}$ from SZ-effect maps,
which is insensitive to dynamical states of clusters
\citep[e.g.,][]{mhbn05}. Here we use a sample of 316 clusters commonly
in the catalogues in \citet{pap+11} and \citet{paa+16} to check the
influence of dynamical state on the mass estimation for clusters.

In the top panels of Figure~\ref{Moff}, we show the distribution of
four dynamical parameters against the ratio of masses estimated from
X-ray and SZ-effect data, and find a clear dependence on the
concentration index $c$, the centroid shift $\omega$ and the
morphology index $\delta$, though the correlation is weaker for the
power ratio $P_3/P_0$. The mass estimated from the X-ray images are
underestimated comparing to masses derived from SZ-effect for
disturbed clusters, but apparently overestimated for relaxed
clusters. The SZ-effect reflects the thermal and non-thermal
components of clusters simultaneously. The underestimation of mass
through X-ray data for disturbed clusters may indicate that disturbed
clusters have higher fraction of energy stored as non-thermal form. In
the middle rank of Figure~\ref{Moff}, the influence of dynamical
parameters is considered for the correction of the mass estimation
from the X-ray images on the $M_{\rm X,500}-M_{\rm SZ,500}$
relation. The Spearman rank-order correlation coefficient increase
when dynamical parameters are considered.

On the other hand, the X-ray luminosity $L_{\rm X,500}$ is widely used
as a mass proxy of clusters. In the bottom panels of
Figure~\ref{Moff}, we show the $L_{\rm X,500}-M_{\rm SZ,500}$
relations for the 316 clusters. Since relaxed and disturbed clusters
usually show different mass scaling relations
\citep[e.g.,][]{omb+06,ldc+06,ldk+09,zjc+13}, we compare the 25\% most
relaxed clusters ($\delta<-0.15$, black circles) with the full
sample. For the 2D correlation (the most left panel), the relaxed
clusters are clearly offset from the best fitted line determined by
the total sample. The offset is significantly reduced when the
dynamical parameters included (the right four panels), which again
means that there is a fundamental plane for the mass, luminosity and
dynamical state for galaxy clusters. The 3D correlations have smaller
intrinsic scatter and larger $r_{\rm s}$ than the 2D correlation.

For the further demonstration, we investigate the influence of
dynamical state to the lensing mass of galaxy clusters which is
believed to be the most reliable. \citet{s15} compiled a catalogue
which contains 485 clusters with mass ($M_{\rm WL,500}$) estimated
through weak lensing analyses. We cross-match our samples with the
catalogues in \citet{s15} and \citet{pap+11}, and get 151 clusters in
common. In Figure~\ref{Moff-wl}, we do similar correlations to
Figure~\ref{Moff} but replace $M_{\rm X,500}$ with $M_{\rm
  WL,500}$. In the upper panels of Figure~\ref{Moff-wl}, the mass
ratio ($M_{\rm X,500}/M_{\rm WL,500}$) is not randomly around but
mostly less than 0, which means that the X-ray mass is
underestimated. In the lower panels, the correlations are enhanced
with slightly larger $r_s$ when dynamical parameters are included.

\section{Summary}
The {\it Chandra} satellite has accumulated X-ray image data for
$\sim$1000 clusters of galaxies, which make it feasible to work out a
catalogue of dynamical parameters for a large sample of clusters.

We collected Chandra data for 964 galaxy clusters. The images for
these clusters are processed following the same procedure, and
smoothed by a Gaussian function with a size of 30 kpc (or 5 arcseconds
for few clusters with no redshift available). Three widely used
dynamical parameters, i.e., the concentration index $c$, the centroid
shift $\omega$ and the power ratio $P_3/P_0$, are calculated from the
X-ray images in a circular region with a radius of 500 kpc. Our
results are consistent with the values available in literatures. We
also derive two adaptive parameters, the profile parameter $\kappa$
and the asymmetry factor $\alpha$ in the best fitted elliptical
region, and define the morphology index $\delta$, which can be
excellent indicator for dynamical state. Finally, we calculate the
profile parameter $\kappa$, the asymmetry parameter $\alpha$ and the
morphology index $\delta$ for 964 clusters (see Table 1).

The four dynamical parameters are well correlated, and can separate
relaxed clusters from disturbed ones. But they are distributed
continuously from the most disturbed status to the relaxed. The
morphology index we derived from the X-ray images and optical
relaxation factor is also significant correlated.

By involving the dynamical parameters, we find that clusters with
diffuse radio giant-halos and mini-halos follow the same fundamental
plane between radio power, cluster mass and morphology index. Mass
estimation for cluster of galaxies if derived from X-ray data is
affected by their dynamical states.

\section*{Acknowledgements}

We thank the referee and Dr. Z. L. Wen for helpful comments. The
authors are supported by the National Natural Science Foundation of
China (11803046, 11988101 and 11833009) and the Young Researcher Grant
of National Astronomical Observatories, Chinese Academy Sciences.
This research has made use of data obtained from the {\it
  Chandra Data Archive} and the {\it Chandra Source Catalog}, and
software provided by the {\it Chandra X-ray Center (CXC)} in the
application packages {\it CIAO}, {\it ChIPS}, and {\it Sherpa}.

\section*{Data availability}

The data underlying this article, the X-ray images for 964 clusters,
and also the code for calculations are available on the web-page
http://zmtt.bao.ac.cn/galaxy\_clusters/dyXimages/.

\bibliographystyle{mnras}


\section{Introduction}
The concept of correlation and regression was originally conceived by Galton when studying how strongly the characteristics of one generation of living things manifested in the following generation \cite{stanton2001galton}. The ideas prompting the development of more mathematically rigorous treatment of correlation
were developed by Karl Pearson in 1896, which yielded the well-known Pearson Product Moment Correlation Coefficient \cite{pearson1896vii} given as
\begin{equation}
\rho(X, Y) = \frac{Cov(X,Y)}{\sigma_X\sigma_Y} = \frac{E(X - \mu_X)(Y - \mu_Y)}{\sqrt{E(X - \mu_X)^2E(Y - \mu_Y)^2}},
\end{equation}
where $X$ and $Y$ are two random variables from a non-degenerative joint distribution $F_{XY}$, $Cov(X,Y)$ denotes the covariance, and $\mu_X$ and $\mu_Y$, $\sigma_X$ and $\sigma_Y$ are the population means and standard deviations, respectively. If we let $(X_1,Y_1)$, $(X_2,Y_2)$, …, $(X_n,Y_n)$ denote $n$ paired i.i.d. observations, then the sample Pearson correlation coefficient is given as
\begin{equation}
\label{eq:sample_pearson}
r(X, Y) = \frac{ \Sigma_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y}) }{\sqrt{\strut\Sigma_{i=1}^n(X_i-\bar{X})^2}\sqrt{\strut\Sigma_{i=1}^n(Y_i-\bar{Y})^2}}.
\end{equation}

Shortly thereafter Pearson's work was published  Spearman introduced the rank correlation coefficient in 1904, with the advantages of being robust to extreme values and disparities of distributions between two variables \cite{spearman1961proof}.  It should be noted however that
K. Pearson, in his biography of Galton, says that the latter "dealt with the correlation of ranks before he even reached the correlation of variates, i.e. about 1875", but Galton apparently published nothing explicitly \cite{kendall1979advanced}.  Mathematically, Spearman’s correlation coefficient is defined as the Pearson correlation coefficient on the ranks of $(X_1,Y_1)$, $(X_2,Y_2)$, …, $(X_n,Y_n)$ and denoted as
\begin{equation}
\rho_s (X,Y)=\rho(a,b),
\end{equation}
where $a_i= Rank(X_i)$ and $b_i=Rank(Y_i)$ are the ranks of $X_i$ and $Y_i$, respectively, $i=1,2, \cdots, n$.
In general, when discussing the Spearman correlation coefficient little attention is given to its population measure.
However,  if we consider that   $a_i/n$ converges to $F(X_{i})$, then one may consider the population measure linked to Spearman's sample correlation coefficient as,
\begin{equation}
\rho_s (X,Y)=\frac{E\{F_X(X)-E[F_X(X)]\}\{F_Y(Y)-E[F_Y(Y)] \}}{E\{F_X(X)-E[F_X(X)]\}^2 E\{F_Y(Y)-E[F_Y(Y)]\}^2}
=\frac{E[F_X(X)-\frac{1}{2}][F_Y(Y)-\frac{1}{2}] }{E[F_X(X)-\frac{1}{2}]^2 E[F_Y(Y)-\frac{1}{2}]^2},
\end{equation}
where $F_X$ and $F_Y$ are the marginal cumulative distribution functions (CDFs) for $X$ and $Y$, respectively.
The sample estimator of $\rho_s$ can be obtained by replacing the original observations with their ranks in Equation \ref{eq:sample_pearson},
\begin{equation}
\label{eqn:spearman}
r_s(X, Y) = r(a, b) = \frac{ \Sigma_{i=1}^n(a_i-\bar{a})(b_i-\bar{b}) }{\sqrt{\Sigma_{i=1}^n(a_i-\bar{a})^2}\sqrt{\Sigma_{i=1}^n(b_i-\bar{b})^2}} = 1-\frac{6\Sigma_{i=1}^n(a_i - b_i)^2}{n(n^2-1)}.
\end{equation}
For samples from a bivariate normal population, there is also a known relation between the Spearman and Pearson correlation coefficients \cite{moran1948rank}, which is
\begin{equation}
E(r_s) = \frac{6}{\pi (n+1)}\{\sin^{-1}\rho + (n-2)\sin^{-1}\frac{1}{2}\rho\}.
\end{equation}
While
Pearson's $\rho$ measures the linear relationship between two random variables it is often described that Spearman's $\rho_s$ measures a monotonic association
between $X$ and $Y$, thus it may be considered more general measure of association, albeit
it does measure the linear association between $F(X)$ and $F(Y)$. Spearman's correlation coefficient is also less sensitive to extreme values because it is rank based. Due to these advantages, it is widely used as a measure of association between two measurements. It is often of interest to test whether two random variables are correlated, i.e. $H_0:\rho_s=0$, for which the common methods include  a $t$-distribution based test, incorrectly based on the approximate
bivariate normality of the ranks, a test based on Fisher's $Z$ transformation, again assuming approximate bivariate normality of the ranks,
and what we term the naive permutation test. The $t$-distribution based test  is commonly used when the sample size is large, with the $t$-statistic defined as
$$t=r_s \sqrt{\frac{n-2}{1-r_s^2} }.$$
Under bivariate normality assumptions this statistic approximately follows student’s $t$ distribution with $n-2$ degrees of freedom under $H_0$. For the test based on Fisher’s $Z$ transformation, the statistic is defined as

$$ Z = \frac{1}{2}\arctan(r_s) = \frac{1}{2}\ln \frac{1+r_s}{1-r_s}. $$

Under bivariate normality assumptions the transformed $Z$ statistic approximately follows normal distribution $N(0, \frac{1.06}{n-3})$ under $H_0$. For small sample size scenarios, naive permutation tests are also often used, where $X$ and $Y$ are randomly shuffled separately to simulate the sample distribution of $r_s$ under $H_0$, which is an exact test for testing the independence between $X$ and $Y$ under the null exchangeability assumptions, but may be
an invalid test for testing $H_0:\rho_s=0$ given $G(F(X),F(Y))= G(F(X))G(F(Y))$ does not imply $\rho_s=0$, where $G(\cdot,\cdot)$ denotes
the joint CDF of $F(X),F(Y)$.

These tests are so widely used that they are often the default options in common statistical software packages such as R \cite{team2013r} and SAS \cite{sas2015base}. However, there is little discussion that these tests relies on the untenable assumption that the underlying sample distribution of the ranks
 has a  bivariate normal distribution, which is in fact an impossibility. Even among those who noted this assumption, there is a misconception that the above tests are robust to such deviations because Spearman’s $\rho_s$ is rank based. This is exemplified in a discussion by Feller et al. in their article \cite{fieller1957tests},
\begin{quote}
Conversely, starting from any bivariate distribution $\phi(X,Y)$ we can always find monotonic transformations $X=f(x)$, $Y=g(y)$ to standardized normal variates $x$ and $y$. The resulting bivariate distribution $\psi(X,Y)$ will not necessarily be bivariate normal, but we think it likely that in practical stations it would not differ greatly from this form. This is a field in which further investigation would be of considerable interest.
\end{quote}
However, as we will show in Section \ref{sec:simulations}, all the commonly used tests about $\rho_s=0$ as discussed above, including
the naive permutation test, are not even asymptotically valid when the non-exchangeabilty assumptions are violated
under $H_0$. In some cases, the type I error can severely drift away from the desired level as the sample size increases! A undesirable feature
that is more notable in the era of ``big-data''.   Another variation of this approach is the
  Fisher-Yates coefficient, which transforms the original $X$ and $Y$ to their corresponding normal quantiles before the testing \cite{fisher1938statistical}. Although the marginal distributions are of the transformed variates take a pseudo normal form, the joint normality of these transformed values is not guaranteed.

In terms of our modified permutation test it is important to note the classic large sample result in Serfling where the “distribution free” large sample normal approximation for the sampling distribution for Pearson's sample correlation coefficient $r$
is derived using the multivariate delta method \cite{kendall1979advanced}. This method guarantees type I error converges to $\alpha$ when $n\to \infty$ given finite fourth moments. A straightforward way to obtain a similar result for Spearman's correlation is given by replacing the i.i.d sample with corresponding ranks. The test is asymptotically valid because the ranks are asymptotically independent, as we will discuss in Section \ref{sec:methods}.
Even though large sample approximations about these estimators are asymptotically valid they tend to suffer inflated type I errors in the small sample setting,
e.g. $n<50$.

To address this issue, we propose a studentized permutation test for  Spearman's correlation $\rho_s$, which extends the work or Diccicio and Romano for Pearson’s correlation coefficient
 $\rho$ \cite{diciccio2017robust}. We will show that the proposed test is asymptotically valid under general assumptions and is exact under
exchangeability assumptions when  $G(F(X),F(Y))= G(F(X))G(F(Y))$, i.e. more simply when  $X$ and $Y$
are independent.
We show that our newly proposed test
has robust Type I error  controls type I error control even when the sample distribution
is dependent (non-exchangeability)  and non-normal. Importantly, the type I error is well controlled when the sample size is small. This will be illustrated by a set of simulation studies. Finally, we will demonstrate the application of this test in real world examples of transcriptomic data of TCGA breast cancer patients, as well as a data set of PSA levels and age.

\section{Methods}
\label{sec:methods}

\subsection{Background information}

In this section, we start by reviewing the robust permutation test for Pearson's correlation coefficient as proposed by Diciccio and Romano\cite{diciccio2017robust}.
Towards this end, we define $\textbf{G}_n$ to be the set of all permutations $\pi$ of $\{1, \dots, n\}$. For testing independence between two random variables $X$ and $Y$, the permutation distribution of any given test statistic $T_n(X^n, Y^n)$ is defined as
\begin{eqnarray}
\label{perm_test}
\hat{R}_n^{T_n}(t) = \frac{1}{n!}\sum_{\pi \in G_n} I\{T_n(X^n, Y^n_{\pi})\le t\}
\end{eqnarray}
where $Y^n_{\pi}$ represents $\{Y_{\pi(1)}, \dots, Y_{\pi(n)}\}$. In this setting, the permutation $\textbf{G}_n$ is all possible pairwise combinations between $X^n$ and $Y^n$. A level $\alpha$ one-sided permutation test rejects if $T_n(X^n, Y^n_{\pi})$ is larger than the $1-\alpha$ quantile of the permutation distribution. The permutation test is exact when exchangeability assumptions hold, that is, the distribution of $ (X^n, Y^n)$ is invariant under the group of transformations $\textbf{G}_n$. The test using the Pearson correlation coefficient $\hat{\rho}$ is exact when used a metric of dependence for testing the null hypothesis of independence given as
$$ H_0: P=P_X \times P_Y, $$
where $P_X$ and $P_Y$ are marginal distributions of $X$ and $Y$, respectively. The null hypothesis of independence is not equivalent to the test about zero correlation given as $H_0: \rho=0$ with the exception of limiting assumptions such as the data are distributed as bivariate normal random variables. In other words, in the general setting two random variables can be dependent but uncorrelated. In such cases, DiCiccio and Romano \cite{diciccio2017robust} have shown that, with finite fourth moments, the permutation distribution of $\hat{\rho}$ converges to $N(0,1)$, but its sampling distribution converges to $N(0, \tau^2)$, where
$$\tau^2 = \frac{\mu_{22}}{\mu_{20}\mu_{02}},$$
and
$$\mu_{pq} = E[(X_1 - \mu_1)^p(Y_1 - \mu_2)^q].$$
Thus the test will not be level $\alpha$ unless $\tau = 1$. In light of this result, DiCiccio and Romano proposed a studentized correlation test statistic, which has been shown to control Type I error asymptotically at $\alpha$ when two random variables are dependent but uncorrelated \cite{diciccio2017robust}. Specifically, the studentized statistic is defined as $S_n = \sqrt{n}\hat{\rho}_n/\hat{\tau}_n$, where
$$\hat{\tau}_n^2 = \frac{\hat{\mu}_{22}}{\hat{\mu}_{20}\hat{\mu}_{02}},$$
$$\hat{\mu}_{pq} = \frac{1}{n} \sum^n_{i=1}(X_{i} - \bar{X})^p(Y_{i} - \bar{Y})^q.$$
The permutation distribution and sampling distribution of $S_n$ both converge to the standard normal distribution asymptotically. It should be noted that even though the results presented in DiCiccio and Romano\cite{diciccio2017robust} are based on large sample approximations, the behavior of this test for small to moderate sample sizes is quite good as born out in their simulation results.

\subsection{ Spearman's permutation correlation test}

 Spearman's coefficient $\rho_s(X, Y)$ is the Pearson correlation coefficient of the ranks of $X$ and $Y$, that is $\rho(a, b)$. When there are ties in the data, their ranks are typically taken as an  average. Unlike Pearson's correlation coefficient $\rho$, which measures the linear relationship between two random variables,  Spearman's correlation coefficient $\rho_s$ measures a monotonic association, thus is far less restrictive. Note that Spearman's correlation coefficient is
also the linear measure between $F(X)$ and $F(Y)$.
 It is also less sensitive to non-normality or extreme values.

Despite the above advantages, it is a misconception that the tests of  $\rho_s=0$ based on the bivariate normality assumptions underlying the
original data
will be robust to the deviation from this assumption. In fact, when the original data are independent, their ranks will will be dependent. Thus, tests of $\rho_s$ typically suffer similar issue as for $\rho$.  We also emphasize that "normality" refers to the \emph{joint} normality
as  opposed to marginal normality, because two random variables that are marginally normal can have a joint non-normal distribution. Therefore, the Fisher-Yates coefficient, which back transforms a variables rank through the normal quantile function does not provide what heuristically one may
consider as a simple correction. In Section \ref{sec:simulations}, we will empirically show that violation of the joint normality assumption will have
severe effect on type I error control. In addition, it is in fact impossible for the joint distribution of the ranks to be bivariate normal.

Our approach is to replace $(X_i, Y_i)$ with their ranks $(a_i, b_i)$ in order to develop a Spearman's correlation permutation test analog to the
Pearson's correlation permutation test, with some subtle differences. The studentized permutation test of the Pearson's $\rho$ only requires finite fourth moments and that observations are i.i.d. Although the pairs of ranks $(a_i, b_i)$ are no longer i.i.d observations, we can show that they asymptotically satisfy this condition.

When $(X_i, Y_i)$ are from paired i.i.d. observations, we have $(F_X(X_i), F_Y(Y_i))$ being i.i.d. as well. Since $a_i=Rank(X_i)=n\hat{F}_X(X_i)$ and $\hat{F}_X(X_i) \to F_X(X_i)$ as $n\to \infty$, we have $\frac{a_i}{n} \to F_X(X_i)$. Therefore, the paired observations $(a_i, b_i)$ are asymptotically i.i.d. Consequently, the exchangeability condition will hold at least asymptotically and the test will be asymptotically exact. Intuitively, when $n$ is sufficiently large, knowing $(a_i, b_i)$ will lend little knowledge on the ranks of another pair $(a_j, b_j)$.  It can also be shown that the correlation between $a_i$ and $a_j (i\ne j)$ is approximately $-\frac{1}{n-1}$. Therefore, for a sample sequence $X^n$, the correlation matrix of the ranks converges to $I_n$ when $n \to \infty$.

Specifically, the one sided studentized permutation test for testing $H_0: \rho_s = 0$ versus $H_1: \rho_s >0$ is performed by the following steps. The test is implemented in the R \verb|perk| (\textbf{per}mutation tests of correlation c\textbf{(k)}oefficients) package, which will be available on CRAN (The Comprehensive R Archive Network, https://cran.r-project.org/) and GitHub (https://github.com/hyu-ub/perk).

\begin{itemize}
	\item For $n$ paired i.i.d. observations $(X_1,Y_1)$, $(X_2,Y_2)$, …, $(X_n,Y_n)$, calculate their ranks within each random variable, $(a_1,b_1)$, $(a_2,b_2)$, …, $(a_n,b_n)$.
	\item Estimate the Spearman's $\rho_s$ using Equation \ref{eqn:spearman} as $r_s$.
	\item Estimate the variance of sample estimates $r_s$ by
	$$\hat{\tau}_n^2 = \frac{\hat{\mu}_{22}}{\hat{\mu}_{20}\hat{\mu}_{02}},$$
	$$\hat{\mu}_{pq} = \frac{1}{n} \sum^n_{i=1}(a_{i} - \bar{a})^p(b_{i} - \bar{b})^q.$$
	\item Calculate the studentized statistic $R_s = r_s/\tau_n$.
	\item Randomly shuffle $(b_1, b_2, ... b_n)$ for $B$ times. For each permutation, calculate the permuted studentized statistic $R_s^k$, $k \in (1, ..., B)$.
	\item Calculate the p-value by $$p=\frac{1}{B}\Sigma^B_{k=1}I(R_s^k>R_s).$$
	\item Reject $H_0$ if $p<\alpha$.
\end{itemize}

\section{Simulations}
\label{sec:simulations}

We examined the Type I error control  across all of the tests introduced above
using distributions commonly found in the literature for these examinations across a wide range of settings \cite{diciccio2017robust, hutson2019robust}. For our simulation study, we focused on testing $H_0: \rho_c = 0$ versus $H_1: \rho_c > 0$, with sample sizes $n = 10, 25, 50, 100, 200$. Each simulation utilized $10,000$ Monte Carlo replications and the number of permutations used is $1,000$. We compared the $t$ test, Fisher's $Z$-transformation (Fisher's $Z$), Fisher-Yates method, Serfling's large sample normal approximation (Asymp Norm), naive permutation test (Permute), and studentized permutation test (Stu Permute). The Type I error control for $\alpha=0.05$ was examined. The simulation scenarios 1 through 5 from DiCiccio and Romano. Two additional distributions were studied as well:

\begin{itemize}
	\item[1.] Multivariate normal (MVN) with mean zero and identity covariance.
	\item[2.] Exponential given as $ (X, Y) = rS^Tu $ where $S = diag(\sqrt{2}, 1)$, $r\sim \exp(1)$, and $u$ is uniformly distributed on the two dimensional unit circle.
	\item[3.] Circular given as the uniform distribution on a two dimensional unit circle.
	\item[4.] $t_{4.1}$ where $X=W+Z$ and $Y=W-Z$, where $W$ and $Z$ are iid $t_{4.1}$ random variables.
	\item[5.] Multivariate $t$-distribution (MVT) with 5 degrees of freedom.
	\item[6.] Mixture of two bivariate normal distributions given as $(X, Y) = WZ_1 + (1-W)Z_2$ where $W \sim Bernoulli(0.5)$, $Z_1 \sim N( \begin{pmatrix} 0 \\0 \end{pmatrix}, \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}) $, $Z_2 \sim N( \begin{pmatrix} 0 \\0 \end{pmatrix}, \begin{pmatrix} 1 & -\rho \\ -\rho & 1 \end{pmatrix}) $. We select a range of $\rho$'s: 0.1, 0.3, 0.6 and 0.9 to simulate different degrees of dependencies between $X$ and $Y$ (MVN 1, MVN 3, MVN6, MVN 9).
	\item[7.] Mixture of four bivariate normal distributions (MVN 45), given as $(X, Y) = \Sigma_{k=1}^4 I(W_i=k)Z_i$ where $P(W=k) = 0.25$, $k=1,2,3,4$. In addition, $Z_k \sim N(\mu_k, I_2) $, where $\mu_1 = (5, 5)^T$, $\mu_2 = (5, -5)^T$, $\mu_3 = (-5, 5)^T$, $\mu_4 = (-5, -5)^T$.
\end{itemize}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.2\textwidth]{./figures/figure1a.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1b.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1c.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1d.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1e.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1f.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1g.png}
	\includegraphics[width=0.2\textwidth]{./figures/figure1h.png}
	\caption{Density plots of the a synthetic ($n=1,000$) for each simulation distributions.}
	\label{fig:heat}
\end{figure}

The results in Table \ref{tab:simu} and Figure \ref{fig:simu} show that the large sample asymptotic normal approximation has
 inflated type I error rates for all distributions when $n \le 50$. The $t$ test, Fisher's $Z$ test, Fisher-Yates, and naive permutation tests tend to be over-conservative for the exponential and circular distributions. While for $t_{4.1}$, the type I error is consistently inflated. Note that, for these tests, such deviation cannot be corrected as sample size increases. Instead, they may converge to an arbitrary level, either lower or higher than $\alpha$.

For MVN 1-9, we simulated a range of dependency among uncorrelated $X$ and $Y$, where MVN 1 has the weakest and MVN 9 has the strongest dependency (Figure \ref{fig:heat}). The above four tests showed the type I error rate inflation becomes increasingly severe as the dependency increases. This demonstrates the failure in controlling type I error results from the data being dependent, which can occur when the underlying distribution is non-normal.

The MVN 45 is a case where the dependency of original data is remedied by using the ranks. In this case, the ranks will distribute as if it comes from a bivariate normal distribution, regardless of the distance between the centers of individual Gaussian sub-populations. Therefore, all four tests show well control of the type I error rate.

On the other hand, the studentized permutation test robustly control type I error for all distributions examined, even when the $n$ is as small as 10. This demonstrates a clear advantage of the proposed test over all other commonly used tests for Spearman's correlation coefficient.

\begin{table}[!htbp] \centering
	\caption{Type I error rate of testing $H_0: \rho_s = 0$ versus $H_0: \rho_s > 0$.}
	\label{tab:simu}
	\begin{tabular}{@{\extracolsep{5pt}} rrrrrrrr}
		\\[-1.8ex]\hline
		\hline \\[-1.8ex]
		Distribution & N & $t$ test & Fisher's $Z$ & Fisher-Yates & Asymp Norm & Permute & Stu Permute \\
		\hline \\[-1.8ex]
		MVN & $10$ & $0.0498$ & $0.0420$ & $0.0479$ & $0.1223$ & $0.0495$ & $0.0457$ \\
		 & $20$ & $0.0501$ & $0.0444$ & $0.0499$ & $0.0811$ & $0.0494$ & $0.0487$ \\
		 & $50$ & $0.0512$ & $0.0473$ & $0.0509$ & $0.0639$ & $0.0508$ & $0.0525$ \\
		 & $100$ & $0.0529$ & $0.0484$ & $0.0527$ & $0.0596$ & $0.0530$ & $0.0518$ \\
		 & $200$ & $0.0479$ & $0.0438$ & $0.0493$ & $0.0525$ & $0.0487$ & $0.0496$ \\
		 \hline \\[-1.8ex]
		Exponential & $10$ & $0.0375$ & $0.0295$ & $0.0303$ & $0.1205$ & $0.0368$ & $0.0498$ \\
		 & $20$ & $0.0363$ & $0.0323$ & $0.0270$ & $0.0776$ & $0.0375$ & $0.0497$ \\
		 & $50$ & $0.0362$ & $0.0321$ & $0.0272$ & $0.0629$ & $0.0371$ & $0.0493$ \\
		 & $100$ & $0.0374$ & $0.0329$ & $0.0255$ & $0.0547$ & $0.0374$ & $0.0495$ \\
		 & $200$ & $0.0362$ & $0.0325$ & $0.0278$ & $0.0507$ & $0.0364$ & $0.0476$ \\
		 \hline \\[-1.8ex]
		$t_{4.1}$ & $10$ & $0.0625$ & $0.0531$ & $0.0693$ & $0.1275$ & $0.0606$ & $0.0489$ \\
		 & $20$ & $0.0631$ & $0.0574$ & $0.0846$ & $0.0768$ & $0.0629$ & $0.0480$ \\
		 & $50$ & $0.0631$ & $0.0579$ & $0.0910$ & $0.0614$ & $0.0628$ & $0.0457$ \\
		 & $100$ & $0.0677$ & $0.0611$ & $0.1082$ & $0.0567$ & $0.0672$ & $0.0491$ \\
		 & $200$ & $0.0728$ & $0.0675$ & $0.1149$ & $0.0572$ & $0.0731$ & $0.0521$ \\
		 \hline \\[-1.8ex]
		Circular & $10$ & $0.0124$ & $0.0095$ & $0.0041$ & $0.0907$ & $0.0124$ & $0.0514$ \\
		 & $20$ & $0.0052$ & $0.0040$ & $0.0004$ & $0.0620$ & $0.0048$ & $0.0474$ \\
		 & $50$ & $0.0031$ & $0.0024$ & $0.0001$ & $0.0575$ & $0.0030$ & $0.0487$ \\
		 & $100$ & $0.0021$ & $0.0015$ & $<0.0001$ & $0.0537$ & $0.0024$ & $0.0497$ \\
		 & $200$ & $0.0016$ & $0.0009$ & $<0.0001$ & $0.0489$ & $0.0020$ & $0.0482$ \\
		 \hline \\[-1.8ex]
		MVT & $10$ & $0.0582$ & $0.0502$ & $0.0626$ & $0.1305$ & $0.0584$ & $0.0497$ \\
		 & $20$ & $0.0596$ & $0.0543$ & $0.0726$ & $0.0825$ & $0.0602$ & $0.0496$ \\
		 & $50$ & $0.0542$ & $0.0489$ & $0.0700$ & $0.0608$ & $0.0553$ & $0.0456$ \\
		 & $100$ & $0.0584$ & $0.0535$ & $0.0740$ & $0.0566$ & $0.0574$ & $0.0496$ \\
		 & $200$ & $0.0590$ & $0.0539$ & $0.0852$ & $0.0527$ & $0.0590$ & $0.0496$ \\
		 \hline \\[-1.8ex]
		MVN 1 & $10$ & $0.0541$ & $0.0458$ & $0.0499$ & $0.1271$ & $0.0531$ & $0.0534$ \\
		 & $20$ & $0.0484$ & $0.0451$ & $0.0490$ & $0.0777$ & $0.0482$ & $0.0498$ \\
		 & $50$ & $0.0481$ & $0.0428$ & $0.0478$ & $0.0631$ & $0.0484$ & $0.0492$ \\
		 & $100$ & $0.0497$ & $0.0452$ & $0.0520$ & $0.0564$ & $0.0498$ & $0.0497$ \\
		 & $200$ & $0.0495$ & $0.0462$ & $0.0523$ & $0.0529$ & $0.0493$ & $0.0488$ \\
		 \hline \\[-1.8ex]
		MVN 3 & $10$ & $0.0547$ & $0.0467$ & $0.0539$ & $0.1341$ & $0.0544$ & $0.0497$ \\
		 & $20$ & $0.0557$ & $0.0514$ & $0.0607$ & $0.0814$ & $0.0551$ & $0.0516$ \\
		 & $50$ & $0.0571$ & $0.0498$ & $0.0607$ & $0.0686$ & $0.0562$ & $0.0508$ \\
		 & $100$ & $0.0530$ & $0.0473$ & $0.0651$ & $0.0561$ & $0.0544$ & $0.0497$ \\
		 & $200$ & $0.0529$ & $0.0468$ & $0.0603$ & $0.0517$ & $0.0508$ & $0.0473$ \\
		 \hline \\[-1.8ex]
		MVN 6 & $10$ & $0.0648$ & $0.0555$ & $0.0740$ & $0.1328$ & $0.0643$ & $0.0461$ \\
		 & $20$ & $0.0694$ & $0.0624$ & $0.0880$ & $0.0832$ & $0.0686$ & $0.0526$ \\
		 & $50$ & $0.0658$ & $0.0606$ & $0.0917$ & $0.0638$ & $0.0650$ & $0.0494$ \\
		 & $100$ & $0.0706$ & $0.0644$ & $0.0988$ & $0.0590$ & $0.0694$ & $0.0498$ \\
		 & $200$ & $0.0678$ & $0.0615$ & $0.1034$ & $0.0527$ & $0.0667$ & $0.0496$ \\
		 \hline \\[-1.8ex]
		MVN 9 & $10$ & $0.0944$ & $0.0843$ & $0.1173$ & $0.1405$ & $0.0924$ & $0.0573$ \\
		 & $20$ & $0.0932$ & $0.0865$ & $0.1276$ & $0.0829$ & $0.0913$ & $0.0502$ \\
		 & $50$ & $0.0965$ & $0.0891$ & $0.1414$ & $0.0662$ & $0.0948$ & $0.0519$ \\
		 & $100$ & $0.0997$ & $0.0928$ & $0.1425$ & $0.0593$ & $0.0994$ & $0.0498$ \\
		 & $200$ & $0.0950$ & $0.0882$ & $0.1531$ & $0.0536$ & $0.0953$ & $0.0494$ \\
		 \hline \\[-1.8ex]
		MVN 45 & $10$ & $0.0564$ & $0.0467$ & $0.0487$ & $0.1301$ & $0.0547$ & $0.0515$ \\
		 & $20$ & $0.0492$ & $0.0444$ & $0.0480$ & $0.0791$ & $0.0491$ & $0.0482$ \\
		 & $50$ & $0.0481$ & $0.0445$ & $0.0507$ & $0.0634$ & $0.0482$ & $0.0512$ \\
		 & $100$ & $0.0492$ & $0.0447$ & $0.0498$ & $0.0557$ & $0.0497$ & $0.0488$ \\
		 & $200$ & $0.0521$ & $0.0468$ & $0.0492$ & $0.0552$ & $0.0509$ & $0.0519$ \\
		\hline \\[-1.8ex]
	\end{tabular}
\end{table}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.95\textwidth]{./figures/figure2.png}
	\caption{Type I error rate of testing $H_0: \rho_s = 0$ versus $H_0: \rho_s > 0$.}
	\label{fig:simu}
\end{figure}

\section{Application}
\label{sec:application}

\subsection{TCGA breast cancer data}
As an illustration of our approach, we tested $H_0: \rho_s = 0$ versus $H_1: \rho_s > 0 $ using The Cancer Genome Atlas (TCGA) breast cancer RNA sequencing (RNA-seq) data. The gene abundance was RSEM normalized \cite{li2011rsem}. Fibroblast growth factor (FGF)2, FGF4, FGF7 and FGF20 are representative paracrine FGFs binding to heparan-sulfate proteoglycan and fibroblast growth factor receptors (FGFRs), whereas FGF19, FGF21 and FGF23 are endocrine FGFs binding to Klotho and FGFRs. FGFR1 is relatively frequently amplified and overexpressed in breast and lung cancer, and FGFR2 in gastric cancer. Moreover, FGF2 activates human dermal fibroblasts through transcriptional downregulation of the TP53 gene \cite{katoh2016fgfr}. In this application, we examine whether the transcriptomic abundance of FGFR1 is correlated with that of TP53. To investigate the performance in small sample settings, we selected 18 samples from 17 mucinous carcinoma patients. The scatter plot of log-transformed TP53 and FGFR1 abundances is shown in Figure \ref{fig:tcga_1} (\textit{left}). The marginal normality of data was examined by Shapiro-Wilk test and the bivariate normality was examined by Henze-Zikler test. The $p$ values of Shapiro-Wilk tests for log-transformed TP53 and FGFR1 abundances are 0.6077 and 0.0644, respectively. The $p$ value of Henze-Zikler test is 0.3478. Although there is no statistical significant, the marginal of FGFR1 abundance likely deviates from normal distribution.

The estimated Spearman's correlation is $\hat{\rho}_s=0.4056$. Table \ref{tab:tcga_1} shows the results of hypothesis testing. Only the result of studentized permutation test is non-significant at $\alpha=0.05$ and suggests there is no evidence of positive correlation between TP53 and FGFR1. In fact the biology does not support a positive correlation either, since FGFR1 mediates negative regulation of TP53 by FGF2 at transcriptional level \cite{katoh2016fgfr}. Indeed, if we include all samples ($n\approx 1200$) from TCGA breast cancer cohort, then all tests will fail to reject $H_1$ with $p$-values over 0.5 except for Fisher-Yates test. Together with the results from the simulations, the result by studentized permutation test is clearly more reliable.

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=0.3\textwidth]{./figures/figure3a.png}
	\includegraphics[width=0.3\textwidth]{./figures/figure3b.png}
	\caption{Scatter plot of log-transformed TP53 versus FGFR1 abundance  for TCGA data (\textit{left}), and age versus $-\log(\text{PSA})$ for the PSA data (\textit{right}).}
	\label{fig:tcga_1}
\end{figure}

\begin{table}[!htbp] \centering
	\caption{Results of testing $H_0: \rho_s = 0$ versus $H_1: \rho_s > 0$ for TCGA breast cancer data and PSA data.}
	\label{tab:tcga_1}
	\begin{tabular}{@{\extracolsep{5pt}} lrr}

		\hline \\[-1.8ex]
		Tests & $p$ value (TCGA) & $p$ value (PSA) \\
		\hline \\[-1.8ex]
		$t$ test & $0.047$ & $<0.001$\\
		Fisher's $Z$ & $<0.001$ & $<0.001$\\
		Fisher-Yates & $<0.001$ & $<0.001$\\
		Asymp Norm & $0.039$ & $<0.001$ \\
		Permute & $0.033$ & $<0.001$ \\
		Stu Permute & $0.081$ & $<0.001$ \\
		\hline \\[-1.8ex]
	\end{tabular}
\end{table}

\subsection{PSA data}
The testing methods were also applied to a data set of age and baseline prostate-specific antigen (PSA) levels \cite{sweeney2015chemohormonal}. The data consists of age and PSA levels of 480 subjects, of which 473 have complete paired observations. The sample Spearman's correlation coefficient between age and PSA is $\hat{\rho}_s=-0.1622$. Since the alternative hypothesis of proposed test is $H_1: \rho_s>0 $, we applied a negative log transformation on PSA levels. Similar as the TCGA example, the marginal normality of data was examined by Shapiro-Wilk test, and the bivariate normality was examined by Henze-Zikler test. The $p$ values of Shapiro-Wilk tests for log-transformed age and PSA levels are 0.0208 and 0.0301, respectively. The $p$ value of Henze-Zikler test is $<0.0001$. The results indicates the distribution is not bivariate normal. Figure \ref{fig:tcga_1} (\textit{right}) shows the scatter plot of age versus $-\log(\text{PSA})$. Table \ref{tab:tcga_1} shows that all tests rejects the $H_1$ and conclude there is a non-zero correlation between age and PSA. The is an example where all tests have consistent results. Although the normality tests are significant, such deviation may have been remedied by using the ranks in this specific example.

\section{Discussion}
\label{sec:discussion}

Conventional tests of the Spearman's correlation rely on normality assumption, including $t$-test, Fisher's $Z$ transformation, and naive permutation test, which fails to control Type I error rates when the assumption is violated. This was illustrated in our simulations studies (Section \ref{sec:simulations}). Such defect cannot be remedied by transforming the marginal distributions such as by Fisher-Yates coefficient. Notably, the deviation from bivariate normality can result in a convergence of type I error rate to an arbitrary level when $n\to \infty$. This indicates that, under scenarios when two random variables are uncorrelated but dependent, the type I error will not be controlled at desired level no matter how large the sample size is.  On the other hand, the Serfling's test based on delta method guarantees that the type I error rate converges to $\alpha$ as long as the fourth order moment is finite. However, it typically suffers an inflated type I error when sample size is under 50.

In this work, we present a robust Spearman's correlation permutation test based on studentized statistic for testing $H_0: \rho_s = 0$ versus $H_1: \rho_s > 0$. The proposed approach is inspired by the work by DiCiccio and Romano \cite{diciccio2017robust}, which was developed for Pearson's correlation. Through extensive simulation studies and real world application, we show the proposed test controls type I error even when sample size is as small as 10 and normality assumption is violated. Therefore, the test is valid in general cases. In addition, the studentized statistic can also be used for bootstrapping tests, so as to test for more general point null hypotheses \cite{hutson2019robust}. In conclusion, the proposed studentized permutation test should be used as a routine for testing non-zero Spearman's correlation coefficient.

\section*{Acknowledgments}
This work was supported by Roswell Park Cancer Institute and National Cancer Institute (NCI) grant P30CA016056,  NRG Oncology Statistical and Data Management Center grant U10CA180822 and IOTN Moonshot grant U24CA232979-01. The results shown here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. The PSA data example is based on research using information obtained from www.projectdatasphere.org, which is maintained
by {\em Project Data Sphere, LLC}. Neither {\em Project Data Sphere, LLC}
nor the owner(s) of any information from the website
have contributed to, approved or are in any way responsible for the contents of this publication.

\bibliographystyle{unsrt}