section

#  A Generalization of the Pearson Correlation
to Riemannian Manifolds

P. Michl  orcid.org/0000-0002-6398-0654

###### Abstract

The increasing application of deep-learning is accompanied by a shift towards
highly non-linear statistical models. In terms of their geometry it is natural
to identify these models with Riemannian manifolds. The further analysis of
the statistical models therefore raises the issue of a correlation measure,
that in the cutting planes of the tangent spaces equals the respective Pearson
correlation and extends to a correlation measure that is normalized with
respect to the underlying manifold. In this purpose the article reconstitutes
elementary properties of the Pearson correlation to successively derive a
linear generalization to multiple dimensions and thereupon a nonlinear
generalization to principal manifolds, given by the Riemann-Pearson
Correlation.

Keywords: Nonlinear Correlation, Riemann-Pearson Correlation, Principal
Manifold

##  1 Introduction

A fundamental issue, that accompanies the analysis of multivariate data,
concerns the quantification of statistically dependency structures by
association measures. Many approaches in this direction can be traced back to
the late 19th century, where the issue was closely related to the task, to
extract laws of nature from two dimensional scatter plots. This in particular
applies to the widespread Pearson correlation coefficient.

######  Definition (Pearson Correlation).

_Let $X\colon\Omega\to\mathbb{R}$ and $Y\colon\Omega\to\mathbb{R}$ be random
variables with finite variances $\sigma_{X}^{2}$ and $\sigma_{Y}^{2}$. Then
the Pearson correlation $\rho_{X,Y}$ is defined by:_

$\rho_{X,Y}\coloneqq\frac{\mathrm{Cov}(X,\,Y)}{\sigma_{X}\sigma_{Y}}$ (1.1)

Due to its popularity and simplicity the Pearson correlation has been
generalized to a variety of different domains of application, including
generic monotonous relationships, relationships between sets of random
variables and asymmetric relationships (Zheng et al. 2010). In the purpose to
provide a generalization to smooth curves and submanifolds, that allow an
incorporation of structural assumptions, some elementary considerations have
to be taken into account, that allow a separation between the pairwise
quantification of dependencies and their global modelling.

##  2 Correlation and Regression Dilution

Pearson’s original motivation, was the regression of a straight line, that
minimizes the averaged Euclidean distance to points, that are scattered about
it (Pearson 1901, p561). Thereby his investigations were preceded by the
observation, that for a measurement series the assumed “direction of
causality” influences the estimate of the slope of the regression line.
Thereby the direction of causality is implicated by the choice of an error
model, that assumes one random variable to be error free and the other to
account for the whole observed error. Pearson empirically observed, that for
$n\in\mathbb{N}$ points, given by i.i.d. realizations
$\boldsymbol{x}\in\mathbb{R}^{n}$ of $X$ and $\boldsymbol{y}\in\mathbb{R}^{n}$
of $Y$, the least squares regression line of $\boldsymbol{y}$ on
$\boldsymbol{x}$ only equals the regression line of $\boldsymbol{x}$ on
$\boldsymbol{y}$, if all points perfectly fit on a straight line. In all other
cases, however, the slopes of the respective regression lines turned out, not
to be reciprocal and their product was found within the interval $[0,\,1)$.
This observation was decisive for Pearson’s definition of the correlation
coefficient. Thereby $\rho_{X,Y}$ is estimated by its empirical counterpart
$\rho_{x,y}$, that replaces variances by sample variances and the covariance
by the sample variance.

######  Lemma 1.

Let $X\colon\Omega\to\mathbb{R}$ and $Y\colon\Omega\to\mathbb{R}$ be random
variables with $n\in\mathbb{N}$ i.i.d. realizations
$\boldsymbol{x}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\in\mathbb{R}^{n}$.
Furthermore let $\beta_{x}\in\mathbb{R}$ denote the slope of the linear
regression of $\boldsymbol{y}$ on $\boldsymbol{x}$ and
$\beta_{y}\in\mathbb{R}$ the slope of the linear regression of
$\boldsymbol{x}$ on $\boldsymbol{y}$. Then:

$\rho_{x,y}^{2}=\beta_{x}\beta_{y}$ (2.1)

###### Proof of Lemma 1.

The following proof is based on (Kenney et al. 1962). The least squares
regression of $\boldsymbol{y}$ on $\boldsymbol{x}$ implicates, that for
regression coefficients $\alpha_{x},\beta_{x}\in\mathbb{R}$ and a normal
distributed random error$\varepsilon\coloneqq Y-(\beta_{x}X+\alpha_{x})$ the
log-likelihood of the realizations is maximized, if and only if the
$\ell^{2}$-norm of the realizations of $\varepsilon$ is minimized, such that:

${\displaystyle
SSE_{y}(\alpha_{x},\,\beta_{x})\coloneqq\sum_{i=1}^{n}(y_{i}-(\beta_{x}x_{i}+\alpha_{x}))^{2}}\to\min$
(2.2)

Since $SSE_{y}$ is a quadratic function of $\alpha_{x}$ and $\beta_{x}$ and
therefore convex, it has a unique global minimum at:

$\displaystyle\frac{\partial}{\partial\alpha_{x}}SSE_{y}$
$\displaystyle=2\sum_{i=1}^{n}(y_{i}-(\beta_{x}x_{i}+\alpha_{x}))(-x_{i})=0$
(2.3) $\displaystyle\frac{\partial}{\partial\beta_{x}}SSE_{y}$
$\displaystyle=2\sum_{i=1}^{n}(y_{i}-(\beta_{x}x_{i}+\alpha_{x}))(-1)=0$ (2.4)

By equating the coefficients, equations 2.3 and 2.4 can be rewritten as a
system of linear equations of $\alpha_{x}$ and $\beta_{x}$:

$\displaystyle\alpha_{x}n+\beta_{x}\sum_{i=1}^{n}x_{i}$
$\displaystyle=\sum_{i=1}^{n}y_{i}$ (2.5)
$\displaystyle\alpha_{x}\sum_{i=1}^{n}x_{i}+\beta_{x}\sum_{i=1}^{n}x_{i}^{2}$
$\displaystyle=\sum_{i=1}^{n}x_{i}y_{i}$ (2.6)

Consequently in matrix notation the vector
$(\alpha_{x},\,\beta_{x})^{\mathrm{T}}$ is determined by:

$\left(\begin{array}[]{c}\alpha_{x}\\\
\beta_{x}\end{array}\right)=\left(\begin{array}[]{cc}n&\sum_{i=1}^{n}x_{i}\\\
\sum_{i=1}^{n}x_{i}&\sum_{i=1}^{n}x_{i}^{2}\end{array}\right)^{-1}\left(\begin{array}[]{c}\sum_{i=1}^{n}y_{i}\\\
\sum_{i=1}^{n}x_{i}y_{i}\end{array}\right)$ (2.7)

Let $\overline{x},\,\overline{y}$ respectively denote the _sample means_. Then
by calculating the matrix inverse, the slope $\beta_{x}$ equates to:

$\beta_{x}=\left(\sum_{i=1}^{n}x_{i}y_{i}-n\overline{x}\overline{y}\right)\left(\sum_{i=1}^{n}x_{i}^{2}-n\overline{x}^{2}\right)^{-1}$
(2.8)

Thereupon by substituting the _sample variance_ :

$\displaystyle\sigma_{x}^{2}$
$\displaystyle\coloneqq\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}$
$\displaystyle=\frac{1}{n}\left(\sum_{i=1}^{n}x_{i}^{2}-n\overline{x}^{2}\right)$
(2.9)

And the _sample covariance_ :

$\displaystyle\mathrm{Cov}(\boldsymbol{x},\,\boldsymbol{y})$
$\displaystyle\coloneqq\frac{1}{n}\sum_{i=1}^{n}(x_{i}-\overline{x})(y_{i}-\overline{y})$
$\displaystyle=\frac{1}{n}\left(\sum_{i=1}^{n}x_{i}y_{i}-n\overline{x}\overline{y}\right)$
(2.10)

It follows from equation 2.8, that:

$\beta_{x}=\frac{\mathrm{Cov}(\boldsymbol{x},\,\boldsymbol{y})}{\sigma_{x}^{2}}$
(2.11)

Conversely the slope $\beta_{y}$ of the linear regression of $\boldsymbol{x}$
on $\boldsymbol{y}$ mutatis mutandis equates to:

$\beta_{y}=\frac{\mathrm{Cov}(\boldsymbol{y},\,\boldsymbol{x})}{\sigma_{y}^{2}}$
(2.12)

By the symmetry of $\mathrm{Cov}$ it the follows, from equations 2.11 and 2.12
that:

$\beta_{x}\beta_{y}=\frac{\mathrm{Cov}(\boldsymbol{x},\,\boldsymbol{y})^{2}}{\sigma_{x}^{2}\sigma_{y}^{2}}=\rho_{xy}^{2}$
(2.13)

∎

Lemma 1 shows, that $\rho_{x,y}$ may be regarded as the geometric mean of the
regression slopes $\beta_{x}$ and $\beta_{y}$, where $\beta_{x}$ and
$\beta_{y}$ respectively describe the causal relationships $X\to Y$ and $Y\to
X$. Thereby $X$ and $Y$ respectively are treated as error free regressor
variables to predict the corresponding response variable, that captures the
overall error. The mutual linear relationship $X\leftrightarrow Y$ is then
described by a regression line, that equally treats errors in both variables.
As an immediate consequence of this symmetry it follows, that this total least
squares regression line is unique, and its slope $\beta_{x}^{\star}$, that
describes $\boldsymbol{y}$ by $\boldsymbol{x}$ is reciprocal to the slope
$\beta_{y}^{\star}$, that describes $\boldsymbol{x}$ by $\boldsymbol{y}$ such
that $\beta_{x}^{\star}\beta_{y}^{\star}=1$. In this sense $\beta_{x}$ and
$\beta_{y}$ may be regarded as biased estimations of $\beta_{x}^{\star}$ and
$\beta_{y}^{\star}$. Thereby the bias generally is known as “regression
dilution” or “regression attenuation”. For the case that both errors are
independent and normal distributed, this bias can be corrected by a prefactor,
that incorporates the error of the respective regressor variable. An
application of this correction to lemma 1 then shows, that $\rho_{X,Y}$ has a
consistent estimations by the sample variances of $\boldsymbol{x}$ and
$\boldsymbol{y}$ and the variances of their respective errors
$\varepsilon_{X}$ and $\varepsilon_{Y}$.

######  Proposition 2.

Let $X\colon\Omega\to\mathbb{R}$ and $Y\colon\Omega\to\mathbb{R}$ be random
variables with $n\in\mathbb{N}$ i.i.d. realizations
$\boldsymbol{x}\in\mathbb{R}^{n}$ and $\boldsymbol{y}\in\mathbb{R}^{n}$ and
random errors $\varepsilon_{X}\sim\mathcal{N}(0,\,\eta_{X}^{2})$ and
$\varepsilon_{Y}\sim\mathcal{N}(0,\,\eta_{Y}^{2})$. Then:

$\rho_{x,y}^{2}\overset{{\scriptstyle
P}}{\to}\left(1-\frac{\eta_{X}^{2}}{\sigma_{x}^{2}}\right)\left(1-\frac{\eta_{Y}^{2}}{\sigma_{y}^{2}}\right),\>\text{for
}n\to\infty$ (2.14)

###### Proof of Proposition 2.

Let $\beta_{x}\in\mathbb{R}$ be the slope of the ordinary least squares (OLS)
regression line of $\boldsymbol{y}$ on $\boldsymbol{x}$, where
$\boldsymbol{x}$ is assumed to realize $X$ with a normal distributed random
error $\varepsilon_{X}\sim\mathcal{N}(0,\,\eta_{X}^{2})$. Then $X$ decomposes
into (i) an unobserved error free regressor variable $X^{\star}$ and (ii) the
random error $\varepsilon_{X}$, such that:

$X\sim X^{\star}+\varepsilon_{X}$ (2.15)

With respect to this decomposition, the slope $\beta_{x}^{\star}$ of the total
least squares (TLS) regression line, that also considers $\varepsilon_{X}$,
then is identified by the slope of the OLS regression of $\boldsymbol{y}$ on
$\boldsymbol{x}^{\star}$, where $\boldsymbol{x}^{\star}$ realizes $X^{\star}$.
Thereupon let $\sigma_{x^{*}}^{2}$ be the empirical variance of
$\boldsymbol{x}^{\star}$, then according to (Snedecor et al. 1967) it follows,
that:

$\beta_{x}\overset{{\scriptstyle
P}}{\to}\frac{\sigma_{x^{*}}^{2}}{\sigma_{x^{*}}^{2}+\eta_{X}^{2}}\beta_{x}^{\star},\>\text{for
}n\to\infty$ (2.16)

Since furthermore $\varepsilon_{X}$ by definition is statistically independent
from $X^{\star}$, it can be concluded, that:

$\displaystyle\sigma_{x}^{2}$
$\displaystyle=\mathrm{Var}(X^{\star}+\varepsilon_{X})$
$\displaystyle=\mathrm{Var}(X^{\star})+\mathrm{Var}(\varepsilon_{X})=\sigma_{x^{*}}^{2}+\eta_{X}^{2}$
(2.17)

Such that:

$\frac{\sigma_{x^{*}}^{2}}{\sigma_{x^{*}}^{2}+\eta_{X}^{2}}\stackrel{{\scriptstyle{\scriptstyle\text{\ref{eq:Correlation
Reliability, proof
3}}}}}{{=}}\frac{\sigma_{x}^{2}-\eta_{X}^{2}}{\sigma_{x}^{2}}=1-\frac{\eta_{X}^{2}}{\sigma_{x}^{2}}$
(2.18)

And therefore by equation 2.16 that:

$\beta_{x}\overset{{\scriptstyle
P}}{\to}\left(1-\frac{\eta_{X}^{2}}{\sigma_{x}^{2}}\right)\beta_{x}^{\star},\>\text{for
}n\to\infty$ (2.19)

Conversely let now $\beta_{y}\in\mathbb{R}$ be the the slope of the OLS
regression of $\boldsymbol{x}$ on $\boldsymbol{y}$, where $\boldsymbol{y}$ is
assumed to realize $Y$ with a random error
$\varepsilon_{X}\sim\mathcal{N}(0,\,\eta_{Y}^{2})$. Then also the corrected
slope $\beta_{y}^{\star}$ mutatis mutandis satisfies the relation given by
equation 2.19 and by the representation of $\rho_{x,y}$, as given by lemma 1,
it then can be concluded, that:

$\displaystyle\rho_{x,y}^{2}$
$\displaystyle\stackrel{{\scriptstyle{\scriptstyle\text{\ref{lem:Correlation
Lin. Regr. Slopes}}}}}{{=}}\beta_{x}\beta_{x}$ (2.20)
$\displaystyle\overset{{\scriptstyle
P}}{\to}\left(1-\frac{\eta_{X}^{2}}{\sigma_{x}^{2}}\right)\left(1-\frac{\eta_{Y}^{2}}{\sigma_{y}^{2}}\right)\beta_{x}^{\star}\beta_{y}^{\star},\>\text{for
}n\to\infty$

The proposition then follows by the uniqueness of the total least squares
regression line for known variances $\eta_{X}^{2}$ and $\eta_{Y}^{2}$, such
that:

$\beta_{y}^{\star}=\frac{1}{\beta_{x}^{\star}}$

∎

##  3 A Generalization to linear Principal Manifolds

Within the same publication, in which Pearson introduced the correlation
coefficient, he also developed a structured approach that determines the
straight line, that minimizes the Euclidean distance (Pearson 1901, p563). His
method, which later received attribution as the method of Principal Component
Analysis (PCA), however, even went further and allowed a canonical
generalization of the problem in the following sense: For $d\in\mathbb{N}$ let
$\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be a multivariate random vector
and for $n\in\mathbb{N}$ let $\boldsymbol{x}\in\mathbb{R}^{n\times d}$ be an
i.i.d. realization of $\boldsymbol{X}$. Then for any given $k\in\mathbb{N}$
with $k\leq d$ the goal is, to determine an affine linear subspace
$L\subseteq\mathbb{R}^{d}$ of dimension $k$, that minimizes the summed
Euclidean distance to $\boldsymbol{x}$. In order to solve this problem, the
fundamental idea of Pearson was, to transfer the principal axis theorem from
ellipsoids to multivariate Gaussian distributed random vectors. Thereupon,
however, the method also can be formulated with respect to generic elliptical
distributions.

######  Definition (Elliptical Distribution).

_For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be a
random vector. Then $\boldsymbol{X}$ is elliptically distributed, iff there
exists a random vector $\boldsymbol{S}\colon\Omega\to\mathbb{R}^{k}$ with
$k\leq d$, which distribution is invariant to rotations, a matrix
$A\in\mathbb{R}^{d\times k}$ of rank $k$ and a vector
$\boldsymbol{b}\in\mathbb{R}^{d}$, such that:_

$\boldsymbol{X}\sim A\boldsymbol{S}+\boldsymbol{b}$ (3.1)

Consequently a random vector $\boldsymbol{X}$ is elliptically distributed, if
it can be represented by an affine transformation of a radial symmetric
distributed random vector $\boldsymbol{S}$. The decisive property, that
underpins the choice of elliptical distributions, lies within their
coincidence of linear and statistical dependencies, which allows to decompose
$\boldsymbol{X}$ in statistically independent components by a linear
decomposition. This property allows, to substantiate the multidimensional
“linear fitting problem” with respect to an orthogonal projection.

######  Proposition 3.

For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be an
elliptically distributed random vector, $L\subseteq\mathbb{R}^{d}$ an affine
linear subspace of $\mathbb{R}^{d}$ and $\pi_{L}$ the orthogonal projection of
$\mathbb{R}^{d}$ onto $L$. Then the following statements are equivalent:

_(i)_ $L$ minimizes the Euclidean distance to $\boldsymbol{X}$

_(ii)_ $\mathbb{E}(\boldsymbol{X})\in L$ and $L$ maximizes the variance
$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$

###### Proof.

Let $\boldsymbol{Y}_{L}\coloneqq\boldsymbol{X}-\pi_{L}(\boldsymbol{X})$, then
the Euclidean distance between $\boldsymbol{X}$ and $L$ can be written:

$\displaystyle d(\boldsymbol{X},\,\pi_{L}(\boldsymbol{X}))^{2}$
$\displaystyle=\mathbb{E}(\|(\boldsymbol{X}-\pi_{L}(\boldsymbol{X}))\|_{2}^{2})$
(3.2) $\displaystyle=\mathbb{E}(\boldsymbol{Y}_{L}^{2})$

This representation can furthermore be decomposed by using the algebraic
formula for the variance:

$\mathbb{E}(\boldsymbol{Y}_{L}^{2})=\mathrm{Var}(\boldsymbol{Y}_{L})+\mathbb{E}(\boldsymbol{Y}_{L})^{2}$
(3.3)

Let now be $\boldsymbol{Y}_{L}^{\perp}\coloneqq\pi_{L}(\boldsymbol{X})$, then
$\boldsymbol{X}=\boldsymbol{Y}_{L}+\boldsymbol{Y}_{L}^{\perp}$ and
$\boldsymbol{Y}_{L}$ and $\boldsymbol{Y}_{L}^{\perp}$ are uncorrelated, such
that:

$\displaystyle\mathrm{Var}(\boldsymbol{X})$
$\displaystyle=\mathrm{Var}(\boldsymbol{Y}_{L}+\boldsymbol{Y}_{L}^{\perp})$
(3.4)
$\displaystyle=\mathrm{Var}(\boldsymbol{Y}_{L})+\mathrm{Var}(\boldsymbol{Y}_{L}^{\perp})$

From equations 3.2, 3.3 and 3.4 it follows, that:

$d(\boldsymbol{X},\,\pi_{L}(\boldsymbol{X}))^{2}=\mathrm{Var}(\boldsymbol{X})-\mathrm{Var}(\boldsymbol{Y}_{L}^{\perp})+\mathbb{E}(\boldsymbol{Y}_{L})^{2}$
(3.5)

Consequentially the Euclidean distance is minimized, if and only if the right
side of equation 3.5 is minimized. The first term
$\mathrm{Var}(\boldsymbol{X})$, however, does not depend on $L$ and since
$\boldsymbol{X}$ is elliptically distributed, the linear independence of
$\boldsymbol{Y}_{L}^{\perp}$ and $\boldsymbol{Y}_{L}$ is sufficient for
statistically independence. It follows, that the Euclidean distance is
minimized, if and only if: (1) The term $\mathbb{E}(\boldsymbol{Y}_{L})^{2}$
is minimized and (2) the term $\mathrm{Var}(\boldsymbol{Y}_{L}^{\perp})$ is
maximized. Concerning (1) it follows, that:

$\displaystyle\mathbb{E}(\boldsymbol{Y}_{L})^{2}$
$\displaystyle=\mathbb{E}(\boldsymbol{X}-\pi_{L}(\boldsymbol{X}))^{2}$
$\displaystyle=(\mathbb{E}(\boldsymbol{X})-\pi_{L}(\mathbb{E}(\boldsymbol{X})))^{2}$

Therefore the term $\mathbb{E}(\boldsymbol{Y}_{L})^{2}$ is minimized, if and
only if $\pi_{L}(\mathbb{E}(\boldsymbol{X}))=\mathbb{E}(\boldsymbol{X})$,
which in turn means that $\mathbb{E}(\boldsymbol{X})\in L$. Concerning (2),
the proposition immediately follows by the definition of
$\boldsymbol{Y}_{L}^{\perp}$. ∎

Let now be $k\leq d$. In order to derive an affine linear subspace
$L\subseteq\mathbb{R}^{d}$ that minimizes the Euclidean distance to
$\boldsymbol{X}$, proposition … states, that is suffices to provide an $L$
which (1) is centred in $\boldsymbol{X}$, such that
$\mathbb{E}(\boldsymbol{X})\in L$, and (2) maximizes the variance of the
projection. In order to maximize $\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$,
however, it is beneficial to give a further representation.

######  Lemma 4.

For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be an
elliptically distributed random vector, $L\subseteq\mathbb{R}^{d}$ an affine
linear subspace of $\mathbb{R}^{d}$, which for an $k\leq d$, a vector
$\boldsymbol{v}\in\mathbb{R}^{d}$ and an orthonormal basis
$\boldsymbol{u}_{1},\,\ldots\,,\,\boldsymbol{u}{}_{k}\in\mathbb{R}^{d}$ is
given by:

$L=\boldsymbol{v}+\bigoplus_{i=1}^{k}\mathbb{R}\boldsymbol{u}_{i}$

Let further be $\pi_{L}\colon\mathbb{R}^{d}\to L$ the orthogonal projection of
$\mathbb{R}^{d}$ onto $L$. Then the variance of the projection is given by:

$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))=\sum_{i=1}^{k}\boldsymbol{u}_{i}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{u}_{i}$

###### Proof.

Let $L^{{}^{\prime}}\coloneqq L+\mathbb{E}(\boldsymbol{X})-\boldsymbol{v}$,
then the orthogonal projection $\pi_{L^{{}^{\prime}}}(\boldsymbol{X})$
decomposes into individual orthogonal projections to the respective basis
vectors, such that:

$\pi_{L^{{}^{\prime}}}(\boldsymbol{X})=\mathbb{E}(\boldsymbol{X})+\sum_{i=1}^{k}\left\langle\boldsymbol{X}-\mathbb{E}(\boldsymbol{X}),\,\boldsymbol{u}_{i}\right\rangle\boldsymbol{u}_{i}$
(3.6)

Let
$\hat{X}_{i}\coloneqq\left\langle\boldsymbol{X},\,\boldsymbol{u}_{i}\right\rangle\boldsymbol{u}_{i}$,
for $i\in\\{1,\,\ldots\,,\,k\\}$. The total variance of this projection is
then given by:

$\displaystyle\mathrm{Var}(\pi_{L^{{}^{\prime}}}(\boldsymbol{X}))$
$\displaystyle\stackrel{{{}_{\ref{eq:L1.2.1}}}}{{=}}\mathrm{Var}\left(\mathbb{E}(\boldsymbol{X})+\sum_{i=1}^{k}\hat{X}_{i}-\sum_{i=1}^{k}\mathbb{E}(\hat{X}_{i})\right)$
(3.7) $\displaystyle=\mathrm{Var}\left(\sum_{i=1}^{k}\hat{X}_{i}\right)$

Since the random variables $\hat{X}_{i}$ by definition are uncorrelated, the
algebraic formula for the variance can be used to decompose the variance:

$\displaystyle\mathrm{Var}\left(\sum_{i=1}^{k}\hat{X}_{i}\right)$
$\displaystyle=\sum_{i=1}^{k}\mathrm{Var}(\hat{X}_{i})$ (3.8)

By equating the term $\mathrm{Var}(\hat{X}_{i})$, for
$i\in\\{1,\,\ldots\,,\,k\\}$ it follows, that:

$\displaystyle\mathrm{Var}(\hat{X}_{i})$
$\displaystyle=\mathrm{Var}\left(\left\langle\boldsymbol{X},\,\boldsymbol{u}_{i}\right\rangle\boldsymbol{u}_{i}\right)$
(3.9)
$\displaystyle=\mathrm{Var}\left(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{u}_{i}\right)\boldsymbol{u}_{i}^{2}$
$\displaystyle=\mathrm{Var}\left(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{u}_{i}\right)$

And furthermore by introducing the covariance matrix
$\mathrm{Cov}(\boldsymbol{X})$:

$\displaystyle\mathrm{Var}\left(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{u}_{i}\right)$
$\displaystyle\stackrel{{{}_{\text{def}}}}{{=}}\mathbb{E}\left((\boldsymbol{X}^{\mathrm{T}}\boldsymbol{u}_{i})^{\mathrm{T}}(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{u}_{i})\right)$
(3.10)
$\displaystyle=\boldsymbol{u}_{i}^{\mathrm{T}}\mathbb{E}\left(\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X}\right)\boldsymbol{u}_{i}$
$\displaystyle\stackrel{{{}_{\text{def}}}}{{=}}\boldsymbol{u}_{i}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{u}_{i}$

Summarized the equations 3.7, 3.8, 3.9 and 3.10 provide a representation for
the variance of the projection to $L^{{}^{\prime}}$:

$\mathrm{Var}(\pi_{L^{{}^{\prime}}}(\boldsymbol{X}))=\sum_{i=1}^{k}\boldsymbol{u}_{i}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{u}_{i}$

Finally the total variance of the projection is invariant under translations
of $L$, such that:

$\displaystyle\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$
$\displaystyle=\mathrm{Var}(\pi_{L}(\boldsymbol{X})+\mathbb{E}(\boldsymbol{X})-\boldsymbol{v})$
(3.11) $\displaystyle=\mathrm{Var}(\pi_{L^{{}^{\prime}}}(\boldsymbol{X}))$
$\displaystyle\stackrel{{\scriptstyle\ref{eq:L1.2.3}}}{{=}}\sum_{i=1}^{k}\boldsymbol{u}_{i}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{u}_{i}$

∎

Lemma 4, shows, that for elliptically distributed random vectors $\mathbf{X}$
the best fitting linear subspaces are completely determined by the expectation
$\mathbb{E}(\boldsymbol{X})$ and the covariance matrix
$\mathrm{Cov}(\boldsymbol{X})$. On this point it is important to notice, that
the covariance matrix is symmetric, which allows its diagonalization with
regard to real valued Eigenvalues.

######  Lemma 5.

For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be an
elliptically distributed random vector, $L\subseteq\mathbb{R}^{d}$ an affine
linear subspace of $\mathbb{R}^{d}$, which for an $k\leq d$, a vector
$\boldsymbol{v}\in\mathbb{R}^{d}$ and an orthonormal basis
$\boldsymbol{u}_{1},\,\ldots\,,\,\boldsymbol{u}{}_{k}\in\mathbb{R}^{d}$ is
given by:

$L=\boldsymbol{v}+\bigoplus_{i=1}^{k}\mathbb{R}\boldsymbol{u}_{i}$

Let further be $\pi_{L}\colon\mathbb{R}^{d}\to L$ the orthogonal projection of
$\mathbb{R}^{d}$ onto $L$, as well as
$\lambda_{1},\,\ldots\,,\,\lambda_{d}\in\mathbb{R}$ the eigenvalues of
$\mathrm{Cov}(\boldsymbol{X})$. Then there exist numbers
$a_{1},\,\ldots\,,\,a_{d}\in[0,\,1]$ with $\sum_{i=1}^{d}a_{i}=k$, such that:

$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))=\sum_{i=1}^{d}\lambda_{i}a_{i}$

###### Proof.

From lemma 4 it follows, that:

$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))=\sum_{j=1}^{k}\boldsymbol{u}_{j}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{u}_{j}$

Since the covariance matrix $\mathrm{Cov}(\boldsymbol{X})$ is a symmetric
matrix, there exists an orthonormal basis transformation matrix
$S\in\mathbb{R}^{d\times d}$ and a diagonal matrix $D\in\mathbb{R}^{d\times
d}$, such that $\mathrm{Cov}(\boldsymbol{X})=S^{T}DS$. Then the variance
$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$ has a decomposition, given by:

$\displaystyle\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$
$\displaystyle=\sum_{i=1}^{k}\boldsymbol{u}_{j}{}^{\mathrm{T}}S^{\mathrm{T}}DS\boldsymbol{u}_{j}$
$\displaystyle=\sum_{i=1}^{k}(S\boldsymbol{u}_{j})^{\mathrm{T}}DS\boldsymbol{u}_{j}$

For $j\in\\{1,\,\ldots\,,\,k\\}$ let now $\boldsymbol{c}_{j}\coloneqq
S\boldsymbol{u}_{j}$ and for $i\in\\{1,\,\ldots\,,\,n\\}$ let the number
$a_{i}\in\mathbb{R}$ be defined by:

$a_{i}\coloneqq\sum_{j=1}^{k}(\boldsymbol{c}_{j}{}_{i})^{2}$

Then according to Lemma 4 the variance $\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$
can be decomposed:

$\displaystyle\sum_{i=1}^{k}\boldsymbol{c}_{j}{}^{T}D\boldsymbol{c}_{j}$
$\displaystyle=\sum_{j=1}^{k}\sum_{i=1}^{d}\boldsymbol{c}_{j}{}_{i}\lambda_{i}\boldsymbol{c}_{j}{}_{i}$
$\displaystyle=\sum_{i=1}^{d}\lambda_{i}\sum_{j=1}^{k}(\boldsymbol{c}_{j}{}_{i})^{2}$
$\displaystyle=\sum_{i=1}^{d}\lambda_{i}a_{i}$

Furthermore since $\boldsymbol{u}_{1},\,\ldots\,,\,\boldsymbol{u}{}_{k}$ is an
orthonormal basis and $S$ an orthonormal matrix it follows that also
$\boldsymbol{c}_{1},\,\ldots\,,\,\boldsymbol{c}{}_{k}$ is an orthonormal
basis. Consequentially for $i\in\\{1,\,\ldots\,,\,d\\}$ it holds, that:

$a_{i}=\sum_{j=1}^{k}(\boldsymbol{c}_{j}{}_{i})^{2}\leq\sum_{j=1}^{k}\left\|\boldsymbol{c}_{j}{}_{i}\right\|_{2}\leq
1$

And furthermore by its definition it follows, that $a_{i}\geq 0$, such that
$a_{i}\in[0,\,1]$. Besides this the sum over all $a_{i}$ equates to:

$\displaystyle\sum_{i=1}^{d}a_{i}$
$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{k}(\boldsymbol{c}_{j}{}_{i})^{2}$
$\displaystyle=\sum_{j=1}^{k}\boldsymbol{c}_{j}^{T}\boldsymbol{c}_{j}=k$

∎

With reference to the principal axis transformation, the eigenvectors of the
covariance matrix are then termed principal components and affine linear
subspaces of the embedding space as linear principal manifolds.

######  Definition (Linear Principal Manifold).

_For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be a
random vector. Then a vector $\boldsymbol{c}\in\mathbb{R}^{d}$ with
$\boldsymbol{c}\neq 0$ is a principal component for $\boldsymbol{X}$, iff
there exists an $\lambda\in\mathrm{\mathbb{R}}$, such that :_

$\mathrm{Cov}(\boldsymbol{X})\cdot\boldsymbol{c}=\lambda\boldsymbol{c}$ (3.12)

_Furthermore let $L\subseteq\mathbb{R}^{d}$ be an affine linear subspace of
$\mathbb{R}^{d}$ with dimension $k\leq d$. Then $L$ is a linear
_$k$-_principal manifold for $\boldsymbol{X}$, if there exists a set set
$\boldsymbol{c}_{1},\,\ldots\,,\,\boldsymbol{c}_{k}$ of linear independent
principal components for $\boldsymbol{X}$, such that:_

$L=\mathbb{E}(\boldsymbol{X})+\bigoplus_{i=1}^{k}\mathbb{R}\boldsymbol{c}_{i}$

_Then $L$ is termed maximal, iff the sum of the Eigenvalues
$\lambda_{1},\,\ldots\,,\,\lambda_{k}$, that correspond to the principal
components $\boldsymbol{c}_{1},\,\ldots\,,\,\boldsymbol{c}_{k}$ is maximal. _

######  Proposition 6.

For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be an
elliptically distributed random vector and $L\subseteq\mathbb{R}^{d}$ an
affine linear subspace. Then the following statements are equivalent:

_(i)_ $L$ minimizes the Euclidean distance to $\boldsymbol{X}$

_(ii)_ $L$ is a maximal linear principal manifold for $\boldsymbol{X}$

###### Proof.

“$\Longrightarrow$” Let $\pi_{L}\colon\mathbb{R}^{d}\hookrightarrow L$ denote
the orthogonal projection of $\mathbb{R}^{d}$ onto $L$. Then according to
proposition 3 $L$ minimizes the averaged Euclidean distance to
$\boldsymbol{X}$, if and only if (i) $\mathbb{E}(\boldsymbol{X})\in L$ and
(ii) $L$ maximizes the variance $\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$. In
particular (i) is satisfied, if and only if an orthonormal basis
$\boldsymbol{u}_{1},\,\ldots\,,\,\boldsymbol{u}{}_{k}\in\mathbb{R}^{d}$ can be
chosen, such that:

$L=\mathbb{E}(\boldsymbol{X})+\bigoplus_{i=1}^{k}\mathbb{R}\boldsymbol{u}_{i}$

Then according to Lemma 5 there exist numbers
$a_{1},\,\ldots\,,\,a_{d}\in[0,\,1]$ with $\sum_{i=1}^{d}a_{i}=k$, such that:

$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))=\sum_{i=1}^{d}\lambda_{i}a_{i}$

Thereupon (ii) is satisfied, if and only if the numbers $a_{i}$ maximize this
sum. Since the covariance matrix $\mathrm{Cov}(\boldsymbol{X})$ is positive
semi-definite, the eigenvalues $\lambda_{i}$ are not negative such that the
sum is maximized for:

$\displaystyle a_{i}$ $\displaystyle=\begin{cases}1&\text{for
}i\in\\{1,\,\ldots\,k\\}\\\ 0&\text{else}\end{cases}$

Such that:

$\displaystyle\sum_{j=1}^{k}\boldsymbol{u}_{j}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{u}_{j}$
$\displaystyle=\sum_{i=1}^{d}\lambda_{i}a_{i}$
$\displaystyle=\sum_{i=1}^{k}\lambda_{i}$
$\displaystyle=\sum_{j=1}^{k}\boldsymbol{c}_{j}^{T}\mathrm{Cov}(\boldsymbol{X})\boldsymbol{c}_{j}$

Accordingly the choice $\boldsymbol{u}_{j}=\boldsymbol{c}_{j}$ for
$i\in\\{1,\,\ldots\,,\,k\\}$ maximizes $\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$
and $L$ has a representation, given by:

$L=\mathbb{E}(\boldsymbol{X})+\bigoplus_{i=1}^{k}\mathbb{R}\boldsymbol{c}_{i}$

”$\Longleftarrow$” Let $L$ have a representation as given by (ii), then (1)
$\mathbb{E}(\boldsymbol{X})\in L$ and (2) the variance
$\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$ is maximized. According to Proposition
3 it follows, that $L$ minimizes the Euclidean distance to $\boldsymbol{X}$. ∎

######  Definition ($L$-Correlation).

_For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be a
random vector and $L$ a maximal linear principal manifold for
$\boldsymbol{X}$. Then for any $i,\,j\in\\{1,\,\ldots\,,\,d\\}$ let the
$L$-Correlation between $X_{i}$ and $X_{j}$ be defined by:_

$\rho_{X_{i},X_{j}\mid L}^{2}\coloneqq R_{i}R_{j}$ (3.13)

_where with the orthogonal projection_ $\pi_{L}\colon\mathbb{R}^{d}\to L$_for
any $i\in\\{1,\,\ldots\,,\,d\\}$ the __reliability_ _of $X_{i}$ with respect
to $L$ is given by:_

$R_{i}\coloneqq
1-\frac{\mathrm{Var}_{i}\left(\boldsymbol{X}-\pi_{L}(\boldsymbol{X})\right)}{\mathrm{Var}_{i}(\boldsymbol{X})}$
(3.14)

######  Proposition 7.

For an elliptically distributed random vector the $L$-Correlation generalizes
the Pearson Correlation to maximal linear principal manifolds.

###### Proof.

Let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{2}$ be an elliptically
distributed random vector and $L$ a maximal linear $1$-principal manifold for
$\boldsymbol{X}$. Then for $i\in\\{1,\,2\\}$ the random error of the variable
$X_{i}$ has a variance:

$\eta_{X_{i}}^{2}=\mathrm{Var}_{X_{i}}\left(\boldsymbol{X}-\pi_{L}(\boldsymbol{X})\right)$

Such that by the definition of the reliability it follows, that:

$R_{i}\stackrel{{\scriptstyle{\scriptstyle\text{\ref{eq:L-Correlation
2}}}}}{{=}}1-\frac{\eta_{X_{i}}^{2}}{\sigma_{X_{i}}^{2}}$

Consequently:

$\rho_{X_{i},X_{j}\mid
L}^{2}=\left(1-\frac{\eta_{X_{i}}^{2}}{\sigma_{X_{i}}^{2}}\right)\left(1-\frac{\eta_{X_{j}}^{2}}{\sigma_{X_{j}}^{2}}\right)$

With $n\in\mathbb{N}$ i.i.d. realizations
$\boldsymbol{x}\in\mathbb{R}^{n\times 2}$ of $\boldsymbol{X}$ an empirical
$L$-Correlation $\rho_{x_{i},x_{j}\mid L}^{2}$ is then given by replacing the
variances by the sample variances. Then by proposition 2 it follows, that:

$\rho_{x_{i},x_{j}}^{2}\overset{{\scriptstyle P}}{\to}\rho_{x_{i},x_{j}\mid
L}^{2},\text{ for }n\to\infty$

∎

##  4 The Riemann-Pearson Correlation

Linear principal manifolds allow the projection of a random vector
$\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ onto a linear subspace
$L\subseteq\mathbb{R}^{d}$, which maximally preserves the linear dependency
structure of $\boldsymbol{X}$ in terms of its covariances. Thereby for the
orthogonal projection $\pi_{L}\colon\mathbb{R}^{d}\hookrightarrow L$, the
variance on $L$, given by $\mathrm{Var}(\pi_{L}(\boldsymbol{X}))$, is referred
as the explained variance and the orthogonal deviation
$\mathrm{Var}(\boldsymbol{X}-\pi_{L}(\boldsymbol{X}))$ as the unexplained
variance. Thereupon by the assumption, that $\boldsymbol{X}$ is elliptically
distributed, it can be concluded, that linear independence coincides with
statistically independence, that that $\pi_{L}(\boldsymbol{X})$ and
$\boldsymbol{X}-\pi_{L}(\boldsymbol{X})$ are statistically independent and
therefore allow the following decomposition:

$\underset{{\scriptstyle\text{total
variance}}}{\underbrace{\mathrm{Var}(\boldsymbol{X})}}=\underset{{\scriptstyle\text{explained
variance}}}{\underbrace{\mathrm{Var}(\pi_{L}(\boldsymbol{X}))}}+\underset{{\scriptstyle\text{unexplained
variance}}}{\underbrace{\mathrm{Var}(\boldsymbol{X}-\pi_{L}(\boldsymbol{X}))}}$

This decomposition, as shown by theorem 7, is of fundamental importance for
the correlation over linear Principal Manifolds, since it determines the
reliabilities of the respective random variable $X$ by the ratio:

$R=1-\frac{{\scriptstyle\text{explained variance}}}{{\scriptstyle\text{total
variance}}}$

$\scriptstyle{\boldsymbol{x}^{1}}$$\gamma$ Figure 4.1: Principal Curve for a
2-dimensional realization

On this point of the discussion it’s just a small step to generalize the
principal components, by a smooth curves
$\gamma\colon[a,\,b]\to\mathbb{R}^{d}$ (figure 4.1). This is particular
appropriate, if the assumption of an elliptically distribution can only hardly
be justified, like for observed dynamical systems. Thereby the evolution
function generates a smooth submanifold $\mathcal{M}\subseteq\mathbb{R}^{d}$
within the observation space $\mathbb{R}^{d}$, and an “error free” observation
can be identified by a random vector $\boldsymbol{X}^{\star}$, with outcomes
on $\mathcal{M}$. Additionally, however, the observation function may be
regarded to be subjected to a measurement error $\boldsymbol{\varepsilon}$. By
the assumption, that $\boldsymbol{\varepsilon}$ has an elliptical
distribution, then the distribution of the observable random vector
$\boldsymbol{X}$ is represented by a elliptical $\mathcal{M}-$distribution.

######  Definition ($\mathcal{M}$-Distribution).

_For $d\in\mathbb{N}$ let
$\boldsymbol{X}^{\star}\colon\Omega\to\mathbb{R}^{d}$ be a random vector and
$\mathcal{M}\subseteq\mathbb{R}^{d}$ a smooth $k-$submanifold of
$\mathbb{R}^{d}$ with $k\leq d$. Then $\boldsymbol{X}^{\star}$ is
$\mathcal{M}-$distributed, iff for the probability density $P$, which is
induced by $\boldsymbol{X}^{\star}$, it holds, that:_

$P(\boldsymbol{X}^{\star}=\boldsymbol{x})>0\Leftrightarrow\boldsymbol{x}\in\mathcal{M}$
(4.1)

_Thereupon a random vector $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ is
elliptically $\mathcal{M}$-distributed, iff there exists an
$\mathcal{M}$-distributed random vector
$\boldsymbol{X}^{\star}\colon\Omega\to\mathbb{R}^{d}$ and an elliptically
distributed random error
$\boldsymbol{\varepsilon}\colon\Omega\to\mathbb{R}^{d}$, such that:_

$\boldsymbol{X}\sim\boldsymbol{X}^{\star}+\boldsymbol{\varepsilon}$ (4.2)

$\mathcal{M}$$\scriptstyle{\sigma_{2}}$$\scriptstyle{\sigma_{3}}$$X_{1}$$X_{2}$$X_{3}$$\scriptstyle{\sigma_{1}}$
Figure 4.2: Elliptical $\mathcal{M}$-distribution in $3$ dimensions

The assumption, that the observed random vector $\boldsymbol{X}$, is
elliptically $\mathcal{M}$-distributed, is very general, but allows an
estimation of $\mathcal{M}$ by minimizing the averaged Euclidean distance to
$\boldsymbol{X}$. Thereby the tangent spaces $T_{x}\mathcal{M}$ have a basis,
given by $k$ principal components of local infinitesimal covariances, such
that the remaining $d-k$ principal components describe the normal space
$N_{x}\mathcal{M}$, which is orthogonal to the tangent space
$T_{x}\mathcal{M}$. Since $T_{x}\mathcal{M}$ and $N_{x}\mathcal{M}$ are
equipped with an induced Riemannian metric, which is simply given by the
standard scalar product, there exists a minimal orthogonal projection
$\pi_{\mathcal{M}}\colon\mathbb{R}^{d}\hookrightarrow\mathcal{M}$, that maps
any realization $\boldsymbol{x}$ of $\boldsymbol{X}$ to a closest point on
$\mathcal{M}$. Then proposition 3 motivates properties for $\mathcal{M}$ to
minimize the averaged Euclidean distance to realizations of $\boldsymbol{X}$.
This provides the definition of smooth $k$-principal manifolds (Hastie et al.
1989, p513).

######  Definition (Principal Manifold).

_For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be a
random vector, $\mathcal{M}\subseteq\mathbb{R}^{d}$ a (smooth) $k$-submanifold
of $\mathbb{R}^{d}$ with $k\leq d$ and
$\pi_{\mathcal{M}}\colon\mathbb{R}^{d}\hookrightarrow\mathcal{M}$ a minimal
orthogonal projection onto $\mathcal{M}$. Then $\mathcal{M}$ is a (smooth)
$k$-principal manifold for $\boldsymbol{X}$, iff
$\forall\boldsymbol{x}\in\mathcal{M}$ it holds, that:_

$\mathbb{E}(\boldsymbol{X}\in\pi_{\mathcal{M}}^{-1}(\boldsymbol{x}))=\boldsymbol{x}$
(4.3)

_Furthermore $\mathcal{M}$ is termed maximal, iff $\mathcal{M}$ maximizes the
explained variance $\mathrm{Var}(\pi_{\mathcal{M}}(\boldsymbol{X}))$._

By extending the local properties of the tangent spaces to the underlying
manifold, by propositions 3 and 6 it can be concluded, that maximal principal
manifolds minimize the Euclidean distance to $\boldsymbol{X}$. Intuitively
this can be understood as follows: The principal manifold property assures,
that:

$\mathrm{Var}(\boldsymbol{X})=\mathrm{Var}(\pi_{\mathcal{M}}(\boldsymbol{X}))+\mathrm{Var}(\boldsymbol{X}-\pi_{\mathcal{M}}(\boldsymbol{X}))$

Consequently the choice of $\mathcal{M}$ maximizes
$\mathrm{Var}(\pi_{\mathcal{M}}(\boldsymbol{X}))$ if and only if it minimizes
$\mathrm{Var}(\boldsymbol{X}-\pi_{\mathcal{M}}(\boldsymbol{X}))$, which equals
the variance of the error and therefore the Euclidean distance. At closer
inspection, however, it turns out, that in difference to linear principal
manifolds, the maximization problem is ill-defined for arbitrary smooth
principal manifolds, since for any finite number of realizations trivial
solutions can be found by smooth principal manifolds, that interpolate the
realizations and therefore provide a perfect explanation. In order to close
this gap, further structural structural assumptions have to be incorporated,
ether by a parametric family $\\{\boldsymbol{f}_{\theta}\\}_{\theta\in\Theta}$
that restricts the possible solutions - or by a regularization, as given in
the elastic map algorithm that penalizes long distances and strong curvature
(Gorban et al. 2008). Due to the complexity of this topic, however, it is left
to the second chapter, where Energy based models are used to overcome this
deficiency. In the following the generalization of the correlation to smooth
principal manifolds for convenience is defined with respect to a principal
manifold $\mathcal{M}$, which is maximal “with respect to appropriate
restrictions”.

######  Definition (Riemann-Pearson Correlation).

_For $d\in\mathbb{N}$ let $\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$ be a
random vector, $\mathcal{M}$ a smooth principal manifold for $\boldsymbol{X}$,
which is maximal “with respect to appropriate restrictions” and
$\pi_{\mathcal{M}}\colon\mathbb{R}^{d}\to\mathcal{M}$ a minimal orthogonal
projection. Then for any $i,\,j\in\\{1,\,\ldots\,,\,d\\}$ the Riemann-Pearson
Correlation between $X_{i}$ and $X_{j}$ is given by:_

$\rho_{X_{i},X_{j}\mid\mathcal{M}}^{2}\coloneqq
R_{i}R_{j}\int_{\mathcal{M}}S_{i,j}(\boldsymbol{x})S_{j,i}(\boldsymbol{x})\,\mathrm{d}P_{\mathcal{M}}$
(4.4)

_where for $i\in\\{1,\,\ldots\,,\,d\\}$ the __reliability_ _of $X_{i}$ with
respect to $\mathcal{M}$ is given by:_

$R_{i}\coloneqq
1-\frac{\mathrm{Var}_{i}\left(\boldsymbol{X}-\pi_{\mathcal{M}}(\boldsymbol{X})\right)}{\mathrm{Var}_{i}(\boldsymbol{X})}$
(4.5)

_and for $i,\,j\in\\{1,\,\ldots\,,\,d\\}$ the __local sensitivity_ _of $X_{i}$
with respect to $X_{j}$ by:_

$S_{i,j}(\boldsymbol{x})\coloneqq\left.\frac{\partial}{\partial
x_{j}}(\boldsymbol{x}-\pi_{\mathcal{M}}(\boldsymbol{x}))\right|_{i}$ (4.6)

######  Proposition 8.

For an elliptical $\mathcal{M}-$distributed random vector $\boldsymbol{X}$ the
Riemann-Pearson Correlation generalizes the $L$-Correlation to smooth
principal manifolds.

###### Proof.

Let $L$ be a maximal linear principal manifold for
$\boldsymbol{X}\colon\Omega\to\mathbb{R}^{d}$, and $\boldsymbol{x}$ a
realization of $\boldsymbol{X}$ then there exists an $\beta\in\mathbb{R}$
with:

$S_{i,j}(\boldsymbol{x})=\left.\frac{\partial}{\partial
x_{j}}(\boldsymbol{x}-\pi_{L}(\boldsymbol{x}))\right|_{i}\equiv\beta$

Furthermore for $c\neq 0$:

$S_{j,i}(\boldsymbol{x})=\left.\frac{\partial}{\partial
x_{i}}(\boldsymbol{x}-\pi_{L}(\boldsymbol{x}))\right|_{j}\equiv\frac{1}{\beta}$

Such that $S_{i,j}(\boldsymbol{x})S_{j,i}(\boldsymbol{x})=1$. Consequently for
$\mathcal{M}=L$ it follows, that:

$\displaystyle\rho_{X_{i},X_{j}\mid\mathcal{M}}^{2}$
$\displaystyle=R_{i}R_{j}\int_{\mathcal{M}}S_{i,j}(\boldsymbol{x})S_{j,i}(\boldsymbol{x})\,\mathrm{d}P_{\mathcal{M}}$
$\displaystyle=R_{i}R_{j}\int_{\mathcal{M}}\,\mathrm{d}P_{\mathcal{M}}$
$\displaystyle=R_{i}R_{j}=\rho_{X_{i},X_{j}\mid L}^{2}$

∎

## References

  * [1] Shurong Zheng, Ning-zhong Shi, and Zhengjun Zhang.  Generalized Measures of Correlation.  Manuscript, pages 1–45, 2010.
  * [2] Karl Pearson.  On lines and planes of closest fit to systems of points in space.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(1):559–572, 1901.
  * [3] J. Kenney and E. Keeping.  Linear regression and correlation.  Mathematics of statistics, 1:252–285, 1962.
  * [4] G. W. Snedecor and W. G. Cochran.  Statistical Methods.  Iowa State University Press, Ames, 6 edition, 1967.
  * [5] Trevor Hastie and Werner Stuetzle.  Principal Curves.  Journal of the American Statistical Association, 84(406):502–516, 1989.
  * [6] Alexander N. Gorban and Andrei Y. Zinovyev.  Principal Graphs and Manifolds.  Arxiv preprint arXiv08090490, page 36, 2008.


#  Screwon spectral statistics and dispersion relation in the quantum Rajeev-
Ranken model

Govind S. Krishnaswamia and T. R. Vishnua,b
aPhysics Department, Chennai Mathematical Institute, SIPCOT IT Park, Siruseri
603103, India
bPhysics Department, Raman Research Institute, Sadashivanagar, Bengaluru
560080, India

Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
(December 20, 2023)

###### Abstract

The Rajeev-Ranken (RR) model is a Hamiltonian system describing screw-type
nonlinear waves (screwons) of wavenumber $k$ in a scalar field theory
pseudodual to the 1+1D SU(2) principal chiral model. Classically, the RR model
based on a quadratic Hamiltonian on a nilpotent/Euclidean Poisson algebra is
Liouville integrable. Upon adopting canonical variables in a slightly extended
phase space, the model was interpreted as a novel 3D cylindrically symmetric
quartic oscillator with a rotational energy. Here, we examine the spectral
statistics and dispersion relation of quantized screwons via numerical
diagonalization validated by variational and perturbative approximations. We
also derive a semiclassical estimate for the cumulative level distribution
which compares favorably with the one from numerical diagonalization. The
spectrum shows level crossings typical of an integrable system. The $i^{\rm
th}$ unfolded nearest neighbor spacings are found to follow Poisson statistics
for small $i$. Nonoverlapping spacing ratios also indicate that successive
spectral gaps are independently distributed. After displaying universal linear
behavior over energy windows of short lengths, the spectral rigidity saturates
at a length and value that scales with the square-root of energy. For strong
coupling $\lambda$ and intermediate $k$, we argue that reduced screwon
energies can depend only on the product $\lambda k$. Numerically, we find
power law dependences on $\lambda$ and $k$ with an approximately common
exponent $2/3$ provided the angular momentum quantum number $l$ is small
compared to the number of nodes $n$ in the radial wavefunction. On the other
hand, for the ground state $n=l=0$, the common exponent becomes 1.

###### Contents

  1. 1 Introduction
  2. 2 Variational approximation to screwon ground state
  3. 3 First order perturbation theory
  4. 4 Screwon spectrum via finite-differences
  5. 5 Energy level statistics
    1. 5.1 Semiclassical cumulative energy level distribution
    2. 5.2 Level spacing distributions
    3. 5.3 Number variance
    4. 5.4 Spectral rigidity
  6. 6 Dispersion relation at strong coupling
  7. 7 Discussion
  8. A Evaluating the spectral rigidity

##  1 Introduction

The Rajeev-Ranken (RR) model [2] is a mechanical system with three degrees of
freedom. It describes nonlinear screw-type continuous waves in a
$1+1$-dimensional nilpotent scalar field theory which is dual to the
SU(2)-principal chiral model [3]. The scalar field $\phi$ is valued in the
$\mathfrak{su}$(2) Lie algebra and classically satisfies the quadratically
nonlinear evolution equations
$\ddot{\phi}=\phi^{\prime\prime}+\lambda[\dot{\phi},\phi^{\prime}]$ for a
dimensionless coupling constant $\lambda$. The screw-type waves are solutions
of the form

$\phi(x,t)=e^{Kx}R(t)e^{-Kx}+mKx,$ (1)

where $R(t)$ is a dynamical $\mathfrak{su}$(2) matrix with the three real
degrees of freedom $R_{a}=i\,{\rm tr}\,(R\sigma_{a})$ for $a=1,2,3$. Moreover,
$K=ik\sigma_{3}/2$ with $k$ a constant with dimensions of a wavenumber and $m$
is a dimensionless real parameter. The resulting equations of motion for $R$
are elegantly expressed as

$\dot{L}=[K,S]\quad\text{and}\quad\dot{S}=\lambda[S,L]\quad\text{where}\quad
L=[K,R]+mK\quad\text{and}\quad S=\dot{R}+\frac{K}{\lambda}.$ (2)

These equations follow from a quadratic Hamiltonian and a nilpotent Poisson
algebra among $L$ and $S$ in a six dimensional phase space with coordinates
$L_{a}=i\,{\rm tr}\,(L\sigma_{a})$ and $S_{a}=i\,{\rm tr}\,(S\sigma_{a})$.
They also admit a Lax pair formulation and classical $r$-matrix leading to a
complete set of independent conserved quantities in involution implying
Liouville integrability [4, 5]. However, $L_{3}=-mk$ is a Casimir invariant of
the nilpotent algebra and $R_{3}$ does not appear in the $L-S$ equations of
motion. The evolution of $R_{3}$ is determined by
$\dot{R_{3}}=S_{3}+k/\lambda$. Thus in this paper, we will treat $L_{3}$ as a
parameter and introduce $R_{3}$ as a coordinate. The equations of motion of
the RR model admit an alternate Hamitonian-Poisson bracket formulation in
terms of canonical variables $R_{1,2,3}$ and their conjugate momenta [4, 5]:

$kP_{1,2}=\dot{R}_{1,2}\pm\frac{1}{2}\lambda kmR_{2,1}\quad\text{and}\quad
kP_{3}=\dot{R}_{3}+\frac{1}{2}\lambda k(R_{1}^{2}+R_{2}^{2}).$ (3)

The Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\left(kP_{1}-\frac{\lambda
kmR_{2}}{2}\right)^{2}+\left(kP_{2}+\frac{\lambda
kmR_{1}}{2}\right)^{2}+\left(kP_{3}-\frac{\lambda
k}{2}(R_{1}^{2}+R_{2}^{2})\right)^{2}\right]$ (5)
$\displaystyle+\frac{k^{2}}{2}(R_{1}^{2}+R_{2}^{2}+m^{2}),$

may be interpreted as that of a cylindrically symmetric anharmonic oscillator
in three dimensions [6]. Indeed, if we regard $R_{1,2,3}$ and $kP_{1,2,3}$ as
the Cartesian components of the position and momentum of a particle of mass
$\mu=1$, then (5) is equivalent to the Hamiltonian

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{p_{x}^{2}+p_{y}^{2}+p_{z}^{2}}{2\mu}+\frac{\lambda
km(xp_{y}-yp_{x})}{2\mu}+\left(\frac{\lambda^{2}k^{2}m^{2}}{8\mu}-\frac{\lambda
kp_{z}}{2\mu}+\frac{k^{2}}{2}\right)(x^{2}+y^{2})$ (7)
$\displaystyle+\frac{\lambda^{2}k^{2}}{8\mu}(x^{2}+y^{2})^{2}+\frac{k^{2}m^{2}}{2}.$

This Hamiltonian describes a particle of mass $\mu$ moving in a cylindrically
symmetric quadratic plus quartic potential with an additional rotational
energy proportional to $L_{z}=xp_{y}-yp_{x}$.

Upon canonical quantization, the cylindrical symmetry of $H$ permits
separation of variables in the Schrödinger equation for the wavefuntion

$\psi(r,\theta,z)=\frac{1}{\sqrt{r}}\varrho(r)\exp{(il\theta)}\exp{\left(\frac{ip_{z}z}{\hbar}\right)}\quad\text{for}\quad
r\geq 0.$ (8)

Here, $l\hbar$ is the eigenvalue of $L_{z}$ and can be any integer while
$p_{z}\in\mathbb{R}$ is the eigenvalue of the operator $p_{z}$, both of which
commute with the Hamiltonian. Division by $\sqrt{r}$ eliminates the
$\varrho^{\prime}$ term in the resulting radial eigenvalue equation

$-\frac{\hbar^{2}}{2\mu}\varrho^{\prime\prime}(r)+\left(U(r)+\frac{\hbar^{2}}{2\mu
r^{2}}\left(l^{2}-\frac{1}{4}\right)+\frac{\hbar l\lambda
km}{2\mu}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}\right)\varrho=E\varrho.$
(9)

Here,

$U(r)=\alpha r^{2}+\beta
r^{4}\quad\text{where}\quad\alpha=\frac{\lambda^{2}k^{2}m^{2}}{8\mu}-\frac{\lambda
kp_{z}}{2\mu}+\frac{k^{2}}{2}\quad\text{and}\quad\beta=\frac{\lambda^{2}k^{2}}{8\mu}.$
(10)

The parameters $m$ and $k^{2}$ have dimensions of length and energy/area,
while $\lambda^{2}$ has dimensions of mass/area. Although $k$ does not have
dimensions of inverse length, we will refer to it as a wavenumber since it
arose as a wavenumber in the screwon ansatz (1). In what follows, we will
often work in units of mass, length and angular momentum where $\mu=1$, $m=1$
and $\hbar=1$. The Hamiltonian in Eqn.(7) is the sum of squares of Hermitian
operators. Thus, the energy eigenvalues must satisfy $E\geq 0$. These are the
energies of the quantized screw-type waves which we refer to as screwons. Due
to translational invariance in the $z$-direction, $E$ includes a free particle
contribution $p_{z}^{2}/2\mu$. We will mostly be interested in the
contribution to the energies coming from the dynamics in the $x$-$y$ plane
which for fixed $p_{z}$ is discrete and labelled by $l$ and a radial quantum
number $n$. Aside from $p_{z}$ these energy levels $E_{n,l}$ depend
parametrically on the coupling $\lambda$ and the screwon wavenumber $k$.
Unlike quantized small oscillations around the vacuum of the scalar field
theory, screwons are new nonlinear degrees of freedom that could play a role
similar to solitons in other field theories. Thus it is interesting to
understand the spectrum of screwons, their shapes, dispersion relations
(dependence on $k$) and their behavior at strong coupling $\lambda$, which
corresponds to high energies in the field theory.

In this article, we study the spectral statistics and dispersion relations of
quantized screwons. In Section 2, we use the (radial) asymptotic behaviour
known from [6] to propose an ansatz for the screwon ground state wavefunction
and obtain a variational estimate for the ground state energy. In Section 3,
we use first order perturbation theory to determine the dependences of
$E_{0,l}$ on $l$ and $\lambda$ at weak coupling. In Section 4, we find the
screwon spectrum by numerical diagonalization of the radial Hamiltonian for
each fixed $l$. While we do this for a range of values of $k$ and $\lambda$,
we only consider the case $p_{z}=1$. In fact, for large coupling and an
appropriate range of moderate $k$, we argue that the shifted spectrum
($E-p_{z}^{2}/2\mu-m^{2}k^{2}/2$) is independent $p_{z}$. Combining the
contributing $l$ sectors we determine the lowest $2\times 10^{4}$ energy
levels. To do this accurately, we optimize the spatial grid (for each value of
$k$ and $\lambda$) to capture the tail as well as all oscillations of the
radial wavefunction. The numerical results are validated using our variational
estimate and perturbative estimate at weak coupling. As expected of a
classically integrable system, the screwon energies display level crossings as
$\lambda$ is varied.

Section 5 is devoted to the level statistics of the screwon spectrum. We
express the cumulative level distribution $n(E)$ as a sum of contributions
$n(E;p_{\theta})$ from all allowed angular momentum sectors. Semiclassical
estimates for $n(E;p_{\theta})$ and $n(E)$ match those obtained from numerical
diagonalization for all but the lowest energies. Asymptotically, we find that
$n(E)$ displays power law growth while
$n(E;p_{\theta})\sim(E-E_{0,l})^{\zeta(|p_{\theta}|)}$, with $\zeta$
increasing monotonically with $|p_{\theta}|$. Here $E_{0,l}$ is the ground
state energy for $p_{\theta}=l\hbar$. We then examine spacing distributions of
the screwon spectrum. The $i^{\rm th}$ unfolded nearest neighbour spacing
distributions follow Poisson statistics although deviations become more
pronounced as $i$ increases. The independence in the distribution of spectral
gaps is also confirmed using the statistics of nonoverlapping spacing ratios.
Both these are as expected of a classically integrable system. Number variance
$\Sigma(E,L)$ and spectral rigidity $\Delta(E,L)$, which encode fluctuations
in the spectrum around a central energy $E$, display universal linear
behaviour over energy windows of short lengths $L$. For larger $L$, they
saturate and display oscillations that should reveal system specific details.
Both the value of $\Delta$ and $L$ at which the saturation occurs show square-
root power law behavior as a function of central energy $E$.

In Section 6, we turn to the screwon dispersion relation, by which we mean the
dependence of energy levels on $\lambda$ and $k$. We focus attention on a
strong coupling regime where the (shifted) energy levels depend on $\lambda$
and $k$ only through their product. Further restricting $k$ to an intermediate
range of wave numbers, ensures that the shifted energies are independent of
$p_{z}$. In this phase, our numerical results show that the shifted energies
depend on $\lambda$ and $k$ through a common power $\eta(l)$. For $l=0$, this
exponent $\eta(0)=1$ for the ground state $n=0$; it rapidly drops to
$\eta(0)=2/3$ with increasing radial quantum number $n\gtrsim 100$. This $2/3$
power law for radially excited screwons continues to hold as long as $|l|\ll
n$ although there are deviations when $|l|$ is comparable to or exceeds $n$.
These results confirm and extend our earlier results on the screwon dispersion
relation obtained via a WKB approximation in [6]. This fractional power law
dispersion relation encodes information on the behavior of quantized screwons
at strong coupling and gives us a glimpse into the possible behavior of the
high energy degrees of freedom of the scalar field theory. We conclude with a
discussion in Section 7, highlighting some interesting outstanding problems
that have arisen from this work.

##  2 Variational approximation to screwon ground state

We may obtain an approximation to the screwon ground state energy using the
Rayleigh-Ritz variational principle. From §4.3 of [6], away from the weak
coupling limit, the asymptotic behaviour of the radial wavefunction (8) is:

$\varrho(r)\sim
r^{-1}\exp{\left(-\frac{\sqrt{2\mu\beta}}{\hbar}\left(\frac{r^{3}}{3}+\frac{\alpha
r}{2\beta}\right)\right)}\quad\text{as}\quad
r\to\infty\quad\text{and}\quad\varrho(r)\to
r^{|l|+\frac{1}{2}}\quad\text{as}\quad r\to 0.$ (11)

An ansatz that incorporates this asymptotic behavior and has no nodes (as
appropriate to the lowest-lying state $n=0$ for any fixed $l$) is

$\varrho_{\rm try}(r)\propto\frac{r^{|l|+\frac{1}{2}}}{1+\zeta
r^{|l|+\frac{3}{2}}}\exp{\left(-\frac{\sqrt{2\mu\beta}}{\hbar}\left(\frac{r^{3}}{3}+\frac{\alpha
r}{2\beta}\right)\right)},$ (12)

where $\zeta\geq 0$ is a variational parameter. The peak in this wavefunction
drops in height and shifts towards $r=0$ as $\zeta$ increases. Note that this
trial wavefunction is not valid in the weak coupling limit since when
$\lambda\to 0$, $\beta\to 0$ and

$\varrho_{\lambda=0}(r)\propto
r^{|l|+1/2}\exp\left(-\frac{\sqrt{\mu}kr^{2}}{2\hbar}\right)L_{p}^{|l|}\left(-\frac{\sqrt{\mu}kr^{2}}{\hbar}\right).$
(13)

Here, $L_{p}^{|l|}$ are generalized Laguerre polynomials [7] and
$p=(n-|l|)/2$. Extremizing the expectation value of the radial Hamiltonian (9)
in the state (12):

$E^{\rm
var}=\min_{\zeta}\frac{\langle\varrho|H|\varrho\rangle}{\langle\varrho|\varrho\rangle}\quad\text{where}\quad\langle\varrho_{1}|\varrho_{2}\rangle=\int
dr\varrho_{1}^{*}(r)\varrho_{2}(r),$ (14)

we find a variational upperbound for the screwon ground state energy. For
instance, when $\lambda=k=l=\mu=m=p_{z}=\hbar=1$, we find the optimal value
$\zeta=0.12$, corresponding to the ground state energy estimate $E^{\rm
var}_{n=0,l=1}=3.416$ which compares favorably with the numerical result
$E_{0,1}=3.415$ from §4. Similarly, for $l=0$, $E^{\rm var}_{n=0,l=0}=1.866$,
which is just above the numerical result $E_{0,0}=1.862$. Moreover, for fixed
$k$, the variational ground state energy increases roughly linearly with
$\lambda$ while the optimal value of $\zeta$ decreases. To examine spectral
statistics of screwons and their dispersion relation, we need to treat excited
states and reliably extract the dependence of energy eigenvalues on
parameters. This requires us to go beyond this variational approximation,
which we do via numerical diagonalization of $H$ in §4.

Figure 1: Comparison of variational (dashed) and finite-difference numerical
(solid) radial wavefunctions $\rho(r)=\sqrt{r}\varrho(r)$ for $n=0$ and
$l=0,1$ with $\lambda=k=1$ and $\hbar=m=\mu=1$. They go to zero as $r$ and
$r^{2}$ as $r\to 0$ and as $e^{-r^{3}/6}$ as $r\to\infty$.

##  3 First order perturbation theory

For weak coupling $\lambda$, we may use perturbation theory to estimate energy
eigenvalues. For fixed $p_{z}$, we write the 2D Hamiltonian (7) as

$\displaystyle H_{2D}$ $\displaystyle=$
$\displaystyle\frac{1}{2\mu}\left(p_{r}^{2}+\frac{p_{\theta}^{2}-\frac{\hbar^{2}}{4}}{r^{2}}\right)+\frac{k^{2}r^{2}}{2}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}+\lambda\left(\frac{kmp_{\theta}}{2\mu}-\frac{kp_{z}r^{2}}{2\mu}\right)$
(16)
$\displaystyle+\lambda^{2}\left(\frac{k^{2}m^{2}r^{2}}{8\mu}+\frac{k^{2}r^{4}}{8\mu}\right)=H_{0}+\lambda
H_{1}+\lambda^{2}H_{2}.$

Note that $H_{0}$ depends on $p_{\theta}=L_{z}$ only through its square. The
eigenvalues of the unperturbed Hamiltonian $H_{0}$ are

$E_{n,l}^{(0)}=(2n+|l|+1)\frac{\hbar|k|}{\sqrt{\mu}}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}.$
(17)

Here $n=0,1,2,\cdots$ and $l$ is an integer. For $n=0$, the corresponding
normalized eigenvectors of $H_{0}$ are

$\psi_{0,l}^{(0)}(r,\theta)=\left[\left(\frac{\sqrt{\mu}k}{\hbar}\right)^{|l|+1}\frac{1}{\pi|l|!}\right]^{\frac{1}{2}}r^{|l|}\exp\left(-\frac{\sqrt{\mu}k}{\hbar}\frac{r^{2}}{2}\right)\exp(il\theta).$
(18)

Thus, to first order in perturbation theory, $E_{0,l}=E_{0,l}^{(0)}+\lambda
E_{0,l}^{(1)}+\cdots$ where

$E_{0,l}^{(1)}=\langle\psi_{0,l}^{(0)}|H_{1}|\psi_{0,l}^{(0)}\rangle=\frac{kml\hbar}{2\mu}-\frac{kp_{z}}{2\mu}(|l|+1)\left(\frac{\hbar}{\sqrt{\mu}k}\right).$
(19)

For instance, if $\hbar=\mu=k=m=p_{z}=l=1$, we find the ground state energy

$E_{0,0}=E_{0,0}^{(0)}+\lambda E_{0,0}^{(1)}+{\cal
O}(\lambda^{2})=2-\frac{\lambda}{2}+{\cal O}(\lambda^{2}).$ (20)

For $\lambda=0.01,0.05,0.08$, the ground state energies $E_{0,0}$ with first
order correction are $1.995,1.975$ and $1.96$ while those obtained numerically
via finite differences are $1.99,1.952$ and $1.925$, showing the loss in
accuracy of first order perturbation theory with increasing $\lambda$.
Nevertheless, for small $\lambda$, perturbation theory predicts a drop in the
energy with increasing $\lambda$. This is visible in the numerical results
shown in Fig. 2.

##  4 Screwon spectrum via finite-differences

Here, we address the question of determining the spectrum of screwons by
numerical diagonalization of the Hamiltonian (7) which (aside from free
particle motion along $z$) describes a cylindrically symmetric anharmonic
oscillator in the $x$-$y$ plane. A direct approach would be to use finite
differences to discretize the 2D Hamiltonian, truncate it and diagonalize the
resulting matrix. An alternative is to numerically determine the spectrum of
the radial Hamiltonian (9) in each angular momentum sector ($l=0,1,2,\ldots$)
and then combine these sectors to obtain the full screwon spectrum. We
implemented both approaches and found the second to be more appropriate for
the following reasons. (a) For comparable precision, diagonalizing the 2D
Hamiltonian was computationally slower. (b) A 2D square lattice does not
possess the circular symmetry of $H$ making it harder to accurately capture
the behavior of the wave function as $r\to 0$ and for large $r$. (c) We wish
to classify energy levels using the quantum numbers $n$ (number of nodes of
$\varrho(r)$) and $l$ to facilitate the search for level crossings. (d) We
need the spectrum of the radial Hamiltonian for fixed $l$ in order to look for
a power law dispersion relation $E_{nl}\sim(\lambda k)^{\eta}$.

To find the spectrum of (9), we discretize the radial coordinate $r=i\delta$
for $i=1,2,\ldots,N_{r}$, with a small spacing $\delta$. This converts (9)
into an eigenvalue problem for an $N_{r}\times N_{r}$ matrix. However, direct
implementation using the wavefunction defined in (8)
$\psi\propto\varrho/\sqrt{r}$ leads to eigenvalues that converge very slowly
(with increasing $N_{r}$) to the known spectrum (17) for $\lambda=0$ when
$l=0$. We traced the numerical difficulty to the behavior of $\varrho(r)$ as
$r\to 0$. From [6] it is known that

$\varrho(r)\to r^{|l|+\frac{1}{2}}\quad\text{as}\quad r\to 0.$ (21)

We then found that the numerical scheme converges much faster if we use (in
place of (8)), the radial wavefunction

$\psi_{a}(r,\theta,z)=\frac{1}{r^{a}}\rho(r)\exp{(il\theta)}\exp{\left(\frac{ip_{z}z}{\hbar}\right)}\quad\text{with}\quad
a=1.$ (22)

This ensures a softer behaviour of the radial wavefunction near the origin:
$\rho(r)\to r^{a+|l|}$ as $r\to 0$. The radial equation corresponding to (22)
is

$-\frac{\hbar^{2}}{2\mu}\left(\rho^{\prime\prime}(r)+\frac{1-2a}{r}\rho^{\prime}+\frac{a^{2}-l^{2}}{r^{2}}\rho\right)+\left(\alpha
r^{2}+\beta r^{4}+\frac{\hbar l\lambda
km}{2\mu}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}\right)\rho=E\rho.$ (23)

As mentioned, we discretize $r$ with equal spacing $r_{i}=i\delta$ for
$i=1,2,\ldots,N_{r}$. To ensure that the discretized Hamiltonian is real
symmetric, we use the centered-difference formulae:

$\rho_{i}^{\prime\prime}=\frac{\rho_{i+1}-2\rho_{i}+\rho_{i-1}}{\delta^{2}}\quad\text{and}\quad\rho_{i}^{\prime}=\frac{\rho_{i+1}-\rho_{i-1}}{2\delta}.$
(24)

The radial equation becomes the system of difference equations:

$\displaystyle-\frac{\hbar^{2}}{2\mu}\left(\frac{\rho_{i+1}-2\rho_{i}+\rho_{i-1}}{\delta^{2}}+\frac{(1-2a)(\rho_{i+1}-\rho_{i-1})}{2i\delta^{2}}\right)$
(25) $\displaystyle+\left(\frac{\hbar^{2}}{2\mu
i^{2}\delta^{2}}\left(l^{2}-a^{2}\right)+\alpha i^{2}\delta^{2}+\beta
i^{4}\delta^{4}+\frac{\hbar l\lambda
km}{2\mu}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}\right)\rho_{i}=E\rho_{i}.$
(26)

We read off the discretized radial Hamiltonian, which is an $N_{r}\times
N_{r}$ tri-diagonal matrix with nontrivial entries

$\displaystyle H^{r}_{ii}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{\mu\delta^{2}}+\alpha i^{2}\delta^{2}+\beta
i^{4}\delta^{4}+\frac{\hbar^{2}}{2\mu
i^{2}\delta^{2}}\left(l^{2}-a^{2}\right)+\frac{\hbar l\lambda
km}{2\mu}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}\quad\text{and}$ (27)
$\displaystyle H^{r}_{i,i\pm 1}$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{2\mu\delta^{2}}\pm\frac{(1-2a)}{2i\delta^{2}}.\quad$
(28)

The ‘boundary’ matrix elements are $H_{1,0}=H_{N_{r},N_{r}+1}=0$. It may also
be written as

$\displaystyle H^{\rm r}_{ij}$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2\mu}\left(-\frac{2}{\delta^{2}}+\frac{a^{2}-l^{2}}{i^{2}\delta^{2}}\right)+\alpha
i^{2}\delta^{2}+\beta i^{4}\delta^{4}+\frac{\hbar l\lambda
km}{2\mu}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}\right]\delta_{ij}$ (30)
$\displaystyle-\frac{\hbar^{2}}{2\mu}\left(\frac{1}{\delta^{2}}+\frac{1-2a}{2i\delta^{2}}\right)\delta_{i,j-1}-\frac{\hbar^{2}}{2\mu}\left(\frac{1}{\delta^{2}}-\frac{1-2a}{2i\delta^{2}}\right)\delta_{i,j+1}\quad\text{for}\quad
1\leq i,j\leq N_{r}.$

Our aim is two-fold: (i) to get the first $N_{l}$ energy levels $E_{n,l}$ for
each of several values of $l$ with $|l|\leq l_{\rm max}$ and (ii) to combine
these spectra to get the first $N$ energy levels accounting for all values of
$l$. The latter is facilitated by the fact that the lowest energy for a given
$l$ increases with $|l|$.

Both (i) and (ii) are to be done for various values of $\lambda$, $k$ and
$p_{z}$ for fixed $\mu=m=\hbar=1$, by a choice of units. To find the lowest
$N_{l}$ levels of the radial Hamiltonian $H^{r}_{l}$, we need to pick a
spacing $\delta$ for the radial grid $r_{i}=i\delta$ for $i=1,\cdots,N_{r}$ to
construct the $N_{r}\times N_{r}$ matrix approximation to $H^{r}$ (30). For
fixed $N_{r}$, $\delta$ must be chosen so that the first $N_{l}$ radial
eigenfunctions $\psi_{n,l}(r)$ are accurately captured. In particular, $r_{\rm
max}=N_{r}\delta$ must be large enough to accommodate all the oscillations in
the most highly excited ($N_{l}^{\rm th}$) state desired, which must decay as
$e^{-\text{const}\times r^{3}}$ as $r\to\infty$ (11). On the other hand,
$\delta$ must be small enough to resolve each of the $N_{l}$ oscillations of
this wavefunction. From (11), the optimal value of $\delta$ is a decreasing
function of both $k$ and $\lambda$. Moreover, due to truncation errors, we may
trust only a small fraction $N_{l}/N_{r}\ll 1$ of the $N_{r}$ energy levels
computed by diagonalizing (30). In practice, for $\lambda=k=p_{z}=\mu=m=1$ we
find that we are able to reliably obtain the lowest $N_{l}=100$ levels by
choosing $N_{r}=5000$ and $r_{\rm max}=25$. In this case, we find that the
lowest $N_{l}$ energy eigenvalues change by less that $0.1\%$ if the size of
the matrix $N_{r}$ is increased from 4000 to 5000. Moreover, for $N_{l}=100$
we find that the highest energy for $l=0$ is $E_{\rm max}=E_{N_{l},0}=803.507$
while the ground state energy in all angular momentum sectors with $|l|>200$
lies above this value. Thus, by merging the spectra for $E\leq E_{\rm max}$
from all sectors with $|l|<l_{\rm max}=200$, we obtain the first $N=21061$
screwon energy levels. We have verified for a few values of $l$ (and
$\lambda=1$) that the ground state energies obtained numerically lie
marginally below the variational upper bounds obtained in §2.

##  5 Energy level statistics

Energy levels $E_{n,l}$ of the radial Hamiltonian $(\ref{e:radial-equation-
dimensionful-varrho})$ may be labelled by a principal quantum number
$n=0,1,2,\ldots$ (the number of nodes of the radial wavefunction $\varrho(r)$
in $(0,\infty)$) and by the angular momentum quantum number $l=0,\pm 1,\pm
2,\ldots$. For $\lambda=0$, the spectrum is given by [6]

$E_{n,l}=(2n+|l|+1)\frac{\hbar|k|}{\sqrt{\mu}}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}.$
(31)

Working in units where $\hbar=m=\mu=1$ and taking for definiteness
$k=p_{z}=1$, gives $E_{n,l}=2n+|l|+2$. The first few levels are as tabulated
below:

$E$ | 2 | 3 | 3 | 4 | 4 | 4 | 5 | 5 | 5 | 5
---|---|---|---|---|---|---|---|---|---|---
$n$ | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 1 | 1
$l$ | 0 | 1 | -1 | 2 | -2 | 0 | 3 | -3 | 1 | -1

For these parameter values, the level with energy $E$ is $(E-1)$-fold
degenerate. More generally, the level with energy $E_{n,l}$ is
$(2n+|l|+1)$-fold degenerate.

For $\lambda>0$, we find that these degeneracies are generally lost.
Heuristically, the rotational energy $\hbar l\lambda km/2\mu$ in (9) suggests
that within a $\lambda=0$ degenerate multiplet, levels with higher $l$ acquire
larger energies. This is illustrated in Fig. 2, which also shows that there
are no avoided level crossings. The resulting accidental degeneracies and the
lack of level repulsion are typical features of integrable systems.

(a)

(b)

Figure 2: Lowest lying screwon energy levels $E_{nl}$ vs coupling $\lambda$
for $k=p_{z}=\mu=m=1$. Degeneracies present when $\lambda=0$ break for
$\lambda>0$ though accidental degeneracies occur at level crossings. The
initial drop in energy with increasing $\lambda$ is explained by 1st order
perturbation theory in §3.

###  5.1 Semiclassical cumulative energy level distribution

Here, we obtain a semiclassical estimate for the cumulative energy level
distribution $n(E,\hbar)$ and compare it with that from numerical
diagonalization. As in §3, we fix $p_{z}$ and restrict attention to the
$x$-$y$ plane so that the energy spectrum $E_{i}$ is discrete. Recall that
$n(E,\hbar)$ is the number of levels with energy less than or equal to $E$:

$n(E,\hbar)=\sum_{i=1}^{N}\theta(E-E_{i}).$ (32)

Here, $\theta(x)$ is the unit step function defined as: $\theta(x)=1$ for
$x\geq 0$ and $\theta(x)=0$ for $x<0$. If we suppose that each state occupies
a volume $h^{2}$ in the 4D phase space, then

$n(E,\hbar)\sim\frac{1}{(2\pi\hbar)^{2}}\int_{H\leq
E}dx\,dy\,dp_{x}\,dp_{y}=\frac{1}{(2\pi\hbar)^{2}}\int_{H\leq
E}dr\,d\theta\,dp_{r}\,dp_{\theta}.$ (33)

The transformation to polar coordinates $r=\sqrt{x^{2}+y^{2}}$,
$\theta=\arccos(x/r)$, $p_{r}=(xp_{x}+yp_{y})/r$ and
$p_{\theta}=xp_{y}-yp_{x}$ is canonical, so the volume element is preserved.
Equation (33) is expected to hold in the semiclassical regime of small $\hbar$
or large $E$. [For the classical $\lambda=0$ model, the combinations of system
parameters (i.e., not including $p_{z}$ and $p_{\theta}$) with dimensions of
action and energy are $km^{2}\sqrt{\mu}$ and $k^{2}m^{2}$, so semiclassical
would mean $\hbar\ll km^{2}\sqrt{\mu}$ or $E\gg k^{2}m^{2}$. For $\lambda\neq
0$, there are other such combinations due to the presence of the dimensionless
parameter $\sqrt{\mu}/m\lambda$]. The classical 2D Hamiltonian (7) (for fixed
$p_{z}$)

$H=\frac{p_{r}^{2}}{2\mu}+\frac{p_{\theta}^{2}}{2\mu r^{2}}+\alpha r^{2}+\beta
r^{4}+\frac{\lambda
kmp_{\theta}}{2\mu}+\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}$ (34)

admits $\theta$ as a cyclic coordinate so that $p_{\theta}$ is conserved. The
integral over $\theta$ gives $2\pi$. We may choose to do the $p_{\theta}$
integral last:

$n(E,\hbar)\sim\frac{1}{(2\pi\hbar)^{2}}\int_{p_{\theta}^{\rm
min}(E)}^{p_{\theta}^{\rm
max}(E)}dp_{\theta}\int_{0}^{2\pi}d\theta\int_{H_{p_{\theta}}(r,p_{r})\leq
E}dr\,dp_{r}=\int_{p_{\theta}^{\rm min}(E)}^{p_{\theta}^{\rm
max}(E)}n(E,\hbar;p_{\theta})\,dp_{\theta}.$ (35)

Here $H_{p_{\theta}}(r,p_{r})$ is the Hamiltonian (34) for fixed $p_{\theta}$
and

$n(E,\hbar;p_{\theta})\sim\frac{1}{2\pi\hbar^{2}}\int_{H_{p_{\theta}}(r,p_{r})\leq
E}dr\,dp_{r}$ (36)

may be interpreted as the semiclassical number of states of given $p_{\theta}$
with $H_{p_{\theta}}\leq E$. The limits of integration $p_{\theta}^{\rm
min,max}(E)$ are obtained by extremizing with respect to $r$ and $p_{r}$ in

$p_{\theta}^{\pm}(r,p_{r};E)=\frac{r}{2}\left[-km\lambda
r\pm\sqrt{-4p_{r}^{2}-(\lambda
kr^{2}-2p_{z})^{2}+8E\mu-4k^{2}\mu(m^{2}+r^{2})}\right].$ (37)

The expressions for $p_{\theta}^{\pm}$ are got by solving the quadratic
equation $H=E$ for $p_{\theta}$. The extrema occur at $p_{r}=0$ and values of
$r$ determined via the roots of a quartic equation in $r^{2}$:

$\displaystyle 9k^{4}\lambda^{4}r^{8}-4k^{3}\lambda^{2}(12\lambda
p_{z}-k\lambda^{2}m^{2}-12k\mu)r^{6}$ (38)
$\displaystyle-8k^{2}[6E\lambda^{2}\mu-11\lambda^{2}p_{z}^{2}+2kp_{z}\lambda(m^{2}\lambda^{2}+8\mu)-k^{2}\mu(5m^{2}\lambda^{2}+8\mu)]r^{4}$
(39)
$\displaystyle+16k(km^{2}\lambda^{2}+4k\mu-4p_{z}\lambda)(p_{z}^{2}-2E\mu+k^{2}m^{2}\mu)r^{2}+16(p_{z}^{2}-2E\mu+k^{2}m^{2}\mu)^{2}=0.\quad$
(40)

We evaluate the integral in (36) numerically for several energies $E$ for
fixed values of $p_{\theta},p_{z}$ and other parameters. The resulting
distributions $n(E,\hbar;p_{\theta})$ are compared in Fig. 3(a) with those
obtained from numerical diagonalization. We find that the semiclassical
estimate for $n(E,\hbar=1;p_{\theta})$ agrees quite well with the numerical
spectrum for each allowed value of $p_{\theta}$ and $E\gg k^{2}m^{2}$. To
obtain $n(E)$, we numerically integrate over $p_{\theta}$ in (35). The
resulting cumulative distribution is compared with that from numerical
diagonalization in Fig. 3(b).

(a)

(b)

Figure 3: (a) Semiclassical cumulative distribution $n(E,\hbar;p_{\theta})$
(squares) nearly matches that from numerical diagonalization (stars) for
$p_{\theta}=0,50,100,150$. (b) Comparison of semiclassical (squares) and
numerical (disks) $n(E,\hbar)$ after combining all allowed values of
$p_{\theta}$. In all cases we take $p_{z}=1$ and $k=\mu=m=\lambda=\hbar=1$.

#### Power law behavior of $n(E;p_{\theta})$ and $n(E)$.

For large $E$, we find that the cumulative level distribution satisfies a
power law $n(E;p_{\theta})\propto(E-E_{0})^{\zeta}$ where $E_{0}(p_{\theta})$
is the ground state energy. Since the term $\lambda kmp_{\theta}/2\mu$ linear
in $p_{\theta}$ is a constant addition to the Hamiltonian (34), it does not
contribute to energy differences. Thus, the power $\zeta(p_{\theta})$ depends
only on $|p_{\theta}|$. We find that it increases monotonically with
$|p_{\theta}|$. Moreover, after combining all angular momentum sectors, $n(E)$
also satisfies a power law for large $E$. For instance, when
$\hbar=k=m=\mu=\lambda=p_{z}=1$, $\zeta(p_{\theta})\approx
0.75,0.78,0.81,0.84,0.88$ for $p_{\theta}=0,\pm 20,\pm 50,\pm 100,\pm 150$
while $n(E)\sim E^{1.5}$ for large $E$.

###  5.2 Level spacing distributions

The quantum energy levels of a classically integrable system are expected to
display certain universal statistical features such as Poissonian level
spacing distributions [8]. To remove the system-specific information and
extract universal properties111To study spectral statistics, we note that
there is no need to ‘purify’ the screwon spectrum obtained in §4. Reflection
through the origin is part of the circular symmetry of (7). of the level
spacings, we ‘unfold’ the energy spectrum so that the mean level spacing is
one [9, 10]. In order to do this, we define unfolded levels
$\xi_{i}=\xi(E_{i})$ for $i=1,2,3,\ldots,N$, where $N$ is the number of levels
considered ($N\approx 20000$ below). The function $\xi$ is chosen to be the
polynomial (typically of degree 5 or 6) that best fits the cumulative energy
distribution function introduced in (32):

$n(E)=\sum_{i=1}^{N}\theta(E-E_{i})=\xi(E)+\eta(E).$ (41)

Thus, $\xi(E)$ is a smoothed version of the ‘staircase function’ $n(E)$, which
satisfies $n(E_{j})=j$ for $j=1,2,\ldots,N$ assuming there are no degenerate
levels. [More generally, if the first $n_{1}$ energies are degenerate
($E_{1}=E_{2}=\cdots=E_{n_{1}}$) and next $n_{2}$ energies are degenerate
($E_{n_{1}+1}=E_{n_{1}+2}=\cdots=E_{n_{1}+n_{2}}$) etc., then
$n(E_{1})=\cdots=n(E_{n_{1}})=n_{1}$ and
$n(E_{n_{1}+1})=\cdots=n(E_{n_{1}+n_{2}})=n_{1}+n_{2}$ and so forth.] Here
$\eta(E)$ contains fluctuations. By construction, the unfolded levels
$\xi_{1},\xi_{2},\ldots$ have a mean spacing of approximately one.

For $i=1,2,3,\ldots$, the $i^{\rm th}$-nearest neighbor (n.n.) spacing
distribution is defined as the normalized distribution of values of
$s=(\xi_{j+i}-\xi_{j})$ for $j=1,2,\ldots N-i$. In Fig. 4, we plot histograms
of the spacing probability distributions for $i=1,2,3$. We find that the n.n.
spacing distribution follows the exponential law $P_{1}(s)=\exp(-s)$, while
the higher order spacing distributions are well approximated by the Poisson
distribution:

$P_{i}(s)=\frac{s^{i-1}\exp(-s)}{(i-1)!}\quad\text{for}\;\;i=1,2,3,\ldots.$
(42)

This Poisson statistics is expected if the n.n. spacing is exponentially
distributed and successive spacings can be treated as independent random
variables.

(a)

(b)

(c)

Figure 4: First, second and third order normalized spacing distributions of
lowest $N\approx 20000$ unfolded energy levels for $\lambda=p_{z}=\mu=k=m=1$.
The solid line is from (42), while the dots represent binheights of the
normalized spacing histograms with binwidths 0.2. The spacing distributions
are seen to obey Poisson statistics as expected of an integrable system. In
fact, we find that Poisson statistics applies for all nonzero coupling
$\lambda$.

Remark: We find that when we use the unfolded spectrum, the n.n. spacing
histograms for all $\lambda>0$ have roughly the same range ($0\leq s\leq 4$).
Moreover, the same binwidth (anything from 0.1 to 1) gives the universal
exponential distribution for all $\lambda>0$. If we do not unfold the spectra,
then the mean spacing varies with $\lambda$ so the spacing distributions have
$\lambda$ dependent ranges. However, it is still possible to extract the
universal exponential spacing distribution provided the binwidth is chosen
appropriately: a different range of binwidths works for different values of
$\lambda$.

Recently, ratios of level spacings have emerged as a useful statistical tool
to study quantum spectra [11, 12]. The $i^{\rm th}$ nonoverlapping spacing
ratio distribution is defined as the normalized distribution of values of

$r_{j}^{(i)}=\frac{E_{j+i}-E_{j}}{E_{j-1}-E_{j-1-i}}\quad\text{for}\quad
j=i+2,\ldots,N-i.$ (43)

If one assumes that the nonoverlapping spacings are independently distributed,
then the probability that the spacing in the numerator is $s^{\prime}$ and
that in the denominator is $s$, is simply the product of individual spacing
probabilities. Thus, the probability that the ratio of spacings is
$s^{\prime}/s=r$ is

$Q_{i}(r)=\iint
P_{i}(s^{\prime})P_{i}(s)\delta\left(\frac{s^{\prime}}{s}-r\right)ds\;ds^{\prime}=\int
P_{i}(rs)P_{i}(s)|s|\>ds=\frac{(2i-1)!}{[(i-1)!]^{2}}\frac{r^{i-1}}{(1+r)^{2i}}.$
(44)

(a)

(b)

(c)

(d)

Figure 5: Histograms of $i^{\rm th}$ nonoverlapping spacing ratio
distributions for $i=1,2,3,4$ for lowest 20000 screwon energy levels with
$\lambda=p_{z}=k=\mu=m=1$. They are seen to match the distributions $Q_{i}(r)$
(44) denoted by solid curves. Unfolding the spectrum does not significantly
affect these histograms.

###  5.3 Number variance

Associated to the unfolded spectrum $\xi_{1},\xi_{2},\ldots\xi_{N}$, we have
the Dirac-comb spectral density function
$d(\xi)=\sum_{i=1}^{N}\delta(\xi-\xi_{i})$. Given a spectral window
$[\xi-L/2,\xi+L/2]$, the number of energy levels it contains is:

$n(\xi,L)=\int_{\xi-L/2}^{\xi+L/2}d(\xi^{\prime})\;d\xi^{\prime}.$ (45)

The ensemble average $\langle n(\xi,L)\rangle$ is defined as the mean value of
$n(\xi,L)$ as $\xi$ ranges over an appropriate portion of the available
spectrum:

$\langle
n(\xi,L)\rangle=(N-{2n_{0}})^{-1}\sum_{j=n_{0}}^{N-n_{0}}n(\xi_{j},L).$ (46)

Here, $n_{0}$ is chosen so that $[\xi-L/2,\xi+L/2]$ always lies within the
available spectrum. A more accurate approach is to define the ensemble average
as $(\xi_{N-n_{0}}-\xi_{n_{0}})^{-1}$ times the integral of $n(\xi,L)$ over
all values of $\xi$ between $\xi_{n_{0}}$ and $\xi_{N-n_{0}}$. However, the
two definitions lead to roughly the same ensemble average as long as $L$ is
not too small and the difference is hardly visible in our plots. If the
unfolded spectrum has approximately unit mean spacing, we expect that $\langle
n(\xi,L)\rangle\approx L$. In Fig. 6(a), we plot $\langle n(\xi,L)\rangle$ for
$1\leq L\leq 320$ by performing an ensemble average over the lowest $N=20,000$
screwon levels, with $n_{0}=500$. $\langle n(\xi,L)\rangle$ is seen to be
$\approx L$ except for very small values of $L$, thus validating the unfolding
procedure.

To study fluctuations in the number of levels around its mean, we consider the
number variance [13]

$\Sigma^{2}(L)=\left<[n(\xi,L)-\left<n(\xi,L)\right>]^{2}\right>,$ (47)

where the ensemble average is performed as for $\langle n(\xi,L)\rangle$. In
Fig 6(b), we plot $\Sigma^{2}(L)$ for the lowest $N=20,000$ screwon levels of
the RR model, with $\lambda=1$.

(a)

(b)

Figure 6: (a) Ensemble average of number of unfolded levels in an interval of
length $L$ follows $\langle n(\xi,L)\rangle\approx L$. (b) Number variance
$\Sigma^{2}(L)$ grows linearly $\approx L$ for $0\leq L\lesssim 20$ and then
oscillates with a wavelength of $\approx 200$ while saturating. Here,
$\lambda=p_{z}=m=k=\mu=1$.

For small $L\lesssim 20$, the number variance $\Sigma^{2}\approx L$. For
larger $L$, $\Sigma^{2}$ saturates and oscillates.

Here we interpret the result that $\Sigma^{2}(L)\approx L$ for small $L$ in
terms of the 2-point correlation function of the spectral density $d(\xi)$.
Indeed, note that

$\Sigma^{2}(L)=\left<(n(\xi,L)-\left<n(\xi,L)\right>)^{2}\right>=\langle
n^{2}\rangle+\langle n\rangle^{2}-2\langle n\langle n\rangle\rangle.$ (48)

Now using $\langle n\rangle\approx L$ and $\Sigma^{2}(L)\approx L$, we deduce
that $\langle n(\xi,L)^{2}\rangle=L^{2}+L$. From the definition

$\langle
n(\xi,L)^{2}\rangle=\int_{\xi-L/2}^{\xi+L/2}d\xi_{1}\int_{\xi-L/2}^{\xi+L/2}d\xi_{2}\langle
d(\xi_{1})d(\xi_{2})\rangle,$ (49)

we see that if

$\langle d(\xi_{1})d(\xi_{2})\rangle=1+\delta(\xi_{1}-\xi_{2})$ (50)

then $\langle n(\xi,L)^{2}\rangle=L^{2}+L$. This may further be interpreted as
saying the connected 2-point correlation function of $d(\xi)$ is a delta
function. To see this, recall that

$\langle d(\xi_{1})d(\xi_{2})\rangle_{c}=\langle
d(\xi_{1})d(\xi_{2})\rangle-\langle d(\xi_{1})\rangle\langle
d(\xi_{2})\rangle.$ (51)

Thus we see that if $\langle
d(\xi_{1})d(\xi_{2})\rangle_{c}=\delta(\xi_{1}-\xi_{2})$ and $\langle
d(\xi_{1})\rangle=1$, then $\langle n(\xi,L)^{2}\rangle=L^{2}+L$ and
consequently $\Sigma^{2}(L)\approx L$. The condition $\langle
d(\xi_{1})\rangle=1$ is expected since the mean level density of the unfolded
spectrum is 1. Of course, this linear growth of $\Sigma^{2}(L)$ is valid only
for small $L$.

It would be interesting to understand the subsequent saturation and
oscillations in $\Sigma^{2}(L)$ possibly in terms of short periodic orbits of
the Rajeev-Ranken model in a semiclassical approximation.

###  5.4 Spectral rigidity

The spectral rigidity $\Delta_{3}(L,E)$ [14] is a measure of fluctuations in
the cumulative spectral distribution (staircase function) $n(E)$. It is
defined as the ‘local average’ of the mean square deviation of the best fit
straight line to $n(E)$ (41) over an energy window
$W=[E-L/2\bar{d}(E),E+L/2\bar{d}(E)]$:

$\displaystyle\Delta_{3}(L,E)$ $\displaystyle=$
$\displaystyle\langle\sigma_{*}(E,L)\rangle\quad\text{where}\quad\sigma_{*}(E,L)={\rm
min}_{A,B}\sigma(E,L,A,B)\quad\text{and}$ (52) $\displaystyle\sigma(E,L,A,B)$
$\displaystyle=$
$\displaystyle(\bar{d}(E)/L)\int_{W}\left[n(E+\epsilon)-(A\epsilon+B)\right]^{2}{\rm
d}\epsilon.$ (53)

Given a central energy $E$, the best fit is computed over an energy range $L$
in units of the local mean spacing $\bar{d}(E)$. We take $\bar{d}(E)$ to be
the smoothed density of states

$\bar{d}(E)=\frac{1}{2\delta}\int_{E-\delta}^{E+\delta}d(E^{\prime}){\rm
d}E^{\prime},$ (54)

which is a measure of the mean spacing between levels in the $\delta$-vicinity
of energy $E$. We comment on the choice of $\delta$ below. Having obtained the
best fit straight line $A_{*}\epsilon+B_{*}$, we calculate the extremal mean
square deviation $\sigma_{*}$ for the chosen $E$ and $L$. Finally, the local
average indicated by $\langle\cdots\rangle$ in (53) is evaluated by averaging
$\sigma_{*}(E,L)$ over a range of energy levels for fixed $L$.

(a)

(b)

(c)

Figure 7: (a) Spectral rigidity $\Delta_{3}(L,E)$ for the screwon spectrum for
central energies $E=300,500,600$ and $750$ (in units where $\mu=m=\hbar=1$)
follow the universal $L/15$ law for small $L$ and then saturate around
$\Delta_{3}^{\rm sat}(E)$ beyond $L=L_{\rm sat}(E)$. Here, for definiteness
$L_{\rm sat}(E)$ is defined as the largest value of $L$ for which
$\Delta_{3}(L,E)$ intersects the $L/15$ straight-line for a given $E$.
Moreover, $\Delta_{3}^{\rm sat}(E)$ is defined as the mean value of
$\Delta_{3}$ over the range $L_{\rm sat}\leq L\leq L_{\rm max}=600$. (b) and
(c) Log-Log plots of $L_{\rm sat}$ and $\Delta_{3}^{\rm sat}$ vs $E$ show that
both display $\sqrt{E}$ behavior. In all plots we have chosen
$\lambda=k=p_{z}=1$.

For the quantum RR model, we have calculated (see Appendix A) the spectral
rigidity $\Delta_{3}(L,E)$ for $L\leq 600=L_{\rm max}$ for several central
energies $E=300,500,600,750$ and displayed the results in Fig. 7. In
estimating the local mean spacing $\bar{d}(E)$, we choose $\delta$ to
accommodate about 80 levels centered at $E$. The value of $\bar{d}(E)$ is
largely insensitive to the choice of $\delta$. On the other hand, the local
average is performed over an energy interval $[E-20,E+20]$. The results are
insensitive to small changes in the energy range for the local average. In
fact, essentially the same results are obtained if the local average is
performed over an ensemble of about 2000 levels centered at $E$. Although
$\Delta_{3}(L,E)$ generally depends on the central energy $E$, for small $L$
it approaches the universal linear shape $\Delta_{3}(L,E)\to L/15$ as expected
from Berry’s arguments [14] for integrable systems based on the semi-classical
trace formula applied to long-periodic orbits. Beyond a critical value of
$L=L_{\rm sat}$, the spectral rigidity oscillates around a saturation value
$\Delta_{3}^{\rm sat}(E)$. Both $L_{\rm sat}$ and $\Delta_{3}^{\rm sat}$
increase when the central energy $E$ is significantly augmented [due to
fluctuations, this monotonicity can fail for energies that are not widely
separated]. In fact, as shown in Figs. 7(b) and 7(c), we find that both
$L_{\rm sat}$ and $\Delta^{\rm sat}_{3}$ are approximately proportional to the
square-root of energy. The latter power law is expected from the work of
Casati, Chirikov and Guarneri [13, 15].

##  6 Dispersion relation at strong coupling

In [6], we had conjectured a $(\lambda k)^{2/3}$ power law dispersion relation
for the energies of highly excited screwons at strong coupling based on a
numerical inversion of the WKB quantization condition. Here, we generalize
this conjecture, make it more precise and provide additional evidence for it
going beyond the WKB approximation.

Recall that the quartic potential $U=\alpha r^{2}+\beta r^{4}$ in the
Hamiltonian (7) involves the coefficients $\alpha$ and $\beta$ (10). At strong
coupling, by which we mean

$\lambda\gg 2\sqrt{\mu}/m\quad\text{and}\quad\lambda\gg\mu k/p_{z},$ (55)

we may neglect the $k^{2}/2$ term in $\alpha$, so that $U$ depends on
$\lambda$ and $k$ only through their product. In fact, in this limit, the
shifted radial Hamiltonian $H-k^{2}m^{2}/2-p_{z}^{2}/2\mu$ (9) also depends on
$\lambda$ and $k$ only through their product. Thus, the shifted energy
eigenvalues $E-k^{2}m^{2}/2-p_{z}^{2}/2\mu$ determined by the radial equation

$-\frac{\hbar^{2}}{2\mu}\left(\frac{{\rm d^{2}}}{{\rm
d}r^{2}}+\frac{1}{r}\frac{{\rm d}}{{\rm
d}r}-\frac{l^{2}}{r^{2}}\right)\rho(r)+\frac{\lambda^{2}k^{2}}{8\mu}\left(\left(m^{2}-\frac{4p_{z}}{\lambda
k}\right)r^{2}+r^{4}\right)\rho+\frac{\hbar l\lambda
km}{2\mu}\rho\approx\left(E-\frac{k^{2}m^{2}}{2}-\frac{p_{z}^{2}}{2\mu}\right)\rho$
(56)

must be a function of $\lambda k$. More precisely, on dimensional grounds,

$E_{n,l}\approx\frac{p_{z}^{2}}{2\mu}+\frac{k^{2}m^{2}}{2}+\frac{\hbar^{2}}{\mu
m^{2}}\epsilon_{n,l}\left(\tilde{\sigma},\frac{mp_{z}}{\hbar}\right)\quad\text{for
large}\quad\lambda\quad\text{satisfying}\quad(\ref{e:strong-coupling-
conditions}).$ (57)

Here, $\hbar^{2}\epsilon_{n,l}/\mu m^{2}$ are the eigenvalues of
$H-k^{2}m^{2}/2-p_{z}^{2}/2\mu$ for large $\lambda$ (satisfying (55)) and
$\epsilon_{n,l}$ is a function of the dimensionless variables
$\tilde{\sigma}=\lambda km^{3}/\hbar$ and $mp_{z}/\hbar$. In the limit (55)
the other independent dimensionless combination $\sqrt{\mu}/m\lambda$ does not
enter as it tends to $0$. Moreover, if

$\frac{4p_{z}}{m^{2}\lambda}\ll k\ll\frac{\lambda p_{z}}{\mu},$ (58)

then $\epsilon_{n,l}$ becomes independent of $p_{z}$ as the $4p_{z}/\lambda k$
term can be ignored relative to $m^{2}$ in (56).

From numerical diagonalization of $H$, we find for fixed $l$ and $p_{z}$ at
strong coupling and moderate $k$ satisfying (55) and (58) that
$\epsilon_{n,l}$ depends on both $\lambda$ and $k$ via power laws. For the
ground state ($n=l=0$) the powers of $\lambda$ and $k$ are approximately equal
to $1$. The powers decrease with increasing level number $n$ (holding $l=0$
fixed) and rapidly approach the common value $\eta=2/3$ for sufficiently
excited states ($n\gtrsim 100$). It would be interesting to analytically
understand the emergence of this $2/3$ power law for highly excited screwons
at strong coupling.

We have not investigated the $l$-dependence of the exponent in detail.
However, we numerically verified that for $n\gg|l|$,
$\epsilon_{n,l}\propto(\lambda k)^{\eta}$ where $\eta\approx 2/3$ is
independent of $l$. For instance, this holds if $n\gtrsim 200$ and $|l|\leq
20$, for $0.08\leq k\leq 5$ holding $\lambda=50$ fixed as well as for
$4\lesssim\lambda\lesssim 35$ holding $k=1$ fixed. However, if $|l|$ is
comparable to or greater than $n$, then $\eta$ depends on $l$. For example, if
$n=0$, then $\eta\approx 1,0.4,0.9$ for $l=0,-20,20$.

(a)

(b)

(c)

(d)

Figure 8: $\log\epsilon_{n,l}$ vs $\log\lambda$ for (a) $n=0$ and (b) n = 499
with $l=0$ when $k=1$. For sufficiently large $\lambda$, the slope of the
fitted straight line decreases from $1$ to $2/3$ as we progress from the
ground state to excited states. $\log\epsilon_{n,l}$ vs $\log k$ for (c) $n=0$
and (d) $n=499$ with $l=0$ for $\lambda=50$. For $k$ satisfying (58) and
sufficiently large $\lambda$, we find that $\epsilon_{n,l}\propto k^{\eta}$
with the power decreasing from $\eta\approx 1$ to $\eta\approx 2/3$ as we
progress from the ground to excited states. In all cases, we take
$p_{z}=m=\mu=\hbar=1$.

##  7 Discussion

There are some interesting directions of research which arise from this work.
Although we have recovered the expected universal behaviour of number variance
and spectral rigidity at small $L$, they both display system-dependent
saturation and oscillations for larger $L$. We would like to understand this
nonuniversal behaviour using Gutzwiller’s trace formula and Berry’s semi-
classical asymptotic theory [14]. For this purpose, we intend to use the exact
solutions of the classical RR model [2, 4, 5] to identify and classify the
shortest periodic orbits. It would also be interesting to understand the power
law behaviour in the saturation of spectral rigidity and the asymptotic
behaviour of the cumulative level distribution function semiclassically. In
another direction, we would like to understand analytically the common power
law behaviour of the dispersion relation for screwons at strong coupling. This
dispersion relation should shed light on the ultraviolet behaviour of screwons
in the field theory.

Acknowledgements: We would like to thank M. V. Berry, M.S. Santhanam and H.
Senapati for useful discussions and references. This work was supported in
part by the Infosys Foundation and grants (MTR/2018/000734, CRG/2018/002040)
from the Science and Engineering Research Board, Govt. of India.

##  Appendix A Evaluating the spectral rigidity

The integrals appearing in the evaluation of spectral rigidity (53) are in
fact finite sums. Writing them as such speeds up numerical calculations. Let
us denote the limits of integration by $E_{\rm min,\rm
max}=E\mp(L/2\bar{d}(E))$. Then $(L/\bar{d}(E))\sigma=I_{1}+I_{2}+I_{3}$,
where

$I_{1}=\int_{E_{\rm min}}^{E_{\rm
max}}(A\epsilon+B)^{2}d\epsilon=\left[\frac{A^{2}\epsilon^{3}}{3}+AB\epsilon^{2}+B^{2}\epsilon\right]_{E_{\rm
min}}^{E_{\rm max}}.$ (59)

Suppose $n_{0}$ is the number of levels with energy $\leq E_{\rm min}$ and $N$
is the number of levels between $E_{\rm min}$ and $E_{\rm max}$, then

$I_{2}=\int_{E_{\rm min}}^{E_{\rm
max}}n(E+\epsilon)^{2}d\epsilon=n_{0}^{2}(E_{n_{0}+1}-E_{\rm
min})+(n_{0}+N)^{2}(E_{\rm
max}-E_{n_{0}+N})+\sum_{k=1}^{N-1}(n_{0}+k)^{2}S_{n_{0}+k}.$ (60)

Here, $S_{n_{0}+k}=E_{n_{0}+k+1}-E_{n_{0}+k}$ is the level spacing. Similarly,
we have

$\displaystyle\frac{I_{3}}{2}$ $\displaystyle=$ $\displaystyle\int_{E_{\rm
min}}^{E_{\rm
max}}n(E+\epsilon)(A\epsilon+B)d\epsilon=B\left[n_{0}(E_{n_{0}+1}-E_{\rm
min})+(n_{0}+N)(E_{\rm
max}-E_{n_{0}+N})+\sum_{k=1}^{N-1}(n_{0}+k)S_{n_{0}+k}\right]$ (62)
$\displaystyle+\frac{A}{2}\left[n_{0}(E_{n_{0}+1}^{2}-E_{\rm
min}^{2})+(n_{0}+N)(E_{\rm
max}^{2}-E_{n_{0}+N}^{2})+\sum_{k=1}^{N-1}(n_{0}+k)E_{n_{0}+k+1}^{2}-E_{n_{0}+k}^{2}\right].\qquad$

## References

  * [1]
  * [2] S. G. Rajeev and E. Ranken, Highly nonlinear wave solutions in a dual to the chiral model, Phys. Rev. D $\mathbf{93}$, 105016 (2016).
  * [3] A. M. Polyakov, Gauge fields and strings, Harwood Academic Publishers, Chur (1987).
  * [4] G. S. Krishnaswami and T. R. Vishnu, On the Hamiltonian formulation and integrability of the Rajeev-Ranken model, J. Phys. Commun. $\mathbf{3}$, 025005 (2019).
  * [5] G. S. Krishnaswami and T. R. Vishnu, Invariant tori, action-angle variables and phase space structure of the Rajeev-Ranken model, J. Math. Phys. $\mathbf{60}$, 082902 (2019).
  * [6] G. S. Krishnaswami and T. R. Vishnu, Quantum Rajeev-Ranken model as an anharmonic oscillator, J. Math. Phys. $\mathbf{63}$, 032101 (2022).
  * [7] W. Pauli, Wave mechanics, Volume 5 of Pauli’s lectures on physics, Dover Publications, New York (2000).
  * [8] M.V. Berry and M. Tabor, Level clustering in the regular spectrum, Proc. R. Soc. A $\mathbf{356}$, 375 (1977).
  * [9] F. Haake, Quantum signatures of chaos, Springer-Verlag, Berlin (2010).
  * [10] T. Guhr, A. Müller-Groeling and H. A. Weidenmüller, Random matrix theories in quantum physics: common concepts, Phys. Rep. $\mathbf{299}$,189 (1998).
  * [11] S. H. Tekur, S. Kumar and M. S. Santhanam, Exact distribution of spacing ratios for random and localized states in quantum chaotic systems, Phys. Rev. E $\mathbf{97}$, 062212 (2018).
  * [12] S. H. Tekur, U. T. Bhosale and M. S. Santhanam, Higher order spacing ratios in random matrix theory and complex quantum systems, Phys. Rev. B $\mathbf{98}$,104305 (2018).
  * [13] H. J. Stöckmann, Quantum chaos, an introduction, Cambridge University Press, Cambridge (1999).
  * [14] M. V. Berry, Semi-classical theory of spectral rigidity, Proc. R. Soc. Lond. A $\mathbf{400}$, 229 (1985).
  * [15] G. Casati, B. V. Chirikov and I. Guarneri, Energy level statistics of integrable quantum systems, Phys. Rev. Lett. $\mathbf{54}$, 1350 (1985).