text:
# A search for radio emission from Galactic supersoft X-ray sources

## 1 Introduction

Supersoft sources are a separate class of X-ray objects. The most popular explanation of these sources is a white dwarf & sub-giant companion with a high accretion rate, 100–1 000 times greater than in cataclysmic variables (van den Heuvel et al. 1992). The large accretion rate creates steady hydrogen burning on the white dwarf surface causing X-ray emission.

Supersoft sources are difficult to detect in the Galaxy as the high column density in the Galactic plane absorbs most of the soft X-ray radiation. Consequently there have been many more detections at high galactic latitudes such as in the LMC, SMC and M31. As of 1999 (Greiner 2000), there were 57 supersoft sources: 10 in the Galaxy, 4 in the SMC, 8 in the LMC and 34 in M31 and 1 in NGC 55. Consequently, even though sources in the LMC and SMC are further away, they have been studied in greater detail due to their larger numbers. For example, in 1997 Fender, Southwell & Tzioumis (1998) searched for radio emission from non-Galactic supersoft sources despite them being prohibitively further away than their Galactic counterparts.

Of the persistent super-soft sources, three have been reported to have outflows, detected by emission lines in their optical and infrared spectra. What is the makeup of these sources, and are they the link between the low velocity jets seen in star-forming systems, and the superluminal jets seen in the microquasars?

### 1.1 Sources with outflow

Three supersoft sources have been observed with weak Doppler-shifted optical emission lines. The first object to be detected was RX J0513.9$``$6951, an LMC object which was discovered during an outburst in 1990 (Schaeidt et al. 1993), and optically identified by Pakull et al. (1993) and Cowley et al. (1993). This source has persistent optical emission with dips in the lightcurve every 100–200 days (Southwell et al. 1996). The evidence for jets in this source comes from strong He ii, H$`\beta `$ and H$`\alpha `$ lines with Doppler-shifted components at $`\pm `$ 4 000 km s-1 (Southwell et al. 1996). These lines are still observed during an optical dip, but with a reduced equivalent width, implying a reduction in the accretion rate.

The second supersoft jet source to be detected was RX J0019.8+2156 (Beuermann et al. 1995), a Galactic object at 2 kpc consisting of a 1 M WD and a 1.5 M donor star (Becker et al. 1998). The optical spectrum shows strong He ii emission and P Cygni absorption in the Balmer lines. This source displays slower red and blue-shifted lines with velocities of $`\pm `$ 815 km s-1, but these lines have a high FWHM of 400 km s-1 (Tomov et al. 1998). Quaintrell & Fender (1998) have also found these features in infrared spectra. Becker et al. (1998) conclude that the jets in this system come from an inclined cone of material blown off from an accretion disc.

The third, and last, supersoft source to be discovered with an outflow is the Galactic source RX J0925.7$``$4758. This was observed by Motch (1998) over two nights in 1997. Spectra showed prominent H$`\alpha `$ blue and red-shifted lines with a velocity of 5 800 km s-1. Motch (1998) modelled the emission from a cone with a large opening angle, at low inclination, similar to RX J0019+2156.

The only radio detection of a supersoft source was that of the symbiotic star AG Draconis (Torbett & Campbell 1987; Seaquist et al. 1993). Observations with the VLA showed two unresolved components separated by $``$ 1″, at fluxes of 200 and 400 Jy, shown in figure 1(a). These authors conclude that the emission is from optically thick thermal free-free emission, based on the majority of symbiotic systems. However, we cannot rule out optically thick emission in this paper.

## 2 Observations

Of the ten Galactic supersoft sources, four were observed by the VLA. These are RX J0019.8+2156, T Pyx, AG Draconis and V1974 Cygni. Three of these sources were observed with the VLA with the exception of AG Dra which was observed using MERLIN for the purpose of resolving any detail.

As we were primarily concerned with initial detections of point sources, a standard continuum setup was used for both the VLA and MERLIN observations. For the VLA we used a 100 MHz bandwidth split into two channels, and for MERLIN we used a 15 MHz bandwidth split into 15 channels. Individual channels were retained in order to reduce the amount of bandwidth smearing of sources offset from the phase-centre.

### 2.1 T Pyxis

T Pyx is a recurrent nova first discovered by Leavitt (1914) at a distance of 1.5 kpc (Shara 1997; Margon & Deutsch 1998). An arcsecond size optical nebula around the source, which shows a large amount of clumpy material on sub-arcsec scales , was detected from this object (Williams 1982; Duerbeck & Sitter 1987). This clumpy material would provide a rich medium for interaction with an outflow, such as the 1 100–1 400 km s-1 jets reported by Shahbaz et al. (1997) (however, this outflow probably originates in the surrounding nebular, and not in jets, see O’Brien & Cohen (1998) and Margon & Deutsch (1998)).

Observations were taken using the VLA in CnB configuration on 2000 March 3, 0415-0615 UT at 8.4 GHz and in D configuration on 2000 September 26 & 28, 1600-1710 & 1400-1630 UT respectively. Data were calibrated using the standard flux calibrator 3C 147, and phase calibrator PKS J0921$``$2618. No emission was observed from either the source, or any of the 7 surrounding sources from the SIMBAD database. We can put a $`1\sigma `$ RMS upper-limit of 22 Jy beam-1 at 4.8 GHz, and 19 Jy beam-1 at 8.4 GHz.

An observing log, together with beam sizes, RMS noise and source details for all sources is given in table 1.

### 2.2 V1974 Cygni ($``$ Nova Cygni 1992)

V1974 Cygni was extensively studied at all wavelengths following an outburst in 1992. Complementary VLA and MERLIN images show that the source initially flared then decayed, as a shell from the outburst expanded (Hjellming 1992; Eyres et al. 1996). They concluded that the emission from the nova was thermal and initially optically thick, then gradually became optically thin. This is expected from an expanding shell of material.

This source was observed by the VLA in CnB configuration on 2000 March 09 & 10 at 1500-1600 and 1330-1430 UT respectively at 8.4 GHz and in DnC configuration on 2000 June 27 at 0800-0945 UT at 4.8 GHz. Data were calibrated using the standard flux calibrator 3C 286 and phase-referenced using Ohio W 538. No emission was detected from Nova Cygni 1992. We can put a $`1\sigma `$ RMS upper limit to any emission of 24.6 Jy beam-1 at 4.8 GHz, and 19.8 Jy beam-1 at 8.4 GHz.

### 2.3 RX J0019.8+2156 ($``$ QR Andromedae)

RX J0019.8+2156 was observed by the VLA in CnB configuration on 2000 March 11 at 1810-2050 UT at 8.4 GHz and in DnC configuration on 2000 June 27 at 1300-1530 UT at 4.8 GHz. The flux calibrator used was 3C 48 and the data were phase referenced using PKS J0010+1724. No emission was detected from RX J0019.8+2156. We can therefore put a $`1\sigma `$ RMS upper limit to any emission of 22.7 Jy beam-1 at 4.8 GHz, and 19.8 Jy beam-1 at 8.4 GHz.

### 2.4 AG Draconis

To confirm the previous VLA detection, shown in figure 1(a), we observed AG Dra on 2000 March 18, 25 & 26 at 1010–1130, 0630–2140 and 0510–2250 UT respectively with MERLIN. Due to the superior resolving power of MERLIN over the VLA we intended to resolve the two unresolved components previously detected, and obtain some information about the source’s structure. The MERLIN image has a beam size of 48.9 $`\times `$ 68.5 mas at an angle of $`27.9\mathrm{°}`$, and a $`1\sigma `$ RMS noise of 103 Jy beam-1 at 4.994 GHz.

The core component is clearly detected in the MERLIN image and is resolved into two point sources at $`\alpha _{\mathrm{J2000}}`$ = 16h 01m 41$`\stackrel{s}{.}`$007, $`\delta _{\mathrm{J2000}}`$ = +66° 48′ 10$`\stackrel{}{.}`$13 with a flux of 570 Jy (designated as N1 in figure 1(b)), and $`\alpha _{\mathrm{J2000}}`$ = 16h 01m 41$`\stackrel{s}{.}`$046, $`\delta _{\mathrm{J2000}}`$ = +66° 48′ 10$`\stackrel{}{.}`$46 and a flux of 420 Jy (designated as N2). These fluxes imply a lower limit of $``$8000 K for the brightness temperature, consistent with the interpretation of the emission being optically thick free-free emission, however, we do not rule out optically thin emission in this paper. We note that the extension to the north of the component N1 is probably due to noise in the image.

The southwest component in the VLA image is not so clearly defined. If the source has not increased in flux between the epochs, and is a point source then it will be indistinguishable from the noise in the MERLIN image (shown in figure 1(c)). We do observe a number of possible point sources with fluxes around 300–350 Jy, but these are ambiguous compared to the noise. We have lowered the contours compared to figure 1(b) to show the noise in more detail.

### 2.5 New objects

All three VLA fields were imaged with a diameter of 10.2 arcmin in order to check for emission from known sources, and to detect new objects in the fields. However, interferometric images all suffer from bandwidth smearing or chromatic aberration. This is particularly important when the observed bandwidth is large and a wide field is imaged. The effects of these are to smear sources in a radial direction from the phase centre, and to reduce the observed peak flux. The consequence of this effect is that the signal to noise ratio increases with off-axis beam angle.

For an observation taken with a bandwidth of 50 MHz at 8.3 GHz, and with a circular beam 2 arcsec in diameter, a point source will suffer a 10 per cent reduction in flux at a distance of 3 arcmin 50 arcsec from the phase centre, and a 50 per cent reduction at a distance of around 11 arcmin, this situation becomes more severe with a smaller beam (see Bridle & Schwab (1989) for details).

We have uploaded to the CDS service for Astronomical Catalogues111http://cdsweb.u-strasbg.fr/cats/Cats.htx a table which lists the sources taken from the Simbad catalogue together with our limits to the radio emission from the source based on the bandwidth smeared $`5\sigma `$ RMS noise. While none of the known sources were detected, we observed 17 new sources with a flux greater than $`5\sigma `$. An uploaded table as above lists these new sources, together with a Gaussian-fitted flux based on the image beam, and true flux due to the bandwidth smearing in the image, (at the bands $`C`$ and $`X`$ given in table 1). A spectral index is also calculated from the bandwidth smeared fluxes and follows the convention of $`S\nu ^\alpha `$.

We can therefore estimate that the density of sources with flux $`>`$ 100 Jy is about 200 sources per square degree at 5 GHz, broadly consistent with limits from the deep radio survey by MERLIN of the Hubble Deep Field (Muxlow et al. 1999). The spectral indices of these sources are generally negative or consistent with a negative value, suggesting that most are non-thermal extragalactic objects, however caution has to been drawn to this conclusion since data at the two frequencies is not simultaneous, and there is a difference in angular resolution between the two frequencies.

## 3 Wind emission

In this section we describe simple analytical models for wind emission and apply them to our sample. Initially we are not interested in where the wind comes from or what it consists of, rather that a generic wind with a particular composition, being ejected at a particular rate from a source at a given distance from us, would produce radio emission. The optically thin flux from a source with these parameters can be calculated and is given by

$$S=23.2\left(\frac{\dot{M}}{\mu v_{\mathrm{km}/\mathrm{s}}}\right)^{4/3}\left(\nu \gamma gZ^2\right)^{2/3}D^2\mathrm{Jy}$$
(1)
where $`\dot{M}`$ is in M yr-1, $`\mu `$ is the mean atomic weight, $`v_{\mathrm{km}/\mathrm{s}}`$ is the wind velocity in km s-1, $`\nu `$ is in Hz, $`\gamma `$ is the ratio of electron to ion number density, $`g`$ is the Gaunt factor, $`Z`$ is the ionic charge and $`D`$ is the distance in kpc (Wright & Barlow 1975). The Gaunt factor here has a slight dependence on temperature and composition and if one does not consider extremely hot winds, then it can be approximated by

$$g=9.77+1.27\mathrm{log}_{10}\left(T^{3/2}/\nu Z\right)$$
(2)
(Leitherer & Robert 1991).

We can now investigate in more detail from where the wind originates, what is its expected composition and what is the expected radio emission from such a region.

### 3.1 Generic secondary wind

To create sustained mass transfer between the binary components either the secondary222In this paper we use the convention that the WD is designated the primary object and the companion, or donor, is designated the secondary. Subscripts 1 and 2 are used in mass and radius relationships to distinguish between these two objects is undergoing Roche lobe overflow, the wind is providing the accretion via Bondi-Hoyle accretion, or both. Which scenario occurs depends on the nature of the secondary, and the binary parameters.

The symbiotic system AG Dra with a $`1.5`$ M K-giant, a $`0.5`$ M white dwarf and a period of 544 days does not normally fill its Roche lobe (the stellar radius of the K-giant is $`30`$ R and the Roche lobe size is $`170`$ R). Since the system must be accreting then the logical conclusion is that the secondary is on the Asymptotic Giant Branch (AGB) and feeding the primary by a slow, dense wind. Out of the four sources considered here, this is the most likely to be feeding its environment directly by a wind. Calculating the expected flux from this wind uses a temperature, velocity, outflow and distance of $`T`$ = 15 000 K, $`v`$ = 30 km s-1, $`\dot{M}`$ = 1.2–1.7 $`\times 10^8`$ M yr-1 and $`D`$ = 1.56–1.81 kpc from Tomova & Tomov (1999), compositions of $`\mu =1.4`$, $`\gamma =1`$ and $`Z=1`$ from Nussbaumer & Vogel (1987) and Spergel, Giuliani & Knapp (1983). These parameters yield a flux at 5 GHz of $`S_C`$ = 12–26 Jy.

In contrast to AG Dra, the three other sources either have a small secondary companion and therefore Roche lobe overflow would be providing the mass-transfer, or as in the case for V1974 Cyg, radio emission would have faded and dissipated from the time of the last outburst in 1992.

### 3.2 Evaporated secondary wind

In the previous section it was shown that only the AGB binary AG Dra would have a strong enough secondary wind to emit significant radio emission. However, if the WD is radiating sufficient soft X-rays then a corona will be formed around the secondary star, which causes a strong stellar wind to be evaporated (Basko & Sunyaev 1973). This type of interaction has been investigated in the context of the supersoft sources by van Teeseling & King (1998, hereafter vTK98) and they show that the mass loss rate in the secondary due to irradiation is

$$\dot{M}_73\varphi \frac{r_2}{a_9}\left(m_2L_{30}\eta _s\right)^{1/2}\mathrm{M}_{}\mathrm{yr}^1$$
(3)
where $`\varphi `$ is an irradiation efficiency factor, $`r_2`$ is the radius of the secondary in solar units, $`a_9`$ is orbital separation in units of $`10^9`$ m, $`m_2`$ is the mass of the secondary in solar units, $`L_{30}`$ is the irradiation luminosity in units of $`10^{30}`$ J s-1 ($`10^{37}`$ erg s-1), $`\eta _s`$ is an ionisation efficiency parameter, and $`\dot{M}_7`$ is given in units of $`10^7`$. It is reasonable, therefore, to assume that since the supersoft sources produce soft X-rays during quiescence from their sustained nuclear burning, that the majority (if not all) of the supersoft sources have mass losses which are wind dominated.

To obtain an estimate for the mass loss from the secondary due to irradiation from the WD, some assumptions to the geometry of the system need to be made. In equation 3, the variable $`\varphi `$ is an efficiency parameter representing the fraction of the companion’s face which is being radiated and the fraction of the wind mass escaping the system. Commentary on the importance of $`\varphi `$ is given in Knigge, King & Patterson (2000), and a value of 1 is expected. The other unknown parameter in equation 3 is the efficiency of the incident radiation in driving a wind, $`\eta _s`$. For a wide range of temperatures, including the expected temperature of the irradiation (few $`\times 10^5`$ K), $`\eta _s1`$. We can therefore calculate an estimation for the mass-loss in the wind from the secondary, which is given in table 2.

Together with the mass-loss rate and equation 1, the flux from the wind can be found. Unfortunately, the expected fluxes are not as precise as one hopes as a large number of important parameters such as composition ($`\mu `$, $`\gamma `$ and $`Z`$), wind temperature and velocity ($`T`$ and $`v_{\mathrm{}}`$) are unknown.

van Teeseling & King (1998) consider He ii to be the main absorber of the WD radiation, which has implications to the composition and temperature of the resulting wind. The temperature of the wind will range from the ionisation temperature of He ii to the temperature of the incident radiation, so we use a wind temperature of $`T5\times 10^5`$ in this paper. The composition of a purely He ii wind would give values for the mean atomic weight $`(\mu =4)`$, the electron to ion number density ratio $`(\gamma =2)`$ and the ionic charge $`(Z=2)`$. It would be surprising if these were realistic values for the wind, but they are useful approximations and only enter the wind-flux equation (1) as $`S\mu ^{4/3}\gamma ^{2/3}Z^{4/3}`$ and so, small changes would not have a large effect on the estimation. The final variable to estimate is the wind velocity, $`v_{\mathrm{}}`$. vTK98 approximate the velocity to be 0.3 times the isothermal sound speed. For sound speeds of the order $`100`$ km s-1 this would give approximate wind velocities of $`30`$ km s-1. Using these approximations to the wind gives flux estimates at 5 GHz ($`S_\mathrm{C}`$) and 8.3 GHz ($`S_\mathrm{X}`$) shown in Table 2.

### 3.3 White dwarf wind

We have tried to estimate the radio flux from a white dwarf wind, but the uncertainties are too large to confidently quote numbers. Both Bondi-Hoyle and Roche lobe overflow accretion cannot be reliably used to predict the mass-transfer rate, and therefore the amount of excess material that the white dwarf can process is also undetermined.

### 3.4 Observational comments to wind emission

From our MERLIN and VLA observations we can provide the flux for AG Draconis and upper limits to the fluxes of the other three supersoft sources in this paper (presented in table 2). Given the spread in the observable system parameters, the expected flux for all three non-detected sources is at, or below the 5 $`\sigma `$ cutoff for detections. The surprise in this analysis is with AG Dra where we detect with MERLIN two unresolved point sources with a combined flux of $``$ 1 mJy. This is much larger than both the wind flux and an evaporated wind at that frequency of around 10-20 Jy and so we can conclude either the assumptions used in these models are unrealistic, or a different emission mechanism is producing the observed radio flux.

It could be significant that the only detected fluxes originate from a K giant secondary rather than F or M dwarfs. The geometry and accretion in AG Dra is not expected to be from Roche lobe overflow, but rather from the wind of the AGB secondary.

## 4 Conclusions

We have searched for quiescent, persistent radio emission from the northern hemisphere supersoft X-ray sources, and taken a high resolution image of the one known source with emission.

We have improved the radio positions, source size and flux from the persistent emitter AG Dra with MERLIN. The core is resolved at the milliarcsec scale into two components with a combined flux of $``$ 1000 Jy. It is possible that the core detected in the VLA image has significant nebulosity and is resolved out by the higher MERLIN resolution. A possible south component detected by the VLA was unconfirmed, although this could be due to it having a flux lower than the noise in the MERLIN image, or also being resolved out.

No new emission was detected from RX J0019.8+2156, T Pyx and V1974 Cygni down to a 1$`\sigma `$ RMS noise level of around 20 Jy beam-1. No emission was detected from any previously known Galactic or extragalactic source in the 3, 104 square arcsec fields imaged, and we place upper-limits to radio emission from these sources.

We have investigated possible causes of radio emission from a wind environment, both directly from the secondary star, and also as a consequence of the high X-ray luminosity from the WD.

A total of 17 new point-sources were imaged with fluxes between 100 and 1500 Jy giving an estimate to the density of sources with a flux $`>`$ 100 Jy of 200 per square degree, at 5 GHz, in this region of the sky.

## 5 acknowledgments

RNO wishes to thank the hospitality of the MERLIN national facility at the Jodrell Bank Observatory, especially Tom Muxlow. SC acknowledges support from grant F/00-180/A from the Leverhulme Trust.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. MERLIN is a national facility operated by the University of Manchester on behalf of PPARC. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France


format: text

text:
# Accurate kinetic energy evaluation in electronic structure calculations with localised functions on real space grids

## 1 Introduction

Density functional theory (DFT) combined with the pseudopotential method has been established as an important theoretical tool for studying a wide range of problems in condensed matter physics . However, the computational cost of performing a total-energy calculation on a system scales asymptotically as the cube of the system size. Consequently, plane-wave pseudopotential DFT can only be used to study systems of up to about one hundred atoms on a single workstation and up to a few hundred atoms on parallel supercomputers. As a result there has been considerable recent effort in the development of methods whose computational cost scales linearly with system size .

A common feature of many of the linear-scaling strategies is the expansion of the single-particle density matrix in terms of a set of localised functions. We refer to these functions as ‘support functions’ . A support function is required to be non-zero only within a spherical region, which we refer to as a ‘support region’, centred on an atomic position. Here we consider a representation of the support functions in terms of a regular real space grid, which constitutes our basis set. If the set of support functions is $`\{\varphi _\alpha \}`$, the single-particle density matrix is expressed as

$$\rho (𝐫^{},𝐫)=\underset{\alpha ,\beta }{}\varphi _\alpha (𝐫^{})K^{\alpha \beta }\varphi _\beta ^{}(𝐫),$$
(1)
where $`K^{\alpha \beta }`$ is the matrix representation of the density matrix in terms of the duals of the support functions. In general the support function set is not orthonormal.

Real space methods have the advantage that they provide a clear spatial partitioning of all quantities encountered in a density functional calculation, a property that is ideal for code parallelisation. As a result, this approach has gained popularity in recent years and a number of such density functional calculations have been reported by different authors . These approaches use finite difference (FD) methods  for the calculation of the kinetic energy. In terms of the support functions the kinetic energy is

$$E_T[\rho ]=\underset{\alpha ,\beta }{}T_{\alpha \beta }K^{\beta \alpha },$$
(2)
where $`T_{\alpha \beta }`$ denotes kinetic energy matrix elements between support functions, given in Hartree atomic units by

$$T_{\alpha \beta }=\frac{1}{2}\varphi _\alpha ^{}(𝐫)^2\varphi _\beta (𝐫)d𝐫.$$
(3)
The evaluation of the kinetic energy matrix elements requires the action of the Laplacian operator on the support functions. Here we will show that in the case of localised support functions, fast Fourier transform (FFT) methods can be adapted for the application of the Laplacian, providing an algorithm with essentially the same computational cost as FD but with higher accuracy and also ready applicability to any grid symmetry.

In the following two sections we present the FD method and our new FFT-based method and compare them both in theory and in practice.

## 2 Theory

For functions represented as values on a regular grid, integrals like the one of equation (3) can be calculated, or rather approximated to increasing accuracy, by a sum over grid points, as long as the value of the integrand is known at every grid point:

$$T_{\alpha \beta }\frac{1}{2}w\underset{𝐫_i}{}\varphi _\alpha ^{}(𝐫_i)\widehat{T}\varphi _\beta (𝐫_i),$$
(4)
where $`\widehat{T}`$ is the Laplacian operator in the discrete representation, $`w`$ is the volume per grid point, and the sum formally goes over all the grid points in the simulation cell.

### 2.1 Finite Differences

The most straightforward approach to the evaluation of the Laplacian operator applied to a function at every grid point is to approximate the second derivative by finite differences of increasing order of accuracy . For example, the $`^2\varphi /x^2`$ part of the Laplacian on a grid of orthorhombic symmetry is

$$\frac{^2\varphi }{x^2}(x_i,y_j,z_k)\frac{1}{h_x^2}\underset{n=A/2}{\overset{A/2}{}}C_n^{(A)}\varphi (x_i+nh_x,y_j,z_k)+O(h_x^A),$$
(5)
where $`h_x`$ is the grid spacing in the $`x`$-direction, $`A`$ is the order of accuracy and is an even integer, and the weights $`C_n^{(A)}`$ are even with respect to $`n`$, i.e. $`C_n^{(A)}=C_n^{(A)}`$. This equation is exact when $`\varphi `$ is a polynomial of degree less than or equal to $`A`$. The leading contribution to the error is of order $`h_x^A`$. The full Laplacian operator for a single grid point in three dimensions consists of a sum of $`(3A+1)`$ terms.

In principle, for well behaved functions, the second order form of equation (5) should converge to the exact Laplacian as $`h0`$. Therefore to increase the accuracy of a calculation one would need to proceed to smaller grid spacings. However, in most cases of interest, this is computationally undesirable and instead, formulae of increasing order are used to improve the accuracy at an affordable cost . Chelikowsky et al. , in their finite difference pseudopotential method, have tested the finite difference expression for up to $`A=18`$ on calculations of a variety of diatomic molecules and have suggested $`A=12`$ as the most appropriate for their purpose, as the higher orders did not provide any significant improvement.

Alternative discretisations of the Laplacian operator are possible, such as the Mehrstellen discretisation of Briggs et al. . This is a fourth order discretisation that includes off-diagonal terms, but only nearest neighbours to the point of interest. It is more costly to compute than the standard fourth order formula of equation (5) and it is still not clear whether its fourth order is sufficient. One may also use FD methods on a grid with variable spatial resolution, such as that of Modine et al.  which is denser near the ionic positions. Such a scheme, however, has the added overhead of a transformation of the Laplacian from Cartesian to curvilinear coordinates. In this paper we use only the FD scheme of equation (5).

The FD approach has desirable properties, both in terms of computational scaling and parallelisation. The Laplacian in the FD representation is a near-local operator, becoming more delocalised with increasing order. Therefore, the cost of applying it to $`N`$ grid points is strictly linear (compared to $`N\mathrm{log}N`$ for Fourier transform methods). Also, as a result of its near-locality, ideal load balancing can be achieved in parallel implementations by partitioning the real space grid into subregions of equal size and distributing them amongst processing elements (PEs) while requiring little communication for applying the Laplacian at the bordering points of the subregions.

If $`N_s`$ represents the size of the system, then the number of support functions will be proportional to $`N_s`$ and so will the total number of grid points in the simulation cell, resulting in a total computational cost proportional to $`N_s^2`$ for the application of the Laplacian on all support functions. More favourable scaling can be achieved by predicting the region in space whithin which the values of a particular function will be of significant magnitude and operating only on this region . Linear-scaling can be achieved by strictly restricting from the outset the support functions to spherical regions centred on atoms . In this case, the cost is $`qN_s`$ with $`q`$ being the cost of applying the Laplacian on the points of a spherical region, which is constant with system size.

FD methods nevertheless have disadvantages that do not appear in the plane-wave formalism. Firstly, there is no *a priori* way of knowing whether a particular order of FD approximation will be sufficient to represent a particular support function accurately. In addition, while plane-wave methods can handle different symmetry groups trivially through the reciprocal lattice vectors of the simulation cell, real space implementations need to consider every symmetry separately and require considerable modifications to the code and higher computational cost. Briggs et al.  have demonstrated this difficulty by performing calculations with hexagonal grids while most common applications of real space methods in the literature are limited to grids of cubic or orthorhombic symmetry .

The computational cost for the calculation of the Laplacian of a single support function with the FD method scales as $`(3A+1)(1+A/D)^3N_{\mathrm{reg}}`$ where $`N_{\mathrm{reg}}`$ is the number of grid points within the support region, and $`D`$ is the number of grid points along the support region diameter and is proportional to $`N_{\mathrm{reg}}^{1/3}`$. This estimate of cost includes all the nonzero values of the Laplacian, which in general occur not only at the grid points inside the support region but also at points outside, up to a distance of $`A/2`$ points from the region’s boundary. It is important to include the contribution to the Laplacian from outside the support region in the sum of equation (4) in order to obtain the best possible accuracy for a given order $`A`$ and also to ensure the Hermiticity of the discretised representation of the Laplacian, $`\widehat{T}`$, and hence of the kinetic energy matrix elements $`T_{\alpha \beta }`$.

### 2.2 Localised discrete Fourier transform

We now present an accurate, linear-scaling method for calculating the kinetic energy matrix elements $`T_{\alpha \beta }`$ of equation (3). We use a mixed space Fourier transform approach that is applicable to any Bravais lattice symmetry. Fourier transformation is a natural method to adopt for this task since in a total-energy calculation one computes other terms, such as the electron density and the Hartree energy, using reciprocal space techniques. This implicitly defines the basis set that we use to be plane-waves and for consistency we should calculate the kinetic energy using the same basis set, i.e. using Fourier transform methods. Thus we calculate the $`^2\varphi `$ term in reciprocal space, where the Laplacian operator is easy to apply, then transform the result back to real space and obtain the matrix elements $`T_{\alpha \beta }`$ by summation over grid points (4). One way to achieve this would be to perform a discrete FFT on each support function $`\varphi `$, using the periodicity of the entire simulation cell. However, unlike the FD algorithm, the FFT is not a local operation and the cost of applying the Laplacian to all the support functions in this way would be proportional to $`N_s^2\mathrm{log}N_s`$, which clearly does not scale linearly with system size.

It is possible to overcome this undesirable scaling without compromising accuracy by performing the FFT over a restricted region of the simulation cell, which we call the ‘FFT box’ (figure 1). Before defining the FFT box, there are two points that should be noted. Firstly, the operator $`\widehat{T}`$ must be Hermitian. This will ensure that the kinetic energy matrix elements, $`T_{\alpha \beta }`$ are Hermitian, and hence the eigenvalues real. Secondly, when calculating two matrix elements such as $`T_{\alpha \beta }`$ and $`T_{\gamma \beta }`$, we require the quantity $`\widehat{T}\varphi _\beta `$ in both cases. To be consistent, our method for calculating the matrix elements must be such that $`\widehat{T}\varphi _\beta `$ is the same in both cases, i.e. we require $`\widehat{T}\varphi _\beta `$ to have a unique and consistent representation throughout the calculation. It is important that both these conditions are satisfied when it comes to optimisation of the support functions during a total-energy calculation, and we shall return to this point later.

In order to fulfil the above requirements, it can be seen that for a given calculation the FFT box must be universal in shape and dimensions. As a result, it must be large enough to enclose any pair of overlapping support functions within the simulation cell. To define a suitable FFT box, we first consider a box with the same unit lattice vectors as the simulation cell, but of dimensions such that it exactly circumscribes the largest support region present in the simulation cell. We then define a box that is commensurate with this, but with sides that are twice as long (and hence a volume eight times as large). This we define to be the FFT box. It is clear that this FFT box is large enough to enclose any pair of support functions exhibiting any degree of overlap.

To calculate a particular matrix element $`T_{\alpha \beta }`$ for two overlapping support functions $`\varphi _\alpha `$ and $`\varphi _\beta `$, we imagine them as being enclosed within the FFT box defined above and we treat this region of real space as a miniature simulation cell. We Fourier transform $`\varphi _\beta `$ using the periodicity of the FFT box and apply the Laplacian at each reciprocal lattice point using standard plane-wave techniques . It is then a simple matter of using one more FFT to back-transform $`^2\varphi _\beta `$ to real space and subsequently calculate $`T_{\alpha \beta }`$ by summation over the grid points of the FFT box, according to equation (4).

The result obtained by this process is equivalent to performing a Fourier transform of $`\varphi _\beta `$ over the whole simulation cell, applying the Laplacian and then interpolating to a coarse, but still regular, reciprocal space grid with only $`N_{box}`$ points, $`N_{box}`$ being the number of grid points in the FFT box, before back-transforming to real space. This coarse sampling in reciprocal space has a negligible influence on the result because each support function is strictly localised in real space and therefore smooth in reciprocal space.

It is worth noting the implicit approximation that we make in calculating the kinetic energy in the way prescribed above. In general, $`^2\varphi _\beta `$ is nonzero outside the support region of $`\varphi _\beta `$ itself, and it is essential to take this into account in the calculations. By construction, we neglect contributions to the kinetic energy from support functions whose support regions do not overlap as we expect them to be negligibly small. This approximation may be controlled via a single parameter, the FFT box size, with respect to which the calculation may be converged if necessary. The same approximation is of course present in the FD method as well.

We expect certain advantages to the FFT box algorithm over FD based methods. Firstly, the FFT box method should be more accurate than any FD scheme since it takes into account information from every single point of the support function and not only locally. However, it is still perfectly local as far as parallelisation is concerned since we only deal with the points within a single FFT box each time, and this constitutes a very small region of the simulation cell. The parallelisation strategy in this case would still consist of partitioning the real space grid of the simulation cell into subregions of equal size and distributing them amongst PEs. Then, FFTs local to each PE are performed on FFT boxes enclosing pairs of overlapping support regions belonging completely to the simulation cell subregion of the given PE. For pairs of overlapping support regions containing grid points common to the subregion of more than one PE, the pair would have to be attributed to one PE and copied as a whole to it for the local FFT to proceed. This would involve some communication overhead, as in the FD case for pairs of overlapping support functions with points in more than one subregion. Another important advantage of the FFT box method is that it is applicable, without any modification, to regular grids of any Bravais lattice symmetry. This is not true of FD methods.

The number of grid points in a cubic FFT box is $`N_{box}`$ (which is related to $`N_{\mathrm{reg}}`$ by $`N_{box}=8\times 6N_{\mathrm{reg}}/\pi 15.3N_{\mathrm{reg}}`$). Therefore the computational cost of applying the FFT method to a single support function in such an FFT box is $`2N_{box}\mathrm{log}N_{box}`$, and thus for all support functions the cost is proportional to $`2N_sN_{box}\mathrm{log}N_{box}`$, where $`N_{box}`$ is independent of $`N_s`$. In other words the cost scales linearly with the number of atoms in the system.

## 3 Tests and discussion

We have performed tests of the FD and FFT box methods for calculating the kinetic energy of localised functions. Choosing a particular type of support function $`\varphi `$ with spherical symmetry, placing one at $`𝐑_\alpha `$ and another at $`𝐑_\beta `$, we rewrite the integral of equation (3) as

$$T(|𝐑_\alpha 𝐑_\beta |)=\frac{1}{2}\varphi ^{}(𝐫𝐑_\alpha )^2\varphi (𝐫𝐑_\beta )d𝐫.$$
(6)
For our first test we calculate the following quantity as a function of the distance $`d`$ between the centres $`𝐑_\alpha `$ and $`𝐑_\beta `$

$$\eta _1(d)=T_{ap}(d)T_{ex}(d),$$
(7)
where $`T_{ex}(d)`$ is the exact value of the integral in the continuous representation of the support functions and $`T_{ap}(d)`$ is its approximation on the real space grid, either by FD or the FFT box method. We chose $`\varphi (𝐫)`$ to be a 2s valence pseudo-orbital for a carbon atom, generated using an atomic norm-conserving carbon pseudopotential  within the local density approximation. The pseudo-orbital is confined in a spherical region of radius $`6.0a_0`$, and vanishes exactly at the region boundary . It is initially generated as a linear combination of spherical Bessel functions, which are the energy eigenfunctions of a free electron inside a spherical box. Our functions are limited up to an energy of 800 eV, resulting in a combination of fourteen Bessel functions. The formula for calculating kinetic energy integrals between Bessel functions is known  and we used it to obtain $`T_{ex}(d)`$ for our valence pseudo-orbital. We then calculated $`\eta _1(d)`$ with a grid spacing of $`0.4a_0`$ (corresponding to a plane-wave cut-off of 839 eV) in an orthorhombic simulation cell, as we are restricted to do so by the FD method. With these parameters $`N_{box}`$ is $`60^3`$, and hence it is trivial to perform the FFT of one support function on a single node. $`\eta _1(d)`$ is plotted for the FFT box method and for various orders of the FD method in the top graph of figure 2.

It can be seen that low order FD methods are inaccurate as compared to the FFT box method, and only when order 28 FD is used does the accuracy approach that of the FFT box method. The A = 12 FD scheme, the highest order that has been used in practice for calculations , gives an error of $`3.97\times 10^5\text{Hartree}`$ at $`d=0`$ as compared to $`1.027\times 10^5\text{Hartree}`$ for the FFT box method. The feature that occurs in the top graph of figure 2, between $`d=5a_0`$ and $`d=7a_0`$, is an artefact of the behaviour of our pseudo-orbitals at the support region boundaries where they vanish exactly, but with a finite first derivative. This causes an enhanced error in all the methods when the edge of one support function falls on the centre of another.

The error in the FFT box method is small, yet non-zero, and we attribute this to the inherent discretisation error associated with representing functions that are not bandwidth limited on a discrete real space grid. Convergence to the exact result is observed as the grid spacing is reduced, as expected.

As our next comparison of the FFT box and FD methods we used the same pseudo-orbitals as before, but considered the quantity

$$\eta _2(d)=T_{ap}(d)T_{PW}(d)$$
(8)
as the measure of the error, where $`T_{PW}`$ is the result obtained by Fourier transforming the support functions using the periodicity of the entire simulation cell. One may think of $`T_{PW}`$ as being the result that would be obtained from a plane-wave code: the support functions may be considered as generalised Wannier functions. Calculating the kinetic energy integrals by performing a discrete Fourier transform on the support functions over the entire simulation cell (an $`O(N_s^2\mathrm{log}N_s)`$ process for all support functions) is equivalent to summing the contributions to the kinetic energy from all of the plane-waves up to the cut-off energy determined by the grid spacing. Thus our FFT box method can be viewed as equivalent to a plane-wave method that uses a contracted basis set (i.e. a coarse sampling in reciprocal space). In some ways $`\eta _2(d)`$ is a better measure of the relative accuracy of the FD and FFT box methods as our goal is to converge to the ‘exact’ result as would be obtained using a plane-wave basis set over the entire simulation cell. $`\eta _2(d)`$ is plotted in the bottom graph of figure 2.

$`T_{PW}`$ was computed using a cell that contained 256 grid points in each dimension. Increasing the cell size further had no effect on $`T_{PW}`$ up to the eleventh decimal place ($`10^{11}\text{Hartree}`$). The plots show that the FFT box method performs significantly better than all orders of FD that were tested. For example at $`d=0`$ the error for A = 28 FD is $`3.49\times 10^6\text{Hartree}`$ as compared to $`1.09\times 10^9\text{Hartree}`$ for the FFT box method. The fact that the FFT box error is so small shows that coarse sampling in reciprocal space has little effect on accuracy, as one would expect for functions localised in real space.

Our implementation can produce similar FFT box results to the above in regular grids of arbitrary symmetry (non-orthogonal lattice vectors) as long as we include roughly the same number of grid points in the support region sphere. As we described earlier the application of the FD method to grids without orthorhombic symmetry is not straightforward.

Furthermore, in our implementation the kinetic energy matrix elements $`T_{\alpha \beta }`$ for both the FFT box method and the FD method (of any order) are Hermitian to machine precision. This is a direct consequence of $`\widehat{T}`$, our representation of the Laplacian operator $`^2`$ on the grid, being Hermitian. As mentioned earlier, this is an important point. The matrix elements $`T_{\alpha \beta }`$ may always be made Hermitian by construction without $`\widehat{T}`$ itself being an Hermitian operator. This would ensure real eigenvalues, as is required. However, when it comes to optimisation of the support functions during a total-energy calculation, we require the derivative of the kinetic energy with respect to the support function values :

$$\frac{E_T}{\varphi _\alpha ^{}(𝐫_i)}=\frac{1}{2}\underset{\beta }{}K^{\beta \alpha }\widehat{T}\varphi _\beta (𝐫_i),\frac{E_T}{\varphi _\alpha (𝐫_i)}=\frac{1}{2}\underset{\beta }{}K^{\alpha \beta }\widehat{T}\varphi _\beta ^{}(𝐫_i),$$
(9)
where the $`𝐫_i`$ are grid points belonging to the support region of $`\varphi _\alpha `$. These relations both hold only if $`\widehat{T}`$ is an Hermitian operator, and support function optimisation can only be performed in a consistent manner if there is one unique representation of $`\widehat{T}\varphi _\beta `$ for each support function $`\varphi _\beta `$. It is also worth noting that the evaluation of these derivatives is the reason why we prefer to perform the sum of equation (4) for the FFT box method in real space, rather than in its equivalent form in reciprocal space. Applying the FFT box method in reciprocal space would be no more costly as far as integral evaluation is concerned but we would require an extra FFT per support function for the subsequent evaluation of equation (9).

For all the methods we describe in this paper we observe variation in the values of the kinetic energy integrals when we translate the system of the two support functions with respect to the real space grid. This is to be expected as the discrete representation of the support functions changes with the position of the support region with respect to the grid. Such variations may have undesirable consequences when it comes to calculating the forces on the atoms. In FFT terminology, they result from irregular aliasing of the high frequency components of our support functions as they are translated in real space. Ideally, in order to avoid this effect, the reciprocal representation of the support functions should contain frequency components only up to the maximum frequency that corresponds to our grid spacing, in other words it should be strictly localised in reciprocal space. Unfortunately this constraint is not simultaneously compatible with strict real space localisation. It should be possible however to achieve a compromise, thus controlling the translation error by making it smaller than some threshold. Such a compromise should involve an increase in the support region radii of our functions by a small factor. This situation is similar to the calculation of the integrals of the nonlocal projectors of pseudopotentials in real space with the method of King-Smith et al.  which requires an increase of the core radii by a factor of 1.5 to 2. For example, if we consider two carbon valence pseudo-orbitals of support radius $`6.0a_0`$ and with $`d=5.0a_0`$ and translate them both in a certain lattice vector direction over a full grid spacing, the maximum variation in the value of the integral with the FFT box method is $`8.28\times 10^6`$ Hartree. If we then do the same with carbon pseudo-orbitals generated with precisely the same parameters but instead with a support radius of $`10.0a_0`$, the maximum variation with respect to translation is reduced to $`2.05\times 10^8`$ Hartree.

## 4 Conclusions

In conclusion, we have presented a new and easy to implement method for calculating kinetic energy matrix elements of localised functions represented on a regular real space grid. This FFT box method is based on a mixed real space – reciprocal space approach. We use well established FFT algorithms to calculate the action of the Laplacian operator on localised support functions, whilst maintaining linear-scaling with system size and near locality of the operation. This makes our FFT box method suitable for implementation in the order-N code that we are developing. We have performed tests of the FFT box method and various orders of FD. Comparing to the exact integrals of the continuous representation, we have demonstrated that our approach is more accurate than low order FD approximations and only when A = 28 FD is used does the accuracy become comparable to that of the FFT box method. We have also highlighted the connection between the FFT box method and plane-wave methods and shown that our approach is up to three orders of magnitude more accurate than A = 28 FD when compared to the ‘exact’ result within the plane-wave basis set of the entire simulation cell. Furthermore, our approach for calculating the kinetic energy is consistent with the way in which other quantities in a total-energy calculation, such as the electron density and the Hartree energy, are computed as these are also calculated using reciprocal space techniques. Finally, we also note that our FFT box method is more versatile than FD as it is applicable to real space grids based on any lattice symmetry whereas FD schemes are usually only applied to orthorhombic grids.

## 5 Acknowledgements

C.-K. S. would like to thank the EPSRC (grant number GR/M75525) for postdoctoral research funding. A. A. M. would like to thank the EPSRC for a Ph.D. studentship. P. D. H. would like to thank Magdalene College, Cambridge for a Reseach Fellowship. C. J. P. would like to thank the EPSRC (grant number GR/N18567) for a ROPA.


text:
# Solution of the one-dimensional Dirac equation with a linear scalar potential

## I Introduction

The linear potential $`V(x)=g|x|`$ is a natural choice for a confining potential in one space dimension. The nonrelativistic Schrödinger equation admits a nearly analytic solution for this potential in terms of an Airy function and the zeros of this function and its derivative. The Dirac equation, on the other hand, appears to be problematic for this potential. If $`V`$ is introduced as the time component of a Lorentz two-vector, no bound-state solutions exist. If it is introduced as a Lorentz scalar, Bhalerao and Ram find only a very limited set of solutions, with no obvious correspondence to the nonrelativistic solutions. Such an outcome in the scalar case is unexpected because the Klein paradox is not a problem; positive and negative energy particles both see a confining potential. A nonrelativistic limit for the positive-energy solutions should reproduce the known nonrelativistic spectrum.

This inconsistency in the scalar case can be resolved. The solution found by Bhalerao and Ram turns out to be over constrained. Here we will construct a more general solution and show that the nonrelativistic results are recovered in an appropriate limit.

To see that the Dirac-equation solution should match on to the nonrelativistic solution, consider the equation (with $`\mathrm{}=1=c`$)

$$[\alpha p+\beta (m+g|x|)]\psi =E\psi $$
(1)
in a representation where

$$\alpha \stackrel{~}{\alpha }\sigma _y=\left(\begin{array}{cc}\hfill 0& \hfill i\\ \hfill i& \hfill 0\end{array}\right),\beta \stackrel{~}{\beta }\sigma _z=\left(\begin{array}{cc}\hfill 1& \hfill 0\\ \hfill 0& \hfill 1\end{array}\right).$$
(2)
As usual, let $`\psi =\stackrel{~}{\psi }\left(\begin{array}{c}\stackrel{~}{u}\\ \stackrel{~}{v}\end{array}\right)`$ and decompose the matrix equation into two coupled equations

$`\stackrel{~}{v}^{}+(m+g|x|)\stackrel{~}{u}`$ $`=`$ $`E\stackrel{~}{u},`$ (3)
$`\stackrel{~}{u}^{}(m+g|x|)\stackrel{~}{v}`$ $`=`$ $`E\stackrel{~}{v}.`$ (4)
For energies $`E=m+\stackrel{~}{\epsilon }`$ near the rest mass $`m`$, with $`\stackrel{~}{\epsilon }m`$, and for weak coupling $`gm^2`$, the second equation yields $`\stackrel{~}{v}\stackrel{~}{u}^{}/2m`$. Substitution into the first equation brings $`\stackrel{~}{u}^{\prime \prime }/2m+g|x|\stackrel{~}{u}\stackrel{~}{\epsilon }\stackrel{~}{u}`$, which is immediately recognized as the nonrelativistic Schrödinger equation.

The solution to the Schrödinger equation is obtained by noting that $`\stackrel{~}{u}^{\prime \prime }=2m(g|x|\stackrel{~}{\epsilon })\stackrel{~}{u}`$ is the differential equation for a shifted Airy function. This yields the normalizable solution $`\stackrel{~}{u}(x)=𝒩\text{Ai}\left((2mg)^{1/3}[|x|\stackrel{~}{\epsilon }/g]\right)`$, with $`𝒩`$ a normalization constant. Continuity of $`\stackrel{~}{u}`$ and $`\stackrel{~}{u}^{}`$ at $`x=0`$ requires that either $`\text{Ai}^{}\left((2mg)^{1/3}\stackrel{~}{\epsilon }/g\right)=0`$ \[for even solutions\] or $`\text{Ai}\left((2mg)^{1/3}\stackrel{~}{\epsilon }/g\right)=0`$ \[for odd solutions\]. Let $`\rho _n`$ and $`\rho _n^{}`$ denote the nth zeros of Ai and $`\text{Ai}^{}`$, respectively. Then the nonrelativistic eigenenergies are $`\stackrel{~}{\epsilon }_n^{}=\rho _n^{}(g^2/2m)^{1/3}`$ for even solutions and $`\stackrel{~}{\epsilon }_n=\rho _n(g^2/2m)^{1/3}`$ for odd. The values of $`\rho _n`$ and $`\rho _n^{}`$ can be obtained from tables in Ref. \[IV\]. The first four of each are as follows: $`\rho _n=2.3381`$, 4.0879, 5.5206, 6.7867 and $`\rho _n^{}=1.0188`$, 3.2482, 4.8201, 6.1633.

We would expect the Dirac equation to yield these same results in the limit of weak coupling. To see that negative energy solutions do not cause any difficulties, we can use the methods of Coutinho, Nogami, and Toyama to prove the following extension of their theorem B: For a scalar potential that is everywhere nonnegative, the positive energy solutions have energy $`Em`$ and the negative energy solutions have $`Em`$. The proof depends on the freedom to pick as real the solutions $`\stackrel{~}{u}`$ and $`\stackrel{~}{v}`$ to the coupled equations (3). Inner products of $`\stackrel{~}{u}`$ and $`\stackrel{~}{v}`$ with the terms of these equations yield

$`{\displaystyle \stackrel{~}{v}\stackrel{~}{u}^{}𝑑x}+{\displaystyle \stackrel{~}{u}(m+V)\stackrel{~}{u}𝑑x}`$ $`=`$ $`E{\displaystyle \stackrel{~}{u}^2𝑑x},`$ (5)
$`{\displaystyle \stackrel{~}{v}\stackrel{~}{u}^{}𝑑x}{\displaystyle \stackrel{~}{v}(m+V)\stackrel{~}{v}𝑑x}`$ $`=`$ $`E{\displaystyle \stackrel{~}{v}^2𝑑x},`$ (6)
where $`g|x|`$ has been replaced by a generic scalar potential $`V`$ and an integration by parts has been performed in the first term of the first equation. For positive $`E`$ and $`V`$, the second equation implies that $`\stackrel{~}{v}\stackrel{~}{u}^{}𝑑x0`$ and then the first equation yields $`Em`$. Analogous steps for negative $`E`$ yield $`\stackrel{~}{v}\stackrel{~}{u}^{}𝑑x0`$ and $`Em`$. Thus the two parts of the spectrum are completely separate, and we are allowed to focus on the positive-energy solutions only.

## II Solution of the Dirac Equation

To solve the Dirac equation directly, we use the same representation as Bhalerao and Ram, that is

$$\alpha \sigma _y=\left(\begin{array}{cc}\hfill 0& \hfill i\\ \hfill i& \hfill 0\end{array}\right),\beta \sigma _x=\left(\begin{array}{cc}\hfill 0& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right).$$
(7)
For $`\psi =\left(\begin{array}{c}u\\ v\end{array}\right)`$ they obtain the coupled equations

$`u^{}+(m+g|x|)u`$ $`=`$ $`Ev,`$ (8)
$`v^{}+(m+g|x|)v`$ $`=`$ $`Eu.`$ (9)
These equations decouple in terms of their variable $`\xi =\sqrt{g}(m/g+|x|)`$, such that for $`x>0`$,

$`\left({\displaystyle \frac{d^2}{d\xi ^2}}+\xi ^2\right)u=(E^2/g+1)u,`$ (10)
$`\left({\displaystyle \frac{d^2}{d\xi ^2}}+\xi ^2\right)v=(E^2/g1)v,`$ (11)
and for $`x<0`$,

$`\left({\displaystyle \frac{d^2}{d\xi ^2}}+\xi ^2\right)u=(E^2/g1)u,`$ (12)
$`\left({\displaystyle \frac{d^2}{d\xi ^2}}+\xi ^2\right)v=(E^2/g+1)v.`$ (13)
Obviously these are harmonic-oscillator-type equations. The normalizable solutions can be constructed from the Hermite functions of order $`\nu `$ and $`\nu +1`$, where $`E^2=2(\nu +1)g`$, as

$`u`$ $`=`$ $`\{\begin{array}{cc}Ce^{\xi ^2/2}H_{\nu +1}(\xi ),\hfill & x>0\hfill \\ C^{}\frac{E}{\sqrt{g}}e^{\xi ^2/2}H_\nu (\xi ),\hfill & x<0,\hfill \end{array}`$ (16)
$`v`$ $`=`$ $`\{\begin{array}{cc}C\frac{E}{\sqrt{g}}e^{\xi ^2/2}H_\nu (\xi ),\hfill & x>0\hfill \\ C^{}e^{\xi ^2/2}H_{\nu +1}(\xi ),\hfill & x<0.\hfill \end{array}`$ (19)
However, because $`\xi `$ is always positive, $`\nu `$ is not restricted to being an integer.

Continuity at $`x=0`$ requires that $`CH_{\nu +1}(\alpha )=C^{}EH_\nu (\alpha )/\sqrt{g}`$ and $`C^{}H_{\nu +1}(\alpha )=CEH_\nu (\alpha )/\sqrt{g}`$, where $`\alpha m/\sqrt{g}`$ is defined as in Ref. \[IV\]. These conditions yield $`C^{}=\pm C`$ and

$$H_{\nu +1}(\alpha )=\pm \frac{E}{\sqrt{g}}H_\nu (\alpha ),$$
(20)
the latter being the eigenvalue condition. This condition has a rich set of solutions when free of the restriction to integer $`\nu `$; there are infinitely many solutions for any positive value of $`\alpha `$.

The sign that appears in the eigenvalue condition (20) corresponds to the parity of the solution. The parity operator is reflection in $`x`$ combined with multiplication by the Dirac matrix $`\beta `$. Because $`\xi `$ is independent of the sign of $`x`$, we find that

$$\beta \left(\begin{array}{c}u(x)\hfill \\ v(x)\hfill \end{array}\right)=\left(\begin{array}{c}v(x)\hfill \\ u(x)\hfill \end{array}\right)=\pm \left(\begin{array}{c}u(x)\hfill \\ v(x)\hfill \end{array}\right).$$
(21)
The presence of such a symmetry is, of course, necessary for the match to the nonrelativistic solution.

## III Nonrelativistic limit

To recover the nonrelativistic solution, we must take an appropriate limit. We write $`E=m(1+\epsilon )`$ and consider small $`\epsilon `$ as well as small $`g`$. The latter corresponds to large $`\alpha `$ and large $`\nu `$. In the limit of large $`\nu `$, we find from Ref. \[IV\] that $`H_\nu `$ has the asymptotic form

$$H_\nu (\xi )2^\nu e^{\xi ^2/2}\mathrm{\Gamma }\left(\frac{\nu +1}{2}\right)\left(\frac{t_\nu }{z_\nu ^21}\right)^{1/4}\text{Ai}(t_\nu ),$$
(22)
where $`z_\nu =\xi /\sqrt{2\nu +1}`$ and, for $`z_\nu 1`$,

$$t_\nu =\left(\frac{3}{4}(2\nu +1)\left[\mathrm{cos}^1z_\nu z_\nu \sqrt{1z_\nu ^2}\right]\right)^{2/3}.$$
(23)
This immediately looks promising because the desired Airy function is present. We next expand in powers of $`\alpha ^1`$ and $`\epsilon `$, with $`y_\pm \epsilon \sqrt{g}|x|/\alpha \pm 1/2\alpha ^2`$ and use of $`\mathrm{cos}^1(1y)\sqrt{2y}+\sqrt{y^3/72}`$, to obtain

$`z_\nu `$ $`=`$ $`{\displaystyle \frac{1+\sqrt{g}|x|/\alpha }{\sqrt{1+2\epsilon +\epsilon ^21/\alpha ^2}}}1y_{},z_{\nu +1}1y_+`$ (24)
$`t_\nu `$ $``$ $`\left({\displaystyle \frac{2m^4}{g^2}}\right)^{1/3}y_{},t_{\nu +1}\left({\displaystyle \frac{2m^4}{g^2}}\right)^{1/3}y_+.`$ (25)
Note that $`\sqrt{g}|x|/\alpha `$ is of order $`\epsilon `$ at the classical turning point, which sets a natural scale for $`x`$, and that the nonrelativistic correspondence implies that $`\epsilon `$ is of order $`(g^2/m^4)^{1/3}\alpha ^{4/3}`$. Thus the $`\epsilon ^2`$ terms can be dropped relative to $`\alpha ^2`$.

At lowest order, these expansions imply

$$e^{\xi ^2/2}H_\nu (\xi ),e^{\xi ^2/2}H_{\nu +1}(\xi )\text{Ai}\left((2m^4/g^2)^{1/3}[|x|/\alpha \epsilon ]\right).$$
(26)
To compare with the original nonrelativistic reduction we must connect the two representations of the Dirac matrices. They are related by a unitary transformation

$$U=\frac{i}{\sqrt{2}}\left(\begin{array}{cc}\hfill 1& \hfill 1\\ \hfill 1& \hfill 1\end{array}\right),$$
(27)
such that

$$\stackrel{~}{\psi }=\left(\begin{array}{cc}\stackrel{~}{u}& \\ \stackrel{~}{v}& \end{array}\right)=U\psi =U\left(\begin{array}{cc}u& \\ v& \end{array}\right).$$
(28)
Therefore, we have $`\stackrel{~}{u}=\frac{i}{\sqrt{2}}(u+v)`$ and $`\stackrel{~}{v}=\frac{i}{\sqrt{2}}(vu)`$, with $`u`$ and $`v`$ given by (16). Thus for large $`\alpha `$, $`\stackrel{~}{u}`$ does indeed reduce to an Airy function with the correct argument, given $`\epsilon =\stackrel{~}{\epsilon }/m`$ and $`\alpha =m/\sqrt{g}`$.

At the next order we obtain, with the aid of Stirling’s formula for $`\mathrm{\Gamma }\left(\frac{\nu +2}{2}\right)`$,

$$H_\nu (\xi )2^{\nu 1/4}e^{\xi ^2/2}\mathrm{\Gamma }\left(\frac{\nu +1}{2}\right)\left(\frac{2m^4}{g^2}\right)^{1/12}(1+y_{}/4)\text{Ai}\left((2m^4/g^2)^{1/3}y_{}\right)$$
(29)
and

$$H_{\nu +1}(\xi )2^{\nu +3/4}e^{\xi ^2/2}\sqrt{\frac{\nu +1}{2}}\mathrm{\Gamma }\left(\frac{\nu +1}{2}\right)\left(\frac{2m^4}{g^2}\right)^{1/12}(1+y_+/4)\text{Ai}\left((2m^4/g^2)^{1/3}y_+\right).$$
(30)
The combination that appears in the eigenvalue condition

$`0=H_{\nu +1}(\alpha )\sqrt{2\nu +2}H_\nu (\alpha )`$

then reduces to

$`0`$ $`\sqrt{2\nu }2^{\nu 1/4}e^{\alpha ^2/2}\mathrm{\Gamma }\left({\displaystyle \frac{\nu +1}{2}}\right)\left({\displaystyle \frac{2m^4}{g^2}}\right)^{1/12}`$ (32)
$`\times \left[(1+y_{}/4)\text{Ai}\left((2m^4/g^2)^{1/3}y_+\right)(1+y_+/4)\text{Ai}\left((2m^4/g^2)^{1/3}y_{}\right)\right]_{x=0}.`$
First-order Taylor expansions of the Airy functions about $`(2m^4/g^2)^{1/3}\epsilon `$ yield

$`0`$ $`\text{Ai}\left((2m^4/g^2)^{1/3}\epsilon \right)(2\alpha )^{2/3}\text{Ai}^{}\left((2m^4/g^2)^{1/3}\epsilon \right)`$ (34)
$`\left[\text{Ai}\left((2m^4/g^2)^{1/3}\epsilon \right)+(2\alpha )^{2/3}\text{Ai}^{}\left((2m^4/g^2)^{1/3}\epsilon \right)\right],`$
where terms of order higher than $`\alpha ^{2/3}`$ have been dropped. For even parity (the upper sign) we have $`\text{Ai}^{}\left((2m^4/g^2)^{1/3}\epsilon \right)0`$ and for odd parity, $`\text{Ai}\left((2m^4/g^2)^{1/3}\epsilon \right)0`$, which are the nonrelativistic eigenvalue conditions. Therefore the eigenvalues will match in the limit of large $`\alpha `$ and small $`\epsilon `$.

## IV Results and conclusions

Comparison of the relativistic and nonrelativistic eigenvalues is made in Fig. 1, where each is plotted as a function of $`1/\alpha `$ for the four lowest levels of each parity. The nonrelativistic values were already obtained above as explicit functions of $`\alpha `$ which are simply plotted as lines in the figure. The relativistic values were computed numerically, with use of Mathematica to solve the eigenvalue condition (20) at selected values of $`\alpha `$. The dimensionless $`\epsilon `$ is related to $`\nu `$ by $`\epsilon =\sqrt{2\nu +2}/\alpha 1`$. For $`1/\alpha `$ near 0, i.e. large $`\alpha `$, the relativistic and nonrelativistic results are indistinguishable. As $`\alpha `$ decreases, they separate smoothly.

From these results we see that the Dirac equation with a scalar linear potential is a well-defined problem in one dimension, with a rich set of solutions and a smooth nonrelativistic limit. As an exercise, one could extend this work to include calculation of the relativistic wave functions and make direct comparisons with the nonrelativistic Airy functions. They will match in the large-$`\alpha `$ limit. A second interesting exercise is a comparison of the ultrarelativistic, strong-coupling limit of the spectrum to the $`\alpha ^{4/3}`$ behavior of the nonrelativistic spectrum. The plots in Fig. 1 appear to imply a $`\alpha ^1`$ behavior for small $`\alpha `$.


text:
# 𝑁=1 Supergravity Chaotic Inflation in the Braneworld Scenario

## I Introduction

Chaotic inflationary models Linde  stand out for their simplicity and fairly natural initial conditions for the onset of inflation. These features are particularly appealing in the context of supergravity and superstring theories, where the natural scale for fields is the Planck scale. In the context of $`N=1`$ Supergravity, however, realizations of chaotic inflation in minimal and in $`SU(1,1)`$ theories are somewhat special as to ensure sufficient inflation Goncharov ; Goncharov1 . Furthermore, chaotic inflation requires super-Planckian field values both to ensure sufficient inflation and the correct amount of Cosmic Microwave Background (CMB) anisotropies, in which case higher order non-renormalizable terms in the superpotential would completely dominate the dynamics since no well motivated symmetry is known to prevent them. Assuming that this problem can be circumvented somehow, one finds that, in order to accommodate in a satisfactory way the bounds on the reheating temperature and energy density fluctuations, a two-scale chaotic inflationary sector is required Bento1 , in contrast with the situation found in $`N=1`$ Supergravity new inflationary type models, where a single scale suffices Ross ; Holman .

In this work, we show that specific features of the braneworld scenario allow for the abovementioned difficulties with higher order non-renormalizable terms to be quite naturally avoided and that current observational contraints can be accounted for with a single scale at the superpotential level.

Higher dimensional superstring motivated cosmological solutions suggest that matter fields that are related to open string modes lie on a lower dimensional brane while gravity propagates in the bulk Polchinski . It is striking that, in these scenarios, extra dimensions are not restricted to be small Antoniadis  and that the fundamental $`D`$-dimensional scale, $`M_D`$, where $`D=4+d`$, can be considerably smaller than the 4-dimensional Planck scale. In this work, we shall consider the $`D=5`$ case. If one assumes that Einstein equations with a negative cosmological constant hold (an anti-De-Sitter space is required) in $`D`$-dimensions and that matter fields are confined to the $`3`$-brane then the $`4`$-dimensional Einstein equation is given by Shiromizu :

$$G_{\mu \nu }=\mathrm{\Lambda }g_{\mu \nu }+\frac{8\pi }{M_P^2}T_{\mu \nu }+\left(\frac{8\pi }{M_5^3}\right)^2S_{\mu \nu }E_{\mu \nu },$$
(1)
where $`T_{\mu \nu }`$ is the energy-momentum on the brane, $`S_{\mu \nu }`$ is a tensor that contains contributions that are quadratic in $`T_{\mu \nu }`$ and $`E_{\mu \nu }`$ corresponds to the projection of the 5-dimensional Weyl tensor on the 3-brane (physically, for a perfect fluid, it is associated to non-local contributions to the pressure and energy flux). The $`4`$-dimensional cosmological constant is related to the $`5`$-dimensional cosmological constant and the $`3`$-brane tension, $`\lambda `$, as

$$\mathrm{\Lambda }=\frac{4\pi }{M_5^3}\left(\mathrm{\Lambda }_5+\frac{4\pi }{3M_5^3}\lambda ^2\right),$$
(2)
while the Planck scale is given by

$$M_P=\sqrt{\frac{3}{4\pi }}\frac{M_5^3}{\sqrt{\lambda }}.$$
(3)
In a cosmological setting, where the $`3`$-brane resembles our Universe and the metric projected onto the brane is an homogeneous and isotropic flat Robertson-Walker metric, the generalized Friedmann equation has the following form Binetruy :

$$H^2=\frac{\mathrm{\Lambda }}{3}+\left(\frac{8\pi }{3M_P^2}\right)\rho +\left(\frac{4\pi }{3M_5^3}\right)^2\rho ^2+\frac{ϵ}{a^4},$$
(4)
where $`ϵ`$ is an integration constant arising from $`E_{\mu \nu }`$. During inflation, the last term in (4), the “dark radiation” term, rapidly vanishes and will be disregarded hereafter. Moreover, observations require that the cosmological constant is negligible in the early Universe, meaning that $`\mathrm{\Lambda }_5`$ and the brane tension were fine tuned, $`\mathrm{\Lambda }_54\pi \lambda ^2/3M_5^3`$. The Friedmann equation can then be written in the following way:

$$H^2=\frac{8\pi }{3M_P^2}\rho \left[1+\frac{\rho }{2\lambda }\right].$$
(5)
Notice that the new term in $`\rho ^2`$ is dominant at high energies, compared to $`\lambda ^{1/4}`$, but quickly decays at lower energies. Requiring the new term to be sub-dominant during nucleosynthesis, implies that $`\lambda \text{ }>(1\text{MeV})^4`$ and, therefore Cline :

$$M_5\text{ }>\left(\frac{1\text{MeV}}{M_P}\right)^{2/3}M_P=10\text{TeV}.$$
(6)
A much more stringent bound can be obtained in models where the fifth dimension is infinite as it yields an extra contribution $`M_5^6\lambda ^2r^2`$ to Newton’s gravitational force, which should be, of course, small beyond scales $`r\text{ }>1mm`$ (see eg. Ref. Floratos  and references therein). From Eq. (3), it follows that $`M_5>10^5\text{TeV}`$. We shall see in the ensuing discussion that cosmological considerations intrinsic to supersymetric cosmology can set even more stringent bounds on $`M_5`$.

In what follows, we shall discuss the effect of the new term in the Friedmann equation on the slow roll-over parameters and show that it allows for a single scale chaotic inflationary model, in the context of $`N=1`$ Supergravity. Moreover, we study how parameters in the superpotential and $`\lambda `$ are affected by constraints on the magnitude of the energy density perturbations required to explain the anisotropies in the CMB radiation observed by COBE as well as on the reheating temperature. In fact, in supergravity cosmological models one has to ensure that the reheating temperature does not exceed $`T_{RH}\text{ }<2.5\times 10^8(100\text{GeV}/m_{3/2})\text{GeV}`$ Ellis , in order not to generate an over-abundance of gravitinos and the ensuing photo-dissociation of light elements at nucleosynthesis. In the context of superstring cosmology, inflationary models have to face further problems such as the fate of the dilaton and moduli fields, the so-called post-modern Polonyi problem Coughlan .

We also show that the observational quantities arising from the model originate in sub-Planckian field regimes and, hence, higher dimensional non-renormalizable operators in the superpotential are not relevant.

## II The Model

We shall assume that the inflaton is the scalar component of a gauge singlet superfield, $`\mathrm{\Phi }`$, in the hidden sector of the theory. We start by splitting the superpotential in a supersymmetry-breaking, a gauge and an inflationary sector:

$$W=P+G+I.$$
(7)
The $`N=1`$ supergravity theory describing the interaction of gauge singlet fields is specified by the Kähler function, in terms of which the scalar potential is given by

$$V=\frac{1}{4}e^K\left[G_a(K^1)_b^aG^b3|W|^2\right]$$
(8)
where $`G_a=K_aW+W_a`$, and the indices $`a,b`$ denote derivatives with respect to the chiral superfields, $`\mathrm{\Phi }`$. We consider the minimal choice for the Kähler function in $`N=1`$ Supergravity:

$$K(\mathrm{\Phi },\mathrm{\Phi }^{})=\mathrm{\Phi }^{}\mathrm{\Phi },$$
(9)
corresponding to canonical kinetic energy for the scalar fields. With this choice, the scalar potential for the inflaton field, $`\varphi `$, is obtained from the superpotential $`I(\mathrm{\Phi })`$ as

$$V(\varphi )=e^{\frac{|\varphi |^2}{M^2}}\left[\left|\frac{I}{\mathrm{\Phi }}+\frac{\mathrm{\Phi }^{}I}{M^2}\right|^23\frac{|I|^2}{M^2}\right]_{\mathrm{\Phi }=\varphi },$$
(10)
where $`M=M_P/\sqrt{8\pi }`$.

Requiring the cosmological constant to vanish and that supersymmetry remains unbroken at the minimum of the potential, $`\mathrm{\Phi }=\mathrm{\Phi }_o`$, leads to the following constraints on the superpotential:

$$I(\mathrm{\Phi }_o)=\frac{I}{\mathrm{\Phi }}(\mathrm{\Phi }_o)=0.$$
(11)
Consider the simplest form for the superpotential $`I`$ which satisfies the above conditions

$$I(\mathrm{\Phi })=\mu (\mathrm{\Phi }\mathrm{\Phi }_0)^2,$$
(12)
where $`\mu `$ is a mass parameter that determines the energy scale for inflation. Hereafter, we set $`\mathrm{\Phi }_0=0`$. The relevant part of the inflaton potential (along the real $`\varphi `$ direction) is then given by:

$$V(\varphi )=\mu ^2\left(4\varphi ^2+\frac{\varphi ^4}{M^2}+\frac{\varphi ^6}{M^4}\right)e^{\frac{\varphi ^2}{M^2}}.$$
(13)
Consistency between the slow-roll approximation and the full evolution equations requires that there are constraints on the slope and curvature of the potential. One can define two slow-roll parameters Maartens

$`ϵ`$ $``$ $`{\displaystyle \frac{M_P^2}{16\pi }}\left({\displaystyle \frac{V^{}}{V}}\right)^2{\displaystyle \frac{1+V/\lambda }{(2+V/\lambda )^2}},`$ (14)
$`\eta `$ $``$ $`{\displaystyle \frac{M_P^2}{8\pi }}\left({\displaystyle \frac{V^{\prime \prime }}{V}}\right){\displaystyle \frac{1}{1+V/2\lambda }}.`$ (15)
Notice that both parameters are suppressed by an extra factor $`\lambda /V`$ at high energies and that, at low energies, $`V\lambda `$, they reduce to the standard form. The end of inflation will take place for a field value, $`\varphi _f`$, such that

$$\mathrm{max}\{ϵ(\varphi _f),|\eta (\varphi _f)|\}=1.$$
(16)
The number of e-folds during inflation is given by $`N=_{t_\mathrm{i}}^{t_\mathrm{f}}H𝑑t`$, which becomes Maartens

$$N\frac{8\pi }{M_P^2}_{\varphi _\mathrm{i}}^{\varphi _\mathrm{f}}\frac{V}{V^{}}\left[1+\frac{V}{2\lambda }\right]𝑑\varphi ,$$
(17)
in the slow-roll approximation. We see that, as a result of the modified Friedmann equation at high energies, the expansion rate is increased by a factor $`V/2\lambda `$, allowing for a smaller initial inflaton field value, $`\varphi _\mathrm{i}`$, for a given number of e-folds, which is crucial for achieving the goal of obtaining sufficient inflation with sub-Planckian field values in our model.

## III Constraints from Reheating

After inflation, the field $`\varphi `$ releases its energy via the coupling to fields in the other sectors in (7), thus reheating the Universe. Since the inflaton is hidden from the other sectors of the theory, it couples to lighter fields with gravitational strength of the order $`\frac{\mu }{M}`$. At minimum, the inflaton field has a mass

$$m_\varphi =2\sqrt{2}\mu ,$$
(18)
leading to a decay width

$$\mathrm{\Gamma }_\varphi \frac{m_\varphi }{(2\pi )^3}\left(\frac{\mu }{M}\right)^3,$$
(19)
and a reheating temperature

$`T_{RH}`$ $``$ $`\left({\displaystyle \frac{30}{\pi ^2g_{RH}}}\right)^{1/4}\sqrt{M\mathrm{\Gamma }}`$ (20)
$``$ $`{\displaystyle \frac{2}{\pi ^2}}\left(\sqrt{{\displaystyle \frac{15}{g_{RH}}}}{\displaystyle \frac{\mu ^3}{M}}\right)^{1/2},`$
where $`g_{RH}`$ is the number of degrees of freedom at $`T_{RH}`$.

As already mentioned, a quite severe upper bound on $`T_{RH}`$ comes from the requirement that gravitinos are not abundantly regenerated in the post-inflationary reheating epoch. Indeed, once regenerated beyond a certain density, stable thermal gravitinos would dominate the energy density of the Universe or, if they decay, have disruptive effects on nucleosynthesis, causing light element photo-dissociation and distortions in the CMB. Avoiding these difficulties implies in the following bounds Ellis :

$$T_{RH}\text{ }<2\times 10^9,6\times 10^9\text{GeV}$$
(21)
for $`m_{3/2}=1,10\text{TeV}`$. For our model, demanding that $`T_{RH}`$ be less than $`2\times 10^9\text{GeV}`$ leads, for $`g_{RH}=150`$, to a limit on parameter $`\mu `$

$$\mu \text{ }<3.7\times 10^6M,$$
(22)
which coincides with the bound obtained in Ref. Bento1  for $`SU(1,1)`$ Supergravity since, in either case, reheating is naturally controlled by the lowest order quadratic term in the superpotential.

We have checked that the above bound also ensures that gravitino production via inflaton decay is sufficiently suppressed, as in Bento1 ; Ross .

Naturally, one could contemplate other more specifical reheating mechanisms. Parametric resonance Kofman  seems to be unimportant given that $`m_\varphi >>\mathrm{\Gamma }_\varphi `$ and the absence of bilinear type couplings in supergravity. On the other hand, anharmonic terms coupling the inflaton with other chiral superfields can be an alternative route for reheating Bertolami . In the braneworld scenario, gravitational particle production has been discussed for an exponentially decaying inflaton potential Copeland .

## IV Constraints from CMB anisotropies

The amplitude of scalar perturbations is given by Maartens

$$A_s^2\left(\frac{512\pi }{75M_P^6}\right)\frac{V^3}{V^2}\left[1+\frac{V}{2\lambda }\right]^3|_{k=aH},$$
(23)
where the right-hand side should be evaluated as the comoving scale equals the Hubble radius during inflation, $`k=aH`$. Thus the amplitude of scalar perturbations is increased relative to the standard result at a fixed value of $`\varphi `$ for a given potential.

In order to obtain the value of $`\varphi `$ when scales corresponding to large-angle CMB anisotropies, as observed by COBE, left the Hubble radius during inflation, we take $`N_{}55`$ and $`\varphi _i=\varphi _{}`$ in Eq. (17). Combining with Eq. (23) and using the fact that the observed value from COBE is $`A_s=2\times 10^5`$, we obtain, after a numerical analysis, a further (stronger) constraint on $`\mu `$

$$\mu \text{ }<2\times 10^8M.$$
(24)
For $`\mu =2\times 10^8M`$, we get $`\varphi _{}=0.09M_P`$ and $`M_5=3\times 10^4M_P`$ (for $`\varphi _f=0.03M_P`$, as obtained from Eq. (16)).

The scale-dependence of the perturbations is described by the spectral tilt Maartens

$$n_s1\frac{d\mathrm{ln}A_s^2}{d\mathrm{ln}k}6ϵ+2\eta ,$$
(25)
where the slow-roll parameters are given in Eqs. (14) and (15). Notice that, as $`V/\lambda \mathrm{}`$, the spectral index is driven towards the Harrison-Zel’dovich spectrum, $`n_s1`$. For our model, Eq. (13), we obtain

$$n_s0.95\text{for}\mu =2\times 10^8M.$$
(26)
The ratio between the amplitude of tensor and scalar perturbations is given by Langlois

$$\frac{A_t^2}{A_s^2}\frac{3M_P^2}{16\pi }\left(\frac{V^{}}{V}\right)^2\frac{2\lambda }{V}.$$
(27)
We get, for our model

$$r4\pi \frac{A_T^2}{A_s^2}0.22\text{for}\mu =2\times 10^8,$$
(28)
which is consistent with current upper limits, $`r<0.4`$ (see eg. DASI ).

For consistency, it should be verified that enough inflation can occur. Indeed, for $`\mu =2\times 10^8`$, we get $`N=65`$ for $`\varphi _i=0.094M_P`$ and a huge amount of e-foldings can be obtained for larger values of $`\varphi _i`$ e.g. $`N=9.1\times 10^6`$ for $`\varphi _i=0.5M_P`$.

For smaller values of $`\mu `$ (see Table 1), we get smaller values of $`\varphi _{}`$ and $`M_5`$, but $`n_s`$ and $`r`$ change very little. The number of e-folds of inflation increases for the same initial value of $`\varphi `$ and $`\varphi _f`$, calculated according to Eq. (16); for instance, for $`\mu =10^{10}M_P`$, $`\varphi _i=0.094M_P`$ and $`\varphi _f=9\times 10^5M_P`$, we get $`N=3\times 10^{11}`$, about five orders of magnitude greater than the value obtained for the case $`\mu =2\times 10^8M`$ (see above). Notice that the values of $`n_s`$ and $`r`$ exhibited in Table 1 are consistent with latest CMB data from DASI DASI , BOOMERANG Boom  and MAXIMA Maxima . However, for $`\mu \text{ }<1\times 10^{10}M`$, the reheating temperature becomes smaller than the typical value for the temperature of the electroweak phase transition, $`T_{EW}300\text{GeV}`$. Thus, $`\mu \text{ }<1\times 10^{10}M`$, implies a premature breaking of the electroweak phase transition, from which we can extract a new bound on the 5-dimensional mass scale (see Table 1)

$$M_5\text{ }>1.3\times 10^6M_P.$$
(29)
Naturally, as field values become smaller, the approximation which consists of expanding the potential in $`\varphi /M`$ and take just the first (quadratic) term becomes better; in this case, it is easy to obtain analytic expressions for the relevant quantities, as it is reduces to the case analysed in Ref. Maartens .

## V Conclusions

Hence, we see that it is possible to sucessfully implemment chaotic inflation in a simple $`N=1`$ Supergravity model within the braneworld scenario, without need to fine tune the parameters of the potential. In fact, a single mass parameter, $`\mu `$, is needed, as in $`N=1`$ Supergravity new inflationary models Ross ; Holman .

An upper bound on $`\mu `$, $`\mu \text{ }<3.7\times 10^6M`$, is obtained from the requirement that sufficiently few gravitinos are regenerated in the post-inflationary reheating epoch. However, a more severe upper bound is obtained from requiring adequate density fluctuations (both slope and amplitude) as observed by COBE, namely $`\mu \text{ }<2\times 10^8M`$ implying, for the 5-dimensional mass scale that $`M_5\text{ }<3\times 10^4M_P`$. Requiring that the reheating temperature be greater than the typical temperature of the electroweak phase transition, $`T_{EW}300\text{GeV}`$, leads to upper bounds on these quantities, namely $`\mu \text{ }>1.06\times 10^{10}M`$ and $`M_5\text{ }>1.3\times 10^6M_P`$.

Finally, we have shown that, remarkably, successful chaotic inflation can be achieved with sub-Planckian field values, thus avoiding the need to invoke hypothetic symmetries that would prevent the presence of higher order non-renormalisable terms in the superpotential. This last feature is important to prevent an overproduction of tensor perturbations in the CMB, a problem that can be particularly acute in chaotic inflation scenarios Enqvist , and that, ultimately due to the quadratic energy density term that appears in the Friedmann equation in the context of the braneworld scenario, is absent in our model.

###### Acknowledgements.

The authors acknowledge the partial support of FCT (Portugal) under the grant POCTI/1999/FIS/36285.