# Exact first moments of the RV coefficient by invariant orthogonal
integration

François Bavaud
University of Lausanne, Switzerland
<EMAIL_ADDRESS>

###### Abstract

The RV coefficient measures the similarity between two multivariate
configurations, and its significance testing has attracted various proposals
in the last decades. We present a new approach, the invariant orthogonal
integration, permitting to obtain the exact first four moments of the RV
coefficient under the null hypothesis. It consists in averaging along the Haar
measure the respective orientations of the two configurations, and can be
applied to any multivariate setting endowed with Euclidean distances between
the observations. Our proposal also covers the weighted setting of
observations of unequal importance, where the exchangeability assumption,
justifying the usual permutation tests, breaks down.

The proposed RV moments express as simple functions of the kernel eigenvalues
occurring in the weighted multidimensional scaling of the two configurations.
The expressions for the third and fourth moments seem original. The first
three moments can be obtained by elementary means, but computing the fourth
moment requires a more sophisticated apparatus, the Weingarten calculus for
orthogonal groups. The central role of standard kernels and their spectral
moments is emphasized.

_Keywords_ RV coefficient $\cdot$ weighted multidimensional scaling $\cdot$
spectral moments $\cdot$ invariant orthogonal integration $\cdot$ Weingarten
calculus

##  1 Introduction

The RV coefficient is a well-known measure of similarity between two datasets,
each consisting of multivariate profiles measured on the same $n$ observations
or objects. This contribution proposes a new approach, the invariant
orthogonal integration, permitting to obtain the exact first four moments of
the RV coefficient under the null hypothesis of absence of relation between
the two datasets. The main results, theorem 1 and corollary 1, are exposed in
section 3.1. The approach is fully nonparametric, and allows the handling of
weighted objets, typically made of aggregates such as regions, documents or
species, which abound in multivariate analysis.

In the present distance-based data-analytic approach, data sets are
constituted by weighted configurations specified by the object weights
together with their pair dissimilarities, assumed to be squared Euclidean.
Factorial coordinates, reproducing the dissimilarities, and permitting a
maximum compression of the configuration inertia, obtain by weighted
multidimensional scaling. The latter, seldom exposed in the literature and
hence briefly recalled in section 2.1, is a direct generalization of classical
scaling. The central step is provided by the spectral decomposition of the
matrix of weighted centered scalar products or kernel. It permits to decompose
the spectral eigenspace into a trivial one-dimensional part, determined by the
object weights, common to both configurations, and a non-trivial part of
dimension $n-1$, orthogonal to the square root of the weights. The weighted RV
coefficient obtains as the normalized scalar product between the kernels of
the two configurations (section 2.2), and turns out to be equivalent to its
original definition expressed by cross-covariances (Escoufier,, 1973; Robert
and Escoufier,, 1976).

After recalling the above preliminaries, somewhat lengthy but necessary, the
heart of this contribution can be uncovered: invariant orthogonal integration
consists in computing the expected null moments of the RV coefficient by
averaging, along the invariant Haar orthogonal measure in the non-trivial
eigenspace, the orientations of one configuration with respect to the other,
by orthogonal transformation of, say, the first eigenspace (section 3.2). It
constitutes a distinct alternative, with different outcomes, to the
traditional permutation approach, whose exchangeability assumption breaks down
for weighted objects: typically, the profile dispersion is expected to be
larger for lighter objects (Bavaud,, 2013) and the $n$ object scores cannot
follow the same distribution. The present approach also yields a novel
significance test for the RV coefficient (equation 16), taking into account
skewness and kurtosis corrections to the usual normal approximation.

Computing the moments of the RV coefficient requires to evaluate the
orthogonal coefficients (23) constituted by Haar expectations of orthogonal
monomials. Low-order moments can be computed, with increasing difficulty, by
elementary means (section 3.3), but the fourth-order moment requires a more
systematic approach (section 3.6), provided by the Weingarten calculus
developed by workers in random matrix theory and free probability. Both
procedures yield the same results for low-order moments (section 3.7), which
is both expected and reassuring.

The first RV moment (11) coincides with all known proposals. The second
centered RV moment (12) is simpler than its permutation analog, and underlines
the effective dimensionality of a configuration. The third centered RV moment
(13) is particularly enlightening: the RV skewness is simply proportional to
the product of the spectral skewness of both configurations, thus elucidating
the often noticed positive skewness of the RV coefficient. The expression for
the fourth centered RV moment (9), (14) is also simple to express and to
compute, yet more difficult to interpret.

##  2 Euclidean configurations in a weighted setting: a concise remainder

###  2.1 Weighted multidimensional scaling and standard kernels

Consider $n$ objects endowed with positive weights $f_{i}>0$ with
$\sum_{i=1}^{n}f_{i}=1$, as well with pairwise dissimilarities
$\mathbf{D}=(D_{ij})$ between pairs of objects. The $n\times n$ matrix
$\mathbf{D}$ is assumed to be squared Euclidean, that is of the form
$D_{ij}=\|\mathbf{x}_{i}-\mathbf{x}_{j}\|^{2}$ for
$\mathbf{x}_{i},\mathbf{x}_{j}\in\mathbb{R}^{r}$, with $r\leq n-1$. The pair
$(\mathbf{f},\mathbf{D})$ constitutes a weighted configuration, with
$f_{i}=1/n$ for unweighted configurations.

Weighted multidimensional scaling aims at determining object coordinates
$\mathbf{X}=(x_{i\alpha})\in\mathbb{R}^{n\times r}$ reproducing the
dissimilarities $\mathbf{D}$ while expressing a maximum amount of dispersion
or inertia $\Delta$ (3) in low dimensions. It is performed by the following
weighted generalization of the well-known Torgerson–Gower scaling procedure
(see e.g. Borg and Groenen,, 2005): first, define
$\bm{\Pi}=\mbox{diag}(\mathbf{f})$, as well as the weighted centering matrix
$\mathbf{H}=\mathbf{I}_{n}-\bm{1}_{n}\mathbf{f}^{\top}$, which obeys
$\mathbf{H}^{2}=\mathbf{H}$. However, $\mathbf{H}^{\top}\neq\mathbf{H}$,
unless $\mathbf{f}$ is uniform.

Second, compute the matrix $\mathbf{B}$ of scalar products by double
centering: $\mathbf{B}=-\frac{1}{2}\mathbf{H}\,\mathbf{D}\,\mathbf{H}^{\top}$.
Third, define the $n\times n$ kernel $\mathbf{K}$ as the matrix of weighted
scalar products :

$\mathbf{K}=\sqrt{\bm{\Pi}}\,\mathbf{B}\sqrt{\bm{\Pi}}\qquad,\qquad\mbox{that
is}\qquad K_{ij}=\sqrt{f_{i}f_{j}}B_{ij}\enspace.$

Fourth, perform the spectral decomposition with $\hat{\mathbf{U}}$ orthogonal
and $\hat{\bm{\Lambda}}$ diagonal

$\mathbf{K}=\hat{\mathbf{U}}\hat{\bm{\Lambda}}\hat{\mathbf{U}}^{\top}\qquad\qquad\hat{\mathbf{U}}\hat{\mathbf{U}}^{\top}=\hat{\mathbf{U}}^{\top}\hat{\mathbf{U}}=\mathbf{I}_{n}\qquad\qquad\hat{\bm{\Lambda}}=\mbox{diag}(\bm{\lambda})\enspace.$
(1)

By construction, $\mathbf{K}$ possesses one trivial eigenvalue $\lambda_{0}=0$
associated to the eigenvector $\sqrt{\mathbf{f}}$ and $n-1$ non-negative
eigenvalues decreasingly ordered as
$\lambda_{1}\,\geq\,\lambda_{2}\,\geq\,\ldots\,\geq\,\lambda_{n-1}\geq 0$,
among which $r=\mbox{rg}(\mathbf{K})$ are strictly positive.

From now on the trivial eigenspace will be discarded: set
$\hat{\mathbf{U}}=(\sqrt{\mathbf{f}}|\mathbf{U})$, where
$\mathbf{U}\in\mathbb{R}^{n\times(n-1)}$ and
$\bm{\Lambda}=\mbox{diag}(\lambda_{1},\ldots,\lambda_{n-1})$. Direct
substitution from (1) yields

$\mathbf{K}=\mathbf{U}\bm{\Lambda}\mathbf{U}^{\top}\quad\qquad\mathbf{U}\mathbf{U}^{\top}=\mathbf{I}_{n}-\sqrt{\mathbf{f}}\sqrt{\mathbf{f}}^{\top}\qquad\quad\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}_{n-1}\qquad\quad\mathbf{U}^{\top}\sqrt{\mathbf{f}}=\mathbf{0}_{n}\enspace.$
(2)

Finally, the searched for coordinates obtain as
$\mathbf{X}=\bm{\Pi}^{-\frac{1}{2}}\mathbf{U}\bm{\Lambda}^{\frac{1}{2}}$, that
is $x_{i\alpha}=u_{i\alpha}\sqrt{\lambda_{\alpha}}/\sqrt{f_{i}}$. One verifies
easily that

$D_{ij}=\sum_{\alpha=1}^{n-1}(x_{i\alpha}-x_{j\alpha})^{2}\qquad\qquad\qquad\Delta=\frac{1}{2}\sum_{i,j=1}^{n}f_{i}f_{j}D_{ij}=\mbox{Tr}(\mathbf{K})=\sum_{\alpha=1}^{n-1}~\lambda_{\alpha}\enspace.$
(3)

The kernels considered here are positive semi-definite and obey in addition
$\mathbf{K}\sqrt{\mathbf{f}}=\mathbf{0}_{n}$. We call them standard kernels.
They can be related to the weighted version of centered kernels of Machine
Learning (see e.g. Cortes et al.,, 2012). To each weighted configuration
$(\mathbf{f},\mathbf{D})$ corresponds a unique standard kernel $\mathbf{K}$,
and conversely.

The matrix
$\mathbf{K}_{0}=\mathbf{I}_{n}-\sqrt{\mathbf{f}}\sqrt{\mathbf{f}}^{\top}$
appearing in (2) constitutes a standard kernel, referred to as the neutral
kernel in view of property
$\mathbf{K}_{0}\mathbf{K}=\mathbf{K}_{0}\mathbf{K}=\mathbf{K}$ for any
standard kernel $\mathbf{K}$. The corresponding dissimilarities are the
weighted discrete distances

$D_{ij}^{0}=\begin{cases}\frac{1}{f_{i}}+\frac{1}{f_{j}}&\text{ for $i\neq
j$}\\\ 0&\text{ otherwise.}\end{cases}$

###  2.2 The RV coefficient

Consider two weighted configurations $(\mathbf{f},\mathbf{D}_{X})$ and
$(\mathbf{f},\mathbf{D}_{Y})$ endowed with the same weights $\mathbf{f}$, or
equivalently two standard kernels $\mathbf{K}_{X}$ and $\mathbf{K}_{Y}$
(Figure 1). Their similarity can be measured by the weighted RV coefficient
defined as

${\tt RV}={\tt
RV}_{XY}=\frac{\mbox{Tr}(\mathbf{K}_{X}\;\mathbf{K}_{Y})}{\sqrt{\mbox{Tr}(\mathbf{K}^{2}_{X})\mbox{Tr}(\mathbf{K}^{2}_{Y})}}$
(4)

which constitutes the cosine similarity between the vectorized matrices
$\mathbf{K}_{X}$ and $\mathbf{K}_{Y}$. As a consequence, ${\tt RV}_{XY}\geq 0$
(since $\mathbf{K}_{X}$ and $\mathbf{K}_{Y}$ are positive semi-definite),
${\tt RV}_{XY}\leq 1$ (by the Cauchy-Schwarz inequality) and ${\tt
RV}_{XX}=1$.

Figure 1: Two weighted configurations $(\mathbf{f},\mathbf{D}_{X})$ (left) and
$(\mathbf{f},\mathbf{D}_{Y})$ (right) embedded in $\mathbb{R}^{n-1}$

Quantity (4) is a straightforward weighted generalization of the RV
coefficient (Escoufier,, 1973; Robert and Escoufier,, 1976): consider
multivariate features $\mathbf{X}\in\mathbb{R}^{n\times p}$ and
$\mathbf{Y}\in\mathbb{R}^{n\times q}$, directly entering into the definition
of $\mathbf{D}_{X}$ and $\mathbf{D}_{Y}$ as coordinates, or equivalently as
$\mathbf{K}_{X}=\sqrt{\bm{\Pi}}\mathbf{X}_{c}\mathbf{X}_{c}^{\top}\sqrt{\bm{\Pi}}$
and
$\mathbf{K}_{Y}=\sqrt{\bm{\Pi}}\mathbf{Y}_{c}\mathbf{Y}_{c}^{\top}\sqrt{\bm{\Pi}}$,
where $\mathbf{X}_{c}=\mathbf{H}\mathbf{X}$ and
$\mathbf{Y}_{c}=\mathbf{H}\mathbf{Y}$ are the centered scores.

The weighted covariances are
$\bm{\Sigma}_{XX}=\mathbf{X}_{c}^{\top}\bm{\Pi}\mathbf{X}_{c}$ and
$\bm{\Sigma}_{YY}=\mathbf{Y}_{c}^{\top}\bm{\Pi}\mathbf{Y}_{c}$. The cross-
covariances are $\bm{\Sigma}_{XY}=\mathbf{X}_{c}^{\top}\bm{\Pi}\mathbf{Y}_{c}$
and
$\bm{\Sigma}_{YX}=\mathbf{Y}_{c}^{\top}\bm{\Pi}\mathbf{X}_{c}=\bm{\Sigma}_{XY}^{\top}$.
The original RV coefficient is defined in the feature space as

${\tt
RV}_{XY}=\frac{\mbox{Tr}(\bm{\Sigma}_{XY}\;\bm{\Sigma}_{YX})}{\sqrt{\mbox{Tr}(\bm{\Sigma}^{2}_{XX})\mbox{Tr}(\bm{\Sigma}^{2}_{YY})}}\enspace.$
(5)

Proving the identity of (4) and (5) is easy.

##  3 Computing the moments of the RV coefficient by invariant orthogonal
integration

###  3.1 Main result and significance testing

Define the CV coefficient by the quantity ${\tt
CV}=\mbox{Tr}(\mathbf{K}_{X}\mathbf{K}_{Y})$.

######  Theorem 1 (Main result).

Under invariant orthogonal integration (section 3.2), the expectation of the
CV coefficient and its centered moments of order 2, 3 and 4 are

$\displaystyle\mathbb{E}({\tt CV})=$
$\displaystyle(n-1)\,\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\,\overline{\bm{\mu}\vphantom{\bm{\lambda}}}$
(6) $\displaystyle\mathbb{E}({\tt CV}_{c}^{2})=$
$\displaystyle\frac{2(n-1)^{2}}{(n-2)(n+1)}\,\overline{\bm{\lambda}_{c}^{2}}\;\overline{\bm{\mu}_{c}^{2}}$
(7) $\displaystyle\mathbb{E}({\tt CV}_{c}^{3})=$
$\displaystyle\frac{8(n-1)^{3}}{(n-3)(n-2)(n+1)(n+3)}\,\overline{\bm{\lambda}_{c}^{3}}\;\overline{\bm{\mu}_{c}^{3}}$
(8) $\displaystyle\mathbb{E}({\tt CV}_{c}^{4})=$
$\displaystyle{\frac{12\,(n-1)^{3}}{(n-4)(n-3)(n-2)n(n+1)(n+3)(n+5)}}\,\Bigl{\\{}4\,(n^{2}-n+2)\,\overline{\bm{\lambda}_{c}^{4}}\>\overline{\bm{\mu}_{c}^{4}}+(n^{4}$
(9) $\displaystyle+n^{3}-15n^{2}$
$\displaystyle-13n+98)\,\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{2}}-4\,(2n^{2}-n-7)\,(\overline{\bm{\lambda}_{c}^{4}}\>\overline{\bm{\mu}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{2}}+\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{4}})\Bigr{\\}}\enspace.$

where ${\tt CV}_{c}={\tt CV}-\mathbb{E}({\tt CV})$. Spectral moments and
centered spectral moments read

$\overline{\bm{\lambda}^{q}}=\frac{1}{n-1}\sum_{\alpha=1}^{n}\lambda_{\alpha}^{q}=\frac{1}{n-1}\sum_{\alpha=1}^{n}\mbox{Tr}(\mathbf{K}_{X}^{q})=\mbox{tr}(\mathbf{K}_{X}^{q})\qquad\qquad\overline{\bm{\lambda}^{q}_{c}}=\frac{1}{n-1}\sum_{\alpha=1}^{n}(\lambda^{c}_{\alpha})^{q}$
(10)

where $\lambda^{c}_{\alpha}=\lambda_{\alpha}-\overline{\bm{\lambda}}$ and
$\mbox{tr}(\mathbf{A})=\mbox{Tr}(\mathbf{A})/(n-1)$ denotes the normalized
trace. Centered spectral moments can be transformed into normalized traces,
and conversely. For instance,
$\overline{\bm{\lambda}_{c}^{3}}=\mbox{tr}(\mathbf{K}_{X}^{3})-3\,\mbox{tr}(\mathbf{K}_{X}^{2})\,\mbox{tr}(\mathbf{K}_{X})+2\,\mbox{tr}^{3}(\mathbf{K}_{X})$.

Identity ${\tt RV}={\tt
CV}/\sqrt{\mbox{Tr}(\mathbf{K}^{2}_{X})\mbox{Tr}(\mathbf{K}^{2}_{Y})}={\tt
CV}/[(n-1)\sqrt{\overline{\bm{\lambda}^{2}}\;\overline{\bm{\mu}^{2}}}\,]$
directly yields:

######  Corollary 1 (First cumulants of the RV coefficient).

Under invariant orthogonal integration, the first cumulants of the RV
coefficient, that is its expectation, variance, skewness and excess kurtosis
are, in order,

$\displaystyle\mathbb{E}({\tt RV})=$
$\displaystyle\frac{1}{n-1}\frac{\mbox{\rm Tr}(\mathbf{K}_{X})\,\mbox{\rm
Tr}(\mathbf{K}_{Y})}{\sqrt{\mbox{\rm Tr}(\mathbf{K}_{X}^{2})\,\mbox{\rm
Tr}(\mathbf{K}_{Y}^{2})}}=\frac{\overline{\bm{\lambda}}\overline{\bm{\mu}\>\vphantom{\lambda}}}{\sqrt{\overline{\bm{\lambda}^{2}}\>\overline{\bm{\mu}^{2}}}}=\frac{\sqrt{\nu({\bm{\lambda}})\nu({\bm{\mu}})}}{n-1}$
(11) $\displaystyle\mbox{$\mathbb{V}$\rm ar}({\tt RV})=$
$\displaystyle\mathbb{E}({\tt
RV}_{c}^{2})=\frac{2(n-1-\nu({\bm{\lambda}}))(n-1-\nu({\bm{\mu}}))}{(n-2)(n-1)^{2}(n+1)}$
(12) $\displaystyle\mbox{$\mathbb{A}$}({\tt RV})=$
$\displaystyle\frac{\mathbb{E}({\tt
RV}_{c}^{3})}{\mathbb{E}^{\frac{3}{2}}({\tt
RV}_{c}^{2})}=\frac{\sqrt{8(n-2)(n+1)}}{(n-3)(n+3)}\>a({\bm{\lambda}})\>a({\bm{\mu}})$
(13)
$\displaystyle{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}({\tt
RV})=$ $\displaystyle\frac{\mathbb{E}({\tt RV}_{c}^{4})}{\mathbb{E}^{2}({\tt
RV}^{2}_{c})}-3=\frac{\mathbb{E}({\tt CV}_{c}^{4})}{\mathbb{E}^{2}({\tt
CV}_{c}^{2})}-3\qquad\quad\mbox{\rm(see (\ref{4MainTheo4}) and
(\ref{4MainTheo2}))}$ (14)

where ${\tt RV}_{c}={\tt RV}-\mathbb{E}({\tt RV})$, and
$a({\bm{\lambda}})=\overline{\bm{\lambda}_{c}^{3}}/(\overline{\bm{\lambda}_{c}^{2}})^{\frac{3}{2}}$
is the spectral skewness. The quantity

$\nu({\bm{\lambda}})=\frac{\mbox{Tr}^{2}(\mathbf{K}_{X})}{\mbox{Tr}(\mathbf{K}_{X}^{2})}=\frac{(\sum_{\alpha\geq
1}\lambda_{\alpha})^{2}}{\sum_{\alpha\geq
1}\lambda^{2}_{\alpha}}=(n-1)\frac{\overline{\bm{\lambda}}^{2}}{\overline{\bm{\lambda}^{2}}}$
(15)

has appeared at times as an adjusted degrees of freedom in multivariate tests
of the general linear model (see e.g. Geisser and Greenhouse,, 1958; Worsley
and Friston,, 1995; Schlich,, 1996; Abdi,, 2010). It provides a measure of
sphericity or effective dimensionality of configuration
$(\mathbf{f},\mathbf{D}_{X})$. Its minimum $\nu({\bm{\lambda}})=1$ is attained
for univariate configurations. Its maximum $\nu({\bm{\lambda}})=n-1$ is
attained for uniform dilatations of the discrete distances
$\mathbf{D}^{0}_{X}$ (section 2.1), in which case $\mbox{$\mathbb{V}$\rm
ar}({\tt RV})=0$ since ${\tt RV}$ is then concentrated on
$\sqrt{\nu({\bm{\mu}})/(n-1)}$.

The second-order Cornish-Fisher cumulant expansion permits to approximatively
redress the normal quantiles by taking into account the skewness and the
"taildeness" of a non-normal distribution (see e.g. Kendall and Stuart,, 1977;
Amédée-Manesme et al.,, 2019). The observed RV is statistically significant at
level $\alpha$ if (one-tailed test)

$\begin{split}&\qquad\qquad\qquad\underbrace{\frac{{\tt RV}-\mathbb{E}({\tt
RV})}{\sqrt{\mbox{$\mathbb{V}$ar}({\tt
RV})}}}_{\mbox{\em$z$-score}}\qquad\>\>>\>\>\underbrace{u_{1-\alpha}}_{\mbox{\em
standard normal quantile}}\qquad\\\
&+\quad\underbrace{\frac{\mbox{$\mathbb{A}$}({\tt
RV})}{6}(u^{2}_{1-\alpha}-1)+\frac{{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}({\tt
RV})}{24}(u^{3}_{1-\alpha}-3u_{1-\alpha})-\frac{\mbox{$\mathbb{A}^{2}$}({\tt
RV})}{36}(2u^{3}_{1-\alpha}-5u_{1-\alpha})}_{\mbox{\em correction to the
normal distribution}}\enspace.\end{split}$ (16)

###  3.2 Invariant orthogonal integration

The rest of the paper is devoted to presenting invariant orthogonal
integration and proving Theorem 1.

Consider two standard kernels
$\mathbf{K}_{X}=\mathbf{U}\bm{\Lambda}\mathbf{U}^{\top}$ and
$\mathbf{K}_{Y}=\mathbf{V}\mathbf{M}\mathbf{V}^{\top}$ with
$\mathbf{M}=\mbox{diag}(\bm{\mu})$, where
$\mathbf{U},\mathbf{V}\in\mathbb{R}^{n\times(n-1)}$ (2). The numerator ${\tt
CV}$ of ${\tt RV}$ in (4) reads

${\tt
CV}=\mbox{Tr}(\mathbf{K}_{X}\mathbf{K}_{Y})=\mbox{Tr}(\mathbf{U}\bm{\Lambda}\mathbf{U}^{\top}\mathbf{V}\mathbf{M}\mathbf{V}^{\top})=\sum_{\alpha=1}^{n-1}\sum_{\beta=1}^{n-1}\lambda_{\alpha}\mu_{\beta}P_{\alpha\beta}$
(17)

where

$P_{\alpha\beta}=\sum_{i,j=1}^{n}u_{i\alpha}u_{j\alpha}v_{i\beta}v_{j\beta}=(\sum_{i=1}^{n}u_{i\alpha}v_{i\beta})^{2}\enspace.$

Identities
$\mathbf{U}\mathbf{U}^{\top}=\mathbf{I}_{n}-\sqrt{\mathbf{f}}\sqrt{\mathbf{f}}^{\top}$,
$\mathbf{V}^{\top}\sqrt{\mathbf{f}}=\mathbf{0}_{n}$ and
$\mathbf{V}^{\top}\mathbf{V}=\mathbf{I}_{n-1}$ from (2) imply the joint
orthogonality property
$\mathbf{V}^{\top}\mathbf{U}\mathbf{U}^{\top}\mathbf{V}=\mathbf{I}_{n-1}$,
yielding

$P_{\bullet\beta}=\sum_{\alpha=1}^{n-1}P_{\alpha\beta}=\sum_{i,j=1}^{n}(\delta_{ij}-\sqrt{f_{i}f_{j}})\,v_{i\beta}v_{j\beta}=\sum_{i=1}^{n}v_{i\beta}^{2}=1\enspace$
(18)

and, similarly, $P_{\alpha\bullet}=1$. Hence, the matrix
$\mathbf{P}=(P_{\alpha\beta})\in\mathbb{R}^{(n-1)\times(n-1)}$ is non-
negative, and doubly stochastic: it expresses as a mixture of permutations of
$\mathcal{S}_{n-1}$ (Birkhoff–von Neumann theorem). In particular, one gets
the crude estimate

$\sum_{\alpha\geq 1}\lambda_{\alpha}\mu_{n-\alpha}\>\leq\>{\tt
CV}\>\leq\>\sum_{\alpha\geq 1}\lambda_{\alpha}\mu_{\alpha}\enspace.$

The null hypothesis $H_{0}$ states that the two configurations
$(\mathbf{f},\mathbf{D}_{X})$ and $(\mathbf{f},\mathbf{D}_{Y})$ are unrelated.
Under $H_{0}$, any relative orientation of a configuration with respect to the
other is equally likely. Hence, the first configuration will be rotated by
replacing $\mathbf{U}=(u_{i\alpha})$ by $\mathbf{U}\mathbf{T}$, where
$\mathbf{T}=(t_{a\alpha})\in\mathbb{O}_{n-1}$, the orthogonal group of
dimension $n-1$. This rotation acts in the non-trivial eigenspace only,
leaving the weights $\mathbf{f}$ unchanged. The term
$\mbox{Tr}(\mathbf{K}^{2}_{X})$ remains the same, and the ${\tt CV}$
coefficient becomes

${\tt
CV}(\mathbf{T})=\mbox{Tr}(\mathbf{U}\mathbf{T}\bm{\Lambda}\mathbf{T}^{\top}\mathbf{U}^{\top}\mathbf{V}\mathbf{M}\mathbf{V}^{\top})=\sum_{\alpha=1}^{n-1}\sum_{\beta=1}^{n-1}\lambda_{\alpha}\mu_{\beta}P_{\alpha\beta}(\mathbf{T})$
(19)

where

$P_{\alpha\beta}(\mathbf{T})=\sum_{a,b=1}^{n-1}t_{a\alpha}t_{b\alpha}\sum_{i,j=1}^{n}u_{ia}u_{jb}v_{i\beta}v_{j\beta}\enspace.$
(20)

The idea of invariant orthogonal integration is to compute the expectation of
the moments

$\mathbb{E}({\tt CV}^{q}):=\int_{\mathbb{O}_{n-1}}{\tt
CV}^{q}(\mathbf{T})\>d\mu(\mathbf{T})\qquad\qquad q=1,2,\ldots$ (21)

by averaging over all possible rotations $\mathbf{T}\in\mathbb{O}_{n-1}$
distributed by the invariant Haar measure $d\mu(\mathbf{T})$ normalized to
$\int_{\mathbb{O}_{n-1}}d\mu(\mathbf{T})=1$. The moment generating function
reads

$\mathbb{E}(\exp(t\,{\tt
CV}))=\int_{\mathbb{O}_{n-1}}\exp(t\>\mbox{Tr}(\mathbf{T}\bm{\Lambda}{\mathbf{T}^{\top}}\mathbf{A}))\>d\mu(\mathbf{T})\quad\quad\mbox{with}~\quad\mathbf{A}=\mathbf{U}^{\top}\mathbf{V}\mathbf{M}\mathbf{V}^{\top}\mathbf{U}\,.$
(22)

Define $[n]=\\{1,2,\ldots,n\\}$. Computing (21) involves the orthogonal
coefficients, defined in whole generality as

$\mathcal{I}_{\mathbf{a}}^{\bm{\omega}}=\int_{\mathbb{O}_{n-1}}d\mu(\mathbf{T})\>t_{a_{1}\omega_{1}}\,t_{a_{2}\omega_{2}}\ldots
t_{a_{2q}\omega_{2q}}$ (23)

where the multi-indices ${\mathbf{a}}=(a_{1}a_{2}\ldots a_{2q})$ and
${\bm{\omega}}=(\omega_{1}\omega_{2}\ldots\omega_{2q})$ are elements of
$[n-1]^{2q}$ that is, ${\mathbf{a}}$ and ${\bm{\omega}}$ are words of length
$2q$ on the alphabet $[n-1]$.

To ease the notations, define $A_{q}=[n-1]^{q}$ and, for
${\bm{\alpha}}=(\alpha_{1}\ldots\alpha_{q})\in A_{q}$, define
${\bm{\alpha}}{\bm{\alpha}}=(\alpha_{1}\alpha_{1}\alpha_{2}\alpha_{2}\ldots\alpha_{q}\alpha_{q})\in
A_{2q}$. Identities (19), (20), (21) and (23) yield

$\displaystyle\mathbb{E}({\tt CV}^{q})$
$\displaystyle=\sum_{\alpha_{1}\ldots\alpha_{q}=1}^{n-1}\underbrace{\lambda_{\alpha_{1}}\cdots\lambda_{\alpha_{q}}}_{\lambda_{\bm{\alpha}}}\sum_{\beta_{1}\ldots\beta_{q}=1}^{n-1}\underbrace{\mu_{\beta_{1}}\cdots\mu_{\beta_{q}}}_{\mu_{\bm{\beta}}}\sum_{i_{1}\ldots
i_{2q}=1}^{n}\underbrace{v_{i_{1}\beta_{1}}v_{i_{2}\beta_{1}}v_{i_{3}\beta_{2}}v_{i_{4}\beta_{2}}\ldots
v_{i_{2q-1}\beta_{q}}v_{i_{2q}\beta_{q}}}_{v_{{\mathbf{i}}{\bm{\beta}}{\bm{\beta}}}}$
$\displaystyle\quad\times\sum_{a_{1}\ldots a_{2q}=1}^{n-1}$
$\displaystyle\underbrace{\mathcal{I}_{a_{1}a_{2}\ldots
a_{2q-1}a_{2q}}^{\alpha_{1}\alpha_{1}\ldots\alpha_{q}\alpha_{q}}}_{\mathcal{I}_{\mathbf{a}}^{{\bm{\alpha}}{\bm{\alpha}}}}\underbrace{u_{i_{1}a_{1}}u_{i_{2}a_{2}}\ldots
u_{i_{2q}a_{2q}}}_{u_{{\mathbf{i}}{\mathbf{a}}}}=\sum_{{\bm{\alpha}}\in
A_{q}}\lambda_{\bm{\alpha}}\sum_{{\bm{\beta}}\in
A_{q}}\mu_{\bm{\beta}}\underbrace{\sum_{{\mathbf{i}}\in[n]^{2q}}v_{{\mathbf{i}}{\bm{\beta}}{\bm{\beta}}}\sum_{{\mathbf{a}}\in
A_{2q}}\mathcal{I}_{\mathbf{a}}^{{\bm{\alpha}}{\bm{\alpha}}}\,u_{{\mathbf{i}}{\mathbf{a}}}}_{\mathbb{E}(P_{{\bm{\alpha}}{\bm{\beta}}})=\mathbb{E}(P_{\alpha_{1}\beta_{1}}P_{\alpha_{2}\beta_{2}}\cdots
P_{\alpha_{q}\beta_{q}})}\,.$ (24)

The knowledge of $\mathcal{I}_{\mathbf{a}}^{{\bm{\alpha}}{\bm{\alpha}}}$
together with joint orthogonality properties will yield exact expressions for
$\mathbb{E}({\tt CV}^{q})$ in terms of spectral moments of ${\bm{\lambda}}$
and ${\bm{\mu}}$, or equivalently in terms of traces of integer powers of
$\mathbf{K}_{X}$ and $\mathbf{K}_{Y}$, as demonstrated in the next sections
for $q=1,2,3,4$.

###  3.3 Computing low-order orthogonal coefficients

Evaluating the orthogonal coefficients (23) is a major topic in random matrix
theory and free probability, and its systematic handling is presented in
section 3.6. Yet, as observed by some authors (see e.g. Aubert and Lam,, 2003;
Braun,, 2006; Yamamoto and Kudo,, 2017), well-inspired invariance
considerations (Lemmas 1 and 2 below) suffice in determining more directly the
values of the orthogonal coefficients of low order.

Since $d\mu(\mathbf{-T})=d\mu(\mathbf{T})$, coefficients (23) are zero unless
each index in ${\mathbf{a}}$ and in ${\bm{\omega}}$ occurs an even number of
times, with a total of $2q$ occurrences, where $q$ defines the order of the
orthogonal coefficient. Also, applying the same permutation on the two multi-
indexes, or exchanging the multi-indexes leaves the coefficients unchanged.
Furthermore, the particular value taken by an index is irrelevant: only
matters its multiplicity. For instance:

$\mathcal{I}_{abcd}^{\alpha\beta\gamma\delta}=\mathcal{I}_{\alpha\beta\gamma\delta}^{abcd}=\mathcal{I}_{dacb}^{\delta\alpha\gamma\beta}\neq\mathcal{I}_{adcb}^{\delta\alpha\gamma\beta}\quad\mbox{in
general ;
}\qquad\qquad\mathcal{I}_{abcc}^{\alpha\alpha\gamma\gamma}=\mathcal{I}_{accb}^{\alpha\alpha\gamma\gamma}=0\quad\mbox{if
$a\neq b$}\enspace.$

Also, for $\alpha\neq\gamma$ and $a\neq b$,

$\mathcal{I}_{aabb}^{\alpha\alpha\gamma\gamma}=\mathcal{I}_{bbaa}^{\alpha\alpha\gamma\gamma}=\mathcal{I}_{1122}^{2211}=\mathcal{I}_{1122}^{1122}\enspace.$

######  Lemma 1 (proved in the Appendix).

Let $\alpha\neq\gamma$ and let ${\bm{\varepsilon}}$ be a multi-index not
containing $\alpha,\gamma$. For any indices $a,b,c,d$ and multi-index
${\bm{e}}$ of the same size as ${\bm{\varepsilon}}$

$\displaystyle\mathcal{I}_{abcd{{\bm{e}}}}^{\alpha\alpha\alpha\alpha{{\bm{\varepsilon}}}}$
$\displaystyle=\mathcal{I}_{abcd{{\bm{e}}}}^{\alpha\alpha\gamma\gamma{{\bm{\varepsilon}}}}+\mathcal{I}_{abcd{{\bm{e}}}}^{\alpha\gamma\alpha\gamma{{\bm{\varepsilon}}}}+\mathcal{I}_{abcd\mathbf{e}}^{\alpha\gamma\gamma\alpha{{\bm{\varepsilon}}}}\enspace.$
(25a) $\displaystyle\mbox{In
particular,}\qquad\mathcal{I}_{aacc{\bm{e}}}^{\alpha\alpha\alpha\alpha{\bm{\varepsilon}}}$
$\displaystyle=\mathcal{I}_{aacc{\bm{e}}}^{\alpha\alpha\gamma\gamma{\bm{\varepsilon}}}+2\mathcal{I}_{aacc{\bm{e}}}^{\alpha\gamma\alpha\gamma{\bm{\varepsilon}}}$
(25b)
$\displaystyle\mbox{and}\qquad\mathcal{I}_{aaaa{\bm{e}}}^{\alpha\alpha\alpha\alpha{\bm{\varepsilon}}}$
$\displaystyle=3\mathcal{I}_{aaaa{\bm{e}}}^{\alpha\alpha\gamma\gamma\vec{{\bm{\varepsilon}}}}\enspace.$
(25c)

######  Lemma 2.

For any multi-index ${\bm{e}}$ not containing $a$, and for any unrestricted
multi-index ${\bm{\varepsilon}}$ of the same size,

$\sum_{a=1}^{n-1}\mathcal{I}^{\alpha\beta{\bm{\varepsilon}}}_{aa{\bm{e}}}=\delta_{\alpha\beta}\>\mathcal{I}^{{\bm{\varepsilon}}}_{\bm{e}}\enspace.$
(26)

###### Proof.

(26) follows directly from $\bm{T}\bm{T}^{\top}=\bm{I}_{n-1}$, that is
$\sum_{a=1}^{n-1}t_{a\alpha}t_{a\beta}=\delta_{\alpha\beta}$. ∎

###  3.4 The first and second moments

The computation of the first moment, which has been derived in the literature
under various strategies, is straightforward :
$\mathcal{I}_{ab}^{\alpha\alpha}=\delta_{ab}\;\mathcal{I}_{aa}^{\alpha\alpha}$,
where $\mathcal{I}_{aa}^{\alpha\alpha}$ is independent of $a$ and $\alpha$. By
(17), (3.2), (26) and joint orthogonality

$\mathcal{I}_{aa}^{\alpha\alpha}=\frac{1}{n-1}\qquad\qquad\mathbb{E}(P_{\alpha\beta})=\frac{1}{n-1}\qquad\qquad\mathbb{E}({\tt
CV})=\frac{1}{n-1}\sum_{\alpha,\beta=1}^{n-1}\lambda_{\alpha}\mu_{\beta}\enspace.$
(27)

The computation of the second moment involves four orthogonal coefficients,
namely (all super- and sub-indices in (28) are distinct)

$E:=\mathcal{I}_{aaaa}^{\alpha\alpha\alpha\alpha}\qquad\quad
F:=\mathcal{I}_{aacc}^{\alpha\alpha\alpha\alpha}=\mathcal{I}_{aaaa}^{\alpha\alpha\gamma\gamma}\qquad\quad
G:=\mathcal{I}_{aacc}^{\alpha\alpha\gamma\gamma}\quad\qquad
H:=\mathcal{I}_{aacc}^{\alpha\gamma\alpha\gamma}\enspace.$ (28)

They satisfy

$F\stackrel{{\scriptstyle(\ref{lemme1})}}{{=}}G+2H\qquad
E\stackrel{{\scriptstyle(\ref{lemme2})}}{{=}}3F\qquad(n-2)G+F\stackrel{{\scriptstyle(\ref{lemme3}),(\ref{firstorder})}}{{=}}\frac{1}{n-1}\qquad(n-2)H+F\stackrel{{\scriptstyle(\ref{lemme3}),(\ref{firstorder})}}{{=}}0$

with solution

$E=3(n-2)\kappa\>\>,\>\>F=(n-2)\kappa\>\>,\>\>G=n\kappa\>\>,\>\>H=-\kappa\>\>,\>\>\kappa=\frac{1}{(n-2)(n-1)(n+1)}$
(29)

Hence (the expression is also valid for four possibly coinciding sub-indices,
since $E=3F$)

$\mathcal{I}_{abcd}^{\alpha\alpha\alpha\alpha}=(\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})\,F$
(30)

and (for $\alpha\neq\gamma$)

$\mathcal{I}_{abcd}^{\alpha\alpha\gamma\gamma}=\delta_{ab}\delta_{cd}\,G+(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})\,H\enspace.$
(31)

As a result, performing $\delta_{\alpha\gamma}$ $\times$ (30) $+$
$(1-\delta_{\alpha\gamma})$ $\times$ (31) yields the general formula, where
the super- and the sub-indices may be distinct or not

$\mathcal{I}_{abcd}^{\alpha\alpha\gamma\gamma}=\kappa\,[n\delta_{ab}\delta_{cd}-(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})-2\delta_{\alpha\gamma}\delta_{ab}\delta_{cd}+(n-1)\delta_{\alpha\gamma}(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})]$

which finally implies, in view of (3.2) and joint orthogonality,

$\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta})=\kappa\,[n-2\delta_{\alpha\gamma}-2\delta_{\beta\delta}+2(n-1)\delta_{\alpha\gamma}\delta_{\beta\delta}]\enspace.$
(32)

Hence, by (3.2) and (10)

$\mathbb{E}({\tt
CV}^{2})=\kappa(n-1)^{3}\,[n(n-1)\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{2}\>\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{2}-2\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{2}\>\overline{\bm{\mu}^{2}}-2\overline{\bm{\lambda}^{2}}\>\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{2}+2\overline{\bm{\lambda}^{2}}\>\overline{\bm{\mu}^{2}}]\enspace.$
(33)

Substracting $\mathbb{E}^{2}({\tt CV})$ obtained in (27), substituting the
value of $\kappa$ in (29) and rearranging terms yields (7).

Expressions for the second centered moments (7) and (12) are simpler than the
corresponding quantities obtained, in the unweighted setting, by averaging
over all permutations of the $n$ objects : the latter contain additional
correction terms, as derived in Kazi-Aoual et al., (1995). See also Heo and
Ruben Gabriel, (1998), Josse et al., (2008) and Abdi, (2010).

###  3.5 The third moment

The third moment reads

$\mathbb{E}({\tt
CV}^{3})=\sum_{\alpha,\beta,\gamma,\delta,\varepsilon,\zeta=1}^{n-1}\lambda_{\alpha}\,\lambda_{\gamma}\,\lambda_{\varepsilon}\,\mu_{\beta}\,\mu_{\delta}\,\mu_{\zeta}\>\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})$
(34)

where

${\textstyle\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})=\sum_{i,j,k,l,s,t=1}^{n}v_{i\beta}v_{j\beta}v_{k\delta}v_{l\delta}v_{s\zeta}v_{t\zeta}\>\sum_{a,b,c,d,e,f=1}^{n-1}\mathcal{I}_{abcdef}^{\alpha\alpha\gamma\gamma\varepsilon\varepsilon}\>u_{ia}u_{jb}u_{kc}u_{ld}u_{se}u_{tf}}$
(35)

and involves eleven third-order orthogonal coefficients, namely (all super-
and sub-indices in (3.5) are distinct)

$\displaystyle L$
$\displaystyle:=\mathcal{I}_{aaaaaa}^{\alpha\alpha\alpha\alpha\alpha\alpha}$
$\displaystyle M$
$\displaystyle:=\mathcal{I}_{aaaacc}^{\alpha\alpha\alpha\alpha\alpha\alpha}$
$\displaystyle N$
$\displaystyle:=\mathcal{I}_{aaccee}^{\alpha\alpha\alpha\alpha\alpha\alpha}$
$\displaystyle P$
$\displaystyle:=\mathcal{I}_{aaaacc}^{\alpha\alpha\alpha\alpha\gamma\gamma}$
$\displaystyle Q$
$\displaystyle:=\mathcal{I}_{aaccaa}^{\alpha\alpha\alpha\alpha\gamma\gamma}$
$\displaystyle R$
$\displaystyle:=\mathcal{I}_{aaacac}^{\alpha\alpha\alpha\alpha\gamma\gamma}$
$\displaystyle S$
$\displaystyle:=\mathcal{I}_{aaccee}^{\alpha\alpha\alpha\alpha\gamma\gamma}$
$\displaystyle T$
$\displaystyle:=\mathcal{I}_{aceeac}^{\alpha\alpha\alpha\alpha\gamma\gamma}$
(36) $\displaystyle U$
$\displaystyle:=\mathcal{I}_{aaccee}^{\alpha\alpha\gamma\gamma\varepsilon\varepsilon}$
$\displaystyle V$
$\displaystyle:=\mathcal{I}_{acacee}^{\alpha\alpha\gamma\gamma\varepsilon\varepsilon}$
$\displaystyle W$
$\displaystyle:=\mathcal{I}_{aeceac}^{\alpha\alpha\gamma\gamma\varepsilon\varepsilon}\enspace.$

Handcrafted computations are a bit awkward, yet feasible, with the result

######  Lemma 3 (proved in the Appendix).

$\begin{split}\frac{1}{\hat{\kappa}}\>\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})&=(n^{2}+n-4)-2(n+1)(\sigma+\tau)+16(\varphi+\psi)+8\sigma\tau\\\
&-8(n-1)(\sigma\psi+\tau\varphi)+8(n-1)^{2}\varphi\psi+2(n-3)(n+3)\omega\end{split}$
(37)

where

$\hat{\kappa}=\frac{\kappa}{(n-3)(n+3)}=\frac{1}{(n-3)(n-2)(n-1)(n+1)(n+3)}$
(38)

and

$\begin{split}&\sigma=\delta_{\alpha\gamma}+\delta_{\alpha\varepsilon}+\delta_{\gamma\varepsilon}\qquad\qquad\quad\tau=\delta_{\beta\delta}+\delta_{\beta\zeta}+\delta_{\delta\zeta}\qquad\qquad\quad\omega=\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\varepsilon}\delta_{\beta\zeta}+\delta_{\gamma\varepsilon}\delta_{\delta\zeta}\\\
&\varphi=\delta_{\alpha\gamma}\delta_{\alpha\varepsilon}\delta_{\gamma\varepsilon}=\delta_{\alpha\gamma}\delta_{\alpha\varepsilon}=\delta_{\alpha\gamma}\delta_{\gamma\varepsilon}=\delta_{\alpha\varepsilon}\delta_{\gamma\varepsilon}\qquad\qquad\psi=\delta_{\beta\delta}\delta_{\beta\zeta}\delta_{\delta\zeta}=\delta_{\beta\delta}\delta_{\beta\zeta}=\delta_{\beta\delta}\delta_{\delta\zeta}=\delta_{\beta\zeta}\delta_{\delta\zeta}\enspace.\end{split}$
(39)

One can check with (32) that

$\sum_{\varepsilon=1}^{n-1}\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})=\sum_{\zeta=1}^{n-1}\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})=\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta})$

as it must. Inserting (37) in (34) and using (10) yields

$\begin{split}\frac{\mathbb{E}({\tt
CV}^{3})}{(n-1)^{4}\hat{\kappa}}=&(n^{2}+n-4)(n-1)^{2}\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{3}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{3}-6(n+1)(n-1)(\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\lambda}^{2}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{3}+\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{3}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}\overline{\bm{\mu}^{2}})+\\\
+16(\overline{\bm{\lambda}^{3}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{3}&+\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{3}\overline{\bm{\mu}^{3}}+6(n^{2}+3)\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\lambda}^{2}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}\overline{\bm{\mu}^{2}}-24(\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\lambda}^{2}}\,\overline{\bm{\mu}^{3}}+\overline{\bm{\lambda}^{3}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}\overline{\bm{\mu}^{2}})+8\overline{\bm{\lambda}^{3}}\,\overline{\bm{\mu}^{3}}\enspace.\end{split}$
(40)

The centered third moment

$\mathbb{E}({\tt CV}_{c}^{3})=\mathbb{E}(({\tt CV}-\mathbb{E}({\tt
CV}))^{3})=\mathbb{E}({\tt CV}^{3})-3\,\mathbb{E}({\tt
CV}^{2})\,\mathbb{E}({\tt CV})+2\,\mathbb{E}^{3}({\tt CV})$ (41)

finally reads, by (27), (33) and (38)

$\begin{split}\frac{\mathbb{E}({\tt
CV}_{c}^{3})}{8(n-1)^{4}\hat{\kappa}}\,=\,&4\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{3}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{3}-6(\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\lambda}^{2}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{3}+\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{3}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}\overline{\bm{\mu}^{2}})+2(\overline{\bm{\lambda}^{3}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}^{3}+\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}^{3}\overline{\bm{\mu}^{3}})+9\,\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\lambda}^{2}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}\overline{\bm{\mu}^{2}}\\\
-3(\overline{\bm{\lambda}}\overline{\bm{\lambda}^{2}}\,\overline{\bm{\mu}^{3}}+&\overline{\bm{\lambda}^{3}}\overline{\bm{\mu}}\overline{\bm{\mu}^{2}})+\overline{\bm{\lambda}^{3}}\,\overline{\bm{\mu}^{3}}=(\overline{\bm{\lambda}^{3}}-3\overline{\bm{\lambda}}\overline{\bm{\lambda}^{2}}+2\overline{\bm{\lambda}}^{3})(\overline{\bm{\mu}^{3}}-3\overline{\bm{\mu}}\overline{\bm{\mu}^{2}}+2\overline{\bm{\mu}}^{3})=\overline{\bm{\lambda}_{c}^{3}}\>\overline{\bm{\mu}_{c}^{3}}\end{split}$

thus proving (8). This exact expression for the third moment seems original,
and is considerably simpler than the corresponding expression derived by
averaging on the $n!$ object permutations (Kazi-Aoual et al.,, 1995). It
depends directly on $n$, but only indirectly on $\mathbf{f}$ through the
eigenvalue spectra. Expression (13) for the RV skewness is particularly
transparent, and elucidates the cause of the marked positive asymmetry of the
RV coefficient, often reported in the literature (see e.g. Mielke,, 1984; Heo
and Ruben Gabriel,, 1998; Josse et al.,, 2008; Zhang et al.,, 2009): plainly,
$a({\bm{\lambda}})>0$ and $a({\bm{\mu}})>0$ for typical scree plots.

###  3.6 The fourth moment

Computing $\mathbb{E}({\tt RV}^{4})$, or equivalently $\mathbb{E}({\tt
CV}^{4})$ is clearly untractable with the former pedestrian approach, and a
more systematic strategy is needed. The latter is provided by the work around
the orthogonal Weingarten functions (see Collins and Śniady,, 2006; Collins
and Matsumoto,, 2009; Matsumoto,, 2012; Collins et al.,, 2013; Mingo and
Popa,, 2013; Mingo and Speicher,, 2017).

Consider $\mathcal{P}_{2q}$, the set of all partitions of
$\\{1,2,\ldots,2q\\}$ whose all blocks are of length two, also called
pairings. There are $(2q-1)!!=(2q-1)(2q-3)\cdots 5\cdot 3$ distinct pairings.
For instance, for $q=4$

$\sigma=(13|25|46|78)\qquad\mbox{and}\qquad\tau=(15|26|34|78)$

constitute such pairings. Their join $\sigma\vee\tau$ (the finest partition
coarser than both $\sigma$ and $\tau$) is $\sigma\vee\tau=(123456|78)$.

In general, the join $\sigma\vee\tau$ of two pairings
$\sigma,\tau\in\mathcal{P}_{2q}$ is a partition made of $N(\sigma\vee\tau)$
blocks of even sizes $2l_{1},2l_{2},2l_{3},\ldots$, with $l_{1}\geq l_{2}\geq
l_{3}\ldots$ and $\sum_{c=1}^{N(\sigma\vee\tau)}l_{c}=q$. The multi-index
$\ell=(l_{1},l_{2},l_{3}\ldots)$ constitutes an integer partition of $q$
(noted $\ell\vdash q$), and defines the type $\ell(\sigma\vee\tau)$ of
$\sigma\vee\tau$.

For $q=4$, five integer partitions or types are possible, namely

$\ell=(1,1,1,1)\equiv(1^{4})\qquad\ell=(2,1,1)\qquad\ell=(2,2)\qquad\ell=(3,1)\quad\quad\ell=(4)\enspace.$

The orthogonal coefficients (23) turn out to express (Collins and Śniady,,
2006) as

$\mathcal{I}^{\bm{\omega}}_{\mathbf{a}}=\sum_{\sigma\in\mathcal{P}_{2q}}\sum_{\tau\in\mathcal{P}_{2q}}\>\delta_{\sigma}({\bm{\omega}})\>\delta_{\tau}({\mathbf{a}})\>\mbox{Wg}(\ell(\sigma\vee\tau))$
(42)

where (considering now $\sigma$ and $\tau$ as permutations exchanging the
indices belonging to the same block of two), the multi-Kronecker symbols
select the pairings $\sigma$ and $\tau$ compatible with the multi-indices, in
the sense

$\delta_{\sigma}({\bm{\omega}})=\prod_{r=1}^{q}\delta_{\omega_{\sigma(2r-1)},\omega_{\sigma(2r)}}\qquad\qquad\delta_{\tau}({\mathbf{a}})=\prod_{r=1}^{q}\delta_{a_{\tau(2r-1)},a_{\tau(2r)}}$
(43)

In other words, $\delta_{\sigma}({\bm{\omega}})=1$ if $\omega_{s}=\omega_{t}$
for each pair $(s,t)$ in $\sigma$ (which implies that all indices in
${\bm{\omega}}$ must occur an even number of times), and
$\delta_{\sigma}({\bm{\omega}})=0$ otherwise.

The quantities $\mbox{Wg}(\ell(\sigma\vee\tau))$ appearing in (42) are the
orthogonal Weingarten functions, and depend upon the dimension $d=n-1$ as
well. They have been computed up to order $q=6$ (Collins and Matsumoto,,
2009). For $q=4$ :

$\displaystyle\mbox{Wg}(1^{4})=$ $\displaystyle\>\phi\>(n-3)(n+2)(n^{2}+4n-4)$
$\displaystyle\mbox{Wg}(2,1,1)=$ $\displaystyle\>\phi\>(-n^{3}-3n^{2}+6n+4)$
$\displaystyle\mbox{Wg}(2,2)=$ $\displaystyle\>\phi\>(n^{2}+3n+14)$ (44)
$\displaystyle\mbox{Wg}(3,1)=$ $\displaystyle\>\phi\>(n-1)(2n+6)$
$\displaystyle\mbox{Wg}(4)=$ $\displaystyle\>-\phi\>(5n+1)$
$\mbox{where}\qquad\phi=\frac{1}{(n-4)(n-3)(n-2)(n-1)n(n+1)(n+3)(n+5)}\enspace.$
(45)

Substituting (42) in (3.2) yields

$\displaystyle\mathbb{E}({\tt CV}^{q})=$
$\displaystyle\sum_{\sigma\in\mathcal{P}_{2q}}\sum_{\tau\in\mathcal{P}_{2q}}\mbox{Wg}(\ell(\sigma\vee\tau))\sum_{{\bm{\alpha}}\in
A_{q}}\delta_{\sigma}({\bm{\alpha}}{\bm{\alpha}})\>\lambda_{\bm{\alpha}}\sum_{{\bm{\beta}}\in
A_{q}}\delta_{\tau}({\bm{\beta}}{\bm{\beta}})\>\mu_{\bm{\beta}}$
$\displaystyle=$
$\displaystyle\sum_{\sigma\in\mathcal{P}_{2q}}\sum_{\tau\in\mathcal{P}_{2q}}\mbox{Wg}(\ell(\sigma\vee\tau))\>\mbox{\rm
Tr}_{\sigma}(\mathbf{K}_{X})\>\mbox{\rm Tr}_{\tau}(\mathbf{K}_{Y})$ (46)
$\displaystyle=$
$\displaystyle\sum_{\sigma\in\mathcal{P}_{2q}}\sum_{\tau\in\mathcal{P}_{2q}}{\scriptstyle(n-1)^{N(\sigma\vee\sigma_{0})+N(\tau\vee\sigma_{0})}}\>\>\mbox{Wg}(\ell(\sigma\vee\tau))\prod_{c=1}^{N(\sigma\vee\sigma_{0})}\overline{{\bm{\lambda}}^{l_{c}}}\>\prod_{\tilde{c}=1}^{N(\tau\vee\sigma_{0})}\overline{{\bm{\mu}}^{l_{\tilde{c}}}}$

where the following lemmas and definitions have been used :

######  Lemma 4 (another variant of joint orthogonality, proved in the
Appendix).

$\sum_{{\mathbf{i}}\in[n]^{2q}}v_{{\mathbf{i}}{\bm{\omega}}}\sum_{{\mathbf{a}}\in
A_{2q}}\delta_{\sigma}({\mathbf{a}})\,u_{{\mathbf{i}}{\mathbf{a}}}=\delta_{\sigma}({\bm{\omega}})$

######  Lemma 5.

Consider the reference pairing
$\sigma_{0}=(12|34|\ldots|2q-1,q)\in\mathcal{P}_{2q}$, and consider the type
$\ell(\sigma\vee\sigma_{0})$, also called coset-type of $\sigma$ (see e.g.
Matsumoto,, 2012; Collins et al.,, 2013). Define

$\mbox{\rm
Tr}_{\sigma}(\mathbf{K})=\prod_{c=1}^{N(\sigma\vee\sigma_{0})}\mbox{\rm
Tr}(\mathbf{K}^{l_{c}})\qquad\qquad\mbox{\rm
tr}_{\sigma}(\mathbf{K})=\prod_{c=1}^{N(\sigma\vee\sigma_{0})}\mbox{\rm
tr}(\mathbf{K}^{l_{c}})$ (47)

Then $\sum_{{\bm{\alpha}}\in
A_{q}}\delta_{\sigma}({\bm{\alpha}}{\bm{\alpha}})\>\lambda_{\bm{\alpha}}=\mbox{\rm
Tr}_{\sigma}(\mathbf{K}_{X})$, which also reads

$\mbox{\rm
Tr}_{\sigma}(\mathbf{K}_{X})=(n-1)^{N(\sigma\vee\sigma_{0})}\;\mbox{\rm
tr}_{\sigma}(\mathbf{K}_{X})=(n-1)^{N(\sigma\vee\sigma_{0})}\;\prod_{c=1}^{N(\sigma\vee\sigma_{0})}\overline{{\bm{\lambda}}^{l_{c}}}\enspace.$
(48)

###### Proof.

Consider ${\bm{\omega}}={\bm{\alpha}}{\bm{\alpha}}\in A_{2q}$. By
construction, $\omega_{(2r-1)}=\omega_{2r}$, that is
$\delta_{\sigma_{0}}({\bm{\omega}})=1$. On the other hand, the term
$\delta_{\sigma}({\bm{\alpha}}{\bm{\alpha}})$ imposes
$\omega_{\sigma(2r-1)}=\omega_{\sigma(2r)}$. Hence all indices of
$\bm{\omega}$ in the blocks of $\sigma\vee\sigma_{0}$ (of sizes
$2l_{1},2l_{2},2l_{3}\ldots,2l_{N(\sigma\vee\sigma_{0})}$) are identical, that
is the sum on ${\bm{\alpha}}$ involves $N(\sigma\vee\sigma_{0})$ unconstrained
indices respectively repeated exactly
$(l_{1},l_{2},l_{3}\ldots,l_{N(\sigma\vee\sigma_{0})})$ times. ∎

Transforming expression (3.6) into an effective formula requires to determine,
among the $((2q-1)!!)^{2}$ pairings $(\sigma,\tau)$ entering into the sum, how
many are jointly of type $\ell(\sigma\vee\tau)$, $\ell(\sigma\vee\sigma_{0})$
and $\ell(\tau\vee\sigma_{0})$.

For $q=4$, Table 1 gives the distribution of joint counts of the
$105^{2}=11025$ pairings $(\sigma,\tau)$, among the $5^{3}=125$ possible
trivariate types. Those counts have, for lack of foreseeable analytical
approach, been mechanically computed with the help of the R package igraph
(Csardi et al.,, 2006), by the functions union() (determining the join of two
pairings coded as binary graphs) and components() (determining the join type).

Table 1: Trivariate type counts for $q=4$: each table refers to the type of
$\sigma\vee\tau$, namely $(1^{4})$, $(2,1,1)$, $(2,2)$, $(3,1)$, $(4)$ from
left to right and top to bottom. Rows refer to the type of
$\sigma\vee\sigma_{0}$, and columns to the type of $\tau\vee\sigma_{0}$.

Working with centered quantities notably simplifies the computations:

######  Lemma 6.

${\tt CV}_{c}={\tt CV}-\mathbb{E}({\tt
CV})=\sum_{\alpha=1}^{n-1}\sum_{\beta=1}^{n-1}\lambda^{c}_{\alpha}\mu^{c}_{\beta}P_{\alpha\beta}\enspace.$

###### Proof.

By (6), (17) and $P_{\alpha\bullet}=P_{\bullet\beta}=1$,

$\displaystyle\sum_{\alpha=1}^{n-1}\sum_{\beta=1}^{n-1}\lambda^{c}_{\alpha}\mu^{c}_{\beta}P_{\alpha\beta}=\sum_{\alpha\beta}(\lambda_{\alpha}-\overline{\bm{\lambda}\vphantom{\bm{\alpha}}})(\mu_{\beta}-\overline{\bm{\mu}\vphantom{\bm{\lambda}}})P_{\alpha\beta}=\sum_{\alpha\beta}\lambda_{\alpha}\mu_{\beta}P_{\alpha\beta}-(n-1)\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}$
$\displaystyle-(n-1)\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}+(n-1)\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\overline{\bm{\mu}\vphantom{\bm{\lambda}}}=\sum_{\alpha\beta}\lambda_{\alpha}\mu_{\beta}P_{\alpha\beta}-(n-1)\overline{\bm{\lambda}\vphantom{\bm{\alpha}}}\,\overline{\bm{\mu}\vphantom{\bm{\lambda}}}={\tt
CV}-\mathbb{E}({\tt CV})={\tt CV}_{c}\enspace.$

∎

Consequently, (3.6) entails

$\mathbb{E}({\tt
CV}_{c}^{q})=\sum_{\sigma\in\mathcal{P}_{2q}}\sum_{\tau\in\mathcal{P}_{2q}}{\scriptstyle(n-1)^{N(\sigma\vee\sigma_{0})+N(\tau\vee\sigma_{0})}}\>\>\mbox{Wg}(\ell(\sigma\vee\tau))\prod_{c=1}^{N(\sigma\vee\sigma_{0})}\overline{{\bm{\lambda}}_{c}^{l_{c}}}\>\prod_{\tilde{c}=1}^{N(\tau\vee\sigma_{0})}\overline{{\bm{\mu}}_{c}^{l_{\tilde{c}}}}$
(49)

in which, for $q=4$, the contributions of $\sigma$ and $\tau$ coset types
$(1^{4})$, $(2,1,1)$ and $(3,1)$, associated to
$\overline{{\bm{\lambda}}_{c}}=0$ or $\overline{{\bm{\mu}}_{c}}=0$, are
zeroed: only $(2,2)$ and $(4)$ survive, with contributions indicated by the
boxed counts in Table 1. Explicitly, the coefficient of
$\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{2}}$
in (49) is

$\displaystyle(n-1)^{4}[12\,\mbox{Wg}(1^{4})+24\,\mbox{Wg}(2,1,1)+60\,\mbox{Wg}(2,2)+48\,\mbox{Wg}(4)]$
$\displaystyle=\>$ $\displaystyle
12\,\phi\,(n-1)^{4}\,(n^{4}+n^{3}-15n^{2}-13n+98)\enspace,$

the coefficient of
$\overline{\bm{\lambda}_{c}^{4}}\>\overline{\bm{\mu}_{c}^{4}}$ is

$\displaystyle(n-1)^{2}[48\,\mbox{Wg}(1^{4})+288\,\mbox{Wg}(2,1,1)+240\,\mbox{Wg}(2,2)+768\,\mbox{Wg}(3,1)+960\,\mbox{Wg}(4)]$
$\displaystyle=\>$ $\displaystyle 48\,\phi\,(n-1)^{4}\,(n^{2}-n+2)\enspace,$

and the coefficient of
$\overline{\bm{\lambda}_{c}^{4}}\>\overline{\bm{\mu}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{2}}$
and
$\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\lambda}_{c}^{2}}\>\overline{\bm{\mu}_{c}^{4}}$
is

$\displaystyle(n-1)^{3}[96\,\mbox{Wg}(2,1,1)+48\,\mbox{Wg}(2,2)+192\,\mbox{Wg}(3,1)+240\,\mbox{Wg}(4)]$
$\displaystyle=\>$ $\displaystyle-48\,\phi\,(n-1)^{4}\,(2n^{2}-n-7)\enspace.$

The final expressions in the above follow from (3.6) and, together with (45),
prove (9). They have been further checked with the software Mathematica.
Expression (9) for the fourth moment is relatively simple, but it lacks
elegance and direct interpretation.

###  3.7 The third moment, revisited

Applying the steps of previous section for $q=3$ reveals that the contribution
of coset-types $(1,1,1)$ and $(2,1)$ for $\sigma$ or $\tau$ is zero by
consequence of centration. Hence, only the coset-types $(3)$ contribute to
(49), which is therefore simply proportional to
$\overline{\bm{\lambda}_{c}^{3}}\,\overline{\bm{\mu}_{c}^{3}}$. The
conciseness of expressions (8) and (13) is thus elucidated. The
proportionality coefficient is determined by the boxed components of Table 2
as

$(n-1)^{2}[8\,\mbox{Wg}(1^{3})+24\,\mbox{Wg}(2,1)+32\,\mbox{Wg}(3)]=\frac{8(n-1)^{3}}{(n-3)(n-2)(n+1)(n+3)}$
(50)

which is exactly expression (8), obtained much more indirectly in section 3.5.
The values of the Weingarten coefficients in (50) were obtained from Collins
and Śniady, (2006). They read with the present notations ($d=n-1$ and (38)) as

$\mbox{Wg}(1^{3})=\hat{\kappa}\>(n^{2}+n-4)\qquad\qquad\mbox{Wg}(2,1)=-\hat{\kappa}\>(n+1)\qquad\qquad\mbox{Wg}(3)=2\hat{\kappa}\enspace.$

Those coefficients coincide, in order, with the values $U$, $V$ and $W$
defined in (3.5) and determined in (51), as they must in view of (42).

In conclusion, the pedestrian approach of sections 3.3, 3.4 and 3.5 exactly
matches the systematic approach of sections 3.6 and 3.7: a circumstance both
expected and relieving, apt to boost confidence in the soundness of the
invariant orthogonal integration approach.

Table 2: Trivariate type counts for $q=3$: each table refers to the type of
$\sigma\vee\tau$, namely $(1^{3})$, $(2,1)$, $(3)$ from left to right. Rows
refer to the type of $\sigma\vee\sigma_{0}$, and columns to the type of
$\tau\vee\sigma_{0}$.

##  4 Discussion and conclusion

The weighted RV coefficient measures the similarity between two weighted
Euclidean configurations, and this contribution proposes exact expressions for
the first four moments of the RV. Considering weighted objects extends the
traditional uniform framework. It also provides precious guidance for
separating the trivial and non-trivial eigenspaces resulting form the spectral
decomposition of the standard kernels occurring in the weighted
multidimensional scaling of both configurations.

Our approach, invariant orthogonal integration, is nonparametric, and consists
in averaging the relative orientation of both configurations by performing
Haar integration on orthogonal matrices $\mathbf{T}\in\mathbb{O}_{n-1}$ acting
in the non-trivial eigenspace only. The resulting expressions are simpler and
easier to interpret than their traditional counterparts obtained by averaging
on permutation matrices $\mathbf{S}$ between $n$ objects. In view of
$\mathbf{S}\mathbf{S}^{\top}=\mathbf{I}_{n}$, permutations also do constitute
rotations, but in $\mathbb{O}_{n}$, and their undiscriminate use is
furthermore questionable in the weighted setting. Comparing the present
approach to the parametric approach, postulating multivariate normal
distribution for the object features, is left open for future investigations.

Also, our approach is object-oriented, as in traditional Data Analysis and
Machine Learning, rather than variable-oriented as in Mathematical Statistics.
Its use requires to dispose of squared Euclidean dissimilarities between
objets, possibly weighted, and some of its numerous applications (including
spatial autocorrelation and network clustering) will be illustrated in
forthcoming publications. This contribution underlines in particular the key
role played by the standard kernel, central to weighted multidimensional
scaling, and whose spectrum governs the values of the RV moments.
Correlatively, it appears that the humble scree plot should deserve more
consideration, beyond its limited role in selecting the number of spectral
dimensions: mentioning and interpreting its effective dimensionality, skewness
and fourth spectral moment could arguably become systematic in practice.

Computing the fourth RV moment did require to recourse to the Weingarten
calculus, whose apparatus, arguably demanding for the neophyte, turned out
decisive for the pursuit of our objective. One may reasonably hope that future
developments along that line will enrich the present results, replacing in
particular the mechanical computation of Tables 1 and 2 by true mathematical
arguments. However, determining the analytical, exact null distribution of
${\tt RV}$, may reveal itself out of reach: as a matter of fact, the moment
generating function (22) is an orthogonal analog of the celebrated Harish-
Chandra trace integral for the unitary group, whose analytical expression has
been determined ever since the fifties (Harish-Chandra,, 1957) (see also e.g.
Tao, (2013) and McSwiggen, (2021)). Yet, discovering a corresponding
expression for the orthogonal case, precisely, has not been achieved so far.

## Proofs

###### Proof of Lemma 1.

Let $\alpha\neq\gamma$ and consider the matrix $\tilde{\bm{T}}$ with
components

$\tilde{t}_{a\alpha}=\cos\xi\,t_{a\alpha}-\sin\xi\,t_{a\gamma}\qquad\quad\tilde{t}_{a\gamma}=\sin\xi\,t_{a\alpha}+\cos\xi\,t_{a\gamma}\qquad\quad\tilde{t}_{a\beta}=t_{a\beta}\quad\mbox{for
}\beta\neq\alpha,\gamma$

for any $a$, where $\xi$ is an arbitrary, fixed angle. Then $\tilde{\bm{T}}$
is an orthogonal matrix, as likely as $\bm{T}$, that is
$d\mu(\tilde{\bm{T}})=d\mu(\bm{T})$. To ease the notations, take
${\bm{\varepsilon}}$ and ${\bm{e}}$ in Lemmas 1 and 2, playing no active role
in what follows, as empty. Then

$\displaystyle\mathcal{I}_{abcd}^{\alpha\alpha\gamma\gamma}$
$\displaystyle=\mathbb{E}([\cos\xi\,t_{a\alpha}-\sin\xi\,t_{a\gamma}][\cos\xi\,t_{b\alpha}-\sin\xi\,t_{b\gamma}][\sin\xi\,t_{c\alpha}+\cos\xi\,t_{c\gamma}][\sin\xi\,t_{d\alpha}+\cos\xi\,t_{d\gamma}])$
$\displaystyle=2\cos^{2}\xi\sin^{2}\xi\mathcal{I}_{abcd}^{\alpha\alpha\alpha\alpha}+(\cos^{4}\xi+\sin^{4}\xi)\mathcal{I}_{abcd}^{\alpha\alpha\gamma\gamma}-2\cos^{2}\xi\sin^{2}\xi\mathcal{I}_{abcd}^{\alpha\gamma\alpha\gamma}-2\cos^{2}\xi\sin^{2}\xi\mathcal{I}_{abcd}^{\alpha\gamma\gamma\alpha}$

Multiplying the l.h.s. by $\cos^{4}\xi+\sin^{4}\xi+2\cos^{2}\xi\sin^{2}\xi=1$
and simplifying yields

$\mathcal{I}_{abcd}^{\alpha\alpha\alpha\alpha}=\mathcal{I}_{abcd}^{\alpha\alpha\gamma\gamma}+\mathcal{I}_{abcd}^{\alpha\gamma\alpha\gamma}+\mathcal{I}_{abcd}^{\alpha\gamma\gamma\alpha}\enspace.$

∎

###### Proof of Lemma 3.

Lemmas (1) and (2) entail the following relations between orthogonal
coefficients (3.5)

$\displaystyle S$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme1})}}{{=}}U+2V$ $\displaystyle
T$ $\displaystyle\stackrel{{\scriptstyle(\ref{lemme1})}}{{=}}V+2W$
$\displaystyle M$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme2})}}{{=}}3N$ $\displaystyle
P$ $\displaystyle\stackrel{{\scriptstyle(\ref{lemme2})}}{{=}}3S$
$\displaystyle R$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme2})}}{{=}}3T$ $\displaystyle
E$ $\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-2)M+L$
$\displaystyle E$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-2)P+M$
$\displaystyle F$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-3)N+2M$
$\displaystyle F$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-3)S+2Q$
$\displaystyle F$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-2)Q+M$
$\displaystyle G$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-3)U+2S$
$\displaystyle H$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-3)V+2T$
$\displaystyle 0$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-2)R+M$
$\displaystyle 0$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-2)T+N$
$\displaystyle 0$
$\displaystyle\stackrel{{\scriptstyle(\ref{lemme3})}}{{=}}(n-3)W+2T$

with solution (recall that $E,F,G,H$ in (29) are already known)

$\begin{split}L&=\frac{15(n-2)}{n+3}\,\kappa\qquad
M=\frac{3(n-2)}{n+3}\,\kappa\qquad N=\frac{n-2}{n+3}\,\kappa\qquad
P=\frac{3(n+2)}{n+3}\,\kappa\\\ Q&=\frac{n}{n+3}\,\kappa\qquad
R=\frac{-3}{n+3}\,\kappa\qquad S=\frac{n+2}{n+3}\,\kappa\qquad
T=\frac{-1}{n+3}\,\kappa\\\ U&=\frac{n^{2}+n-4}{(n-3)(n+3)}\,\kappa\qquad
V=\frac{-(n+1)}{(n-3)(n+3)}\,\kappa\qquad
W=\frac{2}{(n-3)(n+3)}\,\kappa\enspace.\end{split}$ (51)

Consider first $\alpha=\gamma=\varepsilon$, and assume the sub-indices of the
orthogonal coefficients to be matched into three distinct pairs. There are
$5\times 3=15$ such pairings, namely

$\begin{split}\mathcal{I}_{abcdef}^{\alpha\alpha\alpha\alpha\alpha\alpha}&=N\\{{\scriptstyle\delta_{ab}\delta_{cd}\delta_{ef}+\delta_{ac}\delta_{bd}\delta_{ef}+\delta_{ad}\delta_{bc}\delta_{ef}+\delta_{ae}\delta_{cd}\delta_{bf}+\delta_{af}\delta_{cd}\delta_{be}+\delta_{ab}\delta_{ce}\delta_{df}+\delta_{ab}\delta_{cf}\delta_{de}}\\\
+&{\scriptstyle\delta_{ac}\delta_{be}\delta_{df}+\delta_{ac}\delta_{bf}\delta_{de}+\delta_{ad}\delta_{ce}\delta_{bf}+\delta_{ad}\delta_{cf}\delta_{be}+\delta_{ae}\delta_{bc}\delta_{df}+\delta_{ae}\delta_{bd}\delta_{cf}+\delta_{af}\delta_{bc}\delta_{de}+\delta_{af}\delta_{bd}\delta_{ce}}\\}\enspace.\end{split}$
(52)

In (52), the first term preserves the three pairs in the reference partition
$(ab|cd|ef)$, the next six terms preserve one pair only, and the eight
remaining terms mix all pairs. It turns out that (52) also holds for
coinciding pairs in view of $M=3N$ and $L=5M$. By joint orthogonality, the sum
in (35) reads

$\mathbb{E}(P_{\alpha\beta}P_{\alpha\delta}P_{\alpha\zeta})=(1+2\delta_{\beta\delta}+2\delta_{\beta\zeta}+2\delta_{\delta\zeta}+8\delta_{\beta\delta}\delta_{\beta\zeta}\delta_{\delta\zeta})N\enspace.$
(53)

Consider now $\alpha=\gamma\neq\varepsilon$. Distinguishing between cases
preserving or not the pair $(ef)$ yields

$\begin{split}\mathcal{I}_{abcdef}^{\alpha\alpha\alpha\alpha\varepsilon\varepsilon}&=S\\{{\scriptstyle\delta_{ab}\delta_{cd}\delta_{ef}+\delta_{ac}\delta_{bd}\delta_{ef}+\delta_{ad}\delta_{bc}\delta_{ef}}\\}+T\\{{\scriptstyle\delta_{ae}\delta_{cd}\delta_{bf}+\delta_{af}\delta_{cd}\delta_{be}+\delta_{ab}\delta_{ce}\delta_{df}+\delta_{ab}\delta_{cf}\delta_{de}}\\\
&{\scriptstyle+\delta_{ac}\delta_{be}\delta_{df}+\delta_{ac}\delta_{bf}\delta_{de}+\delta_{ad}\delta_{ce}\delta_{bf}+\delta_{ad}\delta_{cf}\delta_{be}+\delta_{ae}\delta_{bc}\delta_{df}+\delta_{ae}\delta_{bd}\delta_{cf}+\delta_{af}\delta_{bc}\delta_{de}+\delta_{af}\delta_{bd}\delta_{ce}}\\}\end{split}$
(54)

which also holds for three preserved pairs since $3S+12T=M$. By joint
orthogonality, (35) reads

$\displaystyle\mbox{for
$\alpha=\gamma\neq\varepsilon$}\qquad\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})$
$\displaystyle=(1+2\delta_{\beta\delta})S+(2\delta_{\beta\zeta}+2\delta_{\delta\zeta}+8\delta_{\beta\delta}\delta_{\beta\zeta}\delta_{\delta\zeta})T$
(55) $\displaystyle\mbox{for
$\alpha=\varepsilon\neq\gamma$}\qquad\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})$
$\displaystyle=(1+2\delta_{\beta\zeta})S+(2\delta_{\beta\delta}+2\delta_{\delta\zeta}+8\delta_{\beta\delta}\delta_{\beta\zeta}\delta_{\delta\zeta})T$
(56) $\displaystyle\mbox{for
$\gamma=\varepsilon\neq\alpha$}\qquad\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})$
$\displaystyle=(1+2\delta_{\delta\zeta})S+(2\delta_{\beta\delta}+2\delta_{\beta\zeta}+8\delta_{\beta\delta}\delta_{\beta\zeta}\delta_{\delta\zeta})T\enspace.$
(57)

In the remaining case $\alpha\neq\gamma\neq\varepsilon$, the same reasoning
yield

$\begin{split}\mathcal{I}_{abcdef}^{\alpha\alpha\gamma\gamma\varepsilon\varepsilon}&=U\\{{\scriptstyle\delta_{ab}\delta_{cd}\delta_{ef}}\\}+V\\{{\scriptstyle\delta_{ac}\delta_{bd}\delta_{ef}+\delta_{ad}\delta_{bc}\delta_{ef}+\delta_{ae}\delta_{cd}\delta_{bf}+\delta_{af}\delta_{cd}\delta_{be}+\delta_{ab}\delta_{ce}\delta_{df}+\delta_{ab}\delta_{cf}\delta_{de}}\\}\\\
+W&\\{{\scriptstyle\delta_{ac}\delta_{be}\delta_{df}+\delta_{ac}\delta_{bf}\delta_{de}+\delta_{ad}\delta_{ce}\delta_{bf}+\delta_{ad}\delta_{cf}\delta_{be}+\delta_{ae}\delta_{bc}\delta_{df}+\delta_{ae}\delta_{bd}\delta_{cf}+\delta_{af}\delta_{bc}\delta_{de}+\delta_{af}\delta_{bd}\delta_{ce}}\\}\end{split}$
(58)

also valid for three preserved pairs since $U+6V+8W=N$, and finally

$\mbox{for
$\alpha\neq\gamma\neq\varepsilon$}\qquad\mathbb{E}(P_{\alpha\beta}P_{\gamma\delta}P_{\varepsilon\zeta})=U+(2\delta_{\beta\delta}+2\delta_{\beta\zeta}+2\delta_{\delta\zeta})V+8\delta_{\beta\delta}\delta_{\beta\zeta}\delta_{\delta\zeta}W\enspace.$
(59)

To ease notations, use definitions (LABEL:NewDiracs), multiply both sides of
(53) by $\varphi$, of (55) by
$\delta_{\alpha\gamma}(1-\delta_{\alpha\varepsilon})(1-\delta_{\gamma\varepsilon})=\delta_{\alpha\gamma}-\varphi$,
of (56) by $\delta_{\alpha\varepsilon}-\varphi$, of (57) by
$\delta_{\gamma\varepsilon}-\varphi$, of (59) by
$(1-\delta_{\alpha\gamma})(1-\delta_{\alpha\varepsilon})(1-\delta_{\gamma\varepsilon})=1-\sigma+2\varphi$,
and add the whole to obtain the unrestricted expression (37). ∎

###### Proof of Lemma 4.

$\displaystyle\sum_{{\mathbf{a}}\in
A_{2q}}\delta_{\sigma}({\mathbf{a}})\,u_{{\mathbf{i}}{\mathbf{a}}}=$
$\displaystyle\sum_{{\mathbf{a}}\in
A_{2q}}\delta_{\sigma}({\mathbf{a}})\,u_{{\mathbf{i}}{\mathbf{a}}}\sum_{{\mathbf{a}}\in
A_{2q}}\prod_{r=1}^{q}\delta_{a_{\sigma(2r-1)},a_{\sigma(2r)}}\,u_{{\mathbf{i}}{\mathbf{a}}}$
$\displaystyle=$
$\displaystyle\prod_{r=1}^{q}\>\sum_{a_{\sigma(2r-1)=1}}^{n-1}u_{a_{\sigma(2r-1)},i_{\sigma(2r-1)}}u_{a_{\sigma(2r-1)},i_{\sigma(2r)}}$
$\displaystyle=$
$\displaystyle\prod_{r=1}^{q}[\delta_{i_{\sigma(2r-1)},i_{\sigma(2r)}}-\sqrt{f_{i_{\sigma(2r-1)}}}\sqrt{f_{i_{\sigma(2r)}}}]\enspace.$

Summing the latter on
$\sum_{{\mathbf{i}}\in[n]^{2q}}v_{{\mathbf{i}}{\bm{\omega}}}$ yields, by the
joint orthogonality property of section 3.2

$\prod_{r=1}^{q}\>\sum_{i_{\sigma(2r-1)=1}}^{n}v_{i_{\sigma(2r-1)},\omega_{\sigma(2r-1)}}v_{i_{\sigma(2r-1)},\omega_{\sigma(2r)}}=\prod_{r=1}^{q}\delta_{\omega_{\sigma(2r-1)},\omega_{\sigma(2r)}}=\delta_{\sigma}({\bm{\omega}})\enspace.$

∎

## References

  * Abdi, (2010) Abdi, H. (2010).  Congruence: Congruence coefficient, RV coefficient, and Mantel coefficient.  Encyclopedia of research design, 3:222–229.
  * Amédée-Manesme et al., (2019) Amédée-Manesme, C.-O., Barthélémy, F., and Maillard, D. (2019).  Computation of the corrected Cornish-Fisher expansion using the response surface methodology: application to var and cvar.  Annals of Operations Research, 281(1):423–453.
  * Aubert and Lam, (2003) Aubert, S. and Lam, C. (2003).  Invariant integration over the unitary group.  Journal of Mathematical Physics, 44(12):6112–6131.
  * Bavaud, (2013) Bavaud, F. (2013).  Testing spatial autocorrelation in weighted networks: the modes permutation test.  Journal of Geographical Systems, 15(3):233–247.
  * Borg and Groenen, (2005) Borg, I. and Groenen, P. J. (2005).  Modern multidimensional scaling: Theory and applications.  Springer Science & Business Media.
  * Braun, (2006) Braun, D. (2006).  Invariant integration over the orthogonal group.  Journal of Physics A: Mathematical and General, 39(47):14581.
  * Collins and Matsumoto, (2009) Collins, B. and Matsumoto, S. (2009).  On some properties of orthogonal Weingarten functions.  Journal of Mathematical Physics, 50(11):113516.
  * Collins et al., (2013) Collins, B., McDonald, D., and Saad, N. (2013).  Compound Wishart matrices and noisy covariance matrices: Risk underestimation.  arXiv preprint arXiv:1306.5510.
  * Collins and Śniady, (2006) Collins, B. and Śniady, P. (2006).  Integration with respect to the Haar measure on unitary, orthogonal and symplectic group.  Communications in Mathematical Physics, 264(3):773–795.
  * Cortes et al., (2012) Cortes, C., Mohri, M., and Rostamizadeh, A. (2012).  Algorithms for learning kernels based on centered alignment.  The Journal of Machine Learning Research, 13:795–828.
  * Csardi et al., (2006) Csardi, G., Nepusz, T., et al. (2006).  The igraph software package for complex network research.  InterJournal, Complex Systems, 1695(5):1–9.
  * Escoufier, (1973) Escoufier, Y. (1973).  Le traitement des variables vectorielles.  Biometrics, 29(4):751–760.
  * Geisser and Greenhouse, (1958) Geisser, S. and Greenhouse, S. W. (1958).  An extension of Box’s results on the use of the $f$ distribution in multivariate analysis.  The Annals of Mathematical Statistics, 29(3):885–891.
  * Harish-Chandra, (1957) Harish-Chandra (1957).  Differential operators on a semisimple Lie algebra.  American Journal of Mathematics, pages 87–120.
  * Heo and Ruben Gabriel, (1998) Heo, M. and Ruben Gabriel, K. (1998).  A permutation test of association between configurations by means of the RV coefficient.  Communications in Statistics-Simulation and Computation, 27(3):843–856.
  * Josse et al., (2008) Josse, J., Pagès, J., and Husson, F. (2008).  Testing the significance of the RV coefficient.  Computational Statistics & Data Analysis, 53(1):82–91.
  * Kazi-Aoual et al., (1995) Kazi-Aoual, F., Hitier, S., Sabatier, R., and Lebreton, J.-D. (1995).  Refined approximations to permutation tests for multivariate inference.  Computational statistics & data analysis, 20(6):643–656.
  * Kendall and Stuart, (1977) Kendall, M. and Stuart, A. (1977).  The advanced theory of statistics. vol. 1: Distribution theory.  London: Griffin.
  * Matsumoto, (2012) Matsumoto, S. (2012).  General moments of the inverse real Wishart distribution and orthogonal Weingarten functions.  Journal of Theoretical Probability, 25(3):798–822.
  * McSwiggen, (2021) McSwiggen, C. (2021).  The Harish-Chandra integral: An introduction with examples.  L’Enseignement Mathématique, 67(3):229–299.
  * Mielke, (1984) Mielke, P. W. (1984).  Meteorological applications of permutation techniques based on distance functions.  Handbook of Statistics, 4:813–830.
  * Mingo and Popa, (2013) Mingo, J. A. and Popa, M. (2013).  Real second order freeness and Haar orthogonal matrices.  Journal of Mathematical Physics, 54(5):051701.
  * Mingo and Speicher, (2017) Mingo, J. A. and Speicher, R. (2017).  Free probability and random matrices, volume 35 of Fields Institute Monographs.  Springer.
  * Robert and Escoufier, (1976) Robert, P. and Escoufier, Y. (1976).  A unifying tool for linear multivariate statistical methods: the RV-coefficient.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 25(3):257–265.
  * Schlich, (1996) Schlich, P. (1996).  Defining and validating assessor compromises about product distances and attribute correlations.  In Data handling in science and technology, volume 16, pages 259–306. Elsevier.
  * Tao, (2013) Tao, T. (2013).  The Harish-Chandra-Itzykson-Zuber integral formula.  What’s new (blog) https://terrytao. wordpress. com/2013/02/08/the-harish-chandra-itzykson-zuber-integral-formula.
  * Worsley and Friston, (1995) Worsley, K. J. and Friston, K. J. (1995).  Analysis of fmri time-series revisited—again.  Neuroimage, 2(3):173–181.
  * Yamamoto and Kudo, (2017) Yamamoto, Y. and Kudo, S. (2017).  Probabilistic analysis of an estimator for the Frobenius norm of a matrix product.  JSIAM Letters, 9:9–12.
  * Zhang et al., (2009) Zhang, H., Tian, J., Li, J., and Zhao, J. (2009).  RV-coefficient and its significance test in mapping brain functional connectivity.  In Medical Imaging 2009: Biomedical Applications in Molecular, Structural, and Functional Imaging, volume 7262, pages 627–635. SPIE.


# Exact solutions of the space time-fractional Klein-Gordon equation with
cubic nonlinearities using some methods

Ayten Özkan  Erdoğan Mehmet Özkan

###### Abstract

Recently, finding exact solutions of nonlinear fractional differential
equations has attracted great interest. In this paper, the space time-
fractional Klein-Gordon equation with cubic nonlinearities is examined.
Firstly, suitable exact soliton solutions are formally extacted by using the
solitary wave ansatz method. Some solutions are also illustrated by the
computer simulations. Besides, the modified Kudryashov method is used to
construct exact solutions of this equation.

###### keywords:

Space time fractional Klein-Gordon equation , ansatz method , modified
Kudryashov method , exact solutions

###### MSC:

[2010] 26A33 , 35R11 , 83C15

††journal: Journal Title

##  1 Introduction

Fractional differential equations are generalization of differential equations
. In recent years, non-linear fractional differential equations (FDEs) have
achieved importance in various disciplines and have become popular. Recently,
the theory and applications of FDEs have been the focus of many studies since
they appear frequently in various applications in mathematics, physics,
biology, engineering, signal processing, systems identification, control
theory, finance, fractional dynamics, and have increasingly fascinated the
attention of many scientists. FDEs have been studied and many researchers
published books and articles in this field [1, 2, 3]. Many methods have been
introduced to obtain exact solutions of FDEs. For instance the first integral
method [4, 5, 6, 7, 8], exp-function method [9, 10, 11], $(G^{\prime}/G)$
expansion method [12, 13, 14], sub-equation method [15, 16], functional
variable method [17, 18], trial equation method [19, 20].

A dependable and powerful method called the ansatz method has been put forward
to search for traveling wave solutions of nonlinear partial differential
equations by Biswas [21, 22] . Although this method has been used by many
authors, the applications of this method are very low in nonlinear FDEs. The
installation of exact and analytical traveling wave solutions of nonlinear
FDEs is one of the most significant and basic duties in nonlinear science,
because they will characterize miscellaneous natural case such as vibrations,
solitons and finite speed distribution. The Ansatz method is one of the
efficient methods used to obtain exact soliton solutions of FDEs.

The solitary wave study has made important progress recently. In mathematics
and physics, a soliton or a solitary wave is a self-reinforcing single wave
that moves at a constant velocity, while maintaining its shape. Solitons
represent solutions of the class of largely weak nonlinear distributive
partial differential equations associated with physical systems. This field of
study has recently made a huge progress [21, 23, 24, 25, 26, 27, 28, 29, 31,
30]. In the present study, FDEs will be converted into integer-order
differential equations by fractional complex transformation, and then various
exact solutions will be obtained to determine singular soliton solutions, dark
soliton solutions and bright soliton solutions [32, 33].

One of the approaches that led to creating exact solutions of fractional
differential equations is a modified version of the Kudryashov method [34].
The modified Kudryashov method is a powerful solution method for finding exact
solutions of nonlinear partial differential equations (PDEs) in mathematical
physics and biology. This method was first applied in fractional differential
equations by Ege and Misirli [35]. Recently, this method has gained
considerable attention due to the ability of PDEs to extract new complete
solutions both in integer order and in fractional order [36, 37, 38].

Nonlinear Klein - Gordon equations have important application areas in science
and engineering such as solid state physics, nonlinear optics and quantum
field theory [39]. This equation is a relativistic field equation for scalar
particles and is a relativistic generalization of the well-known Schrödinger
equation. Despite other relative wave equations, the Klein-Gordon equation is
the most frequently studied equation in quantum field theory, since it is used
to describe particle dynamics [40]. They have been studied by many researchers
and various methods have been used to solve them. Some of these studies can be
listed as follows : Homotopy perturbation method [41], a semi-analytical
method called fractional-reduced differential transformation method with the
appropriate initial condition [42], modified Kudryashov method [43],
fractional complex transformation, $(G^{\prime}/G)$ and $(w/g)$ expansion
methods [44], the well-organized ansatz method [45], a direct analytic method
[46], the modified expanded Tanh method [47]. In this study, ansatz method and
modified Kudryashov method are applied to find out several new explicit exact
solutions of the space time-fractional Klein–Gordon equation with cubic
nonlinearities.

##  2 The modified Riemann-Liouville derivative and methodology of solution

With recent studies, it is well known that the dynamics of many physical
processes are accurately described using FDEs having different kinds of
fractional derivatives. The most popular ones are the Caputo derivative, the
Riemann-Liouville derivative and Grünwald-Letnikov derivative. A different
definition of the fractional derivative is given by Jumarie with a little
modification of the Riemann-Liouville derivative. In [48], $f:R\rightarrow R$,
$\omega\rightarrow f(\omega)$ as a continuous function (not necessarily
differentiable), the modified Riemann-Liouville derivative of order $\alpha$
is given as follows

$D^{\alpha}_{\omega}f(\omega)=\left\\{\begin{array}[]{ll}\frac{1}{\Gamma(1-\alpha)}\frac{d}{d\omega}\int^{\omega}_{0}\frac{f(\tau)-f(0)}{(\omega-\tau)^{\alpha}}d\tau&\textrm{
$,0<\alpha<1$, }\\\ \\\ (f^{(n)}(\omega))^{(\alpha-n)}&\textrm{
$,n\leq\alpha\leq n+1,\ \ n\geq 1$ }\end{array}\right.$ (2.1)

where $\Gamma(.)$ is the Gamma function. In addition, some important
properties of the fractional modified Riemann-Liouville derivative (mRLd) are
listed as follows [49]:

$D^{\alpha}_{\omega}\omega^{\gamma}=\frac{\Gamma(1+\gamma)}{\Gamma(1+\gamma-\alpha)}\textrm{
$,\gamma>1,$ }$ (2.2) $D^{\alpha}_{\omega}(c)=0\textrm{ $(c$ constant) $,$ }$
(2.3)
$D^{\alpha}_{\omega}(af(\omega)+bg(\omega))=aD^{\alpha}_{\omega}f(\omega)+bD^{\alpha}_{\omega}g(\omega),$
(2.4)

where $a\neq 0$ and $b\neq 0$ are constants.

Now, we will take into account the following nonlinear space-time FDE of the
type

$H(u,D^{\alpha}_{t}u,D^{\alpha}_{x}u,D^{2\alpha}_{tt}u,D^{2\alpha}_{xx}u,D^{\alpha}_{t}D^{\alpha}_{x}u...)=0,\
\ \textrm{ $0<\alpha<1$ }$ (2.5)

where $u$ is an unknown functions, $H$ is a polynomial of $u$ and its partial
fractional derivatives, and $\alpha$ is order of the mRLd of the function
$u=u(x,t)$.

The traveling wave transformation is

$\displaystyle u(x,t)=U(\varepsilon),$ (2.6)
$\displaystyle\varepsilon=\frac{kx^{\alpha}}{\Gamma(1+\alpha)}-\frac{ct^{\alpha}}{\Gamma(1+\alpha)},$

with $k\neq 0$ and $c\neq 0$ are constants. We use the chain rule

$\displaystyle D^{\alpha}_{t}u=\sigma_{t}\frac{\partial
U}{\partial\varepsilon}D^{\alpha}_{t}\varepsilon,$ (2.7) $\displaystyle
D^{\alpha}_{x}u=\sigma_{x}\frac{\partial
U}{\partial\varepsilon}D^{\alpha}_{x}\varepsilon,$

with $\sigma_{t},\sigma_{x}$ are sigma indexes [50] and they can be
$\sigma_{t}=\sigma_{x}=L$, where $L$ is a constant.

Substituting (2.6) and applying (2.2) and (2.7) to (2.5), we get following
nonlinear ODE

$N(U,\frac{dU}{d\varepsilon},\frac{d^{2}U}{d\varepsilon^{2}},\frac{d^{3}U}{d\varepsilon^{3}},...)=0.$
(2.8)

##  3 Modified Kudryashov method

Let the exact solution of (2.8) can be showed as follows

$U(\varepsilon)=a_{0}+a_{1}Q(\varepsilon)+...+a_{N}Q(\varepsilon)^{N},$ (3.1)

where $a_{i}$ values ($i=0,1,2,...,N)$ are arbitrary constants to be found
later, but $a_{N}\neq 0$. $Q(\varepsilon)$ has the form

$Q(\varepsilon)=\frac{1}{1+dA^{\varepsilon}}$ (3.2)

which is a solution to the Riccati equation

$Q^{\prime}(\varepsilon)=(Q^{2}(\varepsilon)-Q(\varepsilon))lnA$ (3.3)

where $d$ and $A$ are nonzero constants with $A>0$ and $A\neq 1$. N is
revealed by balancing the highest order derivative and nonlinear terms in
(2.8). Substituting (3.1) into (2.8) and comparing the results of the terms
with a series of nonlinear equations, new exact solutions will be taken for
(2.5).

##  4 Applications

###  4.1 Application of ansatz method to space time fractional Klein-Gordon
equation

We consider the space-time fractional Klein-Gordon equation of the form

$D^{2\alpha}_{tt}u-a^{2}D^{2\alpha}_{xx}+b^{2}u-\lambda u^{3}=0,\ \ (t>0,\ \
0<\alpha\leq 1),\ \ $ (4.1)

where $a,b,\lambda$ are constants [46]. The bright, dark and singular soliton
solutions will be applied to the solitary wave ansatz method. In order to
solve Eq.(4.1), using the traveling wave transformation (2.6), we obtain to an
ODE

$L^{2}(a^{2}k^{2}-c^{2})U^{\prime\prime}-b^{2}U+\lambda U^{3}=0,$ (4.2)

with $U^{\prime}=\frac{dU}{d\varepsilon}$.

####  4.1.1 The bright soliton solution

For the bright soliton solution, we let $A,k$ and, $c$ be abritrary constants.
Then suppose

$U(\varepsilon)=A\textmd{sech}^{p}(\varepsilon),$ (4.3)

where

$\varepsilon=\frac{kx^{\alpha}}{\Gamma(1+\alpha)}-\frac{ct^{\alpha}}{\Gamma(1+\alpha)}.$
(4.4)

It follows from ansatz (4.3) and (4.4) that

$\frac{d^{2}U}{d\varepsilon^{2}}=Ap^{2}\textmd{sech}^{p}(\varepsilon)-Ap(p+1)\textmd{sech}^{p+2}(\varepsilon),$
(4.5)

and

$U^{3}=A^{3}\textmd{sech}^{3p}(\varepsilon).$ (4.6)

Substituting the ansatz (4.3)-(4.6) into (4.2), the following equation is
obtained

$\displaystyle
L^{2}(a^{2}k^{2}-c^{2})Ap^{2}\textmd{sech}^{p}(\varepsilon)-L^{2}(a^{2}k^{2}-c^{2})Ap(p+1)\textmd{sech}^{p+2}(\varepsilon)-b^{2}A\textmd{sech}^{p}(\varepsilon)$
(4.7) $\displaystyle+\lambda A^{3}\textmd{sech}^{3p}(\varepsilon)=0.$

From (4.7), we suppose the exponents $p+2$ and $3p$ are equal and from that
$p$ is determined as 1. When this value is placed in (4.7), it is reduced to
the following equation

$L^{2}(a^{2}k^{2}-c^{2})A\textmd{sech}(\varepsilon)-2L^{2}(a^{2}k^{2}-c^{2})A\textmd{sech}^{3}(\varepsilon)-b^{2}A\textmd{sech}(\varepsilon)+\lambda
A^{3}\textmd{sech}^{3}(\varepsilon)=0.$ (4.8)

From (4.8), we obtain the following system of algebraic equations

$\left\\{\begin{array}[]{ll}\lambda A^{2}-2L^{2}(a^{2}k^{2}-c^{2})=0,\\\
L^{2}(a^{2}k^{2}-c^{2})-b^{2}=0.\end{array}\right.$

Solving this system, we get

$\displaystyle A=\mp\sqrt{\frac{2L^{2}(a^{2}k^{2}-c^{2})}{\lambda}},$ (4.9)
$\displaystyle c=\mp\sqrt{\frac{L^{2}a^{2}k^{2}-b^{2}}{L^{2}}}.$

Finally, we obtain the bright soliton solution for the Fractional Klein-Gordon
as follows

$u(x,t)=\mp\sqrt{\frac{2L^{2}(a^{2}k^{2}-c^{2})}{\lambda}}\textmd{sech}\Big{(}\frac{kx^{\alpha}}{\Gamma(1+\alpha)}\mp\sqrt{\frac{L^{2}a^{2}k^{2}-b^{2}}{L^{2}}}\frac{t^{\alpha}}{\Gamma(1+\alpha)}\Big{)}.$
(4.10)

Figure 1: The solution $u(x,t)$ for equation (4.10) when
$a=2,k=1,b=1,L=1,\lambda=1$.

The solution (4.10) is displayed in Figure 1, in the interval $0<x<10$ and
$0<t<1$.

####  4.1.2 The dark soliton solution

To obtain dark soliton solution , suppose that

$U(\varepsilon)=A\textmd{tanh}^{p}(\varepsilon),$ (4.11)

where

$\varepsilon=\frac{kx^{\alpha}}{\Gamma(1+\alpha)}-\frac{ct^{\alpha}}{\Gamma(1+\alpha)},$
(4.12)

which $k,c$ and $A$ are nonzero constant coefficients. From ansatz (4.11) and
(4.12), we get

$\frac{d^{2}U}{d\varepsilon^{2}}=Ap(p-1)\textmd{tanh}^{p-2}(\varepsilon)-2Ap^{2}\textmd{tanh}^{p}(\varepsilon)+Ap(p+1)\textmd{tanh}^{p+2}(\varepsilon),$
(4.13)

and

$U^{3}=A^{3}\textmd{tanh}^{3p}(\varepsilon).$ (4.14)

Thus, substituting the ansatz (4.11)-(4.14) into (4.2), it is achieved

$\displaystyle
L^{2}(c^{2}-a^{2}k^{2})[Ap(p-1)\textmd{tanh}^{p-2}(\varepsilon)-2Ap^{2}\textmd{tanh}^{p}(\varepsilon)+Ap(p+1)\textmd{tanh}^{p+2}(\varepsilon)]$
(4.15) $\displaystyle+b^{2}A\textmd{tanh}^{p}(\varepsilon)-\lambda
A^{3}\textmd{tanh}^{3p}(\varepsilon)=0.$

From (4.15), equating exponents $p$+2 and 3$p$, that gives rise to $p$=1. By
using this value, Eq. (4.15) reduces to

$L^{2}(c^{2}-a^{2}k^{2})[-2A\textmd{tanh}(\varepsilon)+2A\textmd{tanh}^{3}(\varepsilon)]+b^{2}A\textmd{tanh}(\varepsilon)-\lambda
A^{3}\textmd{tanh}^{3}(\varepsilon)=0.$ (4.16)

From (4.16), we find the algebraic system

$\displaystyle 2L^{2}(c^{2}-a^{2}k^{2})-\lambda A^{2}=0,$ (4.17)
$\displaystyle-2L^{2}(c^{2}-a^{2}k^{2})+b^{2}=0.$

Solving the system (4.17)

$\displaystyle A=\mp\sqrt{\frac{2L^{2}({c^{2}-a^{2}k^{2})}}{\lambda}},$ (4.18)
$\displaystyle c=\mp\sqrt{\frac{{b^{2}+2L^{2}a^{2}k^{2}}}{2L^{2}}}.$

Finally, we get the dark soliton solution for the Fractional Klein-Gordon as
follows:

$u(x,t)=\mp\sqrt{\frac{2L^{2}({c^{2}-a^{2}k^{2})}}{\lambda}}\ \
\textmd{tanh}\Big{(}\frac{kx^{\alpha}}{\Gamma(1+\alpha)}\mp\sqrt{\frac{{b^{2}+2L^{2}a^{2}k^{2}}}{2L^{2}}}\frac{t^{\alpha}}{\Gamma(1+\alpha)}\Big{)}.$
(4.19)

Figure 2: The solution $u(x,t)$ for equation (4.10) when
$a=2,k=1,b=1,L=1,\lambda=1$.

The solution (4.19) is displayed in Figure 2, in the interval $0<x<10$ and
$0<t<1$.

####  4.1.3 The singular soliton solution

In finding singular soliton solution we assume

$U(\varepsilon)=A\textmd{csch}^{p}(\varepsilon),$ (4.20)

with

$\varepsilon=\frac{kx^{\alpha}}{\Gamma(1+\alpha)}-\frac{ct^{\alpha}}{\Gamma(1+\alpha)},$
(4.21)

where $k,c$ and $A$ are nonzero constant coefficients. From ansatz (4.20) and
(4.21), we find

$\frac{d^{2}U}{d\varepsilon^{2}}=Ap^{2}\textmd{csch}^{p}(\varepsilon)+Ap(p+1)\textmd{csch}^{p+2}(\varepsilon),$
(4.22)

and

$U^{3}=A^{3}\textmd{csch}^{3p}(\varepsilon).$ (4.23)

Substituting ansatz (4.20)-(4.23) into (4.2), yields

$\displaystyle
L^{2}(c^{2}-a^{2}k^{2})Ap^{2}\textmd{csch}^{p}(\varepsilon)+L^{2}(c^{2}-a^{2}k^{2})Ap(p+1)\textmd{csch}^{p+2}(\varepsilon)+b^{2}A\textmd{csch}^{p}(\varepsilon)$
(4.24) $\displaystyle-\lambda A^{3}\textmd{csch}^{3p}(\varepsilon)=0.$

In (4.24), when equating exponents $p$+2 and 3$p$, leads $p$=1. Similarly
using $p=1$, equation (4.24) reduces to

$L^{2}(c^{2}-a^{2}k^{2})A\textmd{csch}(\varepsilon)+2L^{2}(c^{2}-a^{2}k^{2})A\textmd{csch}^{3}(\varepsilon)+b^{2}A\textmd{csch}(\varepsilon)-\lambda
A^{3}\textmd{csch}^{3}(\varepsilon)=0.$ (4.25)

From (4.25), we find the algebraic equation system

$\left\\{\begin{array}[]{ll}2L^{2}(c^{2}-a^{2}k^{2})-\lambda A^{2}=0,\\\
L^{2}(c^{2}-a^{2}k^{2})+b^{2}=0.\end{array}\right.$

Solving this system, we get

$\displaystyle A=\mp\sqrt{\frac{2L^{2}(c^{2}-a^{2}k^{2})}{\lambda}}\
(c^{2}-a^{2}k^{2}>0,\lambda<0),$ (4.26) $\displaystyle
c=\mp\sqrt{\frac{L^{2}a^{2}k^{2}-b^{2}}{L^{2}}}\ (L^{2}a^{2}k^{2}-b^{2}>0).$

Finally, we find the singular soliton solution for the Fractional Klein-Gordon
as follows

$u(x,t)=\mp\sqrt{\frac{2L^{2}(c^{2}-a^{2}k^{2})}{\lambda}}\textmd{csch}\Big{(}\frac{kx^{\alpha}}{\Gamma(1+\alpha)}\mp\sqrt{\frac{L^{2}a^{2}k^{2}-b^{2}}{L^{2}}}\frac{t^{\alpha}}{\Gamma(1+\alpha)}\Big{)}.$
(4.27)

Figure 3: The solution $u(x,t)$ for equation (4.27) when
$a=2,k=1,b=1,L=1,\lambda=-1$.

The solution (4.27) is displayed in Figure 3, in the interval $0<x<10$ and
$0<t<1$.

###  4.2 Application of modified Kudryashov method to space time fractional
Klein-Gordon equation

We consider the space-time fractional Klein-Gordon equation of the form (4.1).
In order to solve Eq.(4.1), using the traveling wave transformation (2.6), we
obtain to an ODE

$L^{2}(c^{2}-a^{2}k^{2})U^{\prime\prime}+b^{2}U-\lambda U^{3}=0,$ (4.28)

with $U^{\prime}=\frac{dU}{d\varepsilon}$. The balance of $U^{3}$ and
$U^{\prime\prime}$ gives $N=1$. Therefore, we have

$U(\varepsilon)=a_{0}+a_{1}Q(\varepsilon),\ \ a_{1}\neq 0.$ (4.29)

Substituting the solution (4.29) and its derivative into (4.28) gets

$\displaystyle\Big{(}2a_{1}L^{2}(c^{2}-a^{2}k^{2})(lnA)^{2}-\lambda
a_{1}^{3}\Big{)}Q^{3}(\varepsilon)-3\Big{(}a_{1}L^{2}(c^{2}-a^{2}k^{2})(lnA)^{2}+\lambda
a_{0}a_{1}^{2}\Big{)}Q^{2}(\varepsilon)$ (4.30)
$\displaystyle+\Big{(}a_{1}L^{2}(c^{2}-a^{2}k^{2})(lnA)^{2}+b^{2}a_{1}-3\lambda
a_{0}^{2}a_{1}\Big{)}Q(\varepsilon)+b^{2}a_{0}-\lambda a_{0}^{3}=0.$

Equating the coefficients of each power of $Q(\varepsilon)$ and the constant
term to zero, solving the resulting system of algebraic equations, we get the
following solutions.

Case 1:

$a_{0}=-\frac{b}{\lambda\sqrt{\frac{1}{\lambda}}},\ \
a_{1}=2b\sqrt{\frac{1}{\lambda}},\ \
c=\mp\sqrt{(lnA)^{2}a^{2}k^{2}L^{2}+2b^{2}}$ (4.31)

Substuting (4.31) into (4.29), we have

$U(\varepsilon)=-\frac{b}{\lambda\sqrt{\frac{1}{\lambda}}}+2b\sqrt{\frac{1}{\lambda}}\Big{(}\frac{1}{1+dA^{\varepsilon}}\Big{)},\
\ (\lambda>0).$ (4.32)

Finally, we obtain the exact solution of (4.1)

$u_{1}(x,t)=-\frac{b}{\lambda\sqrt{\frac{1}{\lambda}}}+2b\sqrt{\frac{1}{\lambda}}\Big{(}\frac{1}{1+dA^{\frac{kx^{\alpha}}{\Gamma(1+\alpha)}\mp\frac{\sqrt{(lnA)^{2}a^{2}k^{2}L^{2}+2b^{2}}t^{\alpha}}{\Gamma(1+\alpha)}}}\Big{)},\
\ (\lambda>0).$ (4.33)

Case 2:

$a_{0}=\frac{b}{\lambda\sqrt{\frac{1}{\lambda}}},\ \
a_{1}=-2b\sqrt{\frac{1}{\lambda}},\ \
c=\mp\sqrt{(lnA)^{2}a^{2}k^{2}L^{2}+2b^{2}}$ (4.34)

Substuting (4.34) into (4.29), we get

$U(\varepsilon)=\frac{b}{\lambda\sqrt{\frac{1}{\lambda}}}-2b\sqrt{\frac{1}{\lambda}}\Big{(}\frac{1}{1+dA^{\varepsilon}}\Big{)},\
\ (\lambda>0).$ (4.35)

Finally, we obtain the exact solution of (4.1)

$u_{2}(x,t)=\frac{b}{\lambda\sqrt{\frac{1}{\lambda}}}-2b\sqrt{\frac{1}{\lambda}}\Big{(}\frac{1}{1+dA^{\frac{kx^{\alpha}}{\Gamma(1+\alpha)}\mp\frac{\sqrt{(lnA)^{2}a^{2}k^{2}L^{2}+2b^{2}}t^{\alpha}}{\Gamma(1+\alpha)}}}\Big{)},\
\ (\lambda>0).$ (4.36)

##  5 Conclusion

In this article, the space time-fractional Klein-Gordon equation with cubic
nonlinearities is investigated for soliton and exact solutions. Complex
fractional transformation is utilized to attain the nonlinear ODE from this
equation. Bright, dark and singular soliton solutions are obtained with
solitary wave ansatz method. Moreover, some exact solutions are found with
modified Kudryashov method. The results are proof that these methods are
accurate and effective. In addition, graphs of all soliton solutions are drawn
for the appropriate coefficients.

## References

  * [1] Miller, K. S., Ross, B., An Introduction to the Fractional Calculus and Fractional Differential Equations, Wiley, New York (1993)
  * [2] Podlubny, I., Fractional Differential Equations , Academic Press, California (1999)
  * [3] Kilbas, A. A., Srivastava, H. M., Trujillo, J. J., Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam (2006)
  * [4] Ray, Saha S., New exact solutions of nonlinear fractional acoustic wave equations in ultrasound, Comput Math Appl., 71 (2016) 859-868
  * [5] Taghizadeh, N., Najand, F., M., Soltani Mohammadi V., New exact solutions of the perturbed nonlinear fractional Schrödinger equation using two reliable methods, Appl Math., 10 (2015) 139-148.
  * [6] Mirzazadeh, M., Eslami, M., Biswas, A., Solitons and periodic solutions to a couple of fractional nonlinear evolution equations, Pramana J Phys., 82 (2014) 465-476.
  * [7] Eslami, M., Rezazadeh, H.,The first integral method for Wu-Zhang system with conformable time-fractional derivative, Calcolo., 53 (2016) 475-485.
  * [8] Cenesiz, Y., Baleanu, D., Kurt, A., et al. New exact solutions of Burgers’ type equations with conformable derivative, Waves Random Complex Media., 27 (2017) 103-116.
  * [9] Güner, Ö., Bekir, A. Exact solutions of some fractional differential equations arising in mathematical biology, Int J Biomath., 8 (2015) 1550003-1 - 1550003-17.
  * [10] Bekir, A., Güner, Ö., Çevikel A., C., Fractional complex transform and exp-function methods for fractional differential equations, Abstr Appl Anal. (2013) ,Article ID 426462.
  * [11] Bekir, A., Güner, Ö., Bhrawy, A., H., et al. Solving nonlinear fractional differential equations using expfunction and (G’/G) expansion methods, Rom J Phys., 60 (2015) 360-378.
  * [12] Bin Z., (G’/ G)-expansion method for solving fractional partial differential equations in the theory of mathematical physics, Commun Theor Phys., 58 (2012) 623-630.
  * [13] Bekir, A., Güner, Ö., Exact solutions of nonlinear fractional differential equations by (G’/ G)-expansion method, Chin Phys B., 22 (2013) 110202-1 - 110202-6.
  * [14] Baleanu, D., Uǧurlu, Y., Inc M, et al. Improved (G’/ G)-expansion method for the time-fractional biological population model and Cahn.Hilliard equation, J. Comput Nonlinear Dyn., 10 (2015) 051016.1 - 051016-8.
  * [15] Aksoy, E., Çevikel A.,C., Bekir, A., Soliton solutions of (2+1)-dimensional time-fractional Zoomeron equation, Optik., 127 (2016) 6933-6942.
  * [16] Bekir, A., Aksoy, E., Application of the subequation method to some differential equations of time fractional order, Rom J Phys., 10 (2015) 054503-1 - 054503-5.
  * [17] Matinfar, M., Eslami, M., Kordy, M., The functional variable method for solving the fractional Korteweg de Vries equations and the coupled Korteweg de Vries equations, Pramana J Phys. 85 (2015) 583-592.
  * [18] Güner, Ö., Eser, D., Exact solutions of the space time fractional symmetric regularized long wave equation using different methods, Adv Math Phys., (2014) Article ID 456804.
  * [19] Bulut, H., Baskonus, H.,M., Pandir, Y., The modified trial equation method for fractional wave equation and time fractional generalized Burgers equation, Abstr Appl Anal.,(2013) Article ID 636802.
  * [20] Odabası, M., Mısırlı, E., On the solutions of the nonlinear fractional differential equations via the modified trial equation method, Math Methods Appl Sci. (2015) DOI:10.1002/mma.3533
  * [21] Biswas, A., 1-Soliton solution of the K(m,n) equation with generalized evolution, Phys. Lett. A, 372 (2008), 4601-4602.
  * [22] Krishnan, E.,V., Biswas, A,. Solutions to the Zakharov Kuznetsov equation with higher order nonlinearity by mapping and ansatz methods, Phys. Wave Phenom. 18 (2010) 256-261.
  * [23] Bekir, A., Aksoy, E. and Guner, O., Optical soliton solutions of the Long-Short-Wave interaction system, Journal of Nonlinear Optical Physics and Materials, 22(2) (2013), 1350015(11 pages).
  * [24] Bekir, A. and Guner, O., Bright and dark soliton solutions of the (3 + 1)-dimensional generalized Kadomtsev-Petviashvili equation and generalized Benjamin equation, Pramana-J. Phys., 81(2) (2013), 203-214.
  * [25] Bekir, A. and Guner, O., Topological (dark) soliton solutions for the Camassa-Holm type equations, Ocean Eng., 74 (2013), 276-279.
  * [26] Biswas, A., 1-Soliton solution of the B(m,n) equation with generalized evolution, Commun.Nonlinear Sci. Numer. Simul., 14 (2009), 3226-3229.
  * [27] Biswas, A., Optical solitons with time-dependent dispersion, nonlinearity and attenuation in a power-law media, Commun. Nonlinear Sci. Numer. Simulat., 14 (2009), 1078-1081.
  * [28] Biswas, A., and Milovic,D., Bright and dark solitons of the generalized nonlinear Schrödingers equation, Commun. Nonlinear Sci. Numer. Simulat., 15 (2010), 1473-1484.
  * [29] Biswas, A., Triki,H., Hayat, T., and Aldossary, O. M. 1-Soliton solution of the generalized Burgers equation with generalized evolution, Applied Mathematics and Computation, 217 (2011), 10289- 10294\.
  * [30] Triki,H., and Wazwaz, A. M., Bright and dark soliton solutions for a K(m,n) equation with t-dependent coefficients, Phys. Lett. A, 373 (2009), 2162-2165.
  * [31] Khalique, C. M. and Biswas,A., Optical solitons with parabolic and dual-power law nonlinearity via Lie group analysis, Journal of Electromagnetic Waves and Applications, 23(7) (2009), 963-973.
  * [32] Güner, Ö., Singular and non-topological soliton solutions for nonlinear fractional differential equations, Chinese Physics B, 24 (2015) 100201
  * [33] Mirzazadeh, M., Topological and non-topological soliton solutions to some time-fractional differential equations, Pramana-J. Phys. 85 (2015) 17-29
  * [34] Kudryashov, N. A. One method for finding exact solutions of nonlinear differential equations, Commun. Nonlinear Sci. Numer. Simulat. 17 (6) (2012) 2248–2253.
  * [35] Ege, S.M., Misirli, E. The modified Kudryashov method for solving some fractional-order nonlinear equations, Adv. Difference Equ. 2014 (2014) 135.
  * [36] Korkmaz, A. Exact solitons to (3+1) conformable time fractional jimbo Miwa, Zakharov Kuznetsov and modified Zakharov Kuznetsov equations, Commun.Theor. Phys. 67 (2017) 479–482.
  * [37] Hosseini, K., Ayati, Z., Exact solutions of space-time fractional EW and modified EW equations using Kudryashov method, Nonlinear Sci Lett A., 7 (2016) 58-66.
  * [38] Saha, R.,S., New analytical exact solutions of time fractional KdV KZK equation by Kudryashov methods, Chin Phys B., 25 (2016) 040204-1 - 040204-7.
  * [39] Wazwaz, AM. Compactons, solitons and periodic solutions for some forms of nonlinear Klein– Gordon equations, Chaos, Solitons Fractals., 28 (2006) 1005–1013.
  * [40] Biswas, A. Zony, C., Zerrad, E. Soliton perturbation theory for the quadratic nonlinear Klein–Gordon equation, Applied Mathematics and Computation, 203 (2008) 153.
  * [41] Golmankhaneh A., K., Golmankhaneh, A., Baleanu, D., On nonlinear fractional Klein-Gordon equation. Signal Processing, 91 (2011) 446-51.
  * [42] Tamsir, M., Srivastava, V., Analytical study of time-fractional order Klein-Gordon equation, Alexandria Engineering Journal 55 (2016) 561-567
  * [43] Hosseini, K., Mayeli, P., Ansari, R., Modified Kudryashov method for solving the conformable time-fractional Klein-Gordon equations with quadratic and cubic nonlinearities, Optik, 130 (2017) 737-742.
  * [44] Unsal, O., Guner, O., Bekir, A., Analytical approach for space-time fractional Klein-Gordon equation, Optik 135 (2017) 337-345
  * [45] Hosseini, K., Mayeli, P., Ansari, R., Bright and singular soliton solutions of the conformable time-fractional Klein-Gordon equations with different nonlinearities, Waves in Random and Complex Media 28 (2018) 426-434
  * [46] Çulha, S., Daṣcıog̃lu, A., Analytic solutions of the space time conformable fractional Klein Gordon equation in general form, Waves in Random and Complex Media 29 (2019) 775-790
  * [47] Shallal, M.,A.,Jabbar, H.,N., Ali, K., K., Analytic solution for the space-time fractional Klein-Gordon and coupled conformable Boussinesq equations, Results in Physics 8 (2018) 372-378
  * [48] Jumarie, G., Modified Riemann-Liouville derivative and fractional Taylor series of nondifferentiable functions further results, Comput. Math. Appl., 51 (2006) 1367-1376
  * [49] Jumarie, G., Table of some basic fractional calculus formulae derived from a modified Riemann-Liouville derivative for nondifferentiable functions, Appl. Maths. Lett., 22 (2009) 378-385
  * [50] He, J. H., Elegan,S. K., and Li, Z. B., Geometrical explanation of the fractional complex transform and derivative chain rule for fractional calculus, Physics Letters A, 376 (2012), 257-259.