text:
# Infinite Products of Large Random Matrices and Matrix-valued Diffusion

## 1 Introduction

Applications of random matrix theory (RMT) in physics range from interpretation of complex spectra of energy levels in atomic and nuclear physics , studies of disordered systems , chaotic behavior  to quantum gravity . Use of progressively developed methods of RMT has also proved to be of particular interest in other sciences such as meteorology , image processing , population ecology  or economy .

In this work, we will be concerned with infinite products of random matrices (PRM) sampled from general Gaussian ensembles. Products of that type appear in various branches of physics and interdisciplinary applications  where system dynamics is described in terms of evolution operators with random coefficients. Perhaps the best known examples are those related to thermal properties of disordered magnetic systems, in particular, to the localization of electronic wave functions in random potentials . Attention has also been focused on applications of PRM to the analysis of chaotic dynamical systems , where stability of trajectories and their sensitive dependence on initial conditions are measured in terms of characteristic Lyapunov exponents. Moreover, there are several results exploring the use of PRM formalism in applied physics and interdisciplinary research, ranging from the studies of stability of large eco- and social- systems , adaptive algorithms and the analysis of system performance under the influence of external noises  to image compression  and communication via antenna arrays . However, despite the richness of the potential applications, the number of analytical results available in the PRM theory is still rather limited.

In this paper we develop effective calculational techniques for studying products of random matrices in the large $`N`$ limit, and use these tools to derive the properties of the natural matrix valued generalization of the geometric (multiplicative) diffusion type process. The scalar versions of these processes, leading to log-normal distributions are ubiquitous in various fields  and we expect the matrix valued generalizations to find numerous applications.

Let us briefly review the conventional (scalar) multiplicative (geometric) random walk in one dimension in the presence of some external drift force. The process belongs to the class of (Markovian) Ito diffusion processes, whose stochastic variable $`s`$ undergoes an evolution

$`{\displaystyle \frac{ds}{s}}=\mu (s,t)dt+\sigma (s,t)dW_t`$ (1)
In the above stochastic differential equation (SDE), part of the evolution is driven by a variable $`dW_t`$ that represents the Wiener process (integral of the Gaussian white-noise), respecting

$`dW_t=0dW_t^2=dt`$ (2)
For simplicity we will limit ourselves to the constant drift $`\mu `$ and constant variance $`\sigma `$. Finite time evolution of the system could be viewed as a string of increments

$`s(T)=\underset{M\mathrm{}}{lim}[(1+\mu {\displaystyle \frac{T}{M}}+\sigma \sqrt{{\displaystyle \frac{T}{M}}}x_1)(1+\mu {\displaystyle \frac{T}{M}}+\sigma \sqrt{{\displaystyle \frac{T}{M}}}x_2)\mathrm{}`$
$`\mathrm{}(1+\mu {\displaystyle \frac{T}{M}}+\sigma \sqrt{{\displaystyle \frac{T}{M}}}x_M)]s(0)`$ (3)
where $`M`$ stands for the number of steps in the discretization of the time interval. After averaging $`s(T)`$ over the independent identical distributions (iid) of the Gaussian variables $`x_i`$ and taking the limit $`M\mathrm{}`$, we recover555The log-normal distribution $`s(T)`$ is easy to infer looking at the form of the product (3). Taking the logarithm, expanding and using the central limit theorem we immediately see, that the r.h.s. tends to the Gaussian distribution. the well known solution for the probability density of $`s(T)`$:

$`p(s,T|s_0,0)={\displaystyle \frac{1}{s\sqrt{2\pi \sigma ^2T}}}\mathrm{exp}\left[{\displaystyle \frac{(\mathrm{log}(s/s_0)\mu T+\frac{1}{2}\sigma ^2T)^2}{2\sigma ^2T}}\right]`$ (4)
with $`s`$ restricted to positive values.

The aim of this paper is to study a similar construction in the space of matrices taking the defining equation (3) as a guideline for the generalization. We are therefore interested in properties of the matrix-valued evolution operator defined as

$`Y(T)=\underset{M\mathrm{}}{lim}[(1+\mu {\displaystyle \frac{T}{M}}+\sigma \sqrt{{\displaystyle \frac{T}{M}}}X_1)(1+\mu {\displaystyle \frac{T}{M}}+\sigma \sqrt{{\displaystyle \frac{T}{M}}}X_2)\mathrm{}`$
$`\mathrm{}(1+\mu {\displaystyle \frac{T}{M}}+\sigma \sqrt{{\displaystyle \frac{T}{M}}}X_M)]`$ (5)
where $`\mu `$ is some deterministic “drift” matrix and the stochastic matrices $`X_i`$ belong to identical independent random matrix ensembles. In particular we will be studying the eigenvalue distribution of $`Y(T)`$:

$$\rho _T(z)=\frac{1}{N}\text{tr}\delta ^{(2)}\left(zY(T)\right)$$
(6)
where the average is taken over the stochastic $`N`$ by $`N`$ matrices $`X_i`$ appearing in the definition of $`Y(T)`$.

Contrary to several standard multidimensional extensions of Brownian walks, we concentrate here on studying how the full spectrum of operator $`Y`$ evolves with time $`T`$.

Let us note some key features of the operator $`Y`$. First, since the matrices $`X_i`$ in general do not commute, we are dealing with a ‘path ordered product’. Second, even if the matrices $`X_i`$ are hermitian, their product is not, i.e. the spectrum in general disperses into the complex plane, showing - as pointed out in this paper - some rather unusual feature described as a ”topological phase transition”. Namely the support of the eigenvalue distribution changes from a simply-connected two-dimensional island to a two-dimensional island with a hole.

We note that the time evolution of some initial vector $`|0>`$ under $`Y`$ forms a very general setup for several multivariate stochastic evolutions of the type

$`|\tau =Y(\tau )|0,`$ (7)
Analytical results for such processes are however scarce. In a recent study, Jackson et al.  solved the matrix-valued evolution driven by Gaussian noise of the type considered above in the case of 2 by 2 hermitian matrices. There are also some known results for the diffusion of the norm of vector $`|\tau `$, since this problem can be linked to the Lyapunov exponents of the hermitian matrix $`Y^{}Y`$. We briefly comment on this issue in section 9.

The main result of this work is the derivation of exact (finite $`M`$) and asymptotic ($`M\mathrm{}`$) spectral formulas for the evolution operators $`Y`$, in the limit where the dimension $`N`$ of $`Y`$ tends to infinity. For simplicity, we set to zero the deterministic entries $`\mu `$ – although one can easily include them using our formalism. In this paper we also limit ourselves to ensembles of the Gaussian type. Generalizations to non-Gaussian measures (involving finite and infinite spectral supports) will be investigated elsewhere. We discuss here two cases depending on whether the $`X_i`$’s are (i) hermitian Gaussian matrices (GUE) or (ii) complex (arbitrary) Gaussian matrices (GCE). Similar constructions could be used for real and quaternionic matrices.

The main motivation for this work is to formulate a new and relatively simple theoretical framework for studying spectra of the diffusive-like evolution operators with their further applications to various branches of the complex systems analysis in physics, mathematics and in interdisciplinary science.

In particular, the formalism outlined here could be used as a starting point for matrix versions of super- or sub-diffusions, generated by e.g. matrix Lévy analogs of the stable power-tailed distributions. It will be used also for the study of spectral statistics of infinite products of unitary operators (e.g. Polyakov lines) .

Most of the derivations presented in this paper will be based on diagrammatic methods. In order to make the paper self-contained, in the first two sections we briefly present the formalism, first for hermitian RMM (section 2) and then in section 3, we recall the extension of diagrammatic methods to the non-hermitian case , by using complex Gaussian non-hermitian random matrix model as an example. In section 4 we move to the main part of our paper and introduce our method for the treatment of products of random matrices by performing calculations for a product of two matrices. Then (section 5) we move toward the products of $`M`$ complex matrices, and we perform the limit $`M\mathrm{}`$. We present the final formulas for the spectral density and the equations defining the boundary of the two-dimensional support of the eigenvalues. In sections 6 and 7 we consider the case of hermitian matrices which, surprisingly, is more difficult than the previous one. We analyze first (section 6) the product of two matrices. This is interesting on its own, since it exhibits quite singular behavior – for evolution times smaller than a certain critical value, the a priori complex eigenvalues are ”frozen” to the real axis. In section 7 we move toward the general case. For products of matrices involving more than two hermitian matrices, the eigenvalues always form two-dimensional islands. We solve the spectral problem in the large $`M`$ limit. In section 8 we show the numerical analysis, confirming our analytical predictions and the existence of “topological phase transitions”, corresponding to structural changes of the spectrum.

Finally, in section 9 we briefly discuss the problems of stability of such products and the relation of our results to the standard Lyapunov exponents. Section 10 concludes the paper.

## 2 Hermitian Random Ensembles

A key problem in random matrix theories is to find the distribution of eigenvalues $`\lambda _i`$, in the large $`N`$ (size of the matrix $`H`$) limit, i.e

$`\rho (\lambda )={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\delta (\lambda \lambda _i)`$ (8)
where the averaging $`\mathrm{}`$ is done over the ensemble of $`N\times N`$ random hermitian matrices generated with probability

$`P(H)e^{N\mathrm{Tr}V(H)}.`$ (9)
The eigenvalues of course lie on the real axis. By introducing the resolvent (Green’s function)

$`G(z)={\displaystyle \frac{1}{N}}\mathrm{Tr}{\displaystyle \frac{1}{𝐳H}}.`$ (10)
with $`𝐳=z\mathrm{𝟏}_N`$ and by using the standard relation

$`{\displaystyle \frac{1}{\lambda \pm iϵ}}=\mathrm{P}.\mathrm{V}.{\displaystyle \frac{1}{\lambda }}i\pi \delta (\lambda ).`$ (11)
the spectral function (8) can be derived from the discontinuities of the Green’s function (10)

$`{\displaystyle \frac{1}{2\pi i}}\underset{ϵ0}{lim}\left(G(\lambda iϵ)G(\lambda +iϵ)\right)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}\delta (\lambda H)`$ (12)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\delta (\lambda \lambda _i)=\rho (\lambda ).`$
There are several ways of calculating Green’s functions for HRMM . We will follow the diagrammatic approach introduced by . A starting point of the approach is the expression allowing for the reconstruction of the Green’s function from all the moments $`\mathrm{Tr}H^n`$,

$`G(z)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}{\displaystyle \frac{1}{𝐳H}}={\displaystyle \frac{1}{N}}\mathrm{Tr}\left[{\displaystyle \frac{1}{𝐳}}+{\displaystyle \frac{1}{𝐳}}H{\displaystyle \frac{1}{𝐳}}+{\displaystyle \frac{1}{𝐳}}H{\displaystyle \frac{1}{𝐳}}H{\displaystyle \frac{1}{𝐳}}+\mathrm{}\right]`$ (13)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{z^{n+1}}}\mathrm{Tr}H^n`$
The reason why the above procedure works correctly for hermitian matrix models is the fact that the Green’s function is guaranteed to be holomorphic in the whole complex plane except at most on one or more 1-dimensional intervals. We will use the diagrammatic method to evaluate efficiently the sum of the moments on the right hand side. We will now restrict ourselves to the well known case of a random hermitian ensemble with Gaussian distribution.

The first step is to introduce a generating function with a matrix-valued source $`J`$:

$`Z(J)={\displaystyle 𝑑He^{\frac{N}{2}(\mathrm{Tr}H^2)+\mathrm{Tr}HJ}}.`$ (14)
where we integrate over all $`N^2`$ elements of the matrix $`H`$. All the moments follow directly from $`Z(J)`$ through the relation

$`\mathrm{Tr}H^n={\displaystyle \frac{1}{Z(0)}}\mathrm{Tr}\left({\displaystyle \frac{}{J}}\right)^nZ(J)|_{J=0}`$ (15)
and are straightforward to calculate, since in the Gaussian case the partition function (14) reads $`Z(J)=\mathrm{exp}\frac{1}{2N}\mathrm{Tr}J^2.`$ Accordingly, the lowest nonzero expectation value is

$`H_b^aH_d^c={\displaystyle \frac{^2Z(J)}{J_a^bJ_c^d}}|_{J=0}={\displaystyle \frac{1}{N}}{\displaystyle \frac{J_d^cZ(J)}{J_a^b}}|_{J=0}={\displaystyle \frac{1}{N}}\delta _b^c\delta _d^a.`$ (16)
and the next non-vanishing expectation value reads

$`H_b^aH_d^cH_f^eH_h^g={\displaystyle \frac{1}{N^2}}\left(\delta _b^c\delta _d^a\delta _f^g\delta _h^e+\delta _h^a\delta _g^b\delta _f^c\delta _e^d+\delta _f^a\delta _e^b\delta _h^c\delta _g^d\right)`$ (17)
The key idea in the diagrammatic approach is to associate to the expressions for the moments, like the above, a graphical representation following from a simple set of rules. The power of the approach is that it enables to perform a resummation of the whole power series (13) through the identification of the structure of the relevant graphs.

We depict the “Feynman” rules in Fig. 1, similar to the standard large $`N`$ diagrammatics for QCD . The $`1/𝐳=1/z\delta _a^b`$ in (13) is represented by a horizontal straight line. The propagator (16) is depicted as a double line.

The diagrammatic expansion of the Green’s function is visualized in Fig. 2, where one connects the vertices with the double line propagators in all possible ways. Each “propagator” brings a factor of $`1/N`$, and each loop a factor of $`\delta _a^a=N`$. From the three terms, corresponding to (17) contributing to $`\text{tr}H^4`$ only the first two are presented in Fig. 2 (the third and the fourth diagram). The diagram corresponding to the third term in (17) represents a non-planar contribution which is suppressed as $`1/N^2`$ and hence vanishes when $`N\mathrm{}`$. In general, only planar graphs survive the large $`N`$ limit.

The resummation of (13) is done by introducing the self-energy $`\mathrm{\Sigma }`$ comprising the sum of all one-particle irreducible graphs (rainbow-like). Then the Green’s function reads

$$G(z)=\frac{1}{z\mathrm{\Sigma }(z)}.$$
(18)
In the large $`N`$ limit the equation for the self energy $`\mathrm{\Sigma }`$, follows from resumming the rainbow-like diagrams of Fig. 2. The resulting equation (“Schwinger-Dyson” equation of Fig. 3) encodes pictorially the structure of these graphs and reads

$$\mathrm{\Sigma }=\frac{1}{N}\mathrm{Tr}G\mathrm{𝟏}=\frac{N}{N}G=G.$$
(19)
Equations (18) and (19) give immediately $`G(zG)=1`$ which can be solved to yield

$`G_{}(z)={\displaystyle \frac{1}{2}}(z\sqrt{z^24})`$ (20)
Only the $`G_{}`$ is a normalizable solution, with the proper asymptotic behavior $`G(z)1/z`$ in the $`z\mathrm{}`$ limit. From the discontinuity (cut), using (12), we recover Wigner’s semicircle  for the distribution of the eigenvalues for hermitian random matrices

$`\rho (\lambda )={\displaystyle \frac{1}{2\pi }}\sqrt{4\lambda ^2}.`$ (21)
## 3 Non-hermitian Random Ensembles

The main difficulty in the treatment of non-hermitian random matrix models is the fact that now the eigenvalues accumulate in two-dimensional domains in the complex plane and the Green’s function is no longer holomorphic. Therefore the power series expansion (13) no longer captures the full information about the Green’s function. In particular the eigenvalue distribution is related to the non-analytic (non-holomorphic behaviour of the Green’s function:

$$\frac{1}{\pi }_{\overline{z}}G(z)=\rho (z).$$
(22)
This phenomenon can be easily seen even in the the simplest non-hermitian ensemble — the Ginibre-Girko one , with non-hermitian matrices $`X`$, and measure

$`P(X)=e^{N\mathrm{Tr}XX^{}}.`$ (23)
It is easy to verify that all moments vanish $`\text{tr}X^n=0`$, for $`n>0`$ so the expansion (13) gives the Green’s function to be $`G(z)=1/z`$ (diagrammatically this follows from the fact that the propagator $`<X_b^aX_d^c>`$ vanishes and hence the self-energy $`\mathrm{\Sigma }=0`$). The true answer is, however, different. Only for $`|z|>1`$ one has indeed $`G(z)=1/z`$. For $`|z|<1`$ the Green’s function is nonholomorphic and equals $`G(z)=\overline{z}`$.

In spite of the above difficulties a diagrammatic approach can be used for recovering the full information including the eigenvalue density $`\rho (z)`$. The first step is to define a regularized Green’s function from which one can extract $`\rho (z)`$. This has been done by exploiting the analogy to two-dimensional electrostatics . The resulting regularized Green’s function is however difficult to calculate. The second step is to give an equivalent linearized form which can form a basis for a diagrammatic calculation . We will now briefly review the above developments.

Let us define the “electrostatic potential”

$`F={\displaystyle \frac{1}{N}}\mathrm{Tr}\mathrm{ln}[(𝐳)(\overline{𝐳}^{})+ϵ^2].`$ (24)
Then

$`\underset{ϵ0}{lim}{\displaystyle \frac{^2F(z,\overline{z})}{z\overline{z}}}`$ $`=`$ $`\underset{ϵ0}{lim}{\displaystyle \frac{1}{N}}\mathrm{Tr}{\displaystyle \frac{ϵ^2}{(|𝐳|^2+ϵ^2)^2}}`$ (25)
$`=`$ $`{\displaystyle \frac{\pi }{N}}{\displaystyle \underset{i}{}}\delta ^{(2)}(z\lambda _i)\pi \rho (x,y)`$
represents Gauss law, where $`z=x+iy`$. The last equality follows from the representation of the complex Dirac delta

$`\pi \delta ^{(2)}(z\lambda _i)=\underset{ϵ0}{lim}{\displaystyle \frac{ϵ^2}{(ϵ^2+|z\lambda _i|^2)^2}}`$ (26)
In the spirit of the electrostatic analogy we can define the Green’s function $`G(z,\overline{z})`$, as an “electric field”

$`G{\displaystyle \frac{F}{z}}={\displaystyle \frac{1}{N}}\underset{ϵ0}{lim}\mathrm{Tr}{\displaystyle \frac{\overline{𝐳}X^{}}{(\overline{𝐳}X^{})(𝐳X)+ϵ^2)}}.`$ (27)
Then Gauss law reads

$`_{\overline{z}}G=\pi \rho (x,y)`$ (28)
in agreement with (22).

As mentioned above, instead of working ab initio with the object (27), and in view of applying diagrammatic methods it is much more convenient to proceed differently – as we will now discuss.

Following  we define the matrix-valued resolvent through

$`\widehat{𝒢}`$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}_{\mathrm{B2}}\left(\begin{array}{cc}𝐳\mathrm{X}& \mathrm{i}ϵ\\ \mathrm{i}ϵ& \overline{𝐳}\mathrm{X}^{}\end{array}\right)^1=`$ (31)
$`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{Tr}_{\mathrm{B2}}\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{cc}𝒢_{11}& 𝒢_{1\overline{1}}\\ 𝒢_{\overline{1}1}& 𝒢_{\overline{1}\overline{1}}\end{array}\right)`$ (36)
with

$`A`$ $`=`$ $`{\displaystyle \frac{\overline{𝐳}X^{}}{(\overline{𝐳}X^{})(𝐳X)+ϵ^2}}`$
$`B`$ $`=`$ $`{\displaystyle \frac{iϵ}{(𝐳X)(\overline{𝐳}X^{})+ϵ^2}}`$
$`C`$ $`=`$ $`{\displaystyle \frac{iϵ}{(\overline{𝐳}X^{})(𝐳X)+ϵ^2}}`$
$`D`$ $`=`$ $`{\displaystyle \frac{𝐳X}{(𝐳X)(\overline{𝐳}X^{})+ϵ^2}}`$ (37)
and where we introduced the ‘block trace’ defined as

$`\mathrm{Tr}_{\mathrm{B2}}\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)_{2N\times 2N}\left(\begin{array}{cc}\mathrm{Tr}A& \mathrm{Tr}B\\ \mathrm{Tr}C& \mathrm{Tr}D\end{array}\right)_{2\times 2}.`$ (42)
Then, by definition, the upper-right component $`𝒢_{11}`$, is equal to the Green’s function (27).

The block approach has several advantages. First of all it is linear in the random matrices $`X`$ allowing for a simple diagrammatic calculational procedure. Let us define $`2N`$ by $`2N`$ matrices

$`𝓩=\left(\begin{array}{cc}𝐳& iϵ\mathrm{𝟏}\\ iϵ\mathrm{𝟏}& \overline{𝐳}\end{array}\right),=\left(\begin{array}{cc}X& 0\\ 0& X^{}\end{array}\right).`$ (47)
Then the generalized Green’s function is given formally by the same definition as the usual Green’s function $`G`$,

$`𝒢={\displaystyle \frac{1}{N}}\mathrm{Tr}_{\mathrm{B2}}{\displaystyle \frac{1}{𝓩\text{ }}}.`$ (48)
What is more important, also in this case the Green’s function is completely determined by the knowledge of all matrix-valued moments

$`\mathrm{Tr}_{\mathrm{B2}}𝓩^1𝓩^1\mathrm{}𝓩^1.`$ (49)
This last observation allows for a diagrammatic interpretation. The Feynman rules are analogous to the hermitian ones, only now one has to keep track of the block structure of the matrices, e.g. single straight lines will now be associated with a matrix factor $`𝒵^1`$. We will now demonstrate the above procedure by solving diagrammatically the complex Gaussian Random Matrix Model.

### Example: the Girko-Ginibre ensemble

The ensemble is defined by the measure

$`P(X)e^{N\mathrm{Tr}XX^{}}.`$ (50)
In this case the double line propagators are

$`X_b^aX_d^c`$ $`=`$ $`X_{}^{}{}_{b}{}^{a}X_{}^{}{}_{d}{}^{c}=0`$
$`X_b^aX_{}^{}{}_{d}{}^{c}`$ $`=`$ $`X_{}^{}{}_{b}{}^{a}X_d^c={\displaystyle \frac{1}{N}}\delta _d^a\delta _c^b.`$ (51)
As previously, we introduce the self-energy $`\stackrel{~}{\mathrm{\Sigma }}`$, (but which is now matrix-valued), in terms of which we get a two by two matrix expression

$`𝒢=(𝒵\stackrel{~}{\mathrm{\Sigma }})^1.`$ (52)
where the 2 by 2 matrix $`𝒵`$ reads

$`𝒵=\left(\begin{array}{cc}z& iϵ\\ iϵ& \overline{z}\end{array}\right)`$ (55)
The resummation of the rainbow diagrams for $`\stackrel{~}{\mathrm{\Sigma }}`$ is more subtle, but follows easily from the structure of the propagators (51). The analogue of (19) is now:

$`\stackrel{~}{\mathrm{\Sigma }}`$ $``$ $`\left(\begin{array}{cc}\mathrm{\Sigma }_{11}& \mathrm{\Sigma }_{1\overline{1}}\\ \mathrm{\Sigma }_{\overline{1}1}& \mathrm{\Sigma }_{\overline{1}\overline{1}}\end{array}\right)=\left(\begin{array}{cc}0& 𝒢_{1\overline{1}}\\ 𝒢_{\overline{1}1}& 0\end{array}\right).`$ (60)
The two by two matrix equations (52-60) completely determine the problem of finding the eigenvalue distribution for the Girko-Ginibre ensemble. Inserting (60) into (52) we get:

$$\left(\begin{array}{cc}𝒢_{11}& 𝒢_{1\overline{1}}\\ 𝒢_{\overline{1}1}& 𝒢_{\overline{1}\overline{1}}\end{array}\right)=\frac{1}{|z|^2𝒢_{1\overline{1}}𝒢_{\overline{1}1}}\left(\begin{array}{cc}\overline{z}& 𝒢_{1\overline{1}}\\ 𝒢_{\overline{1}1}& z\end{array}\right).$$
(61)
Note that at this moment we can safely put to zero the regulators $`ϵ`$. Looking at the off-diagonal equation

$`𝒢_{1\overline{1}}={\displaystyle \frac{𝒢_{1\overline{1}}}{|z|^2𝒢_{1\overline{1}}𝒢_{\overline{1}1}}}`$ (62)
we see that there are two solutions: one with $`𝒢_{1\overline{1}}=0`$, and another with $`𝒢_{1\overline{1}}0`$. The first one leads to a holomorphic Green’s function, and a straightforward calculation gives

$`G(z)={\displaystyle \frac{1}{z}}.`$ (63)
The second one is nonholomorhic, imposing the condition

$`|z|^2b^2=1`$ (64)
where we denoted $`𝒢_{1\overline{1}}𝒢_{\overline{1}1}b^2`$, hence

$`G(z,\overline{z})=\overline{z}`$ (65)
which leads, via the Gauss law, to

$`\rho (x,y)={\displaystyle \frac{1}{\pi }}{\displaystyle \frac{}{\overline{z}}}G(z,\overline{z})={\displaystyle \frac{1}{\pi }}.`$ (66)
Both solutions match at the boundary $`b^2=0`$ , which in this case reads $`z\overline{z}=1`$. In such a simple way we recovered the results of Ginibre and Girko for the complex non-hermitian ensemble. The eigenvalues are uniformly distributed on the unit disk $`|z|^2<1`$.

This example illustrates more general properties of the matrix valued generalized Green’s function. Each component of the matrix carries important information about the stochastic properties of the system. There are always two solutions for $`𝒢_{11}`$, one holomorphic, another non-holomorphic. The second one leads, via Gauss law, to the eigenvalue distribution. The shape of the “coastline” bordering the “sea” of complex eigenvalues is determined by the matching conditions for the two solutions, i.e. it is determined by imposing on the non-holomorphic solution for $`b^2`$ the equation $`b^2=0`$. For other important features of the explained above diagrammatic method and more complex examples see .

## 4 Product of Two Complex Matrices

We will now move toward the main goal of this paper, namely the evaluation of the spectral properties of the random matrix products of the form (5) – with zero drift and normalized variance:

$`Y(T)=\underset{M\mathrm{}}{lim}Y_M\underset{M\mathrm{}}{lim}[(1+\sqrt{{\displaystyle \frac{T}{M}}}X_1)(1+\sqrt{{\displaystyle \frac{T}{M}}}X_2)\mathrm{}`$
$`\mathrm{}(1+\sqrt{{\displaystyle \frac{T}{M}}}X_M)]`$ (67)
In principle the methods of section 3 could be used to study such random matrices, however the nonlinear product structure would make the resulting diagrammatic rules exceedingly difficult to control. Fortunately, for the purposes of calculating the Green’s function, and hence the eigenvalue density, the problem may be linearized at the cost of increasing the size of the matrices. In this section we will show how the method works for $`Y_2`$ – the product of two matrices, while in section 5 we will consider the general case including the limit $`M\mathrm{}`$. We therefore have to calculate

$`Y_2=\left(1+\sqrt{{\displaystyle \frac{\tau }{2}}}X_1\right)\left(1+\sqrt{{\displaystyle \frac{\tau }{2}}}X_2\right)A_1A_2`$ (68)
where $`X_1`$ and $`X_2`$ belong to identical independent Gaussian Complex Ensembles (iiGCE=Girko-Ginibre). Since the matrices are non-hermitian we have to use the formalism of section 3 with a conjugate ‘antiholomorphic’ copy. ¿From now on in the rest of the paper we set the regulator $`ϵ`$ to zero. In order to linearize the product structure we use a trick: we introduce a further doubling, thus constructing an auxiliary $`4N`$ by $`4N`$ Green’s function

$`𝒢(w)`$ $`=`$ $`\left[\left(\begin{array}{cccc}𝐰& 0& 0& 0\\ 0& 𝐰& 0& 0\\ 0& 0& \overline{𝐰}& 0\\ 0& 0& 0& \overline{𝐰}\end{array}\right)\left(\begin{array}{cccc}0& A_1& 0& 0\\ A_2& 0& 0& 0\\ 0& 0& 0& A_2^{}\\ 0& 0& A_1^{}& 0\end{array}\right)\right]_{4N\times 4N}^1`$ (77)
$`=`$ $`\left[\left(\begin{array}{cccc}𝐰& \mathrm{𝟏}& 0& 0\\ \mathrm{𝟏}& 𝐰& 0& 0\\ 0& 0& \overline{𝐰}& \mathrm{𝟏}\\ 0& 0& \mathrm{𝟏}& \overline{𝐰}\end{array}\right)\sqrt{{\displaystyle \frac{\tau }{2}}}\left(\begin{array}{cccc}0& X_1& 0& 0\\ X_2& 0& 0& 0\\ 0& 0& 0& X_2^{}\\ 0& 0& X_1^{}& 0\end{array}\right)\right]_{4N\times 4N}^1`$ (86)
$``$ $`\left[𝓦\sqrt{\tau /2}𝒳\right]^1`$ (87)
In the above we first separated the “deterministic” part from the “random one” (second line) and in the last equality we introduced the notation similar to the hermitian and non-hermitian cases analyzed in the previous sections.

In analogy to (42), we define now a block-trace operation $`\text{tr}_{B4}`$, where we trace each $`N`$ by $`N`$ block of the $`4N`$ by $`4N`$ matrix $`𝒢(w)`$ separately. In such a way, we obtain a 4 by 4 auxiliary Green’s function

$`g(w)\left(\begin{array}{cccc}g_{11}& g_{12}& g_{1\overline{1}}& g_{1\overline{2}}\\ g_{21}& g_{22}& g_{2\overline{1}}& g_{2\overline{2}}\\ g_{\overline{1}1}& g_{\overline{1}2}& g_{\overline{1}\overline{1}}& g_{\overline{1}\overline{2}}\\ g_{\overline{2}1}& g_{\overline{2}2}& g_{\overline{2}\overline{1}}& g_{\overline{2}\overline{2}}\end{array}\right)_{4\times 4}={\displaystyle \frac{1}{N}}\text{tr}_{B4}𝒢(w)`$ (92)
The advantage of this function lies in the fact, that the eigenvalues of the product $`A_1A_2`$ coincide with the squares of the eigenvalues of the $`2N\times 2N`$ block matrix

$$\left(\begin{array}{cc}0& A_1\\ A_2& 0\end{array}\right)$$
(93)
This is due to the off-diagonal structure of the sub-blocks in (93). For a general discussion for an arbitrary product see section 5.

As before we define a 4 by 4 matrix of self-energies $`\stackrel{~}{\mathrm{\Sigma }}_{ij}`$

$`g(w)=(𝒲\stackrel{~}{\mathrm{\Sigma }})^1`$ (94)
where $`𝒲`$ is the result of block-tracing ’bold’ $`𝓦`$.

Similarly as in the Girko-Ginibre case, we can now diagrammatically analyze the content of the matrix $`\stackrel{~}{\mathrm{\Sigma }}`$. We see that only the “double line” propagators $`X_1X_1^{}_{X_1}`$ and $`X_2X_2^{}_{X_2}`$ are different from zero. Therefore from all of the 16 elements of matrix $`\stackrel{~}{\mathrm{\Sigma }}`$ only four are nonzero. Moreover, due to identical measures of $`X_1`$ and $`X_2`$ variables, we get

$`\mathrm{\Sigma }_{1\overline{1}}`$ $`=`$ $`\mathrm{\Sigma }_{2\overline{2}}=\alpha g_{\overline{1}1}=\alpha g_{\overline{2}2}`$
$`\mathrm{\Sigma }_{\overline{1}1}`$ $`=`$ $`\mathrm{\Sigma }_{\overline{2}2}=\alpha g_{1\overline{1}}=\alpha g_{2\overline{2}}`$ (95)
where the factor $`\alpha =\tau /2`$ comes from the propagator in the corresponding Schwinger-Dyson equation of the type discussed in the previous chapter.

Filling the matrix $`\stackrel{~}{\mathrm{\Sigma }}`$ with entries (95) we arrive at the 4 by 4 matrix equation for the elements of the Green’s function:

$`\left(\begin{array}{cccc}g_{11}& g_{12}& g_{1\overline{1}}& g_{1\overline{2}}\\ g_{21}& g_{22}& g_{2\overline{1}}& g_{2\overline{2}}\\ g_{\overline{1}1}& g_{\overline{1}2}& g_{\overline{1}\overline{1}}& g_{\overline{1}\overline{2}}\\ g_{\overline{2}1}& g_{\overline{2}2}& g_{\overline{2}\overline{1}}& g_{\overline{2}\overline{2}}\end{array}\right)_{4\times 4}=`$ $`[\left(\begin{array}{cccc}𝐰& \mathrm{𝟏}& 0& 0\\ \mathrm{𝟏}& 𝐰& 0& 0\\ 0& 0& \overline{𝐰}& \mathrm{𝟏}\\ 0& 0& \mathrm{𝟏}& \overline{𝐰}\end{array}\right)`$ (109)
$`\left(\begin{array}{cccc}0& 0& \alpha g_{\overline{1}1}& 0\\ 0& 0& 0& \alpha g_{\overline{1}1}\\ \alpha g_{1\overline{1}}& 0& 0& 0\\ 0& \alpha g_{1\overline{1}}& 0& 0\end{array}\right)]^1`$
As in the Ginibre case, we solve it for $`g_{\overline{1}1}`$ and $`g_{1\overline{1}}`$. Then, using the solution, we calculate $`g_{11}`$. Coming back to original variables with $`w^2=z`$ we set $`G(z,\overline{z})=g_{11}(w)`$. The eigenvalue distribution then follows from

$`\rho (x,y)={\displaystyle \frac{1}{\pi }}{\displaystyle \frac{w}{z}}_{\overline{z}}G(z,\overline{z})`$ (110)
This is a consequence of the relation between the eigenvalues of the block matrix (93) and the eigenvalues of the product $`A_1A_2`$. A general derivation is outlined in the following section.

A little algebra (basically an inversion of the 4 by 4 matrix) shows, that, as before, we get two solutions for $`b^2=g_{1\overline{1}}g_{\overline{1}1}`$$`=g_{2\overline{2}}g_{\overline{2}2}`$: (i) the holomorphic one (corresponding to $`b^2=0`$, i.e. all $`\mathrm{\Sigma }_{ij}`$ equal to zero); (ii) the non-holomorphic one, given by the equation

$`{\displaystyle \frac{\tau }{2}}\left(|w|^2\alpha ^2b^2+1\right)=(\alpha ^2b^21+w^2)(\alpha ^2b^21+\overline{w}^2)\alpha ^2b^2(w+\overline{w})^2`$ (111)
The boundary of the eigenvalue support is given, as before, by the above equation with $`b^2`$ set to 0. It is therefore a fourth order conchoid-like planar curve, which in polar coordinates is given by

$`{\displaystyle \frac{\tau }{2}}(1+r)=r^2+12r\mathrm{cos}\varphi `$ (112)
On Fig. 4 we show the shape of the nonholomorphic island for several sample “evolution times” $`\tau `$.

The Green’s function for the product follows from

$$g_{11}(w)=\frac{w(\overline{w}^21)\alpha ^2b^2\overline{w}}{det}$$
(113)
where

$$det=(\alpha ^2b^21+w^2)(\alpha ^2b^21+\overline{w}^2)\alpha ^2b^2(w+\overline{w})^2$$
(114)
## 5 Product of Arbitrary Number of Complex Matrices

We consider now the product of an arbitrary number of matrices. To see the pattern, let us look briefly at the case of three matrices

$`Y_3`$ $`=`$ $`\left(1+\sqrt{{\displaystyle \frac{\tau }{3}}}X_1\right)\left(1+\sqrt{{\displaystyle \frac{\tau }{3}}}X_2\right)\left(1+\sqrt{{\displaystyle \frac{\tau }{3}}}X_3\right)`$ (115)
$``$ $`A_1A_2A_3`$
where $`X_1,X_2,X_2`$ again belong to the Girko-Ginibre ensemble.

Our approach again follows from a very simple exact relation between the eigenvalues of a product of $`N\times N`$ matrices $`A_1A_2A_3`$ and the eigenvalues of a block matrix

$$_3=\left(\begin{array}{ccc}0& A_1& 0\\ 0& 0& A_2\\ A_3& 0& 0\end{array}\right)_{3N\times 3N}$$
(116)
Indeed if $`\{\lambda _i\}`$ are the eigenvalues of $`A_1A_2A_3`$ then the eigenvalues of the block matrix are $`\{\lambda _i^{1/3},\lambda _i^{1/3}e^{\frac{2\pi i}{3}},\lambda _i^{1/3}e^{\frac{4\pi i}{3}}\}`$. This is an exact relation for any $`N`$ and follows from the relation between the resolvents (here the matrices $`A_i`$ are of finite size and fixed i.e. no averaging)

$$G__3(w)=\frac{1}{3N}\text{tr}\frac{1}{w_3}G_{A_1A_2A_3}(z)=\frac{1}{N}\text{tr}\frac{1}{zA_1A_2A_3}$$
(117)
namely

$$wG__3(w)=zG_{A_1A_2A_3}(w^3z)$$
(118)
This is due to the cyclic structure of the block matrix (116). Obviously, only multiplicities of the cubic powers of $`_3`$ contribute to the trace.

The relation between the eigenvalues now follows from the location of the finite number of poles of both functions.

We can thus safely calculate the eigenvalue density for the block matrix $`_3`$ using the diagrammatic methods and then extract the density for the product through

$$\rho _{A_1A_2\mathrm{}A_M}(z,\overline{z})=\frac{1}{\pi }_{\overline{z}}G_{A_1A_2\mathrm{}A_M}=\frac{1}{M}\frac{w\overline{w}}{z\overline{z}}\rho __M(w,\overline{w})$$
(119)
where $`M=3`$ and $`w^M=z`$ and $`\rho __M(w,\overline{w})=\frac{1}{\pi }_{\overline{w}}G_{}(w,\overline{w})`$.

To extract only the moments involving powers of the triples $`A_1A_2A_3`$, we construct an auxiliary $`6N`$ by $`6N`$ Green’s function, triplicating the $`A_1A_2A_3`$ products by rewriting them as cyclic block matrices:

$`𝒢(w)=`$ $`[\left(\begin{array}{cccccc}𝐰& \mathrm{𝟏}& 0& 0& 0& 0\\ 0& 𝐰& \mathrm{𝟏}& 0& 0& 0\\ \mathrm{𝟏}& 0& 𝐰& 0& 0& 0\\ 0& 0& 0& \overline{𝐰}& 0& \mathrm{𝟏}\\ 0& 0& 0& \mathrm{𝟏}& \overline{𝐰}& 0\\ 0& 0& 0& 0& \mathrm{𝟏}& \overline{𝐰}\end{array}\right)`$ (133)
$`\sqrt{{\displaystyle \frac{\tau }{3}}}\left(\begin{array}{cccccc}0& X_1& 0& 0& 0& 0\\ 0& 0& X_2& 0& 0& 0\\ X_3& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& X_3^{}\\ 0& 0& 0& X_1^{}& 0& 0\\ 0& 0& 0& 0& X_2^{}& 0\end{array}\right)]_{6N\times 6N}^1`$
where again we separated the “deterministic” part from the “random one”. Introducing now the block-trace operation $`\text{tr}_{B6}`$, we obtain a 6 by 6 auxiliary Green’s function $`g(w)=\text{tr}_{B6}𝒢(w)`$.

The generalization for arbitrary $`M`$ is now straightforward. For

$`Y_M=\left(1+\sqrt{{\displaystyle \frac{\tau }{M}}}X_1\right)\left(1+\sqrt{{\displaystyle \frac{\tau }{M}}}X_2\right)\mathrm{}\left(1+\sqrt{{\displaystyle \frac{\tau }{M}}}X_M\right)`$ (134)
we define an auxiliary $`2MN`$ by $`2MN`$ Green’s function of the form

$`𝒢(w)=\left(\begin{array}{cc}𝓦& \sqrt{\frac{\tau }{M}}𝒳\\ \sqrt{\frac{\tau }{M}}𝒳^{}& 𝓦^{}\end{array}\right)^1`$ (137)
where the blocks

$`𝓦=\left(\begin{array}{ccccc}𝐰& \mathrm{𝟏}& 0& \mathrm{}& 0\\ 0& 𝐰& \mathrm{𝟏}& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 𝐰& \mathrm{𝟏}\\ \mathrm{𝟏}& 0& \mathrm{}& 0& 𝐰\end{array}\right)_{MN\times MN}`$ (143)
and

$`𝒳=\left(\begin{array}{ccccc}0& X_1& 0& \mathrm{}& 0\\ 0& 0& X_2& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 0& X_{M1}\\ X_M& 0& \mathrm{}& 0& 0\end{array}\right)_{MN\times MN}`$ (149)
are themselves $`NM`$ by $`NM`$ matrices, i.e. each of the listed elements in $`𝓦`$ and in $`𝒳`$ is itself an $`N`$ by $`N`$ matrix, either diagonal, denoted by a bold symbol, or a random entry $`X_i`$, otherwise an $`N`$ by $`N`$ block of zeroes.

We take now a block-trace operation $`\text{tr}_{BM}`$, where we trace each $`N`$ by $`N`$ block of the $`2MN`$ by $`2MN`$ matrix $`𝒢(w)`$ separately. In such a way, we obtain an $`2M`$ by $`2M`$ auxiliary Green’s function

$`g(w)\left(\begin{array}{cccccc}g_{11}& \mathrm{}& g_{1M}& g_{1\overline{1}}& \mathrm{}& g_{1\overline{M}}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ g_{M1}& \mathrm{}& g_{MM}& g_{M\overline{1}}& \mathrm{}& g_{M\overline{M}}\\ g_{\overline{1}1}& \mathrm{}& g_{\overline{1}M}& g_{\overline{1}\overline{1}}& \mathrm{}& g_{\overline{1}\overline{M}}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ g_{\overline{M}1}& \mathrm{}& g_{\overline{M}M}& g_{\overline{M}\overline{1}}& \mathrm{}& g_{\overline{M}\overline{M}}\end{array}\right)_{2M\times 2M}={\displaystyle \frac{1}{N}}\mathrm{Tr}_{BM}𝒢(w)`$ (156)
As before we define a $`2M`$ by $`2M`$ matrix of self-energies $`\stackrel{~}{\mathrm{\Sigma }}_{ij}`$

$`g(w)=\left[\left(\begin{array}{cc}𝒲& 0\\ 0& 𝒲^{}\end{array}\right)\stackrel{~}{\mathrm{\Sigma }}\right]^1`$ (159)
where $`𝒲`$ is a result of block-tracing $`𝓦`$. We can now diagrammatically analyze the content of the matrix $`\stackrel{~}{\mathrm{\Sigma }}`$. As previously, only “double line” propagators $`X_iX_i^{}_{X_i}`$ for $`i=1,\mathrm{},M`$ are different from zero. Therefore from all of the $`4M^2`$ elements of the matrix $`\stackrel{~}{\mathrm{\Sigma }}`$ only $`2M`$ are different from zero. Due to the symmetries we get

$`\mathrm{\Sigma }_{1\overline{1}}`$ $`=`$ $`\mathrm{}=\mathrm{\Sigma }_{M\overline{M}}=\alpha g_{\overline{1}1}=\alpha g_{\overline{M}M}\alpha g`$
$`\mathrm{\Sigma }_{\overline{1}1}`$ $`=`$ $`\mathrm{}=\mathrm{\Sigma }_{\overline{M}M}=\alpha g_{1\overline{1}}=\alpha g_{1\overline{M}}\alpha \stackrel{~}{g}`$ (160)
where the factor $`\alpha =\tau /M`$ comes from the propagator in the corresponding set of Schwinger-Dyson equations as in the previous chapter.

Inserting now the matrix $`\stackrel{~}{\mathrm{\Sigma }}`$ with entries (160) we arrive at a $`2M`$ by $`2M`$ matrix equation for the elements of the Green’s function. First, we solve it for $`g`$ and $`\stackrel{~}{g}`$. Second, using the solutions, we calculate $`g_{11}`$. Third, using (119) we get an explicit equation for the spectral density

$`\rho (x,y)={\displaystyle \frac{1}{\pi }}{\displaystyle \frac{1}{M}}{\displaystyle \frac{w\overline{w}}{z\overline{z}}}_{\overline{w}}g_{11}(w,\overline{w})`$ (161)
where $`z=x+iy`$. For arbitrary $`M`$, the algebra is rather involved. Luckily, both the block and the cyclic structure of the main entries $`𝒲`$ and $`\stackrel{~}{\mathrm{\Sigma }}`$ make taking an inverse of the $`2M`$ by $`2M`$ matrix possible. The inverse of a cyclic matrix is a cyclic matrix, and its explicit form can be obtained from the solution of an associated recurrence relation, e.g. using the transfer matrix techniques which we will now proceed to do.

At this moment, it is useful to introduce the notation

$`1+|w|^2\alpha ^2b^2`$ $`=`$ $`2|w|\mathrm{cosh}u`$
$`|w|/w`$ $`=`$ $`\delta `$ (162)
The inversion needed for the Green’s function reduces to finding explicit forms of the cyclic $`M`$ by $`M`$ matrices $`𝒞`$, $`\stackrel{~}{𝒞}`$ defined as

$`𝒞=[𝒲𝒲^{}\alpha ^2b^2\mathrm{𝟏}]^1\stackrel{~}{𝒞}=[𝒲^{}𝒲\alpha ^2b^2\mathrm{𝟏}]^1`$ (163)
where $`b^2=g\stackrel{~}{g}`$. $`\stackrel{~}{𝒞}`$ is just the transpose, or equivalently, the complex conjugate of $`𝒞`$. The elements $`c_0,\mathrm{}c_{M1}`$ of the first row of the cyclic matrix $`𝒞`$ fulfill the recurrence relation

$$\left(\begin{array}{c}c_{k+1}\\ c_k\end{array}\right)=\left(\begin{array}{cc}2\delta \mathrm{cosh}u& \delta ^2\\ 1& 0\end{array}\right)\left(\begin{array}{c}c_k\\ c_{k1}\end{array}\right)$$
(164)
with the ‘boundary conditions’

$`\left(\begin{array}{c}c_1\\ c_0\end{array}\right)=\left(\begin{array}{cc}2\delta \mathrm{cosh}u& \delta ^2\\ 1& 0\end{array}\right)\left(\begin{array}{c}c_M\\ c_{M1}\end{array}\right){\displaystyle \frac{1}{w}}\left(\begin{array}{c}1\\ 0\end{array}\right)`$ (173)
The above recurrence can be easily solved using transfer matrix techniques, i.e. by diagonalizing the 2 by 2 transfer matrix in (164) and returning at the end of the calculations to the original basis. An explicit solution for the cyclic string $`c_k`$, for $`k=0,1,..,.M1`$ reads

$`c_k={\displaystyle \frac{1}{w(\mathrm{\Lambda }_+\mathrm{\Lambda }_{})}}\left({\displaystyle \frac{\mathrm{\Lambda }_+^k}{\mathrm{\Lambda }_+^M1}}{\displaystyle \frac{\mathrm{\Lambda }_{}^k}{\mathrm{\Lambda }_{}^M1}}\right)`$ (174)
where $`\mathrm{\Lambda }_\pm =\delta e^{\pm u}`$ denote the eigenvalues of the transfer matrix.

We may now express our key quantities for arbitrary $`M`$ in terms of (174):

$`g_{11}(w)`$ $`=`$ $`[𝒲^{}𝒞]_{11}\overline{w}c_0c_{M1}^{}`$
$`g`$ $`=`$ $`\alpha g[\stackrel{~}{𝒞}]_{11}\alpha gc_0^{}`$
$`\stackrel{~}{g}`$ $`=`$ $`\alpha \stackrel{~}{g}[𝒞]_{11}\alpha \stackrel{~}{g}c_0`$ (175)
where  denotes complex conjugation. As always in the case of the complex spectra, we have two solutions: trivial (holomorphic), corresponding to $`g=\stackrel{~}{g}=0`$, and the nontrivial (non-holomorphic), corresponding to the case $`b^2g\stackrel{~}{g}0`$.

Inserting the nonholomorphic solution into $`g_{11}`$, redefining $`G(z,\overline{z})=g_{11}(w)`$ with $`w^M=z`$, and using (119) yields the eigenvalue density for $`Y_M`$. Since $`_{\overline{w}}\overline{z}=M\overline{z}/\overline{w}`$, the result for the final spectral function can be further simplified

$`\rho (x,y)={\displaystyle \frac{1}{\pi }}{\displaystyle \frac{w}{z}}_{\overline{z}}G(z,\overline{z})`$ (176)
Explicit formulas for arbitrary $`M`$ could be easily constructed e.g. by running any symbolic program using our solution. The shape of the boundary is given by the condition that the nonholomorphic solution meets the holomorphic one on the complex plane.

After solving the problem for arbitrary $`M`$, we can now address the original problem, i.e. the solution for a continuous product of matrices, corresponding to the limit $`M\mathrm{}`$. This means, that we have to perform a careful limiting procedure in the above formulas. Introducing $`U=Mu`$, expanding in $`M`$, and using the polar parametrization $`z=r\mathrm{exp}i\varphi `$, we get the solution

$`G(z,\overline{z})={\displaystyle \frac{1}{\tau }}\mathrm{log}r+{\displaystyle \frac{1}{2}}{\displaystyle \frac{iU}{\tau }}{\displaystyle \frac{\mathrm{sin}\varphi }{\mathrm{sinh}U}}`$ (177)
where $`U=\sqrt{\mathrm{log}^2r\tau ^2b^2}`$ with $`b^2`$ fulfilling the equation

$`\mathrm{cos}\varphi =\mathrm{cosh}U{\displaystyle \frac{\tau }{2U}}\mathrm{sinh}U`$ (178)
The set of equations (177,178) completes the solution of the continuous product of large complex matrices with Gaussian disorder. We would like to note, that due to the limiting procedure $`M\mathrm{}`$ differentiation of the Green’s functions has to be done with care. Taking the large $`M`$ limit of (176) yields

$`\rho (x,y)={\displaystyle \frac{1}{\pi }}{\displaystyle \frac{1}{z}}_{\overline{z}}G(z,\overline{z})`$ (179)
As a cross-check of our calculations we verified that the spectral function obtained in this way is real.

Substituting $`b^2=0`$ in (178), we get the equation for the shape of the support of the eigenvalues for arbitrary $`\tau `$. Note that the shape of the boundary depends on the variable $`\mathrm{log}|z|`$, reflecting in the spectrum the multiplicative (geometric) feature of the stochastic process induced by the multiplication of random matrices. Explicit solution for the boundary reads

$`{\displaystyle \frac{\tau }{2}}\left({\displaystyle \frac{r^21}{\mathrm{ln}r}}\right)=r^2+12r\mathrm{cos}\varphi `$ (180)
It is rewarding to compare solution (180) to the solution of the conchoid-like boundary in the case of two complex matrices (112). In the limit of infinitely many products, the shape of the boundary acquires a new symmetry – (180) is invariant under inversion operation $`r1/r`$. It is visible also from the general form (178), since the radial dependence for $`b^2=0`$ is a function of $`|\mathrm{log}r|`$ only. This symmetry is responsible for the appearance of a structural change of the spectrum – provided $`\tau `$ is sufficiently large, two boundaries, related by inversion in the radial variable, appear. We describe this structural change of the complex spectrum as a ”topological phase transition”. In Section 8, we study numerically the spectral density predicted by our formulas, in particular, the “diffusion” of the shape of the boundary caused by evolution with the “time” parameter $`\tau `$. We also confirm numerically the critical behavior for certain $`\tau _{\mathrm{crit}}`$.

## 6 Product of Two Hermitian Matrices

In this section we will consider the spectral distribution of the product

$`Y_2=\left(1+\sqrt{{\displaystyle \frac{\tau }{2}}}H_1\right)\left(1+\sqrt{{\displaystyle \frac{\tau }{2}}}H_2\right)`$ (181)
where $`H_1,H_2`$ are not arbitrary complex matrices, but they are identical independent hermitian Gaussian matrices (iiGUE). Contrary to naive expectations, the above problem is more subtle then the analogous problem of the two complex Gaussian matrices. First, let us note, that the product of two hermitian matrices is not, in general, a hermitian matrix, so the spectrum of the strings of hermitian matrices develops complex eigenvalues.

Second, the case of only two hermitian matrices of the considered type turns out to be in some sense singular – for certain values of the evolution parameter $`\tau `$, the a priori complex eigenvalues get localized on the real axis. This phenomenon reminds a bit the non-hermitian localization in the Hatano-Nelson model . Third, the problem is algebraically more involved since there will be more nonzero propagators than in the complex case.

Since even a product of two hermitian matrices can develop a priori a complex spectrum, we have to use the same type of auxiliary Green’s function as in the complex Gaussian case. The construction parallels exactly the scheme for two complex matrices considered previously, with the matrices $`X_i`$ and $`X_i^{}`$ replaced both by $`H_i`$. The first new feature appears when we construct the “Schwinger-Dyson” equations for the 4 by 4 self-energy matrix $`\stackrel{~}{\mathrm{\Sigma }}`$. Since the matrices are hermitian ($`H_i=H_i^{}`$), additional nonzero contractions appear in the matrix $`\stackrel{~}{\mathrm{\Sigma }}`$

$`\stackrel{~}{\mathrm{\Sigma }}=\alpha \left(\begin{array}{cccc}0& g_{12}& g_{1\overline{1}}& 0\\ g_{21}& 0& 0& g_{2\overline{2}}\\ g_{\overline{1}1}& 0& 0& g_{\overline{1}\overline{2}}\\ 0& g_{\overline{2}2}& g_{\overline{2}\overline{1}}& 0\end{array}\right)\alpha \left(\begin{array}{cccc}0& h& g& 0\\ h& 0& 0& g\\ \stackrel{~}{g}& 0& 0& \stackrel{~}{h}\\ 0& \stackrel{~}{g}& \stackrel{~}{h}& 0\end{array}\right)`$ (190)
where again the factor $`\alpha =\tau /2`$ comes from the propagator, $`h`$ and $`\stackrel{~}{h}`$ are the new entries and we have exploited the symmetries of the problem (e.g. $`g_{1\overline{1}}=g_{2\overline{2}}`$ etc).

The solution of the 4 by 4 matrix equation for the 4 by 4 Green’s function

$`𝒢(w)=(𝒲\stackrel{~}{\mathrm{\Sigma }})^1`$ (191)
is more involved, since we have to solve it not only for the “gaps” $`g`$ and $`\stackrel{~}{g}`$ but also for additional gaps $`h`$ and $`\stackrel{~}{h}`$. Only then we can calculate $`g_{11}`$, to recover the wanted eigenvalue distribution in the same way as in the complex case (110).

Introducing $`d^2=(1+\alpha h)(1+\alpha \stackrel{~}{h})`$ and $`b^2=g\stackrel{~}{g}`$, we obtain two (modulo their tilded partners) gap equations

$`g`$ $`=`$ $`{\displaystyle \frac{\alpha g}{det}}(\alpha ^2b^2|w|^2d^2)`$
$`h`$ $`=`$ $`{\displaystyle \frac{1}{det}}(\alpha ^2b^2(1+\alpha \stackrel{~}{h})+(1+\alpha h)(\overline{w}^2(1+\alpha \stackrel{~}{h})^2)`$ (192)
coupled to the needed $`g_{11}`$

$`g_{11}={\displaystyle \frac{1}{det}}(\alpha ^2b^2\overline{w}+w(\overline{w}^2(1+\alpha \stackrel{~}{h})^2)`$ (193)
where $`det`$ is the determinant of the matrix $`𝒲\stackrel{~}{\mathrm{\Sigma }}`$ and reads explicitly

$`det=\alpha ^4b^4+(w^2(1+\alpha h)^2)(\overline{w}^2(1+\alpha \stackrel{~}{h})^2)2\alpha ^2b^2(d^2+|w|^2)`$ (194)
Numerical simulations show that the eigenvalue distribution for the product of two matrices is concentrated on an interval on the real axis for $`\tau <1/2`$ and only then when the edge of the interval reaches the origin, do the eigenvalues start to spread into the complex plane (see figures 5 and 6). Let us show how to find the critical value $`\tau =1/2`$ from the above formulas based on the above numerical observation.

We thus have to find the point when the boundary of the two dimensional island hits the origin $`w=0`$. Then the equations (192) and (194) simplify to

$`h`$ $`=`$ $`{\displaystyle \frac{1}{1+\alpha h}}`$ (195)
$`1`$ $`=`$ $`\alpha h\stackrel{~}{h}`$ (196)
where we set $`b^2=g\stackrel{~}{g}=0`$. Then it is easy to verify that this set of equations is satisfied for $`\alpha =1/4`$ which is equivalent to $`\tau =1/2`$.

## 7 Products of Arbitrary Number of Hermitian Matrices

Finally, we consider an arbitrary number of matrices belonging to iiGUE.

$`Y_M=\left(1+\sqrt{{\displaystyle \frac{\tau }{M}}}H_1\right)\left(1+\sqrt{{\displaystyle \frac{\tau }{M}}}H_2\right)\mathrm{}\left(1+\sqrt{{\displaystyle \frac{\tau }{M}}}H_M\right)`$ (197)
As before, the construction of the auxiliary Green’s function parallels the complex case, up to the point when we have to construct explicitly the matrix of self-energies $`\stackrel{~}{\mathrm{\Sigma }}`$, corresponding to the set of $`4M^2`$ Schwinger-Dyson equations. Generalization of the results from the previous sections leads to the explicit form of the auxiliary Green’s function. Introducing

$`𝒢(w)=\left(\begin{array}{cc}𝒲& 𝒮\\ \stackrel{~}{𝒮}& 𝒲^{}\end{array}\right)^1`$ (200)
with blocks

$`𝒲=\left(\begin{array}{ccccc}w& 1\alpha h& 0& \mathrm{}& 0\\ 0& w& 1\alpha h& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& w& 1\alpha h\\ 1\alpha h& 0& \mathrm{}& 0& w\end{array}\right)_{M\times M}`$ (206)
and

$`𝒮=\alpha \left(\begin{array}{ccccc}g& 0& 0& \mathrm{}& 0\\ 0& g& 0& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& g& 0\\ 0& 0& \mathrm{}& 0& g\end{array}\right)_{M\times M}`$ (212)
the consistent set of equations comes as

$`g_{11}`$ $`=`$ $`𝒢_{1,1}`$
$`h`$ $`=`$ $`𝒢_{1,M}`$
$`\overline{h}`$ $`=`$ $`𝒢_{M+1,M+2}`$
$`g`$ $`=`$ $`𝒢_{1,M+1}`$
$`\overline{g}`$ $`=`$ $`𝒢_{M+1,1}`$ (213)
i.e. requires an explicit inversion of the $`2M`$ by $`2M`$ matrix. The algebra is more complicated since, like in the case of two hermitian matrices, we have now two types of gap equations (for $`g`$ and $`h`$ and their tilded partners). Using the solutions, we calculate $`g_{11}`$, then, we get an explicit equation for the spectral density (119). Finally, we perform the limiting procedure. Technically, one has to repeat the transfer matrix technique introduced in section 5 for the complex case. The algebra is lengthy, but due to the cyclic nature of the blocks of matrix $``$ the results form a rather surprisingly simple set of equations. Returning to the original variable $`z=r\mathrm{exp}(i\varphi )`$ we get

$`G(z,\overline{z})={\displaystyle \frac{1}{2\tau }}\mathrm{log}r+{\displaystyle \frac{3}{4}}{\displaystyle \frac{iU}{\tau }}\mathrm{sin}\psi /\mathrm{sinh}U`$ (214)
and the gap equations read

$`h={\displaystyle \frac{1}{2\tau }}\mathrm{log}r{\displaystyle \frac{1}{4}}{\displaystyle \frac{iU}{\tau }}\mathrm{sin}\psi /\mathrm{sinh}U=G1`$ (215)
and

$`\mathrm{cos}\psi =\mathrm{cosh}U{\displaystyle \frac{\tau }{2U}}\mathrm{sinh}U`$ (216)
where $`\psi =\varphi +U\frac{\mathrm{sin}\psi }{\mathrm{sinh}U}`$ and $`U=\sqrt{(\frac{\mathrm{log}r}{2}+\tau /4)^2\tau ^2g\stackrel{~}{g}}`$. The above set of equations completes the solution of the continuous product of matrices with Gaussian hermitian disorder. Spectral density follows from (179). Note, that the dependence on $`\mathrm{log}|z|`$ reflects the multiplicative character of the Brownian motion. Again, the boundary has also an inversion-type symmetry, but now it is of the form $`|z|\mathrm{exp}(\tau )/|z|`$. Due to the non-zero gap $`h`$, the boundary in the hermitian case is given by the pair of coupled transcendental equations and can be presented numerically only. The appearance of the aforementioned symmetry is responsible for the ”topological phase transition” in the spectrum of the product of infinitely many matrices, similar to the one observed for the product of infinitely many complex matrices.

In Section 8, we study numerically the spectral density predicted by our formulas, in particular, the “diffusion” of the shape of the boundary caused by evolution with the “time” parameter $`\tau `$. We confirm the appearance of a “phase transition” for certain $`\tau _{\mathrm{crit}}`$ in the hermitian case as well.

## 8 Numerical Tests

In this section we present the results of numerical simulations, confirming our analytical predictions. We start from the complex case. Figure 7 shows the evolution of the boundary for several sample evolution times $`\tau =0.5,1,1.5,2,3,4,5,6,8,12`$. To present the whole set of the boundaries on the same figure, each boundary is rescaled by a corresponding factor $`\mathrm{exp}(\tau /2)`$. We would like to note here that the effect of rescaling is equivalent to the non-rescaled process with drift $`\mu =\mathrm{𝟏}`$. One can see how the original ellipse-like shape (innermost figure) evolves through a twisted-like shape to the set of double-ring structures. The inner ring, always containing the origin, is so small on the scale of the Fig. 7 that is not visible. At $`\tau =4`$, corresponding to the curve with an inner loop, we observe a topological phase transition. The support of the spectrum is no longer simply connected, for $`\tau 4`$ it is annulus-like, i.e. eigenvalues are expelled from the central region. For even larger times, the outside rim of the annulus approaches the circle, and the inner one shrinks to the point $`z=0`$ in the $`\tau =\mathrm{}`$ limit. Indeed for large $`\tau `$ the radius of the inner ring behaves like

$$r(\tau )e^{\frac{\tau }{2}}$$
(217)
The outer boundary then forms approximately a circle with the radius $`e^{\tau /2}`$. To visualize the repulsion of the eigenvalues from the region around $`z=0`$ we performed higher statistics simulations for $`\tau =4`$ and $`\tau =5`$, corresponding to Figs 9,10, respectively. Note that the eigenvalue distributions are quite high for small $`z`$ and the size of the inner ring is indeed very small. Therefore it is quite difficult to observe numerically the exact exclusion of eigenvalues from the marked circle.

Figure 8 shows the comparison of numerically generated spectra versus the analytical prediction of the shape of the support of eigenvalues of the ”evolution operator”. Again, the same rescaling by $`\mathrm{exp}(\tau /2)`$ was applied.

Figure 11 shows the evolution of the boundary for several sample evolution times $`\tau =0.5,1,1.5,2,3,4,5,6,8,12`$ in the case of GUE. Again, to present the whole set of boundaries on the same figure, each boundary was rescaled by a corresponding factor $`\mathrm{exp}(\tau /2)`$. One can observe how the almost real (for small times $`\tau `$) long-cigar shaped spectrum evolves into a broader shape, developing again singular behavior at $`\tau =4`$ at the leftmost-edge of the spectrum. Again, the spectrum develops a topological transition, although it is not as clearly visible in the figure as in the complex case. The reason is that inversion symmetry has an additional exponential suppression factor, i.e. $`re^\tau /r`$. For very large times, the outside rim of the spectrum approaches the circle, and the inner one shrinks to a point in the $`\tau \mathrm{}`$ limit.

As in the complex case, here we also show a sample simulation (Fig. 12) versus the analytical prediction for the shape of the island.

Finally in figures 5 and 6 above, we demonstrated the ”localization-delocalization” phase transition in the case of the product of two hermitian matrix models of the considered type. Below $`\tau _{\mathrm{crit}}`$, a priori complex spectra condense on the positive real axes. At $`\tau =1/2`$, the eigenvalues start to diffuse onto the $`x0`$ half-plane. For very large times $`\tau `$, eigenvalues start to appear on both sides of the $`x=0`$ axis.

## 9 Stability and Lyapunov Exponents

In this section, we briefly discuss the relation between the spectrum of the “evolution operator” $`Y_M`$ and Lyapunov spectra. In the studies of products of random matrices one of the central issues is the use of limit theorems that would provide an information about the asymptotic limit of the rates of growth and of the spectrum of the product for a sequence of matrices. In particular the limit

$`\underset{M\mathrm{}}{lim}{\displaystyle \frac{1}{M}}\mathrm{ln}|Y_M|\lambda _1`$ (218)
defines the maximum Lyapunov characteristic exponent. The existence of this limit is assured for an infinite product of matrices by the Furstenberg  theorem

$`Y_M={\displaystyle \underset{i=1}{\overset{M}{}}}X(i)`$ (219)
for independent random matrices $`X(i)`$ characterized by the probability distribution function $`d\mu [X]=𝒫[X]d[X]`$. It can be shown under very general assumptions that $`\lambda _1`$ is a non-random, self-averaging quantity, so that the brackets $`\mathrm{}`$ in the above expression can be dropped in almost any realization. The best known example is the use of $`\lambda _1`$ in the description of the DC conductivity of a one dimensional disordered system coupled to two electron reservoirs. The conductivity $`\kappa `$ of the system, given by the Landauer formula yields for large $`M`$

$`\kappa e^{2\lambda _1M}`$ (220)
where $`\lambda _1`$ is the maximal Lyapunov exponent of the product $`Y_M`$. By the Borland conjecture, the inverse of $`\lambda _1`$ gives the localization length measuring the decay of an eigenvector of the system’s Hamiltonian over the chain composed of $`M`$ units.

It turns out that the formalism of generalized Green functions (Section 2 and 3) is particularly suitable for deriving an analytical formula  for $`\lambda _1^1`$ for the one-dimensional disordered wires within the tight-binding Hamiltonian with a site diagonal disorder. Extensions to generalized Lyapunov exponents $`L(q)`$  related to the exponential growth rates of the moments of the matrix product

$`L(q)=\underset{M\mathrm{}}{lim}{\displaystyle \frac{1}{M}}\mathrm{ln}|Y_M|^q`$ (221)
have been also reported for systems described by Markovian evolutionary operators  used in the study of DC conductance statistics of the random dimer model introduced for explaining the exceptionally high electronic transport properties of some conjugated polymer chains.

The eigenvectors of the product $`Y_M^{}Y_M`$ form an orthonormal set of vectors (Lyapunov basis) corresponding to characteristic Lyapunov exponents whose exponentials, $`\mathrm{exp}(\lambda _i)`$ are eigenvalues of the matrix $`|Y_M^{}Y_M|^{1/2M}`$ in the limit of large $`M`$ (Oseledec’s theorem).

The Lyapunov exponents are closely related to the linear stability analysis of the dynamic system $`\dot{x}=F(x)`$ for which small perturbations $`y`$ around a given solution $`x(t)`$ evolve in time according to $`y(t_2)=K(t_2,t_1)y(t_1)`$. They are defined as

$`\lambda _i=\underset{t_2\mathrm{}}{lim}{\displaystyle \frac{1}{t_2t_1}}\mathrm{log}K(t_2,t_1)f_i(t_1)`$ (222)
and for ergodic systems their spectrum describes global properties of the system’s attractor (does not depend on the initial point). In contrast to the Lyapunov spectra, the Lyapunov basis for the ergodic sequence of matrices describes local properties of the PRM system.

Here we would like to underline, that Lyapunov exponents are in general different from the spectrum of eigenvalues obtained in this paper. Since the first Lyapunov exponent is defined as a property of the norm of the operator $`Y`$, it is related to the spectral properties of the $`Y^{}Y`$ operator. This is why the Lyapunov exponents are always real. The spectra studied in this paper correspond to so-called stability exponents, which are directly related to the eigenvalues of the $`Y`$ operator itself. As we have demonstrated above, the eigenvalues of $`Y_M`$ are in general complex numbers and so are the stability exponents:

$`\alpha _i=\underset{M\mathrm{}}{lim}{\displaystyle \frac{1}{M}}\mathrm{ln}\lambda _i`$ (223)
Their real parts coincide with the Lyapunov exponents $`\lambda _i=Re(\alpha _i)`$ only if the matrix $`Y_M`$ can be written in a diagonal form.

The relation between the stability and Lyapunov exponents in a broad class of dynamical systems has been investigated by Goldhirsch et al. . Motivated by mostly numerical results, the authors have conjectured that Lyapunov exponents equal stability ones for the ergodic dynamical systems whose spectrum is non degenerate and the attractor is bounded. The geometrical argument behind that statement is that generically, stretching and shrinking of the relative separation of “trajectories”, starting at a given initial point happen along fixed directions in the phase space. Taken from that perspective further elucidation would be required to clarify the relation between the Lyapunov stability and spectral properties of the generalized stability matrix investigated in this work.

## 10 Summary

We have introduced a natural generalization of the concept of geometric random walk (‘geometric diffusion’) in the space of large, non-commuting matrices. In particular, we proposed a simple diagrammatic method, allowing to calculate spectral properties of the resulting evolution operators as functions of the ”evolution” time $`\tau `$. The construction was explicitly presented in the case of complex Gaussian ensembles and GUE ensembles.

Our method was based on an exact relation between the eigenvalues of a product matrix and the eigenvalues of a certain block matrix. The essential simplification is the fact that the random matrices appear linearly in the block matrix. Although we used the above identification only in the large $`N`$ limit, it could a-priori be used to study the properties of the spectrum of the product matrix also in the microscopic limit, i.e. on the scale of eigenvalue spacing. But then different, non-diagrammatic, methods would have to be employed.

Despite the simplicity of the random matrix models chosen as the building blocks of the matrix-valued diffusion, the spectrum develops a surprising feature, namely a region without eigenvalues appears within the spectrum, thus changing its topological properties. We can thus describe this behavior as a topological phase transition. This points at the appearance of particular repulsion mechanism for sufficiently large evolution times. Our results were based on the simulation based on Gaussian ensembles. We are however tempted to conjecture, that unfolding the spectrum in the vicinity of the point where the topological change occurs may unravel a possible novel universal spectral behavior. We would like to mention, that the critical behavior described here as a ”topological phase transition” appeared only in the limit $`M\mathrm{}`$, as a consequence of the inversion symmetry, acquired in this limit, of the equation describing the boundary of the two-dimensional support of eigenvalues. We would like to mention, that similar, qualitative effect of expulsion of complex eigenvalues from the certain region of complex plane was observed by Feinberg and Zee (first reference in ). They considered however, a single random complex ensemble with non-Gaussian measure and they imposed the rotational symmetry of the spectrum. In our case, we considered the infinite product representing the Ginibre-Girko (therefore Gaussian) multiplicative diffusions, and we obtained the structural change of the spectrum for the general case, with non-trivial dependence on the polar variables. It is an intriguing problem for further studies, to what extent these results reflect the general ”single ring hypothesis”.

We pointed out, in agreement with , that stability exponents considered in our paper carry much richer information about the properties of the system than the Lyapunov exponents.

As a by-product of our analysis, we found a previously unknown (as far as we know) critical behavior in the case of the particular product of two hermitian ensembles. The a-priori complex spectra of the product are nevertheless real for evolution times below a certain critical value $`\tau _{crit}`$, and migrate into the complex plane only for $`\tau >\tau _{crit}`$. A similar in spirit localization-delocalization phase transition was recently observed for a class of non-hermitian Anderson Hamiltonians .

One of the motivations for this work was to propose a general formalism, which can provide a straightforward method of analyzing spectral properties of multivariate diffusion-like processes, with the idea that our method could be used in different branches of theoretical physics and interdisciplinary applications. Some particular applications will be presented elsewhere. Finally, we would like to mention, that the presented formalism allows also to establish a direct link to the diffusion processes based on Free Random Variables techniques , in particular it can demonstrate the emergence of the complex Burgers equations governing the evolution of the spectral generating functions .

### Acknowledgments

This work was partially supported by the Polish State Committee for Scientific Research (KBN) grants 2P03B 09622 (2002-2004), 2P03B08225 (2003-2006) and the EU Center of Excellence in Information Society Technologies ”COPIRA”. JJ and MAN would like to thank the Niels Bohr Institute, where a part of this work has been completed, for the hospitality. Their stay at NBI was supported by “MaPhySto”, the Center of Mathematical Physics and Stochastics, financed by the National Danish Research Foundation.
The authors would like to thank Angelo Vulpiani for interesting discussions. We are also very grateful to Roland Speicher and Piotr Śniady for comments after submitting this paper to math-ph and for valuable remarks in relation to Free Random Variable calculus.


format: text

text:
# Field-induced magnetic anisotropy in La0.7Sr0.3CoO3

## Abstract

Magnetic anisotropy has been measured for the ferromagnetic La0.7Sr0.3CoO3 perovskite from an analysis of the high-field part of the magnetization vs. field curves, i.e., the magnetic saturation regime. These measurements give a magnetic anistropy one order of magnitude higher than that of reference manganites. Surprisingly, the values of the magnetic anisotropy calculated in this way do not coincide with those estimated from measurements of coercive fields which are one order of magnitude smaller. It is proposed that the reason of this anomalous behaviour is a transition of the trivalent Co ions under the external magnetic field from a low-spin to an intermediate-spin state. Such a transition converts the Co3+ ions into Jahn-Teller ions having an only partially quenched orbital angular momentum, which enhances the intra-atomic spin-orbit coupling and magnetic anisotropy.

Attention to cobalt perovskites of formula R1-xAxCoO3 was first paid in the 50’s. At that time, Jonker and van Santen  and Koehler and Wollan  described their basic properties, that were first interpreted by Goodenough . The finding of colossal magnetoresistance effects in Mn perovskites  reactivated the interest on these and related materials, like the title compound. Since the first investigations it was seen that the magnetic and transport properties of cobalt perovskites have a peculiar thermal dependence . To explain it, the existence of a thermally induced change in the electronic configuration of the Co ions was established by Raccah and Goodenough . It is a consequence of the interplay between the intra-atomic exchange and the crystal electric field interactions, which in this particular compound are of comparable magnitude. This fact is highlighted as the source of many of the differences with manganites , where the intra-atomic exchange is rather higher than the crystal electric field interaction. For the particular case of LaCoO3, the electronic state changes at intermediate temperatures. At low temperatures, the predominant state is low-spin ($`t_{2g}^6`$, S=0) and, at intermediate ones, it evolves to an intermediate- ($`t_{2g}^5`$$`e_g^1`$, S=1) or high-spin state ($`t_{2g}^4`$$`e_g^2`$, S=2) , due to the small energy gap ($`0.03`$ eV) between such configurations. Sr-doping leads to a phase segregation into hole-rich ferromagnetic clusters and a hole-poor LaCoO3-like matrix . The hole-rich regions are isolated for low x, but at $`x`$ 0.20 they percolate and the samples become ferromagnetic and metallic. The phase diagram of this very complex system has been given by Señarís Rodríguez and Goodenough .

One of the specific features of La0.7Sr0.3CoO3 is the difficulty of approaching a magnetic saturation state . This has been attributed to the absence of true long-range-order, that leads to a cluster-glass state . Nevertheless, the ultimate origin of the situation is still not clear, and some hypotheses have been launched recently, like the presence of diamagnetic S=0 Co ions that would dilute the magnetic lattice . In any case, after the percolation of the hole-rich ferromagnetic clusters the system can be considered as a ferromagnet.

As a new input to this puzzle, Ibarra et al. reported recently a huge anisotropic magnetostriction  that cannot be explained on the basis of the usual spin-orbit coupling contribution. To explain it, they proposed an orbital instability of Co3+ under an external magnetic field. In this work the hypothesis of the orbital instability is adopted to analyze the high magnetic anisotropy of the La0.7Sr0.3CoO3 system.

The sample was prepared by conventional ceramic methods. X-ray powder diffraction revealed that it was single phase. Initial magnetization data were taken with vibrating-sample and SQUID magnetometers in fields up to 55 kOe. Before each run, the sample was demagnetized heating up to 350 K, well above the Curie point, and cooled down to the test temperature at zero field.

In Fig. 1 the high-field part of the magnetization vs. magnetic field isotherms is shown. As can be seen, the difficulty to achieve saturation is greater than in ferromagnetic manganites. In order to analyze the approach to saturation, magnetization we use the empirical relationship

$$M_H=M_S(T)\frac{a}{H}\frac{b}{H^2}+\chi H$$
(1)
where MH is the component of the magnetization along the field direction, MS is the saturation magnetization and a, b and $`\chi `$ are constants. The term $`\chi H`$ represents the field-induced increase in the spontaneous magnetization of the domains and it is very small at temperatures well below the Curie temperature. The $`a/H`$ term is generally interpreted as due to imperfections such as dislocations or nonmagnetic inclusions, and $`b/H^2`$ is due to crystalline anisotropy. To expedite the fits, the derivative is often used

$$(\frac{dM_H}{dH}\chi )H^3=aH+2b.$$
(2)
Data were fitted to achieve the $`\chi `$ that best linearizes the curves. In Fig. 2(a) we show the representation of the magnetization according to Eq. 2, with the $`\chi `$’s shown in Fig. 2(b). Linear fits were then done for data above 30 kOe in order to obtain $`a`$ and $`b`$. These data, substituted in Eq. 1, allow to fit the magnetization data and to determine the saturation magnetization shown in Fig. 2(c).

From the value of $`b`$, the values of the anisotropy constant can be calculated. For a polycrystalline material with crystallites oriented at random, the approach to saturation may be calculated by averaging the magnetization curves for all crystal orientations . The b constant is then given by

$$b=\beta \frac{K^2}{M_S}$$
(3)
where $`\beta `$ comes from purely geometric considerations. As this study only estimates the orders of magnitude of the magnetic anisotropy, we chose the same value of $`\beta =8/105`$ used for the reference manganese perovskite La2/3Sr1/3MnO3 (Ref. ), understanding that the choice does not affect the final results. The values thus obtained for the anisotropy constant are reported in Fig. 3.

The effects of magnetocrystalline anisotropy are detectable in other features of ferromagnetic materials, like the value of coercive fields. The coercive field can be estimated from the results of Fig. 2 with the well-known relation

$$H_C(T)=A\frac{K}{M_0(T)}$$
(4)
where A is a constant tipically between 0.1 and 1 (Ref. ) and M0 is the spontaneous magnetization. Taking into account that for T $`<`$ 0.8TC M$`{}_{0}{}^{}`$MS (MS can be used instead of M0 at low temperatures), we arrive to the fact that expected coercive fields are of the order of 103 or 104 Oe. The coercive fields of La0.7Sr0.3CoO3 have already been reported by Vázquez et al. ; surprisingly, a clear discrepancy is found with the values expected from Eq. 4, in one or two orders of magnitude, depending on the A value.

A deeper insight is needed, and the value of A is important to quantify this anomaly. In order to determine A, we have analyzed the manganite La0.7Ca0.3MnO3. From magnetization vs. field isotherms at 5 K and fitting the high-field part in the same way as previously, we obtain K= $`3.410^5`$ $`erg/cm^3`$ ($`5.8610^4`$ $`erg/g`$), which is consistent with the results of Perekalina et al. , who give values ranging from $`3.510^4`$ to $`610^4`$ erg/g for the magnetic anisotropy of La0.7(Sr,Pb)0.3MnO3. It is worth mentioning that this value is one order of magnitude smaller than that of the cobaltite. Provided that coercive fields at 5 K are of the order of $``$ 50 Oe and M$`{}_{S}{}^{}650`$ $`emu/cm^3`$, we arrive to the conclusion that A$``$ 0.1, which will be used for the analysis of La0.7Sr0.3CoO3. With this value of A and the data of HC and MS, we can recalculate the anisotropy constant by Eq. 4. The result is shown in Fig. 3. A simple comparison of these data (measured from coercive fields, i.e., at low applied magnetic fields) with Fig. 3 (obtained from the high-field regime) show that they are around one order of magnitude smaller.

The question now is to identify what is happening with the magnetic anisotropy of this compound. For the measurements of the two sets of anisotropy constant values (expected in principle to be similar), the only external parameter that has been changed is the external magnetic field. In one case it is high (saturation regime) and in the other it is low (when measuring coercive fields, that are of the order of $`10^2`$ Oe). Then, it seems that the magnetic anisotropy constant of La0.7Sr0.3CoO3 does not remain the same under external magnetic fields. We see a strong parallelism of this situation with the finding by Ibarra et al.  of a huge anisotropic magnetostriction induced by the external magnetic field. They concluded that the anisotropic magnetostriction cannot be explained on the basis of the usual spin-orbit contribution and proposed that it is due to the orbital instability of the Co ions under the external applied field; i.e., on the transformation of the low-spin state of the Co3+ ions in the insulating matrix (which is diamagnetic and does not contribute to magnetic anisotropy) into the energetically close intermediate-spin configuration, which is a Jahn-Teller ion. As described by Ibarra et al. , under Jahn-Teller distortion, the doubly degenerated $`e_g`$ level of IS-Co3+ splits into two singlets (L=0), and the triplet $`t_{2g}`$ into a singlet and a doublet (L=1). The singlet states of the $`t_{2g}`$ level are occupied by two electrons with opposite spins and the doublet level by three electrons. Due to the degeneracy of the doublet level with nonzero angular momentum, a strong intra-atomic spin-orbit coupling is created which leads to the increase in magnetic anistropy, in contrast to the reference manganese perovskites. Our result can also be interpreted as another evidence in favour of the orbital instability of the Co ions under external magnetic fields in La0.7Sr0.3CoO3.

###### Acknowledgements.

Spanish MCYT, under project FEDER MAT2001-3749, is acknowledged. We also wish to thank L. E. Hueso for data on La0.7Ca0.3MnO3.


text:
# Momentum Imparted by Gravitational Waves

## 1 Introduction

The localization of energy and momentum has been one of the most interesting and important problems in Einstein’s theory of General Relativity (GR). The different attempts at constructing an energy-momentum density do not provide a generally accepted definition. Consequently, different people have different points of view. In one of his papers Cooperstock  argued that in GR, energy and momentum are localized in regions of the non-vanishing energy and momentum tensor. As a result gravitational waves are not carriers of energy and momentum in vacuum. Since the gravitational waves, by definition, have zero stress-energy tensor. Thus the natural question about the existence of these waves was raised. However, the theory of GR indicates the existence of gravitational waves as solutions of Einstein’s field equations . This contrary result arises because energy is not well defined in GR.

Ehlers and Kundt , Pirani  and Weber and Wheeler  resolved this problem for gravitational waves by considering a sphere of test particles in the path of the waves. They showed that when gravitational waves passed through the test particles these acquired a constant momentum from the waves. However, the procedure  does not provide a simple prescription that can be used for arbitrary case. Qadir and Sharif  used an operational procedure, embodying the same principle and found a closed form formula which can be applied to arbitrary spacetimes. This procedure provided the similar results when we applied to plane and cylindrical gravitational waves . In a recent work , this approach has been used to find the momentum imparted by gravitational waves with spherical wavefronts. Interestingly, we obtain physically acceptable results which coincide with results obtained by using M$`\ddot{o}`$ller’s prescription. This paper focuses the problem of finding the momentum imparted by colliding gravitational waves in a closed Friedmann Robertson-Walker (FRW) background and also by gravitational waves with toroidal wavefronts.

We use an operational procedure for expressing the consequences of Relativity in terms of the Newtonian concept of gravitational force. The pseudo-Newtonian gravitational force is defined as the vector whose intrinsic derivative along the separation vector is the maximum tidal force, which is given by the acceleration vector for preferred class of observers. We will not discuss this formalism here as it is available in detail elsewhere . We shall restrict ourselves to briefly mention the essential points of the formalism. In the free fall rest-frame the extended pseudo-Newtonian four-vector force can be written in the form

$$F_0=m[\{\mathrm{ln}(A/\sqrt{g_{00}})\}_{,0}g_{ij,0}g_{,0}^{ij}/4A],F_i=m(\mathrm{ln}\sqrt{g_{00}})_{,i},$$
(1)
where $`A=(\mathrm{ln}\sqrt{g})_{,0},g=det(g_{ij}),i,j=1,2,3`$. This force formula is not uniquely fixed rather than it depends on the choice of frame. We define the quantity whose proper time derivative is $`F_a,(a=0,1,2,3)`$ as the momentum four-vector of the test particle. Its spatial components give the momentum imparted to test particles in the preferred frame (where $`g_{0i}=0`$). Thus the momentum four-vector, denoted by $`p_a`$, is given as

$$p_a=F_a𝑑t,$$
(2)
The rest of the paper is organized as follows. In the next section, we shall describe the colliding gravitational waves in a closed FRW background and toroidal waves. In section three, we use the operational procedure to evaluate momentum imparted to test particles by these waves. Finally, section four is devoted to a brief summary of the results and concluding remarks.

## 2 Gravitational Waves

### 2.1 Colliding Gravitational Waves

The general line element describing the gravitational waves can be written in the form

$$ds^2=e^M(dt^2dx^2)e^U(e^Vdy^2+e^Vdz^2),$$
(3)
where $`U,V`$ and $`M`$ are functions of $`t`$ and $`x`$.

In the presence of a stiff fluid or a scalar field, the field equations imply that $`e^U`$ satisfies the wave equation

$$(e^U)_{tt}(e^U)_{xx}=0.$$
(4)
It is well known that the fluid can be defined in terms of a potential function $`\sigma (t,x)`$ which satisfies the wave equation

$$(e^U\sigma _t)_t(e^U\sigma _x)_x=0.$$
(5)
In terms of this potential, the fluid density is given by

$$32\pi \rho =e^M(\sigma _t^2\sigma _x^2)$$
(6)
and its 4-velocity is proportional to the gradient of $`\sigma `$.

The main gravitational field equation takes the form

$$(e^UV_t)_t(e^UV_x)_x=0$$
(7)
and the remaining equations, which are now automatically integrable, become

$`U_tM_t+U_xM_x+U_{tt}+U_{xx}={\displaystyle \frac{1}{2}}(U_t^2+U_x^2+V_t^2+V_x^2+\sigma _t^2+\sigma _x^2),`$
$`U_xM_t+U_tM_x+2U_{tx}=U_tU_x+V_tV_x+\sigma _t\sigma _x.`$ (8)
In principle, these equations can always be integrated for $`M`$. Techniques for obtaining exact solutions of these equations have been given by Carmeli et al .

It may be noted that the solution describing the closed FRW stiff fluid model can be given by

$`e^U=\mathrm{sin}2t\mathrm{sin}2x=\mathrm{sin}^2(t+x)\mathrm{sin}^2(tx),V=\mathrm{ln}\mathrm{tan}x,`$
$`M=\mathrm{ln}\mathrm{sin}2t\mathrm{ln}\gamma ,\sigma =\sqrt{3}\mathrm{ln}\mathrm{tan}t,`$ (9)
where $`0<t<\pi /2`$ and $`0<x<\pi /2`$. In this spacetime, we can define two null hypersurfaces $`tx+a=0`$ and $`t+xb=0`$, where $`0<a<\pi /4`$ and $`\pi /4<b<\pi /2`$. These can be taken to be wavefronts of approaching gravitational waves which then collide when $`t=\frac{1}{2}(ba)`$ and $`x=\frac{1}{2}(a+b)`$. The FRW solution (9) is taken only as the background region into which these waves propagate.

### 2.2 Toroidal Gravitational Waves

The general line element of the form given by Eq.(3) can be written in cylindrical coordinates as

$$ds^2=e^M(dt^2d\rho ^2)e^U(e^Vd\varphi ^2+e^Vdz^2),$$
(10)
where $`U,V`$ and $`M`$ are functions of $`t`$ and $`\rho `$. As a result, the vacuum field equations imply that $`e^U`$ satisfies the wave equation

$$(e^U)_{tt}(e^U)_{\rho \rho }=0.$$
(11)
The function $`V`$ satisfies the linear equation

$$V_{tt}U_tV_tV_{\rho \rho }+U_\rho V_\rho .$$
(12)
and the remaining equations for $`M`$ are

$`U_{tt}U_{\rho \rho }={\displaystyle \frac{1}{2}}(U_t^2+U_\rho ^2+V_t^2+V_\rho ^2)U_tM_tU_\rho M_\rho ,`$
$`2U_{t\rho }=U_tU_\rho U_tM_\rho U_\rho M_t+V_tV_\rho .`$ (13)
If Eqs.(11) and (12) are satisfied then Eq.(13) can always be integrated for $`M`$.

We can now consider a gravitational wave with toroidal wavefront by putting

$$U=\mathrm{ln}t\mathrm{ln}\rho ,V=\mathrm{ln}t\mathrm{ln}\rho +\stackrel{~}{V}(t,\rho ),$$
(14)
where $`\stackrel{~}{V}`$ takes the form

$$\stackrel{~}{V}(t,\rho )=_{1/2}^{\mathrm{}}\varphi (k)(t\rho )^kH_k(\frac{t^2+\rho ^2a^2}{2t\rho })𝑑k$$
(15)
with an arbitrary function $`\varphi (k)`$. Clearly these solutions reduce to the spherical-fronted waves of  when $`a=0`$.

## 3 Momentum Imparted to Test Particles

In this section, we use the procedure outlined in the first section to calculate momentum imparted to test particles by colliding gravitational waves in the closed FRW spacetime background and also the gravitational waves with toroidal wavefronts. Using the four-vector force formula, given by Eq.(1) in Eq.(3), we have

$$F_0=m[\dot{U}+\frac{\ddot{M}+2\ddot{U}}{\dot{M}+2\dot{U}}\frac{3\dot{U}^2+\dot{V}^2}{\dot{M}+2\dot{U}}],F_1=\frac{m}{2}M^{},F_2=0=F_3.$$
(16)
The corresponding four-vector momentum takes the form

$$p_0=m[U+\mathrm{ln}(\dot{M}+2\dot{U})\frac{3\dot{U}^2+\dot{V}^2}{\dot{M}+2\dot{U}}𝑑t]+f_1(x),$$
(17)
$$p_1=\frac{m}{2}M^{}𝑑t+f_2(x),p_2=constant=p_3,$$
(18)
where $`f_1(x)`$ and $`f_2(x)`$ are arbitrary integration functions of coordinate $`x`$. Here dot and prime indicate derivatives with respect to time $`t`$ and coordinate $`x`$ respectively. It is to be noticed that Eqs.(17) and (18) give the general result of momentum four-vector for colliding and toroidal gravitational waves. Since we are interested to find out the momentum imparted by these gravitational waves, we would use Eq.(18) to calculate the term $`p_1`$. As, it is obvious from this expression, we need to have the value of $`M`$.

### 3.1 Momentum Imparted by Colliding Gravitational Waves

The solution describing the closed FRW stiff fluid model is given by Eq.(9). Substituting these values in Eqs.(17) and (18), we obtain

$$p_0=m\mathrm{ln}(6\mathrm{cot}2t\mathrm{sin}2x)+f_1(x),p_i=constant.$$
(19)
We see that the momentum imparted to test particles becomes constant which can be made zero for a special choice of constant to be zero.

### 3.2 Momentum Imparted by Toroidal Gravitational Waves

The solution in the wave region ($`t\rho a)`$ can be found by solving Eqs.(11) and (12) and is given in the form

$$U=\mathrm{ln}t\mathrm{ln}\rho ,V=\mathrm{ln}t\mathrm{ln}\rho +\stackrel{~}{V}(t,\rho ),$$
(20)
Using Eq.(20), we can write the equations for $`M`$, for $`t\rho a`$, in the form

$`M_t+M_\rho =({\displaystyle \frac{t\rho }{t+\rho }})(\stackrel{~}{V_t}+\stackrel{~}{V_\rho }){\displaystyle \frac{t\rho }{2(t+\rho )}}(\stackrel{~}{V_t}+\stackrel{~}{V_\rho })^2,`$
$`M_tM_\rho =({\displaystyle \frac{t+\rho }{t\rho }})(\stackrel{~}{V_t}\stackrel{~}{V_\rho })+{\displaystyle \frac{t\rho }{2(t\rho )}}(\stackrel{~}{V_t}\stackrel{~}{V_\rho })^2.`$ (21)
We now consider the case of a single component for which $`\stackrel{~}{V}(t,\rho )`$ becomes

$$\stackrel{~}{V}(t,\rho )=a_k(t\rho )^kH_k(\frac{t^2+\rho ^2a^2}{2t\rho }),$$
(22)
where $`k`$ is an arbitrary real parameter such that $`k\frac{1}{2}`$, $`a_k`$ is constant and $`t^2+\rho ^2a^2=2t\rho `$. The parameter $`k`$ can also be complex, but in this case one should consider for $`\stackrel{~}{V}`$ the real or imaginary part of the right hand side of Eq.(22). By substituting Eq.(22) into the transformed form of Eq.(12), it is found that the functions $`H_k(\frac{t^2+\rho ^2a^2}{2t\rho })`$ must satisfy the linear ordinary differential equation which can be reduced to a hypergeometric equation. The condition that $`\stackrel{~}{V}=0`$ can be expressed by the constraint $`H_k(1)=0`$. It can be seen that the outgoing toroidal wave includes an impulsive component if the lowest term in the expansion for $`V`$ has $`k=\frac{1}{2}`$. The front of gravitational wave has a step (or shock) if the lowest term has $`k=\frac{3}{2}`$. The value of $`M`$ for the wave region $`t\rho a`$ turns out to be

$$M=\frac{1}{2k}a_k(t^2\rho ^2)(t\rho )^{k1}H_{k1}\frac{1}{2k}(t\rho )^{2k}a_k^2[k^2H_k^2\frac{(t^2\rho ^2)^2}{4t^2\rho ^2}H_{k1}^2].$$
(23)
It is mentioned here that the dimension of $`a_k`$ is $`L^{2k}`$. To simplify the problem, we choose a particular value of $`k`$, i.e., $`k=1`$ for which the above equation reduces to

$$M=\frac{1}{2}a_1(t^2\rho ^2)H_0\frac{1}{2}(t\rho )^2a_1^2[H_1^2\frac{(t^2\rho ^2)^2}{4t^2\rho ^2}H_0^2],$$
(24)
where

$`H_0=\mathrm{ln}({\displaystyle \frac{t}{\rho a}}),H_1={\displaystyle \frac{1}{2}}[({\displaystyle \frac{t}{\rho a}}+{\displaystyle \frac{\rho a}{t}})ln({\displaystyle \frac{t}{\rho a}})({\displaystyle \frac{t}{\rho a}}{\displaystyle \frac{\rho a}{t}})].`$
Differentiating this value of $`M`$ with respect to $`\rho `$ and substituting in Eq.(18), then after some algebra, we arrive at the following

$`p_1=m{\displaystyle \frac{a_1}{4}}[{\displaystyle \frac{t^3}{3(\rho a)}}({\displaystyle \frac{\rho }{\rho a}}+2)t\rho +2t\rho \mathrm{ln}({\displaystyle \frac{t}{\rho a}})]+m{\displaystyle \frac{a_1^2}{4}}[{\displaystyle \frac{1}{50(\rho a)}}({\displaystyle \frac{13\rho }{\rho a}}`$
$`{\displaystyle \frac{7\rho ^2}{3(\rho a)^2}}1)t^5+{\displaystyle \frac{2}{9}}\rho ({\displaystyle \frac{\rho }{\rho a}}{\displaystyle \frac{5}{6}})t^3+{\displaystyle \frac{1}{2}}\rho (\rho (\rho 1)(\rho a)^2{\displaystyle \frac{\rho ^3}{\rho a}}`$
$`4\rho ^2)t\{{\displaystyle \frac{1}{5(\rho a)}}({\displaystyle \frac{11\rho }{5(\rho a)}}{\displaystyle \frac{7\rho ^2}{10(\rho a)^2}}{\displaystyle \frac{1}{2}})t^5+{\displaystyle \frac{2}{3}}\rho ({\displaystyle \frac{\rho }{\rho a}}+{\displaystyle \frac{2}{3}})t^3`$
$`+{\displaystyle \frac{1}{2}}\rho ^2(\rho a{\displaystyle \frac{\rho ^2}{\rho a}}4\rho )t\}ln({\displaystyle \frac{t}{\rho a}})\{{\displaystyle \frac{\rho }{10(\rho a)^2}}({\displaystyle \frac{\rho }{\rho a}}1)t^5{\displaystyle \frac{2}{3}}t^3\rho `$
$`{\displaystyle \frac{1}{2}}\rho (\rho (\rho a)+(\rho a)^22\rho ^3)t\}(ln({\displaystyle \frac{t}{\rho a}}))^2+f_2(\rho ).`$ (25)
This is the momentum imparted to test particles by gravitational waves with toroidal wavefronts. We see that the momentum term approaches to zero in the limit $`t0`$ for the particular choice of $`f_2=0`$. This result reduces to the momentum imparted by gravitational waves with spherical wavefronts for $`a=0`$ . The interpretation of the term $`p_0`$ of the four-vector momentum has been discussed in detail elsewhere .

The character of the gravitational wave near the wavefront can be determined by the component of $`\stackrel{~}{V}`$ with the minimum value of $`k`$. It can be shown that, near the wavefront as $`t\rho a`$,

$$\stackrel{~}{V}a_k(t+\rho )^1(t\rho +a)^{1+2k},M((1+2k)/2k)a_k(t\rho +a)^{2k}.$$
(26)
Making use of this value of $`M`$ in Eq.(18), we obtain

$$p_1=\frac{m}{4k}(1+2k)a_k(t\rho +a)^{2k}+f_2(\rho ).$$
(27)
This is the momentum imparted to test particles near the toroidal wavefront.

## 4 Discussion

We have evaluated momentum imparted to test particles by colliding gravitational waves in a closed FRW background and also for toroidal gravitational waves by using an entirely different approach. It can be seen that the momentum becomes constant for colliding gravitational waves in the closed FRW background. This is what one would expect from this procedure. It also follows from Eq.(25) that the momentum imparted by toroidal gravitational waves is physically well behaved. This turns out to be zero in the limit $`t0`$ for a particular value of arbitrary function. We also see that this becomes singular at $`\rho =a`$. If we take $`a=0`$, it reduces to the result of spherical wavefronts  and the singularity shifts to $`\rho =0`$ which also acts as a source of the gravitational waves inside the wave region. It is interesting to note that all these results exactly coincide with the results obtained by using M$`\ddot{o}`$ller prescription  for a particular choice of arbitrary function $`f_2`$.

Acknowledgment

The author would like to thank Ministry of Science and Technology for postdoctoral fellowship at University of Aberdeen, Scotland, UK.

References

 Cooperstock, F.I.: Found. Phys. 22(1992)1011; in Topics in Quantum Gravity and Beyond: Papers in Honor of L. Witten eds. Mansouri, F. and Scanio, J.J. (World Scientific, Singapore, 1993) 201; in Relativistic Astrophysics and Cosmology, eds. Buitrago et al (World Scientific, Singapore, 1997)61; Annals Phys. 282(2000)115.

 Misner, C.W., Thorne, K.S. and Wheeler, J.A.: Gravitation (W.H. Freeman, San Francisco, 1973).

 Ehlers, J. and Kundt, W.: Gravitation: An Introduction to Current Research, ed. L. Witten (Wiley, New York, 1962)49.

 Pirani, F.A.E.: Gravitation: An Introduction to Current Research, ed. L. Witten (Wiley, New York, 1962)199.

 Weber, J. and Wheeler, J.A.: Rev. Mod. Phy. 29(1957)509;

 Weber, J.: General Relativity and Gravitational Waves, (Interscience, New York, 1961).

 Qadir, A. and Sharif, M.: Physics Letters A 167(1992)331.

 Sharif, M.: Nuovo Cimento B 116(2001)1311.

 Qadir Asghar and Sharif, M.: Nuovo Cimento B 107(1992)1071;
Sharif, M.: Ph.D. Thesis, Quaid-i-Azam University Islamabad (1991).

 Feinstein, A. and Griffiths, J.B.: Class. Quantum Grav. 11(1994)L109.

 Carmeli, M., Charach, Ch. and Feinstein, A.: Ann. Phys. 150(1983)392.

 Alekseev, G.A. and Griffiths, J.B.: Class. Quantum Grav. 13(1996)2191.

 Alekseev, G.A. and Griffiths, J.B.: Class. Quantum Grav. 12(1995)L-13.

 Qadir, Asghar, Sharif, M. and Shoaib, M.: Nuovo Cimento B 115(2000)419.

 M$`\ddot{o}`$ller, C.: Ann. Phys. (NY) 4(1958)347; 12(1961)118.