text:
# Brane cosmology with a van der Waals equation of state

## Abstract

The evolution of a Universe confined onto a 3-brane embedded in a five-dimensional space-time is investigated where the cosmological fluid on the brane is modeled by the van der Waals equation of state. It is shown that the Universe on the brane evolves in such a manner that three distinct periods concerning its acceleration field are attained: (a) an initial accelerated epoch where the van der Waals fluid behaves like a scalar field with a negative pressure; (b) a past decelerated period which has two contributions, one of them is related to the van der Waals fluid which behaves like a matter field with a positive pressure, whereas the other contribution comes from a term of the Friedmann equation on the brane which is inversely proportional to the scale factor to the fourth power and can be interpreted as a radiation field, and (c) a present accelerated phase due to a cosmological constant on the brane.

Key words: brane cosmology, van der Waals fluid, acceleration field

According to the cosmological observations one can distinguish three distinct periods for the Universe that are related to its acceleration field. The first period refers to an accelerated epoch dominated by a scalar field where a rapid expansion of the Universe characterizes its inflationary phase. The next one is related to a decelerated phase dominated by matter fields. This period is followed by a return to an accelerated epoch dominated by a cosmological constant or dark energy.

Recently several authors have investigated the evolution of the Universe within the framework of the one-brane model of Randall and Sundrum  where the Universe is confined onto a 3-brane, which is a hyper-surface embedded in a five-dimensional space-time called bulk. For recent reviews on this subject one is referred to the works of Brax and van de Bruck  and Langlois  and the references therein.

In the present work we investigate the evolution of a Universe confined onto the 3-brane embedded in a five-dimensional space-time in order to describe the three distinct periods of the Universe beginning with an accelerated phase passing through a decelerated epoch and returning to an accelerated period. For that end we model the cosmological fluid on the brane by the van der Waals equation of state. The use of the van der Waals equation of state in cosmological problems was first proposed by Capozziello and co-workers , who recognize that this equation of state could describe the transition from a scalar field dominated period to a matter dominated epoch without the need of introducing scalar fields. Recently a model for the Universe as a mixture of a van der Waals fluid with dark energy (modeled as quintessence or Chaplygin gas) was proposed in the work  – within the framework of a four-dimensional space-time theory – in order to describe the transition from the accelerated-decelerated-accelerated periods of the Universe.

In this work we show that the Universe on the brane modeled by the van der Waals equation of state evolves in such a manner that the three distinct periods concerning its acceleration field are attained. The initial accelerated epoch is due to the van der Waals fluid which behaves like a scalar field with a negative pressure. The past decelerated period has two contributions: one of them is related to the van der Waals fluid which behaves like a matter field with a positive pressure, whereas the other comes from a term of the Friedmann equation on the brane which is inversely proportional to the scale factor to the fourth power and can be interpreted as a radiation field. The present accelerated period is due to a cosmological constant on the brane.

In the one-brane model of Randall and Sundrum  the Universe is confined onto a hyper-surface – called brane – which is embedded in a five-dimensional space-time with coordinates $`(t,x^1,x^2,x^3,y)`$, called bulk. The brane is located at $`y=0`$ and the line element which describes a spatially flat, homogeneous and isotropic Universe on the brane is given by (see, for example, the reviews )

$$ds^2=g_{MN}dx^Mdx^N=n(t,y)^2dt^2a(t,y)^2\delta _{ij}dx^idx^jdy^2,$$
(1)
where $`g_{MN}`$ is the metric tensor with signature $`(+,,,,)`$. The two functions $`n(t,y)`$ and $`a(t,y)`$ are determined from the five-dimensional Einstein field equations that read

$$R_{MN}\frac{1}{2}Rg_{MN}+\mathrm{\Lambda }g_{MN}=\kappa ^2T_{MN}.$$
(2)
Above, $`R_{MN}`$ is the five-dimensional Ricci tensor, $`R=R_{}^{M}{}_{M}{}^{}`$ its trace, $`\mathrm{\Lambda }`$ denotes a bulk cosmological constant, $`\kappa ^2=8\pi G_5`$ is related to the five-dimensional gravitational constant $`G_5`$, and $`T_{MN}`$ is the five-dimensional energy-momentum tensor.

The five-dimensional energy-momentum tensor $`T_{MN}`$ is decomposed into a sum of two terms: one refers to the bulk whereas the other is related to the brane, i.e.,

$$T_{}^{M}{}_{N}{}^{}=T_{}^{M}{}_{N}{}^{}|_{\mathrm{bulk}}+T_{}^{M}{}_{N}{}^{}|_{\mathrm{brane}}.$$
(3)
The energy-momentum tensor in the bulk is given by

$$T_{}^{M}{}_{N}{}^{}|_{\mathrm{bulk}}=\mathrm{diag}(\rho _B,p_B,p_B,p_B,\rho _B),$$
(4)
where the energy density $`\rho _B`$ and the pressure $`p_B`$ in the bulk do not depend on the $`y`$ coordinate. Furthermore, the energy-momentum tensor on the brane reads

$$T_{MN}|_{\mathrm{brane}}=\delta (y)(\sigma g_{\mu \nu }+𝒯_{\mu \nu })\delta _{}^{\mu }{}_{M}{}^{}\delta _{}^{\nu }{}_{N}{}^{}.$$
(5)
On the brane there exist two contributions: one is related to the constant tension on the brane $`\sigma `$, whereas the other refers to the energy-momentum tensor of the cosmological fluid $`𝒯_{\mu \nu }`$ which is written as

$$𝒯_{}^{\mu }{}_{\nu }{}^{}=\mathrm{diag}(\rho ,p,p,p).$$
(6)
Above, $`\rho `$ and $`p`$ denote the energy density and the pressure of the cosmological fluid on the brane, respectively.

For the line element given by (1), the components of the Einstein field equations (2) become

$$3\left\{\left(\frac{\dot{a}}{a}\right)^2n^2\left[\frac{a^{\prime \prime }}{a}+\left(\frac{a^{}}{a}\right)^2\right]\right\}\mathrm{\Lambda }n^2=\kappa ^2T_{00},$$
(7)
$$\left\{a^2\left[\left(\frac{a^{}}{a}\right)^2+2\frac{a^{}n^{}}{an}+2\frac{a^{\prime \prime }}{a}+\frac{n^{\prime \prime }}{n}\right]\frac{a^2}{n^2}\left[\left(\frac{\dot{a}}{a}\right)^22\frac{\dot{a}\dot{n}}{an}+2\frac{\ddot{a}}{a}\right]\right\}\delta _{ij}$$

$$+\mathrm{\Lambda }a^2\delta _{ij}=\kappa ^2T_{ij},$$
(8)
$$3\left(\frac{\dot{a}}{a}\frac{n^{}}{n}\frac{\dot{a}^{}}{a}\right)=\kappa ^2T_{0y},$$
(9)
$$3\left(\frac{a^{}}{a}\right)^2+3\frac{a^{}n^{}}{an}\frac{3}{n^2}\left[\left(\frac{\dot{a}}{a}\right)^2\frac{\dot{a}\dot{n}}{an}+\frac{\ddot{a}}{a}\right]+\mathrm{\Lambda }=\kappa ^2T_{yy}.$$
(10)
In the above equations the dot and the prime refer to a differentiation with respect to the coordinates $`t`$ and $`y`$, respectively.

We follow the work  and introduce the abbreviations at $`y=0`$: $`a_0(t)a(t,y=0)`$ and $`n_0(t)n(t,y=0)`$. Hence, on the brane the line element is written as

$$ds^2=n_0(t)^2dt^2a_0(t)^2\delta _{ij}dx^idx^j.$$
(11)
By imposing the gauge condition $`n_0(t)=1`$, the time coordinate becomes the proper time on the brane and $`a_0(t)`$ is identified with the scale factor.

If one assumes that there exists no matter flow along the fifth dimension (see, e.g. )– so that $`T_{0y}=0`$ – the integration of equation (9) leads to

$$n(y,t)=\frac{\dot{a}(t,y)}{\dot{a}_0(t)}.$$
(12)
We integrate equations (7) and (8) over $`ϵy+ϵ`$, and by taking the limit $`ϵ0`$, it follows the so-called junction conditions:

$$\frac{\left[\left[a^{}\right]\right]}{a_0}=\frac{\kappa ^2}{3}(\rho +\sigma ),\left[\left[n^{}\right]\right]=\frac{\kappa ^2}{3}\sigma +\frac{2\kappa ^2}{3}\rho +\kappa ^2p,$$
(13)
where $`\left[\left[f\right]\right]=f(0+)f(0)`$ denotes the jump of the function $`f`$ across the brane.

If we take into account in equation (9) the junction conditions (13) together with the condition of no matter flow along the fifth dimension – i.e., $`T_{0y}=0`$ – we get the conservation equation for the energy density on the brane

$$\dot{\rho }+3\frac{\dot{a}_0}{a_0}(\rho +p)=0.$$
(14)
In order to analyze the system of Einstein field equations (7) – (10) in the bulk, we follow the work  and introduce the function

$$F(t,y)(a^{}a)^2\frac{(\dot{a}a)^2}{n^2}.$$
(15)
In terms of the function $`F(t,y)`$ the 00, $`yy`$ and $`ij`$ \- components in the bulk read:

$$F^{}=\frac{2a^{}a^3}{3}(\mathrm{\Lambda }+\kappa ^2\rho _B),\dot{F}=\frac{2\dot{a}a^3}{3}(\mathrm{\Lambda }+\kappa ^2\rho _B),$$
(16)
$$\left(\frac{F^{}}{a^{}}\right)^{}=2\dot{a}a^2(\mathrm{\Lambda }\kappa ^2p_B),$$
(17)
respectively. Now we differentiate (16)1 with respect to time $`t`$ and (16)2 with respect to the coordinate $`y`$, and get that $`\rho _B`$ does not depend on time $`t`$. Furthermore, from the differentiation with respect to time of the expression $`F^{}/a^{}`$ obtained from (16)1 it follows an equation when compared with (17) implies that $`p_B=\rho _B`$. Finally, the integration of (16)1 with respect to the coordinate $`y`$, leads to

$$F=(a^{}a)^2\frac{(\dot{a}a)^2}{n^2}=\frac{a^4}{6}(\mathrm{\Lambda }+\kappa ^2\rho _B)+C.$$
(18)
Above, $`C`$ is a constant since $`\rho _B`$ is time independent.

The Friedmann equation on the brane is obtained from (18) by considering the limit $`y0`$ and the junction conditions (13), yielding

$$\left(\frac{\dot{a}_0}{a_0}\right)^2=\frac{\kappa ^4}{18}\sigma \rho +\frac{\kappa ^4}{36}\rho ^2+\frac{1}{36}(\kappa ^4\sigma ^2+6\mathrm{\Lambda }+6\kappa ^2\rho _B)+\frac{C}{a_0^4}.$$
(19)
From the Friedmann equation one can determine the acceleration equation on the brane. Indeed, the differentiation of (19) with respect to time, yields

$$\frac{\ddot{a}_0}{a_0}=\frac{\kappa ^4}{36}\sigma (\rho +3p)\frac{\kappa ^4}{36}\rho (2\rho +3p)+\frac{1}{36}(\kappa ^4\sigma ^2+6\mathrm{\Lambda }+6\kappa ^2\rho _B)\frac{C}{a_0^4},$$
(20)
thanks to the equation (14) for the energy density on the brane. This equation can be obtained also from equation (10) by using the junction conditions (13).

Now we have a system of differential equations for the determination of the energy density $`\rho (t)`$ and of the scale factor $`a_0(t)`$ on the brane which is composed by the conservation equation for the energy density (14) and by the Friedmann equation (19) – or by the acceleration equation (20). In order to find the time evolution of the energy density and of the scale factor from the system of differential equations one has to prescribe initial conditions for these fields and to close the system by choosing an equation of state which relates the pressure of the cosmological fluid to the energy density on the brane, i.e, $`p=p(\rho )`$.

For the determination of the solution of the system of differential equations described above we introduce the dimensionless quantities

$$t\frac{\kappa ^2\sigma }{\sqrt{18}}t,\rho \frac{\rho }{\sigma },p\frac{p}{\sigma },a_0\frac{a_0(t)}{a_0(0)},$$
(21)
and write the conservation equation for the energy density (14), the Friedmann (19) and the acceleration (20) equations in terms of the dimensionless quantities as

$$\dot{\rho }+3\frac{\dot{a}_0}{a_0}(\rho +p)=0,$$
(22)
$$\left(\frac{\dot{a}_0}{a_0}\right)^2=\rho +\frac{\rho ^2}{2}+\frac{\lambda }{3}+\frac{\chi }{a_0^4},$$
(23)
$$\frac{\ddot{a}_0}{a_0}=\frac{1}{2}(\rho +3p)\frac{\rho }{2}(2\rho +3p)+\frac{\lambda }{3}\frac{\chi }{a_0^4}.$$
(24)
In equations (23) and (24) the following abbreviations were introduced

$$\lambda =9\left[\frac{\mathrm{\Lambda }}{\kappa ^4\sigma ^2}+\frac{\rho _B}{\kappa ^2\sigma ^2}+\frac{1}{6}\right],\chi =\frac{18C}{[a_0(0)\kappa ]^4\sigma ^2}.$$
(25)
The constant $`\lambda `$ in the Friedmann (23) and acceleration (24) equations can be interpreted as a cosmological constant on the brane. Moreover, it is noteworthy to call attention to the fact that equations (23) and (24) can be written as

$$\left(\frac{\dot{a}_0}{a_0}\right)^2=\rho +\rho _{\mathrm{rad}}+\frac{\rho ^2}{2}+\frac{\lambda }{3},$$
(26)
$$\frac{\ddot{a}_0}{a_0}=\frac{1}{2}(\rho +\rho _{\mathrm{rad}}+3p+3p_{\mathrm{rad}})\frac{\rho }{2}(2\rho +3p)+\frac{\lambda }{3},$$
(27)
thanks to the well-known relationships valid in the four-dimensional case where the energy density of the radiation field $`\rho _{\mathrm{rad}}`$ scales as $`1/a_0^4`$ and the radiation pressure is related to the energy density by $`p_{\mathrm{rad}}=\rho _{\mathrm{rad}}/3`$.

From now on we shall analyze the system of differential equations consisting of the conservation equation for the energy density (22) and of the acceleration equation (24). In order to close this system of differential equations we have to choose an equation of state that relates the pressure to the energy density of the cosmological fluid on the brane. In the four-dimensional case one normally uses a barotropic equation of state $`p=w\rho `$ with $`0w1`$ to represent a radiation or matter dominated Universe. The barotropic equation of state with $`1w0`$ is also used to represent the scalar fields: the inflaton in the inflationary period of the Universe, and the quintessence in the dark energy dominated epoch of the Universe. Other equations of state are also used to model the cosmological fluid, namely the Chaplygin equation of state (see e.g. the work  and the references therein) which can describe a transition from a matter dominated period to a cosmological constant dominated epoch and the van der Waals equation of state which can simulate the transition from an inflationary period to a matter field dominated epoch . Here we interested in the description of an accelerated period followed by a decelerated epoch, hence we shall use the van der Waals equation of state

$$p=\frac{8w\rho }{3\rho }3\rho ^2,$$
(28)
where the parameter $`w`$ could be identified with the coefficient of proportionality in the barotropic formula, since for small values of the energy density $`pw\rho `$. In classical thermodynamics (see e.g. ) the equation (28) is a reduced van der Waals equation of state where the free parameter $`w`$ is connected with a dimensionless temperature. Capozziello and co-workers  have used a non-reduced form of the van der Waals equation which is characterized by three free parameters instead of only one. In their works  they also write a reduced form of the van der Waals equation with only one free parameter which is not similar to the reduced van der Waals equation of state normally found in the literature, since the free parameter multiplies also the term proportional to $`\rho ^2`$.

Apart from the initial conditions for the fields of energy density $`\rho (0)`$, scale factor $`a_0(0)`$ and velocity $`\dot{a}_0(0)`$, one has to specify values for the parameters $`\lambda `$, $`\chi `$ and $`w`$ in order to find a solution of the system of differential equations (22) and (24) which is closed by the van der Waals equation of state (28). Here we have specified the initial conditions (by adjusting clocks): $`\rho (0)=1`$ for the energy density, $`a_0(0)=1`$ for the scale factor and $`\dot{a}(0)=\sqrt{3/2+\lambda /3+\chi }`$ for the velocity field, which is a consequence of the Friedmann equation (23). Furthermore, in order to plot the time evolution of the fields in figures 1 and 2 we have chosen the following values for the parameters $`\lambda =0.03`$, $`\chi =0.5`$, and three different values for $`w`$, namely, $`w=0.51`$, $`w=0.501`$, and $`w=0.5001`$. Later on we shall discuss how the changes in the parameters $`w`$, $`\lambda `$ and $`\chi `$ affect the solutions of the system of differential equations.

In figure 1 it is plotted the acceleration field $`\ddot{a}_0(t)`$ as function of time $`t`$ for $`w=0.51`$ (straight line), $`w=0.501`$ (dotted line), and $`w=0.5001`$ (dashed line). We infer from this figure that the acceleration field for these three values of $`w`$ describes the distinct accelerated-decelerated-accelerated phases of the Universe. In the first period, which refers to a scalar field dominated epoch, the positive acceleration grows up to a maximum value followed by a decay towards zero. The next phase is related to a matter dominated period where the acceleration is always negative and decays to a maximum negative value followed by a growth towards zero. The third period is a cosmological constant dominated epoch where the acceleration field assumes a positive value. The energy density $`\rho `$ and pressure $`p`$ fields are plotted in figure 2 as function of time $`t`$ for $`w=0.51`$ (straight line), $`w=0.501`$ (dotted line), and $`w=0.5001`$ (dashed line). We conclude from this figure that for the three values of $`w`$ the energy density field decays with time whereas the pressure field has two distinct behaviors. At the beginning the pressure is negative so that the van der Waals fluid behaves like a scalar field and it is the responsible for the initial accelerated period. At later times the pressure becomes positive so that the van der Waals fluid behaves like a matter field and it is partially responsible for the decelerated period, since the radiation field which scales as $`1/a_0^4`$ also contributes to the decelerated phase of the Universe. When $`w`$ approaches the value of $`w=0.5`$, we can also infer from the figures that the pressure field of the van der Waals fluid at the beginning behaves like an inflaton with an equation of state given by $`p=\rho `$ and it is the responsible for an increase of the early acceleration.

We have chosen three nearby values for $`w`$ in order to show how the acceleration, energy density and pressure fields behave for small changes of $`w`$. We proceed now to discuss how significant changes in the values of $`w`$ have influence on the behavior of these fields for fixed values of $`\lambda `$ and $`\chi `$, here $`\lambda =0.03`$ and $`\chi =0.5`$. For $`w=0.5`$ the energy density remains constant so that the van der Waals fluid behaves like an inflaton with an equation of state given by $`p=\rho `$, the acceleration field grows exponentially and in this case there exists only an accelerated phase for the Universe. For values of $`w<0.5`$ the energy density grows with time, and we infer that this behavior does not represent a physical solution for the evolution of the Universe. For values of $`0.5<w<0.58`$ the initial accelerated period decreases, since the negative part of pressure of the van der Waals fluid decreases, whereas its positive part increases. Hence in the interval $`0.5<w<0.58`$ the three periods accelerated-decelerated-accelerated are present in the evolution of the Universe. These three periods for the Universe can be found for a larger interval for $`w`$ by decreasing the amount of radiation, i.e., by decreasing the value of the constant $`\chi `$. For values of $`w>0.58`$ (and $`\lambda =0.03`$, $`\chi =0.5`$) the acceleration field evolves from a matter dominated Universe where $`\ddot{a}<0`$ to a cosmological constant dominated Universe where $`\ddot{a}>0`$. In this last case there exists no inflationary period.

Let us comment on the behavior of the fields when we change the values of the parameters $`\chi `$ and $`\lambda `$. By increasing the value of $`\chi `$ there exists a more pronounced predominance of the radiation field which scales as $`1/a_0^4`$ and the deceleration period begins at earlier times. By decreasing the value of $`\lambda `$, which is connected to the cosmological constant on the brane, the present accelerated period begins at later times. Moreover, the limit $`\lambda 0`$ leads to an initial accelerated period followed by a decelerated period without a present accelerated period, indicating that $`\lambda `$ plays the role of the dark energy. We could also obtain the present accelerated period – which is characterized by a dark energy dominated Universe – by following the same methodology of the work  and instead of introducing a cosmological constant which represents the dark energy, we could model the dark energy by an equation of state like the Chaplygin or the quintessence equations of state.

As a final comment we have also investigated the solution of the system of differential equations (22) and (24) closed by a barotropic equation of state for the cosmological fluid on the brane, i.e., $`p=w\rho `$ with $`1w1`$. The results we have found for the acceleration field by changing the values of $`w`$ describe either an accelerated phase, or a decelerated epoch followed by an accelerated period of the Universe, i.e., with the barotropic equation of state it was not possible to describe the three phases accelerated-decelerated-accelerated found here by using the van der Waals equation of state.


format: text

text:
# Nitrogen abundances in Planet-harbouring stars

## 1 Introduction

Knowledge of CNO abundances in main sequence stars is clearly important because of the key roles of the former in the chain of nucleosynthesis. They can give us information about the history of element production in different types of stars. Solar-type stars can be very informative as tracers of the chemical and dynamical evolution of the Galaxy (e.g. Pagel & Edmunds Pag81 (1981); Edvardsson et al. Edv93 (1993); Bensby et al. Ben03 (2003)) as they spend between $`10^9`$ to several times $`10^{10}`$ years on the main sequence.

The abundant isotope $`{}_{}{}^{14}\mathrm{N}_{}^{}`$ is synthesized through the CNO cycles in a hydrogen-burning shell of the stellar interior. With regard to the chemistry of planet-harbouring stars, the study of the nitrogen abundances of solar-type dwarfs with and without known planets represents an important key to checking the “self-enrichment” scenario. This hypothesis attributes the origin of the overabundance of metals to the accretion of large amounts of metal-rich, H- and He-depleted rocky planetesimal materials on to the star.

On the one hand, by studying chemical abundances and planetary orbital parameters, several authors reject the possibility of pollution as the main source of the metallicity enhancement observed in planet host stars and propose instead that the metal-rich nature of stars with planets has most probably a “primordial” source in the high metallicity of the primordial cloud (Santos et al. San01 (2001), San02 (2002), San03a (2003)a). The effect on metal abundances of contamination in the outer convective zone by accreted planetary material in FGK main sequence stars has been discussed in Pinsonneault et al. (Pin01 (2001)) and the same conclusion was reached. However, the possibility that pollution may occur has not been excluded. Israelian et al. (Isr01 (2001), 2003a ) have found evidence of the infall of a planet into the planet host star HD 82943. $`{}_{}{}^{6}\mathrm{Li}_{}^{}`$ and $`{}_{}{}^{7}\mathrm{Li}_{}^{}`$ may provide an independent test to confirm this scenario (Israelian et al. 2003a , 2003b ).

On the other hand, the “self-enrichment” scenario would imply that volatile element (with low condensation temperature $`T_\mathrm{C}`$) abundances do not show overabundance in the same way that refractories do, since the elements of the former group are deficient in accreted materials relative to the latter. This hypothesis has to be tested by searching for possible differences among volatile elements, such as CNO, S and Zn, and refractory ones, such as $`\alpha `$s Si, Mg, Ca, Ti and the iron-group elements. Smith et al. (Smi01 (2001)) found out that a small subset of stars with planets bore this accretion signature. They discovered that these stars exhibit a trend of increasing \[$`X`$/H\] with increasing $`T_\mathrm{C}`$ for a given element $`X`$. However, a following study by Takeda et al. (Tak01 (2001)) showed that all the elements, volatiles and refractories alike, of fourteen planet-harbouring stars behave quite similar, suggesting that the enhanced metallicity is not caused by the accretion of rocky material but is rather of primordial origin. Recently, Sadakane et al. (Sad02 (2002)) confirmed this result for nineteen elements, volatiles and refractories, in twelve planet host stars. Gonzalez et al. (Gonz01 (2001)) and Santos et al. (San00 (2000)) did not find significant differences in the \[C/Fe\] and \[O/Fe\] values among planet-harbouring and field stars. Some of these authors derived nitrogen abundances from atomic lines (e.g. Gonzalez & Laws Gonz00 (2000); Gonzalez et al. Gonz01 (2001); Takeda et al. Tak01 (2001); Sadakane et al. Sad02 (2002)). We note that most of them made no mention of their final results and conclusions concerning nitrogen.

Various studies have accumulated evidence that the production of nitrogen at low \[Fe/H\] proceeds principally as a primary rather than a secondary process (Pagel & Edmunds Pag81 (1981); Bessel & Norris Bes82 (1982); Carbon et al. Car87 (1987); Henry et al. Hen00 (2000)). The secondary process dominates at high metallicities. Two possible primary sources of nitrogen have been proposed. The first is intermediate mass (4–8 $`M_{\mathrm{}}`$) stars during their thermally pulsing asymptotic giant branch phases: if temperature conditions are suitable, primary nitrogen can be produced by CN processing in the convective envelope from freshly synthesized dredged up carbon (Marigo Mar01 (2001); van den Hoek & Groenewegen van97 (1997)). The second source is rotating massive stars (Maeder & Meynet Mae00 (2000)). The latter would imply uncoupled N and Fe abundances and a nitrogen overproduction relative to iron.

Recently, Liang et al. (Liang01 (2001)) and Pettini et al. (Pet02 (2002)) concluded that the contribution by massive stars to the nitrogen abundance is generally very low, and that the dominant N contributors are intermediate and low mass stars (ILMS). Although the nitrogen is chiefly primary at low \[Fe/H\], ILMS are also sources of secondary nitrogen. The relative weights of the secondary and primary components depend on the interplay between secondary enrichment caused by dredge-up episodes and the primary contribution given by the CNO cycle during the envelope burning. Shi et al. (Shi02 (2002)) obtained similar results by measuring N i lines in a large sample of disk stars with a metallicity range of $`1.0<`$\[Fe/H\]$`<+`$0.2.

However, the situation is rather confused. Other authors support the origin of primary nitrogen in massive stars (e.g. Izotov & Thuan Izot00 (2000)) because of the low scatter of \[N/O\] ratios in galaxies observed at different stages of their evolution, which would imply no time delay between the injection of nitrogen and oxygen. Meynet & Maeder (Mey02 (2002)) have discussed the effect of rotation in the production of primary and secondary nitrogen in intermediate and high mass stars, introducing axial rotation and rotationally induced mixing of chemical species in their stellar models. They concluded that ILMS are the main source of primary nitrogen at low metallicities.

The main problem that has traditionally made it difficult to measure nitrogen abundances is that there are very few atomic nitrogen lines in the red part of the spectrum. Efficient modern detectors now enable us to obtain high-resolution near-UV spectra, which allows us to carry out a precise and independent study of nitrogen abundance by using the NH band at 3360 Å.

The purpose of this work is to obtain a systematic, uniform and detailed study of nitrogen abundance in two samples of solar-type dwarfs, a large set of planet hosts and a comparison volume-limited sample of stars with no known planetary mass companions. We used two independent indicators in the abundance determination: spectral synthesis of the NH band in near-UV high-resolution spectra and the near-IR N i 7468Å line.

## 2 Observations

Spectra for most of the planet-harbouring stars have been collected and used to derive precise stellar parameters in a series of recent papers (e.g. Santos et al. San01 (2001), San03a (2003)a, San03b (2003)b). For measurement of the N i 7468Å line we used all the spectra of this set obtained with the FEROS spectrograph on the 2.2 m ESO/MPI telescope (La Silla, Chile), the SARG spectrograph on the 3.5 m TNG and the UES spectrograph on the 4.2 m WHT (both at the Roque de los Muchachos Observatory, La Palma), and other spectra from subsequent observational runs using the same instruments.

For NH band synthesis, we used the same spectra as those in Santos et al. (San02 (2002)) to derive beryllium abundances111Observing run 66.C-0116A. New spectra were obtained with UVES spectrograph, at the VLT/UT2 Kueyen telescope (Paranal Observatory, ESO, Chile)222Observing run 68.C-0058A. The observing log for the spectra from campaigns carried out after Santos et al. (San02 (2002)) and Bodaghee et al. (Bod03 (2003)) are listed in Tables 1 and 2. The obtained high S/N spectra have a minimum resolution $`R50000`$.

The data reduction for the SARG and UVES spectra was done using IRAF tools in the echelle package. Standard background correction, flatfield and extraction procedures were used. The wavelength calibration was performed using a ThAr lamp spectrum taken on the same night. The FEROS spectra were reduced using the FEROS pipeline software.

## 3 Analysis

Abundances for all elements were determined using a standard local thermodynamic equilibrium (LTE) analysis with the revised version of the spectral synthesis code MOOG (Sneden Sne73 (1973)) and a grid of Kurucz (Kur93 (1993)) ATLAS9 atmospheres. All the atmospheric parameters, $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$, $`\xi _t`$ and \[Fe/H\], and the corresponding uncertainties were taken from Santos et al. (San03b (2003)b). The adopted solar abundances for nitrogen and iron were $`\mathrm{log}ϵ(N)_{}`$ = 8.05 dex (Anders & Grevesse And89 (1989)) and $`\mathrm{log}ϵ(Fe)_{}`$ = 7.47 dex (Santos et al. San03b (2003)b).

### 3.1 Synthesis of NH band

The NH molecular feature is the strongest one in the spectral region $`\lambda \lambda `$3345–3375 Å. We determined nitrogen abundances by fitting synthetic spectra to data in this wavelength range. The NH spectra were calculated with the dissociation potential $`D_0(NH)=3.37\pm 0.06eV`$ recommended in Grevesse et al. (Gre90 (1990)).

A detailed line list from Yakovina & Pavlenko (Yak98 (1998)) was adopted. These authors obtained a good fit to the Kurucz Solar Atlas (Kurucz et al. Kur84 (1984)) by changing oscillator strength values of the strongest lines and continuum level in a list of atomic and molecular lines from Kurucz’s 1993 database. In our analysis, we slightly modified $`\mathrm{log}`$ gf values of this line list333The full line list is available in www.iac.es (changes smaller than 0.01 dex in most cases) in order to avoid modifying the continuum level when fitting the high-resolution Kurucz Solar Atlas (Kurucz et al. Kur84 (1984)) with a solar model having $`T_{\mathrm{eff}}=5777`$ K, $`\mathrm{log}g=4.44`$ dex and $`\xi _t=`$1.0 km s-1, assuming the solar N abundance of $`\mathrm{log}ϵ(N)_{}`$ = 8.05 dex.

All our targets are solar-type dwarfs with \[Fe/H\] values between $`0.6`$ and $`+0.4`$ dex and $`T_{\mathrm{eff}}`$ between 5000 and $`6300`$ K. Stars with lower $`T_{\mathrm{eff}}`$ were rejected since they would have introduced great uncertainties in the abundance results. The discarded targets were HD 27442, HD 192263, HD 74576 and HD 23484, with $`T_{\mathrm{eff}}`$ of 4825 K, 4947 K, 5000 K and 5176 K, respectively.

The main difficulties were the severe atomic line contamination of the NH features and the crowding of the lines in the UV region. The oscillator strengths of these lines are not known with high accuracy. The crowding of blended lines worsens with the increase of metallicity. This characteristic involves difficulties when fixing the continuum level since the whole spectral region is depressed by line absorption. Continuum was normalized with 5th order polynomia using the CONT task of IRAF. However we could make further improvements in the placement of continuum using the Kurucz Solar Flux Atlas (Kurucz et al. Kur84 (1984)) and DIPSO task of the STARLINK software: continuum level was determined by superimposing observed and synthetic solar spectra and selecting some points of reference. Figure 1 shows significant points in continuum determination, selected in the spectral regions $`\lambda \lambda `$3344–3350 Å , $`\lambda \lambda `$3350–3358 Å and $`\lambda \lambda `$3362–3370 Å for HD 17051. This method was valid for stars with $`T_{\mathrm{eff}}`$ and metallicities around solar values, while in cooler objects we were not able to assure that continuum level could be scaled to the Sun.

In order to find the best synthetic spectrum for each star, we used the program FITTING. This program allows us to create a grid of synthetic spectra for several free parameters, such as abundances, atmospheric parameters, metallicity, etc., and to compare it with the observed spectrum. To generate synthetic spectra, it calls the spectral synthesis code MOOG. The best fit is determined by applying a $`\chi ^2`$ minimization method to the spectral regions considered to be the most significant.

We used a Gaussian function with a FWHM of 0.07 Å for the instrumental broadening, and a rotational broadening function with $`vsini`$= 4 km s-1. No macroturbulence broadening was used. We assumed the same rotational velocity for all the targets since small differences would affect only $`\chi ^2`$ values, not final nitrogen abundances. Only for HD 120136 and HD 19994 we took $`vsini`$= 14.5 km s-1 and $`vsini`$= 8.2 km s-1 from the CORALIE database. We stress that planet searches based on the Doppler technique are biased in favour of old stars with low $`vsini`$ and therefore all our targets are slow rotators. In any case, precise values of $`vsini`$ are not required for the abundance analysis since it is well known that the rotational convolution is not affecting the EW of the line (Gray Gra92 (1992)). Our test with $`\chi ^2`$ analysis confirms this assumption for the stars considered in our study. In the estimate of damping constants, the Unsold aproximation was adopted for all atomic and molecular lines. We have verified that the damping effects do not affect our final abundances.

At first, the abundances of all relevant elements were fitted for a set of stars with different atmospheric parameters. Next, the same fits were performed by scaling the abundance of elements other than nitrogen to the \[Fe/H\] value and by changing only the nitrogen abundance. Same results were obtained in both cases, therefore all the following fittings were carried out by scaling elements other than nitrogen to iron.

The program FITTING introduced the following data in the code MOOG: the atmospheric model, the line list and the range of nitrogen abundance of the grid of synthetic spectra. Then $`\chi ^2`$ values were calculated by comparing each synthetic spectrum with the observed one in the following spectral regions: $`\lambda \lambda `$3344.0–3344.7 Å, $`\lambda \lambda `$3345.9–3346.9 Å, $`\lambda \lambda `$3347.2–3348.2 Å, $`\lambda \lambda `$3353.5–3354.5 Å, $`\lambda \lambda `$3357.0–3358.0 Å, $`\lambda \lambda `$3358.5–3359.0 Å, $`\lambda \lambda `$3361.8–3362.5 Å, $`\lambda \lambda `$3363.0–3363.75 Å, $`\lambda \lambda `$3363.8–3365.5 Å, $`\lambda \lambda `$3368.0–3368.8 Å, $`\lambda \lambda `$3369.1–3369.9 Å, $`\lambda \lambda `$3371.8–3372.8 Å and $`\lambda \lambda `$3374.8–3376.0 Å. $`\chi _{tot}^2`$ was the sum of $`\chi ^2`$ values444$`\chi ^2=_{i=1}^N(F_iS_i)^2`$, where $`F_i`$ and $`S_i`$ are the observed and synthetic spectra fluxes, respectively, at the wavelength i. of the mentioned spectral regions. The nitrogen abundance corresponding to minimum $`\chi _{tot}^2`$ value was extracted and listed in Tables 4 and 5. The fact that $`\chi _{tot}^2`$ is smaller for high $`T_{\mathrm{eff}}`$ can be understood as the consecuence of a suitable placement of continuum in stars more similar to the Sun.

Figure 2 shows fits for a short, significant spectral region. The observed solar spectrum from the Kurucz Solar Atlas (Kurucz et al. Kur84 (1984)) and synthetic spectrum are represented in the left panel. The case of the planet host star HD 17051 is represented in the right panel (the observed spectrum and three synthetic spectra for different values of nitrogen abundance).

Uncertainties in the atmospheric parameters are of the order of 50 K in $`T_{\mathrm{eff}}`$, 0.12 dex in $`\mathrm{log}g`$, 0.08 km s-1 in the microturbulence and 0.05 dex in the metallicity (see Santos et al. San03b (2003)b). In order to know how variations in the atmospheric parameters affect NH abundances, for each parameter a set of three stars having very different parameter values was selected. For each set of stars we then tested \[N/H\] sensitivity to changes in the parameter ($`\pm 100`$ K in the case of $`T_{\mathrm{eff}}`$, $`\pm 0.2`$ dex in $`\mathrm{log}g`$ and \[Fe/H\], and 0.3 km s-1 in $`\xi _t`$). The results are shown in Table 3. No microturbolence effects were included, since an increase of 0.3 km s-1 produced an average decrease of 0.003 dex in nitrogen abundances, which is negligible with regard to the effects of other parameters. Sensitivity of NH abundance was applied in the propagation of each parameter error on abundances. Also, a continuum uncertainty effect of the order of 0.1 dex was considered. All effects were added quadratically to obtain final the uncertainties for the nitrogen abundances.

### 3.2 The N i absorption at 7468.27 Å

Nitrogen abundance analysis based on the feature at 7468.27 Å was carried out diferentially with respect to the Sun. Wavelength and excitation energies of lower level were taken from Moore et al. (Moo66 (1966)) and a solar gf value ($`\mathrm{log}gf=0.15`$) was computed using equivalent width measured in the Kurucz Solar Atlas (Kurucz et al. Kur84 (1984)) (4.2 mÅ) and a solar model with $`T_{\mathrm{eff}}=5777`$ K, $`\mathrm{log}g`$ = 4.44 dex and $`\xi _t=`$1.0 km s-1, assuming the solar N abundance of $`\mathrm{log}ϵ(N)_{}`$ = 8.05 dex. Equivalent widths in the observed stars were determined by Gaussian fitting using the SPLOT task of IRAF, and abundances were computed using the ABFIND driver of MOOG.

The high resolution spectra we used provided us with very reliable measurements in most cases. However, a strong fringing effect around 7468 Å prevented us from measuring high precision EW values in some cases; therefore, the corresponding EW uncertainties are more significant.

The sensitivity of the N i line to variations in atmospheric parameters was estimated in the same way as in the NH case (see Subsection 3.1). The results are shown in Table 3. In order to take into account inaccuracies caused by the continuum determination, equivalent widths for the highest and the lowest continuum levels were measured, and the corresponding abundance errors were added in quadrature to the abundance uncertainties derived from inaccuracies in the atmospheric parameters.

Dependences on $`T_{\mathrm{eff}}`$ and on $`\mathrm{log}g`$ of \[N/Fe\] performed by both methods are presented in Figure 3. We note that no characteristic trends appear in either case. This means that our results are almost free from systematic errors.

## 4 Nitrogen abundances from near-UV spectra

In the literature there are already a few studies that use NH band at 3360 Å to determine nitrogen abundances by spectral synthesis, but only for metal-poor dwarfs (Bessel & Norris Bes82 (1982); Tomkin & Lambert Tom84 (1984); Laird Laird85 (1985); Carbon et al. Car87 (1987)).

We analysed near-UV high-resolution spectra of 28 planet host stars and 25 comparison sample stars with the goal of checking for possible differences in nitrogen abundances between both samples. All the results for planet host and comparison sample stars are presented in Tables 4 and 5, respectively.

Figure 4 shows the \[N/H\] and \[N/Fe\] abundance ratios as functions of \[Fe/H\] for these two samples. These plots indicate that the abundance trend for stars with planets is almost indistinguishable from that of the comparison stars. Because of the iron enhancement that planet-harbouring stars present, their abundance distributions are the high \[Fe/H\] extensions of the curves traced by the comparison sample stars. No peculiar behaviour of the nitrogen abundance seems to be associated with the presence of planets.

On the whole, the \[N/Fe\] vs. \[Fe/H\] trend appears to be almost flat, since the slope value is $`0.10\pm 0.05`$. Nevertheless, a small trend to increase with \[Fe/H\] cannot be excluded. The \[N/H\] vs. \[Fe/H\] plot shows that the behaviour of the N abundance is not significantly different from that of Fe; both elements keep pace with each other. Moreover, no significant nitrogen overproduction relative to iron is observed.

The \[N/H\] distributions for both samples are represented in Fig. 5 (upper panel). An interesting feature of the histogram is that the planet-harbouring stars set does not present a symmetrical distribution, as the comparison sample does. The distribution increases for increasing \[N/H\] values until reaching the highest value, not centred on the \[N/H\] range. However, this peculiarity could be caused by the lack of objects with metallicities higher than $`0.5`$ dex. Santos et al. (San01 (2001), San03a (2003)a, San03b (2003)b) found a similar shape in iron distribution.

## 5 Nitrogen abundances from the N i 7468 Å line

Several studies of nitrogen abundances in metal-rich dwarfs with planets using near-IR N i lines have been published, but data for only a small number of planet host stars have been employed (fourteen planet-harbouring stars in the most extensive case, in Takeda et al. Tak01 (2001)) and most of them do not even mention their nitrogen abundance results (Gonzalez & Laws Gonz00 (2000); Gonzalez et al. Gonz01 (2001); Sadakane et al. Sad02 (2002)). Only Takeda et al. (Tak01 (2001)) have presented their nitrogen abundance results from the measurement of the N i line at 8683 Å, as well as abundance results for another four volatile elements (C, O, S and Zn). They did not observe the tendency predicted by the “self-enrichment” hypothesis (\[$`X`$/H\] $`<`$ \[Fe/H\] for metal-enriched planet host stars), a weak inverse trend (\[$`X`$/H\] $`>`$ \[Fe/H\]) even seemed to appear for N and S.

We obtained nitrogen abundances for 31 planet host stars by equivalent width measurements of the N i line at 7468 Å. All atmospheric parameters, EW values with uncertainties and abundance results are listed in Table 6. Figure 6 shows the \[N/H\] and \[N/Fe\] abundance ratios as functions of \[Fe/H\] for planet host stars, from NH band and N i line sets, and for comparison sample stars. In Figure 5 (top lower panel), the \[N/H\] distributions are shown for both the planet host (with values from both indicators) and comparison star samples. In all the plots the added values follow the same behaviour as those resulting from NH band analysis. The previously commented trends (see Section 4) are confirmed.

## 6 Comparison of the near-UV and 7468 Å line results

Two independent analyses were carried out, using different indicators of nitrogen abundance. We managed thereby to extend our data sample, as well as to check the agreement between both methods and to test the strength of our results.

We had available to us two different sets of data. The first consisted of near-UV spectra of planet host and comparison sample stars. The other was composed of optical and near-IR spectra of stars with planets. The comparison between the results of both indicators was performed using the eight targets that both sets had in common.

Abundances resulting from both indicators are presented in Table 7. We stress the remarkable agreement in most cases. In the case of HD 19994 and HD 217107, both values fall in the uncertainty range resulting from the abundance analysis. Only for HD 46375 do the different indicators give different results.

## 7 Reanalysis of published EW data

Several recent studies have presented measurements of N i lines in targets with known planets (Gonzalez & Laws Gonz00 (2000); Gonzalez et al. Gonz01 (2001); Takeda et al. Tak01 (2001); Sadakane et al. Sad02 (2002)). We have collected all the EW values determined in these papers and used them to calculate nitrogen abundances with atmospheric parameters from Santos et al. (San03b (2003)b).

Gonzalez & Laws (Gonz00 (2000)) and Gonzalez et al. (Gonz01 (2001)) measured N i line at 7468 Å, as we did; therefore, we could compare EW values in some of the targets that both analyses had in common. Takeda et al. (Tak01 (2001)) and Sadakane et al. (Sad02 (2002)) used the N i line at 8683 Å, so we were not able to compare our EW values with theirs, but only \[N/H\]. Besides allowing us to carry out a comparison, the literature values added new abundance measures to our data sample.

This reanalysis was carried out in the same way described in Section 3.2. In the case of the N i line at 8683 Å, the solar gf value computed using the equivalent width (7.0 mÅ) measured in the Solar Atlas (Kurucz et al. Kur84 (1984)) was $`\mathrm{log}gf=0.068`$. All the EW values from the literature, the atmospheric parameters from Santos et al. (San03b (2003)b) and the resulting nitrogen abundances are presented in Table 8. All the results for the targets in common with our data sample are listed in Table 9. We note that the results from different sources are in good agreement. Significant differences exist only for HD 186427, HD 187123 and HD 217107. The \[N/H\] and \[N/Fe\] vs. \[Fe/H\] trends from the literature EW values and from our data are represented in Figure 7. Figure 5 (bottom panel) shows the \[N/H\] distribution for comparison sample and planet host stars from all the sources. We note that the new values confirm the trends discussed in the previous sections (see Section 4).

## 8 Discussion and Conclusions

In the present study we have determined the nitrogen abundances in 51 planet host stars and in 25 stars from a comparison sample with no known planetary mass companion. Among the set of planet hosts, 28 targets have been analysed by spectral synthesis of the NH band at 3360 Å, while equivalent width measurements of the N i 7468.27 Å were carried out for 31 stars. Another 15 planet hosts were subsequently added to our study, with EW values taken from other authors’ works. Moreover, we were able to check that the different indicators are in remarkable agreement for several common targets. This result is an independent and homogeneous analysis of nitrogen abundances in 91 solar-type dwarfs, 66 with planets and 25 from a volume-limited comparison sample.

The behaviour of volatiles in planet-harbouring stars with regard to a comparison sample can be very informative for checking the “self-enrichment” hypothesis. If the accretion of metal-enriched planetary material were a key parameter for the observed enhancement of iron abundances, then volatile abundances in planet host stars would not show as much overabundance as refractories do since elements with a lower condensation temperature are expected to be deficient in accreted materials. In that case, the following trend would be observed in planet host stars: \[N/H\] $`<`$ \[Fe/H\].

Our results suggest that planet hosts stars do not behave in such a way. The \[N/H\] vs. \[Fe/H\] plots show that nitrogen abundance scales with that of the iron. In the \[N/Fe\] vs. \[Fe/H\] plots, no significant trace of \[N/Fe\] $`<`$ 0 appears in stars with planets. Moreover, planet host and comparison sample stars have the same behaviour. Trends from either sample are quite indistinguishable. This emphasizes that the trend predicted from the “self-enrichment” scenario is not observed at all.

Previous studies have obtained similar results for other volatiles in some planet-harbouring stars. Gonzalez et al. (Gonz01 (2001)) corrected the low \[C/Fe\] values found in a previous study (Gonzalez & Laws Gonz00 (2000)) and concluded that \[C/Fe\] and \[O/Fe\] in planet hosts do not display significant differences from those in field dwarfs of the same \[Fe/H\]. Takeda et al. (Tak01 (2001)) found the same trend as we did for N, C, S, O and Zn in 14 planet hosts, and in Sadakane et al. (Sad02 (2002)) this result was confirmed for C and O in another 12 planet-harbouring stars. Our study confirms these results for a significantly larger set of stars with planets, as well as in a comparison sample of stars with no discovered planetary-mass companion. This seems to support a scenario in which the formation of planets is particularly dependent on the high metallicity of the primordial cloud.

In the future, it will be important to pursue further uniform studies of volatile element (C, O, S and Zn) abundances in planet hosts and stars with no known planets. Well-defined \[$`X`$/Fe\] vs. \[Fe/H\] distributions can provide us with precise information in order to confirm or discard hypotheses about planetary formation.

Another important piece of information that emerges from nitrogen abundance trends concerns nitrogen sources. The latest studies propose two sources of primary nitrogen: intermediate mass (4–8 $`M_{\mathrm{}}`$) and massive stars ($`M>8M_{\mathrm{}}`$). If the latter source were the main contributor to nitrogen at high \[Fe/H\], N and Fe abundances would appear to be uncoupled and nitrogen would be overabundant relative to iron. Our study of nitrogen abundances in 91 solar-type dwarfs suggests that nitrogen keeps pace with iron. The \[N/H\] vs. \[Fe/H\] plots show that both elements behave quite similarly, and the \[N/Fe\] vs. \[Fe/H\] plots indicate that no significant nitrogen overabundance exists relative to iron in the metallicity range $`0.6<`$ \[Fe/H\] $`<+0.4`$. We also note that Shi et al. (Shi02 (2002)) recently obtained that \[N/Fe\] is about solar in a sample of 90 disk stars; therefore, they concluded that nitrogen is produced mostly by intermediate mass stars.

Nevertheless, secondary production due to ILMS is the predominant effect on the \[N/Fe\] curve around solar metallicities. Only the \[N/Fe\] trend at \[Fe/H\] $`<1.0`$ can provide conclusive arguments about primary sources. At this time, we are extending our study of \[N/H\] ratios to lower metallicities, by applying the same analysis of the NH band at 3360 Å to a sample of metal-poor dwarfs. The results will be very informative in the framework of investigating the sources of primary nitrogen in the early Galaxy.

###### Acknowledgements.

We would like to thank J. I. González Hernández for kindly providing us with his program FITTING. We wish to thank Dr. G. Meynet (Observatoire de Genève, Switzerland) for many fruitful discussions and his comments on the text. The anonymous referee is thanked for many useful suggestions and comments. IRAF is distributed by the National Optical Astronomy Observatories, operated by the Association of Universities for Research in Astronomy, Inc., under contract with the National Science Foundation, USA.


text:
# Large 𝑁𝑐 Behavior of Light Scalar Meson Nonet Revisited

## 1 Introduction

It has been shown by Peláez  that the complex poles corresponding to the $`\rho `$ and $`K^{}`$ mesons move toward the real axis on the second Riemann sheet as a number of colors, $`N_C`$, becomes large, and that in contrast to the vector mesons those corresponding to the $`\sigma `$ and $`\kappa `$ states move away from the real axis, within the inverse amplitude method (IAM) using full $`O(p^4)`$ amplitudes of the chiral perturbation theory (ChPT).

Stimulated by this work, we have calculated how physical quantities such as phase shifts and cross sections behave on the real axis of the physical sheet as $`N_c`$ increases from 3 to some finite values. In this calculation we adopt an approximate version of the two-channel IAM developed by Oller-Oset-Peláez (OOP), which we call the OOP version. We have obtained the result that while the vector mesons become narrower resonances, the light scalar meson nonet including the $`f_0(980)`$ and $`a_0(980)`$ states fade out as $`N_c`$ exceeds about 6. This result has been confirmed by Peláez in his new paper, excluding an exceptional case of the $`a_0(980)`$ state.

The same issue have been discussed by Oller and Oset using a different model, in which chiral $`O(p^2)`$ amplitudes and possible preexisting tree resonance poles are introduced as ingredients of the model. Their criterion whether a meson is dynamical or not is that the partial wave amplitude of the model can reproduce the experimental data without a nearby preexisting pole. The preexisting poles are assumed to survive in the large $`N_c`$ limit. They have concluded that while the $`\rho `$ and $`K^{}`$ mesons need each preexisting pole, light scalar mesons, possibly except for the $`f_0(980)`$ state, do not necessarily need such poles. We note, however, that if energies to be fitted by the model are restricted to 1.2 GeV, the preexisting pole is not needed for the $`f_0`$ state.

Thus, the above observations are consistent with the common understanding that the members of the vector meson nonet including $`\rho `$ and $`K^{}`$ are typical of $`q\overline{q}`$ mesons in large $`N_c`$ QCD. On the other hand the behavior of scalar mesons is at variance with the nature of the $`q\overline{q}`$ mesons.

In this paper we study again the behavior of two-meson scattering amplitudes when $`N_c`$ increases from 3 under some different conditions from previous work, within the two-channel OOP version with $`O(p^4)`$ amplitudes given in Ref. . The $`O(p^4)`$ amplitudes depend on values of the low energy constants (LEC) of ChPT, denoted as $`L_n`$, which are to be determined phenomenologically so as to reproduce experimental data. It would be difficult, therefore, to discriminate the nature of resonances by studying the $`N_c=3`$ world alone. In order to know the nature of the low mass mesons, it will be very useful to study how light vector and scalar meson states behave and how complex poles of $`f_0(980)`$ and $`a_0(980)`$ move when we increase $`N_c`$ from 3 to some values. Although we do not go far away from the real $`N_c=3`$ world with $`SU(3)\times SU(3)`$ chiral symmetry, we observe that the $`\rho `$ and $`K^{}`$ meson survive as narrow width resonances, but the light scalar nonet fade out as $`N_c`$ exceeds a rather small number about 6. This is the same result as in the previous work.

In the next section the explicit $`N_c`$ dependence of the OOP amplitudes is given, the vector channel and scalar channel are discussed in section 2 and 3, and the conclusions and discussion is given in the last section.

## 2 $`N_c`$ dependence of the OOP amplitudes

In order to carry out the study we have to find the explicit $`N_c`$ dependence of the scattering amplitudes. The amplitudes in ChPT have an explicit $`N_c`$ dependence through the pion decay constant and the LECs. Since the pion decay constant $`f_\pi `$ is of $`O(N_c^{1/2})`$ and the LECs, $`L_1`$, $`L_2`$, $`L_3`$, $`L_5`$ and $`L_8`$ are to be of $`O(N_c)`$, but $`2L_1L_2`$, $`L_4`$, $`L_6`$ and $`L_7`$ are of $`O(1)`$, we put

$`L_n(N_c)`$ $`=`$ $`\widehat{L_n}{\displaystyle \frac{N_c}{3}}+\mathrm{\Delta }L_n,`$ (1)
$`f_\pi ^2(N_c)`$ $`=`$ $`\widehat{f_\pi ^2}{\displaystyle \frac{N_c}{3}},`$ (2)
where $`\widehat{L_n}`$ satisfy the relations, $`2\widehat{L_1}\widehat{L_2}=\widehat{L_4}=\widehat{L_6}=\widehat{L_7}=0`$, and $`\mathrm{\Delta }L_n`$ are of $`O(1)`$. Thus, we have

$$\frac{L_n}{f_\pi ^2}=\frac{\widehat{L_n}}{\widehat{f_\pi ^2}}+\frac{\mathrm{\Delta }L_n}{\widehat{f_\pi ^2}}\frac{3}{N_c}$$
(3)
with $`\widehat{f_\pi }=93`$ MeV. We also assume that the meson decay constants are the same and equal to the pion decay constant $`\widehat{f_\pi }`$ as in Ref. .

The ingredients of the IAM consist of amplitudes of chiral order $`O(p^2)`$ and $`O(p^4)`$ of ChPT. An $`O(p^2)`$ amplitude, denoted by $`T^{(2)}(s,t,u)`$, has a form of a linear function of $`s,t,u`$ divided by $`f_\pi ^2`$, and it is of $`O(N_c^1)`$. A polynomial term of the latter amplitudes, denoted by $`T_{\mathrm{poly}}^{(4)}(s,t,u)`$, is written as a sum of polynomial functions with the LECs as follows:

$$T_{\mathrm{poly}}^{(4)}(s,t,u)=\underset{n=1,8}{}\frac{1}{f_\pi ^2}\left(\frac{L_n}{f_\pi ^2}\right)P_n(s,t,u),$$
(4)
where $`P_n`$ are quadratic functions of $`s,t,u`$ and meson mass squared. The polynomial term $`T_{\mathrm{poly}}^{(4)}`$ is of $`O(N_c^1)`$, because $`L_n/f_\pi ^2`$ scales as $`O(N_c^0)`$ as seen in Eq.(3). An $`s`$-channel loop term given by $`t^{(2)}(s)J(s)t^{(2)}(s)`$ is of $`O(N_c^2)`$, where $`J(s)`$ is the one-loop function regularized as the $`\overline{MS}1`$ scheme at the renormalization scale $`\mu `$, where $`t^{(2)}`$ is a partial wave amplitude derived from $`T^{(2)}(s,t,u)`$. Similarly $`t`$\- and $`u`$-channel loop terms and tadpole terms are of $`O(N_c^2)`$, which are ignored in the OOP version. Thus, the OOP version is expected to be more valid as $`N_c`$ becomes larger. The $`s`$-channel loop terms are indispensable to realize unitarity, although they are of $`O(N_c^2)`$. This difference of the $`N_c`$ dependence produces the different behavior of the amplitudes when $`N_c`$ becomes large.

Our set of the LECs used in this work are determined at the renormalization scale $`\mu =900`$ MeV so as to reproduce experimental phase shifts qualitatively up to about 1.2 GeV. We note that the IAM and OOP amplitudes contain the LECs non-linearly and the fitting region is extended to higher energies. Our sets of $`\widehat{L_n}`$ and $`\mathrm{\Delta }L_n`$ are tabulated in Table I with the set of the large $`N_c`$ model, which correspond to $`\widehat{L_n}`$. The change of the renormalization scale affects values of $`\mathrm{\Delta }L_n`$, but we do not consider the scale change explicitly because the $`\mathrm{\Delta }L_n`$ terms fade out as $`N_c`$ increases.

We emphasize that since the Large $`N_c`$ set can reproduce the low energy scattering behavior rather well even at $`N_c=3`$, except for unessential points, the results are almost the same if we take Large $`N_c`$ set instead of Our $`\widehat{L_n}`$ and put $`\mathrm{\Delta }L_n`$ to the difference between our $`L_n`$ and the Large $`N_c`$ set.

## 3 Vector channels

At first, we discuss the behavior of vector mesons in the single channel calculation. The mass of a vector resonance is controlled by the combination of LEC, $`2L_1L_2+L_3`$,  which is present in the term $`\mathrm{Re}[t^{(2)}t_{\mathrm{poly}}^{(4)}]`$ of $`O(N_c^1)`$, and the loop contribution to the real part is of $`O(N_c^2)`$, so that the mass stays at an almost constant value. This is the IAM expression substituting for the preexisting pole. On the other hand the imaginary part of the loop term contributes to the width. The octet component of the isoscalar vector meson

also has a constant mass and a residue of $`O(N_c^1)`$ below the $`K\overline{K}`$ threshold. Thus, the masses of the vector mesons are of $`O(1)`$, while the widths decrease as $`O(N_c^1)`$.

If we extend the calculation to the multi-channel IAM, the masses are almost unchanged, because the same combination of the LECs dominantly determine the masses. The $`N_c`$ dependence of the phase shifts and that of the cross sections of $`\rho `$ and $`K^{}`$ are shown in Figs. 1 and 2. The same behavior is observed in the case of the Large $`N_c`$ set, though the obtained masses are too small because of too large value of $`|L_3|`$.222The constraint $`2L_1=L_2`$ induces a large cancellation in $`det[t^{(2)}t^{(4)}\mathrm{loop}]`$ even at $`N_c=3`$, so that the resultant amplitude gets a large uncertainty. But the narrowing width with increasing values of $`N_c`$ remains valid.

Thus, we can conclude that the vector mesons described by the OOP version have the nature consistent with the $`q\overline{q}`$ mesons.

## 4 Scalar channels

### 4.1 $`(I,J)=(0,0)`$

This channel contains the controversial $`\sigma (600)`$ and the $`f_0(980)`$ states. Using the two-channel IAM consisting of the $`\pi \pi `$ and $`K\overline{K}`$ channels, we can reproduce experimental data fairly well below 1.2 GeV. The $`N_c`$ dependence of the phase shift and the cross section are shown in Fig. 3, where $`N_c`$ increases from 3 to 15. In contrast to the vector channel we observe that the phase shift becomes flat and the cross section fades out as $`N_c`$ becomes large; the sharp rise of the phase shift near the $`K\overline{K}`$ threshold and the large bump of the cross section near

500 MeV seen at $`N_c=3`$ to 5 disappears at $`N_c=6`$, and then the phase shift and the cross section become almost flat and fade out. Similar drastic change in the $`N_c`$ dependence has also been observed in Ref. , though it is in different context.

The $`f_0(980)`$ pole exists at $`(97522i)`$ MeV at $`N_c=3`$. Where does the pole go as $`N_c`$ increases ? We approximately calculate the pole position by expanding the amplitude in powers of $`k_2`$ up to the first order, where $`k_2`$ is the momentum of the $`K\overline{K}`$ channel. We observe that the pole moves into the upper half plane of the IV sheet from the lower half plane of the II sheet, winding around the branch point at $`K\overline{K}`$ threshold, and goes away from the real axis as shown in the left side of Fig. 4. Pole positions at larger $`N_c`$ cannot be reliable owing to the rough approximation, but this behavior would remain intact.

We briefly comment on effects of adding the $`\eta \eta `$ channel to the OOP amplitude. If we include the $`\eta \eta `$ channel, we find that both real and imaginary part of $`det[t^{(2)}t^{(4)}\mathrm{loopterms}]`$ develop zeros at almost the same point near 770 MeV even at $`N_c=3`$ for a wide range of the LEC sets. Such unreasonable behavior is also seen in the isospinor channel with the $`\pi K`$ and $`\eta K`$ channels as will be noted. These zeros do not violate unitarity, but give very unreasonable behavior to the amplitude, that is too narrow resonant behavior with almost zero width. It should be noted, however, that the behavior of the amplitude, excluding a narrow strip including the zeros, is almost the same with the behavior of the two-channel model except for a shallow dip of the inelasticity at the $`\eta \eta `$ threshold. The unpleasant behavior is not reported in the calculations with full $`T^{(4)}`$. If we eliminate the unreasonable behavior, the fading-out tendency as increasing $`N_c`$ remain valid, though the fading-out occurs at larger $`N_c`$ owing to the additional $`\eta `$ loop contributions.

### 4.2 $`(I,J)=(1,0)`$

This channel contains the $`a_0(980)`$ state and appears as a cusp-like sharp peak as seen in Fig. 5.

The rising phase shift after the cusp bends down and becomes to a flat curve, and the cross section having a sharp peak fades out as $`N_c`$ increases. The pole appears at $`(109117i)`$ MeV on the II sheet at $`N_c=3`$, and it moves from the II sheet to the IV sheet, and leaves rapidly the real axis as shown in the left side of Fig. 4 as $`N_c`$ increases. We note that the real parts of the poles of the $`f_0(980)`$ and $`a_0(980)`$ states are not necessarily degenerate with each other, though both of the peaks of the mass distribution at $`N_c=3`$ appear near the $`K\overline{K}`$ threshold owing to the cusp behavior of the $`a_0(980)`$ state.

### 4.3 $`(I,J)=(1/2,0)`$

The fading-out behavior of the phase shift and cross section with increasing $`N_c`$ is the same as that in the channels discussed above. As stated before there appears an artifact zero near 750 MeV originated from the $`\eta K\eta K`$ component in the OOP version used in this work. So if we eliminate the zero by an interpolation method, we observe that the results by the two channel model are almost the same as those by the single channel calculation by virtue of the weak coupling between the $`\pi K`$ and $`\eta K`$ channels. The calculations making use of the full $`T^{(4)}`$ do not give such an unwanted zero .

## 5 Concluding remarks

We have calculated the $`N_c`$ dependence of the vector and scalar channels stating from $`N_c=3`$ to finite values, 30 for the vector channel and 12 or 15 for the scalar channel within the approximate IAM under the $`N_c`$ dependence of $`L_n/f_\pi ^2`$ given by Eq.(3) and $`f_\pi (N_c)=\sqrt{N_c/3}\times f_\pi (3)`$. And we have observed that the vector mesons survive as sharper resonances at almost the same position, the resonant structures of the scalar channels fade out at rather low values of $`N_c`$ near 5 or 6. By extending the observation we are led to conclude that the vector meson nonet has the nature consistent with the $`q\overline{q}`$ mesons in large $`N_c`$ QCD, but the light scalar meson nonet cannot survive in large $`N_c`$ and then cannot have the nature of the $`q\overline{q}`$ mesons. This conclusion is the same as obtained by Peláez, excluding an exceptional case of the $`a_0(980)`$ state, and it is also consistent with the results obtained in Ref. . Our conclusion supports the arguments that the light scalar mesons are of the $`K\overline{K}`$ molecule, and of $`q^2\overline{q}^2`$ states.

If the $`f_0`$ and $`a_0`$ states are composed of $`(qs)(\overline{q}\overline{s})`$ state, where $`q`$ denotes $`u`$ and/or $`d`$ quark, the similarity between the $`f_0`$ and $`a_0`$ both in mass and generating mechanism is expected. However, the pole positions of the both states can be different from each other by about 100 MeV or more, and the generating mechanism of the $`a_0(980)`$ state would be different from that of $`f_0`$ state: The $`a_0`$ state is generated by the strong channel coupling between $`\pi \eta `$ and $`K\overline{K}`$ channels, but not as a bound state resonance like as the $`f_0`$ state, as seen in the exchange dynamics and in the chiral loop dynamics. There is also the argument that the $`a_0(980)`$ and $`f_0(980)`$ are not elementary particles within the hadronic dynamics. They have observed that the field renormalization constants $`Z`$ of both states are close to 0 using the propagators of existing models, and concluded that a simple $`q\overline{q}`$ or four quark assignment for the $`a_0`$ should be considered with caution and it is certainly questionable for the $`f_0`$. If the light scalar nonet are not $`q\overline{q}`$ mesons, we cannot include them into mass spectra in the low meson dominance hypothesis, because the low mesons are supposed to participate in narrow resonance towers in large $`N_c`$ limit of QCD.

Our conclusion strongly indicates that all of the light scalar nonet are dynamical effects originating from unitarity, chiral symmetry and strong channel couplings. If the mesons in the scalar nonet are dominantly composed of hadronic or four quark component, but include $`|q\overline{q}>_P`$ with a small fraction as in Ref. , we could find out the small $`|q\overline{q}>_P`$ component by increasing $`N_c`$ in theoretical models, because the large hadronic or four quark component fades out and the $`q\overline{q}`$ component remains. At least, our calculation within the two-channel OOP approximation does not indicates that such an intriguing change will occur in larger $`N_c`$ region.


text:
# Front Propagation in Reaction-Superdiffusion Dynamics - Taming L vy Flights with Fluctuations

## Abstract

We investigate front propagation in a reacting particle system in which particles perform scale-free random walks known as L vy flights. The system is described by a fractional generalization of a reaction-diffusion equation. We focus on the effects of fluctuations caused by a finite number of particles. We show that, inspite of superdiffusive particle dispersion and contrary to mean field theoretical predictions, wave fronts propagate at constant velocities, even for very large particle numbers. We show that the asymptotic velocity scales with the particle number and obtain the scaling exponent.

One of the fundamental processes involved in non-equilibrium pattern formation is the spatial propagation of interfaces or fronts. Front propagation usually emerges when a local reaction dynamics interplays with diffusion in space of the reacting agents and has been observed in a wide range of physical, chemical and biological systems Cross and Hohenberg (1993). Prominent examples are patterns of bacterial colonies Golding et al. (1998), spreading phenomena in population genetics Fisher (1937); Aronson and Weinberger (1978) and interface dynamics in superconducters Bartolo and Dorsey (1996). One of the most prominent models which displays propagating fronts is the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation for the spatial concentration $`u(x,t)`$ of a reacting agent,

$$_tu=\gamma u(1u)+\mathrm{\Delta }u,$$
(1)
where diffusive motion of the reacting agents is assumed. However, this assumption cannot be justified for a number of systems. In fact, superdiffusive dispersion in space has been observed in a wide range of physical and biological systems, e.g. intermittent chaotic systems Geisel et al. (1985), bacterial motion Levandowsky et al. (1997), and foraging patterns of albatrosses Viswanathan et al. (1996).

Superdiffusive stochastic motion is usually characterized by a lack of scale in the microscopic single step distribution. One of the most successful theoretical concepts devised for the understanding of superdiffusion is a class of random walks known as L vy flights Shlesinger et al. (1995). A L vy flight consists of random single steps $`\mathrm{\Delta }x`$ which are drawn from an inverse power-law pdf (probability density function) $`p(\mathrm{\Delta }x)|\mathrm{\Delta }x|^{(1+\mu )}`$ characterized by a L vy exponent $`0<\mu <2`$. Due to the heavy tail, the variance in step size is divergent, the process lacks a spatial scale and the position $`X(t)`$ of a L vy flight scales heuristically with time $`t`$ as $`X(t)t^{1/\mu }`$. The associated diffusion equations contain fractional generalizations of ordinary derivatives Metzler and Klafter (2000); Sokolov et al. (2002). These fractional Fokker-Planck equations can exhibit behaviors strikingly different from ordinary ones Brockmann and Geisel (2003a) and have found wide application in physics, e.g. protein motion on folded hetero-polymers Brockmann and Geisel (2003b) and the dynamics of modern epidemics Hufnagel et al. (2004); Brockmann and Hufnagel .

In two recent studies wave front dynamics was shown to be drastically different from ordinary reaction-diffusion dynamics when the reacting agents move superdiffusively Mancinelli et al. (2002); del Castillo-Negrete et al. (2003). The authors considered fractional generalizations of Eq. (1) and showed that the spatio-temporal shape $`u(x,t)`$ of the leading egde of a propagating front has a power-law tail along the spatial coordinate and accelerates exponentially in time as opposed to the constant velocity and exponential decay in space exhibited by ordinarily diffusive system. The predictions made by mean field theory thus indicate that scale-free, superdiffusive dispersion of the reacting agents excludes constant velocity wave fronts and induces an entirely different spatio-temporal behavior.

However, as has been shown in a number of recent studies, the effect of fluctuations can be rather profound in these systems Pechenik and Levine (1999); Brunet and Derrida (2001); Doering et al. (2003). For instance, fluctuations can destabilize homogeneous wave fronts in two-dimensional systems Kessler and Levine (1998). A finite albeit large number $`N`$ of particles or reacting agents leads to a multiplicative noise term in the reaction-diffusion equations and although the variance of the noise is of order $`1/N`$, corrections to macroscopic quantities such as the front speed $`v`$ scale as $`(\mathrm{log}N)^2`$, a significant correction to the mean field approximation even for very large values of $`N`$. Brunet and Derrida Brunet and Derrida (1997); Kessler et al. (1998) extended mean field dynamics by an effective cutoff parameter $`\epsilon `$ for the concentration of particles below which no reaction and hence no exponential growth of the leading edge of a front is possible. Despite the fact that a rigorous equivalence with multiplicative noise is still lacking, the effective cutoff approach is very intuitive and in remarkable agreement with simulations of the full probabilistic dynamics.

Here, we focus on the effect of fluctuations on reaction-superdiffusion kinetics. We show that for arbitrarily small fluctuations (i.e. arbitrarily large particle numbers $`N`$), wave fronts propagate asymptotically at constant velocities. Furthermore, we show that as soon as fluctuations enter the description the algebraic tail along the spatial coordinate of the leading edge disappears and is replaced by an exponential decay. Thus, despite the fact that reacting agents move superdiffusively in space, the wave front patterns are qualitatively the same as in the ordinary diffusion case. We show that a front speed $`v`$ is selected after a transient time and that $`v`$ scales with particle number as $`vN^{1/\mu }`$ for L vy exponents $`\mu <2`$. The results reported here are rather counterintuitive, deviate strongly from the predictions cast by mean field theory, and indicate that fluctuations affect reactions-superdiffusion systems severely and may by no means be neglected.

We begin with a simple two particle type ($`A,B`$) reaction scheme,

$$A_x+B_x\stackrel{k_1}{}2A_x,A_x+B_x\stackrel{k_2}{}2B_x$$
(2)
$$A_x,B_x\stackrel{f(|xy|)}{}A_y,B_y.$$
(3)
Particles of type $`A`$ and $`B`$ either react to produce two particles of type $`A`$ or two particles of type $`B`$ at rate $`k_1`$ and $`k_2`$, respectively. Furthermore, particles of both types may jump from position $`x`$ to position $`y`$ with a probability density rate $`f(|xy|)`$ which we assume to be a decreasing function of distance $`|xy|`$. The dynamic stochastic variables are the particle numbers $`n_A(x,t)`$ and $`n_B(x,t)`$ of particles in an volume of size $`\mathrm{\Omega }`$ around $`x`$ of type $`A`$ and $`B`$, respectively. The volume $`\mathrm{\Omega }`$ is assumed to be large enough to contain a large number of particles but small compared to the spatial extend of the system. The total number of particles in $`\mathrm{\Omega }`$ around $`x`$ is given by $`N(x,t)=n_A(x,t)+n_B(x,t)`$. Without spatial dispersion of particles the local dynamics is governed by the master equation

$$_tp(n,t)=\underset{m}{}w(n|m)p(m,t)w(m|n)p(n,t)$$
(4)
for the probability $`p(n,t)`$ of finding a number $`n=n_A`$ of particles of type $`A`$ at a location $`x`$ with initial condition $`p(n,0)=\delta (nn_0)`$ and the rate

$$w(n|m)=\frac{k_1}{\mathrm{\Omega }}m(Nm)\delta _{n,m+1}+\frac{k_2}{\mathrm{\Omega }}m(Nm)\delta _{n,m1}.$$
(5)
The dynamics of the expectation value of $`n(t)`$ is governed by

$$_tn(t)=\frac{\gamma }{\mathrm{\Omega }}n(t)(Nn(t)),$$
(6)
where $`\gamma =k_1k_2`$. Note that $`n(t)`$ is continuous in $`[0,N]`$. In a spatially extended system the number of particles is a function of position, i.e. $`n=n(x,t)`$. Apart from normalization, $`n(x,t)`$ may be interpreted as the probability of finding an $`A`$ particle in the volume $`\mathrm{\Omega }`$. Dispersion contributes to the change of particles $`_tn(x)`$ according to

$$_tn(x,t)=\text{d}yf(|xy|)\left[n(y,t)n(x,t)\right],$$
(7)
incorporating the probability density rate (Eq. (3)) of jumping from $`y`$ to $`x`$. Eq. (7) defines the operator $``$ acting on the field $`n(x,t)`$. Denoting the spatial density of particles by $`u(x,t)=n(x,t)/\mathrm{\Omega }`$ and combining Eq. (6) with (7) yields

$$_tu=\gamma u(1u)+u,$$
(8)
where we have set without loss of generally the maximum density $`N/\mathrm{\Omega }`$ to unity.

When particles perform ordinary random walks, i.e. if the jump rate $`f(|xy|)`$ is equipped with a length scale $`\sigma `$, the operator $``$ can be approximated by the ordinary Laplacian $`\mathrm{\Delta }`$ on scales larger than $`\sigma `$. However, this description is no longer valid if $`f(x)|x|^{(1+\mu )}`$ with $`0<\mu <2`$. In this case individual jumps lack a scale, particle perform L vy flights, and $``$ is proportional to a non-local singular integral operator,

$$\mathrm{\Delta }^{\mu /2}u(x,t)=C_\mu \text{d}y\frac{\left[u(y,t)u(x,t)\right]}{|xy|^{1+\mu }}.$$
(9)
where $`C_\mu `$ is a constant Brockmann and Sokolov (2002). In Fourier space the operator $`\mathrm{\Delta }^{\mu /2}`$ is equivalent to a multiplication by $`|k|^\mu `$, generalizing the well known $`k^2`$ factor corresponding to the ordinary Laplacian which is why $`\mathrm{\Delta }^{\mu /2}`$ is frequently referred to as a fractional Laplacian.

Since $`u^2>u^2`$ one cannot computed the dynamics of $`n(x,t)`$ from Eq. (8) alone. The crudest approximation is the mean field approach in which fluctuations are ignored ($`u^2u^2`$and thus $`u=u`$) and Eq. (8) reads

$$_tu=\gamma u(1u)+\mathrm{\Delta }^{\mu /2}u,$$
(10)
which is the Fisher-equation (1) for $`\mu =2`$ and its superdiffusive mean field generalization for $`\mu <2`$.

The next better step is to account for fluctuations ($`N1`$) which leads to an additional multiplicative white noise term of variance $`u(1u)/N`$ in Eq. (10). Alternatively one may incorporate the finite particle number as an effective cutoff $`\epsilon 1/N`$ in the logistic growth term retaining a deterministic equation, i.e.

$$_tu=\gamma u(1u)\mathrm{\Theta }(u\epsilon )+\mathrm{\Delta }^{\mu /2}u,$$
(11)
in which $`\mathrm{\Theta }`$ is the Heaviside-function. Qualitatively, this accounts for the fact that no growth on average can occur if the particle concentrations $`u`$ is less than one particle per unit volume. This approach has been applied successfully in the ordinary diffusion scenario Brunet and Derrida (1997). In Eq. (11), mean field theory implies $`\epsilon =0`$.

In the following we investigate front dynamics of the system defined by Eq. (11) and compare its properties to the mean field dynamics ($`\epsilon =0`$) , the limiting case of ordinary diffusion ($`\mu =2`$) and the dynamics of the full probabilistic dynamics (Eqs. (2) and (3)).

Fig. 1 depicts front dynamics obtained from numerical integration of (11) for a concentration $`u(x,t)`$ initially peaked at the origin. For various particle numbers $`N`$ the velocity and shape of the front were computed. Fig. 1(a) shows that even for very large particle numbers wave fronts move asymptotically at constant speeds, a remarkable difference from the mean field limit (dashed line in the inset) which predicts exponentially accelerating front. Fig. 1(b) depicts snapshots of the wavefronts on a double-logarithmic scale. After a transient phase, the shape of the front approaches a steady state with a sharply decreasing boundary at intermediate values for the concentration and an algebraic tail for large $`x`$ with concentrations $`u(x,t)`$ below the cutoff $`\epsilon `$. Fig. 2 displays the steady state front shape on a semilogarithmic scale. As opposed to the algebraic tail predicted by mean field theory, the boundary is an exponential function of relative position $`z=xvt`$, i.e. $`u(z)\mathrm{exp}\lambda z`$, characterized by a spatial scale $`\lambda `$. In summary, the characteristic spatio-temporal wave front solution of Eq. (11) for large times is given by $`u(x,t)\mathrm{exp}\left[\lambda (xvt)\right],`$ for $`1/N<u<1`$ followed by a power-law tail $`u(x,t)(xvt)^{(1+\mu )}`$ for $`u<1/N`$. The decay parameter $`\lambda `$ and the velocity $`v`$ depend on the particle number $`N`$. Qualitatively, the this dependence can be determined in the moving reference frame under the assumption that $`u(x,t)=u(xvt)=u(z)`$. Figs. 1(b) and 2 suggest that $`u(z)=u_1(z)=\mathrm{exp}\lambda z`$ for $`z<z^{}`$ and $`\epsilon <u1`$ and $`u(z)=u_2(z)=A/z^{1+\mu }`$ for $`z>z^{}`$ ($`u<\epsilon `$) where $`z^{}`$ marks the crossover between exponential and algebraic decay. Inserted into Eq. (11) yields

$`v\lambda u_1`$ $``$ $`u_1+\mathrm{\Delta }^{\mu /2}uz<z^{}`$ (12)
$`vu_2^{}`$ $`=`$ $`\mathrm{\Delta }^{\mu /2}uz>z^{}.`$ (13)
Near the crossover the values of the dispersion in both equations are approximately the same, i.e. $`(\mathrm{\Delta }^{\mu /2}u)(z^{}\mathrm{\Delta }z)(\mathrm{\Delta }^{\mu /2}u)(z^{}+\mathrm{\Delta }z)`$ and both equations can be combined to $`v\lambda u_1u_1vu_2^{}`$. With $`u_2^{}(z^{})=(1+\mu )u_2(z^{})/z^{}`$ and $`u_2(z^{})u_1(z^{})`$ we obtain $`v\lambda 1v(1+\mu )/z^{}`$. Since $`z^{}1`$ the second term can be neglected and we find

$$v\lambda ^1.$$
(14)
The dependence of $`\lambda `$ on $`N`$ can be obtained by the short time dynamics. Consider an initial condition $`u(x,0)=\mathrm{exp}(\lambda x)`$ on the half-line $`x0`$. After a short time $`\mathrm{\Delta }t1`$ the wave front is approximately given by $`u(x,\mathrm{\Delta }t)\mathrm{exp}\lambda x+\mathrm{\Delta }t/\lambda x^{1+\mu }.`$ The crossover $`\overline{x}`$ can be defined as the point at which both terms are of the same order of magnitude, implicitly provided by $`\lambda \mathrm{exp}\lambda \overline{x}=\mathrm{\Delta }t/\overline{x}^{1+\mu }`$. Since $`\lambda 1`$ one can approximately solve for the crossover, $`\overline{x}\mathrm{ln}(\lambda ^\mu \mathrm{\Delta }t)/\lambda `$. In order for the exponential to remain invariant under the dynamics, $`\lambda `$ must be chosen such that the crossover coincides with the effective cutoff, i.e. $`u(\overline{x},\mathrm{\Delta }t)=\epsilon `$ which implies the scaling relation $`\lambda \epsilon ^{1/\mu }`$. With Eq. (14) one obtains the scaling law

$$vN^{1/\mu }$$
(15)
for the velocity $`v`$.

Fig. 3 shows the front velocity $`v`$ as a function of particle number $`N`$ obtained by numerical integration of the dynamics (Eq. (11)). The numerics agree well with the scaling law (15) over several orders of magnitude and several choices of the L vy exponent $`\mu `$.

The front propagation characteristics of the macroscopic description in terms of the reaction-superdiffusion equation with effective cutoff (Eq. (11)) coincides with our simulations of the microscopic dynamics defined by Eqs. (2) and (3). The simulation results are summarized in Fig. 4. As in Fig. 1(a) the front velocity $`v`$ is constant after a transient phase, the velocity obeys the scaling law (15) and the wave front decays exponentially in space.

We are convinced that our results are of major importance for the understanding of front propagation in pattern forming systems in which the reactive agents defy the rules of ordinary diffusion. We have shown that constant velocity fronts are typical for pulled front dynamics, contrary to what is expected from mean field approximations and we believe that our results will contribute to the understanding of more complex pattern forming systems such as the geographic spread of human epidemics Brockmann and Hufnagel .

###### Acknowledgements.

This research was supported in part by the National Science Foundation under grant No. PHY99-07949. Discussions with T. Geisel, L. Sander, W. Noyes, V. Hardapple and M. Boone Jr. are greatfully acknowledged.