This version has been removed by arXiv administrators as the submitter did not
have the right to agree to the license at the time of submission [arXiv admin
2024-03-07].


# On Airy Solutions of PII and Complex Cubic Ensemble of Random Matrices, I

Ahmad Barhoumi  Department of Mathematics, Royal Institute of Technology
(KTH), Stockholm, Sweden<EMAIL_ADDRESS>,  Pavel Bleher  Department of
Mathematical Sciences, Indiana University Indianapolis, 402 North Blackford
Street, Indianapolis, IN 46202, USA<EMAIL_ADDRESS>,  Alfredo Deaño
Department of Mathematics, Universidad Carlos III de Madrid, Avda. de la
Universidad 30, 28911 Leganés, Madrid, Spain<EMAIL_ADDRESS>and
Maxim Yattselev  Department of Mathematical Sciences, Indiana University
Indianapolis, 402 North Blackford Street, Indianapolis, IN 46202, USA
<EMAIL_ADDRESS>
###### Abstract.

We show that the one-parameter family of special solutions of PII, the second
Painlevé equation, constructed from the Airy functions, as well as associated
solutions of PXXXIV and SII, can be expressed via the recurrence coefficients
of orthogonal polynomials that appear in the analysis of the Hermitian random
matrix ensemble with a cubic potential. Exploiting this connection we show
that solutions of PII that depend only on the first Airy function
$\mathrm{Ai}$ (but not on $\mathrm{Bi}$) possess a scaling limit in the pole
free region, which includes a disk around the origin whose radius grows with
the parameter. We then use the scaling limit to show that these solutions are
monotone in the parameter on the negative real axis.

###### Key words and phrases:

Airy functions, Painlevé equations, random matrix models

###### 2020 Mathematics Subject Classification:

15B52, 33C10, 33C47

The third author acknowledges financial support from Dirección General de
Investigación e Innovación, Consejería de Educación e Investigación of
Comunidad de Madrid (Spain), and Universidad de Alcalá under grant
CM/JIN/2021-014, as well as from grant PID2021-123969NB-I00, funded by
MCIN/AEI/10.13039/501100011033, and from grant PID2021-122154NB-I00 from
Spanish Agencia Estatal de Investigación. The research of the last author was
supported in part by a grant from the Simons Foundation, CGM-706591.

##  1\. Introduction

At the beginning of the 20th century, it was shown by Painlevé and Gambier
that among the second order differential equations of the form

$q^{\prime\prime}=F(q^{\prime},q,z),$

where $F$ is rational in $q,q^{\prime}$ and analytic in $z$, there are exactly
$50$ canonical equations whose solutions do not possess movable branch points.
Among these, $44$ can either be reduced to linear equations, equations solved
in terms of elliptic functions, or to the remaining $6$ equations. In this
work we are interested in what is now known as the second Painlevé equation,
see (2). It was shown by Gambier [18] that PII has a one-parameter family of
solutions that are expressible in terms of the Airy functions and their
derivatives. These special function solutions appear in many contexts. For
example, in Hermitian random matrix ensembles as the asymptotic analysis of
eigenvalue behavior (in a suitable double scaling limit as the size of the
matrices tends to infinity and near the edge of the spectrum) makes use of a
correlation kernel that in particular cases can be written in terms of Airy
solutions of PII and PXXXIV, see [21]. Another appearance of these functions
in random matrix theory comes from Hermitian ensembles with cubic potential.
The cubic random matrix model, given by the probability distribution

$\tilde{Z}_{N}^{-1}(u)e^{-N\mathrm{Tr}(M^{2}/2-uM^{3})}dM$

on $N\times N$ Hermitian matrices $M$, where $u$ is a parameter, has been
investigated in physical literature by Brézin, Itzykson, Parisi, and Zuber [9]
and Bessis, Itzykson, and Zuber [4]. There it was recognized that the free
energy of this model must possess an expansion in the parameter $N$. The
$N^{-2g}$ coefficient of this expansion, as a function of $u$, possesses a
Taylor series around the origin, whose coefficients count the number of the
three-valent graphs on a Riemann surface of genus $g$. These considerations
were later made rigorous by the middle two authors of the present work [5].
Further investigations of this model were carried out by Huybrechs, Kuijlaars,
and Lejon in [20, 19] and by the authors in [6, 7, 3], see also [1, 2] for the
numerical studies. The main tool of investigation in these works was the
connection to non-Hermitian orthogonal polynomials with respect to the cubic
weight $\exp(-z^{3}/3+tz)$ (parameter $t$ is connected to the parameter $u$),
whose moments can be expressed via the Airy functions $\mathrm{Ai}(z)$ and
$\mathrm{Bi}(z)$ as well as their derivatives. In this work we utilize this
connection to show what the results of [7] can say about the solutions of PII
that depend only on $\mathrm{Ai}(z)$ and its derivatives. The techniques of
[5, 6, 7, 3] can be extended to study all special solutions and we shall do so
in a subsequent publication.

##  2\. Airy Solutions of PII

Each of the Painlevé equations PI–PVI can be written as a Hamiltonian system

$\frac{dq}{dz}=\frac{\partial H}{\partial
p}\quad\text{and}\quad\frac{dp}{dz}=-\frac{\partial H}{\partial q},$

see [24, 23, 22]. In the case of PII, the Hamiltonian $H_{\mathrm{II}}$ is
equal to

$H_{\mathrm{II}}(q,p,z;\alpha)=\frac{1}{2}p^{2}-\left(q^{2}+\frac{z}{2}\right)p-\left(\alpha+\frac{1}{2}\right)q,$

where $\alpha$ is a parameter. Hence, the Hamiltonian system becomes

(1) $q^{\prime}=p-q^{2}-z/2\quad\text{and}\quad p^{\prime}=2pq+\alpha+1/2.$

Eliminating $p$ from (1) gives PII:

(2) $q^{\prime\prime}=2q^{3}+zq+\alpha,$

while eliminating $q$ from (1) yields PXXXIV:

(3) $p^{\prime\prime}=\frac{(p^{\prime})^{2}-(\alpha+1/2)^{2}}{2p}+2p^{2}-zp.$

Notice that the first equation (1) allows one to express $p$ through
$q,q^{\prime}$ and $z$ while the second equation allows to express $q$ through
$p,p^{\prime}$. These connection formulae provide one-to-one correspondence
between solutions of PII and PXXXIV. Morever, if
$p_{\alpha+1/2}(z),q_{\alpha+1/2}(z)$ are solutions of (1) and

(4)
$\sigma_{\alpha+1/2}(z):=H_{\mathrm{II}}\big{(}q_{\alpha+1/2}(z),p_{\alpha+1/2}(z),z;\alpha\big{)},$

then this function solves SII, the Jimbo-Miwa-Okamoto $\sigma$-form of PII:

(5)
$(\sigma^{\prime\prime})^{2}+4(\sigma^{\prime})^{3}+2\sigma^{\prime}(z\sigma^{\prime}-\sigma)=(\alpha/2+1/4)^{2}.$

Conversely, if $\sigma_{\alpha+1/2}(z)$ solves (5), then the functions

(6)
$q_{\alpha+1/2}(z)=\frac{2\sigma_{\alpha+1/2}^{\prime\prime}(z)+\alpha+1/2}{4\sigma_{\alpha+1/2}^{\prime}(z)}\quad\text{and}\quad
p_{\alpha+1/2}(z)=-2\sigma_{\alpha+1/2}^{\prime}(z),$

solve (2) and (3), respectively, see [24, 23, 22, 25].

It is known from the work of Gambier, see [18], that (2), and therefore (3)
and (5) via the correspondence (1), (4), and (6), has a one-parameter family
of solutions expressible through the Airy functions and their derivatives if
and only if $\alpha+1/2$ is an integer. Recall that the standard Airy
functions can be given by their integral representations

(7)
$\mathrm{Ai}(z)=\frac{1}{2\pi\mathrm{i}}\int_{L_{1}-L_{2}}e^{-V(s;z)}ds,\quad
V(s;z):=-\frac{s^{3}}{3}+sz,$

where $L_{k}:=\big{\\{}xe^{\pi\mathrm{i}(-1+2k/3)}:x\in(0,\infty)\big{\\}}$,
$k\in\\{0,1,2\\}$, are rays oriented towards the origin, and

(8)
$\mathrm{Bi}(z)=\frac{1}{2\pi}\left(\int_{L_{0}-L_{1}}+\int_{L_{0}-L_{2}}\right)e^{-V(s;z)}ds.$

When $\alpha=1/2$, it can be shown that (2) has a one-parameter family of
solutions given by

$q_{1}(z;\lambda)=-\frac{d}{dz}\log\left(C_{1}\,\mathrm{Ai}\left(-2^{-1/3}z\right)+C_{2}\,\mathrm{Bi}\left(-2^{-1/3}z\right)\right),$

where $\lambda\in\overline{\mathbb{C}}$ and $C_{1},C_{2}\in\mathbb{C}$ are any
numbers such that $C_{2}/C_{1}=\lambda\in\mathbb{C}$ and $C_{1}=0$ when
$\lambda=\infty$ (the branch of the logarithm is not important as we
immediately take $z$ derivative). The solutions $q_{n}(z;\lambda)$
corresponding to $\alpha=n-1/2$, $n\geq 1$, are constructed recursively via
Bäcklund transformation

$q_{n+1}(z;\lambda)=-q_{n}(z;\lambda)-\frac{2n}{2q_{n}^{2}(z;\lambda)+2q_{n}^{\prime}(z;\lambda)^{2}+z}.$

The solutions for non-positive values of $n$ are then obtained via
$q_{n}(z;\lambda)=-q_{-n+1}(z;\lambda)$. However, the above formula is not the
most convenient way of expressing these solutions. To this end, define
$\tau_{0}(z;\lambda)\equiv 1$ and

(9)
$\tau_{n}(z;\lambda):=\mathrm{det}\left[\frac{d^{j+k}}{dz^{j+k}}\left(C_{1}\,\mathrm{Ai}\left(-2^{-1/3}z\right)+C_{2}\,\mathrm{Bi}\left(-2^{-1/3}z\right)\right)\right]_{j,k=0}^{n-1}$

when $n\in\mathbb{N}$. Then, see [14, 25, 10], the solutions
$q_{n}(z;\lambda)$, $n\geq 1$, and the corresponding functions
$p_{n}(z;\lambda)$ and $\sigma_{n}(z;\lambda)$ solving (3) and (5),
respectively, can be expressed as

(10)
$\begin{cases}q_{n}(z;\lambda)&\displaystyle=\frac{d}{dz}\log\frac{\tau_{n-1}(z;\lambda)}{\tau_{n}(z;\lambda)},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
p_{n}(z;\lambda)&\displaystyle=-2\frac{d^{2}}{dz^{2}}\log\tau_{n}(z;\lambda),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\sigma_{n}(z;\lambda)&\displaystyle=\frac{d}{dz}\log\tau_{n}(z;\lambda).\end{cases}$

It was pointed out by Clarkson in [10], see also [16], based on the numerical
computations, that the functions $q_{n}(z;0)$ seemed to be tronquée solutions
of PII; that is, they have no poles in a sector of the complex plane.
Numerical computations of Fornberg and Weideman in [16] also suggested that
for general $\lambda$ there are three sectors in which the solutions
$q_{n}(z;\lambda)$ are pole free, see Figure 1. These conjectures were proven
for $z$ large by the third author in [12], where large $z$ expansions in the
pole free sectors were obtained. In this note we further provide scaling
limits of the purely Airy solutions. We shall obtain scaling limits of all
special solutions in a subsequent publication. Moreover, it was suggested by
Clarkson in [10], see also [11], again, based on numerical computations, that
$q_{n+1}(z;0)<q_{n}(z;0)$ for $z<0$. We confirm this conjecture for all $n$
large.

Figure 1. Zeros (blue) and poles (red) of $q_{3}(z;0)$ (left panel) and
$q_{3}(z;\infty)$ (right panel).

##  3\. Complex Cubic Ensemble of Random Matrices

Consider the unitary ensemble of random matrices with the cubic potential

$V(M)=-\frac{1}{3}M^{3}+tM,$

where $t$ is a complex parameter. The partition function of this ensemble is
formally defined by the matrix integral over the space of $N\times N$
Hermitian matrices:

$\int_{\mathcal{H}_{N}}e^{-N\mathrm{Tr}\left(-\frac{1}{3}M^{3}+tM\right)}dM.$

Then the formal partition functions of the eigenvalues becomes

$\int_{-\infty}^{\infty}\ldots\int_{-\infty}^{\infty}\prod_{1\leq j<k\leq
N}(s_{j}-s_{k})^{2}\,\prod_{j=1}^{N}e^{-NV(s_{j};t)}ds_{1}\ldots ds_{N},$

where $V(s;t)$ was defined above in (7). This expression is formal because the
integrals are divergent and need regularization. To achieve it, let

$\Gamma=\Gamma(\lambda):=\alpha_{0}L_{0}+\alpha_{1}L_{1}+\alpha_{2}L_{2},$

where the rays $L_{k}$ were defined right after (7) and
$\alpha_{0},\alpha_{1},\alpha_{2}\in\mathbb{C}$ are complex parameters such
that

(11)
$\alpha_{0}=\frac{\lambda}{\pi},\quad\alpha_{1}=-\frac{\lambda}{2\pi}+\frac{1}{2\pi\mathrm{i}},\quad\text{and}\quad\alpha_{2}=-\frac{\lambda}{2\pi}-\frac{1}{2\pi\mathrm{i}},$

when $|\lambda|<\infty$ and $\alpha_{0}=1/\pi$,
$\alpha_{1}=\alpha_{2}=-1/2\pi$ when $\lambda=\infty$. Let

$Z_{N}(t):=\int_{\Gamma}\ldots\int_{\Gamma}\prod_{1\leq j<k\leq
N}(s_{j}-s_{k})^{2}\,\prod_{j=1}^{N}e^{-NV(s_{j};t)}ds_{1}\ldots ds_{N},$

which is well defined (the integrals are convergent) for all values
$t\in\mathbb{C}$, where we understand $\int_{\Gamma}$ as
$\sum\alpha_{k}\int_{\Gamma_{k}}$. As shown in [5], the topological expansion
of the free energy

$F_{N}(t)=\frac{1}{N^{2}}\log Z_{N}(t)$

of the cubic ensemble of random matrices is connected to the enumeration of
regular graphs of degree 3 on Riemann surfaces. The partition function
$Z_{N}(t)$ can also be expressed as

(12) $Z_{N}(t)=N!D_{N-1}(t,N),\quad
D_{n}(t,N):=\left[\int_{\Gamma}s^{i+j}e^{-NV(s;t)}ds\right]_{i,j=0}^{n}.$

The entries of the matrix defining $D_{n}(-t,N)$ can be written as

$\displaystyle\int_{\Gamma}s^{k}e^{-NV(s;-t)}ds$
$\displaystyle=\frac{1}{N^{(k+1)/3}}\int_{\Gamma}s^{k}e^{-V\left(s;-N^{2/3}t\right)}ds$
$\displaystyle=\frac{1}{N^{k+1/3}}\frac{d^{k}}{dt^{k}}\left(\int_{\Gamma}e^{-V\left(s;-N^{2/3}t\right)}ds\right)$
$\displaystyle=\frac{1}{N^{k+1/3}}\frac{d^{k}}{dt^{k}}\left(\alpha_{0}\pi\,\mathrm{Bi}\left(-N^{2/3}t\right)+(\alpha_{1}-\alpha_{2})\pi\mathrm{i}\,\mathrm{Ai}\left(-N^{2/3}t\right)\right),$

where we used expressions (7) and (8) as well as the condition
$\alpha_{0}+\alpha_{1}+\alpha_{2}=0$. That is, upon setting
$z=(\sqrt{2}N)^{2/3}t$, we get that

$\int_{\Gamma}s^{k}e^{-NV(s;-t)}ds=\frac{2^{k/3}}{N^{(k+1)/3}}\frac{d^{k}}{dz^{k}}\left(\alpha_{0}\pi\,\mathrm{Bi}\left(-2^{-1/3}z\right)+(\alpha_{1}-\alpha_{2})\pi\mathrm{i}\,\mathrm{Ai}\left(-2^{-1/3}z\right)\right).$

This formula and the Hankel structure of the determinant $D_{n}(-t,N)$ then
imply that

(13)
$D_{n}(-t,N)=\frac{2^{n(n+1)/3}}{N^{(n+1)(n+2)/3}}\tau_{n+1}(z;\lambda),\quad
z=(\sqrt{2}N)^{2/3}t,$

where $\tau_{n}(z;\lambda)$ was defined in (9). In particular, if we take
$n=N$, the last two representations in (10) can be rewritten in terms of the
free energy $F_{N}(t)$ as

$\begin{cases}p_{N}(z;\lambda)&\displaystyle=-2N^{2}\frac{d^{2}}{dz^{2}}F_{N}\left(-(\sqrt{2}N)^{-2/3}z\right),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\sigma_{N}(z;\lambda)&\displaystyle=N^{2}\frac{d}{dz}F_{N}\left(-(\sqrt{2}N)^{-2/3}z\right).\end{cases}$

The cases of particular interest to us in this note are
$\lambda=0,-\mathrm{i},\mathrm{i}$. These cases are essentially the same as
they correspond to taking the seed function in (9) to be

$\mathrm{Ai}\left(-2^{-1/3}z\right),\quad\mathrm{Ai}\left(-2^{-1/3}ze^{2\pi\mathrm{i}/3}\right),\quad\text{and}\quad\mathrm{Ai}\left(-2^{-1/3}ze^{-2\pi\mathrm{i}/3}\right)$

since $\mathrm{Ai}(z)\pm\mathrm{i}\mathrm{Bi}(z)=2e^{\pm
2\pi\mathrm{i}/3}\mathrm{Ai}(ze^{\mp 2\pi\mathrm{i}/3})$, see [13, Equation
(9.2.1)]. It is worth pointing out that, in these special cases, the functions
$\tau_{n+1}$ also appear in the study of the scaling limits of the largest
eigenvalue in the GUE ensemble, see [17, Proposition 28].

The above representations can be further rewritten utilizing the connection to
orthogonal polynomials. Let $P_{n}(s;t,N)$ be a non-identically zero
polynomial of degree at most $n$ such that

(14) $\int_{\Gamma}s^{k}P_{n}(s;t,N)e^{-NV(s;t)}ds=0,\quad
k\in\\{0,\ldots,n-1\\}.$

Due to the non-Hermitian character of the above relations, it might happen
that a polynomial satisfying (14) is non-unique. In this case we understand by
$P_{n}(z;t,N)$ the monic non-identically zero polynomial of the smallest
degree (such polynomial is always unique). The standard determinantal
representation of orthogonal polynomials yields that

(15)
$h_{n}(t,N):=\int_{\Gamma}P_{n}^{2}(s;t,N)e^{-NV(s;t)}ds=\frac{D_{n}(t,N)}{D_{n-1}(t,N)},$

where $D_{-1}(t,N)\equiv 1$. Observe that $D_{n}(t,N)$ is an entire function
of $t$ and therefore each $h_{n}(t,N)$ is meromorphic in $\mathbb{C}$. Hence,
given $n$, the set of the values $t$ for which there exists
$k\in\\{0,\ldots,n\\}$ such that $h_{k}(t,N)=0$ is countable with no limit
points in the finite plane. Outside of this set the standard argument using
(14) shows that

(16)
$sP_{n}(s;t,N)=P_{n+1}(s;t,N)+\beta_{n}(t,N)P_{n}(s;t,N)+\gamma_{n}^{2}(t,N)P_{n-1}(s;t,N),$

and by analytic continuation this relation extends to those values of $t$ for
which we have that $D_{n-1}(t,N)D_{n}(t,N)\neq 0$ (that is, $n+1$-st and
$n$-th polynomials appearing in (16) have the prescribed degrees), where

(17) $\gamma_{n}^{2}(t,N)=\frac{h_{n}(t,N)}{h_{n-1}(t,N)}.$

Finally, denote by $p_{n,n-1}(t,N)$ the coefficient of $P_{n}(s;t,N)$ next to
$s^{n-1}$, where $p_{0,-1}(t,N):=0$. It easily follows from (16) that

(18) $\beta_{n}(t,N)=p_{n,n-1}(t,N)-p_{n+1,n}(t,N).$

In what follows, when it is important to us to stress that the quantities
appearing in (12), (14), (15), (17), and (18) depend on $\lambda$, we shall
use superscript $(\lambda)$.

######  Theorem 3.1.

Fix $N\geq 1$. Given $\lambda\in\overline{\mathbb{C}}$, it holds for each
$n\geq 1$ that

(19)
$\begin{cases}q_{n}(z;\lambda)&\displaystyle=-(N/2)^{1/3}\,\beta_{n-1}^{(\lambda)}\left(-(\sqrt{2}N)^{-2/3}z,N\right),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
p_{n}(z;\lambda)&\displaystyle=-2(N/2)^{2/3}\,\gamma_{n}^{(\lambda)}\left(-(\sqrt{2}N)^{-2/3}z,N\right)^{2},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\sigma_{n}(z;\lambda)&\displaystyle=-(N/2)^{1/3}\,p_{n,n-1}^{(\lambda)}\left(-(\sqrt{2}N)^{-2/3}z,N\right).\end{cases}$

###### Proof.

Since $\lambda,N$ are fixed, we do not indicate the dependence on them of all
the quantities related to orthogonal polynomials. A combination of (10), (13),
and (15) immediately yields that

$q_{n}(z;\lambda)=-\frac{d}{dz}\log h_{n-1}\left(-(\sqrt{2}N)^{-2/3}z\right).$

Since $P_{n}(s;t)$ is a monic polynomial, its partial $t$ derivative is a
polynomial of degree at most $n-1$, which is, of course, orthogonal to
$P_{n}(s;t)$. Hence, we get from (15), (14), and (18) that

$\displaystyle\frac{d}{dt}h_{n}(-t)$
$\displaystyle=N\int_{\Gamma}sP_{n}^{2}(s;-t)e^{-NV(s;-t)}ds$
$\displaystyle=N\int_{\Gamma}\big{(}s^{n+1}+p_{n,n-1}(-t)s^{n}\big{)}P_{n}(s;-t)e^{-NV(s;-t)}ds$
$\displaystyle=N\int_{\Gamma}P_{n+1}(s;-t)P_{n}(s;-t)e^{-NV(s;-t)}ds+N\beta_{n}(-t)h_{n}(-t)$
(20) $\displaystyle=N\beta_{n}(-t)h_{n}(-t),$

from which the formula for $q_{n}(z;\lambda)$ follows. Now, we get from (15)
that

$\log D_{n}(t)=\sum_{k=0}^{n}\log h_{k}(t).$

The above expression and (20), together with (10) and (13), establish formula
for $\sigma_{n}(z;\lambda)$ since

(21) $\frac{d}{dt}\log
D_{n-1}(-t)=N\sum_{k=0}^{n-1}\beta_{k}(-t)=-Np_{n,n-1}(-t).$

Finally, we get from orthogonality relations (14) that

$\displaystyle 0$
$\displaystyle=\frac{d}{dt}\left(\int_{\Gamma}P_{n}(s;-t)P_{n-1}(s;-t)e^{-NV(s;-t)}ds\right)$
$\displaystyle=\int_{\Gamma}\left(NsP_{n}(s;-t)-p_{n,n-1}^{\prime}(-t)s^{n-1}\right)P_{n-1}(s;-t)e^{-NV(s;-t)}ds$
(22) $\displaystyle=Nh_{n}(-t)-p_{n,n-1}^{\prime}(-t)h_{n-1}(-t),$

which, together with (10), (13), and (17), gives the expression for
$p_{n}(z;\lambda)$. ∎

The results of the previous theorem can equivalently be rewritten in the
following form.

######  Corollary 3.2.

Given $\lambda\in\overline{\mathbb{C}}$, it holds for each $n\geq 1$ that

(23)
$\begin{cases}q_{n}(z;\lambda)&\displaystyle=-2^{-1/3}\,\beta_{n-1}^{(\lambda)}\big{(}-2^{-1/3}z,1\big{)},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
p_{n}(z;\lambda)&\displaystyle=-2^{1/3}\,\gamma_{n}^{(\lambda)}\big{(}-2^{-1/3}z,1\big{)}^{2},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\sigma_{n}(z;\lambda)&\displaystyle=-2^{-1/3}\,p_{n,n-1}^{(\lambda)}\big{(}-2^{-1/3}z,1\big{)}.\end{cases}$

###### Proof.

Again, we do not indicate the dependence on $\lambda$. It readily follows from
(14) that

$0=\int_{\Gamma}s^{k}P_{n}(s;t,N)e^{-NV(s;t)}ds=N^{-(k+1)/3}\int_{\Gamma}s^{k}P_{n}\left(N^{-1/3}s;t,N\right)e^{-V(s;N^{2/3}t)}ds.$

Therefore, it holds that

$P_{n}(s;t,N)=N^{-n/3}P_{n}\big{(}N^{1/3}s;N^{2/3}t,1\big{)}$

when $\deg P_{n}=n$, which happens for the values of $t$ outside of a
countable set without limit points in the finite plane. From this we can
readily deduce that

(24)
$\begin{cases}\beta_{n-1}(t;N)&=N^{-1/3}\beta_{n-1}\left(N^{2/3}t,1\right),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\gamma_{n}^{2}(t;N)&=N^{-2/3}\gamma_{n}^{2}\left(N^{2/3}t,1\right),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
p_{n,n-1}(t,N)&=N^{-1/3}p_{n,n-1}\left(N^{2/3}t,1\right).\end{cases}$

The above relations and (19) easily yield (23). ∎

##  4\. Scaling Limits

Now we can formulate a result on the asymptotic behavior of the Airy solutions
when $\lambda\in\\{0,-\mathrm{i},\mathrm{i}\\}$. To this end, it was shown in
[7, Section 5] that the critical graph of a quadratic differential

$-(1+1/s)^{3}ds^{2},$

say $\mathcal{C}$, consists of 5 critical trajectories emanating from $-1$ at
the angles $2\pi k/5$, $k\in\\{0,1,2,3,4\\}$, one of them being $(-1,0)$,
other two forming a loop crossing the real line approximately at $0.635$, and
the last two approaching infinity along the imaginary axis without changing
the half-plane (upper or lower), see Figure 2 (left panel).

$0$$-1$

${\footnotesize-\sqrt[3]{1/2}}$$\Omega_{(0)}$

$3\cdot 2^{-2/3}e^{-2\pi\mathrm{i}/3}$$3\cdot
2^{-2/3}e^{2\pi\mathrm{i}/3}$$O_{(0)}$

Figure 2. Left panel: the critical graph $\mathcal{C}$; middle panel: the set
$\Delta$ (solid lines) and the domain $\Omega_{(0)}$ (shaded region); right
panel: the domain $O_{(0)}$ (shaded region).

Define

$\Delta:=\big{\\{}x:~{}2x^{3}\in\mathcal{C}\big{\\}}$

and put $\Omega_{(0)}$ to be the subset of the right-half plane bounded by
three smooth subarcs111These subarcs are the one originating at
$e^{\pi\mathrm{i}/3}\sqrt[3]{1/2}$ and extending to infinity in the direction
of the angle $\pi/2$, the one connecting $e^{\pi\mathrm{i}/3}\sqrt[3]{1/2}$
and $e^{-\pi\mathrm{i}/3}\sqrt[3]{1/2}$, and the one originating at
$e^{-\pi\mathrm{i}/3}\sqrt[3]{1/2}$ and extending to infinity in the direction
of the angle $-\pi/2$. of $\Delta$ as on Figure 2 (middle panel). Further put
$\Omega_{(\pm\mathrm{i})}:=e^{\mp 2\pi\mathrm{i}/3}\Omega_{(0)}$. Let

$t(x):=(x^{3}-1)/x.$

The function $t(x)$ is holomorphic in each $\Omega_{(\lambda)}$,
$\lambda\in\\{0,-\mathrm{i},\mathrm{i}\\}$, with non-vanishing derivative
there. Set

$O_{(\lambda)}:=t(\Omega_{(\lambda)}),$

see Figure 2 (right panel). The inverse map $x_{\lambda}(t)$ exists and is
holomorphic in $O_{(\lambda)}$. One can readily check that
$O_{(\pm\mathrm{i})}=e^{\pm 2\pi\mathrm{i}/3}O_{(0)}$ and

(25) $x_{0}(t)=e^{\pm 2\pi\mathrm{i}/3}x_{\pm\mathrm{i}}\left(te^{\pm
2\pi\mathrm{i}/3}\right).$

Clearly, the branch $x_{0}(t)$ is positive on $O_{(0)}\cap(-\infty,\infty)$
(it was shown in [20, Theorem 1.1] that $\partial O_{(0)}$ intersects the real
line at $t_{0}\sim-1.0005424$) and is equal to $1$ at the origin.

######  Theorem 4.1.

Let $\lambda\in\\{0,-\mathrm{i},\mathrm{i}\\}$. For all $n$ large enough it
holds that

$\begin{cases}\displaystyle(2/n)^{1/3}q_{n}\left(-(\sqrt{2}n)^{2/3}t;\lambda\right)&=-x_{\lambda}(t)+\mathcal{O}\big{(}n^{-1}\big{)},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\displaystyle(2/n)^{2/3}p_{n}\left(-(\sqrt{2}n)^{2/3}t;\lambda\right)&=1/x_{\lambda}(t)+\mathcal{O}\big{(}n^{-2}\big{)},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ \displaystyle
2^{1/3}n^{-4/3}\sigma_{n}\left(-(\sqrt{2}n)^{2/3}t;\lambda\right)&=x_{\lambda}(t)-(2x_{\lambda}(t))^{-2}+\mathcal{O}\big{(}n^{-2}\big{)},\end{cases}$

uniformly on closed subset of $\\{t\in
O_{(\lambda)}:|x_{\lambda}(t)|<\sqrt{n}/\delta\\}$ in the first formula and of
$\\{t\in
O_{(\lambda)}:|\mathrm{arg}(t-t_{\lambda})-\theta_{\lambda}|>\delta\\}$ in the
other two formulae, where $\delta>0$ is any, $\theta_{0}=\pi$,
$\theta_{\pm\mathrm{i}}=\mp\pi/3$, and $t_{\lambda}=e^{2\pi\lambda/3}t_{0}$.

###### Proof.

It holds for $\lambda\in\\{-\mathrm{i},\mathrm{i}\\}$ that

$P_{n}^{(0)}(s;t,N)=e^{2\pi\lambda
n/3}P_{n}^{(\lambda)}\left(se^{-2\pi\lambda/3};te^{2\pi\lambda/3},N\right),$

see (14). This implies that

$\begin{cases}\gamma_{n}^{(0)}(t,N)^{2}&=e^{4\pi\lambda/3}\gamma_{n}^{(\lambda)}\big{(}te^{2\pi\lambda/3},N)^{2},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\beta_{n}^{(0)}(t,N)&=e^{2\pi\lambda/3}\beta_{n}^{(\lambda)}\big{(}te^{2\pi\lambda/3},N),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
p_{n,n-1}^{(0)}(t,N)&=e^{2\pi\lambda/3}p_{n,n-1}^{(\lambda)}\big{(}te^{2\pi\lambda/3},N).\end{cases}$

It has been shown in [7, Theorem 4.5] that for $|n-N|\leq N_{0}$ it holds that

(26)
$\begin{cases}\gamma_{n}^{(-\mathrm{i})}(t,N)^{2}&\sim-1/(2x_{-\mathrm{i}}(t))+\sum_{k=1}^{\infty}G_{k}(t;n-N)N^{-k},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ \beta_{n}^{(-\mathrm{i})}(t,N)&\sim
x_{-\mathrm{i}}(t)+\sum_{k=1}^{\infty}B_{k}(t;n-N)N^{-k},\end{cases}$

where the expansions hold uniformly on compact subsets of $O_{(-\mathrm{i})}$
and closed subsets of

$O_{(-\mathrm{i}),\delta}:=\\{t\in
O_{(-\mathrm{i})}:|\mathrm{arg}(t-t_{-\mathrm{i}})-\theta_{-\mathrm{i}}|>\delta\\}$

in the case of $\gamma_{N}^{(-\mathrm{i})}(t,N)$ (the actual condition of
separation from $\partial O_{(-\mathrm{i})}$ is more refined in [7, Theorem
4.5], see [7, Definition 4.1]). The functions $G_{k}(t;n-N)$ and
$B_{k}(t;n-N)$ are analytic in $O_{(-\mathrm{i})}$ and $G_{2k-1}(t;0)\equiv
0$. It is also implicitly contained in [7, Section 8.2] that

(27)
$p_{n,n-1}^{(-\mathrm{i})}(t,N)=-nx_{-\mathrm{i}}(t)+N(2x_{-\mathrm{i}}(t))^{-2}+\mathcal{O}\big{(}N^{-1}\big{)},$

where the error term is again analytic in $O_{(-\mathrm{i})}$, uniform on
compact subsets of $O_{(-\mathrm{i})}$ and closed subsets of
$O_{(-\mathrm{i}),\delta}$ when $n=N$. The second and third claims of the
theorem now follow from Theorem 3.1 by taking $N=n$ in (27) and in the first
formula of (26). The first claim follows by replacing $n$ with $n-1$ in the
second line of (26) and taking $N=n$ there, except this proves it only on
compact subsets of $O_{(\lambda)}$. To overcome this difficulty, write

(28) $\gamma_{n}^{2}(t,n)=-1/(2x_{\lambda}(t))+E_{\gamma_{n}}(t),$

where we suppress the dependence on $\lambda$ of the recurrence coefficients
as emphasizing it is no longer important. The error term $E_{\gamma_{n}}(t)$
is such that

$|E_{\gamma_{n}}(t)|\leq C_{F}/n^{2},\quad t\in F,$

for each closed subset $F\subset O_{(\lambda),\delta}$ and any $\delta>0$.
Notice that $F_{\epsilon}=\\{t:\mathrm{dist}(t,F)\leq\epsilon\\}$ belongs to
$O_{(\lambda),\delta^{\prime}}$ for some appropriate choices of
$\epsilon,\delta^{\prime}>0$. Applying Cauchy integral formula on circles of
radius $\epsilon$ around every point of $F$ yields that
$E_{\gamma_{n}}^{\prime}(t)$ satisfies the same type of bound as
$E_{\gamma_{n}}(t)$, possibly with a different constant. Now, let us write

(29) $\beta_{n-1}(t,n)=x_{\lambda}(t)+E_{\beta_{n-1}}(t).$

It follows from Theorem 3.1 and the second relation in (1) that

$2n\beta_{n-1}(t,n)\gamma_{n}^{2}(t,n)+n=(\gamma_{n}^{2})^{\prime}(t,n).$

Therefore, we get from (28) and (29) that

$E_{\beta_{n-1}}(t)\big{(}1-2x_{\lambda}(t)E_{\gamma_{n}}(t)\big{)}=2x_{\lambda}^{2}(t)E_{\gamma_{n}}(t)-\frac{1}{2n}\left(\frac{x^{\prime}_{\lambda}(t)}{x_{\lambda}(t)}+2x_{\lambda}(t)E_{\gamma_{n}}^{\prime}(t)\right).$

Because $x^{3}-xt-1=0$, we have that $x^{\prime}=x^{2}/(2x^{3}+1)$ for any
branch $x(t)$. Since we already have the estimate of $E_{\beta_{n-1}}(t)$ on
compact subsets of $O_{(\lambda)}$, we can assume without loss of generality
that $|x_{\lambda}(t)|\geq 1$. Thus, it holds that

$E_{\beta_{n-1}}(t)\big{(}1-\mathcal{O}_{\delta}\big{(}n^{-3/2}\big{)}\big{)}=\mathcal{O}_{\delta}\big{(}n^{-1}\big{)}$

uniformly on closed subsets $O_{(\lambda)}$ for which
$1\leq|x_{\lambda}(t)|\leq\sqrt{n}/\delta$. This, of course, finishes the
proof of the theorem. ∎

Theorem 4.1 yields the following corollary, which gives pole and zero free
regions around the origin in the $t$ plane for the Painlevé II functions. The
results in [12] in this direction are asymptotic in $z$ without uniformity in
the parameter $n$, but in our case we can show that the radius of these
domains increases with $n$, something that is consistent with previous
numerical experiments.

######  Corollary 4.2.

For each $t_{*}>0$ such that $\\{|t|\leq t_{*}\\}\subset O_{(\lambda)}$ there
exists $n_{*}$ such that the functions $p_{n}(z;\lambda)$, $q_{n}(z;\lambda)$
and $\sigma_{n}(z;\lambda)$ are analytic and non-vanishing in
$\big{\\{}|z|\leq(\sqrt{2}n)^{2/3}t_{*}\big{\\}}$ for all $n\geq n_{*}$.

The asymptotic behavior of the above quantities in the complement of
$\overline{O}_{(-\mathrm{i})}$ was studied in [3]. However, the results of [3,
Theorem 4.2] on asymptotic behavior of $\beta_{n}$ and $\gamma_{n}^{2}$ are
not complete, not in the least due to the observed presence of their poles in
this region.

##  5\. Monotonicity of Airy Solutions

It was conjectured by Clarkson in [10], see also [11], that

(30) $q_{n+1}(z;0)<q_{n}(z;0)$

is true for every $z\leq 0$ and every $n\geq 1$. In light of Corollary 3.2,
this claim is equivalent to $\beta_{n}(t,1)>\beta_{n-1}(t,1)$ being true for
$t\geq 0$ and $n\geq 1$. Since [11] uses slightly different notation, let us
point out that

$D_{n-1}(t,1)=2^{n(n-1)/3}\mathrm{det}\left[\frac{d^{j+k}}{dt^{j+k}}\mathrm{Ai}(t)\right]_{j,k=0}^{n-1},$

where everything related to orthogonal polynomials corresponds to $\lambda=0$
in (11), see (13) and (9) ($C_{1}=1$ and $C_{2}=0$). This quantity was labeled
by $\Delta_{n}(t)$ in [11, Equation (5.2)]. Hence, it follows from (21), (22),
and (17) as well as (20) and (15) that

(31) $\gamma_{n}^{2}(t,1)=\frac{d^{2}}{dt^{2}}\log
D_{n-1}(t,1)\quad\text{and}\quad\beta_{n}(t,1)=\frac{d}{dt}\log\frac{D_{n-1}(t,1)}{D_{n}(t,1)}.$

These quantities were labeled as $-a_{n}(t)$ and $b_{n-1}(t)$ in [11, Equation
(5.1)]. In [11], these functions were studied since they are solutions of an
_alternative discrete Painlevé I system_ :

$\begin{cases}\gamma_{n}^{2}(t,1)+\gamma_{n+1}^{2}(t,1)+\beta_{n+1}^{2}(t,1)=t\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\gamma_{n}^{2}(t,1)\big{(}\beta_{n+1}(t,1)+\beta_{n}(t,1)\big{)}=-n,\end{cases}$

see [15]. Using rescaling formulae (24), the above relations can be rewritten
as

$\begin{cases}\gamma_{n}^{2}(t,N)+\gamma_{n+1}^{2}(t,N)+\beta_{n+1}^{2}(t,N)=t\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\gamma_{n}^{2}(t,N)\big{(}\beta_{n+1}(t,N)+\beta_{n}(t,N)\big{)}=-n/N,\end{cases}$

which are also known as _discrete string equations_ , see [8]. It was
conjectured in [11, Conjectures 3.2 and 5.3] that

(32)
$0<\beta_{n}^{\prime}(t,1)<1/(2\sqrt{t})\quad\text{and}\quad\sqrt{t}<\beta_{n}(t,1)<\beta_{n+1}(t,1)$

for all $t>0$ and $n\geq 0$. It was also shown that the first inequality in
(32) implies the last inequality there and therefore (30). To see this,
observe that it follows from the difference of the successive top string
equations, (31), and the bottom string equation that

$\beta_{n+1}(t,1)-\beta_{n}(t,1)=\frac{\gamma_{n-1}^{2}(t,1)-\gamma_{n+1}^{2}(t,1)}{\beta_{n+1}(t,1)+\beta_{n}(t,1)}=-(\beta_{n}^{\prime}(t,1)+\beta_{n-1}^{\prime}(t,1))\frac{\gamma_{n}^{2}(t,1)}{n}.$

Now, it was shown in [19, Equation (4.17)] that

$D_{n+1}(t,1)=-(n+1)D_{n-1}(t,1)\int_{t}^{\infty}\frac{D_{n}^{2}(s,1)}{D_{n-1}^{2}(s,1)}ds$

for $n\geq 0$, where $D_{0}(t,1)=\mathrm{Ai}(t)$ and $D_{-1}(t,1)\equiv 1$.
Let $\iota_{1}<0$ be the largest real zero of $\mathrm{Ai}(t)$. Then it
follows from the above relation and (15) that
$\mathrm{sgn}\,h_{n}(t,1)=(-1)^{n}$ for $t\in(\iota_{1},\infty)$. This, of
course, means that $\gamma_{n}^{2}(t,1)<0$ there by (17), which yields the
desired implication on $(\iota_{1},\infty)$.

######  Theorem 5.1.

There exists $n_{0}>0$ such that for all $n\geq n_{0}$ it holds that

$\beta_{n-1}^{\prime}(t,1)>0,\quad t\in[\iota_{1},\infty).$

In particular, (30) holds for $z\leq-2^{1/3}\iota_{1}$ and all $n\geq n_{0}$.

###### Proof.

It was shown in [11, Lemma 3.3] that $\beta_{n}^{\prime}(t,1)>0$ for $t\geq
2^{-2/3}n^{4/3}$. Thus, we are only interested in
$t\in[\iota_{1},2^{-2/3}n^{4/3}]$. It follows from (24) that

$\beta_{n-1}(t,1)=n^{1/3}\beta_{n-1}(n^{-2/3}t,n).$

Hence, it is sufficient for us to show that $\beta_{n}^{\prime}(t,n)>0$ for
$t\in F_{n}:=[-\epsilon,n^{2/3}]$, where $\epsilon>0$ and small. To this end,
recall that $x^{\prime}=x^{2}/(2x^{3}+1)$. Observe also that $x_{0}(t)>0$ for
$t\geq-\epsilon$ (this is the branch that is equal to $1$ at $0$ and is
positive for $t\geq 0$). Hence, $x_{0}^{\prime}(t)>0$ for $t\geq-\epsilon$ and
therefore $x_{0}(t)$ is increasing there. Notice that $x_{0}(t)<2\sqrt{t}$ for
$t\geq 1$ (constant $2$ is by no means optimal). Indeed, otherwise we would
have that

$1=x_{0}(t)\big{(}x_{0}^{2}(t)-t\big{)}>6t^{3/2}\geq 6,$

which is impossible. Thus, $F_{n}$ lies within the set where $|x_{0}(t)|\leq
2n^{1/3}<2n^{1/2}$. Moreover, there exists an open set $U_{n}\supset F_{n}$
satisfying $\mathrm{dist}(\partial U_{n},F_{n})\geq c>0$ that also lies within
$\\{t:|x_{0}(t)|\leq 2n^{1/2}\\}$. As in the proof of Theorem 4.1, we can
write

$\beta_{n-1}(t,n)=x_{0}(t)+E_{\beta_{n-1}}(t),$

where $E_{\beta_{n-1}}(t)$ is analytic in $O_{(0)}$ and satisfies (through the
use of Cauchy integral formula)

$|E_{\beta_{n-1}}^{\prime}(t)|\leq K/n,\quad t\in F_{n},$

for some constant $K>0$. Since we can take $\epsilon$ small enough so that
$2x_{0}^{3}(-\epsilon)>1$, it then holds that

$\beta_{n-1}^{\prime}(t,n)=x_{0}^{\prime}(t)+E_{\beta_{n-1}}^{\prime}(t)>\frac{1}{4x_{0}(t)}-\frac{K}{n}\geq\frac{1}{4x_{0}(n^{2/3})}-\frac{K}{n}>\frac{n^{2/3}-8K}{8n}$

for $t\in F_{n}$, which finishes the proof of the theorem. ∎

## References

  * [1] Gabriel Álvarez, Luis Martínez Alonso, and Elena Medina.  Determination of $S$-curves with applications to the theory of non-Hermitian orthogonal polynomials.  J. Stat. Mech. Theory Exp., (6):P06006, 28, 2013.
  * [2] Gabriel Álvarez, Luis Martínez Alonso, and Elena Medina.  Phase structure and asymptotic zero densities of orthogonal polynomials in the cubic model.  J. Comput. Appl. Math., 284:10–25, 2015.
  * [3] Ahmad Barhoumi, Pavel Bleher, Alfredo Deaño, and Maxim Yattselev.  Investigation of the two-cut phase region in the complex cubic ensemble of random matrices.  J. Math. Phys., 63(6):Paper No. 063303, 40, 2022.
  * [4] D. Bessis, C. Itzykson, and J. B. Zuber.  Quantum field theory techniques in graphical enumeration.  Adv. in Appl. Math., 1(2):109–157, 1980.
  * [5] Pavel Bleher and Alfredo Deaño.  Topological expansion in the cubic random matrix model.  Int. Math. Res. Not. IMRN, (12):2699–2755, 2013.
  * [6] Pavel Bleher and Alfredo Deaño.  Painlevé I double scaling limit in the cubic random matrix model.  Random Matrices Theory Appl., 5(2):1650004, 58, 2016.
  * [7] Pavel Bleher, Alfredo Deaño, and Maxim Yattselev.  Topological expansion in the complex cubic log-gas model: one-cut case.  J. Stat. Phys., 166(3-4):784–827, 2017.
  * [8] Pavel Bleher and Karl Liechty.  Random matrices and the six-vertex model, volume 32 of CRM Monograph Series.  American Mathematical Society, Providence, RI, 2014.
  * [9] E. Brézin, C. Itzykson, G. Parisi, and J. B. Zuber.  Planar diagrams.  Comm. Math. Phys., 59(1):35–51, 1978.
  * [10] Peter A. Clarkson.  On Airy solutions of the second Painlevé equation.  Stud. Appl. Math., 137(1):93–109, 2016.
  * [11] Peter A. Clarkson, Ana F. Loureiro, and Walter Van Assche.  Unique positive solution for an alternative discrete Painlevé I equation.  J. Difference Equ. Appl., 22(5):656–675, 2016.
  * [12] Alfredo Deaño.  Large $z$ asymptotics for special function solutions of Painlevé II in the complex plane.  SIGMA Symmetry Integrability Geom. Methods Appl., 14:Paper No. 107, 19, 2018.
  * [13] F.W.J Olver et al. editors.  NIST digital library of mathematical functions.  http://dlmf.nist.gov.
  * [14] Hermann Flaschka and Alan C. Newell.  Monodromy- and spectrum-preserving deformations. I.  Comm. Math. Phys., 76(1):65–116, 1980.
  * [15] A. S. Fokas, B. Grammaticos, and A. Ramani.  From continuous to discrete Painlevé equations.  J. Math. Anal. Appl., 180(2):342–360, 1993.
  * [16] Bengt Fornberg and J. A. C. Weideman.  A computational exploration of the second Painlevé equation.  Found. Comput. Math., 14(5):985–1016, 2014.
  * [17] P. J. Forrester and N. S. Witte.  Application of the $\tau$-function theory of Painlevé equations to random matrices: PIV, PII and the GUE.  Comm. Math. Phys., 219(2):357–398, 2001.
  * [18] B. Gambier.  Sur les équations différentielles du second ordre et du premier degré dont l’intégrale générale est a points critiques fixes.  Acta Math., 33(1):1–55, 1910.
  * [19] Daan Huybrechs, Arno Kuijlaars, and Nele Lejon.  A numerical method for oscillatory integrals with coalescing saddle points.  SIAM J. Numer. Anal., 57(6):2707–2729, 2019.
  * [20] Daan Huybrechs, Arno B. J. Kuijlaars, and Nele Lejon.  Zero distribution of complex orthogonal polynomials with respect to exponential weights.  J. Approx. Theory, 184:28–54, 2014.
  * [21] Alexander R. Its, Arno B. J. Kuijlaars, and Jörgen Östensson.  Critical edge behavior in unitary random matrix ensembles and the thirty-fourth Painlevé transcendent.  Int. Math. Res. Not. IMRN, (9):Art. ID rnn017, 67, 2008.
  * [22] Michio Jimbo and Tetsuji Miwa.  Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. II.  Phys. D, 2(3):407–448, 1981.
  * [23] Kazuo Okamoto.  Polynomial Hamiltonians associated with Painlevé equations. I.  Proc. Japan Acad. Ser. A Math. Sci., 56(6):264–268, 1980.
  * [24] Kazuo Okamoto.  Polynomial Hamiltonians associated with Painlevé equations. II. Differential equations satisfied by polynomial Hamiltonians.  Proc. Japan Acad. Ser. A Math. Sci., 56(8):367–371, 1980.
  * [25] Kazuo Okamoto.  Studies on the Painlevé equations. III. Second and fourth Painlevé equations, $P_{{\rm II}}$ and $P_{{\rm IV}}$.  Math. Ann., 275(2):221–255, 1986.


# Probabilistic Differentiable Filters Enable Ubiquitous Robot Control with
Smartwatches

Fabian C Weigend, Xiao Liu, Heni Ben Amor All Authors are with the School of
Computing and Augmented Intelligence, Arizona State University $\\{$fweigend,
xliu330<EMAIL_ADDRESS>
###### Abstract

Ubiquitous robot control and human-robot collaboration using smart devices
poses a challenging problem primarily due to strict accuracy requirements and
sparse information. This paper presents a novel approach that incorporates a
probabilistic differentiable filter, specifically the Differentiable Ensemble
Kalman Filter (DEnKF), to facilitate robot control solely using Inertial
Measurement Units (IMUs) from a smartwatch and a smartphone. The implemented
system is cost-effective and achieves accurate estimation of the human pose
state. Experiment results from human-robot handover tasks underscore that
smart devices allow versatile and ubiquitous robot control. The code for this
paper is available at github.com/ir-lab/DEnKF and github.com/wearable-motion-
capture.

##  I Introduction

Examining the human-robot relationship is a central concern in the field of
artificial intelligence and robotics. To facilitate reliable human-robot
collaboration, accurate estimations of the state of humans and robots are
crucial. One challenge is that numerous robotics systems still rely heavily on
costly motion capture systems in estimating the human pose, thereby
restricting their suitability primarily to stationary setups or complex
calibration procedures.

The gold standard for motion capture are cameras [1]. These systems may
feature multiple cameras on multiple base stations or in one device, e.g,
Microsoft Kinect v2[2] or Intel RealSense. These systems usually require a
dedicated setup and suffer from line-of-sight issues. Alternatively, non-
optical systems enable motion capture through Inertial Measurement Units
(IMUs) [3, 4]. For this purpose, the Xsens motion capture system [5] is
commonly used [6]. Also, Sony’s recent Mocopi [7] promises new opportunities
for non-visual real-time human motion capture. A downside of most IMU-based
tracking is that they require users to carefully place multiple specialized
units on their body and a calibration procedure.

Recent research investigates the opportunities of omnipresent IMUs in
smartphones and smartwatches for motion capture. Some even limit themselves to
a single device [8, 9]. These approaches have promising potential for human-
robot collaboration because motion capture through smart devices is ubiquitous
and natural to the user [10]. However, leveraging the data of a single
smartwatch for motion capture of sufficient accuracy for robot control is
challenging, and previous work had to constrain the user to a constant body
forward-facing direction [10].

This work builds upon [10] and incorporates the sensor data of a connected
smartphone to enable ubiquitous robot control without body orientation
constraints. Recursive Bayesian filters, particularly Kalman filters, play a
pivotal role in tasks such as predicting the future movements of human
interaction partners [11], tracking objects over time [12], and ensuring
stability during robot locomotion [13]. We propose that advances in state
estimation, specifically the Differentiable Ensemble Kalman Filter (DEnKF)
[14], allow to enhance the accuracy and stability of motion capture using
ubiquitous smart devices.

As depicted in Figure 1, using a probabilistic differentiable filter provides
us with a distribution of solutions, aiding stability and adding an important
measure of uncertainty for robot control. The proposed filter allows to strike
the balance between less-constrained movements while still achieving stable
and effective pose estimations suitable for human-robot collaboration.

Figure 1: Top: Our Differentiable Ensemble Kalman Filter (DEnKF) achieves
robust body orientation and arm pose estimations from the sensor data of a
single smartwatch and a smartphone. The user can place the phone in any
pocket. Bottom: Orientation predictions are also accurate when the user sits.
We evaluate pose predictions on test data and in a human-robot collaboration
handover task.

##  II Methodology

We present our approach by describing our data collection and defining our
states and observations in Section II-A. Then, we define the used
Differentiable Ensemble Kalman Filter (DEKnF) for human pose state estimations
in Section II-B.

###  II-A Data Collection, Observation and State

For data collection, we develop two apps and make use of a research-grade
motion capture system. The apps for the smartwatch and smartphone stream
sensor data to a remote machine.

Observation: We define that the raw observation $\bf{y}$ consists of the
following values ${\bf y}=[\Delta
t,\boldsymbol{\theta}_{\mathrm{sw}},\boldsymbol{v},\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\phi},\rho,\mathbf{r}_{h}]^{\top}$,
with $\bf{y}\in\mathbb{R}^{22}$, where
$\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\phi}\in\mathbb{R}^{3}$
are the IMU readings. Namely, these are the average values of the linear
acceleration, gravity, and gyroscope since the previous observation $\Delta t$
seconds ago. The averaging is necessary because the apps stream data at around
80 Hz and the accelerometer and gyroscope record measurements at a faster
rate. We also integrate linear acceleration measurements between two
observations and denote them as velocities $\boldsymbol{v}$. In addition, the
apps record the virtual rotation vector sensor ($\bm{\theta}$), which is
provided by Android and Wear OS. We record this orientation in form of a
continuous six-dimensional rotation representation (6DRR), which is well-
suited for training neural networks [15].

Calibration: When the user starts the app on their smartwatch, they usually
hold their arm up, parallel to the chest and hip. We use this start position
to calibrate $\bm{\theta}$ with the first initial orientations
$\lx@glossaries@gls@link{main}{swrot}{\leavevmode\bm{\theta}}_{\text{init}}$,
which gives us the calibrated orientation
$\lx@glossaries@gls@link{main}{swrot}{\leavevmode\bm{\theta}}_{\text{sw}}$.
The value $\rho$ with $\rho\in\mathbb{R}$ is the atmospheric pressure sensor.
We also calibrate this one by subtracting the reading when the user started
the app from all subsequent ones. Finally, the phone provides the body-forward
facing direction. Also here, we record the orientation sensor of the phone and
calibrate it with the first observation in the starting position. Further, we
assume that the user has their forward-facing direction parallel to the watch
at the starting position. This allows us to estimate an offset rotation from
the phone to the forward direction, bringing it into the same global reference
frame as the watch. From this, we denote the body orientation as the sine and
cosine of the calibrated up-axis rotation of the phone, represented as
$\bf{r}_{h}\in\mathbb{R}^{2}$ in our raw observation $\bf{y}$.

Figure 2: Left: Our data collection setup with smart devices and Optitrack
system. Right: The DEnKF model structure. The sensor model projects raw
observations to the observation space, the stochastic transition model
forwards the ensemble one step in $t$, and the KF update step corrects the
state.

State: Our state entails the arm pose and body forward-facing direction of the
human. The ground truth values were recorded with the research-grade optical
motion capture system OptiTrack [1]. We record data from participants who wore
a 25-marker-upper-body suit along with the smartwatch on their left wrist and
a smartphone in their pocket. We collect the upper arm rotation
($\textbf{q}_{\mathrm{u}}$), lower arm rotation ($\textbf{q}_{\mathrm{l}}$)
also in the continuous 6DRR and the body-forward facing direction
$\bf{r}_{h}$. Therefore, the ground truth state for motion capture is denoted
as
$\mathbf{x}=[\text{\lx@glossaries@gls@link{main}{rot_larm}{\leavevmode$\textbf{q}_{\mathrm{l}}$}},\text{\lx@glossaries@gls@link{main}{rot_uarm}{\leavevmode$\textbf{q}_{\mathrm{u}}$}},\mathbf{r}_{h}]$,
where $\mathbf{x}\in\mathbb{R}^{14}$. In our Bayesian filtering framework, we
define the learned observation $\tilde{\mathbf{y}}$ to be equivalent to
$\mathbf{x}$.

###  II-B Differentiable Ensemble Kalman Filter

To model the temporal transition of the human pose $\mathbf{x}$ and map raw
observations $\mathbf{y}$ to the state space, we utilize DEnKF [14, 16, 17] as
shown in Figure 2. This state estimation approach enables us to learn and
infer the dynamics of the human pose over time, while efficiently
incorporating and processing the observed data. We maintain the core
algorithmic steps of an Ensemble Kalman Filter (EnKF) [18] while leveraging
the capabilities of stochastic neural networks (SNNs) [19]. There are two
steps in DEnKF, the _Prediction Step_ propagates the state one step further in
time, and the _Update Step_ corrects the state based on newly collected
observations. Let ${\bf X}_{0:N}$ denote the states of $N$ steps in $t$ with
number of $E$ ensemble members, we initialize the filtering process with ${\bf
X}_{0:N}=[{\bf x}^{1}_{0:N},\dots,{\bf x}^{E}_{0:N}]$, where
$E\in\mathbb{Z}^{+}$.

Prediction Step: In this step, the Transition Model takes the previous states
and predicts the next state. We use a window of $N$ and we leverage the
stochastic forward passes from a trained state transition model to update each
ensemble member:

$\displaystyle{\bf x}^{i}_{t}$ $\displaystyle\thicksim
f_{\boldsymbol{\theta}}({\bf x}^{i}_{t}|{\bf x}^{i}_{t-N:t-1}),\ \forall i\in
E.$ (1)

Matrix ${\bf X}_{t}=[{\bf x}^{1}_{t},\cdots,{\bf x}^{E}_{t}]$ holds the
updated ensemble members which are propagated one step forward through the
state space. Note that sampling from the transition model
$f_{\boldsymbol{\theta}}(\cdot)$ implicitly introduces a process noise.

Update Step: Given the updated ensemble members ${\bf X}_{t}$, a nonlinear
observation model $h_{\boldsymbol{\psi}}(\cdot)$ is applied to transform the
ensemble members from the state space to observation space. The observation
model is realized via a neural network with weights $\boldsymbol{\psi}$:

$\displaystyle{\bf H}_{t}{\bf X}_{t}$
$\displaystyle=\left[h_{\boldsymbol{\psi}}({\bf
x}^{1}_{t}),\cdots,h_{\boldsymbol{\psi}}({\bf x}^{E}_{t})\right],$ (2)
$\displaystyle{\bf H}_{t}{\bf A}_{t}$ $\displaystyle={\bf H}_{t}{\bf
X}_{t}-\left[\frac{1}{E}\sum_{i=1}^{E}h_{\boldsymbol{\psi}}({\bf
x}^{i}_{t}),\cdots,\frac{1}{E}\sum_{i=1}^{E}h_{\boldsymbol{\psi}}({\bf
x}^{i}_{t})\right].$

${\bf H}_{t}{\bf X}_{t}$ is the predicted observation, and ${\bf H}_{t}{\bf
A}_{t}$ is the sample mean of the predicted observation at $t$. EnKF treats
observations as random variables. Hence, the ensemble can incorporate a
measurement perturbed by a small stochastic noise thereby accurately
reflecting the error covariance of the best state estimate [18]. As shown in
Figure 2, we incorporate a Sensor Model that can learn projections between the
learned observation and raw observation space. To this end, we leverage the
methodology of SNN to train a stochastic sensor model that takes N steps of
the raw observation and predicts the current learned observation using
$s_{\boldsymbol{\xi}}(\cdot)$:

$\displaystyle\tilde{{\bf y}}^{i}_{t}$ $\displaystyle\thicksim
s_{\boldsymbol{\xi}}(\tilde{{\bf y}}^{i}_{t}|{\bf y}_{t}),\ \forall i\in E,$
(3)

where ${\bf y}_{t}$ represents the noisy observation. Sampling yields
observations $\tilde{{\bf Y}}_{t}=[\tilde{{\bf y}}^{1}_{t},\cdots,\tilde{{\bf
y}}^{E}_{t}]$ and sample mean $\tilde{{\bf
y}}_{t}=\frac{1}{E}\sum_{i=1}^{i}\tilde{{\bf y}}^{i}_{t}$. The innovation
covariance ${\bf S}_{t}$ can then be calculated as:

$\displaystyle{\bf S}_{t}$ $\displaystyle=\frac{1}{E-1}({\bf H}_{t}{\bf
A}_{t})({\bf H}_{t}{\bf A}_{t})^{T}+r_{\boldsymbol{\zeta}}(\tilde{{\bf
y}_{t}}),$ (4)

where $r_{\boldsymbol{\zeta}}(\cdot)$ is the measurement noise model
implemented using MLP. We use the same way to model the observation noise as
in [20], $r_{\boldsymbol{\zeta}}(\cdot)$ takes an learned observation
$\tilde{{\bf y}_{t}}$ in time $t$ and provides stochastic noise in the
observation space by constructing the diagonal of the noise covariance matrix.
The final estimate of the ensemble ${\bf X}_{t|t}$ can be obtained by
performing the measurement update step:

$\displaystyle{\bf A}_{t}={\bf X}_{t}-\frac{1}{E}\sum_{i=1}^{E}$
$\displaystyle{\bf x}^{i}_{t},\ {\bf K}_{t}=\frac{1}{E-1}{\bf A}_{t}({\bf
H}_{t}{\bf A}_{t})^{T}{\bf S}_{t}^{-1},$ (5) $\displaystyle{\bf X}_{t|t}$
$\displaystyle={\bf X}_{t}+{\bf K}_{t}(\tilde{{\bf Y}}_{t}-{\bf H}_{t}{\bf
X}_{t}),$

where ${\bf K}_{t}$ is the Kalman gain. In inference, the ensemble mean
${\bf\bar{x}}_{t|t}=\frac{1}{E}\sum_{i=1}^{E}{\bf x}^{i}_{t|t}$ is used as the
updated state.

Prediction Targets: Once the estimated state, comprising the rotation values
for both the lower arm and upper arm, is obtained, we utilize forward
kinematics with fixed lower arm length $l_{l}$ and fixed upper arm length
$l_{u}$ to determine the corresponding Cartesian XYZ coordinates of the wrist.
This is done due to the constraint of a fixed sagittal plane orientation
$\mathbf{r}_{h}$, where the elbow position is limited to a sphere around the
shoulder with a radius of $l_{u}$. Additionally, the wrist position must lie
on a manifold defined by spheres with a radius of $l_{l}$ around all possible
elbow positions, as described in prior work [21].

Training: DEnKF contains four sub-modules: a state transition model, an
observation model, an observation noise model, and a sensor model. The entire
framework is trained in an end-to-end manner via a mean squared error (MSE)
loss between the ground truth state $\hat{{\bf x}}_{t|t}$ and the estimated
state ${\bf\bar{x}}_{t|t}$ at every timestep. We also supervise the
intermediate modules via loss gradients
$\mathcal{L}_{f_{\boldsymbol{\theta}}}$ and
$\mathcal{L}_{s_{\boldsymbol{\xi}}}$. Given ground truth at time $t$, we apply
the MSE loss gradient calculated between $\hat{{\bf x}}_{t|t}$ and the output
of the state transition model to $f_{\boldsymbol{\theta}}$ as in Eq. 6. We
apply the intermediate loss gradients computed based on the ground truth
observation $\hat{{\bf y}_{t}}$ and the output of the stochastic sensor model
$\tilde{{\bf y}}_{t}$:

$\displaystyle\mathcal{L}_{f_{\boldsymbol{\theta}}}=\|{\bf\bar{x}}_{t|t-N:t-1}-\hat{{\bf
x}}_{t|t}\|_{2}^{2},\ \ \mathcal{L}_{s_{\boldsymbol{\xi}}}=\|\tilde{{\bf
y}_{t}}-\hat{{\bf y}_{t}}\|_{2}^{2}.$ (6)

All models in the experiments were trained for 50 epochs with batch size 256,
and a learning rate of $\eta=10^{-5}$. We chose the model with the best
performance on a validation set for testing. The ensemble size of the DEnKF
was set to 32 ensemble members.

##  III Evaluation

We discuss the validation accuracy on a separate test dataset. Secondly, we
assess the performance of our model by applying it in a human-robot handover
tasks with a real UR5 robot.

###  III-A Evaluation on Dataset

Figure 3: Prediction error distributions of our DEnKF on the test dataset.

The training dataset comprises data collected from five human subjects.
Written informed consent was obtained and approved by the institutional review
board (IRB) of ASU under the ID STUDY00017558. Each subject was instructed to
wear the smartwatch on their left arm, the smartphone in one of their pockets
and a 25-marker Optitrack suit while performing free random arm motions. The
subjects were also encouraged to change their body forward-facing direction
and to move around in the area covered by the optical tracking system. In
total, we gathered a dataset of 970,493 data points. To further augment the
dataset, we artificially adjusted the calibrated smartwatch and smartphone
data to simulate new body orientations. This is possible because remaining
sensor measurements, e.g., accelerometer or barometer, are in the reference
frame of the watch. In total, our augmented training dataset amounts to
4,259,746 data points.

Figure 4: The handover task system setup is illustrated on the left,
showcasing the workspace divided into three distinct labeled areas. On the
right, picture sequences of the handover results are presented for each of
these designated zones.

We evaluate the DEnKF prediction accuracy on a test dataset completely
separate from the training process. The test dataset comprises 26,688
observations. The participants were asked to perform a series of movements,
including arm swing, arm cross on chest and behind the head, arm raise,
waving, boxing, clapping, walking in a figure-eight pattern, and jogging in a
circular path. Figure 3 summarizes the performance of the DEnKF model on the
test dataset. On average, hand positions are off by 9.94 cm, elbow positions
by 9.27 cm, and body orientation by 7.75 deg. The DEnKF quantifies the state
uncertainty through the distribution of ensemble members, for example,
depicted at the top in Figure 1. Our framework achieves inference speed at
$\sim$62 Hz on a system using an Intel® Xeon(R) W-2125 CPU and NVIDIA GeForce
RTX 2080 Ti.

TABLE I: Model Comparison with related work

|  | Free Forward- | Wrist | Elbow | Hip
---|---|---|---|---|---
Method | Anywhere | Facing Dir. | (cm) | (cm) | (deg)
[22] | $\times$ | $\checkmark$ | $10.93$ | - | -
[23] | $\checkmark$ | $\times$ | $8.50$ | 8.50 | -
[21] | $\checkmark$ | $\times$ | $9.20$ | 7.90 | -
Ours | $\checkmark$ | $\checkmark$ | $9.94$ | 9.27 | 7.75

We compare our results with other related works [22, 23, 21] in Table I. The
method of [22] requires inference in the same environment where the training
data was collected, therefore, it is not applicable _anywhere_. Further, [22]
focuses on pose predictions for the wrist, omitting the rest of the body pose
including elbow or hip orientation. Methods of [23] and [21] demonstrate lower
errors for wrist and elbow but fix the user to a constant forward-facing
direction. In contrast, our proposed method using the DEnKF also provides an
estimate of the Hip pose and allows for ubiquitous pose estimation regardless
of location or changes in body orientation.

###  III-B Handover Task

We demonstrate the efficacy of the trained model in a human-robot handover
task. As depicted in Figure 4, participants sit in a chair and engage in
handover interactions with the robot. Participants are free to rotate on the
chair to experiment with various handover scenarios. The smartphone in their
pocket and the smartwatch on their wrist enable estimating the body
orientation and arm pose to extract the global hand position.

Task Setup: Participants perform six handover interactions. Each handover
interaction is treated as an individual task. At the onset of each task, the
participant holds their hand as shown on the left in Figure 4 and initiates
the task by issuing a voice command. Like in the work of [10], the smartwatch
recognizes the voice command and triggers the robot. Subsequently, the robot
grasps the green cube and moves it towards the tracked hand position of the
participant. When the cube reaches close proximity to the hand, the
participant says “give me the cube” and closes their hand around it. In
response to this voice command, the robot releases its grip on the cube,
signifying the successful completion of the handover.

To evaluate a range of handover positions, we divide the table surface into
three distinct areas labeled as 1, 2, and 3. Each participant is instructed to
perform two handovers in each area, with one at a high position and one at a
low position. In total, every participant performs all six handovers utilizing
smartwatch and smartphone tracking, along with an additional six handovers
conducted with the gold-standard Optitrack for our baseline comparison. The
order of the tasks and tracking modes is randomized to eliminate any potential
biases.

Results: The handover experiment utilizing smartwatches is depicted in Figure
4 (right), illustrating the handover process conducted in zones 1, 2, and 3.
We collected data from five human subjects, with each subject performing a
total of 12 handover tasks.

TABLE II: Handover task results Method | Time (s) | Dist. (cm)
---|---|---
Wearable |  $11.74$${}\pm{}$ | $5.35$ |  $7.74$${}\pm{}$ | $4.68$
Optitrack |  $9.86$${}\pm{}$ | $2.67$ |  $7.77$${}\pm{}$ | $4.71$
Difference |  $1.88$${}\pm{}$ | $4.83$ |  $0.03$${}\pm{}$ | $7.13$

The evaluation metrics in Table II include the task completion time and the
actual 3D distance between hand and cube when the participant gave the voice
command to hand over the cube. As summarized in Table II, we observe a
relatively small disparity in the task completion time, with an average
difference of 1.88 seconds. Further, observed handover 3D distances exhibit a
minimal average difference of 0.03 cm when comparing both methods. Altogether,
the results suggest that handover tasks with smartwatch and smartphone
tracking might take about 1.88 seconds longer but the participant is
comfortable to complete the handover at similar distances. It is important to
note that all handover tasks were accomplished successfully, with users
consistently grasping the cube within the desired zones without any instances
of dropping it. The obtained results indicate that the smartwatch system is
comparable to the Optitrack system for this task, thus establishing its
potential as a cost-efficient alternative.

##  IV Conclusion

This work introduces the integration of a ubiquitous robot control system by
leveraging smart devices and employing differentiable filters. The
experimental results demonstrate that the proposed framework effectively
addresses the estimation challenges related to human arm pose estimation,
especially in scenarios involving diverse human hip orientations.
Additionally, the results from the human-robot handover task showcase that the
proposed system achieves comparable error metrics, highlighting its
effectiveness. With no additional user instrumentation, the proposed framework
offers new and intriguing possibilities for low-cost robot control and human-
robot collaboration applications.

## References

  * [1] G. Nagymáté and R. M. Kiss, “Application of optitrack motion capture systems in human movement analysis: A systematic literature review,” _Recent Innovations in Mechatronics_ , vol. 5, no. 1., p. 1–9., Jul. 2018.
  * [2] C. Zimmermann, T. Welschehold, C. Dornhege, W. Burgard, and T. Brox, “3d human pose estimation in rgbd images for robotic task learning,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 1986–1992.
  * [3] Y. Desmarais, D. Mottet, P. Slangen, and P. Montesinos, “A review of 3d human pose estimation algorithms for markerless motion capture,” _Computer Vision and Image Understanding_ , vol. 212, p. 103275, 2021.
  * [4] T. von Marcard, B. Rosenhahn, M. J. Black, and G. Pons-Moll, “Sparse inertial poser: Automatic 3d human pose estimation from sparse imus,” _CoRR_ , vol. abs/1703.08014, 2017.
  * [5] “Movella xsens motion capture,” https://www.xsens.com/motion-capture/, accessed: 2023-08-24.
  * [6] X. Yi, Y. Zhou, and F. Xu, “Transpose: Real-time 3d human translation and pose estimation with six inertial sensors,” _ACM Transactions on Graphics (TOG)_ , vol. 40, no. 4, pp. 1–13, 2021.
  * [7] “Sony mocopi 3d motion capture system,” https://www.sony.jp/mocopi/, accessed: 2023-08-24.
  * [8] W. Wei, K. Kurita, J. Kuang, and A. Gao, “Real-time limb motion tracking with a single IMU sensor for physical therapy exercises,” in _2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_. IEEE, 2021, pp. 7152–7157.
  * [9] M. Liu, S. Yang, W. Chomsin, and W. Du, “Real-time tracking of smartwatch orientation and location by multitask learning,” in _Proceedings of the 20th ACM Conference on Embedded Networked Sensor Systems_ , ser. SenSys ’22. New York, NY, USA: Association for Computing Machinery, 2023, p. 120–133.
  * [10] F. C. Weigend, S. Sonawani, M. Drolet, and H. B. Amor, “Anytime, anywhere: Human arm pose from smartwatch data for ubiquitous robot control and teleoperation,” 2023.
  * [11] L. Wang, G. Wang, S. Jia, A. Turner, and S. Ratchev, “Imitation learning for coordinated human–robot collaboration based on hidden state-space models,” _Robotics and Computer-Integrated Manufacturing_ , vol. 76, p. 102310, 2022\.
  * [12] S. Chen, “Kalman filter for robot vision: a survey,” _IEEE Transactions on industrial electronics_ , vol. 59, no. 11, pp. 4409–4420, 2011.
  * [13] J. Reher, W.-L. Ma, and A. D. Ames, “Dynamic walking with compliance on a cassie bipedal robot,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 2589–2595.
  * [14] X. Liu, G. Clark, J. Campbell, Y. Zhou, and H. B. Amor, “Enhancing state estimation in robots: A data-driven approach with differentiable ensemble kalman filters,” _arXiv preprint arXiv:2308.09870_ , 2023.
  * [15] Y. Zhou, C. Barnes, J. Lu, J. Yang, and H. Li, “On the continuity of rotation representations in neural networks,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_. IEEE, 2019, pp. 5738–5746.
  * [16] X. Liu, S. Ikemoto, Y. Yoshimitsu, and H. B. Amor, “Learning soft robot dynamics using differentiable kalman filters and spatio-temporal embeddings,” _arXiv preprint arXiv:2308.09868_ , 2023.
  * [17] X. Liu, Y. Zhou, S. Ikemoto, and H. B. Amor, “$\alpha$-mdf: An attention-based multimodal differentiable filter for robot state estimation,” in _7th Annual Conference on Robot Learning_ , 2023.
  * [18] G. Evensen, “The ensemble kalman filter: Theoretical formulation and practical implementation,” _Ocean dynamics_ , vol. 53, no. 4, pp. 343–367, 2003.
  * [19] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_. PMLR, 2016, pp. 1050–1059.
  * [20] A. Kloss, G. Martius, and J. Bohg, “How to train your differentiable filter,” _Autonomous Robots_ , pp. 1–18, 2021.
  * [21] S. Shen, H. Wang, and R. Roy Choudhury, “I am a smartwatch and i can track my user’s arm,” in _Proceedings of the 14th Annual International Conference on Mobile Systems, Applications, and Services_. ACM, 2016, pp. 85–96.
  * [22] W. Wei, K. Kurita, J. Kuang, and A. Gao, “Real-time limb motion tracking with a single imu sensor for physical therapy exercises,” in _2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_. IEEE, 2021, pp. 7152–7157.
  * [23] M. Liu, S. Yang, W. Chomsin, and W. Du, “Real-time tracking of smartwatch orientation and location by multitask learning,” in _Proceedings of the 20th ACM Conference on Embedded Networked Sensor Systems_ , 2022, pp. 120–133.


# A matter dominated navigation Universe in accordance with the Type Ia
supernova data

Xin Li1,3<EMAIL_ADDRESS>1Institute of Theoretical Physics, Chinese Academy
of Sciences, 100190 Beijing, China  Zhe Chang2,3<EMAIL_ADDRESS>Minghua
Li2,3<EMAIL_ADDRESS>2Institute of High Energy Physics, Chinese Academy of
Sciences, 100049 Beijing, China
3Theoretical Physics Center for Science Facilities, Chinese Academy of
Sciences

###### Abstract

We investigate a matter dominated navigation cosmological model. The influence
of a possible drift (wind) in the navigation cosmological model makes the
spacetime geometry change from Riemannian to Finslerian. The evolution of the
Finslerian Universe is governed by the same gravitational field equation with
the familiar Friedmann-Robertson-Walker one. However, the change of space
geometry from Riemannian to Finslerian supplies us a new relation between the
luminosity distant and redshift. It is shown that the Hubble diagram based on
this new relation could account for the observations on distant Type Ia
supernovae.

###### pacs:

02.40.-k,98.80.-k

##  I Introduction

Einstein’s general relativity enables us to come up with a testable theory of
the Universe. HubbleHubble  first found that the galaxies are receding from us
not long after the birth of general relativity. The Hubble’s observation
indicates that our Universe is expanding. Over the past decade, two groups
Riess ; Perlmutter  observing supernovae reported that the luminosity distance
can not be explained by a matter dominated Universe. If one accepts the
convincible assumption of homogeneity and isotropy of the Universe which is
approximately true on large scales, then the general relativity tells us that
we now live in a dark energy dominated Universe and the Universe is
accelerated expanding. Dark energy, which has the property of negative
pressure, is different from the classical matter and the found particles. A
great amount of models have been proposed to study the possible candidates of
dark energy and their dynamics Copeland . The most famous and acceptable
candidate is the cosmological constant. However, the magnitude of cosmological
constant $10^{-3}{\rm eV}^{4}$ is much smaller than the energy density of
vacuum in quantum field theory.

The theory of dark energy Padmanabhan  dominates the modern cosmology in the
past decade. This situation is similar to the raise of the theory of ether,
which has been considered as direct evidence of the aberration of starlight,
an important astronomical effect known since eighteenth century. Although the
rapid progress in technology makes the astronomical observations more and more
accurate, up to now, there is no direct evidence indicate what the dark energy
it is. Since the theory of dark energy is contrived, which requires fine
tuning and apparently cannot be tested in the laboratory or solar system,
several modified theories of general relativity have been developed as the
alternative source of cosmological acceleration Bludman .

Einstein first used the Riemann geometry to describe the theory of
gravitation. In this Letter, we suppose that the spacetime in large scale may
be described by other geometry instead of Riemann geometry. To preserve
redeeming feature of general relativity, this geometry must take Riemann
geometry as its special case. Fortunately, Paul Finsler proposed a natural
generation of Riemann geometry-Finsler geometry.

Finsler geometry as a more general geometry could provide new sight on modern
physics. It is of great interest for physicists to investigated the violation
of Lorentz symmetry Kostelecky . An interesting case of Lorentz violation,
which was proposed by Cohen and GlashowGlashow , is the model of Very Special
Relativity (VSR) characterized by a reduced symmetry SIM(2). In fact, Gibbons,
Gomis and PopeGibbons  showed that the Finslerian line element
$ds=(\eta_{\mu\nu}dx^{\mu}dx^{\nu})^{(1-b)/2}(n_{\rho}dx^{\rho})^{b}$ is
invariant under the transformations of the group DISIM${}_{b}(2)$. In the
framework of Finsler geometry, modified dispersion relation has been
discussedGirelli ; Finsler SR . Also, the model based on Finsler geometry
could explain the recent astronomical observations which Einstein’s gravity
could not. A list includes: the flat rotation curves of spiral galaxies can be
deduced naturally without invoking dark matter Finsler DM ; the anomalous
accelerationAnderson  in solar system observed by Pioneer 10 and 11
spacecrafts should correspond to Finsler-Randers space Finsler PA ; the
secular trend in the astronomical unitKrasinsky ; Standish  and the anomalous
secular eccentricity variation of the Moon’s orbitWilliams  should be
correspond to the effect of the length change of unit circle in Finsler
geometryFinsler AU .

In this Letter, we present a matter dominated navigation model. The influence
of a possible drift (wind) in the navigation cosmolgical model makes the
expanding Universe “accelerated”. It is remarkable that the navigation
cosmological model is described in the framework of Finsler geometry. We find
that the predictions of the navigation cosmological model could account for
the observations of Riess and PerlmutterRiess ; Perlmutter  on distant
supernovae.

##  II formulism

Instead of defining an inner product structure over the tangent bundle in
Riemann geometry, Finsler geometry is based on the so called Finsler structure
$F$ with the property $F(x,\lambda y)=\lambda F(x,y)$ for all $\lambda>0$,
where $x$ represents position and $y$ represents velocity. The Finsler metric
is given asBook by Bao

$g_{\mu\nu}\equiv\frac{\partial}{\partial y^{\mu}}\frac{\partial}{\partial
y^{\nu}}\left(\frac{1}{2}F^{2}\right).$ (1)

In Finsler manifold, there exists a unique linear connection - the Chern
connectionChern . It is torsion freeness and almost metric-compatibility,

$\Gamma^{\alpha}_{\mu\nu}=\gamma^{\alpha}_{\mu\nu}-g^{\alpha\lambda}\left(A_{\lambda\mu\beta}\frac{N^{\beta}_{\nu}}{F}-A_{\mu\nu\beta}\frac{N^{\beta}_{\lambda}}{F}+A_{\nu\lambda\beta}\frac{N^{\beta}_{\mu}}{F}\right),$
(2)

where $\gamma^{\alpha}_{\mu\nu}$ is the formal Christoffel symbols of the
second kind with the same form of Riemannian connection, $N^{\mu}_{\nu}$ is
defined as
$N^{\mu}_{\nu}\equiv\gamma^{\mu}_{\nu\alpha}y^{\alpha}-A^{\mu}_{\nu\lambda}\gamma^{\lambda}_{\alpha\beta}y^{\alpha}y^{\beta}$
and $A_{\lambda\mu\nu}\equiv\frac{F}{4}\frac{\partial}{\partial
y^{\lambda}}\frac{\partial}{\partial y^{\mu}}\frac{\partial}{\partial
y^{\nu}}(F^{2})$ is the Cartan tensor (regarded as a measurement of deviation
from the Riemannian Manifold). In terms of Chern connection, the curvature of
Finsler space is given as

$R^{~{}\lambda}_{\kappa~{}\mu\nu}=\frac{\delta\Gamma^{\lambda}_{\kappa\nu}}{\delta
x^{\mu}}-\frac{\delta\Gamma^{\lambda}_{\kappa\mu}}{\delta
x^{\nu}}+\Gamma^{\lambda}_{\alpha\mu}\Gamma^{\alpha}_{\kappa\nu}-\Gamma^{\lambda}_{\alpha\nu}\Gamma^{\alpha}_{\kappa\mu},$
(3)

where $\frac{\delta}{\delta x^{\mu}}=\frac{\partial}{\partial
x^{\mu}}-N^{\nu}_{\mu}\frac{\partial}{\partial y^{\nu}}$. The notion of Ricci
tensor in Finsler geometry was first introduced by Akbar-ZadehAkbar

$Ric_{\mu\nu}=\frac{\partial^{2}\left(\frac{1}{2}F^{2}R\right)}{\partial
y^{\mu}\partial y^{\nu}},$ (4)

where $R=\frac{y^{\mu}}{F}R^{~{}\kappa}_{\mu~{}\kappa\nu}\frac{y^{\nu}}{F}$.
And the scalar curvature in Finsler geometry is given as
$S=g^{\mu\nu}Ric_{\mu\nu}$.

In standard cosmology, following the cosmological principle, one gets the
Friedmann-Robertson-Walker (FRW) metricWeinberg . In another word, the spatial
part of the Universe is a constant sectional curvature space. Here comes our
major assumption, the gravity in large scale should be described in terms of
Finsler geometry. In light of the cosmological principle, the Finsler
structure of the Universe should be written in such form

$\bar{F}^{2}=dt^{2}-R^{2}(t)F^{2},$ (5)

where $R(t)$ is scale factor with cosmic time $t$, the structure $F$ is a
constant flag curvature space. By making use of the geometrical terms we
mentioned above, one obtains the components of Ricci tensor

$\displaystyle Ric_{00}$ $\displaystyle=$
$\displaystyle-\frac{3\ddot{R}}{R}g_{00},$ (6) $\displaystyle Ric_{ij}$
$\displaystyle=$
$\displaystyle-\left(\frac{\ddot{R}}{R}+\frac{2\dot{R}^{2}}{R^{2}}+\frac{2K}{R^{2}}\right)g_{ij},$
(7)

where the dot denotes a derivative with respect to $t$.

The gravitational field equation in Finsler space should reduce to the
Einstein’s gravitational field equation while the Finsler space reduce to the
Riemannian space. Thus, the symmetrical tensor $G_{\mu\nu}\equiv
Ric_{\mu\nu}-\frac{1}{2}g_{\mu\nu}S$ should involve in the gravitational field
equation in Finsler space. In general, Finsler space is an anisotropic space.
Therefore, it has less Killing vectors than the Riemannian space, and it
breaks the Lorentz symmetryLi . It means the angular momentum in Finsler space
is not a conservative quantity. This fact implies that the energy momentum
tensor is not symmetrical in Finsler space. For these reasons, the
gravitational field equation in Finsler space should be taken in such form

$G_{\mu\nu}+A_{\mu\nu}=8\pi G(T^{s}_{\mu\nu}+T^{a}_{\mu\nu}),$ (8)

where $A_{\mu\nu}$ is an asymmetrical tensor and
$T^{s}_{\mu\nu},T^{a}_{\mu\nu}$ are symmetrical part and asymmetrical part of
energy momentum tensor respectively. This general form is agree with the
result of Asanov, its gravitational field equation contains the asymmetrical
term in Finsler space of Landsberg typeAsanov . Dealing with field equation in
the anisotropic and inhomogeneous cosmology is not a simple staff, and it is
hard to find an exact solution of gravitational field equation while its
energy momentum tensor involves the non diagonal termsMisner . At first
glance, we just deal with the symmetrical part of field equation (8)

$G_{\mu\nu}=8\pi GT_{\mu\nu}.$ (9)

By making use of the equations (6) and (7), the solution of equation (9) is
the same with the Einstein’s field equation deduced by FRW metric. Taking the
energy-momentum tensor $T_{\mu\nu}$ to be the form of perfect fluid, one can
see that the evolution of scale factor $R(t)$ is the same with the Riemannian
case. However, since the luminosity distant is related to the spatial geometry
of the Universe, one may expect that it is different from the Riemannian
luminosity distant.

Unlike Riemann space, a complete classification of the constant flag curvature
spaces remain an unsolved problem. However, by making use of the Zermolo
navigation on Riemannian space, Bao et al. Bao  gave a complete classification
of Finsler-Randers space Randers  of constant flag curvature. The Zermelo
navigation problem Zermelo  aims to find the paths of shortest travel time in
a Riemannian manifold ($M,h$) under the influence of a drift (“wind”)
represented by a vector field $W$. In standard cosmology, our Universe is very
flat now. We may imagine that the Universe is a flat Riemannian manifold with
flat Friedmann metric

$ds^{2}=dt^{2}-R_{h}^{2}(t)(dx^{2}+dy^{2}+dz^{2}),$ (10)

and it influenced by the “wind” $W$. The relation between the Riemannian
manifold ($M,h$) and the Randers metric $F$ is

$F=\frac{1}{\lambda}\left(\sqrt{\lambda
h_{ij}y^{i}y^{j}+(W_{i}y^{i})^{2}}-W_{i}y^{i}\right),$ (11)

where $W_{i}=h_{ij}W^{j}$ and $\lambda=1-h(W,W)$. One should notice that there
is a map between the Riemannian space which influence by the “wind” and the
Randers-Finsler spaceGibbons1 . It means the effect of “wind” already
accounted in the gravitational field equation in Finsler space. The theoremBao
of the classification supplies an interesting case where the Randers metric
$F$ has constant flag curvature $K$: for $K=-\frac{1}{16}\sigma^{2}<0$ and $h$
is flat. And $\sigma$ satisfies the constraint $\mathcal{L}_{W}h=-\sigma h$,
$\mathcal{L}$ denotes Lie differentiation. We set the vector field to be
radial $W=\epsilon(t)(x_{1}\partial x_{1}+x_{2}\partial x_{2}+x_{3}\partial
x_{3})$. Then the flag curvature of Randers metric is
$K=-\frac{1}{4}\epsilon^{2}$. After doing coordinate change and taking
spherical coordinate, we get the Randers metric as

$F=\sqrt{\frac{dr^{2}}{1+\epsilon^{2}r^{2}}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})}-\frac{\epsilon rdr}{1+\epsilon^{2}r^{2}}.$ (12)

After scale change $r\rightarrow 2r/\epsilon$, the metric $F$ changes as

$F=\sqrt{\frac{dr^{2}}{1+4r^{2}}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})}-\frac{2rdr}{1+4r^{2}},$ (13)

and the cosmic scale factor changes as $R(t)=\frac{2R_{h}(t)}{\epsilon(t)}$.
Taking the spatial part of the Finsler structure (5) to be of the form (13),
we get metric of the Universe in the framework of Finsler geometry, and the
space curvature of the Universe is $-1$. One can deduce directly from such
metric that the relation between the redshift $z$ and the scale factor $R(t)$
is the same with the Riemannian case

$1+z=\frac{R_{0}}{R(t)}=\frac{1}{a(t)},$ (14)

where the subscript zero denotes the quantities given at the present epoch.
Since the non-radial part of the metric (13) is the same with FRW metric, the
luminosity distant in the navigation cosmological model is given as
$d_{L}=R_{0}r(1+z)$. However, the proper distant is not the case, and it is
given as

$\displaystyle d_{p}$ $\displaystyle=$
$\displaystyle\int_{0}^{r}\left(\frac{1}{\sqrt{1+4r^{2}}}-\frac{2r}{1+4r^{2}}\right)dr$
(15) $\displaystyle=$
$\displaystyle\frac{\sinh^{-1}2r}{2}-\frac{1}{4}\ln(1+4r^{2}).$

The light traveling along the radial direction satisfies the geodesic equation

$\bar{F}^{2}=dt^{2}-R^{2}(t)\left(\frac{1}{\sqrt{1+4r^{2}}}-\frac{2r}{1+4r^{2}}\right)^{2}dr^{2}=0.$
(16)

Then, we obtain

$d_{p}=\frac{1}{R_{0}}\int^{1}_{(1+z)^{-1}}\frac{da}{a\dot{a}}$ (17)

The equations (9) and (17) can be solved. Supposing the Universe is matter
dominated (no more cosmological constant and dynamical dark energy), we obtain

$d_{p}=\log\frac{(\sqrt{1+\Omega_{m}^{(0)}z}-\sqrt{1-\Omega_{m}^{(0)}})(1+\sqrt{1-\Omega_{m}^{(0)}})}{(\sqrt{1+\Omega_{m}^{(0)}z}+\sqrt{1-\Omega_{m}^{(0)}})(1-\sqrt{1-\Omega_{m}^{(0)}})},$
(18)

where $H_{0}$ is the Hubble constant and $\Omega_{m}$ is the density parameter
for matter and satisfies

$1-\Omega_{m}^{(0)}=\frac{1}{H_{0}^{2}R_{0}^{2}}.$ (19)

Substituting the equation (18) into (15), we obtain the relation between
luminosity distant and redshift in the navigation cosmological model

$H_{0}d_{L}=\frac{1+z}{2\sqrt{1-\Omega_{m}^{(0)}}}\frac{|e^{2d_{p}}-1|}{e^{d_{p}}\sqrt{2-e^{2d_{p}}}}.$
(20)

##  III numerical result

Figure 1: The luminosity distant $Log_{10}(H_{0}d_{L}/c)$ versus the redshift
$z$ for the navigational cosmological model. The data comes from Riess et
al.Riess  and Hubble Space Telescope (HST)Kowalski . And the value of Hubble
constant is set as $H_{0}=67.88km\cdot s^{-1}\cdot Mpc^{-1}$.

Here, we present numerical result for the relation between luminosity distant
and redshift in the framework of the navigation cosmological model. The best
fit curve is shown in Fig.1 with the matter density parameter taken as
$\Omega_{m}^{(0)}=0.92$. And the average value of Chisquare is $1.077588$.
Thus, our prediction could account for experiment data given by Riess et
al.Riess . It should be noticed that the Hubble constant $H_{0}$ we took is
different from the Hubble constant $H_{h0}$ measured in flat space. The
relation for $H$ and $H_{h}$ is

$H=H_{h}-\frac{\dot{\epsilon}}{\epsilon}.$ (21)

By making use of the value of the Hubble constant measured in flat space
$H_{h0}=73\pm 0.3~{}km\cdot s^{-1}\cdot Mpc^{-1}$Spergel , we have

$H_{\epsilon 0}\equiv\frac{\dot{\epsilon}_{0}}{\epsilon_{0}}=5.42km\cdot
s^{-1}\cdot Mpc^{-1}$ (22)

This geometrical parameter $H_{\epsilon 0}$ represents an “accelerated” effect
provided by the vector field $W$.

##  IV conclusion

Our Letter initiates an exploration of the possibility that the empirical
success of the observations of Type Ia supernovae can be regarded as the
influence of a navigational wind. The particles move on Riemnnian manifold and
influenced by a vector field (the “wind”) which proportion to the curvature
$-\frac{\epsilon^{2}}{4}$ of the space of the Universe. Its geodesic indeed is
a Finslerian geodesic. Thus, the observations of Riess et al.Riess  on distant
supernovae may be explained by the effect of the “wind”. Our numerical result
could account for astronomical observations. At last, we point out that the
age of the Universe is about 9.76Gyr in our model. This is contradicted with
the age ($13.5\pm 2Gyr$) of Globular clusters in the MilkyWayJimenez . In the
navigation cosmological model, we only involve the radial “wind”. The non-
radial “wind” should be taken into account in future work, we expect that the
effect of the non-radial “wind” may supply us a reasonable age of the
Universe.

###### Acknowledgements.

We would like to thank Prof. C. J. Zhu and X. H. Mo for useful discussions.
The work was supported by the NSF of China under Grant No. 10525522 and
10875129\.

## References

  * (1) E. Hubble, Proceedings of the National Academy of Sciences of the United states of America 15, 168 (1929).
  * (2) A. G. Riess et al., Astron. J. 116, 1009 (1998); Astron. J. 117, 707 (1999).
  * (3) S. Perlmutter, et al., Astrophys J. 517, 565 (1999).
  * (4) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006).
  * (5) T. Padmanabhan, Phys. Rept. 380, 235 (2003).
  * (6) S. Bludman, arXiv:astro-ph/0605198.
  * (7) D. Colladay and V. A. Kostelecky, Phys. Rev. D55, 6760 (1997); V. A. Kostelecky, Phys. Rev. D 69, 105009 (2004).
  * (8) A. G. Cohen and S.L. Glashow, Phys. Rev. Lett. 97, 021601 (2006).
  * (9) G. W. Gibbons, J. Gomis and C. N. Pope, Phys. Rev. D 76, 081701 (2007).
  * (10) F. Girelli, S. Liberati and L. Sindoni, Phys. Rev. D 75, 064015 (2007).
  * (11) Z. Chang and X. Li, Phys. Lett. B 663, 103 (2008).
  * (12) Z. Chang and X. Li, Phys. Lett. B 668, 453 (2008).
  * (13) J. D. Anderson et al., Phys. Rev. Lett. 81 2858, (1998); Phys. Rev. D 65, 082004 (2002); Mod. Phys. Lett. A 17, 875 (2002).
  * (14) X. Li and Z. Chang, arXiv:gr-qc/0909.3713.
  * (15) G. A. Krasinsky and V. A. Brumberg, Celest. Mech. Dyn. Astrn. 90, 267 (2004).
  * (16) E. M. Standish, Proc. IAU Colloq. 196, 163 (2005).
  * (17) J. G. Williams and D. H. Boggs, in Proceedings of 16th International Workshop on Laser Ranging, ed. S. Schillak, (Space Research Centre, Polish Academy of Sciences), 2009.
  * (18) X. Li and Z. Chang, arXiv:gr-qc/0911.1890.
  * (19) D. Bao, S. S. Chern and Z. Shen, An Introduction to Riemann–Finsler Geometry, Graduate Texts in Mathmatics 200, Springer, New York, 2000.
  * (20) S. S. Chern, Sci. Rep. Nat. Tsing Hua Univ. Ser. A 5, 95 (1948); or Selected Papers, vol. II, 194, Springer 1989.
  * (21) H. Akbar-Zadeh, Acad. Roy. Belg. Bull. Cl. Sci. (5) 74, 281 (1988).
  * (22) S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, Wiley, New York, 1972.
  * (23) X. Li, Z. Chang and X. H. Mo, arXiv:hep-th/1001.2667.
  * (24) G. S. Asanov, Reports on Mathematical Physics 59, 111 (2007).
  * (25) C. W. Misner, K, S, Thorne and J. A. Wheeler, Gravitation, Freeman, San Francisco, 1973.
  * (26) D. Bao, C. Robles and Z. Shen, J. Diff. Geom. 66, 377 (2004); arXiv:math/0311233.
  * (27) G. Randers, Phys. Rev. 59, 195 (1941).
  * (28) E. Zermelo, Z. Angew. Math. Mech. 11(2), 114 (1931).
  * (29) G. W. Gibbons, et al., Phys. Rev. D 79, 044022 (2009).
  * (30) M. Kowalski et al., Astrophys. J. 686, 749 (2008).
  * (31) D. N. Spergel et al., Astrophys. J. Supp. 170, 377 (2007).
  * (32) R. Jimenez, P. Thejll, U. Jorgensen, J. MacDonald and B. Pagel, MNRAS 282, 926 (1996).