text:
# Element Diffusion During Cosmological Structure Formation

## 1 Introduction

Recent attempts to measure the primordial 4He abundance report accuracies better than $`1\%`$ (Izotov & Thuan 1998; Olive, Steigman, & Skillman 1997; Peimbert, Peimbert, & Ruiz 2000). Such accurate measurements are useful for testing the consistency of big bang nucleosynthesis and also for estimating the mean baryon density, $`\mathrm{\Omega }_b`$ (Schramm & Turner 1998). The primordial 4He mass fraction, $`\overline{Y}`$, is often inferred from the emission lines of blue compact galaxies (e.g., Peimbert & Torres-Peimbert 1974; Pagel, Terlevich, & Melnick 1986) under the implicit assumption that it had the same initial value throughout the Universe, including within underdense regions of the intergalactic medium (IGM).

However, a net difference between the element abundances in overdense and underdense regions could arise if diffusion occurs during the gravitational collapse of massive objects. Different atomic species, because of their different masses, will have different thermal speeds and will therefore experience different pressure forces per unit mass. Consequently, different elements would not flow into an overdense region at the same rate; instead, heavier elements would obtain higher relative abundances in galaxies than in voids. Diffusion is generally moderated by a frictional force arising from inter-species collisions. These collisions are most frequent for an ionized gas (Spitzer 1956) in which the characteristic cross–section for ion-proton collisions, $`\pi (Ze^2/k_BT)^2\mathrm{ln}\mathrm{\Lambda }_\mathrm{C}2\times 10^{12}Z(T/10^4\mathrm{K})^2\mathrm{cm}^2`$, is far greater than the collision cross–section for the corresponding atoms in a neutral gas, $`10^{15}\mathrm{cm}^2`$ (here, $`Ze`$ is the ion charge, $`T`$ is the temperature, $`k_B`$ is Boltzmann’s constant, and $`\mathrm{ln}\mathrm{\Lambda }_\mathrm{C}20`$ is the Coulomb Logarithm). We therefore focus our attention on the period between cosmological recombination and reionization during which the IGM was neutral (Loeb & Barkana 2001). During this stage, atomic collision cross sections may be small enough to allow some diffusion to occur.

In this paper we quantify the diffusion of 4He, 7Li, and D relative to 1H that occurs as a spherical, slightly overdense region begins to collapse. We assume a $`\mathrm{\Lambda }`$CDM cosmology, characterized by densities relative to critical of $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_m`$, and $`\mathrm{\Omega }_b=0.045`$ for matter, vacuum, and baryons; a Hubble constant of $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$; a 4He primordial mass fraction of $`\overline{Y}=0.244`$; and primordial 7Li and D number densities relative to 1H of $`\overline{X}_L=1.75\times 10^{10}`$ and $`\overline{X}_D=2.0\times 10^5`$, respectively. The effects of diffusion are most important where the gas pressure matters, namely around the comoving Jeans length scale $`\lambda _J=(2\pi /k_J)`$. For redshifts $`z140`$, when the baryonic temperature is primarily determined by the adiabatic expansion, the Jeans wavenumber is (Barkana & Loeb 2001)

$$k_J=\left[\frac{0.9\left(1+z\right)\mu m_H\mathrm{\Omega }_mH_0^2}{k_BT\left(z\right)}\right]^{1/2}$$
(1)
where $`h=(H_0/100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1)`$, $`\mu `$ is the mean molecular weight, $`m_H`$ is the mass of hydrogen and $`T\left(z\right)`$ the baryonic temperature at a redshift $`z`$.

In § 2, we derive the differential equations used for our analysis. These are linearized equations which are strictly valid only as long as the overdensities remain small. However, we assume that reionization occurs before the collapsing region becomes significantly nonlinear and we stop our calculation at reionization because charged particle cross sections suppress any further diffusion. In § 3, we discuss the initial conditions used for our differential equations. In § 4, we compute both the mean abundance deviation and the typical abundance fluctuation within a collapsing region. We then use a Press-Schechter mass function to average the resulting deviations over all collapsed objects in the Universe. In § 5, we compare our model to recent measurements of 4He abundances. Finally, we present our conclusions in § 6.

## 2 Basic Equations

We first establish our notation and illustrate our assumptions by carrying out a brief derivation of the equations satisfied by density perturbations. Each particle species has a mass density, $`\rho _s`$, and a bulk velocity, $`𝐯_s`$, which are related by the continuity equation,

$$\frac{\rho _s}{t}+_p\left(\rho _s𝐯_s\right)=0,$$
(2)
where we start by using physical (proper) coordinates. Since the Universe is initially homogeneous, $`\rho _s=\overline{\rho }_s\left(t\right)`$ and all particles move with the Hubble flow: $`𝐯_s=\overline{H}\left(t\right)𝐱`$ (Peebles 1993). In this approximation, equation (2) becomes

$$\frac{\overline{\rho }_s}{t}+3\overline{\rho }_s\overline{H}=0.$$
(3)
We next allow for density and velocity perturbations through the definitions

$$\rho _s(𝐱,t)=\overline{\rho }_s\left(t\right)\left[1+\delta _s(𝐱,t)\right],$$
(4)
$$^iv_s^j(𝐱,t)=\sigma _s^{ij}(𝐱,t)+\omega _s^{ij}(𝐱,t)+\left[\overline{H}\left(t\right)+\delta H_s(𝐱,t)\right]\delta _{kronecker}^{ij}.$$
(5)
The velocity gradient is decomposed into a symmetric, traceless tensor ($`\sigma _s^{ij}`$), an antisymmetric tensor ($`\omega _s^{ij}`$), and a trace. We insert these definitions into equation (2) and linearize the perturbations to obtain

$$\frac{\delta _s}{t}+\left(𝐯_s\right)\delta _s+3\delta H_s=0.$$
(6)
Switching from now on from physical to comoving coordinates, and letting an overdot represent a time derivative, equation (6) may be compactly expressed as

$$\dot{\delta _s}=3\delta H_s.$$
(7)
The density and bulk velocity also satisfy a force equation; in comoving coordinates,

$$\dot{𝐯}_s=\rho _s^1\underset{forces}{}𝐅_s.$$
(8)
The sum on the right-hand side includes all forces acting on species $`s`$. We linearize this equation and make use of equation (7) to obtain a new equation which does not explicitly contain the velocity perturbations,

$$\ddot{\delta _s}+2\overline{H}\dot{\delta }_s3\frac{\overline{H}}{t}3\overline{H}^2=\left[\frac{1}{\rho _s}\underset{forces}{}𝐅_s\right].$$
(9)
We now switch independent variables from time to redshift via the Friedmann equation, which leads to

$$H_0t=\frac{2}{3\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }}}\mathrm{sinh}^1\left(\sqrt{\frac{\mathrm{\Omega }_\mathrm{\Lambda }}{\mathrm{\Omega }_m}}\left(1+z\right)^{3/2}\right).$$
(10)
However $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ only becomes significant for $`z1`$; in this paper we consider redshifts $`z1`$ for which

$$H_0t\frac{2}{3\sqrt{\mathrm{\Omega }_m}}\left(1+z\right)^{3/2}.$$
(11)
Denoting a derivative with respect to $`z`$ by a prime, equation (9) becomes

$$\delta _s^{\prime \prime }+\frac{1}{2\left(1+z\right)}\delta _s^{}\frac{3}{\mathrm{\Omega }_mH_0^2\left(1+z\right)^5}\left(\frac{\overline{H}}{t}+\overline{H}^2\right)=\frac{1}{\mathrm{\Omega }_mH_0^2\left(1+z\right)^5}\left[\frac{1}{\rho _s}\underset{forces}{}𝐅_s\right].$$
(12)
We now obtain expressions for the forces. The force due to gravity, when expressed in terms of the potential, $`\varphi `$, is $`𝐅_s=\rho _s\varphi `$ . We insert this into equation (12) and use Poisson’s equation which gives

$$\delta _s^{\prime \prime }+\frac{1}{2\left(1+z\right)}\delta _s^{}=\frac{1}{\mathrm{\Omega }_mH_0^2\left(1+z\right)^5}\left[\frac{1}{\rho _s}_{forces}^{}𝐅_s\right]+\frac{3}{2\mathrm{\Omega }_m}\left(1+z\right)^2\underset{species}{}\mathrm{\Omega }_i\delta _i,$$
(13)
where the primed sum means a sum over all forces except gravity.

Although all types of particles are affected by gravity, we assume that only the baryons ever obtain a non-zero temperature. Each species of baryons experiences a different pressure force which is determined by the negative of its pressure gradient,

$$𝐅_s=P_s=\frac{k_B}{m_s}\left(\rho _sT_s+T_s\rho _s\right),$$
(14)
where we have used $`m_s`$ for the atomic mass and $`T_s`$ for the temperature of species $`s`$. Dividing equation (14) by the mass density and taking the divergence, we see that

$$\frac{𝐅_𝐬}{\rho _s}=\frac{k_B}{m_s}\left(^2T_s+T_s^2\delta _s\right)$$
(15)
to first order in the perturbations. We define $`\overline{T}_s`$ to be the temperature in the absence of perturbations, and express the effects of perturbations in terms of the small quantity $`\tau _s`$: $`T_s=\overline{T}_s\left(1+\tau _s\right)`$. Inserting this into equation (15) and linearizing once more we see that

$$\frac{𝐅_𝐬}{\rho _s}=\frac{k_B\overline{T}_s}{m_s}\left(^2\tau _s+^2\delta _s\right).$$
(16)
The temperature of the baryons is determined by two effects, the coupling of free electrons to the background radiation and the adiabatic expansion (Peebles 1993). This gives

$$\frac{dT_s}{dt}=\frac{x}{1+x}\frac{8\sigma _TaT_\gamma ^4}{3m_ec}\left(T_\gamma T_s\right)+\frac{2T_s}{3\rho _s}\frac{d\rho _s}{dt},$$
(17)
where $`x`$ is the free electron fraction, $`\sigma _T`$ is the Thomson cross section, $`a`$ is the Stefan-Boltzmann constant, $`m_e`$ is the electron mass, $`c`$ is the speed of light, and $`T_\gamma `$ is the temperature of the cosmic microwave background (CMB). Switching the independent variable to $`z`$ and linearizing, equation (17) becomes

$$\frac{dT_s}{dz}=\frac{x}{1+x}\frac{8\sigma _TaT_0^5\left(1+z\right)^{5/2}}{3m_ecH_0\sqrt{\mathrm{\Omega }_m}}\left(1\frac{T_s}{T_0\left(1+z\right)}\right)+\frac{2T_s}{\left(1+z\right)}+\frac{2}{3}T_s\delta _s^{},$$
(18)
where $`T_0=2.73`$ K is the CMB temperature at present.

The mean baryon temperature, $`\overline{T}_s`$, is defined to be a solution of

$$\frac{d\overline{T_s}}{dz}=\frac{x}{1+x}\frac{8\sigma _TaT_0^5\left(1+z\right)^{5/2}}{3m_ecH_0\sqrt{\mathrm{\Omega }_m}}\left(1\frac{\overline{T}_s}{T_0\left(1+z\right)}\right)+\frac{2\overline{T}_s}{\left(1+z\right)}.$$
(19)
This indicates that $`\overline{T}_s`$ is independent of baryonic species and so all baryons would be at the same temperature in the absence of perturbations. We use the computer program RECFAST (Seager, Sasselov, & Scott 1999) to determine the free-electron fraction as a function of redshift. Then we use equation (19) to compute the mean baryon temperature. Both quantities are shown in Figure Element Diffusion During Cosmological Structure Formation.

After linearizing equation (18), we obtain the equations for the temperature perturbations,

$$w_s^{}\frac{x}{1+x}\frac{8\sigma _TaT_0^5\left(1+z\right)^{5/2}}{3m_ecH_0\sqrt{\mathrm{\Omega }_m}}\frac{w_s}{\overline{T}_s}=\frac{x}{1+x}\frac{16\sigma _TaT_0^5\left(1+z\right)^{5/2}}{9m_ecH_0\sqrt{\mathrm{\Omega }_m}}\frac{\delta _s}{\overline{T}_s},$$
(20)
where $`w_s=\tau _s\frac{2}{3}\delta _s`$. With this definition of $`w_s`$, equation (16) becomes

$$\frac{𝐅_𝐬}{\rho _s}=\frac{k_B\overline{T}}{m_s}\left(\frac{5}{3}^2\delta _s+^2w_s\right).$$
(21)
Because of the spatial derivatives, we are forced to work in comoving Fourier space. However, for simplicity of notation we will continue to simply write $`\delta _s`$ and $`w_s`$, although from now on we use these symbols to represent the Fourier components of the perturbations.

Finally we take into account the frictional term arising from inter-species collisions. Physically, this force is non-zero only if two species have different bulk velocities; thus, the force on one species due to a second species may be written as $`𝐅_1=\kappa (𝐱,t)\left(𝐯_2𝐯_1\right)`$ (Burgers 1969). Since the difference in bulk velocities is already linear in the perturbations, we only need to calculate the coefficient $`\kappa `$ to zeroth order for the linearized equations of motion. This allows us to calculate $`\kappa `$ under the assumption that all baryons share the same temperature. In general, the frictional force per unit volume acting on species $`i`$ due to species $`j`$ is obtained by multiplying the number of collisions per unit time per unit volume by the momentum transfer per collision. We multiply this by the probability that species $`i`$ and $`j`$ have initial momentum $`𝐩_i`$ and $`𝐩_j`$ and that the center of mass scattering angle is $`\chi `$. Finally, we integrate over all parameters,

$$𝐅_{ij}=d^3𝐩_id^3𝐩_j𝑑\mathrm{\Omega }n_in_j\sigma _{ij}\left|𝐯_i𝐯_j\right|\mathrm{\Delta }𝐩_i(𝐩_i,𝐩_j,\chi )f_i\left(𝐩_i\right)f_j\left(𝐩_j\right)\frac{1}{\sigma _{ij}}\frac{d\sigma _{ij}}{d\mathrm{\Omega }}(𝐩_i,𝐩_j,\chi ).$$
(22)
We assume that the collisions are elastic and that the atoms can be approximated by hard spheres. We take each species to have a Maxwellian velocity distribution. Following the integration equation (22) becomes

$$𝐅_{ij}=\frac{16m_{ij}}{3\sqrt{2\pi }}\sqrt{\frac{k_B\overline{T}}{m_{ij}}}n_in_j\sigma _{cl}\left(𝐯_i𝐯_j\right)_0^{\mathrm{}}𝑑xx^5e^{x^2}Q_d\left(x\sqrt{\frac{k_B\overline{T}}{m_{ij}}}\right).$$
(23)
Here, $`m_{ij}=m_im_j/(m_i+m_j)`$ is the reduced mass, $`\sigma _{cl}`$ is the classical hard-sphere cross section, $`n_i=\rho _i/m_i`$, and the integral is a quantum mechanical correction (Massey, Burhop, & Gilbody 1969). The function $`Q_d\left(q\right)`$ is given by

$$Q_d=\frac{1}{q^2}\underset{l=0}{\overset{\mathrm{}}{}}\left[4\left(2l+1\right)\mathrm{sin}^2\eta _l2\left(l+1\right)\left(1\mathrm{cos}2\eta _l\mathrm{cos}2\eta _{l+1}+\mathrm{cos}2\left(\eta _l\eta _{l+1}\right)\right)\right].$$
(24)
The integration variable $`x`$ is related to $`q`$ by $`q=xa\sqrt{2m_{ij}k_B\overline{T}/\mathrm{}^2}`$ , and the phase shifts are given by $`\eta _0=q`$ and $`\mathrm{tan}\eta _i=\left(1\right)^{i1}J_{i+1/2}\left(q\right)/J_{i1/2}\left(q\right)`$ , where the functions $`J_n\left(q\right)`$ are Bessel functions of order $`n`$. The distance scale $`a`$ represents the sum of the radii of two colliding atoms; we adopt the values $`r_H=r_D=0.53`$Å , $`r_{He}=0.31`$Å , and $`r_{Li}=1.45`$Å (Clementi, Raimondi, & Reinhardt 1963; Slater 1964). We evaluate the integral in equation (23) numerically and plot the results in Figure Element Diffusion During Cosmological Structure Formation for each type of inter-species collision. Quantum mechanics becomes more important as the temperature decreases but the correction is never greater than a factor of four. For brevity, we henceforth denote this integral by $`I_{ij}\left(z\right)`$.

We now summarize the equations governing each species. For dark matter equation (12) becomes

$$\delta _{dm}^{\prime \prime }+\frac{1}{2}\left(1+z\right)^1\delta _{dm}^{}=\frac{3}{2}\left(1+z\right)^2\underset{species}{}\frac{\mathrm{\Omega }_i}{\mathrm{\Omega }_m}\delta _i.$$
(25)
For each baryonic species there are also pressure and collisional friction terms, and so

$`\delta _i^{\prime \prime }+{\displaystyle \frac{1}{2}}\left(1+z\right)^1\delta _i^{}`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left(1+z\right)^2{\displaystyle \underset{j=allspecies}{}}{\displaystyle \frac{\mathrm{\Omega }_j}{\mathrm{\Omega }_m}}\delta _j`$
$`{\displaystyle \frac{k_B\overline{T}}{m_i\mathrm{\Omega }_mH_0^2}}\left(1+z\right)^3k^2\left({\displaystyle \frac{5}{3}}\delta _i+w_i\right)`$
$`+{\displaystyle \frac{H_0}{G}}\sqrt{{\displaystyle \frac{2k_B\overline{T}}{\mathrm{\Omega }_m\pi ^3}}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{\mathrm{\Omega }_j\sigma _{cl}^{ij}}{\sqrt{m_im_j\left(m_i+m_j\right)}}}\left(1+z\right)^{1/2}I_{ij}\left(z\right)\left(\delta _i^{}\delta _j^{}\right),`$
where $`k`$ is the comoving wave number. From equation (2) we can easily find the differential equations satisfied by the differences relative to hydrogen, $`\mathrm{\Delta }_i=\delta _i\delta _H`$:

$`\mathrm{\Delta }_i^{\prime \prime }+{\displaystyle \frac{1}{2}}\left(1+z\right)^1\mathrm{\Delta }_i^{}`$ $`=`$ $`{\displaystyle \frac{k_B\overline{T}}{\mathrm{\Omega }_mH_0^2}}\left(1+z\right)^3k^2\left[{\displaystyle \frac{5}{3}}\left({\displaystyle \frac{\delta _H}{m_H}}+{\displaystyle \frac{\mathrm{\Delta }_i\delta _H}{m_i}}\right)+\left({\displaystyle \frac{w_H}{m_H}}{\displaystyle \frac{w_i}{m_i}}\right)\right]`$
$`{\displaystyle \frac{H_0}{G}}\sqrt{{\displaystyle \frac{2k_B\overline{T}}{\mathrm{\Omega }_m\pi ^3}}}{\displaystyle \underset{jH}{}}{\displaystyle \frac{\mathrm{\Omega }_j\sigma _{cl}^{Hj}}{\sqrt{m_Hm_j\left(m_H+m_j\right)}}}\left(1+z\right)^{1/2}I_{Hj}\left(z\right)\mathrm{\Delta }_j^{}`$
$`+{\displaystyle \frac{H_0}{G}}\sqrt{{\displaystyle \frac{2k_B\overline{T}}{\mathrm{\Omega }_m\pi ^3}}}{\displaystyle \underset{ji}{}}{\displaystyle \frac{\mathrm{\Omega }_j\sigma _{cl}^{ij}}{\sqrt{m_im_j\left(m_i+m_j\right)}}}\left(1+z\right)^{1/2}I_{ij}\left(z\right)\left(\mathrm{\Delta }_j^{}\mathrm{\Delta }_i^{}\right).`$
## 3 Initial Conditions

We consider dark matter perturbations initially of the spherical top-hat form $`\delta _{dm}(𝐱,z_{ini})=\delta _{ini}\theta \left(R\left|𝐱\right|\right)`$, which has the Fourier components

$$\delta _{dm}(𝐤,z_{ini})=\frac{4\pi \delta _{ini}}{k^3}\left(\mathrm{sin}kRkR\mathrm{cos}kR\right).$$
(28)
Here, $`\delta _{ini}`$ and $`R`$ are constants and $`\theta `$ represents the step function. We take the initial values of $`\mathrm{\Delta }_{He}`$, $`\mathrm{\Delta }_{Li}`$, $`\mathrm{\Delta }_D`$ and their derivatives to be zero, but we would like to choose the initial values of $`\delta _{dm}^{}`$, $`\delta _H`$, and $`\delta _H^{}`$ that correspond to the most rapidly growing mode. To gain some insight into how to make this choice, we consider a simplified problem where baryonic diffusion is negligible. We approximate the baryonic temperature by the temperature of the background radiation (valid for $`z200`$; see Figure Element Diffusion During Cosmological Structure Formation). Equation (2) for the baryon overdensity becomes

$`\delta _b^{\prime \prime }+{\displaystyle \frac{1}{2}}\left(1+z\right)^1\delta _b^{}`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left(1+z\right)^2\left({\displaystyle \frac{\mathrm{\Omega }_{dm}}{\mathrm{\Omega }_m}}\delta _{dm}+{\displaystyle \frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_m}}\delta _b\right)`$
$`{\displaystyle \frac{k_BT_0}{\mu \mathrm{\Omega }_mH_0^2}}\left(1+z\right)^2k^2\delta _b.`$
We use this equation to express $`\delta _{dm}`$, differentiate, and insert the results into equation (25). The outcome is a fourth order differential equation for $`\delta _b`$. It has four independent, power law solutions which we write as $`\delta _b=_{i=1}^4A_i\left(1+z\right)^{p_i}`$. The four $`A_i`$ are the four arbitrary constants associated with a linear fourth order equation. The four $`p_i`$ are given by

$$p_i=\{\begin{array}{c}\frac{1}{4}+\sqrt{\frac{13}{16}\frac{\widehat{k}^2}{2}+\frac{1}{2}\sqrt{\left(\widehat{k}^2\frac{3}{2}\right)^2+6\widehat{k}^2\frac{\mathrm{\Omega }_{dm}}{\mathrm{\Omega }_m}}}\hfill \\ \frac{1}{4}+\sqrt{\frac{13}{16}\frac{\widehat{k}^2}{2}\frac{1}{2}\sqrt{\left(\widehat{k}^2\frac{3}{2}\right)^2+6\widehat{k}^2\frac{\mathrm{\Omega }_{dm}}{\mathrm{\Omega }_m}}}\hfill \\ \frac{1}{4}\sqrt{\frac{13}{16}\frac{\widehat{k}^2}{2}+\frac{1}{2}\sqrt{\left(\widehat{k}^2\frac{3}{2}\right)^2+6\widehat{k}^2\frac{\mathrm{\Omega }_{dm}}{\mathrm{\Omega }_m}}}\hfill \\ \frac{1}{4}\sqrt{\frac{13}{16}\frac{\widehat{k}^2}{2}\frac{1}{2}\sqrt{\left(\widehat{k}^2\frac{3}{2}\right)^2+6\widehat{k}^2\frac{\mathrm{\Omega }_{dm}}{\mathrm{\Omega }_m}}},\hfill \end{array}$$
(30)
where $`\widehat{k}^2`$ is the squared comoving wavenumber normalized by $`(\mu \mathrm{\Omega }_mH_0^2/k_BT_0)`$. The resulting expression for the dark matter overdensity is

$$\delta _{dm}=\underset{i=1}{\overset{4}{}}c_iA_i\left(1+z\right)^{p_i},$$
(31)
with the constants $`c_i`$ given by

$$c_i=p_i^2\left(\frac{2\mathrm{\Omega }_m}{3\mathrm{\Omega }_{dm}}\right)p_i\left(\frac{\mathrm{\Omega }_m}{3\mathrm{\Omega }_{dm}}\right)+\frac{2\widehat{k}^2\mathrm{\Omega }_m3\mathrm{\Omega }_b}{3\mathrm{\Omega }_{dm}}.$$
(32)
Motivated by this solution, we choose our initial conditions so that the perturbation will grow like the most negative $`p_i`$ (i.e., the third one listed in equation 30). We let the three $`A_i`$ corresponding to all but the largest growth rate be zero, and let the final $`A_i`$ be determined by the tophat distribution (equation 28). The remaining initial conditions may now be written in terms of the known $`c_i`$ and $`p_i`$.

## 4 Calculations

### 4.1 Evolution of the k-Modes

We define a new set of variables,

$$\widehat{\delta }_s=\frac{k^3}{4\pi \delta _{ini}\left(\mathrm{sin}kRkR\mathrm{cos}kR\right)}\delta _s,$$
(33)
$$\widehat{\mathrm{\Delta }}_s=\frac{k^3}{4\pi \delta _{ini}\left(\mathrm{sin}kRkR\mathrm{cos}kR\right)}\left(\delta _s\delta _H\right),$$
(34)
$$\widehat{w}_s=\frac{k^3}{4\pi \delta _{ini}\left(\mathrm{sin}kRkR\mathrm{cos}kR\right)}w_s.$$
(35)
These quantities are special because they not only satisfy the same differential equations as the $`\delta `$ and $`w`$, but also satisfy simpler initial conditions, $`\widehat{\delta }_{dm}=1`$ , $`\widehat{\delta }_{dm}^{}=\frac{p_3}{1+z_{ini}}`$ , $`\widehat{\delta }_b=c_3^1`$ , and $`\widehat{\delta }_b^{}=\frac{p_3}{c_3\left(1+z_{ini}\right)}`$ . The parameters $`(\delta _{ini},R)`$ have disappeared from the problem. Therefore we only evolve each normalized $`k`$-mode once; we obtain the physical overdensity for a given $`\delta _{ini}`$ and $`R`$ by multiplying each $`k`$-mode by the appropriate factor from equation (33).

We have set $`z_{ini}=10^3`$ following cosmological recombination. The final redshift is marked by the first of two possible events. Our equations are applicable in the linear regime, so we must stop if the object becomes nonlinear. At a transition to nonlinearity, the density would be greatly enhanced which would result in the suppression of diffusion. Therefore, we expect diffusion to be effective in the linear regime where our equations are valid. The second stopping condition relates to the reionization of the Universe. At that time, large Coulomb cross sections between charged particles are able to effectively halt diffusion. In conclusion, we consider the case where the object becomes nonlinear after reionization and adopt $`z_{reion}=8`$ (Haiman & Loeb 1998, 1999a,b; Gnedin & Ostriker 1997; Chiu & Ostriker 2000; Gnedin 2000) as the stopping point of our diffusion calculation. Since the collapsed fraction of the Universe is small at $`z=8`$ (Barkana & Loeb 2001), most of the material within galaxies that form much later (at $`z3`$) was still in the linear regime at $`z8`$, where our diffusion calculation is adequate.

Several results of having carried out this prescription are shown in Figure Element Diffusion During Cosmological Structure Formation. All length scales have been measured in terms of $`k_J^1`$, defined in equation (1). This expression involves the baryonic temperature which is set by the adiabatic expansion for redshifts less than $`z_t138`$, $`T\left(z\right)=T_0\left(1+z\right)^2\left(1+z_t\right)^1`$; as discussed in § 2. At $`z_{reion}`$, the Jeans wavenumber therefore obtains the value

$$k_J=1040\mathrm{Mpc}^1\left[\left(\frac{9}{1+z}\right)\left(\frac{\mathrm{\Omega }_mh^2}{0.147}\right)\left(\frac{2.73\mathrm{K}}{T_0}\right)\right]^{1/2},$$
(36)
which is the value used in Figure Element Diffusion During Cosmological Structure Formation. We see that the dark matter overdensity has evolved approximately proportionally to $`\left(1+z\right)^1`$, as might have been expected since gravity is the only force applicable to dark matter. On large scales corresponding to $`kk_J`$, the baryons follow the dark matter; however, on small scales corresponding to $`kk_J`$, pressure forces prevent baryonic modes from collapsing. This leads to the slight decline in dark matter amplification near the Jeans scale and the dramatic decline in the amplification of baryonic modes there. In the lower panel of Figure 3, we see that the perturbations of the different baryonic species do not experience exactly the same amplification. At $`kk_J`$, pressure is negligible for all species, and the modes of each are amplified. Since hydrogen is the lightest species, it reaches its Jeans scale at the lowest $`k`$. Thus, all of the $`\widehat{\mathrm{\Delta }}_x`$ become positive at approximately the same wavenumber. As $`k`$ increases, we begin to reach the Jeans scales of the heavier species. Then the growth of $`k`$-modes of this second species is also suppressed, and the $`\widehat{\mathrm{\Delta }}_x`$ approach zero as $`k`$ gets large. However, Figure Element Diffusion During Cosmological Structure Formation indicates that the overall magnitude of diffusion is not solely determined by the Jeans scale of a species, or else the effect would have been largest for 7Li. This is because of the frictional term which serves to suppress diffusion, which is proportional to the square of the radius of an atom and therefore the greatest suppression is associated with the largest atoms.

The normalized overdensities in Figure Element Diffusion During Cosmological Structure Formation were obtained by evolving 1025 $`k`$-modes, equally spaced between with $`0k/k_J100`$. To obtain the physical overdensities, each mode must be multiplied by its normalization factor. The normalization factor depends on the object assembly radius $`R`$, and is fairly rapidly varying compared to the normalized overdensities. To get an accurate determination of the physical overdensity, we interpolate the normalized overdensity so that we sample it at a rate larger than the frequency of the normalization factor. We then use equation (33) to obtain the physical $`k`$-modes. One such sampling is shown in Figure Element Diffusion During Cosmological Structure Formation. The value of $`R`$ used for this plot corresponds to an enclosed mass of $`M=10^9M_{}=\rho _{crit}\mathrm{\Omega }_m\frac{4}{3}\pi R^3`$. The relative errors incurred in interpolation are typically $`10^6`$, and seldom exceed $`10^5`$ with the 1025 initial points. The normalization factor is also proportional to the initial dark matter perturbation, $`\delta _{ini}`$, which we express in terms of the collapse redshift of the object ($`z_{coll}`$) through the approximate relation $`\delta _{ini}\left(1+z_{ini}\right)=\delta _{coll}\left(1+z_{coll}\right)`$ . We take $`\delta _{coll}=1.69`$ for the overdensity of collapse and, unless otherwise stated, $`z_{coll}=2.5`$ for the redshift of collapse. Such objects would be just on the fringe of nonlinearity at $`z_{reion}=8`$. Since diffusion within all objects was shut off at $`z_{reion}`$, the effects of diffusion would be decreased for objects that collapsed later.

### 4.2 Predicted Statistics

#### 4.2.1 Mean Abundances

Next we quantify the extent to which element abundances in collapsed objects differ from their primordial values. Abundances of 4He are commonly reported in terms of the mass fraction,

$$Y=\frac{\overline{\rho }_{He}\left(1+\delta _{He}\right)}{\overline{\rho }_{He}\left(1+\delta _{He}\right)+\overline{\rho }_H\left(1+\delta _H\right)},$$
(37)
where we have neglected the contributions of all species but 4He and 1H to the mass of the Universe. After linearizing this equation, we see that

$$\frac{Y\overline{Y}}{\overline{Y}}=\left(1\overline{Y}\right)\left(\delta _{He}\delta _H\right).$$
(38)
In contrast, the abundances of 7Li and D are typically reported in terms of their number densities relative to hydrogen. In this case,

$$\frac{X_{L,D}\overline{X}_{L,D}}{\overline{X}_{L,D}}=\left(\delta _{Li,D}\delta _H\right).$$
(39)
To see how the element abundances in collapsed objects differ, on average, from their primordial values, we average equations (38)- (39) over the original sphere of comoving radius $`R`$,

$$\mathrm{\Delta }_R=\frac{3}{4\pi R^3}_{sphere}d^3𝐱\frac{d^3𝐤}{\left(2\pi \right)^3}\mathrm{\Delta }_𝐤e^{i𝐤𝐱}.$$
(40)
After integrating, we see that

$$\mathrm{\Delta }_R=\frac{3}{2\pi ^2R^3}_0^{\mathrm{}}\frac{dk}{k}\mathrm{\Delta }_k\left(\mathrm{sin}kRkR\mathrm{cos}kR\right).$$
(41)
We carry out this final integration numerically. As described in § 4.1, we evolve 1025 $`k`$-modes between 0 and 100$`k_J`$ to get the $`\widehat{\mathrm{\Delta }}_x`$ (where $`x`$ could denote any of 4He, 7Li, D). We interpolate to get values for $`\widehat{\mathrm{\Delta }}_x`$ at a fixed number of points, and then multiply each point by its normalization. We then evaluate the integral in equation (41) with these points. Then we repeat the procedure, but with twice as many interpolated points. If the two results differ by less than a fraction of $`5\times 10^3`$, we stop and accept the result.

Figure Element Diffusion During Cosmological Structure Formation illustrates the dependence of these deviations on the size of the collapsed object. Instead of measuring the size of the object by its assembly radius, $`R`$, we instead measure it in terms of its mass defined by

$$M=\rho _{crit}\mathrm{\Omega }_m\frac{4}{3}\pi R^3.$$
(42)
We see that the deviation from the primordial mean shrinks as the collapsed object increases in size approximately according to a power law with slope $`0.33`$. A simple argument yields the same result: assuming that the collapsing region starts out with primordial abundances throughout its volume and that diffusion occurs through the boundary of the region, the deviations from primordial abundances should fall off like $`1/R`$, or $`M^{1/3}`$.

Now that element abundances for particular collapsed objects are known, we proceed by averaging these abundances over all objects in the Universe that have collapsed by a particular redshift. This statistic is approximated by using a Press-Schechter mass function; for 4He,

$$Y=\frac{_{M_{min}}^{\mathrm{}}Y_MM\frac{dn}{dM}𝑑M}{_{M_{min}}^{\mathrm{}}M\frac{dn}{dM}𝑑M},$$
(43)
where $`dn/dM`$ is the number density of objects of mass $`M`$ that have collapsed at a given redshift and $`Y_M`$ is the averaged 4He mass fraction within an object of mass $`M`$. Equation (43) is also appropriate for the 7Li and D abundances if we replace $`Y`$ with $`X_{L,D}`$. Note also that there is a minimum mass scale of galaxies in which star formation is possible that we denote by $`M_{min}`$ (Abel 1995; Tegmark et al. 1997). We use a code developed by Eisenstein & Hu (1999) and Barkana (2000) to calculate $`dn/dM`$ at a redshift of $`2.5`$. The averaged abundances differ from primordial by the amounts shown in Figure Element Diffusion During Cosmological Structure Formation. We also consider a collapse redshift of $`z_{coll}=0`$; in this case Figure Element Diffusion During Cosmological Structure Formation shows that the abundance deviations decrease somewhat. This is not only due to the resulting decrease in $`\delta _{ini}`$, but also due to the dependence of $`dn/dM`$ on collapse redshift. At higher $`z`$, $`dn/dM`$ tends to favor lower mass scales where the abundance deviations are larger (Figure Element Diffusion During Cosmological Structure Formation).

#### 4.2.2 Fluctuations of Abundances Within Objects

The element abundances of collapsed objects can differ from primordial abundances in two ways. First, the abundances of an object could be uniformly different from the primordial values. This possibility was already discussed in the previous section. Second, if the element abundances within an object vary greatly from point to point, differences from the primordial abundances may be observed because of the fluctuations. We consider this latter effect by computing the standard deviation of the abundances within a given object. The lowest-order formula for this quantity is

$$\sigma _Y^2=\overline{Y}^2\left(1\overline{Y}\right)^2\sigma _{\mathrm{\Delta }_Y}^2$$
(44)
for 4He and

$$\sigma _X^2=\overline{X}^2\sigma _{\mathrm{\Delta }_X}^2$$
(45)
for 7Li and D. We have already calculated the mean values of the abundances, so now we only need the mean squares. We proceed as above, but this time we obtain a double integral that must be evaluated numerically,

$$\left(\mathrm{\Delta }^2\right)_M=\frac{3}{8\pi ^4R^3}_0^{\mathrm{}}𝑑kk\mathrm{\Delta }_k_0^{\mathrm{}}𝑑ll\mathrm{\Delta }_l\left(\frac{\mathrm{sin}\left(kl\right)R}{kl}\frac{\mathrm{sin}\left(k+l\right)R}{k+l}\right).$$
(46)
The integration region is the first quadrant of the $`kl`$ plane. Since the integrand is symmetric with respect to the interchange of $`k`$ and $`l`$, it is equal to the integral

$$\left(\mathrm{\Delta }^2\right)_M=\frac{3}{4\pi ^4R^3}_0^{\mathrm{}}𝑑kk\mathrm{\Delta }_k_0^k𝑑ll\mathrm{\Delta }_l\left(\frac{\mathrm{sin}\left(kl\right)R}{kl}\frac{\mathrm{sin}\left(k+l\right)R}{k+l}\right).$$
(47)
We carry out the integrations as with equation (41) and our results are plotted as a function of mass scale in Figure Element Diffusion During Cosmological Structure Formation. Comparing Figures Element Diffusion During Cosmological Structure Formation and Element Diffusion During Cosmological Structure Formation, we see that the relative importance of abundance fluctuations depends on the mass of the object in question. On small mass scales, abundances differ from their primordial values by an approximately uniform amount; but as the object size increases fluctuations from the mean are more likely to produce a difference from the primordial values. We also average the standard deviation over all collapsed objects in the Universe. Our prescription for this is identical to that given by equation (43), with $`\sigma _M`$ taking the place of $`Y_M`$. These results are plotted versus the minimum mass scale for galaxies in Figure Element Diffusion During Cosmological Structure Formation. Comparing this to Figure Element Diffusion During Cosmological Structure Formation, we again see that the relative importance of fluctuations increases with increasing mass scale.

## 5 Observational Data

Though much recent work has focused on uncovering the primordial 7Li and D abundances (e.g., Burles & Tytler 1998a,b; Tytler et al. 1999; Webb et al. 1997; Ryan et al. 2000), the 4He abundance with a reported uncertainty of about $`1\%`$, is still more accurately determined. However, even this accuracy is at least an order of magnitude coarser than the mean effect of diffusion for low-mass objects. But, as discussed in § 4 the 4He abundance is also expected to vary from point to point within an object. In this section, we re-examine the 4He abundance data presented in Izotov & Thuan (1998) to determine whether or not such fluctuations are currently detectable.

In Izotov & Thuan (1998), the helium abundance is assumed to be linear in metallicity for low metallicity objects. Each data point consists of a measurement of the metallicity, $`Z_i`$, a measurement of the helium abundance $`Y_i`$, and estimations for the uncertainties in each, $`ϵ_i`$ and $`\delta _i`$. A straight line is fit to the data, the intercept of which reveals the primordial (zero-metallicity) helium abundance. The model we use is slightly different. Because of diffusion we cannot directly associate 4He abundance measurements of zero-metallicity objects with the primordial 4He abundance. Not only do we expect an object of a particular mass to have a 4He abundance differing from the primordial value as shown in Figure Element Diffusion During Cosmological Structure Formation, but also we expect that different measurements of the single object would display an intrinsic scatter about this mean as shown in Figure Element Diffusion During Cosmological Structure Formation. Therefore, we fit the data to the following model. If an object has metallicity $`Z`$, the probability that it will have a 4He mass fraction $`Y`$ is denoted $`g\left(Y|Z\right)dY`$. We account for measurement uncertainties by defining $`h\left(Y_i|Y\right)`$ to be the probability of measuring $`Y_i`$ if the object actually has abundance $`Y`$, and $`f\left(Z|Z_i\right)`$ to be the probability that an object has metallicity $`Z`$ if $`Z_i`$ is what is observed. Thus, the total probability of a measurement simultaneously obtaining $`Y_i`$ and $`Z_i`$ is

$$P(Y_i|Z_i)=f\left(Z|Z_i\right)g\left(Y|Z\right)h\left(Y_i|Y\right)𝑑Y𝑑Z.$$
(48)
Both $`Z`$ and $`Y`$ are intrinsically positive, but we here assume that the uncertainties in the measurments are small enough so that the probability distributions $`f`$ and $`h`$ may be taken to be Gaussians,

$$f\left(Z|Z_i\right)=\frac{1}{\sqrt{2\pi ϵ_i^2}}\mathrm{exp}\left[\frac{\left(ZZ_i\right)^2}{2ϵ_i^2}\right],$$
(49)
$$h\left(Y_i|Y\right)=\frac{1}{\sqrt{2\pi \delta _i^2}}\mathrm{exp}\left[\frac{\left(YY_i\right)^2}{2\delta _i^2}\right].$$
(50)
As for our model, we assume that

$$g\left(Y|Z\right)=\frac{1}{\sqrt{2\pi \sigma ^2}}\mathrm{exp}\left[\frac{\left(YmZb\right)^2}{2\sigma ^2}\right].$$
(51)
The parameters to be determined by the data are $`m`$, which contains the physics of the metallicity dependence; $`b`$, which is the mean value for the 4He abundance; and $`\sigma `$, which represents the intrinsic scatter about the mean.

We use the integration in equation (48) to obtain

$$P\left(Y_i|Z_i\right)=\frac{1}{\sqrt{2\pi }}\frac{1}{\sqrt{\sigma ^2+m^2ϵ_i^2+\delta _i^2}}\mathrm{exp}\left[\frac{\left(Y_imZ_ib\right)^2}{2\left(\sigma ^2+m^2ϵ_i^2+\delta _i^2\right)}\right].$$
(52)
If we determine the metallicity and helium abundance for $`N`$ objects, the probability for the entire measurement to occur is

$$P_N=\underset{i}{\overset{N}{}}P\left(Y_i|Z_i\right).$$
(53)
To maximize this probability with respect to the parameters $`m`$, $`b`$, and $`\sigma `$, we minimize the $`\chi ^2`$-like function

$$\chi ^2=\underset{i=1}{\overset{N}{}}\left[\mathrm{ln}\left(\sigma ^2+\delta _i^2+m^2ϵ_i^2\right)+\frac{\left(Y_imZ_ib\right)^2}{\sigma ^2+m^2ϵ_i^2+\delta _i^2}\right].$$
(54)
A downhill simplex method and a direction set method yielded almost identical parameters for the minimum $`\chi ^2`$: $`m=4.8\times 10^3`$, $`b=0.244`$, and $`\sigma =1.52\times 10^3`$. We used a Monte Carlo method to set confidence limits on these results. Because of random errors, the observed data points $`Z_i`$ and $`Y_i`$ are not a unique realization of the true values of the quantities. Instead, for a given measurement of $`Z_i`$, there are an infinite number of $`Y_i`$ that could have been measured. The probabilities for obtaining these hypothetical measurements are assumed to be given by equation (52). We have generated $`10^4`$ of these hypothetical data sets by substituting the observed values for $`m`$, $`b`$, and $`\sigma `$ into equation (52). For each hypothetical data set, we again minimized equation (54) to obtain a new triplet of parameters $`m`$, $`b`$, and $`\sigma `$. Figure Element Diffusion During Cosmological Structure Formation shows the resulting intercepts $`b`$ plotted against the scatters $`\sigma `$. The uncertainty in our best-fit sigma and intercept are large, and so the observational situation must improve before we could settle on the true parameters.

## 6 Conclusions

We have explored the effect of diffusion during the linear growth of density perturbations in the neutral IGM before the redshift of reionization and have found that element abundances within collapsed objects differ from their primordial values. Element diffusion leads to two related effects. First, the mean element abundances of a collapsed object could differ from the primordial values; and second, diffusion could result in abundance fluctuations within a single object. For objects with a total mass below $`10^7`$M the former effect dominates; but for larger objects the latter effect is more important. The magnitude of the effect is $`0.005`$$`0.2\%`$, depending on the mass of the object and the element. We have then used a Press-Schechter mass function to average these distributions over all collapsed objects in the Universe, and again found the net effect to be $`<0.1\%`$. This effect is yet too small to be detected since the best observations of element abundances currently have accuracies of $`1\%`$.

If the effect of diffusion on element abundances will eventually be measured with sufficient precision in low-metallicity systems, the results could provide important clues about early structure formation in the IGM. In particular, the mean deviations within objects are sensitive to the minimum mass of a galaxy halo. Also, the overall magnitude of diffusion is sensitive to the redshift of reionization since diffusion is brought to a halt at that time. And as determinations of element abundances get more and more accurate, diffusion will eventually have to be accounted for in order to make comparisons with the predictions of big bang nucleosynthesis.

###### Acknowledgements.

We thank Rennan Barkana for useful comments on the manuscript. This work was supported in part by NASA grants NAG 5-7039, 5-7768, and by NSF grants AST-9900877, AST-0071019.


format: text

text:
# Information loss in an optimal maximum likelihood decoding

## Acknowledgements

I thank Bill Bialek, Anna Montagnini and Alessandro Treves for very useful discussions. This work has been partially supported with a grant of Proff. Treves, of the Human Frontier Science Program, number RG 01101998B.

## References

* Bialek, W., Rieke, F., de Ruyter van Steveninck, R. R., & Warland, D. (1991). Reading a neural code. Science, 252, 1854 - 1857.
* Cover, M. T., & Thomas, J. A., (1991). Elements of Information Theory. New York: Wiley.
* Eskandar, E. N., Richmond, B. J., & Optican, L., M. (1992). Role of inferior temporal neurons in visual memory. I. Temporal encoding of information about visual images, recalled images, and behavioural context. J. Neurophysiol., 68, 1277 - 1295.
* Fano, R. M. (1961) Transmission of Information: A Statistical Theory of Communications. New York: MIT.
* Georgopoulos, A. P., Schwartz, A., & Kettner, R. E. (1986). Neural population coding of movement direction. Science, 233, 1416 - 1419.
* Gochin, P. M., Colombo, M., Dorfman, G. A., Gerstein, G. L., & Gross, C. G. (1994). Neural ensemble encoding in inferior temporal cortex. J. Neurophysiol, 71, 2325 - 2337.
* Golomb, D., Hertz, J., Panzeri, S., Treves, A., & Richmond, B. (1997). How well can we estimate the information carried in neuronal responses from limited samples? Neural Comp., 9, 649 - 655.
* Heller, J., Hertz, J. A., Kjaer, T. W., & Richmond, B. J. (1995). Information flow and temporal coding in primate pattern vision. J. Comput. Neurosci., 2, 175 - 193.
* Kjaer, T. W., Hertz, J. A., & Richmond, B. J. (1994). Decoding cortical neuronal signals: networks models, information estimation and spatial tuning. J. Comput. Neurosci., 1, 109 - 139.
* Kullback, S., (1968). Information theory and statistics. New York: Dover.
* Optican, L. M., & Richmond, B. J. (1987). Temporal encoding of two dimensional patterns by single units in primate inferior temporal cortex: III Information theoretic analysis. J. Neurophysiol., 57, 162 - 178.
* Panzeri, S., & Treves, A. (1996). Analytical estimates of limited sampling biases in different information measures. Network, 7, 87 - 107.
* Panzeri, S., Treves, A., Schultz, S., & Rolls, E. T. (1999). On decoding the responses of a population of neurons from short time windows. Neural Comput., 11, 1553 - 1577.
* Rieke, R., Warland, D., de Ruyter van Steveninck, R. R., & Bialek W., (1996). Spikes: Exporing the Neural Code. Cambridge: MIT Press.
* Rolls, E. T., Critchley, H. D., & Treves, A. (1996). Representation of Olfactory Information in the Primate Orbitofronal Cortex. J. Neurophysiol., 75, (5), 1982 - 1996.
* Rolls, E. T., Treves, A., & Tovée, M. J. (1997). The representational capacity of the distributed encoding of information provided by populations of neurons in primate temporal visual area. Exp. Brain. Res., 114, 149 - 162.
* Rolls, E. T., & Treves, A. (1998). Neural Networks and Brain Function. Oxford: Oxford University Press.
* Rolls, E. T., Treves, A., Robertson, R. G., Georges-Francois, P., & Panzeri, S. (1998). Information About Spatial View in an Ensemble of Primate Hippocampal Cells. J. Neurophysiol., 79, 1797 - 1813.
* de Ruyter van Steveninck, R. R., & Laughlin, S. B. (1996). The rates of information transfer at graded-potential synapses. Nature, 379, 642 - 645.
* Seung, H. S., & Sompolinsky, H. (1993). Simple models for reading neural population codes. Proc. Nac. Ac. Sci. USA, 90, 10749 - 10753.
* Shannon, C. E. (1948). AT&T Bell Laboratories Technical Jounal ,27, 379 - 423.
* Tovée, M. J., Rolls, E. T., Treves, A. & Bellis, R. J. (1993). Information encoding and the responses of single neurons in the primate temporal visual cortex. J. Neurophysiol., 70, 640 - 654.
* Treves, A., & Panzeri, S. (1995). The upward bias in measures of information derived from limited data samples. Neural Comp., 7, 399 - 407.
* Treves, A., Skaggs, W. E., & Barnes, C. A. (1996). How much of the hippocampus can be explained by functional constraints? Hippocampus, 6, 666 - 674.
* Treves, A. (1997). On the perceptual structure of face space. BioSyst., 40, 189 - 196.
* Victor, J. D., & Purpura, K. P. (1996). Nature and precision of temporal coding in visual cortex: a metric space analysis. J. Neurophysiol., 76, 1310 - 1326.
* Victor, J. D., & Purpura, K. P. (1997). Metric-space analysis of spike trains: theory, algorithms and application. Network 8 127 - 164
* Wilson, M. A., & McNaughton, B. L. (1993). Dynamics of the Hippocampal Ensemble Code for Space. Science, 261, 1055 - 1058.


text:
# Low-energy p-d Scattering: High Precision Data, Comparisons with Theory, and Phase-Shift Analyses

## I Introduction

An important question in low-energy nuclear physics is whether a nucleon-nucleon (NN) interaction is sufficient to describe few-nucleon systems or whether a three-nucleon force (3NF) is necessary when the system consists of A$`3`$. To solve for the bound and scattering states of the three-nucleon (3N) system, state-of-the-art theoretical methods like the Faddeev method  (which works well for n-d scattering) have been developed. It was not until the advancement of variational techniques  which could include the Coulomb force rigorously that serious comparisons could be made with p-d scattering data. Current theoretical investigations incorporate phenomenological NN potentials such as CD Bonn , Nijmegen III , and Argonne $`v_{18}`$ (AV18) . These potentials are based on the one-pion-exchange potential for large interparticle separations as well as fits to the present NN database. When calculations using those models are performed in the three-nucleon system, the $`{}_{}{}^{3}\mathrm{He}`$ and $`{}_{}{}^{3}\mathrm{H}`$ binding energies are underestimated by approximately $`0.5`$ to $`1.0`$ MeV depending on which NN potential is employed. The under-prediction of the $`A=3`$ binding energies is the first indication that the two-nucleon interaction alone, as parameterized by the NN potentials in their current forms, is insufficient for an accurate description of the 3N data. To remedy the situation, phenomenological 3N potentials such as the Tuscon-Melbourne (TM) , the Brazil (BR) , and the Urbana IX (UR-IX)  were introduced which are based on the 2$`\pi `$ exchange interaction and adjusted to reproduce the triton binding energies. By applying chiral-perturbation theory to the 3NFs, Friar et al.  showed that the present phenomenological 3NFs are equivalent to first-order in the Lagrangian they constructed except for one term in the TM potential which was found to be of next-to-leading order.

Beyond the binding-energy problem, there is another large discrepancy in the properties of the 3N system, namely the vector analyzing powers (VAPs) $`A_y`$ and $`iT_{11}`$ for N-d scattering. Comparisons between variational calculations  including one of the NN potentials and VAP measurements at $`E_{c.m.}2`$ MeV  show that the calculations under-predict the data by $`30`$% in the maximum of the angular distributions. The disagreement becomes worse at lower energies where the difference increases to $`40`$% at $`E_{c.m.}=431.3`$ keV . This VAP problem was first observed for $`A_y`$ in n-d scattering and has been labelled the “$`A_y`$ puzzle” . The trend of the “$`A_y`$ puzzle” for both $`A_y`$ and $`iT_{11}`$ is a decreasing discrepancy as the energy increases as shown in Fig. 1. The disagreement mostly disappears at $`E_{c.m.}30`$ MeV . The inclusion of a 3NF in the variational calculations lessens the discrepancies with $`A_y`$ and $`iT_{11}`$ data for p-d scattering at the aforementioned energies by no more than $`15`$%. The effect of the magnetic moments was investigated by Stoks  and found not to provide the solution.

To investigate this question further, we have measured with high precision $`A_y`$ and $`iT_{11}`$ as well as the cross section $`\sigma (\theta )`$ and the tensor analyzing powers (TAPs) $`T_{20}`$, $`T_{21}`$, $`T_{22}`$ for p-d elastic scattering at $`E_{c.m.}=667`$ keV. To isolate the deficient phase shifts and mixing parameters in the theoretical calculations, two single-energy phase-shift analyses are performed, one with these data at $`E_{c.m.}=667`$ keV and another with the same observables at $`E_{c.m.}=431.3`$ keV  previously measured by our group. Additional information concerning the experimental details and analysis is available in Ref. . Moreover, the 667 keV cross-section measurement was included in a $`\chi ^2`$ study  with other cross-section data below $`E_{c.m.}=2`$ MeV.

## II Experimental Details

The basic experimental setup and techniques were developed for our measurements of d-p elastic scattering at $`E_{c.m.}=431.3`$ keV and are described in Ref. . The modifications to the setup described in Ref.  for the experiments conducted at $`E_{c.m.}=667\pm 1`$ keV will be discussed. One advantage of making measurements at $`E_{c.m.}=667`$ keV over 431.3 keV is the reduction of the Rutherford scattering of the incident particles from carbon which was present in the targets. The statistical uncertainties in this set of measurements are significantly smaller than in the measurement at 431.3 keV. The only slight disadvantage is the presence of additional proton groups from the $`{}_{}{}^{12}\mathrm{C}(d,p){}_{}{}^{13}\mathrm{C}`$ reaction in the deuteron beam experiments. All six angular distribution measurements employed the FN tandem accelerator at the Triangle Universities Nuclear Laboratory to give a deuteron (proton) beam energy in the center of the target of 2.00 MeV (1.00 MeV).

### A Cross-Section Measurements

The relative cross-section was measured with a beam of deuterons on thin hydrogenated carbon targets  which contained approximately $`(0.5\pm 0.1)\times 10^{18}\mathrm{H}/\mathrm{cm}^2`$ and $`(1.0\pm 0.2)\times 10^{18}\mathrm{C}/\mathrm{cm}^2`$. The rest of the experimental details can be found in Ref. . A sample spectrum for $`{}_{}{}^{1}\mathrm{H}(d,d){}_{}{}^{1}\mathrm{H}`$ scattering at $`\theta _{lab}=26.1^{}`$ is shown in Fig. 2. The absolute cross-section was obtained by normalizing the relative $`\sigma (\theta )`$ for d-p scattering to the well-known $`{}_{}{}^{1}\mathrm{H}(p,p){}_{}{}^{1}\mathrm{H}`$ cross-section determined by the Nijmegen phase-shift analysis . The uncertainties of the absolute cross-section measurement are 0.5% from counting statistics, 0.5% from the normalization to the p-p scattering cross-section, and 0.6% from the beam-current integration. The absolute $`\sigma (\theta )`$ angular distributions are shown in Fig. 3a. The data are divided by the results of the variational calculation to emphasize the differences between the data and the calculations.

### B Analyzing Power Measurements

For the five analyzing power measurements, the techniques used were similar to those described in Refs. . The targets  contained $`(1.0\pm 0.2)\times 10^{18}`$ hydrogen or deuterium atoms and $`(2.0\pm 0.4)\times 10^{18}`$ carbon atoms per $`\mathrm{cm}^2`$. The beam polarization for each experiment was measured online with a polarimeter situated behind the scattering chamber. For the $`T_{20}`$, $`T_{21}`$, and $`T_{22}`$ experiments, a polarimeter  based on the $`{}_{}{}^{3}\mathrm{He}(\stackrel{}{d},p){}_{}{}^{4}\mathrm{He}`$ reaction was utilized. The beam polarizations for the $`T_{20}`$ and $`T_{22}`$ measurements were $`p_{zz}\pm 0.8`$ and $`p_{zz}0`$ for the three spin states with a 3% error. For the $`T_{21}`$ measurement, the beam polarizations were $`p_{zz}\pm 0.8`$ and $`p_{zz}0`$ with a 2% error. For $`iT_{11}`$, a polarimeter  based on the $`{}_{}{}^{2}\mathrm{H}(\stackrel{}{d},p){}_{}{}^{3}\mathrm{H}`$ reaction was used, and the beam polarizations were $`p_z\pm 0.55`$. For the proton beam, the $`{}_{}{}^{6}\mathrm{Li}(\stackrel{}{p},{}_{}{}^{3}\mathrm{He}){}_{}{}^{4}\mathrm{He}`$ reaction  was employed for polarimetry. The beam polarization over the duration of the experiment was $`p_z\pm 0.65`$. The error in the beam polarizations for both experiments was 2%.

The coincidence technique described in Ref.  was used for the VAP experiments. A new addition was 2-$`\mu \mathrm{m}`$-thick mylar foils that were placed in front of the detectors to block any heavy recoil particles such as carbon atoms. By making the measurements at the higher energy, it was possible to cover a larger angular range ($`\theta _{c.m.}50^{}140^{}`$) than obtained at 431.3 keV. A typical timing spectrum for p-d coincidences for a deuteron beam incident on a hydrogen target is shown in Fig. 4. The TAP and VAP angular distributions are given in Figs. 3 and 5, respectively.

## III Comparisons with the Theoretical Calculations

Theoretical calculations to compare with the data have been produced using the AV18 potential and the AV18 plus a 3NF. These calculations employ the Pair-Correlated Hyperspherical Harmonic basis to construct the scattering wave functions and the Kohn variational principle to solve for the scattering-matrix elements . In order to study the sensitivity to different 3NF models, the UR-IX, the TM and its modified version denoted TM interactions have been considered. Originally the TM potential  was given as a sum of four terms characterized by the strengths $`a,b,c,d`$. To assure consistency with chiral symmetry the $`a`$-term in the TM model has been redefined with a new constant $`a^{}=a2m_\pi ^2c`$ (which essentially implies a change of sign for the constant $`a`$) and the $`c`$-term has been set to zero . The TM model, as well as the TM and the UR-IX potentials, have been used recently in an extensive study of 3NF effects in n-d scattering above $`E_{c.m.}=2`$ MeV . For completeness, here we will present results using the same set of constants as Ref.  for the TM and TM models.

Accompanying the data presented in Figs. 3 and 5 are several separate variational calculations, one including the AV18 potential and the others with the AV18 and either the UR-IX or TM potentials. It is clear that the calculations with or without the 3N potential reproduce $`\sigma (\theta )`$ and the TAP data fairly well. On the other hand, the agreement with the VAP data is much poorer. These results are qualitatively similar to those observed at higher energies which were mentioned in the Introduction. For the present case, the differences between the NN-potential-only calculations and the $`A_y`$ and $`iT_{11}`$ measurements are approximately 40% and 37%, respectively at the maximum of the angular distributions. With the inclusion of the UR-IX potential, the agreement for all of the observables improves; however, the change in the VAPs is only marginal and a large discrepancy still remains on the order of 36% and 29% for $`A_y`$ and $`iT_{11}`$, respectively. When the TM 3N-potential is included, the VAP calculations show slightly worse agreement with the data (39% and 34% differences) than in the case of the UR-IX potential. However, when the TM potential is used, the calculations are closer to the VAP data than for the case of the UR-IX potential; the difference between $`A_y`$($`iT_{11}`$) data and the calculations is 32%(25%). The same trend has been observed in Ref.  at $`E_{c.m.}=2`$ MeV. The difference between the calculations with the TM potential and those with UR-IX is about 5%. While each of the 3NFs mentioned here provide some improvement to the $`A_y`$ and $`iT_{11}`$ calculations, none of the potentials eliminates the descrepancy. To obtain a more quantitative comparison, the reduced $`\chi ^2`$, $`\chi _N^2`$, was calculated for each observable, and the results are described in Sec. IV.

Besides the TM potential, other modified NN and 3N force models have recently become available. For example, a new NN interaction constructed from chiral perturbation theory seems to give a better description of the VAP data at low energy . A different possibility may be in the construction of a new 3N force. Hüber and Friar  have suggested that a possible candidate should be a spin-dependent 3NF since the VAP is a difference between polarized cross-sections. Following this direction, a phenomenological 3N $`\stackrel{}{L}\stackrel{}{S}`$ force has recently been proposed . This new force essentially modifies the scalar function $`v_{11}^{ls}(r)`$ already present in the $`\stackrel{}{L}\stackrel{}{S}`$ term of the AV18 potential in triplet spin and isospin channels. Explicitly the following form has been proposed

$$V_{3N}^{ls}=\underset{i<j}{}v_{11}^{ls}(r_{ij})𝐋_{ij}𝐒_{ij}P_{11}(ij)+W_0\mathrm{e}^{\alpha \rho }\underset{i<j}{}𝐋_{ij}𝐒_{ij}P_{11}(ij).$$
(1)
where $`P_{11}`$ is a projector in channels with spin $`S_{ij}=1`$ and isospin $`T_{ij}=1`$. The hyperradius $`\rho `$ is

$$\rho =\sqrt{\frac{2}{3}(r_{12}^2+r_{23}^2+r_{31}^2)}$$
(2)
and $`W_0`$ and $`\alpha `$ are parameters characterizing the strength and range of the three-body term. This force influences the splitting of the $`{}_{}{}^{4}P_{J}^{}`$ phases and the magnitude of $`ϵ_{3/2^{}}`$ without appreciably affecting the other parameters. Moreover, when the deuteron and spectator nucleon are separated by a large distance, the proposed interaction reduces to the original NN $`\stackrel{}{L}\stackrel{}{S}`$ potential $`v_{11}^{ls}(r_{ij})`$. Fig. 6 shows $`A_y`$ and $`iT_{11}`$ calculated using this new potential for three different strength and range parameters. The values for $`W_0`$ and $`\alpha `$ have been taken from Ref.  and have been fixed in order to reproduce the VAP data at $`E_{c.m.}=2`$ MeV. The addition of the 3N spin-orbit force has provided much better agreement with our VAP data than previous calculations. The various phase-shift parameters will be discussed in detail in Sec. V. Moreover the calculated differential cross sections and TAPs remain essentially unchanged  by the new term.

## IV Calculations of Reduced $`\chi ^2`$

To investigate quantitatively the effects of the UR-IX potential, reduced $`\chi ^2`$ calculations were completed with both the AV18 potential and the AV18+UR-IX potentials. More details for the $`\sigma (\theta )`$ comparison are presented in Ref. , and the results are listed in the first two lines of Table I. None of the data were renormalized in order to find a minimal $`\chi _N^2`$. From $`\chi _N^2`$ calculations, it is possible to draw conclusions similar to those from the visual comparisons: the $`\sigma (\theta )`$ and the TAP data are reproduced very well while VAP data show poorer agreement with theory. However, all of the $`\chi _N^2`$ are improved when a 3NF is included. For $`\sigma (\theta )`$, $`\chi _N^2`$ improves by an order of magnitude. The $`\chi _N^2`$ is reduced significantly for each TAP. The calculations of $`A_y`$ and $`iT_{11}`$ with UR-IX show improved agreement with the data; however, the $`\chi _N^2`$ remains in the hundreds.

## V Phase-Shift Analysis

The technique of Phase Shift Analysis (PSA) is useful for determining the dependences of the observables on the various partial waves, and also for studying the sources of disagreements between theory and experiment. Here we consider data at two energies: $`E_{c.m.}=667`$ keV, the energy of the present data, and also $`E_{c.m.}=431.3`$ keV, the energy of previous measurements in our laboratory . The PSAs reported here are carried out in the framework of Seyler , where the scattering matrix is parameterized in terms of phase shifts and mixing parameters. The phase shift parameters in general depend on orbital angular momentum $`L`$, total angular momentum $`J`$, and the channel spin $`S`$. The mixing parameters connect partial waves with different $`L`$ and/or $`S`$ for a given $`J`$ and parity.

Before embarking upon a detailed PSA, it is important to understand the convergence of the observables versus the angular momentum cutoff. Studies were performed using the AV18+UR-IX calculations for the phase shifts and mixing parameters. The dependence on the angular momentum cutoff $`L_{max}`$ was examined by setting to zero all phase shifts and mixing parameters with $`L>L_{max}`$ for different values of $`L_{max}`$. Note that the effect of the Coulomb potential is retained in all partial waves, so it is expected that the observables will converge quickly. We have quantified the convergence of the observables by using $`\chi _N^2`$ comparisons of the calculations with different $`L_{max}`$ values to the experimental data at $`E_{c.m.}=431.3`$ and 667 keV. The results are shown in Table II. First of all we find very little change in $`\chi _N^2`$ between $`L_{max}=3`$ and $`L_{max}=4`$, indicating, as expected, that the results are already well-converged with $`L_{max}=3`$ at these low energies. It is also interesting to note that the cross section and TAPs are fairly well described with $`L_{max}=2`$, while only $`L_{max}=1`$ is required to produce a similar level of convergence for the VAPs. Smaller, but statistically-significant improvements are obtained by extending to $`L_{max}=3`$. We have also observed in these studies that the observables $`T_{20}`$ and $`T_{22}`$ are highly dependent on the $`L=2`$ phase shifts and mixing parameters.

Single-energy PSAs have been performed by using the phase shifts and mixing parameters produced by the AV18+UR-IX potential as starting values. Some of the parameters are allowed to vary freely; the best-fit values for the variable parameters are determined by $`\chi ^2`$ minimization with attention given to the fitting procedure to avoid local minima. Phase shifts were included up to a maximum angular momentum of $`L=4`$, but only those with $`L2`$ were allowed to vary.

Previous phase-shift analyses  have linked the VAP problem to the $`P`$-wave phase-shifts. To test the effects of the different $`P`$-wave phase-shifts and their associated mixing parameters, two-parameter fits where made. Tables I and III show the $`\chi _N^2`$ at the two energies for sample trials where only the $`P`$-wave phase-shifts were varied. By varying the $`{}_{}{}^{2}P_{J}^{}`$ phase-shifts and/or the mixing parameter $`ϵ_{1/2^{}}`$ freely (trials 1 and 2), there is very little change in the $`\chi _N^2`$ with the exception of $`iT_{11}`$ at 431.3 keV. With the combination of $`ϵ_{3/2^{}}`$ and a $`{}_{}{}^{4}P_{J}^{}`$ phase-shift (trials 3, 4, and 5), good fits were produced with the most dramatic improvement occurring in each VAP. At $`E_{c.m.}=667`$ keV, the $`\chi _N^2`$ reduced by a factor of approximately 300 for $`A_y`$ and a factor of 30 for $`iT_{11}`$. For the $`E_{c.m.}=431.3`$ keV case, the fits produced reductions in $`\chi _N^2`$ by factors of approximately 20 and 10 for $`A_y`$ and $`iT_{11}`$, respectively, with the exception of the fit to $`iT_{11}`$ with the free parameters $`{}_{}{}^{4}P_{3/2}^{}`$ and $`ϵ_{3/2^{}}`$. A comparison of the three $`{}_{}{}^{4}P_{J}^{}`$ trials indicates that the two-parameter fit with the $`{}_{}{}^{4}P_{1/2}^{}`$ phase-shift and $`ϵ_{3/2^{}}`$ gives the best result although the differences with the other two trials are very small. These trials so far illustrate the influence of the $`P`$-wave phase-shifts and $`ϵ_{3/2^{}}`$; however, a more realistic fit (trial 6) is to allow all of the $`S`$\- and $`P`$-wave phase-shifts with $`ϵ_{1/2^{}}`$ and $`ϵ_{3/2^{}}`$ to vary freely (9 parameters in all). The remaining phases and mixing parameters were taken from the AV18+UR-IX calculations . The results of the fits are listed in Table IV.

While there were improvements with the full $`S`$\- and $`P`$-wave fits over the two parameter $`{}_{}{}^{4}P_{J}^{}`$ fits, the reduction in $`\chi _N^2`$ is small. To make a more quantitative comparison between trials, errors were determined for specific parameters by the change in a specific parameter necessary to increase the minimal $`\chi ^2`$ by 1 . Allowing more parameters to vary led to a larger error in each parameter which underlies our attempts to find the best fit with the least number of parameters. The uncertainties in phase shifts and mixing parameters when all $`S`$\- and $`P`$-wave parameters were varied are given in Table IV. These errors reflect the contribution of statistical uncertainties in the data only.

The PSA presented here produced very good fits to the data sets at $`E_{c.m.}=431.3`$ keV and 667 keV with sensible results. The greatest influence on the $`A_y`$ and $`iT_{11}`$ fits came from the $`{}_{}{}^{4}P_{J}^{}`$ phase-shifts and $`ϵ_{3/2^{}}`$ parameter. This result was noticed in the earlier analyses  at higher energies. It is difficult to provide a more precise value from the PSA for each phase-shift with the present data sets at 431.3 keV and 667 keV. What these analyses do confirm is the need to modify the $`{}_{}{}^{4}P_{J}^{}`$ phase-shifts and the mixing between the $`{}_{}{}^{4}P_{3/2}^{}`$ and the $`{}_{}{}^{2}P_{3/2}^{}`$ phase-shifts (as evident from the change in $`ϵ_{3/2^{}}`$). In fact, $`ϵ_{3/2^{}}`$ was increased (in absolute value) by about $`19\%`$ from the 667 keV PSA and by about $`28\%`$ from the 431.3 keV PSA over the predicted value even when the UR-IX potential was included in the calculations. The reader should be cautioned that the errors from both analyses allow much wider ranges of values for $`ϵ_{3/2^{}}`$. The theoretical calculations also appear to underestimate the splitting $`\mathrm{\Delta }P={}_{}{}^{4}P_{5/2}^{}{}_{}{}^{4}P_{1/2}^{}`$. For example at 667 keV the calculated splitting is $`\mathrm{\Delta }P=0.87^{}`$ compared to $`0.96^{}`$ obtained in the fit. At 431.3 keV the difference is even bigger; the calculated value is $`\mathrm{\Delta }P=0.48^{}`$ compared to $`0.63^{}`$ from the fit. The small, but nevertheless important, underestimation of $`\mathrm{\Delta }P`$ has been found also in the calculations at $`E_{c.m.}=2`$ MeV .

## VI Conclusions

In conclusion, we have presented high-precision measurements of $`\sigma (\theta )`$, $`A_y`$, $`iT_{11}`$, $`T_{20}`$, $`T_{21}`$, and $`T_{22}`$ for p-d elastic scattering at a center-of-mass energy of $`667`$ keV. The cross section and TAP data show excellent agreement with the variational calculations using the AV18 and UR-IX potentials. The “$`A_y`$ puzzle” persists at this energy with the inclusion of the UR-IX in the calculations giving a small improvement in the agreement with the VAP data. In order to extend the analysis to different forms of the 3NFs, the AV18+TM and AV18+TM potential models have also been considered in the comparisons to the data. The calculations with either TM potential under-predict the VAP data by a similar amount as in the case of the AV18+UR-IX results. From the reduced $`\chi ^2`$ analysis, it is clear that adding a 3NF into the calculations improves the agreement with the data for all of the observables.

From a single-energy PSA at very low energies ($`E_{c.m.}=431.3`$ keV and 667 keV), the under-prediction of the VAP data by the calculations appears to be caused by certain $`P`$-wave components. By varying two parameters, $`ϵ_{3/2^{}}`$ and one of the $`{}_{}{}^{4}P_{J}^{}`$ phase-shifts, and maintaining all of the other parameters as determined from theory, an almost perfect fit was achieved for $`A_y`$. For $`iT_{11}`$, the $`\chi _N^2`$ was reduced by more than an order of magnitude to approximately 1.3 for the 431.3 keV PSA and 1.8 for the 667 keV analysis. To obtain an even better fit, all of the $`S`$\- and $`P`$-wave phase-shifts along with $`ϵ_{1/2^{}}`$ and $`ϵ_{3/2^{}}`$ were varied freely to obtain $`\chi _N^2`$ of $`1.0`$ for each observable. The results of this analysis are the same within the errors as the PSA with $`ϵ_{3/2^{}}`$ and a $`{}_{}{}^{4}P_{J}^{}`$ phase-shift as the only free parameters. From comparisons of the PSA results to calculations it is evident that the theory underestimates the splitting of the $`{}_{}{}^{4}P_{J}^{}`$ phase shifts as well as the mixing between $`{}_{}{}^{2}P_{3/2}^{}`$ and $`{}_{}{}^{4}P_{3/2}^{}`$.

The inclusion of a phenomenological 3N spin-orbit force, which fixes the $`P`$-wave mixing in the calculations, provides better agreement with the VAP data without disturbing the agreement obtained for the cross-section and TAP data. While this fact does not determine the existence of a 3NF, it indicates a need to develop and refine the current 3N potential models . In addition, one can not rule out the possibility that currently employed NN $`{}_{}{}^{3}P_{J}^{}`$ interactions are incorrect .

Finally, the $`A_y`$ problem can be studied with information from p-$`{}_{}{}^{3}\mathrm{He}`$ elastic scattering data. The variational calculations employed in this paper have recently been extended to the p-$`{}_{}{}^{3}\mathrm{He}`$ scattering observables . In this case there is an underprediction by approximately 40% of new $`A_y`$ data taken at an $`E_{c.m.}`$ below 2 MeV. Moreover, $`A_y`$ for p-$`{}_{}{}^{3}\mathrm{He}`$ scattering is an order of magnitude larger than for the p-d scattering case. This fact makes it an ideal candidate for future investigations, and more p-$`{}_{}{}^{3}\mathrm{He}`$ elastic scattering data at low energies is needed.

###### Acknowledgements.

The authors would like to thank Dr. K. D.Veal for his assistance in the data taking. This work was supported in part by the US DOE under Grant No. DE-FG02-97ER41041.


text:
# 𝜔 Centauri as a Disrupted Dwarf Galaxy: Evidence from Multiple Stellar Populations

## 1. Introduction

As part of our investigation of the luminosity-metallicity relation of the RR Lyrae stars in $`\omega `$ Cen, we have obtained $`BV`$ and $`Ca`$ & Strömgren CCD frames with the CTIO 0.9m telescope that cover 40 $`\times `$ 40 $`arcmin^2`$ in a 3 $`\times `$ 3 grid centered on the cluster. As a byproduct of this investigation, we obtained high-quality homogeneous $`BV`$ color-magnitude (CM) data for more than 130,000 stars in the field toward $`\omega `$ Cen, which represents the most extensive photometric survey to date for this cluster as far as the stars brighter than the main-sequence (MS) turnoff are concerned. This study (Lee et al. 1999) has shown, for the first time, that $`\omega `$ Cen has several distinct stellar populations, which is reminiscent of the Sagittarius dwarf galaxy (Layden & Sarajedini 1997). This feature was later confirmed by Pancino et al. (2000) from their independent photometry. In this paper, we present our progress in more detailed analysis of the data along with the population models.

## 2. Multiple Stellar Populations in $`\omega `$ Centauri

Figure 1 presents $`BV`$ CM diagram for all stars in our program field including RR Lyrae variable stars. We were able to remove some foreground Galactic disk star contamination, as they are well discriminated from more metal-poor member stars of $`\omega `$ Cen in the analysis based on our $`Ca`$ & Strömgren $`by`$ photometry (see Rey et al., this volume). Here we confirm our previous discovery that the red-giant-branch (RGB) of $`\omega `$ Cen has four distinct populations: the most metal-poor sequence, metal-poor sequence, metal-rich sequence, and the most metal-rich sequence well separated from others. We believe this is a clear evidence for the multiple stellar populations in $`\omega `$ Cen.

The horizontal-branch (HB) distribution is also consistent with the discrete nature of RGB. We have the blue HB and RR Lyrae variables mainly associated with the two metal-poor RGBs, some red HB stars associated with the metal-rich population, and finally the red clump superimposed on the RGB, which must be associated with the most metal-rich population.

## 3. Population Models

The presence of multiple stellar populations in $`\omega `$ Cen is confirmed by our population models (Figure 2; see also Ree et al., this volume). Our models are constructed based on the $`Y^2`$ isochrones (Yi et al. 2001) and corresponding HB evolutionary tracks (Yi, Demarque, & Kim 1997) calculated with updated input physics. The reader is referred to Lee, Demarque, & Zinn (1994) and Park & Lee (1997) for the details of model constructions. From the population models, we confirm four distinct populations with different metallicities can reproduce the observed discrete nature of the RGB. The models also confirm that the RR Lyrae variables are mainly produced by two metal-poor populations, with some contribution from the 2nd most metal-rich population.

For the RR Lyrae stars in $`\omega `$ Cen, we have new metallicity measurements from our $`Ca`$ & Strömgren $`by`$ photometry (see Rey et al., this volume), and therefore more detailed comparison with the models is possible. From the luminosity-metallicity relation and the period-shift – metallicity correlation of the RR Lyrae variables in $`\omega `$ Cen, we confrim that our models reproduce even the fine details within the instability strip, including the sudden upturn of the RR Lyrae luminosity and corresponding increase in period-shift at \[Fe/H\] $``$ -1.5, which is a result of redward evolution from the blue HB (see Rey et al., this volume; Yoon & Lee, this volume). This suggests we have a very good understanding of what is happening in the instability strip of $`\omega `$ Cen.

Our models also reproduce the red clump associated with the most metal-rich RGB, but only when the age of the most metal-rich population is some 4 Gyr younger than the most metal-poor population. This is illustrated in Figure 3, where we can see the variation of the red clump location in the CM diagram under different assumptions regarding the age difference betwen the most metal-rich and the most metal-poor populations. By comparing these models with the observed color of the red clump, which is estimated from the analysis of the luminosity functions of the RGBs superimposed on the red clump (see Rey et al., this volume), we conclude that the $`\mathrm{\Delta }`$t of $``$ 4 Gyrs is our best estimate for the age spread within the $`\omega `$ Cen, in the sense that the metal-rich population is younger. This is in qualitative agreement with the results obtained directly from the Strömgren photometry of main-sequence stars (Hughes & Wallerstein 2000; Hilker & Richtler 2000), although these results are subject to large errors associated with the photometry and metallicity groupings.

These unusual characteristics of $`\omega `$ Cen are very similar to the case of Sagittarius dwarf galaxy, which is in the process of disruption with the Milky Way (Ibata, Gilmore, & Irwin 1994). In Figure 4, we have compared our population models with the observations centered on the M54 (Layden & Sarajedini 2000), which is believed to be the nucleus of the Sagittarius dwarf galaxy. Note again that the blue HB is from the most metal-poor (and older) population, the red HB is from the intermediate metallicity population, and finally the red clump is from the most metal-rich (and younger) population. We confrim here that the Sagittarius dwarf galaxy also has several distinct stellar populations with the internal age-metallicity relation that spans $``$ 7 Gyrs. Other Local Group dwarf galaxies are also known to have similarly complex star formation histories (Mateo 1998).

## 4. Discussion

### 4.1. The Origin of $`\omega `$ Centauri

The multiple populations and the internal age-metallicity relation found in our work (and also suggested by others in this workshop), which resemble the characteristics of Sagittarius dwarf galaxy, suggest $`\omega `$ Cen was massive enough for self-enrichment and several bursts of star formation. The relatively extended enrichment period of $``$ 4 Gyrs then indicates that the initial evolution of $`\omega `$ Cen occured away from the dense central regions of young Milky Way, because, near the Galactic center like the current location of $`\omega `$ Cen, we would expect the gaseous materials to have been stripped from the proto $`\omega `$ Cen on a much shorter timescale as a result of disk shocking and/or similar processes. This view is not inconsistent with $`\omega `$ Cen’s rather unusual kinematics and orbit (Dinescu, this volume).

From the discrete nature of RGB, some may suggest that the $`\omega `$ Cen could be a merger of several globular clusters. But mergers of (at least) four clusters are very unlikely, if not impossible, in the Milky Way, considering the rather high velocities between the clusters in the halo. We can not rule out this scenario in dwarf galaxies, however, because internal velocity dispersion is much lower in dwarf galaxies (van den Bergh 1996).

All of these information strongly suggest that the $`\omega `$ Cen is a relict or a nucleus of a disrupted dwarf galaxy. The case of $`\omega `$ Cen and that of the Sagittarius dwarf system therefore provide direct evidence for past and continuing accretion of protogalactic fragments, which suggest that similar accretion events may have continued throughout the history of Milky Way formation.

### 4.2. Other Globular Clusters with Multiple Populations?

It is interesting to note that two globular clusters, now identified as the nuclei or parts of disrupted dwarf galaxies, are among the most massive globular clusters in the Milky Way. $`\omega `$ Cen is, of course, the most massive globular cluster, and M54, the nucleus of the Sagittarius dwarf, is the 2nd most massive globular cluster. Then, how about the other massive globular clusters? There are several Galactic globular clusters that resemble $`\omega `$ Cen, but among them, we found the cases of NGC 6388 and NGC 6441 are the most interesting. They are 3rd and 5th most massive globular clusters, respectively, and they all have very peculair bimodal HB distributions. Our population models for these unusual clusters show that the adoption of two distinct populations within the systems and very mild internal age-metallicity relations between the two populations can reproduce the observed features on the HB and RGB (see Ree et al., this volume). This conclusion is also supported by the fact that the mean period of RR Lyrae variables in these clusters are too long for their high metallicities, which is understood in two populations secnario, where the RR Lyraes are highly evolved stars from the older and metal-poor blue HB population (see Ree et al., this volume).

This suggests some massive globular clusters in the Milky Way are probably not genuine globular clusters, but infact relicts of disrupted dwarf galaxies like $`\omega `$ Cen and M54.

### 4.3. Classification of Galactic Globular Cluster System

Our conclusion on the origin of $`\omega `$ Cen and other massive globular clusters, when combined with other recent findings on the origin of Galactic globular cluster system, suggest that the present day Galactic globular clusters are infact subdivided into three different types:

(a) Clusters formed in a collapsing proto-Galaxy (e.g., Eggen, Lyndel-Bell, & Sandage 1962). They are mostly in the inner halo and are only genuine Galactic globular clusters.

(b) Clusters originally formed in satellite dwarf galaxies later accreted to the Milky Way. They include globular clusters (Pal 12, Ter 7, Ter 8, & Arp 2) belong to Sagittarius dwarf galaxy, and young outer halo clusters with retrograde motion (Zinn 1993; van den Bergh 1993). Yoon & Lee (this volume) also reported very compelling evidence that metal-poor (\[Fe/H\] $`<`$ -2.0) Oosterhoff II clusters have the positional and orbital characteristics fully consistent with the hypothesis that they were originated from a satellite galaxy.

(c) Nuclei (or relicts) of disrupted dwarf galaxies. $`\omega `$ Cen and M54, the nucleus of the Sagittarius dwarf, are undoubtedly belong to this type. Our population models (Ree et al., this volume) suggest other massive globular clusters with bimodal HBs (e.g., NGC 6388 and NGC 6441) may also belong to this category.

## References

Eggen, O. J., Lynden-Bell, D., & Sandage, A. 1962, ApJ, 136, 748

Ibata, R. A., Gilmore, G. & Irwin, M. J. 1994, Nature, 370, 194

Hilker, M., & Richtler, T. 2000, A&A, 362, 895

Hughes, J., & Wallerstein, G. 2000, AJ, 119, 1225

Layden, A. C., & Sarajedini, A. 1997, ApJ, 486, L107

. 2000, AJ, 119, 1760

Lee, Y. -W., Demarque, P., & Zinn, R. 1994, ApJ, 423, 248

Lee, Y. -W., Joo, J. -M., Sohn, Y. -J., Rey, S. -C., Lee, H. -c., & Walker, A. R. 1999, Nature, 402, 55

Mateo, M. 1998, ARA&A, 36, 435

Pancino, E., Ferraro, F. R., Bellazzini, M., Piotto, G., & Zoccali, M. 2000, ApJ, 534, L83

Park, J. -H., & Lee, Y. -W. 1997, ApJ, 476, 28

van den Bergh, S. 1993, ApJ, 411, 178

. 1996, ApJ, 471, L31

Yi, S., Demarque, P., & Kim, Y. -C. 1997, ApJ, 482, 677

Yi, S., Demarque, P., Kim, Y. -C., Lee, Y. -W., Ree, C. H., Lejeune, T., & Barnes, S. 2001, ApJS, 136, 417

Zinn, R. J. 1993, in The Globular Cluster–Galaxy Connection, ed. G. Smith & J. Brodie (San Francisco: ASP), 38


text:
# Self-assembled guanine ribbons as wide-bandgap semiconductors

## Abstract

We present a first principle study about the stability and the electronic properties of a new biomolecular solid-state material, obtained by the self-assembling of guanine (G) molecules. We consider hydrogen-bonded planar ribbons in isolated and stacked configurations. These aggregates present electronic properties similar to inorganic wide-bandgap semiconductors. The formation of Bloch-type orbitals is observed along the stacking direction, while it is negligible in the ribbon plane. Global band-like conduction may be affected by a dipole-field which spontaneously arises along the ribbon axis. Our results indicate that G-ribbon assemblies are promising materials for biomolecular nanodevices, consistently with recent experimental results.

The possibility that biomolecules may conduct current is intriguing in the development of molecular electronics. In particular, because of their sequence-specific properties, DNA molecules and their derivatives are attractive for nanometer scale electronics, where their self-assembling capability might be used to wire electronic-materials in a programmable way braun , as well as to form novel conducting materials danny ; rinaldi ; rinaldi1 .

Indeed, in recent experiments rinaldi ; rinaldi1 , a semiconducting behavior was displayed by a device constituted of a metallic nanogate filled with a dried solution of a lipophilic deoxyguanosine derivative. This material is composed of hydrogen-bonded planar ribbons, which tend to stack and form solid-state fibers. On a scale of $`100`$ nm, the ribbons self-organize in ordered structures, giving rise to a film of nanocrystalline character rinaldi ; rinaldi1 ; gottarelli . The results of transport and photocurrent experiments identify different regimes, depending on the width ($`L`$) of the nanogate: when $`L100200`$ nm, few nanoscrystals are present between the electrodes and the I-V characteristics present the typical metal-semiconductor-metal behavior; on the contrary, when the gap is reduced ($`L<100`$ nm), a single crystal is trapped between electrodes, and the device acts as a rectifying diode. Two concurrent mechanisms have been proposed to explain electron transport in such devices: band-like conduction within the single crystals, and incoherent hopping among neighboring grains. Band conduction and hopping mechanism have been also invoked to describe charge migration in other guanine-rich systems (e.g. G-quartets) or poly(dG)-poly(dC) DNA danny ; barton ; schuster ; giese ; jortner ; however, no conclusive agreement was reached so far in this respect.

In this paper, we present a first principle study of the electronic and conduction properties of self-assembled guanine ribbons: using an extended-solid approach, we analyze the effects of H-bonding and $`\pi `$-$`\pi `$ coupling, that control the microscopic interactions along the ribbon axis and along the stacking direction, respectively. Our method is based on ab-initio molecular dynamics simulations: we performed total energy and force calculations, optimizing at the same time the atomic positions and the single-particle electronic wave-functions, according to the Car-Parrinello approach fhi . The electronic structure is described within the Density Functional Theory (DFT) dft ; Generalized Gradient Approximation (BLYP-GGA) blyp  is applied in exchange-correlation functional to allow an accurate description of H-bonding. The atomic potentials are described by ab-initio soft pseudopotentials in the Troullier-Martins mt  formulation. The electronic wave functions are expanded in a plane-wave basis with an energy cutoff of 50 Ry.

With respect to real deoxyguanosine ribbons, our model leaves out the sugar (deoxyribose) and the attached paraffinic chains, because they are not involved in transport phenomena barton1 ; ladik . Since guanine costitutes the building block of the ribbon aggregates, we first study the geometry and the electronic structure of the isolated molecule (see Figure 1, inset). Our method provides an accurate description of the structural parameters: bond lengths and angles are reproduced with an error smaller than 2% with respect to X-ray experimental data saenger . The total charge density has a highly asymmetric distribution, which induces a net dipole ($`|\overline{\mu }|=7`$ D) oriented as shown in the inset of Figure 1. The analysis of the electronic structure reveals $`\sigma `$ and $`\pi `$ electron states. The out-of-plane nature of $`\pi `$ orbitals, in particular of the HOMO (Fig. 2a) and of the LUMO, may control the interactions in the stacking direction.

To simulate extended ribbons, we assembled coplanar G molecules according to the experimental geometry gottarelli . By optimizing the structure of an isolated G dimer, we obtained the intermolecular distance. We used this value to fix the periodicity of the extended structures, whose geometry was also optimized. The periodicity of the ribbons (1.1nm) is in good agreement with the experimental data (1.2nm). The model ribbons (Fig. 1) are 1D planar wires of guanines connected by a H-bond network. The donor-acceptor distances (NH$`\mathrm{}`$ N=2.9 Å, NH$`\mathrm{}`$ O=2.8 Å) are consistent with previous theoretical results for similar systems (e.g. DNA base-pairs) sponer .

H-bonding is a weak interaction, not capable of modifing the internal configuration of the single guanine; we find that its main effect is to stabilize the structure, as proved by large calculated adsorption energy ($`\mathrm{\Delta }E=810`$ meV/G). On the other hand, H-bonding does not favor the lateral coupling of electron states between neighboring molecules. As shown in Figure 2b, where the charge density plot of the HOMO is presented, the ribbon states are localized around the single guanines, with no Bloch-type delocalization through the hydrogen-bonds. The density of states (DOS, top of Fig. 3a) shows isolated sharp peaks corresponding to the energy levels of G, and the energy bands (not shown here) are dispersionless. We thus conclude that band-like conductivity, through extended states, is not a viable mechanism for electron transport along the ribbon axis.

The picture changes drastically when we consider stacked configurations. In a previous work prb , we demonstrated that the electronic properties of vertical stacks of G molecules are very sensitive to the $`\pi \pi `$ coupling: A strong $`\pi \pi `$ superposition is mandatory to have a Bloch-like conduction along the stacking direction. Here, we present a stacked configuration of planar ribbons, where the stacking direction is orthogonal to the plane of the ribbon. Even though experimental results prove that the ribbons self-organize in fibers rinaldi , the X-ray structure of single crystal is still missing, therefore the real stacking geometry is unknown. We report our results for a configuration in which the ribbons are perfectly eclipsed one above each other and the band dispersion is expected to be maximized. In this configuration, we investigated the concurrent effects of H-bonding in the ribbon plane and $`\pi \pi `$ coupling along the stacking direction: a reduced value of the adsorption energy ($`\mathrm{\Delta }E=310`$ meV/G), with respect the isolated wire, is the index of electrostatic repulsion of $`\pi \pi `$ superposition.
We then examined the electron states near the gap edge, that affect the conduction properties; these are $`\pi `$ orbitals very sensitive to the vertical superposition. As an example, the LUMO of the eclipsed structure is shown in Fig. 2c: the $`\pi \pi `$ coupling induces the formation of Bloch-type orbitals delocalized along stacking direction. On the other hand, as it happens for isolated ribbons, the electronic states are not extended in the plane, indicating that H-bonding gives negligible contribution to band conductivity. This is consistent with the bandstructure (Fig. 3b). The HOMO- and LUMO-derived bands are dispersive along the stacking direction ($`\mathrm{\Gamma }A`$); while they are almost flat along the ribbon axis ($`\mathrm{\Gamma }X`$). The calculated effective masses for electrons and holes are $`m_e=1.24`$ and $`m_h=1.07`$ (in units of the free electron mass) along the stacking direction, whereas they are very high along the axis of the ribbons. We conclude that the $`\pi \pi `$ interaction enhances band properties, improving the mobility of the charge carriers only along the stacking direction. The DOS (bottom of Fig. 3a) for the eclipsed ribbons confirms this behavior, and shows wider ranges of allowed energies with respect to isolated ribbons, expecially in the region around the gap. Band properties of fibers may be similar to those of inorganic wide-bandgap semiconductors (e.g. GaN, AlN pugh ) for which Bloch conduction has been proved in presence of doping. Our results are in agreement with the experimental evidence of a semiconducting behavior of guanine ribbons.

In order to explain the rectifying properties of single nanocrystals, we must consider the polar nature of the G-ribbons. Due to the relative orientation of the single G’s in the ribbons, molecular dipoles arrange to form a net dipole oriented along the ribbon axis (Fig. 1). The calculated values of dipoles, for finite sequences of 2 ($`|\overline{\mu }|=14.6`$ D) and 4 ($`|\overline{\mu }|=32.8`$ D) guanines, indicate that a dipole-dipole interaction arises along ribbons. The associated electric field may couple with proper Bloch-momentum of charge carriers, defining a different effective direction for electron transport. The presence of a favorite direction justify the diode-like character of experimental I-V curves.

In summary, we have demonstrated that supramolecular aggregates of guanines may form solid-state materials with interesting conduction properties. Ribbon fibers present wide-bandgap semiconducting behavior, compatible with Bloch-like conduction along stacking direction. A strong dipole field influences band properties of such nanocrystals, inducing preferential direction in movement of electron and holes, as confirmed by recent experimental results. rinaldi .

Acknowledgements
Funding was provided by INFM through Progetto Calcolo Parallelo, and PRA-SINPROT.


text:
# Critical behaviour of the Ginzburg-Landau model in the type II region

## Abstract

We study the critical behaviour of the three-dimensional U(1) gauge+Higgs theory (Ginzburg-Landau model) at large scalar self-coupling $`\lambda `$ (“type II region”) by measuring various correlation lengths as well as the Abrikosov-Nielsen-Olesen vortex tension. We identify different scaling regions as the transition is approached from below, and carry out detailed comparisons with the criticality of the 3d O(2) symmetric scalar theory.

The three-dimensional (3d) U(1) gauge+Higgs theory (Ginzburg-Landau (GL) model, scalar electrodynamics) is an effective theory for phase transitions in superconductors, liquid crystals, and possibly cosmology. The Lagrangian is

$$_{GL}=\frac{1}{4}F_{ij}^2+|D_i\varphi |^2+m_3^2|\varphi |^2+\lambda _3|\varphi |^4.$$
(1)
Here $`\varphi `$ is a complex scalar field, $`F_{ij}=_iA_j_jA_i`$, and $`D_i=_i+ie_3A_i`$. The parameters $`m_3`$, $`e_3^2`$, $`\lambda _3`$ have the dimension GeV. Physics then depends on the dimensionless ratios

$$x\lambda _3/e_3^2,ym_3^2/e_3^4.$$

The GL model has no local order parameters. However, there are non-local quantities, like the photon mass $`m_\gamma `$, and the vortex tension $`T`$ , which vanish in the symmetric and are non-zero in the broken phase. The schematic phase diagram is shown in Fig. 1. At small $`x`$ (type I region) the transition is of the first order, while at large $`x`$ (type II region) the transition is believed to be of the second order, with a (presumably) tricritical point in between, at $`x0.3`$ .

In this work we study the phase transition deep in the type II region, at $`x=\lambda _3/e_3^2=16`$, shown as a thick vertical line in Fig. 1, and compare the scaling behaviour with that of the 3d XY model. The universality class of the transition in the type II domain has not been previously resolved unambiguously with lattice simulations.

The scaling behaviour of the system is determined by the longest correlation lengths, or the lightest “masses”. We shall consider the photon and scalar masses $`m_\gamma `$ and $`m_s`$ as functions of $`\tau (yy_c)`$. The relation between masses varies strongly as $`\tau `$ is varied; consequently, in the broken phase, one can argue that there are three distinct possibilities:

1. Mean field: The textbook case, valid if $`\tau 0`$, leading to $`m_\gamma (|\tau |/x)^{1/2}m_s(2|\tau |)^{1/2}`$.

2. XY scaling: Here scaling is dominated by the dynamics of the scalar field; indeed, since $`\lambda _3e_3^2`$, it is not unreasonable for this to happen. Such scaling is supported by RG arguments , and has been observed experimentally . In this domain $`m_s|\tau |^{\nu _{\mathrm{xy}}},m_\gamma |\tau |^{\nu _{\mathrm{xy}}/2},`$ with $`\nu _{\mathrm{xy}}0.67`$, and we expect $`m_s<2m_\gamma `$, otherwise the scalar would decay into two photons.

However, we do not expect this scaling to remain valid very close to the critical point: the XY model has a massless Goldstone at $`\tau 0`$, whereas the GL model has a massless photon only when $`\tau 0`$. Moreover, due to logarithmic confinement in 3d U(1) gauge theory, one can expect the scalar to remain massive at the transition point.

3. Inverted XY scaling: Here the GL “temperature” parameter $`\tau `$ is supposed to map to the inverted temperature ($`\tau `$) of the 3d XY model, and the broken/symmetric phases of the GL model are supposed to correspond to the symmetric/broken phases of the XY model. This scenario is supported by duality arguments . When $`\tau <0`$ the GL model in (2+1)d has two massive photon polarisation states, which become critical at $`\tau =0`$. At $`\tau >0`$, these d.o.f’s become the massless photon, and a massive vector “resonance”. These counts match exactly the two massive scalars in the symmetric phase of the XY model, which, in the broken phase, become the massless Goldstone and the massive radial mode.

Identification of the photon with the critical degree of freedom of the XY model suggests $`m_\gamma |\tau |^{\nu _{\mathrm{xy}}}`$; however, other estimates in the literature suggest $`m_\gamma |\tau |^{0.5\mathrm{}0.67}`$ . The scalar remains massive, since all diverging d.o.f’s have been accounted for. (Alternatively, the scalar mass could vanish, but with an exponent smaller than the photon one).

Duality arguments also relate the Abrikosov-Nielsen-Olesen vortex tension $`T`$ to the scalar mass in the dual theory, giving $`T|\tau |^{\nu _{\mathrm{xy}}}.`$

We use the standard lattice discretization of Eq. (1), with a non-compact U(1) gauge action. For details, see . We choose a fixed lattice spacing ($`e_3^2a=1`$); observing the universal behaviour does not require taking the continuum limit.

Unfortunately, in the GL model there is no gauge invariant local observable which would correspond to magnetization. Vortex percolation has often been used to probe the critical behaviour ; however, on the lattice this is not rigorous . We shall here only measure quantities such as the specific heat, scalar and photon correlation lengths, and the vortex tension.

In this paper we report on correlation lengths. In the GL model, these are measured using the following operators (for concreteness, measurements taken along the $`x_3`$-direction):

| Scalar | $`S=\varphi ^{}\varphi `$, |
| --- | --- |
| Photon (plaquette) | $`B_3=ϵ_{123}F_{12}=\mathrm{}_{12}`$, |
| Vector | $`V_i=\text{Im }\varphi ^{}D_i\varphi `$ . |

The vector and photon operators are measured at finite “transverse momentum” to, say, the $`x_1`$-direction. Additionally, each channel contains states with multiple photons, complicating the scalar channel measurement in particular. In the symmetric phase the scalar is always a resonance, decaying into two or more photons. For each operator we use several levels of blocking + smearing to reduce noise, and we diagonalize the full cross-correlation matrix between the blocked operators, as described, e.g., in .

In Fig. 2 the photon mass squared is shown. The critical exponent is $`2\nu 0.59`$. Inverted XY scaling is usually thought to imply a larger exponent, $`1.0<2\nu <2\nu _{\mathrm{xy}}1.34`$ (see above). Thus, we conclude that inverted XY scaling, at least in this form, is not valid for these $`\tau =yy_c`$.

In the scalar channel we include up to four lowest two-photon states at zero total transverse momentum in the cross-correlation analysis, requiring diagonalization of matrices up to $`14^2`$, and improve thus on earlier work at $`x=2`$ . In Fig. 3 we show measurements on $`64^3`$ lattices. Three of the states clearly correspond to two photons, and the scalar ($`\varphi ^{}\varphi `$) is a resonance in the symmetric phase.

The scalar (resonance) mass is shown separately in Fig. 4. Approaching $`y_c`$ from the broken phase, the mass appears to vanish with a critical exponent $`m_s|\tau |^{0.48}`$. The value does not agree well with “standard” XY scaling, $`\nu _{\mathrm{xy}}0.67`$, but it is marginally compatible with the inverted XY scaling scenario, where the scalar mass should either remain finite at $`y_c`$, or at least have a critical exponent smaller than $`\nu _{\mathrm{xy}}`$.

In conclusion, the photon correlation length does not display (at least with the signatures commonly suggested) inverted XY scaling, as predicted by duality arguments. In addition, the critical exponent for the vortex tension $`T`$ is several standard deviations off from the inverted XY scaling prediction (to be reported elsewhere).

Of course, one possible reason for the discrepancy is that at $`x=16`$ the scaling window may be an extremely narrow band around $`y_c`$, rendering it almost impossible to observe the true scaling in Monte Carlo simulations.