text:
# Detection of OH Towards the Extreme Carbon Star IRC+10216

## 1 Introduction

IRC+10216 is the prototypical late-stage, carbon-rich asymptotic giant branch (AGB) star. This star is extremely well-studied, due to its close proximity and large mass loss rate; the star is losing mass at a rate of $`3\times 10^5`$M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$ (Glassgold 1996), producing a dense, dusty circumstellar envelope, well shielded from interstellar ultraviolet (ISUV) radiation. The dense, shielded environment of IRC+10216’s envelope is home to a rich circumstellar chemistry. To date, more than 50 molecules have been detected around IRC+10216. Due to the extreme carbon-rich nature of the source (C/O $`1.4`$) (Glassgold 1996), the detection of any oxygen-bearing molecules other than CO and small amounts of SiO and HCO+ was entirely unexpected. Recently, however, Melnick et al. (2001) reported the detection of an emission feature at $`556.936`$GHz, attributed to the $`1_{10}1_{01}`$ transition of water vapor. The detection of water vapor in such a carbon-rich circumstellar envelope was interpreted by Melnick et al. (2001) as evidence for the existence of an extrasolar cometary system, analogous to the Solar System’s Kuiper Belt, in orbit around IRC+10216. In this system, the luminosity of the central star has increased dramatically due to the later stages of post-main sequence stellar evolution, causing the icy bodies in the Kuiper Belt analog to vaporize (Ford & Neufeld 2001) and produce the water vapor observed by Melnick et al. (2001).

Of possible concern for the cometary hypothesis is the existence of a peculiar subclass of carbon stars which exhibit OH and H2O maser emission, as well as silicate emission features (e.g. Engels 1994). These features would normally be unexpected in carbon stars, for the reasons explained above. This subclass also shows a number of other features which distinguish peculiar carbon stars from normal carbon stars, including relatively high 13C/12C ratios, IRAS colors compatible with low mass loss rate oxygen stars, and depleted s-process element abundances relative to other carbon stars. These peculiar carbon stars may represent objects in the transition stage between oxygen star and carbon star, and their photospheric chemistry is likely to be out of thermochemical equilibrium, or they may be members of binary pairs of one oxygen star and one carbon star (Wallerstein & Knapp 1998). In any case, IRC+10216 displays none of the ancillary characteristics associated with these peculiar carbon stars, aside from OH and H2O emission. Due to the large C/O ratio in IRC+10216 and the fact that it is not a member of a binary system, as well as the lack of ancillary characteristics, we exclude the possibility that OH or H2O emission from IRC+10216 arises from membership in this peculiar class of carbon stars.

Since the Melnick et al. (2001) result was based on the detection of a single line of water vapor, we became interested in finding a method of verifying the identification of water vapor with the emission feature at $`556.936`$GHz. To this end, we realized that the water vapor emitted by icy bodies in the cometary system would be carried by the circumstellar outflow and shielded by the dust, until the outflow became diffuse enough to allow penetration by the ISUV field. Once the ISUV field could penetrate the outflow, the water vapor would be photodissociated into OH and atomic hydrogen. The OH could be detected, and such a detection would yield support for the Melnick et al. (2001) identification. Based on our model for the production of OH, we proposed a deep search for OH around IRC+10216 (a previous search by Andersson, Wannier & Lewis (1992) resulted in a non-detection); we were awarded $`30`$ nights at the Arecibo Observatory111The Arecibo Observatory is part of the National Astronomy and Ionosphere Center, which is operated by Cornell University under a cooperative agreement with the National Science Foundation. for our search. We present our observations in §2, discuss the implications of our results in §3 and present our conclusions in §4.

## 2 Observations

The radio telescope at Arecibo Observatory (AO) is an altitude-azimuth telescope with a fixed $`304.8`$m primary spherical reflector and a movable Gregorian dome which contains secondary and tertiary reflectors, as well as multiple receivers. For the observations discussed in this paper, the L-wide receiver was used with dual linear polarizations with each of four boards of the correlator in a 9-level mode, centered on the frequency of interest and covering a bandwidth of $`3.125`$MHz divided into $`1024`$channels, resulting in an unsmoothed velocity resolution of $`0.55\mathrm{km}\mathrm{s}^1`$ at $`1665`$ and $`1667`$MHz. Observations were taken 2001 November 27 – December 10; 2002 February 2 – 3, 12, 14, 23; 2002 March 2 – 4, 24, 27 – 31; and 2002 April 1, which resulted in a total on-source integration time of $`1613`$minutes. At least one flux calibration was done during each observing session, on one, two or three of the following sources: B0838+133, J0935+086, B1117+146, B1040+123 or B1005+077, depending upon the amount of time available. Our calibration sources were observed at a variety of zenith angles, so we may take an average gain from these measurements as representative of the telescope gain for our target observations. We found the telescope gain to be $`8.2`$K/Jy at $`1666`$MHz. All observations were conducted using standard position switching, with equal duration $`5`$minute on-source and off-source scans, except on 2001 November 28, when the scan duration varied between $`3`$ and $`5`$minutes in an effort to maximize observing time. The transitions of primary interest, the $`1665`$ and $`1667`$MHz main lines of OH, were observed during all sessions. Other transitions which were only observed during some of the sessions were the $`1612`$ and $`1721`$MHz satellite lines of OH, the $`1420`$MHz hyperfine line of HI, the $`1347`$MHz $`2_12_1`$ line of $`\nu _2=1`$ HCN and the $`1743`$MHz $`32`$ line of HC9N.

Data reduction was accomplished by converting the raw AO data to Continuum and Line Analysis Single-dish Software (CLASS) format. To prevent contamination of the data by RFI signals, each scan was examined by eye, and any scans with apparent RFI (i.e. those with prominent, rapidly varying signals) were eliminated. The remaining scans were summed to produce a total spectrum for each observed frequency. The total spectrum was then smoothed, using the CLASS “smooth” macro, resulting in a smoothed velocity resolution of $`1.1\mathrm{km}\mathrm{s}^1`$. A first order polynomial baseline was subtracted from the smoothed data. In the case of the OH main lines, an emission feature was readily apparent in the final data, and a Gaussian profile was fit to each feature. The calibrated main line spectra, together with the Gaussian fits, are displayed in Figure 1. The fit parameters for these lines are listed in Table 1. The errors listed in the table are $`1\sigma `$; thus the $`1665`$ and $`1667`$MHz lines are detected at a significance of $`4.7\sigma `$ and $`8.8\sigma `$, respectively. We determine an average fractional abundance for OH of $`x(\mathrm{OH})4\times 10^8`$ using the $`1667`$MHz line flux, $`F_{1667}=39\mathrm{mK}\mathrm{km}\mathrm{s}^1=4.8\mathrm{mJy}\mathrm{km}\mathrm{s}^1`$, and the relationship between $`F_{1667}`$ and the number of OH molecules, $`N(\mathrm{OH})=4\pi D^2F_{1667}/f_uA_{ul}h\nu `$, where $`D=170`$pc is the distance to IRC+10216, $`f_u=5/16`$ is the fraction of molecules in the upper state (the value of 5/16 is plausible for the most likely excitation conditions), $`A_{ul}=7.8\times 10^{11}\mathrm{s}^1`$ is the Einstein A coefficient for the transition and $`h\nu `$ is the energy of the transition. Once we have found $`N(\mathrm{OH})`$, we may divide by the number of H nuclei in the Arecibo beam to find $`x(\mathrm{OH})`$. Based on a density of H nuclei of $`n_H=3\times 10^{37}(\dot{M}_5/v_{exp,6})/r^2`$cm-3 (Glassgold, 1996), where $`\dot{M}_5=3`$ is the stellar mass loss rate in units of $`10^5\mathrm{M}_{}\mathrm{yr}^1`$, and $`v_{exp,6}=1.4`$ is the expansion velocity222The expansion velocity determined by Melnick et al. (2001) from the water vapor line in IRC+10216 is $`v_{exp,6}=1.7`$. We use $`v_{exp,6}=1.4`$ instead, since this is the value found previously for IRC+10216 (e.g. Skinner et al. 1999), and because there are large errors associated with the Melnick et al. measurement. of the circumstellar envelope in units of $`10^6\mathrm{cm}\mathrm{s}^1`$, we find $`7.8\times 10^{55}`$ hydrogen nuclei, assuming a maximum radius of $`r_{max}=10^{17}`$cm. We choose $`r_{max}=10^{17}`$cm because we expect nearly all of the OH to lie inside of this radius, and we do not correct for beam sensitivity since this is within the half-power beamwidth of AO and any such correction would be at the level of a few percent at most. This yields a rough determination of $`x(\mathrm{OH})`$ which is relatively model independent but somewhat misleading. Since we believe the OH to be in a photodissociated shell or ring around IRC+10216, the OH fractional abundance in the outer regions of the circumstellar envelope around IRC+10216 may reach values several times larger than our minimum.

We failed to detect any signal originating in the envelope of IRC+10216 at the five other frequencies we observed. Severe radio frequency interference (RFI) across all observations prevents us from placing meaningful upper limits on the strength of the $`1721`$MHz satellite line of OH or the $`1743`$MHz line of HC9N. Strong contamination from interstellar gas at the velocity of IRC+10216 prevents us from making any strong statements about stellar emission at $`1420`$MHz, the hyperfine line of HI. Our results suggest the necessity of interferometric observations to detect and map the contribution to HI emission from the star. For the remaining lines, we were able to determine upper limits on the observed flux. These are also listed in Table 1.

## 3 Discussion

In conjunction with the observations discussed above, we constructed several models which attempted to predict the strength and shape of the OH lines in IRC+10216. First, we used a photodissociation model for H2O and OH, assuming a dusty spherical outflow, illuminated from the exterior by the interstellar ultraviolet (ISUV) field. We used photodissociation rates and dust opacities based on Roberge et al. (1991). The rates in Roberge et al. (1991) were calculated for a plane-parallel slab, so a correction for the spherical geometry of the circumstellar envelope around IRC+10216 was applied. This correction is described in detail in Appendix A. Our model assumed spherical shielding of the molecules due to the spherical distribution of dust around IRC+10216, but did not assume a spherical molecular distribution, since molecular self-shielding was neglected, as is appropriate under the condition of small molecular abundances. A given model produces the relative abundance of H2O and OH as a function of radial distance. We generated a number of abundance distributions by varying the unshielded molecular photodissociation rate, the stellar mass loss rate, $`\dot{M}`$, and the dust opacity. The unshielded molecular photodissociation rate yields the strength of the ISUV field, the stellar mass loss rate determines the density of the shielding dust, and the dust opacity indicates the shielding strength of the dust. Photodissociation models were parameterized by assuming a photodissociation rate, $`\mathrm{\Gamma }_i`$, of the form $`\mathrm{\Gamma }_i(A_v)=2C_iG_0\mathrm{exp}(\alpha _iA_v\beta _iA_v^2)\mathrm{s}^1`$, where $`i`$ represents either OH or H2O, $`A_v`$ is the depth in visual magnitudes and $`C_i`$, $`\alpha _i`$ and $`\beta _i`$ are factors derived from Roberge et al. (1991) using the procedure described in Appendix A. The factor $`G_0`$ represents the relative strength of the ISUV field and $`G_0=1`$ defines the average ISUV field (Draine 1978). We also considered variations in the ratio of $`A_v`$ to $`N_H`$, where $`N_H`$ is the column density of hydrogen nuclei. This ratio, $`A_v/N_H`$, represents how much extinction, and therefore molecular shielding, is produced per unit column density of hydrogen nuclei. For standard interstellar grains considered by Roberge et al. (1991), $`A_v/N_H=(A_v/N_H)_05.6\times 10^{22}\mathrm{mag}\mathrm{cm}^2`$. The photodissociation model was then used as an input to a second program, which calculated the expected strength and shape of the OH emission lines, including instrumental effects (ie: convolution of beam sensitivity with target intensity) and Gaussian microturbulence. For the beam convolution, we assumed the antenna response is a Gaussian, with $`\theta _e=104^{\prime \prime }`$ where $`\theta _e`$ is the angular scale on which the sensitivity is reduced by a factor $`e`$ relative to its on-axis value. This corresponds to a half-power beamwidth, $`\mathrm{HPBW}=2.89^{}`$ (Heiles 1999).

We have chosen five representative photodissociation models, with parameters listed in Table 2, labeled by the letters A through E. Model A assumes the usual ISUV field strength, $`G_0=1`$, normal interstellar grains, $`A_v/N_H=(A_v/N_H)_0`$, and a mass loss rate, $`\dot{M}=3\times 10^5\mathrm{M}_{}\mathrm{yr}^1`$, which is in the middle of the range of mass loss rates measured for IRC+10216. Models B through E demonstrate the effect of departures from the values of Model A. The results of our photodissociation code for Models A through E are plotted in Figure 2. The solid line represents the fractional abundance of H2O relative to molecular hydrogen, $`x(\mathrm{H}_2\mathrm{O})`$, while the dashed line is the fractional abundance of OH relative to molecular hydrogen, $`x(\mathrm{OH})`$, both as a function of astrocentric radius. All models are normalized to a peak fractional water abundance of $`x(\mathrm{H}_2\mathrm{O})=1.1\times 10^6`$, consistent with the middle of the range of values for $`x(\mathrm{H}_2\mathrm{O})`$ determined by Melnick et al. (2001) from observations of the water vapor around IRC+10216. Notice that for Models B through D, even relatively large changes (factors of 2 or 3) in the values of $`G_0`$, $`A_v/N_H`$ and $`\dot{M}`$ have a relatively small effect on the radial location of the peak in $`x(\mathrm{OH})`$. The small value of $`G_0`$ in Model C and the large value of $`A_v/N_H`$ in Model D are fairly unrealistic (Jura 1983), but the plots of Models A through D demonstrate the range of possible sizes of the OH shell with even marginally realistic parameters. Indeed, the only photodissociation model which produces a shell significantly larger than Model A is Model E, which requires an implausibly large $`A_v/N_H`$ and a relatively improbable $`G_0`$.

Initially, we consider the line shape produced by OH distributed in a spherical shell. Based on photodissociation Model A, our line profile program predicted that the OH emission would come from an almost entirely unresolved spherical shell, and the line profile would thus be very nearly a “top hat” profile. The predicted profile is displayed in Figure 3 as line model S1, and the parameters for that line model as well as several other line models are listed in Table 3. After initial observations indicated the absence of detectable emission at $`26.5\mathrm{km}\mathrm{s}^1`$, the systemic velocity of the source (Cernicharo, Guélin & Kahane 2000), other models were discussed. One possibility we considered was that the OH shell was larger than we originally predicted, due either to a larger than expected dust opacity or to a weaker than expected ISUV field. If the OH shell were sufficiently enlarged, it would be partially resolved by AO and would produce a double-peaked line profile. Unfortunately, as discussed above, extremely unreasonable parameters are required to produce a substantially enlarged OH shell. Even using photodissociation model C, displayed in Figure 3 as line model S2, the line profile remains essentially a “top hat”. Only photodissociation Model E produces a double peaked line (line model S3 in Figure 3).

The final possibility we considered was an edge-on ring, which would also give rise to a double-peaked profile for the optically thin OH, even if the ring were unresolved. A ring structure for the OH would arise naturally if the water vapor were distributed in a disk, as might be expected for a cometary system analogous to the Kuiper Belt. Several ring profiles are plotted in Figure 4; the input parameters for the ring line models are listed in Table 3 together with the spherical line models. The rings are defined in spherical coordinates as extending from $`\theta =\frac{\pi }{2}\frac{\psi }{2}`$ to $`\theta =\frac{\pi }{2}+\frac{\psi }{2}`$, where $`\psi `$ is the opening angle of the disk/ring, taken as an adjustable parameter in our model; $`\psi =180`$degrees is a sphere. Note that changing the opening angle from $`10`$degrees (model R1) to $`1`$degree (model R2) has little effect on the line shape. We note that an unresolved expanding disk would lead to an elliptical profile for the optically thick water line. The water line profile obtained by Melnick et al. (2001) is plotted in Figure 1. It does not have a sufficiently large signal to noise ratio to distinguish between an elliptical profile and the parabolic profile expected for an expanding sphere. Finally, in model R3, we demonstrate the effect of changing the assumed microturbulent velocity, $`b`$, from the value measured previously for several other lines in IRC+10216, $`b=0.65\mathrm{km}\mathrm{s}^1`$ (Skinner et al. 1995), to $`b=2.0\mathrm{km}\mathrm{s}^1`$. Such a change increases the width of the peaks and tends to decrease the ratio of height of the line peak to the height of the line center. The microturbulent velocity of the OH gas could be increased by the injection of energy due to the photodissociation process itself.

Although the predicted line shape varies substantially depending on the parameters of the photodissociation model and geometry, the predicted line strength is primarily a function of the assumed fractional water abundance. For a given photodissociation model and geometry, the line strength varies linearly with assumed fractional water abundance. Unfortunately the fractional water abundance is uncertain by more than a factor of two, with possible values in the range $`x(\mathrm{H}_2\mathrm{O})=424\times 10^7`$ and a preferred value of $`x(\mathrm{H}_2\mathrm{O})=1.1\times 10^6`$ (Melnick et al. 2001). The large range of possible values for the fractional water abundance is due primarily to the uncertainty in the distance to IRC+10216, which is unknown to a factor of $`2`$. Various authors report reasonable distance estimates between $`300`$pc (Doty & Leung 1997) and $`110`$pc (Groenewegen, Van Der Veen & Matthews 1998). We adopt a distance of $`170`$pc (e.g. Winters, Dominik & Sedlmayr 1994; LeBertre 1997; Skinner et al. 1999) throughout. For line model S1 (based on photodissociation model A) and $`x(\mathrm{H}_2\mathrm{O})=1.1\times 10^6`$, we predict a line strength of $`10.5\mathrm{mJy}\mathrm{km}\mathrm{s}^1`$ (or $`86.1\mathrm{mK}\mathrm{km}\mathrm{s}^1`$) for the $`1667`$MHz line, However, due to the uncertainty in the fractional water abundance, this corresponds to possible line strengths ranging from $`3.82`$ to $`22.9\mathrm{mJy}\mathrm{km}\mathrm{s}^1`$ ($`31.3`$ to $`188\mathrm{mK}\mathrm{km}\mathrm{s}^1`$). The predicted line fluxes listed in Table 3 are for the $`1667`$MHz line and are calculated assuming $`x(\mathrm{H}_2\mathrm{O})=1.1\times 10^6`$. The $`1665`$MHz line is $`5/9`$ the strength of the $`1667`$MHz line in LTE for optically thin emission.

While the observed strengths of the $`1667`$ and $`1665`$MHz lines agree well with our line strength predictions, the observed OH line shapes are in very poor agreement with our predictions. The observed line profiles are neither “top hats” nor double-peaked; indeed, they are not even symmetric about $`v_{LSR}=26.5\mathrm{km}\mathrm{s}^1`$, the systemic velocity of IRC+10216. This asymmetry is unique333Observations by Lucas & Cernicharo (1989) and Lucas & Guilloteau (1992) found asymmetric line profiles for HCN in IRC+10216. However, this emission originates in the inner regions of the circumstellar outflow (inside the acceleration zone). The observed asymmetry was attributed to masing. \- every other molecule observed in the outer envelope of IRC+10216 has a profile which is symmetric about the systemic velocity (e.g. Cernicharo et al. 2000). Rather, the $`1667`$ and $`1665`$MHz emission features are blueshifted with respect to the systemic velocity, and they have central velocities which differ slightly (see Table 1). Based on the errors quoted in Table 1, we might conclude that the velocities of the two features differ at a statistically significant level. However, we note that the errors are based on fitting a Gaussian profile to a line which is not likely to be a true Gaussian. A more appropriate (and conservative) estimate of the error associated with the line center determination is simply the velocity resolution of our spectrum. For the $`1665`$ and $`1667`$MHz lines, the size of a resolution element after smoothing is $`1.1\mathrm{km}\mathrm{s}^1`$. Given this error estimate, we conclude that the $`1665`$ and $`1667`$MHz lines may have the same central velocity.

For comparison with our pre-observation models, we plot our “best fit” model over the data for the $`1667`$MHz line in Figure 5. Our “best fit” is model R3, except that we must assume a fractional water abundance of $`2.4\times 10^6`$ in order to match the line fluxes for velocities between $`40`$ and $`30\mathrm{km}\mathrm{s}^1`$. Also, in order to match the velocity of the observed line with the blueshifted peak of the model, we assume a $`v_{LSR}=23.5\mathrm{km}\mathrm{s}^1`$, which represents a redshift of $`3\mathrm{km}\mathrm{s}^1`$ with respect to prior measurements of the $`v_{LSR}`$ using other species (e.g. Cernicharo et al. 2000). It is obvious from the plot that the blueshifted peak of the model fits the data reasonably well, but the redshifted peak of the model is entirely missing from the data.

The observed shapes of the $`1667`$ and $`1665`$MHz OH lines are quite unexpected. There are only three general mechanisms which might create the observed line profiles; these mechanisms are: 1) absorption, 2) masing or 3) an intrinsic asymmetry in the spatial distribution of the OH. In the following sections, we examine each mechanism in detail and discuss its plausibility.

### 3.1 Absorption

It seems reasonable to begin with the assumption that the OH emission from IRC+10216 is intrinsically symmetric. If that were the case, then in order to explain the OH line shapes, there must exist some absorbing medium within the circumstellar envelope of IRC+10216 which prevents emission from the far (redshifted) side of the envelope from reaching us, while emission from the near (blueshifted) side escapes and is detected. Since these observations were done at wavelengths which are large, relative to previous observations of this star, there is one source of opacity that could have previously gone undetected - electron free-free absorption. However, this absorption would be accompanied by thermal continuum emission. In order to be an effective absorber of the redshifted half of the OH line, the absorber would need to be optically thick at a radius of $`10^{17}`$cm, where the OH distribution peaks (though the peak radius of the OH distribution varies somewhat depending on the parameters of the photodissociation model, $`10^{17}`$cm is a reasonable value for most models). At such a large radius, the envelope gas temperature is quite cold, about $`11`$K, and the thermal emission from the absorber would correspond to a blackbody at the same temperature. Unfortunately, our observations are not well suited to determining the continuum brightness at $`1665`$ and $`1667`$MHz, but Sahai, Claussen, & Masson (1989) have observed the continuum brightness at slightly higher frequencies, specifically at $`15`$ and $`20`$GHz. The predicted flux densities due to free-free emission at $`15`$ and $`20`$GHz are $`250`$ and $`140`$mJy, respectively, for an object which is optically thick at $`1665`$ and $`1667`$MHz. Sahai et al. (1989) find the flux densities to be only $`6`$ and $`1.4`$mJy at $`15`$ and $`20`$GHz, respectively. The absence of strong continuum emission indicates that at these frequencies, the envelope of IRC+10216 cannot be optically thick at radii large compared to that of the AO beam, effectively ruling out the absorption scenario. In addition, the electron density required to make the envelope of IRC+10216 optically thick would be substantially in excess of the densities presently expected based on theoretical models of the envelope chemistry.

### 3.2 Masing

A second possible explanation for the observed line profiles is maser amplification of the $`1667`$ and $`1665`$MHz lines. This possibility is attractive, because it is the only one which could potentially be responsible for a velocity difference between the $`1667`$ and $`1665`$MHz lines. While the modeling of circumstellar maser emission is difficult and a detailed model of possible OH masers around IRC+10216 is beyond the scope of this paper, we can compare our observations to other asymptotic giant branch (AGB) stars with known OH maser emission. However, we should be cautious about drawing strong conclusions based on this comparison, due to the fact that all stars with confirmed circumstellar OH masers are oxygen-rich stars with much higher fractional OH abundances or carbon-rich stars which have recently made the transition from the oxygen-rich state and also have higher fractional OH abundances than IRC+10216.

There is at least one star which is known to have exhibited OH emission with velocity characteristics similar to IRC+10216, albeit at much larger flux densities. In observations reported by Etoka et al. (2001), the star R Crt displayed maser emission in the $`1667`$ and $`1665`$MHz main lines, but not in the $`1612`$MHz satellite line; this in itself is not particularly unusual among stars which exhibit circumstellar OH masers. More importantly, in observations on 1990 November 17, the $`1667`$ and $`1665`$MHz emission appeared only at velocities blueshifted with respect to the systemic velocity of R Crt, but with velocities which were offset from each other by $`1\mathrm{km}\mathrm{s}^1`$. Significantly, the ratio of the $`1667`$ to $`1665`$MHz line strengths was $`1.7`$ during this observation, not far from the ratio of $`1.8`$ expected in thermal equilibrium. As the OH emission from R Crt is known to be maser amplified, this example would seem to argue for at least the possibility of masing in IRC+10216.

Several objections may be raised against the masing scenario, however, which would need to be resolved before masing could be accepted as the source of the line profile asymmetry in IRC+10216. The first objection is that the predicted fluxes for IRC+10216 of the $`1667`$ and $`1665`$MHz lines in LTE from line model S1 are $`1.3`$ and $`0.73\mathrm{mJy}\mathrm{km}\mathrm{s}^1`$ respectively, for velocities between $`40`$ and $`30\mathrm{km}\mathrm{s}^1`$ and an $`x(\mathrm{H}_2\mathrm{O})=4\times 10^7`$. These predictions are only a factor of $`4`$ less than the observed line strengths, even assuming the minimum possible $`x(\mathrm{H}_2\mathrm{O})`$, and using line model S1, which predicts the smallest flux density for velocities between $`40`$ and $`30\mathrm{km}\mathrm{s}^1`$. The substantial blue/red asymmetry of the line flux requires strong masing. The asymmetry of the $`1667`$MHz line can be quantified by fitting a Gaussian with fixed parameters to the redshifted side of the line and taking a ratio of the $`3\sigma `$ upper limit on redshifted flux, $`F_{red}`$, with the measured flux from the blueshifted side of the line, $`F_{blue}=39.0\mathrm{mK}\mathrm{km}\mathrm{s}^1`$. For a Gaussian with a line center fixed at $`v_{LSR}=10.9\mathrm{km}\mathrm{s}^1`$ and a fixed width of $`5.8\mathrm{km}\mathrm{s}^1`$, we find a flux of $`4.38\pm 11.4\mathrm{mK}\mathrm{km}\mathrm{s}^1`$, $`3\sigma `$. We chose the center of the red feature so that the blue and red sides of the line would be symmetric about $`v_{LSR}=23.5\mathrm{km}\mathrm{s}^1`$. Thus, we find a $`3\sigma `$ upper limit $`F_{red}/F_{blue}0.4`$, or an amplification of the blue side of the line by a factor of at least $`2.5`$, clearly requiring significant amplification. If the OH emission is substantially amplified by maser action, the $`1667`$ and $`1665`$MHz lines ought to have a brightness temperature comparable to or larger than the temperature of the microwave background ($`3`$K, or about $`3700\mathrm{mJy}\mathrm{km}\mathrm{s}^1`$ between $`40`$ and $`30\mathrm{km}\mathrm{s}^1`$). This might lead one to conclude, based on our measurements, that the observed $`1667`$ and $`1665`$MHz lines are too dim to be strongly maser amplified. Although the OH emission is not that bright, if maser amplification were confined to a small spot, the beam dilution in AO’s large field of view could mean that our observations are still consistent with amplification. Specifically, the maser spot would need to have a diameter smaller than $`10`$arcsec or $`1700`$AU at a distance of $`170`$pc, assuming an intrinsic brightness of $`3`$K or larger. This is not particularly small, considering that OH maser spots around other evolved stars can have diameters smaller than $`25`$AU (Chapman, Cohen & Saikia 1991), so the dim lines do not exclude masing.

A far stronger objection stems from theoretical models of circumstellar OH masers by Elitzur (1978). Elitzur found that OH masers in oxygen-rich stars were created by pumping from the far-IR emission of warm dust, and that for the $`1667`$ and $`1665`$MHz main lines, the maser gain is proportional to the density of OH molecules (below a threshold density). Although Elitzur (1978) did not model OH densities or gas temperatures as low as those expected in IRC+10216, a reasonable extrapolation of his models indicates that $`|\tau |`$ would be $`1`$, where $`\tau `$ is the line center optical depth.

### 3.3 Intrinsic Asymmetry

The only remaining general class of explanations for the OH line asymmetry is an intrinsic asymmetry in the spatial distribution of the OH around IRC+10216. There are two possible mechanisms for producing such a spatially asymmetric distribution. One mechanism is an asymmetric ISUV field in the vicinity of IRC+10216; this would produce a larger fractional abundance of photodissociation products on one side of the circumstellar envelope relative to the other. Since OH is a photodissociation product, this scenario could produce an asymmetric line profile if the viewing geometry were properly aligned. However, all other molecular photodissociation products would then have the same line profile as OH. This is not observed (e.g. Cernicharo et al. 2000).

The other way to produce an asymmetric distribution of OH would be an asymmetric distribution of its parent molecule, H2O. The observed water vapor emission is essentially symmetric about the systemic velocity of IRC+10216. However, the H2O emission and the OH emission do not arise from precisely identical spatial regions (see Figure 2). The time required for a parcel of gas to travel a given radial distance can be found from the simple relationship $`R=v_{out}t`$, where $`v_{out}`$ is the terminal velocity of the stellar outflow. Due to the direct relationship between the radial coordinate and time ($`R=v_{out}t`$), the fact that the OH and H2O are not spatially coincident means that the present distributions of OH and H2O were formed at two different epochs, and the OH and H2O need not have identical geometries at the present time. If the spatial distribution of water were asymmetric a few hundred years ago, but is symmetric now, the OH and H2O observations could be reconciled with the intrinsic asymmetry scenario.

One possible explanation for a change in the distribution of water vapor over such a relatively short timescale involves the time dependence of the stellar mass loss rate. The circumstellar envelope of IRC+10216 contains multiple dust shells which are the result of episodic mass loss. Detailed images of these shells by Mauron & Huggins (1999, 2000), show that the shells are incomplete - they are arcs rather than rings. Additionally, the dust distribution in the innermost regions of the envelope ($`r15^{\prime \prime }`$) is an extremely asymmetric bipolar outflow. The timescale on which these dust shells are produced is of order a few hundred years. A new mass loss study by Fong, Meixner & Shah (2003) demonstrates that mass loss in IRC+10216 has been ongoing for at least $`7000`$years and that mass loss occurred in a clumpy and asymmetric fashion. Rapid variations in the mass loss rate are also possible; recent high-resolution IR images show dust distribution asymmetries on size scales as small as tens of milliarcseconds and imply mass loss variations on timescales of years (Tuthill et al. 2000; Weigelt et al. 2002; Men’shchikov, Hofmann & Weigelt 2002). Ford & Neufeld (2001) showed that the vaporization of water from icy bodies orbiting IRC+10216 is sensitively dependent on the stellar flux impinging on an icy body. Since the dust distribution is highly asymmetric, the luminosity of the central star will “leak” out in an asymmetric and probably quite patchy pattern, after multiple scatterings and dust absorption and re-emission. Thus, the flux impinging on different icy bodies at the same astrocentric radius will not be necessarily be the same; it will depend on the particular geometry of radiation leakage through the circumstellar dust. Therefore, the water vaporization rate from the surface of a particular icy body (and the distribution of water vapor in the circumstellar outflow) will depend very substantially on the time varying dust distribution in the circumstellar outflow. Note that this scenario is still capable of producing the water line profile observed by Melnick et al. (2001). If the water vapor is more or less off to one side of the star (from a terrestrial observer’s perspective), there will be roughly equal amounts of material coming toward us and moving away from us, still yielding a roughly symmetric profile. Additionally, since the water line is optically thick, it is probable that some clumpiness or patchiness in the water distribution would be smoothed out in the line profile.

## 4 Summary

We have detected the $`1665`$ and $`1667`$MHz lines of OH in emission towards the carbon-rich AGB star IRC+10216. This detection supports the identification of the $`556.936`$GHz emission detected toward IRC+10216 by Melnick et al. (2001) with the $`1_{10}1_{01}`$ transition of water vapor. The $`1665`$ and $`1667`$MHz line profiles are unexpectedly narrow and blueshifted with respect to the systemic velocity of IRC+10216. We believe that the unanticipated line shape is the result of either maser amplification or some intrinsic spatial asymmetry in the OH distribution around IRC+10216, or possibly both. While we favor the spatial asymmetry explanation, due to early theoretical work by Elitzur (1978), further modeling of the excitation conditions of OH around IRC+10216 using a code which includes strong masing would help clarify the possible role of maser amplification in this source. Additional theoretical work may also be required to understand the potential spatial asymmetry of OH, despite the apparent symmetry of its parent molecule, H2O. Finally, while it would be most desirable to spatially resolve the OH emission in IRC+10216, such observations cannot presently be completed within an acceptable amount of telescope time.

We are grateful to the staff at AO, particularly Chris Salter, Phil Perillat, Karen O’Neil and Mike Nolan. We are especially grateful for the outstanding service observing work performed by Hector Hernandez. We would also like to thank Jennifer Wiseman, David Hollenbach, Al Glassgold and Barry McKernan for helpful discussions. This work was supported by subcontract SV252005 from the Smithsonian Institution; K.E.S.F. was also supported by an American Dissertation Fellowship from the American Association of University Women (AAUW). The National Astronomy and Ionosphere Center is operated by Cornell University under a cooperative agreement with the National Science Foundation.

## Appendix A Correction to Photodissociation Rates due to Spherical Symmetry

For the purposes of work described in this paper, we want a simple expression for the photodissociation rates of H2O and OH as a function of depth, measured in magnitudes of visual extinction, $`A_v`$, in a sphere which is illuminated isotropically from the exterior by the ISUV field. Calculations done by Roberge et al. (1991) (R91 hereafter) provide rates for a plane-parallel slab geometry. We describe here a method which can be applied to the rates calculated with plane-parallel geometry which will give approximate rates for the spherical case.

R91 calculated photodissociation and photoionization rates for a number of species, as a function of visual extinction into plane-parallel clouds illuminated isotropically from both sides by the average ISUV field (Draine 1978). In general, the photodestruction rate, $`\mathrm{\Gamma }`$, for a species, $`i`$, is

$$\mathrm{\Gamma }_i(A_v)=4\pi _{\lambda _H}^{\lambda _i}J_\lambda (A_v)\sigma _i(\lambda )𝑑\lambda .$$
(A1)
The rate, $`\mathrm{\Gamma }_i`$, depends on the extinction, $`A_v`$, measured inward from the surface of the cloud in units of visual magnitudes. Here $`\sigma _i(\lambda )`$ is the wavelength dependent cross section for the process of interest, and the integration runs from $`\lambda _H=91.2`$nm to the threshold wavelength $`\lambda _i`$ for the process of interest. The quantity $`J_\lambda (A_v)`$ is the mean intensity of radiation at depth $`A_v`$ in photons cm-2 s-1 sr-1 nm-1. R91 used the spherical harmonics method to solve numerically the radiative transfer equation for $`J_\lambda (A_v)`$ for clouds of $`A_v^{tot}=1,10`$ and $`100`$mag. The solution for $`J_\lambda (A_v)`$ allowed R91 to calculate the photodestruction rates for a large number of species at various depths $`A_v`$ into each of their model clouds ($`A_v^{tot}=1,10`$ and $`100`$mag). In addition, R91 found that a bi-exponential function, of the form

$$\mathrm{\Gamma }_i(A_v)=C_i\mathrm{exp}(\alpha _iA_v\beta _iA_v^2)\mathrm{s}^1$$
(A2)
provided a reasonable fit to their numerical results, where $`C_i`$ is the unattenuated photodestruction rate and $`\alpha _i`$ and $`\beta _i`$ are fit parameters, for a given species $`i`$.

A point $`P`$ in the interior of an R91 plane-parallel slab is illuminated by radiation coming from many directions. If we consider a slab with $`A_v=100`$, we can assume that near one surface of the slab, essentially all photodissociating radiation originates at the nearer surface. Under such conditions, the point $`P`$ “sees” rays diminished by an “effective visual extinction”, $`A_v^{}=A_v/\mathrm{cos}\theta `$ coming from a direction $`\theta `$, where the angle is measured from the normal to the slab and $`A_v`$ is the depth measured from the surface of the cloud to the point $`P`$ along the normal ray. In order to calculate photodissociation rates valid for a point $`Q`$ interior to a sphere, we realized that we should decompose the photodissociating radiation arriving at point $`P`$ inside a slab into constituent rays coming from a direction $`\theta `$. This would give us the photodissociation rate due to a single ray. Then, we can use the result for a single ray to integrate over all rays arriving at the point $`Q`$ in the interior of the sphere. Such a point $`Q`$ “sees” rays diminished by a different “effective visual extinction”, $`A_v^{\prime \prime }=\theta A_v/\mathrm{sin}\theta `$. We can then integrate over all rays arriving at point $`Q`$, to determine the total photodissociation rate for a point internal to a sphere.

So, given a formula for the photodissociation rate $`\mathrm{\Gamma }_i^{pp}=C_i\mathrm{exp}(\alpha _i^{pp}A_v\beta _i^{pp}A_v^2)`$ in a plane-parallel slab, we should be able to find best fit values for the parameters $`\alpha _i^r`$ and $`\beta _i^r`$ subject to the condition that

$$\frac{\mathrm{\Gamma }_i^{pp}\mathrm{\Gamma }_i^{pr}}{\mathrm{\Gamma }_i^{pp}}0.2$$
(A3)
for all $`A_v`$ of interest, where

$$\mathrm{\Gamma }_i^{pr}=C_i_0^{\pi /2}\mathrm{exp}\left(\alpha _i^rA_v^{}\beta _i^r\left(A_v^{}\right)^2\right)\mathrm{sin}\theta d\theta .$$
(A4)
Here $`\mathrm{\Gamma }_i^{pr}`$ is the best fit photodissociation rate at a point $`P`$ interior to a plane-parallel slab, found by integrating over each individual ray terminating at point $`P`$. The parameters $`\alpha _i^r`$ and $`\beta _i^r`$ are simply the attenuation parameters along any single ray. Once we have $`\alpha _i^r`$ and $`\beta _i^r`$, we can repeat a similar procedure to find a best fit $`\alpha _i^{sp}`$ and $`\beta _i^{sp}`$. Specifically, we wish to satisfy the condition

$$\frac{\mathrm{\Gamma }_i^{sr}\mathrm{\Gamma }_i^{sp}}{\mathrm{\Gamma }_i^{sr}}0.2$$
(A5)
for all $`A_v`$ of interest, where

$$\mathrm{\Gamma }_i^{sr}=C_i_0^\pi \mathrm{exp}\left(\alpha _i^rA_v^{\prime \prime }\beta _i^r\left(A_v^{\prime \prime }\right)^2\right)\mathrm{sin}\theta d\theta $$
(A6)
and

$$\mathrm{\Gamma }_i^{sp}=2C_i\mathrm{exp}(\alpha _i^{sp}A_v\beta _i^{sp}A_v^2).$$
(A7)
Here $`\mathrm{\Gamma }_i^{sr}`$ is the photodissociation rate at a point $`Q`$ interior to a sphere, which is found by integrating over all single rays arriving at point $`Q`$ (calculated using $`\alpha _i^r`$ and $`\beta _i^r`$ which were found previously to fit the attenuation along a single ray). $`\mathrm{\Gamma }_i^{sp}`$ is a simple, though approximate, expression for the photodissociation rate at a point $`Q`$ interior to a sphere, and $`\alpha _i^{sp}`$ and $`\beta _i^{sp}`$ are the best fit attenuation parameters for points interior to a sphere.

In the particular situation of a circumstellar envelope, we note that most of the photodissociation of H2O and OH occurs at low $`A_v`$, near the surface of the sphere, and thus it is most important that $`\mathrm{\Gamma }_i^{sp}`$ be accurate for relatively small $`A_v`$. To that end, we fit our own bi-exponential to the tabulated photodissociation rates at small $`A_v`$ from R91, which yielded the $`\alpha _i^{pp}`$ and $`\beta _i^{pp}`$, listed in Table 4. We then applied the fitting procedure described above and found values for $`\alpha _i^r`$, $`\beta _i^r`$, $`\alpha _i^{sp}`$ and $`\beta _i^{sp}`$, also listed in Table 4.


format: text

text:
# Curvaton reheating: an application to braneworld inflation

## I Introduction

The curvaton lyth1 ; lyth2  is an interesting new proposal for explaining the observed large-scale adiabatic density perturbations in the context of inflation. Perturbations are acquired by a field other than the inflaton during inflation, whose energy density is subdominant during inflation but which comes to dominate sometime afterwards, at which point its initial isocurvature perturbations convert to adiabatic moll ; lyth1  (see also TT ). Decay of the curvaton to conventional matter locks the perturbations in as adiabatic, though there can be complex phenomenology if some species have already decoupled from the thermal bath by the time this happens lyth2 ; ES .

An aspect of this process which has been implicit in most papers is that this process also offers a new mechanism of reheating. The first detailed discussion of this aspect appeared recently in Ref. FL  by Feng and Li, who named it curvaton reheating. The curvaton is a type of matter whose energy density is not diluted during inflation, and some or all of the present material in the Universe may have survived inflation via the curvaton. This reheating mechanism complements the other two known mechanisms, the standard one being decay of the inflaton energy density into conventional matter and the rarely-used alternative being gravitational particle production at the end of inflation gravpart .

In this paper we will explore an application of curvaton reheating to a particular model of braneworld inflation known as steep inflation cope . This models brings the inflationary epoch to an end in a novel way. The potential chosen (for example a steep exponential) is so steep that it would not drive inflation in the standard cosmology, but it may do so in the presence of high-energy braneworld corrections to the Friedmann equation, for example as found in the Randall–Sundrum Type II scenario. As the energy density decreases, those corrections become unimportant, and the scalar field makes a transition into a kinetic energy dominated regime bringing inflation to an end. As the inflaton survives this process without decay, an alternative reheating mechanism is required.

In Ref. cope , it was proposed that reheating took place via gravitational particle production gravpart ; MJ1 . However this brings problems to the scenario. Because of the inefficiency of the process, there is a long kinetic energy dominated regime, which allows short-wavelength gravitational waves to reach excessive amplitudes sahni . There are likely also to be problems with the low reheat temperature achieved, which is barely sufficient to allow standard nucleosynthesis to proceed. The scenario also suffers from a separate problem, which is that it tends to predict perturbation spectra which are too far from scale-invariance, and too high an amplitude of long-wavelength gravitational waves; for the exponential potential the predicted spectra are in conflict with the most recent observations.

In this paper, we propose that implementation of curvaton reheating in place of gravitational particle production can redress those difficulties. Curvaton reheating can be much more efficient, allowing a high reheat temperature and preventing short-wavelength gravitational wave dominance. Further, adiabatic perturbations produced via the curvaton may be closer to scale-invariance, and in particular the contribution of gravitational waves to the microwave anisotropies may be considerably suppressed. Two related papers appeared as we were completing this work; Ref. FL  discuss curvaton reheating in a quintessential inflation model PV  in the standard cosmology, while Ref. Dim  considers adding the curvaton to the steep inflation model, but focussing on particle physics model building aspects rather than observational phenomenology and constraints.

A summary of this paper is as follows. In Section II, we review the main features of steep inflation in the braneworld scenario. We present only the main results of these kind of models, including the production of primordial gravitational waves. In Section III, we study the evolution of a curvaton field and its perturbations in steep inflation, and we focus our attention on all the constraints to be applied upon the free parameters of the model. We consider separately the cases that arise depending on whether the curvaton field dominates after or before it decays into radiation. In Section IV we put all constraints together and analyze the parameter space of the model of curvaton field $`+`$ steep inflation which satisfies all the requirements. In the final section we draw some conclusions and topics to be investigated in future publications.

## II The Steep Inflation model

In this section we give a brief description of an inflationary model in the braneworld scenario, assuming that the inflaton field is a minimally coupled real scalar field $`\varphi `$ endowed with an exponential potential of the form

$$V(\varphi )=V_0e^{\alpha \sqrt{\kappa _0}\varphi },$$
(1)
where $`V_0,\alpha `$ are free parameters and $`\kappa _0=8\pi G=8\pi m_{\mathrm{Pl}}^2`$ with $`m_{\mathrm{Pl}}`$ the Planck mass. More details of the model can be found in Refs. cope ; sahni . A somewhat related idea has been discussed in Ref. guendel .

### II.1 Inflationary dynamics

In the $`(4+1)`$-dimensional brane scenario inspired by the Randall–Sundrum model, the standard Friedman equation is modified to

$$H^2=\frac{\kappa _0}{3}\rho \left(1+\frac{\rho }{2\lambda _\mathrm{b}}\right),$$
(2)
where $`\lambda _\mathrm{b}`$ is the brane tension and, taking the inflaton field to be the dominant component, $`\rho =\dot{\varphi }^2/2+V(\varphi )`$. The scalar field is confined to the brane so that its equation of motion is the usual Klein–Gordon one

$$\ddot{\varphi }+3H\dot{\varphi }+\frac{dV}{d\varphi }=0,$$
(3)
where dots denote temporal derivatives.

For large energies $`\rho \lambda _\mathrm{b}`$, the inflaton field experiences an increased friction as it rolls down its potential, and then one can expect inflation to occur even for steep potentials (i.e. even for $`\alpha 1`$). This can be seen from the slow-roll parameters for the exponential potential MWBH , which for large energies read

$$ϵ=\eta =2\alpha ^2\frac{\lambda _\mathrm{b}}{V}.$$
(4)
and thus, inflation can occur whenever $`V>2\alpha ^2\lambda _\mathrm{b}`$.

The equations of motion can be solved exactly in the slow-roll approximation, from which we find the following expressions for the scale factor $`a`$, the Hubble factor $`H`$ and the inflaton field $`\varphi `$, respectively,

$`\mathrm{ln}{\displaystyle \frac{a}{a_\mathrm{i}}}`$ $`=`$ $`(𝒩+1)e^{H_\mathrm{f}(tt_\mathrm{f})},`$ (5)
$`H`$ $`=`$ $`H_\mathrm{f}e^{H_\mathrm{f}(tt_\mathrm{f})},`$ (6)
$`\alpha \sqrt{\kappa _0}\varphi (t)`$ $`=`$ $`H_\mathrm{f}(tt_\mathrm{f})\mathrm{ln}{\displaystyle \frac{V_\mathrm{f}}{V_0}}.`$ (7)
Here ‘i(f)’ indicates the initial (final) values of the different quantities during inflation, and $`𝒩`$ is the number of $`e`$-foldings as can be seen from Eq. (5). Observe that, in the braneworld scenario, the slow-roll approximation for the exponential potential is a much better approximation than in standard cosmology.

Inflation ends by violation of the slow-roll conditions, when $`ϵ\eta 1`$, and so the value of the potential energy at the end of inflation is $`V_\mathrm{f}2\alpha ^2\lambda _\mathrm{b}`$. Using Eqs. (5)–(7), we find that the initial and final values of the Hubble factor and the potential energy, respectively, are related through

$$\frac{H_\mathrm{i}}{H_\mathrm{f}}=\frac{V_\mathrm{i}}{V_\mathrm{f}}𝒩+1,$$
(8)
where

$$H_\mathrm{f}=\alpha ^2\sqrt{\frac{2\lambda _\mathrm{b}\kappa _0}{3}}.$$
(9)
So far the model has three free parameters, but only two of them are relevant for inflation, the brane tension $`\lambda _\mathrm{b}`$ and the self-interaction of the inflaton field $`\alpha `$.

The amplitude of scalar density perturbations can be calculated for the exponential potential, and it is given by cope ; sahni

$$A_\mathrm{s}^2\frac{8}{75}\frac{V_\mathrm{i}^4}{\alpha ^2\lambda _\mathrm{b}m_{\mathrm{Pl}}^4}=\frac{128}{75}\alpha ^6(𝒩+1)^4\frac{\lambda _\mathrm{b}}{m_{\mathrm{Pl}}^4}.$$
(10)
In the usual scenario of steep inflation, the values of the inflation parameters can be adjusted so that we get the observed value $`A_\mathrm{s}^24\times 10^{10}`$. In particular, we can determine the value of the brane tension in terms of $`\alpha `$ as

$$\lambda _\mathrm{b}\frac{2.3\times 10^{10}}{\alpha ^6}\frac{m_{\mathrm{Pl}}^4}{(𝒩+1)^4}.$$
(11)
Notice that any deviation of $`\lambda _\mathrm{b}`$ from the value obtained in Eq. (11) would imply either a larger or smaller amplitude for the scalar spectrum.

As was pointed out in Ref. sahni , the number of $`e`$-foldings can be unambiguously determined for the exponential potential, for which we obtain $`𝒩70`$. This result is based on Eq. (10) and is independent of the value of $`A_\mathrm{s}^2`$.

### II.2 After inflation

Soon after inflation ends, the brane term becomes unimportant and we recover the standard Friedman equation. Hence, the friction on the inflaton field diminishes and the inflaton energy is dominated by the kinetic term and the scalar energy density falls as stiff matter, $`\rho _\varphi =\rho _\varphi ^{(\mathrm{kin})}(a_{\mathrm{kin}}/a)^6`$. The epoch that follows is called the ‘kinetic epoch’ or ‘kination’ MJ2 , and we shall use ‘kin’ to label the value of the different quantities at the beginning of the kinetic epoch.

Since we are now working within standard cosmology (at energies $`\rho <\lambda _\mathrm{b}`$), the Hubble factor evolves as

$`H=H_{\mathrm{kin}}\left({\displaystyle \frac{a_{\mathrm{kin}}}{a}}\right)^3,`$ $`H_{\mathrm{kin}}^2{\displaystyle \frac{\kappa _0}{3}}\rho _\varphi ^{(\mathrm{kin})}.`$ (12)
The kinetic regime does not commence immediately after the end of inflation, but the values of the Hubble parameter at the end of inflation and at the beginning of the kinetic regime are related by the formula sahni

$$\frac{H_{\mathrm{kin}}}{H_\mathrm{f}}=0.085\frac{0.688}{\alpha ^2}.$$
(13)
A small amount of radiation has been produced quantum-mechanically during inflation with an energy density $`\rho _\mathrm{r}0.01g_pH_\mathrm{f}^4`$, where $`g_p=10100`$ is the number of different species created from the vacuum. If the inflaton potential does not decay just after inflation, then this is the only mechanism to reheat the Universe.

With the density of scalar matter falling more rapidly than that of radiation, the Universe continues evolving according to Eq. (12) until the small amount of radiation comes to dominate. Since the ratio of scalar matter to radiation at the end of inflation is large, $`(\rho _\varphi /\rho _\mathrm{r})_\mathrm{f}1`$, radiation comes to dominate the expansion of the Universe only after a very long kinetic epoch. This is troublesome since the energy density of gravitational waves produced quantum-mechanically can exceed the nucleosynthesis constraints, as shown in Ref. sahni .

In braneworld inflation, the amplitude of gravitational waves $`h_{\mathrm{GW}}^2`$ is enhanced with respect to the standard scenario, and its value at the end of inflation is given by

$$h_{\mathrm{GW}}^23\alpha ^2(𝒩+1)\frac{H_\mathrm{i}^2}{m_{\mathrm{Pl}}^2}.$$
(14)
Gravitational waves behave as massless scalar fields and then their amplitude remains constant during inflation. During the kinetic epoch, the energy density of gravitational waves evolves as sahni

$$\rho _g=\frac{32}{3\pi }h_{\mathrm{GW}}^2\rho _\varphi \left(\frac{a}{a_{\mathrm{kin}}}\right)^2.$$
(15)
From this, the energy density of gravitational waves at the time of scalar stiff matter–radiation equality $`(\rho _\varphi =\rho _\mathrm{r})`$ is

$$\frac{\rho _g}{\rho _\mathrm{r}}|_{a=a_{\mathrm{eq}}}=\frac{64}{3\pi }h_{\mathrm{GW}}^2\left(\frac{a_{\mathrm{eq}}}{a_{\mathrm{kin}}}\right)^2\frac{64}{3\pi }h_{\mathrm{GW}}^2\frac{\rho _\varphi }{\rho _\mathrm{r}}|_\mathrm{f}1.$$
(16)
Thus, the energy density of gravitational waves on short wavelength scales sahni  exceeds the contribution of radiation. It happens then that gravitational waves come to dominate the stiff scalar matter well before radiation, which in fact never dominates, and then nucleosynthesis constraints are not satisfied at all. Therefore, the original steep inflation model seems to be ruled out sahni .

However, the production of gravitational waves can be kept under control if the kinetic epoch is shortened. The latter can be done either by a sufficiently early decay of the inflaton field or by having a post-inflationary equation of state softer than that of stiff matter. The latter picture is quite appealing since the surviving inflaton field could become part of the dark matter at late times sahni ; lidsey ; FJ .

## III The curvaton field

In this section we explore the possibility that the inflaton field is not responsible for providing the primordial fluctuations required for the formation of structure in the late Universe, as in the case where the value of the brane tension is lower than that in Eq. (11). Instead, the primordial fluctuations are generated by a second scalar field $`\sigma `$, usually called the curvaton field, through initially isocurvature perturbations lyth1 . For simplicity, we assume here that the curvaton field has a quadratic scalar potential $`U(\sigma )=m^2\sigma ^2/2`$, and hence obeys the Klein–Gordon equation

$$\ddot{\sigma }+3H\dot{\sigma }+m^2\sigma =0,$$
(17)
where $`m`$ is the curvaton mass.

We will now follow the evolution of the curvaton field through different stages and find the constraints upon the free parameters of the model in order to have a viable curvaton scenario together with steep inflation. In doing this, we shall follow a similar procedure to that in Ref. nicola . We begin with a brief overview of the dynamics.

First of all, it is assumed that the curvaton field coexists with the inflaton field during inflation, during which the inflaton energy density is the dominant component. The curvaton field should survive the rapid expansion of the Universe, and for that it has to be effectively massless. As we shall see below, this condition is translated into a constraint of the curvaton mass $`m`$. If this constraint is satisfied, the curvaton field will remain constant to its initial value.

The following stage is that when the curvaton field becomes effectively massive, and this will generally happen during the kinetic epoch. In order to prevent a period of curvaton-driven inflation, the Universe must remain inflaton-dominated until this time. The latter condition imposes a constraint on the initial values the curvaton field can take, $`\sigma _\mathrm{i}`$. Once effectively massive, the curvaton field oscillates around the minimum of its potential and its energy density evolves as non-relativistic matter.

The final stage is that of curvaton decay into radiation, and the standard big bang cosmology is recovered afterwards. In the general case, curvaton decay should occur before nucleosynthesis, but other constraints arise depending on epoch of the decay, governed by the decay parameter $`\mathrm{\Gamma }_\sigma `$. There are two scenarios to be considered, depending whether the curvaton field decays before or after it becomes the dominant component of the Universe. Whatever the case, the energy from gravitational waves should be maintained under control, which will actually become an extra constraint involving different parameters of the model. Also, we calculate the curvaton perturbations for each case which will help us to reduce the number of free parameters. At this point, the gaussianity conditions on the perturbations will be taken into account too.

The introduction of the curvaton field adds three free parameters to the steep inflation model: the curvaton mass $`m`$, the initial value of the curvaton field $`\sigma _\mathrm{i}`$ and the curvaton decay parameter $`\mathrm{\Gamma }_\sigma `$. However, a successful model needs to satisfy the conditions outlined above which are strong enough to restrict the allowed curvaton parameters.

### III.1 Inflation

We first study the evolution of the curvaton field during braneworld inflation with an exponential inflaton potential as described in Section II. We change the independent variable to $`z3H/H_\mathrm{f}`$, for which the equation of motion of the curvaton field Eq. (17) reads

$$z^2\sigma ^{\prime \prime }+z(1z)\sigma ^{}+\frac{m^2}{H_\mathrm{f}^2}\sigma =0,$$
(18)
where primes denote derivatives with respect to $`z`$.

Assuming that the curvaton field is effectively massless during inflation $`mH_\mathrm{f}`$, the last term in Eq. (18) is small compared to the others, so that the curvaton equation of motion is

$$z^2\sigma ^{\prime \prime }+z(1z)\sigma ^{}0,$$
(19)
whose exact solutions are

$`\sigma (z)`$ $`=`$ $`\sigma _\mathrm{i}+\sigma _\mathrm{i}^{}z_\mathrm{i}e^{z_\mathrm{i}}\left[\mathrm{Ei}(z)\mathrm{Ei}(z_\mathrm{i})\right],`$ (20)
$`\sigma ^{}`$ $`=`$ $`\sigma _\mathrm{i}^{}{\displaystyle \frac{z_\mathrm{i}}{z}}\mathrm{exp}\left(zz_\mathrm{i}\right).`$ (21)
Here, $`\sigma _\mathrm{i}`$ and $`\sigma _\mathrm{i}^{}`$ are the initial values of the field and its derivative, and $`\mathrm{Ei}(z)`$ is the exponential integral. In the case that $`𝒩=70`$, so that $`z`$ runs in the range $`z_\mathrm{i}=213z_\mathrm{f}=3`$ (see Eq. (8)), the ratio of the final value of the curvaton field to its initial value is

$$\frac{\sigma _\mathrm{f}}{\sigma _\mathrm{i}}=11.005\frac{\sigma _\mathrm{i}^{}}{\sigma _\mathrm{i}}.$$
(22)
We see that the final value $`\sigma _\mathrm{f}`$ depends on the initial conditions. The interesting case is that in which the curvaton energy is equally distributed into kinetic and potential at the beginning of inflation, so that we can take $`\dot{\sigma _\mathrm{i}}^2m^2\sigma _\mathrm{i}^2`$. Thus, from $`\dot{\sigma _\mathrm{i}}=H_\mathrm{f}z_\mathrm{i}\sigma _\mathrm{i}^{}`$ we obtain

$$\frac{\sigma _\mathrm{i}^{}(z)}{\sigma _\mathrm{i}(z)}\frac{m}{H_\mathrm{f}z_\mathrm{i}}=\frac{1}{3(𝒩+1)}\frac{m}{H_\mathrm{f}}.$$
(23)
Therefore, the condition $`mH_\mathrm{f}`$ makes the curvaton field remain constant under general initial conditions, and we can write $`\sigma _\mathrm{f}\sigma _\mathrm{i}`$. Notice also that $`\dot{\sigma }_\mathrm{f}=0`$.

### III.2 Kinetic epoch

As we anticipated at the beginning of this section, the dynamics of the curvaton field after inflation are different to the scenarios considered before in the literature. This difference is mainly because the curvaton field coexists with the surviving inflaton field, whose energy density is by far the dominant one at the end of inflation. The consequences of this will be explored in this section and our immediate concern will be to work out the corresponding gravitational wave constraint.

According to our hypothesis, it is during the kinetic epoch that the curvaton field becomes effectively massive. This happens at a time when $`mH`$, from which we obtain (see Eq. (12))

$$\frac{m}{H_{\mathrm{kin}}}=\frac{a_{\mathrm{kin}}^3}{a_\mathrm{m}^3},$$
(24)
where ‘m’ labels the quantities at the time when the curvaton becomes massive. Up to this point, the curvaton field remained effectively massless and we can consider that $`\sigma _\mathrm{m}\sigma _\mathrm{i}`$.

In order to prevent a period of curvaton-driven inflation, the Universe must still be dominated by the scalar stiff-matter at this point. Hence, using Eqs. (12) and (24), we arrive to the restriction

$$\frac{m^2\sigma _\mathrm{i}^2}{2\rho _\varphi ^{(\mathrm{m})}}=\frac{4\pi }{3}\frac{m^2\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2m^2}1\sigma _\mathrm{i}^2\frac{3}{4\pi }m_{\mathrm{Pl}}^2.$$
(25)
Eq. (25) is sufficient to achieve subdominance at the time when the curvaton field becomes massive. But the curvaton energy should also be subdominant at the end of inflation, which implies

$$\frac{U_\mathrm{f}}{V_\mathrm{f}}=\frac{m^2\sigma _\mathrm{i}^2}{4\alpha ^2\lambda _\mathrm{b}}\frac{3m^2m_{\mathrm{Pl}}^2}{16\pi \alpha ^2\lambda _\mathrm{b}}=\alpha ^2\frac{m^2}{H_\mathrm{f}^2},$$
(26)
where we have used Eqs. (9) and (25). Thus, the curvaton mass should obey the stronger constraint

$$m\alpha ^1H_\mathrm{f}.$$
(27)
After the curvaton field becomes effectively massive, its energy decays as non-relativistic matter in the form

$$\rho _\sigma =\frac{m^2\sigma _\mathrm{i}^2}{2}\frac{a_\mathrm{m}^3}{a^3}.$$
(28)
### III.3 Curvaton domination

We calculate here the constraints to be applied on the curvaton model for the case in which the curvaton field comes to dominate the cosmic expansion, dealing first with the gravitational wave constraint. Since the curvaton energy redshifts as in Eq. (28), we find at the time of stiff scalar matter–curvaton matter equality $`(\rho _\sigma =\rho _\varphi )`$ that<sup>1</sup><sup>1</sup>1Note that the label $`a_{\mathrm{eq}}`$ has different meanings in Eqs. (16), (29) and (35).

$$\frac{\rho _\sigma }{\rho _\varphi }|_{a=a_{\mathrm{eq}}}=\frac{4\pi }{3}\frac{m^2\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2H_{\mathrm{kin}}^2}\frac{a_\mathrm{m}^3}{a_{\mathrm{kin}}^3}\frac{a_{\mathrm{eq}}^3}{a_{\mathrm{kin}}^3}=1,$$
(29)
Using Eqs. (24) and (29) in Eq. (15), we can calculate the energy density of gravitational waves at this time, which should obey the constraint

$$\frac{\rho _g}{\rho _\sigma }|_{a=a_{\mathrm{eq}}}=\frac{64}{3\pi }h_{\mathrm{GW}}^2\left(\frac{3}{4\pi }\frac{m_{\mathrm{Pl}}^2}{\sigma _\mathrm{i}^2}\frac{H_{\mathrm{kin}}}{m}\right)^{2/3}1.$$
(30)
To finish with, we derive the conditions to be applied on the decay parameter $`\mathrm{\Gamma }_\sigma `$. On the one hand, we require the curvaton field to decay before nucleosynthesis, and then $`H_{\mathrm{nucl}}=10^{40}m_{\mathrm{Pl}}<\mathrm{\Gamma }_\sigma `$. On the other hand, we also require the curvaton to decay after domination, and then $`\mathrm{\Gamma }_\sigma <H_{\mathrm{eq}}`$, where $`H_{\mathrm{eq}}`$ is found by using Eqs. (12), (24) and (29):

$$H_{\mathrm{eq}}=H_{\mathrm{kin}}\frac{a_{\mathrm{kin}}^3}{a_{\mathrm{eq}}^3}=\frac{4\pi }{3}\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}m.$$
(31)
Hence, the constraint upon the decay parameter is

$$10^{40}m_{\mathrm{Pl}}<\mathrm{\Gamma }_\sigma <\frac{4\pi }{3}\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}m.$$
(32)
### III.4 Curvaton decay before domination

We now turn our attention to the case in which the curvaton field decays before it dominates the cosmological expansion, but after it becomes massive so that we can use Eq. (28). As before, we calculate first the gravitational wave constraint. The curvaton field decays at a time when $`\mathrm{\Gamma }_\sigma =H`$ and then

$$\frac{\mathrm{\Gamma }_\sigma }{H_{\mathrm{kin}}}=\frac{a_{\mathrm{kin}}^3}{a_\mathrm{d}^3},$$
(33)
where ‘d’ labels the different quantities at the time of curvaton decay. The radiation produced from curvaton decay has an energy density

$$\rho _r^{(\sigma )}=\frac{m^2\sigma _\mathrm{i}^2}{2}\frac{a_\mathrm{m}^3}{a_\mathrm{d}^3}\frac{a_\mathrm{d}^4}{a^4}.$$
(34)
We will now assume that this radiation energy density is much larger than that produced by inflation, a necessary requirement if we are to maintain the gravitational waves under control. From Eq. (34), radiation equals the stiff scalar matter $`(\rho _\mathrm{r}^{(\sigma )}=\rho _\varphi )`$ at a time given by

$$\frac{\rho _\mathrm{r}^{(\sigma )}}{\rho _\varphi }|_{a=a_{\mathrm{eq}}}=\frac{4\pi }{3}\frac{m^2\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2H_{\mathrm{kin}}^2}\frac{a_\mathrm{m}^3}{a_{\mathrm{kin}}^3}\frac{a_{\mathrm{eq}}^2}{a_{\mathrm{kin}}^2}\frac{a_\mathrm{d}}{a_{\mathrm{kin}}}=1.$$
(35)
This equation defines the ratio $`a_{\mathrm{kin}}/a_{\mathrm{eq}}`$ to be used in Eq. (15). Hence, the constraint from gravitational waves now reads

$$\frac{\rho _g}{\rho _\mathrm{r}}|_{a=a_{\mathrm{eq}}}=\frac{16}{\pi ^2}h_{\mathrm{GW}}^2\sqrt[3]{\frac{\mathrm{\Gamma }_\sigma }{H_{\mathrm{kin}}}}\frac{m}{H_{\mathrm{kin}}}\frac{m_{\mathrm{Pl}}^2}{\sigma _\mathrm{i}^2}1,$$
(36)
where we have also used Eqs. (24) and (33).

Finally, we derive the new constraints on $`\mathrm{\Gamma }_\sigma `$. The curvaton field should decay after it becomes massive, so that $`\mathrm{\Gamma }_\sigma <m`$; and before it dominates the expansion of the Universe, $`\mathrm{\Gamma }_\sigma >H_{\mathrm{eq}}`$ (see Eq. (31)). Therefore, the constraints on the decay parameter are

$$\frac{4\pi }{3}\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}m<\mathrm{\Gamma }_\sigma <m.$$
(37)
### III.5 The surviving inflaton field

For completeness, we should add a small note here about the evolution of the inflaton field after curvaton domination and/or decay. It is well known that a steep exponential potential has a tracker behaviour where it evolves with the same equation of state as the dominant component scal . However, during the kination regime, the field is fast-rolling, and this typically takes it to very low energy densities, so that it can join the tracker only late in the cosmological evolution. Detailed analysis, following Ref. FJ , would be needed to determine its precise evolution, but it is unlikely to have any significant cosmological consequences unless the exponential form of the potential becomes modified at low energies.

### III.6 Curvaton fluctuations

Eqs. (30) and (37) are the constraints on the free parameters of the model coming from gravitational waves, but their actual form is not useful yet because of the term $`h_{\mathrm{GW}}^2`$. The latter is indeed related to other parameters of the model through the curvaton fluctuations, as we shall show in this section.

The curvaton fluctuations $`\delta \sigma _k`$ obey the linearized Klein–Gordon equation

$$\ddot{\delta \sigma _k}+3H\dot{\delta \sigma _k}+\left[\frac{k^2}{a^2}+m^2\right]\delta \sigma _k=0,$$
(38)
where $`k`$ is the comoving wave number.

During the time the fluctuations are inside the horizon, they obey the same differential equation as the inflaton fluctuations do, for which we conclude that they acquire the amplitude $`\delta \sigma _\mathrm{i}H_\mathrm{i}/2\pi `$. Once the fluctuations are out of the horizon, they obey the same differential equation as the unperturbed curvaton field does and then we expect that they remain constant during inflation, under quite general initial conditions.

Up to this point, the evolution of the curvaton fluctuations resembles that of previous scenarios already present in the literature. But, in a similar manner to our analysis of the homogeneous curvaton field in Sections III.3 and III.4, we should now take into account that the inflaton field survives the inflationary period, and study how the final spectrum of perturbations could be modified by this.

Following previous work, we can calculate the spectrum of the Bardeen parameter $`𝒫_\zeta `$, whose observed value is about $`2\times 10^9`$. When decay occurs after curvaton domination, the produced perturbation is lyth1 ; lyth2 ; nicola

$$𝒫_\zeta \frac{1}{9\pi ^2}\frac{H_\mathrm{i}^2}{\sigma _\mathrm{i}^2}.$$
(39)
This case is by far the simplest one. The spectrum of fluctuations is automatically gaussian since $`\sigma _\mathrm{i}^2H_\mathrm{i}^2/4\pi ^2`$, and is independent of $`\mathrm{\Gamma }_\sigma `$, a feature that will simplify the analysis of the parameter space. Moreover, even though the curvaton field coexists with the inflaton field up to this epoch, the spectrum of fluctuations is the same as in the standard scenario.

In the limit in which the curvaton decays when subdominant, the Bardeen parameter is given by lyth1 ; lyth2

$$𝒫_\zeta \frac{r_\mathrm{d}^2}{36\pi ^2}\frac{H_\mathrm{i}^2}{\sigma _\mathrm{i}^2},$$
(40)
where the normalization takes into account that the dominant component at that time is the scalar stiff matter. Here $`r_\mathrm{d}`$ is the ratio of curvaton energy density to stiff scalar matter at curvaton decay, which is found by using Eqs. (24) and (33) in

$$r_\mathrm{d}=\frac{\rho _\sigma }{\rho _\varphi }|_{a=a_\mathrm{d}}=\frac{4\pi }{3}\frac{m^2\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2H_{\mathrm{kin}}^2}\frac{a_\mathrm{m}^3}{a_{\mathrm{kin}}^3}\frac{a_\mathrm{d}^3}{a_{\mathrm{kin}}^3}=\frac{4\pi }{3}\frac{m}{\mathrm{\Gamma }_\sigma }\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}.$$
(41)
Contrary to the previous case, Eqs. (40) and (41) show explicitly the influence of the surviving inflaton field, and they are markedly different to the standard scenario in which the curvaton field coexists only with radiation.

## IV Model constraints

In this section, we use the different constraints on the curvaton field to study under which conditions it can be used together with steep inflation. The two cases that were worked out before will be treated separately.

### IV.1 Curvaton decay after domination

In this case, Eqs. (8), (9), and (39) allow us to fix $`\lambda _\mathrm{b}`$ in terms of $`\sigma _\mathrm{i}`$ and $`\alpha `$ as

$$\frac{\lambda _\mathrm{b}}{m_{\mathrm{Pl}}^4}=\frac{27\pi }{16}\frac{𝒫_\zeta }{(𝒩+1)^2\alpha ^4}\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2},$$
(42)
and the amplitude of primordial gravitational waves Eq. (14) can be written in terms of $`𝒫_\zeta `$, $`\sigma _\mathrm{i}`$ and $`\alpha `$ as

$$h_{\mathrm{GW}}^2=27\pi ^2\alpha ^2(𝒩+1)𝒫_\zeta \frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}.$$
(43)
Thus, Eqs. (30) and (27) are transformed into the following constraints on the curvaton mass:

$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $``$ $`10^3\sqrt{\pi ^3(𝒩+1)}\alpha \left(2.6\alpha ^221\right)𝒫_\zeta ^2{\displaystyle \frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}},`$ (44)
$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $``$ $`{\displaystyle \frac{3\pi 𝒫_\zeta ^{1/2}}{\alpha (𝒩+1)}}{\displaystyle \frac{\sigma _\mathrm{i}}{m_{Pl}}},`$ (45)
where we have used Eq. (13).

The initial value of the curvaton field $`\sigma _\mathrm{i}`$ is restricted by Eq. (25) and by another restriction coming from the assumption that the inflaton fluctuations are not relevant. The latter is obtained by using Eqs. (11) and (42), which give

$$\frac{\sigma _\mathrm{i}}{m_{\mathrm{Pl}}}<\frac{5}{\sqrt{72\pi }}\frac{1}{\alpha (𝒩+1)}\frac{A_\mathrm{s}}{\sqrt{𝒫_\zeta }}.$$
(46)
In other words, Eq. (46) just means that the brane tension has a lower value than that of Eq. (11).

Finally, Eq. (32) restricts the value of the decay parameter $`\mathrm{\Gamma }_\sigma `$, which can be transformed into another constraint upon $`m`$ and $`\sigma _\mathrm{i}`$ as

$$\frac{m}{m_{\mathrm{Pl}}}\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}\frac{3}{4\pi }\times 10^{40}.$$
(47)
The different conditions on the model of steep inflation $`+`$ curvaton field constrain the values of $`\lambda _\mathrm{b}`$, $`\sigma _\mathrm{i}`$ and $`m`$, but not of $`\alpha `$ and $`𝒩`$ which remain free parameters of the model. Therefore, we can draw the allowed region for the former in a plot $`m`$ versus $`\sigma _\mathrm{i}`$ by fixing the value of $`\alpha `$. Such a plot is shown in Fig. 1 for $`\alpha =15`$ and $`𝒩=70`$, in which we can see that there are viable models satisfying all the requirements.

However, the allowed region of the parameter space is reduced for larger values of the self-interaction of the inflaton field $`\alpha 1`$, as we can see from Eqs. (44) and (45). Thus, it is not always possible to find a successful model of curvaton field $`+`$ steep inflation, and this is mainly because of the gravitational wave constraint Eq. (36).

### IV.2 Curvaton decay before domination

For this case, it is convenient to use Eqs. (40) and (41) to write $`\sigma _\mathrm{i}`$ in terms of the other quantities as

$$\frac{\sigma _\mathrm{i}^2}{m_{\mathrm{Pl}}^2}=\frac{81}{4}\frac{𝒫_\zeta }{(𝒩+1)^2}\frac{H_{\mathrm{kin}}^2}{H_\mathrm{f}^2}\frac{m_{\mathrm{Pl}}^2}{H_{\mathrm{kin}}^2}\frac{\mathrm{\Gamma }_\sigma ^2}{m^2}.$$
(48)
On the other hand, the amplitude of primordial gravitational waves Eq. (14) is rewritten as

$$h_{\mathrm{GW}}^2=3\alpha ^2(𝒩+1)^3\frac{H_\mathrm{f}^2}{H_{\mathrm{kin}}^2}\frac{H_{\mathrm{kin}}^2}{m_{\mathrm{Pl}}^2},$$
(49)
and then the gravitational wave constraint Eq. (36) reads

$$\frac{64\alpha ^2}{27\pi ^2𝒫_\zeta }(𝒩+1)^5\frac{H_\mathrm{f}^4}{H_{\mathrm{kin}}^4}\frac{H_{\mathrm{kin}}^4}{m_{\mathrm{Pl}}^4}\left(\frac{H_{\mathrm{kin}}}{\mathrm{\Gamma }_\sigma }\right)^{5/3}\frac{m^3}{H_{\mathrm{kin}}^3}1.$$
(50)
In the same manner, Eq. (37) is now written as

$$\frac{27\pi 𝒫_\zeta }{(𝒩+1)^2}\frac{H_{\mathrm{kin}}^2}{H_\mathrm{f}^2}\frac{m_{\mathrm{Pl}}^2}{H_{\mathrm{kin}}^2}\frac{\mathrm{\Gamma }_\sigma ^2}{m}<\mathrm{\Gamma }_\sigma <m.$$
(51)
Notice that, from Eqs. (25) and (51) respectively, there are three similar constraints on the curvaton mass of the form

$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $`>`$ $`\left({\displaystyle \frac{\sqrt{27\pi 𝒫_\zeta }}{𝒩+1}}{\displaystyle \frac{H_{\mathrm{kin}}}{H_\mathrm{f}}}{\displaystyle \frac{m_{\mathrm{Pl}}}{H_{\mathrm{kin}}}}\right){\displaystyle \frac{\mathrm{\Gamma }_\sigma }{m_{\mathrm{Pl}}}},`$ (52)
$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $`>`$ $`\left({\displaystyle \frac{\sqrt{27\pi 𝒫_\zeta }}{𝒩+1}}{\displaystyle \frac{H_{\mathrm{kin}}}{H_\mathrm{f}}}{\displaystyle \frac{m_{\mathrm{Pl}}}{H_{\mathrm{kin}}}}\right)^2{\displaystyle \frac{\mathrm{\Gamma }_\sigma }{m_{\mathrm{Pl}}}},`$ (53)
$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $`>`$ $`{\displaystyle \frac{\mathrm{\Gamma }_\sigma }{m_{\mathrm{Pl}}}},`$ (54)
while Eqs. (27) and (50) give the following constraints, respectively,

$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $``$ $`{\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{H_f}{H_{\mathrm{kin}}}}{\displaystyle \frac{H_{\mathrm{kin}}}{m_{\mathrm{Pl}}}},`$ (55)
$`{\displaystyle \frac{m}{m_{\mathrm{Pl}}}}`$ $``$ $`\left[{\displaystyle \frac{27\pi ^2𝒫_\zeta }{64\alpha ^2(𝒩+1)^5}}{\displaystyle \frac{H_{\mathrm{kin}}^4}{H_\mathrm{f}^4}}\right]^{1/3}\left({\displaystyle \frac{m_{\mathrm{Pl}}^8}{H_{\mathrm{kin}}^8}}{\displaystyle \frac{\mathrm{\Gamma }_\sigma ^5}{m_{\mathrm{Pl}}^5}}\right)^{1/9}.`$ (56)
At this point, we should take into account the gaussianity condition upon the curvaton perturbations, since it is not automatically satisfied as in Section IV.1. The gaussianity condition implies that the perturbations must be small compared to the mean value of the curvaton field, and then $`\sigma _\mathrm{i}^2H_\mathrm{i}^2/4\pi `$. Using Eqs. (8), (13) and (48), we obtain

$$\frac{m}{m_{\mathrm{Pl}}}\frac{9\sqrt{\pi 𝒫_\zeta }}{(𝒩+1)^2}\frac{H_{\mathrm{kin}}^2}{H_\mathrm{f}^2}\frac{m_{\mathrm{Pl}}^2}{H_{\mathrm{kin}}^2}\frac{\mathrm{\Gamma }_\sigma }{m_{\mathrm{Pl}}}.$$
(57)
All these constraints can be fulfilled if $`m,\mathrm{\Gamma }_\sigma 0`$ independently of the value of $`H_{\mathrm{kin}}`$, but we know that this should not be possible since we still expect to recover the standard big bang scenario at nucleosynthesis. Hence, we will also impose the restriction $`\mathrm{\Gamma }_\sigma >10^{40}m_{\mathrm{Pl}}`$.

The final constraint has to do with the inflaton fluctuations, which should be negligible compared to the curvaton ones. In principle, this is a constraint upon $`\lambda _b`$ through Eq. (11), but Eqs. (9) and (13) allow us to transform it into a constraint on $`H_{\mathrm{kin}}`$ in the form

$$\frac{H_{\mathrm{kin}}}{m_{\mathrm{Pl}}}\frac{6.21\times 10^5}{\alpha ^3(𝒩+1)^2}\left(0.085\alpha ^20.688\right).$$
(58)
To unravel the parameter space allowed by the different constraints, we proceed as follows. We will take as free parameters the inflationary ones, $`\alpha ,𝒩`$ and $`H_{\mathrm{kin}}`$, but the latter will take only values permitted by Eq. (58). Once we fix the values of $`\alpha `$, $`𝒩`$ and $`H_{\mathrm{kin}}`$ by hand, we plot only the strongest constraint of Eqs. (52)–(54) (the other two will be automatically satisfied) and the constraints in Eqs. (55), (56) and (57). An example of this is shown in Fig. 2 for the values $`\alpha =15`$ and $`𝒩=70`$, where we also assumed that inflaton fluctuations contribute with as much as the $`10\%`$ of the COBE signal.

The different constraints permit some room for viable models, though the allowed region is smaller for smaller values of $`H_{\mathrm{kin}}`$, i.e. smaller values of $`\lambda _\mathrm{b}`$, and for larger values of $`\alpha `$. Then, we conclude that a realistic scenario in which the curvaton field decays before domination is very restricted by the different constraints.

## V Conclusions

The steep inflation model is of interest as it ends inflation in a novel manner, through reduction of inflation-sustaining braneworld corrections, and because it generates inflation from the steep exponential potentials often found in supergravity phenomenology. Unfortunately, the original steep inflation model appears unable to match observations, due to an excessive amplitude of short-scale gravitational waves and because its perturbations are not close enough to scale-invariance.

In this paper we have shown that curvaton reheating can alleviate these difficulties, by dramatically raising the reheat temperature (as compared to gravitational particle production reheating), while simultaneously suppressing the large-angle gravitational wave contribution. The model is subject to a complex network of constraints, but viable models do exist, provided the brane tension $`\lambda _\mathrm{b}`$ and the exponential slope $`\alpha `$ are not too large. We should stress that our analysis applies only to the case where the inflaton survives until late times, since an early decay of the inflaton field would lead to the standard scenario already studied in the literature lyth1 ; nicola ; lyth2 .

An interesting feature of the steep inflation model is that the inflaton field may survive to the present, and with suitable modification of the potential from its exponential form at low energies it could act as either dark matter lidsey ; FJ  or quintessence cope ; PV ; HL . This resembles the scenarios given in Ref. lyth2  concerning residual isocurvature matter perturbations. For instance, if cold dark matter is created before curvaton decay, then there is a maximal correlation between curvature and CDM-isocurvature perturbations which is at variance with current observations. We expect that the inflaton perturbations would be affected in the same way if the inflaton field is to be the dark matter at late times. The case of quintessence should not be troublesome, since the quintessence perturbations should decay upon horizon entry. However these situations merit more detailed investigation.

###### Acknowledgements.

L.A.U.-L. was supported by CONACyT, México under grant 010385, and A.R.L. in part by the Leverhulme Trust.


text:
# Third-harmonic exponent in three-dimensional 𝑁-vector models

## Abstract

We compute the crossover exponent associated with the spin-3 operator in three-dimensional $`O(N)`$ models. A six-loop field-theoretical calculation in the fixed-dimension approach gives $`\varphi _3=0.601(10)`$ for the experimentally relevant case $`N=2`$ (XY model). The corresponding exponent $`\beta _3=1.413(10)`$ is compared with the experimental estimates obtained in materials undergoing a normal-incommensurate structural transition and in liquid crystals at the smectic-A–hexatic-B phase transition, finding good agreement.

In nature many physical systems undergo phase transitions belonging to the universality classes of the O($`N`$) vector models. In particular, the XY model, corresponding to $`N=2`$, describes the $`\lambda `$-transition in 4He, (anti)ferromagnets with easy-plane anisotropy, density-wave systems, etc.; see Ref.  for a review. The critical exponents associated with the order parameter have been accurately measured both experimentally and theoretically . Moreover, in some XY systems it is also possible to measure experimentally the critical exponents associated with secondary order parameters. This is the case of liquid crystals , of normal-incommensurate transitions , and of graphite-intercalation compounds .

The most relevant exponent is the second-harmonic one that has been recently computed to high-precision using field-theoretical methods in Refs. . The $`ϵ`$-expansion gives $`\varphi _2=1.174(12)`$, while the fixed-dimension expansion gives $`\varphi _2=1.184(12)`$. The fourth-harmonic crossover exponent was reported in Refs. : $`\varphi _4=0.077(3)`$ ($`ϵ`$ expansion) and $`\varphi _4=0.069(5)`$ (fixed-dimension expansion). Here, we wish to determine the third-harmonic exponent by means of a six-loop perturbative calculation in the fixed-dimension approach, extending previous three-loop determinations . Such a calculation is also relevant for some crossover phenomena, in which the XY symmetry is reduced to that of the three-state Potts model, as it happens in cubic magnets in the presence of stress or of appropriate magnetic fields .

In the field-theoretical approach one starts from the usual $`\varphi ^4`$ Hamiltonian

$$=d^dx\left[\frac{1}{2}(_\mu \varphi )^2+\frac{1}{2}r\varphi ^2+\frac{1}{4!}u(\varphi ^2)^2\right],$$
(1)
where $`\varphi _a(x)`$ is an $`N`$-component real field. The XY model corresponds to $`N=2`$, but here we will keep $`N`$ generic. Secondary order parameters are associated with operators $`𝒪^{(l)}`$ that are polynomials of order $`l`$ in the fields and that transform irreducibly under the spin-$`l`$ representation of the $`O(N)`$ group. In particular, the third-harmonic operator is

$$𝒪_{abc}^{(3)}=\varphi _a\varphi _b\varphi _c\frac{\varphi ^2}{N+2}\left(\varphi _a\delta _{bc}+\varphi _b\delta _{ac}+\varphi _c\delta _{ab}\right).$$
(2)
We wish now to compute the crossover exponent $`\varphi _3`$ associated to $`𝒪_{abc}^{(3)}`$ and the corresponding exponents $`\beta _3`$ and $`\gamma _3`$ given by

$`\beta _3`$ $`=`$ $`2\alpha \varphi _3,`$ (3)
$`\gamma _3`$ $`=`$ $`2+\alpha +2\varphi _3.`$ (4)
The exponents $`\beta _3`$ and $`\gamma _3`$ describe respectively the critical (singular) behavior of the average $`𝒪^{(3)}(x)|t|^{\beta _3}`$ and of the susceptibility $`\chi _𝒪_x𝒪^{(3)}(0)𝒪^{(3)}(x)_c|t|^{\gamma _3}`$.

For this purpose we determine the renormalization function $`Z_3(g)`$ from the one-particle irreducible three-point function $`\mathrm{\Gamma }_3^{(3)}(0)`$ with an insertion of the operator $`𝒪_{abc}^{(3)}`$ at zero external momenta, i.e. we set

$$𝒪_{abc}^{(3)}\varphi _a\varphi _b\varphi _c^{1PI}=AZ_3^1(g),$$
(5)
where $`A`$ is a numerical coefficient that ensures that $`Z_3(0)=1`$, and $`g`$ is the four-point renormalized coupling. Then, we compute the renormalization-group function

$$\eta _3(g)\frac{\mathrm{ln}Z_3}{\mathrm{ln}m}|_u=\beta (g)\frac{d\mathrm{ln}Z_3}{dg},$$
(6)
and $`\eta _3=\eta _3(g^{})`$, where $`g^{}`$ is the fixed-point value of $`g`$. Finally, the renormalization-group scaling relation (valid in three dimensions)

$$\varphi _3=\left(\eta _3+\frac{3}{2}\frac{3}{2}\eta \right)\nu $$
(7)
allows us to determine $`\varphi _3`$.

We computed $`\mathrm{\Gamma }_3^{(3)}(0)`$ to six loops. The calculation is rather cumbersome since it requires the evaluation of a few thousand Feynman diagrams. We handled it with a symbolic manipulation program, which generates the diagrams and computes the symmetry and group factors of each of them. We used the numerical results compiled in Ref.  for the integrals associated with each diagram. We obtained

$`\eta _3(\overline{g})`$ $`=`$ $`{\displaystyle \frac{6}{N+8}}\overline{g}+{\displaystyle \frac{2(N+10)}{(N+8)^2}}\overline{g}^2{\displaystyle \frac{128.736+15.4900N0.650238N^2}{(N+8)^3}}\overline{g}^3`$ (12)
$`+{\displaystyle \frac{1148.68+191.005N+1.82163N^2+0.283028N^3}{(N+8)^4}}\overline{g}^4`$
$`{\displaystyle \frac{12606.9+2550.46N+64.4818N^22.34060N^30.152501N^4}{(N+8)^5}}\overline{g}^5`$
$`+{\displaystyle \frac{161373.+38736.8N+1874.23N^25.98451N^3+1.88168N^4+0.094179N^5}{(N+8)^6}}\overline{g}^6`$
$`+O(\overline{g}^7),`$
where, as usual, we have introduced the rescaled coupling $`\overline{g}`$ defined by

$$g=\frac{48\pi }{8+N}\overline{g}.$$
(13)
Here $`g`$ is the usual four-point renormalized coupling normalized so that $`g=u/m`$ ($`m`$ is the renormalized mass) at tree level.

Field-theoretical perturbative expansions are divergent, and thus, in order to obtain accurate results, an appropriate resummation is required. We use the method of Ref.  that takes into account the large-order behavior of the perturbative expansion, see, e.g., Ref. . Mean values and error bars are computed using the algorithm of Ref. .

Given the expansion of $`\eta _3(\overline{g})`$, we determine the perturbative expansion of $`\varphi _3(\overline{g})`$, $`\beta _3(\overline{g})`$, and $`\gamma _3(\overline{g})`$ using relations (4) and (7). Then, we resum the perturbative series and compute them at $`\overline{g}=\overline{g}^{}`$ . For $`N=2`$ we obtain $`\varphi _3=0.5963(21)`$, $`0.5968(2)`$, $`\beta _3=1.398(8)`$, $`1.405(3)`$, $`\gamma _3=0.800(7)`$, $`0.808(13)`$, where for each exponent we report the estimate obtained from the direct analysis and from the analysis of the series of the inverse. The two estimates obtained for each exponent agree within error bars, but, with the quoted errors, the scaling relations (4) are not well satisfied. For instance, using $`\nu =0.67155(27)`$ (Ref. ) and $`\beta _3=1.403(8)`$ we obtain $`\varphi _3=0.611(8)`$, while using the same value of $`\nu `$ and $`\gamma _3=0.803(13)`$ we have $`\varphi _3=0.606(6)`$. These two estimates are slightly higher than those obtained from the analysis of $`\varphi _3(g)`$ and $`1/\varphi _3(g)`$. Clearly, the errors are somewhat underestimated, a phenomenon that is probably connected with the nonanalyticity  of the renormalization-group functions at the fixed point $`\overline{g}^{}`$.

In order to obtain a conservative estimate, we have thus decided to take as estimate of $`\varphi _3`$ the weighted average of the direct estimates and of the estimates obtained using $`\beta _3`$ and $`\gamma _3`$ together with the scaling relations (4). The error is such to include all estimates. The other exponents are dealt with analogously . The final results for several values of $`N`$ are reported in Table I.

Let us compare these results with previous ones for $`N=2`$. Ref.  reports $`\beta _3=3\beta +3\nu x_3`$, where $`x_30.174`$ from the $`ϵ`$-expansion and $`x_30.276`$ from the fixed-dimension expansion. It follows $`\beta _31.397`$ and $`\beta _31.602`$ in the two cases. These results are reasonably close to ours. The exponent $`\beta _3`$ has been determined at the smectic-A–hexatic-B phase transition in liquid crystals, obtaining  $`\beta _34.8\beta 1.66`$, using $`\beta =0.3485(4)`$ (Ref. ). Such an exponent has also been measured in some materials exhibiting a structural normal-incommensurate phase transition. From the analysis of x-ray scattering data in erbium Ref.  obtains $`\beta _3=1.8(3)`$, while two different experiments in Rb2ZnCl4 give respectively $`\beta _3=1.80(5)`$ (Ref. ) and $`\beta _3=1.50(4)`$ (Ref. ). Keeping into account that the experimental errors seems to be underestimated, there is reasonable agreement with our results. The exponent $`\delta ^{}\beta /\varphi _3`$ was measured at the trigonal-to-pseudotetragonal transition in -stressed SrTiO3 obtaining  $`\delta ^{}=0.62(10)`$. Such an estimate is in good agreement with our prediction $`\delta ^{}=0.58(1)`$. Finally, we should mention that our results are also relevant for polymer physics. Indeed, we can derive from the estimates obtained for $`N=0`$ the partition-function exponent $`p`$ for nonuniform  star polymers with three arms: $`p=3(\gamma +\nu )/2+\varphi _3=3.06(1)`$, where we used $`\gamma =1.1575(6)`$ (Ref. ).