text:
# Generation of Primordial Cosmological Perturbations from Statistical Mechanical Models

## I Introduction

A central problem in contemporary cosmology is the quantitative explanation of the inhomogeneity observed in the Universe at large scales (see e.g. pee93 ; padm ). Such inhomogeneity is probed both indirectly for early times via the fluctuations in temperature in the microwave sky, and directly at present in the distribution of matter in space. All currently standard cosmological models work within a paradigm in which these fluctuations are the result of (mostly) gravitational evolution operating on some initial very small fluctuations. Common to all such models is the assumption - supported in particular by observations of the microwave background cobe \- of a very specific form of these initial perturbations at large scale, i.e. Gaussian fluctuations with a Harrison-Zeldovich (HZ, often referred to as “scale-invariant”) power spectrum har ; zel . In a recent paper some of us hz  have discussed how, in a simple classification of correlated systems, this spectrum corresponds to highly ordered glass-like or lattice-like distributions. In particular such systems are characterised in real space by the striking behaviour of the variance in the mass, $`\mathrm{\Delta }M^2`$, contained in a given volume, e.g. in spheres of radius $`R`$, $`\mathrm{\Delta }M^2R^2`$, or $`R^2\mathrm{log}R`$, i.e. it is essentially proportional to the surface area of the sphere. This is to be contrasted with the Poisson behaviour (proportional to the volume $`R^3`$) characteristic of many equilibrium systems with short range interactions away from phase transitions or that characteristic of long-range correlations at a critical point in many systems (for which $`\mathrm{\Delta }M^2R^\alpha `$ with $`\alpha >3`$). Actually $`R^2`$ is the slowest growth possible for any isotropic translationally invariant distribution of points beck .

In this paper we examine the nature of the correlations in such super-homogeneous systems and compare them with equilibrium correlations in systems with long range interactions studied in statistical mechanics martin-yalcin ; levesque-weis-lebo ; review-martin . More precisely we consider the One Component Plasma (OCP), which is simply a system of charged point particles interacting through a repulsive $`1/r`$ potential, in a uniform background which gives overall charge neutrality (for a review, see OCP-review ). We discuss how the OCP can be modified so that its equilibrium correlations are precisely those considered in cosmology. This involves an appropriately designed attractive short range potential as well as a repulsive $`1/r^2`$ potential at large scales. The latter corresponds to a four dimensional Coulomb potential with the particles confined to three dimensions.

Our analysis of these parallels has more than theoretical interest. It also provides a new method for generating initial conditions (IC) for numerical studies of the formation of structure in cosmology HE ; white93 . We explain that this method avoids the problems associated with currently used algorithms which involve superimposing small perturbations on a lattice or glass-like distribution. The latter is understood to be a “sufficiently uniform” discretization of a constant density background, not as representing the perturbed distribution. The problem with this choice of initial state is that the unperturbed distribution has its own inherent fluctuations/correlations, which can play an important role in the evolution of the system melott ; bsl02 ; bjsl02 ; slbj02 .

## II The HZ spectrum in cosmology

Let us recall first the necessary essentials of modeling perturbations from a uniform density used in cosmology. These are described as a continuous stationary stochastic process (SSP), with the property that the volume average of observables, i.e. empirical averages, approach the ensemble averages in the large volume limit. The stationarity refers in this context to the statistical invariance under spatial translation. Moreover one assumes also statistical isotropy of the stochastic process, i.e. invariance under rotation of the ensemble quantities. Thus if $`\rho (\stackrel{}{r})`$ is the density field, we have

$$\rho (\stackrel{}{r})=\rho _0>0,$$
(1)
and we can define a reduced 2-point correlation function $`\stackrel{~}{\xi }(r)`$ by

$$\rho (\stackrel{}{r}_1)\rho (\stackrel{}{r}_2)\rho _0^2\left[1+\stackrel{~}{\xi }(r)\right]$$
(2)
where $`r=|\stackrel{}{r}_1\stackrel{}{r}_2|`$. Corresponding isotropies are assumed to hold for the higher order correlation functions. Further the probability distribution for the fluctuations in the initial uniform density are generally assumed to be Gaussian.

Much more used in cosmology than $`\stackrel{~}{\xi }(r)`$ is the equivalent $`k`$-space quantity, the power spectrum $`P(\stackrel{}{k})`$ which is defined as

$$P(\stackrel{}{k})=\underset{V\mathrm{}}{lim}\frac{|\delta _\rho (\stackrel{}{k})|^2}{V}$$
(3)
where $`\delta _\rho (\stackrel{}{k})=_Vd^3re^{i\stackrel{}{k}\stackrel{}{r}}\delta (\stackrel{}{r})`$ is the Fourier integral in the volume $`V`$ of the normalized fluctuation field $`\delta (\stackrel{}{r})=(\rho (\stackrel{}{r})\rho _0)/\rho _0`$. In a statistically isotropic SSP this depends only on $`|\stackrel{}{k}|`$ with

$$P(\stackrel{}{k})P(k)=_0^{\mathrm{}}\stackrel{~}{\xi }(r)\frac{\mathrm{sin}(kr)}{kr}4\pi r^2𝑑r.$$
(4)
It follows from its definition that $`P(k)0`$, and, in the case we are considering, its integral, which is equal to $`\stackrel{~}{\xi }(0)`$, is assumed to be finite.

In current cosmological models it is generally assumed that the power spectrum $`P(k)k`$ at small $`k`$. This is known as the HZ or “scale-invariant” power spectrum. It is believed to describe the “primordial” fluctuations at very early times, the putative remnants of a period of “inflation” prior to the ordinary Big Bang phase pee93 ; padm . The linear $`k`$ behaviour is cut-off by a short distance scale (needed to ensure integrability) i.e. for $`k`$ larger than some cut-off, $`k_c`$, $`P(k)`$ decreases faster than $`k^3`$. The reason for the appearance of this particular spectrum is tied to considerations about the cosmological model (see pee93 ; padm ). Thus in the homogeneous and isotropic Friedmann-Robertson-Walker (FRW) cosmology, the only characteristic length scale is the size of the horizon $`R_H(t)`$ (the region causally connected at a time $`t`$). The HZ spectrum then corresponds to the choice which gives (as discussed in hz  this condition is in fact only satisfied for a spherical Gaussian window function, and not in a real physical sphere)

$$\sigma _M^2(R_H(t))=\mathrm{constant},$$
(5)
independent of $`t`$,where $`\sigma _M^2(R)`$ is the normalized mass variance

$$\sigma _M^2(R)=\mathrm{\Delta }M^2/M^2.$$
(6)
The HZ spectrum then follows from a consistency criterion for the treatment of the perturbations in the model: any other spectrum will give mass fluctuations which dominate over the homogeneous background either at some time in the future or past.

All current cosmological models share the HZ spectrum (or something very close to it) as the “primordial” form of their perturbations. This spectrum evolves, at larger and larger scales, until the characteristic time at which the densities of matter and radiation energy are equal. This evolution which depends strongly on the particular model used may change completely the form of the spectrum for distances smaller than a characteristic scale $`k_{eq}^1`$ (i.e. $`k>k_{eq}`$). For larger physical scales ($`k<k_{eq}`$), however, the evolution is almost exactly the same in all models, and it leaves the primordial HZ form intact.

Viewed in the general framework of correlated processes a very crucial property of the HZ type power spectrum is that

$$\underset{k0}{lim}P(k)=0$$
(7)
which in real space implies that

$$\underset{R\mathrm{}}{lim}\frac{\mathrm{\Delta }M^2}{V(R)}=\underset{R\mathrm{}}{lim}_{V(R)}d^3r\stackrel{~}{\xi }(r)=0$$
(8)
where $`V(R)`$ is the volume of a sphere of radius $`R`$ (with arbitrary origin). Eq. (8) is to be contrasted with the behaviour of this integral in a Poisson distribution, when it yields a finite positive constant, and with that associated with power-law correlations in critical systems hz  for which the same integral of $`\stackrel{~}{\xi }(r)`$ diverges.

Systems satisfying Eq. (8) are thus “more homogeneous” than a Poisson type system, as can be seen when one considers the behaviour of the mass variance in spheres. For a spectrum such that $`P(k)k^n`$ for $`k0`$ and appropriately cut-off at large $`k`$ we have for large $`R`$

$$\sigma _M^2(R)\{\begin{array}{ccc}1/R^{3+n}\text{if}n<1\hfill & & \\ \mathrm{log}(R)/R^4\text{if}n=1\hfill & & \\ 1/R^4\text{if}n>1\hfill & & \end{array}$$
(9)
In terms of the non-normalized quantity $`(\mathrm{\Delta }M)^2\sigma _M^2(R)R^6`$, Eq. (8) says that for $`n>0`$ we have a slower increase of the variance as a function of $`R`$ than for Poisson fluctuations, corresponding to $`n=0`$, with the limiting behaviour corresponding to a variance which is proportional to the surface area of the sphere. These systems are thus characterised by surface fluctuations, ordered (or homogeneous enough, one could say) to give this very particular behaviour. In Fig. 1 are shown, in two dimensions, a Poissonian distribution and a super-homogeneous distribution. The relatively greater uniformity of the latter is clearly identifiable.

## III Distributions of points with surface growth of the squared variance

A simple cubic lattice is the simplest example of a discrete set of points which shows this limiting behaviour ($`\sigma _M^21/R^4`$) of the variance kendall . The result is not hard to understand: The standard deviation $`\sigma _M`$ is that measured by averaging over spheres of radius $`R`$ centered at a randomly chosen point of the unit cell. It is proportional to the typical variation of points in the volume, that in the case of a lattice is the square root of the average number of points in the last spherical shell of thickness equal to the lattice unit. A better (more statistically uniform) example of the same kind is the so-called shuffled lattice hz ; this is a lattice whose sites are independently randomly displaced by a distance $`x`$ in all directions from their initial position according to some distribution $`p(x)`$ which has a finite second moment. In this case we find $`P(k)k^2`$ at small $`k`$ and, consequently, again $`\sigma _M^2(R)1/R^4`$ at large $`R`$. The simple lattice, however, is not a SSP, and even the shuffled lattice, though it can be defined as a SSP, is not statistically isotropic because the underlying lattice structure is not completely erased by the shuffling gacs .

To construct a statistically isotropic and homogeneous particle distribution with such a behaviour of $`\sigma ^2(R)`$ is non-trivial. A particular example is the so-called “pinwheel” tiling of the plane radin1 ; radin2 . The generation algorithm for it is defined by taking a right angled triangle with sides of respective length one and two (and hypotenuse $`\sqrt{5}`$) and, at the first step, forming five similar square triangle of sides $`1/\sqrt{5}`$ and $`2/\sqrt{5}`$ respectively as shown in Fig. 2. At the second step we expand these new triangles, to the size of the original triangle, and repeat the procedure ad infinitum, so that they cover the plane completely. Finally, placing a point randomly inside each elementary triangle will give the superhomogeneous point distribution which is statistically isotropic (with a continuous power spectrum).

## IV A physical example: The One-Component Plasma

In the study of the one component plasma (OCP) one considers the equilibrium statistical mechanics (at inverse temperature $`\beta `$) of a system of charged point particles interacting through a Coulomb potential, in a continuous uniform background giving overall charge neutrality OCP-review . Taking the particles to carry unit mass and charge the microscopic mass density of the particles is given by

$$\rho (\stackrel{}{r})=\underset{i}{}\delta (\stackrel{}{r}\stackrel{}{r}_i)$$
(10)
where $`\stackrel{}{r}_i`$ is the position of the $`i`$th particle. The total microscopic charge density is

$$\rho _c(\stackrel{}{r})=\underset{i}{}\delta (\stackrel{}{r}\stackrel{}{r}_i)n$$
(11)
where

$$n=\rho (\stackrel{}{r}).$$
(12)
The two-point correlation function is defined as

$`\stackrel{~}{\xi }(\stackrel{}{r}\stackrel{}{r}^{})`$ $`=`$ $`{\displaystyle \frac{\rho (\stackrel{}{r})\rho (\stackrel{}{r}^{})}{\rho ^2}}1`$ (13)
$`=`$ $`{\displaystyle \frac{\delta (\stackrel{}{r}\stackrel{}{r}^{})}{n}}+h(\stackrel{}{r}\stackrel{}{r}^{}).`$
The first term comes from the diagonal terms $`\stackrel{}{r}=\stackrel{}{r}^{}`$ in $`\rho (\stackrel{}{r})\rho (\stackrel{}{r}^{})`$ and is characteristic of any point distribution.

It can be proven quite generally review-martin  that at high temperatures the OCP has a translation invariant isotropic distribution with $`\stackrel{~}{\xi }(\stackrel{}{r})`$ satisfying Eq. (8) with $`r=|\stackrel{}{r}|`$. In terms of $`h(r)`$ this gives

$$nh(r)d^3\stackrel{}{r}=1$$
(14)
Eq.(14) is usually referred to as a sum rule which originates from the perfect screening of each charge caused by the long-range nature of the Coulomb potential. To see this let us suppose that an external infinitesimal charge density $`\rho _{\mathrm{ext}}=ϵe^{i\stackrel{}{k}.\stackrel{}{r}}`$ of very long wavelength is applied to the system. It creates an external electric potential

$$\mathrm{\Phi }(\stackrel{}{r})=\frac{\rho _{\mathrm{ext}}(\stackrel{}{r}^{})}{|\stackrel{}{r}\stackrel{}{r}^{}|}d^3\stackrel{}{r}^{}=\frac{4\pi }{k^2}ϵe^{i\stackrel{}{k}.\stackrel{}{r}}$$
(15)
and a perturbation to the Hamiltonian

$$H_{\mathrm{ext}}=\rho _c(\stackrel{}{r})\mathrm{\Phi }(\stackrel{}{r})d^3\stackrel{}{r}=ϵ\frac{4\pi }{k^2}\rho _c(\stackrel{}{r})e^{i\stackrel{}{k}.\stackrel{}{r}}d^3\stackrel{}{r}.$$
(16)
Using linear response theory the induced charge in the system is given by

$$\rho _{\mathrm{ind}}(\stackrel{}{r}^{})=\beta \rho _c(\stackrel{}{r}^{})H_{\mathrm{ext}}$$
(17)
where the average is over the unperturbed statistical distribution. Thus, assuming that the applied charge is perfectly screened i.e. the system responds with an induced charge density $`\rho _{\mathrm{ind}}=\rho _{\mathrm{ext}}`$ we have, in the limit, $`k0`$,

$$ϵe^{i\stackrel{}{k}.\stackrel{}{r}^{}}\beta ϵ\frac{4\pi }{k^2}\rho _c(\stackrel{}{r}^{})\rho _c(\stackrel{}{r})e^{i\stackrel{}{k}.\stackrel{}{r}}d^3\stackrel{}{r}$$
(18)
and therefore, since

$$P(k)=\frac{1}{n}+\stackrel{~}{h}(k)$$
(19)
where $`\stackrel{~}{h}(k)`$ is the Fourier transform of $`h`$, we find that

$$P(k)\frac{k^2}{4\pi n^2\beta }$$
(20)
for small $`k`$.

The behaviour of the power spectrum at small $`k`$ is traceable through Eq. (15) as being simply, up to a factor $`\beta `$, the inverse of the Fourier transform, (in the sense of distributions) of the repulsive $`1/r`$ potential. It is evident that what one needs to obtain the same kind of behaviour but with $`P(k)k`$ at small $`k`$ is a repulsive potential of which the Fourier transform behaves as $`1/k`$, i.e. a repulsive $`1/r^2`$ potential. Before considering this in more detail we discuss some issues relevant to the cosmological context which we have passed over without comment so far.

## V Discrete vs. Continuous

An essential difference between the case of cosmological perturbations and the system just discussed is that the former refers to a continuous density field while the latter describes the correlation properties of a set of discrete points. To interpret the latter as giving a realization of the former we obviously need to specify explicitly how to relate the two. There is no unique prescription to pass from a discrete to a continuous field. A prescription simply corresponds to a regularisation of the Dirac delta with a function $`W_L(\stackrel{}{r})`$ with the property

$`W_L(\stackrel{}{r})=L^3W_o({\displaystyle \frac{\stackrel{}{r}}{L}}),{\displaystyle W_o(\stackrel{}{r})d^3\stackrel{}{r}}=1`$ (21)
where $`L`$ is the characteristic scale introduced by the regularisation e.g. the Gaussian

$$W_L(\stackrel{}{r})=\left(\frac{1}{\sqrt{2\pi }L}\right)^3\mathrm{exp}\left(\frac{r^2}{2L^2}\right).$$
(22)
For any finite value of $`L`$ we can define a continuous density field $`\rho _L(\stackrel{}{r})`$ as the convolution of this function with the discrete density field. The pair correlation function of the continuous field can then be written also as a double convolution integral of $`W_L(\stackrel{}{r})`$ and the correlation function in the discrete case. The singularity in the latter is thus also removed by the regularisation. The power spectrum of this continuous field is then simply given as

$$P_L(\stackrel{}{k})=|W_L(\stackrel{}{k})|^2P_D(\stackrel{}{k})$$
(23)
where $`W_L(\stackrel{}{k})`$ is the Fourier transform of the regularisation $`W_L(\stackrel{}{r})`$ and $`P_D(\stackrel{}{k})`$ is the power spectrum of the discrete field. In particular for the Gaussian smoothing function we have

$$P_L(\stackrel{}{k})=\mathrm{exp}(k^2L^2)P_D(\stackrel{}{k}).$$
(24)
For the OCP discussed above the correlation function, in the weak coupling limit, is approximately the Debye-Hückel formula

$$h(r)=\beta \frac{e^{\kappa r}}{r}\kappa ^2=4\pi \beta n$$
(25)
so that

$$P_D(\stackrel{}{k})=\frac{1}{n}\left[\frac{k^2}{\kappa ^2+k^2}\right]$$
(26)
which is a non-integrable function of $`k`$. Integrability is ensured by the cut-off imposed at large $`k`$ in the power spectrum through the smoothing function. An assumed property of the primordial density field in standard cosmological models is that it is Gaussian distributed. Since the density field is inherently positive this assumption of Gaussianity is more properly attributed to fluctuations (expected to be small) around the uniform density. Small fluctuations in the discrete OCP are in fact Gaussian martin-yalcin .

## VI Initial Conditions in N-body simulations

A context in which it is necessary to have a concrete realization of the density $`\rho (\stackrel{}{r})`$ in cosmological models is that of the N-body approach to the problem of structure formation HE ; white93 . In these numerical simulations an initial configuration (IC), which should represent the universe at early times, is evolved under its self-gravity. The goal of these simulations is an understanding of how structures grow in such models, and whether they are compatible with observations of the distribution of matter (as probed, primarily, by the distribution of galaxies and galaxy clusters). The initial continuous density field is thus necessarily represented by a discrete set of points, which, because the discretization scale does not represent the real underlying physical one, must be considered just as we have discussed above as representing the continuous density field via some appropriate smoothing.

In practice such a discretization is always generated in a very particular way EDFW ; white93 : points of equal mass are displaced from a “pre-initial” configuration in a way prescribed by the correlation field one wants to set up. This prescription can be understood as follows. Superimposing an infinitesimal displacement field $`\stackrel{}{u}(\stackrel{}{r})`$ on a uniform density $`\rho _o`$ one has

$$\frac{\rho (\stackrel{}{r})\rho _o}{\rho _o}\stackrel{}{}.\stackrel{}{u}(\stackrel{}{r})$$
(27)
and therefore in $`k`$-space the power spectrum of the correlated density field $`\rho (\stackrel{}{r})`$ produced is roughly proportional to $`k^2P_u(k)`$, where $`P_u(k)`$ is that of the displacement field. The “pre-initial” configuration used, which gives the initial unperturbed positions of the particles, is a discretization of the initial uniform density field $`\rho _o`$ and the displacement field is specified by the power spectrum of the continuous field which is desired as IC in the simulation. In practice the former is taken either to be an exact simple lattice or a “glassy” configuration white93 . When this latter configuration is used it is generated by evolving the system with the N-body code, but with the sign of gravity reversed. This corresponds to a system essentially like the OCP (at low temperature), and the long-time evolution brings it to a highly uniform configuration with normalized variance decaying at large scales as $`1/R^4`$ bsl02 .

In the generation of the IC in this way the fluctuations associated with the “pre-initial” configuration are neglected, or rather implicitly assumed to be negligible with respect to those which are introduced by the perturbations to it. The problem with this procedure is that it neglects precisely these “pre-initial” fluctuations. (If we assume the initial continuous field to have power-spectrum $`P_i(k)`$ what results from the infinitesimal shuffling is a power spectrum $`P(k)=P_i(k)+k^2P_u(k)`$). Such a neglect is justified if these fluctuations play no role in the evolution of the system, as one wants to see only the growth of structure coming from the perturbations superimposed on the uniform density field. Indeed the reason the lattice or “glassy” configuration is used is that they are configurations in which the gravitational force is (effectively) zero, whereas in a Poisson distribution for example (which is a priori as good a discretization of the uniform background) it is not. While the Poisson distribution will thus evolve and form structures even without additional perturbation, the perfect lattice will not. That the discretized lattice or glass solves this problem is however far from clear: they both represent unstable point configurations with respect to gravity. A small applied perturbation will thus in general give an evolution which depends on the correlations and/or fluctuations in the initial unperturbed configurations. In principle such effects can be kept under control by making the discretization of the original density sufficiently fine. In practice in N-body simulations the “non-linear regime” in which structure formation takes place is comparable to the discretization scale, and the role of such effects seems, at the very least, problematic melott ; bsl02 ; bjsl02 ; slbj02 .

One way in which this difficulty could be got around is by generating directly a discrete distribution whose regularisation is exactly the desired continuous field. Our observations on the OCP provides in principle a way of realizing this possibility.

## VII IC as equilibrium of a modified OCP

We have seen that the OCP equilibrium correlations give surface fluctuations ($`\mathrm{\Delta }M^2R^2`$), but with a power spectrum at small $`k`$ which goes like $`k^2`$. By considering instead a repulsive $`1/r^2`$ potential, whose Fourier transform is $`2\pi ^2/k`$, we obtain from Eq. (19)

$$P(k)\frac{k}{2\pi ^2n^2\beta },\mathrm{for}k0$$
(28)
and

$$h(r)\frac{1}{2\pi ^4\beta n^2}\frac{1}{r^4}$$
(29)
for $`r(2\pi ^2\beta m)^1`$. The change from the exponential decay of Eq. (25) to a power-law decay is a result of the different analyticity properties of the two power-spectra: the $`k^2`$ behaviour is analytic at the origin, guaranteeing a rapidly decaying behaviour of its Fourier transform, while the $`k`$ spectrum is not.

In the context of cosmological N-body simulations what one needs is not simply the primordial HZ power spectrum with some appropriate small scale cut-off: what is simulated is only a part of the cosmological evolution, starting from a time at which the initial spectrum of fluctuations is already significantly modified from its primordial form. While purely gravitational evolution at these early times does not modify the HZ spectrum, non-gravitational effects, present until the time when the universe becomes dominated by matter, do so significantly. The nature of these modifications depends on the details of the cosmological model, but in all cases it affects only $`k`$ larger than a characteristic $`k_{eq}`$ corresponding to the causal horizon at the “time of equality” (when matter starts to dominate over radiation). One then has a power spectrum of the form

$$P(k)=k^nf(k)$$
(30)
with $`n`$ exactly (or very close to) unity and $`f(k)`$ a model- dependent form for the spectrum with $`f(0)=a>0`$. In Figure 3

are shown the power spectra of two representative cosmological models pee93 ; padm . Both are CDM (“cold dark matter”) type. HDM (“hot dark matter”) models, now strongly disfavored, have a more abrupt (typically exponential) cut-off at larger $`k`$, due to the “wash-out” of small scale structure associated with their higher velocity dispersions.

To produce such a spectrum as the equilibrium one of an OCP like system requires further modification of the form of the potential. Just like the standard OCP, an unmodified $`1/r^2`$ potential will give, in the weak-coupling limit, a spectrum which becomes flat (i.e. Poisson like) at large $`k`$. This is simply due to the fact that in this limit the thermal fluctuations dominate over the potential at small scales, effectively randomizing the particle positions up to some scale. A crude guess of what potential would produce the behaviour of a typical cosmological model can be obtained by supposing that at small scales

$$1+h(r)e^{\beta V(r)}$$
(31)
which corresponds to completely neglecting collective effects ($`1+h`$ represents the relative probability compared to random of finding a particle at distance $`r`$ from a given one). Given that the desired fluctuations always have small amplitudes ($`|h|1`$), we would then need to be in the regime of temperature and scales such that $`|\beta V(r)|1`$, so that $`h(r)\beta V(r)`$. The potential should thus be attractive at smaller scales, as the system is positively correlated at those scales. A $`k^2`$ behaviour at large $`k`$, which is often used (see Figure 3) as an initial condition approximating cosmological models (CDM type) in this regime (beyond the “turn-over”), may be obtained from an attractive $`1/r`$ potential, e.g.

$$V(r)=\frac{1}{r^2}\eta \frac{e^{\mu r}}{r}.$$
(32)
By modifying the parameters $`\mu `$ and $`\eta `$, as well as the temperature, both the amplitude of the $`P(k)`$ and the location of a change from a $`P(k)k^2`$ to a behaviour $`P(k)k`$ for small $`k`$ can be controlled. The potential (32) is repulsive at short distances, but it may be necessary to make it more strongly repulsive in order to to ensure the system is not unstable to collapse ruelle . Such issues, as well as the practical numerical implementation of these observations will be investigated in forthcoming work. Once this method has been implemented to set up the initial conditions, it will naturally be interesting to understand the effect of this change in initial conditions on the dynamical evolution of gravitational N-body simulations.

## VIII Acknowledgments

We thank Dominique Levesque for useful discussions about the OCP and Charles Radin for information about pinwheel tilings. FSL thanks the Swiss National Science Foundation for post-doctoral support. The work of JLL was supported by NSF grant DMR 98-13268 and AFOSR grant AF 49620-01-1-0154.


format: text

text:
# External and intrinsic anchoring in nematic liquid crystals: A Monte Carlo study

## I Introduction

In confined nematic liquid crystals with a large surface–to–volume ratio the aligning effects of the confining surfaces are of great importance in determining the equilibrium director configuration Zumer . There are two major contributions to these surface aligning effects, the first one originating from direct interactions between nematic molecules and the solid substrate (external anchoring), while the second one is due to incomplete anisotropic nematic-nematic interactions in the vicinity of the sample surface (intrinsic anchoring). An understanding of these confinement related aligning mechanisms is of great importance not only from the fundamental point of view, but also, e.g., for the design and construction of liquid crystal–based optical devices. Within phenomenological approaches anchoring effects are usually characterized by two parameters: the preferred molecular alignment direction at the sample surface (the easy axis) and the free energy coefficient $`W`$ penalizing any deviation from this direction (the anchoring strength) Rapini . Here the $`W`$ coefficient depends on $`S`$, the nematic order parameter. This dependence seems to be strongly related to the specific properties of a given interface: for example, for a system of hard rods confined between hard walls one finds $`WS`$ Okano , while in experiments measuring anchoring inside polycarbonate membranes even $`WS^4`$ can be obtained Alenka . Given, moreover, $`KS^2`$ deGennes ; Cleaver  ($`K`$ denoting the Frank elastic constant), with decreasing $`S`$ the extrapolation length $`b=K/W`$ deGennes  may either increase (as seen experimentally in thermotropics Alenka ), or decrease (as obtained from simulations with hard particles Andrienko ). In the bulk, the value of $`S`$ is primarily determined by temperature, while close to an interface it may also be affected by incomplete bulk interactions, as well as by (dis)ordering effects of the (possibly rough) confining substrate. In the latter case $`W`$ is often assumed to be simply proportional to $`𝒲S_0S(0)`$, where $`𝒲`$ represents the surface coupling constant, $`S_0`$ the surface–imposed value of $`S`$ (determined by substrate roughness), and $`S(0)`$ its actual surface value Nobili92 .

While both the anchoring strength and the easy axis can be determined experimentally Blinov , it is also possible to deduce them from simulations based on pairwise intermolecular interactions. For example, intrinsic anchoring has been analyzed in a pseudomolecular continuum approach with ellipsoidal molecules Silvia , using a lattice approximation in the zero–temperature limit hex , or analyzing different types of interfaces in Gay–Berne systems Emerson ; Mills ; Bates97 . Furthermore, surface anchoring strength has also been measured in a system of hard ellipsoids in contact with a hard wall Andrienko ; Allen . The anchoring strength reported in most analyses shows that anchoring is rather strong and that the corresponding extrapolation length is of the order of a few molecular dimensions, while its experimental values are typically above 100 nm Blinov . On the other hand, a recent analysis of external anchoring — also based on a lattice model, but for nonzero temperature — was presented in Ref. Nikolai , showing that the extrapolation length can increase significantly when the nematic–to–isotropic (NI) transition is approached, as also seen experimentally Alenka . Moreover, a lattice gas approach has been adapted recently to study nematic interfaces Bates02 .

Motivated by these developments, in this paper we extend the analyses performed in Refs. hex ; Nikolai  to nonzero temperatures (employing Monte Carlo simulations) and repeat the measurement of the extrapolation length both for external and intrinsic anchoring, the former in presence of rough solid substrates and the latter for a free nematic interface. The term “rough” here refers to flat but partially disordered substrates. We will first briefly recall the features of the lattice model used in Ref. hex , then discuss the modifications needed to perform the present analysis, and finally, present and discuss the results.

## II Simulation model and anchoring measurement

We model a liquid crystalline slab of thickness $`d`$, using a modification of the well–known Lebwohl–Lasher (LL) model Lebwohl . In the LL model elongated nematic molecules are represented by freely-rotating spin-like particles that are attached to lattice points of a cubic lattice. Since the pairwise interaction energy used in the LL model is spatially isotropic and thus does not depend on the relative position of the particles, it cannot produce any orienting effects at a free nematic surface and is therefore not suitable for studies of intrinsic anchoring. Therefore, we choose a more general potential to model the interparticle interaction energy resulting from anisotropic van der Waals forces. These forces can give rise to the nematic phase Maier  even without hard rod repulsion. For two neighboring nematic particles $`i`$ and $`j`$ with orientations given by unit vectors $`𝐮_i`$ and $`𝐮_j`$, separated by $`𝐫_{ij}`$, we model the interaction potential as Nehring :

$`U_{ij}=ϵ\left[𝐮_i𝐮_j3\nu (𝐮_i{\displaystyle \frac{𝐫_{ij}}{r_{ij}}})(𝐮_j{\displaystyle \frac{𝐫_{ij}}{r_{ij}}})\right]^2`$ (1)
with $`ϵ>0`$, which for $`\nu =0`$ reduces to the standard LL potential, while for $`\nu =1`$ one has the induced dipole-induced dipole interaction. The spatial anisotropy parameter $`\nu `$ enables us to continuously vary the relative importance of the spatially anisotropic ($`𝐫_{ij}`$-dependent) contribution to the interaction law (1). Although the range of van der Waals potential is proportional to $`r_{ij}^6`$, for computational efficiency we consider only interactions between nearest neighbors. Thereby the intrinsic anchoring energy $`W`$ and the elastic constant $`K`$ are underestimated, but we expect the errors in the estimation of $`W`$ and of the extrapolation length $`K/W`$ not to exceed 20% and 30%, respectively.

In a nematic slab, boundary conditions are typically fixed by the interaction with solid walls, or, alternatively, through orienting effects near free nematic interfaces. In the simulation, each of the solid walls is represented by a layer of fixed particles $`𝐩_i`$ ($`|𝐩_i|=1`$), either all perfectly aligned or somewhat disordered with some residual order. The nematic–wall interaction is modeled via

$`U_{ij}^s=ϵ_s\left[𝐩_i𝐮_j\right]^2,`$ (2)
as promoted, e.g., by short–range steric forces, where, again, $`𝐩_i`$ and $`𝐮_j`$ are nearest neighbors on the lattice. A dimensionless anchoring strength parameter can be defined as $`w=ϵ_s/ϵ`$. In the case of a perfectly aligned surface one has $`𝐩_i=𝚷`$, i.e., all particles are aligned along the easy axis $`𝚷`$, while for disordered surfaces $`𝚷`$ represents the average orientation of $`𝐩_i`$. On the other hand, a free nematic interface is simply modeled through the missing–neighbor effect.

Instead of the cubic lattice used in the LL model, in present simulations we use a simple hexagonal lattice to model the liquid crystal, which is necessary to avoid unphysical bulk easy axes hex  for $`\nu 0`$. These axes arise as a direct consequence of using a lattice approximation for modeling an anisotropic liquid and are present for any spatially anisotropic potential in the cubic lattice. They are absent, however, in the hexagonal lattice, provided that the spins $`𝐮_i`$ are assumed to be two–dimensional vectors confined to hexagonal planes. The geometry of our sample is shown in Fig. 1. The $`z`$ axis is normal to the confining interfaces and the hexagonal planes are parallel to the $`xz`$ plane. We stress that our analysis must be restricted to low values of the anisotropy parameter $`\nu `$ in order to avoid solid-like periodic director profile solutions that are stable for $`\nu >0.3`$ hex .

The strength of any anchoring can be measured by imposing an elastic distortion so that the average surface molecular orientation deviates from the easy axis defined by the anchoring. The magnitude of this deviation can then be used to estimate the anchoring strength and the corresponding extrapolation length deGennes . The elastic distortion in a nematic slab can be imposed either by applying a magnetic field whose orientation must not coincide with the direction of the easy axis, or by antagonistic anchoring conditions at opposing surfaces. In the zero-temperature analysis of Ref. hex  the magnetic field method was used. On the other hand, for nonzero temperatures a strong enough magnetic field can enhance the degree of nematic order and even shift the NI phase transition, which can present additional difficulties in interpreting the results. Therefore, we decided not to use the magnetic field approach; instead, we consider a hybrid cell-like sample with antagonistic boundary conditions.

Consider the $`\nu =0`$ case with pure external anchoring first. The left ($`z=0`$) interface — where the anchoring will be measured — is chosen to be a solid wall, either perfectly ordered or somewhat disordered, promoting homeotropic alignment through weak external anchoring ($`w=0.1`$). At the same time, the right ($`z=d`$) interface is also taken to be a solid wall, yet with perfect planar alignment and strong external anchoring ($`w=1`$). For the case $`\nu 0`$, we replace the solid wall at $`z=0`$ with a free nematic interface, enabling us to study intrinsic anchoring alone. Note that for $`\nu 0.3`$ the easy axis of the free interface is homeotropic hex , thus providing the same confinement type as with external anchoring. As a result, in both cases a combined bend and splay elastic deformation is expected to appear in the sample. The deformation should be present as long as the sample thickness $`d`$ exceeds $`d_c=|(K/W)_0(K/W)_d|`$, where $`(K/W)_0`$ and $`(K/W)_d`$ are the extrapolation lengths corresponding to the effective anchoring on the left and the right wall, respectively, while for $`d<d_c`$ a uniform director structure is observed Barbero . In a hybrid cell close to NI transition a non-bent biaxial structure is also possible, consisting of two strata of uniform alignment, separated by a biaxially ordered layer Andreja . Recall that unlike in the original LL model in the present study spins $`𝐮_i`$ are two–dimensional vectors and hence unable to reproduce biaxiality.

For hybrid boundary conditions in the one–constant approximation the Frank elastic theory (assuming a constant degree of nematic order throughout the slab) predicts a perfectly linear director tilt angle profile $`\varphi (z)`$, where $`\varphi `$ is measured, e.g., with respect to the slab normal $`z`$. The torque balance condition at, e.g., the left surface ($`z=0`$) can be written as $`(d\varphi /dz)_0=\frac{1}{2}(W/K)_0\mathrm{sin}2\varphi (0)`$ and enables us to deduce $`(K/W)_0`$ from the measured values of $`(d\varphi /dz)_0`$ and $`\varphi (0)`$. Note that for small $`\varphi (0)`$ (i.e., strong enough anchoring) the above condition simplifies to $`(d\varphi /dz)_0=(W/K)_0\varphi (0)`$ and allows us to determine $`(K/W)_0`$ simply by extrapolating the profile $`\varphi (z)`$ graphically across the sample boundary to $`\varphi =0`$ corresponding to the homeotropic easy axis. Note also that if the degree of nematic order is subject to variations — which is usually the case near sample boundaries in a layer of thickness $`\xi `$ ($`\xi `$ denoting the nematic correlation length) — the profile $`\varphi (z)`$ may deviate from the above predicted linear behavior. In this case the extrapolation of the profile towards the surface must be performed from far enough in the bulk where the order parameter profile is constant, that is from $`\xi <z<d\xi `$. Note also that whenever $`(K/W)_0`$ approaches $`d`$ (for weak anchoring or in a thin sample), $`(d\varphi /dz)_0`$ and, consequently, $`(K/W)_0`$ are accompanied by a significant systematic error.

Our Monte Carlo (MC) simulations are now performed as follows. The size of the simulation box size was set to $`48^3`$ for a total of 105984 nematic particles, excluding the boundary particles $`𝐩_i`$ in the layers at $`z=0`$ and $`z=d`$. In case of disordered substrates particle orientations in the confining layers are generated following a probability distribution $`f(\varphi )\mathrm{exp}(𝒫\mathrm{cos}^2\varphi )`$, which gives a homeotropic easy axis ($`\varphi =0`$) and a $`𝒫`$-dependent degree of order $`S_0`$. For example, $`𝒫0`$ corresponds to a distribution close to isotropic, while $`𝒫\mathrm{}`$ yields a perfectly aligned substrate. Then, the hexagonal lattice is divided into three sublattices, as shown in Fig. 1, ensuring that the bonds between neighboring particles on the lattice never connect two particles from the same sublattice. Considering the simple hexagonal lattice as tripartite enables us to vectorize the simulation algorithm, which provides a significant speed-up in calculations. Further, in the $`x`$ and $`y`$ directions periodic boundary conditions are assumed. We start either from a random configuration in two dimensions (recall that $`𝐮_i`$ are restricted to lie in hexagonal planes), or from an equilibrated configuration at a temperature slightly higher than the simulated one, if this is available. Then we apply the standard Metropolis algorithm Metropolis . For our vectorized algorithm to work correctly, in each MC cycle we first attempt (and accept/reject) trial moves involving particles in the first sublattice and only then proceed to the second one, and after that to the third one. In generating a new trial configuration each time only a single particle is involved. We have typically performed $`10^5`$ MC cycles for equilibration, followed by $`10^5`$ production cycles to accumulate averages of interest.

We measure the extrapolation length by analyzing the director profile $`\varphi (z)`$. The $`\varphi (z)`$ dependence is calculated by accumulating the two independent components of the two–dimensional ordering matrix $`Q_{\alpha \beta }(z)=2u_i^\alpha u_i^\beta _z\delta _{\alpha \beta }`$, where $`\alpha `$ and $`\beta `$ can be either $`x`$ or $`z`$, and $`u_i^\alpha `$ represents the $`\alpha `$ component of the unit vector $`𝐮_i`$. The average $`\mathrm{}_z`$ is performed both over all particles in the layer centered at $`z`$ and over the production MC cycles. Then, the averaged ordering matrix is diagonalized and the positive eigenvalue is identified as the two–dimensional scalar order parameter $`s(z)`$. Accordingly, the corresponding eigenvector is the director. Note that in Ref. Nikolai  the extrapolation length in the LL model with three–dimensional spins was studied, using a different method to obtain the director profile. Rather than computing the ordering matrix, the polar angle averaged over all particles in a layer and over production MC cycles was computed; i.e., no information was kept about the azimuthal orientation of the particles. This method overestimates the actual average director tilt angle for small values of this angle. The strong temperature dependence of the extrapolation length found in Ref. Nikolai  is most likely an artifact of the incorrect method used in that paper.

## III External anchoring

We analyze the spatially isotropic intermolecular interaction first by setting $`\nu =0`$ in Eq. (1). Such a model is now similar to the standard LL model, yet it is characterized by a different coordination number (eight in the former model, six in the latter) and the dimensionality of nematic spins (two in the former model, three in the latter). These distinctions lead to a different balance between the decrease of internal energy and loss of orientational entropy upon going from the isotropic into the nematic phase, and thus shift the NI transition to a temperature higher than in the LL model. Monitoring temperature scans of internal energy, specific heat, and the order parameter $`s`$ (lowercase $`s`$ will be used throughout the text for the two–dimensional order parameter), for $`\nu =0`$ we estimate the dimensionless transition temperature in a $`48^3`$ bulk sample (with full periodic boundary conditions) to be $`T_{NI}^{}=k_BT_{NI}/ϵ1.52\pm 0.01`$. Here a dimensionless temperature scale, $`T^{}=k_BT/ϵ`$, has been introduced.

Consider the director and order parameter profiles \[$`\varphi (z)`$ and $`s(z)`$, respectively\] in a hybrid nematic cell where both substrates are smooth and hence impose perfect nematic order with $`s_0=1`$. Recall that for $`\nu =0`$ one is dealing exclusively with external anchoring. Fig. 2(a) shows $`\varphi (z)`$ profiles for different temperatures. One can readily observe that far enough from the NI transition the profiles approach a linear function (predicted also from Frank elastic theory deGennes ), with minor changes in slope $`|d\varphi /dz|`$ only close to substrates where the degree of order can vary. In particular, wherever $`s(z)`$ exceeds its bulk value $`s_b`$, the nematic becomes more difficult to deform, which results in a reduction in slope (and vice versa). Moreover, note that far enough from the NI transition \[in our simulations for values of reduced temperature $`\tau =(T_{NI}^{}T^{})/T_{NI}^{}0.07`$\] the $`\varphi (z)`$ profiles are essentially insensitive to changing temperature. Given $`w=0.1`$, as chosen above for the homeotropically anchored smooth left wall, the extrapolation length is estimated as $`b=11.6a(1\pm 7\%)`$, where $`a`$ is the lattice spacing. Significant deviations from the linear $`\varphi (z)`$ profile can be observed only close to $`T_{NI}^{}`$ \[for $`\tau 0.02`$\], when the nematic far enough from the walls melts and thereby avoids elastic distortion — see Fig. 2 (b). Then molecular alignment becomes homeotropic in the vicinity of the left surface, followed by a region of (nearly) isotropic liquid in the slab center, and by a region of planar alignment close to the right surface. Note that one can avoid this nematic “meltdown” by reducing the deformation strength, e.g., by setting $`0<\varphi (d)\pi /2`$ or by reducing $`w`$ at the $`z=d`$ wall. In this case $`b`$ is observed to be essentially $`T^{}`$-independent even up to $`\tau 0.02`$. Note, however, that thereby the immediate vicinity of the NI transition actually has not been probed: for a realistic liquid crystal with $`T_{NI}300`$ K, $`\tau 0.02`$ corresponds to temperatures as much as 6 K below the transition. Alternatively, the deformation strength can be reduced also by significantly increasing the system size. Further, in Fig. 2 (b) one can always observe an increase of $`s`$ when approaching the substrate at $`z=d`$ promoting planar alignment. Similarly, $`s`$ decreases close to the $`z=0`$ substrate because of a weaker nematic–substrate coupling ($`w=0.1`$ as opposed to $`w=1`$ for the $`z=d`$ substrate).

Turn now to cases with a rough homeotropic substrate, simulated by generating an ensemble of fixed particles $`𝐩_i`$ at $`z=0`$ with $`s_0<1`$. Qualitatively, the behavior of $`\varphi (z)`$ and $`s(z)`$ does not change, including the insensitivity of $`b`$ to variations of $`T^{}`$ in the temperature range presently accessible. However, as shown in Fig. 3 for $`T^{}=1.26`$, a smaller $`s_0`$ reflects in a reduced anchoring strength at $`z=0`$ and hence in a reduced slope $`|d\varphi /dz|`$, which, in turn, results in a larger $`b`$: for $`s_00.78`$ we find $`b=14.9a(1\pm 7\%)`$, for $`s_00.41`$, $`b=31a(1\pm 7\%)`$, and for $`s_00.15`$, $`b=50a(1\pm 48\%)`$ (for a summary see Table 1). Note that in the latter case $`b`$ already approaches the sample thickness $`d`$ and is already close to violating the stability condition $`db`$ for the deformed director structure Barbero . Moreover, assuming $`Ks_b^2`$ and $`Ws_0s(0)`$, the extrapolation length should scale as $`bs_b^2/s(0)s_0`$. The full quantitative agreement, however, turns out to be rather poor (even within error bars).

Note that in the vicinity of fairly disordered substrates the nematic becomes elastically softer than in the bulk, which increases the local slope of the director profile $`\varphi (z)`$ in comparison with the bulk. As the NI transition is approached, the thickness of this disordered layer ($`\xi `$) starts to grow, thereby indirectly affecting the bulk slope of $`\varphi (z)`$, which — in principle — could result in a more pronounced temperature dependence for $`b`$ that is extrapolated from the bulk. Note, however, that a significantly thicker sample than the present one is required to actually observe this scenario.

We also studied the influence of the hexagonal lattice structure and the confinement of $`𝐮_i`$ to the $`xz`$ plane on our results, by analyzing external anchoring in the original LL model where the spins $`𝐮_i`$ are three–dimensional vectors. Both the LL and hexagonal lattice models yield the same qualitative behavior: a nearly temperature-independent $`b`$ which increases in value as the surface roughness is increased.

## IV Intrinsic anchoring

Setting $`\nu 0`$, intrinsic anchoring appears at the interfaces, in addition to the external contribution present already for $`\nu =0`$. Note that changing $`\nu `$ affects not only anchoring, but also changes the elastic softness of the nematic Barbero2 : for $`\nu `$ not too large ($`\nu 0.78`$), the Frank elastic constant decreases with increasing $`\nu `$. Consequently, upon increasing $`\nu `$, the NI transition temperature is found to decrease (see Table 1). Therefore, when comparing $`b`$ for different values of $`\nu `$, it is appropriate to rescale the temperature $`T^{}`$ and compare measurements with comparable values of the reduced temperature $`\tau `$ (resulting in similar values of $`s_b`$).

To facilitate the analysis, we decided to remove the wall at $`z=0`$ allowing us to deal with intrinsic anchoring alone. When $`\nu 0.3`$ (but nonzero) the intrinsic easy axis remains homeotropic hex , and the hybrid cell-like geometry studied in the previous section is maintained. As in the case of pure external anchoring, we find that the intrinsic extrapolation length $`b`$ shows no temperature dependence (within estimated error; see Fig. 4). However, $`b`$ does depend on the anisotropy of the nematic-nematic interparticle potential: for $`\nu =0.05`$ and $`\nu =0.1`$ one finds $`b=16.3a(1\pm 13\%)`$ and $`b=4.6a(1\pm 10\%)`$, respectively (see Table 1 and Fig. 5). This trend can be attributed to an increase of the anchoring energy $`W`$, as well as to a simultaneous decrease of the elastic constant $`K`$ upon increasing $`\nu `$. Note that in the analysis of Ref. hex  $`b=8a(1\pm 13\%)`$ was found for $`\nu =0.1`$ at zero temperature. The disagreement with the current largely temperature–independent estimate can be attributed to considering only nearest neighbors in the present analysis, which — as already stated — results in an underestimation of $`b`$.

Note, moreover, that for $`\nu >0.1`$ extrapolation lengths are in the microscopic range (of the order of a few — up to 5 — lattice spacings $`a`$). On the other hand, experimental values of the extrapolation length are typically of the order of 100 nm or greater Blinov . We can obtain quantitative agreement between our results and experiments by a significant decrease of the $`\nu `$ parameter, as also suggested in Ref. hex . A small value of $`\nu `$ in Eq. (1) promotes parallel alignment, as is favored, e.g., by steric repulsions in a system of hard rods. A decrease in $`\nu `$ might therefore be regarded as an effective inclusion of steric repulsions which are absent in our lattice model.

## V Conclusions

In this paper we have studied the anchoring of a nematic liquid crystal to a solid substrate and to a free interface, using numerical simulations of a simple hexagonal lattice model of two–dimensional spins interacting through a spatially anisotropic potential. We focused on the roles of substrate roughness and of the spatial anisotropy of the interparticle potential. An elastic deformation was imposed on the simulation cell, allowing us to extract the surface extrapolation length $`b`$ from the director profile. Setting the anisotropy of the potential to zero allowed us to probe the strength of the external anchoring (the anchoring arising from the direct interaction of the nematic with the surface). For temperatures not too close to the nematic–isotropic transition the extrapolation length was found to be largely temperature independent, and to grow with increasing surface roughness (as characterized by the distribution of the local preferred axis on the surface). This increase is physically reasonable, reflecting the decrease in the overall anchoring of the nematic as the surface becomes more disordered. Qualitatively similar behavior was found when we simulated the original Lebwohl–Lasher model where the spins are three–dimensional.

With the full anisotropic interaction potential present, intrinsic anchoring arises due to the incomplete spin–spin interactions at the sample surface. For a free nematic interface (and external anchoring absent), we extracted the intrinsic surface extrapolation length. Like the extrapolation length in the presence of only external anchoring, this length is temperature–independent. We found that the strength of the intrinsic anchoring grows together with the elastic softness of the nematic as the interaction potential anisotropy is increased, leading to smaller values of the extrapolation length. Thus, obtaining agreement with the experimentally measured values of $`b`$, on the order of 100 nm, requires a relatively small value of the interaction anisotropy parameter $`\nu `$ approaching $`10^2`$.

## Acknowledgments

Computational work in support of this research was performed at the Theoretical Physics Computing Facility at Brown University. Financial support from the National Science Foundation under Grant Nos. INT–9815313 and DMR–0131573, and from the Slovenian Office of Science (Programmes No. P0-0503-1554 and 0524-0106) is gratefully acknowledged.


text:
# Theoretical models of ferromagnetic III-V semiconductors

\[

## Abstract

Recent materials research has advanced the maximum ferromagnetic transition temperature in semiconductors containing magnetic elements toward room temperature. Reaching this goal would make information technology applications of these materials likely. In this article we briefly review the status of work over the past five years which has attempted to achieve a theoretical understanding of these complex magnetic systems. The basic microscopic origins of ferromagnetism in the (III,Mn)V compounds that have the highest transition temperatures appear to be well understood, and efficient computation methods have been developed which are able to model their magnetic, transport, and optical properties. However many questions remain.

\]

It is hoped that the emerging field of semiconductor spintronics, which studies the controlled flow of charge and spin in a semiconductor, will lead to the development of new non-volatile, high-density, and high-speed information technologies. An important milestone in this field was the discovery five years ago of ferromagnetism in Mn-doped, p-type GaAs observed at temperatures in excess of 100 Kelvin . Ferromagnetism at room temperature with full participation of itinerant carriers would be a major breakthrough in semiconductor spintronics. With this aim, intensive material research is currently in progress on transition metal doped III-V semiconductors. In this brief review we present a snap shot of the theoretical part of this endeavor, which has progressed hand in hand with experimental developments.

A qualitative picture of the electronic structure of III-V diluted magnetic semiconductors (DMSs) was proposed by Dietl and coworkers . Their model is based on the internal reference rule  which states that energy levels derived from the d-shell of the magnetic ion are constant across semiconductor compound families with properly aligned bands. The application of this rule to III-V materials is illustrated in Fig 1 for ionized magnetic impurities substituted on cation sites . The position of the A2 level for Mn in GaAs, deep in the valence band, suggests that Mn in GaAs is 2+, that its d-shell is occupied by five electrons, and that incorporation of Mn in this and many other (III,Mn)V compounds will result a large density of valence band holes that can mediate ferromagnetic coupling between the $`S=5/2`$ Mn local moments. For low carrier densities, these valence band holes will be bound to the Mn ion, leading to shallow acceptor levels. This model of carrier-induced ferromagnetism is now fairly well established for (Ga,Mn)As as a result of electron paramagnetic resonance experiments , x-ray magnetic circular dichroism measurements , and magneto-transport data .

The acceptor impurity levels of Fe and Co  are unlikely to lead to high valence band hole concentrations in arsenides or antimonides, as shown in Fig. 1. The possible coexistence of acceptor and neutral magnetic impurities suggests that a ferromagnetic double-exchange mechanism can dominate in this case . The scenario, however, has not yet been confirmed experimentally. For nitrides and phosphides the ionization energy of the acceptor levels is even higher and the A2 Mn impurity level is not deep in the valence band, making it likely that d-charge fluctuations cannot be neglected. Experimentally, ferromagnetism with critical temperatures exceeding room-temperature has been claimed in (Ga,Mn)N  on the basis of hysteresis in magnetization measurements, however, the electronic configuration of Mn impurity in this compound has not yet been conclusively established, and these claims have not yet been widely confirmed. The observation of the $`p`$-type conduction achieved in some samples  seems to suggest that the simple picture of A1 Mn acceptor level according to Fig. 1 may not apply in this case.

Guided by the complex phenomenology very briefly discussed above and by available experimental data on III-V DMSs, several theoretical groups have performed detailed studies of these new ferromagnetic materials using both microscopic and effective Hamiltonian approaches. The results of some of this theoretical research are reviewed in the following paragraphs.

Ab-initio theories are an invaluable tool for studying the microscopic origins of ferromagnetism and for predicting electronic, magnetic, and structural ground-state properties of new compounds. The local density approximation (LDA) combined with disorder-averaging coherent-potential approximation (CPA) was used  to study GaAs and GaN DMSs doped with V, Cr, Mn, Fe, Co, and Ni. A ferromagnetic ground-state was predicted for V, Cr, and Mn doping while Fe, Co, and Ni impurities were predicted to lead to a spin-glass state. The structural properties of (Ga,Mn)As have been calculated in the dilute limit and for hypothetical zinc-blende MnAs using pseudoatomic-orbital-based LDA algorithms . The calculated lattice constant of MnAs is similar to the lattice constant of GaAs. This result appears to suggest that the relatively large compressive strains observed in (Ga,Mn)As epilayers grown on GaAs are caused by intrinsic defects rather than by the substitutional Mn.

The electronic structure of (Ga,Mn)As supercells with Mn in the substitutional position, with As-antisites, or Mn-interstitial defects have been studied by various methods within the LDA . Typical spin-polarized band and Mn 3d projected density-of-states results are shown in Fig. 2  for substitutional Mn. It has been conclusively established by the LDA calculations that in the presence of substitutional Mn, the Fermi energy crosses the valence band near the $`\mathrm{\Gamma }`$ point, introducing one hole per Mn. Clear theoretical evidence for As 4p-like band holes has been obtained in LDA+U calculations  and the Mn 3d density of states calculated by this method is in a good agreement with photoemission data . As-antisite and Mn-interstitial defects act as a double donor, according to the microscopic calculations , and, therefore, are the likely sources of carrier compensation observed in experiment. In summary, first-principles calculations support the picture of carrier induced ferromagnetism in some (III,Mn)V ferromagnets, but detailed predictions depend on the approximation used for exchange and correlations and also perhaps on technical details of the calculation. This activity is on-going and further progress should be expected.

Semi-phenomenological effective Hamiltonian theories in which the low-energy degrees of freedom are Mn ion S=5/2 spins and valence band holes, discussed in the following paragraphs, have been developed to make quantitative predictions for zero and finite temperature magnetic, transport and optical properties of bulk samples  and heterostructures . In (Ga,Mn)As, e.g., these models are partially justified by the negligible changes in the valence band structure of the host semiconductor observed in angle-resolved photoemission experiments . The key term in the effective Hamiltonian description is a Kondo-like exchange interaction, $`J_{pd}_I𝐒_I𝐬(𝐑_I)`$, which arises microscopically from sp-d hybridization and strong d-orbital Coulomb interactions. Here $`𝐒_I`$ is the $`S=5/2`$ local moment operator on Mn site $`I`$ and $`𝐬(𝐫)`$ is the multi-band envelope function hole spin-operator . The simplest version of this model combines mean-field theory with a virtual crystal approximation for the random distribution of Mn moments . Illustrative results of the model are presented in Fig. 3 for the critical temperature in (Ga,Mn)As and (In,Mn)As. These calculations and the model predictions for the strain-engineered magnetic anisotropy, anomalous Hall conductivity , domain structure , and magnetic circular dichroism  are in good agreement with experiment.

The reliability of the standard mean-field approximation has been examined by accounting for the Coulomb interactions among holes , which enhance the critical temperature, and for correlations in Mn moment orientations which reduce the energetic cost of incompletely aligned spin configurations  and therefore lower the critical temperature. The calculations show that heavy-hole - light-hole mixing is responsible for a relatively large spin stiffness, resulting in only $`20`$ % reduction of the critical temperature due to spin-wave fluctuations in arsenides and antimonides . Moreover, the reduction is practically canceled by the enhancement of the critical temperature due to hole-hole interactions (see inset of Fig. 3) which explains the success of the mean-field theory in these ferromagnets .

Theoretical models beyond the virtual crystal approximation have been used to study the effects of disorder on transport properties, and magnetic properties in general, of the DMSs. The Boltzmann equation with Born approximation scattering rates have provided estimates of the anisotropic magnetoresistance effect (AMR) of order $`110`$ % , in good agreement with experiment . AMR of a similar magnitude is predicted by Kubo formula for the ac-conductivity, as illustrated in Fig. 4. The key for understanding the transport and magnetic anisotropy effects is a strong spin-orbit coupling in the host semiconductor valence-band. The absolute dc-conductivity estimates (see the left inset of Fig. 4)  are typically only 2-10 times larger than experimental values . The discrepancy is likely partly due to inaccuracies in the model scattering amplitudes and multiple-scattering or localization effects omitted by this approximation. However experimental conductivities are still increasing as growth and annealing procedures are optimized so the discrepancy could be due mainly to unintended extrinsic disorder. Kubo formula calculations of infrared conductivities, illustrated in the right inset of Fig. 4, of metallic (Ga,Mn)As DMSs demonstrate that transverse f-sum rule measurements can be used to extract accurate values for the free carrier density . The hole-fluid effective Hamiltonian theory discussed above uses no free parameters and is expected to be reliable in samples with high critical temperatures where the holes are metallic and their interaction with Mn and other impurities is effectively screened.

More detailed studies of disorder effects in DMSs have combined Kondo description of the spin interactions with Monte Carlo techniques applied to both metallic  (high carrier density) and insulating  (low carrier density) regimes, and with the CPA  or spin-wave theory  for the metallic samples. The CPA calculations, which start from a tight-binding approximation band model and use RKKY theory of the effective Mn-Mn interaction, predict a factor of 3 enhancement of the critical temperature due to disorder for typical Mn and hole densities in (Ga,Mn)As DMSs . Other studies using dynamical mean field theory  suggest similar trends with some additional features when considering the effective $`J_{pd}`$ coupling as a variable. On the other hand, the random Mn distribution and spatial fluctuations of the magnetization suppress the ferromagnetic transition temperature in Monte Carlo simulations . The non-collinear nature of the ground state due to disorder has been pointed out in the spin-wave theory calculations . It was also suggested  that highly frustrated magnetic correlations may result from the strong spin-orbit interaction in the valence band and the spatial disorder.

The metal-insulator transition regime was explored by exact diagonalization techniques . Consistent with experiment, the studies predict a metal-insulator transition for 1% Mn doping in GaAs DMSs . It has been emphasized  that defect correlations play an important role in the transition and the calculations  have provided evidence in favor of the applicability of hole-fluid and impurity-band magnetic-polaron models on metallic and insulating sides of the transition respectively. Monte Carlo studies of the insulating phase have found spatially inhomogeneous magnetizations bellow the ferromagnetic critical temperature and an enhancement of the critical temperature due to disorder .

We acknowledge helpful discussions with Tomasz Dietl, and Jan Mašek who also provided figures presented in this review. The work was supported by the Grant Agency of the Czech Republic, DARPA, DoE, and the Welch Foundation.


text:
# X-Ray Spectroscopy and Atomic Data

## 1 Introduction

With the advent of the high-resolution grating spectrometers on board the Chandra and XMM-Newton x-ray observatories, spectroscopy has taken central stage in the rejuvenated field of X-Ray Astronomy. The increasing role of spectroscopy as a tool for astrophysical measurements has naturally drawn attention to the relevant atomic physics and atomic data, as those are directly linked with the ability to draw meaningful conclusions from observed spectra. The x-ray wavelength band between 1 and 100 Å, which is covered complementally by the Chandra and XMM-Newton spectrometers, contains a rich forest of spectral lines emitted by highly charged ions that form at electron temperatures of $`kT_e`$ = 0.1 - 3 keV. At these high temperatures, most of the cosmically abundant elements from C to Ni are stripped down to their K shell (i.e., $`n`$ = 1, $`n`$ being the principal quantum number) and emit relatively few spectral lines. Additionally, many strong lines of L-shell ($`n`$ = 2) ions of Si to Ni fall in this wavelength range. The eight ionization stages in the L shell can provide more precise information on the temperature structure of the source, independent of elemental abundances, than the two K-shell ionization stages are capable of.

## 2 The Central Role of Fe-L

The most prominent L-shell x-ray lines in astrophysical spectra are those of iron. Fe-L line emission has been used recently to probe the hot temperature structure in stellar coronal sources (Brinkman et al., 2001), supernova remnants (van der Heyden et al., 2002) and galaxy clusters (Peterson et al., 2001). In general, transitions for both K- and L- shell ions can be calculated with available atomic codes. Nevertheless, in the past, the Fe L-shell lines have been considered highly uncertain, presumably because of the increasing complexity of multi-electron atoms.

To address these issues, an x-ray spectroscopy laboratory-astrophysics program was initiated 11 years ago and is still active today (PI’s Kahn and Beiersdorfer). The program is built around the unique capabilities of the LLNL EBIT to measure electron ion interactions. Driven by the astrophysical motivation, these efforts have, so far, focused primarily on the Fe L-shell complex. In particular, high-precision measurements for wavelengths and collisional excitation cross-sections (relative and absolute) have been published. Also, peculiar line ratios of Fe16+ that have been puzzling astronomers for years have been investigated and much of the atomic uncertainty has been disentangled from the real astrophysical issues. For more details see the contribution of Brown et al. in these proceedings.

As the first astrophysical grating spectra became available, our team has made a systematic attempt to test the ability of existing models to reproduce the observed emission line intensities. It was found (Behar, Cottam, & Kahn, 2001a; Brinkman et al., 2001) that, generally, the observed line intensities could be fairly well reproduced by state-of-the-art distorted wave calculations. For instance, the HULLAC code (Bar-Shalom, Klapisch, & Oreg, 2001) was used in those works. In particular, the Fe-L line intensities in the model including all of the high-$`n`$ lines were found to fare quite well, implying that the calculated excitation rates and ensuing radiative cascades are fairly adequate. Conversely, accurate wavelengths still needed to be incorporated from laboratory measurements. A more detailed confrontation of the atomic calculations with Fe-L spectra of stellar coronae can be found in Behar et al. (2001b), where a comparison of calculations with the spectra of bright stellar coronal sources such as Capella and HR 1099 confirms that, just like the K-shell atomic data, the Fe L-shell data, when calculated correctly, are highly reliable and therefore very useful. The latest versions of the widely used databases MEKAL (Mewe, Kaastra, & Liedahl, 1995) and APEC (Smith et al., 2001) now incorporate similar HULLAC data calculated by D. Liedahl with an earlier version of the code.

## 3 Remaining Atomic Data Issues

Although in general state-of-the-art models are doing well, several discrete, but nonetheless important, discrepancies still remain. In particular, the ratios of the 3s - 2p line intensities relative to those of the 3d - 2p transitions of the same charge state were found to be anomalously high for both Fe16+ and Fe17+. Similar effects have also been found in other late-type stars, and even in elliptical galaxies (Xu et al., 2002), which suggests that they are not associated with the astrophysical conditions. The origin of this discrepancy has been recently studied in many theoretical and experimental works (Laming et al., 2001; Doron & Behar, 2002; Beiersdorfer et al., 2002), not all in agreement with each other. Additionally, the L-shells of other cosmically abundant elements remain largely unexplored. Major efforts in this direction are being conducted by Lepson et al. (see these proceedings).

This is the place to note that the atomic processes associated with the x-ray emission depend on the type of plasma at the source. So far, we have focused on hot, collisional plasmas that are governed by electron impact ionization and excitation. This type of plasma is relevant to x-ray observations of stellar coronae, supernova remnants, the hot ISM in old galaxies, the intergalactic haloes of galaxy clusters, and to the overwhelming majority of laboratory plasmas. On the other hand, x-ray sources that are ionized and excited by an external radiation field, e.g., active galactic nuclei (AGN) and x-ray binaries, require atomic data for photon impact processes. Modeling line emission from photoionized sources involves radiative recombination and photoexcitation rates (i.e., oscillator strengths). In a recent paper on the x-ray spectrum of the type II AGN NGC 1068 (Kikhabwala et al., 2002), we have shown that the available data for the K-shell ions are very good for reproducing the x-ray spectra from photoionized sources. Work is in progress to test the status of the Fe-L data for photoionized sources. Similar data are routinely used to model absorption spectra. In absorption, however, one also observes inner-shell transitions, for which until recently there was an enormous lack of atomic data.

## 4 The Urgent Need for Inner-shell Absorption Measurements

The grating spectrometers on board Chandra and XMM-Newton enable us for the first time to detect x-ray absorption lines due to inner-shell photoexcitation. Since gas under almost any conditions absorbs x-rays, these lines are ubiquitous to x-ray absorption observations. Detections of inner-shell absorption lines have been reported mostly in the ionized outflows of AGN, but also for absorption by neutral ISM. Inner-shell absorption can probe the entire range of ionization states from neutral up to highly-ionized Li-like species. Consequently, these lines impose unprecedented, strict, constraints on the ionization structure in the absorbing medium. Since there were extremely few relevant atomic data for these features in the literature, in order to analyze the spectra, we had to embark on a tedious endeavor of calculating numerous lines (Behar, Sako, & Kahn, 2001c; Behar & Netzer, 2002). Calculations by the Ohio State team have also contributed to this effort (Pradhan, 2000; Nahar, Pradhan, & Zhang, 2001; Pradhan et al., 2002). None of these calculations have been benchmarked in the laboratory.

One case where x-ray absorption lines (including inner-shell excited) are particularly useful is for measuring velocities in AGN outflows (e.g., Kaspi et al., 2001). For these measurements, the rest frame wavelengths of the lines need to be known to very high accuracies. The case of inner-shell K$`\alpha `$ absorption by oxygen is particularly interesting because it could potentially relate the traditional x-ray absorber (O6+ and O7+) with absorbers of other wavebands (e.g., O5+ in the UV). Whether these absorbers represent the same kinematic systems or not has been debated in the AGN community for some time now. The correct diagnostics of the inner-shell absorption lines could potentially provide a conclusive answer to this interesting astrophysical question. In order to demonstrate the large uncertainties of the currently available atomic data, in Table 1 we present three different calculations for the wavelengths of the strongest K$`\alpha `$ inner-shell absorption lines of O1+ through O5+ and also compare them with the deduced wavelengths from the Chandra/LETGS observation of NGC 5548, courtesy of J. Kaastra. The wavelengths from NGC 5548 were obtained assuming that the outflow velocities of all ions are similar to that of the well-calibrated O6+ velocity. Actually, this assumption may not be valid, but with the lack of laboratory measurements, it provides a rough idea of where to expect these lines.

Although one might have expected the R-matrix method to be the most rigorous, it is clear from the table that at the current state of the atomic data, it is virtually impossible to determine what are the correct rest frame positions of these lines. The discrepancies among the various methods reach 50 mÅ, which corresponds here to uncertainties in the measured velocities of 700 km/s. This uncertainty is of the same order as typical outflow velocities in nearby active galaxies, implying that these lines are practically useless for this purpose, despite their great potential. Laboratory measurements are desperately needed.

## 5 Suggested Z Pinch Measurements

The currently best available facility for producing inner-shell absorption lines and measuring their wavelengths and optical depths is the z pinch at Sandia National Laboratory. The powerful z pinch experiments (x-ray fluxes reaching 1019 erg s-1 cm-2) produce long-lived (6 ns), steady- state, photoionized gas (Bailey et al., 2001), which can be studied with high resolution spectrometers, both in emission and in absorption simultaneously. The ionizing spectrum can be characterized rather accurately by a blackbody spectrum. The control over the position and density of the targets in these experiments provides a sensitive handle on the ionization state and column density in the absorbing medium. A demonstration of these capabilities was given by Bailey et al. (2001), where absorption by photoionized Ne has been measured. The spectrum obtained with a crystal spectrometer in that experiment shows many individual lines that are nicely resolved, which allows for accurate wavelength and equivalent-width measurements.

## 6 Summary and Prospects

Many atomic data needs for hot, collisional plasmas relevant to x-ray spectroscopic observations, now regularly obtained with the gratings on board Chandra and XMM-Newton, have been provided by the ongoing LLNL EBIT programs of Kahn and Beiersdorfer. Particularly, these programs have provided the most important data for Fe-L emission by hot gas and work is in progress to measure many more high-quality data for other L-shell systems. Now that we actually have high-resolution cosmic spectra, we can determine better than before which atomic data are the most crucial. Thus, the most urgent needs of the x-ray astronomy community for collisional plasmas will continue to be addressed with the EBIT Laboratory Astrophysics Program. Photoionized gases have received less attention as they are more rare in nature and very few laboratory experiments have sufficient x-ray flux to produce them. Particularly missing are measurements for inner-shell absorption lines. The current uncertainties in the positions of these lines considerably limit our ability to use them for diagnostics. In the future, we intend to try to use z pinch experiments to remedy this deficiency.


text:
# 𝐾⁻-Nucleus Scattering at Low and Intermediate Energies.

## I Introduction

The works on meson baryon dynamics of Ref. KSW95  showed how Chiral symmetry constraints could be accommodated within a unitarity approach, able to describe resonances. This proved to be crucial to disentangle the intricate interaction between antikaons and nucleons at low energies RO97 Ga02 . The model of Ref. RO97 , was employed in Ref. RO00  to microscopically derive an optical potential for the $`K^{}`$ in nuclear matter in a self-consistent manner. Self-consistency turned out to be a crucial ingredient to derive the $`K^{}`$-nucleus potential and led to an optical potential considerably more shallow than those found in Refs. BF97 FG99 .

In Refs. zaki  and baca , the predictions of the chirally inspired potential of Ref. RO00  for measured shifts and widths of $`K^{}`$ atoms were evaluated, and it was found that this potential provides an acceptable description of the observed kaonic atom states, through the whole periodic table. Despite of having both real and imaginary parts of quite different depth, some other empirical optical potentials (FG93 FG99 ) also examined in Ref. baca , led to acceptable descriptions of the experimentally available $`K^{}`$ atom data as well. However, there were appreciable differences among the predicted widths for deeply bound antikaon nuclear states, not detected yet, when different potentials were used. The aim of this paper is to explore the possibility of differentiating between several $`K^{}`$ nucleus optical potentials by means of the scattering data. The extrapolation to finite $`K^{}`$ kinetic energies of the potential of Ref. RO00  requires at least the inclusion of the $`p`$-wave part of the $`K^{}`$ selfenergy. This was performed in Ref. pwave , and tested for $`K^{}p`$ scattering in Ref. Ji02 . However, even after having included $`p`$-wave contributions, one cannot expect reliable predictions from the theoretical potential of Refs. RO00  and pwave , at the lowest energy for which there exist experimental data ($`800`$ MeV for the $`K^{}`$ momentum), where $`d`$ and $`f`$ waves contributions are relevant. Besides, as we will show, for this relatively high energy, the impulse approximation works reasonably well, which is a clear indication that these data do not have much information on the details of the $`K^{}`$-nucleus dynamics. Thus, we have also focused our attention at the typical momentum of the $`K^{}`$ after the $`\varphi `$-meson decay ($``$ 127 MeV) with the hope that the scattering experience could be performed at DA$`\mathrm{\Phi }`$NE or at KEK or in the future Japanese Hadron Collider (JHC).

## II $`K^{}`$-Nucleus Optical Potentials

We solve the Klein Gordon equation

$$\left(\stackrel{}{}^2+\mu ^2+2\omega V_{\mathrm{opt}}\right)\mathrm{\Psi }=\left(\omega V_C\right)^2\mathrm{\Psi },$$
(1)
where $`\omega `$ is the Center of Mass (CM) $`K^{}`$–nucleus energy, $`V_C`$ and $`V_{\mathrm{opt}}`$ are the finite-size Coulomb and optical $`K^{}`$–nucleus potentials and $`\mu `$ the reduced $`K^{}`$-nucleus mass. At large distances and for a CM scattering angle $`\theta `$, the $`K^{}`$ wave-function $`\mathrm{\Psi }(\stackrel{}{r})`$ behaves as $`I(r)+f(\theta )S(r)`$, with $`I(r)`$ and $`S(r)`$ the standard wave functions for Coulomb scattering from a punctual charge $`Z`$ and $`f(\theta )`$ the scattering amplitude, which normalization is determined by its relation to the CM differential elastic cross section $`d\sigma _\mathrm{e}/d\mathrm{\Omega }=\left|f(\theta )\right|^2`$. The integrated cross sections read:

$`\sigma _\mathrm{e}`$ $`=`$ $`{\displaystyle \frac{\pi }{q^2}}{\displaystyle \underset{l}{}}(2l+1)\left|1\eta _le^{2\mathrm{i}\left(\sigma _l+\delta _l\right)}\right|^2,`$ (2)
$`\sigma _\mathrm{t}`$ $`=`$ $`{\displaystyle \frac{2\pi }{q^2}}{\displaystyle \underset{l}{}}(2l+1)\left[1\eta _l\mathrm{cos}\left(2\left(\sigma _l+\delta _l\right)\right)\right],`$ (3)
$`\sigma _{\mathrm{re}}`$ $`=`$ $`{\displaystyle \frac{\pi }{q^2}}{\displaystyle \underset{l}{}}(2l+1)\left[1\eta _l^2\right]`$ (4)
with $`q`$ the CM $`K^{}`$ momentum and $`\sigma _l`$, $`\delta _l`$ and $`\eta _l`$ the standard Coulomb phase shifts, the additional phase shifts due to strong interaction and the inelasticities appearing in the standard partial wave decomposition of $`f(\theta )`$ (see Ref. NOG93 ). While the elastic ($`\sigma _\mathrm{e}`$) and total ($`\sigma _\mathrm{t}`$) cross sections are infinite, the reaction ($`\sigma _{\mathrm{re}}`$) cross section is finite because of the short-range of the nuclear interaction.

The $`K^{}`$-nucleus optical potential, $`V_{\mathrm{opt}}`$, is related to the $`K^{}`$-selfenergy, $`\mathrm{\Pi }(q^0,|\stackrel{}{q}|),`$ inside of a nuclear medium, neglecting isovector effects, by

$$2\omega V_{\mathrm{opt}}(r)=\mathrm{\Pi }(m+T,(q_{}^{0}{}_{}{}^{2}m^2)^{\frac{1}{2}};\rho (r))$$
(5)
where $`m`$ and $`T`$ are the $`K^{}`$ mass and laboratory kinetic energy and $`\rho `$ is the sum of proton and neutron densities.

From the antikaon-selfenergy, as determined by, Refs. RO00  and pwave , we define the first selfenergy used in this work ($`\mathrm{\Pi }^{\mathrm{TH}}`$). This selfenergy does not have any free parameters, all the needed input is fixed either from studies of meson-baryon scattering in the vacuum or from previous studies of pion-nucleus dynamics NOG93 . It provides an acceptable ($`\chi ^2/dof`$ of 2.9) description of the set of 63 shifts and widths of $`K^{}`$ atom levels used in Ref. baca . We have neglected all type of non-localities, since they lead to changes in the results presented here of 3% at most.

As in Ref. baca , we also construct a modified selfenergy, which we call $`\mathrm{\Pi }^{\mathrm{THPH}}`$, by adding to $`\mathrm{\Pi }^{\mathrm{TH}}`$ a phenomenological part linear in density. This phenomenological part is determined by a constant $`\delta b_0`$ which we fix to the value $`(0.12\mathrm{i}0.38)\mathrm{fm}`$, obtained in Ref. baca  from a $`\chi ^2`$-fit to the kaonic atom data. The new selfenergy reads:

$`\mathrm{\Pi }^{\mathrm{THPH}}(r)`$ $`=`$ $`\mathrm{\Pi }^{\mathrm{TH}}(r)4\pi \delta b_0\rho (r)`$ (6)
The third selfenergy considered in this work is just obtained from the Impulse Approximation (IA), i.e., $`t\rho `$ form for the $`K^{}`$ selfenergy, and it neglects all orders higher than the leading one, in the density expansion. It reads:

$`\mathrm{\Pi }^{\mathrm{IA}}(r)`$ $`=`$ $`4\pi {\displaystyle \frac{\sqrt{s}}{M}}b_0^{\mathrm{IA}}(\sqrt{s})\rho (r)`$ (7)
with $`b_0^{\mathrm{IA}}(\sqrt{s})=\frac{1}{4}\left(3{}_{1}{}^{}f(\sqrt{s})+{}_{0}{}^{}f(\sqrt{s})\right)`$, $`M`$ the nucleon mass, $`\sqrt{s}`$ the total CM $`K^{}N`$ energy and $`{}_{I=1,0}{}^{}f`$ the isoscalar and isovector forward antikaon-nucleon scattering amplitudes, which partial wave decomposition reads:

$${}_{I}{}^{}f(\sqrt{s})=\underset{l}{}\left((l+1){}_{I}{}^{}f_{l}^{j_+}(\sqrt{s})+l{}_{I}{}^{}f_{l}^{j_{}}(\sqrt{s})\right)$$
(8)
with $`j_\pm =l\pm 1/2`$ the total angular momentum. At threshold, $`b_0^{\mathrm{IA}_{\mathrm{thr}}}b_0^{\mathrm{IA}}(m+M)=(0.15+\mathrm{i}0.62)\mathrm{fm}`$ Ma81 . The IA leads to extremely poor results for kaonic atoms BF97 ,baca . This is a clear indication that higher density corrections, not taken into account within the IA, are extremely important for kaonic atoms.

Finally, we have also considered two other antikaon selfenergies fitted to the kaonic atom data and energy independent (FG93  and FG99 ):

$`\mathrm{\Pi }^{2\mathrm{D}\mathrm{D}}`$ $`=`$ $`4\pi \left(1+{\displaystyle \frac{\mu }{M}}\right)\left(b_0^{\mathrm{IA}_{\mathrm{thr}}}+B_0\left({\displaystyle \frac{\rho }{\rho _0}}\right)^\alpha \right)\rho `$ (9)
$`\mathrm{\Pi }^{\mathrm{IAPH}}`$ $`=`$ $`4\pi \left(1+{\displaystyle \frac{\mu }{M}}\right)\stackrel{~}{b}_0\rho `$ (10)
with $`\stackrel{~}{b}_0=(0.52+\mathrm{i}0.80)\mathrm{fm}`$, $`B_0=(1.62\mathrm{i}0.028)\mathrm{fm}`$ and $`\alpha =0.273`$ as determined from $`\chi ^2`$-fits to $`K^{}`$-atom data in Ref. baca . Note that, though both $`\mathrm{\Pi }^{\mathrm{IAPH}}`$ and $`\mathrm{\Pi }^{\mathrm{IA}}`$ are linear in density selfenergies, they lead to substantially different potentials, since the real parts of the coefficients $`\stackrel{~}{b}_0`$ and $`b_0^{\mathrm{IA}_{\mathrm{thr}}}`$ differ both in sign and in size.

## III Results and Concluding Remarks

Since the $`K^{}`$ lifetime is relatively small, in practical terms it is experimentally difficult to count with low energetic $`K^{}`$ beams. However, all selfenergies described in the previous section, except for that obtained in the IA, are valid only near threshold. Thus, we have studied the case $`q_{\mathrm{lab}}=127`$ MeV, since this is the $`K^{}`$ momentum after the $`\varphi `$-meson decay.

In Fig. 1 we present results obtained with the $`K^{}`$ selfenergies fitted to the kaonic atom data. We also show results obtained by using the IA, where we have approximated the IA selfenergy at $`q_{\mathrm{lab}}=127`$ MeV by its threshold value quoted above. Strong interaction integrated elastic, reaction and total cross sections are also given in the top part of Table 1. We obtain these cross sections after having got rid of the Coulombian interaction, otherwise the total and elastic cross sections would diverge, i.e., we compute strong phase-shifts and inelasticities ($`\delta _l`$ and $`\eta _l`$) in presence of the Coulomb interaction, and afterwards we set to zero the Coulombian phase shifts, $`\sigma _l`$, in the formulae of Eqs. (2)-(4).

As can be seen in the figure, the non-linear density dependent $`K^{}`$ selfenergy, $`\mathrm{\Pi }^{2\mathrm{D}\mathrm{D}}`$, and the linear density dependent, threshold IA selfenergy, $`\mathrm{\Pi }(\sqrt{s})=\mathrm{\Pi }^{\mathrm{IA}}|_{\sqrt{s}=m+M}`$, provide extraordinarily similar results. Since both models have the same linear term in density, this is a clear indication that the reaction takes place in the surface of the nuclei, because of the big imaginary part of the potentials. The semiphenomenological $`\mathrm{\Pi }^{\mathrm{THPH}}`$ selfenergy, has a stronger departure from a linear behaviour in density than $`\mathrm{\Pi }^{2\mathrm{D}\mathrm{D}}`$, it has a smaller imaginary part (see Fig. 1 of Ref. baca ) and all of these explain the bigger differences with the IA model. Results, in particular position of the minima, obtained with $`\mathrm{\Pi }^{\mathrm{THPH}}`$ are clearly distinguishable from those obtained with any of the other three models also plotted in the figure, pointing out to a clear different density behavior likely due to the selfconsistent derivation of it. As a matter of example, for 12C in the region around $`\theta =60^\mathrm{o}`$, 2DD, IA and IAPH give similar elastic cross sections of about 33 mb/sr, whereas THPH gives about 52 mb/sr. This difference is appreciable and the size of the cross sections, tens of mb/sr, might allow to measure such a difference at DA$`\mathrm{\Phi }`$NE or KEK or in the future at the JHC. The differences are even bigger for larger angles, around the minimum of the THPH cross section (region 110-$`130^\mathrm{o}`$), but there, the cross sections are smaller, which makes harder to get the required statistics to see the effect. Besides, theoretical results in the neighborhood of a minimum are subject to more uncertainties. Similar conclusions can be drawn from the 40Ca and 208Pb results. In what respects to the integrated cross sections of Table 1, 2DD and IA give similar cross sections, though the IA reaction cross section is always slightly bigger, because the imaginary part of the $`B_0`$ parameter in Eq. (9) is negative. The IAPH model always provides the biggest reaction cross sections, because its selfenergy has also the largest imaginary part among all models considered (see Fig. 1 of Ref. baca ). Thus one can differentiate two sets of models, i.e, 2DD and IA selfenergies from THPH and IAPH ones. Besides, measurements, with precisions of about 10%, of the reaction cross sections would disentangle between THPH and IAPH models.

Let us look now to the experimental data. There only exist data Ma82  on $`K^{}`$ differential elastic cross sections for $`q_{\mathrm{lab}}=800`$ MeV and from 12C and 40Ca. In Fig. 2 we compare the IA predictions (solid line), including up to $`f`$waves (from Ref. Go77 ), to data. There also exist some data on total cross sections (mb), $`338\pm 8`$ \[$`306\pm 8`$\], Bu68  from 12C at $`q_{\mathrm{lab}}=800`$ \[$`655`$\] MeV. The IAs+p+d+f model provides again an acceptable description: 345 mb at 800 MeV and 304 mb at 655 MeV. Thus, this region is less sensitive to in nuclear medium effects than the $`K^{}`$ atom one. Indeed, the Glauber approximation also describes the 800 MeV scattering data, as discussed in SC99 . This work also corroborates that the imaginary part of the $`K^{}N`$ amplitude, obtained from the $`K^{}`$-nucleus scattering data, is close to that deduced in the vacuum. Besides, the contribution of $`d`$ and $`f`$ waves, not included in the THPH model, turn out to be important (compare the solid line to the IA results obtained when only $`s`$ and $`p`$ waves –dot-dashed line– are considered, in Fig. 2). In addition, the models of Refs. RO97  and pwave  for the $`s`$ and $`p`$ $`K^{}N`$ waves, though realistic near threshold, can not be safely extrapolated to momenta as high as 800 MeV. Thus, one expects the poor description of data provided by the THPH model. It is however surprising, that the IAPH predictions turn out to be almost indistinguishable from the IAs+p+d+f ones. This is merely a coincidence and, it occurs since accidentally at this momentum, $`b_0^{\mathrm{IA}}(\sqrt{s})`$ is approximately equal to $`\stackrel{~}{b}_0`$.

To finish, we also present results at an intermediate $`K^{}`$ momentum ($`q_{\mathrm{lab}}=300`$ MeV), despite the fact that there exist no data. For this momentum, calculations based on the IA shows that higher waves than the $`p`$ one have a small/moderate contribution and therefore can be neglected in some approximations. Thus in Fig. 3 and bottom part of Table 1, we compare again the THPH, 2DD and IAPH models for the $`K^{}`$ selfenergy inside the nuclear medium, together with the IA results including up to the $`p`$-wave, or up to the $`d`$-wave (partial waves are taken from Ref. Go77 ). The first observation is that the 2DD model differs now more than for the $`127`$ MeV momentum case, from the IA models. This is mainly due to the effect of $`p`$-wave in the latter ones. The second observation is that the semiphenomenological model THPH leads to a pattern clearly different than the rest of selfenergies, not only for the elastic differential cross section but also in the integrated ones compiled in Table 1. This is in principle good news, because then a scattering measurement in this region of $`K^{}`$ momentum will be definitive to disentangle between this approach and the others considered in Fig. 3 and bottom of Table 1. However a word of caution must be said here, the $`s`$-wave part of the antikaon selfenergy of Ref. RO00  is based on a model for the $`K^{}N`$ scattering in the free space that, though it is quite successful near threshold, predicts amplitudes for the isoscalar channel around $`q_{\mathrm{lab}}=300`$ MeV, with real parts which are in total disagreement (in sign and in size, see Ref. ORB02 ) with the analysis of Ref. Go77 . Thus, most probably one cannot trust the THPH model to describe the $`K^{}`$-dynamics at this momentum. Indeed, there is no reason either to believe more in the 2DD and IAPH models, and we believe that the more reliable predictions for $`q_{\mathrm{lab}}=300`$ MeV are those based on the IA.

## Acknowledgments

We thank to E.Oset and A.Ramos for useful discussions. A.J.M. acknowledges a fellowship from the University of Granada. This work was supported by Junta de Andalucía, DGI and FEDER funds (BFM2002-03218).