text:
# Soft gluons and gauge-invariant subtractions in NLO parton-shower Monte Carlo event generators

## 1 Introduction

Current parton-shower Monte Carlo event generators are essentially leading-order (LO) QCD tools. For precision phenomenology at present and future high-energy colliders, it is valuable to be able to go beyond this level of approximation . Although there are a number of treatments of various next-to-leading-order (NLO) effects in event generators (see  and references in ), there is as yet no method for going beyond the leading approximation systematically. This implies that event generators are not able to incorporate fully the known NLO (and NNLO) calculations of hard scattering cross sections.

Recently, systematic subtractive procedures have been proposed  to correct this situation. An important step in the implementation of this program is to show how to decompose Feynman graphs into sums of terms over different regions; the terms are to be arranged to correspond to factors in a factorization theorem that is in a form suitable for the Monte Carlo application. In Ref.  an implementation was presented for a simple, but phenomenologically important case: the photon-gluon fusion process in leptoproduction. This case was simple because soft gluons do not enter at the leading power, so that leading regions do not overlap. The purpose of the present work is to show how to extend the method to decompose graphs with soft gluons and hence with overlapping leading regions. We will treat the simplest such case in leptoproduction, i.e., the photon-quark process to NLO.

A Monte Carlo event generator treats observables that are not infrared and collinear safe. So the kind of subtraction method used in Refs.  and  is not directly applicable, since it relies on the calculated observable being infrared and collinear safe. In particular we cannot use a cancellation of soft gluon contributions between real and virtual graphs. More general factorization theorems are to be used in which the soft contributions factor instead of canceling , and all the factors are defined gauge-invariantly in terms of Wilson-line operators . We will require that the following properties be satisfied:

* The integrand for the hard scattering coefficient is to be an integrable function, and not merely a distribution. Thus the hard scattering coefficient is usable even when one is calculating an infrared and collinear unsafe observable. See Eq. (12) for an example.
* The terms in the expansion of each Feynman graph should arise from matrix elements of gauge-invariant operators.
* In particular, the necessary cut-offs on rapidity integrations are to be defined gauge-invariantly. This involves the use of Wilson-line operators along non-lightlike directions.
* The evolution equations  with respect to these cut-offs should be simple, in the sense that there should be no power-law remainder terms. That is, the equations are homogeneous, like the renormalization-group equations, rather than the Callan-Symanzik equations. This is achieved by the technique proposed in Ref. .

In Sec. 2 we briefly describe the framework provided by the subtraction method for Monte Carlo event generators. In Sec. 3 we treat gauge-invariant subtractions and construct the decomposition of the partonic cross section. In Sec. 4 we comment on graph-by-graph subtractions. In Sec. 5 we give conclusions.

## 2 Subtraction scheme

To put our result in the Monte Carlo context, we schematically represent the cross section in an event generator as

$$\sigma [W]=\underset{\mathrm{final}\mathrm{states}X}{}W(X)\mathrm{PS}\widehat{H}.$$
(1)
Here $`W`$ is a weight function that specifies the definition of the particular cross section under consideration. The symbol $`\mathrm{PS}`$ denotes the parton shower and the symbol $``$ denotes its action on the initial and final partons in the hard scattering, whose cross section is denoted by $`\widehat{H}`$. In a standard Monte Carlo, the hard scattering is taken to the leading order, $`H^{(\mathrm{LO})}`$, and $`\mathrm{PS}`$ denotes showering from the partons in $`H^{(\mathrm{LO})}`$. In an NLO Monte Carlo, the cross section involves a structure of the form

$$\mathrm{PS}\left[H^{(\mathrm{LO})}+\alpha _s\left(H^{(\mathrm{NLO})}\mathrm{PS}_I(1)H^{(\mathrm{LO})}\mathrm{PS}_F(1)H^{(\mathrm{LO})}\right)\right].$$
(2)
Here the first term in the square brackets is the LO hard-scattering function, and the second term is the subtracted NLO hard-scattering function. $`H^{(\mathrm{NLO})}`$ is the result of computing the partonic cross section from the NLO graphs, while $`\mathrm{PS}_I(1)`$ and $`\mathrm{PS}_F(1)`$ are the order $`\alpha _s`$ approximations to the initial-state and final-state showering. The subtraction terms avoid double counting of events already included by showering from $`H^{(\mathrm{LO})}`$.

We will construct a decomposition of $`H^{(\mathrm{NLO})}`$ into a sum of pieces, one for each of the leading regions $`R`$,

$$H^{(\mathrm{NLO})}=\underset{\mathrm{regions}R}{}A_H(R)+\text{nonleading power},$$
(3)
that holds uniformly over the whole of the phase space. Each of the pieces in (3) will contain counterterms that prevent double counting and provide the suppression for going outside the region in which that particular piece was originally supposed to give a good approximation to the matrix element. This is to be contrasted with approaches based on splitting the phase space in different domains and using different approximations to the matrix element in these different domains (see, e.g., Refs. ).

Once we have established a result of the type (3), the term corresponding to the ultraviolet region will give precisely the subtracted hard-scattering function to be used in Eq. (2). The collinear terms will correspond to factors in the cross section that are associated with showering; they therefore imply the evolution kernels to be used in the showering. The soft term would correspond to a new element in the Monte Carlo, but, at the order to which we work, we will show that the soft term can be eliminated by a suitable choice of the cut-offs. This is a result analogous to one in . However the results of  indicate that this is unlikely to be an all-order result.

Our decomposition of graphs of order $`\alpha _s`$ entails a specific definition of the collinear factors, which will not necessarily coincide with the definitions used in any current event generator. However, we will not address in this paper the issue of how to obtain a showering algorithm that corresponds to our collinear terms. A correct solution of this problem will encompass known results about coherent emission of gluons and angular ordering of gluon emission .

## 3 Decomposition of partonic cross section

Let us consider deep inelastic scattering $`\gamma ^{}+PX`$, in which we consider a generic cross section or observable associated with the reaction $`\gamma ^{}(q)+q(p)g(k)+q(k^{})`$ — see Fig. 1. We denote this by $`\mathrm{\Sigma }[\phi ]`$, where $`\phi `$ is a weight function that is the product of the weight function $`W(X)`$ concerned with the final state $`X`$ in Eq. (1) and the factors in the cross section associated with the showering, including the parton density. Thus $`\phi `$ contains all the infrared sensitive and nonperturbative parts of the observable.

We work in the $`\gamma ^{}+\text{hadron}`$ reference frame, and we use light-front coordinates $`v^\mu =(v^+,v^{},𝐯_T)`$ with $`v^\pm =(v^0\pm v^3)/\sqrt{2}`$. The hadron and photon momenta are $`P^\mu =(P^+,m^2/2P^+,\mathrm{𝟎}_T)`$ and $`q^\mu =(xP^+,Q^2/2xP^+,\mathrm{𝟎}_T)`$. Then we parameterize the gluon momentum $`k`$ as

$$k^\mu =(\alpha xP^+,\beta \frac{Q^2}{2xP^+},\left|𝐤_T\right|\widehat{\mathit{\varphi }}),$$
(4)
where $`\widehat{\mathit{\varphi }}`$ is a unit transverse vector at azimuthal angle $`\varphi `$.

For our calculation, the external partons are on-shell, and the incoming quark $`p`$ has zero transverse momentum, so that $`\mathrm{\Sigma }`$ can be written as follows:

$$\mathrm{\Sigma }[\phi ]=_0^{\mathrm{}}𝑑\alpha _0^{\mathrm{}}𝑑\beta _0^{2\pi }\frac{d\varphi }{2\pi }\phi (x,Q^2,\alpha ,\beta ,\varphi )J(x,\alpha ,\beta )(\alpha ,\beta ).$$
(5)
Here, $`J`$ is the Jacobian factor

$$J(x,\alpha ,\beta )=\frac{1}{16\pi ^2}\frac{1}{1+\alpha \beta }\mathrm{\Theta }\left(\frac{1x}{x}\alpha \right)\mathrm{\Theta }\left(1\frac{x}{1x}\alpha \beta \right),$$
(6)
and $``$ is the next-to-leading-order matrix element for $`\gamma ^{}q`$ obtained by contracting the photon Lorentz indices with the projector corresponding to the structure function $`F_2`$

$$=4e_q^2g_s^2C_FM(\alpha ,\beta ),M(\alpha ,\beta )=(1\beta )^2\frac{1+(1+\alpha \beta )^2}{\alpha \beta (1+\alpha \beta )}+2+6\frac{(1\beta )^2}{1+\alpha \beta }.$$
(7)
The physical region for $`\alpha ,\beta `$ is the interior of the triangle in Fig. 2.

Standard arguments  determine the infrared sensitive regions contributing to the leading power behavior of $`\mathrm{\Sigma }[\phi ]`$, which are located on Fig. 2 as follows: The region in which the gluon is collinear to the initial state is a neighborhood of the axis $`\beta =0`$, the region in which the gluon is collinear to the final state is a neighborhood of the axis $`\alpha =0`$, and the soft region is a neighborhood of the origin $`\alpha =0`$, $`\beta =0`$. The truly hard region lies away from the $`\alpha =0`$ and $`\beta =0`$ axes.

To obtain a decomposition for $`\mathrm{\Sigma }`$ of the type of Eq. (3), we now employ the technique of Ref. . This generalizes the R-operation techniques of renormalization. (See Ref.  for a related approach.) To ensure that the procedure is gauge-invariant, each of the terms in the right hand side of Eq. (3) is constructed from matrix elements involving Wilson line operators,

$$V_I(n)=𝒫\mathrm{exp}\left(ig_{\mathrm{}}^0𝑑ynA(yn)\right),V_F(n)=𝒫\mathrm{exp}\left(ig_0^+\mathrm{}𝑑ynA(yn)\right),$$
(8)
with suitable directions $`n`$ for the lines. Evolution equations in $`n`$ enable one to connect the results corresponding to different directions . We define light-like vectors $`\widehat{p}=(1,0,\mathrm{𝟎}_T)`$, $`\widehat{p}^{}=(0,1,\mathrm{𝟎}_T)`$. We will also use non-lightlike vectors $`u=(u^+,u^{},\mathrm{𝟎}_T)`$, $`u^{}=(u^+,u^{},\mathrm{𝟎}_T)`$, all of whose components are positive. It is convenient to define $`\eta =(2x^2P_{}^{+}{}_{}{}^{2}/Q^2)u^{}/u^+`$, and $`\eta ^{}=(Q^2/2x^2P_{}^{+}{}_{}{}^{2})u_{}^{}{}_{}{}^{+}/u_{}^{}{}_{}{}^{}`$.

As in , we start with the smallest region, the soft region $`\alpha ,\beta 0`$, and determine the corresponding contribution to the matrix element (7):

$$M_S(\alpha ,\beta )=\frac{2}{\alpha \beta }\frac{2}{(\alpha +\eta ^{}\beta )\beta }\frac{2}{\alpha (\beta +\eta \alpha )}.$$
(9)
Observe that the first term in the right hand side of this formula is just obtained by taking the soft approximation to Eq. (7). It can be thought of as the one-loop contribution to the square of a vacuum–to–gluon matrix element of a product of eikonal Wilson lines taken along lightlike directions $`\widehat{p}`$, $`\widehat{p}^{}`$ . This first term reproduces the behavior of the matrix element $`M`$ when $`\alpha `$ and $`\beta `$ simultaneously approach zero. But there are also logarithms in its integral associated with the collinear regions where $`\alpha /\beta `$ or $`\beta /\alpha `$ go to zero. The subtractions provided by the other two terms conveniently cancel these regions. They can be derived from operators analogous to those for the first term, except for replacing one of the lightlike eikonal lines by a line along a non-lightlike direction. In particular, the second term subtracts the divergence from a region collinear to the initial state, i.e., from the region $`\beta /\alpha 0`$. At the same time, the non-lightlike vector $`u^{}`$ in this term provides a cut-off on the region of small $`\alpha `$. Similarly, the third term in Eq. (9) subtracts the divergence from the region collinear to the final state, i.e., from the region $`\alpha /\beta 0`$.

Next we construct terms for the collinear regions. By applying a treatment  analogous to that for the soft region, we arrive at

$$M_I(\alpha ,\beta )=\frac{1}{\beta }\frac{1+(1+\alpha )^2}{\alpha (1+\alpha )}\frac{2}{\alpha \beta }+\frac{2}{(\alpha +\eta ^{}\beta )\beta }=\frac{1}{\beta }\left(\frac{\alpha }{1+\alpha }+\frac{2}{\alpha +\eta ^{}\beta }\right)$$
(10)
for the region collinear to the initial state, and

$`M_F(\alpha ,\beta )`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{(1\beta )+(1\beta )^3}{\beta }}{\displaystyle \frac{2}{\alpha \beta }}+{\displaystyle \frac{2}{\alpha (\beta +\eta \alpha )}}`$ (11)
$`=`$ $`{\displaystyle \frac{1}{\alpha }}(4+3\beta \beta ^2+{\displaystyle \frac{2}{\beta +\eta \alpha }})`$
for the region collinear to the final state. The first term of the expression in the middle in each of these equations is the unsubtracted collinear approximation to the original matrix element (i.e., the $`\beta 0`$ or $`\alpha 0`$ limit of $`M`$). We will comment below on the subtraction terms.

The fully subtracted matrix element, associated with the hard region, is then given by

$`M_H(\alpha ,\beta )`$ $`=`$ $`MM_SM_IM_F`$ (12)
$`=`$ $`(1\beta )^2{\displaystyle \frac{1+(1+\alpha \beta )^2}{\alpha \beta (1+\alpha \beta )}}+2+6{\displaystyle \frac{(1\beta )^2}{1+\alpha \beta }}`$
$`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{1+(1+\alpha )^2}{\alpha (1+\alpha )}}{\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{1+(1\beta )^2}{\beta }}(1\beta )+{\displaystyle \frac{2}{\alpha \beta }}`$
$`=`$ $`\beta +{\displaystyle \frac{\alpha }{(1+\alpha )(1+\alpha \beta )}}+{\displaystyle \frac{6(1\beta )^2}{(1+\alpha \beta )}}.`$
This matrix element is finite in all of the infrared regions. It can be safely integrated down to $`\alpha =0`$ or $`\beta =0`$. Moreover, it is independent of the choice of the non-lightlike directions $`u`$, $`u^{}`$: the dependence on $`\eta `$, $`\eta ^{}`$ has canceled in Eq. (12).

Eqs. (9)–(12) achieve a decomposition of the type (3) for the original matrix element. There is one term for each of the leading regions — in particular, a soft term. We now ask: can we reorganize it in a way that is suited for a parton-shower algorithm, such as, e.g., the algorithm  used in the event generators ?

The soft term can be eliminated by choosing the vectors $`u`$ and $`u^{}`$ so that $`\eta \eta ^{}=1`$; then there are only collinear terms, as is appropriate to match the structure of the parton-shower Monte Carlo algorithm, which has only initial-state or final-state branchings. The symmetric choice $`\eta =\eta ^{}=1`$ gives

$$M_I^{(\mathrm{MC})}(\alpha ,\beta )=\frac{1}{\beta }\frac{1+(1+\alpha )^2}{\alpha (1+\alpha )}\frac{1}{\alpha }\frac{2}{\alpha +\beta },$$
(13)
$$M_F^{(\mathrm{MC})}(\alpha ,\beta )=\frac{1}{\alpha }\frac{1+(1\beta )^2}{\beta }(1\beta )\frac{1}{\beta }\frac{2}{\alpha +\beta }.$$
(14)
All of the infrared contributions are now associated with configurations that are either collinear to the initial state or collinear to the final state. We have inserted superscripts in the left hand sides of Eqs. (13) and (14) to indicate that this particular choice of the non-lightlike directions gives rise to a structure that corresponds to that of Monte Carlo algorithms.

Eq. (12) determines the form of the counterterms to be used in the matrix element; Eqs. (13) and (14) determine the form of the counterterms to be used in the showering. The subtracted NLO matrix element (12) only receives contributions from the truly ultraviolet region. As for the modification to the showering, consider Eq. (13). $`M_I^{(\mathrm{MC})}`$ is associated with the first branching from the initial state. As noted below Eq. (11), the first term in the right hand side is obtained from taking the collinear approximation $`\beta 0`$ to the original matrix element. The coefficient of $`1/\beta `$ is the standard quark $``$ quark splitting kernel:

$$\frac{1+(1+\alpha )^2}{\alpha (1+\alpha )}=P_{qq}\left(\frac{1}{1+\alpha }\right).$$
(15)
This first term corresponds to the standard form of the showering. It gives a good approximation in the initial-state collinear region, i.e., $`\beta 0`$. The second term in the right hand side of Eq. (13), on the other hand, is non-standard. It effectively provides a cut-off when $`\alpha 0`$. Note that the second term is suppressed in the collinear region $`\beta 0`$ at fixed $`\alpha `$. That is, the modified showering coincides with the usual one in the collinear limit and differs from it away from the collinear limit. Analogous remarks can be made based on the formula (14) associated with the final state.

Observe that if one regulated the $`\alpha 0`$ behavior of the first term in Eq. (13) by subtracting its $`\alpha 0`$ limit, given by $`2/(\alpha \beta )`$, this would bring about an extra $`\beta 0`$ singularity. This would not be suited for our application in a Monte Carlo algorithm. In contrast, the second term in Eq. (13) represents precisely what is, from our point of view, a better choice of a counterterm: it subtracts the $`\alpha 0`$ singularity without introducing any extra singular behavior at $`\beta 0`$.

Note also that this counterterm cuts off the integration over the region of small $`\alpha `$ at a value of order $`\beta `$: the leading behavior of $`M_I^{(\mathrm{MC})}`$ for small $`\alpha `$ is of the type

$$\frac{2}{\beta (\alpha +\beta )}+\text{regular}\text{terms}\text{in}\alpha .$$
(16)
Then we see that the procedure based on gauge-invariant subtractions that we have just applied, compared to the cut-off method, tells us precisely where the cut-off is to be placed. The position of the cut-off on $`\alpha `$ turns out to be $`\beta `$-dependent. In more physical terms, this indicates that the cut-off to be applied in the initial-state shower and the cut-off to be applied in the final-state shower are not to be set independently, but they are related.

## 4 Graph-by-graph subtractions

In our calculational example, the form of the four terms was determined by requiring them to be obtainable from matrix elements of gauge-invariant operators. However, the procedure can be applied graph-by-graph, as we will now explain. The basis of this procedure is the derivation of factorization for soft factors given in Refs. .

Consider the connection of a gluon to a subgraph that consists of lines that are all collinear in the $`+`$ direction. This factor we denote by $`J^\mu (k)`$, where $`\mu `$ is the Lorentz index that couples to the gluon. In a region where the gluon is either soft or collinear in the opposite direction, it is a good approximation to replace $`J^\mu `$ by its $`+`$ component, and to replace the momentum $`k`$ in $`J`$ by its $``$ component. Multiplying and dividing by $`k^{}`$ gives

$`J^\mu (k)`$ $``$ $`g_+^\mu J^+(0,k^{},\mathrm{𝟎}_T)`$ (17)
$`=`$ $`g_+^\mu {\displaystyle \frac{1}{k^{}}}\left[k^{}J^+(0,k^{},\mathrm{𝟎}_T)\right].`$
Since the last factor is of the form of a gluon Green function contracted with the gluon’s momentum, a Ward identity can be applied. Then after a sum over all relevant graphs we obtain factors corresponding to matrix elements of operators that include path-ordered exponentials .

When this replacement is made, together with the corresponding replacement for gluonic connections to the opposite “jet” factor, some of the terms have new divergences when the gluon’s rapidity goes to $`+\mathrm{}`$ or $`\mathrm{}`$, as explained earlier and in Ref. . These divergences are cancelled by further subtractions that are constructed by manipulations on the $`1/k^{}`$ and $`1/k^+`$ factors, to give a result like Eq. (9).

## 5 Conclusions

In conclusion, the above calculation tells us how to organize the subtractions in situations in which soft gluons are present and both initial-state and final-state branchings contribute. This is one of the issues that have to be dealt with to construct NLO Monte Carlo event generators. The subtractive technique  used for this calculation is based on graph-by-graph R-operation methods and enforces gauge invariance by relating the counterterms to matrix elements of Wilson-line operators. The calculations of this paper do not address the issue of how to generate the whole shower corresponding to the subtracted collinear terms. The investigation of this is left to future work.

###### Acknowledgments.

This research is supported in part by the US Department of Energy under grant No. DE-FG02-90ER-40577.


format: text

text:
# Antiflow of kaons in relativistic heavy ion collisions

## Abstract

We compare relativistic transport model calculations to recent data on the sideward flow of neutral strange $`K_s^0`$ mesons for Au+Au collisions at 6 AGeV. A soft nuclear equation of state is found to describe very well the positive proton flow data measured in the same experiment. In the absence of kaon potential, the $`K^0`$ flow pattern is similar to that of protons. The kaon flow becomes negative if a repulsive kaon potential determined from the impulse approximation is introduced. However, this potential underestimates the data which exhibits larger antiflow. An excellent agreement with the data is obtained when a relativistic scalar-vector kaon potential, that has stronger density dependence, is used. We further find that the transverse momentum dependence of directed and elliptic flow is quite sensitive to the kaon potential in dense matter.

The study of hadron properties in hot and dense medium created in relativistic heavy ion collision is of considerable interest. Of particular importance is the in-medium modification of kaon properties as it is related to chiral symmetry restoration  and neutron star properties . Since the work of Kaplan and Nelson  on the possibility of kaon condensation in dense matter, considerable effort has been devoted to the understanding of kaon properties in dense matter from heavy ion collisions \[4-10\]. All these investigations have lead to a general consensus that kaons feel a weak repulsive potential while antikaons feel a strong attractive potential ; the latter is responsible for $`K^{}`$ condensation in neutron stars. Since kaons are primarily produced at an early hot and dense stage in heavy ion collision, a promising tool to search for kaon in-medium properties is via pressure-induced collective flow effects. The strength of the flow depends on the nuclear equation of state (EOS), the pressure gradient developed, and the kaon potential or dispersion relation in dense matter. Indeed, Li et al. , based on a relativistic transport model, have found that the repulsive kaon potential in nuclear medium tends to deflect kaons away from nucleons leading to an anticorrelated flow. On the other hand, antikaons were found  to have a similar flow as that of nucleons due to their attractive potential.

The theoretical predictions of exploiting the transverse flow to study the kaon self-energy in the medium have attracted much attention experimentally. The FOPI Collaboration at SIS/GSI has found  for the reaction Ni+Ni at $`E_{\mathrm{beam}}=1.93`$ AGeV a very small antiflow for kaon, i.e., the kaon average transverse momentum in the reaction plane as a function of rapidity is almost zero around midrapidity. The kaon sideward flow for this reaction, and especially for a heavier system Ru(1.69 AGeV)+Ru, was found to be anticorrelated (correlated) with that of protons at low (high) transverse momenta. The azimuthal angular distributions of kaons measured by the KaoS Collaboration  in Au+Au collisions at 1 AGeV indicate that kaons are emitted preferentially perpendicular to the reaction plane. Also large $`K^{}`$ production cross-section was observed by this group in Ni+Ni collisions at $`E_{\mathrm{beam}}=(0.81.8)`$ AGeV . All these experimental findings at beam energies of $`12`$ AGeV clearly suggest the existence of an in-medium repulsive kaon-nucleon potential and a strong attractive antikaon-nucleon potential; the latter is also supported by $`K^{}`$ atomic data .

Flow analysis for $`K^+`$, $`K^{}`$, and $`K_s^0`$ have been performed by several collaborations at AGS/BNL in Au+Au collision at beam energies of $`212`$ AGeV . Most of these experiments are confined to near-central collisions where nearly vanishing kaon transverse flow versus rapidity was observed. Several suggestions put forward about the origin of vanishing flow are, for example, isotropic production of kaons in hadron-hadron collisions ; the colliding hadrons have opposite rapidities ; cancellation between the negative flow at low transverse momenta and the positive flow at high transverse momenta .

Recently, the E895 Collaboration at AGS/BNL has measured  the directed flow of neutral strange $`K_s^0`$ mesons in 6 AGeV Au+Au collisions for central and mid-central events. The $`K_s^0`$ was found to have a considerable antiflow relative to that of the proton observed in the same experiment . Using a relativistic transport model (ART)  for heavy ion collisions we investigate in this paper if the appreciable kaon antiflow can be explained by the repulsive kaon potential in dense nuclear matter. We shall demonstrate that the typical relativistic scalar-vector mean field potential can indeed reproduce the data. We also predict that the flow anisotropies is sensitive to the kaon potential and is thus a useful probe of the EOS in high density matter.

In the ART model employed here, the nucleon mean field is parameterized by the usual Skyrme-type with a soft EOS corresponding to an incompressibility of $`K=210`$ MeV and a stiff EOS with $`K=380`$ MeV at normal nuclear matter density of $`\rho _0=0.16`$ fm-3. In this model, the imaginary part of the kaon self-energy is approximately treated by kaon-hadron scatterings and its real part is given by the mean-field potential. The model has been used to study many aspects of heavy ion collisions at the AGS energies . In the original ART model, kaon production in hadron-hadron scatterings is instantaneous. In the present study, we introduce a formation time of $`\tau _{\mathrm{fo}}^K=1.2`$ fm/c.

Various approaches have been adopted to evaluate the kaon potential in dense matter . In the present study, we shall use the two commonly used forms that have been quite successful in explaining several experimental data. One of these is based on the kaon dispersion relation determined from the kaon-nucleon scattering length using the impulse approximation:

$$w_K(𝐩,\rho _b)=\left[m_K^2+𝐩^24\pi \left(1+\frac{m_K}{m_N}\right)\overline{a}_{KN}\rho _b\right]^{1/2},$$
(1)
where $`m_K`$ and $`m_N`$ are the kaon and nucleon bare masses, $`\rho _b`$ is the baryon density and $`\overline{a}_{KN}0.255`$ fm is the isospin-averaged kaon-nucleon scattering length. The kaon potential in nuclear medium can then be defined as

$$U_K(𝐩,\rho _b)=w_K(𝐩,\rho _b)\left(m_K^2+𝐩^2\right)^{1/2}.$$
(2)
This yields a repulsive potential of 30 MeV at $`\rho _0`$ for kaon at zero momentum. The other form of kaon potential used here is the scalar-vector potential determined from the chiral Lagrangian and can be written as

$$w_K(𝐩,\rho _b)=\left[m_K^2+𝐩^2a_K\rho _S+(b_K\rho _b)^2\right]^{1/2}+b_K\rho _b,$$
(3)
where $`b_K=3/(8f_\pi ^2)0.333`$ GeV fm3, $`\rho _S`$ is the scalar density. If only the Kaplan-Nelson term, which is related to the $`KN`$ sigma term due to explicit chiral symmetry breaking, is considered, then the parameter $`a_K`$ which determines the strength of the attractive kaon potential is $`a_K=\mathrm{\Sigma }_{KN}/f_\pi ^2`$. Since the exact value of $`\mathrm{\Sigma }_{KN}`$ and the magnitude of the higher order corrections are poorly known, we take a different approach where $`a_K0.173`$ GeV2 fm3 is determined by fitting the kaon repulsive potential of $`U_K(\rho _0)=+30`$ MeV as obtained in the impulse approximation.

Let us now make a systematic study of the directed flow of particles measured by the E895 Collaboration for the Au+Au collisions at $`E_{\mathrm{beam}}=6`$ AGeV . Before we confront with the $`K_s^0`$ flow data, it is instructive to find how far the ART model can describe the observed proton flow . Calculations are performed for impact parameters $`b=57`$ fm, since the events analyzed for the proton flow originate from this domain . Figure 1 shows for protons the mean transverse momentum in the reaction plane $`p_x`$ as a function of rapidity $`y`$ (scaled to the beam rapidity). The data reveals a large proton flow with a slope of $`F=dp_x/d(y/y_{c.m.})180`$ MeV/c at midrapidity. The ART model with a soft nuclear equation of state (solid line) provides a good agreement with the data over a large rapidity range. On the other hand,

FIG. 1. The proton directed transverse flow $`p_x`$ as a function of normalized rapidity $`y/y_{cm}`$ in Au+Au collisions at $`E_{\mathrm{beam}}=6`$ AGeV. The ART model calculations for a soft (solid line) and stiff (dashed line) nuclear EOS at an impact parameter $`b=57`$ fm are compared with the E895 data (filled circles).

the stiff EOS (dashed line) results in a much stronger proton flow compared to the data. The enhanced pressure gradient for the stiff EOS provides a larger driving force for the expansion of the system and thereby a stronger collective flow. Henceforth, we restrict only to the soft nuclear EOS when comparing with the $`K_s^0`$ antiflow data measured in the same experiment.

In Fig. 2, the ART model predictions of the sideward kaon flow as a function of rapidity are compared with the E895 data of $`K_s^0`$ mesons reconstructed from charged pion decays . Since the antikaon $`\overline{K}^0`$ has a contribution of $`\stackrel{<}{}10\%`$ in the data sample at AGS energies, the observed $`K_s^0`$’s are therefore primarily from neutral kaons $`K^0`$’s. It is thus justified to study the collective flow of $`K^0`$ mesons in our model. In contrast to the protons, a pronounced in-plane antiflow for the $`K_s^0`$ is observed in the data with a slope of $`F127\pm 20`$ MeV/c. The $`K_s^0`$ flow data corresponds to central and mid-central events with $`b\stackrel{<}{}7`$ fm and for transverse momentum $`p_t700`$ MeV/c . We therefore restrict our calculations to these cut-off values. As is evident from Fig. 2, in the absence of kaon mean field, the $`K^0`$’s have a flow pattern similar to that of the nucleons. This is not surprising as most of the $`K^0`$’s are produced in the early compression stage of the reaction at time $`t\stackrel{<}{}5`$ fm/c and the rescattering of the kaons with the nucleons in the dense matter thus causes them to flow in the direction of the nucleons.

We shall now demonstrate the effect of inclusion of kaon potential on the $`K^0`$ flow. Using the kaon potential determined from the impulse approximation, the $`K^0`$’s are repelled from the nucleons resulting in antiflow with respect to the nucleons (see Fig. 2). However, the $`p_x`$ for kaon is found to underestimate the data as in Ref. .

FIG. 2. The kaon directed transverse flow $`p_x`$ as a function of normalized rapidity $`y/y_{cm}`$ in Au+Au collisions at $`E_{\mathrm{beam}}=6`$ AGeV. The filled circles are the E895 data while the other curves correspond to ART model calculations ($`b\stackrel{<}{}7`$ fm and $`p_t700`$ MeV/c) without kaon mean field potential and with kaon dispersion relations obtained from impulse approximation and scalar-vector potential. For the latter kaon potential, results with zero kaon formation time $`\tau _{\mathrm{fo}}^K`$ is also shown. The inset shows the density dependence of the kaon potential $`U_K`$ at zero momentum obtained in the impulse approximation and scalar-vector potential.

In fact, theoretical calculations  have shown that the simple impulse approximation also underestimates experimental kaon-nucleus scattering data. The scalar-vector potential has clearly a stronger density dependence compared to the impulse approximation (see inset of Fig. 2) since the vector potential dominates over the scalar potential at densities above the normal nuclear matter value. The kaon flow predicted from the scalar-vector potential is found to have a strikingly good agreement with the $`K_s^0`$ flow data as seen in Fig. 2. Although the final kaon flow is opposite to the proton flow, the primordial kaons flow with the nucleons up to the maximum compression stage of the collision, and their $`p_x`$ values are nearly identical irrespective of the kaon potential employed. This stems from the fact that at the early stage of the collision the density gradient inside the participant matter is rather small and as a consequence the force acting on the kaon during its propagation is negligible. The potential (or baryonic density) gradient is substantial only during the expansion of the system so that the kaon potential could then effectively repel the $`K^0`$’s away from the nucleons. This suggests that at the AGS energies, kaon flow may provide information about the kaon dispersion relation for densities below $`4\rho _0`$. Since the kaons produced in the ultradense stage mostly travel with nucleons, the kaon flow is expected to be sensitive to its formation time. This effect is illustrated in Fig. 2 for the scalar-vector potential where the kaon formation time $`\tau _{\mathrm{fo}}^K`$ is set to zero. The early production of $`K^0`$ then leads to larger rescattering in the initial stage of the reaction. Consequently, the antiflow of $`K^0`$’s is suppressed especially near the midrapidity. It is worth mentioning that the scalar-vector potential has been successful in reproducing several data at GSI energies of $`12`$ AGeV .

Additional insight into the origin of the transverse flow may be gained by considering the azimuthal angular distribution $`dN/d\varphi `$ of the kaons at midrapidity by means of the Fourier expansion

$$\frac{dN}{d\varphi }=v_0\left[1+2v_1\mathrm{cos}\varphi +2v_2\mathrm{cos}(2\varphi )+\mathrm{}\right].$$
(4)
Here $`v_0`$ is the normalization constant, and the Fourier coefficients $`v_1=\mathrm{cos}\varphi =p_x/p_t`$ represents the directed flow and $`v_2=\mathrm{cos}(2\varphi )=(p_x^2p_y^2)/p_t^2`$, referred to as the elliptic flow , reflects the azimuthal asymmetry of the emitted particles. While $`v_2>0`$ indicates in-plane enhancement, $`v_2<0`$ suggests a squeeze-out orthogonal to the reaction plane.

In Fig. 3 we depict the $`p_t`$-dependence of $`v_1`$, dubbed as differential flow, for midrapidity $`K^0`$ ($`|y/y_{c.m.}|<0.3`$) without any kaon mean field and with the inclusion of scalar-vector kaon potential. The results are for semicentral ($`b=6`$ fm) Au+Au collisions at 6 AGeV which can describe equally well the flow data. In absence of kaon potential, $`v_1`$ is found to be almost zero at low $`p_t`$ but has a large positive value at high $`p_t`$. The low-$`p_t`$ kaons are produced mostly in the relatively dilute region of the central expanding hadronic matter and are preferentially emitted in the antiflow direction away from the spectators. Consequently, they freeze out quite early in the reaction leading to a small in-plane flow. In contrast, kaons with large final $`p_t`$ originate mainly inside the baryon-rich participant zone. These kaons are likely to suffer multiple rescatterings with the participant baryons and with the spectators. Hence the transverse momenta of these kaons are pushed to higher values with the final flow direction (and magnitude) similar to those of the nucleons. The differential flow for midrapidity nucleons seen here is similar to that observed experimentally .

In presence of repulsive scalar-vector potential for $`K^0`$, the differential flow of kaons with $`p_t1.2`$ GeV/c becomes negative while kaons with higher $`p_t`$ remains positive. The large negative $`v_1`$ for low-$`p_t`$ kaons follows from the fact that these kaons are formed near the surface of the participant matter and thus feel a pronounced repulsive (kaon) potential gradient away from the baryons. Also the force acting on a kaon being inversely proportional to energy, the low energy kaons are strongly repelled from the baryons. Since the high-$`p_t`$ kaons are produced deep inside the baryon-rich medium they are less repelled by the kaon potential both due to smaller density gradient and eventually because of their larger energy. These kaons thus remain to flow in the same direction as the baryons. Note that the $`p_t`$-cut in the E895 data for $`K_s^0`$ flow ($`p_t0.7`$ GeV/c) leads to enhanced

FIG. 3. The transverse momentum dependence of the directed flow $`v_1`$ for midrapidity protons and kaons without and with a kaon scalar-vector potential. The results are for Au(6 AGeV)+Au collisions at an impact parameter of 6 fm.

antiflow for the $`K_s^0`$ than that when integrated over the entire $`p_t`$ spectra.

In Au(6 AGeV)+Au collisions the passage time of projectile and target spectator is given by $`t2R/(\gamma v)6.4`$ fm/c. During this time, the spectators prevent the participating hadrons from developing an in-plane flow. Therefore the hadronic matter is initially squeezed-out preferentially orthogonal to the reaction plane. The sign and magnitude of elliptic flow $`v_2`$, apart from the passage time, is also determined by the transverse expansion time of the participant matter. In the later stages of the reaction, however, the geometry of the participant region is such that in-plane emission is favored resulting in a positive $`v_2`$ . Most of the kaons at the AGS energies are produced during the compression stage at times less than $`5`$ fm/c and thus would be shadowed by the spectators. Therefore kaons may serve as a sensitive probe to the EOS in the ultradense matter. Indeed, in peripheral and semicentral Au+Au collisions at $`1`$ AGeV, preferential emission of $`K^+`$ meson perpendicular to the reaction plane has been observed . This may be interpreted as an evidence of repulsive kaon potential in the medium .

In Fig. 4, the elliptic flow $`v_2`$ as a function of $`p_t`$ is shown for protons and kaons using the same cuts as in Fig. 3. The proton $`v_2`$ gradually increases with $`p_t`$ and have positive values even for small $`p_t`$. The low-$`p_t`$ particles generally freezes out in the compression stage orthogonal to the reaction plane due to spectator shadowing, while the high-$`p_t`$ particles undergo multiple rescattering. The experimental observation  for proton $`v_2`$ indicates a transition from squeeze-out to a preferential in-plane emission with increasing beam energies; the transition occurs at $`E_{\mathrm{beam}}6`$ AGeV. The ART model  and the UrQMD model  can describe the data consistently with a stiff nuclear EOS ($`K=380`$ MeV); the $`v_2`$ for protons

FIG. 4. Same as Fig. 3, but for the elliptic flow $`v_2`$.

in the ART model is however insensitive to the EOS for $`E_{\mathrm{beam}}6`$ AGeV. This is in contrast to the prediction of Ref.  where a softening of the EOS with incident energy is deduced.

The elliptic flow for $`K^0`$ in absence of any kaon potential is found to be slightly larger than that of the protons. This is caused by smaller shadowing by the spectators because of the time delay in their production and due to small $`KN`$ total cross-section of about 10 mb. With the inclusion of the scalar-vector potential, the $`v_2`$ for kaons is reduced. An explanation to this is as follows. The geometrical shape of the participant region at the initial stages entails a kaon repulsive potential gradient that is larger in the direction perpendicular to the reaction plane than parallel to the plane. Consequently, kaons in its mean field escape more freely perpendicular to the reaction plane. This effect is more pronounced for low-$`p_t`$ kaons produced near the surface of the participant region. We however find that at the AGS energies considered here, the sensitivity of the squeeze-out to the kaon dispersion relation is smaller compared to that observed at GSI energies .

In summary, recent data on the sideward flow of neutral strange $`K_s^0`$ meson in 6 AGeV Au+Au collisions have been compared to ART model calculations. The directed flow data for $`K_s^0`$ is opposite to that of protons and reveals a large antiflow. A soft nuclear EOS is found to reproduce very well the in-plane proton flow data at the same beam energy. We have used two different kaon dispersion relations: the impulse approximation where the mean field has a linear density dependence, and the typical relativistic scalar-vector potential. The $`K^0`$ mesons which flow with the nucleons in absence of a kaon mean-field potential exhibit an antiflow for the kaon potential derived from impulse approximation. However, the antiflow is found to underestimate the data. An overall good agreement with the data is achieved for the scalar-vector potential that has been equally successful in explaining the kaon data at SIS/GSI energies. Since the vector potential dominates at high densities, the scalar-vector potential becomes more repulsive than that in the impulse approximation. The $`K^0`$ flow at AGS energies is sensitive to the kaon potential at densities smaller than about four times the normal nuclear matter value as the kaon flow is essentially determined during the expansion stage of the collision. At higher densities reached at the early stages, the kaons generally flow with the nucleons and their flow magnitude is thus determined by the nuclear EOS. Kaons with small transverse momentum have a large negative differential flow due to stronger repulsion by the kaon potential, while positive kaon differential flow is a consequence of strong positive baryonic flow. In the semicentral collisions about midrapidity, the kaons are found to have an in-plane elliptic flow. The magnitude of this flow is slightly larger than that of the protons since the kaons suffer smaller shadowing. In presence of scalar-vector kaon potential, the kaon elliptic flow is decreased. This suggests that analysis of directed flow in conjunction with elliptic flow of kaons may provide useful information on the kaon dispersion relation in dense matter.

This work was supported in part by the National Science Foundation under Grant No. PHY-9870038, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081.


text:
# Breathing and randomly walking pulses in a semilinear Ginzburg-Landau system

## I Introduction

The important role played by localized pulses (sometimes called “autosolitons” ) in models of the Ginzburg-Landau (GL) type  is well known. They find applications to plasmas and plasma-like media , -, nonequilibrium semiconductors , hydrodynamics (Poiseuille flow , electrohydrodynamic convection ), nonlinear optics - , etc.

The simplest model which gives rise to localized pulses is the complex cubic GL equation,

$$u_t=\gamma _0u\left(gi\sigma \right)|u|^2u+(\gamma _1+i\gamma _2)u_{xx},$$
(1)
where $`\gamma _0`$, $`g`$, and $`\gamma _1`$ are positive coefficients of the linear gain, nonlinear losses, and effective diffusion (dispersive losses), while $`\sigma `$ and $`\gamma _2`$ (which may have either sign) account for the nonlinear frequency shift and dispersion, respectively. Equation (1) always has a single exact solitary-pulse (SP) solution , which is, however, unstable, as its background, i.e., the trivial solution $`u=0`$, is obviously unstable due to the presence of the linear gain.

The simplest possibility to modify the model so that to produce stable SPs is to convert it into the quintic GL equation , setting $`\gamma _0<0`$ and $`g<0`$, and adding a nonlinear quintic dissipative term in order to provide for the overall stabilization. In the quintic GL equation, the trivial solution is stable since $`\gamma _0<0`$, and two solitary-pulse solutions, one unstable and one stable, may coexist at fixed values of parameters .

In this connection, it is necessary to mention a modification of the usual cubic GL equation proposed by Schöpf and Kramer  (see also the works ), in which $`\gamma _0=1`$ and $`g`$ is negative, while no higher-order stabilizing nonlinearity is added. The corresponding model lacks global stability, as a very large perturbation with a small wavenumber will provoke blowup. Nevertheless, numerical simulations have shown that finite-amplitude solutions in the form of both chaotic and stationary arrays of pulses may exist in this model, due to the stabilizing effect exerted by the dispersion and diffusion.

The quintic GL equation is, generally, a phenomenological model, as the derivation from the first principles can scarcely stop at the fifth-order nonlinearity (provided that it is a crucially important term rather than a small correction). Although this model can sometimes describe experimentally observed waves with a fairly high accuracy, an example being the so-called dispersive chaos in the traveling-wave convection in binary fluids , it is desirable to find more models that allow for the existence of stable SPs and can be derived from the first principles. A possibility is to consider a model of a dual-core nonlinear optical fiber, in which the linear gain, dispersion, effective diffusion (spectral filtering), and Kerr (nondissipative cubic) nonlinearity are present in one (active) core, while the other (passive) one, linearly coupled to the active core, has only linear loss . Thus, the model consists of a cubic GL equation linearly coupled to the second, purely linear, ordinary differential equation for the local amplitudes $`u`$ and $`v`$ of the waves in the two subsystems:

$`u_t`$ $`=`$ $`\gamma _0u+i\sigma |u|^2u+(\gamma _1+i\gamma _2)u_{xx}+iv,`$ (2)
$`v_t`$ $`=`$ $`\mathrm{\Gamma }v+iu`$ (3)
(note that, in the application to optical fibers, the evolutional variable $`t`$ is actually the distance along the fiber, while $`x`$ is the so-called local time ). The parameters $`\gamma _0`$, $`\gamma _1`$, $`\gamma _2`$ and $`\sigma `$ of the active subsystem have the same meaning as in the cubic GL equation (1), there is no saturating nonlinear term, $`\mathrm{\Gamma }`$ is the loss factor in the passive subsystem, and the coefficient of the linear coupling between the subsystems is set equal to $`1`$.

Note that an essentially more general version of the model, with the dispersion, Kerr, and filtering terms present in the equation for the passive core was also studied in detail . However, both physical arguments and numerical results presented in  clearly demonstrate that the version of the model with the very simple form of the additional equation displayed above is quite sufficient to grasp all the essential dynamical features of the dual-core optical fiber with one active and one passive cores .

Although the model (3), as well as those considered in refs.  and , has only cubic nonlinearity, it is very different from them, having another stabilization mechanism. In particular, it has been demonstrated analytically and numerically that there is no global instability in the present model (i.e., it never gives rise to blowup) .

Parameters in eqs. (3) can be chosen so that to provide for the stability of the zero solution, which lends localized pulses, if any, a chance to be stable , see below. Then, two exact zero-velocity SP solutions to eqs. (3) can be found, following the pattern of the exact solution to the cubic GL equation, in the form

$$u=u_0[\mathrm{sech}(kx)]^{(1+i\mu )}e^{i\omega t},v=v_0[\mathrm{sech}(kx)]^{(1+i\mu )}e^{i\omega t}$$
(4)
. The two exact solutions have different values of $`u_0,v_0,k,\mu `$ and $`\omega `$. Direct numerical simulations have clearly shown that one of the two exact SPs may be stable . The simulations have revealed nontrivial stability borders for this SP in the parameter space of the model (“nontrivial” implies that the stability borders are different from stability conditions for the zero background). Interactions between stable SPs were also simulated in detail in this model ; in particular, it has been found that stable bound states of two pulses exist. Three-pulse states exist too, but they are unstable against symmetry-breaking perturbations.

The objective of this work is to search for stable localized states in the model (3) beyond instability borders of the stationary pulses. In section 2, we demonstrate that effectively stable pulses and sets of pulses exist in the region where the zero background is unstable. It should be mentioned that a possibility for a pulse to survive above the background-instability threshold is known in a GL model of another type , but in that case the stabilization is provided by the fact that the pulse is moving, due an extra symmetry-breaking term added to the GL equation, and, in a system with periodic boundary conditions, it may perform a round trip quickly enough to suppress the growing local perturbations. In our model, the situation is different: immediately above the background-instability threshold, we observe stable pulses resting on top of a background standing-wave pattern with a small amplitude. With further increase of the overcriticality, the background pattern keeps a rather small amplitude but becomes chaotic. In the course of this transition, SPs remain stable, although they become randomly walking ones, rather than keeping zero velocity. While we did not try to elaborate in detail how the randomization of the background gives rise to the random walk of the pulse, the effect seems to be quite similar to the Brownian motion, resembling the Gordon-Haus effect (random walk of solitons interacting with random noise) in optical fibers . Moreover, additional pulses are spontaneously generated at larger values of the overcriticality, and spectacular correlations in the random walk of far separated pulses are observed in this case, due to their mutual repulsion.

Another finding, presented in section 3, is that one of generic modes of the instability of stationary SPs (in the case when the zero background is stable) does not destroy them, but instead transforms them into stable breathers (vibrating SPs). Stable bound states of two breathers were found too, while bound states with a larger number of breathers were not found. All these types of stable nonstationary SPs are fairly novel to the present model, and may be of considerable interest for the general analysis of models of the GL type.

## II Solitary pulses above the background-instability threshold

### A Stability of the zero background

As it was explained above, the first necessary condition for the unequivocal stability of SPs is stability of the zero solution. To investigate this issue, one should linearize eqs. (3) and look for a perturbation eigenmode $`\mathrm{exp}\left(\lambda t+iqx\right)`$, where $`q`$ is an arbitrary real perturbation wavenumber, and $`\lambda `$ is the corresponding instability growth rate. This yields a dispersion equation for $`\lambda (q)`$ :

$$\lambda ^2+(\mathrm{\Gamma }\gamma _0+\gamma _1q^2+i\gamma _2q^2)\lambda +\mathrm{\Gamma }(\gamma _0+\gamma _1q^2+i\gamma _2q^2)=0.$$
(5)
Onset of instability takes place when $`\mathrm{Re}\lambda `$ changes its sign at some (critical) value of $`q`$. As it follows from eq. (5), this critical condition gives rise to an equation for the critical value:

$$\frac{\gamma _2^2q^4\mathrm{\Gamma }(\gamma _0+\gamma _1q^2)}{(\mathrm{\Gamma }\gamma _0+\gamma _1q^2)^2}+\mathrm{\Gamma }(\gamma _0+\gamma _1q^2)+1=0.$$
(6)
In the further analysis, we will display typical numerical results obtained at fixed values

$$\sigma =1,\gamma _2=5,\gamma _1=0.9,$$
(7)
while $`\gamma _0`$ will be varied at some fixed values of $`\mathrm{\Gamma }`$. Note that the variation of the gain parameter $`\gamma _0`$ is the most physically meaningful way to scan the results, as it can be easily changed in the optical experiment by adjusting the pump power, while the other parameters are fixed for a given experimental setup.

Figure 1 displays the line $`\gamma _0(q)`$, as obtained from eq. (6) (neutral stability curve) at the fixed values (7) and $`\mathrm{\Gamma }=1.35`$. Above the line, one has $`\mathrm{Re}\lambda (q)>0`$ with $`\mathrm{Im}\lambda (q)0`$, i.e., oscillatory instability takes place. The curve has a minimum $`\left(\gamma _0\right)_{\mathrm{min}}=0.527`$ at $`q=0.54`$. That is, the zero solution is stable if $`\gamma _0\left(\gamma _0\right)_{\mathrm{min}}`$, and when the critical gain $`\left(\gamma _0\right)_{\mathrm{min}}`$ is reached, the zero background becomes unstable against finite-wavenumber perturbations. In view of the importance of the critical value of the gain, in fig. 2 we display it as a function of $`\mathrm{\Gamma }`$. The cusp at $`\mathrm{\Gamma }=1.82`$ is related to a switch of the instability mode at this point: at $`\mathrm{\Gamma }<1.82`$, the finite-wavenumber instability takes place, while at $`\mathrm{\Gamma }>1.82`$ the value $`\gamma _0=\left(\gamma _0\right)_{\mathrm{min}}`$ is attained at $`q=0`$.

### B Standing and walking solitary pulses

Simulations of the dynamics of SPs were performed by means of a pseudospectral method, assuming periodic boundary condition in $`x`$, with the period $`60`$ or $`120`$ (in some cases, the period was $`480`$, see below). As the initial configuration, we took the exact SP solutions (4), borrowing expressions for its parameters from . In those cases when SPs were found to be stable in the simulations reported in Ref. , we also saw that they were stable. In this section, we focus on a possibility of existence of robust pulses above the critical gain $`\left(\gamma _0\right)_{\mathrm{min}}`$, when it appears that SPs cannot be stable.

Figure 3 displays SP produced by the simulations at $`\gamma _0=0.54`$ and $`\mathrm{\Gamma }=1.35`$. As it is seen from fig. 1, this point is located slightly above the neutral stability curve, so that the corresponding overcriticality parameter is $`ϵ\left[\gamma _0\left(\gamma _0\right)_{\mathrm{min}}\right]/\left(\gamma _0\right)_{\mathrm{min}}0.025`$. As is seen in fig. 3a, the weak background instability generates a small-amplitude standing-wave pattern, the stable localized pulse existing on top of it. To further illustrate the properties of this solution, in fig. 3b we specially display the time dependence of $`\mathrm{Re}u`$ at two points: in the center of the localized pulse, and in the background. Both dependences are periodic, but their frequencies are very different (the background oscillations are much slower), i.e., SP and the background are not synchronized.

The background oscillations are regular at $`\gamma _0=0.54`$; however, they become chaotic for $`\gamma _0>0.55`$, at the same value $`\mathrm{\Gamma }=1.35`$. We use the term ”chaotic”, since the time evolution of the background is irregular and has no periodicity. Figure 4a displays a localized pulse existing on top of the chaotic background at $`\gamma _0=0.61`$, the corresponding overcriticality being $`ϵ=0.157`$. Figure 4(b) displays the time evolution of $`|u(x,t)|`$. Interaction of the localized pulse with the chaotic background apparently gives rise to a random walk of the pulse; however, the pulse is not destroyed. The temporal evolution of the SP’s central coordinate $`X_p`$ is shown in fig. 5a. To analyze the character of the random walk of the localized pulse, we display in fig. 5(b) the time-averaged squared displacement $`\left(\mathrm{\Delta }X_p\right)^2`$, where $`\mathrm{\Delta }X_p\mathrm{\Delta }X_p(t+\mathrm{\Delta }t)\mathrm{\Delta }X_p(t)`$, vs. the temporal interval $`\mathrm{\Delta }t`$. As it is obvious from fig. 5b, the random walk approximately obeys the diffusive law, $`\left(\mathrm{\Delta }X_p\right)^2\mathrm{\Delta }t`$.

Spontaneous formation of new pulses takes place for larger $`\gamma _0`$. Note that spontaneous generation of effectively stable pulses by the unstable background was observed in simulations of another GL model in Ref. . It is not quite clear whether there is a definite border for the spontaneous generation of SPs from the chaotic background. The pattern shown in fig. 6 was generated from a nearly zero initial state (transient evolution is not shown) at $`\gamma _0=0.64`$, corresponding to $`ϵ=0.214`$. Three localized pulses are spontaneously generated in this case. The background is chaotic, and the three SPs are not identical. It is clearly observed that the pulses repel each other, so that they keep nearly equal separations between themselves (recall we simulated the system (3) with periodic boundary conditions). Figure 6b displays the time evolution of $`|u(x,t)|`$. A noteworthy feature of the established dynamical state is that the random walk of all the three pulses is highly synchronized, so that the pulses keep nearly equal separations between themselves in the course of the random walk. The number of the spontaneously generated pulses grows with the system’s size. In fig. 7 we show the established state at the same values of the parameters as in fig. 6, but with the size of the system $`480`$ instead of $`120`$. The separations between SPs remain nearly equal, while the synchronization of the random walk of SPs is weaker in a large system with the chaotic background. The synchronization does take place between adjacent pulses, but the full global synchronization does not occur, as is seen in fig. 7.

The number of the pulses in an established state increases with $`\gamma _0`$ but it is not a uniquely defined function of the system’s parameters and size. More detailed simulations (not shown here) demonstrate that extra pulses can be added, and some may be removed, i.e., their exact number depends on the way the state was made.

Finally, at large values of the overcriticality, many pulses are created, and the system goes over into a more irregular “turbulent” state (in which there is no spatially regular structure of localized pulses). The transition occurs close to at $`\gamma _0=1.0`$; however, this depends on the number of pulses. Examples are shown in fig. 8 for $`\gamma _0=1.05`$ (a) and $`\gamma _0=2`$ (b) (again, with $`\mathrm{\Gamma }=1.35`$). Note that these states have been obtained starting from nearly-zero initial conditions. Generation of pulse sets and the turbulent behavior are seen in these pictures. The time evolution is very fast, and the final number of pulses is larger for $`\gamma _0=2.`$

## III Breathers and their bound states

The exact pulse solutions to eqs. (3) found in  are stable only in a part of the parameter region where they exist . Besides the loss of the stability of the zero background, other, more nontrivial destabilization mechanisms occur as well . In particular, it seems very plausible that the stationary pulse can be destabilized through a Hopf bifurcation.

In fig. 9, we show a stable breather replacing the stationary pulse at $`\gamma _0=0.18`$ and $`\mathrm{\Gamma }=0.2`$. Figure 9a displays the time evolution of the profile $`|u(x,t)|`$, and fig. 9(b) displays the time evolution of $`|u(x_p,t)|`$ at the peak position. Limit-cycle intrinsic oscillations of the pulse are clearly seen. Breathing localized pulses were earlier found in a reaction-diffusion model  and in the quintic GL equation  (note that the breathing solution in fig. 9a has a form similar to the well-known two-soliton solution of the nonlinear Schrödinger equation, $`u=4\eta [\mathrm{cosh}(3kx)+3\mathrm{exp}(4i\eta ^2t)\mathrm{cosh}(kx)]/[\mathrm{cosh}(4kx)+4\mathrm{cosh}(2kx)+3\mathrm{cos}(4\eta ^2t)]`$). To present a whole family of the breathers, in fig. 10a we show the amplitude of the oscillations of $`|u(x_p,t)|`$ at the breather’s central point as a function of $`\gamma _0`$ for $`\mathrm{\Gamma }=0.2`$. At the value $`\gamma _00.156`$, at which the oscillation amplitude vanishes, the breather changes into the stationary SP. Figure 10b displays the square amplitude $`A^2`$ as a function of $`\gamma _0`$ near the critical point. The linear dependence of $`A^2`$ on $`\gamma _0`$ implies $`A\sqrt{\gamma _0\gamma _{0c}}`$ near the critical point $`\gamma _{0c}`$, that is, the bifurcation is the supercritical Hopf bifurcation.

If the gain parameter $`\gamma _0`$ becomes too small, no nontrivial solution can exist in the system. Thus, there is a minimum value of $`\gamma _0`$ confining the existence region of the localized pulses. Collecting results from many runs of the simulations performed at different values of $`\gamma _0`$ and $`\mathrm{\Gamma }`$, while the other parameters took the fixed values (7), in fig. 11 we have drawn the phase diagrams of the localized solutions existing below the threshold of the zero-background’s instability. The zero background is stable beneath the line marked by squares, the breathing pulses change into stationary ones at the line marked by crosses, and the stationary pulses are changed by the zero solution at the line marked by diamonds. In accord with what was said above, the border between the stationary pulses and breathers (crosses in fig. 11) is sharp, i.e., there is no overlapping between them. That is a natural consequence of the supercritical character of the Hopf bifurcation (see above), which makes coexistence of stable stationary and oscillating solutions very implausible.

The Hopf-bifurcation line collides with the zero-background’s instability line near $`\mathrm{\Gamma }=0.3`$ and $`\gamma _0=0.27`$. For $`\mathrm{\Gamma }>0.3`$, stable breathing pulses do not exist. We stress that, unlike the effectively stable nonbreathing pulses considered in section 2, we have not observed stable breathing pulses on top of the small-amplitude background above the zero-background’s instability threshold in the region $`\mathrm{\Gamma }<0.3`$ .

Bound states of stable localized pulses also play an important role in the dynamics of GL-type models (see, e.g.,  and ). Our simulations have revealed a possibility of forming a stable bound state of two breathers. Figure 12 displays a bound state found at $`\gamma _0=0.18`$ and $`\mathrm{\Gamma }=0.2`$. A noteworthy feature of this bound state is that the phase shift between internal vibrations in the two breathers is $`\pi /2`$, i.e., one breather has a maximum amplitude when the amplitude of its mate is minimal, and vice versa. Attempts to generate a bound state of three breathers have failed. In this connection, it is relevant to mention that fully stable bound states of two stationary pulses can be readily generated in the same model, while complexes consisting of three pulses exist, but they are destroyed by perturbations breaking their symmetry .

## IV Conclusion

In this work, we have considered a model based on the cubic complex Ginzburg-Landau equation which is linearly coupled to an additional linear equation, which was introduced in  as a simplification of a more general model originally proposed in . The model describes a dual-core nonlinear optical fiber with one active and one passive core. By means of systematic simulations, we have found new types of stable localized excitations in this model, existing in addition to the earlier found stationary pulses. The newly found excitations include pulses existing on top of a small-amplitude background (which is regular or chaotic) above the threshold of instability of the zero solution, and breathers into which stationary pulses are transformed by a Hopf bifurcation below the zero-solution instability threshold. A sharp border between the stable stationary pulses and breathers was identified. Stable bound states of two breathers with a phase shift $`\pi /2`$ between their internal vibrations have been found too. Above the threshold, the pulses are, respectively, standing if the background oscillations are regular, and randomly walking if the background is chaotic. With the increase of the system’s size, additional randomly walking pulses are spontaneously generated. The random walk in a multi-pulse state is synchronized (but not completely) due to repulsion between the pulses. At a large overcriticality, the multi-pulse state goes over into a turbulent one. Lastly, breathers do not exist on top of the small-amplitude background.

One of the authors (B.A.M.) appreciates financial support from the Japan Society for Promotion of Science, and hospitality of the Department of Applied Physics at the Faculty of Engineering, Kyushu University (Fukuoka, Japan).

## FIGURE CAPTIONS

Fig. 1. The neutral stability curve $`\gamma _0(q)`$ at $`\sigma =1`$, $`\gamma _1=0.9`$, $`\gamma _2=5`$, and $`\mathrm{\Gamma }=1.35`$.

Fig. 2. The minimum value of the gain giving rise to instability of the zero background vs. $`\mathrm{\Gamma }`$ at $`\sigma =1`$, $`\gamma _1=0.9`$, and $`\gamma _2=5`$.

Fig. 3. A stable solitary pulse existing on top of the small-amplitude background at $`\gamma _0=0.54`$ and $`\mathrm{\Gamma }=1.35`$: the wave field $`|u(x)|`$ (a), and the time dependence of $`\mathrm{Re}u`$ at two fixed spatial points, one at the center of the pulse and one in the background (b).

Fig. 4. Snapshot of the wave field $`|u(x,t)|`$ in the localized pulse existing on top of the chaotic background in the case $`\gamma _0=0.61`$ and $`\mathrm{\Gamma }=1.35`$ (a), and the time evolution of $`|u(x,t)|`$ (b).

Fig. 5. Details of the random walk of the coordinate $`X_p(t)`$ (a), and diffusive growth of the mean-square displacement $`\left(X_p(t+\mathrm{\Delta }t)X_p(t)\right)^2`$ (b) of the solitary pulse from fig. 4.

Fig. 6. Instantaneous snapshot (a) and the time evolution (b) of the wave field $`|u|`$ at $`\gamma _0=0.61`$ and $`\mathrm{\Gamma }=1.35`$. Synchronization of the random walk of three spontaneously formed pulses is seen.

Fig. 7. The same as in fig. 6(b), but for the system four times as long. As an initial condition, we use four copies of the snapshot pattern of fig. 6(a), that is, $`u(x,0)`$ for $`0<x<120`$ takes the randomly perturbed value of fig. 6(a), and $`u(x,0)`$ for $`120<x<240`$ also takes the value of fig. 6(a) perturbed randomly, and so on. So initially there are 12 pulses. It is clearly seen that the synchronization of the random walk of different pulses is imperfect.

Fig. 8. A “turbulent” configuration of $`|u(x,t)|`$ at $`\gamma =1.05`$ (a) and $`\gamma =2`$ (b) for $`\mathrm{\Gamma }=1.35`$.

Fig. 9. A stable breather observed at $`\gamma _0=0.18`$ and $`\mathrm{\Gamma }=0.2`$: the field $`|u(x,t)|`$ (a), and its value at the central point of the breather vs. time (b).

Fig. 10. The amplitude $`A`$ of the oscillations of the breather’s field at its central point (see fig. 9b) vs. the gain parameter $`\gamma _0`$ (a), and $`A^2`$ vs. $`\gamma _0`$ near the critical point (b). The value $`\gamma _00.156`$, at which the oscillation amplitude vanishes, is a point of the transition (Hopf bifurcation) between the stationary pulse and breather.

Fig. 11. The phase diagram in the parametric plane ($`\mathrm{\Gamma }`$,$`\gamma _0`$), the other parameters taking the fixed values (7). The zero background is stable beneath the line marked by squares; localized pulses disappear (decaying to zero) beneath the line marked by diamonds; crosses mark the border between the stable stationary solitary pulses (above the border) and stable breathers (below the border).

Fig. 12. A stable bound state of two breathers in the same case as in fig. 9. Two pulses are set up apart from each other as an initial condition.


text:
# Filling up the Two-Body Partial Decay Width of P-Wave Charmonia with the Color Octet Contribution

Quarkonia are special hadronic systems in which the constituent quarks because of their large mass move very slowly therefore they are open to non-relativistic treatment. So quarkonium production or decay, once supplemented with non-perturbative wavefunctions or matrix elements can be calculated using the small velocity as an expansion parameter . Recently some exclusive $`\chi _J`$ decay channels have been re-measured at the Beijing Spectrometer by the BES collaboration . An attempt to compare with existing calculations of the channels $`\chi \pi \pi `$ and $`\chi p\overline{p}`$ revealed that these were out-of-date. Various parameters have changed since the last serious attempts at calculating the partial widths and the understanding on some of the hadronic wavefunctions has also improved. It is therefore necessary to re-examine the decay calculations.

The decay processes of interest are best calculated using the perturbative QCD based scheme of Brodsky and Lepage (BL) . In this scheme, decay amplitudes are given by a convolution of non-perturbative hadronic distribution amplitudes and the associated perturbative hard part through which the hard momentum flows. Using this scheme with the latest parameters such as $`\mathrm{\Lambda }_{\mathrm{QCD}}`$, decay constants and distribution amplitudes of pion and nucleon, it was found that in both channels the theoretical width of the $`c\overline{c}`$ system fell far short of the experimental results. The widths for the decay into $`\pi \pi `$ are at most 23% for $`\chi _0`$ and 28% for $`\chi _2`$ of the experimental data. Those for the decay into $`p\overline{p}`$ are even worse at only 8.3% for $`\chi _1`$ and 10.4% for $`\chi _2`$. The next logical step to find out why there are such large discrepancies is to try the calculations again using the improved BL scheme of Sterman et al  hoping that the advantages in the improved scheme could compensate for whatever in the original scheme that might have caused the deficiencies . Unfortunately this is not to be. The new scheme actually yielded widths of the same magnitude as obtained from the old. Varying all parameters within their uncertainties to enlarge the widths did not help either. The results were still at least a factor of two or more below. The widths at the size of the experimental data remained unattainable. So the small theoretical results did not arise because of some artifact in the scheme or the poor choice of parameters used.

While searching for an explanation, it became obvious that color octet might be the answer as this was shown to be a necessary component in inclusive P-wave charmonium decay to ensure infrared finite answers . Strangely color octet has never been considered within the context of exclusive reactions. At first sight it cannot be the answer because higher Fock state contributions are suppressed by large momentum flow within the BL scheme. On the other hand, a contradiction would arise if none of the exclusive $`\chi _J`$ decay channels, many-body decays included, requires the octet contribution. This is because the octet state must be there to ensure infrared finite inclusive width, which is just the sum of partial widths of all exclusive channels.

An examination of the charmonium wavefunctions show that in momentum space the P-wave wavefunction is down by $`k/M`$ in comparison to the S-wave due to angular momentum, where $`k`$ is the internal momentum of the heavy quarks within the charmonium and $`M`$ is the mass of the latter. Extending this consideration to the decay amplitude of the $`\pi \pi `$ and $`p\overline{p}`$ channels, it can be shown that by examining the large mass $`M`$ dependence within the BL scheme that the singlet and octet amplitudes for the decay into $`\pi \pi `$ go like $`1/M^3`$ and those for $`p\overline{p}`$ channel both also have the same dependency of $`1/M^5`$ . So despite the fact that the octet is a higher state, it is not suppressed relative to the valence state. The suppression of the higher state in this case can no longer be thought of as one when the P-wave valence state is brought down to the same level by the orbital angular momentum. Once this theoretical hurdle has been overcome, it remains to work out new ingredients to be used in the calculations. One is the yet unknown wavefunction of the color octet state and the other is more about the practicality of what to do with the constituent gluon. Neither of these was ever needed or encountered in the calculations of exclusive processes. Nevertheless the former can be modeled and fit to data, and the latter argued to yield a dominant contribution if the constituent gluon line ends on the hard part in the Feynman graphs . With these details sorted out, it was found that indeed the color octet contributions were capable of bringing the theoretical widths up to the size of experimental data . Thus confirming the theoretical expectations and resolving the problem of too small theoretical widths.


text:
# Reconstructing the inflationary power spectrum from CMBR data

## I introduction

Fluctuations in the Cosmic Microwave Background Radiation (CMBR) were detected for the first time by the COBE satellite in 1992 . Subsequently it was realised that precision measurements of the fluctuation spectrum can be used to infer almost all of the fundamental cosmological parameters . The method most commonly used for determining parameters from the data is to maximise the likelihood function, $``$, over the space of model parameters to be determined, $`\theta =\{\mathrm{\Omega },\mathrm{\Omega }_m,\mathrm{\Omega }_b,H_0,\tau ,\mathrm{}\}`$. A large number of papers have dealt with with this issue in great detail. However, when calculating the theoretical CMBR power spectrum needed for the likelihood function, a necessary input is the initial power spectrum of fluctuations, usually assumed to have been produced during the inflationary epoch.

There exists a plethora of different inflationary models, each having a specific prediction for the produced fluctuation spectrum. In the simplest slow-roll models, the spectrum is described by a simple featureless power-law $`P(k)k^{n_s1}`$, where $`n_s1`$ . However, the power spectrum could easily behave differently, depending both on the specific form of the inflaton potential, $`V(\varphi )`$, and on the presence of different physical phenomena during inflation. A plethora of different models predicting such behaviour exists . Examples of such phenomena are the resonant production of particles during inflation, proposed by Chung et al. , and the multiple inflation model discussed by Adams, Ross and Sarkar .

In almost all existing likelihood analyses the initial power spectrum is assumed to have the power-law form $`P(k)k^{n_s1}`$, where $`n_s`$ is a constant . This type of analysis restricts the parameter estimation in such a way that no non-power law features can be detected. Since there are so many different inflationary models, each with unique predictions, it is highly desirable to have a more model-independent way of estimating the initial power spectrum.

Such a possibility has been discussed by Souradeep et al.  and by Wang, Spergel and Strauss . In both treatments, the power spectrum was described as a featureless spectrum, but binned in $`k`$-space. Using a Fisher matrix analysis, it was shown by Wang, Spergel and Strauss that it is possible to determine the power in each bin to reasonable precision, so that some features could be detectable. In the next section we discuss the prospects for detecting features in more detail, also using the Fisher matrix technique. Particularly, it is shown that there is a trade-off between the precision to which the power in each bin can be measured and the width of each bin.

In Section III we go on to discuss in detail how such a likelihood analysis with a binned power spectrum is performed in practise. This analysis highlights some of the difficulties which can be expected. In general it is quite difficult to reconstruct a power spectrum with smooth features using a set of discontinuous bin amplitudes. If the binning is too coarse, the individual bins being comparable in size or broader than the power spectrum features, the spectrum reconstruction can easily give misleading results. For a reasonable reconstruction it is necessary to choose a binning which is significantly finer than the features to be detected. However, maximising the likelihood function is a highly non-linear optimization problem. If the number of free parameters in the fit, $`N=N_\theta +N_{\mathrm{bins}}`$, is too large the computation time to retrieve the maximum likelihood becomes very long. Nevertheless, it is demonstrated that it is computationally feasible to use at least 20 bins in $`k`$-space.

Finally, Section IV contains a discussion of the results, with emphasis on how the techniques can be used on results from high precision satellite experiments, like MAP and PLANCK .

## II Fisher Matrix Analysis

CMBR temperature fluctuations are usually expressed in terms of spherical harmonics as

$$\frac{\mathrm{\Delta }T}{T}(\theta ,\varphi )=\underset{lm}{}a_{lm}Y_{lm}(\theta ,\varphi ).$$
(1)
From these $`a_{lm}`$ coefficients it is possible to construct the power spectrum as

$$C_l=|a_{lm}|^2,$$
(2)
where the average in principle is an ensemble average. Such an ensemble average can obviously not be performed since we have access to only one realisation of the underlying distribution. However, in a universe which is isotropic the ensemble average can be replaced by an average over $`m`$-values for a given $`l`$ .

It is possible to estimate the precision with which the cosmological model parameters can be extracted from a given hypothetical data set. The starting point for any parameter extraction is the vector of data points, $`x`$. This can be in the form of the raw data, or in compressed form, either as the $`a_{lm}`$ coefficients or the power spectrum, $`C_l`$. Each data point has contributions from both signal and noise, $`x=x_{\mathrm{CMBR}}+x_{\mathrm{noise}}`$. If both signal and noise are Gaussian distributed it is possible to build a likelihood function from the measured data which has the following form

$$(\mathrm{\Theta })\mathrm{exp}\left(\frac{1}{2}x^{}[C(\mathrm{\Theta })^1]x\right),$$
(3)
where $`\mathrm{\Theta }=(\mathrm{\Omega },\mathrm{\Omega }_b,H_0,n,\tau ,\mathrm{})`$ is a vector describing the given point in model parameter space and $`C(\mathrm{\Theta })=xx^T`$ is the data covariance matrix. In the following we shall always work with data in the form of a set of power spectrum coefficients, $`C_l`$.

If the data points are uncorrelated so that the data covariance matrix is diagonal, the likelihood function can be reduced to $`e^{\chi ^2/2}`$, where

$$\chi ^2=\underset{l=2}{\overset{N_{\mathrm{max}}}{}}\frac{(C_{l,\mathrm{obs}}C_{l,\mathrm{theory}})^2}{\sigma (C_l)^2},$$
(4)
is a $`\chi ^2`$-statistics and $`N_{\mathrm{max}}`$ is the number of power spectrum data points .

The maximum likelihood is an unbiased estimator, which means that

$$\mathrm{\Theta }=\mathrm{\Theta }_0.$$
(5)
Here $`\mathrm{\Theta }_0`$ indicates the true parameter vector of the underlying cosmological model and $`\mathrm{\Theta }`$ is the average estimate of parameters from maximising the likelihood function.

The likelihood function should thus peak at $`\mathrm{\Theta }\mathrm{\Theta }_0`$, and we can expand it to second order around this value. The first order derivatives are zero, and the expression is thus

$$\chi ^2=\chi _{\mathrm{min}}^2+\underset{i,j}{}(\theta _i\theta )\left(\underset{l=2}{\overset{N_{\mathrm{max}}}{}}\frac{1}{\sigma (C_l)^2}\left[\frac{C_l}{\theta _i}\frac{C_l}{\theta _j}(C_{l,\mathrm{obs}}C_l)\frac{^2C_l}{\theta _i\theta _j}\right]\right)(\theta _j\theta ),$$
(6)
where $`i,j`$ indicate elements in the parameter vector $`\mathrm{\Theta }`$. The second term in the second derivative can be expected to be very small because $`(C_{l,\mathrm{obs}}C_l)`$ is in essence just a random measurement error which should average out. The remaining term is usually referred to as the Fisher information matrix

$$F_{ij}=\frac{^2\chi ^2}{\theta _i\theta _j}=\underset{l=2}{\overset{N_{\mathrm{max}}}{}}\frac{1}{\sigma (C_l)^2}\frac{C_l}{\theta _i}\frac{C_l}{\theta _j}.$$
(7)
The Fisher matrix is closely related to the precision with which the parameters, $`\theta _i`$, can be determined. If all free parameters are to be determined from the data alone without any priors then it follows from the Cramer-Rao inequality  that

$$\sigma (\theta _i)=(F^1)_{ii}$$
(8)
for an optimal unbiased estimator, such as the maximum likelihood .

In general, the standard error on the $`C_l`$ coefficients can be written as

$$\sigma (C_l)=\left[\frac{2}{(2l+1)f_{\mathrm{sky}}}\right]^{1/2}(C_l+\mathrm{\Delta }_{\mathrm{exp}}).$$
(9)
Here, $`f_{\mathrm{sky}}`$ is the sky coverage and $`\mathrm{\Delta }_{\mathrm{exp}}`$ is a (Gaussian) experimental error. Notice that even for $`\mathrm{\Delta }_{\mathrm{exp}}=0`$, $`\sigma (C_l)0`$. This derives from the fact that the ensemble average in Eq. (2) has been replaced by an average over $`m`$-values, and is usually referred to as “Cosmic Variance” .

The CMBR is also predicted to be polarized, and this polarization power spectrum can in principle be measured. For scalar perturbations there are only three independent quantities: The temperature, $`C_{l,T}`$, the $`E`$-field polarization, $`C_{l,E}`$, and the temperature-polarization cross-correlation, $`C_{l,C}`$. Then the Fisher matrix of Eq. (7) is instead

$$F_{ij}=\underset{l=2}{\overset{N_{\mathrm{max}}}{}}\underset{X,Y}{}\frac{C_{l,X}}{\theta _i}\mathrm{Cov}^1(C_{l,X},C_{l,Y})\frac{C_{l,X}}{\theta _j}.$$
(10)
In this case the covariance matrix, $`\mathrm{Cov}(C_{l,X},C_{l,Y})`$, is a symmetric $`3\times 3`$ matrix with the following elements

$`\mathrm{Cov}(C_{l,T},C_{l,T})`$ $`=`$ $`\left[{\displaystyle \frac{2}{(2l+1)f_{\mathrm{sky}}}}\right](C_{l,T}+\mathrm{\Delta }_{\mathrm{exp},1})^2`$ (11)
$`\mathrm{Cov}(C_{l,E},C_{l,E})`$ $`=`$ $`\left[{\displaystyle \frac{2}{(2l+1)f_{\mathrm{sky}}}}\right](C_{l,E}+\mathrm{\Delta }_{\mathrm{exp},2})^2`$ (12)
$`\mathrm{Cov}(C_{l,C},C_{l,C})`$ $`=`$ $`\left[{\displaystyle \frac{2}{(2l+1)f_{\mathrm{sky}}}}\right](C_{l,C}^2+(C_{l,T}+\mathrm{\Delta }_{\mathrm{exp},1})(C_{l,E}+\mathrm{\Delta }_{\mathrm{exp},2}))`$ (13)
$`\mathrm{Cov}(C_{l,T},C_{l,E})`$ $`=`$ $`\left[{\displaystyle \frac{2}{(2l+1)f_{\mathrm{sky}}}}\right]C_{l,C}^2`$ (14)
$`\mathrm{Cov}(C_{l,T},C_{l,C})`$ $`=`$ $`\left[{\displaystyle \frac{2}{(2l+1)f_{\mathrm{sky}}}}\right]C_{l,C}(C_{l,T}+\mathrm{\Delta }_{\mathrm{exp},1})`$ (15)
$`\mathrm{Cov}(C_{l,E},C_{l,C})`$ $`=`$ $`\left[{\displaystyle \frac{2}{(2l+1)f_{\mathrm{sky}}}}\right]C_{l,C}(C_{l,E}+\mathrm{\Delta }_{\mathrm{exp},2}),`$ (16)
where $`\mathrm{\Delta }_{\mathrm{exp},i}`$ are again experimental experimental errors related to pixel noise and beamwidth .

For the purposes of the present paper we assume that $`f_{\mathrm{sky}}=1`$, $`\mathrm{\Delta }_{\mathrm{exp},i}=0`$ and $`l_{\mathrm{max}}=1500`$, corresponding to a full-sky survey up to $`l_{\mathrm{max}}=1500`$, limited only by cosmic variance. If the temperature power spectrum alone is considered this is not too far from what can be expected with MAP, whereas PLANCK will likely measure both temperature and polarization to this accuracy.

As the free cosmological parameters we use the matter density, $`\mathrm{\Omega }_m`$, the cosmological constant, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, the baryon density, $`\mathrm{\Omega }_b`$, the Hubble parameter, $`H_0`$, the optical depth to reionization, $`\tau `$, and the overall normalization, $`Q`$. Instead of using the spectral index of a featureless power-law spectrum, $`n_s`$, as the last free parameter, we bin the spectral index in $`k`$-space. In practise we use $`N`$ bins of equal size in $`\mathrm{log}(k)`$, from $`k=1.27\times 10^5\mathrm{Mpc}^1`$ to $`k=0.25\mathrm{Mpc}^1`$ (this covers the entire range of $`k`$-space which is visible in the CMBR). In each bin the power spectrum is assumed to follow a power-law, $`P(k)_ik^{n_s(i)1}`$. Thus, the parameter vector is

$$\mathrm{\Theta }=\{\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_b,H_0,\tau ,Q,n_s(1),\mathrm{},n_s(N)\}$$
(17)
As the reference model around which to calculate the Fisher matrix, we take the standard CDM model with, $`\mathrm{\Omega }=\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_b=0.05`$, $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and $`\tau =0`$. The reference model has a power-law initial spectrum with $`n_s(i)=1`$ for all $`i`$.

In Fig. 1 the precision with which $`n_s(i)`$ can be measured is shown for different numbers of bins, for the case where only the temperature power spectrum can be measured. Not surprisingly, the precision with which the effective spectral index in each bin can be measured depends strongly on the number of bins. For only one bin, $`\mathrm{\Delta }n_s/n_s=8.21\times 10^3`$, whereas for 20 bins, the smallest $`\mathrm{\Delta }n_s/n_s`$ is $`3.49\times 10^2`$. Thus, there is a trade-off between the resolution in $`k`$-space and the resolution in $`n_s`$. Detecting a low-amplitude narrow feature is not possible, whereas either a broad low-amplitude or a narrow high-amplitude feature should be detectable. In agreement with Wang, Spergel and Strauss  we find that the CMBR data is most sensitive in the range $`k10^310^1\mathrm{Mpc}^1`$. If polarization can also be measured, the precision with which the $`n_s(i)`$ can be measured increases. Fig. 2 shows the same as Fig. 1, but with the inclusion of polarization. Indeed, the precision on the individual bin amplitudes increases by a large factor. For only one bin, $`\mathrm{\Delta }n_s/n_s=2.90\times 10^3`$, whereas for 20 bins, the smallest $`\mathrm{\Delta }n_s/n_s`$ is now $`9.49\times 10^3`$.

Next, a fundamental question is how much the ability to determine the other fundamental parameters depends on the number of bins in $`k`$-space. Figs. 3 and 4 show the precision in measuring the other cosmological parameters as a function of $`N`$, the number of bins. Fig. 2 is for the case where only the temperature anisotropy can be measured. For all the cosmological parameters, except $`\tau `$, there little degeneracy between these parameters and the power spectrum indices, meaning that there is little loss of ability to pin down these parameters. From $`N=10`$, the resolution begins to decrease, although quite slowly, and at $`N=30`$ it is for instance down by about 40% for $`\mathrm{\Omega }_b`$. Clearly, the spectral indices, $`n_s(i)`$, are not very degenerate with the other cosmological parameters. The reason is that the discontinuous binning in the power spectrum is very hard to mimic by continuous changes in other parameters. This finding is in agreement with Ref. , where such a binning was also used. However, there is a marked difference compared with the findings of Ref. . In this work, smoothness of the power spectrum was enforced by using smooth shape functions peaked at a given point in $`k`$-space instead of discontinuous binning. These authors find that the uncertainty in the other parameters increases very rapidly with $`N`$, rising by more than a factor of 10 at $`N=30`$ . However, using smooth functions to describe the power spectrum could indeed be expected to be able to mimic changes in other parameters to a much higher degree than our approach. So it is not surprising that a much higher degree of degeneracy is found. On the other hand, as will be seen in the next section, using a discontinuous binning has the very clear disadvantage that it can be difficult to achieve reasonable likelihood fits if the number of bins is too small.

In line with this discussion, it was shown by Kinney  that smooth power spectrum features can mimic variations in the cosmological parameters almost exactly (to within cosmic variance). In this work 75 bins in $`k`$-space were used and the power spectrum smoothened by cubic spline interpolation, leading to a very large degree of degeneracy between the $`k`$-space power spectrum and the other cosmological parameters.

The one exception to the above is the optical depth to reionization, $`\tau `$. Here there is a very large degree of degeneracy, which increases rapidly with the number of bins. This effect was also described by Wang, Spergel and Strauss . It happens because the effect of reionization is to suppress the amplitude at small scales by a factor $`e^{2\tau }`$. This can be quite easily mimicked by suppressing the $`k`$-space power spectrum at small scales.

However, Fig. 4 shows the case where polarization can also be measured. In this case, the degeneracy between $`\tau `$ and the $`k`$-space spectrum is broken, and the ability to determine the value of $`\tau `$ is increased by a large factor. For the other cosmological parameters, the precision is also increased, but not as dramatically.

## III Numerical Spectrum Reconstruction

Having done this initial estimate of how precisely the initial power spectrum can be estimated, it is important to see how such a power spectrum reconstruction works in practise. In order to investigate this, we have produced 10 synthetic spectra based on a underlying theoretical model. These spectra are calculated assuming Gaussian errors given by cosmic variance only. In order to simulate a possible feature in the power spectrum we introduce a Gaussian “bump”, so the spectral index has the form

$$n_s=1A\mathrm{exp}\left(\frac{(\mathrm{log}(k_{})\mathrm{log}(k_{0,}))^2}{\alpha }\right),$$
(18)
where $`k_{}=k/1\mathrm{Mpc}^1`$. This bump is characterised by three parameters: $`A`$, the amplitude, $`k_{0,}`$, the position in $`k`$-space, and $`\alpha `$, the width. For the underlying cosmological model we have chosen the same standard CDM model as in the last section, so the model parameters are: $`\mathrm{\Omega }=\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_b=0.05`$, $`H_0=50\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`A=0.7`$, $`\alpha =0.1`$ and $`\mathrm{log}k_{0,}=2`$. Thus, the bump has been placed closed to the region where the CMBR data is most sensitive. Note that in this section we use only information related to the temperature power spectrum, not the polarisation spectrum.

On each of these synthetic spectra we try reconstructing the power spectrum using the logarithmic binning method introduced in the previous section. The likelihood optimization algorithm is based on the simulated annealing principle described in Ref. . This method has the advantage that it is very fast for highly non-linear optimization over many-dimensional parameter spaces, exactly the problem that is faced here. In Fig. 3 we show the values extracted by the parameter reconstruction for three different numbers of bins, $`N=5,10`$ and 20.

In the $`N=5`$ case the individual bins are broader than the bump to be reconstructed. This shows up in the fact that no good fit is achieved ($`\chi _{\mathrm{min}}^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=12.7`$), and the spectrum reconstruction is quite poor. The reconstructed spectrum does show some evidence of a lower spectral index around the position of the bump, but the magnitude is not correct. For $`N=10`$ the reconstruction is already much better and reproduces the actual slope of the underlying spectrum. However, $`\chi _{\mathrm{min}}^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=8.06`$, so the fit is not very good in this case either. In the last case where $`N=20`$ the reconstructed spectrum again follows the underlying one reasonably well. However, at $`N=20`$ the uncertainty on the individual bin amplitudes is already significantly larger than for $`N=10`$. By going to an even higher number of individual bins, the spectrum bump will be smeared out beyond recognition. The $`\chi ^2`$-fit is significantly better than for $`N=10`$ ($`\chi _{\mathrm{min}}^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=2.84`$), although far from what should be expected from a “good” fit ($`\chi _{\mathrm{min}}^2/\mathrm{d}.\mathrm{o}.\mathrm{f}1`$).

The reason for these very high $`\chi ^2`$ values is that it is quite hard to mimic a smooth feature in the power spectrum with a sequence of discontinuous bins. Thus, even though it is indeed possible to map out the shape of the power spectrum, it is not possible to achieve a good fit in terms of the likelihood without letting $`N\mathrm{d}.\mathrm{o}.\mathrm{f}`$. However, going much beyond $`N=20`$ increases the uncertainty in the individual bin amplitudes and doing the likelihood maximization is already quite demanding numerically at $`N=20`$.

Finally, with a large number of bins (10 or 20) the reconstruction at scales smaller than the bump (large $`k`$) is consistently predicting a too low spectral index. This effect is worrying and to investigate whether it is a generic feature of any discontinuous binning method, or specific to our choice, where each bin has a “tilt”, $`P(k)_ik^{n_s(i)1}`$, we also do the reconstruction using the same method as Wang, Spergel and Strauss . Here, the power in each bin is assumed to have constant amplitude over the entire bin so that $`P(k)_i=A_i`$. We show the reconstruction for the case of $`N=20`$ in Fig. 6. For this case, the reconstruction of the power spectrum bump is actually better than for the “tilt” spectrum reconstruction. There is no trace of the spuriously low spectral index found at high $`k`$. However, the $`\chi ^2`$ is substantially worse for this method, with $`\chi _{\mathrm{min}}^2/\mathrm{d}.\mathrm{o}.\mathrm{f}=3.45`$. From this, a robust conclusion is that for discontinuous binning the obtainable fits are quite poor in terms of $`\chi ^2`$, regardless of the binning method. However, it is possible to detect features in the spectrum by this method. At least the tilted binning method can show spurious effects, so one should be careful about making strong conclusions without going through several different methods of spectrum reconstruction.

Just for comparison we have also done a spectrum reconstruction where the free parameters are those describing a smooth scale invariant power-law spectrum ($`n_s=1`$) with a Gaussian bump: $`A,k_{0,}`$ and $`\alpha `$, instead of using $`n_s(i)`$. The values extracted from likelihood maximisation on the 10 synthetic spectra are given in Table I. Clearly, in this case, the true underlying spectrum is recovered to very good accuracy because the same functional form is used in the fit as in the underlying power spectrum. Also the $`\chi ^2`$ values are completely consistent with what is expected from a good fit.

So using smooth fitting functions could have the advantage of being better able to fit smooth features in the underlying initial power spectrum. However, the discontinuous binning method has the advantage of being very model independent. So a possible way to proceed would be to initially try the discontinuous binning method, and only for refinements go to more complicated fitting functions.

## IV Discussion

We have discussed the possibility for measuring the shape of the primordial power spectrum by using precision measurements of CMBR fluctuations. The method employed was to parametrise the spectrum as a power-law, binned in $`k`$-space. An analysis of the obtainable precision in measuring the power spectrum, using a Fisher matrix analysis, showed that it should be possible to detect features that are sufficiently broad by using CMBR observations.

However, there are some practical problems involved in this spectrum reconstruction. If the power spectrum is parametrised by binning it in $`k`$-space, it is very difficult to mimic any smooth features in the underlying power spectrum. Especially if the feature is comparable in width to the individual bins the reconstructed spectrum can show spurious features. Nevertheless, by going to a sufficiently fine binning it is possible to map out the general shape of the power spectrum, although still very difficult to get a good fit in terms of the likelihood.

The fact that it is possible to reconstruct the underlying $`k`$-space power spectrum from CMBR observations opens up the possibility of probing physics at the time of spectrum formation by observing the universe at CMBR formation ($`T4000`$K). As was for instance discussed in Refs. , if the primordial fluctuations are produced during the inflationary epoch, it is possible to reconstruct the inflationary potential, $`V(\varphi )`$, from the $`k`$-space spectrum.

However, in addition to that it will be possible to detect non-power-law features in the power spectrum using this discrete binning method. Such features are for instance predicted to occur in the multiple inflation model of Adams, Ross and Sarkar . Here, the effect of symmetry breaking during the inflationary period was discussed. If such spontaneous symmetry breaking occurs along flat directions, short inflationary periods will result, leaving distinct non power-law features in the power spectrum. Another possible mechanism for producing features is the resonant production of particles during inflation discussed by Chung et al. .

Note that the CMBR is mainly sensitive to power spectrum at $`k10^310^1\mathrm{Mpc}^1`$. Any features outside this region in $`k`$-space would not be detectable, so the CMBR can at most probe a very small region of the inflationary period. However, as discussed for instance in Ref. , by also using data from large scale structure (LSS) surveys such as the Sloan Digital Sky Survey it is possible to increase the region of $`k`$-space that can be probed. LSS surveys are sensitive to smaller scales than the CMBR (for instance the Sloan survey should probe the region $`k10^21\mathrm{Mpc}^1`$), and by combining MAP/PLANCK with such surveys it should be possible to have sensitivity in the region $`k10^31\mathrm{Mpc}^1`$. In addition to this, including data from large scale surveys can break some of the cosmological parameter degeneracies . For instance, using CMBR data alone it is impossible to determine either $`\mathrm{\Omega }_m`$ or $`\mathrm{\Omega }_\mathrm{\Lambda }`$ separately with high precision. The parameter that can be measured is the total energy density $`\mathrm{\Omega }=\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }`$. However, the LSS surveys are sensitive to a different combination of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and by adding this information the degeneracy between them can be broken .

###### Acknowledgements.

Use of the CMBFAST code developed by Seljak and Zaldarriaga  is acknowledged.