text:
# On the temperature dependence of 2D “metallic” conductivity in Si inversion layers at intermediate temperatures

\[

## Abstract

We show that the recent experimental claim {Pudalov et al. Phys. Rev. Lett.91, 126403 (2003) } of observing “interaction effects in the conductivity of Si inversion layers at intermediate temperatures” is incorrect and misleading. In particular, the temperature dependent conductivity $`\sigma `$, in contrast to the resistivity (which is what is shown in the paper), does not have a linear temperature regime, rendering the extraction of the slope $`d\sigma /dT`$ completely arbitrary. We also show that, at least for higher densities, the standard semiclassical transport theory, which includes realistic disorder effects such as scattering by screened charged impurity and surface roughness, gives essentially quantitative agreement with the experimental data.

PACS Number : 71.30.+h; 71.27.+a; 73.40.Qv

\]

In a recent Letter , here after referred to as I, Pudalov et al. have presented experimental results on the temperature dependent resistivity $`\rho (T)`$ of Si inversion layers comparing the data to the so-called interaction theory  of Zala et al. Although claims of “rigorous experimental test” and “excellent agreement” are made rather uncritically in I, the purpose of the current Comment is to point out the misleading and essentially incorrect nature of the main claims in I. In particular, we point out that (1) the comparison to the theory of ref.  carried out in I is inappropriate and arbitrary; and (2) the well-established semiclassical Boltzmann transport theory  for Si inversion layers using screened charged impurity scattering and interface roughness scattering provides quantitative agreement with much of the higher density data ($`n5\times 10^{11}cm^2`$) presented in I with the agreement becoming qualitative at lower densities.

The interaction theory , with which Pudalov et al. compare their data, predicts a linear T-dependence of the conductivity $`\sigma (T)/\sigma _0=1+F\frac{T}{T_F}`$ in the intermediate temperature ballistic regime where $`\sigma _0`$ is the $`T=0`$ “Drude” conductivity and $`T_F`$ the Fermi temperature. The slope $`FF(r_s)`$, which depends on the dimensionless $`r_s`$ ($`n^{1/2}`$) parameter, cannot be calculated within the theory and is predicted to change its sign at high density (which has not yet been observed experimentally). The absolute necessary condition (but, by no means sufficient) for verifying the theory of ref.  is obviously an observation of a linear temperature dependent conductivity over a reasonable temperature range. This is particularly so in view of the fact that the interaction theory  can only predict the qualitative leading order temperature dependence, but not the quantitative magnitude of the slope $`F`$ for the $`\sigma (T)`$ curves. A rudimentary analysis of the data presented in I clearly shows that $`\sigma (T)`$ results of I are nonlinear, and thus do not satisfy the minimal necessary condition needed for a comparison with the interaction theory. The authors of I have taken the misleading step of presenting their data for the resistivity $`\rho (T)`$, rather than the conductivity $`\sigma (T)=[\rho (T)]^1`$. The $`\rho (T)`$ data in I superficially appear relatively more linear than the corresponding $`\sigma (T)`$ results, and quite trivially even a linear $`\rho (T)`$ is not equivalent to a linear $`\sigma (T)`$ unless the temperature correction is small, which is manifestly not the case in I (i.e. the temperature dependence in not weak except at the highest densities where the Boltzmann theory gives quantitatively accurate results).

We show in Fig. 1(a) our best estimate for $`d\sigma /dT`$ extracted numerically for the data of I, and it is obvious that a linear-T approximation to $`\sigma (T)`$, i.e. a constant $`d\sigma /dT`$ over any reasonable range of temperature, simply does not apply anywhere in the data making the whole exercise of the comparison to the interaction theory completely arbitrary. In Fig. 1(b) we show an essentially parameter-free quantitative comparison between the high-density data in I and the standard semiclassical Boltzmann transport theory including realistic effects of screened interface charged impurity and surface roughness scattering. This quantitative agreement remains good down to about $`5\times 10^{11}cm^2`$ below which the experimental $`\rho (T)`$ shows stronger nonlinear temperature dependence at higher ($`>3K`$) temperature, but the intermediate temperature ($`1K`$) slope $`d\sigma /dT`$ continues to agree well with the prediction of the realistic Boltzmann theory. The quantitative agreement between the Boltzmann theory and the data of I can be further improved (down to lower densities) by including scattering by bulk charged impurity centers and by adjusting the carrier density and/or carrier effective mass (as well as an improved screening function including effects of scattering and local field correction).

We conclude by pointing out that the $`\rho (T)`$ data presented in I shows smooth evolution from high to low density with only the quantitative temperature dependence becoming stronger in a continuous manner with decreasing carrier density. This indicates that the temperature dependent disorder (e.g. screened Coulomb scattering) is playing an important role since at high and intermediate densities one gets quantitative agreement between the realistic Boltzmann theory and the experimental $`\rho (T)`$ — qualitatively, lowering density is weakening screening leading to stronger temperature dependent disorder and consequently larger resistivity.


format: text

text:
# Revealing polarons with high pressure on low electron-doped manganites

## I INTRODUCTION

The study of high pressure effects on the magnetic transitions and the electrical transport properties of manganite compounds can provide useful information about the relevant mechanism in magnetic ordering and its relation to transport properties. Indeed, this was the case for the pressure dependence of the Curie temperature ($`T_c`$($`P`$)), mostly studied, up to now, in the hole-doped part of the phase diagram of these compounds.

It was experimentally shown Laukhin et al. (1997); Moritomo et al. (1997) that the pressure coeficient (dln($`T_c`$)/d$`P`$) as a function of $`T_c`$ follows a sort of universal curve, independent of the mean A ionic radius RA (1.124 Å$`<`$ RA $`<`$ 1.147 Å) for the A0.67B0.33MnO3 compounds (A = Pr, Sm, Nd, Y, La; B = Ca, Sr), or the doping level $`x`$ (0.2 $`x`$ 0.4) for La1-xSrxMnO3 and Nd1-xSrxMnO3.

It has also been shown Laukhin et al. (1997) that a qualitative understanding of this curve can be obtained in terms of the double exchange (DE) model Zener (1951); Anderson and Hasagawa (1955). Within this model, pressure increases the transfer integral of the eg electron, which is hopping from Mn3+ to Mn4+. This results in a broadening of the bandwidth (W0) which yields an increase of $`T_c`$. This can be understood by considering the pressure dependence of two geometric or steric factors, both of which control W0: The Mn-O distance (dMn-O) and the bending angle ($`\mathrm{\Theta }`$) of the Mn-O-Mn bond. However, calculations based on a simple model related to pressure variations of the steric factors do not give a good quantitative agreement with the experimental data. A better quantitative prediction can be obtained by considering the polaronic modifications of W0 due to the Jahn-Teller (JT) cooperative effect, which was clearly established for some of the manganite members Zhao et al. (1996). It has been shown that a positive contribution to the pressure dependence of $`T_c`$ arises from the negative pressure derivative of the isotope exponent, $`\alpha =1/2\gamma E_{JT}/h\omega `$, where $`\gamma `$ is a positive constant ($``$ 1), EJT is the Jahn-Teller energy and $`\omega `$ an appropriate optical frequency. Even though, the experimental results could not be suitably fitted by this model, as for example, in the case of the studies of the oxygen-isotope effects under pressure in the La0.65Ca0.35MnO3 samples Wang et al. (1999), the dlnTc/d$`P`$ value observed for the 18O sample (23%) is higher than the one predicted by this model(6%). It was shown Lorenz et al. (2001) that this controversy can be solved considering the polaron theory in the intermediate electron-phonon coupling region ($`\lambda `$ 1) and in the adiabatic approximation Alexandrov et al. (1994), which also gives an adequate description of the electrical conductivity for temperatures above $`T_c`$.

The study of electron-doped manganites can shed light on the issue whether polarons are relevant or not to the $`T_c`$($`P`$) dependence in manganites. As the low eg level occupancy of the Mn4+ ions no longer favors the JT distortion, JT polaronic effects are not expected to contribute to transport properties or to magnetic ordering in the same way they are suspected to participate in the hole-doped part of the phase diagram. Then, a small $`T_c`$($`P`$) dependence for the electron-doped manganites would reveal the important role played by JT polarons in the hole-doped samples and may also indicate, on the other hand, an active role of other type of charge carriers, like lattice or magnetic polarons, associated with the electron-doped manganites in other studies Batista et al. (1998); Chen and Allen (2001).

For low doping levels (0.05 $`x`$ 0.15), Ca1-xYxMnO3 (CYMO) and Ca1-xSmxMnO3 (CSMO) compounds are typical electron-doped manganites. For this doping regime a competition between antiferromagnetic (AF) superexchange and ferromagnetic (FM) double exchange (DE) interactions has been revealed by dc-magnetization and electrical transport measurements as a function of temperature and magnetic field for the CYMO system Garbarino et al. (2001); Aliaga et al. (2001, 2002). Neutron diffraction studies Martin et al. (1999a, 2000) for CSMO ($`x`$=0.10 and $`x`$=0.15) showed that, besides the FM phase, a Pnma structure associated with a G-type AF is also present in the paramagnetic state. At 10 K, the sample with $`x`$=0.10 exhibit a G-type AF phase, but for $`x`$=0.15 the C type magnetic order (characterized by a P 21/m space group) dominates, coexisting with smaller G AF clusters.

In this paper we study the high pressure effects on the resistivity and ac susceptibility as a function of temperature of the low electron-doped CYMO and CSMO manganese perovskites. Although an expected weak pressure dependence of $`T_c`$($`P`$) was observed, a quantitative agreement could not be obtained by solely considering the pressure sensitivity of the steric factors within a DE scenario. Our results indicate that the characteristic hopping conduction and the value of the pressure coefficient d(ln$`T_c`$)/d$`P`$ can be ascribed to the polaronic nature of carriers, which can be associated with small Fröhlich polarons in a weak to intermediate coupling regime.

## II EXPERIMENTAL

Single phase and well oxygenated Ca1-xYxMnO3 (CYMO; $`x=`$0.06; 0.07; 0.08; 0.10) and Ca1-xSmxMnO3 (CSMO; $`x=`$0.10; 0.15) ceramic samples were prepared by solid state reaction. This samples had been previously studied and the details of their synthesis and characterization have been published elsewhere Garbarino et al. (2001); Aliaga et al. (2001); Maignan et al. (1998); Martin et al. (1999a); Mahendiran et al. (2000); Respaud et al. (2001); Martin et al. (1999b, 2000); Vega et al. (2002).

Resistivity ($`\rho `$($`T`$)) and ac susceptibility ($`\chi _{ac}`$($`T`$)) were measured as a function of temperature ( 4 K $`T`$ 300 K) applying high hydrostatic pressures up to 1 GPa. A self-clamping cell was used with a 50-50 mixture of kerosene and transformer oil as the pressure transmitting medium. $`\rho `$($`T`$) was measured following a standard 4 terminal DC technique, while $`\chi _{ac}`$($`T`$) was evaluated using an applied ac magnetic field of 1 Oe at an excitation frequency of $``$ 1 kHz.

## III RESULTS AND DISCUSSION

Resistivity as a function of temperature for different pressures is displayed in Fig. 1 for CYMO ($`x`$=0.08; 0.10) and CSMO ($`x`$=0.10). A metallic conductivity at room temperature followed by a semiconducting-like behavior when decreasing temperature can be observed for both compounds. The low temperature resistivity shows in all cases a divergence characteristic of an insulator behavior.

The ferromagnetic transition temperature ($`T_c`$) can be associated with an inflexion or a small drop in the resistivity, depending on doping level, as well as with a peak in the ac susceptibility, as can be noticed in the inset of Fig. 1. The low temperature part of $`\chi _{ac}`$($`T`$), which can be associated with the behavior of a cluster glass Maignan et al. (1998), will not be considered here, as we want to focus our study on the appearance of the ferromagnetic ordering and its sensitivity to the applied pressure.

Pressure decreases resistivity and enhances the resistivity drop at $`T_c`$. $`T_c`$ was determined using the peak of the logarithmic temperature derivative of $`\rho `$($`T`$) and $`\chi _{ac}`$($`T`$) curves, as illustrated in the inset of Fig.2. The obtained pressure sensitivity of $`T_c`$, depicted in Fig. 2 and in Fig. 3, depends on the technique used to characterize the sample. A non-monotonic dependence is observed for the resistivity data: $`T_c`$ remains nearly constant for low pressures and then softly increases for higher pressures with a slope $``$ 5-7 K GPa-1. On the other hand, susceptibility data shows essentially a linear increase of $`T_c`$ with pressure with a similar slope.

These differences can be associated with the fact that $`\rho `$(T) is related not only to the magnetic ordering of the sample but also to a percolation problem. The percolation scenario can be generated by the coexistence of different conducting phases and by their pressure-dependent relative distribution. This was already demonstrated by neutron diffraction studies of these samples, as the coexistence of FM clusters in an AFM matrix. Although the pressure dependence of $`T_c`$ obtained from $`\chi _{ac}`$($`T`$) should then be more reliable, both methods produce similar results if only the high pressure portion of the curves is analyzed. Considering this criterion, the pressure coefficient, dln$`T_c`$/d$`P`$, can be determined from the $`T_c`$($`P`$) curves obtained using both techniques. As shown in Fig.4, this coefficient is for all the samples (0.05 $`\pm `$ 0.01) GPa-1, which is far from the usual values Moritomo et al. (1997) (universal curve) obtained for hole-doped manganites in this $`T_c`$ range ($``$ 0.20 GPa-1).

In a DE scenario the FM transition temperature ($`T_c`$) is proportional to the bare bandwidth ($`W_0`$), which depends on the steric parameters, dMn-O and $`\mathrm{\Theta }`$($`P`$Medarde et al. (1995). It can be shown that the pressure coefficient can then be estimated as Laukhin et al. (1997):

$$\frac{dlnT_c}{dP}=3.5\kappa (d_{MnO})2tan(\mathrm{\Theta })\mathrm{\Theta }\kappa (\mathrm{\Theta })$$
(1)
where $`\kappa (d_{MnO})`$ and $`\kappa (\mathrm{\Theta })`$ are the bond length and bond angle compressibilities, respectively. Using equation 1 and the values already measured for other manganese perovskites Radaelli et al. (1997), we can evaluate a contribution, at the most, of 0.01 GPa-1 to the pressure coefficient produced by the variation of these steric factors. This small value, compared with our experimental results in Fig.4, reveals the existence of other contributions, which may be based on the polaronic nature of electrical carriers. Though JT collective polarons are not expected in this samples, local lattice distortions or magnetic interactions may play an important role, which may be noticed in the transport properties.

In this sense, the temperature dependence of the resistivity of these samples seems to be a characteristic behavior of the electron-doped manganites, as previously reported in the literature Maignan et al. (1998); Neumeier and Cohn (2000); Garbarino et al. (2001); Aliaga et al. (2001). The observed dependencies (see Fig.1) are typical of degenerate semiconductors, where a metallic-like conductivity can be observed at room temperature, while a thermally activated hoppping of carriers dominates the conduction mechanism at low temperatures. The possibility of having a Variable Range Hopping conduction Mott and Davis (1971) was also considered, but an examination of the data rules out this possibility.

We analyze our data considering that the electrical transport at intermediate temperatures is effectively related to the hopping of polarons. As it was shown for the hole-doped manganites Lorenz et al. (2001), an adequate description should be performed in the adiabatic limit, which must be maintained for the electron-doped samples, as similar characteristic times are expected to be involved. In this case, the resistivity can be expressed as:

$$\rho (T)=\rho _0\mathrm{exp}(\frac{E_a}{k_BT})$$
(2)
with an activation energy, $`E_a`$, that can be approximated by

$$E_a=\frac{1}{2}E_pJ\frac{1}{2}E_p\frac{W_0}{2z}$$
(3)
where $`E_p`$ is the binding energy of the polaron, $`J`$ the electron transfer integral between nearest neighbor sites, $`W_0`$ is the bare bandwidth, $`z`$ the lattice coordination number and $`k_B`$ the Boltzmann constant.

In the case that carriers could be associated with a polaronic nature, the bare bandwidth $`W_0`$ should be narrowed as a consequence of the electron-phonon interaction, producing an effective smaller $`W`$. $`W`$ and $`W_0`$ can be related by the following expression:

$$W=W_0F^\lambda (E_p,\lambda ,\omega _0)$$
(4)
where $`E_p`$ is the polaron binding energy, $`\lambda =\frac{E_p}{W_0/2}`$ the dimensionless electron-phonon coupling constant and $`\omega _0`$ a characteristic optical phonon frequency. The $`F^\lambda `$ are different functions, defined according to the value of $`\lambda `$, which determines the coupling regime Alexandrov et al. (1994); Alexandrov and Mott (1995). For example, in the strong coupling regime ($`\lambda 1`$), $`F^\lambda \mathrm{exp}(E_p/\mathrm{}\omega _0)`$, but this function is not valid anymore for the intermediate coupling regime as pointed out by Alexandrov and Mott Alexandrov and Mott (1995).

From Eqs. (3) and (4) we obtain that

$$E_a=\frac{1}{2}E_p\frac{CT_c}{F^\lambda }$$
(5)
Following the description of Lorenz et al Lorenz et al. (2001), we have plotted in Fig. 5 the activation energy $`E_a`$ as a function of $`T_c`$ obtained for different pressures for the CYMO and CSMO with different doping levels. Except for the lowest $`T_c`$ values, where, as already mentioned, their determination is influenced by the percolation problem, an almost linear dependence is obtained. This indicates that $`E_p`$ and $`F^\lambda `$ are essentially constants within our experimental pressure range, independently of the pressure variation of all relevant parameters.

A linear fit of the curves of Fig. 5 yields an estimation of $`E_p`$ 120-180 meV. The increasing value with increasing doping may point to a tendency to localization, produced by the size-mismatch of replacing $`Ca`$ with small ions like $`Y`$ or $`Sm`$. In this case, considering that $`W_0`$ 1 eV, the coupling constant can be estimated as $`\lambda 0.3`$, which determines a weak to intermediate coupling regime. This value is over 4 times smaller than the one typically obtained for the hole-doped manganites Dagotto et al. (2001).

Based on the fact that the electron-doped manganites are very poor conductors, particularly for the low doping range here analyzed, we can assume the framework introduced for doped manganites and other oxides by Alexandrov et al Alexandrov and Bratkovsky (1999); Alexandrov and Kornilovitch (1999), where small or lattice polarons with a long range Fröhlich interaction are the quasiparticles involved in the electrical transport properties of these materials. Thus, for these small Fröhlich polarons we can consider the Fröhlich electron-phonon coupling constant

$$\alpha =\frac{e^2}{K_p}(\frac{m}{2\mathrm{}^3\omega _0})^{1/2}(\frac{1}{4\pi ϵ_0})$$
(6)
where $`K_p^1=ϵ_{\mathrm{}}^1+ϵ_s^1`$ is the effective dielectric constant and $`ϵ_s`$ and $`ϵ_{\mathrm{}}`$ the static and the high frequency dielectric constants, respectively, $`e`$ the electron charge, $`m`$ the electron band mass and $`ϵ_0`$ the vacuum permitivity.

Following Feynman’s approach for the weak to intermediate coupling regimeFeynman (1955); Mitra et al. (1987), we can express $`E_p`$ and the polaronic effective mass, $`m^{}`$ as a function of $`\alpha `$:

$$E_p=(\alpha +0.0123\alpha ^2)\mathrm{}\omega _0$$
(7)
$$\frac{m^{}}{m}=1+\frac{\alpha }{6}+0.025\alpha ^2$$
(8)
Considering a typical phonon frequency ($`\omega _0`$) of the order of 400 to 500 cm-1, estimated from Raman spectroscopy studies of CaMnO3 Abrashev et al. (2002), and taking a mean value of $`E_p=`$ 0.15 eV from our data, we obtain from Eqs. 6, 7 and 8 that $`K_p=`$ 11 $`\pm `$ 1 , $`\alpha =`$ 1.5 $`\pm `$ 0.1 and $`m^{}/m`$ 1.3, respectively. The last two values are in accordance with the assumption of being in the intermediate coupling regime Mitra et al. (1987) (1 $`\alpha 6`$; m/m $`<`$ 3). Also, assuming that the electron band mass is similar to the electron mass ($`mm_e`$), the polaron localization radius $`R_p=`$ ($`\mathrm{}`$/m$`\omega _0`$)1/2 $``$ 8 $`\AA `$ 2 dMn-O-Mn determines a polaron size in the limit of being considered small. While the radius of the lattice distortion produced by the long range Fröhlich interaction can be estimated Mott and Davis (1971) as $`R_d=5\mathrm{}K_p/(me^2)`$ 30 $`\AA `$. These values are in good accordance with the assumption of an electrical transport mediated by small Fröhlich polarons, i.e. a small size for the electron localization and a large size for the lattice distortion. To our knowledge, there is at the time no experimental determination of $`ϵ_{\mathrm{}}`$ for any electron-doped manganite, though a very high low-frequency dielectric constant $`ϵ`$ up to 108 was measured at room temperature in Ca1-xLaxMnO3 ($`x`$ 0.1) Zawilski et al. (2002). As $`K_p`$ 3 was estimated for LaMnO3 Alexandrov and Bratkovsky (1999), a higher value for the electron-doped CaMnO3 can be expected, as a consequence derived from Eq. 6 and its lower coupling constant, when compared with the former manganite.

At this point, we can try to make a rough approximation to determine what is the contribution to the pressure coefficient produced by the dependence of the bandwidth on the mass enhancement, related to the polaronic nature of carriers. To do so, we can consider the approximation of slow polarons ($`kq_p`$, where $`k`$ and $`q_p`$ are the electron and the phonon momenta, respectively)Alexandrov and Mott (1995) where the bandwidth $`W`$ is essentially modified as

$$W=\frac{W_0}{m^{}/m}$$
(9)
Using Eqs. 8 and 9, the pressure coefficient can then be expressed as

$$\frac{dlnTc}{dP}\frac{dlnW_0}{dP}\frac{m}{m^{}}\frac{d(m^{}/m)}{d\alpha }\frac{d\alpha }{d\omega }\frac{d\omega }{dP}$$
(10)
From Eqs.6 and 7 and assuming that $`d\omega /dP`$ = $`\delta \omega _0`$ (with $`\delta `$ $``$ 0.01 GPa-1) from reference Postorino et al., 2002, we can evaluate the second term of Eq. 10, which gives a positive contribution of $``$ 0.02 GPa-1, which is in good agreement with the results presented in Fig.4. Therefore, these results reveal the necessity of considering an additional contribution to the pressure coefficient besides the one which arises out of the steric factors. The additional contribution can be consistently related to the pressure sensitivity of intermediate-coupling Fröhlich polaronic effects on the bandwidth.

## IV CONCLUSIONS

We have studied the pressure sensitivity of the resistivity and of the ferromagnetic ordering of electron-doped manganese perovskites (CYMO and CSMO) in the low doping regime. The obtained values of dln$`T_c`$/d$`P`$ are considerably smaller than the ones measured for hole-doped manganites in the same $`T_c`$ range. Within a DE scenario this values can not be explained as a consequence of the pressure variation of steric factors. A better quantitative agreement can be obtained if a scenario of small intermediate-coupled Fröhlich polarons is considered.

## V ACKNOWLEDGEMENTS

We would like to acknowledge financial support by ANPCyT Grant PICT97-03-1700, Fundación Antorchas A-13622/1-31, and UBACyT JW11. We are indebted to P. Levy, F. Parisi and R. Zysler for fruitful discussions and to V. Bekeris and P. Tamborenea for a critical reading of the manuscript. We also acknowledge technical assistance from C. Chiliotte, D. Giménez, E. Pérez Wodtke and D. Rodríguez Melgarejo.


text:
# In-silico folding of a three helix protein and characterization of its free-energy landscape in an all-atom forcefield

(June 18, 2003)

## Abstract

We report the reproducible first-principles folding of the 40 amino acid, three-helix headpiece of the HIV accessory protein in a recently developed all-atom free-energy forcefield. Six of twenty simulations using an adapted basin-hopping method converged to better than 3 Å backbone RMS deviation to the experimental structure. Using over 60,000 low-energy conformations of this protein, we constructed a decoy tree that completely characterizes its folding funnel.

Available genomic and sequence information for proteins contains a wealth of biomedical information that becomes accessible when translated into three-dimensional structure baker01 . While theoretical models for protein structure prediction schonbrunn02 ; liwo02  that partially rely on experimental information have shown consistent progressmoult01 , the assessment of de-novo strategies that rely on sequence information alone has been much less favorable bonneau01 . The development of such techniques, in particular of transferable first-principle all-atom folding methods, would significantly benefit the understanding of protein families where little experimental information is available, the prediction of novel folds as well as the investigation of protein association and dynamics which are presently difficult to probe experimentally. Recent progress for small peptides liwo02 ; snow02 ; simmerling02 ; schug03  documents both the feasibility of this approach as well as its limitations hansmann02 ; lin03 , in particualr those associated with the direct simulation of the folding process through molecular dynamics duan98 .

Overwhelming experimental evidence supports the thermodynamic hypothesis anfinsen73  that many proteins are in thermodynamic equilibrium with their environment: their native state thus corresponds to the global minimum of their free energy landscape onuchic97 . The free energy of the system is accessible either indirectly by explicit ensemble averaging of the combined internal energy of protein and solvent, or directly in a free-energy forcefield where an implicit solvation model approximates direct interactions with the solvent as well as most of the entropic contributions. We developed an all-atom protein forcefield (PFF01) herges03b ; schug03 ; herges04a  with an area-based implicit solvent model that approximates the free energy of peptide conformations in the natural solvent. Using a rational decoy approach this forcefield was optimized to correctly predict the native structure of the 36-amino acid headgroup of villin duan98 ; hansmann02 ; lin03 . Without further parameter adjustment we then simulated the structurally conserved 40 amino-acid headpiece of the autonomously folding HIV accessory protein (1F4I-40) withers00  with a modified basin hopping technique nayeem91 ; wales96 . Out of twenty simulations the five energetically lowest correctly reproduced the NMR structure of this three-helix protein with a backbone RMS deviation of less than 3 Å. The combination of decoy based model development for the free energy with efficient stochastic optimization methods suggests a viable route for protein structure prediction at the all atom level with present day computational resources.

Model: We have recently developed an all-atom protein forcefield (PFF01), which was used to reproducibly fold the 20 amino acid trp-cage protein schug03 . PFF01 herges04a  comprises an atomically resolved electrostatic model with group specific dielectric constants and a Lennard Jones parameterization that was adapted to the experimental distance distributions from crystal structures of 138 proteins abagyan94 . Interaction with the fictitious solvent are modeled in a simple solvent accessible surface approach eisenberg86 , where the solvation free-energies per unit surface were fitted to the enthalpies of solvation of the Gly-X-Gly series of peptides sharp91 . The only low-energy degrees of freedom available to the peptide during the folding process are rotations of the dihedral angles of the backbone and the sidechains, these are the only moves considered during the simulation. There are two move-classes, small random rotations about a single angle and library moves, which set a particular backbone dihedral to a permitted value in the Ramachandran plot.

If an accurate model for the free energy of the protein in its environment is available, stochastic optimization methods can be used to locate the global optimum of the free-energy landscape orders of magnitude faster than traditional simulation techniques. We adapted the basin hopping technique nayeem91  (BHT) for protein simulations by replacing a single minimization step with a simulated annealing run kirkpatrick83  with self-adapting cooling cycle and length. At the end of one annealing step the new conformation was accepted if its energy difference to the current configuration was no higher than a given threshold energy $`\epsilon _T`$, an approach recently proven optimal for certain optimization problems schneider98 . Each simulation was performed in three separate steps: First we used a high temperature bracket of 800/300 K ($`\epsilon _T`$ =15 K) for the annealing window and reduced the strength of the solvent terms in the forcefield by 20%. The second step started from the final configurations of the first run, used the same annealing window but the full solvent interactions. In the third step the resulting structures were further annealed in a low temperature bracket of 600/3 K ($`\epsilon _T`$=1K). Within each annealing run the temperature was geometrically decreased, also the number of steps per annealing run was gradually increased to ensure better convergence. In total each simulation comprised 107 energy evaluations with a computational effort roughly corresponding to a 10 ns MD simulation (in vacuum with 1 fs timestep). After this time no significant energy fluctuations occurred in the simulations, indicating that each had settled into a metastable configuration.

Results: Using PFF01 we performed 20 independent modified basin hopping simulations of the structurally conserved 40 amino-acid headpiece of the HIV accessory protein (pdb-code 1F4I, sequence: QEKEAIERLK ALGFEESLVI QAYFACEKNE NLAANFLLSQ). The best structures found in each run were ranked according to their energy and the backbone RMS deviation (RMSB) to the NMR structure was computed. Table (1) demonstrates that the five lowest structures had to good accuracy converged to the NMR structure of the protein. The first non-native decoy appears in position six, with an energy deviation of 5 kcal/mol (in our model) and a significant RMSB deviation. The table demonstrates that all low-energy structures have essentially the same secondary structure, i.e. position and length of the helices are always correctly predicted, even if the protein did not fold correctly. The degree of secondary structure content and similarity decreases for the decoys with higher energy (data not shown) in good correlation with their energy. From the standpoint of the optimization approach (but not necessarily for the physical folding scenario) this suggests that successful folding requires the formation of the correct secondary structure content as a prerequisite.

The good agreement between the folded and the experimental structure is also evident from Figure (1), which shows the secondary structure alignment of the native and the folded conformations. The good physical alignment of the helices illustrates the importance of hydrophobic contacts to correctly fold this protein. An independent measure to assess the quality of these contacts is to compare the Cβ-Cβ distances (which correspond to the NOE constraints of the NMR experiments that determine tertiary structure) in the folded structure to those of the native structure. The color coded Cβ-Cβ distance in Figure (2) demonstrates a 66 % (80 %) coincidence of the Cβ-Cβ distance distances to within one (1.5) standard deviations of the experimental resolution. The dark diagonal block indicate intra-helical contacts, which are, perhaps not too surprisingly, resolved to very good accuracy. The off-diagonal dark blocks, however, indicate that also a large fraction of long range native contacts is reproduced correctly.

Starting from intermediate structures of the folding simulations we generated over 60,000 low-energy conformations (decoys). Decoys with a root mean square deviation of the backbone (RMSB) deviation of less than 3 Å were grouped into families with free energy brackets of 2 kcal/mol. We then resolved the topological hierarchywales96 ; brooks01  of the associated potential energy surface through the construction of a decoy tree (Figure (3)) that illustrates the low-energy structure of the free energy surface. Beginning from the best conformation, we draw a vertical line for each decoy family in this window. Moving upward in energy the number of decoys in each family grows almost exponentially in the low energy region which we can resolve well. As a result the diversity of each family grows until different families unite. Family membership is associative, i.e. as soon as two decoys in different branches have an RMSB deviation of less than 3 Å, all members of both families belong to one superfamily. Pictorially this representation results in an inverted tree-like structure that characterizes the topology of the metastable states of the free energy surface.

Trees with very short stems and many low-energy branches are characteristic of glassy potential energy surfaces, which are associated with Levinthals paradox levinthal68 ; honig99  in the context of protein folding. Well structured trees with few terminal branches suggest the existence of a folding funnel onuchic97 , consistent with the “new” paradigm for protein folding karplus94 . From this perspective the structure in Figure (3) this tree is indicative of the existence of a very broad folding funnel onuchic97  with pronounced competing secondary metastable conformations, which are depicted as the non-native terminal minima of the tree. The discretization of the energy axis in intervals of 2 kcal/mol starting from the native conformation results in a smoothing of the free-energy surface. Each line in the figure corresponds to a family containing hundreds tothousands of structures, which are all associated with the same low-lying metastable conformation (the terminal point of the branch). Simulations trapped in such a metastable state must overcome a potential energy barrier of the order of the energy difference to the next highest branching point of the tree to visit another structure. The branching points of the tree were constructed only from structures of the decoy set and not through independent transition state search among the terminal structures. In addition, main-chain entropy is neglected in this analysis, which results in an overestimation of the barrier. Further investigations to more accurately determine the transition states are presently under way.

The lowest competing terminal branch (branch C), associated with decoy D06 in Table (1), is less than 5 kcal/mol above the best native decoy. Decoy D06 has comparable energy to competing decoys but much larger RMSB deviation and has few long-range native contacts. This structure (see Figure (4)) has also three helices of comparable length and sequence location, but differs from the native structure in the relative alignment of the helices with respect to each. The RMSB deviation of the low-lying terminal structures to the NMR structures is large (i.e. comparable with the RMSB deviation to unfolded structures), indicating that conserved secondary structure elements, rather than distance constraints characterize the folding funnel. The low lying local minima are thus characterized by varying spatial arrangements of similar secondary structure elements, a property that is ill represented by either RMSB deviation or the number of native contacts.

Discussion: With this work we have demonstrated that accurate free-energy forcefields can predict the native structure of proteins with nontrivial tertiary structure with present day computational resources. This result represents an in-silico realization of the thirty year-old thermodynamic hypothesis that proteins are in thermodynamic equilibrium with the environment under physiological conditions anfinsen73 . Under this hypothesis, the native structure of the protein can be predicted using stochastic optimization methods orders of magnitude faster than by direct simulation. Our results demonstrate that the important influence of the solvent can be modeled with a relatively simple solvent accessible surface approach.

The analysis of the free energy surface supports the funnel paradigm of protein folding for a nontrivial protein with a significantly larger hydrophobic core than was previously possible. The relatively small number of terminal branches of the decoy tree offers the first glimpse on the experimentally inaccesible structure of the folding funnel. It suggest the exsitence of a broad folding funnel with well defined secondary metastable states which may constitute important folding intermediates. The free energy optimization approach used here permitted the characterization of these low-lying states, which surprisingly share very similar secondary structure with the native configurations. Investigations of other proteins must show, whether this pattern persist also for other proteins.

It should be noted that the success of the optimization approach depends strongly on the ability of the optimization technique to differentiate between the low-lying minima of the FES in a realistic forcefield. The performance of optimization methods is often strongly problem dependent, but with 1F4I, 1VII and 1L2Y three nontrivial model systems exist on which different optimization methods can be evaluated. The decoy sets generated and insights regarding low-lying metastable states can also serve as a test-bed for the development of coarse-grained protein models. PFF01 is presently be validated for other peptides and proteins and rationally evolved to correctly predict the structure of larger fold families.

Acknowledgment: We are grateful to S. Gregurick, J. Moult and J. Pedersen for discussions and portions of the code used in these simulations. This work was supported by the Deutsche Forschungsgemeinschaft, the Bode Foundation and the BMWF.


text:
# 1 Introduction

## 1 Introduction

Extremely asymmetrical scattering (EAS) is a radically new type of Bragg scattering in slanted, periodic, strip-like wide gratings, when the first diffracted order satisfying the Bragg condition (scattered wave) propagates parallel to the front grating boundary . The main characteristic feature of EAS is the strong resonant increase of the scattered wave amplitude, compared to the amplitude of the incident wave at the front boundary . Other unique features of EAS include additional resonances in non-uniform gratings , detuned gratings , in wide gratings when the scattered wave propagates at a grazing angle with respect to the boundaries , and the unusually strong sensitivity of EAS to small variations of mean structural parameters at the grating boundaries . The additional resonances may result in amplitudes of the scattered and incident waves in the grating that can be dozens or even hundreds of times larger than that of the incident wave at the front boundary .

One of the main physical reasons for all these unusual features of EAS is the diffractional divergence of the scattered wave inside and outside the grating . Indeed, the scattered wave results from scattering of the incident wave inside the grating, and propagates parallel to the grating. Thus, it is equivalent to a beam located within the grating, of an aperture equal to the grating width. Such a beam will diverge outside the grating due to diffraction. Therefore, steady-state EAS is formed by the two physical effects – scattering and diffractional divergence. Based on the equilibrium between these processes, an approximate analytical method of analysis of EAS, applicable to all types of waves, has been developed .

A reasonable question is which of the EAS resonances are practically achievable? Any resonance is characterized by a time of relaxation, and if this time is too large, the corresponding resonance cannot be achieved in practice. In the case of EAS, large relaxation times may result in excessively large apertures of the incident beam being required for the steady-state regime to be realized . It is obvious that knowledge of relaxation times and relaxation (non-steady-state) processes during EAS is essential for the successful development of practical applications of this highly unusual type of scattering.

Until recently, only estimates of relaxation times for EAS in uniform gratings have been made . The analysis was based on physical assumptions and speculations , rather than direct treatment of non-steady-state regime of EAS. Simple analytical equations for the relaxation times were derived . However, the accuracy of these equations is questionable, since their derivation did not take into account re-scattering processes in the grating. Moreover, these equations are not applicable in the presence of the additional resonances .

Therefore, the first aim of this paper is to present an efficient numerical algorithm for the rigorous analysis of non-steady-state EAS of bulk electromagnetic waves in wide uniform and non-uniform holographic gratings. The second aim is to investigate non-steady-state EAS and accurately determine relaxation times in narrow and wide uniform periodic gratings. In particular, amplitudes of the incident wave (0th diffracted order) and scattered wave (+1 diffracted order) inside and outside the grating will be determined as functions of time and coordinates after the incident wave is “switched on.”

## 2 Structure and methods of analysis

Consider an isotropic medium with a slab-like, uniform, holographic grating with a sinusoidal modulation of the dielectric permittivity (see fig 1):

$$\begin{array}{cccc}\hfill ϵ_\mathrm{s}& =& ϵ+ϵ_\mathrm{g}\mathrm{exp}(\mathrm{i}q_xx+\mathrm{i}q_yy)+ϵ_\mathrm{g}^{}\mathrm{exp}(\mathrm{i}q_xx\mathrm{i}q_yy),\hfill & \text{if}\mathrm{\hspace{0.33em}0}<x<L,\hfill \\ \hfill ϵ_\mathrm{s}& =& ϵ,\hfill & \text{if}x<0\text{or}x>L,\hfill \end{array}$$
(1)
where $`L`$ is the grating width, $`ϵ_\mathrm{g}`$ is the grating amplitude, the mean dielectric permittivity $`ϵ`$ is the same inside and outside the grating, $`q_x`$ and $`q_y`$ are the $`x`$ and $`y`$ components of the reciprocal lattice vector $`𝐪`$, $`q=2\pi /\mathrm{\Lambda }`$, $`\mathrm{\Lambda }`$ is the period of the grating, and the coordinate system is shown in figure 1. There is no dissipation in the medium ($`ϵ`$ is real and positive), and the structure is infinite along the $`y`$ and $`z`$ axes.

Non-steady-state EAS in this structure occurs when the incident wave is switched on at some moment of time. Then, both the incident and scattered wave amplitudes inside and outside the grating evolve in time and gradually relax to their steady-state values at $`t=+\mathrm{}`$. This occurs when the incident pulse, having an infinite aperture along the $`y`$ and $`z`$ axes, is the product of a sinusoid and a step function of time.

However, the numerical analysis of this infinitely long (in time) pulse is inconvenient, since its Fourier transform contains a $`\delta `$-function. Therefore, in order to calculate non-steady-state amplitudes of the incident and scattered waves in the structure at an arbitrary moment of time $`t=t_0`$, we consider a rectangular (in time) sinusoidal incident pulse of time length $`2t_0`$, amplitude $`E_{00}`$, and angle of incidence $`\theta _0`$ (fig. 1), which has a simple analytical Fourier transform. We also assume that the incident pulse is switched on instantaneously (with the amplitude $`E_{00}`$) everywhere in the grating. That is, we ignore the process of propagation of the pulse through the grating. This assumption allows a substantial reduction in computational time. It is correct only if the time $`t_0`$, at which non-steady-state amplitudes are calculated, is significantly larger than $`t_\mathrm{g}=(L/\mathrm{cos}\theta _0)/(cϵ^{1/2})`$ – the time of propagation of the incident pulse across the grating. We will see below that this condition is satisfied (for most times) since typical relaxation times $`\tau t_\mathrm{g}`$.

Thus, at $`t=0`$, the front of the pulse is at the rear grating boundary, i.e. the incident field is zero behind the grating (at $`x>L`$). This means that if the angle of incidence $`\theta _00`$ (fig. 1), then the amplitude of the incident pulse is not constant along a wave front. Therefore, the incident pulse experiences diffractional divergence, and the spatial Fourier transform should be used. However, this divergence can be ignored in our further considerations, since it may be significant only within several wavelengths near the front of the pulse, i.e. for times $`10^{14}`$ s or smaller (this is the time for the pulse to travel a distance of a few wavelengths). This time is several orders of magnitude less than typical relaxation times (see below). Therefore, only for very short time intervals after switching an incident pulse on can noticeable errors result from the above approximations.

The frequency spectrum of this input is determined from its analytical Fourier transform. As a result, the incident pulse is represented by a superposition of an infinite number of plane waves having different frequencies and amplitudes, all incident at $`\theta _0`$. Note that the spectrum is the incident pulse depends on the pulse width $`2t_0`$, and is different for every time $`t_0`$ at which we calculate the fields. The steady-state response of the grating to each plane wave is determined by means of the rigorous theory of steady-state EAS , based on the enhanced T-matrix algorithm , or the approximate theory  (if applicable). Thus the frequency spectrum of the incident and scattered waves inside and outside the grating is obtained, and their amplitudes as functions of the $`x`$-coordinate at any moment time can be found using the inverse Fourier transform. Due to the geometry of the problem, the non-steady-state incident and scattered wave amplitudes do not depend on the $`y`$-coordinate.

Note that the inverse Fourier transform is taken at $`t=t_0`$, i.e. at the middle of the incident pulse, in order to minimize numerical errors. The inverse Fourier transform is found by direct integration, to allow a non-uniform set of frequency points to be used. The rapid variation of the frequency spectra at certain points, and the wide frequency spectrum of the input for small $`t_0`$ make it infeasible to use the fast Fourier transform. The calculations are carried out separately for each moment of time. Therefore, there is no accumulation of numerical errors, and the errors that are noticeable at small times $`10^{14}`$ s (see above) disappear at larger times.

This approach is directly applicable for all shapes of the incident pulse, as well as for an incident beam of finite aperture. However, for beams of finite aperture, we should also use the spatial Fourier integral.

Since the approximate theory is directly applicable for all types of waves (including surface and guided waves in periodic groove arrays ), the developed approach to non-steady-state scattering is also applicable for all types of waves (if used together with the approximate theory).

## 3 Numerical Results

Using the described numerical algorithm, non-steady-state EAS of bulk TE electromagnetic waves in uniform holographic gratings given by eqn (1) has been analyzed. The grating parameters are as follows: $`ϵ=5`$, $`ϵ_\mathrm{g}=5\times 10^3`$, $`\theta _0=45^{}`$, and the wavelength in vacuum (corresponding to the central frequency $`\omega _0`$ of the spectrum of the incident pulse) $`\lambda _0=1\mu `$m. The Bragg condition is assumed to be satisfied precisely for the first diffracted order at $`\omega =\omega _0`$:

$$𝐤_1(\omega _0)=𝐤_0(\omega _0)𝐪,$$
(2)
where $`𝐤_0(\omega )`$ are the frequency dependent wave vectors of the plane waves in the Fourier integral of the incident pulse, $`𝐤_1(\omega )`$ are the wave vectors of the corresponding first diffracted orders (scattered waves), $`k_1(\omega _0)=k_0(\omega _0)=\omega _0ϵ^{1/2}/c`$, $`𝐤_1(\omega _0)`$ is parallel to the grating boundaries (fig. 1), and $`c`$ is the speed of light. Note that if $`\omega \omega _0`$, $`𝐤_1(\omega )`$ is not parallel to the grating boundaries .

Typical relaxation of amplitudes of the scattered wave (+1 diffracted order) and incident wave (0th diffracted order) over time inside and outside the gratings of widths $`L=10\mu `$m, $`28\mu `$m, and $`80\mu `$m is shown in the animations in fig. 2. The time dependencies of non-steady-state amplitudes of the scattered wave (+1 diffracted order) at $`x=0`$, $`L/2`$, and $`L`$, and the transmitted wave (0th diffracted order at $`x=L`$) are shown in fig. 3 for the same grating widths. Note that $`L=28\mu `$m is approximately equal to the critical width $`L_\mathrm{c}`$ . Physically, $`L_\mathrm{c}/2`$ is the distance within which the scattered wave can be spread across the grating by means of the diffractional divergence, before being re-scattered by the grating . All the curves in figures 2 and 3 can equally be regarded as approximate or rigorous, since both the theories in the considered structure give practically indistinguishable results .

If $`L<L_\mathrm{c}`$, then the relaxation at the front and rear grating boundaries occurs practically simultaneously (see figures 2(a) and 3(a)). In the middle of the grating, the scattered wave amplitude grows slightly faster at small $`t`$. This is due to energy losses from the scattered wave, caused by diffractional divergence of the wave near the boundaries (the edges of the scattered beam). This effect becomes more obvious with increasing grating width, and is especially strong if $`LL_\mathrm{c}`$ (fig. 3(b)). This is because the effect of the diffractional divergence in the middle of the grating (in the middle of the beam) becomes weaker with increasing grating width (i.e. beam aperture). At the same time, at the edges of the beam (at the grating boundaries) the divergence is strong, resulting in a significant reduction of the rate of change of the non-steady-state scattered wave amplitudes (compare curves (1) and (2) with curve (3) in fig. 3(b)). However, in wide gratings (with $`L>L_\mathrm{c}`$), the relaxation occurs first near the front boundary, and then the steady-state scattered wave amplitudes start spreading towards the rear boundary (fig. 2(c)). Therefore, the relaxation process in the middle of the grating and especially at the rear boundary tends to slow down compared to that near the front boundary – compare curves (1) and (2) in figs. 3(b,c).

The relaxation near the front boundary takes place more or less smoothly, except for some not very significant oscillations in wide gratings near the end of the relaxation process (curves (1) in figs. 3(a–c)). The same happens in the middle of the grating and at the rear boundary in narrow gratings ($`L<L_\mathrm{c}`$). However, if $`LL_\mathrm{c}`$, the relaxation curves in the middle of the grating and at the rear boundary display a complex and unusual behavior at small and large time intervals (see curves (2) and (3) in figs. 3(b,c)).

The unusual and complex character of relaxation processes in wide gratings is especially obvious from the time dependencies of non-steady-state incident (transmitted) wave amplitude at the rear grating boundary (figs. 2(c), 3(d)). In wide gratings, these dependencies are characterized by significant oscillations with minima that are close to zero – see fig. 2(c) and curve (iii) in fig. 3(d). The typical number of these oscillations increases with increasing grating width (compare curves (ii) and (iii) in fig. 3(d)). The minima in these oscillations tend to zero with increasing $`L`$. When the amplitude of the transmitted wave at the rear grating boundary is close to zero, almost all energy of the incident wave is transferred to the scattered wave. In wide gratings this may happen several times during the relaxation process (figs. 2(c), 3(d)).

The relaxation times at the grating boundaries and at $`x=L/2`$, determined as the times at which the amplitudes reach $`(11/\mathrm{e})`$ of their steady-state values, are:

| $`L`$ | $`\tau |_{x=0}`$ | $`\tau |_{x=L/2}`$ | $`\tau |_{x=L}`$ |
| --- | --- | --- | --- |
| $`10\mu `$m | $`4\times 10^{11}`$ s | $`4\times 10^{11}`$ s | $`4\times 10^{11}`$ s |
| $`28\mu `$m | $`6\times 10^{12}`$ s | $`1.7\times 10^{12}`$ s | $`2.4\times 10^{11}`$ s |
| $`80\mu `$m | $`7.4\times 10^{12}`$ s | $`1.2\times 10^{11}`$ s | $`3.4\times 10^{11}`$ s |

The relaxation times for narrow gratings, determined by means of the developed algorithm, are about three times smaller than those previously estimated . The significant overestimation of relaxation times in paper  is due to not taking into account the effects of re-scattering of the scattered wave in the grating. Re-scattering reduces the transmitted wave amplitude during the relaxation process (fig. 3(d)). Thus, the energy flow into the scattered wave is increased, and the relaxation times are reduced (for more detailed discussion of re-scattering see ).

During the process of relaxation, the scattered wave propagates a particular distance along the $`y`$-axis. Therefore, the relaxation times determine critical apertures $`l_\mathrm{c}`$ (along the $`y`$-axis) of the incident beam, that are required for steady-state EAS to be achieved (see also ). For example, for fig. 3(a) (with the largest values of $`\tau `$), the critical aperture of the incident beam is $`l_\mathrm{c}=c\tau ϵ^{1/2}\mathrm{cos}\theta _00.4`$ cm, which does not present any problem in practice.

## 4 Conclusions

This paper has developed an efficient numerical algorithm for the approximate and rigorous numerical analysis of non-steady-state EAS of waves in uniform slanted gratings. An unusual type of relaxation with strong oscillations of the incident and scattered wave amplitudes has been predicted for gratings wider than the critical width.

If used in conjunction with the approximate theory of steady-state EAS , the developed algorithm is immediately applicable for the analysis of non-steady-state EAS of all types of waves, including surface and guided modes in periodic groove arrays.

Typical relaxation times have been calculated for narrow and wide gratings. It has been shown that these times are significantly smaller than previous estimates . The corresponding critical apertures of the incident beam that are required for achieving steady-state EAS have also been determined. The obtained results demonstrate that steady-state EAS can readily be achieved in practice for reasonable beam apertures and not very long gratings.


text:
# Bounce Models in Brane Cosmology and a Gravitational Stability Condition

## I INTRODUCTION

In the brane-world scenarios, our conventional universe is a 3-brane embedded in a higher dimensional space. While gravity can freely propagate in all dimensions, the standard matter particles and forces are confined to the 3-brane only . In recent years, five-dimensional (5D) brane cosmological models have received extensive studies \[2-9\]. It is noticed that one of the many interesting features of brane models concerns the big bang singularity: Due to the existence of extra dimensions, the conventional big bang singularity could be removed. So the big bang is perhaps not the beginning of time but a transition from a pre-existing phase of the universe to the present expanding phase, and our universe may have existed for an infinite time prior to the putative big bang. In the ekpyrotic model , it was suggested that the universe was produced from a collision between our brane and a bulk brane. In the cyclic model , the universe undergoes an endless sequence of cosmic epochs that begin with a big bang and end in a big crunch. In the bounce models \[7-9,12-14\], it was shown that the scale factor could evolve across a finite (but non-zero) minimum which represents a bounce (as opposed to a bang). The purpose of this paper is to study the bounce property of brane models more generally.

In a previous paper , a five-dimensional big bounce cosmological solution with a non-compact fifth dimension was presented. By using this solution as valid in the bulk, a global brane model is derived in Ref. , in which the model has two 3-branes with the extra dimension compactified on an $`S_1/Z_2`$ orbifold. This brane model is of the type of Binetruy, Deffayet and Langlois  in which no cosmological constant was introduced in the bulk. In this paper, we are going to generalize it by adding a cosmological constant in the bulk, and then to study evolutions of the models.

The plan of this paper is as follows. In Section II, we look for general solutions for cosmological models with two 3-branes and for three cases with $`\mathrm{\Lambda }=0`$, $`\mathrm{\Lambda }>0`$ and $`\mathrm{\Lambda }<0`$, respectively. In Section III, we study the gravitational force acted on a test particle in the vicinity of a brane and derive a stability condition. In Section IV, we give several simple exact solutions as an illustration, and study the global evolutions and the stability of the brane models as well as the big bounce singularity. Section V is a short discussion.

## II GENERAL BRANE SOLUTIONS WITH A COSMOLOGICAL CONSTANT

Let the 5D metric being

$$dS^2=B^2dt^2A^2\left(\frac{dr^2}{1kr^2}+r^2d\mathrm{\Omega }^2\right)dy^2,$$
(1)
where $`B=B(t_,y)`$ and $`A=A(t_,y)`$ are two scale factors, $`k`$($`=\pm 1`$ or $`0`$) is the 3D curvature index, and $`d\mathrm{\Omega }^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)`$. This metric describes 5D cosmological models with spherical symmetry in 3D and with a static fifth dimension. Now we consider brane models with the extra fifth dimension $`y`$ compactified on a small circle for which the first brane is at $`y=y_1=0`$ and the second is at $`y=y_2`$ with $`y_2>0`$. The 5D Einstein equations read

$`G_{AB}`$ $`=`$ $`\kappa ^2\left(\mathrm{\Lambda }g_{AB}+T_{AB}\right),`$
$`T_{AB}`$ $`=`$ $`{\displaystyle \underset{i=1,2}{}}\left[\left(\rho _i+p_i\right)u^\alpha u^\beta p_ig^{\alpha \beta }\right]g_{\alpha A}g_{\beta B}\delta \left(yy_i\right),`$ (2)
where upper case Latin letters denote 5D indices (0,1,2,3;5), lower case Greek letters denote 4D indices (0,1,2,3), $`\kappa ^2`$ is the 5D gravitational constant, $`u^\alpha dx^\alpha /ds`$ is the 4D velocity in comoving coordinates with $`ds`$ being the 4D line-element and $`u^\alpha (u^0,0,0,0)`$, and $`\rho _i`$ and $`p_i`$ are energy density and pressure on the $`i`$-th brane, respectively. With use of the metric (1), the non-vanishing equations in (2) are found to be

$`G_{00}`$ $`=`$ $`3\left({\displaystyle \frac{\dot{A}^2}{A^2}}+k{\displaystyle \frac{B^2}{A^2}}\right)3B^2\left({\displaystyle \frac{A^{\prime \prime }}{A}}+{\displaystyle \frac{A^2}{A^2}}\right)`$ (3)
$`=`$ $`\kappa ^2B^2\left[\mathrm{\Lambda }+{\displaystyle \underset{i=1,2}{}}\rho _i\delta \left(yy_i\right)\right],`$
$$G_{05}=3\left(\frac{\dot{A}^{}}{A}\frac{\dot{A}}{A}\frac{B^{}}{B}\right)=0,$$
(4)
$$G_{55}=\frac{3}{B^2}\left[\frac{\ddot{A}}{A}+\frac{\dot{A}}{A}\left(\frac{\dot{A}}{A}\frac{\dot{B}}{B}\right)+k\frac{B^2}{A^2}\right]+\frac{3A^{}}{A}\left(\frac{A^{}}{A}+\frac{B^{}}{B}\right)=\kappa ^2\mathrm{\Lambda },$$
(5)
$`\left(1kr^2\right)G_{11}`$ $`=`$ $`r^2G_{22}=r^2\mathrm{sin}^2\theta G_{33}`$ (6)
$`=`$ $`{\displaystyle \frac{A^2}{B^2}}\left[{\displaystyle \frac{2}{A}}+{\displaystyle \frac{\dot{A}}{A}}\left({\displaystyle \frac{\dot{A}}{A}}{\displaystyle \frac{2\dot{B}}{B}}\right)+k{\displaystyle \frac{B^2}{A^2}}\right]`$
$`+A^2\left[{\displaystyle \frac{B^{\prime \prime }}{B}}+{\displaystyle \frac{2A^{\prime \prime }}{A}}+{\displaystyle \frac{A^{}}{A}}\left({\displaystyle \frac{A^{}}{A}}+{\displaystyle \frac{2B^{}}{B}}\right)\right]`$
$`=`$ $`\kappa ^2A^2\left[\mathrm{\Lambda }{\displaystyle \underset{i=1,2}{}}p_i\delta \left(yy_i\right)\right],`$
where an overdot and a prime denote partial derivatives with respect to $`t`$ and $`y`$, respectively.

Note that there are different ways to express the solutions of above equations (see, for example, ). Here we follow our previous work  for convenience. We see that equation (4) can be integrated once, giving

$$B=\frac{\dot{A}}{\mu (t)},$$
(7)
where $`\mu (t)`$ is an arbitrary function. Using this relation to eliminate $`B`$ from equation (3), we obtain

$$\left(A^2\right)^{\prime \prime }+\frac{2}{3}\kappa ^2\left[\mathrm{\Lambda }+\underset{i=1,2}{}\rho _i\delta \left(yy_i\right)\right]A^2=2\left(\mu ^2+k\right).$$
(8)
Also, with use of the relation (7), equation (5) can then be integrated, giving

$$\left(\mu ^2+kA^2\right)A^2=\frac{1}{6}\kappa ^2\mathrm{\Lambda }A^4+F,$$
(9)
where $`F`$ is an integration “constant” which actually could be a function of $`y`$.

In the bulk, the two equations (8) and (9) require that $`F=K=const`$. So the Einstein equation in the bulk is

$$\left(A^2\right)^2+\frac{2}{3}\kappa ^2\mathrm{\Lambda }A^4=4A^2\left[\mu ^2(t)+k\right]4K.$$
(10)
Meanwhile, we can verify that equation (6) is then satisfied identically.

From a known exact bulk solution of the equation (10) we can easily extend it from bulk to branes and obtain a global two-brane solution. The $`S_1/Z_2`$ symmetry requires that we firstly should write a bulk solution from $`g_{AB}(t,y)`$ to $`g_{AB}(t,\left|y\right|)`$. Then, according to Israel’s jump conditions, the two scale factors $`A(t,\left|y\right|)`$ and $`B(t,\left|y\right|)`$ are required to be continuous across the branes localized in $`y=y_i`$. Their first derivatives with respect to $`y`$ can be discontinuous across the branes. And their second derivatives with respect to $`y`$ can give a Dirac delta function. Thus the resulting terms with a delta function appearing in the LHS of equations (3) and (6) must be matched with the corresponding terms containing a delta function in the RHS of them in order to satisfy the field equations. For $`A_i`$ we have, for instance,

$$A_i^{\prime \prime }=\widehat{A_i^{\prime \prime }}+\left[A^{}\right]_i\delta (yy_i),$$
(11)
where $`\widehat{A_i^{\prime \prime }}`$ represents terms in $`A_i^{\prime \prime }`$ that are not contributed to the delta function $`\delta (yy_i)`$, and $`\left[A^{}\right]_i`$ is the jump in the first derivative across $`y=y_i`$, defined by $`\left[A^{}\right]_iA^{}(y_i^+)A^{}(y_i^{})`$. The $`Z_2`$ reflection symmetry leads to

$$A^{}(y_i^{})=A^{}(y_i^+),\left[A^{}\right]_i2A^{}(y_i^+).$$
(12)
Thus the Einstein equations (3) and (6) on the branes give

$`\kappa ^2\rho _i`$ $`=`$ $`{\displaystyle \frac{3}{A_i}}\left[A^{}\right]_i={\displaystyle \frac{6A^{}(y_i^+)}{A_i}},`$
$`\kappa ^2p_i`$ $`=`$ $`{\displaystyle \frac{1}{B_i}}\left[B^{}\right]_i+{\displaystyle \frac{2}{A_i}}\left[A^{}\right]_i={\displaystyle \frac{2B^{}(y_i^+)}{B_i}}+{\displaystyle \frac{4A^{}(y_i^+)}{A_i}}.`$ (13)
Meanwhile, the 5D conservation law $`T_{A;B}^B=0`$ gives the usual 4D equation of conservation on the branes:

$$\dot{\rho }_i+3\left(\rho _i+p_i\right)\frac{\dot{A}_i}{A_i}=0.$$
(14)
We also need to define the Hubble and deceleration parameters on the branes appropriately. Be aware that the proper time on a given $`y=`$constant hypersurface is defined by $`d\tau =Bdt`$, so, with use of the relation (7), the Hubble and deceleration parameters are defined as (see )

$$H(t,y)\frac{1}{B}\frac{\dot{A}}{A}=\frac{\mu }{A},q(t,y)=\left[\frac{A}{B}\frac{}{t}\left(\frac{\dot{A}}{B}\right)\right]\left(\frac{\dot{A}}{B}\right)^2=\frac{A\dot{\mu }}{\mu \dot{A}}.$$
(15)
So, on the branes, we have

$$H_i=\frac{\mu }{A_i},q_i=\frac{A_i\dot{\mu }}{\mu \dot{A}_i}.$$
(16)
Using (16) and (13) in (10), we obtain

$$H_i^2+\frac{k}{A_i^2}=\frac{\kappa ^4}{36}\rho _i^2+\frac{\kappa ^2}{6}\mathrm{\Lambda }+\frac{K}{A_i^4}.$$
(17)
Thus we recover the induced Friedmann equation on the branes that was discussed and studied widely in literature \[2-6,9,17\]. In what follows we consider the three cases for the brane models with $`\mathrm{\Lambda }=0`$, $`\mathrm{\Lambda }>0`$, and $`\mathrm{\Lambda }<0`$, respectively.

### II.1  TYPE I BRANE MODELS: $`\mathrm{\Lambda }=0`$

For $`\mathrm{\Lambda }=0`$, the general solutions of (10) are found to be

$$A^2=\left(\mu ^2+k\right)y^22\nu \left|y\right|+\frac{\nu ^2+K}{\mu ^2+k}.$$
(18)
This type of solutions contains two arbitrary functions $`\mu (t)`$ and $`\nu (t)`$ and two constants $`k`$ and $`K`$. Differentiating this $`A^2`$ with respect to $`y`$, we obtain

$$AA^{}=\left(\mu ^2+k\right)y\nu \frac{\left|y\right|}{y}.$$
(19)
Therefore we get

$`A^{}(0^+)`$ $`=`$ $`A^{}(0^{})={\displaystyle \frac{\nu }{A_1}},`$
$`A^{}(y_2^+)`$ $`=`$ $`A^{}(y_2^{})={\displaystyle \frac{(\mu ^2+k)y_2\nu }{A_2}},`$ (20)
where the first brane is at $`y=y_1=0`$ and the second is at $`y=y_2>0`$. Then, from (7) we get

$`B^{}(0^+)`$ $`=`$ $`B^{}(0^{})={\displaystyle \frac{1}{\mu }}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\nu }{A_1}}\right),`$
$`B^{}(y_2^+)`$ $`=`$ $`B^{}(y_2^{})={\displaystyle \frac{1}{\mu }}{\displaystyle \frac{}{t}}\left[{\displaystyle \frac{(\mu ^2+k)y_2\nu }{A_2}}\right].`$ (21)
So equation (13) gives

$`\kappa ^2\rho _1`$ $`=`$ $`{\displaystyle \frac{6\nu }{A_1^2}},`$
$`\kappa ^2p_1`$ $`=`$ $`{\displaystyle \frac{2}{\dot{A}_1}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\nu }{A_1}}\right){\displaystyle \frac{4\nu }{A_1^2}},`$ (22)
and

$`\kappa ^2\rho _2`$ $`=`$ $`{\displaystyle \frac{6}{A_2}}\left({\displaystyle \frac{(\mu ^2+k)y_2\nu }{A_2}}\right),`$
$`\kappa ^2p_2`$ $`=`$ $`{\displaystyle \frac{2}{\dot{A}_2}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{(\mu ^2+k)y_2\nu }{A_2}}\right){\displaystyle \frac{4}{A_2}}\left({\displaystyle \frac{(\mu ^2+k)y_2\nu }{A_2}}\right).`$ (23)
### II.2 TYPE II BRANE MODELS: $`\mathrm{\Lambda }>0`$

For $`\mathrm{\Lambda }>0`$, the general solutions of (10) can easily be found. Then, by changing $`y`$ to $`\left|y\right|`$, we obtain

$`A^2`$ $`=`$ $`a(t)\mathrm{cos}\left(\theta (t)+L^1\left|y\right|\right)+2L^2\left(\mu ^2+k\right),`$
$`a(t)`$ $`=`$ $`2L\sqrt{L^2\left(\mu ^2+k\right)^2K},L^1\sqrt{{\displaystyle \frac{2\kappa ^2\mathrm{\Lambda }}{3}}},`$ (24)
where $`\mu (t)`$ and $`\theta (t)`$ are two arbitrary functions. By differentiating $`A^2`$ with respect to $`y`$, we obtain

$$2AA^{}=L^1a(t)\mathrm{sin}\left(\theta +L^1\left|y\right|\right)\frac{\left|y\right|}{y}.$$
(25)
So we have

$`A^{}(0^+)`$ $`=`$ $`A^{}(0^{})={\displaystyle \frac{a\mathrm{sin}\theta }{2LA_1}},`$
$`A^{}(y_2^+)`$ $`=`$ $`A^{}(y_2^{})={\displaystyle \frac{a\mathrm{sin}\left(\theta +L^1y_2\right)}{2LA_2}},`$ (26)
Using these and the relation $`B=\dot{A}/\mu `$, we get

$`B^{}(0^+)`$ $`=`$ $`B^{}(0^{})={\displaystyle \frac{1}{2L\mu }}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sin}\theta }{A_1}}\right),`$
$`B^{}(y_2^+)`$ $`=`$ $`B^{}(y_2^{})={\displaystyle \frac{1}{2L\mu }}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sin}\left(\theta +L^1y_2\right)}{A_2}}\right),`$ (27)
So equation (13) gives

$`\kappa ^2\rho _1`$ $`=`$ $`{\displaystyle \frac{3a\mathrm{sin}\theta }{LA_1^2}},`$
$`\kappa ^2p_1`$ $`=`$ $`{\displaystyle \frac{1}{L\dot{A}_1}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sin}\theta }{A_1}}\right){\displaystyle \frac{2a\mathrm{sin}\theta }{LA_1^2}},`$ (28)
and

$`\kappa ^2\rho _2`$ $`=`$ $`{\displaystyle \frac{3a\mathrm{sin}\left(\theta +L^1y_2\right)}{LA_2^2}},`$
$`\kappa ^2p_2`$ $`=`$ $`{\displaystyle \frac{1}{L\dot{A}_2}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sin}\left(\theta +L^1y_2\right)}{A_2}}\right)+{\displaystyle \frac{2a\mathrm{sin}\left(\theta +L^1y_2\right)}{LA_2^2}}.`$ (29)
### II.3 TYPE III BRANE MODELS: $`\mathrm{\Lambda }<0`$

For $`\mathrm{\Lambda }<0`$, the general solutions of (10) are found to be

$`A^2`$ $`=`$ $`a(t)\mathrm{cosh}\left(\xi (t)+L^1\left|y\right|\right)2L^2\left(\mu ^2+k\right),`$
$`a(t)`$ $`=`$ $`2L\sqrt{L^2\left(\mu ^2+k\right)^2+K},L^1\sqrt{{\displaystyle \frac{2}{3}}\kappa ^2\mathrm{\Lambda }},`$ (30)
where $`k`$ and $`K`$ and $`L`$ are constants, $`\mu =\mu (t)`$ and $`\xi =\xi (t)`$ are two arbitrary functions. Differentiating $`A^2`$ with respect to $`y`$, we obtain

$$2AA^{}=L^1a(t)\mathrm{sinh}\left(\xi +L^1\left|y\right|\right)\frac{\left|y\right|}{y}.$$
(31)
So we find

$`A^{}(0^+)`$ $`=`$ $`A^{}(0^{})={\displaystyle \frac{a\mathrm{sinh}\xi }{2LA_1}},`$
$`A^{}(y_2^+)`$ $`=`$ $`A^{}(y_2^{})={\displaystyle \frac{a\mathrm{sinh}\left(\xi +L^1y_2\right)}{2LA_2}},`$ (32)
Then, using these and the relation $`B=\dot{A}/\mu `$, we find

$`B^{}(0^+)`$ $`=`$ $`B^{}(0^{})={\displaystyle \frac{1}{2L\mu }}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sinh}\xi }{A_1}}\right),`$
$`B^{}(y_2^+)`$ $`=`$ $`B^{}(y_2^{})={\displaystyle \frac{1}{2L\mu }}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sinh}\left(\xi +L^1y_2\right)}{A_2}}\right).`$ (33)
So the equation (13) gives

$`\kappa ^2\rho _1`$ $`=`$ $`{\displaystyle \frac{3a\mathrm{sinh}\xi }{LA_1^2}},`$
$`\kappa ^2p_1`$ $`=`$ $`{\displaystyle \frac{1}{L\dot{A}_1}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sinh}\xi }{A_1}}\right)+{\displaystyle \frac{2a\mathrm{sinh}\xi }{LA_1^2}},`$ (34)
and

$`\kappa ^2\rho _2`$ $`=`$ $`{\displaystyle \frac{3a\mathrm{sinh}\left(\xi +L^1y_2\right)}{LA_2^2}},`$
$`\kappa ^2p_2`$ $`=`$ $`{\displaystyle \frac{1}{L\dot{A}_2}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{a\mathrm{sinh}\left(\xi +L^1y_2\right)}{A_2}}\right){\displaystyle \frac{2a\mathrm{sinh}\left(\xi +L^1y_2\right)}{LA_2^2}}.`$ (35)
## III A GRAVITATIONAL STABILITY CONDITION FOR PARTICLES ON BRANES

Generally speaking, brane models require the fifth dimension to be compactified to a small size $`y_2`$. So a stable brane model should have a stable size. The well-known Goldberger-Wise mechanism  was to use a scalar field in the bulk to stabilize the size of the fifth dimension. In this way, the model has the tendency to recover it’s size after a perturbation. Here, a perturbation means a small change for the size.

Here we wish to say that even if the size of the fifth dimension is somehow stabilized, one still can ask question whether particles on the brane may leave the brane and escape into the bulk. Arkani-Hamed et al  pointed out that in sufficiently hard collisions, particles on the brane can acquire momentum in the extra dimensions and escape from our 4D world, carrying away energy. If this happens continuously, it will cause another kind of instability problem. Note that here the instability just means brane particles are not in a stable position along the fifth dimension. Now a natural question is that once a particle left the brane, will it escape in the bulk forever or return to our brane again?

To answer this question we noticed that the bulk discussed in this paper only contains a cosmological constant term $`\mathrm{\Lambda }`$ and so is empty. Therefore, it is reasonable to expect that particles inside the bulk should obey the 5D geodesic equation

$$\frac{d^2x^A}{dS^2}+\mathrm{\Gamma }_{BC}^A\frac{dx^B}{dS}\frac{dx^C}{dS}=0,$$
(36)
which was used previously  for a similar purpose as here. From this equation, a 5D gravitational force acted on the bulk test particle can be defined as

$$F^A=\mathrm{\Gamma }_{BC}^A\frac{dx^B}{dS}\frac{dx^C}{dS}.$$
(37)
For simplicity, we assume that the particle is temporary at rest in the bulk. Then, along the fifth dimension, (37) gives

$$F^5=\frac{B^{}}{B}.$$
(38)
In the vicinity of the $`i`$-th brane, we use the general results (13) and (12) in (38). Then we obtain

$$F^5(y_i^+)=F^5(y_i^{})=\frac{\kappa ^2}{3}\left(\rho _i+\frac{3}{2}p_i\right).$$
(39)
If this force $`F^5`$ is attractive, i.e., $`F^5(y_i^+)<0`$ and $`F^5(y_i^{})>0`$, the particle would be “pulled” back to our brane. Thus we obtain a new kind of stability condition for the $`i`$-th brane as

$$\rho _i+\frac{3}{2}p_i>0.$$
(40)
This is a reasonable condition which holds for ordinary matters including dark matter. However, if the brane matter contains a dark energy component such as a cosmological constant or quintessence, the condition (40) may or may not hold, depending on how much the dark energy is contained on the brane. For example, we let

$$\rho _i=\overline{\rho }_i+\overline{\lambda }_i,p_i=\overline{p}_i\overline{\lambda }_i,$$
(41)
where $`\overline{\lambda }_i`$ is a cosmological term on the $`i`$-th brane. Then the condition (40) becomes

$$\overline{\rho }_i+\frac{3}{2}\overline{p}_i>\frac{1}{2}\overline{\lambda }_i.$$
(42)
Suppose initially condition (42) holds, then, as the universe expands, both $`\overline{\rho }_i`$ and $`\overline{p}_i`$ decrease while $`\overline{\lambda }_i`$ keeps unchanged. So gradually the universe will enter in an unstable stage in which particles and energy on the brane may escape in the bulk and the 4D conservation law of energy (14) may not hold anymore. Recent observations  reveal that presently we are living in an accelerating stage dominated by a dark energy term such as $`\overline{\lambda }_i`$ in (42). Thus we see that both the acceleration of our universe and the unstable nature of the brane particles could be explained as due to the same repulsive force of dark energy. In the following section we will show that some brane models do not satisfy the condition (40).

## IV SOME SIMPLE BIG BOUNCE BRANE MODELS

In Section II we have obtained three types of global brane solutions corresponding to $`\mathrm{\Lambda }=0`$, $`\mathrm{\Lambda }>0`$ and $`\mathrm{\Lambda }<0`$, respectively. Each type contains two arbitrary functions: $`\mu (t)`$ and $`\nu (t)`$ for type I, $`\mu (t)`$ and $`\theta (t)`$ for type II, and $`\mu (t)`$ and $`\xi (t)`$ for type III. All the two scale factors $`A(t,y)`$ and $`B(t,y)`$ and the densities $`\rho _i`$ and $`p_i`$ on the branes are expressed via the two functions. From the relation $`B=\dot{A}/\mu `$ and metric (1) we see that the form of $`Bdt`$ is invariant under an arbitrary coordinate transformation $`t\stackrel{~}{t}(t)`$. This freedom could be used to fix one of the two arbitrary functions. Another freedom may correspond, as is in the standard general relativity, to the unspecified equation of state of matter. So, generally speaking, if the matter content on the first brane is known, then the two arbitrary functions could be fixed. Then the whole solutions could be fixed too. Then we will know the matter content on the second brane. Therefore, the brane solutions obtained in Sec. II are quite general.

To compare these cosmological solutions with observations, we need know clearly the matter content on our brane. This may need careful analysis and might be complicated mathematically, and we are not going to do it here in this paper. However, as an illustration, we will pick up several simple exact models in the following and to exhibit typical features of brane cosmology.

### IV.1 A SIMPLE TWO-BRANE MODEL WITH $`\mathrm{\Lambda }=0`$

Consider the $`\mathrm{\Lambda }=0`$ brane solutions (18), for which we choose

$$k=0,K=1,\nu (t)=t_b/t,\mu (t)=\left(2t\right)^{1/2},$$
(43)
where $`t_b>0`$ is a constant. In this way, the solution (18) becomes

$`A^2`$ $`=`$ $`2t\left[1+\left({\displaystyle \frac{\left|y\right|2t_b}{2t}}\right)^2\right],`$
$`B^2`$ $`=`$ $`\left[1\left({\displaystyle \frac{\left|y\right|2t_b}{2t}}\right)^2\right]^2\left[1+\left({\displaystyle \frac{\left|y\right|2t_b}{2t}}\right)^2\right]^1.`$ (44)
So on the first brane we have

$`A_1^2`$ $`=`$ $`2t\left[1+\left({\displaystyle \frac{t_b}{t}}\right)^2\right],B_1^2=\left[1\left({\displaystyle \frac{t_b}{t}}\right)^2\right]^2\left[1+\left({\displaystyle \frac{t_b}{t}}\right)^2\right]^1,`$
$`\kappa ^2\rho _1`$ $`=`$ $`{\displaystyle \frac{3t_b}{t^2+t_b^2}},\kappa ^2p_1={\displaystyle \frac{2t_b}{t^2t_b^2}}{\displaystyle \frac{t_b}{t^2+t_b^2}}.`$ (45)
On the second brane we have

$`A_2^2`$ $`=`$ $`2t\left[1+\left({\displaystyle \frac{y_22t_b}{2t}}\right)^2\right],B_2^2=\left[1\left({\displaystyle \frac{y_22t_b}{2t}}\right)^2\right]^2\left[1+\left({\displaystyle \frac{y_22t_b}{2t}}\right)^2\right]^1`$
$`\kappa ^2\rho _2`$ $`=`$ $`{\displaystyle \frac{6\left(y_22t_b\right)}{4t^2+\left(y_22t_b\right)^2}},\kappa ^2p_2={\displaystyle \frac{4\left(y_22t_b\right)}{4t^2\left(y_22t_b\right)^2}}{\displaystyle \frac{2\left(y_22t_b\right)}{4t^2+\left(y_22t_b\right)^2}}.`$ (46)
It is easy to see for the first brane that there is a critical time $`t=t_b>0`$ at which the scale factor $`A_1(t)`$ reaches to a non-zero minimum $`2\sqrt{t_b}`$. After this critical time, $`A_1(t)`$ increases, implying the universe is expanding. For $`tt_b`$ , we have $`A_1(t)\sqrt{2t}`$, $`B_1(t)1`$, $`\rho _13p_1`$, and the universe evolves as if in the radiation-dominated standard Friedmann model. Before this critical time, the universe contracts from infinity. It was shown  that as the coordinate time $`t`$ varies from zero to $`t_b`$ , the proper time $`\tau `$ varies from $`\mathrm{}`$ to $`\tau _b`$ , where $`\tau _b`$ is a finite constant and corresponds to $`t_b`$ . Thus we can call $`t_b`$ as to represent a big bounce (as opposed to a big bang), and before this bounce the universe has existed for an infinitely long time. Meanwhile, rather than the big bang singularity, which should correspond to $`A_1(t)=0`$, this big bounce singularity corresponds to $`B_1(t)=0`$. From the solution (45) we can also see that as $`t`$ varies across $`t_b`$ , the energy density $`\rho _1`$ remains finite while the pressure $`p_1`$ changes from $`\mathrm{}`$ to $`+\mathrm{}`$ . This implies a phase transition happened at the bounce. Here we plot the evolution of the scale factor $`A_1(t)`$ with $`t_b=1`$ in Figure 1.

The evolution of the second brane is similar as the first one with the bounce point at $`t=t_{b2}\left|t_by_2/2\right|`$. From equation (46) we see that whether the energy density on the second brane is positive or negative depends on the “size” of the fifth dimension $`y_2`$. If $`y_2>2t_b`$, we have $`\rho _2>0`$. If $`y_2<2t_b`$, we have $`\rho _2<0`$.

Now let us look at the stability condition (40). Using (45) and (46) in (40), we obtain

$$\kappa ^2\left(\rho _1+\frac{3}{2}p_1\right)=\frac{3t_b}{t^2t_b^2}+\frac{3t_b}{2\left(t^2+t_b^2\right)},$$
(47)
$$\kappa ^2\left(\rho _2+\frac{3}{2}p_2\right)=\frac{6\left(y_22t_b\right)}{4t^2\left(y_22t_b\right)^2}+\frac{3\left(y_22t_b\right)}{4t^2+\left(y_22t_b\right)^2}.$$
(48)
For the first brane we find that $`\rho _1+\frac{3}{2}p_1>0`$ for $`t>t_b`$. So, in the expanding stage (after the bounce), the first brane is stable. There are two cases for the second brane. If $`y_2>2t_b`$, we have $`\rho _2>0`$ and $`\rho _2+\frac{3}{2}p_2>0`$ in the stage $`t>t_{b2}`$ where $`t_{b2}\left|y_22t_b\right|`$ represents the bounce time. So after the bounce this brane is stable. If $`y_2<2t_b`$, we have $`\rho _2<0`$ and $`\rho _2+\frac{3}{2}p_2<0`$ in the stage $`t>t_{b2}`$. So after the bounce this brane is unstable. Thus, to obtain a model with two stable branes, we must have $`y_2>2t_b`$. An interesting special case is $`y_2=4t_b`$ for which we have $`\rho _2=\rho _1>0`$ and the model is completely symmetric and particles on the branes are stable.

### IV.2 AN OSCILLATING UNIVERSE MODEL WITH $`\mathrm{\Lambda }>0`$

The notion of an oscillatory universe can be traced back to 1930’s and has continued to attract interest \[23-25,11\]. Tolman  discussed it within the framework of general relativity assuming a closed universe ($`k=+1`$) in which the universe undergoes a sequence of cycles of expansion and contraction. Dicke and Peebles  restudied Tolman’s model and pointed out that an oscillating universe could provide an escape from some of the cosmological problems such as the horizon and homogeneity problems. “As the universe ages”, they wrote, “more and more of it becomes visible that earlier was beyond the horizon and presumably causally disconnected from us …. How then are we to understand the remarkable familiarity of the objects just appearing on the horizon? Perhaps by tracing the evolution back through the big bang to an earlier collapsing phase”. However, Tolman’s oscillatory models were constrained by having to pass through the big bang singularity in which the energy density and temperature diverge. Therefore, Dicke and Peebles expected that “some future and better theory might show that the collapse of the universe would lead to a ‘bounce’ instead of a singularity”.

Now let us consider the type II ($`\mathrm{\Lambda }>0`$) brane solutions. As we discussed at the beginning of this section that the metric form of $`Bdt`$ is invariant under an arbitrary coordinate transformation $`t\stackrel{~}{t}(t)`$. Meanwhile, the scale factor $`A(t,y)`$ in (24) contains two arbitrary functions $`\mu (t)`$ and $`\theta (t)`$. So, by choosing the time coordinate properly, we can set $`\theta (t)=\pi t`$ without loss of generality. So the cosine term in the solution leads generally to an oscillating universe model. For illustration, we choose

$`\mu ^2(t)`$ $`=`$ $`t,\theta (t)=\pi t,`$
$`L`$ $`=`$ $`K=1,k=0.`$ (49)
Then

$`A^2`$ $`=`$ $`2t+2\sqrt{t^21}\mathrm{cos}\left(\pi t\left|y\right|\right)`$
$`a(t)`$ $`=`$ $`2\sqrt{t^21},{\displaystyle \frac{2\kappa ^2\mathrm{\Lambda }}{3}}=1,`$ (50)
So on the branes we have

$$A_i^2=2t+2\sqrt{t^21}\mathrm{cos}\left(\pi ty_i\right).$$
(51)
We plot the evolution of the scale factor $`A_1(t)`$ in Figure 2. From this figure we see clearly that the big bang spacetime singularity of the standard cosmology is replaced here by a series of smooth big bounces at $`t=t_b1`$, $`3`$, $`5`$, $`7`$, ….

Another typical feature of Figure 2 is that there is a “beginning of time” at $`t=1`$ in the model and the oscillating amplitudes of the followed cycles increase monotonously. This reminds us of the well-known entropy problem of the old oscillatory universe model. Tolman  pointed out that if the total entropy of the universe can only increase, then, in an oscillatory model, the entropy produced during one cycle would add to the entropy produced in the next cycle, causing the oscillating amplitude of each cycle to be larger than the one before it. Extrapolating backward in time, the universe would have to be finite cycles old. More discussions about the entropy problem can be found in literature . We see that our oscillating universe model described by equation (51) and Figure 2 coincides with Tolman’s description perfectly and, therefore, could provide a realized framework to discuss the entropy problem.

The two equations (28) and (29) become

$`\kappa ^2\rho _1`$ $`=`$ $`{\displaystyle \frac{6\sqrt{t^21}\mathrm{sin}\pi t}{A_1^2}},`$
$`\kappa ^2p_1`$ $`=`$ $`{\displaystyle \frac{2}{\dot{A}_1}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\sqrt{t^21}\mathrm{sin}\pi t}{A_1}}\right)+{\displaystyle \frac{4\sqrt{t^21}\mathrm{sin}\pi t}{A_1^2}},`$ (52)
and

$`\kappa ^2\rho _2`$ $`=`$ $`{\displaystyle \frac{6\sqrt{t^21}\mathrm{sin}\left(\pi ty_2\right)}{A_2^2}},`$
$`\kappa ^2p_2`$ $`=`$ $`{\displaystyle \frac{2}{\dot{A}_2}}{\displaystyle \frac{}{t}}\left[{\displaystyle \frac{\sqrt{t^21}\mathrm{sin}\left(\pi ty_2\right)}{A_2}}\right]{\displaystyle \frac{4\sqrt{t^21}\mathrm{sin}\left(\pi ty_2\right)}{A_2^2}}.`$ (53)
So the energy density on the branes changes sign periodically with time, implying that the universe may have a negative energy density. This is an unusual feature of the oscillatory brane model. It is reasonable to assume that the size the fifth dimension is much smaller than the period of each cycle of the universe, i.e., $`y_22\pi `$. Then from equations (51) - (53) we find $`A_2(t)`$ $`A_1(t)`$, $`\rho _2\rho _1`$ and $`p_2p_1`$. So if our brane has a positive energy density at present stage of the universe, the other brane would have a negative energy density. Applying this in the stability condition (40), we find that $`\left(\rho _2+\frac{3}{2}p_2\right)\left(\rho _1+\frac{3}{2}p_1\right)`$. This means that if one of the two branes is at a stable stage, the other one must be at an unstable stage. As a whole, we conclude that brane particles in this oscillating model are not stable.

### IV.3 A SIMPLE UNIVERSE MODEL WITH $`\mathrm{\Lambda }<0`$

For this type of brane solutions (30), we take a similar choice as in equations (49):

$`\mu ^2(t)`$ $`=`$ $`t,\xi (t)=t,`$
$`L`$ $`=`$ $`K=1,k=0.`$ (54)
Then we obtain

$`A^2`$ $`=`$ $`2\sqrt{t^2+1}\mathrm{cosh}\left(t\left|y\right|\right)2t,`$
$`A_i^2`$ $`=`$ $`2\sqrt{t^2+1}\mathrm{cosh}\left(ty_i\right)2t.`$ (55)
We plot the evolution of $`A_1(t)`$ of the first brane in Figure 3, from which we see that $`A_1(t)`$ reaches to its minimum at the bounce point $`t=t_b0.5`$. Before the bounce, the universe contracts from infinity; after the bounce, the universe expands to infinity again.

From equations (34) and (35) we obtain

$`\kappa ^2\rho _1`$ $`=`$ $`{\displaystyle \frac{6\sqrt{t^2+1}\mathrm{sinh}t}{A_1^2}},`$
$`\kappa ^2p_1`$ $`=`$ $`{\displaystyle \frac{2}{\dot{A}_1}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\sqrt{t^2+1}\mathrm{sinh}t}{A_1}}\right){\displaystyle \frac{4\sqrt{t^2+1}\mathrm{sinh}t}{A_1^2}},`$ (56)
and

$`\kappa ^2\rho _2`$ $`=`$ $`{\displaystyle \frac{6\sqrt{t^2+1}\mathrm{sinh}\left(ty_2\right)}{A_2^2}},`$
$`\kappa ^2p_2`$ $`=`$ $`{\displaystyle \frac{2}{\dot{A}_2}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{\sqrt{t^2+1}\mathrm{sinh}\left(ty_2\right)}{A_2}}\right)+{\displaystyle \frac{4\sqrt{t^2+1}\mathrm{sinh}\left(ty_2\right)}{A_2^2}}.`$ (57)
Thus the energy density is positive on the first brane and negative on the second brane. If we assume $`\left|t\right|y_2`$, then we have $`\rho _2\rho _1`$ and $`p_2p_1`$. So probably we are living on the first brane. As for the stability condition (40), we also have $`\left(\rho _2+\frac{3}{2}p_2\right)\left(\rho _1+\frac{3}{2}p_1\right)`$ as in the above oscillating model. So particles in this model are also unstable.

## V DISCUSSION

In this paper we have derived, in Sec. II, a class of five-dimensional cosmological solutions with two 3-branes and with the fifth dimension being static and compactified on a small circle. The bulk contains only a cosmological constant $`\mathrm{\Lambda }`$. It is found that for all the three cases of $`\mathrm{\Lambda }`$ ($`\mathrm{\Lambda }=0`$, $`\mathrm{\Lambda }>0`$, $`\mathrm{\Lambda }<0`$) the solutions contain two arbitrary functions of time. One of these two freedoms might be explained as due to the unspecified time coordinate in the 5D metric (1), leaving another to account for various contents of the cosmic matter. In Section III we have used the 5D geodesic equations to study the gravitational force acted on a test particle in the vicinity of a brane. This force could be interpreted as generated by matters on the brane. By requiring this force being attractive and so to grip particles from escaping into the bulk, we have derived a physically reasonable gravitational stability condition as given in equation (40). For illustration and for simplicity we presented three simple exact models in Section IV by choosing the two arbitrary functions properly. From these simple models we found that the conventional space-time singularity “big bang” could be replaced in brane models by a matter singularity “big bounce” at which the “size” of the 3D space is finite and the energy density does not diverge, while the pressure diverges. This enable us to expect that in brane cosmological models the “history” of our universe could be traced back across the big bounce to a pre-existing phase. This is clearly of great interest and deserve more studies. The stability of brane particles of these simple models are also studied.

Here we want to discuss more about the oscillating universe models given in equations (50) - (53). As pointed out by Tolman that an oscillatory model could resolve the horizon and the homogeneity problems. However, the main difficulty of Tolman’s oscillatory model is having to pass through the big bang space-time singularity. Now brane models could remove the big bang singularity in a satisfactory way and thus rescued Tolman’s old model. Meanwhile, Tolman’s entropy problem also get resolved.

We should emphasis that the present work is exploratory. Be aware that the general brane solutions given in this paper contain two arbitrary functions. This would enable us to discuss more observations such as the acceleration of the universe . We leave these studies in the future.

###### Acknowledgements.

The author thanks Guowen Peng and Lixin Xu for useful discussions. This work was supported by NSF of P. R. China under Grant 10273004.


text:
# Does merger-induced core-collapse produce 𝛾-ray bursts in type Ib & Ic supernovae?

### Acknowledgments.

I thank Kaiyou Chen, Stirling Colgate and Chris Fryer for their support, guidance, conversations and insight, and Alexander Heger for timely comments. The work was performed under the auspices of the Department of Energy.

## References

1. Klebesadel, R. W., Strong, I. B. & Olson, R. A. Observations of Gamma-Ray Bursts of Cosmic Origin. Astrophys. J. 182, L85-L88 (1973).

2. Galama, T. J. et al. An unusual supernova in the error box of the $`\gamma `$-ray burst of 25 April 1998. Nature 395, 670-672 (1998).

3. Hjorth, J. et al. A very energetic supernova associated with the $`\gamma `$-ray burst of 29 March 2003. Nature 423, 423-426 (2003).

4. Saio, H. & Nomoto, K. Inward Propagation of Nuclear-burning Shells in Merging C-O and He White Dwarfs. Astrophys. J. 500, 388-397 (1998).

5. Lyne, A. G. et al. The discovery of a millisecond pulsar in the globular cluster M28. Nature 328, 399-400 (1987).

6. Chen, K., Middleditch, J. & Ruderman, M. Low mass X-ray binaries and millisecond pulsars in globular clusters. Astrophys. J. 408, L17-L20 (1993).

7. Alpar, M. A., Cheng, A. F., Ruderman, M. A. & Shaham, J. A new class of radio pulsars. Nature 300, 728-730 (1982).

8. Bailyn, C. D. Blue Stragglers and Other Stellar Anomalies: Implications for the Dynamics of Globular Clusters. Annu. Rev. Astron. & Astrophys. 33, 133-162 (1995).

9. Chen, K. & Leonard, P. J. T. Does the coalescence of white dwarfs produce millisecond pulsars in globular clusters? Astrophys. J. 411, L75-L78 (1993).

10. Nisenson, P., Papaliolios, C., Karovska, M. & Noyes, R. Detection of a very bright source close to the LMC supernova SN 1987A. Astrophys. J. 320, L15-L18 (1987).

11. Meikle, W. P. S., Matcher, S. J. & Morgan, B. L. Speckle interferometric observations of supernova 1987A and of a bright associated source. Nature 329, 608-611 (1987).

12. Nisenson, P. & Papaliolios, C. A Second Bright Source Detected near SN1987A. Astrophys J. 518, L29-L32 (1999).

13. Wang, L. et al. The Axisymmetric Ejecta of Supernova 1987A. Astrophys. J. 579, 671-677 (2002).

14. Plait, P. C., Lundqvist, P., Chevalier, R. A. & Kirshner, R. P. HST observations of the ring around SN 1987A. Astrophys. J. 439, 730-751 (1995).

15. Burrows, C. J. et al. Hubble Space Telescope Observations of the SN 1987A Triple Ring Nebula. Astrophys. J. 452, 680-684 (1995).

16. Chen, K. & Colgate. S. A. in Conf. Proc. Millisecond Pulsars. A Decade of Surprise (eds Fruchter, A. S., Tavani, M., Backer, D. C.) 72, 205 (Astron. Soc. Pacif., San Francisco, 1995).

17. Lubow, S. H. & Shu. F. H. Gas dynamics of semidetached binaries. II - The vertical structure of the stream. Astrophys. J. 207, L53-L55 (1976).

18. Crotts, A. P. S. & Heathcote, S. R. SN 1987A’s Circumstellar Envelope. II. Kinematics of the Three Rings and the Diffuse Nebula. Astrophys. J. 528, 426-435 (2000).

19. Chanan, G. A., Middleditch, J. & Nelson, J. E. The Geometry of the Eclipse of a Pointlike Star by a Roche-Lobe Companion. Astrophys. J. 208, 512-517 (1976).

20. Middleditch, J. et al. Rapid photometry of supernova 1987A: a 2.14 ms pulsar? New Astron. 5, 243-283 (2000).

21. Cen, R. A Possible Lateral Gamma-Ray Burst Jet from SN 1987A. Astrophys. J. 524, L51-L54 (1999).

22. Wijnands, R. An observational review of accretion-driven millisecond X-ray pulsars. Preprint at http://arXiv.org/astro-ph0309347 and references therein (2003).

23. Iwamoto, K. et al. A hypernova model for the supernova associated with the $`\gamma `$-ray burst of 25 April 1998. Nature 395, 672-674 (1998).

24. MacFadyen, A. I. & Woosley, S. E. Collapsars: Gamma-Ray Bursts and Explosions in ‘Failed Supernovae’. Astrophys. J. 524, 262-289 (1999).

25. Berger, E., Kulkarni, S. R., Frail, D. A. & Soderberg, A. M. A Radio Survey of Type Ib and Ic Supernovae: Searching for Engine Driven Supernovae. Preprint at http://arXiv.org/astro-ph/0307228, (2003).

26. Kulkarni, S. R. et al. Radio emission from the unusual supernova 1998bw and its association with the $`\gamma `$-ray burst of 25 April 1998. Nature 395, 663-669 (1998).

27. Frail, D. A. et al. Beaming in Gamma-Ray Bursts: Evidence for a Standard Energy Reservoir. Astrophys. J. 562, L55-L58 (2001).

28. Chakrabarty, D., Pivovaroff, M. J., Hernquist, L. E., Heyl, J. S. & Narayan, R. The Central X-Ray Point Source in Cassiopeia A. Astrophys. J. 548, 800-810 (2001).

29. Middleditch, J. & Kristian, J. A. A search for young luminous optical pulsars in extragalactic supernova remnants. Astrophys. J. 279, 157-161 (1984).

30. Hankins, T. H. Kern, J. S., Weatherall, J. C. & Eilek, J. A. Nanosecond radio bursts from strong plasma turbulence in the Crab pulsar. Nature 422, 141-143 (2003).