text:
# An efficient shock-capturing central-type scheme for multidimensional relativistic flows

## 1 Introduction

Relativistic flows and shocks play an essential role in modern high energy astrophysics, both for the interpretation of various observed features and for the description of the physical processes they give rise to, as, for example, the acceleration of highly energetic particles. Among the various astrophysical objects in which relativistic flows have been invoked to explain the observed properties, the best studied are probably:

1. Active galactic nuclei (AGNs), for which the presence of relativistic bulk motions (up to Lorentz factors of order 10) was soon suggested (Rees 1967). Associated to AGNs are often highly collimated relativistic jets (Begelman et al. 1984; and, for a review, Ferrari 1998), seen as apparent super-luminal lobes in some radio-loud sources and whose shock fronts are among the most promising candidates for high-energy particle acceleration. Mildly relativistic jets appear to be also associated to galactic X-ray compact sources (generally called microquasars, see Mirabel & Rodriguez 1999, for a review).
2. Gamma-ray bursts (GRBs), supposed to be originated during the collapse of the iron core of massive stars and the subsequent fireball explosion (Piran 1999; Kobayashi et al. 1999), which give rise to an expanding blast wave of pairs and hadrons with typical Lorentz factors of $`10^210^3`$ (Meszaros & Rees 1992), whose kinetic energy is then believed to be converted into gamma rays via cyclotron radiation and/or inverse Compton scattering.
3. Pulsar wind nebulae, assumed to be bubbles of relativistic particles and magnetic fields emitted by a pulsar as a relativistic wind with Lorentz factor ranging from $`10^4`$ to $`10^7`$, depending on the model of pair production in the pulsar magnetosphere (Michel & Li 1999). The wind region may be confined by a termination shock, generated by the interaction with outer supernova matter, as is the case for the synchrotron emitting nebulae called plerions (Rees & Gunn 1974; Kennel & Coroniti 1984); or, if the pulsar is moving with supersonic speed through the interstellar medium, by the resulting ram pressure, in this case giving rise to a bow-shock pulsar wind nebula that may be detected in H$`\alpha `$ (Bucciantini & Bandiera 2001).

As we can see from these few but important examples, there is a strong interest in the astrophysical community to the development of computational codes for the numerical modeling and simulation of relativistic flows.

Over the last decade, high resolution shock-capturing methods of Godunov type, successfully applied in classical fluid dynamics, have started to be employed for the case of relativistic hydrodynamics as well (Marquina et al. 1992; Schneider et al. 1993; Balsara 1994; Duncan & Hughes 1994; Eulderink & Mellema 1994; Font et al. 1994; Dolezal & Wong 1995; Falle & Komissarov 1996; Donat et al. 1998; Aloy et al. 1999). These schemes are characterized by the following main features: a conservative form of the discretized equations, in order to capture weak solutions and satisfy jump relations; a reconstruction phase, to recover variables at inter-cell locations where fluxes have to be computed; an upwind phase, in which an exact or approximate solution to the local Riemann problem is found. The simulated flows achieve high accuracy in smooth regions and, at the same time, shock profiles are stable and sharply defined. This is the reason for the great success of upwind differencing over central (or spectral) differencing. In the latter case, artificial viscosity must be introduced in order to damp the spurious oscillations that always form near discontinuities (Gibbs phenomena), thus leading to artificial heating and unwanted damping of physical waves.

When discontinuous solutions are of main interest, second order (both in time and space) total variation diminishing (TVD) schemes coupled with an accurate Riemann solver are probably the best choice, in terms of sharp resolution of discontinuities and overall stability. However, in the general case smooth wave-like or even turbulent fields appear together with discontinuous solutions, thus second order schemes are no longer able to resolve both features with enough precision. By relaxing the stringent TVD condition, a class of essentially non-oscillatory (ENO) higher-order schemes were first proposed by Harten et al. (1987), later modified in a more efficient implementation by Shu & Osher (1989). The ENO philosophy is to use adaptive stencils to reconstruct variables and fluxes at cell interfaces. Thus, in smooth regions symmetric stencils will be used, whereas near discontinuities the stencil will shift to the left or to the right, selecting the smoother part of the flow and thus achieving the same high resolution (typically from third to fifth order in the case of weighted ENO schemes, see Jiang & Shu 1996) everywhere, without the need of expensive adaptive mesh refinement (AMR) techniques. The price to pay is the presence, near shocks, of small oscillations of the order of the truncation error (Gibbs oscillations given by centered stencils that cross a discontinuity would be much higher, since they are proportional to the jump itself, independently on the resolution).

Like for most TVD Godunov-type schemes, also for ENO algorithms the local spectral decomposition in the building of numerical fluxes is commonly adopted. It was strongly advocated already in the original paper, where a version of the ENO scheme with the (much) simpler component-wise reconstructions was found to fail, giving significant oscillations, in some Riemann problems. However, due to the high resolution of ENO methods, characteristic decomposition on every point of the interpolation stencil becomes prohibitive when moving to three-dimensional simulations, especially for relativistic flows where Jacobian matrices are more complex than in the Eulerian case. The same kind of problem has to be faced in magnetohydrodynamics (MHD, see Londrillo & Del Zanna 2000, from now on LD), because of the increasing number of variables, equations and eigenmodes. The worst possible case is obviously that of relativistic MHD, for which only second order 1-D (Balsara 2001) and 2-D (Komissarov 1999) Godunov schemes have been presented so far.

Following the scheme proposed in LD, in the present series of two papers a shock-capturing scheme that avoids the expensive characteristic decomposition, still retaining non-oscillatory properties, will be suggested, first for ideal relativistic hydrodynamics (RHD, in this paper) and then for ideal relativistic magnetohydrodynamics (RMHD, in the next paper of the series). Our method is based on the so-called central schemes, that extend the first-order local Lax-Friedrichs (LLF) flux splitting and assume an average over the Riemann fan at every cell interface (see Nessyahu & Tadmor 1990, for the original paper and Kurganov et al. 2001, for the latest developments in this field). The basic properties of our method can be found in Liu & Osher (1998), where a third order multidimensional central scheme was presented, with the following main features:

1. Point values rather than cell averages are used (that is to say finite differences instead of finite volumes, see Sect. 3 for details), thus removing the need of staggered grids (employed in previous central schemes) and making the extension to multiple dimensions a trivial task.
2. The semi-discrete form of the equations is solved, so that time integration can be achieved with any solver for ordinary differential equations (ODEs), for example TVD Runge-Kutta methods.
3. No characteristic decomposition and Riemann solvers are required: fluxes are reconstructed and derived component-wise, thus achieving a great simplicity in the programming and, above all, efficiency. The only spectral pieces of information needed are the local highest characteristic speeds (related to the Courant coefficient).
4. A new third order reconstruction algorithm is introduced and tested, called Convex-ENO (CENO), which has the fundamental property to stay as close as possible to the lower order TVD limited linear reconstruction (by using the minmod limiter in this case), thus reducing even to first order where needed, while retaining high accuracy in smooth region. This was recognized to be the key point of the surprising success of central schemes (Tadmor, cited in their acknowledgements), confirmed for the MHD case by LD.

Here, like in LD, we take advantage of all these positive features and we modify the original scheme by splitting the CENO reconstruction into two separate routines, so that reconstruction can be applied to primitive variables rather than to fluxes (the resulting scheme appears to be more robust, especially in the relativistic case). Moreover, different limiters (e.g. the monotonized centered) and solvers (the upwind HLL, from Harten, Lax, and van Leer) are introduced and tested in our framework.

We will show that our scheme, by providing an overall high resolution, is able to compensate for the larger smearing of discontinuities, especially contact discontinuities (due to the use of solvers based on just one or two characteristic speeds) and the results are comparable with those obtained by more sophisticated (but much more computationally expensive, especially in 3-D) relativistic codes. Concluding, notice that the proposed scheme can be applied unchanged to Euler HD too, where only the definition of conservative variables, fluxes and characteristic speeds are different.

## 2 Ideal RHD equations

The covariant equations in special relativistic hydrodynamics (RHD) are (e.g. Landau & Lifshitz 1959):

$$_\alpha (\rho u^\alpha )=0,$$
(1)
$$_\alpha (wu^\alpha u^\beta +pg^{\alpha \beta })=0,$$
(2)
where velocities are normalized against the speed of light ($`c=1`$) and Greek letters indicate four-vectors, while Latin indexes will be devoted to spatial 3-D vectors. The metric tensor defines space-time properties and here a Minkowski flat space with $`g^{\alpha \beta }=\mathrm{diag}\{1,1,1,1\}`$ will be assumed throughout for simplicity, with coordinates $`x^\alpha =(t,x^j)`$. The physical quantities involved in the conservation laws are the rest mass density $`\rho `$, the kinetic pressure $`p`$, the relativistic enthalpy $`w=e+p`$ ($`e`$ is the energy per unit volume, including rest mass energy) and the four-velocity $`u^\alpha =(\gamma ,\gamma v^j)`$, where $`\gamma =(1v^2)^{1/2}`$ is the Lorentz factor. The system has to be closed with an equation of state $`p=p(\rho ,e)`$, and here the relation for an ideal gas

$$p=(\mathrm{\Gamma }1)(e\rho )e=\rho +p/(\mathrm{\Gamma }1)$$
(3)
will be considered, where the adiabatic index should be taken as $`\mathrm{\Gamma }=5/3`$ for the mildly relativistic case and as $`\mathrm{\Gamma }=4/3`$ for the ultrarelativistic case $`e\rho `$.

Godunov-type shock-capturing numerical methods developed for classical Euler equations can actually be applied to any multidimensional system of hyperbolic conservation laws of the form

$$\frac{u}{t}+\underset{i=1}{\overset{d}{}}\frac{f^i(u)}{x^i}=0,$$
(4)
where $`u`$ is the vector of conserved variables and $`f^i`$ are their corresponding fluxes, along each direction ($`d`$ is the number of spatial dimensions). Equations (1) and (2) are automatically cast in this form by just defining

$$u(v)=[\rho \gamma ,w\gamma ^2v^j,w\gamma ^2p]^T,$$
(5)
$$f^i(v)=[\rho \gamma v^i,w\gamma ^2v^iv^j+p\delta ^{ij},w\gamma ^2v^i]^T,$$
(6)
where $`v=[\rho ,v^j,p]^T`$ are called primitive variables, and therefore the numerical techniques largely used for Euler HD equations can be applied to RHD too. The hyperbolic nature of Eq. (4) is guaranteed provided the local sound velocity $`c_s`$ is sub-luminal, that is for causal equations of state (Anile 1989). For $`\mathrm{\Gamma }`$-law gases where Eq. (3) holds, the relativistic sound speed is given by

$$c_s^2=\left(\frac{p}{e}\right)_s=\left(\frac{e}{p}\right)_s^1=a^2\left(1+\frac{a^2}{\mathrm{\Gamma }1}\right)^1\mathrm{\Gamma }p/w,$$
(7)
which is obviously always less than unity, where $`sp\rho ^\mathrm{\Gamma }`$ is the specific entropy and $`a^2=\mathrm{\Gamma }p\rho ^1`$ defines the classical sound speed.

In any numerical time advancing routine, primitive variables $`v`$ have to be derived from the conservative ones at least once per time step, and if this is trivial for Euler equations it is not so in the RHD case, in which a numerical nonlinear root-finding technique must be employed. If we define $`W=w\gamma ^2`$ and

$$u=[D,Q^j,E]^T,$$
(8)
where $`D=\rho \gamma `$, $`Q^j=Wv^j`$, and $`E=Wp`$, the system to be inverted can be cast in the single equation for $`\gamma `$:

$$W(\gamma )^2(1\gamma ^2)Q^2=0,$$
(9)
with $`w=\rho +\mathrm{\Gamma }_1pW=D\gamma +\mathrm{\Gamma }_1\gamma ^2p`$, here $`\mathrm{\Gamma }_1=\mathrm{\Gamma }/(\mathrm{\Gamma }1)`$, and $`p=WE`$, so that

$$W(\gamma )=\frac{E\mathrm{\Gamma }_1\gamma ^2D\gamma }{\mathrm{\Gamma }_1\gamma ^21}.$$
(10)
Once the Lorentz factor is found numerically with the requested accuracy, the primitive variables are easily recovered thanks to the relations above.

## 3 A novel ENO-based central scheme

In this section ENO methods and their implementation for component-wise central schemes, together with our specific modifications, will be presented. For a general introduction and review of ENO methods for hyperbolic conservation laws see Shu (1997).

Consider a numerical discretization of Eq. (4), in the one-dimensional case $`d=1`$ to begin with. Given an interval $`[a,b]`$, $`N`$ numerical cells $`I_i=[x_{i1/2},x_{i+1/2}]`$ of equal length $`\mathrm{\Delta }x=(ba)/N`$ can be defined, with cell centers (grid nodes) given by

$$x_i=a+(i1/2)\mathrm{\Delta }x;i=1,\mathrm{},N.$$
(11)
To an order $`r`$ of spatial accuracy, the numerical value of any quantity $`v(x)`$ (at a given time) will be denoted as $`v_i=v(x_i)+O((\mathrm{\Delta }x)^r)`$ at grid points (the so-called point values) or $`v_{i\pm 1/2}=v(x_{i\pm 1/2})+O((\mathrm{\Delta }x)^r)`$ at cell boundaries, where the order of accuracy refers only to smooth regions, in which the larger stencil can be used for interpolation. Moreover, cell averages are defined as

$$\overline{v}_i=\frac{1}{\mathrm{\Delta }x}_{x_{i1/2}}^{x_{i+1/2}}v(x)dx,$$
(12)
and only for schemes up to second order, $`r2`$, they do coincide with point values $`v_i`$.

Most of the shock-capturing schemes evolve the cell averaged conserved quantities $`\overline{u}_i`$ in time, as obtained by integrating Eq. (4) over the cell $`I_i`$. However, in the multidimensional case, the resulting numerical fluxes $`\widehat{f}_i`$, that discretize the physical ones along one direction, have to be averaged along the transverse directions. This implies a truly multidimensional numerical interpolation for high order schemes, which is usually complex and computationally expensive. This approach is generally referred to as the finite volume approximation. Like in LD and in many other ENO schemes, from the works by Shu & Osher (1988; 1989) onwards, we will adopt here finite difference approximations, based on point values. In the semi-discrete formalism (that is retaining continuous time dependency in the spatially discretized quantities), Eq. (4) becomes ($`d=1`$):

$$\frac{\mathrm{d}u_i}{\mathrm{d}t}=\frac{\widehat{f}_{i+1/2}\widehat{f}_{i1/2}}{\mathrm{\Delta }x},$$
(13)
where $`\widehat{f}_{i\pm 1/2}`$ are high-order approximations to the primitives of physical fluxes, that is to say that cell averages of the $`\widehat{f}(x)`$ function must coincide with point values $`f_i`$ of the flux function $`f(x)`$, to the given accuracy ($`\widehat{f}_{i\pm 1/2}f_{i\pm 1/2}`$ up to second order). The extension to the multidimensional case is now straightforward, since interpolations from cell centers to cell boundaries are made separately, dimension by dimension, and the other derivatives are just subtracted from the right hand side of Eq. (13) exactly in the same fashion as for the discretized $`x`$ derivative.

As specified by Shu & Osher (1988), Eq. (13) must be integrated in time by using proper multi-level Runge-Kutta methods corresponding to the high order of spatial accuracy, so here we will always employ the optimal third order TVD algorithm:

$`u^{(1)}`$ $`=`$ $`u^n+\mathrm{\Delta }t[u^n]`$
$`u^{(2)}`$ $`=`$ $`{\displaystyle \frac{3}{4}}u^n+{\displaystyle \frac{1}{4}}u^{(1)}+{\displaystyle \frac{1}{4}}\mathrm{\Delta }t[u^{(1)}]`$ (14)
$`u^{n+1}`$ $`=`$ $`{\displaystyle \frac{1}{3}}u^n+{\displaystyle \frac{2}{3}}u^{(2)}+{\displaystyle \frac{2}{3}}\mathrm{\Delta }t[u^{(2)}]`$
where the superscript $`n`$ indicates the time stepping discretization, $`u^{(1)}`$ and $`u^{(2)}`$ refer to intermediate integration stages, and $`[u]`$ is simply the right hand side of Eq. (13). The explicit time advancing scheme above is stable under the CFL (Courant-Friedrichs-Lewy) condition $`c<1`$, where the Courant coefficient appears in the definition of the maximum time step allowed:

$$\mathrm{\Delta }t=\frac{c}{\mathrm{max}_i(\alpha ^i/\mathrm{\Delta }x^i)},$$
(15)
with $`\alpha ^i`$ being the largest speed of propagation of characteristic waves in the direction $`i`$. In the relativistic case, there is clearly a lower limit for the time step, given in Eq. (15) by simply taking $`\alpha ^i=1`$ along all directions.

In order to complete the description of the scheme, the interpolation techniques and the approximate Riemann solver that defines the numerical inter-cell fluxes $`\widehat{f}_{i+1/2}`$ must be given. The following steps are taken, for every Runge-Kutta sub-cycle and for every direction in the multidimensional case:

1. Primitive variables are recovered from the conservative ones, according to the recipe already given in the previous section (this step is actually taken just once for each sub-cycle):

$$\{u_i\}\{v_i\};i=1,\mathrm{},N.$$
(16)
2. Primitive variables are reconstructed at cell interfaces to give two states, called left and right inter-cell states:

$$\{v_i\}\{v_{i+1/2}^L\},\{v_{i+1/2}^R\};i=0,\mathrm{},N.$$
(17)
The two reconstructions are based on polynomial interpolation over different sets of stencils, one centered on $`x_i`$ to approximate the left state $`v_{i+1/2}^L`$ and the other centered on $`x_{i+1}`$ to approximate the right state $`v_{i+1/2}^R`$. The interpolation is performed separately on each variable $`v(x)`$ according to a Rec routine based on the CENO (Convex ENO, Liu & Osher 1998) technique, as described in Appendix A. Contrary to common ENO schemes, our method uses point values $`v_i`$ to yield point value reconstructed quantities $`v_{i+1/2}`$, instead of starting from cell averages $`\overline{v}_i`$. The order of the reconstruction is $`r=3`$ in smooth regions, but it reduces to linear reconstruction or even to first order (by using minmod-type limiters) near discontinuities. In this latter case $`v_{i+1/2}^L=v_i`$ and $`v_{i+1/2}^R=v_{i+1}`$, for each discontinuous field.
3. At each inter-cell point $`x_{i+1/2}`$ the local Riemann problem must be solved, in some approximate way. The most accurate solution would be given by an exact solver that computes the evolved state from $`v^L`$ and $`v^R`$ and then the correspondent flux, to be finally identified with $`f_{i+1/2}`$. Another possibility is to define an average state $`\stackrel{~}{v}_{i+1/2}`$ and to decompose fluxes according to a linearized problem with the Jacobian $`f/u`$ calculated there (Roe matrix approach). Our choice is to avoid spectral decomposition, which is computationally expensive, and to take an average over the local Riemann fan as in central-type schemes, already mentioned in the introduction. A simple two-speeds Riemann solver is the HLL one, first used in central schemes by Kurganov et al. (2001), that retains an upwind nature in the sense that coincides with $`f^L`$ or $`f^R`$ if the Riemann fan does not cross the inter-cell itself (all the eigenvalues have the same sign). This may be written in the form

$$f^{\mathrm{HLL}}=\frac{\alpha ^+f^L+\alpha ^{}f^R\alpha ^+\alpha ^{}(u^Ru^L)}{\alpha ^++\alpha ^{}},$$
(18)
where all quantities are calculated from the reconstructed values of Eq. (17) by using Eqs. (5) and (6). Here the $`\alpha ^\pm `$ coefficients take into account the highest speeds at the two sides of the Riemann fan, which can be estimated from the maximum and minimum eigenvalue $`\lambda ^\pm `$ of the Jacobians at the left and right states:

$$\alpha ^\pm =\mathrm{max}\{0,\pm \lambda ^\pm (v^L),\pm \lambda ^\pm (v^R)\}.$$
(19)
For relativistic flows, the required eigenvalues are (after splitting the velocity along and perpendicular to the direction of spatial derivation, according to the relativistic rule for addition of velocity vectors):

$$\lambda ^\pm =\frac{v_{}(1c_s^2)\pm c_s\sqrt{(1v^2)(1v_{}^2v_{}^2c_s^2)}}{1v^2c_s^2},$$
(20)
that reduce simply to $`\lambda ^\pm =(v\pm c_s)/(1\pm vc_s)`$ in the one-dimensional case. Note that the maximum and minimum eigenvalues at the two reconstructed states are the only spectral pieces of information required. In this way shocks are handled correctly, whereas contact discontinuities and shear waves, corresponding to intermediate eigenvalues $`\lambda =v_{}`$, can appear somehow smeared, when compared with the results from proper Riemann solvers.

The simplest, smoothest, but also most dissipative numerical flux is the (local) Lax-Friedrichs one, given by

$$f^{\mathrm{LLF}}=\frac{1}{2}[f^L+f^R\alpha (u^Ru^L)],$$
(21)
in which $`\alpha =\mathrm{max}\{\alpha ^+,\alpha ^{}\}`$ and therefore the averaging region is always symmetric with respect to $`x_{i+1/2}`$ (the LLF numerical flux is the prototype of central schemes). Note that the maximum of the set of values $`\{\alpha _{i+1/2}\}`$ yields the CFL-related $`\alpha `$ coefficient appearing in Eq. (15), and actually they coincide for the so-called global LF numerical flux. The calculation of eigenvalues may be avoided completely in a global LF scheme that uses $`\alpha =1`$ everywhere, although this is the most smearing case. Note that Eq. (18) reduces to Eq. (21) when $`\alpha ^+=\alpha ^{}=\alpha `$, that is for symmetric Riemann fans in which $`v_{}=0\lambda ^+=\lambda ^{}`$.
4. From the point values of numerical fluxes, the approximations of their derivatives must be finally calculated (for lower than third order schemes this step can be avoided):

$$\{f_{i+1/2}\}\{\widehat{f}_{i+1/2}\};i=0,\mathrm{},N.$$
(22)
These are the numerical fluxes that actually enter Eq. (13). This last step usually does not appear in other ENO high-order schemes, since the reconstruction is made directly on fluxes and steps 2 and 4 are taken simultaneously (the numerical flux is calculated before reconstruction, using flux splitting methods). However, we have noticed that for central-type schemes like the one presented here, which avoid characteristic decomposition, the proposed approach is more robust and less oscillatory. This final interpolation is again performed separately on each variable with another CENO-based routine (Der, see Appendix A).

## 4 Numerical results

In this section the proposed scheme is validated against typical tests available in the current literature, separated here in one-dimensional tests, essentially shock-tube Riemann problems, and a few 2-D and 3-D experiments. All the simulations have been run on a single PC-Linux (AMD Athlon 1GHz CPU, 512Mb RAM, Pacific-Sierra Vast/f90 compiler, –O2 optimization and single precision calculations) in order to demonstrate that our scheme is not particularly demanding in terms of computational resources. In 3-D, with a grid of $`100^3`$ nodes, typical speeds are of about 2 mins per Runge-Kutta iteration (three internal sub-cycles for the third order scheme), so a simulation like the spherical explosion takes a few hours of computational time. However, the code is fully parallelized with Message Passing Interface (MPI) directives and has been successfully tested on a variety of supercomputers, like CrayT3E, IBM SP3/4 and Beowulf clusters.

In the following tests, we will always use the proposed third order CENO3 as a base scheme, unless specified otherwise. The approximate Riemann solvers employed are the HLL or the LLF, as described in the previous section, whereas the slope limiters used in the linear part of the reconstruction (see Appendix A) are the Monotonized Centered (MC) or the MinMod (MM). All the simulations presented here use the classical value $`\mathrm{\Gamma }=5/3`$ for the adiabatic index.

### 4.1 One-dimensional tests

Shock-capturing numerical schemes, as in the classical hydrodynamical case, must be able to reproduce the discontinuous profiles involved in the so-called shock-tube problems. Given a unitary numerical 1-D pipe of $`N`$ grid points, two constant states $`(\rho ,v,p)`$ are taken on the left ($`0x0.5`$) and on the right ($`0.5<x1`$) with respect of a diaphragm, placed initially at $`x=0.5`$ and then removed. Typical patterns seen in the subsequent evolution are shocks, contact discontinuities (characterized by density jumps not accompanied by discontinuities in normal velocity and pressure) and rarefaction waves. In the relativistic regime these features are qualitatively unchanged, since the structure of the characteristics is the same, but density jumps are not limited by any function of the adiabatic index and rarefaction waves do not yield straight profiles, due to the nonlinear Lorentz transformation formulae.

First we present two relativistic blast wave explosion problems, characterized by an initial static state with temperature and pressure much higher in the region on the left:

$$\text{Blast wave 1:}\{\begin{array}{cc}(\rho ,v,p)^L=\hfill & (10,0,13.3),\hfill \\ (\rho ,v,p)^R=\hfill & (1,0,10^6),\hfill \end{array}$$

$$\text{Blast wave 2:}\{\begin{array}{cc}(\rho ,v,p)^L=\hfill & (1,0,1000),\hfill \\ (\rho ,v,p)^R=\hfill & (1,0,0.01),\hfill \end{array}$$

as in the Donat et al. (1998) paper. The first test is only mildly relativistic, while the second is more severe, with a shock speed corresponding to $`\gamma 6`$. In Fig. 1 and Fig. 2 the two tests are simulated with $`N=400`$ and Courant number $`c=0.5`$, with the base scheme CENO3-HLL-MC. Note the total absence of oscillations and the accuracy in the definition of shocks and rarefaction waves. The contact discontinuity is more smeared, due to the use of a solver that does not take into account that intermediate wave and above all to the reconstruction method, which is the same for all the quantities, thus we cannot steepen (e.g. with the superbee limiter) the contact discontinuity alone, as usually done in characteristics based schemes. In particular, the height of the density peak in Fig. 3 provides a measure of the numerical viscosity of the scheme: the exact value should be around $`10.5`$, while here we find about $`7.3`$, slightly more than what shown in Donat’s paper (in their third order simulation). A smaller value, and a larger spreading, is apparent in Fig. 3, where the most smearing case is tested, that is minmod limiter and global Lax-Friedrichs solver with $`\alpha =1`$, corresponding to $`\lambda ^\pm \pm 1`$ (speed of light) in Eq. (20). However, given the extreme simplicity of the scheme, even this case should not be regarded as completely unusable.

The next example considered is the notorious relativistic shock reflection problem, where an ultrarelativistic cold wind hits a wall, a shock propagates backwards and a static region of relativistically hot gas ($`e\rho c_s^2\mathrm{\Gamma }1`$) is left behind. The numerical box is again $`[0,1]`$ (we use $`N=250`$) and the reflecting wall is placed at $`x=1`$. The physical values employed are:

$$\text{Shock reflection:}(\rho ,v,p)=(1,0.99999,0.01),$$

that corresponds to a Lorentz factor as high as $`\gamma 223`$, about the highest allowed in single precision calculations. This is a very severe test, for the high velocities involved and for the so-called wall heating problem, visible as a dip in the density profile near the reflecting wall. This is a classical numerical artifact, due basically to the implicit numerical viscosity present in every scheme. In Fig. 4 a simulation with second order accuracy (both in time and space), HLL solver and MM limiter is shown for time $`t=0.75`$ ($`c=0.5`$). This is the only case where the third order CENO3 has failed, giving significant postshock oscillations, that can only be reduced by lowering the Courant number or by enhancing numerical viscosity, though they never completely disappear. However, the overheating error in the density is here only $`2.3\%`$, in spite of the second order of accuracy and of our simplified scheme, to be compared with the $`2.5\%`$ value of Donat’s paper in their third order implementation of the celebrated Marquina’s scheme.

As a last 1-D application we present a test proposed in the Dolezal & Wong (1995) paper, which nicely shows the capacity of ENO-based schemes of treating accurately both discontinuous and smooth features occurring close together and at the same time. A shock tube is perturbed in its right-hand state ($`0.5<x1`$) with a density sinusoidal profile:

$$\text{Density perturbation:}\{\begin{array}{cc}(\rho ,v,p)^L=\hfill & (5,0,50),\hfill \\ (\rho ,v,p)^R=\hfill & (2+0.3\mathrm{sin}50x,0,5).\hfill \end{array}$$

The subsequent evolution in time shows the interaction of the blast wave with the density wave (the values are not exactly the same as in Dolezal’s paper, where a nuclear equation of state is used), which enters the expanding heated region and modulates its density plateau. In Fig. 5 the extreme accuracy of our third order CENO3 scheme is shown, by comparing the resulting density profiles at time $`t=0.4`$ in two runs, one with just $`N=200`$ points (diamonds), tested against one with $`N=2000`$ grid points (solid line). Again, by comparing our results to those obtained by the third order characteristics based ENO scheme in Dolezal’s paper, it is apparent that no significant differences arise, in spite of the much less effort involved in our scheme.

### 4.2 Multidimensional tests

Multidimensional relativistic simulations are more difficult than one-dimensional ones because the velocity components are spatially interpolated separately, possibly causing the condition $`v^2<1`$ to fail in ultrarelativistic regimes due to numerical errors in the reconstruction. For this reason in some cases we had to reduce to first order reconstruction, namely when the Lorentz factor reaches $`\gamma =10`$. Note that this does not mean in any sense that this threshold cannot be exceeded, but only that the accuracy is lower in these regions, which usually are located at shock fronts where the order would have been lowered anyway. We have also tested reconstruction on four-velocity components, which are not bounded, but in this way the problem is just shifted to the primitive variables finding routine, thus with no improvement at all.

While it easy to test one-dimensional codes, since Riemann problems can be solved exactly through iterative algorithms, it is not so in more than one dimension, where it is rare to see really quantitative numerical scheme validations. Here we will propose a 2-D Riemann problem, 2-D and 3-D explosions, compared with the correspondent 1-D cylindrical and spherical solutions, and finally one simulation of a relativistic jet, which is not a quantitative test but it is a true astrophysical application.

The two-dimensional counterpart of shock tubes is a square domain divided in four quadrants of constant values at initial time and free evolution for $`t>0`$. In Lax & Liu (1998) all the possible different configurations involving 1-D shocks, 1-D rarefactions waves and 2-D slip lines (contact discontinuities) were studied in detail. In Fig. 6 we show the output (contours of the density logarithm for time $`t=0.4`$) of a situation similar to their configuration 12, where the four boundaries defines two contact discontinuities and two 1-D shocks (symmetric with respect to the main diagonal). Here a relativistic version of this test is proposed, with the following initial settings:

$$\text{Riemann 2-D:}\{\begin{array}{cc}(\rho ,v_x,v_y,p)^{NE}=\hfill & (0.1,0,0,0.01),\hfill \\ (\rho ,v_x,v_y,p)^{NW}=\hfill & (0.1,0.99,0,1),\hfill \\ (\rho ,v_x,v_y,p)^{SW}=\hfill & (0.5,0,0,1),\hfill \\ (\rho ,v_x,v_y,p)^{SE}=\hfill & (0.1,0,0.99,1).\hfill \end{array}$$

Note that we have not taken exact 1-D shocks across the N and E interfaces, and this may be recognized by observing the evolved discontinuities in Fig. 6, converging towards the NE corner with their complete Riemann fans. In the rest of the domain the structure evolves with curved shock fronts and a complicated pattern in the SW quadrant, reminiscent of an oblique jet with bow shock and converging internal shock fronts. The lines in the SW direction with respect to the bow shock are actually due to spurious waves created at the W and S interfaces by the numerical diffusion term in the energy equation (left and right states have a jump in kinetic energy), and cannot be removed if not by using a Roe-type solver. Note that the most diffusive case (LLF solver and MM limiter) was used in the simulation, since other cases resulted more unstable at the same high resolution and with the same Courant number $`c=0.5`$.

Fig. 7 and Fig. 8 refer respectively to the 2-D cylindrically symmetric and 3-D spherically symmetric blast wave explosion problem below:

$$\text{Radial blast wave:}\{\begin{array}{ccc}(\rho ,v_r,p)=\hfill & (1,0,1000)\hfill & r0.4,\hfill \\ (\rho ,v_r,p)=\hfill & (1,0,1)\hfill & r>0.4,\hfill \end{array}$$

where in both cases the runs have been performed in a Cartesian unit box (with reflective right hand side boundary conditions along each direction) and compared with the correspondent 1-D radial solution, obtained with $`N=800`$ grid points in cylindrical or spherical coordinates. The results are satisfactory, since the radial symmetry is well preserved (here the results along the diagonal are shown) and high Lorentz factors are reached without particular problems ($`\gamma >25`$ in the spherically symmetric case). Here the scheme used is CENO3 with HLL solver and MC limiter, with $`c=0.5`$ and $`c=0.3`$ in the 2-D and 3-D cases, respectively.

Finally, as a typical astrophysical test, we simulate the propagation of an axisymmetric jet in 2-D cylindrical coordinates $`(z,r)`$. Note that jet simulations are a very hard test for codes not based on characteristics decomposition, because of usually stronger numerical viscosity at shear layers. The domain is $`0<r<8`$ and $`0<z<20`$, with reflective boundary conditions on the axis $`r=0`$ and simple extrapolation at the other boundaries (except at $`z=0`$ within the jet radius, where initial values are kept constant). At $`t=0`$ a relativistic jet with $`v_z=0.99`$ and density 100 times less than the surroundings (but same pressure) is located at $`r1`$ and $`z1`$

$$\text{Jet:}\{\begin{array}{ccc}(\rho ,v_z,v_r,p)=\hfill & (0.1,0.99,0,0.01)\hfill & r1,z1,\hfill \\ (\rho ,v_z,v_r,p)=\hfill & (10,0,0,0.01)\hfill & \text{outside};\hfill \end{array}$$

density and velocity jumps are actually smoothed in order to reduce spurious transverse waves that appear due to numerical viscosity at the shear layer. The material is injected with a Lorentz factor $`\gamma 7.1`$, corresponding to a relativistic Mach number (e.g. Duncan & Hughes 1994) of $`=\gamma v/\gamma _{c_s}c_s17.9`$. The jet evolution is followed until $`t=40`$, as shown in Fig. 9, where density contours and gray shades in logarithmic scale are presented. The code settings are CENO3-HLL-MC with Courant coefficient $`c=0.25`$, while the resolution employed is $`400\times 160`$, corresponding to 20 grid points per jet radius. Note that the smearing of contact discontinuities, unavoidable in methods not based on characteristics decomposition, is actually small and vortices due to Kelvin-Helmoltz instabilities are nicely defined, as well as the external bow shock, the internal Mach disk and other shocks reflected off the axis. Moreover, notice the absence of the so-called carbuncle problem, usually manifesting as an extended nose in the front of the jet on the axis (e.g. Quirk 1994).

## 5 Conclusions

The central shock-capturing scheme of Londrillo & Del Zanna (2000) is extended to the relativistic case, here for ideal RHD flows and to RMHD in the next paper of the series. Compared to other schemes proposed for relativistic astrophysical problems over the last decade in the literature, our method is extremely simple and efficient, since no eigenvector decomposition and Riemann solvers are involved. In spite of this, due to the high order accuracy achieved in smooth regions and to TVD limiting near discontinuities, provided by the CENO reconstructions employed, our results compete with those obtained by more sophisticated (but less computationally efficient) codes, even in 2-D and 3-D highly relativistic test problems. Different geometries (cylindrical and spherical coordinates) and boundary conditions have also been tested. All the numerical experiments have been run on a simple PC-Linux machine, even the most demanding 3-D test with one million grid points, in order to demonstrate the efficiency of our code. We believe that our simple scheme can be successfully employed in many relativistic simulations of astrophysical interest, either in the fluid case and above all in the most demanding magnetic case, which will be the subject of a forthcoming paper.

###### Acknowledgements.

The authors thank P. Londrillo and R. Bandiera for valuable discussions. This work has been partly supported by the Italian Ministry for University and Research (MIUR) under grants Cofin2000–02–27 and Cofin2001–02–10.

## Appendix A Convex ENO interpolation routines

Given a function $`v(x)`$ with grid point values $`\{v_i\}`$, for $`i=1,\mathrm{},N`$, let us see how non-oscillatory reconstructions (Rec) and second derivatives (Der) are defined in our scheme. These procedures are based on the original CENO reconstruction from Liu & Osher (1998) and they have been described also in LD (we repeat them here for completeness).

In the range $`x_{i1/2}x_ix_{i+1/2}`$, we first consider the two linear polynomials $`L_i(x)`$, based on the stencils $`[x_{i1},x_i]`$ or $`[x_i,x_{i+1}]`$, and we choose their convex combination which is closest to $`v_i`$, that happens to coincide with the TVD limited reconstruction:

$$\stackrel{~}{L}_i(x)=v_i+v_i^{}\left(\frac{xx_i}{\mathrm{\Delta }x}\right).$$
(23)
Here the non-oscillatory undivided first derivative $`v_i^{}`$ is given by the minmod function:

$$v_i^{}=\mathrm{mm}[\mathrm{\Delta }_{}v_i,\mathrm{\Delta }_+v_i],$$
(24)
and $`\mathrm{\Delta }_\pm v_i=\pm (v_{i\pm 1}\pm 1)`$. Thus, to second order of accuracy the CENO technique simply reduces to the usual TVD monotone slope limiting. Although not in the philosophy of the original CENO approach, we have tested other limiters. The most compressive ones, like superbee are too oscillatory in our component-wise framework, while a good compromise between sharpening and stability appears to be provided by the so-called monotonized centered (TVD for $`1\theta 2`$):

$$v_i^{}=\mathrm{mm}[\theta \mathrm{\Delta }_{}v_i,\mathrm{\Delta }_0v_i,\theta \mathrm{\Delta }_+v_i],$$
(25)
where $`\mathrm{\Delta }_0v_i=(\mathrm{\Delta }_+v_i+\mathrm{\Delta }_{}v_i)/2=(v_{i+1}v_{i1})/2`$ and $`\theta =2`$ in order to bias towards central interpolation (since twice the one-sided derivatives might be chosen near discontinuities, the Courant coefficient should be lowered correspondingly for overall stability). In both relations the function $`\mathrm{mm}`$ is defined as

$$\mathrm{mm}(a_1,a_2,\mathrm{})=\{\begin{array}{cc}\mathrm{min}_j\{a_j\}\hfill & a_j0j,\hfill \\ \mathrm{max}_j\{a_j\}\hfill & a_j0j,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}$$
(26)
To achieve an accuracy of third order in smooth regions, we need to calculate the three quadratic polynomials ($`k=1,0,+1`$, $`j=i+k`$):

$$Q_i^k(x)=v_j+\mathrm{\Delta }_0v_j\left(\frac{xx_j}{\mathrm{\Delta }x}\right)+\frac{1}{2}\mathrm{\Delta }_0^2v_j\left(\frac{xx_j}{\mathrm{\Delta }x}\right)^2,$$
(27)
where $`\mathrm{\Delta }_0^2v_j=\mathrm{\Delta }_+\mathrm{\Delta }_{}v_j=\mathrm{\Delta }_{}\mathrm{\Delta }_+v_j=v_{j+1}2v_j+v_{j1}`$, together with the corresponding weighted differences

$$d_i^k=\alpha ^k(Q_i^k\stackrel{~}{L}_i).$$
(28)
Then, we take again their convex combination which is closest to $`\stackrel{~}{L}_i`$ at a given point $`x`$ (this procedure is the generalization of minmod limiting to higher orders) and finally, in smooth regions where all $`d_i^k`$ have the same sign, we take

$$\stackrel{~}{Q}_i=Q_i^{k_0},|d_i^{k_0}|=\mathrm{min}_k(|d_i^k|),$$
(29)
otherwise we exit from the interpolation routine with $`\stackrel{~}{Q}_i=\stackrel{~}{L}_i`$. As for the coefficient $`\theta `$ in the limiter of the lower order, the weights $`\alpha ^k`$ are chosen in order to bias towards central interpolation, and following the original recipe for component-wise CENO schemes we take $`\alpha ^{\pm 1}=1`$ and $`\alpha ^0=0.7`$. In our Rec routine, the left and right states are calculated at the same time and the output is finally

$$v_{i+1/2}^L=\stackrel{~}{Q}_i(x_{i+1/2}),v_{i1/2}^R=\stackrel{~}{Q}_i(x_{i1/2}).$$
(30)
The other CENO routine used in our scheme is Der, which allows one to calculate non-oscillatory second derivatives $`\overline{v}_i^{\prime \prime }`$, needed to transform cell averages into point values, calculated at the same point:

$$\{\overline{v}_i\}\{v_i\},v_i=\overline{v}_i\frac{1}{24}\overline{v}_i^{\prime \prime }.$$
(31)
In this case, the above CENO selection criterion applies with null linear polynomials, $`\stackrel{~}{L}_i=0`$, and quadratic ones given by ($`k=1,0,+1`$, $`j=i+k`$):

$$Q_i^k=\mathrm{\Delta }_0^2\overline{v}_j\overline{v}_{j+1}2\overline{v}_j+\overline{v}_{j1},$$
(32)
then finally $`\overline{v}_i^{\prime \prime }=\stackrel{~}{Q}_i`$ (which is zero, so $`\overline{v}_iv_i`$ as for second order approximation, in the case different signs of the $`Q_i^k`$ terms).


format: text

text:
# Vacuum selection by inflation as the origin of the dark energy

## Abstract

I propose a new mechanism to account for the observed tiny but finite dark energy in terms of a non-Abelian Higgs theory, which has infinitely many perturbative vacua characterized by a winding number, in the framework of inflationary cosmology. Inflation homogenizes field configuration and practically realizes a perturbative vacuum with vanishing winding number, which is expressed by a superposition of eigenstates of the Hamiltonian with different vacuum energy density. As a result, we naturally find a nonvanishing vacuum energy density with fairly large probability, under the assumption that the cosmological constant vanishes in some vacuum state. Since the predicted magnitude of dark energy is exponentially suppressed by the instanton action, we can fit observation without introducing any tiny parameters.

preprint: OU-TAP 178
One of the greatest mysteries of contemporary cosmology is the origin of dark energy or vacuum-like component of cosmic energy indicated by the analysis of type Ia supernovae observations1-3, whose magnitude is some 120 orders smaller than its natural value, the Planck density4. In this essay I propose a novel mechanism to account for the observed tiny but finite dark energy without introducing any small numbers in the framework of inflationary cosmology5 making use of a non-Abelian Higgs theory, which is a part of typical unified theories of elementary interactions, with infinitely many perturbative vacua characterized by winding numbers. In the proposed mechanism inflationary expansion in the early universe plays an essential role in realizing the appropriate quantum state. So let us first recall its cosmological implications.

Inflation stretches preexistent inhomogeneities and yields an extremely large smooth domain in which the current Hubble volume is contained. It also provides a generation mechanism of density fluctuations which originate in the quantum nature of the inflation-driving scalar field6-8. In calculating perturbations, we consider the field to be in the vacuum state and calculate its zero-point quantum fluctuations. This treatment should be an extremely good approximation, if not exact, as long as we focus on perturbation within currently observable Hubble volume. This is because, even if the scalar field was in some excited state initially, subsequent inflation would exponentially dilute the preexistent quanta. Indeed the trace of such an initial excitation could be seen only if an observation was made across an exponentially large region far beyond the current Hubble horizon, which is of course infeasible. Thus for all the practical purposes we can regard that the quantum state of the scalar field is no different from the vacuum state in this context. That is to say, the vacuum state is dynamically selected as a result of cosmological inflation.

Next let us consider a non-Abelian gauge theory with the gauge group $`G`$ which has infinitely many perturbative vacuum states classified by the winding number $`n`$. Most gauge groups used in grand unified theories such as $`G=`$SU($`N`$), SO($`N`$), and Sp($`N`$) with $`N2`$ possess this property. For definiteness, however, we concentrate on the simplest theory with this property, $`G=`$SU(2), below. In these theories perturbative vacuum states with different winding numbers cannot be transformed from each other by a continuous gauge transformation9 and there is an energy barrier between them. Quantum transition, however, is possible and described by an instanton solution which is an Euclidean solution connecting two adjacent perturbative vacua with a finite Euclidean action $`S_0=8\pi ^2/g^2`$ with $`g`$ being the gauge coupling constant10-11. As a result the true vacuum state is expressed by a superposition of perturbative vacua, $`|n`$, as

$$|\theta =\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e^{in\theta }|n,$$
(1)
where $`\theta `$ is a real parameter9,12. One can easily find that this state is a real eigenstate of the Hamiltonian $``$ by calculating the amplitude $`\theta ^{}\left|e^T\right|\theta `$ with the dilute instanton approximation12, namely, adding contributions of $`m`$ instantons and $`\overline{m}`$ anti-instantons keeping in mind that each (anti-)instanton changes the winding number by $`(1)`$1.

$`\theta ^{}\left|e^T\right|\theta ={\displaystyle \underset{n,n^{}}{}}n^{}\left|e^T\right|ne^{in\theta in^{}\theta ^{}}={\displaystyle \underset{n,n^{}}{}}{\displaystyle \underset{m,\overline{m}}{}}{\displaystyle \frac{1}{m!}}{\displaystyle \frac{1}{\overline{m}!}}\left(KVTe^{S_0}\right)^{m+\overline{m}}\delta _{m\overline{m},nn^{}}e^{in\theta in^{}\theta ^{}}`$ (2)
$`={\displaystyle \underset{n,n^{}}{}}e^{in(\theta \theta ^{})}{\displaystyle \underset{m,\overline{m}}{}}{\displaystyle \frac{1}{m!}}{\displaystyle \frac{1}{\overline{m}!}}\left(KVTe^{S_0}\right)^{m+\overline{m}}e^{i(m\overline{m})\theta }=\mathrm{exp}\left(2KVTe^{S_0}\mathrm{cos}\theta \right)\delta \left(\theta \theta ^{}\right),`$ (3)
Here $`KVTe^{S_0}`$ represents contribution of a single instanton or anti-instanton where $`K`$ is a positive constant and $`VT`$ represents spacetime volume.

Then the above equality (2) clearly shows that each $`\theta `$-vacuum $`|\theta `$ has a different energy density than the perturbative vacuum $`|n`$ by $`\mathrm{\Delta }\rho =2Ke^{S_0}\mathrm{cos}\theta .`$ Apparently the $`\theta =0`$ vacuum has the lowest energy, but one cannot conclude that this is the only vacuum state, because the other $`\theta `$-vacua are also stable against gauge-invariant perturbations9,12. In fact, the factor $`K`$ is divergent due to the contribution of arbitrary large instantons in pure gauge theory13. In order to obtain a physical cutoff scale, let us introduce an SU(2) doublet scalar field $`\mathrm{\Phi }`$ with a potential $`V[\mathrm{\Phi }]=\lambda (|\mathrm{\Phi }|^2M^2/2)^2/2`$ following ’t Hooft13. For $`M0`$, instead of the instanton solution of pure gauge theory we should use a constrained instanton solution which also connects adjacent perturbative vacua with a finite Euclidean action14,15. Then the prefactor is given by16 $`K(8\pi /g^2)^4M^4/2`$ and the expectation value of the vacuum energy density in perturbative vacuum $`|n`$ and that in $`\theta `$-vacuum $`|\theta `$ are given by

$`n||n`$ $`=`$ $`\left({\displaystyle \frac{8\pi }{g^2}}\right)^4M^4e^{\frac{8\pi ^2}{g^2}}+\rho _\mathrm{v}(\theta =0),`$ (4)
$`\rho _\mathrm{v}(\theta )\theta ||\theta `$ $`=`$ $`\left({\displaystyle \frac{8\pi }{g^2}}\right)^4M^4e^{\frac{8\pi ^2}{g^2}}(1\mathrm{cos}\theta )+\rho _\mathrm{v}(\theta =0)\rho _\mathrm{d}(1\mathrm{cos}\theta )+\rho _{\mathrm{v0}},`$ (5)
respectively. Here $`\rho _\mathrm{d}(8\pi /g^2)^4M^4e^{\frac{8\pi ^2}{g^2}}`$, and $`\rho _\mathrm{v}(\theta =0)\rho _{\mathrm{v0}}`$ is the vacuum energy density in the state $`|\theta =0`$. Since different vacuum states have different vacuum energy density, they behave differently in the presence of gravity. Hence any good solution to the conventional cosmological constant problem, namely why the vacuum energy density vanishes, should specify at which vacuum state the cosmological constant vanishes. Unfortunately, since there is no completely satisfactory solution to this long-standing problem yet, let us simply normalize the vacuum energy density to vanish in the CP-symmetric $`\theta =0`$ vacuum state and set $`\rho _{\mathrm{v0}}=0`$ for the moment just for definiteness222We point out that the wormhole mechanism17 is an example of proposed solutions to the conventional cosmological constant problem which predicts the vacuum energy density vanishes at $`\theta =0`$ vacuum18.. As will be seen below, however, this choice is not essential and our mechanism can work even if the vacuum energy density vanishes at some other vacuum except at $`|\theta =\pi `$.

Now let us consider what quantum state is chosen by cosmological inflation in this non-Abelian-Higgs system. Let us assume the energy scale of the Higgs field, $`M`$, is much larger than the Hubble parameter during inflation, $`H`$, which is shown to be the case later. As in the case of singlet scalar field, sufficiently long inflation would homogenize the Higgs field configuration over an exponentially large domain with $`|\mathrm{\Phi }|=M/\sqrt{2}`$. This in turn means that the state with vanishing winding number $`|n=0`$ is practically realized as long as we concentrate on the scales within the current Hubble volume, because this is the only state with homogeneous scalar field configuration among many possible vacuum states, perturbative $`|n`$ or real $`|\theta `$. In terms of the energy eigenstates this state is expanded as

$$|n=0=_0^{2\pi }\frac{d\theta }{2\pi }|\theta ,$$
(6)
which consists of superposition of all possible real vacuum state $`|\theta `$ with equal weight.

While we can calculate the expectation value of vacuum energy density in this state as (4), it is more appropriate to discuss the probability distribution function (PDF) of vacuum energy because this state consists of superposition of energy eigenstates with different energy density. Using (6) we can easily calculate the PDF of vacuum energy density, $`P(\rho _\mathrm{v})`$, in this state as

$`P(\rho _\mathrm{v})`$ $`=`$ $`\delta \left(\rho _\mathrm{v}\rho _\mathrm{v}(\theta )\right)_\theta ={\displaystyle _0^{2\pi }}{\displaystyle \frac{d\theta }{2\pi }}\delta \left(\rho _\mathrm{v}\rho _\mathrm{d}(1\mathrm{cos}\theta )\right)`$ (7)
$`=`$ $`\{\begin{array}{cc}\frac{1}{\pi \left(2\rho _\mathrm{v}\rho _\mathrm{d}\rho _\mathrm{v}^2\right)^{1/2}}& \mathrm{for}0\rho _\mathrm{v}2\rho _\mathrm{d},\hfill \\ 0& \mathrm{otherwise}.\hfill \end{array}`$ (10)
We can also calculate the cumulative probability distribution function to find $`\rho _\mathrm{v}\rho `$ as

$$P(\rho _\mathrm{v}\rho )=_\rho ^{2\rho _\mathrm{d}}P(\rho )𝑑\rho =\frac{1}{2}\frac{1}{\pi }\mathrm{arcsin}\left(\frac{\rho }{\rho _\mathrm{d}}1\right),$$
(11)
for $`0\rho 2\rho _\mathrm{d}`$. Figures 1 and 2 depict these PDFs (10) and (11), respectively. As is seen there $`P(\rho _\mathrm{v})`$ is sharply peaked at $`\rho _\mathrm{v}=0`$ corresponding to $`\theta =0`$ and $`\rho _\mathrm{v}=2\rho _\mathrm{d}`$ corresponding to $`\theta =\pi `$ and the probability to find $`\rho _\mathrm{v}2\rho _\mathrm{d}`$ is fairly large as seen in Figure 2.

Thus, if the perturbative vacuum $`|n=0`$ is stable in the cosmological time scale, the PDF of the vacuum energy density is still given by (10) today in our universe which experienced inflation. Hence it is by no means surprising that we observe a finite dark energy around $`2\rho _\mathrm{d}`$ today.

We can easily match the predicted value, $`2\rho _\mathrm{d}`$, with the observed one, $`10^{120}M_G^4`$, by appropriately choosing $`M`$ and $`g`$, where $`M_G`$ is the reduced Planck scale. We also demand that the tunneling rate from $`|n=0`$, $`\mathrm{\Gamma }Ke^{\frac{16\pi ^2}{g^2}}`$, is small enough to guarantee that there is no transition within the current Hubble volume, $`H_0^3`$, in the cosmic age, $`\mathrm{\Gamma }H_0^41`$. Then we obtain

$$\pi /\alpha +2\mathrm{ln}\alpha =60\mathrm{ln}10+2\mathrm{ln}\left(M/M_G\right),\text{and}M\alpha M_G,$$
(12)
where $`\alpha g^2/(4\pi )`$ is the coupling strength at the energy scale $`M/\sqrt{2}`$. If the above inequality is marginally satisfied, we find $`\alpha =1/44.4`$ and $`M=5\times 10^{16}`$GeV. If, on the other hand, we take $`M=M_G`$ so that the cutoff scale of instanton is identical to the presumed field-theory cutoff, we find $`\alpha =1/47.`$

Several comments are in order. (i) Since the scale of the Higgs field $`M`$ turns out to be much larger than the Hubble parameter during inflation, which is constrained as19-21 $`H/(2\pi )4\times 10^{13}`$GeV, quantum fluctuation does not affect the realization of the state $`|n=0`$. (ii) Since the reheat temperature after inflation is much smaller than $`M`$, there is no thermal transition to destabilize the state $`|n=0`$ after inflation. (iii) We may find CP nonconserving value, $`\theta 0,\pi `$. But this would not cause phenomenological troubles because our $`\theta `$ is not directly related with that of strong interaction. (iv) It is not mandatory to assume that the vacuum energy density vanishes at the state $`|\theta =0`$. One may take other vacuum state to normalize the origin of the vacuum energy density. For example, if one assumes that $`\rho _\mathrm{v}`$ vanishes at $`|n`$ or $`|\theta =\pi /2`$, $`P(\rho _\mathrm{v})`$ will be peaked at $`\rho _\mathrm{v}=\pm \rho _\mathrm{d}`$ instead, which we may interpret that we live in a universe with $`\rho _\mathrm{v}\rho _\mathrm{d}`$. Then the values of the model parameters should be slightly changed accordingly.

In conclusion, we have reached a striking conclusion that the observed tiny dark energy can be realized by the interplay of high energy physics between grand unification and Planck scales and inflationary cosmology. Furthermore by virtue of the exponential dependence of the dark energy on the action of instanton, we can reproduce its observed value without introducing any tiny numbers.

References

1. S. Perlmutter et al., Astrophys. J. 517, 565 (1999).
2. B.P. Schmidt et al., Astrophys. J. 507, 46 (1998).
3. A.G. Riess et al., Astron. J. 116, 1009 (1998).
4. S. Weinberg, Rev. Mod. Phys. 61, 1 (1988).
5. For a review of inflation, see, e.g. A.D. Linde, Particle Physics and Inflationary Cosmology (Harwood, Chur, Switzerland, 1990).
6. S.W. Hawking, Phys. Lett. 115B, 295(1982).
7. A.A. Starobinsky, Phys. Lett. 117B, 175(1982).
8. A.H. Guth and S-Y. Pi, Phys. Rev. Lett. 49, 1110(1982).
9. R. Jackiw and C. Rebbi, Phys. Rev. Lett. 37, 172 (1976).
10. A.M. Polyakov, Phys. Lett. 59B, 82 (1975).
11. A.A. Belavin, A.M. Polyakov, A.S. Schvartz, and Yu.S. Tyupkin, Phys. Lett. 59B, 85 (1975).
12. C.G. Callan Jr, R.F. Dashen, and D.J. Gross, Phys. Lett. 63B, 334 (1976).
13. G. ’t Hooft, Phys. Rev. Lett. 37, 8 (1976); Phys. Rev. D14, 3432 (1976); D18, 2199 (1978)(E).
14. I. Affleck, Nucl. Phys. B 191, 429 (1981).
15. M. Nielsen and N.K. Nielsen, Phys. Rev. D61, 105020 (2000).
16. J. Yokoyama, Phys. Rev. Lett. 88, 151302 (2002).
17. S.R. Coleman, Nucl. Phys. B 310, 643 (1988).
18. H.B. Nielsen and M. Ninomiya, Phys. Rev. Lett. 62, 1429 (1989).
19. V.A. Rubakov, M.V. Sazin, and A.V. Veryaskin, Phys. Lett. 115B, 189(1982)
20. L.F. Abbott and M.B. Wise, Nucl. Phys. B 244, 541 (1984).
21. C.L. Bennet etal., Astrophys. J. Lett. 464, L1(1996).


text:
# 1 Introduction

## 1 Introduction

Recently, classical time dependent solutions in string theory describing the decay process of an unstable D-brane was constructed . The decay of an unstable D-brane is caused by the tachyon field on the brane which is rolling down toward the minimum of the potential. In particular, it was shown in  that the rolling of tachyon makes the unstable D-brane toward the tachyon matter which is pressureless gas with nonzero energy density. These solutions are described by boundary states of the unstable D-brane including boundary perturbation associated with the rolling tachyon. In , it was also shown that the effective action proposed by  with a specific tachyon potential reproduces the following two properties of the tachyon matter: 1) absence of plane wave solutions, and 2) exponential fall off of the pressure. However, since this action with the potential has not been derived from any off-shell string theory, the validity of the analysis based on this action is not clear.

On the other hand, the action of boundary string field theory (BSFT) for bosonic string theory  successfully describes some dynamics of the tachyon . In the superstring case, the BSFT actions were obtained for non-BPS D-branes  and D-brane- anti D-brane pairs . These actions exactly describe lower dimensional D-branes as solitons and are expected to satisfy the first property above because of their forms of tachyon potentials. Moreover, we know the coupling of massless closed string modes, including the RR-fields , to the fields in the unstable D-brane in BSFT. Therefore it is both reasonable and interesting to consider the rolling of the tachyon in BSFT.

In this paper, we analyse the classical decay process of the unstable D-brane using the BSFT action. More precisely, we will use the usual BSFT action which is obtained by neglecting $`T`$ and higher derivatives. This action is exact as long as we use the linear tachyon profile $`T=v+u_\mu x^\mu `$. 333In this paper, we only consider the superstring.

At first sight, one might think that the rolling tachyon will hit the singularity within finite time, since the geodesic distance between the points $`T=0`$ and $`T=\mathrm{}`$ in the space of tachyon field $`T`$ is finite . However, as we will show in the following section, the time dependent solutions of the action approach $`T\pm x^0`$ asymptotically and keep rolling forever. The solution with the asymptotic behavior such that $`Tx^0`$ as $`x^0\pm \mathrm{}`$ can be thought of as the space-like brane considered in .

We will also show that the pressure of the system will rapidly converge to zero while the energy density is conserved. Hence the system will become the tachyon matter which is consistent with the Sen’s result.

Unfortunately, since our solution is slightly away from the linear profile, where the BSFT action is reliable, some of our results may receive corrections. We will discuss this issue in section 4.

One of the interesting subject for the rolling tachyon of unstable D-brane systems is the relevance in cosmology . We will also consider the system coupled with gravity, and examine cosmological evolution using the BSFT action.

The rest of the paper is organized as follows. In section 2, we analyse the BSFT action for a non-BPS D-brane, and obtain the asymptotic solution of the rolling tachyon. We show that the system will become the tachyon matter in the large time limit. In section 3, we couple the system to gravity and examine the cosmological evolution. Section 4 is devoted to conclusion and discussion.

## 2 Tachyon matter in BSFT

In this paper we consider the non-BPS D$`p`$-brane in type II string theory. The effective theory of non-BPS D$`p`$-brane is described by the BSFT action

$`S=\stackrel{~}{T}_p{\displaystyle d^{p+1}xe^{\frac{T^2}{4}}\left(\frac{\alpha ^{}}{2}\eta _{\mu \nu }_\mu T_\nu T\right)},`$ (2.1)
where

$`(z)=z{\displaystyle \frac{4^z}{2}}{\displaystyle \frac{\mathrm{\Gamma }(z)^2}{\mathrm{\Gamma }(2z)}}={\displaystyle \frac{\sqrt{\pi }\mathrm{\Gamma }(z+1)}{\mathrm{\Gamma }(z+1/2)}}`$ (2.2)
and $`\stackrel{~}{T}_p`$ is the tension of the non-BPS D$`p`$-brane. Here we only consider the tachyon field and set the other fields to zero for simplicity. This action is considered to be exact if we set $`T=v+u_\mu x^\mu `$, where $`v`$ and $`u_\mu `$ are constants. Using the expansion

$`(z)=1+2\mathrm{log}2z+𝒪(z^2),`$ (2.3)
we can read the potential and kinetic term for the tachyon

$`S=\stackrel{~}{T}_p{\displaystyle d^{p+1}xe^{\frac{T^2}{4}}\left(1+\mathrm{log}2\alpha ^{}_\mu T^\mu T+𝒪((_\mu T^\mu T)^2)\right)}.`$ (2.4)
Next, we compute the Hamiltonian as in . In this paper, we focus on the time dependence of the tachyon field, and mainly assume $`_iT=0`$ for $`i=1,\mathrm{},p`$. Then the momentum conjugate to $`T`$ is

$`\mathrm{\Pi }={\displaystyle \frac{\delta S}{\delta \dot{T}}}=\stackrel{~}{T}_pe^{\frac{T^2}{4}}\alpha ^{}\dot{T}^{}(z),`$ (2.5)
where $`\dot{T}`$ stands for time derivative of $`T`$ and $`z=\frac{\alpha ^{}}{2}\dot{T}^2`$. The Hamiltonian density is

$`=T_{00}=\stackrel{~}{T}_pe^{\frac{T^2}{4}}D(z),`$ (2.6)
where

$`D(z)=(z)2z^{}(z)=(z)\left(1+2\mathrm{log}4z+4z(\psi (z)\psi (2z))\right).`$ (2.7)
where $`\psi (z)`$ is the poly-gamma function defined as

$`\psi (z){\displaystyle \frac{d}{dz}}\mathrm{log}\mathrm{\Gamma }(z)={\displaystyle \frac{\mathrm{\Gamma }^{}(z)}{\mathrm{\Gamma }(z)}}.`$ (2.8)
The behavior of the functions $`D(z)`$ and $`(z)`$ for $`1<z0`$ is depicted in Figure 2 and Figure 2

We would like to consider a solution of the equations of motion which represent the tachyon rolling down to the bottom of the potential. So, we require $`T^2\mathrm{}`$ at $`x^0\mathrm{}`$. Then, (2.6) implies that $`D(z)`$ should become infinity in the limit $`x^0\mathrm{}`$, since the Hamiltonian $`H=𝑑x^p`$ is conserved. Therefore, it is important to examine the singularity of $`D(z)`$ to know the asymptotic behavior of the solution.

Now, we summarize some properties of $`\mathrm{\Gamma }(z)`$. It is holomorphic except $`z=n`$ with $`n=0,1,\mathrm{}`$, and has a pole $`\mathrm{\Gamma }(z)\frac{(1)^n}{n!}\frac{1}{z+n}`$ near $`z=n`$. And $`1/\mathrm{\Gamma }(z)`$ is holomorphic in $`|z|<\mathrm{}`$. The poly-gamma function $`\psi (z)`$ has singularities $`\psi (z)\frac{1}{z+n}`$ at $`z=n`$ with $`n=0,1,\mathrm{}`$. Using these properties, it is easy to see that $`(z)`$ and $`D(z)`$ behaves as

$`(z)`$ $``$ $`n4^n{\displaystyle \frac{(2n)!}{(n!)^2}}{\displaystyle \frac{1}{z+n}},`$ (2.9)
$`D(z)`$ $``$ $`{\displaystyle \frac{2n^2}{4^n}}{\displaystyle \frac{(2n)!}{(n!)^2}}{\displaystyle \frac{1}{(z+n)^2}},`$ (2.10)
near $`z=n`$ where $`n=1,2,\mathrm{}`$444 Note that $`z=0`$ is a regular point for $`(z)`$ and $`D(z)`$, as we have already seen in (2.3)..

The fact that $`D(z)`$ has singularities for finite $`z`$ is crucial for our consideration. Indeed, if $`D(z)`$ had no singularities except $`z=\mathrm{}`$, the tachyon might become infinite within finite time. Actually, we can show that it does occur for the effective action of the Minahan-Zwiebach model  which has $`(z)=1+2\mathrm{log}2z`$, i.e.

$`S=\stackrel{~}{T}_p{\displaystyle d^{p+1}xe^{\frac{T^2}{4}}\left(1+\mathrm{log}2\alpha ^{}_\mu T^\mu T\right)},`$ (2.11)
which is obtained by dropping $`𝒪((_\mu T^\mu T)^2)`$ terms in (2.4). In this case, $`H=\stackrel{~}{T}_pe^{\frac{T^2}{4}}(1+\alpha ^{}\mathrm{log}2\dot{T}^2)`$, which means $`\dot{T}\sqrt{\frac{H}{\stackrel{~}{T}_p}e^{\frac{T^2}{4}}1}`$. Then we can see that $`t=^T𝑑T^{}(\frac{H}{\stackrel{~}{T}_p}e^{\frac{T_{}^{}{}_{}{}^{2}}{4}}1)^{\frac{1}{2}}^T𝑑T^{}\sqrt{\frac{\stackrel{~}{T}_p}{H}}e^{\frac{T_{}^{}{}_{}{}^{2}}{8}}`$ and then $`t`$ is finite at $`T\mathrm{}`$.

Let us consider the behavior of $`D(z)`$ near $`z=1`$, which is the nearest singularity from the origin $`z=0`$. From (2.10), we know that the function $`D(z)`$ behaves as

$`D(z){\displaystyle \frac{1}{(z+1)^2}},`$ (2.12)
near $`z=1`$. As explained above, if we set the initial condition for $`z`$ to be in $`1<z0`$, we expect the asymptotic behavior of the solution is $`T\pm x^0`$ in the limit $`x^0\mathrm{}`$ and $`\mathrm{}`$ (we set $`\alpha ^{}=2`$ for notational simplicity).

Then, suppose that the asymptotic form of the solution at $`x^0\mathrm{}`$ is of the form

$`T=a+x^0+ϵ(x^0),`$ (2.13)
where $`ϵ(x^0)`$ is considered to be a small perturbation. The energy conservation condition implies

$`\dot{ϵ}(x^0)e^{\frac{(x^0)^2}{8}},`$ (2.14)
where $`\dot{ϵ}`$ denotes the time derivative of $`ϵ`$. Then, the asymptotic behavior of the pressure can be computed as

$`P=`$ $`=`$ $`\stackrel{~}{T}_pe^{\frac{T^2}{4}}(z)`$ (2.15)
$``$ $`e^{\frac{(x^0)^2}{8}}.`$ (2.16)
Therefore, this analysis also suggest that the system behaves as a pressureless gas for large time, which is consistent with the result in . So, we conclude that the asymptotic solution (2.13) represents the tachyon matter in BSFT.

Note, however, that there is a subtle difference between the two analyses. The pressure near the end of the rolling given in  using the boundary state is given by

$`Pe^{\alpha x^0}`$ (2.17)
which can also be reproduced by using the effective action of  if we take a specific form of the tachyon potential . The sign of the pressure near the end of the evolution in our model is also different from Sen’s one. Our analysis suggests that the pressure approaches zero from the positive side, while it is negative in . There are several possible reasons for the discrepancy. First, the boundary state of  is different from the solution considered here. Then, the asymptotic form of pressure might be depend on the specific forms of the solutions. Second, the massive modes ignored in the analysis might change the asymptotic form since the difference between the solution and the linear profile could be affected by them. There is another possibility that the higher derivative term could change the asymptotic form.

In the above consideration, we only considered $`z`$ which satisfies $`1<z0`$. We can also consider an asymptotic solution such that $`z`$ approaches $`1`$ from the $`z<1`$ side. 555 We thank J. Minahan for pointing out an error in the argument about this solution in the first version of the paper.  This solution has a peculiar property. Since the function $`D(z)`$ increases when $`z`$ approaches $`1`$ from the $`z<1`$ side, the tachyon slows down while it rolls down the potential. The situation is the same for the other singularities $`z=n`$ ($`n=1,2,\mathrm{}`$) of the function $`D(z)`$, and there are asymptotic solutions approaching $`z=n`$. In any cases the pressure converges to zero, and the system end up with the tachyon matter at late time.

In principle, we can construct a full solution connecting its asymptotic behavior $`T\pm x^0`$. There are basically two types of solutions distinguished by the sign of the asymptotic values of $`T`$, namely, time-like kink solutions, which approach $`Tx^0`$ (or $`Tx^0`$) as $`x^0\pm \mathrm{}`$. and bounce solutions $`T|x^0|`$ (or $`T|x^0|`$) as $`x^0\pm \mathrm{}`$. This can be seen from the equation of motion,

$`{\displaystyle \frac{1}{2}}TD(z)+2\ddot{T}D^{}(z)=0.`$ (2.18)
Note that the left hand side is proportional to $`\dot{}\dot{T}^1`$, and hence it is equivalent to the energy conservation condition as long as $`\dot{T}0`$. Since we know $`D(z)>0`$ and $`D^{}(z)<0`$ from Figure 2, (2.18) implies that $`T`$ is accelerated toward the bottom of the the potential. Thus, once the tachyon starts rolling down the potential, it will never stop. Then, suppose that we set the initial condition, such that the tachyon rolls up the potential. Since $`D(z)1`$ and the equality holds only for $`z=0`$, the energy conservation condition implies that the tachyon cannot stop if the energy density is bigger than the tension $`\stackrel{~}{T}_p`$ and get over the potential barrier, while it will stop at the point where $`\stackrel{~}{T}_pe^{\frac{T^2}{4}}`$ is equal to the energy and turn back if the energy density is smaller than the tension $`\stackrel{~}{T}_p`$.

Note that both of these two types of solutions are localized in the time direction. The time-like kink solutions has “charge” 666 This “charge” is not the usual conserved charge. It is defined as the integral of the RR-field over a large sphere surrounding the finite time region of the non-BPS D-brane in the ten dimensional space-time. , which is calculated from the Chern-Simons term for the non-BPS D-brane , $`S_{CS}=\stackrel{~}{T}_p\mathrm{Str}e^{2\pi \alpha ^{}}C`$, where $``$ is the superconnection. Then, we can identify this type of solutions as the space-like branes considered in .

Note that this time-like kink solution cannot be obtained by simple Wick rotation of usual kink solution for the BSFT action. The kink solution is of the form $`T=ux`$ where $`u`$ is a parameter and we have to take the limit $`u\mathrm{}`$ to minimize the action. The Wick rotation trick $`uiu`$ and $`xix^0`$ implies a configuration with $`zu^2\mathrm{}`$. However, there is an upper bound for the speed of the rolling of the tachyon in the time-like kink solutions, which is not compatible with this limit.

## 3 Cosmological Evolution using the BSFT action

Here, we sketch the cosmological evolution by the BSFT action (2.1) following . Although it is difficult to make a precise analysis in the string field theory point of view, since the higher derivative terms are truncated in the action, it is still interesting to see the qualitative nature of the model and compare it with other models for the cosmological evolution driven by the rolling tachyon . (See also .)

We consider the BSFT action (2.1) in 4 dimension, turning on the minimal coupling to gravity, with the usual Einstein-Hilbert action. We assume that the tachyon only depends on time and the metric is the Robertson-Walker metric

$`ds^2=dt^2+a(t)\left({\displaystyle \frac{dr^2}{1Kr^2}}+r^2d\mathrm{\Omega }_2^2\right).`$ (3.1)
The equations of motion are

$`\left({\displaystyle \frac{\dot{a}}{a}}\right)^2`$ $`=`$ $`{\displaystyle \frac{8\pi G}{3}}\rho {\displaystyle \frac{K}{a^2}},`$ (3.2)
$`{\displaystyle \frac{\ddot{a}}{a}}`$ $`=`$ $`{\displaystyle \frac{4\pi G}{3}}(\rho +3P),`$ (3.3)
where $`\rho `$ and $`P`$ are the energy density and pressure given in (2.6) and (2.15), respectively.

$`\rho `$ $`=`$ $`\stackrel{~}{T}_pe^{\frac{T^2}{4}}D(z),`$ (3.4)
$`P`$ $`=`$ $`\stackrel{~}{T}_pe^{\frac{T^2}{4}}(z),`$ (3.5)
where $`z=\dot{T}^2`$.

From (3.2), we can immediately see that $`\dot{a}`$ is always positive for $`K0`$ cases, describing the expanding universe, since the function $`D(z)`$ is positive the weak energy condition $`\rho >0`$ is satisfied.

The equations of motion also imply

$`\dot{\rho }=3{\displaystyle \frac{\dot{a}}{a}}(\rho +P),`$ (3.6)
which is equivalent to the conservation equation for the energy momentum tensor. Inserting (3.4) and (3.5), this equation can be rewritten as

$`{\displaystyle \frac{1}{2}}TD(z)+2\ddot{T}D^{}(z)=3{\displaystyle \frac{\dot{a}}{a}}(D(z)(z))\dot{T}^1.`$ (3.7)
Suppose that we set the initial condition $`T(\overline{t})>0`$ and $`\dot{T}(\overline{t})=0`$ at $`x^0=\overline{t}`$. We can easily see that the right hand side converges to zero in the limit $`\dot{T}0`$. Then, $`D(0)=1`$ and $`D^{}(0)=F^{}(0)2\mathrm{log}2<0`$ imply $`\ddot{T}>0`$ which shows that the tachyon starts rolling down the potential. This argument also implies that the tachyon will not stop rolling once it starts, as we observed in the previous section for the model without gravity.

The ratio of the pressure and the energy density is given by

$`\omega ={\displaystyle \frac{P}{\rho }}={\displaystyle \frac{1}{1+2\mathrm{log}4z+4z(\psi (z)\psi (2z))}},`$ (3.8)
which has maximum near $`z0.7`$ with $`\omega (0.7)0.1`$ and $`\omega (0)=1,\omega (1)=0`$ as depicted in Figure 3.

Hence, the dominant energy condition $`|\omega |1`$ is satisfied. However the strong energy condition

$`\rho +3P0,`$ (3.9)
does not hold for $`|z|0.278\mathrm{}`$. As we can see from the equation (3.3), this implies the acceleration of the scale factor in the initial stage with small $`|z|`$.

In summary, the qualitative behavior of the model is roughly the same as that given in . It would be interesting to examine further to see if it can be considered as a realistic model of inflation or dark matter. (See )

## 4 Conclusion and Discussion

In this paper, we have analysed the classical decay process of unstable D-branes in superstring theory using the BSFT action. We have shown that the solutions of the equations of motion for the tachyon field asymptotically approach to $`T=\pm x^0`$ producing the tachyon matter. We also considered the cosmological evolution driven by the rolling tachyon using the BSFT action.

We used the BSFT action obtained by dropping terms which include second or higher derivatives of the tachyon field. However, in the time dependent solution of the equations of motion for this effective action, the second or higher derivatives of the tachyon field is not precisely zero. Therefore, there could be some corrections for the analysis if we include the terms we dropped.

Let us make some comments about the possible corrections. First, recall that the BSFT action that we used can be trusted for the tachyon of the form $`T=a+ux^0`$. We showed that the tachyon approaches the linear profile $`T\pm x^0`$ in the asymptotic region. This comes from the fact that the function $`D(z)`$ has a pole at $`z=1`$. The function $`D(z)`$ could receive some corrections. However, we expect that even if we include these corrections, it still has a singularity for the tachyon configuration $`T=a\pm x^0`$. This is because the action (2.1) is exact for the linear profile and we know that the action has a singularity for this tachyon configuration. It is quite natural to expect that the singularity coincides for the action and the energy density. Therefore, we strongly believe that there are solutions which approach $`T\pm x^0`$ in the limit $`|x^0|\mathrm{}`$ even if we include the corrections.

On the other hand, we cannot exclude the possibility that there are solutions which does not approach $`T\pm x^0`$ when we include the higher derivative corrections.

The precise form of the solution for the tachyon field seems to be difficult to obtain in this approach. In addition to the complexity of the action, there are some more fundamental difficulties. If we include the higher derivative corrections, the massive modes can no longer be ignored . The BSFT action including these modes has not been obtained because it corresponds to the non-renormalizable terms in the boundary interaction. Therefore, our solution cannot be trusted away from the asymptotic region. However, most of our qualitative arguments given in section 2 and section 3 only require the broad behavior of the functions $`(z)`$ and $`D(z)`$ and so we hope that these analyses successfully catch the essential features of the theory.

We have observed that the asymptotic behavior of the pressure is slightly different from that obtained in . These ambiguities might be the origin of this discrepancy. In order to clarify these issues, it would be important to study the relations between the BSFT approach and the boundary state approach given in .

In this paper we only turned on the tachyon field. The action including both the gauge field strength $`F`$ and the tachyon for a non-BPS D-brane was obtained in a closed form  as

$`S`$ $`=`$ $`\sqrt{2}T_9{\displaystyle d^{10}xe^{\frac{1}{4}T^2}\sqrt{det(\eta _{\mu \nu }+F_{\mu \nu })}\left(\frac{1}{4\pi }𝒢^{\mu \nu }_\mu T_\nu T\right)},`$ (4.1)
where

$`𝒢^{\mu \nu }\left({\displaystyle \frac{1}{1+F}}\right)^{(\mu \nu )}.`$ (4.2)
We can also obtain the action including scalar fields, representing the transverse fluctuation of the D-brane, by T-dualizing the above action. It would be interesting to include the effect of the gauge field and scalar fields in the consideration of the rolling tachyon. See  for the related analysis. In particular, the time evolution of the system with the electric flux  is interesting to study.

As argued in , all the D$`p`$-branes can be constructed from infinitely many unstable D-particles, i.e. non-BPS D0-branes in type IIB theory or D0-brane - anti D0-brane pairs in type IIA theory. Therefore, it is possible to analyse the decay process of the non-BPS D$`p`$-brane in terms of the unstable D-particles. Then, the tachyon matter produced by the decay of the non-BPS D$`p`$-brane can be thought of as the pressureless gas consists of the decay products of the unstable D-particles. It might be interesting to consider the tachyon matter in this point of view.

### Note added

After this work was completed, a related paper  appeared.

### Acknowledgments

We would like to dedicate this article to the memory of Sung-Kil Yang from whom we learned a lot of things. We will not forget his warm encouragements.

S.T. would like to thank J. de Boer for helpful discussions and T. Takayanagi and J. Minahan for useful comments and discussions. S.S. is grateful to T. Okuda and K. Okuyama for useful discussions.