text:
# 1 Introduction

## 1 Introduction

Type IIn supernovae (SN IIn) with narrow H$`\alpha `$ emission, introduced as a separate family by Schlegel (1990), explode in a very dense CS environment. This is demonstrated by the presence of a strong CS H$`\alpha `$, high bolometric luminosity and strong broad H$`\alpha `$ emission powered by the CS interaction (Chugai 1990, 1992). An analysis of the optical effects of the CS interaction provides an efficient diagnostic tool for the CS density around SN IIn. The use of this probe led to the detection of the unusually dense CS environment around SN 1987F, (Chugai 1990), SN 1997ab (Salamanca et al. 1998), SN 1997cy (Turatto et al. 2000) and some other SN IIn. For SN 1997ab with the velocity of the CS gas of $`90`$ km s-1 the derived mass loss is enormous, $`10^2M_{}`$ yr-1 (Salamanca et al. 1998). The mechanism for such a tremendous mass loss rate is unclear; it exceeds the mass loss rate of the most extreme cases of a red supergiant superwind by at least a factor of ten.

The problem of the mechanism of a powerful mass loss rate by SN IIn presupernovae is becoming urgent in view of the recent results of the SN 1994W study, which show that the CS shell in this case was created by mass loss with an average rate of $`0.2M_{}`$ yr-1 and the enormous kinetic luminosity, two orders of magnitude greater than the radiative luminosity of a massive presupernova (Chugai et al. 2003). It was suggested there that such a powerful mass loss by presupernova was related to an explosive event about 1.5 yrs before the SN 1994W outburst. A specific feature of this supenova is the relatively high velocity of the CS gas ($`u10^3`$ km s-1) that eventually leads to the energy problem for the superwind mechanism.

Originally, the idea of explosive mass loss several years before the supernova explosion was proposed by Weaver and Woosley (1979) in connection with a possible strong Ne flash in the degenerate O/Ne/Mg core. Grasberg and Nadyozhin (1986) suggested explosive ejection of the presupernova envelope roughly 50 days prior to supernova explosion in order to account for the narrow lines in SN 1983K. However at present the explosive mass ejection by presupernovae is just a working hypothesis, especially keeping in mind that Woosley et al (2002) recently expressed doubts about the reality of their mechanism. In this respect, the study of signatures of the explosive mass ejection by presupernovae (large mass and energy of the CS shell and small age) can provide us with interesting information about poorly undersood phenomena occuring in presupernovae on the eve of a supernova explosion.

It was already noted that explosive mass ejection by presupernovae might occur in those SN IIn which show high CS velocity ($`1000`$ km s-1) and CS subordinate hydrogen and metal lines (Chugai et al. 2003). Apart from SN 1994W this type of SN IIn includes another well observed supernova SN 1995G (Pastorello et al. 2002). This supernova has a specific light curve interpreted as a result of CS interaction (Pastorello et al. 2002).

In the present paper we use the model of the bolometric light curve in order to extract information about the CS gas density. This will permit us to derive mass, energy, and age of the CS envelope and thus a conclusion concerning its origin. In section 2 we give a brief description of the model, in section 3 we explore model sensitivity to parameters and demonstrate the uncertainty of the parameters recovered in the case of SN 1997cy. In section 4 we model the light curve of SN 1995G and derive the mass of the CS envelope. The discussion of results and their relationship to the H$`\alpha `$ intensity and profile is presented in section 5.

The paper is based upon the photometry and spectra presented by Pastorello et al. (2002). We use the Hubble constant of $`70`$ km s-1 Mpc-1.

## 2 Model

The light curve model applied below was described earlier (Chugai 2001). Here we only briefly recapitulate its essential features. We consider the expansion of the supernova in the CS enevelope with a given density distribution and velocity in the thin shell approximation (Chevalier 1982a). We neglect the structure of the interaction region which consists of the two shock waves, forward and reverse, with a density peak in between. This structure may be described by self-similar solution under proper conditions (Nadyozhin 1981, 1985; Chevalier 1982b). We are interested in the more general case of density distributions with non-zero CS velocity when the self-similar solution is not applicable, so the dynamics of the thin shell will be calculated numerically.

The interaction of a supernova with its CS environment on a time scale greater than several days weakly depends on the initial epoch of the interaction. We assume the following: the interaction begins at the presupernova radius $`R_0`$ and the supernova expands freely ($`v=r/t`$), and its density distribution is an inner plateau with the outer power law density drop ($`\rho v^9`$). The numerical solution of the equation of motion with the final velocity of the CS gas provides the thin shell radius $`R(t)`$, and the relative velocities of the supernova and CS gas flows. These permit us to calculate the kinetic luminosity of the forward and reverse shock, which may be transformed into X-ray luminosities of both shocks and eventually into optical bolometric luminosity (Chugai 1992). The contribution of the luminosity supplied by the internal energy stored in the supernova during the explosion is calculated according to an analytical approximation (Arnett 1980, 1982). In our model the full light curve is a linear superposition of this luminosity and the interaction luminosity. This approach permits us to take into acount the radiation of the initial internal energy at the early epoch in a straighforward way.

The CS density is described by the density $`\rho _0`$ at the radius $`10^{15}`$ cm, or density parameter $`w=4\pi r^2\rho `$ at this radius, and by the power index ($`s`$) in the power law ($`\rho r^s`$). The extent of the CS envelope is characterized by the outer radius $`R_\mathrm{b}`$. With the mass ($`M`$) and energy ($`E`$) of the supernova envelope we have five parameters, which are constrained by the light curve, photospheric radius, supernova expansion velocity and phase of the late rapid light curve decline.

## 3 Parameter sensitivity and SN 1997cy

The sensitivity of the model to parameter variations is demonstrated by the models A1, A2, B, C, D (Fig. 1) with parameters shown in Table 1. The table presents beginning from the second column: supernova mass ($`M`$), energy ($`E`$), power index of the CS density distribution ($`s`$), density parameter ($`w_0`$) and density ($`\rho _0`$) at the radius $`10^{15}`$ cm, outer radius of the CS envelope ($`R_\mathrm{b}`$), and the mass of the CS envelope ($`M_{\mathrm{cs}}`$). The latter is the output value and not a free parameter. We adopt the CS velocity 1000 km s-1, presupernova radius $`R_0=1000R_{}`$, and 56Ni mass $`0.003M_{}`$ following the estimate for SN 1994W (Sollerman et al. 1998). The model A1 is adopted as a standard, with which all the other models are compared. For all the cases Fig. 1 shows the bolometric light curves and the thin shell velocities.

We begin with a brief discussion of the general properties of models. First, they all show a shoulder, which reflects the overtaking of the CS shell boundary by the forward shock. Second, after that time the luminosity is still high. The source of this radiation is the inner shock which is driven by the supersonic velocity jump between the outer supernova material and the thin shell. As the thin shell is accelerated, the inner shock luminosity decreases. Note that the luminosity of the inner shock is lower in the case of a lower mass of the supernova envelope (model B). Finally, for most models, the contribution of the internal energy of the supernova to the luminosity is relatively small, except for the model A2 with a noticeable early bump during the first $`50100`$ days.

Now let us consider the effects of parameter variations. The reduction of the supernova kinetic energy by an order of magnitude reduces the early luminosity by more than an order of magnitude with a less pronounced luminosity decrease at the late epoch (Fig. 1a). The five-fold mass reduction (model B) results in the substantially higher early luminosity due to the higher velocity of the outer shock (Fig. 1b). Still the model B shows a faster decay of the early luminosity because of the fast crossing of the CS envelope. The model C with the steeper density distribution ($`s=2`$) for a similar mass of the CS shell gives, as expected, a higher early luminosity and faster decay (Fig. 1c). The model D with twice as higher mass of the CS envelope results in the higher luminosity and lower velocity of the thin shell (Fig. 1d).

The above analysis implies that a combination of low mass and low energy is able to produce the rather bright SN IIn phenomenon as a result of the efficient deceleration of the bulk of the supernova envelope. As an ilustration of this statement let us consider the light curve of SN 1997cy (type IIn), studied by Turatto et al. (2000). In the cited paper the proposed interaction model suggests the large supernova energy ($`3\times 10^{52}`$ erg) and mass ($`20M_{}`$). Although this model is plausible, one cannot rule out an alternative model with moderate values of both energy and mass. This is demonstrated by models cy1 and cy2 (Fig. 2 and Table 1). The adopted velocity of the CS envelope is 10 km s-1. Model cy1 with $`M=5M_{}`$ and $`E=2\times 10^{51}`$ erg provides a satisfactory agreement with data, while the model cy2 with $`M=1.5M_{}`$ and $`E=10^{51}`$ erg shows an even slighly better fit since it better reproduces the luminosity drop at $`t>600`$ days. The velocity of the thin shell is different in these models and such a difference may be crucial in discarding inappropriate models on the basis of the velocity information provided by the line profiles. Unfortunately, there is no straighforward procedure for determining the thin shell velocity in SN 1997cy. So the uncertainty of the choice of the mass and energy remains.

## 4 SN 1995G light curve and the mass of CS envelope

Let us briefly consider additional observational constraints on the model apart from the bolometric light curve of SN 1995G. The energy distribution in early spectra of SN 1995G provides the estimates of temperature and radius of the photosphere on days 2 and 36 (Pastorello et al. 2002). We believe that the photospheric radius at the early epoch is approximately equal to the radius of the thin shell. The arguments are based upon the result that the thin dense shell is opaque in the optical band in the case of a very dense CS environment ($`w10^{17}`$ g cm-1) for about one-two months (Chugai 2001). Thus, the first additional observational constraint on the model suggests that the early photospheric radius should approximately coincide with the radius of the thin shell, if the CS density is rather high, i.e. $`w10^{17}`$ g cm-1.

The next constraint of the model is provided by the expansion velocity of the thin shell. The observational information about this velocity should be extracted from the line profiles, particularly from the maximal velocity of the broad component. Since at the early epoch Thomson scattering can contribute in the broad wings (Chugai et al. 2003), in order to estimate the thin shell velocity we rely on the late nebular spectra of SN 1995G on days 265 and 561. In Fig. 3 the Fe II 5018 Å line is shown for both epochs (Pastorello et al. 2002). This line is free from blending which makes it a reliable indicator of the broad component velocity. The broad component is identified with the dense thin shell, possibly partially fragmented, at the boundary between the supernova and the CS gas as in SN 1994W (Chugai et al. 2003). The maximal velocity in the blue wing of the broad component is estimated using a simple procedure of the linear approximation of the line and continuum flux (Fig. 3). The derived maximal velocities are 3000 km s-1 and 2700 km s-1 on days 265 and 561, respectively, with a possible uncertainty of 10%. We attribute these velocities to the thin shell. Similar values ($`3000`$ km s-1) are shown by the maximal velocity in the H$`\alpha `$ blue wing. However, in this case Fe II emission lines may contribute to the blue wing flux, which may hamper the reliability of the estimate of the maximal velocity.

The expansion velocity of the CS envelope estimated from narrow absorption lines is $`u=750`$ km s-1 according to Pastorello et al. (2002). However, a somewhat larger value $`u850`$ km s-1 is demonstrated by the emission lines of the infrared triplet of Ca II (Pastorello et al. 2002, their Fig. 10). We adopt here $`u=800`$ km s-1.

Finally, yet another constraint on the model is the observation that after day 700 the light curve shows a more rapid decay which is interpreted as the result of the overtaking of the outer boundary of the CS envelope by the forward shock (Pastorello et al. 2002). This fact will be used to estimate the outer radius of the CS envelope.

Preliminary computations of an extended set of the expansion dynamics and bolometric light curves for SN 1995G reveal the following important feature: it emerges that within empirical constraints the model is not sensitive to either of the two guiding supernova parameters, mass and energy. Selecting mass as a primary parameter, we found that a fit is achievable for a wide range of mass values. Two models with masses $`2M_{}`$ (model G1) and $`10M_{}`$ (model G2), with other parameters given in Table 1, show an acceptable fit of the bolometric light curve, photospheric radius and thin shell velocity (Fig. 4). Note that to reach agreement with the photospheric radius in the first epoch we added 20 days to the age of the supernova given by Pastorello et al. (2002) which means that the explosion of SN 1995G is assumed to occur 20 days earlier than the zero point accepted in the cited paper. However, when some observational phase is mentioned in the text, we retain formally its day according to Pastorello et al. (2002). The models G1 and G2 confirm that the density and mass of the CS envelope do not depend on the adopted supernova mass (Table 1). The density at the radius of $`10^{15}`$ cm is $`9\times 10^{15}`$ g cm-3, or in terms of the hydrogen concentration for a normal abundance, $`n4\times 10^9`$ cm-3.

The independence of the CS gas density on the adopted supernova mass has a simple explanation. For the high CS density required in case of SN 1995G both shock waves (forward and reverse) are essentially radiative during a long period ($`700`$ d) which means that for a given total radiated energy the amount of the dissipated kinetic energy in shocks is invariant. Since the forward shock dominates in the luminosity, the latter suggests that the overall radiated energy should be of the order of the average value of $`0.5M_{\mathrm{cs}}(vu)^2`$, which thus must be invariant. Given observational constraints on the shell velocity ($`v`$) and CS velocity ($`u`$), the total mass of the CS envelope, thus, must also be invariant in different models.

In both models the kinetic energy is lower than the typical energy of core collapse supernovae ($`10^{51}`$ erg). Although the question of the typical value of the energy for SN IIn is open, it would be instructive to consider a case of a ”standard” energy. The model G3 with the energy $`E=10^{51}`$ erg shows a sensible fit and again requires the similar CS shell density and mass (Table 1) thus supporting the independence of these values on the adopted supernova mass at least in the range of $`220M_{}`$. Note, the latter statement can be reformulated in terms of supernova energy as a guiding parameter. In that case the derived CS density around SN 1995G is independent of the supernova energy at least in the range of $`(0.241)\times 10^{51}`$ erg.

The fact of the weak dependence of the CS gas density on the adopted supernova mass in the interaction model is of importance for the diagnostics of the CS density around SN IIn. The model of the bolometric light curve in combination with the velocity of the thin shell thus essentially provides us with a confident estimate of the CS density in the limit of the radiative forward shock wave. This tool was already used earlier for SN 1987F in an assumption of the standard supernova energy (Chugai 1992). Now it has become clear that in the limit of high CS density ($`w10^{17}`$ g cm-1) the derived CS density practically does not depend on the adopted supernova mass, at least in the range of $`2M20M_{}`$.

## 5 Discussion

The CS gas density distribution found above has direct implications for the interpretation of the spectrum of SN 1995G. The similarity of early spectra of SN 1995G and SN 1994W, particularly the presence of the strong effect of Thomson scattering in the H$`\alpha `$ profile (Chugai et al. 2003) suggests that a substantial fraction of H$`\alpha `$ is emitted by the CS enevlope. Is this picture consistent with the above CS density estimate?

The H$`\alpha `$ luminosity of the CS enevelope in the range of $`r_1<r<r_2`$ in our model ($`\rho r^2`$) is

$$L(\mathrm{H}\alpha )=\frac{1}{4\pi r_1}\alpha _{32}h\nu _{23}(xXwN_\mathrm{A})^2\left(1\frac{r_1}{r_2}\right),$$
(1)
wher $`\alpha _{32}`$ is the effective recombination coefficient of H$`\alpha `$ emissivity, $`h\nu _{23}`$ is the energy of the H$`\alpha `$ photon, $`x`$ is the ionization degree, $`X`$ is the hydrogen abundance, $`N_\mathrm{A}`$ is the Avogadro number. Substituting in equation (1) the values for day 2, namely, the inner radius $`r_1`$ equal to the radius of the photosphere $`1.1\times 10^{15}`$ cm (Pastorello et al. 2002), outer radius $`r_2=R_\mathrm{b}=2\times 10^{16}`$ cm, $`\alpha _{32}=1.2\times 10^{13}`$ cm3 s-1 (for electron temperature $`10^4`$ K), $`w=1.2\times 10^{17}`$ g cm-1, one obtains $`L(\mathrm{H}\alpha )=6.2\times 10^{40}x^2`$ erg s-1 assuming $`X=0.7`$. On the other hand, the observed H$`\alpha `$ luminosity on day 2 is $`L(\mathrm{H}\alpha )=6.3\times 10^{40}`$ erg s-1 (Pastorello et al. 2002). Following the SN 1994W case we assume that the contribution of the CS component in H$`\alpha `$ is at least half of the total line luminosity (Chugai et al. 2003). Comparing the model and observational luminosity of the CS component we thus conclude that both values agree, if the average ionization degree of the model CS envelope is $`0.7<x1`$.

The optical depth of the CS envelope to Thomson scattering for the same epoch is $`\tau _\mathrm{T}=k_\mathrm{T}wx/4\pi r_12.6x`$. For the above range of the ionization degree this gives $`1.8<\tau _\mathrm{T}<2.6`$. The presence of the strong red wing in H$`\alpha `$ on day 2 (Pastorello et al. 2002, their Fig. 10) suggests significant Thomson scattering, so one expects $`\tau _\mathrm{T}>1`$. Moreover, close similarity of this profile with the H$`\alpha `$ in SN 1994W on day 30 (Chugai et al. 2003, their Fig. 10) suggests that the optical depth of the CS envelope in SN 1995G is conceivably as large as $`\tau _\mathrm{T}2`$, which is consistent with the above range of Thomson optical depth.

To summarize, the density of the CS envelope derived from the interaction model agrees with both the H$`\alpha `$ line luminosity and the presence of strong effects of Thomson scattering in this line.

Generally, analysis of CS absorption lines might provide us additional information about the CS density. However this approach would require a rather complicated model of ionization and excitation in the CS envelope. Simple considerations based, for instance, on Fe II absorption, provide a rough lower limit (Sollerman et al. 1998). The requirement that absorption line with the optical depth $`\tau `$ is present in the spectrum is $`\tau >1`$. For an envelope of a size $`r`$ with the average velocity dispersion on the line of sight of the order of the expansion velocity $`v`$ (greater than thermal velocity), the latter condition is

$$\tau \sigma (\nu )n_1r\sigma _0n_1r\lambda _{12}/v>1,$$
(2)
where $`\sigma _0=(\pi e^2/m_ec)f_{12}`$ is the integrated over frequency cross-section $`\sigma (\nu )`$, $`f_{12}`$ is the oscillator strength, $`n_1`$ is the concentration at the lower level, $`\lambda _{12}`$ is the wavelength. For the Fe II 5018 Å absorption ($`f_{12}=0.01`$) adopting on day 2 the excitation temperature of the lower level equal to the photospheric one, 8800 K (Pastorello et al. 2002), and taking $`r`$ equal to the photospheric radius ($`1.1\times 10^{15}`$ cm) one gets the lower limit of the hydrogen concentration assuming a solar Fe abundance $`n>3\times 10^7`$ cm-3 in qualitative agreement with the density found from the light curve analysis.

The outer boundary of the CS enevelope ($`2\times 10^{16}`$ cm) combined with the CS expansion velocity 800 km s-1 implies the age of the CS envelope $`t_{\mathrm{cs}}8`$ yr which is close to the estimate of the starting time for the strong mass loss, $`12`$ yr before the supernova explosion found by Pastorello et al. (2002). The mass of the CS envelope ($`1M_{}`$) combined with the age thus suggests the average mass loss rate $`\dot{M}0.1M_{}`$ yr-1 which is an enormous value. The estimated total kinetic energy of the CS enevelope is $`E_{\mathrm{cs}}6\times 10^{48}`$ erg and the average kinetic luminosity of the mass loss is then $`E_{\mathrm{cs}}/t_{\mathrm{cs}}2.4\times 10^{40}`$ erg s-1. This value is almost a factor of two in excess of the typical radiative luminosity of a massive presupernova ($`10^5L_{}`$). Thus the mass loss certainly cannot be attributed to the superwind.

We propose, therefore, that the mass ejection in the SN 1995G presupernova was initiated by some powerful energy release in the hydrodynamic time scale approximately 8 years before the major supernova explosion. In fact we are reproducing here arguments used in the case of SN 1994W to conclude that the CS enevelope around SN 1994W was lost as a result of an explosive event $`1.5`$ yr before the supernova explosion (Chugai et al. 2003). It was proposed there that the explosive mass ejection was initiated by the flash of nuclear burning of Ne in the degenerate O/Ne/Mg core. This assumption follows the original hypothesis of Weaver and Woosley (1979) concerning the behavior of presupernovae with initial masses of $`11M_{}`$. A similar possibility might occur also in the case of SN 1995G. Note that the age of the CS envelope around SN 1995G ($`8`$ yr) lies within the range of the phase for the Ne burning in massive star cores $`110`$ years before the supernova explosion (Heger 1998).

If the initial mass of the SN 1995G presupernova was actually close to $`11M_{}`$ then given a neutron star of $`1.4M_{}`$ the supernova ejecta cannot exceed $`10M_{}`$. In that case our interaction model predicts that the supernova kinetic energy $`6\times 10^{50}`$ erg (Table 1). This is lower than the value $`1.5\times 10^{51}`$ erg adopted for SN 1994W (Chugai et al. 2003). It may well be that the differences of supernova energies and ages of CS shells of SN 1995G and SN 1994W are possibly related to slight differences in their initial masses or in slightly different evolutionary histories.

If the mass ejection of presupernova SN 1995G had an explosive nature then the CS envelope expansion regime must be close to a free expansion ($`ur`$) law, at least in the outer layers. In this respect the increase of velocity derived from CS absorption lines of Fe II between days 330 and 560 (Pastorello et al. 2002) is qualitatively consistent with the possible free expansion CS kinematics.

The envelope ejection with the mass of $`1M_{}`$ and kinetic energy of $`6\times 10^{48}`$ erg should be accompanied by the optical flash eight years before the explosion of SN 1995G. In a simple analytical model of the light curve (Chugai 1991) assuming presupernova radius $`100<R_0<1000R_{}`$ we estimate the absolute magnitude of flash maximum as $`12.5>M_V>13.5`$ mag. with the duration of the light curve of $`80120`$ days. For the distance of 63 Mpc the apparent magnitude of maximum for this flash is $`21.8>V>20.8`$ mag. The only available image of the host galaxy NGC1643 close to the suggested time of the presupernova flash is a UK Schmidt plate taken 12 Dec. 1982, 13 years before the SN 1995G explosion. Therefore this image does not constrain the presupernova flash history. Inspection of this plate reveals no object brighter than 20.5 mag. in J band.

## 6 Conclusion

We performed the modeling of the bolometric light curve and expansion dynamics of SN 1995G in the dense CS environment. As a result we obtained the density and mass of the CS envelope which do not depend on the adopted mass of the supernova envelope. The derived mass of the CS envelope combined with the velocity of the CS gas leads us to conclude that the CS envelope was ejected as a result of energetic hydrodynamical process eight years before the explosion of SN 1995G. We speculate that this mass ejection was initiated by the powerful flash of nuclear burning in the O/Ne/Mg core of presupernova following the earlier conjecture proposed in the case of SN 1994W.

We thank the UK Schmidt Unit for providing scanned images of Schmidt plates around NGC1643. The research by N. Chugai is supported by RFBR grant 01-02 16295.

References

1. Arnett W.D., Supernovae: A survey of current research, Ed.
M.J. Rees, R.J. Stoneham (Reidel, 1982), p. 221.
2. Arnett W.D., Astrophys. J. 237, 541 (1980).
3. Chevalier R.A., Astrophys. J. 259, 302 (1982a).
4. Chevalier R.A., Astrophys. J. 258, 790 (1982b).
5. Chugai N.N., Sov. Astr. Lett. 16, 457 (1990).
6. Chugai N.N., Sov. Astr. Lett. 17, 210 (1991)
7. Chugai N.N., Sov. Astr. 36, 63 (1992)
8. Chugai N.N., Mon. Not. R. astron. Soc. 326, 1448 (2001).
9. Chugai N.N., Blinnikov S.I., Lundqvist P. et al., Mon. Not. R. astron.
Soc. (in preparation) (2003).
10. Falk S.W. and Arnett D.W., Astrophys. J. Suppl. 33, 515 (1977).
11. Grasberg E.K., Nadyozhin D.K., Sov. Astr. Lett. 12, 68 (1986).
12. Heger A., The presupernova evolution of rotating massive stars,
Ph.D dissertation, MPA 1120 (1998).
13. Nadyozhin D.K., Preprint ITEP No. 1 (1981).
14. Nadyozhin D.K., Astrophys. Space. Sci. 112, 225 (1985).
15. Pastorello A., M. Turatto, S. Benetti et al., Mon. Not. R. astron. Soc.
333, 27 (2002).
16. Salamanca I., Cid-Fernandes R., Tenorio-Tagle G. et al., Mon. Not.
R. astron. Soc. 300, L17 (1998).
17. Schlegel E., Mon. Not. R. astron. Soc. 244, 269 (1990).
18. Sollerman J., Cumming R.J. and Lundqvist P., Astrophys. J. 493,
933 (1998).
19. Turatto M., Suzuki T., Mazzali P.et al., Astrophys. J. 534, L57 (2000).
20. Weaver T.A. and Woosley S.E., BAAS 11, 724 (1979).
21. Woosley S.E., Heger A. and Weaver T.A., Rev. Mod. Phys. 74,
1015 (2002).


format: text

text:
# Equilibration of Long Chain Polymer Melts in Computer Simulations

## I Introduction

A systematic investigation of the structure-property relations for polymeric systems by computer simulations requires the preparation of equilibrated melts of long, entangled chains. For temperatures well above the glass transition, this can, in principle, be achieved using sufficiently long molecular dynamics (MD) or Monte Carlo (MC) simulations. However since the longest relaxation of an entangled polymer melt of length $`N`$ scales at least as $`N^3`$, giving at least $`N^4`$ in cpu time, this method is only feasible for relatively short chain lengths. Depending on the detail of the model the longest chains which can presently be equilibrated by brute-force MD or MC simulations are on the order of 2-7 entanglements lengths $`N_e`$. Using united atom potentials for polyethylene in which one treats the carbon and its bonded hydrogen as a single united atom, it is possible to equilibrate high temperature $`(T450K)`$ melts for chains of the order of approximately $`2N_e`$ or about 200 monomers. For coarse-grained bead-spring models like the one employed here it is possible to study longer chain lengths, on the order of up to 500 monomers or approximately $`7N_e`$. However this reaches the very limits of present day fastest computers. An increase of cpu speed by an order of magnitude does even not allow for a doubling of the chain length. While the present chain lengths are sufficient to follow the dynamics well into the reptation regime, they are not long enough to study structure-property relations, such as the plateau modulus $`G_N^o`$. This would require chains of order $`1020N_e`$. The situation is even worse for branched polymer melts or polymers at interfaces, where equilibration times even for relatively short chains are prohibitively large to use brute force simulations to produce equilibrated melts.

Fortunately to produce equilibrated samples, there is no need to follow the slow physical dynamics of entangled polymer melts. One possibility is to temporarily turn off the excluded volume interactions and to allow the chains to pass through each other. In order to obtain a well-defined Monte Carlo algorithm, it is useful to combine this idea with parallel tempering. While it has been demonstrated that this is feasible in principle, the method has so far not proven particularly efficient, the main reason being the large amount of computer time spent on unphysical Hamiltonians. More conventional MC algorithms which can be used to equilibrate polymer melts include reptation moves (generalized slithering snake algorithms), configuration bias algorithms, and concerted rotation algorithms. Being relatively local, these methods work best for moderate chain lengths and densities. The complete equilibration of very long chain melts still requires long runs. The fastest method, the generalized slithering snake algorithm, theoretically scales as $`N^2`$. An alternative are algorithms, which are able to move large sections of a chain at once. The prototype for such methods is the pivot algorithm in which a monomer is chosen at random and one of the two segments formed by that partitioning is pivoted rigidly in a random direction about the unit. A Boltzmann weight is used to determine if the move is accepted or not. This method is highly efficient for single chains in a continuum, implicit solvent. Unfortunately, a direct application to dense melts is not feasible since any large scale conformational change is bound to violate the packing constraints imposed by the chain excluded volume. However it is possible to introduce a double bond crossing algorithm in a way that two new bonds or bridges are formed across a pair of chains, creating two new chains each substantially different from the original. For the bead spring model concerned here this can simply be done by cutting two bonds and introducing two new bonds as discussed in more detail in Sec. VII. For atomistic models, the bond length is much shorter then the diameter of a monomer, and it is necessary to construct new bridges between the chain segments. Karayiannis et al. have shown that by constructing two trimer bridges between the two properly chosen dimers along the backbone one can quickly equilibrate long chain polyethylene melts. For the two chains involved this is essentially a double-pivot move, where each chain experiences such a pivot rotation.

In addition to improved equilibration algorithms there is an obvious interest in methods for generating initial melt configurations which are as close to equilibrium as possible. In our previous work on polymer melts, we prepared most melts by first generating an ensemble of chains with the correct end-to-end distance and randomly placing them in the simulation cell. We then introduced a weak, non-diverging excluded volume potential in the form of e.g. a cosine potential, $`A(1+\mathrm{cos}(\pi r/2^{1/6}\sigma ))`$, between non-bonded monomers, where $`r`$ is the distance between two monomers and $`\sigma `$ is the bead diameter. The amplitude $`A`$ was then increased over a short time interval until the minimum distance between monomers was sufficient to switch on the Lennard-Jones potential without creating numerical instabilities. In the first part of this paper we confirm observations by Brown et al. that this method deforms the (longer) chains. As a consequence, long chain melts are not as well equilibrated as originally believed. The deformation turns out to be strongest at short length scales and completely relaxes away only after a time $`\tau _{max}(N)`$, which is the time for a chain to move its own size. This long time is needed because the local chain-chain topology can only relax by the slow physical dynamics. For linear chains, $`\tau _{max}(N)\tau _e(N/N_e)^3`$, where the entanglement time $`\tau _e=\tau _{\mathrm{Rouse}}(N_e)`$ is the Rouse time of a chain of length $`N_e`$ (experimentally the measured exponent is closer to 3.4 than 3 except for extremely long chains). This however was exactly what we tried to avoid by that approach. Not surprisingly, the longer the chain length $`N`$, the more severe the effect is. For relatively short chains, such as those investigated in our previous studies of the crossover from Rouse to reptation dynamics $`(N350)`$, the simulations were run long enough that the results are independent of the starting state. However for longer chains, which cannot be run long enough compared to the longest relaxation time, this simple procedure for generating starting states is inadequate.

In the present paper we report results from an extensive effort to prepare well-equilibrated melts of bead-spring polymers with chain lengths ranging from $`N=350`$ to $`N=7000`$. For this purpose we use both approaches outlined above. First we show how to avoid the local stretching with only minimal computational effort by a suitable modification of our standard method. Then we demonstrate the capacity of the double-pivot algorithm to dynamically equilibrate melts of medium sized chains of length up to $`N=700`$. The paper is organized as follows: In Sec. II we define the model. In Sec. III, we present some simple theoretical estimates of the structure of bead-spring chains with different intrinsic bending stiffness in dense melts. In Sec. IV, we discuss ways to characterize the single chain structure and extract suitable target functions for long chain melts from long simulations of short chain melts. In Sec. V we examine carefully the standard procedure used in the past to prepare melt configurations. We show that for long chains this method deforms the chains at short length scales and that these deformations relax completely only after the chains have moved their own size. In Sec. VI, we present the new pre-packing procedure, which significantly reduces the density fluctuations particularly for long chains. This combined with a gradual introduction of the excluded volume, results in well equilibrated long chain melts in a reasonable amount of cpu time. In Sec. VI we describe the double-pivot algorithm and compare results for the internal dimensions of the chain with those produced from pre-packing and a gradual introduction of the excluded volume. In Sec. VIII we test the two methods for chains with local bending rigidity. A brief summary of our main conclusions are given in Sec. IX.

## II Model

We use a coarse grained model in which the polymer is treated as a string of beads of mass $`m`$ connected by a spring. The beads interact with a purely repulsive Lennard-Jones excluded volume interaction,

$$\mathrm{U}_{\mathrm{LJ}}(r)=\{\begin{array}{ccc}4ϵ\left\{(\sigma /r)^{12}(\sigma /r)^6+\frac{1}{4}\right\}\hfill & & rr_c\hfill \\ 0\hfill & & rr_c\hfill \end{array}$$
(1)
cutoff at $`r_c=2^{1/6}\sigma `$. The beads are connected by a finite extensible non-linear elastic potential (FENE),

$$\mathrm{U}_{\mathrm{FENE}}(\mathrm{r})=\{\begin{array}{ccc}0.5\mathrm{kR}_\mathrm{o}^2\mathrm{ln}\left(1(\mathrm{r}/\mathrm{R}_\mathrm{o})^2\right)\hfill & \mathrm{r}\mathrm{R}_\mathrm{o}& \\ \mathrm{}\hfill & \mathrm{r}>\mathrm{R}_\mathrm{o}& \end{array}$$

in addition to the Lennard-Jones interaction. The model parameters are the same as in ref. , namely $`k=30ϵ/\sigma ^2`$ and $`R_o=1.5\sigma `$ unless otherwise noted. The temperature $`T=ϵ/k_B`$. The basic unit of time is $`\tau =\sigma (m/ϵ)^{1/2}`$. We performed constant volume simulations of monodisperse melts at a segment density $`\rho =0.85\sigma ^3`$. The temperature was kept constant by coupling the motion of each bead weakly to a heat bath with a local friction $`\mathrm{\Gamma }`$. Unless otherwise specified $`\mathrm{\Gamma }=0.5m/\tau `$. The equations of motion were integrated using a velocity Verlet algorithm with a time step $`\mathrm{\Delta }t=0.012\tau `$. The average bond length is $`<b^2>^{1/2}=0.97\sigma `$. The polymer melts studied consisted of $`M`$ chains of length $`N`$ beads. The chain lengths studied varied from $`N=25`$ to $`7000`$. The number of chains are each system is specified in Table I.

Equilibration algorithms for polymer melts which follow the physical dynamics encounter a strong growth of the necessary equilibration times with chain length. Per unit of time $`\tau `$, the simulation of a melt of $`M`$ chains of length $`N`$ requires a computational effort of $`10^2M\times N`$ particle updates. Since the longest relaxation time is of the order $`\tau _{max}\tau _e(N/N_e)^3`$, equilibration of our largest system with $`M=80`$ and $`N=7000`$ would require $`10^2M\times N(\tau _e/\tau )(N/N_e)^38\times 10^{17}`$ particle updates. Assuming a typical performance of $`2.5\times 10^5`$ particle updates per second, which is typical for this model on a single processor, the required cpu time for the equilibration is about $`10^5`$ years. On a modern parallel cluster such as the DEC alpha CPlant cluster at Sandia, the times are somewhat less. On 64 processors, we get about 20 steps per second for our largest system, which means that the total time for the chains to move on average their own size is approximately 1600 years. Increasing the numbers of processors helps somewhat but since the efficiency of the computation decreases as the number of processors is increased for fixed system size, this type of brute force approach for long chain melts is clearly not feasible.

In our previous studies, we considered only fully flexible chains. However to demonstrate the effectiveness of the methodology, we also include a nonzero bending stiffness, which is modelled by

$$U_{bend}(\theta )=k_\theta \left(1\mathrm{cos}\theta \right)$$
(2)
where $`\mathrm{cos}\theta _i=(\widehat{𝐫}_i\widehat{𝐫}_{i1})(\widehat{𝐫}_{i+1}\widehat{𝐫}_i)`$, where $`(\widehat{𝐫}_i\widehat{𝐫}_{i1})`$ is the unit vector in the bond direction.

## III Chain characteristic in the melt

In dense polymer melts excluded volume interactions between different parts of a given polymer chain are screened. Since in our model there are no explicit intrachain interactions beyond those between neighboring bonds, the single chain structure should be well characterized by the expectation value $`\mathrm{cos}\theta `$ of the bond angle $`\theta `$. For example, the knowledge of $`\mathrm{cos}\theta `$ is sufficient to calculate an estimate for the mean-square end-to-end extension $`R^2`$ of chain segments of length $`N`$ monomers and $`n=N1`$ bonds,

$$R^2(n)=nb^2\left(\frac{1+\mathrm{cos}\theta }{1\mathrm{cos}\theta }\frac{1}{n}\frac{2\mathrm{cos}\theta \left(1\mathrm{cos}\theta ^n\right)}{\left(1\mathrm{cos}\theta \right)^2}\right)$$
(3)
where the asymptotic prefactor is referred to as

$$c_{\mathrm{}}=\frac{1+\mathrm{cos}\theta }{1\mathrm{cos}\theta }.$$
(4)
The simplest estimate of $`\mathrm{cos}(\theta )`$ only takes the bare bending energy Eq. (2) into account. Neglecting small variations in the bond length around the mean value $`b=0.97\sigma `$, one finds:

$`\mathrm{cos}\theta `$ $`=`$ $`{\displaystyle \frac{_0^\pi \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{exp}\left(\beta U_{bend}(\theta )\right)d\theta }{_0^\pi \mathrm{sin}\theta \mathrm{exp}\left(\beta U_{bend}(\theta )\right)d\theta }}`$ (5)
$`=`$ $`{\displaystyle \frac{1+e^{2\beta k_\theta }\left(\beta k_\theta 1\right)+\beta k_\theta }{\left(e^{2\beta k_\theta }1\right)\beta k_\theta }},`$ (6)
where $`\beta =1/k_BT`$. Eq. (5) is represented as a solid line in Fig. 1. However this result underestimates the effective chain stiffness, since the chain cannot fold back. This effect can be accounted for approximately by modelling the chains as non-reversal-random-walks (NRRWs) with excluded volume interactions between next-nearest neighbors along the backbone. As a function of the bond angle $`\theta `$ their distance is given by $`l_\theta =|\stackrel{}{r}_{i+1}\stackrel{}{r}_{i1}|=b[2(1+\mathrm{cos}\theta _i)]^{1/2}`$ so that

$`\mathrm{cos}\theta `$ $`=`$ $`{\displaystyle \frac{_0^\pi \mathrm{sin}\theta \mathrm{cos}\theta \mathrm{exp}\left(\beta (U_{bend}(\theta )+U_{LJ}(l_\theta ))\right)d\theta }{_0^\pi \mathrm{sin}\theta \mathrm{exp}\left(\beta (U_{bend}(\theta )+U_{LJ}(l_\theta ))\right)d\theta }}`$ (7)
which has to be evaluated numerically (see the black dots in Fig. 1). For $`k_\theta =0`$ this gives $`c_{\mathrm{}}=1.76`$ compared to $`c_{\mathrm{}}=1`$ obtained from Eq.( 3) and in good agreement with simulation results for chains of length $`20N350`$. This value is slightly smaller than the value $`c_{\mathrm{}}=2.06`$ obtained at the $`\mathrm{\Theta }`$-point for the same Lennard-Jones interaction truncated at $`2.5\sigma `$ in an implicit, continuum solvent.

Ensembles of single chains which obey Eq. (3) can be generated by simple sampling. This is particularly simple, if the allowed bond vectors are distributed evenly over a solid angle with $`\theta \theta _{max}`$. In this case

$$\frac{c_{\mathrm{}}1}{c_{\mathrm{}}+1}=\mathrm{cos}\theta _{\theta _{max}}=\left(\mathrm{cos}\theta _{max}/2\right)^2$$
(8)
so that the appropriate value for $`\theta _{max}`$ can be easily determined for any desired effective stiffness.

## IV Target functions and suitable ways of characterizing single chain structure

The estimates of the chain structure presented in the previous section were obtained within an effective single chain picture. In reality, the formal characterization of a polymer melt is a complicated many-body problem. For example, we have completely ignored the influence of packing effects on the chain conformations. As was recently shown the local melt structure very sensitively depends on the ratio of bond length to effective excluded volume of the beads. This effect is not well represented by the value of the overall chain extension. Here we use our standard bond length/diameter ratio and focus on the single chain properties, however the methods discussed below apply to systems with different parameters as well. In the present case, sufficiently long MD or MC simulations represent the only “ab initio” method which allows to take into account all interactions.

In order to characterize the chain conformations in the melt we mainly rely on mean-square internal distances $`R^2(|ij|,N)`$ averaged over all segments of size $`n=|ij|`$ along the chains, where $`i<j[1,N]`$ are monomer indices. In addition we characterize deviations from Gaussian statistics on short length scales using the ratio $`R^4/R^2^2`$. Figures 1 and 2 show a comparison of simulation results for equilibrated melts and our naive estimates. The results from the melts have been run sufficiently long that the chains have moved several times their own size. For $`N=350`$, $`k_\theta =0`$, the run time $`T_t=2.6\times 10^6\tau `$. Clearly, the description of chains in a melt as simple freely jointed chains is too naive. Nevertheless, the deviations are not large and the measured and estimated values for $`c_{\mathrm{}}`$ agree quite well. In particular, the simulation data show neither unexpected features nor significant finite chain length effects.

The rest of the present paper deals with algorithms which try to circumvent the slow physical equilibration path. These algorithms will be judged according to two criteria: their capacity to reproduce the target functions and their performance. We mostly concentrate on the fully flexible case. In the end, we come back to the stiffer systems as a kind of “blind test”.

## V Standard procedure for Preparing Polymer Melts

Our standard method for preparing melt conformations is as follows:

1. Start from an ensemble of chains with the correct end-to-end distance $`R^2(N)=b^2c_{\mathrm{}}(N1)`$ on large length scales. For a known $`c_{\mathrm{}}`$ this step is trivial.
2. Arrange the chains randomly in the simulation box.
3. Introduce a weak, non-diverging excluded volume potential. A convenient form for the soft potential is a cosine potential,

$$\mathrm{U}_{\mathrm{soft}}(r)=\{\begin{array}{cc}A(1+\mathrm{cos}(\pi r/r_c))\hfill & rr_c\hfill \\ 0\hfill & rr_c\hfill \end{array}$$

between non-connected monomers with $`r_c=2^{1/6}\sigma `$. The initial amplitude of $`A=4ϵ`$ was linearly increased to a final values of $`A100ϵ`$ (“push-off”) over a short time interval of 10-20 $`\tau `$. In the final state the inter monomer distances are large enough to allow one to switch to the LJ potential without creating numerical instabilities. During this phase, we often use a stronger coupling to the thermostat by increasing the friction constant to $`\mathrm{\Gamma }=2.0m/\tau `$. In addition, the velocities of all monomers are set to zero every $`50`$ time steps. This facilitates monomers that are strongly overlapping to separate. We refer to this method in which the excluded volume interactions are introduced fairly rapidly as “fast push-off.”
4. Relax the system with the full LJ-potential by a short MD run.

Figure 3 shows that this procedure actually deforms the chains so that the melt is not equilibrated after step 3. The deformations are strongest on short lengths. With increasing $`N`$ only the amplitude but not the position of the deformation along the chain changes. The monomer displacements during step 3 are too small to affect the conformations on large scales, allowing us to tune them to preselected values. Since the largest length scales also take the longest time to relax, one might hope that the MD relaxation in step 4 can be significantly shorter than for a starting conformation of, say, completely stretched chains. To be more specific, a naive expectation is that chain segments of length $`nN`$ equilibrate on time scales comparable to the Rouse/reptation time of chains of length $`n`$ independent of $`N`$. Figure 4 compares the result of equilibration runs of length $`T_t=10^4\tau \tau _{max}(N=90)`$ to the target function in Fig. 2. As expected, the internal distances measured for chains of length $`N=25`$ and $`N=50`$ coincide. However segments with $`n50`$ of longer chains for $`N=350`$ are still far from being equilibrated after $`10^4\tau `$. Even after a time $`t\tau _{\mathrm{Rouse}}(N)\tau _e(N/N_e)^210^5\tau `$ (not shown) the chains of length $`N=350`$ were not fully equilibrated. Rather the local equilibration is completed only after the chains have moved their own size, $`t>\tau _{max}(N)`$. In lattice simulations it is sometimes possible to avoid the slow reptation dynamics by allowing the chains to cut through each other. In this case $`\tau _{max}`$ is given by a Rouse-time and scales as $`N^2`$. In off-lattice bead-spring models topology conservation is the result of energetic barriers which prevent chain crossing. Since these barriers are a result of the microscopic interaction potentials, they are difficult to manipulate without affecting the local chain structure. Hence Rouse-like relaxation mechanisms due to barrier crossing dominate reptation only for extremely long (and therefore inaccessible) chain lengths. Attempts to circumvent these barriers using parallel tempering  have met with limited success at least as far as efficiency gains are concerned.

## VI Optimized Procedure for Polymer Melt Conformations

The methods presented in this section aim at a more careful implementation of the idea underlying our standard procedure: the local build-up of the characteristic melt packing for chains with the correct large length scale statistics. In this section we show that we can achieve this goal by first pre-packing the phantom chains and a subsequent slow and improved push-off scheme for the excluded volume. Both steps are needed to produce well equilibrated melts in a reasonable amount of cpu time. Applying only one does not achieve this objective. All results in this section are averaged over five independent realizations of the packing and push-off procedures. The system sizes studied are listed in Table I.

### A Prepacking of phantom chains

The chain deformations created by the standard procedure are due in part to the large local density fluctuations in a “melt” of randomly overlapping NRRW chains. As the chain length increases, these fluctuations increase in size as shown in Fig. 5. In this section we describe a dynamic Monte-Carlo-like pre-packing procedure for phantom chains, which significantly reduces these density fluctuations.

As a merit function we have chosen the fluctuations $`E=n^2n^2`$ in the number $`n_i`$ of particles found in a sphere of radius $`d`$ around particle $`i`$. In a homogeneously dense system $`E=0`$. We perform a zero temperature Monte Carlo (MC) optimization in which all moves which decrease the density fluctuations are accepted and all moves which increase the density fluctuations are rejected. In the course of the packing run $`d`$ is decreased from $`d=4\sigma `$ to $`2\sigma `$. We use a linked cell structure in order to calculate the $`n_i`$ efficiently. Choosing $`d>\sigma `$ effectively increases the range of the excluded volume correlations.

All MC moves change the positions or orientations of entire chains which are treated as rigid. This pre-packing procedure therefore does not affect the single chain statistics, which by construction already have the correct end-to-end distance $`R`$ as well as the targeted internal distance distributions. We consider five types of MC moves:

of individual chains in a random directions.

of individual chains by random angles around random axes through their centers of mass.

of individual chains at random mirror planes going through the center of mass.

of individual chains at their centers of mass.

of two chains preserving the center of mass positions.

Typical run times were of the order of days on individual workstations and therefore negligible.

The amplitude of the density fluctuations is drastically reduced as the packing proceeds. The final result is best assessed by the $`q0`$ limit of the system structure function $`S(q)`$. Figure 5 shows that the pre-packing reduces the density fluctuations to one percent of the initial value. Applying the fast push-off procedure outlined in Sec. V to these pre-packed states, we see by comparing Figs. 3 and 6a that the local deformations are eliminated for short chains $`(N50)`$ but not for long chains. For long chains, there is a significant reduction in the local stretching but not enough to avoid the necessity of subsequent long MD runs, which again would need to be of the order of the disentanglement time $`\tau _{max}(N)`$.

### B Slow push-off

In order to further reduce the chain deformations we have modified the push-off procedure in three ways:

1. Replace the cos-potential Eq. (3) with a force-capped-Lennard-Jones-potential, .
2. Keep the full excluded volume between next-nearest neighbor monomers to preserve to non-reversal-random-walk structure (see sec. III).
3. Increase the push-off time.

Force-capping  can either be defined directly via a maximum force $`F_{\mathrm{m}ax}`$ or implicitly via a distance $`r_{fc}`$ where $`\mathrm{U}_{\mathrm{LJ}}^{}{}_{}{}^{}(r_{fc})F_{\mathrm{m}ax}`$. At larger separations, the interaction is defined by Eq. (1), at smaller distances the force becomes constant so that the interaction potential grows only linearly with $`(rr_{fc})`$:

$$\mathrm{U}_{\mathrm{FCLJ}}(r)=\{\begin{array}{cc}(rr_{fc})\mathrm{U}_{\mathrm{LJ}}^{}{}_{}{}^{}(r_{fc})+\mathrm{U}_{\mathrm{LJ}}(r_{fc})& r<r_{fc}\\ \mathrm{U}_{\mathrm{LJ}}(r)& rr_{fc}\end{array}$$

In the present case, $`r_{fc}`$ is gradually reduced from the Lennard-Jones cut-off radius $`r_c=2^{1/6}\sigma `$ to $`0.8\sigma `$ which is significantly smaller than the relevant interparticle distances.

Force-capping has the advantage that the soft potential systematically approaches the true potential. In contrast, the cos-potential overcompensates the missing singularity at the origin by a fast rise of the repulsive interactions close to the cut-off distance. This leads to slight differences in the monomer packing. As a result, the chain conformations are slightly perturbed after switching to the full LJ-potential.

We have increased the push-off time in order to introduce the excluded volume interactions quasi-statically. In practice, one has in this context to worry about two issues: (i) How slow is slow enough? (ii) What is the equilibrium statistics for polymers with force-capped interactions? Concerning the first point, we have chosen a push-off time of $`5000\tau `$, which is of the order $`\tau _{\mathrm{Rouse}}(50)`$. This covers the typical subchain regime, up to which the intrachain distances still were disturbed by a fast push-off. The second question is more difficult to answer. If we simply switch off the repulsive monomer-monomer interaction, we find $`R^2=0.94\sigma ^2N`$, while the effective stiffness $`c_{\mathrm{}}1.7`$ of our chains in the melt is due to local bead packing effects. Thus a slow push off starting from a potential for which the end-to-end distance is an ideal random walk is also dangerous. For bead-spring chains these effects can easily be accounted for by treating the chains as NRRWs by always keeping the full LJ excluded volume between next-nearest neighbors along the chain. For atomistic polymer models, it is necessary to extend this to larger distances along the chain (“pentane effect”). Once this adjustement is made the push off procedure cannot be “too slow”and the chain size varies only very weakly with the maximum excluded volume force.

Figure 6b shows that the new procedure works quite well. Independent of the total chain length $`N`$ the chains are no longer stretched on intermediate scales and the large scales remain unaffected by the introduction of excluded volume interactions. As an additional check, we have varied the stiffness of the initially generated random walks over a range of about 20% (Fig. 7a). Fig. 7b compares the internal distances after pre-packing and slow push-off and confirms our previous observations. The chains are fully equilibrated on short scales, where the results are independent of the initial conditions. In contrast, internal distances beyond $`n=100`$ remain practically unaffected. We therefore conclude that our new method allows the preparation of well-equilibrated melts of extremely long chains within reasonable computational effort, provided the correct asymptotic chain stiffness $`c_{\mathrm{}}`$ is known from a careful extrapolation of simulation results for well-equilibrated short and medium chain length melts. For our largest system of $`M=80`$ chains of length $`N=7000`$, an equilibration time of $`5000\tau `$ takes about one month on a single processor or about 8 hours on the 100 processor DEC alpha cluster. Thus even with moderate computational means it is possible to equilibrate very large systems using a combined pre-packing and slow push-off procedure.

## VII Double-Pivot-MD Hybrid Algorithm

The present section deals with the problem of accelarating the dynamic equilibration of a long chain polymer melt. We describe a double-pivot-MD hybrid algorithm which performs this task significantly faster than simple MD relaxation. In contrast to the methods discussed in the previous section, the dynamic equilibration of a dense polymer melt does not require independent knowledge of the effective $`c_{\mathrm{}}`$. However, the preparation of locally equilibrated samples can significantly reduce the required relaxation times.

In the original pivot algorithm one choses a monomer at random and one of the two segments formed by that partitioning is pivoted rigidly in a random direction about the selected monomer. The pivot algorithm is extremely efficient when applied to isolated chains in an implicit solvent. However, in a melt or even in a solution at moderate density, pivot moves applied to a single chain are bound to be rejected, since they lead to strong excluded volume interactions with other chains.

An attractive alternative are Monte Carlo moves involving at least two chains, which change the connectivity within the melt in a way that the overall packing of beads remains (almost) unaffected. Such algorithms were first introduced for lattice models and recently extended to off-lattice models. The efficient equilibration usually comes at the expense of a certain amount of polydispersity which is introduced into the samples, though recently Karayiannis have overcome this limitation using double-bridging Monte Carlo moves. For the bead-spring polymers studied here, the algorithm is straightforward to implement since the bond length $`l`$ is comparable to the excluded volume parameter $`\sigma `$. This certainly is a unique situation, as for many bead spring coarse grained models of specific chemical species, the above ratio is not as close to one. For example in the case of a united atom model for polyethylene where the bond length $`l=1.54\AA `$ is significantly smaller then $`\sigma 4.0\AA `$, it is necessary to construct two trimer bridges between the two properly chosen dimers along the backbone. For the bead spring model studied here, complex bridging moves are not necessary. For this reason we refer to the method as the double-pivot algorithm since it is much closer to the pivot algorithm used for single chains in an implict, continuum solvent.

In order to maintain a monodisperse melt we search for a pair of spatial neighbor monomers $`i`$ and $`j`$ on different chains or on the same chain, which happen to also be the same distance from one end of their respective chain. Then one can replace two bonds along the original chains with two new bonds or bridges across the pair of chains. The total change in energy $`\mathrm{\Delta }E`$ is local and consists of the sum of the energies of the two new bonds minus the energies of the two previous bonds as well as the difference in the bending energy, which involves the sum of the energies of four new angles minus the energies of the four previous angles for the case $`k_\theta 0`$. The move is accepted by a standard Metropolis criterion, namely with a Boltzmann weight $`\mathrm{exp}(\mathrm{\Delta }E/k_BT)`$ if $`\mathrm{\Delta }E>0`$ and $`1`$ if $`\mathrm{\Delta }E0`$. As in the original pivot algorithm, the double-pivot algorithm generates new chains which are substantially different from the original two chains since as many as $`N/2`$ monomers of a chain are replaced by monomers from the other chain. This results in immediate large changes in the end-to-end distance and radius of gyration.

We have implemented the double-pivot algorithm into both our shared memory and distributed memory MD codes. The details of the implementation and a discussion of its computational efficiency will be published elsewhere. Briefly in our shared memory code, every $`35`$ time steps, we randomly chose a monomer $`i`$ and check its non-bonded nearest neighbors to determine if any satisfy the condition that they are equidistant from the end of their chain. If so, the energy change $`\mathrm{\Delta }E`$ of the double-pivot move is determined and the move accepted or rejected on the basis of the Metropolis criterion. If the move is accepted, monomers and/or the connectivity table are re-labeled depending on which of the two codes is used and the MD simulation is continued. If the move is rejected or no suitable pair is identified, a new monomer is chosen at random and the process repeated. If no suitable pair is generated after searching a specified fraction of monomers (usually $`25\%`$ depending on the chain length $`N`$), we continue with the MD simulations for another $`35`$ steps and repeat the procedure. On the distributed memory code, the procedure is similar except that each processor randomly choses a monomer $`i`$ and checks its non-bonded nearest neighbors on the same processor to see if they satisfy the condition that they are equidistant from the end of their chain. If no suitable non-bonded neighbor is found, another monomer $`i`$ is randomly selected. The process is continued until a specified fraction of the momomers, typically 50%, are tested. To facilitate determining the distance from the free end, each chain carries an extra label starting from $`1`$ at either end to $`N/2`$ at the center. Because each processor searches a unique region of space independently, it is possible to have multiple successful moves each time the procedure is applied. The search is restricted to sets of monomers on the same processor to avoid the need to communicate changes in chain topologies between processors. This restriction also prevents two or more processors from performing simultaneous swaps that could energetically conflict with each other. While this restriction means (slightly) less swaps are considered, this is more than outweighed by the parallelism, e.g. up to $`P`$ swaps take place in one time step, where $`P`$ is the number of processors. With the standard parameters of the FENE potential, $`k=30ϵ/\sigma ^2`$, the acceptance rate is quite low. A way to improve the acceptance rate is to start with a somewhat lower value, $`k=1015ϵ/\sigma ^2`$ and gradually increase $`k`$ during the course of the simulation. This change has little effect on the overall structure of the chains but increases significantly the acceptance rate. Using our distributed memory code on 27 DEC alpha processors for a system of $`M=500`$ chains of length $`N=500`$, the number of successful exchanges for $`10^5\mathrm{\Delta }t`$ with $`k_\theta =0`$ was approximately $`1750`$ for $`k=30ϵ/\sigma ^2`$ compared to $`8500`$ for $`k=20ϵ/\sigma ^2`$ and $`49500`$ for $`k=10ϵ/\sigma ^2`$.

To demonstrate the effectiveness of the double-pivot algorithm, we have applied it to a melt of $`M=500`$ chains of length $`N=500`$ generated using the standard fast push-off procedure outlined in Sec. V. Figure 8 shows the internal distances as a function of time during a double-pivot/MD relaxation. Note that the characteristic “bump” in the internal distance function does not simply decay. Rather, Figure 8 shows that the (seemingly equilibrated) largest scales are first swollen, before the chain dimensions decay homogeneously on all scales. Equilibration for the reduced spring constant $`k=10ϵ/\sigma ^2`$ was achieved after about $`6\times 10^4\tau `$ or 10 accepted pivot moves per monomer. The penalty for reducing the spring constant to $`k=10`$ is an increase of the mean-square bond length by $`7\%`$, while the chain stiffness $`c_{\mathrm{}}`$ remains unaffected. In subsequent runs, the spring constant $`k`$ is increased to slowly “reel in” in the chains. This required an additional $`6\times 10^4\tau `$. We made no attempt to optimize the number of steps for each value of $`k`$. Nevertheless, the required equilibration times are considerably shorter than the time $`\tau _{max}(500)4\times 10^6\tau `$ for the chains to move their own size in a pure MD simulation.

These observations illustrate the close analogy between the pivot algorithm for single chains and the double-pivot algorithm for dense polymer melts. As pointed out by Madras and Sokal, it is useful to distinguish between decorrelation and equilibration in discussing the performance of such algorithms. Pivoting is extremely efficient in decorrelating large length scales, provided the initial conformation is properly equilibrated. Up to small corrections, the number of pivot moves needed to decorrelate the end-to-end distance of isolated chains in an implicit solvent is $`𝒪(10)`$ independent of chain length. Since the same holds for subchains of arbitrary length, correspondingly more moves are required to decorrelate higher order (Rouse) modes. In other words, the decorrelation proceeds from large to small scales. However, the same is not correct for the equilibration of a perturbed initial configuration. In this case it is necessary to run the system up to the decorrelation time of the shortest perturbed length scale in order to equilibrate the chains on all length scales. For single chain studies aiming at global properties Madras and Sokal  therefore recommend to apply the pivot algorithm to equilibrated initial configurations generated by other techniques. In the present case, the combination of pre-packing with a slow push-off at least partially fulfills these requirements, since the chains are locally equilibrated.

Since equilibration requires that each monomer comes close enough to suitable exchange partners, one may ask whether of $`𝒪(10)`$ exchanges per monomer independent of $`N`$ is sufficient or whether the exchanges are simply being made back and forth between a limited number of states. To estimate whether enough independent configurations are being sampled, consider the time for a monomer to diffusive a typical distance between monomers with the same index $`(\rho /(N/2))^{1/3}`$ between successful pivot moves. If we assume that monomers diffuse with the typical $`\delta r^2=\sqrt{k_BTb^2t/\zeta }`$ Rouse dynamics, this corresponds to

$$t_{piv}=\left(\frac{(N/2)}{\rho }\right)^{4/3}\frac{\zeta }{k_BTb^2}$$
(9)
which is equivalent to the Rouse time of a chain of length $`N_{eff}=(N/2)^{2/3}`$. For our longest chains $`N=7000`$, $`N_{eff}236`$. The Rouse time for a chain of length $`N=236`$ is approximately $`7\times 10^4\tau `$, which is still much less than the time needed to make $`𝒪(10)`$ exchanges per monomer even if we used the reduced spring constant $`k=10`$. Thus only for much longer chains or more complex architectures does one have to worry about sufficient independent configurations being sampled by this hybrid double-pivot-MD algorithm.

## VIII Local bending rigidity

As a last point we present results for chains with intrinsic stiffness Eq. (2). In section IV we referred to this case as a “blind test,” because the starting states were generated based on our simple estimate Eq. (7) for the effective chain stiffness in the melt and before the brute-force MD runs for the target functions shown in Fig. 2 were completed. As can be seen in Fig. 1, our original estimates were slightly too large. A comparison with the target functions after the pre-packing and slow push-off phase (Fig. 9a) shows that as for fully flexible chains the correct chain statistics is reproduced on short length scales. However, we are now in a situation comparable to the situtaion presented in Fig. 7 where the large length scales are not fully equilibrated, because the initial chains were generated with a slightly incorrect effective stiffness.

Subsequently, we ran the combined double-pivot-MD simulation for $`1.08\times 10^5\tau `$ with $`k=10`$ for the first 4 millions steps, $`k=20`$ for the second $`4`$ million steps and $`k=30`$ for the last million steps for a system of $`200`$ chains of $`N=350`$ for $`k_\theta =1.5`$ and $`2.0`$ and $`N=700`$ for $`k_\theta =0.75`$. The resulting mean squared end-to-end distance were in excellent agreement with the target functions generated by brute force MD simulations for shorter chains. Table II shows an example of the number of successful exchanges per $`10^5`$ steps for the system of $`200`$ chains of length $`N=350`$ monomers for three values of the spring constant $`k`$ for four values of $`k_\theta `$ for our distributed memory code run on 27 DEC alpha processors. As seen from this table, the procedure used gives approximately $`912`$ successful exchanges for $`k=10`$, $`2`$ successful exhanges for $`k=10`$ and $`0.1`$ for $`k=30`$. Fortunately there is little difference in the local structure of the chain for $`k=20`$ compared to $`k=30`$, though as seen from Fig. 8, the chain is expanded for $`k=10`$. Thus there is a definite tradeoff between reducing the bond spring constant $`k`$, which increases the number of exchanges while at the same time increasing the mean squared bond length. Note that for longer chains lengths, the length of the run and/or the number of processors would have to be increased so that the number of successful exchanges is on the order of $`𝒪(10)`$. This limits the double-pivot/MD hydrid method to moderate chain lengths.

## IX Conclusions

In this paper, we have discussed the preparation of equilibrated melts of long chain polymers in computer simulations. The interest in reliable algorithms for this task is due to two problems: (i) the prohibitively long relaxation times encountered in brute force MD simulations and (ii) the difficulty to predict the chain statistics on large scales (i.e. the effective chain stiffness $`c_{\mathrm{}}`$) on the basis of microscopic intra- and interchain interactions.

Our results can be summarized as follows:

* It is insufficient to judge the quality of the equilibration from the statistics of the chain end-to-end distances or radius of gyration. The first just measures one length, while the latter gives a non specific average over all internal distances. As a much more sensitive criterion we measure mean square internal distances on all length scales from the monomer size up to the contour length of the chains.
* Our standard method to rapidly introduce excluded volume interactions between randomly assembled Gaussian chains with the correct overall statistics (“fast push-off”) fails, when judged according to this improved criterion. The fast push-off introduces significant chain distortions on length scales of the order of the tube diameter. They can only vanish when the local melt topology is equilibrated. In MD simulations this is only possible via reptation dynamics so that the proper equilibration requires runs whose length exceeds the disentanglement time of the chains. To overcome this problem we have introduced two different methods.
* We have modified our standard approach by reducing the density fluctuations in the assembly of Gaussian chains (“pre-packing”) and by introducing the excluded volume interactions in a quasi-static manner (“slow push-off”). As before, this method requires prior knowledge of $`c_{\mathrm{}}`$. Tests for bead-spring models with chain lengths up to $`N=7000`$ $`(N/N_e=𝒪(100))`$ show that suitable push-off times are the Rouse times of the characteristic chain lengths of the overshoot. In the present case this time is of the order of a few entanglement times, independent of $`N`$.
* We have applied a “double-pivot” algorithm which is inspired by the highly efficient pivot algorithm for single chains and the double-bridging method for dense systems. For intermediate chain lengths, this is a viable method to equilibrate melts in a reasonable amount of cpu time, particularly when one has no prior knowledge of $`c_{\mathrm{}}`$.

The combination of MD relaxation, double-pivot and slow push-off allows the efficient and controlled preparation of equilibrated melts of short, medium and long chains respectively. While our results were obtained for an off-lattice bead spring model with chain lengths up to $`N=7000`$ beads, the methods should work as efficiently for lattice polymer models such as the bond fluctuation model. The pre-packing and gradual introduction of the excluded volume are also applicable to united and explicit atom simulations. Furthermore, the methods are applicable to polydisperse melts and to branched polymers, including long chain branching, comb and star polymers. The double-pivot and double bridging algorithms has also been used to equilibrate long end-tethered, polymer brushes in contact with a melt of long chains. For these more complex chain architectures, in which the chains are not necessarily Gaussian and there is no a priori way to know the conformation of the chains, algorithms like the double-pivot and double-bridging algorithms are presently the only possible way to equilibrate systems faster than extremely long MD or MC runs.


text:
# References

## Abstract

In cosmological backgrounds, there can be ’partially massless’ higher spin fields which have fewer degrees of freedom than their massive partners. The equations for the partially massless spin-2 fields are usually taken to be the linearized Einstein equations augmented with a ’tuned’ Pauli-Fierz mass. Here, we add more powers of curvatures and show that for the string-generated Einstein-Gauss-Bonnet model, partially massless spin-2 fields have real mass in AdS, in contrast to the Einstein level result. We discuss the implication of this for the AdS/CFT applications and briefly study the $`C^4`$-corrected $`AdS_5\times S^5`$ solution in type IIB SUGRA.

BRX-TH-508

Partially Massless Spin-2 Fields in String Generated Models
Bayram Tekin 111e-mail: tekin@brandeis.edu
Department of Physics, Brandeis University, Waltham, MA 02454, USA

In Anti-de-Sitter and de-Sitter spacetimes, massive higher spin fields ( $`s2`$ ) exhibit a novel phenomenon: partial masslessness (or the appearance of extra gauge invariances), if the masses of the fields are ‘tuned’ to the background cosmological constant . A massless spin-2 field (which we shall exclusively employ with here) in a flat $`D`$-dimensional background has $`D(D3)/2`$ degrees of freedom ( DOF) and a massive one has $`(D+1)(D2)/2`$ DOF. On the other hand, in a constant curvature background there is a third possibility, for which the spin-2 field has one fewer than the generic massive case: Hence the name partially massless. A second derivative, scalar gauge invariance appears at the partially massless point. For these fields, there is a crucial difference between dS and the AdS spaces: The tuned $`m^2`$ is positive for the former and negative for the latter.

In the literature, the original analysis involves the equations of the linearized gravity augmented with the linear Pauli-Fierz mass term. In this note, we shall slightly generalize by considering also (the linearizations of ) higher curvature terms that are generated in certain string models and study their effects on the partially massless fields.

The motivation for this work comes from the possible applications of partially massless spin-2 fields to AdS/CFT duality . As is well-known, bulk fields in AdS correspond to operators in the boundary conformal field theory. According to the AdS/CFT dictionary, the masses of the bulk fields determine the conformal dimensions of the boundary fields. In particular, it was shown by Dolan et. al. , that partially massless spin-2 fields $`(h_{\mu \nu })`$ in the bulk of AdS space correspond, on the boundary, to fields $`(L_{ij})`$, which obey a specific conformally-invariant differential equation. If the boundary of the AdS is flat, the differential operator reduces to partial conservation equation of a tensor field (operator) : $`_i_jL^{ij}=0`$. As demonstrated in , due to $`m^2`$ being negative for partially massless fields in AdS, the first descendant of the corresponding partially conserved operator has a negative norm, rendering the boundary theory unphysical. Therefore, at the linear level there seems to be no room for the partially massless fields in AdS. For dS, on the other hand, we do not have a proper dS/CFT correspondence to make use of these fields. What happens beyond the (already studied ) quadratic Einstein level, in AdS, is our main aim in this work.

Here we suggest that, in certain cases ( though, not in the most famous $`AdS_5\times S^5`$ example), stringy higher curvature corrections flip the sign of $`m^2`$ in AdS. In particular, we will show that if the equations of motion for the massive spin-2 field in an AdS background are taken to be the linear version of the Einstein plus Gauss-Bonnet theory, the partially massless fields have positive mass.

We also study the $`(\text{Weyl})^4`$ corrected spin-2 action in the maximally supersymmetric $`AdS_5\times S^5`$ case, and show that partially massless fields are ruled out here, supporting the conclusion of . \[ Strictly speaking, by imposing the Pauli-Fierz term, which does not arise in string theory, we are changing the background geometry. \]

We first review the main ingredients of the partially massless spin-2 fields in a cosmological background, and then consider the spin-2 field equations derived from higher curvature theories. Finally, we study the $`AdS_5\times S^5`$ background as a solution to Type IIB String/SUGRA equations with $`R^4`$ quantum corrections.

Our conventions are: signature $`(,+,+,\mathrm{}+)`$, $`[_\mu ,_\nu ]V_\lambda =R_{\mu \nu \lambda }^\sigma V_\sigma `$, $`R_{\mu \nu }R_{\mu \lambda \nu }^\lambda `$. We derive the equations of motions of a spin 2-field by linearizing the gravity part of the low energy string-generated gravity equations around the relevant background. Let us start with the lowest order Einstein equations

$`R_{\mu \nu }{\displaystyle \frac{1}{2}}g_{\mu \nu }R+\mathrm{\Lambda }g_{\mu \nu }=0,`$ (1)
linearized around $`\overline{g}_{\mu \nu }`$, for which the Ricci tensor is $`\overline{R}_{\mu \nu }=\frac{2}{D2}\mathrm{\Lambda }\overline{g}_{\mu \nu }`$. We obtain the following covariantly conserved equation for the massless spin-2 field $`h_{\mu \nu }g_{\mu \nu }\overline{g}_{\mu \nu }`$

$`𝒢_{\mu \nu }^LR_{\mu \nu }^L{\displaystyle \frac{1}{2}}\overline{g}_{\mu \nu }R^L{\displaystyle \frac{2}{D2}}\mathrm{\Lambda }h_{\mu \nu }=0,`$ (2)
where $`R^L=(g^{\mu \nu }R_{\mu \nu })_L`$ and the linear part of the Ricci tensor is defined as

$`R_{\mu \nu }^LR_{\mu \nu }\overline{R}_{\mu \nu }={\displaystyle \frac{1}{2}}(\overline{\mathrm{}}h_{\mu \nu }\overline{}_\mu \overline{}_\nu h+\overline{}^\sigma \overline{}_\nu h_{\sigma \mu }+\overline{}^\sigma \overline{}_\mu h_{\sigma \nu }).`$ (3)
All derivatives are with respect to the background metric, which also raises and lowers the indices. Here: $`h=\overline{g}^{\mu \nu }h_{\mu \nu }`$ and $`\overline{\mathrm{}}=\overline{g}^{\mu \nu }\overline{}_\mu \overline{}_\nu `$. To get a massive spin-2 field, we add the usual ghost-free Pauli-Fierz mass, with the correct relative sign:

$`𝒢_{\mu \nu }^L+{\displaystyle \frac{1}{2}}m^2(h_{\mu \nu }\overline{g}_{\mu \nu }h)=0`$ (4)
For $`m^20`$, taking the double divergence and the trace of this equation, one obtains, respectively:

$`\overline{}_\mu \overline{}_\nu h^{\mu \nu }\mathrm{}h=0,`$ (5)
$`\left[\mathrm{\Lambda }{\displaystyle \frac{(D1)}{2}}m^2\right]h=0.`$ (6)
Unlike the generic $`(m^2,\mathrm{\Lambda })`$ theory, $`h_{\mu \nu }`$ does not have to be traceless (i.e $`h0`$) at the partially massless point for which the mass is tuned as

$`m^2={\displaystyle \frac{2}{D1}}\mathrm{\Lambda }.`$ (7)
One condition on the field drops out, yielding instead the following (higher derivative) scalar gauge invariance

$`\delta h_{\mu \nu }=\overline{}_\mu \overline{}_\nu \xi +{\displaystyle \frac{2}{(D1)(D2)}}\overline{g}_{\mu \nu }\mathrm{\Lambda }\xi `$ (8)
Needless to say (7) assigns a negative $`m^2`$ for the partially massless fields in AdS. This by itself does not pose a threat since AdS allows negative $`m^2`$ for various fields, including massive spin-2 field as long as the corresponding Brietenlohner-Freedman  bounds are satisfied.

In , whose notations and results we follow, partially massless spin-2 field was employed in the context of AdS/CFT. The authors showed that such a field in the bulk of $`AdS_D`$ corresponds to a symmetric tensor $`L_{ij}`$ on the boundary, which satisfies the following conformally invariant differential equation

$`\left[\overline{}_i\overline{}_j+{\displaystyle \frac{1}{D3}}\stackrel{~}{R}_{ij}\right]L^{ij}=0,`$ (9)
where $`\stackrel{~}{g}_{ij}`$ is the conformal metric on the boundary. This differential equation was studied before by Eastwood and collaborators . As in the simplest examples of AdS/CFT, let us assume that the boundary is conformally flat, then (9) reduces to the more transparent partial conservation law

$`_i_jL^{ij}=0.`$ (10)
Suppose now $`|L^{ij}>`$ is a symmetric primary state on which the $`D1`$ dimensional conformal group acts. Let $`P_i`$ be the momenta; as a a result of (10), we require $`P_iP_j|L^{ij}>`$ to be a null state. This then implies the existence of a negative norm state: $`||P_iL^{ij}>|^2<0`$. In the language of the bulk fields, this negative norm state corresponds to $`m^2`$ being negative for the partially massless AdS fields (7).

At the linear level, the above result signals an impasse for the partially massless fields in the context of AdS/CFT. But here, I will argue that higher curvature corrections can drastically change the linear level results in certain theories. In fact, though rare, it is not unheard-of that higher derivative terms alter the sign of the linear level theory. \[ To warm up, let us recall the non-trivial example of the Topologically Massive Gravity (TMG) in $`2+1`$ dimensions . This model is constructed by adding a third derivative Cotton tensor (or Chern-Simons term ) to the Einstein one $`R_{\mu \nu }\frac{1}{2}g_{\mu \nu }R+\frac{1}{M}C_{\mu \nu }=0.`$ Unlike pure 3D Einstein gravity, TMG is a dynamical theory: Linearized equations show that there is a massive scalar degree of freedom with mass $`|M|`$. For this scalar mode to be non-ghost with positive Hamiltonian, the sign of the Einstein action has to be flipped from its usual one once the higher derivative topological term is introduced . Even though TMG is quite distinct from the theories we are about to consider, it nevertheless provides us with a concrete example of how higher derivative terms can change linear level signs.\]

Let us now turn to Einstein-Gauss-Bonnet model , which is expected to arise in some string theories \[ or string theory compactifications\].

$`I={\displaystyle d^Dx\sqrt{g}\left\{\frac{R}{\kappa }+\gamma (R_{\mu \nu \rho \sigma }^24R_{\mu \nu }^2+R^2)\right\}}.`$ (11)
String theory dictates the sign of $`\gamma `$ to be positive. Even without an explicit cosmological constant, AdS \[not dS!\] is a solution  to the equations of motion derived from (11)

$`0={\displaystyle \frac{1}{\kappa }}(R_{\mu \nu }{\displaystyle \frac{1}{2}}g_{\mu \nu }R)`$
$`+2\gamma \left\{RR_{\mu \nu }2R_{\mu \sigma \nu \rho }R^{\sigma \rho }+R_{\mu \sigma \rho \tau }R_\nu ^{\sigma \rho \tau }2R_{\mu \sigma }R_\nu ^\sigma {\displaystyle \frac{1}{4}}g_{\mu \nu }(R_{\tau \lambda \rho \sigma }^24R_{\sigma \rho }^2+R^2)\right\}.`$ (12)
The cosmological constant of the AdS space that solves the equations is determined by ($`\gamma ,\kappa `$)

$`\mathrm{\Lambda }={\displaystyle \frac{2(D1)(D2)}{\kappa \gamma (D3)(D4)}}.`$ (13)
Linearizing (12) about AdS, one finds  how the presence of $`\gamma `$ modifies the massive spin-2 equation,

$`𝒢_{\mu \nu }^L\left[{\displaystyle \frac{1}{\kappa }}+{\displaystyle \frac{4\mathrm{\Lambda }\gamma (D4)(D3)}{(D2)(D1)}}\right]+{\displaystyle \frac{1}{2\kappa }}m^2(h_{\mu \nu }\overline{g}_{\mu \nu }h)=0`$ (14)
subject to the condition (13). Thus one actually obtains the sign-flipped version of (4)

$`𝒢_{\mu \nu }^L{\displaystyle \frac{1}{2}}m^2(h_{\mu \nu }\overline{g}_{\mu \nu }h)=0`$ (15)
The crucial and non-trivial point here is that both (4) and (14) come from proper, non-ghost Lagrangians, despite ‘contradicting’ each other. It is , by now, clear that the partially massless states in AdS will be non-tachyonic for Einstein-Gauss-Bonnet theory. Namely for (15): $`m^2=2\mathrm{\Lambda }/(D1)>0`$. It also follows that if this model describes spin-2 fields in the bulk of $`AdS_D`$, for $`D5`$, partially massless fields correspond to conformally invariant boundary operators with positive descendants. Thus adding the Gauss-Bonnet term would circumvent the objection of  ( at least beyond four dimensions ).

On the other hand, if we include an explicit cosmological constant ( $`\mathrm{\Lambda }_0`$, whose sign is arbitrary ), the above analysis branches into several directions. For given $`(\mathrm{\Lambda }_0,\kappa ,\gamma )`$ there are two different AdS spacetimes which solve the equations. \[ See the related work \]. The cosmological constant of these spaces are

$`\mathrm{\Lambda }_\pm ={\displaystyle \frac{1\pm \sqrt{1+8\kappa \stackrel{~}{\gamma }\mathrm{\Lambda }_0}}{4\kappa \stackrel{~}{\gamma }}}`$ (16)
where $`\stackrel{~}{\gamma }=\gamma (D4)(D3)/[(D2)(D1)]`$ and $`8\kappa \mathrm{\Lambda }_0\stackrel{~}{\gamma }1`$. The equation (15) becomes

$`𝒢_{\mu \nu }^L\pm {\displaystyle \frac{1}{2}}{\displaystyle \frac{m^2}{\sqrt{1+8\kappa \mathrm{\Lambda }_0\stackrel{~}{\gamma }}}}(h_{\mu \nu }\overline{g}_{\mu \nu }h)=0`$ (17)
and the partially massless fields have the mass

$`m^2=\pm {\displaystyle \frac{2\mathrm{\Lambda }_\pm }{D1}}\sqrt{1+8\kappa \mathrm{\Lambda }_0\stackrel{~}{\gamma }}`$ (18)
One, now, has various choices, depending on the sign of $`\mathrm{\Lambda }_0`$ and on the choice of $`\mathrm{\Lambda }_+`$ or $`\mathrm{\Lambda }_{}`$. For $`\mathrm{\Lambda }_0<0`$, we have $`\mathrm{\Lambda }_\pm <0`$, but only $`\mathrm{\Lambda }_{}`$ gives real mass to partially massless fields.

Having pointed out that certain higher curvature spin-2 theories have non-tachyonic masses for the partially massless fields, can we now make use of this result in AdS/CFT? Unfortunately, I do not know any examples of latter where there are Gauss-Bonnet corrections to bulk gravity. One would expect that less supersymmetric forms of AdS/CFT duality might allow Gauss-Bonnet terms \[See  and the references therein.\] But, for now, let us briefly study the well-understood maximally supersymmetric case. In low energy Type IIB SUGRA the next to leading order corrections (in the $`\alpha ^{}`$ expansion ) around the $`AdS_5\times S^5`$ vacuum are of the $`C^4`$ type( here $`C^4`$ schematically represents the two inequivalent contractions of four Weyl tensors $`C_{\mu \nu \alpha \beta }`$ . Switching off the dilaton and all the fields but spin-2, the action reads

$`I={\displaystyle \frac{1}{\alpha _{}^{}{}_{}{}^{4}}}{\displaystyle d^{10}x\sqrt{g}\left(R+k\alpha _{}^{}{}_{}{}^{3}f_4(\rho ,\overline{\rho })C^4\right)}`$ (19)
where $`\rho `$ is the usual complex coupling and $`f_4`$ is given by an Eisenstein series, whose details will not be relevant to us, but can be found in . $`k`$ is a constant. $`AdS_5\times S^5`$ is a background solution. ( Due to the flux from the self-dual five form $`F_5`$, a cosmological constant is generated.) $`AdS_5\times S^5`$, with the same radius on $`AdS_5`$ and $`S^5`$ is conformally flat. This vacuum is expected to be a solution to full string theory:Thus the appearance of $`C^4`$, which vanishes for the conformally flat geometries, is no surprise. To see if the higher order terms can affect the partially massless fields, one can linearize the equations derived from (19) around $`AdS_5\times S^5`$ and add a Pauli-Fierz mass term. But it is immediately clear that there will be no contributions from the $`C^4`$ terms since, at $`O(h^2)`$ level one has $`\overline{C}^2\delta (\sqrt{g}C^2)`$, which vanishes for all conformally flat spaces. Thus $`C^4`$ corrections cannot change the structure of the partially massless fields. With negative $`m^2`$, since the latter predict the existence of negative norm states in the boundary $`𝒩=4`$ $`SU(N)`$ Yang-Mills theory, they ought to be ruled out in the bulk. It is important to note that if the corrections were of the generic $`R^n`$-type, where $`R`$ represents scalar,Ricci, or Riemann curvatures, the equations for the spin-2 field and its partially massless point would have changed drastically. On the other hand $`C^3`$, and $`C^2`$ corrections do not change the picture nor do arbitrary power of $`C`$. Up to $`O(h^2)`$, $`C^3`$ terms vanish just like the $`C^4`$ ones for AdS. Even though $`C^2`$ does not vanish at this order, it does not contribute to the trace of the equations, and thus does not affect the partially massless point.

To summarize, we have shown that, unlike the linear spin-2 theory, partially massless fields have positive $`m^2`$ in AdS, if Gauss-Bonnet corrections are added for $`D5`$. In the context of AdS/CFT, these partially massless fields in the bulk correspond to operators with positive norm descendants in the boundary.

This paper is dedicated to the memory of Ian Kogan (1958-2003). I would like to thank S. Deser and A. Waldron for useful discussions. I also would like to thank A. Lawrence, H. Schnitzer and N. Wyllard for directing me to the references regarding the $`R^4`$ corrections in string theory. This research was supported by NSF Grant 99-73935.