Average Spectral Density of Multiparametric Gaussian Ensembles of Complex Matrices

A change in system conditions can affect the matrix elements $H_{kl}$, resulting in evolution of $\rho$ in matrix space. But as ensemble parameters are functions of system parameters, they also vary, leading to an evolution of $\rho$ in ensemble space too. Here the matrix space is a $2N^{2}$ dimensional space (consisting of $N^{2}$ real and $N^{2}$ imaginary components of the matrix elements, with each point in the space corresponds to a $N\times N$ non-Hermitian matrix. The ensemble space is a $2N^{2}$ dimensional space consisting of independent ensemble parameters $x_{kl},y_{kl}$ with each point representing the state of the ensemble.

For the evolving ensemble to continue representing the system, it is desirable that the evolution of $\rho$ in matrix space is exactly mimicked by that in ensemble space. This motivates us to consider a first order variation of the ensemble parameters $y_{kl;s}\rightarrow y_{kl;s}+\delta y_{kl;s}$ and $x_{kl;s}\rightarrow x_{kl}+\delta x_{kl;s}$ over time and seek (i) what type of matrix space dynamics they correspond to, and, (ii) can it lead to an equilibrium ensemble with stationary correlations? More explicitly, we consider following combination of the first order parametric derivatives

$\displaystyle T\;\rho\equiv\sum_{k,l}^{N}\sum_{s=1}^{\beta}\left[A_{kl;s}{\partial\rho\over\partial y_{kl;s}}+B_{kl;s}{\partial\rho\over\partial x_{kl;s}}\right]$

(2)

with

	$\displaystyle A_{kl;s}$	$\displaystyle=$	$\displaystyle 2y_{kl;s}\left(\gamma-2y_{kl;s}\right)-x^{2}_{kl;s}$
	$\displaystyle B_{kl;s}$	$\displaystyle=$	$\displaystyle 2x_{kl;s}\,\left(\gamma-4\,y_{kl;s}\right)$		(3)

Here $\gamma$ is an arbitrary parameter that gives the freedom to choose the end point of the evolution. For brevity of notations, hereafter $\sum_{k,l}^{N}\sum_{s=1}^{\beta}$ will be written just as $\sum_{k,l;s}$.

With help of its Gaussian form, the multi-parametric evolution of $\rho$ can be described exactly in terms of a diffusion, with finite drift, in $H$-matrix space,

$\displaystyle T\rho+C_{1}\,\rho=\sum_{k,l;s}{\partial\over\partial H_{kl;s}}\left[{\partial\over\partial H_{kl;s}}+\gamma\;H_{kl;s}\right]\rho$

(4)

with $C_{1}=\sum_{k,l;s}\left(\gamma-2\;y_{kl;s}\right)$.

As eq.(4) is difficult to solve for generic parametric values, we seek a transformation from the set of $M$ parameters $\{y_{kl;s},x_{kl;s}\}$ to another set $\{t_{1},\ldots,t_{M}\}$, with $M=2N^{2}$, such that only $t_{1}$ varies under the evolution governed by the operator $T$ and rest of them i.e $t_{2},\ldots,t_{M}$ remain constant: $T\rho+C_{1}\,\rho\equiv\frac{\partial\rho}{\partial t_{1}}$ or $T\rho_{1}\equiv{\partial\rho_{1}\over\partial t_{1}}$ where $\rho_{1}=C_{2}\;\rho$ with $C_{2}$ defined by the condition $TC_{2}-C_{1}C_{2}={\partial C_{2}\over\partial Y}=0$ (discussed in appendix A). This implies

$\displaystyle{\partial\rho_{1}\over\partial t_{1}}$

$\displaystyle=$

$\displaystyle\sum_{k,l;s}{\partial\over\partial H_{kl;s}}\left[{\partial\over\partial H_{kl;s}}+\gamma\;H_{kl;s}\right]\rho_{1}$

(5)

Clearly $t_{k}$ should satisfy the condition that

$\displaystyle\sum_{k,l;s}\left[A_{kl;s}{\partial t_{k}\over\partial y_{lk;s}}+B_{kl;s}{\partial t_{k}\over\partial x_{kl;s}}\right]=\delta_{k1}$

(6)

As discussed in appendix A, the above set of $M$ equations can be solved by standard method of characteristics. The solution for $k=1$ is

$\displaystyle t_{1}$

$\displaystyle=$

$\displaystyle\pm{2\over M}\sum_{k,l;s}\int\frac{{\rm d}y_{kl;s}}{y_{kl;s}\,w}+c_{0}$

(7)

with $c_{0}$ as the constant of integration, $w={2\gamma-4y_{kl;s}-x_{kl;s}^{2}}$ and $x_{kl;s}$ expressed as function of $y_{kl;s}$. Here the sign for $t_{1}$ is chosen so as to keep it increasing above the initial value (following standard diffusion equation approach where $t_{1}$ plays the role of time and is always positive).

We note that the form of $t_{1}$ in eq.(7) is obtained by assuming all $y_{kl;s}\not=0$. For $y_{kl;s}=0$, eq.(3) gives $A_{kl;s}=0$ and the corresponding derivative do not contribute to eq.(2). Indeed, for the cases in which a specific parameter $y_{kl;s}=0$ throughout the evolution (with $k,l,s$ arbitrary), the number $M$ of evolving parameters is less than $2N^{2}$.

For the case $x_{kl;s}\not=0$, the constants $c_{kl;s}$ and $\tilde{c}_{kl;s}$ are given by the relations $x_{kl;s}=\tilde{c}_{kl;s}(y_{kl;s}-y_{lk;s})(\gamma-2y_{kl;s}-2y_{lk;s})$ and $y_{lk;s}=c_{kl;s}\;y_{kl;s}$. Integration in eq.(7) further gives

$\displaystyle t_{1}=\pm{2\over\gamma M}\sum_{k,l;s}\int\frac{{\rm d}y_{kl;s}}{y_{kl;s}\,(2\gamma-a_{1}y_{kl;s}+a_{2}y_{kl;s}^{2}-a_{3}y_{kl;s}^{3})}+c_{0}.$

(8)

with $a_{1}=4+\gamma^{2}\tilde{c}_{kl;s}^{2}(1-c_{kl;s})^{2}$, $a_{2}=4\gamma\tilde{c}_{kl;s}^{2}(1+c_{kl;s})(1-c_{kl;s})^{2}$, $a_{3}=4\tilde{c}_{kl;s}^{2}(1-c_{kl;s}^{2})^{2}$. Here the symbol $\sum_{k,l;s}^{\prime}$ indicates the summation over non-zero $y_{kl;s}$ only.

Similarly solution of eq.(6) for $k>1$ leads to $M-1$ constants of integrations which can be determined by the initial conditions as well as global constraints on the dynamics. (As a similar approach for Hermitian operators has been discussed in a series of studies psall , and in psnh for non-Hermitian operators, the details of the derivation are not included in this work). To distinguish it from $t_{k}$ with $k>1$, hereafter $t_{1}$ will be referred as $Y$.

As a simple example (also used in our numerics discussed in section V), here we consider the case with $x_{kl;s}=0$. The latter corresponds to $\tilde{c}_{kl;s}=0$ and thereby $w=2(\gamma-2y_{kl;s})$. Following from eq.(7), we then have

$\displaystyle Y\equiv t_{1}=\pm{2\over M}\sum_{k,l,s}^{\prime}{\rm ln}\;{y_{kl;s}\over|\gamma-2y_{kl;s}|}+c_{0},$

(9)

and $t_{k}=c_{k}$ for $k>1$ with $c_{k}$ as arbitrary constants of integration. (This can also be seen directly from eq.(6): for $x_{kl;s}=0$, it reduces to the condition $\sum_{k,l,s}^{\prime}y_{kl;s}\;{\partial t_{k}\over\partial y_{kl;s}}=\delta_{k1}$ which on solving gives eq.(9) for $Y$).

An important insight in eq.(5) can be gained by noting the similarity of its form with the standard Fokker Planck (FP) equation, describing the diffusive dynamics of $2N^{2}$ ”particles” $H_{kl;s}$ in ”time” $t_{1}\equiv Y$ in a harmonic external potential and subjected to white noise. As the standard FP formulation is based on consideration of positive time like variable, hereafter we will consider $Y\geq 0$; the latter is ensured by choosing the appropriate sign in eq.(8) or eq.(9). Similar to FP equation, $\rho_{1}$ approaches the steady state in the limit ${\partial\rho\over\partial Y}\to 0$ of eq.(5) or as $Y\to\infty$. For example, a choice of almost all $y_{kl;s}\rightarrow N/(1-\tau^{2})$ and $x_{kl;s}\rightarrow(-1)^{s-1}\tau N/(1-\tau^{2})$ with $\gamma=1$ fulfills the condition for $\tau\rightarrow 0,\;\pm 1$ and can therefore lead to three different types of steady states, namely, Ginibre ($\tau=0$), GUE ($\tau=1$) and GASE ($\tau=-1$). As expected, the solution of eq.(5) with its left side set to zero is consistent with already known distributions for Ginibre, GUE and GASE (also follows from eq.(1) on substitution of the above values of $x_{kl;s},y_{kl;s}$). As mentioned previously, for a random matrix ensemble to represent a complex system, an appropriate choice of the ensemble parameters must take into account all system-specifics, thereby rendering them to be functions of the system parameters. This in turn makes $Y$ as the functional of the latter and thereby a measure of the complexity of the system.

As mentioned below eq.(3), $\gamma$ is an arbitrary parameter and can be set to unit value too. Nonetheless it is useful to retain $\gamma$ in the formualtion for two reasons: (i) it gives the freedom to analyze the statistical behavior of the system in $\gamma=0$ or $\infty$ limit (e.g. the dynamics in absence off the confining harmonic potential part or in reducing the role of diffusion relative to drift term), (ii) $\gamma$ also plays the role of a spectral unit in numerical analysis (e.g. while analyzing the approach of different ensembles to Ginibre universality class).

We emphasize that the form of eq.(5) is based on the specific forms for $A_{kl;s}$ and $B_{kl;s}$ given in eq.(3). Indeed it is motivated by our search for a combination of the parametric derivatives that would lead to a Dyson Brownian dynamics model in non-Hermitian matrix space (introduced by Dyson for Hermitian case and discussed in detail e.g. in haak ; meta ), starting from an arbitrary initial condition and with a Ginibre ensemble as the equilibrium limit. The search is motivated in turn by a hope to connect and utilise already existing technical tools developed for Dyson Brownian ensembles of Hermitian matrices haak ; meta ; apbe ; psrev ; circ ; psall ; psall1 ; psco .

The above mentioned connection motivates a natural query as to what type of dynamics $\rho$ (eq.(1)) undergoes if subjected to an evolution governed by an arbitrary combinations of parametric derivatives? Alternatively, if we consider a Brownian dynamics model for a non-Gaussian ensemble (i.e the combination of matrix element derivatives in right side of eq.(5)), can it still be described by a combination of the first order or higher order parametric derivatives? Some insights in this context are discussed in previous studies on Hermitian matrix ensembles psalt ; psco . As discussed in section V of psalt , not all parametric combinations would lead to an evolution in the Hermitian matrix space with an equilibrium limit. Besides, as discussed in section II.C or II. D of psco , not all combinations seem to be easily reducible to a single parametric dependence of the evolution (although a detailed investigation is needed in this context).

Contrary to a spectral degeneracy of a Hermitian matrix, a nonhermitian matrix is non-diagonalisable at the exceptional points i.e the points in the parameter space at which some of the eigenvalues become degenerate and the corresponding eigenvectors coalesce into one ir2 . But as the number of exceptional points is at most finite, the non-Hermitian matrices are considered to be diagonalisable. We hereafter refer the parametric conditions avoiding exceptional points as non-exceptional.

Assuming non-exceptional parametric conditions, a non-Hermitian matrix can be diagonalized by a transformation of the type $Z=UHV$ with $Z$ as the matrix of eigenvalues $z_{j}$ and $U$ and $V$ as the left and right eigenvector matrices respectively. Here we consider the case of complex matrices ($\beta=2$) only, with its eigenvalues $z_{j}\equiv\sum_{s=1}^{2}(i)^{s-1}z_{js}$, in general, distributed over an area in the complex plane and $U,V$ as unitary matrices; here subscript $s$ refers to the real or imaginary components of the eigenvalues). Let $\tilde{P}(z,y,x)$ be the jpdf of the eigenvalues at arbitrary values of $y_{kl;s}$ and $x_{kl;s}$ with $z\equiv\{z_{i}\}$. An integration of eq.(5) over all eigenfunctions components leads to an evolution of $\tilde{P}(z,y,x)$ with $t_{1}$; for symbols consistency with earlier works, hereafter $t_{1}$ will be referred as $Y$. Let $P_{1}(z_{1},\ldots,z_{N})$ be related to the normalized distribution by $\tilde{P}=C_{2}P_{1}/C$, its $Y$ governed evolution can then be described as (details discussed in appendix B)

$\displaystyle{\partial P_{1}\over\partial Y}$

$\displaystyle=$

$\displaystyle\sum_{s=1}^{2}\sum_{n=1}^{N}{\partial\over\partial z_{ns}}\left[{\partial\over\partial z_{ns}}-2\,{\partial{\rm ln}|\Delta_{N}(z)|\over\partial z_{ns}}+\gamma z_{ns}\right]P_{1}$

(10)

with $\beta=2$ and

$\displaystyle\Delta_{N}(z)\equiv\prod_{k<l}^{N}(z_{k}-z_{l})$

(11)

As discussed in appendix B, the derivation of eq.(10) is based on the set of eq.(109) describing the response of the eigenvalues as well as the eigenvectors to a small change in the matrix element $H_{kl}$; the assumption of $U,V$ as unitary matrices is necessary for these relations. As discussed in psnh , an assumption of $UV=I$ only leads to different dynamics of eigenvalues and eigenfunctions with varying matrix elements; this in turn leads to a n evolution equation for $P_{1}$ different from eq.(10).

In the limit ${\partial P_{1}\over\partial Y}\rightarrow 0$ or $Y\rightarrow\infty$, the evolution described by eq.(10) reaches the equilibrium steady state of the Ginibre ensemble with $P_{1}(z;\infty)\,\equiv|Q_{N}|^{2}=\prod_{j<k}|\Delta_{N}(z)|^{2}{\rm e}^{-{\gamma\over 2}\sum_{k}|z_{k}|^{2}}$; the steady state limit occurs for the system conditions resulting in almost all $y_{kl;s}\rightarrow N$ and $x_{kl;s}\rightarrow 0$ with $\gamma=1$. Thus eq.(10) describes a cross-over from a given initial ensemble (with $Y=Y_{0}$) to Ginibre ensemble as the equilibrium limit with $Y-Y_{0}$ as the crossover parameter. The non-equilibrium states of the crossover, given by non-zero finite values of $Y-Y_{0}$, are various ensembles of the complex matrices corresponding to varying values of $y_{kl}$’s and $x_{kl}$’s thus modelling different complex systems. A point worth emphasizing here is following: while the ”equilibrium” and ”non-equilibrium” states mentioned above occur on the ensemble space, such ensembles may well-represent complex systems in e.g. thermal non-equilibrium stages too.

Experimental results are often based on a single spectrum analysis; this leaves the unfolding by $\rho_{sm}(z)$ as the only option. But theoretical modelling of a single complex system by an ensemble of its replicas is often based on the assumptions of local ergodicity i.e the averages of its properties over a given spectral range, say $\Delta^{2}z\equiv\Delta z\Delta z^{*}$ around $z$, for a single matrix (referred as the spectral average) are same as the averages over the ensemble at a fixed $z$ (referred as the ensemble average). Mathematically this implies

$\displaystyle\langle\rho_{sm}(z)^{2}\rangle-\langle\rho_{sm}(z)\rangle^{2}\rightarrow 0,\qquad\qquad\langle\rho_{sm}(z)\rangle=\langle\rho(z)\rangle$

(14)

with $\langle.\rangle$ implying an ensemble average for a fixed $z$ and thereby permits it to be used for a theoretical unfolding of the spectrum. Prior to application of theoretical results to experimental data, it is therefore necessary to check the ergodicity of the spectrum. (An important point worth emphasizing here is as follows. The concept of ergodicity, a standard assumption in conventional statistical mechanics, implies analogy of the phase space averaging (i.e averaging over an ensemble of micro states) of a measure with time averaging. In case of a random matrix ensemble, the phase space averages are replaced by the averages over an ensemble of matrices and time averages by the spectral averages; the ergodicity of a measure now implies the analogy of ensemble averaging with spectral averaging. )

The fluctuations in local spectral density can in principle be measured in terms of the $n^{\rm th}$ order spectral density correlation $R_{n}(z_{1},\ldots,z_{n};Y)$, defined as psnh

$\displaystyle R_{n}={N!\over{(N-n)!}}\int P(z_{1},\ldots,z_{N};Y)\,{\rm D}\Omega_{n+1}$

(15)

with ${\rm D}\Omega_{n+1}\equiv\prod_{k=n+1}^{N}\;{\rm d}z_{k}\,{\rm d}z^{*}_{k}$. The above definition is valid in general for an arbitrary non-Hermitian ensemble. As clear from a comparison with eq.(13), the case $n=1$ corresponds to the ensemble averaged spectral density $R_{1}=N\,\langle\rho\rangle$ and is subjected to following normalization condition: $\int{\rm d}z\;{\rm d}z^{*}\,R_{1}(z)=N$.

An ergodic level density does not by itself implies the ergodicity of the local density fluctuations. Here we generalise the definition for the non-Hermitian case: the ergodicity of a local property say ”$f$” can be defined as follows: let $P(z,s)$ be the distribution function of $f$ in the neighbourhood ”s” of an arbitrary region ”$z$” over the ensemble and ${\tilde{P}}_{k}(z,s)$ be the distribution function of $f$ over a single realisation, say ”$k^{th}$” of the ensemble: $f$ is defined as ergodic if it satisfies the following conditions: (i) the distribution for almost each realization is same i.e ${\tilde{P}}_{k}(z,s)$ is independent of $k$ for almost all $k$, implying ${\tilde{P}}_{k}(z,s)={\tilde{P}}(z,s)$, (ii) the distribution over the ensemble is analogous to that of single realization $P(z,s)={\tilde{P}}(z,s)$, (iii) $P(z,s)$ is stationary on the spectral plane i.e independent of $z$: $P(z,s)=P(z^{\prime},s)$ for arbitrary points $z,z^{\prime}$ on the complex plane, However if condition (iii) is not fulfilled, $f$ is then referred as locally ergodic. More clearly, an ensemble is locally ergodic in a spectral range say $\Delta z$ around $z$, if the averages of its local properties over the range for a single matrix (referred as the spectral average) are same as the averages over the ensemble at a fixed $z$ (referred as the ensemble average). It is worth re-emphasizing here the relevance of the translational/ rotational invariance in the bulk of a spectrum: as clear from the above, it implies an invariance of the local fluctuation measures from one spectral point on the complex plane to another. A non-invariant spectrum can at best be locally ergodic and therefore a complete information about the spectrum requires its analysis at various spectral ranges.

A generic random matrix ensemble need not contain an explicit time dependence and it is necessary to first define ”stationarity” in this context. This refers to a dynamics of the ensemble, in ensemble space, due to variation of the ensemble parameters; this in turn leads to leads to variation of the statistical properties. The stationarity (equilibrium) limit here corresponds to the ensemble parameter values for which ensemble becomes invariant under a set of transformations related to global constraints. The intermediate states of the ensemble as its parameters vary from an arbitrary initial state to stationary limit, are then referred as ”non-stationary” or “non-equilibrium” ensembles.

As clear from the above, a non-stationarity of the spectrum in present context refers to its pseudo time-dependence. For example, as mentioned above, the solution of eq.(10) corresponds to stationary Ginibre ensemble in the equilibrium limit $Y\to\infty$. The intermediate states of the evolution for finite $Y$ are then referred as non-stationary ensembles. It is important to reemphasize the difference between non-ergodicity and non stationarity of the spectrum: the former refers to a non homogeneous statistics on the complex spectral plane at a fixed $Y$, the latter refers to change of statistics at a specific spectral point of interest with changing $Y$.

As discussed in previous section, the ensemble density (1) is both non-ergodic as well as non-stationary. In the next section, we analyze the dynamics of its spectral density in complex plane as $Y$ varies.

By expanding eq.(20) and eq.(21) in powers of $r$, we have $C(r,\theta,r^{\prime},\theta^{\prime})=-\,\sum_{m=1}^{\infty}\left({r\over r^{\prime}}\right)^{m}\,\cos[m(\theta-\theta^{\prime})]$ and $D(r,\theta,r^{\prime},\theta^{\prime})=\sum_{m=1}^{\infty}\left({r\over r^{\prime}}\right)^{m}\,\sin[m(\theta-\theta^{\prime})]$. Substitution of the above in eq.(18) and eq.(19) leads to

	$\displaystyle I_{c}(r,\theta)$	$\displaystyle=$	$\displaystyle 2\,\sum_{m=1}^{\infty}r^{m-1}\,{\tilde{Q}}_{m}(\theta)$		(25)
	$\displaystyle I_{d}(r,\theta)$	$\displaystyle=$	$\displaystyle-2\,\sum_{m=1}^{\infty}r^{m-1}\,{\tilde{S}}_{m}(\theta)$		(26)

with

	$\displaystyle{\tilde{Q}}_{m}(r,\theta)$	$\displaystyle=$	$\displaystyle P\int{\rm d}\theta^{\prime}\,\cos[m(\theta-\theta^{\prime})]\,\langle{\left(1\over r^{\prime}\right)^{m}}\rangle_{\theta^{\prime}}$		(27)
	$\displaystyle{\tilde{S}}_{m}(r,\theta)$	$\displaystyle=$	$\displaystyle P\int{\rm d}\theta^{\prime}\,\sin[m(\theta-\theta^{\prime})]\,\langle{\left(1\over r^{\prime}\right)^{m}}\rangle_{\theta^{\prime}}$		(28)

with $\langle{1\over r^{\prime m}}\rangle_{\theta^{\prime}}=\int{{\rm d}r^{\prime}\,r^{\prime}\,{\mathcal{R}}_{2}\over r^{\prime m}}$. Using eq.(25) and eq.(26), we have

$\displaystyle{1\over r}\,{\partial(r\,I_{c})\over\partial r}+{1\over r}\,{\partial I_{d}\over\partial\theta}={2\over r}\,\sum_{m=1}^{\infty}\left({\partial(r^{m}\;{\tilde{Q}}_{m})\over\partial r}-r^{m-1}\,{\partial{\tilde{S}}_{m}\over\partial\theta}\right)$

(29)

As both ${\tilde{Q}}_{m}$ and ${\tilde{S}}_{m}$ are $\theta$-dependent, substitution of the above leaves eq.(23) non-separable in $r,\theta$ variables. In the regime $r\sim{1\over\sqrt{N}}$, however, the series in eq.(29) can contribute significantly only if ${\tilde{Q}}_{m}$ and ${\tilde{S}}_{m}$ vary as $N^{(m-1)/2}$ or faster.

The latter not only requires $\langle\left({1\over r^{\prime}}\right)^{m}\rangle_{\theta^{\prime}}\sim N^{(m-1)/2}$ (such that $\int{{\rm d}r^{\prime}\,{\mathcal{R}}_{2}\over r^{\prime m-1}}\sim N^{(m-1)/2}$) but also its dependence on $\theta$ (to keep non-zero contributions from the integrals in eq.(29) for all $m\geq 1$). But as $\lim_{r,r^{\prime}\to 0}{\mathcal{R}}_{2}(r,r^{\prime},\theta,\theta^{\prime})$ is not very sensitive to angular dependence, this in turn suggests a $\theta^{\prime}$-independence of $\langle{1\over r^{\prime m}}\rangle_{\theta^{\prime}}$ near origin ($r^{\prime}\leq{1\over\sqrt{N}}$) and is also confirmed by our numerical analysis of three different ensembles, discussed later in section V. As a consequence, for almost all cases (e.g except for those with phase singularities at the origin psps ), the contribution of the series in eq.(29) to eq.(23) is negligible with respect to other terms and we can write

$\displaystyle\left[{\partial^{2}\over\partial r^{2}}+{1\over r^{2}}\;{\partial^{2}\over\partial\theta^{2}}+{2\over r}\,{\partial\over\partial r}\right]\,{\mathcal{R}}_{1}+E\,{\mathcal{R}}_{1}=0$

(30)

Substituting eq.(24) in eq.(30), we now have

	$\displaystyle r^{2}\,{\partial^{2}U\over\partial r^{2}}+2\,r\;{\partial U\over\partial r}+\left(E\,r^{2}-E_{1}\right)\,U$	$\displaystyle=$	$\displaystyle 0\,$		(31)
	$\displaystyle{\partial^{2}T(\theta)\over\partial\theta^{2}}+E_{1}\,T$	$\displaystyle=$	$\displaystyle 0$		(32)

Using separation of variables, the solution of above equations can be given as

	$\displaystyle U_{\mu}(r)$	$\displaystyle=$	$\displaystyle r^{-1/2}\,\left[a_{\mu}\,J_{\mu}\left({r\sqrt{E}}\right)+b_{\mu}\,N_{\mu}\left(r\sqrt{E}\right)\right]$		(33)
	$\displaystyle T_{\nu}(\theta,Y)$	$\displaystyle=$	$\displaystyle c_{\nu}\,\cos\left[\sqrt{E_{1}}\;\theta+d_{\nu}\right]$		(34)

with $J_{\mu},N_{\mu}$ as the Bessel’s functions, $E_{1}={1\over 4}-{\mu}^{2}$, $a_{\mu},b_{\mu},c_{\nu},d_{\nu}$ as arbitrary constants, determined by the boundary conditions on ${\mathcal{R}}_{1}(r,\theta)$; the condition that latter remains finite at $r=0$ then implies $b_{\mu}=0$. The single valuedness of $T(\theta)$ as $\theta\to\theta+2\pi$ requires $d_{\nu}=2\pi k$ and $\sqrt{E_{1}}=n$ with $k,n$ as integers and thereby $\mu=\pm\sqrt{1/4-n^{2}}$. To keep $\mu$ as a real number, only allowed values of $n=0$. The later in turn gives $\mu^{2}=1/4$. Noting that $J_{-\mu}(x)=(-1)^{n}\,J_{\mu}(x)$, hereafter we use $\mu=1/2$.

As the arbitrary constant $E$ is not subjected to any additional known constraints except semi-positive definite one, a valid particular solution can occur for continuous range of $E$ from $0\to\infty$. By substitution of the solutions in eq.(33) and eq.(34) for $n=0$ (equivalently $E_{1}=0$) in eq.(22), the general solution of eq.(17) for arbitrary initial condition (but finite at $r=0$) and for $0\leq r\leq{1\over\sqrt{N}}$ can now be given as

$\displaystyle R_{1}(r,\theta;Y)=r^{-1/2}\,\int_{0}^{\infty}{\rm d}E\;a(E)\,J_{1/2}\left(r\sqrt{E}\right)\,{\rm exp}\left[-E\,(Y-Y_{0})\right]$

(35)

with constant $a(E)$ determined by the initial condition of $R_{1}$ (examples discussed in appendix E). A circular symmetric behavior of $R_{1}$ is also confirmed by our numerics discussed in section V, thus lending credence to our separability conjecture in eq.(24).

The existence of a stationary solution in limit $Y-Y_{0}\to\infty$ requires $E=0$. The solution in the limit becomes

$\displaystyle R_{1}(r,\theta;\infty)=\lim_{E\to 0}\;{a(E)\over\sqrt{r}}\,J_{1/2}\left(r\sqrt{E}\right)={1\over\sqrt{\pi}}\hskip 21.68121pt{\rm for}\;r\leq{1\over\sqrt{N}}$

(36)

Here the condition,that $\lim_{r\to 0}R_{1}(r,\theta)$ remains finite requires $a(E)\propto E^{-1/4}$; this in turn gives a constant $R_{1}(r,\theta;\infty)$ for $r\leq{1\over\sqrt{N}}$.

We note that the solution in eq.(35) is obtained for smooth behavior of the density at $r=0$. The solution for the cases in which initial spectral density is singular at $r=0$, we must also consider $b_{\mu}\not=0$.

Expanding eq.(20) and eq.(21) in powers of ${1\over r}$, we have $C(r,\theta,r^{\prime},\theta^{\prime})=\sum_{m=0}^{\infty}\left({r^{\prime}\over r}\right)^{m+1}\,\cos[m(\theta-\theta^{\prime})]$, $D(r,\theta,r^{\prime},\theta^{\prime})=\sum_{m=0}^{\infty}\left({r^{\prime}\over r}\right)^{m+1}\,\sin[m(\theta-\theta^{\prime})]$. Substitution of the above in eq.(18) and eq.(19) leads to

	$\displaystyle I_{c}(r,\theta)$	$\displaystyle=$	$\displaystyle-2\,\sum_{m=0}^{\infty}\left({1\over r}\right)^{m+1}\,Q_{m}(r,\theta)$		(37)
	$\displaystyle I_{d}(r,\theta)$	$\displaystyle=$	$\displaystyle-2\,\sum_{m=0}^{\infty}\left({1\over r}\right)^{m+1}\,S_{m}(r,\theta)$		(38)

with

	$\displaystyle Q_{m}(r,\theta)$	$\displaystyle=$	$\displaystyle{\bf P}\int{\rm d}\theta^{\prime}\,\cos[m(\theta-\theta^{\prime})]\;\langle{(r^{\prime})^{m}}\rangle$		(39)
	$\displaystyle S_{m}(r,\theta)$	$\displaystyle=$	$\displaystyle{\bf P}\int{\rm d}\theta^{\prime}\,\sin[m(\theta-\theta^{\prime})]\;\langle{(r^{\prime})^{m}}\rangle.$		(40)

with $\langle(r^{\prime})^{m}\rangle=\int{\rm d}r^{\prime}\,r^{\prime m+1}\,{\mathcal{R}}_{2}(r,\theta,r^{\prime},\theta^{\prime})$ as a function of $r,\theta,\theta^{\prime}$. Using the above, we have