text:
# Matching Conditions in Atomistic-Continuum Modeling of Materials

## Abstract

A new class of matching condition between the atomistic and continuum regions is presented for the multi-scale modeling of crystals. They ensure the accurate passage of large scale information between the atomistic and continuum regions and at the same time minimize the reflection of phonons at the interface. These matching conditions can be made adaptive if we choose appropriate weight functions. Applications to dislocation dynamics and friction between two-dimensional atomically flat crystal surfaces are described.

Traditionally two apparently separate approaches have been used to model a continuous medium. The first is the continuum theory, in the form of partial differential equations describing the conservation laws and constitutive relations. This approach has been impressively successful in a number of areas such as solid and fluid mechanics. It is very efficient, simple and often involves very few material parameters. But it becomes inaccurate for problems in which the detailed atomistic processes affect the macroscopic behavior of the medium, or when the scale of the medium is small enough that the continuum approximation becomes questionable. Such situations are often found in studies of properties and defects of micro- or nano- systems and devices. The second approach is atomistic, aiming at knowing the detailed behavior of each individual atom using molecular dynamics or quantum mechanics. This approach can in principle accurately model the underlying physical processes. But it is often times prohibitively expensive.

Recently an alternative approach has been explored that couples the atomistic and continuum approaches. . The main idea is to use atomistic modeling at places where the displacement field varies on an atomic scale, and the continuum approach elsewhere. Two representative examples of such a coupled approach are the quasi-continuum method  and the coarse-grained molecular dynamics . Both methods have been successfully applied to quasi-static examples, but extension to dynamic problems has not been straightforward. The main difficulty lies in the proper matching between the atomistic and continuum regions. Since the details of lattice vibrations, the phonons, which are an intrinsic part of the atomistic model, cannot be represented at the continuum level, conditions must be met that the phonons are not reflected at the atomistic-continuum interface. Since the atomistic region is expected to be a small part of the computational domain, violation of this condition quickly leads to local heating of the atomistic region and destroys the simulation result. In addition, the matching between the atomistic-continuum interface has to be such that large scale information is accurately transmitted in both directions.

The main purpose of this paper is to introduce a new class of matching conditions between atomistic and continuum regions. These matching conditions have the property that they allow accurate passage of large scale (scales that are represented by the continuum model) information between the atomistic and continuum regions and no reflection of phonon energy to the atomistic region. These conditions can also be used in pure molecular dynamics simulations as the border conditions to ensure no reflection of phonons at the boundary of the simulation.

For the sake of clarity we will first explain the main issues and ideas on a simple problem: a one-dimensional chain of particles coupled by springs:

$$\ddot{u}_j=u_{j+1}2u_j+u_{j1}$$
(1)
The spring constant is set to be 1. After discretization in time, we have

$$\frac{u_j^{n+1}2u_j^n+u_j^{n1}}{\mathrm{\Delta }t^2}=u_{j+1}^n2u_j^n+u_{j1}^n$$
(2)
where $`u_j^n`$ is the displacement of the $`j`$-th particle at time $`t=n\mathrm{\Delta }t`$.

Equation (2) is supposed to be solved for all integers $`j`$. Now let us assume that we will truncate the computational domain and only compute $`u_j^n`$ for $`j0`$. At $`j=0`$, we will impose a new boundary condition to make sure that the phonons arriving from $`j>0`$ are not reflected back at $`j=0`$.

The phonon spectrum for (2) is obtained by looking for solutions of the type $`u_j^n=e^{i(n\omega \mathrm{\Delta }t+j\xi )}`$. This gives us the relation

$$\frac{1}{\mathrm{\Delta }t}\mathrm{sin}\frac{\omega \mathrm{\Delta }t}{2}=\mathrm{sin}\xi $$
(3)
At $`j=0`$, we replace (2) by

$$u_0^n=\underset{k,j0}{}a_{k,j}u_j^{nk}a_{0,0}=0$$
(4)
We would like to determine the coefficients $`\{a_{k,j}\}`$. For the simple problem at hand, it is possible to obtain analytical formulas of $`\{a_{k,j}\}`$ such that the imposition of (4) together with the solution of (2) for $`j>0`$ reproduces exactly the solution of (2) if it was solved for all integer values of $`j`$, i.e. an exact reflectionless boundary condition can be found. The details of this is given in . This exact boundary condition should be the same as the one found numerically in . It represents the exact Green’s function for (2) which is nonlocal. However, this procedure appears to be impractical for realistic models, particularly when the atomistic region moves with time which is the case that interests us. Therefore we will not pursue this direction here.

A practical solution is to restrict (4) to a finite number of terms and look for the coefficients $`\{a_{k,j}\}`$ that minimize reflection. In order to do this, let us look for solutions of the type

$$u_j^n=e^{i(n\omega \mathrm{\Delta }t+j\xi )}+R(\xi )e^{i(n\omega \mathrm{\Delta }tj\xi )}$$
(5)
where $`\omega `$ is given by (3). $`R(\xi )`$ is the reflection coefficient at wavenumber $`\xi `$. Inserting (5) into (4), we obtain

$$R(\xi )=\frac{a_{k,j}e^{i(j\xi k\omega \mathrm{\Delta }t)}1}{a_{k,j}e^{i(j\xi +k\omega \mathrm{\Delta }t)}1}$$
(6)
The optimal coefficients $`\{a_{k,j}\}`$ are obtained by

$$\mathrm{min}_0^\pi W(\xi )|R(\xi )|^2𝑑\xi $$
(7)
subject to the constraint

$$R(0)=0,R^{}(0)=0,\text{ etc.}$$
(8)
Here $`W(\xi )`$ is a weight function, which is chosen to be $`W(\xi )=1`$ in the examples below.

Condition (8) guarantees that large scale information is transmitted accurately, whereas (7) guarantees that the total amount of reflection is minimized. This procedure offers a lot of flexibility. For example, instead of $`_0^\pi |R(\xi )|^2𝑑\xi `$, we can minimize the total reflection over certain carefully selected interval. Another possibility is to choose the weight function to be the (empirically computed) energy spectrum. The coefficients $`\{a_{k,j}\}`$ may then change in time to reflect the change of the nature of the small scales. In practice, we found it preferable to use $`_0^{\pi \delta }|R(\xi )|^2𝑑\xi `$ with some small $`\delta `$, instead of $`_0^\pi |R(\xi )|^2𝑑\xi `$, in order to minimize the influence of $`\xi =\pi `$ for which we always have $`R(\pi )=1`$.

Let us look at a few examples. If in (4) we only keep the terms involving $`a_{1,0}`$ and $`a_{1,1}`$, then imposing the condition $`R(0)=0`$ gives

$$u_0^n=(1\mathrm{\Delta }t)u_0^{n1}+\mathrm{\Delta }tu_1^{n1}$$
(9)
If instead we keep terms involving $`a_{0,1},a_{1,0}`$ and $`a_{1,1}`$, we can then impose both $`R(0)=0`$ and $`R^{}(0)=0`$. This gives us

$$u_0^n=u_1^{n1}+\frac{1\mathrm{\Delta }t}{1+\mathrm{\Delta }t}(u_0^{n1}u_1^n)$$
(10)
Conditions of the type (9) and (10) are intimately related to the absorbing boundary conditions proposed and analyzed in  for the computation of waves. These conditions perform well for low wavenumbers but are less satisfactory at high wavenumbers.

To improve the performance at high wavenumbers let us consider a case that include terms with $`k2,j3`$ and minimize $`_0^{\pi \delta }|R(\xi )|^2𝑑\xi `$ (with $`\delta =0.125\pi `$) subject to the condition $`R(0)=0`$, the optimal coefficients can be easily found numerically and are given by

$$(a_{k,j})=\left(\begin{array}{ccc}1.95264069& 7.4207\times 10^2& 1.4903\times 10^2\\ 0.95405524& 7.4904\times 10^2& 1.5621\times 10^2\end{array}\right)$$
(11)
If instead we only include terms such that $`k3,j2`$, then

$$(a_{k,j})=\left(\begin{array}{cc}2.9524& 1.5150\times 10^2\\ 2.9065& 3.0741\times 10^2\\ 0.9541& 1.5624\times 10^2\end{array}\right)$$
(12)
The resulting reflection coefficients $`R`$ are displayed in Figure 1.

We have applied this method to a number of problems. As the simplest model that encompasses most of the issues in a coupled atomistic/continuum simulation, we consider the Frenkel-Kontorova model

$$\ddot{x}_j=x_{j+1}2x_j+x_{j1}U^{}(x_j)+f$$
(13)
where $`U`$ is a periodic function with period 1, $`f`$ is an external forcing. The continuum limit of this equation is simply the Klein-Gordan equation

$$u_{tt}=u_{xx}Ku+f$$
(14)
where $`K=U^{\prime \prime }(0)`$. We consider the case when there is a dislocation and study its dynamics under a constant applied forcing. We take $`U(x)=(x[x])^2`$ where $`[x]`$ is the integer part of $`x`$. In this example we take (13) as our atomistic model, and (14) as our continuum model. For the coupled atomistic-continuum method, we use a standard second order finite difference method for (14) in the region away from the dislocation, and we use (13) in the region around the dislocation. However, we also place finite difference grid points in the atomistic region. At these points, the values are obtained through averaging the values from the atomistic model. At the interface between the atomistic and continuum regions, we decompose the displacement into a large scale and a small scale part. The large scale part is computed on the finite difference grid, using (12). The small scale part is computed using the reflectionless boundary conditions described earlier. The interfacial position between the MD and continuum regions is moved adaptively according to an analysis of the wavelet coefficients or the local stress. Both strategies lead to similar results. Care has to be exercised in order to restrict the size of the atomistic region. For example, when wavelet coefficients are used in the criteria to move the atomistic region, we found it more efficient to use the intermediate levels of the wavelet coefficients rather than the finest level.

We first consider the case when a sharp transition is made between the atomistic and continuum regions with a 1:16 ratio for the size of the grids. Figure 2 is a comparison of the displacement and velocity fields computed using the full atomistic model and the coupled atomistic/continuum model, with $`f=0.04`$. The atomistic region has 32 atoms. The full atomistic simulation has $`10^4`$. Dislocation appears as a kink in the displacement field. Notice that at the atomistic/continuum interface, there is still substantial phonon energy which is then suppressed by the reflectionless boundary condition. No reflection of phonons back to the atomistic region is observed.

We next consider a case with $`f=0.02`$, which alone is too weak to move the dislocation, but to the left of the dislocation, we add a sinusoidal wave to the initial data. The dislocation moves as a consequence of the combined effect of the force and the interaction with the wave. Yet in this case the same atomistic/continuum method predicts an incorrect position for the dislocation, as shown in Figure 3. The discrepancy seems to grow linearly in time. Improving the matching conditions does not seem to lead to significant improvement.

The difference between this case and the case shown in Figure 2 is that there is substantially more energy at the intermediate scales. This is clearly shown in the energy spectrum that we computed for the two cases but it can also be seen in Figure 3 where an appreciable amount of small scale waves are present in front of the dislocation. Such intermediate scales are suppressed in a method that uses a sharp transition between the atomistic and continuum regions. We therefore consider the next alternative in which the atomistic/continuum transition is made gradually in a 1:2 or 1:4 ratio between neighboring grids. The right column in Figure 3 shows the results of such a method that uses a gradual 1:2 transition. We see that the correct dislocation position is now recovered.

Our second example is the friction between atomically flat crystal surfaces. To model this process atomistically, we use standard molecular dynamics with the Lennard-Jone potential . The two crystals are separated by a horizontal interface. The atoms in the bottom crystal are assumed to be much heavier than the atoms on top. A constant shear stress is applied at the top surface. The main issue here is how dissipation takes place. Physically the mean kinetic energy is dissipated through conversion into phonons which then convert into heat and exit the system. A standard practice in modeling such a process is to add a friction term to the molecular dynamics in order to control the temperature of the system . In contrast, we ensure the proper dissipation of phonons to the environment by imposing the reflectionless boundary conditions for the phonons. As a result we obtain a linear relationship between the mean displacement of the atoms in the top crystal as function of time, see Figure 4. The temperature of the system also saturates (Figure 4). Also plotted in Figure 4 is the result of the mean displacement computed using the combined atomistic/continuum method. Here the continuum model is the linear elastic wave equation with Lame coefficients computed from the Lennard-Jones potential. The agreement between the full atomistic and the atomistic/continuum simulation is quite satisfactory.

In conclusion, we presented a new strategy for the matching condition at the atomistic/continuum interface in multiscale modeling of crystals. These conditions are adaptive if we choose the weight functions in (7) to reflect the evolving nature of the small scales. They minimize the reflection of phonons and at the same time ensure accurate passage of large scale information.

This work is supported by NSF through a PECASE award and by ONR grant N00014-01-1-0674.


format: text

text:
# On a possible GRB-supernova time sequence

## Abstract

The data from the Chandra satellite on the iron emission lines in the afterglow of GRB 991216 are used to give further support for the EMBH theory, which links the origin of the energy of GRBs to the extractable energy of electromagnetic black holes (EMBHs), leading to an interpretation of the GRB-supernova correlation. Following the relative space-time transformation (RSTT) paradigm and the interpretation of the burst structure (IBS) paradigm, we introduce a paradigm for the correlation between GRBs and supernovae. The following sequence of events is shown as kinematically possible and consistent with the available data: a) the GRB-progenitor star $`P_1`$ first collapses to an EMBH, b) the proper GRB (P-GRB) and the peak of the afterglow (E-APE) propagate in interstellar space until the impact on a supernova-progenitor star $`P_2`$ at a distance $`2.69\times 10^{17}`$ cm, and they induce the supernova explosion, c) the accelerated baryonic matter (ABM) pulse, originating the afterglow, reaches the supernova remnants $`18.5`$ hours after the supernova explosion and gives rise to the iron emission lines. Some considerations on the dynamical implementation of the paradigm are presented. The concept of induced supernova explosion introduced here specifically for the GRB-supernova correlation may have more general application in relativistic astrophysics.

black holes, gamma ray bursts, supernovae  11affiliationtext: ruffini@icra.it

We have seen in the previous two letters how the fit of the data from the RXTE (Corbet & Smith, 2000) and Chandra (Piro et al., 2000) satellites on the afterglow of the GRB 991216 offers a tool to determine the only two free parameters of the EMBH theory. This theory links the energy source of GRBs to the electromagnetic mass energy of black holes (Christodoulou & Ruffini, 1971). We have also seen how this theory has consequences for the interpretation of the structure of GRBs (the IBS paradigm, see Ruffini et al., 2001b). The same analysis also yields information about the density and overall distribution of the baryonic matter in the remnant left over from the gravitational collapse of the GRB-progenitor star to an EMBH, see Ruffini et al. (2001a) and Ruffini et al. (2001b), and see also Ruffini et al. (2001g). Similarly it allows one to probe the density distribution of the interstellar matter around the newly formed EMBH (see e.g. Ruffini (2001); Ruffini et al. (2001e)).

The aim of this letter is to point out that the data on the iron lines from the Chandra satellite on the GRB 991216 (Piro et al., 2000) and similar observations from other sources (Piro et al., 1999; Amati et al., 2000; Piro et al., 2000) make it possible to extend this analysis to a larger distance scale, possibly all the way out to a few light years, and consequently probe the distribution of stars in the surroundings of the newly formed EMBH. These considerations lead to a new paradigm for the interpretation of the supernova-GRB correlation.

That indeed a correlation between the occurrence of GRBs and supernova events exists has been established by the works of Bloom et al. (1999); Galama et al. (1998a, b, 2000); Kulkarni et al. (1998); Piro et al. (1998); Pian et al. (1999); Reichart (1999); van Paradijs et al. (2000). Such an association has been assumed to indicate that GRBs are generated by supernovae explosions (see e.g. Kulkarni et al., 1998).

We propose that, if relativistic effects are properly taken into account, then an alternative, kinematically viable, explanation can be given of the supernova-GRB association. We again use GRB 991216 as a prototypical case. The same theoretical considerations have been also applied to other cases, including GRB 980425 and SN 1998BW (Ruffini et al., 2001d).

We focus on the detailed kinematical description of this GRB-supernova time-sequence process and outline a possible dynamical scenario. We introduce a process by which a massive GRB-progenitor star $`P_1`$ of mass $`M_1`$ undergoes gravitational collapse to an EMBH. During this process a dyadosphere is formed and subsequently the P-GRB and the E-APE of Ruffini et al. (2001b) are generated in sequence (see also Ruffini et al., 2001a). They propagate and impact, with their photon and neutrino components, on a second supernova-progenitor star $`P_2`$ of mass $`M_2`$. Assuming that both stars were generated approximately at the same time, we expect to have $`M_2<M_1`$. For a wide range of parameters, such a collision will not affect the star $`M_2`$ (Chardonnet & Ruffini, 2001). Under some special conditions of the thermonuclear evolution of the supernova-progenitor star $`P_2`$, the collision can induce a supernova explosion.

We assume that the star $`P_2`$ is close to the line of sight of the EMBH. We will see in the following that this gives an upper limit to the distance $`D_{P_2}=2.69\times 10^{17}`$ cm of the supernova-progenitor star $`P_2`$ from the EMBH. The location of the star $`P_2`$ will then be constrained between the transparency point of the P-GRB, $`9.692\times 10^{13}`$ cm (see Ruffini et al., 2001a), and the above upper limit, and will be a function of the angle subtended by the line of sight and the star $`P_2`$, as seen from the EMBH. The energy-momentum deposited by the GRB in the collision with the star is in the range $`10^{39}10^{45}`$ ergs (see e.g. Ruffini et al., 2001h).

Especially relevant to our model are the following data from the Chandra satellite (see Piro et al., 2000):
1) At the arrival time of 37 hr after the initial burst there is evidence of iron emission lines for GRB 991216.
2) The emission lines are present during the entire observation period of $`10^4`$ sec. The iron lines could also have been produced earlier, before Chandra was observing. Thus the times used in these calculations are not unique: they do serve to provide an example of the scenario.
3) The emission lines appear to have a peak at an energy of $`3.49\pm 0.06`$ keV which, at a redshift $`z=1.00\pm 0.02`$ corresponds to an hydrogen-like iron line at 6.97 keV at rest. This source does not appear to have any significant motion departing from the cosmological flow. The iron lines have a width of 0.23 keV consistent with a radial velocity field of $`0.1c`$.

From the theoretical slope of the afterglow, presented in Ruffini et al. (2001b), we see that the flux of the afterglow observed by Chandra is in excellent agreement with the general afterglow slopes. Clearly the iron lines are only a small fraction of the observed flux.

We assume the laboratory frame as an inertial system of reference in which both stars $`P_1`$ and $`P_2`$ are at rest. A second asymptotic inertial reference frame is assumed in which the detector is at rest. Therefore a detector arrival time $`t_a^d`$ is defined related to the laboratory time $`t`$ by

$$t_a^d=\left(1+z\right)t_a=\left(1+z\right)\left(t\frac{r\left(t\right)}{c}\right),$$
(1)
where $`z`$ is the cosmological redshift, which in the case of GRB 991216 is equal to $`z=1.0`$, and $`r\left(t\right)`$ is the radius of the expanding pulse at laboratory time $`t`$.

To explain the above observations we propose that the expansion of the accelerated baryonic matter (ABM pulse see Ruffini et al., 2001a) relativistically expanding away from newly formed EMBH reaches $`P_2`$ with a delay in arrival time of 18.5 hr (the details of the computation are given in Bianco et al., 2001b). The associated afterglow then illuminates the expanding supernova shell, producing the observed iron emission lines. The Chandra satellite observations then offer the first data on such an induced-supernova-explosion (Piro et al., 2000).

On the basis of the explicit computations of the different eras presented in Ruffini et al. (2001a), we make three key points:
1) An arrival time of 37 hr in the detector frame corresponds to a radial distance from the EMBH traveled by the ABM pulse of $`2.69\times 10^{17}`$ cm in the laboratory frame (see Ruffini et al., 2001a).
2) It is likely that a few stars are present within that radius as members of a cluster. It has become evident from observations of dense clusters of star-forming regions that a stellar average density of typically $`10^2\text{pc}^3`$ (Beck et al., 2000) should be expected. There is also the distinct possibility for this case and other systems that the stars $`P_1`$ and $`P_2`$ are members of a detached binary system.
3) The possible observations at different wavelengths (Ruffini et al., 2001d) crucially depend on the relative intensities between the GRB and the supernova as well as on the value of the distance and the redshift of the source.

In order to reach an intuitive understanding of these complex computations we present a schematic very simplified diagram (not to scale) in Fig. 1.

We now describe the specific data of this GRB-supernova time-sequence (GSTS) paradigm:

1) The two stars $`P_1`$ and $`P_2`$ are separated by a distance $`D_{P_2}=2.69\times 10^{17}`$ cm in the laboratory frame, see Fig. 1. Both stars are at rest in the inertial laboratory frame. At laboratory time $`t=0`$ and at comoving time $`\tau =0`$, the gravitational collapse of the GRB-progenitor star $`P_1`$ occurs. The initial emission of gravitational radiation or a neutrino burst from the event then synchronizes the arrival times $`t_a=0`$ for the supernova-progenitor star $`P_2`$ and $`t_a^d=0`$ for the distant observer at rest with the detector. The electromagnetic radiation emitted by the gravitational collapse process is instead practically zero, by comparison, due to the optical thickness of the material at this stage (Bianco et al., 2001a).

2) From the determination of the parameters obtained in Ruffini et al. (2001b) and the computations in Ruffini et al. (1999, 2000), at laboratory time $`t_1=3.295\times 10^3`$ s and at a distance from the EMBH of $`D_1=9.692\times 10^{13}`$ cm, the condition of transparency for the pair-electromagnetic-baryon (PEMB) pulse is reached and the P-GRB is emitted, see Fig. 1, Ruffini et al. (2001a) and Ruffini et al. (2001b). This time is recorded in arrival time at the detector $`t_{a_1}^d=0.1361`$ sec, and, at $`P_2`$, at $`t_{a_1}=6.805\times 10^2`$ sec.

The fact that the PEMB pulse in an arrival time of 0.1361 sec covers a distance of $`9.692\times 10^{13}`$ cm gives rise to an apparent superluminal effect. In order to clarify this apparent inconsistency we have introduced, from the computed values of $`t_1`$ and $`t_{a_1}^d`$, an “effective” Lorentz gamma factor by Eq.(2)

$$\gamma \sqrt{\frac{t}{2t_a^d}},$$
(2)
which gives $`\gamma _1110.0`$. Then we can straightforwardly explain the difference between the two times as $`t_a^d=t/\left(2\gamma ^2\right)`$, see Fig. 1.

This introduction of an “effective” Lorentz gamma factor has no predictive power, it can only be introduced a posteriori, as an heuristic tool in order to draw the qualitative diagram in Fig. 1. In practice the entire integration of the equations must be accomplished taking into account all changes of the time varying Lorentz gamma factors and the corresponding space and time variables throughout each era (see for details Bianco et al., 2001b; Ruffini et al., 2001f).

3) At laboratory time $`t=1.653\times 10^6`$ s and at a distance from the EMBH of $`4.863\times 10^{16}`$ cm in the laboratory frame, the peak of the E-APE is reached which is recorded at the arrival time $`t_a=11.86`$ sec at $`P_2`$ and $`t_a^d=23.72`$ s at the detector. This also gives rise to an apparent superluminal effect which can be also explained following the same arguments given above in point 2). This event has not been represented in Fig. 1 in order not to confuse the image.

4) At a distance $`D_{P_2}=2.69\times 10^{17}`$ cm, the two bursts described in the above points 2) and 3) collide with the supernova-progenitor star $`P_2`$ at arrival times $`t_{a_1}=6.805\times 10^2`$ s and $`t_a=11.86`$ s respectively. They can then induce the supernova explosion of the massive star $`P_2`$.

5) The associated supernova shell expands with velocity $`0.1c`$.

6) The expanding supernova shell is reached by the ABM pulse generating the afterglow with a delay of $`t_{a_2}=18.5`$ hr in arrival time following the arrival of the P-GRB and the E-APE. This time delay coincides with the interval of laboratory time separating the two events, since the $`P_2`$ is at rest in the inertial laboratory frame (see Fig. 2).

Again as explained above in point 2), this time delay can be interpreted a posteriori by introducing (Bianco et al., 2001b) an “effective” Lorentz gamma factor defined by

$$\gamma \sqrt{\frac{D_{P_2}}{2ct_{a_2}}},$$
(3)
which gives $`\gamma 8.21`$. Then we can heuristically visualize the time delay as $`t_{a_2}^d=D_{P_2}/\left(2c\gamma ^2\right)`$. Clearly the results presented in Fig. 2 do follow the complete integration of the equations of motion of the system through each different era defined in Ruffini et al. (2001a).

The ABM pulse will have travelled in the laboratory frame a distance $`D_{P_2}D_1D_{P_2}=2.69\times 10^{17}`$ cm in a laboratory time $`t_2t_1t_2=9.02\times 10^6`$ s (neglecting the supernova expansion).

This again gives rise to an apparent superluminal effect which can be interpreted heuristically, as in point 2), as a relativistic motion of the ABM pulse with an “effective” Lorentz factor of $`\gamma _25.82`$, see Fig. 1. The era IV extends from the point of transparency (point 4 in Fig. 1 of Ruffini et al., 2001a) all the way to the collision of the pulse with the supernova shell, which occurs at $`\gamma 3.38`$. By this time the supernova shell has reached a dimension of $`1.997\times 10^{14}`$ cm, which is consistent with the observations from the Chandra satellite.

In the above considerations of GRB 991216 the supernova remnant has been assumed to be close to, but not exactly along, the line of sight extending from the EMBH to the distant observer. However, such a case should exist for other GRBs and would lead to an observation of iron absorption lines as well as to an increase in the radiation observed in the afterglow corresponding to the crossing of the supernova shell by the ABM pulse. In fact, as the ABM pulse engulfs the baryonic matter of the remnant, above and beyond the normal interstellar medium baryonic matter, the conservation of energy and momentum implies that a larger amount of internal energy is available and radiated in the process (Ruffini et al., 2001e). This increased energy-momentum loss will generally affect the slope of the afterglow decay, approaching more rapidly a nonrelativistic expansion phase (details are give in Ruffini et al., 2001f).

If we now turn to the possibility of dynamically implementing the scenario, there are, at least, three different possibilities:
1) Particularly attractive is the possibility that a massive star $`P_2`$ has rapidly evolved during its thermonuclear evolution to a white dwarf (see e.g. Chandrasekhar, 1978). It it then sufficient that the P-GRB and the E-APE implode the star sufficiently as to reach a central density above the critical density for the ignition of thermonuclear burning. Consequently, the explosion of the star $`P_2`$ occurs, and a significant fraction of a solar mass of iron is generated. These configurations are currently generally considered precursors of some type I supernovae (see e.g. Filippenko, 1997, and references therein).
2) Alternatively, the massive star $`P_2`$ can have evolved to the condition of being close to the point of gravitational collapse, having developed the formation of an iron-silicon core, type II supernovae. The above transfer of energy momentum from the P-GRB and the E-APE may enhance the capture of the electrons on the iron nuclei and consequentely decrease the Fermi energy of the core, leading to the onset of gravitational instability (see e.g. Bethe, 1991, p. 270 and followings). Since the time for the final evolution of a massive star with an iron-silicon core is short, this event will require a perhaps unlikely coincidence.
3) The pressure wave may trigger massive and instantaneous nuclear burning process, with corresponding changes in the chemical composition of the star, leading to the collapse.

The GSTS paradigm has been applied to the case of the correlation between SN 1998bw and GRB 980425, which, with a redshift of 0.0083, is one of the closest and weakest GRBs observed. In this case the EMBH appears to have a significantly lower value of the parameter $`\xi `$ but the validity of the GSTS paradigm presented here is fully confirmed (see Ruffini et al., 2001d).

The GSTS paradigm and the concept of induced supernova explosion, which we have introduced for the collapse to an EMBH, may play a role also in the case of a collapse of a white dwarf core to a neutron star in a binary system. It may solve the long lasting problem of the almost equality of neutron star masses observed in some binary pulsars (see e.g. Taylor & Weisberg, 1989).

We thank three anonymous referees and J. Wilson for their remarks, which have improved the presentation of this letter


text:
# Modified Breit-Wigner formula for mesonic resonances describing OZI decays of confined 𝑞𝑞̄ states and the light scalar mesons

## 1 Introduction

Lattice QCD in principle offers the most direct way to link to experiment what we believe to be the fundamental theory of strong interactions. However, for the time being only quenched calculations are capable of making, with reasonable accuracy, predictions for e.g. masses of mesons, as a result of the confinement sector of QCD . Relating such predictions to experimental data is quite a controversial issue though, since the states obtained in quenched lattice QCD (qlQCD) are manifestly stable, a circumstance which is clearly not the case in experiment, at least for most mesons. The reason of course is the difficulty of incorporating quark-pair creation in qlQCD, thus impeding the process of OZI-allowed strong decay, which is responsible for the large widths of many mesonic resonances. But also the central positions of such resonances may be quite displaced due to the effects of strong decay , when compared with the stable qlQCD states. Moreover, the contributions of virtual hadronic decays through quark loops owing to the presence of closed OZI-allowed thresholds, which should also lead to real mass shifts, are not fully accounted for in qlQCD . An additional complication is the observation that not even the number of $`q\overline{q}`$ states in the $`J^P=0^+`$ sector of the qlQCD predictions seems to agree with experiment , thus contributing to the general confusion concerning especially the light scalar mesons.

It is evidently a very unsatisfactory state of affairs that, notwithstanding the ever improving accuracy of the numerical predictions in qlQCD, these cannot be trusted to unmistakably confirm possible signals of new physics. In the present paper, we shall propose an alternative method of data analysis, which does allow for an accurate link between qlQCD, or any other, model-dependent formulation of confinement, and experiment. In the process we are also going to find a reliable approach to the light scalar mesons ( $`<1`$ GeV).

The idea is simple and amounts to the observation that valence-quark-pair creation connects the states of qlQCD to the resonances which in experiment are seen in elastic and inelastic meson-meson cross sections. Hence, if we model quark-pair creation in such a way that it theoretically can be turned on and off, then we are capable to predict the qlQCD states in the model limit of no-quark-pair creation. Such a philosophy already underlied an elaborate coupled-channel quark model , designed to simultaneously describe mesonic bound-state spectra, resonances, and meson-meson scattering. However, in spite of the model’s success to reproduce a host of experimental data with a very limited number of parameters, it is clearly not suited as a tool for data analysis, owing to the specific model choice of the confining $`q\overline{q}`$ potential, and furthermore the rather complicated matrix expressions needed to obtain $`S`$-matrix-related observables. Therefore, in this work our strategy will be the following. By replacing the specific confinement part of the Hamiltonian of Refs.  by a more general one, we allow here for any arbitrary discrete spectrum of “bare” $`q\overline{q}`$ states. At the same time, by using a simpler form for the coupling potential, describing the transitions between the $`q\overline{q}`$ and meson-meson sectors through $`{}_{}{}^{3}P_{0}^{}`$ quark-pair creation, the complexity of the model’s scattering solutions is substantially reduced. The resulting exact, closed-form formula for the $`S`$-matrix is subsequently fitted to the experimental data, i.e., partial-wave cross sections or phase shifts for meson-meson scattering, by adjusting the model parameters. Afterwards, we may study the theoretical limit of vanishing valence-quark-pair creation, described by a sole parameter, which decouples the bare states from the meson-meson continuum. The latter are supposed to be equivalent to the states found in qlQCD. Consequently, the masses of the thus resulting spectrum could then be compared with those from qlQCD calculations.

The organisation of the present paper is as follows. In section 2 we develop a model-independent partial-wave $`S`$-matrix for the description of elastic and inelastic two-meson scattering. A so-called “Resonance-Spectrum Expansion” (RSE) is discussed in section (2.2). In section 3 the RSE is compared to the data for elastic $`K\pi `$ $`S`$-wave phase shifts and $`P`$-wave cross sections in the one-threshold limit. The complex singularities of the corresponding $`S`$-matrices are shown to allow for an easy relation to the $`q\overline{q}`$ confinement spectrum. The conclusions are presented in section 4.

## 2 Breit-Wigner-like scattering amplitudes

A suitable model in which the communication between the confinement sector of strong interactions with the two-meson continuum is enabled through valence-quark-pair creation has been proposed in a series of articles . It allows for the determination of partial-wave two-meson elastic and inelastic scattering cross sections, as well as for the search of singularities in the partial-wave coupled-channel two-meson scattering matrix for all possible valence flavours. Hence, it enables the calculation of resonances above and bound states below the lowest threshold, i.e., the meson spectrum. In the limit of no valence-quark-pair creation one obtains the confinement spectrum, or bare states, which may be compared to the states of qlQCD. The model has four parameters and four constituent quark masses. No distinction is made between up and down quarks. One of the model parameters parametrises the employed confinement mechanism (harmonic oscillator), whereas the other three parametrise the communication between the two distinct sectors of the model, that is, the permanently closed confinement channels and the meson-meson continuum channels. The model yields, with one set of four parameters and one set of four constituent quark masses, good results for heavy quarkonia  and light-meson spectra , including the scalars , as well as for two-meson scattering data .

Let us first study the generic form of the scattering matrix for permanently closed channels coupled to several meson-meson scattering channels, in a simplified version of the above model.

### 2.1 Scattering matrix for several coupled channels

When we describe quarkonia by wave functions $`\psi _c`$ and two-meson systems by wave functions $`\psi _f`$, then we obtain for their time evolution the wave equation

$$\left(EH_c\right)\psi _c\left(\stackrel{}{r}\right)=V_t\psi _f\left(\stackrel{}{r}\right)\text{and}\left(EH_f\right)\psi _f\left(\stackrel{}{r}\right)=\left[V_t\right]^T\psi _c\left(\stackrel{}{r}\right).$$
(1)
Here, $`H_c`$ describes the dynamics of confinement in the interaction region, $`H_f`$ the dynamics of the scattered particles, and $`V_t`$ the communication between the two different sectors.

For the dynamics of confinement we understand here that, as a function of the interquark distance $`r`$, the resulting quark-antiquark binding forces grow rapidly outside the interaction region. Consequently, we must eliminate $`\psi _c`$ from Eqs. (1), since it is vanishing at large distances and thus unobservable. Formally, it is easy to do so. We then obtain the relation

$$\left(EH_f\right)\psi _f\left(\stackrel{}{r}\right)=\left[V_t\right]^T\left(EH_c\right)^1V_t\psi _f\left(\stackrel{}{r}\right).$$
(2)
By comparison of Eq. (2) with the usual expressions for the scattering wave equations, we must conclude that the generalised potential $`V`$, which results from the set of coupled equations (1), is here given by

$$V=\left[V_t\right]^T\left(EH_c\right)^1V_t.$$
(3)
In the momentum representation Eq. (3) takes the form

$$\stackrel{}{p}|V|\stackrel{}{p}^{}=\stackrel{}{p}|\left[V_t\right]^T\left(EH_c\right)^1V_t|\stackrel{}{p}^{}.$$
(4)
Let us denote the configuration-space representation of the properly normalised eigensolutions of the operator $`H_c`$ of Eqs. (1), corresponding to the energy eigenvalue $`E_{n\mathrm{}_c}`$, by

$$\stackrel{}{r}|n\mathrm{}_cm=Y_m^{(\mathrm{}_c)}\left(\widehat{r}\right)_{n\mathrm{}_c}\left(r\right),\text{where}\{\begin{array}{c}n=0\text{}1\text{}2\text{}\mathrm{},\hfill \\ \mathrm{}_c=0\text{}1\text{}2\text{}\mathrm{},\hfill \\ m=\mathrm{}_c\text{}\mathrm{}\text{}+\mathrm{}_c.\hfill \end{array}$$
(5)
Here, $`n`$ and $`\mathrm{}_c`$ represent the orbital radial and angular excitations of the quark-antiquark system, respectively. Hence, when we let the self-adjoint operator $`H_c`$ act to the left, we obtain for Eq. (4) the expression

$`\stackrel{}{p}|V|\stackrel{}{p}^{}`$ $`=`$ $`{\displaystyle \underset{n\mathrm{}_cm}{}}\stackrel{}{p}|\left[V_t\right]^T|n\mathrm{}_cmn\mathrm{}_cm|\left(EH_c\right)^1V_t|\stackrel{}{p}^{}`$ (6)
$`=`$ $`{\displaystyle \underset{n\mathrm{}_cm}{}}\stackrel{}{p}|\left[V_t\right]^T{\displaystyle \frac{|n\mathrm{}_cmn\mathrm{}_cm|}{EE_{n\mathrm{}_c}}}V_t|\stackrel{}{p}^{}.`$
Evaluation of this equation leads to the Born term of the transition amplitude.

However, now we find it opportune to select the operators $`H_f`$ and $`V_t`$ such that it becomes possible to determine all higher-order terms of the transition amplitude. By doing so, we leave no doubt about the analyticity and unitarity properties of the resulting scattering matrix. In configuration space we define these operators by

$$H_f=\frac{1}{2}\mu ^1_r^2+M_1+M_2,\text{and}V_t=\frac{\lambda }{a^{3/2}}\overline{V}_t\delta \left(ra\right),$$
(7)
where $`\mu `$ represents the reduced-mass matrix of the meson-meson system, and $`M_{1,2}`$ stand for matrices that contain the masses of the two mesons in each scattering channel. We limit ourselves to diagonal mass matrices in this work.

The transition potential $`V_t`$ in Eqs. (7) is, as we shall see below, a reasonable approximation to the quark-pair-creation transition potentials described in Ref. . It is parametrised by $`\lambda `$, which determines its intensity, and by $`a`$, which stands for the average distance between the interacting particles (either a quark and an antiquark, or two mesons) where transitions from one sector to the other take place. In practice, $`a`$ should come out of the order of 1 fm, which is confirmed by adjusting the model parameters to the experimental data. $`\overline{V}_t`$ is the matrix that contains the relative intensities for transitions between the meson-meson and quark-antiquark sector(s). Notice that we assume here spherical symmetry for all interactions.

With the choices of Eqs. (7), we obtain for the full (to all orders in $`\lambda `$) partial-wave scattering matrix the exact expression

$`S_{\mathrm{}}\left(E\right)=\left[1\mathrm{\hspace{0.33em}2}i{\displaystyle \frac{\lambda ^2}{a}}p^{1/2}H^{(2)}(a)\mu ^{1/2}\left[\overline{V}_t\right]^T{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left|_{n\mathrm{}_c}(a)\right|^2}{EE_{n\mathrm{}_c}}}\overline{V}_t\mu ^{1/2}J(a)p^{1/2}\right]`$
$`\times \left[1+\mathrm{\hspace{0.33em}2}i{\displaystyle \frac{\lambda ^2}{a}}p^{1/2}H^{(1)}(a)\mu ^{1/2}\left[\overline{V}_t\right]^T{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left|_{n\mathrm{}_c}(a)\right|^2}{EE_{n\mathrm{}_c}}}\overline{V}_t\mu ^{1/2}J(a)p^{1/2}\right]^1,`$
where $`p`$, $`\mu `$, $`J(a)`$ and $`H^{(1,2)}(a)`$ are diagonal matrices throughout this work, with as many diagonal elements as meson-meson channels considered. The non-vanishing matrix element

$$p_i=\left[p\right]_{ii},$$
(9)
represents the relative linear momentum in the centre-of-mass (CM) system of the $`i`$-th scattering channel. The diagonal elements of $`J(a)`$ and $`H^{(1,2)}(a)`$ are related to the usual spherical Bessel and Hankel functions by

$$\left[J(a)\right]_{ii}=p_iaj_\mathrm{}_i\left(p_ia\right)\text{and}\left[H^{(1,2)}(a)\right]_{ii}=p_iah_\mathrm{}_i^{(1,2)}\left(p_ia\right),$$
(10)
where $`\mathrm{}_i`$ stands for the relative angular momentum in the $`i`$-th scattering channel.

The matrix $`\overline{V}_t`$ contains the coupling constants which are worked out in Ref. . In case only one $`q\overline{q}`$ channel is considered ($`S`$-wave meson-meson scattering for isodoublet and isovector), the matrix $`\overline{V}_t`$ is just a row vector. Then the expression

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left|_{n\mathrm{}_c}(a)\right|^2}{EE_{n\mathrm{}_c}}$$
(11)
is just a real number, that is, a function of the total CM energy $`E`$.

In case one considers more $`q\overline{q}`$ channels ($`S`$-wave meson-meson scattering for isoscalar coupled $`n\overline{n}`$ and $`s\overline{s}`$ channels, or $`P`$\- and higher-wave scattering), $`\overline{V}_t`$ has as many rows as $`q\overline{q}`$ channels. Then the resonance sum (11) is a real matrix of the size of the number of $`q\overline{q}`$ channels.

### 2.2 The Resonance-Spectrum Expansion

Expression (2.1) contains very little model dependence. It combines simple kinematics with the experimental observation that resonances occur in non-exotic scattering processes of mesons. Not even assumptions are made about possible final-state interactions, which is expressed by the choice for $`H_f`$ in Eqs. (7). Consequently, since it is not contaminated with model-dependent prejudices, expression (2.1) is extremely suitable for the analysis of experimental results in two-meson scattering. Precise determination of the experimental values for $`E_{n\mathrm{}_c}`$ and $`\left|_{n\mathrm{}_c}(a)\right|^2`$ in the resonance sum (11) will give support to the study of the confinement dynamics and the mechanism of hadronic decay. Below we shall show how well the procedure works for data analysis.

In Refs.  all possible pseudoscalar-pseudoscalar, pseudoscalar-vector, and vector-vector scattering channels are coupled, through $`{}_{}{}^{3}P_{0}^{}`$ nonstrange and strange quark-pair creation, to the relevant valence quark-antiquark channels. For isovector and isodoublet flavours one then has one (for $`S`$-wave scattering) or two (for $`P`$ and higher waves) permanently closed channels coupled to many scattering channels. For the light isoscalars the number of permanently closed channels is doubled, one channel for $`n\overline{n}`$ and one for $`s\overline{s}`$. The intensities of the relative couplings for the transitions between permanently closed channels and the various scattering channels are controlled by flavour independence, which is an observed property of strong interactions . This has been worked out in Ref. , and in some more detail in Refs. . As a result, also in Refs.  only one overall parameter is left for all possible transition intensities, which can be switched on and off. The behaviour of this more complex model is, especially near the lowest threshold, very similar to the behaviour of the scattering matrix given in formula (2.1).

In Refs. , a harmonic oscillator with flavour-independent frequency was chosen for the description of the confinement dynamics in the permanently closed channels. Hence, by switching off the overall transition parameter, one obtains the harmonic-oscillator spectrum. On the other hand, by switching it on, the experimental data for two-meson scattering are reproduced, as well as bound states like the Kaon and the $`J/\mathrm{\Psi }`$. Here, we do not intend to specify the confinement Hamiltonian $`H_c`$ of Eqs. (1), but shall follow a different strategy. We observe that the expressions $`_{n\mathrm{}_c}\left(a\right)`$ in formula (2.1), i.e., the values of the eigenfunctions of the confinement operator $`H_c`$ at distance $`a`$, are $`c`$-numbers independent of the total CM energy $`E`$. In model , $`_{n\mathrm{}_c}(a)`$ and $`E_{n\mathrm{}_c}`$ represent the harmonic-oscillator eigenstates and eigenvalues, respectively, which need only one free parameter, the oscillator frequency. However, since the confinement mechanism is here supposed to be unknown, we might as well substitute $`\left|_{n\mathrm{}_c}\left(a\right)\right|^2`$ and $`E_{n\mathrm{}_c}`$ by arbitrary non-negative real constants, $`B_{n\mathrm{}_c}`$ and $`E_{n\mathrm{}_c}`$, to be adjusted to the experimental data, i.e.,

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left|_{n\mathrm{}_c}(a)\right|^2}{EE_{n\mathrm{}_c}}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{B_{n\mathrm{}_c}}{EE_{n\mathrm{}_c}}.$$
(12)
In practical calculations, one may limit the sum in this “Resonance-Spectrum Expansion” (RSE) to a few ($`=N`$) resonances only, and approximate the sum of the remaining terms by a constant, assuming $`EE_{n\mathrm{}_c}`$ for $`n>N`$. This way one obtains Breit-Wigner-like expressions.

By redefining $`\lambda `$ and the $`B_{n\mathrm{}_c}`$’s one might take the above-referred constant equal to $`1`$, according to

$$\lambda ^2\underset{n=0}{\overset{\mathrm{}}{}}\frac{B_{n\mathrm{}_c}}{EE_{n\mathrm{}_c}}\lambda ^2\left\{\underset{n=0}{\overset{N}{}}\frac{B_{n\mathrm{}_c}}{EE_{n\mathrm{}_c}}\mathrm{\hspace{0.33em}1}\right\}.$$
(13)
An alternative approach is to absorb $`\lambda ^2`$ into the $`B_{n\mathrm{}_c}`$’s and then separate the relevant terms and the remaining sum. We shall not follow this strategy, since we want to keep explicit the dependence on the parameter which provides the communication between the scattering and confinement sectors.

At this point it is opportune to discuss the model dependence of our procedure. By the substitution (12), any relation to the quantum numbers of the confinement sector is lost. We just continue to label the $`B_{n\mathrm{}_c}`$’s and $`E_{n\mathrm{}_c}`$’s in order to distinguish them properly. It also implies that a direct reference to the $`{}_{}{}^{3}P_{0}^{}`$ mechanism is lost. What is left is just the choice (7) for the remaining operators, where $`H_f`$ does not even contain final-state interactions, and $`V_t`$ is a spherically symmetric local approximation to almost any possible transition mechanism that provides the experimentally observed OZI-allowed communication between the confinement and meson-meson scattering sectors. In summary, the model only assumes that non-exotic meson-meson scattering is dominated by the coupling to $`s`$-channel resonances.

## 3 One threshold

Let us consider the case of one permanently closed channel coupled to one meson-meson scattering channel. Using formulae (2.1) and (12), one deduces the partial-wave scattering phase shift $`\delta _{\mathrm{}}(p)`$ for elastic meson-meson scattering, which is a function of the relative momentum in the CM frame,

$$\text{cotg}\left(\delta _{\mathrm{}}(p)\right)=\frac{n_{\mathrm{}}(pa)}{j_{\mathrm{}}(pa)}\left[2\lambda ^2\mu paj_{\mathrm{}}^2(pa)\underset{n=0}{\overset{\mathrm{}}{}}\frac{B_{n\mathrm{}_c}}{EE_{n\mathrm{}_c}}\right]^1,$$
(14)
with $`\overline{V}_t`$ absorbed in $`\lambda `$.

Notice that the partial-wave phase shifts vanish for $`\lambda 0`$, unless one takes at the same time the limit $`EE_{n\mathrm{}_c}`$, which represents therefore the no-interaction limit in the two-meson system. It is indeed the parameter $`\lambda `$ that switches valence-quark-pair creation on and off.

Formula (14) has similar features as standard Breit-Wigner approximations  for resonant phenomena. However, from the values of $`E_{n\mathrm{}_c}`$ and $`B_{n\mathrm{}_c}`$, one cannot read off the positions of the singularities. At most, one might determine an approximate formula which is good for small values of $`\lambda `$, namely

$$E_{n\mathrm{}_c}B_{n\mathrm{}_c}\left[\underset{n^{}n}{}\frac{B_{n^{}\mathrm{}_c}}{E_{n\mathrm{}_c}E_{n^{}\mathrm{}_c}}\frac{i}{2\lambda ^2\mu paj_{\mathrm{}}(pa)h_{\mathrm{}}^{(1)}(pa)}\right]^1,$$
(15)
concerning singularities in the complex-energy plane. Their precise locations can be determined numerically, starting from the approximate values (15).

The values of $`E_{n\mathrm{}_c}`$ correspond to the confinement spectrum. These are the quantities of interest in this work. Let us consider next the case of $`K\pi `$ scattering, as an example.

### 3.1 $`𝑲𝝅`$ $`𝑷`$-wave scattering

Since the pion is a rather light particle, we prefer to employ relativistic kinematics for the relation between the relative linear momentum $`p`$ in the $`K\pi `$ system and the total CM energy $`E`$, i.e.,

$$E=\sqrt{p^2+m_\pi ^2}+\sqrt{p^2+m_K^2}\text{or}p=\frac{E}{2}\left[\left\{1\left(\frac{m_\pi +m_K}{E}\right)^2\right\}\left\{1\left(\frac{m_\pi m_K}{E}\right)^2\right\}\right]^{1/2}.$$
(16)
Correspondingly, for the reduced mass of the $`K\pi `$ system we define

$$\mu (E)=\frac{1}{2}\frac{dp^2}{dE}=\frac{E}{4}\left[1\left(\frac{m_K^2m_\pi ^2}{E^2}\right)^2\right].$$
(17)
In Fig. 1 we show the result of formula (14) for the cross sections of $`I=\frac{1}{2}`$ elastic $`P`$-wave scattering, for the parameter values $`\lambda =0.75`$ GeV<sup>-3/2</sup> and $`a=5`$ GeV<sup>-1</sup>, and with the substitution

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{B_{n,0}}{EE_{n,0}}\frac{0.5}{E0.945}\mathrm{\hspace{0.33em}1}\text{GeV}^{\mathrm{\hspace{0.17em}2}},$$
(18)
where we neglect possible $`\mathrm{}_c=2`$ contributions.

The theoretical curve agrees well with the data. For the $`P`$-wave $`K\pi `$ scattering length, we find here the somewhat too low result $`a_1^{1/2}=0.0085m_\pi ^1`$, to be compared to the experimental values taken from Ref.  (in units of $`m_\pi ^1`$), namely $`0.017`$ , $`0.018\pm 0.002`$ , and $`0.018`$ , or to the chiral-perturbation-theory result $`0.013\pm 0.003`$ . Nevertheless, the procedure of substitution (18) leads to a more than satisfactory description in the relevant domain of CM energy, thus allowing to read off the value for the bare $`K^{}`$(892) mass, i.e., 0.945 GeV.

Moreover, one may inspect the scattering matrix, which follows from expression (14) after substitution (18), for its singularities in the complex-energy plane. One finds one pole at

$$0.887\mathrm{\hspace{0.33em}0.027}i\text{GeV}.$$
(19)
The relation between the position of the singularity and the Breit-Wigner-like parameters is lost. However, one gains a simple relation to the confinement spectrum. Moreover, in the substitution (18) one may take an arbitrary number of resonances into account.

In order to verify that singularity (19) stems from the bare state at 945 MeV, we may stepwise switch off the model parameter $`\lambda `$ and inspect the theoretical positions of the corresponding singularities. This procedure is shown in Fig. (2).

It demonstrates beyond any doubt the relation between the singularity (19) and the bare state at 945 MeV. Notice that the motion of the singularity for small values of $`\lambda `$ is perturbative and quadratic in $`\lambda `$, as indicated by expression (15). However, for larger values of the coupling constant the singularity positions become more and more nonperturbative.

### 3.2 $`𝑲𝝅`$ $`𝑺`$-wave scattering

In Fig. 3 we show the result of formula (14) for the phase shifts of $`I=\frac{1}{2}`$ elastic $`S`$-wave scattering, for the values $`\lambda =0.75`$ GeV<sup>-3/2</sup> and $`a=3.2`$ GeV<sup>-1</sup>, and with the substitution

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{B_{n,1}}{EE_{n,1}}\frac{1.0}{E1.31}+\frac{0.2}{E1.69}\mathrm{\hspace{0.33em}1}\text{GeV}^{\mathrm{\hspace{0.25em}2}}.$$
(20)
Also for $`S`$-wave scattering we find good agreement between the theoretical curve and the data. The difference between $`P`$ and $`S`$ wave for the radius $`a`$ of the transition potential agrees well with the corresponding difference for the full pseudoscalar-pseudoscalar transition potentials depicted in Figs. 7 and 8 of Ref. , respectively. For the $`S`$-wave $`K\pi `$ scattering length we find here $`a_0^{1/2}=0.22m_\pi ^1`$, to be compared to the experimental values taken from Ref.  (in units of $`m_\pi ^1`$), i.e., $`0.33\pm 0.01`$ , $`0.237`$ , $`0.240\pm 0.002`$ , $`0.13\pm 0.09`$ , and $`0.22\pm 0.04`$ , or to the chiral-perturbation-theory result $`0.17\pm 0.02`$ .

From Eq. (20) one reads for the lowest $`J^P=0^+`$ isodoublet eigenstates of confinement the bare masses

$$E_{0,1}=\mathrm{\hspace{0.33em}1.31}\text{GeV}\text{and}E_{1,1}=\mathrm{\hspace{0.33em}1.69}\text{GeV}.$$
(21)
However, when we search formula (14), after substitution (20), for singularities, then we find, besides the two corresponding singularities at

$$(1.458\mathrm{\hspace{0.33em}0.118}i)\text{GeV}\text{and}(1.713\mathrm{\hspace{0.33em}0.019}i)\text{GeV},$$
(22)
also an additional one at

$$(0.714\mathrm{\hspace{0.33em}0.228}i)\text{GeV}.$$
(23)
In Fig. 4 we depict the transformation of the theoretical and experimental phase shifts of Fig. 3 into the $`I=\frac{1}{2}`$ elastic partial $`S`$-wave cross sections, which show a clear signal at about 830 MeV, with a width of some 500 MeV, as possibly also seen in a very recent experiment by the E791 collaboration . Notice, however, that neither the theoretical, nor the experimental cross sections exhibit a dip, characteristic for Breit-Wigner distributions.

The singularity (23), which in model  comes out at $`(0.727\mathrm{\hspace{0.33em}0.263}i)\text{GeV}`$, is interpreted as the isodoublet partner $`K_0^{}(727)`$ of the $`f_0(4001200)`$. A corresponding resonance has also been reported in several other works  in more recent years. The preliminary experimental result of the E791 collaboration reported in Ref.  awaits further confirmation.

The singularities (22) correspond to the ground state ($`n=0`$, $`\mathrm{}_c=1`$) at $`E_{0,1}`$ and the first radially excited state ($`n=1`$, $`\mathrm{}_c=1`$) at $`E_{1,1}`$ of the confinement spectrum, respectively. By reducing the value of the coupling constant $`\lambda `$, the singularities (22) move towards $`E_{0,1}`$ and $`E_{1,1}`$, respectively, which can be most clearly understood from expression (15). The complex-energy singularities of the partial $`S`$-wave scattering matrix (2.1), resulting from stepwise reducing the theoretical coupling constant $`\lambda `$ in formula (14), are depicted in Fig. (5).

It clearly demonstrates the relation between the singularities (22) and the bare states (21). Nonperturbative effects for larger values of the coupling constant can in particular be observed for the lower of the two resonances. First- (or second-) order perturbative calculations would result in completely different positions for the singularities corresponding to the model value of the coupling constant. Especially the real part of the mass shift is strongly affected by higher-order corrections. In a recent $`K`$-matrix analysis , as well as in a covariant quarkonium model , it is indeed also found that the bare state might be appreciably lower than the central resonance position for the $`K_0^{}(1430)`$.

The $`K_0^{}(727)`$ singularity (23) does not have a direct relation to any of the states stemming from the confinement mechanism. This fact is also most clearly demonstrated by the theoretical process of decoupling the $`\pi K`$ sector from the strange-nonstrange quarkonium sector, as depicted in Fig. (6).

The singularity (23) acquires a larger imaginary part when $`\lambda `$ is reduced, thus describing a state with ever increasing width, which is a highly nonperturbative phenomenon. For vanishing coupling constant, one observes that the $`K_0^{}(727)`$ disappears into the scattering background. This is a very important observation, since it implies that the $`K_0^{}(727)`$ singularity is a consequence of the transition mechanism that provides the communication between the $`\pi K`$ sector and the strange-nonstrange quarkonium sector. No such phenomena are observed for $`P`$\- and higher-wave meson-meson scattering. So we may conclude that a corresponding effect is screened from observation by the centrifugal barrier in Eq. (2). Consequently, the mechanism of valence-quark-pair creation appears to be more important in $`S`$-wave meson-meson scattering.

The absence of a direct relation between the $`K_0^{}(727)`$ singularity and the nonstrange-strange quarkonium sector has for the first time been demonstrated in Ref. , where similar singularities are shown to describe the properties of the $`f_0(4001200)`$, $`f_0(980)`$, and $`a_0(980)`$ resonances, thereby resolving two important issues: the nature of the light scalar mesons and the completion of the light scalar nonet.

The positions of the various singularities for $`\pi K`$ $`S`$-wave scattering in the complex-energy plane as a function of $`\lambda `$ can of course be obtained from the analyticity properties of the scattering matrix $`S(E)`$ (2.1), or, as no approximations are made in the determination of $`S(E)`$, directly from the dynamical equations (1) and (7). In an extensive study on multichannel scattering with permanently closed channels , resonances like the here described $`K_0^{}(727)`$ are distinguished from the fundamental resonances related to the bound states of $`H_c`$, and referred to as hadron molecular states. The idea of meson molecules has been worked out in several papers , as a possible explanation for the light scalar mesons $`a_0`$(980) and $`f_0`$(980). However, the term molecule gives the wrong impression that the $`q\overline{q}`$ component has no importance for these states. From the dynamical equations (1) one may determine the contribution of the strange-nonstrange quarkonium sector to the wave functions for any of the states under the $`K_0^{}(727)`$ resonance.

The theoretical distribution of Fig. 4 could well be used for experimental analysis, by optimising the adjustment of the here proposed parameters to the data. Especially the absence of the dip in the cross section is well accounted for in formula (14) with substitution (20). It has moreover the advantage that a direct relation exists between the theoretical distribution and the phase shifts and scattering matrices, which enables the precise location of the singularity associated with the $`K_0^{}(727)`$.

## 4 Conclusions

We have shown that the RSE parameters $`E_{n\mathrm{}_c}`$ of formula (12) relate experiment to qlQCD calculations of hadron masses better than the usual central Breit-Wigner positions of resonances do. This is not only a consequence of the potentially large mass shifts due to hadronic decay, but also owing to the nice feature of the RSE procedure that $`S`$-matrix singularities not originating from genuine confinement can easily be isolated. As an important application of the RSE, it is shown here, through the example of the $`K_0^{}(727)`$, how the light scalar mesons can be described by $`S`$-matrix singularities which are not directly related to the ground states of the $`J^P=0^+`$ confinement spectrum.

The RSE parameters $`B_{n\mathrm{}_c}`$ of formula (12) incorporate the unknowns of hadronic decay processes. Empirical knowledge of these parameters will certainly give a substantial contribution to the detailed study of hadronic decay at low and intermediate energies.

With respect to the quantitative conclusions of the present work, a word of caution is in place. Since expression (2.1) allows for it, we may also inspect the effects of higher thresholds. This has been carried out for harmonic-oscillator confinement in the model of Refs. , with the following results: for $`P`$-wave non-exotic two-meson scattering, both with $`J^{PC}=0^+`$ and $`J^{PC}=1^{}`$, the real parts of the singularities which correspond to the ground states at $`E_{0,0}`$ of the confinement spectrum come out some 300–400 MeV below $`E_{0,0}`$. For the higher radial excitations these shifts are considerably smaller. For $`S`$-waves the shifts are also smaller for the ground states and, moreover, in the positive direction. When we compare these findings with the one-threshold results shown above, we must conclude that higher thresholds should be taken into account for a more quantitative determination of the confinement and decay mechanisms. In particular, the second singularity of formula (22) might come out rather displaced, if higher thresholds are accounted for. This will be investigated in future work.

Acknowledgements. We are indebted to Jeffrey Appel for drawing our attention to the need for a more phenomenological approach to data analysis in coupled-channel models for confinement and hadronic decay. We also wish to thank Alexander A. Osipov, Brigitte Hiller, Alex Blin and Hans Walliser for useful discussions and suggestions. This work is partly supported by the Fundação para a Ciência e a Tecnologia of the Ministério da Ciência e da Tecnologia of Portugal, under contract numbers PESO/P/PRO/15127/99, POCTI/35304/FIS/2000 and CERN/P/FIS/40119/2000.


text:
# Near action degeneracy of periodic orbits in systems with non-conventional time reversal

## I Introduction

The form factor $`K(\tau )`$ is defined as the Fourier transformation of the two point energy correlation function of the quantum system where $`\tau `$ is time measured in units of the Heisenberg time $`T_H`$. If the Gutzwiller trace formula is used for the density of energy levels the form factor becomes a double sum over the classical periodic orbits $`\gamma `$ with the period $`T=\tau T_H`$:

$$K(\tau )=\underset{\mathrm{}0}{lim}\frac{1}{T_H}\underset{\gamma ,\gamma ^{}}{}A_\gamma A_\gamma ^{}^{}\mathrm{e}^{i(S_\gamma S_\gamma ^{})/\mathrm{}}\delta \left(T\frac{T_\gamma +T_\gamma ^{}}{2}\right).$$
(1)
Here, $`A_\gamma `$ are the stability coefficients of the orbits; the angle brackets signify average over a small interval of time .

Only pairs of orbits whose action difference is not large compared with $`\mathrm{}`$ can make a contribution surviving the time averaging. The diagonal approximation  takes into account the diagonal terms $`\gamma =\gamma ^{}`$ and pairs of mutually time reversed orbits in the case that the dynamics of the system is invariant with respect to time reversal (TR). The result of the diagonal approximation is then

$$K_{\mathrm{diag}}(\tau )=2\tau .$$
(2)
which is to be compared with the form factor of the Gaussian ensemble of random orthogonal matrices for $`0<\tau <1`$,

$$K^{GOE}(\tau )=2\tau \tau \mathrm{log}(1+2\tau ).$$
(3)
Obviously, the diagonal approximation reproduces only the first term of the Taylor expansion of the random matrix form factor. It may be expected that taking into account less obvious pairs of periodic orbits with small action difference higher terms of the random matrix form factor (3) will be recovered, in line with the conjecture stated by Bohigas, Giannoni and Schmidt . In the theory of disordered systems the so-called weak diagonal approximation was used by Smith et al  who stressed the importance of pairs of trajectories with the multiple-loop topology. The pairs of orbits whose contribution seems to be responsible for the higher order terms of the form factor in the case of clean chaos have been discovered only recently by Sieber and Richter .

The Sieber-Richter pair of orbits is schematically shown in Fig.1. One of its members contains a self-intersection with a small crossing angle $`ϵ`$ and consists of two loops, one of which is passed clockwise and another one counterclockwise. Playing with small deformations of such an orbit it can be demonstrated that there exists a partner periodic orbit which is almost everywhere exponentially close to the original one; however, at one place in configuration space the partner orbit has an avoided crossing at the place of a self-intersection.

The partner can be regarded as the same two-loop construction as the orbit with crossing, however, now both loops are passed in the same sense of traversal. The time reversal symmetry is essential for the existence of the partner; otherwise retracing the motion along one of the loops would not be allowed. The difference of actions of the orbit and its partner is

$$\mathrm{\Delta }S\frac{p^2ϵ^2}{2m\lambda }$$
(4)
where $`\lambda `$ is the Lyapunov constant of the orbit. For small $`ϵ`$ it can be of the order of $`\mathrm{}`$.

The stumbling block is to evaluate the number of the Sieber-Richter pairs and to sum over their contributions. Up to now it has been done only for the billiards in the space of constant negative curvature where all orbits have the same Lyapunov constant and the Maslov index equals zero. Contribution of all pairs Fig.1 have been found to be $`2\tau ^2`$  which coincides with the term of the order $`\tau ^2`$ in the expansion of the random matrix form factor (3).

One of the open problems in the Sieber-Richter theory is connected with the so-called non-conventional time reversal (NCTR) symmetry. It is often encountered in systems in external magnetic field provided that there exists a suitable spatial symmetry . Their Hamiltonian is invariant under the conventional time reversal combined with an appropriate spatial symmetry operation. Each periodic orbit of a NCTR-symmetric system has a twin with the same action. However, unlike the conventional TR dynamics, the twin is not the same orbit retraced backwards.

The statistical properties of the energy levels of the systems with NCTR and the conventional TR symmetry are identical. In particular, the form factor of the systems with NCTR is given by (3). Therefore, the diagonal approximation is not sufficient, and pairs of orbits with very close action must also exist and contribute to $`K(\tau )`$.

At the first glance it seems that the Sieber-Richter arguments collapse in the case of NCTR. Even a very weak magnetic field destroys the closeness of action in the pair because of the totally different magnetic flux (due to the fact that one of the loops is passed in the opposite sense by the members of the pair). In a stronger field whose impact on the trajectory cannot be neglected the pairs in Fig. 1 simply cannot exist because passing the same loop in the opposite direction would contradict the equations of motion. The spatial symmetry implied by NCTR does not help. Since the correct semiclassical explanation of the form factor (3) must be essentially the same, be it systems with TR or NCTR, the inability to produce the contributing pairs of orbits in the NCTR-symmetric systems might compromise the whole Sieber-Richter theory. Below we show that these apprehensions are groundless, and the case of NCTR does not present any new difficulties.

Consider a two dimensional motion in the orthogonal uniform magnetic field $`𝐁=B𝐞_z`$ and suppose that the potential energy has the symmetry

$$V(x,y)=V(x,y).$$
(5)
We shall use the gauge $`A_x=By,A_y=A_z=0`$. Then the classical Hamiltonian of the system

$$H=\frac{1}{2m}\left(p_x+\frac{eBy}{c}\right)^2+\frac{p_y^2}{2m}+V(x,y)$$
(6)
will be invariant with respect to NCTR consisting of the conventional time reversal TR (changing the sign of the canonical momenta) followed by the reflection in the $`x`$ axis of the plane (replacement $`yy,p_yp_y`$). An arbitrary periodic orbit of our system with the trajectory

$$x=x(t),y=y(t),$$
(7)
and the momentum $`p_y(t)=m\dot{y}(t)`$ along $`y`$ has a NCTR twin,

$$\stackrel{~}{x}(t)=x(t),\stackrel{~}{y}(t)=y(t),\stackrel{~}{p}_y(t)=p_y(t).$$
(8)
The twin has the same magnetic flux and action. Its trajectory is obtained from the original one by reflection in the $`x`$ axis while the sense of traversal on both orbits is the same. Let us canonically transform the variables as

$$x^{}=x,y^{}=p_y,p_x^{}=p_x,p_y^{}=y$$
(9)
and consider how our two orbits project on the new configuration space $`x^{}y^{}`$. The original orbit will be described by the equations

$$x^{}(t)=x(t),y^{}(t)=p_y(t)$$
(10)
with $`x,p_y`$ given by Eqs. (7). Its NCTR twin obeys

$$\stackrel{~}{x}^{}(t)=x(t),\stackrel{~}{y}^{}(t)=p_y(t).$$
(11)
The only difference between (10) and (11) is the change of sign of $`t`$. Consequently, in the $`xp_y`$ plane the two orbits are depicted by the same closed curve traversed in opposite directions: NCTR acts on a periodic orbit in the coordinates $`xp_y`$ exactly like the usual time reversal in the ordinary configuration space.

Now it is easy to see that in the case of NCTR the Sieber-Richter arguments for the existence of the two-loop pairs of orbits Fig.1 remain fully valid. However, these pairs may only exist in the $`xp_y`$ plane (Fig.2a) which is the configuration plane in the new coordinates. Only in this projection of the phase space to a two-dimensional submanifold switching the sense of traversal of a loop caused by the replacement of crossing by an avoided crossing is compatible with the equations of motion. Therefore, only in the $`xp_y`$ plane pairs as depicted in Fig.1 exist and have close actions in strong magnetic fields. Evaluation of the number of Sieber-Richter pairs for the system with NCTR and their contribution to the form factor therefore does not lead to new difficulties as compared to the problem of conventional TR-invariant systems.

It is instructive to see how the pairs Fig. 2a look like when the respective motion is projected onto the usual configuration space $`xy`$ (Fig.2b). They have little in common with the double-loop Sieber-Richter pairs of Fig.1. However, the general idea of building a new orbit with practically the same action by gently reconnecting parts of the original orbit and its NCTR twin is still obvious. Depending on the projection chosen, the criterion for finding Sieber-Richter pairs changes. In systems with NCTR one may either look for the two loop orbits with small opening angle in the $`xp_y`$ plane or search for the orbits like in Figure 2b. The latter may be preferable for systems like billiards in the magnetic field whose trajectory in the $`xp_y`$ plane is discontinuous.

It may be somewhat puzzling that a jump in representation is needed for the Sieber-Richter treatment when the magnetic field is switched on: instead of $`xy`$ space we must shift to $`xp_y`$. However such a jump is only natural in view of the change of the universality class of the dynamics. A chaotic system with the NCTR symmetry belongs to the Gaussian orthogonal ensemble only in the presence of the magnetic field. When the field is switched off the spectrum splits into two independent subspectra (even and odd with respect to $`yy`$). Superposition of two such spectra whose levels may cross each other creates a specific ensemble obviously different from GOE, usually called GOE$`\times `$GOE.

The systems with the plane of symmetry in the uniform magnetic field constitute the most important but not the only example of NCTR. Consider e.g. a two-dimensional system with a center of symmetry

$$V(x,y)=V(x,y)$$
(12)
in an extremely non-uniform magnetic field

$$𝐁=y𝐞_z.$$
(13)
In the gauge $`A_x=By^2/2,A_y=A_z=0`$ the Hamiltonian of the system is invariant with respect to the NCTR composed of the time and spatial inversion which means that the momenta are unchanged. The partner obtained from a periodic orbit by this symmetry operation coincides with the original orbit if we draw it in the plane of momenta $`p_xp_y`$, hence it is in this plane that the Sieber-Richter pairs are described by small intersection angle $`ϵ`$.

To summarize, the Sieber-Richter double loop pairs Fig. 1 or Fig. 2a may be observed only in the configuration space whose generalized coordinates are unchanged by the particular NCTR. A periodic orbit and its NCTR twin are depicted in this space by the same curve traversed in opposite direction. This may be the usual coordinate space $`xy`$ in systems with conventional time reversal symmetry, the $`xp_y`$ space in systems with the NCTR $`tt,yy`$, and the $`p_xp_y`$ space if NCTR is described by $`tt,xx,yy`$. An attempt to break out of this symmetry-dictated space by a canonical transformation mixing the coordinates and momenta will immediately strip the Sieber-Richter pairs of their double-loop, intersection/avoided-crossing appearance. However, using appropriate citeria to describe the close action partners can be recognized in principle in any two dimensional projection of the phase space.

Appendix

Calculations in  are based on the geometric evaluation of the shifts of momenta using the fact that velocities and momenta in the conventional coordinate space are practically identical. After our canonical transform the connection between the new “velocities” and the momenta becomes more complicated. In particular the direction of the momentum is no longer tangent to the trajectory in the new configuration space. That means that the Sieber-Richter result has to be rederived.

We assume the existence of a self-crossing periodic orbit with a small opening angle like the one shown in Fig. 1, 2a. A coordinate frame is introduced with its origin at the point of crossing, and the $`x^{}`$ axis along the bisector of the small angle. Consider the Poincaré section at $`x^{}=0`$ with the coordinates $`y^{}`$ and $`p_y^{}`$ on the crossing plane. The true Poincaré map (Fig.3) is obtained when passages of the $`x^{}=0`$ plane with a certain sign of $`\dot{x^{}}`$, say, $`\dot{x^{}}>0`$ are marked. The self-crossing orbit will then be depicted by a periodic point $`O`$ on the $`p_y^{}`$ axis with $`p_y^{}`$ positive and small. The self-crossing TR (NCTR) twin of this orbit will produce another periodic point $`O^{}`$ symmetrical with respect to the $`y^{}`$ axis: $`y^{}=0,p_y^{}<0`$.

We shall concentrate, however, on the submaps $`R`$, $`L`$ of the Poincaré map describing the transform of $`y^{},p_y^{}`$ generated by the right and left loop of the orbit. We shall mark the crossing of the $`x^{}=0`$ plane both for $`\dot{x^{}}>0`$ and $`\dot{x^{}}<0`$; such a break of the rules is needed since $`R`$ and $`L^{}`$ are not true Poincaré maps. We shall also be interested in the TR (NCTR) submaps obtained by passing the loops of the orbit in the direction opposite to Fig. 3; these will be denoted $`R^{}`$ and $`L`$, respectively. The periodic point $`O`$ of the total Poincaré map is simultaneously the periodic point of the submaps $`R`$ and $`L`$ whereas $`O^{}`$ is the periodic point of the submaps $`L^{}`$ and $`R^{}`$.

Each submap can be characterized by its stability matrix connecting the initial and final deflections of the $`y^{}`$ coordinate and momentum from the periodic point of the respective submap. We shall need $`M_R`$ (the stability matrix of the right loop passed as it is shown in Fig. 3) and $`M_L^{}`$ (the one for the left loop followed in the direction opposite to Fig.3). The large eigenvalues $`\mathrm{\Lambda }_R`$ and $`\mathrm{\Lambda }_L^{}`$ of these two matrices can be evaluated as $`\mathrm{exp}(\lambda T)`$ where $`\lambda `$ is the Lyapunov constant and $`T=T_R,T_L^{}`$ is the period of the respective loop. The periods of the orbits in the sum for the form factor (1) are of the order $`\mathrm{}^1\mathrm{}`$, therefore the larger eigenvalues of the stability matrices are exponentially large whereas the smaller ones ($`1/\mathrm{\Lambda }`$) are exponentially close to zero. The respective eigenvectors determine the unstable and stable directions of each of the submaps.

Consider Fig.4 where periodic points $`O,O^{}`$ of the submaps $`R,L^{}`$ and their stable and unstable directions are shown. Let us investigate the application of the submap $`R`$ to an initial point $`P`$ chosen in the vicinity of the crossing $`A`$ of the stable direction of $`R`$ and unstable direction of $`L^{}`$. We shall represent the initial radius-vector by an expansion in powers of $`\mathrm{\Lambda }_R^1`$,

$$𝐫_P=𝐞_s^R\left(l_{OA}+\frac{c_1}{\mathrm{\Lambda }_R}+\mathrm{}\right)+𝐞_u^R\left(\frac{l_{OB}}{\mathrm{\Lambda }_R}+\frac{d_2}{\mathrm{\Lambda }_R^2}+\mathrm{}\right).$$
(14)
Here $`𝐞_s^R,𝐞_u^R`$ are the eigenvectors of the stability matrix of $`R`$ along the stable and unstable directions, $`l_{OA}`$ and $`l_{OB}`$ are distances from $`A`$ and $`B`$ to the periodic point $`O`$ ($`B`$ is the crossing of the unstable direction of $`R`$ and stable direction of $`L^{}`$); the coefficients $`c_1,d_2,\mathrm{}`$ are so far undetermined.

After the loop $`R`$ has been completed the point $`P`$ will be mapped to the point $`Q`$ obtained by squeezing along the stable and stretching along the unstable direction with the coefficient $`\mathrm{\Lambda }_R`$:

$$𝐫_Q=𝐞_s^R\left(\frac{l_{OA}}{\mathrm{\Lambda }_R}+\frac{c_1}{\mathrm{\Lambda }_R^2}+\mathrm{}\right)+𝐞_u^R\left(l_{OB}+\frac{d_2}{\mathrm{\Lambda }_R}+\mathrm{}\right).$$
(15)
It is seen that $`Q`$ is infinitely close to the crossing point $`B`$ of the unstable direction of $`R`$ and the stable direction of $`L^{}`$; thus the distance of $`Q`$ from the stable manifold of $`L^{}`$ is exponentially small and depends on the coefficients $`c_1,d_2,\mathrm{}`$.

Now consider the loop $`L^{}`$ taking $`Q`$ as its initial point. The loop will practically annihilate the stable component $`𝐞_u^Rl_{OB}`$ and place the final point somewhere on its unstable manifold. The exact position of the final point on the unstable manifold of $`L^{}`$ depends on $`c_1,d_2,\mathrm{}`$; these can be finetuned so that the final point will coincide with the initial point $`P`$ of the loop $`R`$. ***More accurately, demanding that the final point of the second loop coincides with the initial point $`P`$ of the first loop we obtain a set of equations for consecutive definition of $`c_1,d_2,c_2,d_3\mathrm{}`$. The coefficients in the expansion will not be growing and convergence for $`𝐫_P`$ guaranteed provided $`\mathrm{\Lambda }_R<\mathrm{\Lambda }_L`$ which can be assumed wihout loss of generality. But that would mean that $`P`$ is a periodic point of the composition of the submaps $`R`$ and $`L^{}`$. It corresponds to a new periodic orbit composed of the deformed loops $`R`$ and $`L^{}`$. The new orbit crosses the “true” Poincaré map at the point infinitely close to the crossing of the stable manifold of $`R`$ and unstable manifold of $`L^{}`$. The TR (NCTR) twin of the new periodic orbit can be found by considering the sequence of the submaps $`R^{}`$ and $`L`$. The new orbit and its twin are of course the Sieber-Richter partners with avoided crossing of the orbit with self-intersection.

As can be seen we have not used the connection between the velocities and momenta of the usual configuration space. Note also that the effect of the pair formation seems to remain in force even when the nonlinear corrections to the submaps $`R,L^{}`$ are taken into account, with their stable and unstable manifolds depicted by curves rather than straight lines, as long as the stable and unstable manifold intersect only once. For small angle $`ϵ`$, this condition should be fulfilled.


text:
# Rapidity particle spectra in sudden hadronization of QGP

\[

## Abstract

We show that the remaining internal longitudinal flow of colliding quarks in nuclei offers a natural explanation for the diversity of rapidity spectral shapes observed in Pb–Pb 158$`A`$ GeV nuclear collisions. Thus QGP sudden hadronization reaction picture is a suitable approach to explain the rapidity spectra of hadrons produced.

PACS: 12.38.Mh, 25.75.-q

\]

If a quark-gluon plasma deconfined state super cools and decomposes in a sudden fashion into free streaming hadrons , we should be able to see in the hadron spectra certain features of quark statistical distributions. Such a QGP fingerprint on $`m_{}`$-spectra of strange baryons and antibaryons which are seen to be identical is one of key indicators for the sudden hadronization picture .

A question which has remained open is how sudden hadronization model can be reconciled with the diverse shapes of rapidity spectra of observed hadrons. Great differences are known to be present: while hadrons made of quarks brought into the collision, such as, e.g., protons, show a pronounced double hump structure, progressively more peaked distributions are seen for particles which are more and more composed of newly made quarks .

We present a schematic description of this phenomenon based on the assumption that the quark momentum distribution in QGP is imprinted on hadrons produced, and not erased in a post hadronization period. The model we explore is a direct continuation of the study of $`m_{}`$ shape and abundance of hadrons produced .

We note that many light quarks are brought into the collision fireball by colliding nuclei — in 158$`A`$ GeV Pb–Pb central interactions there are 1000 + such valance quark participants and there is a similar total number of hadrons produced. In general the two incoming quark fluids can retain memory of the early stage in the collision: in the center of momentum frame two quark fluids collide with initial rapidities $`y_i=\pm 2.92`$. In order to reproduce the proton data, we assume in what follows, that the average remaining longitudinal quark flow is $`y_f=\pm 1`$. This average includes quarks from $`q\overline{q}`$ pairs produced, which are mostly found at central rapidity. We further assume that the average flow of $`\overline{s},s,\overline{q}`$ quarks is zero.

If quarks were not confined, their emission spectrum would be computed in a thermal emission model with flow. We proceed to establish this spectrum. Each emitting volume element is assumed to be in motion, with a four-vector $`u^\mu `$, and in the local frame at rest the temperature in the volume element is $`T`$, and the local quark distribution is of Boltzmann type. In order to count the states, we apply Touschek invariant phase space measure for the level density :

$$\frac{V_0d^3p}{2\pi ^3}e^{E/T}\frac{V_\mu p^\mu }{(2\pi )^3}d^4p\mathrm{\hspace{0.17em}2}\delta _0(p^2m^2)e^{p_\mu u^\mu /T}.$$
(1)
Here, $`V_\mu =V_0u_\mu `$ is the flowing volume element as observed in the laboratory frame. $`V_0`$ is the comoving volume element in the local rest frame. $`\delta _0`$ is the Dirac-delta function for the positive roots only. This is an invariant function for all proper Lorentz transformations. Using the rapidity and transverse momentum $`p_{}`$ (or mass $`m_{}=\sqrt{m^2+p_{}^2}`$) variables of a particle in a cylindrical symmetric decomposition with respect to the collision axis

$$p^\mu =(m_{}\mathrm{cosh}y,p_{}\mathrm{cos}\varphi _p,p_{}\mathrm{sin}\varphi _p,m_{}\mathrm{sinh}y),$$
(2)
we recover the well known invariant measure in Eq. (1):

$`2\delta _0(p^2m^2)d^4p={\displaystyle \frac{d^3p}{E}}`$ $`=`$ $`m_{}dm_{}dyd\varphi _p`$ (3)
$`=`$ $`p_{}dp_{}dyd\varphi _p.`$ (4)
The four-velocity $`u^\mu `$ allowing for both transverse and longitudinal flows is:

$`u^\mu `$ $`=`$ $`(u_0,u_x,u_y,u_z),u^2=1;`$ (5)
$`u_0`$ $`=`$ $`\mathrm{cosh}y_{}\mathrm{cosh}y_{};u_x=\mathrm{sinh}y_{}\mathrm{cos}\varphi _u;`$ (6)
$`u_y`$ $`=`$ $`\mathrm{sinh}y_{}\mathrm{sin}\varphi _u;u_z=\mathrm{cosh}y_{}\mathrm{sinh}y_{}.`$ (7)
We recall the usual relations:

$$\frac{m_{}}{m}=\gamma _{}=\mathrm{cosh}y_{},\frac{p_{}}{m}=v_{}\gamma _{}=\mathrm{sinh}y_{},$$

which illustrate the structural identity of Eq. (2) with Eq. (5).

It is straightforward now to obtain the quantity required to construct spectra Eq. (1) ,

$$u_\mu p^\mu =\gamma _{}\left[m_{}\mathrm{cosh}(yy_{})p_{}v_{}\mathrm{cos}\varphi \right],$$
(8)
where $`\varphi =\varphi _p\varphi _u`$. For a suitable choice of coordinate system in which the $`x`$-axis is pointing in the direction of transverse flow vector, $`\varphi _u=0`$, and we can use the particle emission angle as the azimuthal angle of integration for radiation emitted from the volume element: $`\varphi =\varphi _p`$. This has the advantage that we can easily restrict the counting of particles to those emitted outwards including $`0<\varphi _p\pi /2`$ instead of the full range $`0<\varphi _p\pi `$

In Eq. (8), the emitted particle variables are as usual $`y,m_{}`$ and the flow vector variables used are $`y_{}`$ and $`v_{}`$ with $`\gamma _{}=1/\sqrt{1v_{}^2}`$. The invariant spectrum written explicitly is:

$`{\displaystyle \frac{d^2N}{m_{}dm_{}dy}}`$ $`=`$ $`{\displaystyle \frac{d\varphi \gamma _{}}{(2\pi )^3}\left[m_{}\mathrm{cosh}(yy_{})p_{}v_{}\mathrm{cos}\varphi \right]}`$ (10)
$`\times e^{\gamma _{}\left[m_{}\mathrm{cosh}(yy_{})p_{}v_{}\mathrm{cos}\varphi \right]/T}.`$
There are some differences in detail between this traditional expression  and recent proposals based on consideration of the freeze-out surface geometric properties . Certain features we find are the same in both approaches: the rapidity $`y`$ distribution is not appreciably influenced by transverse dynamics (i.e., transverse velocity and temperature), which effectively separates. The $`m_{}`$ spectra are only dependent on $`T,v_{}`$, but practically not on the longitudinal flow $`y_{}`$. We note that Eq. (10) retains absolute number of particles independent of the flow parameters, if integral over $`\varphi `$ includes the full range $`0\varphi \pi `$.

The rapidity spectra of massless quanta in the deconfined phase are shown in figure 1. We see how the thermal spectrum ($`m=0`$, $`T=145`$ MeV, $`v_{}=0.52`$) gradually ‘flows’ apart as $`y_{}`$ is increased from $`y_{}=0`$ (dotted) to $`y_{}=1.5`$ (solid) in step of 0.5 . Included in the double numerical integral (over $`\varphi `$ and $`m_{}`$) are particles emitted outwards with reference to the direction of the flow vector of the volume element. Model studies of kinetic dynamics show that inside the QGP matter we find a mix of fluids, the incoming quarks are still retaining some of the original flow while all newly made quark pairs have no memory of the initial condition of colliding matter. In particular, strange quark pairs made in the plasma do not flow.

Final state hadrons are born from deconfined quarks with such spectra. In the previously studied case of $`m_{}`$ spectra  it was found that the transverse flow parameters of the hadronizing quark matter were in agreement with results of chemical particle production parameters. Thus the shape of $`m_{}`$ spectra of hyperons, antihyperons and kaons was as if these particles were directly made of hadronizing QGP quanta. Addressing here the rapidity spectra we will take this point of view but proceed more carefully. As hyperon are formed at the fireball breakup, any remaining longitudinal flow present among fireball constituents will be imposed on the product particle, thus $`\mathrm{\Lambda }`$-spectra containing potentially two original valence quarks are stretched in $`y`$, which $`\overline{\mathrm{\Lambda }}`$-$`y`$-spectra are not, as they are made from newly formed particles. One would thus expect that anti-hyperons comprising only newly made quarks will appear without flow, while $`\mathrm{\Lambda }`$-spectra will show flow, if such was still present.

We thus implement a particle-staged turn-on of longitudinal flow at SPS where in central Pb–Pb collisions the incoming valance quark number and total number of hadrons produced is about the same. For protons we assume that effectively all quarks are remembering some of the original flow. For $`\mathrm{\Lambda }`$ with one strange quark we take a mix of 2/3 weight in the spectrum with flow and 1/3 without. For particles like $`K^+(u\overline{s})`$ we take a 50-50% mix and for all newly made particles like $`\overline{\mathrm{\Lambda }},\mathrm{\Phi }`$ we assume that only no flow components contribute.

To describe proton rapidity spectra we assume after inspection of  that the longitudinal flow is $`y_{}\pm 1`$, and choose this value without a further attempt to fit the spectra. With $`m_q=0`$, $`T=145`$ MeV and $`v_{}=0.52c`$ (both parameters could be chosen very differently without altering the result, these were taken in view of results of ), we obtain the proton rapidity spectrum we present in figure 2, bottom panel. We note that we have to average over positive and negative flows $`y_{}=\pm 1`$, since the collision in the CM frame involves both.

The strange quark content of the $`\mathrm{\Lambda }`$ contributes with relative strength 33% in the central region, which results in a flat distribution — see the second panel from the bottom in figure 2. Next we take a 50-50% mix which gives us the result shown in the third panel from the bottom, corresponding to the $`K^+(u\overline{s})`$ rapidity distribution and when resonance decays are considered, also pion distribution. Finally, in the top panel we show prototype of the rapidity distribution arising for all completely newly made particles such as $`\overline{\mathrm{\Lambda }},\mathrm{\Phi }`$. Though far from perfect, we see in this schematic study practically exactly the behavior observed.

The same schematic rapidity behavior is seen in S–S collisions at 200$`A`$ GeV for the particles here considered. A direct comparison of the charged particle distributions between Pb–Pb and S–S collision systems seen in Ref. suggests that $`y_{}`$ is about 0.4 units of rapidity larger in the lighter collision system. Though seemingly a small change, this opens by 50% the gap between the fluids as we saw in figure 1 and a much more pronounced central rapidity reduction in certain particle abundances ensues.

Our schematic method can help understand the particle rapidity distributions at RHIC. The PHOBOS collaboration has presented a nearly flat distribution of charged hadrons with small peaks near to pseudo rapidity $`\eta =1.5`$ at the edge . Trying our model we finds that rather than to have a 3-fluid quark source, we need now a continuous distribution of sources residing in the range of rapidity $`1.6y1.6`$. In other words if the hadron rapidity distribution is as before offering an image of the quark flow, we have a flow plateau as expected in the hydrodynamical scaling solution proposed by Bjørken .

In summary of this brief contribution, we have shown, within a schematic dynamic model, that quark flow in QGP can explain the diversity of particle rapidity spectra observed at $`\sqrt{s_{NN}}=17.2`$ GeV, where clearly we do not see yet a high energy scaling behavior. What we have shown is that a three component fluid, two made of quarks still flowing against each other, and a third one made of newly made components (antiquarks, strange quarks), is capable to imprint on the rapidity spectra just the right systematic behavior. This transfer of quark dynamics to hadron dynamics is similar to what was found in the study of the $`m_{}`$ spectra. For hadrons to emerge with momentum spectra of quarks we need a hadronization with sudden breakup of the (super cooled) QGP plasma phase into free-streaming hadrons. Thus this study shows that such a picture is capable to explain rapidity spectra.

Acknowledgments: Work supported in part by a grant from the U.S. Department of Energy, DE-FG03-95ER40937. Laboratoire de Physique Théorique et Hautes Energies, University Paris 6 and 7, is supported by CNRS as Unité Mixte de Recherche, UMR7589.


text:
# INFLATON FIELD AND PRIMORDIAL BLACKHOLE

## 1 INTRODUCTION :

In astrophysics, it is known that the black holes are the ultimate corpse of a collapsing massive star when the mass of a given star exceeds about twice that of the solar mass. It is probable in cosmology that other kind of black holes having mass smaller than the solar mass may have produced. These are lighter so are hotter and important in the context of determination of Hawking radiation.1 The later types of black holes are popularly known as primordial black holes (henceforth, PBH). Such kind of blackholes are likely to be formed during quantum fluctuation of the matter fields in the early universe during inflation. The inflationary scenario of the early universe are derived out of an inflaton field guided by a potential during its slow rolling approximation. The condition that PBH may produce requires a cosmological constant was obtained by this slow roll approximation as the potential can be treated as a constant (flat potential) during the short inflationary regime.2 Although the chance of getting PBH’s are small today may have been more in the very early universe. They might have played a crucial role in the structure formation of the late universe.3 Recently, Meng et al.4 studied the creation of PBH from inflaton field starting from the static spherically symmetrical de Sitter black hole spacetimes with charge and mass. They examined the PBH properties, for both charged and uncharged cases and obtained a relationship between the primordial black hole pair formation and the energy transfer of the inflaton field.

In this work we consider a special form of Reissner-Nordstr$`\ddot{o}`$m-de Sitter model with metric ansatz5 given by

$$ds^2=A(r,t)dt^2+A(r,t)^1dr^2+r^2d\mathrm{\Omega }_2^2$$
(1)
with

$$A(r,t)=\left(1\frac{a}{r}\right)^2\left[1\frac{\mathrm{\Lambda }}{3}\left(r^2+2ar+3a^2\right)\right]$$

where $`a`$ is a constant and $`\mathrm{\Lambda }`$ is the cosmological term. The mass and charge of the corresponding black hole solution are given by

$$𝐌=a\left(1\frac{2}{3}\mathrm{\Lambda }a^2\right),$$
(2)
$$𝐐^2=a^2\left(1\mathrm{\Lambda }a^2\right).$$
(3)
Here, the cosmological parameter $`\mathrm{\Lambda }`$, is taken to be a function of the homogeneous scalar field $`\varphi (t)`$ which causes inflation. We choose for $`\mathrm{\Lambda }`$ the the form of potential used by Peebles and Vilenkin6 as

$$\mathrm{\Lambda }(\varphi )=b\left(\varphi ^4+\varphi _o^4\right).$$
(4)
It is to be noted that in the pre-inflationary era, $`\varphi `$ dominates the dynamics and hence $`\mathrm{\Lambda }`$ is considered as a variable quantity. But as the inflationary expansions proceeds, $`\varphi `$ gradually becomes insignificant and at present we may take $`\varphi `$ to be a constant.

Further, in order to solve the flatness problem during the inflation we must have

$$\frac{1}{2}\left(\frac{\mathrm{\Lambda }^{}}{\mathrm{\Lambda }}\right)^2=\frac{8\varphi ^6}{(\varphi _o^4+\varphi ^4)^2}<<1$$
(5)
and

$$\frac{\mathrm{\Lambda }^{\prime \prime }}{\mathrm{\Lambda }}=\frac{12\varphi ^2}{(\varphi _o^4+\varphi ^4)}<<1.$$
(6)
If we introduce, $`\rho =r/𝐌`$, $`z=\mathrm{\Lambda }𝐌^2`$, then we get from $`A(r,t)=0`$

$$z=\frac{3\left(1+\frac{2𝐐^2}{a^2}\right)^2}{\left[\rho ^2\left(1+\frac{2𝐐^2}{a^2}\right)^2+6\rho \left(1+\frac{2𝐐^2}{a^2}\right)+27\right]}.$$
(7)
We note that $`A(r,t)=0`$ has two non-degenarate and a degenerate roots so

$$\frac{dz}{d\rho }=0$$
(8)
gives us a unique root

$$\rho =3/\left(1+2𝐐^2/a^2\right)=\rho _{}$$
(9)
with $`z(\rho _{})=\frac{1}{6}\left(1+2𝐐^2/a^2\right)^2`$.

Thus we have the following two situations

(i) if $`z(\rho _{})=\mathrm{\Lambda }𝐌^2`$, then there is only one horizon

(ii) if $`z(\rho _{})>\mathrm{\Lambda }𝐌^2`$, then no horizon exists.

The first situation is known as the critical solution and can be written in explicit form

$$b\left(\varphi _o^4+\varphi ^4\right)𝐌^2=\frac{1}{6}\left[1+\frac{2𝐐^2}{9𝐌^2}\left(1+\frac{2𝐐^2}{a^2}\right)^2\right]^2$$
(10)
provided

$$8𝐐^2<9𝐌^2.$$
(11)
However for an uncharged case we have

$$b\left(\varphi _o^4+\varphi ^4\right)𝐌^2=\frac{1}{6}$$
(12)
To conclude we obtain the criterion for a primordial black hole formation using a variable cosmological constant as suggested be Peebles and Vilenkin as the potential for the inflaton field.

Acknowledgements :

The authors (BCP and SC) like to thank the Inter-University University Centre for Astronomy and Astrophysics ( IUCAA ), Pune for warm hospitality and providing facilities to complete this work during a visit.