# Spin-0 bosons near rotating stars

E. O. Pinho<EMAIL_ADDRESS>Departamento de Física, CFM -
Universidade Federal de
Santa Catarina; C.P. 476, CEP 88.040-900, Florianópolis, SC, Brazil  C. C.
Barros Jr<EMAIL_ADDRESS>Departamento de Física, CFM - Universidade
Federal de
Santa Catarina; C.P. 476, CEP 88.040-900, Florianópolis, SC, Brazil

###### Abstract

In this work we study the effect of rotating stars in the behavior of bosons
close to their surfaces. For this task, metrics determined by the rotation of
these stars will be taken into account. We will consider the Klein-Gordon
equation in the Hartle-Thorne metric and in the one proposed by Berti et al.
by considering some kinds of stars. Pions and Higgs bosons will be
investigated as examples. By solving the equations, it may be observed that
the main effect that may be observed is a rotational phase, $\delta_{R}$, that
depends on the the angular momentum of the star, $J$, in the wave function of
the boson. Corrections of higher orders in $J$ are also investigated.

##  I Introduction

The observation of the effects of the general theory of relativity in physical
systems always has been a very challenging task. Since the first observation
of light deflection by the Sun or the effect of the frame dragging due to the
rotation of the earth Schiff (1960), Lense and Thirring (1918), the Lense-
Thirring effect, that took decades to be observed by the Gravity probe B
Everitt  _et al._ (2011) or even the recent detection of gravitational waves
Abbot  _et al._ (2016), we have witnessed great experimental efforts in order
to obtain reliable results.

A fundamental question in physics is how quantum mechanics and the general
relativity are related, or how the general relativity affects quantum systems.
At the moment a large number of this kind of systems have been studied. Since
the early studies from Parker Parker (1980) many other have been done, as for
example quantum oscillators Ahmed (2022), Ahmed (2023), Santos  _et al._
(2019), Soares  _et al._ (2021), Yang  _et al._ (2021), Rouabhia  _et al._
(2023), magnetic fields in the Melvin metric Santos and Barros Jr. (2016),
cosmic strings Santos and Barros Jr. (2017), Santos and Barros Jr. (2018a),
Vitória and Bakke (2018a), the Casimir effect Bezerra  _et al._ (2016), Santos
and Barros Jr. (2018b) and many other kinds of systems Sedaghatnia  _et al._
(2019) Guvendi  _et al._ (2022), Vitória and Bakke (2018b) Barros Jr. (2005),
Barros Jr. (2006).

In this work we will study the effect of the dragging of the spacetime due to
the rotation of stars in quantum systems, and for this purpose we will take
spin-0 bosons as test particles. A description of the spacetime of “slowly”
rotating relativistic stars was constructed during the 1960s by James B.
Hartle and Kip S. Thorne Hartle (1967), Hartle and Thorne (1968), for both the
exterior and interior regions of a star. The limits of the approximation, or
in other words, the “slowness” of the stars, in fact are not so narrow when
considering physical stars, and this criterion is not valid only in the slow-
rotation realm, yet it seems as if the solution is accurate for many systems,
even for some that would not have been considered slow then.

So in this paper, we will use the external Hartle-Thorne metric as given by
Hartle and Thorne (1968) and also the corrections proposed in Berti  _et al._
(2005), and then write a Klein-Gordon equation. We will solve the equations
and then show numerical results for some typical systems.

This paper is divided into five sections beyond the introduction. In section
II we will describe the line elements we will be using in our calculations and
in Sec. III we will sketch the process by which we take the line element and
expand it up to first order in the star angular momentum $J$. Then we will
write the Klein-Gordon equation and solve it in various ways, in terms of
several different special functions. In section IV, we will show an expansion
that contains terms up to second order of the angular momentum $J$ and the
resulting power series solution to the Klein-Gordon equation obtained in this
spacetime. The resulting solution corresponds to the wave function for a free
particle plus additional mass and angular momentum corrections. We will also
compare this solution with the one obtained in section III, thereby showing
that they have intrinsic similarities, despite having been obtained in
conceptually different ways. In Sec. V the solution will then be studied for
some physical systems with different values of the mass and angular momentum
and in Sec. VI we will draw the conclusions of this work.

##  II The Two Metrics

In this section we will study the spacetime structure by means of metrics that
take into account the rotation of stars. For this purpose we will consider the
Hartle-Thorne metric Hartle (1967), Hartle and Thorne (1968), and the metric
proposed by Berti et al. Berti  _et al._ (2005). As far as we are interested
in studying observable effects of bosons near rotating stars, we will only
consider the solutions of the equations in a external region.

The Hartle-Thorne metric is a well-known spacetime metric that describes the
geometry determined by slowly rotating stars. It was first obtained in a
sequence of papers in the late 1960s and it has been used ever since in a
variety of important works. The metric was initially constructed out of a very
general canvas, that of an axially symmetric, stationary system that rotates
uniformly and “slowly”. This “slow” rotation refers to the fluid’s angular
velocity $\Omega$, which must be slow enough for the fractional changes in
pressure, energy density, and gravitational field due to rotation to be small.
These considerations lead to the criterion presented by the authors

$\Omega^{2}\ll\left(\frac{c}{R}\right)^{2}\frac{GM}{Rc^{2}}$ (1)

for a star of mass $M$ and raduis $R$. We note, however, that this criterion
is actually wide-reaching, unlike it may be thought in a first view it and
does not necessarily correspond to a slow rotation if compared with the values
of real stars Krastev  _et al._ (2008).

So, a general line element may be proposed in order to study a rotating object
by considering it as an axially symmetric stationary system in terms of
spherical coordinates

$ds^{2}=-H^{2}dt^{2}+E^{2}dr^{2}+r^{2}K^{2}[d\theta^{2}+\sin^{2}{\theta}(d\phi-
Ldt)^{2}],$ (2)

where $H$, $E$, $K$, and $L$ are functions of $r$ and $\theta$. The rotation
is introduced as a perturbation and then the line element can be written as
Hartle and Thorne (1968)

$\begin{split}ds^{2}=&-\left(1-\frac{2M}{r}-\frac{2J^{2}}{r^{4}}\right)\times\\\
&\times\Bigg{\\{}1+2\Bigg{[}\frac{J^{2}}{Mr^{3}}\left(1+\frac{M}{r}\right)+\frac{5}{8}\frac{Q-\frac{J^{2}}{M}}{M^{3}}{Q_{2}}^{2}\left(\frac{r}{M}-1\right)\Bigg{]}P_{2}(\cos{\theta})\Bigg{\\}}dt^{2}\\\
&+\left(1-\frac{2M}{r}-\frac{2J^{2}}{r^{4}}\right)^{-1}\times\\\
&\times\Bigg{\\{}1-2\Bigg{[}\frac{J^{2}}{Mr^{3}}\left(1-\frac{5M}{r}\right)+\frac{5}{8}\frac{Q-\frac{J^{2}}{M}}{M^{3}}{Q_{2}}^{2}\left(\frac{r}{M}-1\right)\Bigg{]}P_{2}(\cos{\theta})\Bigg{\\}}dr^{2}\\\
&+r^{2}\Bigg{\langle}\Bigg{\langle}1+2\Bigg{\langle}-\frac{J^{2}}{Mr^{3}}\left(1+\frac{M}{r}\right)+\frac{5}{8}\frac{Q-\frac{J^{2}}{M}}{M^{3}}\Bigg{\\{}\frac{2M}{\sqrt{r(r-2M)}}{Q_{2}}^{1}\left(\frac{r}{M}-1\right)+\\\
&-{Q_{2}}^{2}\left(\frac{r}{M}-1\right)\Bigg{\\}}\Bigg{\rangle}P_{2}(\cos{\theta})\Bigg{\rangle}\Bigg{\rangle}\left\\{d\theta^{2}+\sin^{2}{\theta}\left[d\phi-\frac{2J}{r^{3}}dt\right]^{2}\right\\}\end{split}$
(3)

where $Q$ is the star’s mass quadrupole moment, $J$ its angular momentum, and

$\begin{split}{Q_{2}}^{1}(\zeta)&=\sqrt{\zeta^{2}-1}\left[\frac{3\zeta^{2}-2}{\zeta^{2}-1}-\frac{3}{2}\zeta\ln{\frac{\zeta+1}{\zeta-1}}\right]\\\
{Q_{2}}^{2}(\zeta)&=-\frac{3\zeta^{3}-5\zeta}{\zeta^{2}-1}+\frac{3}{2}(\zeta^{2}-1)\ln{\frac{\zeta+1}{\zeta-1}},\end{split}$
(4)

are associated Legendre functions of the second kind. In Abramowicz  _et al._
(2003) it was pointed that this metric is not consistently truncated, and for
this reason it presents some errors. Then a corrected metric has been
proposed, but Berti et al. Berti  _et al._ (2005) studied this metric and
found some minor sign errors. They proposed the solution that will be the
second metric we will be using in this paper and it may be considered as a
corrected version of the original Hartle-Thorne metric. From Berti  _et al._
(2005), the corresponding line element is given by

$\begin{split}g_{rr}&=\left(1-\frac{2M}{r}\right)^{-1}[1+j^{2}G_{1}+qF_{2}]+\mathcal{O}(\epsilon^{3})\\\
g_{tt}&=-\left(1-\frac{2M}{r}\right)[1+j^{2}F_{1}+qF_{2}]+\mathcal{O}(\epsilon^{3})\\\
g_{\theta\theta}&=r^{2}[1+j^{2}H_{1}-qH_{2}]+\mathcal{O}(\epsilon^{3})\\\
g_{\phi\phi}&=g_{\theta\theta}\sin^{2}{\theta}+\mathcal{O}(\epsilon^{3})\\\
g_{t\phi}&=\left(\frac{2jM^{2}}{r}\right)\sin^{2}{\theta}+\mathcal{O}(\epsilon^{3})\end{split}$
(5)

where

$\begin{split}F_{1}&=-pW+A_{1}\\\
F_{2}&=5r^{3}p(3u^{2}-1)(r-M)(2M^{2}+6Mr-3r^{2})-A_{1}\\\
A_{1}&=\frac{15r(r-2M)(1-3u^{2})}{16M^{2}}\ln{\frac{r}{r-2M}}\\\
A_{2}&=\frac{15r(r-2M^{2})(3u^{2}-1)}{16M^{2}}\ln{\frac{r}{r-2M}}\\\
G_{1}&=p[(L-72M^{5}r)-3u^{2}(L-56M^{5}r)]-A_{1}\\\
H_{1}&=A_{2}+\frac{(1-3u^{2})(16M^{5}+8M^{4}r-10M^{2}r^{3}+15Mr^{4}+15r^{5})}{8Mr^{4}}\\\
H_{2}&=-A_{2}+\frac{5(1-3u^{2})(2M^{2}-3Mr-3r^{2})}{8Mr}\\\
L&=80M^{6}+8M^{4}r^{2}+10M^{3}r^{3}+20M^{2}r^{4}-45Mr^{5}+15r^{6}\\\
p&=\frac{1}{8Mr^{r}(r-2M)}\\\
W&=(r-M)(16M^{5}+8M^{4}r-10M^{2}r^{3}-30Mr^{4}+15r^{5})\\\
&+u^{2}(48M^{6}-8M^{5}r-24M^{4}r^{2}-30M^{3}r^{3}-60M^{2}r^{4}+135Mr^{5}-45r^{6})\\\
j&=\frac{J}{M^{2}}\\\ q&=\frac{Q}{M^{3}}\\\ u&=\cos{\theta}\end{split}$ (6)

and $\epsilon=\Omega/\Omega^{*}$ with $\Omega^{*}=\sqrt{M/R^{3}}$. Observing
these results it seems to be reasonable to consider both metrics in the study
of quantum systems. As far as the corrections proposed in Berti  _et al._
(2005) are very small, we will verify if they are important in the behavior of
spin-0 bosons near some kinds of stars.

##  III Klein-Gordon equation in the Hartle-Thorne metric

In this section we will study the Klein-Gordon equation in the metrics
presented in Sec. II up to the first order in the angular momentum ratio
$j_{1}=J/r^{2}$. We must remark that it is a good quantity for the analysis
that we intend to do in this work since it presents small numerical values for
the stars. For the sun, taking $r=R$ for example, $j_{1}=9.7\times 10^{-13}$,
for a typical white dwarf, $j_{1}\sim 10^{-10}$ and for a typical neutron star
$j_{1}\sim 10^{-5}$. As we can see in a first approximation we may neglect
terms of the order $j_{1}^{2}$ as far it is a very small number in a region
outside the star. So, in this section we will work out the solutions up to the
first order in $j_{1}$, and in next section we will study the corrections up
to $j_{1}^{2}$.

In Hartle and Thorne (1968) the line element (3) has been expanded in powers
of $M/r$, that for the external region of the Sun for example, is smaller than
$M/R\sim 2\times 10^{-6}$ and when using numerical data relative to the Sun,
an approximate form of this line element with an accuracy of the order of 1
part in $10^{15}$ has been proposed

$\begin{split}ds^{2}=&-\left[1-\frac{2M}{r}+\frac{2Q}{r^{3}}P_{2}(\cos{\theta})\right]dt^{2}+\left[1-\frac{2M}{r}+\frac{2Q}{r^{3}}P_{2}(\cos{\theta})\right]^{-1}dr^{2}\\\
&+r^{2}\left[1-\frac{2Q}{r^{3}}P_{2}(\cos{\theta})\right]\left[d\theta^{2}+\sin^{2}{\theta}\left(d\phi-\frac{2J}{r^{3}}dt\right)^{2}\right]\end{split}$
(7)

We will study stars with nearly symmetric mass distributions by taking $Q=0$,
and as far as $J^{2}/r^{6}<<M/r$ for the external line element we may write

$\begin{split}ds^{2}=&-\left(1-\frac{2M}{r}\right)dt^{2}-\frac{4J}{r}\sin^{2}{\theta}dtd\phi\\\
&+\left(1-\frac{2M}{r}\right)^{-1}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}{\theta}d\phi^{2}\right),\end{split}$
(8)

that is the metric that will be considered in the calculations.

The Klein-Gordon equation for an arbitrary spacetime may be obtained by
replacing the regular derivatives in the usual Klein-Gordon equation by
covariant derivatives, and then we get Birrell and Davies (1982), Parker and
Toms (2009)

$\left[-\frac{1}{\sqrt{-g}}\partial_{\mu}\big{(}g^{\mu\nu}\sqrt{-g}\partial_{\nu}\big{)}+\mu^{2}\right]\psi=0.$
(9)

Since the components of our metric are not functions of $t$ or $\phi$, we can
do the separation of variables in the equation and obtain

$\psi(t,r,\theta,\phi)=e^{im\phi}e^{-i\omega t}\Theta(\theta)R(r),$ (10)

where $m=0,\pm 1,\pm 2...$ and we still have to solve the equations for $R(r)$
and $\Theta(\theta)$. The equation for $\Theta(\theta)$ is

$\left[\frac{d^{2}}{d\theta^{2}}+\frac{\cos\theta}{\sin\theta}\frac{d}{d\theta}+\lambda^{\prime}-\frac{m^{2}}{\sin^{2}\theta}\right]\Theta(\theta)=0,$
(11)

which solution is simply the set of associated Legendre functions
$P_{l}^{m}(\theta)$ with $\lambda^{\prime}=l(l+1)$ where $l$ is an integer.
The radial equation is

$\left(\mu^{2}-\omega^{2}A+\frac{l(l+1)}{r^{2}}+\frac{4m\omega
JA}{r^{3}}\right)R-\frac{1}{r^{2}}\frac{d}{dr}\left(\frac{r^{2}}{A}\frac{dR}{dr}\right)=0,$
(12)

where $A=\left(1-\frac{2M}{r}\right)^{-1}$. By substituting
$R(r)=\frac{\sqrt{A}}{r}u(r)$, making the change of variables
$x=\frac{r}{2M}$, and rearranging fractions, we arrive at the following
equation

$\frac{d^{2}u}{dx^{2}}+\left[a+\frac{b}{x^{2}}+\frac{c}{x}+\frac{d}{(x-1)^{2}}+\frac{e}{x-1}\right]u=0,$
(13)

where

$\begin{split}a&=4M^{2}(\omega^{2}-\mu^{2})\\\ b&=\frac{1}{4}\\\
c&=\frac{1}{2}+l(l+1)-\frac{4m\omega J}{M}\\\
d&=\frac{1}{4}+4M^{2}\omega^{2}-\frac{4m\omega J}{M}\\\
e&=-\frac{1}{2}-4M^{2}\mu^{2}-l(l+1)+8M^{2}\omega^{2}+\frac{4m\omega J}{M}\
.\end{split}$ (14)

This equation is the normal form of the confluent Heun equation, and its
canonical form is given by

$\frac{d^{2}H}{dx^{2}}+\left(\alpha+\frac{\beta+1}{x}+\frac{\gamma+1}{x+1}\right)\frac{dH}{dx}+\left(\frac{\sigma}{x}+\frac{\tau}{x-1}\right)H=0$
(15)

where $H(x)=HeunC(\alpha,\beta,\gamma,\delta,\eta,x)$ are the so-called
confluent Heun functions with

$\begin{split}\sigma&=\frac{1}{2}(\alpha-\beta-\gamma-\alpha\beta-\beta\gamma)-\eta\\\
\tau&=\frac{1}{2}(\alpha+\beta+\gamma+\alpha\gamma+\beta\gamma)+\delta+\eta\end{split}$
(16)

By means of an integrating factor, it is possible to relate the equations (13)
and (15), and thus each of their coefficients,

$\begin{split}\alpha&=\pm 4M\sqrt{\mu^{2}-\omega^{2}}\\\ \beta&=0\\\
\gamma&=\pm 4\sqrt{\frac{m\omega J}{M}-M^{2}\omega^{2}}\\\
\delta&=8M^{2}(2\omega^{2}-\mu^{2})\\\ \eta&=\frac{4m\omega J}{M}-l(l+1)\
.\end{split}$ (17)

With these expressions we determined the exact analytical solution of the
Klein-Gordon equation in the Hartle-Thorne metric given by eq. (8). By
observing that some terms of eq. (12) are still very small, we may obtain a
simpler solution with the purpose to have a better understanding of the
physics of this system relating it to more common special functions,

$\frac{d^{2}u}{dr^{2}}+\left[k^{2}+\frac{a^{\prime}}{r^{2}}+\frac{b^{\prime}}{r}+\frac{c^{\prime}}{(r-2M)^{2}}+\frac{d^{\prime}}{r-2M}\right]u=0,$
(18)

where

$\begin{split}a^{\prime}&=\frac{1}{4}\\\
b^{\prime}&=\frac{1}{4M}+\frac{l(l+1)}{2M}-\frac{Jm\omega}{M^{2}}\\\
c^{\prime}&=\frac{1}{4}+4M^{2}\omega^{2}-\frac{2Jm\omega}{M}\\\
d^{\prime}&=-\frac{1}{4M}-2M\mu^{2}-\frac{l(l+1)}{2M}+4M\omega^{2}+\frac{Jm\omega}{M^{2}}\\\
k^{2}&=\omega^{2}-\mu^{2}.\end{split}$ (19)

Observing that $a^{\prime}/(r-2M)^{2}$ and $b^{\prime}/(r-2M)$, for fixed
parameters and $r$ values, are, numerically, several orders of magnitude
smaller than the other three terms that multiply the function $u(r)$ directly
in eq. (18), a reasonable approximation is to neglect these terms. By
performing a change of variables $z=2ik(r-2M)$, the equation becomes

$u^{\prime\prime}(z)+\left[-\frac{1}{4}+\frac{1/4-\lambda^{2}}{z^{2}}+\frac{\varepsilon}{z}\right]u(z)=0,$
(20)

where $\varepsilon=d^{\prime}/(2ik)$ e $\lambda=\sqrt{1/4-c^{\prime}}$. The
equation above is the so-called Whittaker equation, it is a modified version
of the confluent hypergeometric equation and it can be written in the form of
the Kummer’s equation through a transformation of the function. This equation
possesses two solutions, the Whittaker functions $M_{\varepsilon,\lambda}(z)$
and $W_{\varepsilon,\lambda}(z)$, which can be written, respectively, in terms
of a Kummer function $M(A,B,z)$ and a Tricomi function $U(A,B,z)$ that are
solutions to the confluent hypergeometric equation Abramowitz and Stegun
(1964),

$\displaystyle M_{\varepsilon,\lambda}(z)$
$\displaystyle=e^{-\frac{1}{2}z}z^{\frac{1}{2}+\lambda}M(\frac{1}{2}+\lambda-\varepsilon,1+2\lambda,z),$
(21) $\displaystyle W_{\varepsilon,\lambda}(z)$
$\displaystyle=e^{-\frac{1}{2}z}z^{\frac{1}{2}+\lambda}U(\frac{1}{2}+\lambda-\varepsilon,1+2\lambda,z).$
(22)

As far as the term proportional to $c^{\prime}$ is still much smaller than the
remaining terms of eq. (18) and in the region exterior to the stars $r>>2M$ we
have $r\simeq r-2m$ and then the equation becomes, with a good accuracy

$\frac{d^{2}u}{dr^{2}}+\left[k^{2}+\frac{d^{\prime}}{r}\right]u=0.$ (23)

Supposing a solution with the form

$u(r)=re^{ikr}w(r)$ (24)

and a change of variables $z=2ikr$, we obtain the equation

$zw^{\prime\prime}(z)+(2-z)w^{\prime}(z)-\left(1-\frac{d^{\prime}}{2ik}\right)w(z)=0,$
(25)

that has the solution

$w(z)=c_{1}\,M\left(1-\frac{d^{\prime}}{2ik},2,z\right)+c_{2}\,U\left(1-\frac{d^{\prime}}{2ik},2,z\right),$
(26)

and for the physical problem that we are studying $c_{2}=0$. The Kummer
function may defined as

$M\left(1-\frac{d^{\prime}}{2ik},2,z\right)=\sum_{s=0}^{\infty}\frac{\left(1-\frac{d^{\prime}}{2ik}\right)_{s}}{(2)_{s}\,s!}z^{s},$
(27)

that asymptotically has the behavior of a Bessel function, then for large
values of $z$ (that is exactly the region where this calculation is
considered), the solution may be written in the form

$R(r)=\frac{e^{i\left[kr+\frac{d^{\prime}}{2k}\ln{r}\right]}}{r},$ (28)

or

$R(r)=\Psi(r)e^{i\frac{jm\omega}{2k}\ln r}=\Psi(r)e^{i\delta_{R}}$ (29)

that shows explicitly the phase originated by the effect of the rotation of
the star in the lowest order, $\delta_{R}$, that may be interpreted as a
rotational phase. In the next section we will investigate the solution for
higher orders of $j_{1}$.

##  IV Expansion of the Hartle-Thorne metric to second order in J

As far as some mistakes have been observed in the Hartle-Thorne metric given
by eq. (3), if we are interested in studying terms of the order of $j_{1}^{2}$
and higher in the solution we have to use the Berti et al. metric (5) in our
calculations. It is interesting to note that if we keep the terms up to first
order in $j_{1}$ and $Q=0$, eq. (7) is recovered, that means that the results
presented in the last section are correct and the differences occur for higher
terms in $j_{1}$. So, in this section we will derive an equation taking into
account the terms that appear in the corrected metric and write down a
solution based upon an infinite series ansatz. In order to observe the
behavior of the corrections in the solution an easier way is to solve the
equations for fixed values of the angle $\theta$. Considering these
assumptions, the resulting radial equation will become

$\left[g^{rr}\partial_{r}^{2}+\tilde{g}\partial_{r}\left(-\omega^{2}g^{tt}-m^{2}g^{\phi\phi}+m\omega
g^{t\phi}-\mu^{2}\right)\right]f(r)=0$ (30)

where

$\tilde{g}=\partial_{r}g^{rr}+\frac{1}{2g}g^{rr}\partial_{r}g.$ (31)

After some manipulations the equation may be written in the form

$\left(M_{0}+\frac{M_{1}}{r}+\frac{M_{2}}{r^{2}}+\frac{M_{3}}{r^{3}}+\frac{M_{4}}{r^{4}}\right)g^{\prime\prime}(r)+\left(N_{0}+\frac{N_{1}}{r}+\frac{N_{2}}{r^{2}}+\frac{N_{3}}{r^{3}}+\frac{N_{4}}{r^{4}}\right)g(r)=0$
(32)

where

$f(r)=\frac{A}{r-2M}g(r)$ (33)

and

$\begin{split}M_{0}&=1\\\ M_{1}&=-6M\\\ M_{2}&=12M^{2}\\\ M_{3}&=-8M^{3}\\\
M_{4}&=2J^{2}(5-12u^{2})\\\ N_{0}&=\omega^{2}-\mu^{2}\\\
N_{1}&=4M\mu^{2}-2M\omega^{2}\\\ N_{2}&=m^{2}/(u^{2}-1)-4M^{2}\mu^{2}\\\
N_{3}&=4m^{2}M/(u^{2}-1)+2Jm\omega+2M\\\
N_{4}&=-4M^{2}+4m^{2}M^{2}/(u^{2}-1)-4JmM\omega+2J^{2}\omega^{2}(3u^{2}-2).\end{split}$
(34)

This equation may be solved by the method presented in Elizalde (1987) where
the solution has the form

$\displaystyle g(r)$ $\displaystyle=\alpha\mathrm{e}^{\pm iF(r)}$ (35)
$\displaystyle F(r)$
$\displaystyle=kr+a_{0}\ln{\frac{r}{2M}}+\sum_{n=1}^{\infty}a_{n}\left(\frac{2M}{r}\right)^{n}.$
(36)

and the exponent $F(r)$ is a sum of functions and a series of powers of $M/r$
where the $a_{n}$ coefficients have to be determined. It is easy to find the
first few coefficients by simply substituting this ansatz back into equation
(32), however, as $n$ increases the expressions become larger and larger. The
first five coefficients are given by

$\begin{split}k&=\sqrt{\frac{N_{0}}{M_{0}}}\\\
a_{0}&=\frac{-k^{2}M_{1}+N_{1}}{2kM_{0}}\\\
a_{1}&=\frac{1}{4MkM_{0}}\Bigg{(}a_{0}^{2}+2a_{0}kM_{1}+k^{2}M_{2}-N_{2}\pm
ia_{0}M_{0}\Bigg{)}\\\
a_{2}&=\frac{1}{4(2M)^{2}kM_{0}}\Bigg{(}a_{0}^{2}M_{1}+2a_{0}kM_{2}+k^{2}M_{3}-N_{3}-4Ma_{0}M_{0}a_{1}\\\
&-4MkM_{1}a_{1}\pm ia_{0}M_{1}\mp i4MM_{0}a_{1}\Bigg{)}\\\
a_{3}&=\frac{1}{6(2M)^{3}kM_{0}}\Bigg{(}a_{0}^{2}M_{2}+2a_{0}kM_{3}+k^{2}M_{4}-N4-4Ma_{0}M_{1}a_{1}\\\
&-4MkM^{2}a_{1}+4M^{2}M_{0}a_{1}^{2}-16M^{2}a_{0}M_{0}a_{2}-16M^{2}kM_{1}a_{2}+\\\
&\pm ia_{0}M_{2}\mp i4MM_{1}a_{1}\mp i24M^{2}M_{0}a_{2}\Bigg{)}.\end{split}$
(37)

For large values of $r$ we have

$R(r)\sim\frac{e^{i\left[kr+a_{0}\ln{r}\right]}}{r}.$ (38)

that is the same result obtained in Sec. III, in a very different way, shown
in eq.(28).

##  V Analyzing the solutions

In this section we will investigate the solution (35) for several physical
configurations and observe how it changes depending on the mass, angular
velocity and energy of the star or of the particle whose motion the Klein-
Gordon equation purports to describe. It is also interesting to study high
energy and low energy test particles. We will start by listing all physical
configurations we chose in order to verify the size of the mass and angular
momentum corrections. We selected three celestial bodies: the Sun Workman  _et
al._ (2022), the neutron star PSR B1257+12 Zheng  _et al._ (2015),Wolszczan
and Frail (1992), and the white dwarf PG 2131+066 Reed  _et al._ (2000). All
numerical values used in this section are given in Planck units, where
$c=G=\hbar=1$. For the sun we have

$\begin{split}M_{\odot}&=9.136\times 10^{37}\\\ J_{\odot}&=1.82\times
10^{75}\\\ R_{\odot}&=4.30835481\times 10^{43},\end{split}$ (39)

for the neutron star,

$\begin{split}M_{ns}&=1.4M_{\odot}\\\ J_{ns}&=0.00297J_{\odot}\\\
R_{ns}&=0.000015R_{\odot},\end{split}$ (40)

and for the white dwarf the corresponding values are

$\begin{split}M_{wd}&=0.608M_{\odot}\\\ J_{wd}&=0.1432J_{\odot}\\\
R_{wd}&=0.0186R_{\odot}.\end{split}$ (41)

We also chose two test particles with very different masses, the pion with
mass $\mu_{\pi}$=0.139 GeV Workman  _et al._ (2022), and the Higgs boson with
mass $\mu_{Higgs}$= 125.1 GeV, that in Planck units are

$\begin{split}\mu_{\pi}&=1.1056\times 10^{-20}\\\ \mu_{Higgs}&=1.0245\times
10^{-17}\end{split}$ (42)

We separate eq. (36) in four different parts: $\mathrm{A}$ is the real part of
$F$ containing only mass terms; $\mathrm{B}$ is the imaginary part of $F$
containing only mass terms; $\Gamma$ is the real part of $F$ containing the
angular momentum corrections, i.e., terms that are proportional to the angular
momentum $J$ or to its square $J^{2}$; $\Delta$ is the imaginary part of $F$
containing the angular momentum contributions.

$F(r)=\mathrm{A}(r)+\mathrm{B}(r)+\Gamma(r)+\Delta(r)$ (43)

These functions have the same structure of $F(r)$ given by eq. (36) and may be
calculated by the sum of five terms determined by different dependencies on
the variable $r$. In the Tables, each column is related with one kind of
dependence on the variable $r$. We will show the results for the six
combinations of stars and particles given by eq. (39)–(42) for particles with
energies $\omega=2\mu$, $r=R$, that means particles near the surfaces, $m=1$
and $\theta=\pi/4$. These results are displayed in Tab. I-VI.

| $r$ | $\ln{r}$ | $1/r$ | $1/r^{2}$ | $1/r^{3}$
---|---|---|---|---|---
$\mathrm{A}$ | $8.25031\times 10^{23}$ | $5.04993\times 10^{19}$ | $-1.71066\times 10^{13}$ | $-3.6348\times 10^{7}$ | $-1.02727\times 10^{2}$
$\mathrm{B}$ | 0 | 0 | $2.47395\times 10^{-6}$ | $4.24684\times 10^{-12}$ | $1.29492\times 10^{-17}$
$\Gamma$ | 0 | 0 | 0 | $-5.66093\times 10^{-13}$ | $-8.81302\times 10^{-2}$
$\Delta$ | 0 | 0 | 0 | 0 | $6.86148\times 10^{-37}$

Table 1: Celestial object: Sun. Particle: Pion.

| $r$ | $\ln{r}$ | $1/r$ | $1/r^{2}$ | $1/r^{3}$
---|---|---|---|---|---
$\mathrm{A}$ | $1.23755\times 10^{19}$ | $5.2965\times 10^{18}$ | $-2.23527\times 10^{18}$ | $-4.43284\times 10^{17}$ | $-1.16929\times 10^{17}$
$\mathrm{B}$ | 0 | 0 | $2.30902\times 10^{-1}$ | $3.69947\times 10^{-2}$ | $1.05281\times 10^{-2}$
$\Gamma$ | 0 | 0 | 0 | $-7.48238\times 10^{-6}$ | $-2.30951\times 10^{8}$
$\Delta$ | 0 | 0 | 0 | 0 | $6.04614\times 10^{-25}$

Table 2: Celestial object: Neutron Star. Particle: Pion.

| $r$ | $\ln{r}$ | $1/r$ | $1/r^{2}$ | $1/r^{3}$
---|---|---|---|---|---
$\mathrm{A}$ | $1.53456\times 10^{22}$ | $2.2049\times 10^{19}$ | $-3.39985\times 10^{14}$ | $-2.36138\times 10^{10}$ | $-2.18152\times 10^{6}$
$\mathrm{B}$ | 0 | 0 | $8.0869\times 10^{-5}$ | $4.53782\times 10^{-9}$ | $4.52286\times 10^{-13}$
$\Gamma$ | 0 | 0 | 0 | $-2.34386\times 10^{-10}$ | $-2.81012\times 10^{2}$
$\Delta$ | 0 | 0 | 0 | 0 | $1.52739\times 10^{-32}$

Table 3: Celestial object: White Dwarf. Particle: Pion.

| $r$ | $\ln{r}$ | $1/r$ | $1/r^{2}$ | $1/r^{3}$
---|---|---|---|---|---
$\mathrm{A}$ | $7.64512\times 10^{26}$ | $4.6795\times 10^{22}$ | $-1.58518\times 10^{16}$ | $-3.36818\times 10^{10}$ | $-9.51912\times 10^{4}$
$\mathrm{B}$ | 0 | 0 | $2.47395\times 10^{-6}$ | $4.24684\times 10^{-12}$ | $1.29492\times 10^{-17}$
$\Gamma$ | 0 | 0 | 0 | $-5.66093\times 10^{-13}$ | $-8.16655\times 10$
$\Delta$ | 0 | 0 | 0 | 0 | $7.40464\times 10^{-40}$

Table 4: Celestial object: Sun. Particle: Higgs boson.

| $r$ | $\ln{r}$ | $1/r$ | $1/r^{2}$ | $1/r^{3}$
---|---|---|---|---|---
$\mathrm{A}$ | $1.14677\times 10^{22}$ | $4.90798\times 10^{21}$ | $-2.0713\times 10^{21}$ | $-4.10768\times 10^{20}$ | $-1.08352\times 10^{20}$
$\mathrm{B}$ | 0 | 0 | $2.30902\times 10^{-1}$ | $3.69947\times 10^{-2}$ | $1.05281\times 10^{-2}$
$\Gamma$ | 0 | 0 | 0 | $-7.48238\times 10^{-6}$ | $-2.1401\times 10^{11}$
$\Delta$ | 0 | 0 | 0 | 0 | $6.52476\times 10^{-28}$

Table 5: Celestial object: Neutron star. Particle: Higgs boson.

| $r$ | $\ln{r}$ | $1/r$ | $1/r^{2}$ | $1/r^{3}$
---|---|---|---|---|---
$\mathrm{A}$ | $1.42199\times 10^{25}$ | 2.04316 $\times 10^{22}$ | $-3.15-45\times 10^{17}$ | $-2.18816\times 10^{13}$ | $-2.02149\times 10^{9}$
$\mathrm{B}$ | 0 | 0 | $8.0869\times 10^{-5}$ | $4.53782\times 10^{-9}$ | $4.52286\times 10^{-13}$
$\Gamma$ | 0 | 0 | 0 | $-2.34386\times 10^{-10}$ | $-2.60399\times 10^{2}$
$\Delta$ | 0 | 0 | 0 | 0 | $1.64829\times 10^{-35}$

Table 6: Celestial object: White dwarf. Particle: Higgs boson.

As we can see in the tables, the star to which the corrections are more
significant is the neutron star. However, even the largest angular momentum
corrections are small if compared with the other terms. For this reason we
will continue this analysis considering neutron stars. By varying the physical
quantities in question (such as the energy of the particle, and the angular
momentum of the star) we will further explore our solution to the radial
equation (32). We begin by defining three functions related to the solution
which we will plot and which might give us some information on its status as a
solution. These functions are

$\displaystyle f_{1}$ $\displaystyle=iF$ (44) $\displaystyle f_{2}$
$\displaystyle=\frac{e^{iF}}{r-2M}$ (45) $\displaystyle f_{3}$
$\displaystyle=\frac{\Gamma+\Delta}{|\mathrm{A}+\mathrm{B}+\Gamma+\Delta|}$
(46)

The function $f_{1}$ is precisely the exponent in the solution (35). The
function $f_{2}$ is eq. (35), along with the integrating factor defined in
(33); the end result is that $f_{2}$ represents the radial solution. The
function $f_{3}$ is defined as the difference between the exponent (36) and
the mass contributions as defined in (43), divided by the absolute value of
$F$. Since all that is left in the numerator of $f_{3}$ is the sum of angular
momentum contributions $\Gamma+\Delta$ as defined in (43), this function tells
us the approximate order of magnitude of the angular momentum contributions
when compared to the mass contributions.

Fig. 1 and Fig. 2 represent the function $f_{1}$ in two regions, in Fig. 2 we
show the results to the higher values of the energy. As we can see, the major
difference between figures 1 and 2 is that the value of the function
stabilizes as the energy increases, and even big differences in energy account
for little changes in $f_{1}$. For lower energies, even with small changes,
the value of $f_{1}$ changes considerably.

Figure 1: Absolute value of the real part of $f_{1}$ for some values of the
energy. Figure 2: Absolute value of the real part of $f_{1}$ for higher values
of the energy. Figure 3: Real part of $f_{2}$ as function of $r$. Figure 4:
Real part of $f_{2}$ as function of $r$.

In Fig. 3 and 4 we represent the function $f_{2}$, the radial wave function,
in two regions according to the values of the energies of the pions. As we can
see in Fig. 3, for large distances, far away from the surface of the star, the
curves seem to coalesce, and this effect is especially visible for higher
energies.

The last two graphs, Fig. 5 and 6, represent the function $f_{3}$, which shows
the size of the angular momentum contribution as a variation in its value is
done, while maintaining the star’s original mass and radius fixed. We observe
that the largest angular momentum selected, $J_{6}=10,000J_{ns}$, where
$J_{ns}$ is the angular momentum of the neutron star as given in (40),
probably cannot be reached for a real star as it exceeds the maximum possible
value given in the literature Krastev  _et al._ (2008) for the angular
velocity (and therefore angular momentum) of a neutron star. We used it
regardless in order to show what happens to the solution for extreme
rotations.

Figure 5: Imaginary part of $f_{3}$ as a function of $r$ for the actual
angular momentum $J_{ns}$ of the neutron star and for variations of $J$.
Figure 6: Imaginary part of $f_{3}$ as a function of $r$.

Fig. 5 and 6 show that the shape of the curves does not change considerably
when varying the angular momentum, however, the absolute value of the function
does. This is the expected result: the larger the angular momentum of the
neutron star, the larger will be its contribution to the solution. We can also
investigate what happens with variations of the energy $\omega$ in the case
where the angular momentum used was the neutron star’s own. The results are
shown in Fig. 7.

Figure 7: Absolute value of the imaginary part of $f_{3}$ as a function of $r$
for the different energies $\omega$ and fixed angular momentum $J=J_{ns}$.

As we can see in figure 7, the angular momentum contribution is larger for
higher energies. Indeed, our solution seems to take the angular momentum into
account more significantly when the star is rotating quickly (though still
under the “slow” rotation regime as dictated by the original paper) and the
particle is a high-energy zero-spin boson. We must remark that in the results
presented in this section, the velocity of the particle $v_{p}>>v_{R}$, where
$v_{R}$ is the rotation velocity of the star for $r=R$. If we have small
energies, $v_{p}\sim v_{R}$, more terms will be needed in eq. (36).

##  VI Conclusions

In this work we studied the behavior of spin-0 bosons near rotating stars. For
this purpose we considered the Klein-Gordon equation in the region external to
the star with the spacetime structure determined by the Hartle-Thorne metric
(8) up to first order in $j_{1}$, that is a very good approximation. The
corrections for the metric proposed in Berti  _et al._ (2005) appears for
terms of order of $j_{1}^{2}$, and in general are very small. By solving the
Klein-Gordon equation derived in the metric (8) a general solution have been
found, with the radial part given in terms of a confluent Heun function. By
analyzing the terms of the radial equation for $r>R$ it is possible to neglect
some of then and still have an accurate result as far as they are very small.
With this procedure we obtained a simpler solution that shows explicitly the
main effect of the rotation of the star in the wave function that is given by
the rotational phase defined from eq. (29)

$\delta_{R}=\frac{jm\omega}{2k}\ln r\ .$ (47)

Then an analysis for higher powers of $j_{1}$ have been performed considering
the Klein-Gordon equation in the Berti metric, eq. (5), and solving it by the
method presented in Elizalde (1987). We presented the calculations up to
$a_{3}$ in eq. (36), and the numerical calculations have shown that for bosons
with relativistic energies the solutions converge easily and the effect of the
rotation of the stars is very small. Among the considered physical systems,
the contribution of the terms depending on $j_{1}$ is more important for the
neutron stars. For lower values of the energies ($v_{p}\sim v_{R}$) this
effect increases and more terms are needed in the solution. So, from the
results we can observe that an important aspect of this effect is that it
increases as the rotation velocity of the star increases relatively to the
velocity of the test particle. For these reasons, for stars this effect must
always be small as far as $v_{R}$ cannot reach high values and we expect that
near black holes or galaxies it becomes more important. This kind of behavior
is similar to the pions and to the Higgs bosons, that is for light or heavy
bosons.

The inclusion of mass asymmetries of the stars is straightforward, we just
have to take $Q\neq 0$ in eq. (7) and follow the procedure presented in this
work, and then, it is possible to find this effect in the wave function.

If we think in terms of experimental results, the best way is to study this
effect is to look for observables for which the rotational phase (47) is
important, as for example in the analysis of the final state interactions of
decays or in the study of interacting particles. We must also remark that even
if this kind of effect is small, an experimental verification is a proof of
the frame dragging by rotating objects in and consequently of the general
theory of relativity influence in quantum systems.

##  VII Acknowledgments

We would like to thank CNPq for the financial support.

## References

  * Schiff (1960) L. I. Schiff, Physical Review Letters 4, 215 (1960).
  * Lense and Thirring (1918) J. Lense and H. Thirring, Physikalische Zeitschrift 19, 156 (1918).
  * Everitt  _et al._ (2011) C. F. Everitt, D. DeBra, B. Parkinson, J. Turneaure, J. Conklin, M. Heifetz, G. Keiser, A. Silbergleit, T. Holmes, J. Kolodziejczak, _et al._ , Physical Review Letters 106, 221101 (2011).
  * Abbot  _et al._ (2016) B. P. Abbot _et al._ , Physical Review Letters 116 (2016), 10.1103/physrevlett.116.061102.
  * Parker (1980) L. Parker, Phys. Rev. Lett. 44, 1559 (1980).
  * Ahmed (2022) F. Ahmed, Int. J. Mod. Phys. A 37, 2250186 (2022).
  * Ahmed (2023) F. Ahmed, Commun. Theor. Phys. 75, 025202 (2023).
  * Santos  _et al._ (2019) L. C. N. Santos, C. E. Mota, and C. C. Barros Jr., Adv. High Energy Phys. 2019, 2729352 (2019), arXiv:1901.04104 [hep-th] .
  * Soares  _et al._ (2021) A. R. Soares, R. L. L. Vitória, and H. Aounallah, Eur. Phys. J. Plus 136, 966 (2021).
  * Yang  _et al._ (2021) Y. Yang, Z.-W. Long, Q.-K. Ran, H. Chen, Z.-L. Zhao, and C.-Y. Long, Int. J. Mod. Phys. A 36, 2150023 (2021), arXiv:2106.12225 [quant-ph] .
  * Rouabhia  _et al._ (2023) T. I. Rouabhia, A. Boumali, and H. Hassanabadi, Phys. Part. Nucl. Lett. 20, 112 (2023).
  * Santos and Barros Jr. (2016) L. C. N. Santos and C. C. Barros Jr., Eur. Phys. J. C 76, 560 (2016), arXiv:1508.07307 [hep-th] .
  * Santos and Barros Jr. (2017) L. C. N. Santos and C. C. Barros Jr., Eur. Phys. J. C 77, 186 (2017), arXiv:1608.06692 [hep-th] .
  * Santos and Barros Jr. (2018a) L. C. N. Santos and C. C. Barros Jr., Eur. Phys. J. C 78, 13 (2018a), arXiv:1801.01024 [hep-th] .
  * Vitória and Bakke (2018a) R. L. L. Vitória and K. Bakke, Eur. Phys. J. C 78, 175 (2018a), arXiv:1802.07536 [gr-qc] .
  * Bezerra  _et al._ (2016) V. B. Bezerra, M. S. Cunha, L. F. F. Freitas, C. R. Muniz, and M. O. Tahim, Modern Physics Letters A 32, 1750005 (2016).
  * Santos and Barros Jr. (2018b) L. C. N. Santos and C. C. Barros Jr., Int. J. Mod. Phys. A 33, 1850122 (2018b).
  * Sedaghatnia  _et al._ (2019) P. Sedaghatnia, H. Hassanabadi, and F. Ahmed, Eur. Phys. J. C 79, 541 (2019).
  * Guvendi  _et al._ (2022) A. Guvendi, S. Zare, and H. Hassanabadi, Phys. Dark Univ. 38, 101133 (2022).
  * Vitória and Bakke (2018b) R. L. L. Vitória and K. Bakke, Eur. Phys. J. Plus 133, 490 (2018b).
  * Barros Jr. (2005) C. C. Barros Jr., Eur. Phys. J. C 42, 119 (2005), arXiv:physics/0409064 .
  * Barros Jr. (2006) C. C. Barros Jr., Eur. Phys. J. C 45, 421 (2006), arXiv:hep-ph/0504179 .
  * Hartle (1967) J. B. Hartle, Astrophys. J. 150, 1005 (1967).
  * Hartle and Thorne (1968) J. B. Hartle and K. S. Thorne, Astrophys. J. 153, 807 (1968).
  * Berti  _et al._ (2005) E. Berti, F. White, A. Maniopoulou, and M. Bruni, Monthly Notices of the Royal Astronomical Society 358, 923 (2005).
  * Krastev  _et al._ (2008) P. G. Krastev, B.-A. Li, and A. Worley, The Astrophysical Journal 676, 1170 (2008).
  * Abramowicz  _et al._ (2003) M. A. Abramowicz, G. J. E. Almergren, W. Kluzniak, and A. V. Thampan, “The hartle-thorne circular geodesics,” (2003), arXiv:gr-qc/0312070 [gr-qc] .
  * Birrell and Davies (1982) N. D. Birrell and P. C. W. Davies,  _Quantum Fields in Curved Space_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1982).
  * Parker and Toms (2009) L. E. Parker and D. Toms,  _Quantum Field Theory in Curved Spacetime: Quantized Field and Gravity_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2009).
  * Abramowitz and Stegun (1964) M. Abramowitz and I. Stegun,  _Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables_, Applied mathematics series (Dover Publications, 1964).
  * Elizalde (1987) E. Elizalde, Physical review D: Particles and fields 36, 1269 (1987).
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Zheng  _et al._ (2015) H. Zheng, K.-J. Sun, and L.-W. Chen, The Astrophysical Journal 800, 141 (2015).
  * Wolszczan and Frail (1992) A. Wolszczan and D. A. Frail, Nature 355, 145 (1992).
  * Reed  _et al._ (2000) M. D. Reed, S. D. Kawaler, and M. S. O’Brien, The Astrophysical Journal 545, 429 (2000).


#  Asymptotic distribution for the proportional covariance model

Myung Geun Kim111 E-mail<EMAIL_ADDRESS>
###### Abstract

Asymptotic distribution for the proportional covariance model under
multivariate normal distributions is derived. To this end, the parametrization
of the common covariance matrix by its Cholesky root is adopted. The
derivations are made in three steps. First, the asymptotic distribution of the
maximum likelihood estimators of the proportionality coefficients and the
Cholesky inverse root of the common covariance matrix is derived by finding
the information matrix and its inverse. Next, the asymptotic distributions for
the case of the Cholesky root of the common covariance matrix and finally for
the case of the common covariance matrix itself are derived using the
multivariate $\delta$-method. As an application of the asymptotic distribution
derived here, a hypothesis for homogeneity of covariance matrices is
considered.

Keywords: Asymptotic distribution, Cholesky decomposition, maximum likelihood
estimators, proportionality of covariance matrices.

##  1 Introduction

Let ${\bf X}_{k}$ be a generic $p$-variate random vector having a multivariate
normal distribution $N_{p}(\mbox{\boldmath$\mu$}_{k},{\bf\Sigma}_{k})$ for
each $k=1,...,K,$ where the $\mbox{\boldmath$\mu$}_{k}$ are column vectors of
the $p$-dimensional Euclidean space and the ${\bf\Sigma}_{k}$ are $p$ by $p$
positive definite covariance matrices. Let the ${\bf S}_{k}$ be the mutually
independent unbiased estimators of the ${\bf\Sigma}_{k}$, each for a sample of
size $N_{k}$ and let $n_{k}=N_{k}-1$. Then the $n_{k}{\bf S}_{k}$ are
independently distributed according to Wishart distribution
$W_{p}({\bf\Sigma}_{k},n_{k})$.

The hypothesis that $K$ covariance matrices are proportional to each other can
be expressed as

$H_{p}:{\bf\Sigma}_{k}=c_{k}{\bf\Sigma}_{1}\mbox{~{}~{}~{}~{}for all k = 2,
... ,K}$ (1)

where the $c_{k}$ are unknown positive proportionality coefficients. In the
two sample case, estimation for proportional covariance models was studied by
Khatri(1967), Guttman, Kim and Olkin (1985) and Rao (1983). The proportional
covariance model was adopted in classification by 0wen (1984). An algorithm to
find the maximum likelihood estimators was given independently by 0wen (1984),
Manly and Rainer (1987) and Eriksen (1987). The convergence of the algorithm
and the uniqueness of maximum likelihood estimates were proved by Eriksen
(1987). Using a parametrization of the spectral decomposition of the common
covariance matrix, Flury (1986) considered maximum likelihood estimation of
the proportionality coefficients and the common covariance matrix, and
suggested an iterative method for performing computations. Under the same
parametrization as in Flury (1986), Boente et al. (2007, 2009) studied
influence diagnostics and robust estimation for proportional covariance
models.

The Cholesky decomposition theorem states that there exists a lower triangular
matrix ${\bf A}$ with positive diagonal elements such that

${\bf\Sigma}_{1}={\bf A}{\bf A}^{T},$

and this expression is unique. We use $T$ to denote a transpose of a matrix.
Then the model (1) is parametrized as

${\bf\Sigma}_{k}=c_{k}{\bf A}{\bf A}^{T}~{}~{}~{}(1\leq k\leq K).$ (2)

For notational convenience, it is assumed that $c_{1}$ is identical with 1.
The parametrization (2) can be considered as the special case of the model
treated by Jöreskog (1971) and the easily accessible LISREL (linear structural
relationships) program is available for getting the maximum likelihood
estimates of the model parameters (Jöreskog and Sörbom, 1996).

Let ${\bf B}^{T}={\bf A}^{-1}$. The matrix ${\bf B}$ is then an upper
triangular matrix with positive diagonal elements. Then the parametrization
(2) can be written as

${\bf\Sigma}_{k}^{-1}=\frac{1}{c_{k}}{\bf B}{\bf B}^{T}~{}~{}~{}(1\leq k\leq
K).$ (3)

In this work, we derive the asymptotic distribution of the maximum likelihood
estimators of the proportionality coefficients and the common covariance
matrix, using the well-known properties of the maximum likelihood estimators
under regularity conditions (Rao, 1973, Chapter 6). The procedure for doing so
is as follows. Firstly, the information matrix for the parameters $c_{k}$ and
${\bf B}$ in (3) is derived. Secondly, the asymptotic covariance matrix for
the maximum likelihood estimators of the parameters $c_{k}$ and ${\bf B}$ in
(2) is found by inverting the corresponding information matrix. Thirdly, we
derive the asymptotic covariance matrix for the maximum likelihood estimators
of the parameters $c_{k}$ and ${\bf A}$ in (2), using the multivariate
$\delta$-method. Finally, the asymptotic covariance matrix for the maximum
likelihood estimators of the proportionality coefficients and the common
covariance matrix is derived using the multivariate $\delta$-method. Finally,
a hypothesis for homogeneity of covariance matrices is considered as an
application of the asymptotic distribution derived here.

##  2 The information matrix for the parameters $c_{k}$ and ${\bf B}$

In this section, we derive the information matrix for the $c_{k}$ and ${\bf
B}$. The likelihood function of the ${\bf\Sigma}_{k}$ given the ${\bf S}_{k}$
is

$L({\bf\Sigma}_{1},...,{\bf\Sigma}_{K})=C\times\prod_{k=1}^{K}[\mid{\bf\Sigma}_{k}\mid^{-\frac{n_{k}}{2}}\mbox{exp}\\{-\frac{n_{k}}{2}\mbox{tr}({\bf\Sigma}_{k}^{-1}{\bf
S}_{k})\\}],$

where $C$ is a constant not depending on the parameters. Since we are
interested in covariance matrices and a covariance matrix is location
invariant, there is no restriction on mean vectors.

Under the parametrization (3), we have

$\displaystyle\mid{\bf\Sigma}_{k}^{-1}\mid$ $\displaystyle=$ $\displaystyle
c_{k}^{-p}\prod_{i=1}^{p}b_{ii}^{2}$
$\displaystyle\mbox{tr}({\bf\Sigma}_{k}^{-1}{\bf S}_{k})$ $\displaystyle=$
$\displaystyle\frac{1}{c_{k}}\sum_{i=1}^{p}{\bf b}_{i}^{T}{\bf S}_{k}{\bf
b}_{i},$

where $b_{ij}$ denotes the $(i,j)$th element of ${\bf B}$ and ${\bf b}_{i}$
the $i$th column of ${\bf B}$. Thus the log-likelihood function, ignoring
constant term, becomes

$l({\bf c},{\bf
B})=n_{+}\sum_{i=1}^{p}\mbox{log}(b_{ii})-\sum_{k=1}^{K}\frac{n_{k}}{2}\\{p\mbox{log}(c_{k})+\frac{1}{c_{k}}\sum_{i=1}^{p}{\bf
b}_{i}^{T}{\bf S}_{k}{\bf b}_{i}\\},$

where $n_{+}=\sum_{k=1}^{K}n_{k}$ and ${\bf c}=(c_{2},...,c_{K})^{T}$. It is
easily shown that all the expectations of first derivatives of $l({\bf c},{\bf
B})$ with respect to the parameters are zero, and the expectations of second
derivatives of the log-likelihood function are as follows:

$\displaystyle-E\\{\frac{\partial^{2}l({\bf c},{\bf B})}{\partial
c_{k}^{2}}\\}$ $\displaystyle=$
$\displaystyle\frac{pn_{k}}{2c_{k}^{2}}~{}~{}~{}~{}(2\leq k\leq K)$
$\displaystyle-E\\{\frac{\partial^{2}l({\bf c},{\bf B})}{\partial
b_{ii}^{2}}\\}$ $\displaystyle=$ $\displaystyle n_{+}(a_{ii}^{2}+{\bf
a}_{(i)}^{T}{\bf a}_{(i)})~{}~{}~{}~{}(1\leq i\leq p)$
$\displaystyle-E\\{\frac{\partial^{2}l({\bf c},{\bf B})}{\partial
b_{ii}\partial c_{k}}\\}$ $\displaystyle=$
$\displaystyle-\frac{n_{k}a_{ii}}{c_{k}}~{}~{}~{}~{}(2\leq k\leq K,~{}1\leq
i\leq p)$ $\displaystyle-E\\{\frac{\partial^{2}l({\bf c},{\bf B})}{\partial
b_{ji}\partial c_{k}}\\}$ $\displaystyle=$ $\displaystyle 0~{}~{}~{}~{}(2\leq
k\leq K,~{}j<i)$ $\displaystyle-E\\{\frac{\partial^{2}l({\bf c},{\bf
B})}{\partial b_{ji}\partial b_{j^{\prime}i}}\\}$ $\displaystyle=$
$\displaystyle n_{+}{\bf a}_{(j)}^{T}{\bf a}_{(j^{\prime})}~{}~{}~{}~{}(1\leq
j,j^{\prime}<i),$

where $a_{ij}$ denotes the $(i,j)$th element of ${\bf A}$ and ${\bf a}_{(i)}$
the column vector formed by the elements in the $i$th row of ${\bf A}$, and
all the other cases have zero expectations.

Let ${\bf I}({\bf c},{\bf B})$ be the information matrix for the parameters
${\bf c}$ and ${\bf B}$, and let it be partitioned as

${\bf I}({\bf c},{\bf B})=\pmatrix{{\bf I}_{11}({\bf c},{\bf B})&{\bf
I}_{12}({\bf c},{\bf B})\cr{\bf I}_{21}({\bf c},{\bf B})&{\bf I}_{22}({\bf
c},{\bf B})\cr},$

such that ${\bf I}_{11}({\bf c},{\bf B})$ is the information matrix for ${\bf
c}$ and ${\bf I}_{22}({\bf c},{\bf B})$ for the nonzero elements of ${\bf B}$.
Let $r_{k}=n_{k}/n_{+}.$ Then ${\bf I}_{11}({\bf c},{\bf B})$ is a diagonal
matrix of order $K-1$:

$\frac{p}{2}\mbox{diag}(\frac{r_{2}}{c_{2}^{2}},...,\frac{r_{K}}{c_{K}^{2}}).$

We are interested in nonzero elements of ${\bf B}$, $b_{ij}~{}(i<j)$ and the
number of nonzero elements of ${\bf B}$ is $p(p+1)/2$. Let ${\bf b}_{i.j}$
$=(b_{1i},...,b_{ji})^{T}.$ Then we can notice that ${\bf b}_{i.1}=b_{1i}$ and
${\bf b}_{i.p}={\bf b}_{i}.$ Let $\alpha$
$=(r_{2}/c_{2},...,r_{K}/c_{K})^{T}.$ We write as ${\bf 1}_{p.i}$ the unit
column vector of dimension $p$ having one in the $i$th position and zeroes in
the other positions. Note that ${\bf 1}_{1.1}=1.$ Then the $(K-1)$ by
$p(p+1)/2$ matrix ${\bf I}_{12}({\bf c},{\bf B})$ is

$\bordermatrix{&{\bf b}_{1.1}^{T}&{\bf b}_{2.2}^{T}&\ldots&{\bf
b}_{p.p}^{T}\cr{\bf
c}&-a_{11}\mbox{\boldmath$\alpha$}&-a_{22}\mbox{\boldmath$\alpha$}{\bf
1}_{2.2}^{T}&\ldots&-a_{pp}\mbox{\boldmath$\alpha$}{\bf 1}_{p.p}^{T}}.$

Let ${\bf A}_{i}$ be the leading principal submatrix of ${\bf A}$ having order
$i$. Note that ${\bf A}_{1}=a_{11}$ and ${\bf A}_{p}={\bf A}.$ Then the
$p(p+1)/2$ by $p(p+1)/2$ matrix ${\bf I}_{22}({\bf c},{\bf B})$ is a block
diagonal matrix whose $i$th block element for ${\bf b}_{i.i}$ is given by

$a_{ii}^{2}{\bf 1}_{i.i}{\bf 1}_{i.i}^{T}+{\bf A}_{i}{\bf A}_{i}^{T}.$

This $i$th block element is of order $i$.

##  3 Asymptotic distribution of the maximum likelihood estimators of the
$c_{k}$ and ${\bf B}$

The asymptotic distribution of the maximum likelihood estimators of ${\bf c}$
and ${\bf B}$ is characterized by its asymptotic covariance matrix, and the
asymptotic covariance matrix is just the inverse of the corresponding
information matrix.

Let ${\bf V}({\bf c},{\bf B})$ be the asymptotic covariance matrix for the
maximum likelihood estimators of the parameters ${\bf c}$ and the nonzero
elements of ${\bf B}$, and let it be partitioned as

${\bf V}({\bf c},{\bf B})=\pmatrix{{\bf V}_{11}({\bf c},{\bf B})&{\bf
V}_{12}({\bf c},{\bf B})\cr{\bf V}_{21}({\bf c},{\bf B})&{\bf V}_{22}({\bf
c},{\bf B})\cr},$

such that ${\bf V}_{11}({\bf c},{\bf B})$ is the covariance matrix for the
maximum likelihood estimator of ${\bf c}$ and ${\bf V}_{22}({\bf c},{\bf B})$
for the maximum likelihood estimators of the nonzero elements of ${\bf B}$.

The following lemma gives the inverse of ${\bf I}_{22}({\bf c},{\bf B})$ and
enables us to get ${\bf V}_{11}({\bf c},{\bf B})$,

Lemma 3.1. For each $i=1,...,p,$ the inverse of $a_{ii}^{2}{\bf 1}_{i.i}{\bf
1}_{i.i}^{T}+$ ${\bf A}_{i}{\bf A}_{i}^{T}$ is

${\bf B}_{i}{\bf B}_{i}^{T}-\frac{1}{2}{\bf b}_{i.i}{\bf b}_{i.i}^{T}$

where ${\bf B}_{i}$ is the leading principal submatrix of ${\bf B}$ having
dimension $i$.

Lemma 3.1 is easily proved using the result (Sherman-Morrison formula) by
Bartlett (1951) and the fact that ${\bf B}_{i}^{T}={\bf A}_{i}^{-1}.$

Lemma 3.2. The matrix ${\bf I}_{11}({\bf c},{\bf B})~{}-$ ${\bf I}_{12}({\bf
c},{\bf B}){\bf I}_{22}({\bf c},{\bf B})^{-1}{\bf I}_{21}({\bf c},{\bf B})$,
${\bf U}_{11}({\bf c},{\bf B})$ say, of dimension $K-1$ is given by

$\frac{p}{2}\\{\mbox{diag}(\frac{r_{2}}{c_{2}^{2}},...,\frac{r_{K}}{c_{K}^{2}})-\mbox{\boldmath$\alpha$}\mbox{\boldmath$\alpha$}^{T}\\}.$

Proof of Lemma 3.2. The $i$th block element of ${\bf I}_{12}({\bf c},{\bf
B}){\bf I}_{22}({\bf c},{\bf B})^{-1}$ is a $(K-1)$ by $i$ matrix for the pair
${\bf c}$ and ${\bf b}_{i.i}$ and is given by

$-\frac{1}{2}\mbox{\boldmath$\alpha$}{\bf b}_{i.i}^{T}.$

Then ${\bf I}_{12}({\bf c},{\bf B}){\bf I}_{22}({\bf c},{\bf B})^{-1}{\bf
I}_{21}({\bf c},{\bf B})$ is given by

$\frac{p}{2}\mbox{\boldmath$\alpha$}\mbox{\boldmath$\alpha$}^{T}$

and therefore the proof is completed.

The inverse of ${\bf U}_{11}({\bf c},{\bf B})$ is just ${\bf V}_{11}({\bf
c},{\bf B})$ and is stated in the following lemma.

Lemma 3.3. The $(K-1)$ by $(K-1)$ matrix ${\bf V}_{11}({\bf c},{\bf B})$ is
given by

$\frac{2}{p}\\{\mbox{diag}(\frac{c_{2}^{2}}{r_{2}},...,\frac{c_{K}^{2}}{r_{K}})+\frac{1}{r_{1}}{\bf
c}{\bf c}^{T}\\}.$

The matrix ${\bf V}_{11}({\bf c},{\bf B})$ is easily derived using the result
by Bartlett (1951).

The matrix ${\bf I}_{12}({\bf c},{\bf B}){\bf I}_{22}({\bf c},{\bf B})^{-1}$
is derived in the proof of Lemma 3.2 and ${\bf V}_{12}({\bf c},{\bf B})~{}=$
$-{\bf V}_{11}({\bf c},{\bf B}){\bf I}_{12}({\bf c},{\bf B}){\bf I}_{22}({\bf
c},{\bf B})^{-1}$. Hence we get easily the following lemma.

Lemma 3.4. The $(K-1)$ by $p(p+1)/2$ matrix ${\bf V}_{12}({\bf c},{\bf B})$ is
given by the following matrix multiplied by $(1/pr_{1})$:

$\bordermatrix{&{\bf b}_{1.1}^{T}&{\bf b}_{2.2}^{T}&\ldots&{\bf
b}_{p.p}^{T}\cr{\bf c}&{\bf c}{\bf b}_{1.1}^{T}&{\bf c}{\bf
b}_{2.2}^{T}&\ldots&{\bf c}{\bf b}_{p.p}^{T}}.$

Since ${\bf V}_{22}({\bf c},{\bf B})~{}=$ ${\bf I}_{22}({\bf c},{\bf
B})^{-1}-$ ${\bf I}_{22}({\bf c},{\bf B})^{-1}$${\bf I}_{21}({\bf c},{\bf
B})$${\bf V}_{12}({\bf c},{\bf B})$, Lemma 3.1, the proof of Lemma 3.2 and
Lemma 3.4 yield ${\bf V}_{22}({\bf c},{\bf B})$.

Lemma 3.5. The $p(p+1)/2$ by $p(p+1)/2$ matrix ${\bf V}_{22}({\bf c},{\bf B})$
is given by

$\bordermatrix{&{\bf b}_{1.1}^{T}&{\bf b}_{2.2}^{T}&\ldots&{\bf
b}_{p.p}^{T}\cr{\bf b}_{1.1}&{\bf B}_{1}{\bf B}_{1}^{T}+d{\bf b}_{1.1}{\bf
b}_{1.1}^{T}&e{\bf b}_{1.1}{\bf b}_{2.2}^{T}&\ldots&e{\bf b}_{1.1}{\bf
b}_{p.p}^{T}\cr{\bf b}_{2.2}&e{\bf b}_{2.2}{\bf b}_{1.1}^{T}&{\bf B}_{2}{\bf
B}_{2}^{T}+d{\bf b}_{2.2}{\bf b}_{2.2}^{T}&\ldots&e{\bf b}_{2.2}{\bf
b}_{p.p}^{T}\cr\vdots&\vdots&\vdots&\ddots&\vdots\cr{\bf b}_{p.p}&e{\bf
b}_{p.p}{\bf b}_{1.1}^{T}&e{\bf b}_{p.p}{\bf b}_{2.2}^{T}&\ldots&{\bf
B}_{p}{\bf B}_{p}^{T}+d{\bf b}_{p.p}{\bf b}_{p.p}^{T}},$

where $d=\\{1-(1+p)r_{1}\\}/2pr_{1}$ and $e=(1-r_{1})/2pr_{1}$.

We use hereafter the hat notation to denote the maximum likelihood estimator.

The general properties of the maximum likelihood estimators yield the
asymptotic distribution of $\hat{{\bf c}}$ and $\hat{{\bf B}}$ under the
proportional covariance model.

Theorem 3.1. As $n_{+}$ tends to infinity with each $r_{k}$ held fixed, the
asymptotic distribution of

$\sqrt{n_{+}}(\hat{{\bf c}}-{\bf c},\hat{{\bf B}}-{\bf B})$

is a multivariate normal with mean zero and covariance matrix ${\bf V}({\bf
c},{\bf B})$.

Lemma 3.3 shows that the asymptotic distribution of
$\sqrt{n_{+}}(\hat{c}_{2}-c_{2},...,\hat{c}_{K}-c_{K})$ is free of
${\bf\Sigma}_{1}$. Guttman et al. (1985) derived the asymptotic distribution
of $\sqrt{n_{+}}(\hat{c}_{2}-c_{2})$ when $K=2$ and Lemma 3.3 reduces to
Theorem 4.1 of Guttman et al. (1985).

When $r_{1}=1$, that is, the single-population case, the coefficients $d$ and
$e$ become $-1/2$ and zero, respectively. Hence Lemma 3.5 shows that the
$\hat{{\bf b}}_{i.i}$ $(i=1,...,p)$ are asymptotically independent. By the
invariance property of the maximum likelihood estimator, the Cholesky inverse
root of ${\bf S}_{1}$ is just $\hat{{\bf B}}$, that is, ${\bf
S}_{1}^{-1}=\hat{{\bf B}}\hat{{\bf B}}^{T}.$ Thus Lemma 3.5 gives the
asymptotic distribution of the Cholesky inverse root of the sample covariance
matrix.

##  4 Asymptotic distribution of the maximum likelihood estimators of the
$c_{k}$ and ${\bf A}$

The asymptotic distribution of the maximum likelihood estimators of the
$c_{k}$ and ${\bf A}$ is derived from that of the $c_{k}$ and ${\bf B}$, using
the multivariate $\delta$-method.

Let ${\bf V}({\bf c},{\bf A})$ be the asymptotic covariance matrix for the
maximum likelihood estimators of the parameters ${\bf c}$ and ${\bf A}$. We
write ${\bf V}({\bf c},{\bf A})$ as a partitioned form:

${\bf V}({\bf c},{\bf A})=\pmatrix{{\bf V}_{11}({\bf c},{\bf A})&{\bf
V}_{12}({\bf c},{\bf A})\cr{\bf V}_{21}({\bf c},{\bf A})&{\bf V}_{22}({\bf
c},{\bf A})\cr},$

such that ${\bf V}_{11}({\bf c},{\bf A})$ is the covariance matrix for the
maximum likelihood estimator of ${\bf c}$ and ${\bf V}_{22}({\bf c},{\bf A})$
for the maximum likelihood estimators of the nonzero elements of ${\bf A}$.

Since ${\bf B}^{T}={\bf A}^{-1}$, we have the following lemma.

Lemma 4.1. The partial derivative of $a_{ji}$ with respect to $b_{hg}$ is

$\frac{\partial a_{ji}}{\partial
b_{hg}}=\left\\{\begin{array}[]{cl}-a_{jg}a_{hi}&(i\leq h\leq g\leq j)\\\
0&\mbox{otherwise.}\end{array}\right.$

Let ${\bf a}_{i.j}=(a_{1i},...,a_{ji})^{T}$ and ${\bf a}_{i(-j)}$
$=(a_{j+1,i},...,a_{pi})^{T}$. Note that ${\bf a}_{i(-0)}$ and ${\bf a}_{i.p}$
are equivalent to ${\bf a}_{i}$, the $i$th column of ${\bf A}$ and that ${\bf
a}_{i.1}=a_{1i}$. Then the Jacobian, $J({\bf A},{\bf B})$, of ${\bf A}$ with
respect to ${\bf B}$ is

$\bordermatrix{&{\bf b}_{1.1}^{T}&{\bf b}_{2.2}^{T}&{\bf
b}_{3.3}^{T}&\ldots&{\bf b}_{p.p}^{T}\cr{\bf a}_{1}&-{\bf a}_{1}{\bf
a}_{1.1}^{T}&-{\bf a}_{2}{\bf a}_{1.2}^{T}&-{\bf a}_{3}{\bf
a}_{1.3}^{T}&\ldots&-{\bf a}_{p}{\bf a}_{1.p}^{T}\cr{\bf a}_{2(-1)}&{\bf
0}&-{\bf a}_{2(-1)}{\bf a}_{2.2}^{T}&-{\bf a}_{3(-1)}{\bf
a}_{2.3}^{T}&\ldots&-{\bf a}_{p(-1)}{\bf a}_{2.p}^{T}\cr{\bf a}_{3(-2)}&{\bf
0}&{\bf 0}&-{\bf a}_{3(-2)}{\bf a}_{3.3}^{T}&\ldots&-{\bf a}_{p(-2)}{\bf
a}_{3.p}^{T}\cr\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\cr{\bf
a}_{p(-p+1)}&{\bf 0}&{\bf 0}&{\bf 0}&\ldots&-{\bf a}_{p(-p+1)}{\bf
a}_{p.p}^{T}}.$

Hence the Jacobian $J({\bf A},{\bf B})$ is a block upper triangular matrix.

Under the transformation ${\bf B}^{T}={\bf A}^{-1}$, the asymptotic
distribution of $\hat{{\bf c}}$ is not changed and in fact ${\bf V}_{11}({\bf
c},{\bf A})=$ ${\bf V}_{11}({\bf c},{\bf B})$.

Since ${\bf V}_{12}({\bf c},{\bf A})$ $={\bf V}_{12}({\bf c},{\bf B})J({\bf
A},{\bf B})^{T}$, we get the following lemma.

Lemma 4.2. The $(K-1)$ by $p(p+1)/2$ matrix ${\bf V}_{12}({\bf c},{\bf A})$ is
given by the following matrix multiplied by $-1/pr_{1}$:

$\bordermatrix{&{\bf a}_{1}^{T}&{\bf a}_{2(-1)}^{T}&\ldots&{\bf
a}_{p(-p+1)}^{T}\cr{\bf c}&{\bf c}{\bf a}_{1}^{T}&{\bf c}{\bf
a}_{2(-1)}^{T}&\ldots&{\bf c}{\bf a}_{p(-p+1)}^{T}}.$

Proof of Lemma 4.2. By Lemma 3.4, ${\bf V}_{12}({\bf c},{\bf A})$ has as its
$i$th block element for the covariance matrix of $\hat{{\bf c}}$ and
$\hat{{\bf a}}_{i(-i+1)}$

$-\frac{{\bf c}}{pr_{1}}\sum_{j=i}^{p}{\bf b}_{j.j}^{T}{\bf a}_{i.j}{\bf
a}_{j(-j+1)}^{T}~{}=~{}-\frac{1}{pr_{1}}{\bf c}{\bf a}_{i(-i+1)}^{T}$

because ${\bf B}^{T}{\bf A}={\bf I}_{p}$.

Let ${\bf A}_{-i}$ be the matrix formed by deleting the first $i$ columns and
the first $i$ rows of ${\bf A}$ simultaneously. Note that ${\bf A}_{-0}={\bf
A}$ and ${\bf A}_{-p+1}={\bf a}_{p(-p+1)}=a_{pp}$.

Lemma 4.3. The $p(p+1)/2$ by $p(p+1)/2$ matrix ${\bf V}_{22}({\bf c},{\bf A})$
is given by

$\bordermatrix{&{\bf a}_{1}^{T}&{\bf a}_{2(-1)}^{T}&\ldots&{\bf
a}_{p(-p+1)}^{T}\cr{\bf a}_{1}&{\bf A}{\bf A}^{T}+d{\bf a}_{1}{\bf
a}_{1}^{T}&e{\bf a}_{1}{\bf a}_{2(-1)}^{T}&\ldots&e{\bf a}_{1}{\bf
a}_{p(-p+1)}^{T}\cr{\bf a}_{2(-1)}&e{\bf a}_{2(-1)}{\bf a}_{1}^{T}&{\bf
A}_{-1}{\bf A}_{-1}^{T}+d{\bf a}_{2(-1)}{\bf a}_{2(-1)}^{T}&\ldots&e{\bf
a}_{2(-1)}{\bf a}_{p(-p+1)}^{T}\cr\vdots&\vdots&\vdots&\ddots&\vdots\cr{\bf
a}_{p(-p+1)}&e{\bf a}_{p(-p+1)}{\bf a}_{1}^{T}&e{\bf a}_{p(-p+1)}{\bf
a}_{2(-1)}^{T}&\ldots&(1+d){\bf a}_{p(-p+1)}{\bf a}_{p(-p+1)}^{T}},$

where $d$ and $e$ are defined as in Lemma 3.5.

Proof of Lemma 4.3. The $(i,j)$th block element of ${\bf V}_{22}({\bf c},{\bf
A})$ is the covariance matrix of $\hat{{\bf a}}_{i(-i+1)}$ and $\hat{{\bf
a}}_{j(-j+1)}$, and is given by

$\sum_{l=i}^{p}\sum_{m=j}^{p}{\bf a}_{l(-i+1)}{\bf a}_{i.l}^{T}{\bf
V}_{22}[l,m]{\bf a}_{j.m}{\bf a}_{m(-j+1)}^{T}~{}~{}~{}~{}(i\leq j),$

where ${\bf V}_{22}[l,m]$ denotes the $(l,m)$th block element of ${\bf
V}_{22}({\bf c},{\bf B})$ which is the covariance matrix of $\hat{{\bf
b}}_{l.l}$ and $\hat{{\bf b}}_{m.m}$. Let $\delta_{ij}$ be the Kronecker
delta. Since ${\bf a}_{i.l}^{T}{\bf b}_{l.l}~{}=$ $\delta_{il}$ and ${\bf
a}_{i.l}^{T}{\bf B}_{l}~{}=$ ${\bf 1}_{l.i}^{T}$, the $(i,j)$th block element
of ${\bf V}_{22}({\bf c},{\bf A})$ becomes

$e(1-\delta_{ij}){\bf a}_{i(-i+1)}{\bf a}_{j(-j+1)}^{T}+\delta_{ij}({\bf
A}_{-i+1}{\bf A}_{-i+1}^{T}+d{\bf a}_{i(-i+1)}{\bf a}_{i(-i+1)}^{T}),$

which completes the proof.

Lemmas 3.3, 4.2 and 4.3 yield the asymptotic distribution of $\hat{{\bf c}}$
and $\hat{{\bf A}}$ under the proportional covariance model.

Theorem 4.1. The asymptotic distribution of

$\sqrt{n_{+}}(\hat{{\bf c}}-{\bf c},\hat{{\bf A}}-{\bf A})$

is a multivariate normal with mean zero and covariance matrix ${\bf V}({\bf
c},{\bf A})$, as $n_{+}$ tends to infinity with each $r_{k}$ kept fixed.

Consider the single-population case, $r_{1}=1$. Since the coefficients $d$ and
$e$ become $-1/2$ and zero, respectively, Lemma 4.3 shows that the $\hat{{\bf
a}}_{i(-i+1)}$ $(1\leq i\leq p)$ are asymptotically independent. In this
single-population case, ${\bf S}_{1}=\hat{{\bf A}}\hat{{\bf A}}^{T}$ and
therefore Lemma 4.3 gives the asymptotic distribution of the Cholesky root of
the sample covariance matrix when the underlying distribution is a
multivariate one. The asymptotic distribution of the sample generalized
variance being equal to the product of the $\hat{a}_{ii}^{2}$ is easily
derived from Lemma 4.3 and it can be found also in Theorem 7.5.4 of Anderson
(1984).

##  5 Asymptotic distribution of the maximum likelihood estimators of the
$c_{k}$ and ${\bf\Sigma}_{1}$

Now, we are in a position to derive the asymptotic distribution of the maximum
likelihood estimators of the proportionality coefficients and the common
covariance matrix, and it will be done using the multivariate $\delta$-method.

Let $\sigma_{ij}$ be the $(i,j)$th element of ${\bf\Sigma}_{1}$. The common
covariance matrix ${\bf\Sigma}_{1}$ has $p(p+1)/2$ independent parameters. Let
${\bf I}_{p.i}$ be the $(p-i)$ by $p$ matrix in which the first $i$ columns
are zero vectors and the remaining $p-i$ columns form the unit matrix of order
$p-i$. When $i=0$, ${\bf I}_{p.0}$ becomes the unit matrix of order $p$ and is
denoted by ${\bf I}_{p}$. For the $i$th column $\mbox{\boldmath$\sigma$}_{i}$
of ${\bf\Sigma}_{1}$, $\mbox{\boldmath$\sigma$}_{i(-j)}$ is defined similarly
to ${\bf a}_{i(-j)}$.

Since $\sigma_{ij}$ $={\bf a}_{(i)}^{T}{\bf a}_{(j)}$, we get the Jacobian,
$J({\bf\Sigma}_{1},{\bf A})$, of ${\bf\Sigma}_{1}$ with respect to ${\bf A}$
as follows:

$\bordermatrix{&{\bf a}_{1}^{T}&{\bf a}_{2(-1)}^{T}&{\bf
a}_{3(-2)}^{T}&\ldots&{\bf a}_{p(-p+1)}^{T}\cr&a_{11}{\bf
I}_{p}&&&&\cr\mbox{\boldmath$\sigma$}_{1}&+&{\bf 0}&{\bf 0}&\ldots&{\bf
0}\cr&{\bf a}_{1}{\bf 1}_{p.1}^{T}&&&&\cr&a_{21}{\bf I}_{p.1}&a_{22}{\bf
I}_{p-1}&&&\cr\mbox{\boldmath$\sigma$}_{2(-1)}&+&+&{\bf 0}&\ldots&{\bf
0}\cr&{\bf a}_{1(-1)}{\bf 1}_{p.2}^{T}&{\bf a}_{2(-1)}{\bf
1}_{p-1.1}^{T}&&&\cr&a_{31}{\bf I}_{p.2}&a_{32}{\bf I}_{p-1.1}&a_{33}{\bf
I}_{p-2}&&\cr\mbox{\boldmath$\sigma$}_{3(-2)}&+&+&+&\ldots&{\bf 0}\cr&{\bf
a}_{1(-2)}{\bf 1}_{p.3}^{T}&{\bf a}_{2(-2)}{\bf 1}_{p-1.2}^{T}&{\bf
a}_{3(-2)}{\bf
1}_{p-2.1}^{T}&&\cr\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\cr&a_{p1}{\bf
I}_{p.p-1}&a_{p2}{\bf I}_{p-1.p-2}&a_{p3}{\bf
I}_{p-2.p-3}&&a_{pp}\cr\mbox{\boldmath$\sigma$}_{p(-p+1)}&+&+&+&\ldots&+\cr&{\bf
a}_{1(-p+1)}{\bf 1}_{p.p}^{T}&{\bf a}_{2(-p+1)}{\bf 1}_{p-1.p-1}^{T}&{\bf
a}_{3(-p+1)}{\bf 1}_{p-2.p-2}^{T}&&{\bf a}_{p(-p+1)}}.$

Hence the Jacobian $J({\bf\Sigma}_{1},{\bf A})$ is a block lower triangular
matrix.

Let ${\bf V}({\bf c},{\bf\Sigma}_{1})$ be the asymptotic covariance matrix for
the maximum likelihood estimators of the parameters ${\bf c}$ and
${\bf\Sigma}_{1}$, and let it be partitioned as

${\bf V}({\bf c},{\bf\Sigma}_{1})=\pmatrix{{\bf V}_{11}({\bf
c},{\bf\Sigma}_{1})&{\bf V}_{12}({\bf c},{\bf\Sigma}_{1})\cr{\bf V}_{21}({\bf
c},{\bf\Sigma}_{1})&{\bf V}_{22}({\bf c},{\bf\Sigma}_{1})\cr},$

such that ${\bf V}_{11}({\bf c},{\bf\Sigma}_{1})$ is the covariance matrix for
the maximum likelihood estimators of ${\bf c}$ and ${\bf V}_{22}({\bf
c},{\bf\Sigma}_{1})$ for the maximum likelihood estimator of the independent
parameters of ${\bf\Sigma}_{1}$.

The covariance matrix of $\hat{{\bf c}}$, ${\bf V}_{11}({\bf
c},{\bf\Sigma}_{1})$, is equivalent to ${\bf V}_{11}({\bf c},{\bf B})$ given
in Lemma 3.3.

We need some useful formulas to get ${\bf V}_{12}({\bf c},{\bf\Sigma}_{1})$
which are stated in the following lemma which is easily proved.

Lemma 5.1. For $j\leq i$, we have
(i) ${\bf I}_{p-j+1.i-j}{\bf a}_{j(-j+1)}={\bf a}_{j(-i+1)}$
(ii) ${\bf a}_{j(-j+1)}^{T}{\bf 1}_{p-j+1.i-j+1}=a_{ij}$
(iii) $\mbox{\boldmath$\sigma$}_{i(-i+1)}=\sum_{j=1}^{i}a_{ij}{\bf
a}_{j(-i+1)}.$

The covariance matrix ${\bf V}_{12}({\bf c},{\bf\Sigma}_{1})$ of $\hat{{\bf
c}}$ and $\hat{{\bf\Sigma}}_{1}$ is given by ${\bf V}_{12}({\bf c},{\bf
A})J({\bf\Sigma}_{1},{\bf A})^{T}$ as in the following lemma.

Lemma 5.2. The $(K-1)$ by $p(p+1)/2$ matrix ${\bf V}_{12}({\bf
c},{\bf\Sigma}_{1})$ is given by the following matrix multiplied by
$-2/pr_{1}$:

$\bordermatrix{&\mbox{\boldmath$\sigma$}_{1}^{T}&\mbox{\boldmath$\sigma$}_{2(-1)}^{T}&\ldots&\mbox{\boldmath$\sigma$}_{p(-p+1)}^{T}\cr{\bf
c}&{\bf c}\mbox{\boldmath$\sigma$}_{1}^{T}&{\bf
c}\mbox{\boldmath$\sigma$}_{2(-1)}^{T}&\ldots&{\bf
c}\mbox{\boldmath$\sigma$}_{p(-p+1)}^{T}}.$

Proof of Lemma 5.2. The covariance matrix of $\hat{{\bf c}}$ and
$\hat{\mbox{\boldmath$\sigma$}}_{i(-i+1)}$ which is the $i$th block element of
${\bf V}_{12}({\bf c},{\bf\Sigma}_{1})$ is given by

$-\frac{1}{pr_{1}}{\bf c}\sum_{j=1}^{i}{\bf a}_{j(-j+1)}^{T}(a_{ij}{\bf
I}_{p-j+1.i-j}^{T}+{\bf 1}_{p-j+1.i-j+1}{\bf a}_{j(-i+1)}^{T})$

which becomes by Lemma 5.1

$-\frac{2}{pr_{1}}{\bf c}\mbox{\boldmath$\sigma$}_{i(-i+1)}^{T}.$

Let
${\bf\Sigma}_{(i,j)}=(\mbox{\boldmath$\sigma$}_{j(-i+1)},...,\mbox{\boldmath$\sigma$}_{p(-i+1)})$.
Then ${\bf\Sigma}_{(i,j)}$ is a $(p-i+1)$ by $(p-j+1)$ matrix consisting of
the elements in the lower right corner of ${\bf\Sigma}_{1}$ starting at
$\sigma_{ij}$. The $(p-i+1)$ by $(p-j+1)$ matrix ${\bf A}_{(i,j)}$ is defined
similarly to ${\bf\Sigma}_{(i,j)}$.

Lemma 5.3. The $p(p+1)/2$ by $p(p+1)/2$ matrix ${\bf V}_{22}({\bf
c},{\bf\Sigma}_{1})$ has as its $(i,j)$ block element

$\sigma_{ij}{\bf\Sigma}_{(i,j)}+4e\mbox{\boldmath$\sigma$}_{i(-i+1)}\mbox{\boldmath$\sigma$}_{j(-j+1)}^{T}+\mbox{\boldmath$\sigma$}_{j(-i+1)}\mbox{\boldmath$\sigma$}_{i(-j+1)}^{T}~{}~{}~{}~{}(i\leq
j)$

which is the covariance matrix of $\hat{\mbox{\boldmath$\sigma$}}_{i(-i+1)}$
and $\hat{\mbox{\boldmath$\sigma$}}_{j(-j+1)}$, where $d$ and $e$ are defined
as in Lemma 3.5.

Since $3+4d=4e+1$, the covariance matrix of
$\hat{\mbox{\boldmath$\sigma$}}_{i(-i+1)}$ becomes

$\sigma_{ii}{\bf\Sigma}_{(i,i)}+(3+4d)\mbox{\boldmath$\sigma$}_{i(-i+1)}\mbox{\boldmath$\sigma$}_{i(-i+1)}^{T}.$

More formulas are needed to prove Lemma 5.3. Let ${\bf a}_{(i)(-j)}~{}=$
$(a_{i,j+1},$…$,a_{ip})^{T}.$ Then we get the following lemma together with
Lemma 5.1 useful for proving Lemma 5.3.

Lemma 5.4. For $i\leq j$, $1\leq m\leq i$ and $1\leq l\leq j$, we have

(i) ${\bf 1}_{p-m+1.i-m+1}^{T}{\bf A}_{-m+1}={\bf a}_{(i)(-m+1)}^{T}$
(ii) ${\bf I}_{p-m+1.i-m}{\bf A}_{-m+1}={\bf A}_{(i,m)}$
(iii) ${\bf a}_{(i)(-m+1)}^{T}{\bf A}_{(j,m)}^{T}=\sum_{u=m}^{i}a_{iu}{\bf
a}_{u(-j+1)}^{T}$
(iv) ${\bf a}_{(i)(-m+1)}^{T}{\bf a}_{(j)(-m+1)}=\sum_{u=m}^{i}a_{iu}a_{ju}$
(v)
$\mbox{\boldmath$\sigma$}_{i(-i+1)}\mbox{\boldmath$\sigma$}_{j(-j+1)}^{T}=\sum_{l=1}^{j}\sum_{m=1}^{i}a_{im}a_{jl}{\bf
a}_{m(-i+1)}{\bf a}_{l(-j+1)}^{T}$

(vi) $\mbox{\boldmath$\sigma$}_{j(-i+1)}\mbox{\boldmath$\sigma$}_{i(-j+1)}^{T}$ | $=\sum_{l=1}^{j}\sum_{m=1}^{i}a_{im}a_{jl}{\bf a}_{l(-i+1)}{\bf a}_{m(-j+1)}^{T}$
---|---
| $=\sum_{m=1}^{i}a_{im}a_{jm}{\bf a}_{m(-i+1)}{\bf a}_{m(-j+1)}^{T}$
| $+\sum_{m=1}^{i}\\{a_{jm}{\bf a}_{m(-i+1)}(\sum_{l=m+1}^{i}a_{il}{\bf
a}_{l(-j+1)}^{T})\\}$
| $+\sum_{m=1}^{i}\sum_{l=m+1}^{j}a_{im}a_{jl}{\bf a}_{l(-i+1)}{\bf
a}_{m(-j+1)}^{T}$

(vii) ${\bf\Sigma}_{(i,j)}=\sum_{u=1}^{m-1}{\bf a}_{u(-i+1)}{\bf
a}_{u(-j+1)}^{T}+{\bf A}_{(i,m)}{\bf A}_{(j,m)}^{T}$

(viii) $\sigma_{ij}{\bf\Sigma}_{(i,j)}$ | $=(\sum_{m=1}^{i}a_{im}a_{jm})(\sum_{l=1}^{p}{\bf a}_{l(-i+1)}{\bf a}_{l(-j+1)}^{T})$
---|---
| $=\sum_{m=1}^{i}a_{im}a_{jm}{\bf A}_{(i,m)}{\bf A}_{(j,m)}^{T}$
| $+\sum_{l=1}^{i-1}(\sum_{m=l+1}^{i}a_{im}a_{jm}){\bf a}_{l(-i+1)}{\bf
a}_{l(-j+1)}^{T}.$

Proof of Lemma 5.3. For $i\leq j$, the $(i,j)$th block element of ${\bf
V}_{22}({\bf c},{\bf\Sigma}_{1})$ which is the covariance matrix of
$\hat{\mbox{\boldmath$\sigma$}}_{i(-i+1)}$ and
$\hat{\mbox{\boldmath$\sigma$}}_{j(-j+1)}$ is given by

$\sum_{l=1}^{j}[\sum_{m=1}^{i}\\{(a_{im}{\bf I}_{p-m+1.i-m}+{\bf
a}_{m(-i+1)}{\bf 1}_{p-m+1.i-m+1}^{T}){\bf V}_{22}[m,l]\\}$ $(a_{jl}{\bf
I}_{p-l+1.j-l}^{T}+{\bf 1}_{p-l+1.j-l+1}{\bf a}_{l(-j+1)}^{T})],$ (4)

where ${\bf V}_{22}[m,l]$ denotes the $(m,l)$th block element of ${\bf
V}_{22}({\bf c},{\bf A})$. By Lemmas 5.1 and 5.4, summing the summands in (4)
over $1\leq l=m\leq i$ gives

$\displaystyle\sum_{m=1}^{i}\\{a_{im}a_{jm}{\bf A}_{(i,m)}{\bf
A}_{(j,m)}^{T}+4da_{im}a_{jm}{\bf a}_{m(-i+1)}{\bf a}_{m(-j+1)}^{T}$
$\displaystyle+$ $\displaystyle a_{jm}{\bf a}_{m(-i+1)}{\bf
a}_{(i)(-m+1)}^{T}{\bf A}_{(j,m)}^{T}+a_{im}{\bf A}_{(i,m)}{\bf
a}_{(j)(-m+1)}{\bf a}_{m(-j+1)}^{T}$ $\displaystyle+$
$\displaystyle(\sum_{u=m}^{i}a_{iu}a_{ju}){\bf a}_{m(-i+1)}{\bf
a}_{m(-j+1)}^{T}\\}$ (5)

and taking a summation of the summands in (4) over $l\neq m$ yields

$4e\mbox{\boldmath$\sigma$}_{i(-i+1)}\mbox{\boldmath$\sigma$}_{j(-j+1)}^{T}-4e\sum_{m=1}^{i}a_{im}a_{jm}{\bf
a}_{m(-i+1)}{\bf a}_{m(-j+1)}^{T}.$ (6)

Applying Lemma 5.4 again to (5) and (6) completes the proof since $4d+3=1+4e$.

From Lemmas 3.3, 5.2 and 5.3, we get the asymptotic distribution of $\hat{{\bf
c}}$ and $\hat{{\bf\Sigma}}_{1}$ under the proportional covariance model (1).

Theorem 5.1. The statistic

$\sqrt{n_{+}}(\hat{{\bf c}}-{\bf c},\hat{{\bf\Sigma}}_{1}-{\bf\Sigma}_{1})$

is asymptotically distributed according to a multivariate normal distribution
with mean zero and covariance matrix ${\bf V}({\bf c},{\bf\Sigma}_{1})$ as
$n_{+}$ tends to infinity with each $r_{k}$ held fixed.

The asymptotic distribution of $\hat{{\bf c}}$ is free of ${\bf\Sigma}_{1}$ as
can be seen in Lemma 3.3, and that of $\hat{{\bf\Sigma}}_{1}$ is independent
of ${\bf c}$. However, $\hat{{\bf c}}$ and $\hat{{\bf\Sigma}}_{1}$ are not
asymptotically independent.

When $r_{1}=1$, the $(i,j)$th block element of ${\bf V}_{22}({\bf
c},{\bf\Sigma}_{1})$ becomes

$\sigma_{ij}{\bf\Sigma}_{(i,j)}+\mbox{\boldmath$\sigma$}_{j(-i+1)}\mbox{\boldmath$\sigma$}_{i(-j+1)}^{T}$
(7)

because $e$ is zero. In this case, the maximum likelihood estimator of
${\bf\Sigma}_{1}$ becomes ${\bf S}_{1}$, derived from the Wishart
distribution, which is an unbiased estimator of ${\bf\Sigma}_{1}$. However,
the maximum likelihood estimator of ${\bf\Sigma}_{1}$ based on the
multivariate normal distribution is in fact $(n_{1}/N_{1}){\bf S}_{1}$. Since
${\bf S}_{1}$ and $(n_{1}/N_{1}){\bf S}_{1}$ have the same limiting
distribution, this case $(r_{1}=1)$ yields the asymptotic distribution of the
maximum likelihood estimator of ${\bf\Sigma}_{1}$. We can easily check if (7)
is identical with (15) in Anderson (1984, p.82).

##  6 A hypothesis for homogeneity of covariance matrices

As an application of Theorem 5.1, consider a hypothesis for homogeneity of
covariance matrices

$H_{0}:c_{2}=...=c_{K}=1.$ (8)

The inverse ${\bf U}_{11}({\bf c},{\bf B})$ of ${\bf V}_{11}({\bf
c},{\bf\Sigma}_{1})$ is needed to test the hypothesis (8), which is given in
Lemma 3.2.

The asymptotic distribution of

$n_{+}(\hat{c}_{2}-1,...,\hat{c}_{K}-1){\bf U}_{11}({\bf c},{\bf
B})(\hat{c}_{2}-1,...,\hat{c}_{K}-1)^{T}$

under (8) is a chi-squared distribution with degrees of freedom $K-1$.
However, ${\bf U}_{11}({\bf c},{\bf B})$ is still unknown and can be estimated
by replacing ${\bf c}$ with the corresponding maximum likelihood estimator. We
denote the estimated ${\bf U}_{11}({\bf c},{\bf B})$ by $\hat{{\bf
U}}_{11}({\bf c},{\bf B})$. Since $\hat{{\bf c}}$ is a consistent estimator of
${\bf c}$, $\hat{{\bf U}}_{11}({\bf c},{\bf B})$ converges in probability to
${\bf U}_{11}({\bf c},{\bf B})$ and therefore a test statistic for testing the
hypothesis (8) is given by

$n_{+}(\hat{c}_{2}-1,...,\hat{c}_{K}-1)\hat{{\bf U}}_{11}({\bf c},{\bf
B})(\hat{c}_{2}-1,...,\hat{c}_{K}-1)^{T}$ (9)

whose asymptotic distribution under (8) is a chi-squared distribution with
$K-1$ degrees of freedom (Rao, 1973, Chapter 6). The test statistic (9) can be
expressed in a simple form

$n_{+}\\{\frac{pr_{1}(1-r_{1})}{2}-pr_{1}\sum_{k=2}^{K}\frac{r_{k}}{\hat{c}_{k}}+\frac{p}{2}[\sum_{k=2}^{K}\frac{r_{k}}{\hat{c}_{k}^{2}}-(\sum_{k=2}^{K}\frac{r_{k}}{\hat{c}_{k}})^{2}]\\}$
$=\frac{n_{+}p}{2}[\sum_{k=1}^{K}\frac{r_{k}}{\hat{c}_{k}^{2}}-(\sum_{k=1}^{K}\frac{r_{k}}{\hat{c}_{k}})^{2}],$

where $\hat{c}_{1}\equiv 1$.

References

Anderson, T.W. (1984). An Introduction to Multivariate Statistical Analysis.
2nd ed. Wiley, New York.

Bartlett, M.S. (1951). An inverse matrix adjustment arising in discriminant
analysis. Ann. Math. Statist. 22, 107–111.

Boente, G., Critchley, F. and Orellana, L. (2007). Influence functions of two
families of robust estimators under proportional scatter matrices, Statistical
Methods and Applications, 15, 295–327

Boente, G., Critchley, F. and Orellana, L. (2009). Robust tests for the common
principal components model, Journal of Statistical Planning and Inference,
139, 1332–1347

Eriksen, P.S. (1987). Proportionality of covariance matrices. Ann. Statist.
15, 732–748.

Flury, B.N. (1986). Proportionality of k covariance matrices. Statist. Prob.
Letters 4, 29–33.

Guttman, I., D.Y. Kim and I. Olkin (1985). Statistical inference for constants
of proportionality. In Multivariate Analysis, Vol VI (P.R. Krishnaiah, ed.),
Elsevier Science Publishers B.V.

Jöreskog, K.G. (1971). Simultaneous factor analysis in several populations.
Psychometrika 36, 409–426.

Jöreskog, K.G. and D. Sörbom (1996). LISREL 8: User’s reference guide, Hove
and London: Scientific Software International.

Khatri, C. G. (1967). Some distribution problems connected with the
characteristic roots of $S_{1}S_{2}^{-1}$. Annals of Mathematical Statistics
38, 944–948.

Manly, B. F. and Rayner, J. C. W. (1987). The Comparison of Sample Covariance
Matrices Using Likelihood Ratio Tests. Biometrika, 74, 841–847.

Owen, A. (1984). A neighborhood.based classifier for landsat data. Canadian
Journal of Statistics, 12, 191–200.

Rao, C.R. (1973). Linear Statistical Inference and its Applications. 2nd ed.
Wiley, New York.

Rao, C. R. (1983). Likelihood tests for relationships between covariance
matrices. In Studies in Econometrics, Time Series and Multivariate Statistics.
S. Karlim, T. Ameniya, L. A. Goodman, editors. pp. 529–543. Academic Press,
New York.