text:
# N-Site approximations and CAM analysis for a stochastic sandpile

## I Introduction

Sandpile models are the prime example of self-organized criticality (SOC) , in which a system with an absorbing-state phase transition is forced to its critical point , leading to scale-invariance in the apparent absence of parameters . The absorbing-state phase transition, which depends, as is usual, on the fine-tuning of one or more control parameters, is evident in sandpiles with a fixed number of particles , models that have come to be called fixed-energy sandpiles (FES). While most studies of sandpiles have probed the driven case , there is great interest in understanding the scaling properties of FES models as well .

Previous studies of FES reveal that they exhibit a phase transition between an absorbing and an active state as the particle density $`\zeta `$ (which is the temperaturelike control parameter) is increased beyond a critical value . Until now, all quantitative results on FES have been obtained from simulations. It is therefore of interest to apply theoretical methods to such models. One such approach is Suzuki’s coherent anomaly method (CAM) for analyzing a series of cluster approximations. It has been shown to yield good estimates for critical properties both in  and out of equilibrium . In this work I develop $`n`$-site approximations for a one-dimensional sandpile model, and analyze the results using the CAM. This represents the first application of the CAM to a model representative of the class of absorbing-state phase transitions in systems with a conserved density .

In this paper I study a FES with a height restriction. From the theoretical viewpoint, an inconvenient feature of sandpile models is the unbounded number of particles that may occupy the same site; this complicates attempts to derive cluster approximations. In Manna’s stochastic sandpile , sites with height $`z2`$ are active. If we restrict the height (or number of sand grains per site), to be $`2`$, the effect on critical properties should be minimal, aside from a possible shift in the critical density $`\zeta _c`$. This expectation was recently verified numerically: the restricted-height stochastic sandpile belongs to the same universality class as its unrestricted counterpart . I study the restricted-height model for calculational convenience, as a representative of a broader universality class that includes Manna’s stochastic sandpile, the conserved lattice gas, and the conserved threshold transfer process .

The balance of the paper is organized as follows. The model is defined in Sec. II, followed by a discussion of $`n`$-site approximations in Sec. III. Numerical results are presented in Sec. IV. The CAM analysis is discussed in Sec. V, and in Sec. VI I present a brief summary.

## II Model

The model is defined on a ring of $`L`$ sites with periodic boundaries. (The cluster approximations effectively study the $`L\mathrm{}`$ limit.) The configuration is specified by the number of particles $`z_i=0,1`$, or 2 at each site; sites with $`z_i=2`$ are said to be active, and have a toppling rate of unity. The continuous-time (sequential), Markovian dynamics consists of a series of toppling events at individual sites. When site $`i`$ topples, two particles attempt move to randomly chosen nearest neighbors $`j`$ and $`j^{}`$ of $`i`$. The new position of each particle is accepted if and only if the target site has fewer than two particles. I consider a stochastic toppling rule in which the two particles move independently. Any particle attempting to move to a site already harboring two particles is sent back to the toppling site. (Thus an attempt to send two particles from site $`j`$ to site $`k`$, with $`z_k=1`$, results in $`z_k=2`$ and $`z_j=1`$.) Transition probabilities are listed in Table I.

This model, and a closely related one (with a cooperative toppling rule), were studied via simulation in Ref. , which showed that the critical exponents $`\beta `$ and $`\nu _{}`$ are the same as for the unrestricted Manna sandpile. A similar conclusion was reached in Ref.  for a two-dimensional restricted-height sandpile (called the conserved threshold transfer process in that work).

## III Cluster approximations

I have derived dynamic $`n`$-site cluster approximations for the one-dimensional restricted sandpile model. Such approximations often yield qualitatively correct phase diagrams . The procedure parallels that used by Ferreira and Mendirata to study the one-dimensional contact process . The $`n`$-site approximation consists of a set of coupled differential equations for the probabilities $`P_𝒞^{(n)}`$ of each $`n`$-site configuration, $`𝒞`$. (There are $`3^n`$ such configurations, but the number of independent probabilities is $`3^n/2`$, due to symmetries.) The system is assumed homogeneous, so that the $`P_𝒞^{(n)}`$ are independent of position.

Since transitions in a set of $`n`$ contiguous sites generally depend on sites outside the cluster, the $`n`$-site probabilities are coupled to those for $`n+1`$ and so on, generating an infinite hierarchy of equations. The $`n`$-site approximation truncates this hierarchy by approximating $`m`$-site probabilities (for $`m>n`$) in terms of $`n`$-site conditional probabilitites. In the $`n`$-site approximation, the joint probability for a sequence of $`n+1`$ sites is approximated so :

$`P^{(n+1)}(z_1,\mathrm{},z_{n+1})`$ $``$ $`P^{(n)}(z_{n+1}|z_n,\mathrm{},z_2)P^{(n)}(z_n,\mathrm{},z_1)`$ (1)
$`=`$ $`{\displaystyle \frac{P^{(n)}(z_{n+1},\mathrm{},z_2)P^{(n)}(z_n,\mathrm{},z_1)}{P^{(n1)}(z_n,\mathrm{},z_2)}}.`$ (2)
The equations for one- and two-site approximations are relatively simple to derive, and are described in Ref. . I have developed a computational algorithm capable of generating the approximation for arbitrary $`n`$. Each configuration $`𝒞=(z_n,\mathrm{},z_1)`$ is represented by an integer

$$I(𝒞)=\underset{k=1}{\overset{n}{}}z_k3^{k1}.$$
(3)
The calculation begins with the generation of all configurations, corresponding to each integer from zero (all sites empty) up to the maximum, $`3^n1`$ (all sites doubly occupied); the symmetry (under inversion) of each configuration is determined. If $`𝒞`$ is not symmetric, then it and its mirror image $`𝒞_R`$ must have the same probability, and only the smaller of $`𝒞`$ and $`𝒞_R`$ is treated explicitly, reducing the number of variables by roughly half.

Next, a list of all possible transitions is constructed. Here it is useful to distinguish between central transitions (involving a toppling at one of the sites 2,…,$`n1`$) and boundary transitions, in which either site 1 or site $`n`$, or one of the peripheral sites (0 or $`n+1`$) topples. The rate of a central transition $`𝒞𝒞^{}`$ is the product of a branching probability $`p_b`$ (for the particles to be redistributed in a particular manner, as in Table I), and the intrinsic toppling rate, which is unity. Consider, for example, configuration $`𝒞=(21120)`$. The transition rates associated with a toppling at the second site (counting from the right), are

$`w[(21120)(21102)]={\displaystyle \frac{1}{4}},`$

$`w[(21120)(21201)]={\displaystyle \frac{1}{2}},`$

$`w[(21120)(21210)]={\displaystyle \frac{1}{4}}.`$

For each configuration $`𝒞`$, the set of allowed transitions to other states, and the associated rates, are stored. In the case of a central transition from $`𝒞`$ to $`𝒞^{}`$, the contribution to the time-derivative of the probability has the form:

$`{\displaystyle \frac{dP(𝒞^{})}{dt}}=w[𝒞𝒞^{}]P(𝒞).`$There is of course a corresponding loss term for $`P(𝒞)`$:

$`{\displaystyle \frac{dP(𝒞)}{dt}}=w[𝒞𝒞^{}]P(𝒞).`$Thus central transitions contribute to the evolution of the probability distribution precisely as in the master equation.

For boundary transitions, one does not have access to the $`n+1`$-site probabilitites required to mount a complete description, and so must resort to the truncation scheme embodied in Eq. (2). For example, the contribution to $`dP(2,z_2,\mathrm{},z_n)/dt`$ due to the transition $`(2,z_2,\mathrm{},z_n)(1,z_2,\mathrm{},z_n)`$ is

$`{\displaystyle \frac{1}{2}}P^{(n+1)}(0,2,z_2,\mathrm{},z_n)+{\displaystyle \frac{3}{4}}P^{(n+1)}(1,2,z_2,\mathrm{},z_n).`$The $`P^{(n+1)}`$ are estimated using Eq. (2). For boundary transitions one stores not only the rate, but the two configurations (aside from the original one, $`𝒞`$) whose probabilities are needed to evaluate $`dP(𝒞)/dt`$. With this information available, one can evaluate the derivatives $`dP(𝒞)/dt`$ for all possible configurations, given the probability distribution.

The evolution of the probability distribution is found via numerical integration, using a fourth-order Runge-Kutta scheme . The integration is halted when a stationary distribution is attained, that is, when the time-derivatives $`dP(𝒞)/dt`$ all have an absolute value smaller than $`\delta `$ (typically, $`\delta =10^{13}`$). An interesting technical point concerns the evaluation of the $`n1`$-site marginal distribution. There are evidently two equivalent expressions that may be used:

$`P^{(n1)}(z_{n1},\mathrm{},z_1)={\displaystyle \underset{z_n=0}{\overset{2}{}}}P^{(n)}(z_n,z_{n1},\mathrm{},z_1),`$and

$`P^{(n1)}(z_{n1},\mathrm{},z_1)={\displaystyle \underset{z_0=0}{\overset{2}{}}}P^{(n)}(z_{n1},\mathrm{},z_1,z_0).`$Numerical stability is greatly enhanced using the mean of the two expressions given above.

For sizes $`n7`$, very near the critical point, relaxation to the stationary distribution is very slow, and the following procedure proves advantageous. Let $`D=\mathrm{max}_𝒞|dP(𝒞)/dt|`$ be the largest derivative (in absolute value). The properties of interest (principally, the active-site density) are recorded as a function of $`D`$, and the integration halted when $`D<10^{10}`$. Fig. 1 shows the result for the active site density, for $`\zeta `$ slightly above $`\zeta _c`$. The stationary value is obtained via extrapolation to $`D=0`$, usually via a quadratic fit to the four data points for smallest $`D`$. (The resulting correction is typically less than 1% of the value at $`D=10^{10}`$.) I also studied the order parameter relaxation rate $`\gamma =|\dot{\rho }/\rho |`$ for each $`n`$ at a series of $`\zeta `$ values near, but below, $`\zeta _c`$. These data are used to estimate the critical exponent $`\nu _{||}`$ in Sec. V.

## IV Numerical results

I derived cluster approximations for $`n11`$ sites, yielding the stationary active-site density $`\rho _a`$ as a function of $`\zeta `$. The $`n`$-site approximation predictions for $`\rho _a(\zeta )`$ ($`n=3`$ to 11) are compared against simulation in Fig. 2; the theoretical curves appear to approach the simulation result systematically.

For each $`n`$, the active-site density is zero below a certain critical value, $`\zeta _{c,n}`$. Since the phase transition in the stochastic sandpile is continuous , one expects the same to be true of the cluster approximations. This is indeed the case for $`n4`$, but for $`n=5`$, 6 and 7 there is a very small discontinuity in $`\rho _a`$ ($`\stackrel{<}{}10^3`$, invisible on the scale of Fig. 2), as we decrease $`\zeta `$. Since the same procedure is used for all $`n`$, the discontinuity is unlikely to be artefact of the numerical method. On the other hand, I do not regard the discontinuity as physically significant; it appears to represent an unphysical feature of the cluster approximations for certain $`n`$ values, very near the critical density. In the CAM analysis I disregard the behavior of $`\rho _a`$ in the immediate vicinity of $`\zeta _{c,n}`$, and instead analyze its properties at points somewhat removed from the transition.

There remains, naturally, the problem of estimating the critical density, $`\zeta _{c,n}`$. For each $`n`$ I determine the critical density by fitting the four or five data points nearest the transition, where $`\rho 10^3`$ or less. In each case, I plot $`\mathrm{ln}\rho _a`$ versus $`\mathrm{ln}(\zeta \zeta _{c,n})`$, varying $`\zeta _{c,n}`$ to obtain the best power-law fit. The associated slopes vary between 1 (for $`n=3`$) and about 0.25 (for $`n=6`$), but these, again, are regarded as unimportant details of the approximation in question. It is important to stress that, in the cases where the transition is apparently discontinuous, the difference between the location of the discontinuity and the extrapolated value of $`\zeta _{c,n}`$ is less than one part in 105, and that the estimates for $`\zeta _c`$ and critical exponents are insensitive to these tiny differences. The values of $`\zeta _{c,n}`$ obtained in this manner are listed in Table II.

Using the results for $`\zeta _{c,n}`$, I estimate $`\zeta _c=lim_n\mathrm{}\zeta _{c,n}`$ by plotting $`\mathrm{\Delta }_n\zeta _c\zeta _{c,n}`$ versus $`n`$ in a double-logarithmic plot, varying $`\zeta _c`$ to obtain the best power-law fit. The latter is obtained using $`\zeta _c`$ in the range 0.929 - 0.931, yielding $`\zeta _c=0.930(1)`$, in good agreement with the simulation result of 0.92965 . The finite-size scaling prediction for the critical point shift is : $`\mathrm{\Delta }_n^\nu _{}1/n`$. I obtain a good fit to the data (see Fig. 3) using $`\nu _{}=1.66`$ (as found in simulations ), including a correction to scaling term:

$`\mathrm{\Delta }_n^\nu _{}{\displaystyle \frac{A}{n}}+{\displaystyle \frac{B}{n^2}}.`$The numerical data are consistent with the simulation estimate for $`\nu _{}`$, but not sufficient to furnish an independent estimate of the exponent.

## V Coherent anomaly analysis

A detailed explanation of the CAM procedure is given in Ref. ; it may be understood on the basis of finite-size scaling . The approach here parallels that used by Tomé and de Oliveira in their study of the Domany-Kinzel model . To begin, one argues that the cluster size $`n`$ plays the role of an effective system size $`L`$ as regards scaling properties. This is because the $`n`$-site approximation effectively cuts off correlations of range $`>n`$ (notwithstanding the fact that cluster approximations nominally treat an infinite system). Thus, as noted above, one expects a critical point shift $`\mathrm{\Delta }_nn^{1/\nu _{}}`$. Finite-size scaling theory also yields the relation $`\rho _{a,n}(\zeta _c)n^{\beta /\nu _{}}\mathrm{\Delta }_n^\beta `$ for the order parameter in a finite system, at the (true) critical point. For $`\zeta >\zeta _{c,n}`$, $`\rho _{a,n}(\zeta )`$ is a smooth function. Thus we are led to a scaling hypothesis for the order parameter :

$$\rho _{a,n}(\zeta )=\mathrm{\Delta }_n^\beta f(\frac{\zeta \zeta _{c,n}}{\mathrm{\Delta }_n}),$$
(4)
where $`f(x)`$ is a scaling function with $`f(0)=0`$. If we suppose that $`f(x)x^{\beta _{MF}}`$ for $`0x1`$, then $`\rho _{a,n}(\zeta )=A_n(\zeta \zeta _{c,n})^{\beta _{MF}}`$, where the amplitude $`A_n`$ diverges as $`n\mathrm{}`$:

$$A_n\mathrm{\Delta }_n^{(\beta _{MF}\beta )}.$$
(5)
This is the usual CAM relation. On the other hand, the hypothesis that $`n`$ is equivalent to a finite system size leads directly to:

$$\rho _{a,n}(\zeta _c)\mathrm{\Delta }_n^\beta .$$
(6)
This expression involves the behavior of the $`n`$-site approximation at the critical point $`\zeta _c`$ not $`\zeta _{c,n}`$. It is interesting to note that the hypothesis of an effective system size directly proportional to $`n`$ is not strictly necessary. The scaling relations involving $`\mathrm{\Delta }_n`$ follow from the more general hypothesis of an effective system size $`L_{eff}=L_{eff}(n)`$, for example $`L_{eff}n^\varphi `$ with $`\varphi >0`$.

### A CAM analysis for $`\beta `$

As noted above, the $`n`$-site approximations for the order parameter $`\rho _{a,n}`$ are not all well behaved in the vicinity of $`\zeta _{c,n}`$. For this reason, analysis of $`\rho _{a,n}`$ at $`\zeta _{c,n}`$ will not yield a consistent set of well defined amplitudes $`A_n`$. But since $`\rho _{a,n}(\zeta )`$ is well behaved for $`\zeta >\zeta _{c,n}`$, we can study its scaling at some point intermediate between $`\zeta _{c,n}`$ and $`\zeta _c`$. In particular, the scaling hypothesis Eq. (4) implies that if we fix $`x=(\zeta \zeta _{c,n})/\mathrm{\Delta }_n`$, then

$$\frac{d\rho _{a,n}}{d\zeta }=f(x)\mathrm{\Delta }_n^{\beta \beta _{MF}}.$$
(7)
Our strategy is to analyze the order parameter data reasonably near the $`n`$-site critical value, but away for $`\zeta _{c,n}`$ itself, where $`\rho _{a,n}`$ is singular. A crucial point in this analysis is the postulate that the mean-field exponent $`\beta _{MF}=1`$, regardless of the behavior of $`\rho _{a,n}`$ in the immediate vicinity of $`\zeta _{c,n}`$. The motivation for this assumption is, firstly, that $`\beta _{MF}`$ is clearly unity for $`n=1`$, 2, or 3; secondly, that a critical exponent such as $`\beta `$ is determined, in mean-field theory, by symmetry properties of the order parameter, and hence should not vary with $`n`$; and thirdly, that $`\beta _{MF}=1`$ generically for phase transitions to an absorbing state . (The basis for this last assertion is that the mean-field equation for the order parameter will have the form $`d\rho /dt=A\rho B\rho ^2`$, barring some coincidence or a symmetry that renders $`A`$ and/or $`B`$ zero .)

I evaluate $`d\rho _{a,n}/d\zeta `$ (numerically, using an interval $`\mathrm{\Delta }=0.0005`$), for fixed $`x=1/4`$; the results are shown in Fig. 4. Least-squares linear fits to the data for $`n=8`$ \- 11 yield, via Eq. (7), the value $`\beta =0.408(6)`$, where the figure in parenthesis denotes the uncertainty. Using $`\zeta _c=0.929`$ instead of the best estimate, 0.930, I find $`\beta =0.421(5)`$. Thus a reasonable estimate for $`\beta `$ is 0.41(1). (A similar analysis, but evaluating the derivatives at $`x=1/2`$, yields $`\beta =0.42`$.)

The above analysis is complemented with a study of $`\rho _{a,n}(\zeta _c)`$, as suggested by Eq. (6). The graph of $`\rho _{a,n}(\zeta _c)`$ versus $`\mathrm{\Delta }_n`$ shows (on log scales, see Fig. 5), a fair amount of curvature, making determination of $`\beta `$ more difficult in this case. Linear fits to the data for $`n`$ = 7-9, 8-10, and 9-11 yield, respectively, $`\beta `$ = 0.471, 0.460, and 0.448, consistent with an approach to the value of $`0.41`$ for large $`n`$. Verification of convergence must naturally await the evaluation of approximations for larger clusters.

A further point of interest is the validity of the scaling hypothesis, Eq. (4). The data collapse shown in Fig. 6, a plot of $`\rho ^{}=\mathrm{\Delta }_n^\beta \rho _{a,n}`$ versus $`x=(\zeta \zeta _{c,n})/\mathrm{\Delta }_n`$ provides support for the hypothesis. (In Ref.  a similar collapse is demonstrated for the Domany-Kinzel model.)

### B CAM analysis for $`\nu _{||}`$

As shown in Refs. , the CAM approach is readily extended to dynamics. Let

$$\gamma _n(\zeta )=\frac{1}{\rho }\frac{d\rho }{dt}$$
(8)
be the relaxation rate in the $`n`$-site approximation, and let $`\gamma (\zeta )`$ be the true relaxation rate. Then we expect $`\gamma |\zeta \zeta _c|^{\nu _{||}}`$, while $`\gamma _n|\zeta \zeta _{c,n}|^{\nu _{||,MF}}`$ in the $`n`$-site approximation, where the mean-field exponent is $`\nu _{||,MF}=1`$ for models with an absorbing-state phase transition . A scaling hypothesis, analogous to Eq. (4), for the relaxation rate, is

$$\gamma _n(\zeta )=\mathrm{\Delta }_n^{\nu _{||}}g(\frac{\zeta \zeta _{c,n}}{\mathrm{\Delta }_n}),$$
(9)
where the scaling function $`g`$ vanishes when its argument is zero. Supposing that $`g(x)|x|^{\nu _{||,MF}}`$, we see that $`\gamma _n(\zeta )\overline{\gamma }_n|\zeta \zeta _{c,n}|`$, where the amplitude follows

$$\overline{\gamma }_n\mathrm{\Delta }_n^{\nu _{||}\nu _{||,MF}}.$$
(10)
I determine the relaxation rate numerically for $`\zeta \stackrel{<}{}\zeta _{c,n}`$, and from these data extract the amplitudes $`\overline{\gamma }_n`$. The results, shown in Fig. 7, display substantial curvature on a log-log plot, so that direct determination of the critical exponent $`\nu _{||}`$ is not feasible. Simulations  yield estimates for $`\nu _{||}`$ in the range 2.3 - 2.6. The CAM results are consistent with values in this range, if we include a correction to scaling. The solid line in Fig. 7 is given by

$$\mathrm{ln}\overline{\gamma }_n=1.5\mathrm{ln}\mathrm{\Delta }_n+A\mathrm{\Delta }_nB,$$
(11)
with fit parameters $`A=7.435`$ and $`B=2.125`$, consistent with a correction to scaling expression $`\overline{\gamma }_n\mathrm{\Delta }_n^{\nu _{||}1}(1+A\mathrm{\Delta }_n)`$, with $`\nu _{||}=2.5`$. While results for larger clusters will be needed to determine $`\nu _{||}`$ with precision, one can at least assert that the present results are consistent with the rather imprecise estimates from simulations.

## VI Discussion

I have devised a computational algorithm for generating $`n`$-site cluster approximations for a one-dimensional stochastic sandpile model with a fixed particle density. To facilitate the analysis I impose the height restriction $`z_i2`$. Analyzing the results for $`n11`$, I obtain the estimates $`\zeta _c=0.930(1)`$, $`\beta =0.41(1)`$ and $`\nu _{||}2.5`$, all in agreement with simulation . While the results are not very precise, they provide significant independent support for the simulational findings, showing that FES models with a strict activity threshold belong to a universality class distinct from that of directed percolation (DP) . Sandpile models in which sites with an above-threshold height can remain stable (so-called “sticky grains”), have recently been shown to belong to the DP class , but such is not the case for the model studied here. Application of the methods used in this work to other sandpile models should prove illuminating.

Acknowledgments

I thank Mário de Oliveira, Tânia Tomé, Attila Szolnoki, Gyorgy Szabó, Géza Ódor, and Paulo Alfredo Gonçalves Penido for helpful comments. This work was supported by CNPq, and CAPES, Brazil.

e-mail: dickman@fisica.ufmg.br

FIGURE CAPTIONS

FIG. 1. Active-site density $`\rho `$ versus $`D=\mathrm{max}_𝒞|dP(𝒞)/dt|`$, for $`n=10`$, $`\zeta =0.8860`$.

FIG. 2. Stationary active-site density $`\rho `$ versus particle density $`\zeta `$. Solid curves: $`n`$-site approximations for $`n`$ = 3 - 11; points: simulation results for a system of 5000 sites.

FIG. 3. Critical point shift $`\mathrm{\Delta }_n^\nu _{}`$ versus$`1/n`$. The solid line is a fit including a correction term as described in the text.

FIG. 4. $`\rho _n^{}d\rho _{a,n}/d\zeta `$ versus $`\mathrm{\Delta }_n`$. The derivative is evaluated at $`x=1/4`$. The slope of the solid line is -0.592, corresponding to $`\beta =0.408`$.

FIG. 5. Active-site density $`\rho _{a,n}(\zeta _c)`$ versus $`\delta _n`$. The slope of the solid line is 0.408.

FIG. 6. Scaled active-site density $`\rho ^{}=\mathrm{\Delta }_n^\beta \rho _{a,n}`$ versus scaled particle density $`x=(\zeta \zeta _{c,n})/\mathrm{\Delta }_n`$, for $`n=6`$ \- 11.

FIG. 7. Relaxation rate amplitude $`\overline{\gamma }_n`$ versus $`\mathrm{\Delta }_n`$. The solid curve is a fit including a correction to scaling term as described in the text.


format: text

text:
# 1 Introduction

## 1 Introduction

In this note we will address the issue of universality of the probability density function (p.d.f.) of the product $`B^1X`$ of real and complex random matrices.

In order to motivate our discussion, before delving into random matrix theory, let us discuss a simpler problem. Thus, consider the random variables $`x`$ and $`y`$ drawn from the normal distribution

$$G(x,y)=\frac{1}{2\pi \sigma ^2}e^{\frac{x^2+y^2}{2\sigma ^2}}.$$
(1.1)
Define the random variable $`z=\frac{x}{y}`$. Obviously, its p.d.f. is independent of the width $`\sigma `$ of $`G(x,y)`$, and it is a straightforward exercise to show that

$$P(z)=\frac{1}{\pi }\frac{1}{1+z^2},$$
(1.2)
i.e., the standard Cauchy distribution.

A slightly more interesting generalization of (1.1) is to consider the family of joint probability density (j.p.d.) functions of the form

$$G(x,y)=f(x^2+y^2),$$
(1.3)
where $`f(u)`$ is a given appropriate p.d.f., subjected to the normalization condition

$$\underset{0}{\overset{\mathrm{}}{}}f(u)𝑑u=\frac{1}{\pi }.$$
(1.4)
A straightforward calculation of the p.d.f. of $`z=\frac{x}{y}`$ leads again to (1.2). Thus, the random variable $`z=\frac{x}{y}`$ is distributed according to (1.2), independently of the function $`f(u)`$. In other words, (1.2) is a universal probability density function.111We can generalize (1.3) somewhat further, by considering circularly asymmetric distributions $`G(x,y)=f(ax^2+by^2)`$ (with $`a,b>0`$ of course, and the r.h.s. of (1.4) changed to $`\sqrt{ab}/\pi `$), rendering (1.2) a Cauchy distribution of width $`\sqrt{b/a}`$, independently of the function $`f(u)`$. $`P(z)`$ is universal, essentially, due to rotational invariance of (1.3). More generally, $`P(z)`$ must be independent, of course, of any common scale of the distribution functions of $`x`$ and $`y`$.

We will now show that an analog of this universal behavior exists in random matrix theory. Our interest in this problem stems from the recent application of random matrix theory made in  to calculate the complexity of an analog computation process , which solves linear programming problems.

## 2 The universal probability distribution of the product $`B^1X`$ of real random matrices

Consider a real $`m\times (m+n)`$ random matrix $`A`$ with entries $`A_{i\alpha }(i=1,\mathrm{}m;\alpha =1,\mathrm{}m+n)`$. We take the j.p.d. for the $`m(m+n)`$ entries of $`A`$ as

$$G(A)=f(\mathrm{tr}AA^T)=f\left(\underset{i,\alpha }{}A_{i\alpha }^2\right),$$
(2.1)
with $`f(u)`$ a given appropriate p.d.f.. From222We use the ordinary Cartesian measure $`dA=d^{m(m+n)}A=_{i\alpha }dA_{i\alpha }`$. Similarly, $`dB=d^{m^2}B`$ and $`dX=d^{mn}X`$ for the matrices $`B`$ and $`X`$ in (2.6) and (2.18).

$$G(A)𝑑A=1$$
(2.2)
we see that $`f(u)`$ is subjected to the normalization condition

$$\underset{0}{\overset{\mathrm{}}{}}u^{\frac{m(m+n)}{2}1}f(u)𝑑u=\frac{2}{S_{m(m+n)}},$$
(2.3)
where

$$S_d=\frac{2\pi ^{\frac{d}{2}}}{\mathrm{\Gamma }\left(\frac{d}{2}\right)}$$
(2.4)
is the surface area of the unit sphere embedded in $`d`$ dimensions. This implies, in particular, that $`f(u)`$ must decay faster than $`u^{m(m+n)/2}`$ as $`u\mathrm{}`$, and also, that if $`f(u)`$ blows up as $`u0+`$, its singularity must be weaker than $`u^{m(m+n)/2}`$. In other words, $`f(u)`$ must be subjected to the asymptotic behavior

$$u^{m(m+n)/2}f(u)0$$
(2.5)
both as $`u0`$ and $`u\mathrm{}`$.
We now choose $`m`$ columns out of the $`m+n`$ columns of $`A`$, and pack them into an $`m\times m`$ matrix $`B`$ (with entries $`B_{ij}`$). Similarly, we pack the remaining $`n`$ columns of $`A`$ into an $`m\times n`$ matrix $`X`$ (with entries $`X_{ip}`$). This defines a partition

$$A(B,X)$$
(2.6)
of the columns of $`A`$.

The index conventions throughout this paper are such that indices

$$i,j,\mathrm{}\mathrm{range}\mathrm{over}1,2,\mathrm{},m,$$

$$p,q,\mathrm{}\mathrm{range}\mathrm{over}1,2,\mathrm{},n,$$
(2.7)
and $`\alpha `$ ranges over $`1,2,\mathrm{},m+n.`$

In this notation we have $`\mathrm{tr}AA^T=_{i,j}B_{ij}^2+_{i,p}X_{ip}^2=\mathrm{tr}BB^T+\mathrm{tr}XX^T`$, and thus (2.1) reads

$$G(B,X)=f(\mathrm{tr}BB^T+\mathrm{tr}XX^T).$$
(2.8)
We now define the random matrix $`Z=B^1X`$. Our goal is to calculate the j.p.d. $`P(Z)`$ for the $`mn`$ entries of $`Z`$. $`P(Z)`$ is clearly independent of the particular partitioning (2.6) of $`A`$, since $`G(B,X)`$ is manifestly independent of that partitioning. The main result in this section is stated as follows:

###### Theorem 2.1

The j.p.d. for the $`mn`$ entries of the real random matrix $`Z=B^1X`$ is independent of the function $`f(u)`$ and is given by the universal function

$$P(Z)=\frac{C}{[det(11+ZZ^T)]^{\frac{m+n}{2}}},$$
(2.9)
where $`C`$ is a normalization constant.

###### Remark 2.1

The probability density function (2.9) is a special (spherically symmetric) case of the so-called333Our notations in Remark (2.1) are slightly different from the notations used in . In particular, we interchanged their $`\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$, and also denoted their $`(TM)^T`$ by $`ZM`$ here. Finally, we applied the identity $`det(11+AB)=det(11+BA)`$ to arrive after all these interchanges from their equation (4.2.1) to (2.10). matrix variate $`t`$-distributions : The $`m\times n`$ random matrix $`Z`$ is said to have a matrix variate $`t`$-distribution with parameters $`M,\mathrm{\Sigma },\mathrm{\Omega }`$ and $`q`$ (a fact we denote by $`ZT_{n,m}(q,M,\mathrm{\Sigma },\mathrm{\Omega })`$) if its p.d.f. is given by

$$D(det\mathrm{\Sigma })^{\frac{n}{2}}(det\mathrm{\Omega })^{\frac{m}{2}}\left[det\left(11_m+\mathrm{\Sigma }^1(ZM)\mathrm{\Omega }^1(ZM)^T\right)\right]^{\frac{1}{2}(m+n+q1)},$$
(2.10)
where $`M,\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$ are fixed real matrices of dimensions $`m\times n`$, $`m\times m`$ and $`n\times n`$, respectively. $`\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$ are positive definite, and $`q>0`$. The normalization coefficient is

$$D=\frac{1}{\pi ^{\frac{mn}{2}}}\frac{\underset{j=1}{\overset{n}{}}\mathrm{\Gamma }\left(\frac{m+n+qj}{2}\right)}{_{j=1}^n\mathrm{\Gamma }\left(\frac{n+qj}{2}\right)}.$$
(2.11)
It arises in the theory of matrix variate distributions as the p.d.f. of a random matrix which is the product of the inverse sqaure root of a certain Wishart-distributed matrix and a matrix taken from a normal distribution, and by shifting this product by $`M`$, as described in . Our universal distribution (2.9) corresponds to setting $`M=0,\mathrm{\Sigma }=11_m,\mathrm{\Omega }=11_n`$ and $`q=1`$ in (2.10) and (2.11).

###### Remark 2.2

It would be interesting to distort the parent j.p.d. (2.1) into a non-isotropic distribution and see if the generic matrix variate $`t`$-distribution (2.10) arises as the corresponding universal probability distribution function in this case.

To prove Theorem (2.1), we need

###### Lemma 2.1

Given a function $`f(u)`$, subjected to (2.3), the integral

$$I=𝑑Bf(\mathrm{tr}BB^T)|\text{det}B|^n$$
(2.12)
converges, and is independent of the particular function $`f(u)`$.

###### Remark 2.3

A qualitative and simple argument, showing the convergence of (2.12), is that the measure $`d\mu (B)=dB|\text{det}B|^n`$ scales as $`d\mu (tB)=t^{m(m+n)}d\mu (B)`$, and thus has the same scaling property as $`dA`$ in (2.2), indicating that the integral (2.12) converges, in view of (2.5). To see that I is independent of $`f(u)`$ one has to work harder.

Proof. We would like first to integrate over the rotational degrees of freedom in $`dB`$. Any real $`m\times m`$ matrix $`B`$ may be decomposed as

$$B=𝒪_1\mathrm{\Omega }𝒪_2$$
(2.13)
where $`𝒪_{1,2}𝒪(m)`$, the group of $`m\times m`$ orthogonal matrices, and $`\mathrm{\Omega }=\mathrm{Diag}(\omega _1,\mathrm{},\omega _m)`$, where $`\omega _1,\mathrm{},\omega _m`$ are the singular values of $`B`$. Under this decomposition we may write the measure $`dB`$ as

$$dB=d\mu (𝒪_1)d\mu (𝒪_2)\underset{i<j}{}|\omega _i^2\omega _j^2|d^m\omega ,$$
(2.14)
where $`d\mu (𝒪_{1,2})`$ are Haar measures over the $`𝒪(m)`$ group manifold. The measure $`dB`$ is manifestly invariant under actions of the orthogonal group $`𝒪(m)`$

$$dB=d(B𝒪)=d(𝒪^{}B),𝒪,𝒪^{}𝒪(m),$$
(2.15)
as should have been expected to begin with.

###### Remark 2.4

Note that the decomposition (2.13) is not unique, since $`𝒪_1𝒟`$ and $`𝒟𝒪_2`$, with $`𝒟`$ being any of the $`2^m`$ diagonal matrices $`\mathrm{Diag}(\pm 1,\mathrm{},\pm 1)`$, is an equally good pair of orthogonal matrices to be used in (2.13). Thus, as $`𝒪_1`$ and $`𝒪_2`$ sweep independently over the group $`𝒪(m)`$, the measure (2.14) over counts $`B`$ matrices. This problem can be easily rectified by appropriately normalizing the volume $`𝒱_m=𝑑\mu (𝒪_1)𝑑\mu (𝒪_2)`$. One can show444One simple way to establish (2.16), is to calculate $`𝑑B\mathrm{exp}\frac{1}{2}\mathrm{tr}B^TB=(2\pi )^{\frac{m^2}{2}}=𝒱_m\underset{\mathrm{}}{\overset{\mathrm{}}{}}d^m\omega _{i<j}|\omega _i^2\omega _j^2|\mathrm{exp}\frac{1}{2}_i\omega _i^2`$ . The last integral is a known Selberg type integral . that the correct normalization of the volume is

$$𝒱_m=\frac{\pi ^{\frac{m(m+1)}{2}}}{2^m_{j=1}^m\mathrm{\Gamma }\left(1+\frac{j}{2}\right)\mathrm{\Gamma }\left(\frac{j}{2}\right)}.$$
(2.16)
Let us now turn to (2.12). The integrals over the orthogonal group in (2.12) clearly factor out, and we obtain

$$I=𝒱_m\underset{\mathrm{}}{\overset{\mathrm{}}{}}\underset{i=1}{\overset{m}{}}d\omega _i\underset{j<k}{}|\omega _j^2\omega _k^2|\left(\underset{i=1}{\overset{m}{}}\omega _i\right)^nf\left(\underset{i=1}{\overset{m}{}}\omega _i^2\right).$$
(2.17)
Finally, we change the integration variables in (2.17) to the polar coordinates associated with the $`\omega _i`$. The angular part of that integral is fixed only by dimensionality and by the factor $`_{j<k}|\omega _j^2\omega _k^2|\left(_{i=1}^m\omega _i\right)^n`$, and is thus independent of the function $`f(u)`$.

To prove that $`I<\mathrm{}`$ we need only consider integration over the radius $`r^2=_{i=1}^m\omega _i^2`$, since integration over the angles obviously produces a finite result. Using (2.3), we find that the radial integral in question is

$$\underset{0}{\overset{\mathrm{}}{}}𝑑rr^{m^2+nm1}f(r^2)=\frac{1}{2}\underset{0}{\overset{\mathrm{}}{}}u^{\frac{m(m+n)}{2}1}f(u)𝑑u=\frac{2}{S_{m(m+n)}},$$

independently of $`f(u)`$.

We are ready now to prove Theorem (2.1):

Proof. By definition555Our notation is such that $`\delta (X)=_{i=1}^m_{p=1}^n\delta (X_{ip})=_{p=1}^n\delta ^{(m)}(X_p)`$, $`X_p`$ being the $`p`$-th column of $`X`$.,

$`P(Z)`$ $`=`$ $`{\displaystyle 𝑑B𝑑Xf(\mathrm{tr}BB^T+\mathrm{tr}XX^T)\delta (ZB^1X)}`$ (2.18)
$`=`$ $`{\displaystyle 𝑑B𝑑Xf(\mathrm{tr}BB^T+\mathrm{tr}XX^T)|\text{det}B|^n\delta (XBZ)}.`$
Integration over $`X`$ gives:

$`P(Z)`$ $`=`$ $`{\displaystyle 𝑑Bf(\mathrm{tr}BB^T+\mathrm{tr}BZZ^TB^T)|\text{det}B|^n}`$ (2.19)
$`=`$ $`{\displaystyle 𝑑Bf[\mathrm{tr}B(11+ZZ^T)B^T]|\text{det}B|^n}`$
The $`m\times m`$ symmetric matrix $`11+ZZ^T`$ can be diagonalized as $`11+ZZ^T=𝒪\mathrm{\Lambda }𝒪^T`$, where $`𝒪`$ is an orthogonal matrix, and $`\mathrm{\Lambda }=\text{Diag}(\lambda _1,\mathrm{},\lambda _m)`$ is the corresponding diagonal form. Obviously, all $`\lambda _i1`$, since $`ZZ^T`$ is positive definite. Substituting this diagonal form into (2.19) we obtain

$`P(Z)={\displaystyle 𝑑Bf(\mathrm{tr}B𝒪\mathrm{\Lambda }𝒪^TB^T)|\text{det}B|^n}.`$
From the invariance of the determinant $`|\text{det}B𝒪|=|\text{det}B|`$ and of the volume element $`d(B𝒪)=dB`$ under orthogonal transformations we have:

$$P(Z)=𝑑Bf(\mathrm{tr}B\mathrm{\Lambda }B^T)|\text{det}B|^n.$$
(2.20)
Let us now rescale $`B`$ as $`\stackrel{~}{B}=B\sqrt{\mathrm{\Lambda }}`$. Thus,

$$det\stackrel{~}{B}=\sqrt{det\mathrm{\Lambda }}detB,\mathrm{and}d\stackrel{~}{B}=(det\mathrm{\Lambda })^{\frac{m}{2}}dB.$$
(2.21)
Finally, substituting (2.21) in (2.20) we obtain

$`P(Z)`$ $`=`$ $`{\displaystyle \frac{d\stackrel{~}{B}}{(det\mathrm{\Lambda })^{\frac{m}{2}}}f(\mathrm{tr}\stackrel{~}{B}\stackrel{~}{B}^T)\left(\frac{|\text{det}\stackrel{~}{B}|}{(det\mathrm{\Lambda })^{\frac{1}{2}}}\right)^n}`$ (2.22)
$`=`$ $`{\displaystyle \frac{C}{(det\mathrm{\Lambda })^{(m+n)/2}}}={\displaystyle \frac{C}{[det(11+ZZ^T)]^{(m+n)/2}}},`$
where $`C`$ is the normalization constant

$$C=𝑑Bf(\mathrm{tr}BB^T)|\text{det}B|^n$$
(2.23)
rendering

$$P(Z)𝑑Z=1.$$
(2.24)
$`C`$ is nothing but the integral (2.12). Thus, according to Lemma (2.1), $`C<\mathrm{}`$ and is also independent of the function $`f(u)`$.

###### Remark 2.5

The j.p.d. $`P(Z)`$ in (2.9) is manifestly a symmetric function only of the eigenvalues of $`ZZ^T`$, and thus, a symmetric function only of the singular values of $`Z`$.

###### Remark 2.6

From the normalization condition (2.24) we obtain an alternative expression for the normalization constant (2.23) as

$$\frac{1}{C}=\frac{1}{𝑑Bf(\mathrm{tr}BB^T)|\text{det}B|^n}=\frac{dZ}{[det(11+ZZ^T)]^{(m+n)/2}},$$
(2.25)
which is manifestly independent of the particular function $`f(u)`$, in accordance with Lemma (2.1). The integral over the matrix $`Z`$ can be reduced to a multiple integral of the Selberg type  over the singular values of the matrix $`Z`$, which can be carried out explicitly:

$$\frac{dZ}{[det(11+ZZ^T)]^{(m+n)/2}}=\pi ^{\frac{mn}{2}}\underset{j=1}{\overset{n}{}}\frac{\mathrm{\Gamma }\left(\frac{j}{2}\right)}{\mathrm{\Gamma }\left(\frac{m+j}{2}\right)}.$$
(2.26)
For particular choices of the function $`f(u)`$, we can use (2.25) to derive explicit integration formulas. For example, the function

$$f(u)=\frac{e^u}{\pi ^{m(m+n)/2}}$$
(2.27)
(i.e., the entries $`A_{i\alpha }`$ in (2.1) are i.i.d. according to a normal distribution of variance $`1/2`$) satisfies (2.3). Thus, we obtain from (2.25) that

$$𝑑Be^{\mathrm{tr}BB^T}|\text{det}B|^n=\pi ^{\frac{m^2}{2}}\underset{j=1}{\overset{n}{}}\frac{\mathrm{\Gamma }\left(\frac{m+j}{2}\right)}{\mathrm{\Gamma }\left(\frac{j}{2}\right)}.$$
(2.28)
Note that the integral on the left-hand side of (2.28) can also be reduced to a multiple integral of the Selberg type (this time, over the singular values of $`B`$), which can be carried out explicitly. The result is

$$\pi ^{\frac{m^2}{2}}\underset{j=1}{\overset{m}{}}\frac{\mathrm{\Gamma }\left(\frac{n+j}{2}\right)}{\mathrm{\Gamma }\left(\frac{j}{2}\right)}.$$

Since this must coincide with (2.28), we obtain the identity

$$\underset{j=1}{\overset{m}{}}\frac{\mathrm{\Gamma }\left(\frac{n+j}{2}\right)}{\mathrm{\Gamma }\left(\frac{j}{2}\right)}=\underset{j=1}{\overset{n}{}}\frac{\mathrm{\Gamma }\left(\frac{m+j}{2}\right)}{\mathrm{\Gamma }\left(\frac{j}{2}\right)}.$$
(2.29)
###### Example 1

For $`n=1`$, i.e., the case where $`X`$ and $`Z`$ are $`m`$ dimensional vectors, (2.9) simplifies into the $`m`$ dimensional Cauchy distribution

$$P(Z)=\frac{C}{(1+Z^TZ)^{(m+1)/2}}.$$
(2.30)
This is so because for $`n=1`$, the matrix $`ZZ^T`$ has $`m1`$ eigenvalues equal to 0, that correspond to the $`m1`$ dimensional subspace of vectors orthogonal to $`Z`$, and one eigenvalue equal to $`Z^TZ`$. Thus, $`det(11+ZZ^T)=1+Z^TZ`$. Eq. (2.30) then follows by substituting this determinant into (2.22).

## 3 The universal probability distribution of the product $`B^1X`$ of complex random matrices

The results of the previous section are readily generalized to complex random matrices. One has only to count the number of independent real integration variables correctly. In what follows we will use the notations defined in the previous section (unless specified otherwise explicitly). Thus, consider a complex $`m\times (m+n)`$ random matrix $`A`$ with entries $`A_{i\alpha }(i=1,\mathrm{}m;\alpha =1,\mathrm{}m+n)`$. We take the j.p.d. for the $`m(m+n)`$ entries of $`A`$ as

$$G(A)=f(\mathrm{tr}AA^{})=f\left(\underset{i,\alpha }{}|A_{i\alpha }|^2\right),$$
(3.1)
with $`f(u)`$ a given appropriate p.d.f.. From666We use the Cartesian measure $`dA=d^{2m(m+n)}A=_{i\alpha }d\mathrm{Re}A_{i\alpha }d\mathrm{Im}A_{i\alpha }`$, with analogous definitions for $`dB`$ and $`dX`$ below.

$$G(A)𝑑A=1$$
(3.2)
we see that $`f(u)`$ is subjected to the normalization condition

$$\underset{0}{\overset{\mathrm{}}{}}u^{m(m+n)1}f(u)𝑑u=\frac{2}{S_{2m(m+n)}}.$$
(3.3)
This implies that $`f(u)`$ must be subjected to the asymptotic behavior

$$u^{m(m+n)}f(u)0$$
(3.4)
both as $`u0`$ and $`u\mathrm{}`$.
As in the previous section, we choose a partition

$$A(B,X)$$
(3.5)
of the columns of $`A`$. Thus, $`\mathrm{tr}AA^{}=_{i,j}|B_{ij}|^2+_{i,p}|X_{ip}|^2=\mathrm{tr}BB^{}+\mathrm{tr}XX^{}`$, and thus (3.1) reads

$$G(B,X)=f(\mathrm{tr}BB^{}+\mathrm{tr}XX^{}).$$
(3.6)
We now define the random matrix $`Z=B^1X`$. Our goal is to calculate the j.p.d. $`P(Z)`$ for the $`mn`$ entries of $`Z`$. The main result in this section is stated as follows:

###### Theorem 3.1

The j.p.d. for the $`mn`$ entries of the complex random matrix $`Z=B^1X`$ is independent of $`f(u)`$ and is given by the universal function

$$P(Z)=\frac{C}{[det(1+ZZ^{})]^{m+n}},$$
(3.7)
where $`C`$ is the normalization constant (3.12).

Proof. The proof proceeds in a similar manner to the proof of Theorem (2.1). The only important difference is that now $`\delta (X)=_{i=1}^m_{p=1}^n\delta (\mathrm{Re}X_{ip})\delta (\mathrm{Im}X_{ip})=_{p=1}^n\delta ^{(2m)}(X_p)`$, $`X_p`$ being the $`p`$-th column of $`X`$. One obtains

$`P(Z)`$ $`=`$ $`{\displaystyle 𝑑B𝑑Xf(\mathrm{tr}BB^{}+\mathrm{tr}XX^{})\delta (ZB^1X)}`$ (3.8)
$`=`$ $`{\displaystyle 𝑑Bf[\mathrm{tr}B(11+ZZ^{})B^{}]|\text{det}B|^{2n}}`$
where we have integrated over $`X`$.

The $`m\times m`$ complex hermitean matrix $`11+ZZ^{}`$ can be diagonalized as $`11+ZZ^{}=𝒰\mathrm{\Lambda }𝒰^{}`$, where $`𝒰`$ is a unitary matrix, and $`\mathrm{\Lambda }=\text{Diag}(\lambda _1,\mathrm{},\lambda _m)`$ is the corresponding diagonal form. Obviously, all $`\lambda _i1`$, since $`ZZ^{}`$ is positive definite. Substituting this diagonal form into (3.8) we obtain

$$P(Z)=𝑑Bf(\mathrm{tr}B𝒰\mathrm{\Lambda }𝒰^{}B^{})|\text{det}B|^{2n}=𝑑Bf(\mathrm{tr}B\mathrm{\Lambda }B^{})|\text{det}B|^{2n},$$
(3.9)
where we used the invariance of the determinant $`|\text{det}B𝒰|=|\text{det}B|`$ and the invariance of the volume element $`d(B𝒰)=dB`$ under unitary transformations.

As in the previous section we now rescale $`B`$ as $`\stackrel{~}{B}=B\sqrt{\mathrm{\Lambda }}`$. Thus,

$$det\stackrel{~}{B}=\sqrt{det\mathrm{\Lambda }}detB,\mathrm{and}d\stackrel{~}{B}=(det\mathrm{\Lambda })^mdB.$$
(3.10)
Finally, substituting (3.10) in (3.9) we obtain

$`P(Z)`$ $`=`$ $`{\displaystyle \frac{d\stackrel{~}{B}}{(det\mathrm{\Lambda })^m}f(\mathrm{tr}\stackrel{~}{B}\stackrel{~}{B}^{})\left(\frac{|\text{det}\stackrel{~}{B}|}{(det\mathrm{\Lambda })^{\frac{1}{2}}}\right)^{2n}}`$ (3.11)
$`=`$ $`{\displaystyle \frac{C}{(det\mathrm{\Lambda })^{m+n}}}={\displaystyle \frac{C}{[det(11+ZZ^{})]^{m+n}}},`$
where $`C`$ is the normalization constant

$$C=𝑑Bf(\mathrm{tr}BB^{})|\text{det}B|^{2n}$$
(3.12)
rendering

$$P(Z)𝑑Z=1.$$
(3.13)
Finally, one can show, in a manner analogous to Lemma (2.1), that $`C<\mathrm{}`$ and that it is independent of the particular function $`f(u)`$.

There are obvious analogs to the remarks made in the previous section, which follow from Theorem (3.1), which we will not write down explicitly.

## 4 More on the distribution of solutions of systems of linear equations with random coefficients: extension of a result due to Girko

The methods of the previous sections may be applied in studying the distribution of solutions of systems of linear equations with random coefficients. For concreteness, let us concentrate on real linear systems in real variables.

Consider a system of $`m`$ real linear equations

$$\underset{\alpha =1}{\overset{m+n}{}}A_{i\alpha }\xi _\alpha =b_i,i=1,\mathrm{}m$$
(4.1)
in the $`m+n`$ real variables $`\xi _\alpha `$. With no loss of generality, we will treat the first $`m`$ components of the vector $`\xi `$ as unknowns, and the remaining $`n`$ components of $`\xi `$ as given parameters. Thus, we split

$$\xi =\left(\begin{array}{c}z\\ u\end{array}\right)$$
(4.2)
where $`z`$ is the vector of unknowns $`z_i=\xi _i`$ ($`i=1,\mathrm{}m`$) and $`u`$ is the vector of parameters $`u_p=\xi _{p+m}`$ ($`p=1,\mathrm{}n`$). Similarly, we split the matrix of coefficients

$$A=(B,X),$$

where the $`m\times m`$ matrix $`B`$ (with entries $`B_{ij}`$) and the $`m\times n`$ matrix $`X`$ (with entries $`X_{ip}`$) were defined in (2.6). Thus, we may rewrite (4.1) explicitly as a system for the $`z_i`$:

$$Bz=bXu.$$
(4.3)
If we consider an ensemble of systems (4.1), in which $`A`$ and $`b`$ are drawn according to some probability law, the unknowns $`z_i`$ become random variables, which depend on the parameters $`u_p`$. Girko proved the following theorem for a particular family of such ensembles:

###### Theorem 4.1

(Girko)
If the random variables $`A_{i\alpha }`$ and $`b_i`$ ($`i=1,\mathrm{}m;\alpha =1,\mathrm{}m+n`$) are independent, identically distributed variables, having a stable distribution law with the characteristic function

$$g(t;\alpha ,c)=e^{c|t|^\alpha },0<\alpha 2;c>0,$$
(4.4)
then the random variables $`z_i`$ ($`i=1,\mathrm{}m`$) are identically distributed with the probability density function

$$p(\zeta ;\alpha ,\beta )=\frac{2}{\beta }\underset{0}{\overset{\mathrm{}}{}}r\rho (\frac{r\zeta }{\beta };\alpha )\rho (r;\alpha )𝑑r,$$
(4.5)
where $`\rho (r;\alpha )`$ is the probability density of the postulated stable distribution, and

$$\beta =\left(1+\underset{p=1}{\overset{n}{}}|u_p|^\alpha \right)^{\frac{1}{\alpha }}.$$
(4.6)
The ratios $`z_i/z_j`$ ($`ij;i,j=1,\mathrm{}m`$) have the density $`p(r;\alpha ,1)`$.

In the special case $`\alpha =2`$ in (4.4), the random variables $`A_{i\alpha }`$ and $`b_i`$ are normally distributed. For this case Girko obtained

###### Corollary 4.1

If, under the conditions of Theorem (4.1), $`\alpha =2`$, then

$$p(\zeta ;2,\beta )=\frac{\beta }{\pi (\zeta ^2+\beta ^2)},$$
(4.7)
i.e., $`\zeta `$ follows a Cauchy distribution of width $`\beta `$.

When $`\alpha =2`$, the j.p.d. of the $`A_{i\alpha }`$ and the $`b_i`$ is

$$G(A,b)=(2\pi \sigma ^2)^{\frac{m(m+n+1)}{2}}e^{\frac{1}{2\sigma ^2}(\mathrm{tr}A^TA+b^Tb)},$$
(4.8)
which is a special case of the j.p.d.’s we have discussed in the previous sections. Thus, in the spirit of the discussion in the previous sections, we will study systems of linear equations (4.1) with random coefficients $`A`$ and inhomogeneous terms $`b`$ with j.p.d.’s of the form

$$G(A,b)=f(\mathrm{tr}A^TA+b^Tb),$$
(4.9)
with $`f(u)`$ a given appropriate p.d.f. subjected to the normalization condition

$$\underset{0}{\overset{\mathrm{}}{}}u^{\frac{m(m+n+1)}{2}1}f(u)𝑑u=\frac{2}{S_{m(m+n+1)}}.$$
(4.10)
Our goal is to calculate the j.p.d. $`P(z;u)`$ for the $`m`$ unknowns $`z_i`$. We summarize our main result in this section as

###### Theorem 4.2

If the random variables $`A_{i\alpha }`$ and $`b_i`$ ($`i=1,\mathrm{}m;\alpha =1,\mathrm{}m+n`$) are distributed with a j.p.d. given by (4.9), with $`f(u)`$ being any appropriate probability density function subjected to (4.10), then the random variables $`z_i`$ ($`i=1,\mathrm{}m`$) are distributed with the universal j.p.d. function

$$P(z;u)=C\frac{\beta }{(\beta ^2+z^Tz)^{\frac{m+1}{2}}},$$
(4.11)
independently of the function $`f(u)`$, where

$$\beta =\sqrt{1+u^Tu}$$
(4.12)
and $`C`$ is a normalization constant given by

$$C=𝑑A|\mathrm{det}B|f(\mathrm{tr}A^TA).$$
(4.13)
###### Remark 4.1

Note that Theorem (4.2) generalizes the case $`\alpha =2`$ of Girko’s result, Theorem (4.1), from the particular $`f(u)e^u`$ to a whole class of probability densities $`f(u)`$, and moreover, it determines for this class of distributions the (universal) j.p.d. of the $`z_i`$’s, and not only the distribution of a single component. Thus, it is an interesting question whether Girko’s result could be generalized also to other ensembles of systems of linear equations as well.

Proof. By definition, from (4.3),

$`P(z;u)`$ $`=`$ $`{\displaystyle 𝑑A𝑑bf(\mathrm{tr}A^TA+b^Tb)\delta (zB^1(bXu))}`$ (4.14)
$`=`$ $`{\displaystyle 𝑑A𝑑bf(\mathrm{tr}A^TA+b^Tb)|\text{det}B|\delta (Bz+Xub)}`$
$`=`$ $`{\displaystyle 𝑑Af[\mathrm{tr}A^TA+(A\xi )^T(A\xi )]|\text{det}B|},`$
where in the last step we integrated over the $`m`$ dimensional vector $`b`$ and used $`Bz+Xu=A\xi `$.

The last expression in (4.14) is manifestly invariant under $`𝒪(m)\times 𝒪(n)`$ orthogonal transformations

$$P(𝒪_1z;𝒪_2u)=P(z;u),𝒪_1𝒪(m),𝒪_2𝒪(n),$$
(4.15)
due to the invariance of the measure $`dA=dBdX=d(B𝒪_1)d(X𝒪_2)`$, the invariance of the determinant $`|\mathrm{det}B|=|\mathrm{det}B𝒪_1|`$, and the invariance of the trace $`\mathrm{tr}A^TA=\mathrm{tr}(B𝒪_1)^T(B𝒪_1)+\mathrm{tr}(X𝒪_2)^T(X𝒪_2).`$ With this symmetry at our disposal, we may simplify the calculation of $`P(z;u)`$ by rotating the vectors $`z`$ and $`u`$ into fixed convenient directions, e.g., into the directions in which only $`z_1`$ and $`u_1`$ do not vanish:

$$z_i^{(0)}=z_0\delta _{i1},u_p^{(0)}=u_0\delta _{p1}.$$
(4.16)
with

$$z_0=(z^Tz)^{\frac{1}{2}},u_0=(u^Tu)^{\frac{1}{2}}.$$
(4.17)
Thus, we obtain

$`P(z;u)`$ $`=`$ $`{\displaystyle 𝑑Af[\mathrm{tr}A^TA+(Bz^{(0)}+Xu^{(0)})^T(Bz^{(0)}+Xu^{(0)})]|\text{det}B|}`$ (4.18)
$`=`$ $`{\displaystyle 𝑑B𝑑Xf(S)|\text{det}B|},`$
where

$$S=\underset{j=2}{\overset{m}{}}B_j^TB_j+\underset{p=2}{\overset{n}{}}X_p^TX_p+(1+z_0^2)B_1^TB_1+2u_0z_0B_1^TX_1+(1+u_0^2)X_1^TX_1,$$
(4.19)
in which $`B_i`$ is the $`i`$-th column of $`B`$, and $`X_p`$ is the $`p`$-th column of $`X`$.

The bilinear form involving $`B_1`$ and $`X_1`$ in (4.19) may be diagonalized as

$$(1+\xi ^T\xi )\left(\frac{z_0B_1+u_0X_1}{\sqrt{\xi ^T\xi }}\right)^2+\left(\frac{u_0B_1z_0X_1}{\sqrt{\xi ^T\xi }}\right)^2,$$
(4.20)
where we have used $`z_0^2+u_0^2=\xi ^T\xi `$. We now perform a rotation in the $`B_1X_1`$ plane, followed by a scale transformation of the first term in (4.20), thus defining

$$B_1^{}=(1+\xi ^T\xi )^{\frac{1}{2}}\frac{z_0B_1+u_0X_1}{\sqrt{\xi ^T\xi }},X_1^{}=\frac{u_0B_1z_0X_1}{\sqrt{\xi ^T\xi }},$$
(4.21)
such that $`d^mB_1^{}d^mX_1^{}=(1+\xi ^T\xi )^{\frac{m}{2}}d^mB_1d^mX_1`$. We will also need the inverse transformation for $`B_1`$

$$B_1(B_1^{},X_1^{})=\frac{1}{\sqrt{\xi ^T\xi }}\left(\frac{z_0B_1^{}}{(1+\xi ^T\xi )^{\frac{1}{2}}}+u_0X_1^{}\right)$$
(4.22)
in order to express the matrix $`B`$ in terms of the primed column vectors:

$$\stackrel{~}{B}=(B_1(B_1^{},X_1^{}),B_2,\mathrm{},B_m).$$
(4.23)
Thus, using (4.20) - (4.23) and the trivial fact that $`dA=_{i=1}^md^mB_i_{i=p}^nd^mX_p`$, we obtain

$$P(z;u)=\frac{1}{(1+\xi ^T\xi )^{\frac{m}{2}}}𝑑Af(\mathrm{tr}A^TA)|\text{det}\stackrel{~}{B}|,$$
(4.24)
where we have removed the primes from the integration variables. We are not done yet, since $`\stackrel{~}{B}_1`$, the first column of $`\stackrel{~}{B}`$, depends on $`z_0`$ and $`u_0`$. To rectify this problem, we note from (4.22) that

$$\stackrel{~}{B}_1(B_1,X_1)=\sqrt{\frac{1+u_0^2}{1+\xi ^T\xi }}\left(B_1\text{cos}\theta +X_1\text{sin}\theta \right),\text{cos}\theta =\frac{z_0}{\sqrt{\xi ^T\xi (1+u_0^2)}}.$$
(4.25)
Thus, performing one final rotation by an angle $`\theta `$ in the $`B_1X_1`$ plane, which leaves, of course, $`dA`$ and $`\mathrm{tr}A^TA`$ invariant, we see that in terms of the rotated columns $`|\text{det}\stackrel{~}{B}|=\sqrt{\frac{1+u_0^2}{1+\xi ^T\xi }}|\text{det}B|`$, and thus, finally, we obtain that

$$P(z;u)=\frac{\sqrt{1+u_0^2}}{(1+\xi ^T\xi )^{\frac{m+1}{2}}}𝑑Af(\mathrm{tr}A^TA)|\text{det}B|,$$
(4.26)
which coincides with (4.11), due to (4.12) and (4.13).

###### Remark 4.2

The fact that the integral $`C=𝑑A|\mathrm{det}B|f(\mathrm{tr}A^TA)`$ is convergent and independent of the function $`f(u)`$ can be proved by decomposing $`A`$ into its singular values , essentially in a manner similar to our proof of Lemma (2.1), but with slight modifications in (2.13) and (2.14) due to the fact that $`A`$ is a rectangular matrix rather than a square matrix. We shall not get into these technicalities here, which the reader may find in . Note, however, that $`C`$ may be determined from the normalization of $`P(z;u)`$:

$$d^mzP(z;u)=1=C\beta \frac{d^mz}{(\beta ^2+z^Tz)^{\frac{m+1}{2}}}=\frac{C\pi ^{\frac{m+1}{2}}}{\mathrm{\Gamma }\left(\frac{m+1}{2}\right)}.$$

Thus,

$$C=\frac{\mathrm{\Gamma }\left(\frac{m+1}{2}\right)}{\pi ^{\frac{m+1}{2}}}=\frac{2}{S_{m+1}}.$$
(4.27)
###### Remark 4.3

We note that for $`n=0`$ (i.e., when $`A=B`$) and $`u=0`$, (4.3) degenerates into $`Bz=b`$, which is precisely the case $`n=1`$ (and $`zZ`$) in the conditions for Theorem (2.1), which we analyzed in Example (1). Thus, (4.11), evaluated at $`n=1`$ and $`u=0`$ must coincide with (2.30), as one can easily check it does.

Since Theorem (4.2) states the explicit form (4.11) of $`P(z;u)`$, we can now use it to derive, e.g., the probability density of the distribution of a single component $`z_i`$ and that of the ratio of two different components, mentioned in Girko’s Theorem (4.1):

###### Corollary 4.2

The $`m`$ components $`z_i`$ are identically distributed, with the probability density of any one of the components $`z_i=\zeta `$ given by (4.7) of Corollary (4.1).

Proof. That the $`z_i`$ are identically distributed is an immediate consequence of the rotational invariance of $`P(z;u)`$ in (4.11). The proof is completed by performing the necessary integrals:

$`p(\zeta ;\beta )`$ $`=`$ $`{\displaystyle 𝑑z_2\mathrm{}𝑑z_mP(z;u)_{|_{z_1=\zeta }}}=C\beta {\displaystyle \frac{dz_2\mathrm{}dz_m}{(\beta ^2+\zeta ^2+_{i=2}^mz_i^2)^{\frac{m+1}{2}}}}`$ (4.28)
$`=`$ $`{\displaystyle \frac{C\pi ^{\frac{m1}{2}}}{\mathrm{\Gamma }\left(\frac{m+1}{2}\right)}}{\displaystyle \frac{\beta }{\beta ^2+\zeta ^2}}.`$
Thus, from (4.27) we obtain the desired result that $`p(\zeta ;\beta )=\frac{\beta }{\pi (\beta ^2+\zeta ^2)}.`$ This result, should have been anticipated, since the universal formula (4.11) holds, in particular, for the the Gaussian distribution (4.8).

Finally, we have

###### Corollary 4.3

The ratios $`z_i/z_j`$ ($`ij;i,j=1,\mathrm{}m`$) have the density $`P(r)=p(r;1)`$.

Proof. The ratio $`z_i/z_j`$ is dimensionless, and thus its distribution cannot depend on the width $`\beta `$, which is the only dimensionful quantity in (4.11). The proof amounts to performing the necessary integrals, e.g., for the random variable $`z_1/z_2`$:

$`P(r)`$ $`=`$ $`{\displaystyle d^mzP(z;u)\delta \left(r\frac{z_1}{z_2}\right)}=C\beta {\displaystyle \frac{|z_2|dz_2\mathrm{}dz_m}{[\beta ^2+(r^2+1)z_2^2+_{i=3}^mz_i^2]^{\frac{m+1}{2}}}}`$ (4.29)
$`=`$ $`{\displaystyle \frac{2CS_{m2}}{r^2+1}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{m2}{2}\right)\mathrm{\Gamma }\left(\frac{3}{2}\right)}{2\mathrm{\Gamma }\left(\frac{m+1}{2}\right)}}={\displaystyle \frac{1}{\pi (r^2+1)}}=p(r;1).`$
## Acknowledgments

I thank Asa Ben-Hur, Shmuel Fishman and Hava Siegelmann for a stimulating collaboration on complexity in analog computation, which inspired me to consider the problems studied in this paper. I also thank Ofer Zeitouni for making some comments. Shortly after the first version of this work appeared, its connection with the theory of matrix variate distributions was pointed out to me by Peter Forrester. I am indebted to him for enlighting me on that connection and for suggesting references . This research was supported in part by the Israeli Science Foundation.


text:
# Shadowing in photo-absorption : role of in-medium hadrons

## Abstract

We study the effects of in-medium hadronic properties on shadowing in photon-nucleus interactions in Glauber model as well as in the multiple scattering approach. A reasonable agreement with the experimental data is obtained in a scenario of downward spectral shift of the hadrons. Shadowing is found to be insensitive to the broadening of the spectral functions. An impact parameter dependent analysis of shadowing might shed more light on the role of in-medium properties of hadrons.

PACS Nos. : 25.20.Dc,25.40.Ep,24.10.Ht,12.40.Vv

Theoretical studies based on various models of hadronic interaction predict a reduction in the mass of hadrons at or above nuclear matter density . The observation of enhanced dilepton production in the low invariant mass domain in the heavy ion collision experiments  does seem to indicate a non-trivial modification of the properties of light vector mesons, particularly the $`\rho `$ meson, in hot and/dense medium . However, the complicated dynamics both in the initial as well as final states in these experiments inhibits a firm conclusion about nuclear medium effects at present.

On the other hand, the experiments like photo-absorption on nuclei provide much cleaner systems for the study of in-medium properties of mesons . With the availability of better photon beams, there has been a renewed interest in the photo-absorption processes.

The phenomena of shadowing plays an important role in the photo-nuclear reactions. The photo-nuclear data at lower energies, for different nuclei, seem to indicate an early onset of shadowing . In ref. this feature has been interpreted as a signature for a lighter $`\rho `$\- meson in the medium, where the shadowing effect was evaluated within a Glauber- Gribov multiple scattering theory  and generalized vector meson dominance (VMD). In contrast the authors in ref. have claimed that the early onset of shadowing can be understood within simple Glauber theory  if one takes the negative real part of the $`\rho `$N scattering amplitude into account which corresponds to a higher effective in-medium $`\rho `$ meson mass. In a subsequent paper , the authors have concluded that the enhancement of shadowing at low energies occurs due to lighter $`\rho `$ mesons as well as intermediate $`\pi ^0`$ produced in non-forward scattering.

In the backdrop of these different inferences, we have made an attempt to understand the role of in-medium properties of hadrons in the phenomenon of shadowing in photo-nuclear reactions. The shadowing in photon-nucleus reactions can be written as,

$`{\displaystyle \frac{A_{eff}}{A}}={\displaystyle \frac{\sigma _{\gamma A}}{A\sigma _{\gamma N}}}`$ $`=`$ $`1+{\displaystyle \frac{\delta \sigma _{\gamma A}}{A\sigma _{\gamma N}}}`$ (1)
where $`\sigma _{\gamma A}=A\sigma _{\gamma N}+\delta \sigma _{\gamma A}`$ consists of the incoherent scattering of the photon from individual nucleons and a correction due to the coherent interaction with several nucleons.The later has been evaluated using both multiple scattering approach  as well as Glauber’s formula along with VMD . Resonance contribution to $`\gamma A`$ cross-section for photon energy $`1.2`$ GeV has been estimated by using the prescription given in .

Let us first consider the in-medium effects on the kinematics of $`\gamma A`$ collisions. For a photon incident on a nucleus with energy $`E_{\gamma L}`$, the energy in the rest frame of the nucleon is given by

$`E_\gamma =\gamma _FE_{\gamma L}(1\beta _F\mathrm{cos}\theta _L),`$ (2)
where $`\beta _F=p_F/E_F`$ and $`\theta _L`$ is the angle between the incident photon and the Fermi momentum of the nucleon which now depends on the space co-ordinate through the density $`n(r)`$. The square of the centre of mass energy $`s`$, of the $`\gamma N`$ system can then be written as,

$`s`$ $`=`$ $`(p_\gamma +p_F)^2`$ (3)
$`=`$ $`m_{N}^{}{}_{}{}^{2}+2\gamma _Fm_N^{}E_{\gamma L}(1\beta _F\mathrm{cos}\theta _L),`$ (4)
where $`m_N^{}`$ is the effective nucleon mass inside the nucleus. The modification of vector meson masses in nuclear environment has been studied in different models . Here, we have used two different models namely universal scaling scenario (USS)  and Quantum Hadrodynamical model (QHD) . In USS, the effective hadronic masses ($`m_H^{}`$) vary with nuclear density as

$$\frac{m_H^{}}{m_H}=10.2x,$$
(5)
where $`x=n(r)/n_0(r)`$, $`n_0(r)`$ being the normal nuclear matter density. In QHD the effective masses of nucleons and vector mesons are calculated using standard techniques of thermal field theory  and is parametrized as :

$`{\displaystyle \frac{m_H^{}}{m_H}}=1+{\displaystyle \underset{j=1}{}}a_jx^j.`$ (6)
For nucleons $`a_1=0.351277`$ and $`a_2=0.0766239`$; in case of $`\rho `$, $`a_1=1.30966`$, $`a_2=1.78784`$, $`a_3=1.17524`$ and $`a_4=0.294456`$ and finally for $`\omega `$, $`a_1=0.470454`$, $`a_2=0.313825`$ and $`a_3=0.0731274`$. No medium effect on the $`\varphi `$ meson is considered here as it is expected to be small .

At lower energies $`\gamma A`$ interaction is known to be dominated by resonance production . Beyond this region we have used VMD for a better description. The vector meson produced inside the nucleus will have an effective mass depending on the density of the nuclear medium as seen by the meson. This change in mass would affect the coherence length $`\lambda `$ which corresponds to the time scale of the fluctuation between the bare photon and the hadronic component of the physical photon. For small $`\lambda `$, the hadron mediated interaction may become indistinguishable from the bare photon interaction and there will not be any shadowing. In the present case, $`\lambda `$ becomes a function of the radial distance inside the nucleus. For a vector meson with effective mass $`m_V^{}`$, one gets,

$`\lambda ={\displaystyle \frac{1}{E_\gamma \sqrt{E_\gamma ^2m_{V}^{}{}_{}{}^{2}}}}{\displaystyle \frac{2E_\gamma }{m_V^2}},`$ (7)
where $`E_\gamma `$ itself depends on the position of the struck nucleon through eq. (2). Multiple scattering with the nucleons in the nucleus results in a modified mass of this meson. The correction to the nuclear photo-absorption cross section due to multiple scattering can be written in terms of $`n`$-fold multiple scattering amplitude $`𝒜^{(n)}`$ as ,

$$\delta \sigma _{\gamma A}=\frac{1}{2m_Nk}Im\underset{n=2}{\overset{A}{}}𝒜^{(n)}$$
(8)
$`k`$ is the wave vector of the photon and the $`n=1`$ term corresponds to the incoherent part. In the present work we show results upto double scattering only (see ref.  for details).

In the high energy limit, under the eikonal approximation, the summation of the multiple scattering series goes over to the Glauber’s formula if one neglects the width of the vector mesons. The photon entering the nucleus at an impact parameter $`b`$ produces a vector meson at position $`z_1`$. Inside the nucleus, the coherence length $`\lambda `$, in general, would be different at $`z_1`$ ($`\lambda _1`$) and $`z_2`$ ($`\lambda _2`$) as the different densities will yield different masses. The expression for the shadowing part of the cross section is then given by ,

$`\delta \sigma _{VA}={\displaystyle \frac{g_V^2}{4\pi \alpha }}\delta \sigma _{\gamma A}=`$ (10)
$`{\displaystyle \frac{1}{2kk_V}}{\displaystyle d^2b𝑑z_1𝑑z_2\mathrm{exp}\left[\frac{1}{2}\sigma _{VN}(z^{})n(b,z^{})𝑑z^{}\right]}`$
$`\times `$ $`k_V(z_1)\sigma _{VN}(z_1)k_V(z_2)\sigma _{VN}(z_2)n^{(2)}(b,z_1,z_2)`$ (11)
$`\times `$ $`[(\alpha _V(z_1)\alpha _V(z_2)1)\mathrm{cos}(({\displaystyle \frac{z_1}{\lambda _1}}{\displaystyle \frac{z_2}{\lambda _2}})`$ (12)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{z_1}^{z_2}}\alpha _V(z^{})\sigma _{VN}(z^{})n(b,z^{})dz^{})(\alpha _V(z_1)+\alpha _V(z_2))`$ (13)
$`\times `$ $`\mathrm{sin}(({\displaystyle \frac{z_1}{\lambda _1}}{\displaystyle \frac{z_2}{\lambda _2}})+{\displaystyle \frac{1}{2}}{\displaystyle _{z_1}^{z_2}}\alpha _V(z^{})\sigma _{VN}(z^{})n(b,z^{})dz^{})],`$ (14)
where $`\alpha _V=Ref_{VV}/Imf_{VV}`$ is the ratio of the real and imaginary part of the $`VN`$ forward scattering amplitude . $`\sigma _{VN}`$ is the $`VN`$ scattering cross section  and $`k_V`$ is the wave vector of the vector meson. The attenuation of the vector meson amplitude is described by the exponential factor. We have included 2-body correlation in the two-particle density as , $`n^{(2)}(b,z_1,z_2)=n(b,z_1)n(b,z_2)[1j_0(q_c|z_1z_2|)]`$
where $`q_c=780`$ MeV and $`j_0`$ is the spherical Bessel function.

The authors of ref.  have indicated that the scattering of the vector meson with the nucleons in the nucleus leads to a change in its mass ($`\mathrm{\Delta }m2\pi n(r)Ref/m`$ and concluded that using an external mass would mean an overcounting of the medium effects. This observation may not be valid entirely. Let us consider the QHD model to discuss this point. It is well known that the small increase in the vector meson mass due to its interaction with the Fermi sea is overwhelmed by the large decrease due to Dirac sea interaction. This statement, though model dependent, does point out that the vacuum fluctuation (VF) which in free space renormalizes the particle to its physical mass, may have a different role in the medium. Presently, this phenomena can not be described from first principles. The increase in $`\rho `$ mass due to negative real part of the $`\rho N`$ scattering amplitude ($``$ 10 - 100 MeV  depending on the parameterization) may not be effective enough due to the larger drop from VF corrections and the net decrease might show up in the experimental data. So, while considering the vector meson in the medium, one should consider both the effects. Furthermore, in QHD, a drop in nucleon mass causes a larger drop in $`\rho `$ mass. Hence, for the present study it is also necessary to consider the effective nucleon mass inside the nucleus which was ignored in the previous studies. We should mention here that any increase in mass will lead to a reduction in the shadowing because of the decrease in $`\lambda `$. Moreover, the experimental data from other sources, e.g. heavy ion collisions  and proton-nucleus collisions  seem to indicate a softening of the vector meson spectral function.

Before presenting the theoretical results we discuss the available experimental data. To get the experimental numbers for $`A_{eff}`$, we have used $`\sigma _{\gamma A}`$ from ref . and $`\gamma `$-proton cross section from refs .. The $`\gamma `$-neutron cross section is obtained as $`\sigma _{\gamma n}=\sigma _{\gamma d}\sigma _{\gamma p}+\sigma _{\gamma G}`$, where $`\sigma _{\gamma d}`$ is taken from ref. and $`\sigma _{\gamma G}`$ is the Glauber correction which is known to be small at lower energies . The data for $`\gamma p`$ and $`\gamma n`$ are interpolated for the relevant energies corresponding to given $`\sigma _{\gamma A}`$ . The average photon-nucleon cross section for a nucleus with mass number $`A`$ is given as

$`\sigma _{\gamma N}={\displaystyle \frac{Z\sigma _{\gamma p}+(AZ)\sigma _{\gamma n}}{A}}`$ (15)
from which the experimental numbers are obtained as,

$`{\displaystyle \frac{A_{eff}}{A}}={\displaystyle \frac{\sigma _{\gamma A}}{A\sigma _{\gamma N}}}.`$ (16)
We discuss the results now. Depending on the size of the nucleus we have used two different density distributions; for $`A<16`$ the shell model density profile of Ref.  and for heavier nuclei ($`A>16`$) the density profile from Ref.  has been used. According to eq. (2), $`E_\gamma `$ is a function of angle, $`\theta _L`$ for non-zero $`p_F`$. The results which are presented below have been averaged over all the angles. We find that the effect of Fermi momentum in the kinematics (through eq. (2)) is negligibly small.

The variation of $`\lambda `$ with $`E_\gamma `$ for $`\rho `$ mesons is plotted in Fig. 1 for Pb and Al nuclei. In order to take both the in-medium mass and width into account we have folded the coherence length with the spectral function (as in eq. (18) below). We observe that $`\lambda `$ is larger at lower values of $`E_\gamma `$ in the QHD scenario than in the universal scaling approach. Again the heavier nuclei seem to be affected earlier. These observations are crucial in the understanding of shadowing effects.

Fig. 2 shows the variation of $`A_{eff}/A`$ with $`E_\gamma `$ for different nuclei. The region below 1.2 GeV can be well described by the contribution of the baryonic resonances alone. For the region $`1.2<E_\gamma <3`$ GeV we have used VMD with (Glauber) and without (multiple scattering) eikonal approximation. We find that the USS gives a better description of the data both in the multiple scattering approach and Glauber model than the scenario with vacuum properties of hadrons. In QHD the drop of mass being larger, in general, the data is underestimated. It is also necessary to point out that though the resonance and VMD regions have been differentiated by a vertical line at $`E_\gamma =1.2`$ GeV, there might as well be an overlap of both the pictures around this point. On an average, our results show that the experimental data over the entire range of photon energy under consideration are reasonably well reproduced by the downward shift of the spectral function within the USS.

As mentioned before the shift in the hadronic spectral function in the nuclear medium is an unsettled issue. The experimental data on dilepton production from heavy ion collisions at CERN super-proton synchrotron energies can be explained either by shifting the pole mass ($`m_V`$) to a lower value or by increasing the width ($`\mathrm{\Gamma }_V(M)`$) of the spectral function. We would like to demonstrate here how these kind of medium effects (pole mass shift or broadening) affect the shadowing. For this we have considered the quantity,

$$A_{eff}=\frac{\rho (M)A_{eff}(M)𝑑M}{\rho (M)𝑑M}$$
(17)
where

$$\rho (M)=\frac{1}{\pi }\frac{M\mathrm{\Gamma }_V(M)}{(M^2m_V^2)^2+M^2\mathrm{\Gamma }_V^2(M)}.$$
(18)
$`A_{eff}(M)`$ gets maximum weight at the peak of the spectral function, i.e. from the point $`M^2=m_V^2`$ and the contribution from either side of this point is approximately averaged out. Therefore, the results become sensitive to the pole mass and is largely insensitive to the broadening of the spectral function. In fig. 3 we show the quantity $`A_{eff}/A`$ as a function of photon energy for $`\rho `$-meson only. As explained above the results for vacuum mass 770 MeV and widths 150 MeV (solid line) and 230 MeV  (dotted line) are indistinguishable. However, a pole mass shift in USS shows a larger shadowing as $`\lambda `$ increases substantially in this case. Moreover, an increase ($``$ 40 MeV ) in $`\rho `$ mass results in the decrease in shadowing (dot-dashed line).

In all the results shown above, the shadowing is evaluated by integrating over all the values of impact parameter. On the other hand, an incident photon passing through the nucleus peripherally would see less shadowing due to lower densities. To visualize it we have plotted $`d(A_{eff}/A)/d^2b`$ with impact parameter in fig. 4. While going from C to Pb, we observe that at the lower impact parameter shadowing is larger for lighter nuclei. This phenomena is a reflection of the nuclear density profile, which for lighter nuclei is larger in the core region compared to the heavier nuclei. It will be interesting to know whether one can define a centrality parameter for $`\gamma A`$ collisions as is usually done for heavy ion collisions ( percentage minimum bias etc.).

To conclude, we have studied the effects of in-medium properties of hadrons on shadowing in photo-absorption processes both in the framework of Glauber model and multiple scattering approach. The general pattern of experimental data seem to prefer a dropping vector meson mass scenario. The universal scaling appears to be closer to the data. The shadowing effect is insensitive to the spectral broadening of the vector meson in the nuclear medium. In contrast to the previous works the spatial dependence of the masses of both vector mesons and nucleon are considered here. The effect of Fermi motion is found to be small in the kinematics of the process. However, the effect of two-body correlation is important as its absence overestimates the data. We would also like to comment on QHD. The simple Walecka model, which we have used here has its own limitations (e.g. large incompressibility etc.). In this model, the reduction in the nucleon and vector meson masses is substantially larger than other models, which leads to large amount of shadowing.

Finally, we emphasize that our understanding of the shadowing phenomena vis-a-vis Glauber model may be improved through corrections to the approximations inherent in the model as well as by examining the model parameters critically. The leading correction to the Glauber model due to deviation from eikonal propagation gives rise to a correction $`(2P_{cm}^AR)^1`$ relative to the Glauber scattering amplitude . Here $`P_{cm}^A`$ is the centre of mass momentum of the nucleus and $`R`$ is its charge radius. In the present case the correction is rather small for $`E_\gamma >1`$ GeV. Moreover, a refinement of the Glauber model (multiple scattering) parameters, e.g., vector meson-nucleon scattering amplitude, two nucleon correlation etc. might give a good agreement with the data even with the vacuum properties of the hadrons. Hence a better estimate of these quantities is essential for a definitive statement regarding the role of medium effects on shadowing in photo-absorption processes. Experimental data with better statistics would certainly help us to resolve these uncertainties.


text:
# Do dusty A stars exhibit accretion signatures in their photospheres?

## 1 Introduction

IRAS observations have shown that about 20% of all nearby A stars are surrounded by dust (Cheng et al. Cheng92 (1992)); the most prominent prototypes of these dusty A stars are Vega and $`\beta `$ Pictoris. In order to understand the nature of these dusty A stars, Holweger & Rentzsch-Holm (HRH95 (1995)) and Holweger et al. (HHK99 (1999), in the following HHK99) searched for gas in the circumstellar (CS) environment of a sample of dusty and dust-free A stars. They found narrow absorption features in the Ca ii K and the Na i D lines in about 30% of the stars. The occurrence of these features shows a strong correlation with the projected rotational velocity, which suggests the presence of gas concentrated in a disk-like structure, and therefore proves that these imprints on the photospheric spectrum originate in the circumstellar rather than in the interstellar medium (Holweger & Rentzsch-Holm HRH95 (1995)).

Since main-sequence A stars possess extremely shallow surface convection zones, their atmospheres are stable enough to allow diffusion processes to occur (Michaud Mich70 (1970)). This leads to a variety of peculiarities like the AmFm-phenomenon, which has been specifically investigated in open clusters with well-known ages (see e.g. Alecian Ale96 (1996), Hui-Bon-Hoa & Alecian Hui98 (1998) and references therein). Moreover, the shallow convection zones allow any contamination from circumstellar or interstellar matter to show up in the photospheres of these stars. The metal underabundances of a small subgroup of A stars, such as the $`\lambda `$ Bootis stars, can be explained in terms of accretion of metal-deficient gas from their circumstellar environment (Venn & Lambert VL90 (1990)). In this subgroup of stars the dust phase is enriched in condensable elements whereas the gas phase is heavily metal-depleted, except for those elements which have a low condensation temperature, like C, N, O, and S. Conversely, in the framework of the ”planet migration“ model, the metal overabundances of solar-type stars with planets are attributed to the accretion of rocky planets onto the star (Smith et al. Smith (2001)). A recent study of non-interacting solar-type binary stars revealed the existence of pairs of stars, where one component was more metal-rich than the other (Gratton et al. Gratton (2001)). There is a trend towards increasing abundance difference with condensation temperature, pointing towards the accretion of dust-rich or rocky material onto the stellar surface from either the inner part of a protoplanetary disk or from rocky planets. Despite the wealth of literature dicussing the circumstellar disks and peculiar abundance patterns in main-sequence A stars, a definite connection between these phenomena has been discussed but never really proven.

A study on the surface composition of $`\beta `$ Pictoris carried out by Holweger et al. (Holweger97 (1997)) showed that this example of a star with CS matter, has solar abundances ’within a factor two or less’ and therefore reveals no signs of accretion. On the other hand Vega, another example of a star with CS matter, is often called a ”mild” $`\lambda `$ Bootis star (Venn & Lambert VL90 (1990); Lemke & Venn Lemke3 (1996)), because it shows moderate metal underabundances that can be explained by the accretion of gas depleted in condensable elements. In general this raises some questions when applied to normal A stars (i.e. stars without known peculiarities): Is the result obtained for $`\beta `$ Pic or Vega outstanding?: Can the presence of CS matter affect the abundances?: Is it possible to detect these signs of accretion?

These questions can be answered by an extended search for the differences between the abundance pattern of dusty and dust-free normal A stars. Therefore, we have carried out a detailed abundance analysis of oxygen, calcium, barium, yttrium, and iron in a sample of nearby dusty and dust-free normal A stars. These elements are selected as representatives of the light elements with low condensation temperatures (O), and the heavy elements with higher condensation temperatures (Ca, Ba, Y, and Fe). Ba and Ca are strongly affected by diffusion whereas Fe is not. We scrutinized our sample to detect systematic differences in surface composition between dusty and dust-free stars.

## 2 The stellar sample

We have compiled a sample of 21 northern sky object to compare the abundances of dusty and dust-free A stars. It contains 15 normal A stars taken from Cheng et al. (Cheng92 (1992)), who cross-correlated the catalogue of A stars within 25 pc (Woolley et al. W70 (1970)) with the IRAS Faint Source Survey. For comparsion, we have also observed 6 $`\lambda `$ Bootis stars where we expect to see typical accretion patterns in the atmospheres. Some of the objects that are accessible from ESO/Chile have been observed before in Ca ii K (Lemke Lem89 (1989), Stürenburg S93 (1993), Holweger & Rentzsch-Holm HRH95 (1995), HHK99, but nothing is known about their surface composition except the Ca abundance. We have added one star from the literature to our sample, namely $`\beta `$ Pic. Further details on the program stars and the spectral regions observed are given in Table 1.

## 3 Observations

We obtained high-resolution spectra with the 1.52m telescope at the Observatoire de Haute-Provence equipped with the Aurélie spectrograph (Gillet et al. Gil94 (1994)) on 4 nights (2000 January 18 – January 22). In the first two nights, we observed the Ca ii K line, using grating Nr. 6, centered at 3934 Å, at a resolution of R=110 000. In addition to studying the presence of narrow absorption features in the Ca ii K line, it can also be used to make a reliable determination of the projected rotational velocities. An exact knowledge of this parameter is indispensable in our abundance analysis, especially in spectral regions with many blend lines and for fast rotators. Furthermore, Ca has proven to be a sensitive indicator of diffusion processes. We observed Ba for the same reason. Sodium (Na i D1 and D2) observations were carried out for the two stars where we found narrow absorption features in Ca ii K during the first two nights. In addition, we observed two comparison stars with pure stellar profiles to correct for telluric lines. On the third and fourth night we observed Ba ii and the O i triplet at 7771-4 Å, using grating Nr. 5, centered at $`\lambda _c`$ = 4934 Å, and $`\lambda _c`$ = 7800 Å, with R=60 000. A target list, including the elements observed, is given in Table 1.

For eight stars (HR 804, HR 1483, HR 1666, HR 2491, HR 2763, HR 5531, HR 6378, and HR 8728) we have spectra of the Ca ii K line region taken in 1996 and 1997 with the ESO CES system, at a resolution of R=70 000 (see HHK99).

## 4 Data reduction

Data reduction was carried out using ESO MIDAS standard routines. We used Th–Ar spectra for the wavelength calibration of the Ca ii K and Ba ii lines, and Th–Ne spectra for the O i data (Odorico et al. OGP87 (1987)).

We have constructed a template to correct for telluric lines in the Na i D region by dividing the pure stellar spectrum of HR 4534 by its synthetic spectrum. This template contains only the telluric lines, and was used to clean up the spectra of the two stars HR 1989 and HR 6556.

## 5 Fundamental parameters

We obtained values of Teff and $`\mathrm{log}g`$ from uvby$`\beta `$ photometry, using an updated version of the UVBYBETA code (Moon & Dworetsky Moon (1986); Napiwotzki et al. Napiwotzki (1993)). Typical errors for these parameters were $`\mathrm{\Delta }`$ T$`{}_{\mathrm{eff}}{}^{}=200`$ K and $`\mathrm{\Delta }\mathrm{log}g=0.1`$ dex. The final parameters are given in Table 2, columns 3 and 4. Using these parameters and a solar composition (Anders & Grevesse Anders (1989)), we used the ATLAS9 code (Kurucz Kurucz1 (1992)) to calculate plane-parallel atmospheric models. For all stars, we assume a microturbulence, $`\xi `$, of $`3.0`$ km s-1, consistent with previous work (Holweger & Rentzsch-Holm HRH95 (1995); HHK99).

## 6 Abundance analysis

Abundance analysis is done by synthesizing a spectrum using the Kiel line formation code, LINFOR. We derive abundances for five elements, O, Ca, Ba, Fe, and Y, in the three spectral windows. Here, we describe the input line data and the non-LTE calculations for O i.

### 6.1 Line data

The line data used in this analysis is taken from the VALD database (Kupka (1999)) and is summarized in Table 3. We use a total of 176 lines to fit the spectra in the spectral range of the Ba and Y lines, but in most cases, due to heavy blending, the Fe abundance is derived from a strong blend of Fe lines at 4958 Å. We then use this Fe abundance to derive a Ba abundance from the Ba/Fe blend at 4934 Å, and derive the Y abundance from the Y/Ba blend at 4900 Å.

### 6.2 Non-LTE effects

We discuss here the non-LTE effects of the individual elements, although a detailed non-LTE abundance determination was only performed for the O triplet, where we expect large departures from LTE (Paunzen et al. Paunzen (1999)). We derived level populations for all the relevant energy levels of neutral O, using the Kiel NLTE code, the model atmospheres described in Sect. 5, and the O i model atom described in Paunzen et al. (Paunzen (1999)). Non-LTE corrections were typically of the order of $`0.5`$ dex (see Table 4).

For Fe non-LTE corrections have been shown to be well below $`+0.25`$ dex for all program stars (Rentzsch-Holm Rentzsch (1996)). We did not perform NLTE calculations in this case since the typical error from our abundance analysis was around 0.2 dex.

Non-LTE corrections of Ba have been shown to depend strongly on the Ba abundance itself (Lemke Lemke2 (1990)). In Lemke’s sample of normal A stars, the abundance corrections range from $`0.07`$ to $`+0.30`$ dex. The overabundances we derived in this work are up to 2.0 dex, so can never be removed by non-LTE corrections.

Since we observe Y ii and Ca ii, we do not expect significant departures from LTE, as these are the dominant ionization states of the atom. This has been justified by Stürenburg (S93 (1993)), who found a mean correction of 0.03 dex for the Ca ii K resonance line.

### 6.3 Abundances

The spectral abundances derived from this work are summarized in Table 4. It proved impossible to derive reliable Ba and Y abundances in the star HR 5793 due to its large rotational velocity and the poor S/N ratio in the spectral range of the Ba lines. Given the large number of Fe lines we were able to derive an Fe abundance. In HR 5351, we could only deduce the Ba abundance. The typical error in the abundance analysis is 0.02–0.05 dex, depending on the spectral quality. In addition, the uncertainty in the stellar parameters may introduce an error of the order of 0.2 dex.

In Vega, our Ca abundance, $`\mathrm{log}ϵ(\mathrm{Ca})=5.83`$, compares well with the value of 5.82, obtained by Lemke (Lemke2 (1990)). The non-LTE O abundance, \[O\]$`=0.17`$, derived from the triplet at 7771 Å, also agrees well with the LTE abundance derived from visible lines that are barely affected by non-LTE, \[O\]$`=0.18`$, (Qiu et al. Qiu (2001)). In Sirius, our Ca ii abundance of 5.97 could indicate that Lemke (Lemke2 (1990)) slightly overestimated the non-LTE effects in neutral Ca when deriving Ca abundances, $`\mathrm{log}ϵ(\mathrm{Ca},\mathrm{LTE})=5.65`$, $`\mathrm{log}ϵ(\mathrm{Ca},\mathrm{NLTE})=6.26`$.

We completed our data using the following abundance analyses; Lemke (Lem89 (1989), Lemke2 (1990)) for HR 2491 (Ba, Fe): Gigas (Gigas1 (1986), Gigas2 (1988)) for HR 7001 (Ba, Fe): Holweger et al. (Holweger97 (1997)) for HR 2020 (Ca, Ba, Fe): Dunkin et al. (Dunkin (1997)) for HR 8728 (Fe): Cowley & Aikman (Cowley1 (1980)) for HR 6378 (Y, Fe).

## 7 Circumstellar lines

Our analysis of the Ca ii K data gave evidence for narrow absorption features in the rotationally broadened stellar line profile of two stars, HR 1989 and HR 6556.

### 7.1 HR 1989

Narrow absorption lines had already been detected in HR 1989 by Stürenburg (S93 (1993)), Bohlender & Walker (Bohlender (1994)), and Hauck et al. (Hauck (1998)). The latter gives an equivalent width of 72.2 mÅ, at a redshift of 11.1 km s-1. Our data gives a similar equivalent width, of 65 mÅ, but at a blueshift of 1 km s-1 (Fig. 1a). Therefore we conclude that the absorption feature varies with radial velocity with respect to the star, proving that in this case it is of circumstellar origin.

Despite the poor S/N ratio of the spectrum, the Na i D lines in HR 1989 show strong narrow absorptions at the line centers. Their equivalent width are 72.2 and 33 mÅ for the D2 and D1 line respectively, but the uncertainty on these values is large since the stellar line fit is so poor. Sfeir et al. (Sfeir (1999)) obtained data for lines-of-sight towards 143 stars, one being HR 1989, to map the local bubble. They derived equivalent widths of 34.1 and 16.4 mÅ for the Na i D2 and D1 lines, respectively. Using our Ca observations and Sfeir et al. (Sfeir (1999)) Na observations, we obtained a Na i/Ca ii ratio that is lower than 1. This also hints at a circumstellar rather than interstellar origin of the features provided that we ignore possible variabilities between the different observing periods.

### 7.2 HR 6556

In HR 6556, HHK99 found a circumstellar Ca absorption with an equivalent width of 22 mÅ, at a blueshift of 35 km s-1. Our data gives a slightly smaller equivalent width of 14 mÅ, at a blueshift of 32 km s-1 (Fig. 1b). Therefore this feature appears to be stable within the error limits. The Na spectrum of HR 6556 is so noisy that we are not able to deduce the presence or absence of narrow circumstellar lines.

A study of the interstellar Na i density distribution in the solar neighbourhood revealed that the observed Na column density in HR 6556 did not agree with neighbouring values (Vergely et al. Vergely (2001)). In fact, the discrepancy is about a factor of 200. We interpret this as evidence for circumstellar gas around HR 6556.

## 8 Dusty versus dust-free A stars

To ensure that our sample was homogeneous, we removed all the newly discovered spectroscopic binaries (see Table 1), since the abundance determination in these objects underestimates the real abundances. We also excluded the $`\lambda `$ Bootis stars, because these belong to a subgroup of A stars where the mechanism giving rise to the metal-poor abundance pattern is not yet fully understood (Paunzen P99 (1999)). In addition some of the $`\lambda `$ Bootis stars are beyond 50 pc, and IRAS measurements give only upper limits on the infrared fluxes, so that we do not know whether they are dusty or not. The remaining stars were split into two groups, those with (5) and without (7) infrared excess.

Even small amounts of accreted material will show up in the photospheric composition of A stars, since they have shallow convection zones ($`10^9`$ M, Turcotte & Charbonneau Turcotte (1993)). Besides accretion, meridional mixing and diffusion are the most important processes affecting the stellar abundance pattern. While meridional mixing simply wipes out any inhomogeneities, diffusion acts selectively on certain elements, imprinting its own abundance pattern which differs from the accretion pattern.

Surprisingly, there is no difference between the abundance patterns of dusty A stars and those of dust-free A stars. Table 5 shows that the mean abundance patterns of both groups of stars agree within the abundance analysis error. Only for Y do we have a discrepancy, which can be explained by the fact that this element is only determined in very few stars.

We do not find any correlation between abundances and effective temperature or gravity. However there is a pronounced correlation between the Ba abundances and the rotational velocity, with the overabundances diminishing with increasing $`v\mathrm{sin}i`$ (Fig. 2b). According to Michaud (Mich70 (1970)), the overabundances of Ba are due to diffusion processes, and their disappearence at higher $`v\mathrm{sin}i`$ is a result of efficient meridional mixing at rotational velocities in excess of 100 km s-1 (Turcotte & Charbonneau Turcotte (1993)). The Y abundances show also a marginal increase towards lower $`v\mathrm{sin}i`$. Cowley (Cowley0 (1976)) noted that diffusion theory does not predict large Y overabundances, which is in good agreement with our findings. We expect mild underabundances of Ca from diffusion theory, a trend which is qualitatively seen for stars with low rotational velocities in Fig. 2a. In all three cases, no distinction can be made between the dusty and dust-free stars in our sample.

A plot of the specific element abundance versus \[Ba/Fe\] ratio is used as a diagnostic tool to detect diffusion signatures. The use of the ratio ensures that any star-to-star metallicity variation cancels, and the influence of accretion likewise is excluded — in the standard accretion model Fe and Ba are both entirely condensed in dust grains — so that we are left only with diffusion effects. Fig. 3a and 3b illustrate how Ca and Ba are affected by diffusion: Ba shows large overabundances, but there is a trend for Ca to be depleted in the photosphere in the presence of strong diffusion. This is in perfect qualitative agreement with the theoretical predictions of Michaud (Mich70 (1970)). As expected, we find that in our program stars, O and Fe are not affected by diffusion due to their high abundance and large number of lines (Fig. 3c).

We derive least-square fits for the correlations between the Ca, Ba and Fe abundances with the \[Ba/Fe\] ratio, using the dusty and dust-free stars (dashed-dotted lines in Fig. 3). The parameters of the individual fits are calculated assuming a typical abundance error of 0.2 dex for all stars. The results are summarized in Table 6. Based on this analysis, we can definitely exclude a pure accretion signature for the dusty A stars.

## 9 Discussion

We showed that dusty and dust-free A stars possess the same abundance patterns and fit both into the diffusion scenario. We do not detect any specific accretion signature in the stars with circumstellar dust, even though there are a few slow rotators amongst them. This result agrees with the findings of Holweger et al. (Holweger97 (1997)) for $`\beta `$ Pictoris, which, despite of its conspicuous circumstellar disk, shows no signs of accretion. Typical disk masses derived for the dusty stars of our sample range from $`10^8`$ to $`10^4`$ M. So in all cases the disk contains enough material to contaminate the stellar photosphere. In their study Cheng et al. (Cheng92 (1992)) excluded any confusion of the infrared-excess with a background source, which leaves us with a question: Why do the dusty A stars seem to be unaffected by the surrounding circumstellar dust?

There are two possibilities which may explain the non-detection of any peculiar abundance pattern in dusty A stars. First, the stars may accrete circumstellar matter at an extremely low rate, e.g. $`\stackrel{}{M}<10^{14}`$ M yr-1. This would mean that the accretion rate is too low to overcome diffusion and therefore cannot significantly contaminate the convection zone of the star. The second possibility is that there is no interaction between the star and the disk, implying the existence of a gap, similar to that which is apparant in the case of $`\beta `$ Pictoris and HR 4796.

###### Acknowledgements.

We thank the DFG for subsidizing this project by a travel and observing grant (KA 1581/1-1) and Helen Fraser for a careful reading of the manuscript.