text:
# Recovering the Topology of the Initial Density Fluctuations Using the IRAS Point Source Catalogue Redshift Survey

## INTRODUCTION

The IRAS Point Source Catalogue Redshift Survey  lends itself to topological studies because of its high sampling densities and large volume. It countains approximately $`15000`$ galaxies to the full depth of the PSC ($`0.6Jy`$). Its sky coverage is $`84.1`$ per cent, where only the zone of avoidance is excluded, here defined as an infrared background exceeding $`25\mathrm{M}\mathrm{J}\mathrm{y}\mathrm{s}\mathrm{r}^1`$ at $`100\mu \mathrm{m}`$, and a few unobserved or contaminated patches at higher latitude. The excluded regions are coded in an angular mask, as shown in Fig. 1. Canavezes et al.  analysed the topology of PSCz and showed that it is consistent with the topology of CDM models: The genus curves retain the w-shape characteristic of random-phase density fields even at small smoothing lengths, where non-linear evolution has already generated significant skewness of the 1-point PDF. These non-linearities are, however, detected by a depressed amplitude of the genus curves - amplitude drops - consistent with those detected for the CDM models. Above $`10h^1\mathrm{Mpc}`$ strong phase correlations were not detected in PSCz and those that were detected at smaller scales are expected in the framework of mildly non-linear gravitational evolution. This supports the hypothesis that structure grew from random-phase initial conditions.

In this paper we attempt to subtract these phase correlations by reversing gravity. If it is true that structure does indeed originate in random-phase fluctuations, then reversing gravitational evolution should provide us with an initial density fluctuation field free from phase correlations on any scale.

There are several methods to achieve this goal. Nusser & Dekel , showed how to express the Zel’dovich approximation in a set of Eulerian coordinates and reverse it to obtain a reconstructed density field. They tested the method with N-body simulations and found satisfactory results. Recently, some improvements have been made on the original method (e.g. Gramann ). Nusser, Dekel & Yahil  applied the Nusser & Dekel  approximation to the $`1.2`$ Jy $`IRAS`$ redshift survey to recover the 1-point probability distribution function of the initial density field. Their results were consistent with Gaussian initial conditions. The PSCz presents us with an unprecedented number of resolution elements and should constitute therefore the best available data set on which to apply a reconstruction technique.

The method we choose to follow here is the quasi-linear method proposed by Nusser & Dekel . In section 1 we describe this method; In section 2 we describe the construction of the density maps used as the input of the time-machine and describe our error estimates; Section 3 is devoted to testing the method using N-body simulations of a CDM model and in section 4 we apply the method to the PSCz and present our results.

## 1 THE RECONSTRUCTION METHOD

The reconstruction method we follow is based on two premises: The velocity field is, when smoothed on a scale of a few Mpc, irrotational, and the Zel’dovich approximation is accurate over the mildly non-linear regime.

Our first step is to try to express the Zel’dovich approximation in a set of Eulerian coordinates.

Given an initial comoving position $`𝐪`$ for a given particle, a final position $`𝐱(𝐪,𝐭)`$ will have the form

$$𝐱(𝐪,t)=𝐪+P(𝐪,t).$$
(1)
The Zel’dovich approximation states that the displacement term $`P(𝐪,t)`$ can be written as a product of two functions, each a function only of one of the variables $`𝐪`$ or $`t`$, i.e., we can separate the variables $`q`$ and $`t`$:

$$𝐱(𝐪,t)=𝐪+D(t)\mathrm{\Psi }(𝐪).$$
(2)
This is simply a linear approximation with respect to the particle displacements rather than density. Particles in the Zel’dovich approximation follow straight lines:

$$𝐯(𝐪,t)=a(t)\frac{d𝐱}{dt}=a(t)\dot{D}(t)\mathrm{\Psi }(𝐪).$$
(3)
The density fluctuation is given by

$$\delta (𝐪,t)=\overline{\rho }/J(𝐪,t)1,$$
(4)
where $`J(𝐪,t)`$ is the Jacobian of the coordinate transformation $`𝐪𝐱`$.

As long as second order terms in $`\delta `$ and $`𝐯`$ are negligible, $`\delta (t)D(t)`$, where $`D(t)`$ is the growing mode solution in linear perturbation theory (see e.g. Peebles ).

Unlike linear theory, where the density at a given position evolves according to the linear growth rate, under the Zel’dovich approximation infinite density can develop in a finite time as a result of the convergence of particle trajectories into a pancake. in order to integrate a density field back in time, Nusser & Dekel  found the differential equation in Eulerian space which contains the Zel’dovich approximation. This is derived from the standard equation of motion of dust particles in an expanding Universe (e.g. Peebles ):

$$\frac{d𝐯}{dt}+H𝐯=\frac{1}{a}\mathrm{\Phi }_g.$$
(5)
Here $`𝐯`$ stands for the peculiar velocity of the particle, i.e., $`𝐯=a(t)\frac{d𝐱}{dt}`$, where $`𝐱`$ are the comoving positions; $`=(/x,/y,/z)`$ and $`\mathrm{\Phi }_g`$ is the gravitational potential which is related to the local density fluctuation via the Poisson equation

$$^2\mathrm{\Phi }_g=\frac{3}{2}H^2\mathrm{\Omega }a^2\delta .$$
(6)
Using the normalized variables $`\theta `$ and $`\phi _g`$ defined thus:

$`\theta (𝐱,t)`$ $``$ $`{\displaystyle \frac{𝐯(𝐱,t)}{a\dot{D}}}=\psi (𝐪)`$ (7)
$`\varphi _g(𝐱,t)`$ $``$ $`{\displaystyle \frac{\mathrm{\Phi }_g(𝐱,t)}{a^2\dot{D}}},`$ (8)
equations 5 and 6 reduce to

$$\dot{\theta }+\frac{3H\mathrm{\Omega }}{2f(\mathrm{\Omega })}\theta =\phi _g,$$
(9)
where $`f(\mathrm{\Omega })\dot{D}/HD\mathrm{\Omega }^{0.6}`$ .

Notice that the Zel’dovich approximation is contained in 7, which is to say

$$\dot{\theta }=0$$
(10)
along the trajectory of a given particle, i.e., along a line of constant $`𝐪`$. Hence, 9 is reduced, under the Zel’dovich approximation, to

$$\frac{3H\mathrm{\Omega }}{2f(\mathrm{\Omega })}\theta =\phi _g.$$
(11)
$`\theta `$ is therefore the gradient of a potential:

$$\theta =\phi _v,$$
(12)
with

$$\varphi _v=\frac{2f(\mathrm{\Omega })}{3H\mathrm{\Omega }}\phi _g+F(t).$$
(13)
We can obviously set $`F(t)0`$ since $`\theta `$ is the only physical (measurable) quantity.

Equation 10 can now be expanded in Eulerian coordinates:

$`{\displaystyle \frac{d\theta }{dt}}=0`$ $``$ $`\left[{\displaystyle \frac{}{t}}+{\displaystyle \frac{d𝐱}{dt}}\right]\theta (𝐱,t)=0`$ (14)
$``$ $`{\displaystyle \frac{\theta }{t}}+\dot{D}(t)\psi (𝐪)\theta =0`$
$``$ $`{\displaystyle \frac{\theta }{t}}+\dot{D}\left(\theta \right)\theta =0.`$
Making use of the identity $`\left(\theta \right)\theta =1/2\left(\theta \right)^2\theta \left(\theta \right)`$, taking into account the irrotationality of $`\theta `$, and substituting 12 we arrive at

$`{\displaystyle \frac{\phi _v}{t}}+\dot{D}{\displaystyle \frac{1}{2}}\left(\theta \right)^2=0`$ $``$ $`\left[{\displaystyle \frac{\phi _v}{t}}+{\displaystyle \frac{\dot{D}}{2}}\left(\phi \right)^2\right]=0`$ (15)
$``$ $`{\displaystyle \frac{\phi _v}{t}}+{\displaystyle \frac{\dot{D}}{2}}\left(\phi _v\right)^2=F(t).`$
Again, because only gradients of $`\phi _v`$ have a physical meaning, we can take $`F(t)`$ to be zero. We finally arrive at the equation:

$$\frac{\phi _v}{D}+\frac{1}{2}\left(\phi _v\right)^2=0.$$
(16)
This is the differential equation which expresses the Zel’dovich approximation in Eulerian coordinates. Knowing $`\phi _v`$ at any time enables us to compute $`\phi _v`$ at any other time simply by integrating 16 forwards or backwards.

The first step is to calculate $`\phi _v`$ at the present time. That can be easily achieved using Poisson’s equation, which in Fourier space reads:

$$k^2\stackrel{~}{\mathrm{\Phi }}_g=\frac{3}{2}H^2\mathrm{\Omega }a^2\stackrel{~}{\delta },$$
(17)
where the tildes denote the Fourier transforms. Using 8 we arrive at:

$$Dk^2\stackrel{~}{\phi }_v=\stackrel{~}{\delta }.$$
(18)
So, in order to obtain the velocity potential of a given smooth density field with periodic boundary conditions, we first transform it to Fourier space using some FFT code, we then divide the obtained field by $`k^2`$ (normalizing D to unity at present), and then transform back to real space to obtain $`\phi _v`$. $`\phi _v`$ is then integrated back in time in the most trivial way: At every step, we calculate $`1/2\left(\phi _v\right)^2`$ from the potential field; this is then used as a first order Taylor correction to predict the value of $`\phi _v`$ at the next step. Once a suitable number of integration steps has been performed, we use equation 18 again to obtain the value of $`\delta `$ at the corresponding value of $`D`$.

For our application - topology, we expect the exact number of integrations to be irrelevant. After a certain number of iterations we expect to reach the linear regime in all scales, which means that the genus curves of the density fields will no longer change. This particular characteristic of the genus can in fact be useful to test the convergence of the method. We obtained several reconstructed density fields for different values of the number of integration steps, and found the topologies to converge. This result is shown is Fig. 2.

## 2 CONSTRUCTION OF THE DENSITY MAPS

The construction of the density maps we follow is, to a large extent, equivalent to the method followed by Canavezes et al. . We employ the same fit for the PSCz selection function $`s(z)`$:

$$s(z)=\frac{\psi }{z^\alpha \left(1+\left(\frac{z}{z^{}}\right)^\gamma \right)^{\beta /\gamma }},$$
(19)
with the parameters shown in Table 1 , and estimate the density $`\rho (r)`$ by

$$\rho (r)\frac{m(r)}{s(r)},$$
(20)
where $`m(r)`$ is the discrete point distribution. However, the angular mask needs to be treated with care. Because the structure we can see behind the mask will have a significant effect on the evolution of structure on the whole observable area, we need to resort to some sort of filling, prior to using the time-machine.

We consider two different types of filling: One, a random filling, in which fake objects are placed randomly over the masked region with an average number density equal to the average number density of the whole survey (weighted by the selection function); And another, a cloning, in which fake objects are placed randomly in each bin, but with an average number density equal to the average number density of the neighbouring observable bins, again weighted by the selection function.

On both cases, we create a box containing a sphere of radius $`R_{max}`$ as defined in Table 2. $`R_{max}`$ is the maximum distance up to which the average distance between two neighbouring galaxies in PSCz is smaller than the adopted smoothing length. The size of the box will then be $`\frac{2}{\sqrt{3}}Rmax`$. Although boundary conditions are not periodic, we smooth the density field in Fourier space assuming periodic boundary conditions. We consider this to be preferable to zero padding, as this would create an artificial boundary that would eventually affect the gravitational time-machine.

We then apply the Zel’dovich time-machine on the box thus obtained for both filling techniques. In the subsequent topological analysis, we limit ourselves to the sphere inside the box with radius $`R_{max}`$.

### Smoothing Procedure

Ideally, one would want to apply the operator time-reverse directly upon the unsmoothed density field obtained from PSCz, and only then smooth it on a range of scales to be able to calculate its genus curve. However, in order to obtain meaningful results, we need a smooth field a priori. This poses a very important problem: These operators (time-reverse and smoothing) do not in general commute. In other words, the topology of the final density field can be significantly different whether the smoothing operation is performed before or after applying the time-machine to our original field.

To illustrate this point, let us consider a Gaussian random density field smoothed on some scale $`\lambda `$. After applying the time-machine this density field will still be Gaussian. However, if we are to naively calculate its genus curve choosing for isodensity contours the same isodensity contours defined prior to applying te time-machine, we would wrongly conclude that the field is not Gaussian. In order to circumvent this problem we would be forced to smooth the final field once again on a larger scale, thus reducing considerably the statistical significance of our results.

There is, however, an alternative way to solve this dilemma: One can try to find a smoothing operator that comutes with our reconstruction operator.

How can we look for such a smoothing operator? One way of ensuring this is to find an operator that does not change the density field if applied once again at a latter stage. By other words, if we ensure that our final density field remains the same regardless of how many times we perform the smoothing operation throughout the time-reversing process, then it is because these two operators do indeed comute. A class of such operators are the adaptive smoothing operators. By adaptive smoothing, we understand a local smoothing of variable smoothing length according to the local structure. There are several types of adaptive smoothing. The ideal adaptive smoothing operator ensures that the total mass enclosed by isodensity contours remains constant throughout the time-reversing process. A simplified version assumes spherical symmetry. At each point we employ a smoothing length $`\lambda `$ such that $`\lambda ^31/\rho =1/\left(\rho _0+\delta \right)`$. The proportionality constant defines the characteristic smoothing length, which is the smoothing length employed when the overdensity $`\delta `$ vanishes. For highly dense regions, a smaller value of $`\lambda `$ will be chosen, whereas voids will be smoothed with a larger value of $`\lambda `$.

When we apply the time-machine operator on a map smoothed in this way, the isodensity contours around a cluster will move, but the total mass enclosed by them will remain approximately constant, as long as the spherical symmetry hypothesis is a good approximation.

In our application to the PSCz we use a spherically symmetric adaptive smoothing algorithm. We start by creating a set of maps obtained from the original map (this might be either the PSCz data or simulation data) by smoothing it on a range of scales around a given characteristic scale $`l_0`$, using Gaussian filters of the form

$$G_\lambda (x)=\frac{1}{\pi ^{3/2}\lambda ^3}e^{x^2/\lambda ^2}$$
(21)
where $`\lambda `$ varies around $`l_0`$. We then look for the appropriate value of $`\lambda `$ at each position $`x`$ by enforcing the equation

$$(1+\delta _\lambda )\lambda ^3=l_0^3,$$
(22)
where $`\delta _\lambda `$ is the density contrast when the original map is smoothed on a scale $`\lambda `$. This ensures that the mass enclosed on a sphere of radius $`\lambda `$ is, to first order in a Taylor expansion, constant throughout the whole map.

### Error Estimates

There are three different sources of error that enter our results: Shot noise, cosmic variance and the errors associated with the Zel’dovich approximation itself.

The most accurate and realistic way of estimating these errors, which ultimately enter the genus curves of the reconstructed PSCz density fields, is to analyse the topology of fake PSCz catalogues drawn from N-body simulations of some standard CDM model. Since we intend to test the validity of the Zel’dovich time-machine with N-body simulations, it seems reasonable to extent this philosophy to the calculation of the statistical errors themselves. This is achieved by seeking a galaxy number density field with a Poisson distribution, whose expectation value is identical to the density field of the N-body simulation multiplied by the PSCz selection function. This is equivalent to observing the density field of the N-body simulation in a similar way to the way PSCz observes the density field of the real Universe. This galaxy number density field is then divided by the PSCz selection function again to obtain a ”PSCZ-noisy” distance-independent estimate of the real galaxy number density. Each of the mock PSCz catalogues is adaptively smoothed using the algorithm described above and then used as input in the Zel’dovich time-machine. The genus curves of the reconstructed fields are calculated and the variance obtained over 10 mock PSCz catalogues is used as the statistical error estimate.

## 3 TESTING THE METHOD WITH N-BODY SIMULATIONS OF CDM MODELS.

As we mentioned previously, we intend to test the regime of validity of the Zel’dovich time-machine by means of an N-body simulation of a CDM model. We apply the time-machine on the present density field drawn from the simulation and compare the density field thus obtained with the density field drawn from the original test-mass positions used as input in the simulation. More specifically, we need only to compare the topologies of the reconstructed and original fields through their genus curves.

For this purpose, we use two N-body simulations corresponding to two cold dark matter models, kindly provided by the Virgo consortium .The simulations have been performed with an AP3M-SPH code named HYDRA . Here we consider the SCDM model and the $`\tau `$CDM model, which parameters are shown in Table 3. Because this simulations contain CDM on periodic boxes of size $`239.5h^1\mathrm{Mpc}`$ and use such a large number of particles, they constitute an ideal ground for this test.

We start by binning the particles in cells of side length $`2h^1\mathrm{Mpc}`$ and smooth further on some characteristic smoothing length $`\lambda `$, using the adaptive smoothing algorithm described in section 2, in an analogous procedure to the way we treat the PSCz galaxies. This ”double” smoothing is required in the real Universe in order to minimize shot noise effects, but it is also necessary in order to eliminate severe non-linearities since the Zel’dovich time-machine is expected not to work over such regimes. It is also needed to smooth out regions of orbit-crossing.

The next step is to calculate the genus curves of both the original density field that is used as input in the simulation, smoothed in the same manner as described above, and the reconstructed density field obtained after applying the Zel’dovich time-machine to the present density field smoothed in the same way. In order to compare both genus curves we compute their amplitudes and their amplitude drops. It is very important to note that if we are to obtain gaussianized versions of the density fields at the present time, they will not have, in principle, the same genus amplitudes as the density fields themselves, even when these are in the Gaussian regime. This is so because of the very nature of adaptive smoothing. When a given density field is smoothed adaptively its genus curve will have a different shape and amplitude than the shape and amplitude of the genus curve of the same field when this is smoothed with a constant Gaussian window. We need to be extremely careful not to draw the wrong conclusions about the Gaussian or non-Gaussian nature of our density fields. After applying the time-machine, however, the amplitude of fluctuations will be reduced significantly. In fact they will be reduced to such an extent as to make the field look almost homogeneous. This means that the isodensity contours on adaptively smoothed maps will be very close to the isodensity contours on maps smoothed using a constant Gaussian window. Hence, it is possible to calculate amplitude drops and determine the Gaussian (or non-Gaussian) nature of our density maps after applying the Zel’dovich time-machine, i.e., to determine the Gaussian (or non-Gaussian) nature of the reconstructed fields.

In order to test the convergence of the reconstruction method, we obtained several reconstructed density fields for diffent values of the total number of integration steps, as mentioned in 1. Fig. 2 shows a particular slice of the reconstructed density fields at different integration steps when the method is applied to the maps at $`z=0`$ obtained from the $`\tau `$CDM simulation by smoothing on characteristic lengths of $`8h^1\mathrm{Mpc}`$ and $`10h^1\mathrm{Mpc}`$. As it is evident from this Figure, the topology of the density fields is undistinguishable for values of $`z`$ greater than $`4`$ on smoothing scales of both $`8h^1\mathrm{Mpc}`$ and $`10h^1\mathrm{Mpc}`$. Only the amplitude of the density fluctuations changes, indicating that we are now in the linear regime. In Fig. 3 this is made even clearer. Here we show a point-by-point comparison of the reconstructed $`\tau `$CDM density fields at redshifts $`z=4`$ and $`z=9`$, for the characteristic smoothing lengths of $`8h^1\mathrm{Mpc}`$ and $`10h^1\mathrm{Mpc}`$. We only show one in eight of all points, chosen randomly. It is obvious from this plot that the shape of fluctuations did not change from $`z=4`$ to $`z=9`$, for either of the smoothing lengths adopted. It is also obvious that the amplitude of fluctuations was reduced by a factor of $`2`$, as it is expected in the linear regime (notice that in the linear regime the growing mode of fluctuations varies as $`1/(1+z)`$). Hence, our time-machine shows the correct assimptotic behaviour.

### Results

Figure 4 shows the isodensity contours obtained for the SCDM simulation and the $`\tau `$CDM simulation, when the fields are adaptively smoothed on a characteristic scale of $`5h^1\mathrm{Mpc}`$ in a comoving box of size $`240h^1\mathrm{Mpc}`$. The first row shows the fields at the present time. The second row shows the fields after applying the Zel’dovich time-machine back to a redshift of 10 and the third row shows the original fields used as input in the simulations. In each panel the thick line represents the contour where the density contrast is zero. Dashed lines represent contours of negative density contrast whereas solid lines represent positive density contrasts. The contour spacing is $`0.2/(1+z)`$. Because this is normalized to the linear reconstruction case where $`\delta 1/(1+z)`$, the linear theory reconstruction density maps look exactly like the present day maps (first row), albeit their much lower fluctuation amplitude. From the contours on both the reconstructed maps and the original maps we notice that their fluctuation amplitudes are larger than the fluctuation amplitude of a linearly reconstructed map. This result is expected because fluctuations will grow faster at the later stages of evolution, in the mildly non-linear regime. It is obvious from Fig. 4 that the Zel’dovich time-machine is able to change the rank order of isodensity contours. This means that the genus curves will also be changed. In fact we can predict, just by looking at Fig. 4 that the amplitude of the genus curves will increase after applying the Zel’dovich time-machine, as it is expected.

Figures 5 and 6 show the genus curves obtained for the SCDM model and the $`\tau `$CDM model respectively, at selected characteristic smoothing lengths, namely at $`5h^1\mathrm{Mpc}`$, $`8h^1\mathrm{Mpc}`$ and $`14h^1\mathrm{Mpc}`$. We restrict our analysis to this range of smoothing lengths, as we know a priori that the Zel’dovich time-machine is only relevant in the mildly non-linear regime. In both figures, the first column shows the genus curves of the reconstructed field and its randomized counterpart, whereas the second column shows the genus curves of the original fields, used as input for the simulations. We deliberately do not show the genus curves of the present day adaptively smoothed density fields, since, as we argued in the previous chapter, the isodensity contours are, at this stage, different from the isodensity contours of the fields smoothed using a constant kernel.

It is striking in Fig. 5 and Fig. 6 the consistency between the genus curves of the reconstructed field and of its randomized counterpart. This is particularly true for the genus amplitudes, which indicates that there is very little phase correlation in the reconstructed fields on all scales, as it is expected from Gaussian initial conditions. With regard to the genus amplitudes per se, we find consistency with the genus amplitudes of the original density fields, although we notice a tendency for these to be slightly higher in the reconstructed fields. These means that the second moment of the power spectrum is recovered slightly in excess. This becomes more important as we approach the smaller scales, which is to say as the system becomes more non-linear.

It is also interesting to notice a depression at high values of $`\nu `$ on the genus curves of the reconstructed fields. We find this depression to be present in all cases, both for the Virgo simulations and the PSCz data. We expect this to be a particular feature of the reconstruction technique, i.e., ultimately dependent on the Zel’dovich approximation itself. At high values of $`\nu `$ the regime is considerably non-linear and so we expect the Zel’dovich approximation to perform rather poorly. This problem is enhanced by the fact that the density fields have been smoothed adaptively. Nevertheless, we do not expect this feature to affect our results considerably as the calculation of the amplitude of the genus curves is restricted to the range $`1<\nu <1`$.

Note that there is no sampling noise in these density fields. The tremble in the genus curves is ultimately due to cosmic variance.

In Fig. 7 we plot the amplitude drops obtained for the reconstructed SCDM model and for the reconstructed $`\tau `$CDM model as a function of smoothing length. Also shown are the amplitude drops obtained for the models at $`z=0`$ . In all scales, the values obtained for the reconstructed fields are closer to unity than those of the present-day fields. This indicates that the Zel’dovich time-machine is a good tool in recovering Gaussian initial density fields. However, at small scales, we detect a slight departure from the expected value of unity, even for the reconstructed fields. This is not surprising since at this scales non-linearities become strong. The degree of strength of the non-linearities depends on the particular model. As we see from Fig. 7 they are more important in the $`\tau `$CDM model than in the SCDM model, in agreement to what was found by Canavezes et al. .

Fig. 8 shows the genus amplitudes of the reconstructed fields for both the SCDM model and the $`\tau `$CDM model, together with the genus amplitudes of the original fields. In the case of the $`\tau `$CDM model the agreement between the original genus amplitude and the genus amplitude of the reconstructed field, is striking, in particular when we consider scales above $`7h^1\mathrm{Mpc}`$. For the SCDM model the agreement is not as good. The reconstructed genus amplitudes appear to be slightly higher than the true original ones.

These results indicate that the Zel’dovich time-machine is effective in recovering the right genus amplitude drops on scales larger than $`8h^1\mathrm{Mpc}`$, although recovered amplitudes cannot be considered reliable in the strict sense of the word.

## 4 RECOVERING THE INITIAL DENSITY FIELD FROM PSCZ

As mentioned previously, one of the difficulties in using the PSCz in an attempt to recover the initial density fluctuations in the Universe, is the fact that the masked region can make up to $`20\%`$ of the whole observable area. Since this region will have a definite gravitational influence on the observable area when we attempt to reconstruct to original density field, some sort of filling is essential. As mentioned in 2 we employ two different techniques: The random filling and the cloning. Fig. 10 shows how the two different fillings appear to the eye. On the upper panel we show the map where the mask has been filled randomly and on the lower panel we show the map where the mask has been cloned. Although it is impossible to the eye to recognize any significant difference, we will be able to detect some noticeable differences on the topologies of the reconstructed density fields.

### Results & Discussion

Fig. 9 shows the genus curves obtained for the reconstructed PSCz fields. On the first column we plot the genus curves of the fields for which the mask has been filled using a cloning technique, and on the second column we plot the genus curves of the fields where the mask has been filled randomly. The thick solid lines refer to the reconstructed PSCz fields proper, whereas the thin solid lines refer to the randomized versions of those fields,i.e., the fields obtained from the reconstructed density fields by randomizing phases in Fourier space subject to the reality constraint $`\delta _𝐤=\delta _𝐤^{}`$ and keeping the same power spectrum. The dotted and dashed lines are the best fitting random-phase curves to both the reconstructed PSCz density fields and their randomized versions.

It is striking to notice the proximity between genus curves, even at small smoothing lengths. The statistical errors drawn from the N-body simulation are not shown in Fig. 9 because of the high degree of correlation between the points in the genus curves.

Independently of whether the mask has been randomly filled or cloned, the recovered genus curves seem consistent with random-phase Gaussian fluctuations because the amplitude drops appear small. However, the amplitude themselves seem to depend on the mask filling technique.

In Fig. 12 we plot the amplitude drops obtained for the reconstructed PSCz density fields, as a function of smoothing length, as well as the error bars drawn from the $`\tau `$CDM model following the method outlined in chapter 2. The amplitude drops of the reconstructed N-body simulations are also shown for comparison.

In contrast with Fig. 11, where phase correlations were found for PSCz at small smoothing lengths, the reconstructed PSCz fields do not show any significant phase correlations on scales ranging from $`5h^1\mathrm{Mpc}`$ to $`14h^1\mathrm{Mpc}`$. This reinforces the hypothesis that density fluctuation originate from random-phase Gaussian initial conditions.

## AKNOWLEDGEMENTS

We are grateful to the Virgo consortium (J. Colberg, H. Couchman, G. Efstathiou, C. S. Frenk, A. Jenkins, A. Nelson, J. Peacock, F. Pearce, P. Thomas and S. D. M. White) for provinding simulation data ahead of publication. AC aknowledges the support of FCT (Portugal).


format: text

text:
# Can the inflaton and the quintessence scalar be the same field?

## 1 Introduction

The very last years have witnessed growing interest in cosmological models with $`\mathrm{\Omega }_m1/3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }2/3`$, following the most recent observational data (see for example the discussion in the paper by Bahcall et al.  and references therein). A very promising candidate for a dynamical cosmological constant is a rolling scalar field, named “quintessence”.  The main motivation for constructing such dynamical schemes resides in the hope of weakening the fine tuning issue implied by the smallness of $`\mathrm{\Lambda }`$. In this respect, a very suitable class of models is provided by inverse power scalar potentials which admit attractor solutions characterized by a negative equation of state.

Consider the cosmological evolution of a scalar field $`Q`$, with potential $`V(Q)=M^{4+\alpha }Q^\alpha `$, $`\alpha >0`$, in a regime in which the scalar energy density is subdominant with respect to the background. Then it can be shown  that the solution $`Qt^{1n/m}`$, with $`n=3(w_Q+1)`$ and $`m=3(w_B+1)`$, is an attractor in phase space. We have defined $`w_Q`$ to be the equation of state of the scalar field $`Q`$, and $`w_B`$ that of the background ($`=1/3`$ for radiation and $`=0`$ for matter). The equation of state of the scalar field on the attractor is found to be $`w_Q=(\alpha w_B2)/(\alpha +2)`$, which is always negative during matter domination. As a consequence, the ratio of the scalar to background energy density is not constant but scales as $`\rho _Q/\rho _Ba^{mn}`$, thus growing during the cosmological evolution, since $`n`$ $`<m`$. The behaviour of these solutions is determined by the cosmological background and for this reason they have been named “trackers” in the literature.  A good feature of these models is that for a very wide range of the initial conditions the scalar field will reach the tracking attractor before the present epoch. Depending on the initial time, you can have several tens of orders of magnitude of allowed initial values for the scalar energy density. This fact, together with the negative equation of state, makes the trackers feasible candidates for explaining the cosmological observation of a presently accelerating universe. The point at which the scalar and matter energy densities are of the same order depends on the mass scale in the potential. This is fixed by requiring that $`\mathrm{\Omega }_Q=𝒪(1)`$ today.

An interesting question, then, is whether the ‘quintessence’ scalar and the inflaton field, which dominate the expansion of the universe at very different times, could indeed be the same field. If this is the case, it should also be possible to uniquely fix the initial conditions for the ‘quintessential’ rolling from the end of inflation. Models in which one single scalar field drives both inflation and the late time cosmological accelerated expansion are named “quintessential inflation” models.

## 2 A particle physics model: Supersymmetric QCD

As first noted by Binètruy , supersymmetric QCD theories with $`N_c`$ colors and $`N_f<N_c`$ flavors may give an explicit realization of a quintessence model with an inverse power law scalar potential. The matter content of the theory is given by the chiral superfields $`Q_i`$ and $`\overline{Q}_i`$ ($`i=1\mathrm{}N_f`$) transforming according to the $`N_c`$ and $`\overline{N}_c`$ representations of $`SU(N_c)`$, respectively. In the following, the same symbols will be used for the superfields $`Q_i`$, $`\overline{Q}_i`$, and their scalar components. Supersymmetry and anomaly-free global symmetries constrain the superpotential to the unique exact form

$$W=(N_cN_f)\left[\mathrm{\Lambda }^{(3N_cN_f)}/\mathrm{det}T\right]^{\frac{1}{N_cN_f}}$$
(1)
where the gauge-invariant matrix superfield $`T_{ij}=Q_i\overline{Q}_j`$ appears. $`\mathrm{\Lambda }`$ is the only mass scale of the theory.

We consider  the general case in which different initial conditions are assigned to the different scalar VEV’s $`Q_i=\overline{Q}_i^{}q_i`$, and the system is described by $`N_f`$ coupled differential equations. In analogy with the one-scalar case, we look for power-law solutions of the form

$$q_{tr,i}=C_it^{p_i},i=1,\mathrm{},N_f.$$
(2)
It is straightforward  to verify that for fixed $`N_f`$ (and when $`\rho _Q\rho _B`$), a solution exists with $`p_ip=p(N_c)`$ and $`C_iC=C(N_c,\mathrm{\Lambda })`$ and that it is the same for all the $`N_f`$ flavors. The equation of state of the tracker is given by

$$w_Q=\frac{1+r}{2}w_B\frac{1r}{2},$$
(3)
where we have defined $`rN_f/N_c`$. Then, even if the $`q_i`$’s start with different initial conditions, there is a region in field configuration space such that the system evolves towards the equal fields solutions (2), and the late-time behavior is indistinguishable from the case considered by Binetrúy  where equal initial conditions for the $`N_f`$ flavors were chosen. In spite of this, the multi-field dynamics introduces some new interesting features. For example, we have found that (in the two–field case) for any given initial energy density such that, for $`q_1^{in}/q_2^{in}=1`$, the tracker is joined before today, there exists always a limiting value for the fields’ difference above which the attractor is not reached in time. A more detailed discussion and numerical results about the two-field dynamics can be found in Masiero et al.

## 3 Quintessential inflation

As already discussed, the range of initial conditions which allows $`\rho _Q`$ to join the tracker before the present epoch is very wide. Nevertheless, it should be noted that in principle we do not have any mechanism to prevent $`\rho _Q^{in}`$ from being outside the desired interval. In this respect, an early universe mechanism which could naturally set $`\rho _Q^{in}`$ in the allowed range for a late time tracking, would be highly welcome. Moreover, if we require the quintessence scalar to be identified with the inflaton, we would at the same time obtain a tool for handling the initial conditions and a simple unified picture of the early and late time universe dynamics, which we call “quintessential inflation”.

The basic idea is to consider inflaton potentials which, as it is typical in quintessence, go to zero at infinity like inverse powers. In this way it is possible to obtain a late time quintessential behaviour from the same scalar that in the early universe drives inflation. The key point resides in finding a potential which satisfies the condition that inflation and late time tracking both occur, and that they occur at the right times (see Peloso et al. ).

Two models have been shown to fulfill these two requirements. One example  is a first-order inflation model with potential going to zero at infinity like $`\varphi ^\alpha `$. A bump at $`\varphi M_p`$ allows for an early stage of inflation while the scalar field gets “hung up” in the metastable vacuum of the theory. Nucleation of bubbles of true vacuum through the potential barrier sets the end of the accelerated expansion and starts the reheating phase. After the reheating process is completed, the quintessential rolling of the scalar $`\varphi `$ starts and its initial conditions (uniquely fixed by the end of inflation) can be shown to naturally be within the range which leads to a present day tracking. As an alternative, we considered  the model of hybrid inflation proposed by Kinney et al.  This is shown to naturally include a late–time quintessential behavior. As typical of hybrid schemes, the potential at early times (that is until the inflaton field is smaller than a critical value $`\varphi _c`$) is dominated by a constant term and inflation takes place. Eventually the inflaton rolls above $`\varphi _c`$, rendering unstable the second scalar of the model, $`\chi `$. This auxiliary field starts oscillating about its minimum and in this stage the universe is reheated. After $`\chi `$ has settled down, the inflaton continues its slow roll down the residual potential, which goes to zero at infinity like $`\varphi ^2`$, thus allowing for a quintessential tracking solution. Also in this case the initial conditions for the quintessential part of the model are not set by hand, but depend uniquely on the value of the inflaton field at the end of reheating.

## Acknowledgments

I would like to thank Antonio Masiero, Massimo Pietroni and Marco Peloso with whom I obtained the results presented in this talk.


text:
# Quasi-Localization of Gravity by Resonant Modes

## Abstract:

We examine the behaviour of gravity in brane theories with extra dimensions in a non-factorizable background geometry. We find that for metrics which are asymptotically flat far from the brane there is a resonant graviton mode at zero energy. The presence of this resonance ensures quasi-localization of gravity, whereby at intermediate scales the gravitational laws on the brane are approximately four dimensional. However, for scales larger than the lifetime of the graviton resonance the five dimensional laws of gravity will be reproduced due to the decay of the four-dimensional graviton. We also give a simple classification of the possible types of effective gravity theories on the brane that can appear for general non-factorizable background geometries.

preprint: hep-th/0002161

The past two years have produced several surprising results in the field of gravity in extra dimensions. Firstly, Arkani-Hamed, Dimopoulos and Dvali  showed that the size of compact extra dimensions could be as large as a millimeter, with the fundamental Planck scale as low as a TeV, without running into contradiction with the current short-distance gravitational measurements, if the standard model fields are localized on a 4D “brane”. Subsequently Randall and Sundrum (RS)  found that using a non-factorizable “warped” geometry for the extra dimension one could in fact have an infinitely large extra dimension, and still reproduce Newton’s Law at large distances on the brane. The key observation of RS is that in their scenario there is a localized zero-energy graviton bound-state in 5D which should be interpreted as the ordinary 4D graviton. In this scenario the geometry of the extra dimension plays a crucial rôle and localization of gravity on the brane is impossible if the geometry far from the brane is asymptotically flat. In view of this fact, it is even more surprising that, as Gregory, Rubakov and Sibiryakov (GRS)  recently showed, even if the geometry is asymptotically flat far from the brane, it is still possible to find a phenomenologically viable model if one does not insist on the Newton potential being valid at arbitrarily large scales, but instead only requires it to hold over a range of intermediate scales (a related proposal can be found in ).

The aim of this paper is to present a physical explanation of the results of GRS and to provide a universal description of warped gravitational theories. We will show that the reason behind the GRS result is that even though there is no zero energy bound-state graviton in their model, there is a resonant mode—a “quasi bound-state”—at zero energy, which plays the rôle of the 4D graviton. The existence of this resonance at zero energy implies the “quasi-localization” of gravity—that is, it tends to produce a region of intermediate scales on which the gravitational laws appear to be four dimensional. We find that the long-distance scale at which gravity appears to be five dimensional again is inversely related to the width of the graviton resonance. The physics behind the appearance of this new scale is that at very large time scales the graviton decays into plane waves away from the brane, and thus reproduces 5D gravity at large distances. As the width of the resonance approaches zero, the lifetime becomes large and in the limit one regains the RS model. On the other hand, if the resonance becomes very wide, it is basically washed out from the spectrum, its effects become unimportant, and there is no longer a region where gravity is effectively 4D.

We show that the existence of such a resonance is expected in these kinds of theories when the geometry is asymptotically flat space far from the brane. Thus we find a simple way of classifying warped gravitational theories: if the ground state wave function is normalizable, one has localization of gravity à la Randall and Sundrum. If the ground state wavefunction is not normalizable, but the geometry does not asymptote to flat space, then there is simply no effective 4D gravity; however, if the ground state wavefunction is non-normalizable and the geometry asymptotes to flat space, we have quasi-localization of gravity via the resonance à la GRS.

The most general 5D metric with 4D Poincaré symmetry can be written

$$ds^2g_{\mu \nu }dx^\mu dx^\nu =e^{A(z)}\left(\eta _{ab}dx^adx^bdz^2\right).$$
(1)
We will assume that the “warp factor” $`A(z)`$ is symmetric and, for simplicity, a non-decreasing function of $`z`$ for $`z>0`$. Furthermore we will also assume that the matter is localized on a brane at $`z=0`$. We now consider fluctuations around the 4D Minkowski metric of the form $`h_{ab}(x,z)=e^{3A(z)/4}\psi (z)\stackrel{ˇ}{h}_{ab}(x)`$, with $`\eta ^{cd}_c_d\stackrel{ˇ}{h}_{ab}(x)=m^2\stackrel{ˇ}{h}_{ab}(x)`$, where $`m`$ is the four-dimensional Kaluza-Klein mass of the fluctuation. The behaviour of the fluctuation in the transverse space is governed by $`\psi (z)`$ which satisfies a Schrödinger-like equation:

$$\frac{d^2\psi (z)}{dz^2}+V(z)\psi (z)=m^2\psi (z),V(z)=\frac{9}{16}A^{}(z)^2\frac{3}{4}A^{\prime \prime }(z).$$
(2)
Notice that (2) always admits a (not necessarily normalizable) zero-energy wavefunction

$$\widehat{\psi }_0(z)=\mathrm{exp}\left[\frac{3}{4}A(z)\right],$$
(3)
which potentially describes the 4D graviton.

The auxiliary quantum system described by (2) encodes all the properties that we need in order to establish whether under favourable conditions there exists an effective 4D Newton potential on the brane. The relevant quantity to consider is the induced gravitational potential between two unit masses on the brane. A discrete (normalized) eigenfunction $`\psi _m(z)`$ of the quantum system contributes

$$\frac{\psi _m(0)^2}{M_{}^3}\frac{e^{mr}}{r},$$
(4)
to the potential, where $`M_{}`$ is the fundamental Plank scale in 5D. On the other hand, continuum modes (normalized as plane waves asymptotically) contribute

$$_{m_0}^{\mathrm{}}𝑑m\frac{\psi _m(0)^2}{M_{}^3}\frac{e^{mr}}{r}.$$
(5)
In the cases which we consider the potential $`V(z)`$ will always approach 0 as $`|z|\mathrm{}`$ and so $`\widehat{\psi }_0(z)`$ is the only possible bound-state of the system and the continuum begins at $`m_0=0`$. The potential has the characteristic volcano shape with a central well surrounded by barriers that decay to zero . We can distinguish 3 classes depending on the asymptotic behaviour of $`\widehat{\psi }_0(z)`$: (a) $`\widehat{\psi }_0(z)`$ is normalizable and consequently falls off faster than $`|z|^{1/2}`$; (b) $`\widehat{\psi }_0(z)`$ is non-normalizable and falls off as a power $`|z|^\alpha `$ ($`\alpha \frac{1}{2}`$); and (c) $`\widehat{\psi }_0(z)`$ is non-normalizable and asymptotes to a constant so the 5D spacetime is asymptotically flat. (Note that this case requires a region where $`A^{\prime \prime }(z)<0`$, and therefore cannot arise as a scalar field domain wall in an otherwise flat background, for example.)

For case (a), $`\widehat{\psi }_0(z)`$ gives rise to the usual 4D Newton potential, as is apparent from (4). In  it was shown that the effects of the continuum modes generically give small corrections. The intuitive reason for this is that in order for $`\widehat{\psi }_0(z)`$ to be normalizable, the tunneling probability for continuum modes of small $`m`$ through the barriers of $`V(z)`$ must vanish as $`m0`$; this implies that $`\psi _m(0)^2`$ must be less singular than $`m^1`$. Hence the integral over the continuum modes in (5) always yields higher order corrections in $`r^1`$ relative to the leading $`r^1`$ piece.

For case (b) there are only continuum modes. In these cases $`\psi _m(0)^2m^{2\alpha }`$, for small $`m`$, and so the potential (5) behaves as $`r^{2\alpha 2}`$ . Consequently there is no effective 4D gravity. Note that in this case gravity does not necessarily appear 5D at large distances, either.

Case (c) is the main focus of interest. As in case (b) there is only a continuum contribution to the potential (5); however, we will argue that under favourable conditions an effective 4D Newton’s Law is recovered as in  in some intermediate regime $`r_1rr_2`$, where the parameters $`r_{1,2}`$ depend on details of the warp factor $`A(z)`$. This region of “quasi-localization” of gravity is intimately connected with the quasi-bound-state $`\widehat{\psi }_0(z)`$ which causes a resonance in $`\psi _m(0)^2`$ at $`m=0`$. The new long distance scale $`r_2`$ depends upon the width of the resonance: the narrower the resonance the larger the scale $`r_2`$. In the limit in which the width goes to zero $`\widehat{\psi }_0(z)`$ becomes normalizable, $`r_2\mathrm{}`$ and 4D gravity is obtained for all distance scales $`r_1`$, as in the original RS model .

In order to prime our intuition it is useful to consider the “volcano box” potential (Figure 1), for which the continuum modes can be calculated exactly.

The solution for a symmetric continuum wavefunction has the form:

$$\psi _m(z)=\frac{1}{\sqrt{c(m)^2+d(m)^2}}\{\begin{array}{cc}\mathrm{cos}k_1z\hfill & |z|z_1\hfill \\ a(m)e^{k_2(zz_1)}+b(m)e^{k_2(zz_1)}\hfill & z_1|z|z_2\hfill \\ c(m)\mathrm{cos}k_3(zz_2)+d(m)\mathrm{sin}k_3(zz_2)\hfill & |z|z_2,\hfill \end{array}$$
(6)
where

$$k_1=\sqrt{V_1+m^2},k_2=\sqrt{V_2m^2},k_3=\sqrt{m^2},$$
(7)
and

$$\begin{array}{cc}\hfill c(m)& =\mathrm{cos}k_1z_1\mathrm{cosh}k_2(z_2z_1)\frac{k_1}{k_2}\mathrm{sin}k_1z_1\mathrm{sinh}k_2(z_2z_1),\hfill \\ \hfill d(m)& =\frac{k_2}{k_3}\mathrm{cos}k_1z_1\mathrm{sinh}k_2(z_2z_1)\frac{k_1}{k_3}\mathrm{sin}k_1z_1\mathrm{cosh}k_2(z_2z_1).\hfill \end{array}$$
(8)
In order to have a quasi-bound-state at $`m=0`$ we require a solution of the form (6) with $`\psi _0(z)=c(0)`$ for $`|z|z_2`$. This happens when the parameters of the potential satisfy

$$\sqrt{V_2}\mathrm{tanh}\sqrt{V_2}(z_2z_1)=\sqrt{V_1}\mathrm{tan}\sqrt{V_1}z_1.$$
(9)
The quantity we are after is

$$\psi _m(0)^2=\frac{1}{c(m)^2+d(m)^2}.$$
(10)
For small $`m`$, we can expand $`c(m)`$ and $`d(m)`$ in powers of $`m`$. The important point is that $`c(m)`$ has an expansion in even powers of $`m`$ while $`d(m)`$ has a expansion in odd powers of $`m`$:

$$c(m)=c_0+c_1m^2+\mathrm{},d(m)=d_0m+d_1m^3+\mathrm{}.$$
(11)
Hence for small $`m`$, $`\psi _m(0)^2`$ has the Breit-Wigner form indicative of a resonance

$$\psi _m(0)^2=\frac{𝒜}{m^2+\mathrm{\Delta }m^2}+𝒪(m^4),$$
(12)
where the width of the resonance is $`\mathrm{\Delta }m=|c_0|/\sqrt{d_0^2+2c_0c_1}`$. The width depends in a complicated way on the parameters of the potential $`V(z)`$. A narrow resonance can be achieved by having $`c_0`$ small, which requires $`\sqrt{V_2}(z_2z_1)1`$. In this limit the width is approximately

$$\mathrm{\Delta }m\frac{8}{(1+V_2/V_1)(z_1+1/\sqrt{V_2})}e^{2\sqrt{V_2}(z_2z_1)}.$$
(13)
Intuitively, the behaviour of (13) can be deduced by the following reasoning: in order to get a narrow resonance we require that the tunneling probability for modes of small $`m`$ through the barriers of the potential be very small, which is achieved by having $`\sqrt{V_2}(z_2z_1)1`$. In fact the exponential factor in (13) can be deduced from a simple WKB analysis. In this approximation the tunneling probability for the eigenfunction $`\psi _m(z)`$ is

$$T(m)\mathrm{exp}\left[2_{z_1}^{z_2}𝑑z\sqrt{V(z)m^2}\right]=\mathrm{exp}\left[2\sqrt{V_2m^2}(z_2z_1)\right].$$
(14)
The width of the quasi-bound-state is $`\mathrm{\Delta }mT(0),`$ giving the exponential dependence in (13).

We expect the existence of a resonance to be generic, since it is caused by the non-normalizable mode $`\widehat{\psi }_0(z)`$. Let us suppose that the resonance is sufficiently narrow that we can approximate $`\psi _m(0)^2`$ by

$$\psi _m(0)^2=\frac{𝒜}{m^2+\mathrm{\Delta }m^2}+f(m),$$
(15)
where $`f(m)`$ is some underlying function which rises from 0 to 1 as the energy goes from 0 to just over the height of the barrier. If we assume that $`f(m)m^\beta `$ ($`\beta >0`$), for small $`m`$, the gravitational potential (5) is

$$U(r)=\frac{r_2M_{}^3𝒜}{r}_0^{\mathrm{}}𝑑x\frac{e^{xr/r_2}}{x^2+1}+𝒪(1/r^{\beta +2}),$$
(16)
where $`r_2=1/\mathrm{\Delta }m`$. The contribution from the resonance gives the 4D Newton’s Law for $`rr_2`$ with Newton’s constant

$$G_N=\frac{\pi r_2𝒜}{2M_{}^3},$$
(17)
whereas for $`rr_2`$ the contribution goes as $`1/r^2`$ and so the 5D Newton potential is recovered at very large distances. The rise of the underlying part of the continuum $`f(m)`$ gives the short distance corrections to the 4D Newton’s Law and sets the lower scale $`r_1`$.

The question now is under what conditions is the resonance narrow so that $`r_2`$ can be large. To investigate this we can, following our analysis of the volcano box potential, use the WKB approximation as a guide. The point is that the width of the resonance is proportional to the tunneling probability $`T(m)`$ for the continuum modes of zero energy, through the barrier of the potential, evaluated at $`m=0`$. The WKB approximation gives

$$T(m)\mathrm{exp}\left[2_{z_1}^{\mathrm{}}𝑑z\sqrt{V(z)m^2}\right],$$
(18)
where the integral is over the barrier region of the potential. Notice that the integral in (18) is convergent only when $`V(z)`$ falls off faster than $`z^2`$; precisely the situation for case (c) above. In this case the limit $`T(0)`$ is finite and since we expect $`\mathrm{\Delta }mT(0)`$ the requirement for a narrow resonance is

$$_{z_1}^{\mathrm{}}𝑑z\sqrt{V(z)}1.$$
(19)
In other words, the barriers of the potential must be sufficiently powerful. For case (a) above, on the other hand, the integral diverges and $`T(0)=0`$, as expected since $`\widehat{\psi }_0(z)`$ is normalizable.

Now we consider some concrete examples. To recap what we need in order to get quasi-localized gravity on the brane is some geometry which is asymptotically flat in the transverse dimension. This was realized in the simple model of GRS  by patching together 5D AdS space onto 5D Minkowski space at some point $`z=z_0`$:

$$A(z)=\{\begin{array}{cc}2\mathrm{log}(k|z|+1)\hfill & |z|z_0\hfill \\ 2\mathrm{log}(k|z_0|+1)\hfill & |z|z_0.\hfill \end{array}$$
(20)
This case can be exactly solved and one can show that quasi-localization occurs for $`k^1rk^2z_0^3`$. This allows us to check our simple WKB analysis; in this case for $`z_0k^1`$ the WKB integral is $`𝑑z\sqrt{V(z)}\sqrt{15/4}\mathrm{ln}z_0`$ giving $`r_2z_0^{\sqrt{15}}`$ (to compare with the exact behaviour $`r_2z_0^3`$). In the GRS scenario the patching of AdS to flat space requires the existence of new branes at $`z_0`$. However, this is not a necessary feature and one can invent scenarios which smoothly interpolate between the AdS geometry and flat space; for instance by taking

$$A(z)=2\mathrm{log}\left(\frac{1}{k|z|+1}+a\right),$$
(21)
which can also be smoothed at $`z=0`$ by taking, for example,

$$A(z)=\mathrm{log}\left(\frac{1}{k^2z^2+1}+a^2\right).$$
(22)
In this case, the crossover occurs smoothly at $`z_01/(ka)`$. For this example, it is a simple matter to numerically solve the differential equation (2) and find $`\psi _m(0)^2`$ as a function of $`m`$. Figure 2 illustrates this function for some values of the parameters giving rise to a fairly broad resonance so that the rise of $`\psi _m(0)^2`$ to 1 can also be seen. For small $`a`$ the resonance becomes very narrow and one can verify numerically that it approximates the form in (15).

Figure 3 shows a close-up of the resonance to illustrate the Breit-Wigner form.

We can find the height of the resonance from our knowledge of $`\widehat{\psi }_0(z)`$:

$$\frac{𝒜}{\mathrm{\Delta }m^2}=\frac{2}{\pi }\widehat{\psi }_0(0)^2=e^{3(A(0)A(\mathrm{}))/2}\frac{2}{\pi }a^3.$$
(23)
In addition, when the resonance is very narrow, we can find the width to leading order in $`a`$ by using the fact that, as $`a0`$, $`\widehat{\psi }_0(z)`$ becomes normalizable and the effect of the resonance should appproximate a delta function. This gives

$$\frac{𝒜}{\mathrm{\Delta }m}=\frac{4k}{\pi }.$$
(24)
Hence, $`\mathrm{\Delta }m2ka^32k^2z_0^3=r_2^1`$.

We have introduced the idea of “quasi localization” of gravity on the brane as the general notion lying behind the scenario of GRS  whereby over a large range of intermediate distance scales the gravitational potential is, to a good accuracy, Newtonian. We argued that the relevant geometry for quasi localization is when the transverse geometry becomes asymptotically flat. The new large scale above which the gravitational potential becomes five-dimensional depends on the scale at which the crossover to flat space occurs. Physically we can describe the onset of 5D gravity on the brane at this new scale by saying that the effective 4D graviton is unstable and decays into the KK continuum. Obviously in order to be phenomenologically viable the lifetime must be very long. Although we have only considered the case with one extra dimension we expect that the same phenomenon to occur with any number of extra dimensions. In the same vein, we note that quasi-localization may be of relevance in string theory where the transverse geometry to a $`p`$-brane soliton is indeed asymptotically flat.

We would like to thank Yuri Shirman for useful discussions and fruitful collaboration. C.C. is an Oppenheimer Fellow at the Los Alamos National Laboratory. C.C., J.E. and T.J.H. are supported by the US Department of Energy under contract W-7405-ENG-36.

Note added: After this paper has been completed a revised version of Ref.  (as well as Ref. ) appeared, which also independently noted the presence of the resonance in these theories.


text:
# From Historical Perspectives to Some Modern Possibilities

## References

Aller, L. H. 1956, “Gaseous Nebulae” (London: Chapman and Hall)

Aller, L. H., Ufford, C. W. & Van Vleck, J. H. 1949, ApJ, 109, 42

Aller, L. H., Hyung, S. & Feibelman, W. A. 1999, Proc. Nat. Acad. Sci., 96, 5371

Campbell, W. A. & Moore, J. H. 1918, Lick Obs. Publ., 13, 82.

Curtis, H. D. 1918, Lick Obs. Publ., 13, 57

Feibelman, W. A., Aller, L. H., Keyes, C. D. & Czyzak, S. J. 1985, Proc. Nat. Acad. Sci., 84, 2292.

Green, W. K. 1917, Lick Obs. Bull., 9, 72

Keenan, F. P., Aller, L. H., Bell, K. L. Hyung, S., McKenna, F., & Ramsbottom, C. 1996, MNRAS 287. 1073

Keenan, F. C., McKenna, F., Bell, K. L., Ramsbottom, C., Wickstead, A., Aller, L. H., & Hyung, S. 1997 ApJ, 487, 457

Keenan, F. C., Aller, L. H., Bell, K., Espey, Feibelman W. A., Hyung, S., McKenna, F., & Ramsbottom D. 1998, MNRAS, 295, 688

Keenan, D. C., Aller, L. H., Bell, K. L., Hyung, S., Crawford, F. L., Feibelman, W. A., McKenna, F., & McLaughlin, B. M. 1999, MNRAS, 304, 27

Menzel, D. H., & Aller, L. H. 1941, ApJ, 93, 230

Minkowski, R. 1964, PASP, 76. 197

Minkowski, R., & Osterbrock, D. E. 1960, ApJ, 131, 547

Perrine, C. D., 1929, Astr. Nach., 237, 89

Pradhan, A. K. 1976, MNRAS, 177, 31

Rowlands, N. Houck, J. U., & Herter, E. 1994, ApJ, 427, 867

Sabbadin, F., & Hamazaoglu, G., 1982, A&A, 109, 134

Seaton M. J. 1954a MNRAS 114, 154

Seaton M. J. 1954b, Annales d’Astrophysique, 17, 396

Seaton, M. J., & Osterbrock, D.E. 1957, ApJ, 125, 66

Wilson, O. C. 1950, ApJ, 111, 279

Zanstra, H. 1931, ZAp, 2, 329


text:
# A-Tract Induced DNA Bending is a Local Non-Electrostatic Effect

## Results and Discussion

Although the macroscopic curvature of DNA induced by adenine-tracts (A-tracts) was discovered almost two decades ago  structural basis for this phenomenon remains unclear. A few models considered originally  suggested that it is caused by intrinsic conformational preferences of certain sequences, but all these and similar theories failed to explain experimental data obtained later.  Calculations show that the B-DNA duplex is mechanically anisotropic,  that bending towards minor grooves of some A-tracts is strongly facilitated,  and that the macroscopic curvature becomes energetically preferable once the characteristic A-tract structure is maintained by freezing or imposing constraints.  However, the static curvature never appears spontaneously in calculations unbiased a priori and these results leave all doors open for the possible physical origin of the effect. In the recent years the main attention has been shifted to specific interactions between DNA and solvent counterions that can bend the double helix by specifically neutralizing some phosphate groups.  The possibility of such mechanism is often evident in protein-DNA complexes, and it has also been demonstrated by direct chemical modification of a duplex DNA.  In the case of the free DNA in solution, however, the available experimental observations are controversial.  Molecular dynamics simulations of a B-DNA in an explicit counterion shell could neither confirm nor disprove this hypothesis.  Here we report the first example where stable static curvature emerges spontaneously in molecular dynamics simulations. Its direction is in striking agreement with expectations based upon experimental data. However, we use a minimal B-DNA model without counterions, which strongly suggests that they hardly play a key role in this effect.

Figure 1 exhibits results of a 10 ns simulation of dynamics of a 25-mer B-DNA fragment including three A-tracts separated by one helical turn. This sequence has been constructed after many preliminary tests with shorter sequence motives. Our general strategy came out from the following considerations. Although the A-tract sequences that induce the strongest bends are known from experiments, probably not all of them would work in simulations. There are natural limitations, such as the precision of the model, and, in addition, the limited duration of trajectories may be insufficient for some A-tracts to adopt their specific conformation. Also, we can study only short DNA fragments, therefore, it is preferable to place A-tracts at both ends in order to maximize the possible bend. There is, however, little experimental evidence of static curvature in short DNA fragments, and one may well expect the specific A-tract structure to be unstable near the ends. That is why we did not simply take the strongest experimental “benders”, but looked for sequence motives that in calculations readily adopt the characteristic local structure, with a narrow minor groove profile and high propeller twist, both in the middle and near the ends of the duplex. The complementary duplex $`\mathrm{AAAATAGGCTATTTTAGGCTATTTT}`$ has been constructed by repeating and inverting one such motive.

The upper trace in plate (a) shows the time dependence of rmsd from the canonical B-DNA model. It fluctuates below 4 Å sometimes falling down to 2 Å, which is very low for the double helix of this length indicating that all helical parameters are well within the range of the B-DNA family. The lower surface plot shows the time evolution of the minor DNA groove. The surface is formed by 75 ps time-averaged successive minor groove profiles, with that on the front face corresponding to the final DNA conformation. The groove width is evaluated by using space traces of C5’ atoms as described elsewhere . Its value is given in angströms and the corresponding canonical B-DNA level of 7.7 Å is marked by the straight dotted lines on the faces of the box. It is seen that the overall groove shape has established after 2 ns and remained stable later, with noticeable local fluctuations. In all A-tracts the groove strongly narrows towards 3’ ends and widens significantly at the boundaries. There are two less significant relative narrowings inside non A-tract sequences as well.

Dynamics of $`\mathrm{B}_\mathrm{I}\mathrm{B}_{\mathrm{II}}`$ backbone transitions are shown in plate (b). The BI and BII conformations are distinguished by the values of two consecutive backbone torsions, $`\epsilon `$ and $`\zeta `$. In a transition they change concertedly from (t,g-) to (g-,t). The difference $`\zeta \epsilon `$ is, therefore, positive in BI state and negative in BII, and it is used in Fig. (d) as a monitoring indicator, with the corresponding gray scale levels shown on the right. Each base pair step is characterized by a column consisting of two sub-columns, with the left sub-columns referring to the sequence written at the top in 5’-3’ direction from left to right. The right sub-columns refer to the complementary sequence shown at the bottom. It is seen that, in A-tracts, the BII conformation is preferably found in ApA steps and that $`\mathrm{B}_\mathrm{I}\mathrm{B}_{\mathrm{II}}`$ transitions in neighboring steps often occur concertedly so that along a single A-strand $`\mathrm{B}_\mathrm{I}`$ and $`\mathrm{B}_{\mathrm{II}}`$ conformations tend to alternate. The pattern of these transitions reveals rather slow dynamics and suggests that MD trajectories in the 10 ns time scale are still not long enough to sample all relevant conformations. Note, for instance, a very stable $`\mathrm{B}_{\mathrm{II}}`$ conformation in both strands at one of the GpG steps.

Plate (c) shows the time evolution of the overall shape of the helical axis. The optimal curved axes of all DNA conformations saved during dynamics were rotated with the two ends fixed at the OZ axis to put the middle point at the OX axis. The axis is next characterized by two perpendicular projections labeled X and Y. Any time section of the surfaces shown in the figure gives the corresponding axis projection averaged over a time window of 75 ps. The horizontal deviation is given in angströms and, for clarity, its relative scale is two times increased with respect to the true DNA length. Shown on the right are two perpendicular views of the last one-nanosecond-average conformation. Its orientation is chosen to correspond approximately that of the helical axis in the surface plots.

It is seen that the molecule maintained a planar bent shape during a substantial part of the trajectory, and that at the end the bending plane was passing through the three A-tracts. The X-surface clearly shows an increase in bending during the second half of the trajectory. In the perpendicular Y-projection the helical axis is locally wound, but straight on average. The fluctuating pattern in Y-projection sometimes reveals two local maxima between A-tracts, which corresponds to two independent bends with slightly divergent directions. One may note also that there were at least two relatively long periods when the axis was almost straight, namely, around 3 ns and during the fifth nanosecond. At the same time, straightening of only one of the two bending points is a more frequent event observed several times in the surface plots.

Finally, plate (d) shows the time fluctuations of the bending direction and angle. The bending direction is characterized by the angle between the X-projection plane in plate (c) and the $`xz`$ plane of the local DNA coordinate frame constructed in the center of the duplex. According to the Cambridge convention  the local $`x`$ direction points to the major DNA groove along the short axis of the base-pair, while the local $`z`$ axis direction is adjacent to the optimal helicoidal axis. Thus, a zero angle between the two planes corresponds to the overall bend to the minor groove exactly at the central base pair. In both plots, short time scale fluctuations are smoothed by averaging with a window of 15 ps. The total angle measured between the opposite axis ends fluctuates around 10-15 in the least bent states and raises to average 40-50 during periods of strong bending. The maximal instantaneous bend of 58 was observed at around 8 ns.

The bending direction was much more stable during the last few nanoseconds, however, it fluctuated at a roughly constant value of 50starting from the second nanosecond. This value means that the center of the observed planar bend is shifted by approximately two steps from the middle base pair so that its preferred direction is to the minor groove at the two ATT triplets, which is well distinguished in plate (c) as well, and corresponds to the local minima in the minor groove profiles in plate (a). During the periods when the molecule straightened the bending direction strongly fluctuates. This effect is due to the fact that when the axis becomes straight the bending plane is not defined, which in our case appears when the central point of the curved axis passes close to the line between its ends. It is very interesting, however, that after the straightening, the bending is resumed in approximately the same direction.

Figure 2 exhibits similar data for another 10 ns trajectory of the same DNA fragment, computed in order to check reproducibility of the results. A straight DNA conformation was taken from the initial phase of the previous trajectory, energy minimized, and restarted with random initial velocities. It shows surprisingly similar results as regards the bending direction and dynamics in spite of a somewhat different minor groove profile and significantly different distribution of $`\mathrm{B}_\mathrm{I}`$ and $`\mathrm{B}_{\mathrm{II}}`$ conformers along the backbone. Note that in this case the helical axis was initially S-shaped in X-projection, with one of the A-tracts exhibiting a completely opposite bending direction. Fluctuations of the bending direction are reduced and are similar to the final part of the first trajectory, which apparently results from the additional re-equilibration. In this case the maximal instantaneous bend of 71 was observed at around 4 ns.

Comparison of traces in plates (a) and (d) in Figs. 1 and 2 clearly shows that large scale slow fluctuations of rmsd are caused by bending. The rmsd drops down to 2 Å when the duplex is straight and raises beyond 6 Å in strongly bent conformations. In both trajectories the molecule experienced many temporary transitions to straight conformations which usually are very short living. These observations suggest that the bent state is relatively more stable than the straight one and, therefore, the observed behavior corresponds to static curvature. In conformations averaged over successive one nanosecond intervals the overall bending angle is 35-45 except for a few periods in the first trajectory. Figure 3 shows a snapshot from around 8.5 ns of the second trajectory where the rmsd from the straight canonical B-DNA reached its maximum of 6.5 Å. The strong smooth bent towards the minor grooves of the three A-tracts is evident, with the overall bending angle around 61.

All transformations exhibited in Figs. 1 and 2 are isoenergetic, with the total energy fluctuating around the same level established during the first nanosecond already, and the same is true for the average helicoidal parameters. Plates (b), however, indicate that there are much slower motions in the system, and this observation precludes any conclusions concerning the global stability of the observed conformations. Moreover, we have computed yet another trajectory for the same molecule starting from the canonical A-DNA form. During 10 ns it converged to a similarly good B-DNA structure with the same average total energy, but the bending pattern was not reproduced. It appears, therefore, that the conformational space is divided into distinct domains, with transitions between them probably occurring in much longer time scales. However, the very fact that the stable curvature in good agreement with experimental data emerges in trajectories starting from a featureless straight canonical B-DNA conformation strongly suggests that the true molecular mechanism of the A-tract induced bending is reproduced. Therefore, it cannot depend upon the components discarded in our calculations, notably, specific interactions with solvent counterions and long-range electrostatic effects.

We are not yet ready to present a detailed molecular mechanism responsible for the observed curvature because even in this relatively small system it is difficult to distinguish the cause and the consequences. We believe, however, that all sorts of bending of the double helical DNA, including that produced by ligands and that due to intrinsic sequence effects, have its limited, but high flexibility as a common origin. Its own conformational energy has the global minimum in a straight form, but this minimum is very broad and flat, and DNA responds by distinguishable bending to even small perturbations. The results reported here prove that in the case of A-tracts these perturbations are produced by DNA-water interactions in the minor groove. Neither long range phosphate repulsion nor counterions are essential. The curvature is certainly connected with the specific A-tract structure and modulations of the minor groove width, but it does not seem to be strictly bound to them. In dynamics, conformations, both smoothly bent and kinked at the two insertions between the A-tracts, are observed periodically. Note also, that the minor groove profile somewhat differs between the two trajectories and that it does not change when the molecule straightens. We strongly believe, however, the experimental data already available will finally allow one to solve this problem by theoretical means, including the approach described here, and we continue these attempts.

## Methods

Molecular dynamics simulations have been performed with the internal coordinate method (ICMD)  including special technique for flexible sugar rings . The so-called “minimal B-DNA” model was used  which consists of a double helix with the minor groove filled with explicit water. Unlike the more widely used models, it does not involve explicit counterions and damps long range electrostatic interactions in a semi-empirical way by using distance scaling of the electrostatic constant and reduction of phosphate charges. The DNA model was same as in earlier reports,  namely, all torsions were free as well as bond angles centered at sugar atoms, while other bonds and angles were fixed, and the bases held rigid. AMBER94  force field and atom parameters were used with TIP3P water  and no cut off schemes. With a time step of 10 fs, these simulation conditions require around 75 hours of cpu per nanosecond on a Pentium II-200 microprocessor.

The initial conformations were prepared by vacuum energy minimization starting from the fiber B-DNA model constructed from the published atom coordinates.  The subsequent hydration protocol to fill up the minor groove  normally adds around 16 water molecules per base pair. The heating and equilibration protocols were same as before . During the runs, after every 200 ps, water positions were checked in order to identify those penetrating into the major groove and those completely separated. These molecules, if found, were removed and next re-introduced in simulations by putting them with zero velocities at random positions around the hydrated duplex, so that they could readily re-join the core system. This procedure assures stable conditions, notably, a constant number of molecules in the minor groove hydration cloud and the absence of water in the major groove, which is necessary for fast sampling . The interval of 200 ps between the checks is small enough to assure that on average less then one molecule is repositioned and, therefore, the perturbation introduced is considered negligible.

## Acknowledgements

I thank R. Lavery for useful discussions as well as critical comments and suggestions to the paper.

## Appendix

This section contains comments from anonymous referees of peer-review journals where the manuscript has been considered for publication, but rejected.

### A Journal of Molecular Biology

#### 1 First referee

Dr. Mazur describes molecular dynamics simulations where a correct static curvature of DNA with phased A-tracts emerges spontaneously in conditions where any role of counterions or long range electrostatic effects can be excluded.

I have several problems with this manuscript:

1) The observed curvature is dependent on the starting model. In fact the manuscript uses the phrase ‘stable static curvature’ incorrectly to describe what is probably a trapped metastable state. The observed curve is neither stable nor static.

2) The choice of DNA sequence seems to be biased toward that which gives an altered structure in simulations, ad is not that which gives the most pronounced bend in solution. I would suggest a comparison of (CAAAATTTTTG)n and (CTTTTAAAAG)n.

3) The result is not consistent with solution results. See for example:

Prodin, F., Cocchione, S., Savino, M., & Tuffillaro, A. “Different Interactions of Spermine With a Curved and a Normal DNA Duplex - (Ca(4)T(4)G)(N) and (Ct(4)a(4)G)(N) - Gel -Electrophoresis and Circular-Dichroism Studies” (1992) Biochemistry International 27, 291-901.

Brukner, l, Sucis, S., Dlakic, M., Savic, A., & Pongor, S. “Physiological concentrations of magnesium ions induces a strong macroscopic curvature in GGGCCC - containing DNA” (1994) J. Mol. Biol. 236, 26-32.

Diekmann, S., & Wang, J. C. “On the sequence determinants and flexibility of the kinetoplast DNA fragment with abnormal gel electrophoretic mobilities” (1985) J. Mol. Biol. 186, 1-11.

Llaudnon, C. H., & Griffith, J. D. “Cationic metals promote sequence-directed DNA bending” (1987) Biochemistry 26, 3759-3762.

4) The result is not consistent with other simulations. See for example:

Feig, M., & Pettitt, B. M. “Sodium and Chlorine ions as part of the DNA solvation shell” (1999) Biophys. J. 77, 1769-81.

5) The results should be given by objective statistical descriptions rather than a series of spot examples, as in “sometimes reveals two independent bends”.

#### 2 Second referee

This manuscript describes the modeling of a 25-residue DNA duplex using molecular dynamics simulations. The DNA sequence in question contains 3 A/T tracts arranged in-phase with the helix screw and thus is expected to manifest intrinsic bending. Unlike previous MD studies of intrinsically bent DNA sequences, these calculations omit explicit consideration of the role of counterions. Because recent crystallographic studies of A-tract-like DNA sequence have attributed intrinsic bending to the localization of counterions in the minor groove, the present finding that intrinsic bending occurs in the absence of explicit counterions is important for understanding the underlying basis of A-tract-dependent bending.

Overall, the MD procedure appears sound and the calculations were carried out with obvious care and attention to detail. There are two specific issues raised by this study that should be addressed in revision, however.

1. Although the sequence chosen for this study was based on a canonical, intrinsically-bent motif consisting of three A tracts, it is unclear to what extent intrinsic bending has been experimentally shown for this particular sequence. There are known sequence-context effects that modulate A-tract-dependent bending and thus the author should refer the reader to data in the literature or show experimentally that intrinsic bending of the expected magnitude occurs for this particular sequence. Moreover, one A tract is out-of-phase with respect to the others and it is therefore not clear how this contributes to the overall bend. The author is understandably concerned about end effect with short sequences; this problem can be ameliorated by examining DNA fragments that constrain multiple copies of the chosen motif or by extending the ends of the motif with mixed-sequence DNA.

2. Notwithstanding the authors remark bout separating the cause and the effects with respect to intrinsic bending some comments about the underlying mechanism of bending seem appropriate. It would be particularly useful to know whether average values of any specific conformational variables are unusual or whether strongly bent states are consistent with narrowing of the minor groove within A-tracts, for example.


text:
# Shape anisotropy and Voids

## Abstract

Numerical simulations on a 2-dimensional model system showed that voids are induced primarily due to shape anisotropy in binary mixtures of interacting disks. The results of such a simple model account for the key features seen in a variety of flux experiments using liposomes and biological membranes Sita .

A variety of lipid molecules contribute to the structure and barrier function of biological membranes. Part of each molecule is hydrophilic (head) and part hydrophobic (tail). This amphipathic nature, in the presence of water, leads to a bilayer structure in which hydrophilic head groups have maximum contact with water and hydrophobic tails have minimum contact with water. The process of membrane formation is one of minimizing the free energy and maximizing the stability of the structure Stein .

A major question has been whether the bilayer is to be viewed as an isotropic homogeneous phase or as a heterogeneous phase Lee . If osmotic contraction of the bilayer vesicles leads to an altered hydraulic conductivity (water flux coefficient), one obviously favors a heterogeneous membrane model. Otherwise a homogeneous, isotropic phase model would be adequate, obviating the need to look for a fine structure within the bilayer. Using the erythrocyte as an experimental system ( in which the area of the biconcave cell does not change when it is osmotically expanded to a spherical shape), it was concluded that hydraulic conductivity was stretch independent, i.e., in support of the isotropic model Sita . An alternative way to assess hydraulic conductivity is to use hydrogen peroxide as an analog of water: since many experimental systems have catalase (an enzyme that degrades hydrogen peroxide to molecular oxygen and water) within the vesicle/cell, an assay of this occluded catalase directly permits one to measure the conductivity to exogenous hydrogen peroxide. Under equilibrium conditions of assay, the rate of degradation would be same as the rate of permeation of the peroxide into the vesicle. Thus, one can directly assess the stretch sensitivity of the membrane by osmotic titrations with osmolites, using non-electrolytes like hydrogen peroxide as probes of flux. In the course of these experiments it was found that Sita : (i) among all the lipid combinations tested, the phosphatidylcholine (PC) vesicles and intact erythrocytes, both, did not show a decrease in occluded catalase activity on osmotic compression of the membrane (ii) on the other hand, all other membrane systems, such as peroxisomes, E.coli, macrophages showed stretch (osmotic) sensitivity (iii) so did liposomes made from these cells and organelles (iv) further, when binary mixtures were investigated, only cardiolipin and cerebrosides when added to PC (5 to 10$`\%`$ of PC) conferred stretch sensitivity in liposomes (v) these binary mixtures also exhibited enhanced activation volume(osmotic sensitivity) and diminished activation energy for hydrogen peroxide flux (vi) further, glucose was readily permeable across these membranes of binary mixtures (vii) addition of cholesterol, which is abundantly found in erythrocytes, inhibited the stretch sensitivity to peroxide permeation (viii) evidence was also seen that this diffusion increases with decrease in temperature, i.e., the process has a negative temperature coefficient.

These studies on biological membranes and liposomes, in which composition as well as dynamics considerably vary, prompted us to question as to what constitutes a minimal description to account for variable permeability induced by doping across a liposomal membrane. For instance, cardiolipin enhanced permeation of hydrogen peroxide and molecules as large as glucose Sita . Though possible descriptors could be many, (composition, structure, dynamics in terms of inter and intra molecular potentials), we adopted a bare-bone approach to resolve this complex issue to arrive at a minimal description adequate to account for the observations.

Structural changes in the membrane are best identified by non-interactive molecules and therefore leaks across bilayers are commonly studied using non-electrolytes Stein . The diagnostic for non-specific permeation is size dependence such that, these hydrated solutes intercalate, penetrate and navigate through such interstices, spaces or voids, stochastically or in files to reach the other side of the membrane Lee . In order to capture the diverse features in a parsimonious manner, we restrict to a two dimensional cross-section of a three dimensional system. Such a restriction is reasonable since the probe particle permeating across the membrane at any instant of time experiences the effective cross section rather than the three dimensional obstruction. The permeation across the membrane depends primarily on the availability of free space or voids. Thus the problem reduces to the study of packing of 2-dimensional objects at the first instance. Then one needs to determine which factor(s) determine the appearance and size-distribution of voids in such a 2-dimensional system.

The configuration space of this model system (membrane) is a 2-dimensional box with toroidal boundary conditions. The constituents of this two dimensional box are the circular disks (and/or the rigid combinations of the circular disks as dopants) of unit radii. A circular disk simply represents the hard core scattering cross section, seen by the passing particle (a non-electrolyte which acts as a probe), across the thickness spanned by two lipid molecules, viewed somewhat as cylinders. A typical dopant is two or more circular disks rigidly joined in a prespecified geometry. These constituents are identified by the position coordinates of their centers, the angle made by the major axes with the side of the box ( in case of dopants) and the radii of the circular disks. It is reasonable that the disks (which represent molecules, with long range attractive interaction and hard core repulsion near center, contained within a structure) interact pair wise via Lennard-Jones potential (a measure of the interaction energy), which has the form

$`V_{LJ}(r_{ij})=4ϵ{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=i+1}{\overset{N}{}}}\left(({\displaystyle \frac{\sigma }{r_{ij}}})^{12}({\displaystyle \frac{\sigma }{r_{ij}}})^6\right)`$
where, $`r_{ij}`$ is the distance between the centers of the $`i^{\mathrm{th}}`$ and $`j^{\mathrm{th}}`$ disks, $`\sigma `$ determines the range of hard core part in the potential and the $`ϵ`$ signifies the depth of attractive part. While studying the binary mixtures, we consider different shape anisotropic combinations (impurities or dopants) of $`\kappa `$ number of circular disks. We treat these combinations as one unit. e.g. rodn denotes a single dopant made of $`n`$ unit circular disks rigidly joined one after another in a straight line. The impurities interact with constituent circular disks via potential,

$`V(r_{ij})={\displaystyle \underset{\alpha =1}{\overset{\kappa }{}}}V_{LJ}(r_{i_\alpha j})`$
where, $`r_{i_\alpha j}`$ is the distance between the centers of $`\alpha ^{\mathrm{th}}`$ disk in $`i^{\mathrm{th}}`$ impurity and the $`j^{\mathrm{th}}`$ circular disk, and among themselves interact via

$`V(r_{ij})={\displaystyle \underset{\alpha =1}{\overset{\kappa _1}{}}}{\displaystyle \underset{\beta =1}{\overset{\kappa _2}{}}}V_{LJ}(r_{i_\alpha j_\beta })`$
where, $`r_{i_\alpha j_\beta }`$ is the distance between the centers of $`\alpha ^{\mathrm{th}}`$ disk in $`i^{\mathrm{th}}`$ impurity and the $`\beta ^{\mathrm{th}}`$ disk in $`j^{\mathrm{th}}`$ impurity.

An $`r`$-void is defined as a closed area in a membrane devoid of disks or impurities, and sufficient to accommodate a circular disk of radius $`r`$ Gauri . Clearly, larger voids also accommodate smaller probes, i.e., an $`r`$-void is also an $`r^{}`$-void if $`r^{}<r`$. Similarly, the voids for the particle of size zero are the voids defined in the conventional sense, i.e., a measure of the net space unoccupied by the disks.

Equilibrium configurations of the model system are obtained by a Monte Carlo method (using the Metropolis algorithm) starting with a random placement of the disks Gauri ; Binder 111The equilibrium configurations thus obtained are further confirmed by simulated annealing NR .. The box was filled with disks such that they occupy 70% area of the box,i.e.,loosely packed to facilitate the formation of voids. The temperature parameter, $`T`$, was so chosen that the quantity $`k_BT<4ϵ`$. This ensured an approximate hexagonal arrangement of the disks and the presence of very few large voids in the absence of dopants. (Fig. 1a; far too less r-voids for $`r0.5`$222It may be recalled that glucose offers approximately half the radius of the PC cross section, yielding a relevant definition for a larger void of interest.). Similar numerical simulations are performed on the model system with dopants. The number of dopants is chosen to be 10% Sita  of the number of circular disks with a constraint that the total occupied area of the box is still 70% as the focus was on the redistribution of void sizes. Fig. 1b illustrates the formation of larger voids in the vicinity of the rod2 impurities.

The variation in the number of $`r`$-voids as a function of the size of the permeating particle(using the digitization algorithm described in Gauri ) is shown in Fig. 2. When only circular disks are present (dotted curve), hardy any large $`r`$-voids are seen. When mixed with the anisotropic impurities, say, rod2, a distinct increase in the number of large $`r`$-voids is seen with appropriate redistribution of smaller $`r`$-voids (solid curve). This result is consistent with the unexpected permeation of large molecules such as glucose through the doped membrane, observed experimentally Sita .The difference curve showed the formation of a significant extent ($``$ 30%) of $`r`$-voids of size $`0.5`$ and above.

Is the induction of large voids due to the anisotropy in potential of the impurities and, should the large voids form around the rods, the centers of anisotropy? Firstly, we carried out simulations with large circular disks in place of rod2 as impurities. The radius of large disk was chosen in such a way that the area occupied by each of the large disk is same as that of a rod2. Fig. 3 shows the result of such simulations. The curve (a) in Fig. 3 represents the difference curve of $`r`$-void distribution of pure membrane and that of membrane doped with rod2 impurity. The curve (b) represents the same when the membrane is doped with large circular disks. It can be clearly seen that the number of larger $`r`$-voids is always less in the latter case, thus confirming the role of shape anisotropy in the induction of large $`r`$-voids. Further, simulations are carried out with rod2 of smaller size and rod4 type impurities. The curves (c), (d) in Fig. 3 respectively shows the corresponding difference curves. It shows an interesting feature that the peak of the difference curve shifts with the change in the type of anisotropy. This suggests a possible way of constructing membrane with selective permeability properties. The simulation regime adopted here, limited the exploration of ternary mixtures in yielding statistically significant results on transport. However, by using rod2 of 0.5 size (Fig.3, curve (c)) (an oval approximation of the small dopant cholesterol), we could demonstrate a shift in void sizes to left in binary mixtures, consistent with our experimental results in ternary mixtures with cholesterol Sita .

Further, we considered rodn type impurities. Fig. 4 represents the relation between the length dimension of rodn and the number of $`r`$-voids (for r=0.55). The anisotropy in the potential of rodn increases with $`n`$, such that the number of large $`r`$-voids should increase with $`n`$. Fig. 4 indeed shows a jump when the rod2 impurities are added and afterwards, it shows a slow and almost linear increase with increase in $`n`$.

Since dopants induced voids, their influence is most likely to be seen in their own vicinity, enhancing the “local transport”. As the dopants exhibit different potential in different directions, certain positions of the constituents are preferred from the point of view of energy minimization, eventually giving rise to voids in the vicinity of impurities. To verify this, we calculate the local permeation probability for particle of size $`r`$, which is a ratio of the area of $`r`$-voids and the area of the local neighborhood. Fig. 5 shows the local permeation probability around ten randomly chosen impurities and ten randomly chosen circular disks. The higher local permeation probability is seen to be associated with rod2.

The model is realistic in that, one can compute elastic properties (surface tension like attribute) by stretching the membrane from one side and computing the energy change, which yielded a change of $``$ 13 dyn/cm which is of the same order as the observed surface tension and changes there of in bilayers Jan . The model is limited in relation to an investigation of temperature effects which requires incorporation of multiple time scales. The model is really general because it not restricted by the parameter space of the components and therefore it is extendible to a variety of phenomena including transport in weakly bound granular media, voids in polymers, modeling of zeolites which may act like a sponge absorbing only desired species.

In summary, we proposed a two dimensional computational model system comprising a mixture of objects interacting via Lennard-Jones potential to explain anomalous permeation seen in bilayers. The significance of these observations on permeation, in what is essentially a granular medium (with long range attraction as an added feature), relates to development of large voids seeded by impurities. Unlike shape anisotropy, the change in composition (via a change in $`\sigma `$), lateral pressure/density, and presumably temperature (though dynamic simulations would need to be done) would all simply produce voids within the bounds of a hexagonal array. Even among various dopants in shapes, X, L, Y, Z, T symmetrical or otherwise, all other factors being equal, size per size, rodn produce the largest voids Gauri .

Thus the largest $`r`$-voids induced by anisotropy yield biological observations of interest, without manifesting as ordered geometric structures. In this sense these $`r`$-voids differ from the conventional results obtained in studies on granular media (which also do not usually incorporate long range interactions) Granular . It is increasingly becoming clear that voids play a pivotal role in relating the dynamics of biopolymers to specific functional states Raj .

We thank C. N. Madhavarao for discussions. The author (GRP) is grateful to CSIR (India) for fellowship and (ADG) is grateful to DBT (India) for research grant.

FIGURE CAPTIONS

Typical equilibrium configurations for interacting disks. The parameter, $`\sigma `$ in Lennard-Jones potential is 2 units. The size of the box is 50$`\times `$50 (a) A pure membrane with 556 unit circular disks. (b) A doped membrane with 464 circular disks and 46 rod2(shown in gray).

Distribution of $`r`$-voids in two different configurations. The main graph shows the number of $`r`$-voids as a function of the relative size ($`r`$) of the probe particle. The vertical bars represent the error margins at the corresponding points. the dotted curve gives the distribution in pure membrane while the solid curve shows the same in a membrane doped with rod2 (10:1). Difference curve clearly demonstrates the presence of large voids in the doped membrane.

Difference curves of distribution of $`r`$-voids. Curves are obtained by treating the void distribution for pure membrane as the base. Difference curves for membranes(a) doped with rod2 (b) doped with large circular disks (of radii $`\sqrt{2}`$) which occupy the same area as that of rod2. These induce smaller number of large voids as compared to rod2. (c) doped with small rods. Irrespective of their small size, they induce large voids, but the peak shifts towards left. (d) doped with rod4. Significantly large voids are induced. The peak shifts towards the right.

Dependence of the number of $`r`$-voids on the length of the rod-shaped impurities. The graph shows a steady increase in the number of $`r`$-voids (for $`r=0.55`$) with $`n`$. The first expected large jump in the number of voids because of the shape anisotropy is seen clearly when the configuration consists of molecules in the shape of circular disks and rod2.

Local permeation probability in a doped model system. The points show the local permeation probability around ten randomly chosen unit disks and ten randomly chosen rod2s. Further, as a guide line, averages are shown by the heights of the boxes, clearly indicating significantly more permeation in the neighborhood of rod2.