text:
# Hole-doping effects on a frustrated spin ladder

## I Introduction

Quantum spin systems with geometrical frustration have attracted considerable attention recently. There are a number of typical materials, such as the spin-gap compound SrCu2(BO3)2 Kageyama  with the orthogonal-dimer structure, Shastry ; Miyahara ; Koga\_let  the pyrochlore compounds with the lattice of corner-sharing tetrahedral network of magnetic ions, Harris ; Canals ; Isoda ; Koga ; Tsunetsugu ; Ueda ; Fujimoto  etc.

Besides frustrated spin systems studied intensively so far, metallic systems with geometrical frustration have attracted much interest. For instance, it has been claimed that frustration due to a tetrahedral network of V ions may be important to understand heavy-fermion behavior of the metallic pyrochlore compound $`\mathrm{LiV}_2\mathrm{O}_4`$. LiV2O4  Also, another prototypical pyrochlore compound $`\mathrm{Y}(\mathrm{Sc})\mathrm{Mn}_2`$ Y(Sc)Mn2  shows spin-liquid behavior in a metallic phase.

Stimulated by these experimental findings, we study here the hole-doping effects on fully frustrated spin systems. To this end, we consider a simplified model in one dimension (1D), which still possesses the essence of strong frustration in the above-mentioned systems. To be precise, we employ the t-J ladder model with the orthogonal dimer structure, which also has linked-tetrahedra as a key structure. The advantage to employ the 1D model is that this simplification allow us to accurately calculate physical quantities with use of such a reliable numerical method as the density matrix renormalization group (DMRG),White  and thus precisely study the effects of frustration on hole-doped systems.

This paper is organized as follows. In Sec. II, we introduce the frustrated t-J ladder model, and then in Sec. III study a metal-insulator transition upon hole doping. By computing spin correlation functions, we clarify the origin of the first-order phase transition found there. In Sec. IV, spin excitations in a metallic phase are discussed in detail. It is clarified to what extent a spin-gap metallic phase persists around the insulating phase. In particular, we show that a spin triplet excitation in the spin-gap phase can be described in terms of two independent quasiparticles. We obtain the phase diagram in Sec. V, and summary and discussions are given in Sec. VI.

## II Frustrated Ladder Model

The model we treat here is a frustrated electron ladder which is an extension of the spin-1/2 Heisenberg model with the orthogonal-dimer structure introduced by Gelfand. Gelfand  The model is illustrated schematically in Fig. 1. It has two competing antiferromagnetic interactions $`J_1`$ and $`J_2`$, giving rise to strong geometrical frustration. In a hole-doped case, we use the t-J model where electron hopping is assumed to have two distinct values, as shown in Fig. 1. The Hamiltonian is thus given by

$`H`$ $`=`$ $`{\displaystyle \underset{(i,j),\sigma }{}}t_\alpha c_{i\sigma }^{}c_{j\sigma }+{\displaystyle \underset{(i,j)}{}}J_\alpha \left(𝐒_i𝐒_j{\displaystyle \frac{1}{4}}n_in_j\right),`$ (1)
where the summation is taken over nearest neighbor bonds ($`i`$,$`j`$) indicated by the thick and thin lines in Fig. 1. Here, $`c_{i\sigma }`$ is the annihilation operator of an electron with spin $`\sigma `$ at the $`i`$-th site and $`n_i`$ ($`S_i`$) is the electron number (spin-1/2) operator. As mentioned above, hopping $`t_i`$ (exchange coupling $`J_i`$) takes either $`t_1`$ or $`t_2`$ ($`J_1`$ or $`J_2`$). Note that doubly occupied state at each site are implicitly forbidden in the t-J model.

We note here that the above model Gelfand ; Honecker ; Brenig ; Kim  has several remarkable properties in addition to those for ordinary 1D frustrated spin systems: (i) it is a 1D analogue of the 2D orthogonal-dimer model (Shastry-Sutherland model Shastry ; Miyahara ) relevant for SrCu2(BO3)2, for which the ground state is exactly given by a direct product of decoupled dimers; (ii) it is also regarded as a 1D variant of the pyrochlore system (so-called linked-tetrahedra chain), which has frustrated tetrahedra as a key structure. Gelfand  A remarkable point related to (ii) is that, for noninteracting electrons on this lattice, the model has a dispersionless flat-band mode,

$`\epsilon (k)=\{\begin{array}{c}t_1\hfill \\ t_14t_2\mathrm{cos}(k),\hfill \end{array}`$ (4)
which reflects special geometry of the system, as is the case for the 3D pyrochlore lattice. Therefore, using this simplified model we can also explore the role of a flat-band mode, which is important to understand properties of the pyrochlore system near half filling.

## III Metal-Insulator transition

Let us start our discussions with the ground-state properties of the t-J model in a metallic phase. We particularly focus on the nature of a metal-insulator transition induced by hole doping. Here, we use a variant of the infinite DMRG method, so-called product-wave-function renormalization group method Nishino , to compute several quantities in a bulk limit.

In the undoped spin model, the system is either in the exact dimer-singlet phase ($`J_2/J_1<0.71`$) or the plaquette-singlet phase ($`J_2/J_1>0.71`$).Gelfand  We thus discuss the hole-doping effects on these phases separately. In the following analyses, we set $`t_2/t_1=J_2/J_1\gamma `$, and mainly show the results of $`J_1/t_1=0.25`$ ($`t_1=1.0`$) for simplicity.

### III.1 dimer phase

We begin with the dimer phase, which is stabilized at half filling for the ratio of the exchange couplings, $`\gamma <0.71`$. The calculated electron density $`n`$ is shown in Fig. 2(a) as a function of the chemical potential $`\mu `$.

It is seen that hole-doping smoothly drives the system to a metallic phase. At quarter filling ($`n=1/2`$), there appears a plateau in $`n`$, which is due to formation of a CDW insulating state. It is instructive to compare this smooth transition with first-order transitions found in the magnetization curve at half filling in the dimer phase: Honecker  the magnetization exhibits jumps among the plateaus at zero, half, and full magnetization. The discontinuity in the magnetization reflect the fact that a triplet excited state in the dimer phase is completely localized due to the orthogonal-dimer structure. This is not the case for hole doping: doped holes are indeed heavy but still mobile, allowing a smooth transition to the metallic phase.

We next compute spin correlation functions, which are defined in the rung ($``$) and chain ($``$) directions as,

$`\overline{S}_{}`$ $`=`$ $`S_1^zS_2^z,`$
$`\overline{S}_{}`$ $`=`$ $`S_1^zS_3^z(=S_1^zS_4^z).`$
Shown in Fig. 2(b) are the computed spin correlation functions. As the system approaches half filling from a metallic side, $`\overline{S}_{}`$ along the rung rapidly decreases down to $`1/4`$ whereas $`\overline{S}_{}`$ along the chain becomes almost zero. Recall that the values of $`\overline{S}_{}=1/4`$ and $`\overline{S}_{}=0`$ are expected for a system with isolated dimers. Therefore, the above results of correlation functions imply that the system near half filling is regarded as a resonating metallic state of decoupled dimers.

### III.2 plaquette phase

When holes are doped into the plaquette-singlet phase, the nature of a metal-insulator transition is completely changed. Shown in Fig. 3(a) is the electron density as a function of the chemical potential. It is remarkable that the system is driven to a metallic phase via a first-order phase transition, which is characterized by a jump in the electron density.

We can see the origin of the first-order metal-insulator transition by checking short-range spin correlation functions in the rung as well as leg directions, shown in Fig. 3(b). At half filling, $`\overline{S}_{}>0`$ and $`\overline{S}_{}<0`$, being consistent with the plaquette-singlet ground state. However, once the system is driven to a metallic phase via the first-order transition, both of $`\overline{S}_{}`$ and $`\overline{S}_{}`$ become negative. Furthermore, when the system approaches half-filling from a metallic side, $`\overline{S}_{}`$ decreases with negative values while $`\overline{S}_{}`$ gets very small. Namely, in a metallic phase close to half filling dimer correlations along the rung are largely enhanced. Hence, in the vicinity of half filling, this dimer-dominant state strongly competes with the plaquette-singlet state energetically, giving rise to the first-order phase transition observed above.

We note here that such a first-order metal-insulator transition only occurs when holes are doped in the plaquette phase, but not in the dimer phase. It is also instructive to compare this with the results known for the magnetization curve at half filling, where first-order phase transitions occur only in the dimer phase ($`\gamma <0.71`$). Honecker  These apparently different aspects in phase transitions come from the fact that hole doping induces a competition between the dimer- and plaquette-singlet states, whereas a magnetic field induces a competition between the dimer- (or plaquette-) singlet state and the triplet state.

## IV Spin excitations in a metallic phase

We now discuss spin excitations in a hole-doped system near half filling. We start with low-energy excitations in a system with two holes. As mentioned above, the system is driven to a metallic phase immediately upon hole doping for the dimer phase ($`0<\gamma <0.71`$). On the other hand, for the plaquette phase ($`\gamma >0.71`$), doped holes may be in a localized state, as a result of a first-order transition. Namely, holes are accommodated in a degenerate dispersionless level up to a critical hole density. In the following, we thus pay our attention to the former region, $`0<\gamma <0.7`$, which is more interesting as far as the system with a small number of holes is concerned.

In order to discuss excitations in a doped case, we consider two quantities, i.e. the spin gap $`\mathrm{\Delta }_S(N_h)`$ and the binding energy $`\mathrm{\Delta }_B(N_h)`$ respectively defined for $`N_h`$ (even number) holes on $`2\times L`$ lattice sites as,

$`\mathrm{\Delta }_S(N_h)`$ $`=`$ $`E_0(N_h,S^z=1)E_0(N_h,S^z=0),`$ (5)
$`\mathrm{\Delta }_B(N_h)`$ $`=`$ $`2E_0(N_h1,S^z=1/2)E_0(N_h,S^z=0)`$ (6)
$`E_0(N_h2,S^z=0).`$
In Fig. 4, we show the spin gap $`\mathrm{\Delta }_S`$ and the binding energy $`\mathrm{\Delta }_B`$ computed for two doped holes ($`N_h=2`$). The inset shows a finite-size scaling of $`\mathrm{\Delta }_S`$ and $`\mathrm{\Delta }_B`$ for $`\gamma =0.6`$. It should be noticed that the spin gap and the binding energy monotonically increase as $`\gamma `$ increases from $`\gamma 0.35`$. It is seen that the spin gap $`\mathrm{\Delta }_S`$ just coincides with the binding energy $`\mathrm{\Delta }_B`$, implying that a pair-breaking excitation in a metallic phase may be described by two independent quasiparticles. It is to be noted that the above spin gap in a doped system is totally different from the spin gap defined at half filling, $`\mathrm{\Delta }_S=J_1`$, as we will see below.

Shown in Fig. 5 is a spatial distribution of the hole density $`n_i^h`$ ($`n_i^h1n_i`$) and the spin density $`S_i^z`$ for two holes on a $`L=64`$ system. Here filled circles and open triangles show the results for the singlet ground state and for a triplet excited state, respectively. For $`\gamma =0.6`$, as shown by filled circles in Fig. 5(c), a hole-pair exists in the ground state, implying the existence of an attractive interaction between two holes in the ground state. For a triplet excitation, the bound state is broken and two quasiparticles appear, as seen from open triangles in Fig. 5(c) and (d). On the other hand, the situation is somewhat different for $`\gamma =0.3`$, where the repulsive interaction seems to exist even in the ground state, as shown by filled circles in Fig. 5(a) and (b). We have checked that a repulsive interaction in the latter case is mainly caused by finite-size effects. Therefore, it is expected that such a repulsive interaction becomes small for large $`L`$, and in the thermodynamic limit, doped holes may behave as two independent quasiparticles in the case of $`\gamma =0.3`$.

We now turn to the case with a finite density of holes. In Fig. 6, we show the spin gap as a function of the electron density $`n`$ for $`\gamma =0.6`$. The dotted lines with filled circles are the results of $`\mathrm{\Delta }_S`$ for finite-size systems. It is seen that the value of the spin gap dramatically changes between half-filling and two-hole case, and then smoothly changes in the region $`N_h2`$. We have performed a finite-size scaling analysis of the spin gap from the data shown by dotted lines, obtaining the curve shown by the solid line. Note that the spin gap has a discontinuity at half filling, implying that the spin gap in a metallic phase has a different origin from the half-filled case. Also shown by the broken line is the binding energy $`\mathrm{\Delta }_B`$ computed by a similar finite-size scaling analysis. Since $`\mathrm{\Delta }_S`$ is nearly equal to $`\mathrm{\Delta }_B`$, we can say that excited states may be described by two independent quasiparticles, as we have already seen in a system with two holes.

In order to see the nature of the attractive interaction between holes more clearly, we compute short-range hole-hole correlation functions $`\overline{H}_\alpha `$ in a bulk limit, which are defined in the rung ($``$) and chain ($``$) directions,

$`\overline{H}_{}`$ $`=`$ $`n_1^hn_2^hn^h^2,`$
$`\overline{H}_{}`$ $`=`$ $`n_1^hn_3^hn^h^2.`$
In Fig. 7, we show $`\overline{H}_\alpha `$ as a function of $`n`$ for $`\gamma =0.3`$ and $`0.6`$. For $`\gamma =0.6`$, there is the region with both of $`\overline{H}_{}`$ and $`\overline{H}_{}`$ being positive, where an effective interaction between holes is attractive in both directions. On the other hand, for $`\gamma =0.3`$, both of $`\overline{H}_{}`$ and $`\overline{H}_{}`$ are negative (i.e. repulsive interaction) in the whole range shown in the figure. By comparing these results with Fig. 4, we see that the region, where holes attract each other, roughly corresponds to the spin-gap metallic region. This means that the mechanism of spin-gap formation in a metallic phase, which is different from that in the insulating phase, is closely related to the attractive force between holes.

## V Phase diagram

Based on the above results, we now obtain the ground-state phase diagram, which is shown in Fig. 8. As found by Gelfand,Gelfand  the ground state at half filling is either in the dimer-singlet (or plaquette-singlet) phase with spin gap at half filling for $`\gamma <0.71`$ ($`\gamma >0.71`$). These two phases are separated via a first-order phase transition when $`\gamma `$ is changed.

We have found that the transition from the insulating plaquette phase to a metallic phase is of first-order, being contrasted to the dimer phase showing a second-order transition. Furthermore, a spin-gap metallic phase, which is indicated by the shaded region in Fig. 8, persists in the vicinity of half-filling, which is indeed caused by strong frustration. For the value $`J_1/t_1=0.25`$ shown in Fig. 8, we have found that a spin-gap metallic phase appears in the region of $`\gamma >0.35(\gamma _S)`$.

When the electron density is further decreased, we encounter a CDW insulating phase at quarter filling, characterized by the plateau in $`n`$-$`\mu `$ curve. The region of a CDW phase becomes narrower as $`\gamma `$ increases, and finally vanishes beyond a certain critical value. As is seen from Fig. 2(b) and Fig. 3(b), a metallic phase continuously spreads in the whole region of $`\gamma `$, where short-range spin correlations along the rung, $`\overline{S}_{}`$, always take negative values. This implies that dimer correlations are dominant for the ground state in a metallic phase, as is usually the case for a doped antiferromagnetic ladder. However, it should be noticed that in a metallic phase close to the plaquette-singlet insulator, a competition between the dimer- and plaquette- singlet states may give rise to dual properties in physical quantities.

We have thus far shown the calculated results by choosing a specific ratio of hopping and exchange. It is necessary to check what is changed when we take other choices of the parameters. For instance, we show the spin gap of a system with two holes in Fig. 9 when the ratio of $`J_1/t_1`$ is changed. As $`J_1/t_1`$ decreases, the critical point, $`\gamma _S`$, where the spin gap starts to develop, is shifted to a larger value. Note, however, that $`\gamma _S0.5`$ in the limit of $`J_1t_1`$. Namely, even if $`J_1/t_1`$ is small, a spin-gap metallic phase appears at least in the region of $`0.5<\gamma <0.71`$. This feature is different from the case of the t-J ladder without frustration, Tune\_lad ; Dagotto ; WhitScala ; Siller ; Poilblanc ; Hayward  which has a spin gap only for relatively large $`J_1/t_1`$.

## VI Summary and Discussions

We have studied the hole-doping effects on a spin ladder system with the orthogonal-dimer structure, which also possesses linked-tetrahedra as a key structure. We have investigated how a metal-insulator transition is affected by geometrical frustration, by computing the electron density and several correlation functions by means of the DMRG method.

It has been found that when the system is in the dimer phase, hole-doping smoothly drives the system to a metallic phase. However, when the system is in the plaquette phase, a first-order transition occurs upon hole doping. The latter first-order transition occurs as a consequence of a strong competition between the dimer- and plaquette-singlet states. This is different from ordinary metal-insulator transitions for a frustrated t-J chain with next-nearest neighbor exchange couplings studied so far.Ogata

The first-order phase transition found here is closely related to unique geometry of the linked-tetrahedra lattice, which causes a dispersionless (spatially isolated) mode, as mentioned in the beginning of the paper. It is to be noted that such a dispersionless mode caused by special geometry exists for the 3D pyrochlore lattice and also for the 2D orthogonal-dimer lattice. In this sense, the present results may capture some of characteristic properties common in this class of fully frustrated systems with doped holes. For instance, a pyrochlore antiferromagnetic system, indeed shows a competition between the dimer and plaquette phases, being similar to the present results. Koga  Therefore, we think that unusual properties may appear in a metallic phase of the doped pyrochlore system around the insulating phase. This is consistent with the results deduced by preservative analyses. Isoda ; Fujimoto  We also expect that such peculiar properties may emerge even in 2D orthogonal-dimer systems when holes are doped into an insulating phase around the boundary between the dimer and plaquette phases. Koga\_let  Detailed theoretical studies on hole-doped pyrochlore systems as well as 2D orthogonal-dimer systems are left as an interesting open problem.

## Acknowledgements

This work was partly supported by a Grant-in-Aid from the Ministry of Education, Science, Sports and Culture of Japan. A part of computations was done at the Supercomputer Center at the Institute for Solid State Physics, University of Tokyo and Yukawa Institute Computer Facility. A. Kawaguchi was supported by Japan Society for the Promotion of Science.


format: text

text:
# X-ray Timing of the Enigmatic Neutron Star RX J0720.4-3125

## 1 Introduction

RX J0720.4$``$3125 was discovered by Haberl et al. (1997) as a soft ($`kT80`$ eV), bright X-ray source in the ROSAT All-Sky Survey. Given its very low hydrogen column density ($`N_H1\times 10^{20}\text{ cm}^2`$), nearly sinusoidal 8.39-s pulsations, relatively constant X-ray flux, and very faint ($`B=26.6`$ mag), blue optical counterpart (Kulkarni & van Kerkwijk, 1998; Motch & Haberl, 1998), it was classified as a nearby, isolated neutron star.

As one of the closest ($`d300`$ pc; Kaplan, van Kerkwijk, & Anderson 2002) neutron stars, RX J0720.4$``$3125 occupies a central position in our study of these objects. However, the long period is puzzling, and has led to three models: an old, weakly magnetized neutron star accreting matter from the interstellar medium (Wang, 1997; Konenkov & Popov, 1997); a middle-aged pulsar with $`10^{12}`$ G magnetic field whose radio beams are directed away from the Earth (Kulkarni & van Kerkwijk, 1998); or an old magnetar (neutron star with magnetic field $`>10^{14}`$ G; Duncan & Thompson 1992) that is kept warm by the decay of its strong magnetic field (Heyl & Hernquist, 1998; Heyl & Kulkarni, 1998). These models predict different period derivatives: $`\dot{P}<5\times 10^{15}\text{ s s}^1`$, $`\dot{P}10^{15}`$$`10^{13}\text{ s s}^1`$, and $`\dot{P}\mathrm{few}\times 10^{13}\text{ s s}^1`$, respectively.

Motivated thus, we undertook timing observations of RX J0720.4$``$3125 using the Chandra X-ray Observatory (CXO), supplemented with analysis of archival data from ROSAT and BeppoSAX.

After submission of this paper, we became aware of the work of Zane et al. (2002) reporting a timing analysis of RX J0720.4$``$3125. A “Notes added in manuscript” section regarding that analysis can be found at the end of the manuscript.

## 2 Observations

The primary data consist of two sets of observations obtained from Chandra: one using the HRC in the spectroscopic mode (HRC-S) with the Low Energy Transmission Grating (LETG), and one using ACIS in the continuous clocking (CC) mode. The primary and archival datasets are summarized in Table 1.

We processed the HRC-S data using the standard pipeline111http://asc.harvard.edu/ciao/threads/spectra\_letghrcs/ and extracted $`0^{\mathrm{th}}`$ order events from a circle with radius 10 pixels ($`1\stackrel{}{\mathrm{.}}3`$). For the $`\pm 1^{\mathrm{st}}`$ orders, we extracted events from a region $`0.0006\mathrm{°}`$ wide in the cross-dispersion direction (the tg\_d coordinate) and from $`0.08\mathrm{°}`$ to $`0.35\mathrm{°}`$ along the dispersion direction (the tg\_r coordinate). We extracted events from the ACIS data within $`\pm 1\mathrm{}`$ of the source. We then used the axBary program to barycenter the events in both these datasets.

The best fit position for RX J0720.4$``$3125, found by averaging the $`0^{\mathrm{th}}`$ order data from the three HRC-S datasets, is (J2000) $`\alpha =07^\mathrm{h}20^\mathrm{m}24\stackrel{\mathrm{s}}{\mathrm{.}}96`$, $`\delta =31\mathrm{°}25\mathrm{}49\stackrel{}{\mathrm{.}}6`$, with rms uncertainty of $`0\stackrel{}{\mathrm{.}}6`$ in each coordinate due to CXO aspect uncertainties. This is consistent ($`1\stackrel{}{\mathrm{.}}4`$ away) with the optical position (Kulkarni & van Kerkwijk, 1998). The X-ray source appears unresolved and its profile is consistent with that of a point source (half-power radius of $`0\stackrel{}{\mathrm{.}}5`$).

For the ROSAT HRI data, we extracted the events within a circle of radius 45 pixels ($`22\stackrel{}{\mathrm{.}}5`$) centered on the source. We used a circle of radius 200 pixels ($`100\mathrm{}`$) for the PSPC data. These events were barycentered using the ftools programs abc and bct and corrected to Barycentric Dynamical Time (TDB) according to Cox (2000, p. 14).

We extracted the BeppoSAX LECS events within a circle with radius of 25 pixels ($`200\mathrm{}`$) centered on the source and restricted to those with pulse-invariant (PI) amplitudes that were less than 90 (energies $`<0.95`$ keV), in order to maximize the signal-to-noise. Finally, we barycentered the events with the SAXDAS tool baryconv.

## 3 Timing Analysis

For each dataset, we computed $`Z_1^2`$ power spectra around the known 8.39-s period. Specifically, we explored the period range from 8.376 s to 8.405 s in steps of $`7\text{ }\mu \mathrm{s}`$ (oversampling by factors of 20–800 relative the nominal step-size of $`P_0^2/\mathrm{\Delta }T`$, where $`P_0=8.39`$ s is the approximate period and $`\mathrm{\Delta }T`$ is the span of the dataset from Table 1). As can be seen from Figure 1 all but the Chandra HRC-S and ROSAT HRI-2 datasets yielded unambiguous period estimates. For the HRC-S and HRI-2 sets the period estimates are ambiguous because the large gaps in the observations result in strong side-lobes. In Figure 2, we display the best-fit periods for the unambiguous determinations as well as possible periods for the HRC-S and HRI-2 datasets.

As can be seen from Figure 2, the ambiguity of the HRC-S and HRI-2 datasets can be resolved provided we assume (reasonably) that the period evolves smoothly with time. Our choice of period (for HRC-S and HRI-2) and the best fit periods (for the other datasets) are shown in Table 1. The errors on the periods were determined using the analytical expression from Ransom (2001). While that expression was derived for FFT power spectra, $`Z_1^2`$ power spectra have the same statistics (both are exponentially distributed) so the same relations should apply (we have verified this with numerical simulations). We also show in Table 1 times-of-arrival (TOAs) for each of the datasets.

The data in Table 1 are consistent with there being no measurable $`\dot{P}`$: fitting for a linear spin-down gives $`P=8.391115(8)`$ s at MJD 51633 and $`\dot{P}=(1\pm 12)\times 10^{14}\text{ s s}^1`$, with $`\chi ^2=1.54`$ for 4 degrees-of-freedom (DOF). If, instead, we fit only for a constant period, we find $`P=8.391115(8)`$ s, with $`\chi ^2=1.64`$ for 5 DOF. Therefore we can constrain the secular period derivative to be $`|\dot{P}|<3.6\times 10^{13}\text{ s s}^1`$ (3-$`\sigma `$).

The folded lightcurve is largely sinusoidal, with an rms pulsed-fraction (the rms of the lightcurve divided by the mean) of 8% for both Chandra datasets. However, this pulsed-fraction is energy-dependent: the fraction rises with decreasing energy (see Figure 3), in agreement with the XMM analysis (Paerels et al., 2001).

## 4 Discussion

Our upper limit of $`|\dot{P}|<3.6\times 10^{13}\text{ s s}^1`$ is sufficiently high that we cannot meaningfully constrain the accretion model, for which we expect $`\dot{P}<5\times 10^{15}\text{ s s}^1`$ (the limit is the case where all of the required $`\dot{M}`$ of $`10^{12}\text{ g s}^1`$ couples to the neutron star at the corotation radius, giving maximum torque per unit mass). We note that the accretion model is no longer viable for another similar but better studied isolated neutron star, RX J1856.5$``$3754 (van Kerkwijk & Kulkarni, 2001). Regardless, the accretion model is best confronted by measuring the proper motion and distance, and looking for evidence of sufficiently dense ambient gas (deep Keck H$`\alpha `$ imaging and HST astrometric observations are in progress).

However, we can constrain the pulsar and magnetar models. We can draw four inferences common to both models. First, the spin-down luminosity, $`\dot{E}=I\dot{\omega }\omega <2.4\times 10^{31}`$ erg s-1; here, $`\omega =2\pi /P`$. Second, in the framework of a simple (vacuum magnetic dipole radiation) pulsar model, the physical age is roughly approximated (provided the current spin period is much larger than that at birth and that the magnetic field does not decay significantly) by the so-called characteristic age: $`\tau _cP/(2\dot{P})>4\times 10^5`$ yr. Third, the strength of the dipole field is $`B=3.2\times 10^{19}(P\dot{P})^{1/2}<6\times 10^{13}`$ G. Fourth, we assume that the X-ray emission (well described by a blackbody; Haberl et al. 1997; Paerels et al. 2001) is cooling flux from the surface. The bolometric cooling luminosity is $`L_{\mathrm{cool}}2\times 10^{32}d_{300}^2\text{ ergs s}^1`$ (Haberl et al., 1997), using the distance estimate of $`300d_{300}`$ pc derived by scaling from RX J1856.5$``$3754 (Kaplan et al., 2002).

Knowledge of $`L_{\mathrm{cool}}`$ enables us to estimate the cooling age, $`t_{\mathrm{cool}}`$, of RX J0720.4$``$3125. Magnetic fields, especially intense $`B`$ fields such as those proposed for magnetars, can profoundly influence the cooling of neutron stars. To this end, we use the curves of $`L`$ vs. $`t`$ from Heyl & Kulkarni (1998) and find $`t_{\mathrm{cool}}(510)\times 10^5`$ yr, assuming a 50% uncertainty in the distance and with only a slight dependence on $`B`$. This age is consistent with the characteristic age derived above.

In the magnetar model (we assume that the $`B`$-field decay is dominated by the slower irrotational mode; see Heyl & Kulkarni 1998) the expected $`B`$ field at about $`10^6`$ yr is, $`2\times 10^{14}`$ G, well above the upper limit obtained from our $`\dot{P}`$ limit. Models that are consistent with both $`L_{\mathrm{cool}}`$ and our limit on $`\dot{P}`$ (and thus an upper limit on the current value of $`B`$) are those with initial $`B10^{14}`$ G. Based on this, we conclude that RX J0720.4$``$3125 is not a magnetar, motivating us to consider the pulsar model.

Earlier, Kulkarni & van Kerkwijk (1998) did not accept the radio pulsar model because in 1998 there were no radio pulsars with such long periods. However, over the past four years we have come to appreciate the existence of pulsars with $`B>10^{13}`$ G (Camilo et al. 2000; see also Figure 4). In particular, the parameters of RX J0720.4$``$3125 are not too dissimilar to those of PSR J1814$``$1744, which has $`P4`$ s and $`\dot{P}7.4\times 10^{13}\text{ s s}^1`$. Thus, the past objections against the radio pulsar model are no longer tenable, and RX J0720.4$``$3125 seems fully compatible with being an off-beam high-$`B`$ pulsar. If that is the case, then we expect $`\dot{P}10^{13}\text{ s s}^1`$, a value that we should be able to measure in the near future.

A separate possibility is that RX J0720.4$``$3125 is an off-beam pulsar with age compatible with $`t_{\mathrm{cool}}`$, but with a conventional ($`10^{12}`$ G) magnetic field and $`\dot{P}10^{15}\text{ s s}^1`$ (and therefore $`\tau _c10^8`$ yr). RX J0720.4$``$3125 could then be similar to the 8.5-s, $`2\times 10^{12}`$-G pulsar PSR J2144$``$3944 (Young, Manchester, & Johnston 1999; see Figure 4). With a braking index of 3, the age of a pulsar is $`\tau =\tau _c\left(1(P_0/P)^2\right)`$, where $`P_0`$ is the initial spin-period. If RX J0720.4$``$3125 does have $`\tau t_{\mathrm{cool}}`$ and $`B10^{12}`$ G, we find and $`P_08.3`$ s, very close to $`P`$. Such as a pulsar would be an example of the “injection” hypothesis (Vivekanand & Narayan, 1981), where pulsars are born with initial spin periods $`P_010`$ ms (as for the Crab). Such long initial periods are allowed and perhaps expected in some models of neutron-star formation (e.g., Spruit & Phinney, 1998), where the precise initial period depends very sensitively on the details of the formation mechanism and may range over four orders of magnitude. While there are a few pulsars whose characteristic ages are factors of 10–100 times the ages derived from supernova remnant associations (Pavlov et al., 2002), this would be the first case for a source with $`P_0>1`$ s.

We make the following parenthetical observation: for most known pulsars the X-ray pulsed-fraction (largely) increases with photon energy (Perna, Heyl, & Hernquist, 2001), whereas for RX J0720.4$``$3125 we see the opposite effect. However, for the pulsars (e.g., PSR B0656+14, PSR B1055$``$52), the X-ray luminosity has a strong, highly-pulsed, non-thermal component with $`L_{X,\text{non-th}}10^3\dot{E}`$ (Becker & Trümper, 1997). Furthermore, in such objects heating of the polar caps by pulsar activity (probably dependent on $`\dot{E}`$) is likely significant. The interplay of these components with the viewing geometry can result in the large range of observed phenomena (e.g., Perna et al., 2001). For RX J0720.4$``$3125, though, with its small $`\dot{E}`$ there is little reason to expect a strong non-thermal contribution or a hot polar cap (although there must be some inhomogeneities to give the observed pulsations). We conjecture that the increase in the pulse fraction with decreasing photon energy is primarily due to the absence the additional components.

## 5 Conclusions

In this Letter, based on X-ray timing data and cooling models, we argue that the nearby soft X-ray source RX J0720.4$``$3125 is not a middle-aged magnetar but is likely a $`10^6`$-yr off-beam pulsar. To accommodate its age and long period we speculate that it either has $`B10^{13}`$ G or was born with $`P_08.3`$ s, a very surprising result as both source types are, at present, considered to be rare. We now consider the larger ramifications of our conclusions.

A volume-limited sample of neutron stars offers us an opportunity to sample the diversity of such sources. In this respect, soft X-ray surveys provide the best such samples since all neutron stars — normal radio pulsars, high-$`B`$ pulsars, magnetars, and the mysterious Cas-A-like neutron stars — will cool through soft X-ray emission well into their middle ages. Indeed, this expectation is borne out by the local sample: pulsars such as PSR B0656+14, Geminga222Presumably a standard pulsar that is not beamed toward us., RX J0720.4$``$3125, a youngish magnetar (see below), and finally the very mysterious RX J1856.5$``$3754, of which we know nothing other than it is a cooling neutron star333It is further worth noting that the sample of soft X-ray neutron stars has at least three long-period objects (Haberl, Pietsch, & Motch, 1999)..

The number of neutron stars belonging to a given class depends not only on the sensitivity of the X-ray survey but also on the product of the birth rate and the cooling age. Thus, for example, magnetars with their longer-lasting cooling radiation may dominate the local population despite a lower birthrate (Heyl & Kulkarni, 1998). This bias and the long period led us to speculate that RX J0720.4$``$3125 was an old magnetar, a conclusion we have now refuted. In contrast, the soft thermal X-ray source RX J1308.8+2127, with $`P=5.2`$s and $`\dot{P}10^{11}\text{ s s}^1`$ (Hambaryan et al., 2002) appears to have a magnetar-strength field.

The proximity of RX J0720.4$``$3125 argues for a substantial Galactic population of similar sources, but very few such radio pulsars are known. The cause of this paucity is that radio surveys select against long-period pulsars, especially those with $`B>10^{13}`$ G, in several ways. (1) The beaming factor is known to decrease with increasing period, reaching 3% at $`P10`$ s (Tauris & Manchester, 1998). (2) As can be seen from Figure 4, the lifetime of a radio pulsar decreases with increasing $`B`$: a $`B10^{12}`$ G neutron star crosses the the radio death line at $`10^8`$ yr whereas a $`B10^{13}`$ G pulsar dies at $`2\times 10^7`$ yr. The loss of throughput of a pulsar survey for a 5-s pulsar relative to a 1-s pulsar from these two effects alone is nearly one order of magnitude. (3) The true loss is even greater since long-period signals are frequently classified as interference (we note that population models do not constrain the population of long-period pulsars \[Hartman et al. 1997\], mainly due to reasons 1 and 2). Young high-$`B`$ pulsars in supernova remnants would almost certainly create visible plerion nebulae due to their high $`\dot{E}`$’s, while long-period injected pulsars of similar ages would be invisible except for their cooling radiation (without assuming that the radio beams are directed toward the Earth). Thus, injected pulsars detectable only via X-ray emission may be present in many “hollow” supernova remnants (i.e. those without visible plerions). Radio pulsar searches better tuned to long periods and very deep radio and X-ray searches for young pulsars in supernova remnants may uncover the postulated class of long-period sources.

We thank M. Cropper and S. Zane for discussions on the TOAs and for alerting us to the different time systems used by different programs. We thank D. Frail and J. Heyl for useful discussions, and we thank an anonymous referee for constructive comments. We have used the NASA-maintained HEASARC web site for archival data retrieval and subsequent analysis. D.L.K. holds a fellowship from the Fannie and John Hertz Foundation, and his and S.R.K.’s research are supported by NSF and NASA. M.H.v.K. is supported by a fellowship from the Royal Netherlands Academy of Arts and Sciences. Notes Added In Manuscript After we submitted our paper we became aware of a paper by Zane et al. (2002) reporting timing analysis of RX J0720.4$``$3125. Our period determinations of the archival data (PSPC, HRI-1, LECS, HRI-2, HRC-S) are in excellent agreement with those of Zane et al. Both papers also report new determinations, which are: Chandra ACIS-S3 (our paper) and XMMa (2000 May 13) and XMMb (2000 November 21; both from Zane et al.). We restricted our analysis to an incoherent combination of the various datasets, i.e. we looked for secular evolution of the period determined from each observation separately. We did not attempt to phase connect the datasets. Zane et al. do present a coherent analysis, using the archival and XMM data. However, in our opinion such an analysis is premature and not robust. First, it is premature, because the XMM derived periods of known pulsars have fractional errors $`|\mathrm{\Delta }P/P|`$ ranging from $`1.9\times 10^7`$ to $`1.2\times 10^5`$ (as reported by the XMM calibration team; Kuster et al. 2002). This error alone may result in systematic uncertainties as high as $`\dot{P}`$ of $`6\times 10^{12}`$ over the 6-month duration of the XMM datasets and $`5\times 10^{13}`$s s-1 over the entire span of the observations. Second, it is not robust, as phase connection without any ambiguity requires that the datasets be separated by time intervals less than “coherence” timescale, $`P^2/\sigma _P`$ (where $`\sigma _P`$ is the uncertainty in the measurement of $`P`$) and none of the datasets (including the XMM datasets) satisfy this condition. We note that neither of the two primary solutions from Zane et al. fits the TOAs in Table 1.


text:
# Comparison of space-time evolutions of hot/dense matter in √𝑠_{𝑁𝑁}=17 and 130 GeV relativistic heavy ion collisions based on a hydrodynamical model

## I Introduction

One of the main purposes of ultra-relativistic heavy ion collision experiments is to explore the property of the hot and dense matter QM2 (2002). Recently the new experiment has begun to work at BNL-RHIC of which higher collision energy than other experiments up to now provides us chances to produce a new state of matter, quark-gluon plasma (QGP) with distinct possibility. However, the complicated collision processes composed of multiparticle productions and many-body interactions make it difficult to understand the properties of the hot matter. Therefore, a simple dynamical description of the system as a basis for the deeper understanding is indispensable.

Relativistic hydrodynamical models are well-established phenomenological tools for describing high energy nucleus-nucleus collisions and subsequent multiparticle production Bjorken (1983); von Gersdorff et al. (1986); Akase et al. (1991); Alam et al. (1996); Ornik et al. (1996); Hung and Shuryak (1998); Kolb et al. (2000). In this paper, we use a (3+1)-dimensional hydrodynamical model Ishii and Muroya (1992) with a QCD phase transition. We assume a cylindrical symmetry to the collision dynamics. Thus, our discussion is limited to the central collisions only. By virtue of the simple picture of our model, we can easily analyze both SPS and RHIC data with the same numerical code. Most of hydrodynamical calculations for RHIC data use Bjorken’s scaling solution Bjorken (1983) for the longitudinal direction. For example, Kolb et al. analyzed hadronic tranverse mass spectra and anisotropic flow Kolb et al. (2001). Zschiesche et al. Zschiesche et al.  investigated the HBT radii. These calculations assume the longitudinal boost-invariant infinite source. Though recently some hybrid models are used Bass and Dumitru (2000); Teaney et al. (2001, ) for the description of the hadronic phase, we here use a conventional description in which the hadronic phase is in local equilibrium.

In this paper, concentrating our discussion on the central collisions, we reproduce the single-partcle spectra of hadrons at the beginning. In the hydrodynamical model, single-particle distributions are used as inputs in order to determine initial parameters rather than outputs. However, it is not trivial whether we succeed to reproduce experiments with “natural parameters” or not. Based on the solutions of hydrodynamical equations, we discuss the physical properties and the space-time evolution of the fluids in SPS and RHIC. We also evaluate the two-pion correlation functions and analyze the HBT radii. As a subsequent work of Ref. Hirano et al. , we focus our discussion on comparison of the RHIC results and the SPS results based on the same numerical code.

As is well-known, the two-particle correlation function gives us information on the size of the particle source Wiedemann and Heinz (1999); Weiner (2000). In the cases of the relativisitic heavy ion collisions, the correlation function tells us about the freeze-out which should be far from the static source. Thus, dynamical models such as hydrodynamical models are indispensable for understanding the relation between observed correlation functions and the space-time history of the system. However, up to now, any dynamical model assuming QGP failed to explain the experimental HBT radii in RHIC consistently with the single-particle spectra Zschiesche et al. ; Soff et al. (2001) and elliptic flow Heinz and Kolb , as being known as “HBT puzzle”.111In Ref. Humanic , a hadronic rescattering model is shown to reproduce these quantities. We study the HBT radii in the framework of a hydrodynamical model which takes account of both transverse and longitudinal flow appropriately with a simple initial condition.

In the next section, we explain our model. In Sec. III, we discuss the space-time evolution of the fluid. In Sec. IV, we present the result of two-particle correlation. Section V is devoted for the concluding remarks.

## II Hydrodynamical Model

Let the system achieve the local thermal and chemical equilibrium shortly after a collision of two incident nuclei. This relaxation process cannot be described by the hydrodynamical model. The hydrodynamical model starts at initial time, $`\tau _0`$, at which thermal and chemical equilibrium are established at least locally. The hydrodynamical equations are given as

$$_\mu T^{\mu \nu }=0.$$
(1)
We assume the perfect fluid for simplicity. Hence, energy-momentum tensor is given as

$$T^{\mu \nu }=(ϵ+P)U^\mu U^\nu Pg^{\mu \nu },$$
(2)
with $`U^\mu ,ϵ`$ and $`P`$ being four velocities of a fluid element, energy density and pressure, respectively. These are treated as local quantities. We numerically solve the above equations together with the net baryon number conservation law,

$$_\mu (n_\text{B}U^\mu )=0,$$
(3)
where $`n_\text{B}`$ is the net baryon number density and is also treated as a local quantity. Putting the $`z`$ axis as a collision axis, we use a cylindrical coordinate system as follows;

$`t`$ $`=\tau \mathrm{cosh}\eta ,`$ (4)
$`x`$ $`=r\mathrm{cos}\varphi ,`$ (5)
$`y`$ $`=r\mathrm{sin}\varphi ,`$ (6)
$`z`$ $`=\tau \mathrm{sinh}\eta .`$ (7)
Focusing our discussion on central collisions, we may assume the cylindrical symmetry on the system. Therefore, by virtue of an identity $`U_\mu U^\mu =1`$, the four velocity can be expressed by two rapidity-like varibles $`Y_\text{L}`$ and $`Y_\text{T}`$;

$`U^\tau `$ $`=\mathrm{cosh}(Y_\text{L}\eta )\mathrm{cosh}Y_\text{T},`$ (8)
$`U^\eta `$ $`=\mathrm{sinh}(Y_\text{L}\eta )\mathrm{cosh}Y_\text{T},`$ (9)
$`U^r`$ $`=\mathrm{sinh}Y_\text{T}.`$ (10)
Most of hydrodynamical calculations which analyze RHIC data use Bjorken’s scaling solution $`Y_\text{L}=\eta `$. Putting the solution as an ansatz reduces numerical tasks very much but the analyses are limited to midrapidity region only. We solve not only transverse expansion but also the longitudinal expansion explicitly. Numerical procedure for solving the coupled equations (1) and (3) is explained in Ishii and Muroya (1992). In this algorhythm, we solve the entropy and baryon number conservation law explicitly. Throughout our calculation, the total energy, entropy and baryon number are conserved within 5% of accuracy at the time step $`\delta \tau =0.01`$ fm/c.

In order to solve the hydrodynamical equations, we must fix the equation of state (EoS). We adopt a bag model EoS with a first order phase transition. The QGP phase is composed of a free gas of massless u, d, s quarks and gluons. Hadronic phase is also assumed to be a free gas but with excluded volume correction. All hadrons are included up to 2 GeV/$`c^2`$ of mass except for hyperons. Putting the critical temperature as $`T_\text{c}=160`$ MeV at vanishing baryon density, we get the bag constant $`B^{1/4}=233`$ MeV. We display the pressure as a function of temperature and baryonic chemical potential in Fig. 1. See Nonaka et al. (2000) for the further detail of the EoS and numerical treatment of the first order phase transition in solving the hydrodynamical equations.

We assume that the system achieves local equilibrium and begins to expand hydrodynamically at $`\tau =\tau _0=1.0`$ fm/$`c`$. We put the initial conditions on this hyperbola. Bjorken’s scaling solution is used as the initial condition of the longitudinal flow. Transverse flow is simply neglected at the initial time. We parameterize the initial energy density distribution $`ϵ(\tau _0,\eta ,r)`$ and net baryon number density distribution $`n_\text{B}(\tau _0,\eta ,r)`$ as,

$`ϵ(\tau _0,\eta ,r)`$ $`=ϵ_{\text{max}}\mathrm{exp}\left[{\displaystyle \frac{(|\eta |\eta _0)^2}{2\sigma _\eta ^2}}\theta (|\eta |\eta _0){\displaystyle \frac{(rr_0)^2}{2\sigma _r^2}}\theta (rr_0)\right],`$ (11)
$`n_\text{B}(\tau _0,\eta ,r)`$ $`=n_{\text{B0}}\left\{\mathrm{exp}\left[{\displaystyle \frac{(\eta \eta _\text{D})^2}{2\sigma _\text{D}^2}}\right]+\mathrm{exp}\left[{\displaystyle \frac{(\eta +\eta _\text{D})^2}{2\sigma _\text{D}^2}}\right]\right\}\mathrm{exp}\left[{\displaystyle \frac{(rr_0)^2}{2\sigma _r^2}}\theta (rr_0)\right].`$ (12)
The energy density distribution of the longitudinal direction (11) has a central plateau characterized by $`\eta _0`$ and a gaussian tail whose width is given by $`\sigma _\eta `$ (Fig. 2), while the net baryon number distribution is a superposition of the two gaussians of which peaks exist at $`\pm \eta _\text{D}`$. For the transverse direction, both are parametrized by a flat region with gaussian smearing near the edge (Fig. 3). For a nucleus with mass number $`A`$, the relation among these quantities is given by $`\sigma _\text{r}+r_0=1.2A^{1/3}`$. 222We adopt the initial condition as a natural and the simplest extension of the (1+1)-dimensional Bjorken’s picture Bjorken (1983) and as a basis for the further improvement.  Once these parameters are fixed, we can solve the hydrodynamical equations and pursue the space-time evolution of the fluid. These initial parameters are so chosen that the model reproduces the single-particle spectra measured in the experiments. The single-particle spectra can be calculated by making use of of the Cooper-Frye formula Cooper and Frye (1974)

$$E_𝐤\frac{dN_i}{d^3𝐤}=\frac{g_i}{(2\pi )^3}_\mathrm{\Sigma }k_\mu 𝑑\sigma ^\mu \frac{1}{\mathrm{exp}[(U_\nu k^\nu \mu _\text{B})/T_\text{f}]1},$$
(13)
where $`g_i`$ is a degeneracy of the hadrons and $`T_\text{f}`$ is a freeze-out temperature. The sign is plus for fermions and minus for bosons. Integration is performed on 3-dimensional freeze-out hypersurface $`\mathrm{\Sigma }`$. By virtue of the Lagrangian hydrodynamics, contribution from the time-like hypersurface is expected to be small and the space-like hypersurface dominates the particle emission at freeze-out; we employ the non-covariant prescription $`k_\mu d\sigma ^\mu k_\tau d\sigma ^\tau `$ for the sake of simplicity in the numerical treatment. In this approximation, total counted energy evaluated from Eq.(13) is slightly larger than 90% of the total energy of the initial fluid; we regard the approximation works well enough. At SPS, we assume that the freeze-out occurs at a energy density $`ϵ_\text{f}`$ and at a temperature $`T_\text{f}`$, at the RHIC energy. We also assume that the thermal and the chemical freeze-out are taken to happen simultaneouly. We show the freeze-out lines and the phase boundary on $`T\mu _\text{B}`$ plane in Fig. 4. Note that two freeze-out lines in the figure do not differ at low baryonic chemical potential (Fig. 4).

We take account of the particles emitted from resonance decay as well as the direct emission from the freeze-out hypersurface. We include the decay processes $`\rho 2\pi `$, $`\omega 3\pi `$, $`\eta 3\pi `$, $`K^{}\pi K`$, and $`\mathrm{\Delta }N\pi `$ Sollfrank et al. (1991); Hirano (2001). These resonances are also assumed to be thermally emitted from the freeze-out hypersurface.

Two sets of initial parameter are summarized in Table 1. Figures 5-7 show single-particle spectra in 17 $`A`$GeV Pb+Pb collisions at SPS. Our model well reproduces the experimental data with parameters in Table 1. Also in 130 $`A`$GeV Au+Au collisions, our model shows good agreement with the data as in Figs. 8-11. However, we note that our model fails to produce enough number of anti-protons and overestimates the kaon yield in Fig. 10, where we multiply factors of 0.6 for kaons and 3.5 for anti-protons for clear comparison of the slopes Hirano et al. . This discrepancy may indicate the need for more sophisticated freeze-out mechanism.

## III Space-time evolution

In this section, we present the numerical solution of the hydrodynamical equations and discuss differences in the space-time evolutions of the fluids between RHIC and SPS. Figures 12 and 13 show space-time evolution of the fluid on the transverse plane. We also display the space-time evolution on the longitudinal plane in Fig. 14. Longitudinal flow and transverse flow are shown in Figs. 15 and 16, respectively. From Table 1, the maximum energy density in RHIC is only 5% higher than the one in SPS. Only 5% higher energy density for the almost 50% larger $`dN/dY`$ seems surprising result. We show the number density of the thermal negative pions emitted into midrapidty as a function of the space-time rapidity $`\eta `$ of the freeze-out point (Fig. 17). This figure informs us that thermal contribution of the volume element at $`\eta =0`$ to the particles into midrapidity is only 9% larger in RHIC than in SPS. However, the wider region of the freeze-out hypersurface in $`\eta `$ contributes to the midrapidity particle distribution in RHIC more than in SPS. As a result, factor 1.5 times larger number of particles obtained in the midrapidity region after summing up particles emitted at different $`\eta `$. The difference between RHIC and SPS in Fig. 17 originates in the longitudinal extent $`\eta _0+\sigma _\eta `$ (see also Fig. 2) and is direct consequence of longitudinal dynamics. We also plot the entropy per unit flow rapidity $`dS/dY_\text{L}`$ in Fig. 18. This is a conserved quantity if the boost-invariance is kept. In both RHIC and SPS, reflecting the deviation from the scaling solution shown in Fig. 15, entropy is shifted to the larger flow rapidity. Reduction of entropy at $`Y_\text{L}=0`$, where $`Y_\text{L}=\eta =0`$ always holds, comes from $`dY_\text{L}/d\eta `$ which is larger than unity Eskola et al. (1998). Thus, the shift at RHIC is smaller than the one at SPS since the deviation from the boost-invariant solution is small (Fig. 15). At SPS, the difference of $`dS/dY_\text{L}`$ between the initial stage and the final stage is larger than the case at RHIC.

Though our maximum energy density 6.0 GeV/fm3 at initial is also so smaller than other calculation Kolb et al. (2001); Zschiesche et al. ; Hirano (2002), this is due to the difference of initial time and transverse energy density profile of which nuclear thickness is considered. As for the average energy density at midrapidity, we get $`ϵ_{\text{RHIC}}=3.9`$ GeV/fm3 and $`ϵ_{\text{SPS}}=3.77`$ GeV/fm3. $`ϵ_{\text{RHIC}}`$ is a little smaller than an estimation of Ref. Adcox et al. (2001). As a result of such a little difference in energy density, the space-time evolutions of two cases do not alter much in Figs. 12 and 13. The most different point is a longitudinal extension of the fluid, $`\eta _0+\sigma _\eta `$. In RHIC, it is twice as large as in SPS. This is a consequence of much higher collision energy at RHIC. Indeed, the total energy of the fluid is 25290 GeV at RHIC, which is 99% of total collision energy. Hence, higher collision energy does not lead to higher energy density but is used to produce the matter with large volume at $`\tau _0=`$1 fm/$`c`$.

The output from the fluids is summarized in Table 2. Total net baryon number of the fluid is much smaller in RHIC than in SPS, as well as the mean chemical potential on the freeze-out hypersurface. This difference can be seen in the space-time evolution of temperature on $`\eta \tau `$ plane (Fig. 14). As shown in Fig. 4, phase boundary can no longer be specified by temperature only but depends on both temperature and chemical potential in high net baryon density. For example, $`T=158`$ MeV corresponds to the hadronic phase at vanishing baryon density. However, it can be in QGP phase at $`\mu _\text{B}=400`$ MeV and be in mixed phase at some $`\mu _\text{B}`$. This behavior is seen $`T=158`$ MeV contour in Fig. 14, where the baryonic chemical potential becomes higher as $`\eta `$ increases and a fluid element near $`\eta 0.5`$ stays at mixed phase for a long time. Such behavior does not appear in the RHIC case where chemical potential is small enough to characterize the mixed phase by the almost constant temperature.

In Fig. 15 where deviation from Bjorken’s scaling solution $`Y_\text{L}\eta `$ is plotted, acceleration is larger in SPS than in RHIC because of the steeper pressure gradient of $`\eta `$-direction. Finally, the lifetime of each phase is also shown in Table 2.

## IV Two-particle correlation

In this section, we present the result of the two-pion correlation function and HBT radii based on the numerical solution of the relativistic hydrodynamical equation. For simplicity, we assume that all the pions are emitted from a chaotic source and neglect the resonance contribution. Then, the two-particle correlation function is easily calculated through

$$C_2(q^\mu ,K^\mu )=1+\frac{|I(q^\mu ,K^\mu )|^2}{I(0,k_1^\mu )I(0,k_2^\mu )},$$
(14)
where $`K^\mu =(k_1^\mu +k_2^\mu )/2`$, $`q^\mu =k_1^\mu k_2^\mu `$, respectively Shuryak (1973); Hama and Padula (1988). Here $`k_i^\mu `$ is on-shell mometum of $`i`$-th pion. We put

$$I(q^\mu ,K^\mu )=K_\tau 𝑑\sigma ^\tau (x)\sqrt{f(k_1,x)f(k_2,x)}e^{iq_\nu x^\nu },$$
(15)
so that $`I(0,k^\mu )`$ reduces to the Cooper-Frye formula with $`f(k,x)`$ being the Bose-Einstein distribution function. Considering the experimental momentum acceptance, we integrate the correlation function with respect to the average momentum in region $`\mathrm{\Omega }`$ as

$$C_2(q^\mu )|_\mathrm{\Omega }=1+\frac{_\mathrm{\Omega }K_\text{T}𝑑K_\text{T}𝑑Y|I(q^\mu ,K^\mu )|^2}{_\mathrm{\Omega }K_\text{T}𝑑K_\text{T}𝑑YI(0,k_1^\mu )I(0,k_2^\mu )}.$$
(16)
The HBT radii are obtained by fitting the calculated correlation function (16) to gaussian fitting function

$`C_{\text{2fit}}(q^\mu )`$ $`=1+\mathrm{exp}(R_{\text{side}}^2q_{\text{side}}^2R_{\text{out}}^2q_{\text{out}}^2R_{\text{long}}^2q_{\text{long}}^2R_{\text{ol}}^2q_{\text{out}}q_{\text{long}}).`$ (17)
For RHIC data, in which rapidity acceptance $`|Y|0.5`$, the out-long cross term $`R_{\text{ol}}`$ Chapman et al. (1995a) can be ignored. According to the azimuthal symmetry, we can put $`K_\text{T}=K_x`$, $`q_{\text{side}}=q_y`$ and $`q_{\text{out}}=q_x`$.

Results of $`K_\text{T}`$ dependence of the HBT radii are presented in Fig. 19 and 20, where we show the transverse momentum dependence HBT radii of the SPS Pb+Pb collisions and $`M_\text{T}\sqrt{K_\text{T}^2+m_\pi ^2}`$ dependence of the RHIC Au+Au collisions, respectively. In addition to the three radius parameters, we also present the ratio of $`R_{\text{out}}`$ to $`R_{\text{side}}`$ for better comparison between two collisions Rischke and Gyulassy (1996).

Sideward HBT radii (upper figures in Figs. 19 and 20) are consistent with the experiments in both RHIC and SPS. Larger radii than other calculations come from initial large transverse size of the fluid. Outward HBT radii show quantitative agreement for the SPS data. However, qualitative behavior shows some deviation from the experimental data; our result takes the maximum value at $`K_\text{T}0.3`$ GeV/$`c`$ while the experiment data seem to monotonically decrease except for the highest $`K_\text{T}`$ bin. For the RHIC data, experimental data show steep decrease with $`M_\text{T}`$. On the other hand, our results are similar to the ones of SPS because of the similarity in the space-time evolution of both fluids. As for the longitudinal HBT radii, our model reproduces the qualitative behavior of the results of both experiments but shows a little overestimate at low $`M_\text{T}`$ of RHIC result. Our result suggests that the longitudinal finite size effect is essential for understanding the behavior of $`R_{\text{long}}`$ even at RHIC because other calculations assuming infinite boost invariant source show larger $`R_{\text{long}}`$ Zschiesche et al. ; Heinz and Kolb . HBT radius in the longitudinal direction has two kinds of origin; spatial extent of the fluid Chapman et al. (1995a) and thermal suppression caused by rapid expansion. The emission region is roughly characterized by a product of the Boltzmann factor $`\mathrm{exp}[m_\text{T}\mathrm{cosh}(Y_\text{L}Y)/T]`$ and a geometrical factor (e.g., $`\mathrm{exp}(\eta ^2/2(\mathrm{\Delta }\eta )^2)`$ ). The deviation from the scaling solution causes stronger thermal suppression. As a result, our solution provides smaller longitudinal HBT radius. Our $`R_{\text{out}}/R_{\text{side}}`$ moderately increases with $`K_\text{T}`$ in both SPS and RHIC. This tendency is also seen in a hybrid calculation Soff et al. (2001) in spite of the quite different description of the hadronic phase. Our results show very good agreement with the SPS result, while the RHIC data clearly shows different behavior. As concerns the behaviors of $`R_{\text{out}}`$ and $`R_{\text{out}}/R_{\text{side}}`$, opaque source Heiselberg and Vischer (1998) is a possible explanation if the current formalism of the two-particle correlation is correct. Though the opaque property appears in the hydrodynamical model by virtue of transverse flow Morita et al. (2000), it is still insufficient to reproduce the small $`R_{\text{out}}`$ in the RHIC experiment. More theoretical investigation will be required for solving the problem Humanic ; McLerran\_nucl0205028 .

Figure 21 shows the correlation functions projected onto each components of relative momenta. Transverse momentum of emitted pair is $`0.125K_\text{T}0.225`$ GeV/$`c`$, which corresponds to the lowest momentum bin in Fig. 20. In each figure, calculated correlation functions are corrected by a common factor as $`\lambda =0.6`$. The other kind of the reduction of the correlation function at $`q_i=0`$ ($`i=`$side, out, long) is caused by an integration with respect to other components of the relative momenta. The resultant factor is in proportional to $`1/R_i`$. Our correlation function seems to be consistent with the experimental data for the outward and the longitudinal direction in spite of the overestimation of the HBT radii. Because the HBT radii correspond to the inverse width of the correlation function, difference in width between $`R=7`$ fm and $`R=6`$ fm is only about 5 MeV in relative momentum. The small value of sideward correlation function at $`q_{\text{side}}=0`$ indicates that outward and longitudinal HBT radii are larger than the experimental results.

## V Concluding Remarks

In this paper, we investigate single-particle distributions and two-particle correlation functions in SPS Pb+Pb 17 $`A`$GeV collisions and RHIC Au+Au 130 $`A`$GeV collisions based on a hydrodynamical model in which both longitudinal and transverse expansion are taken into account. As long as the single-particle spectra, the hydrodynamical model well describes both SPS and RHIC data. The initial parameter set in the model for both collisions indicates that initial energy density in RHIC is only slightly higher than the one in SPS and much larger extent of hot matter is produced in RHIC, if we compare them at the same initial time and by similar initialization (Fig. 3 and Eqs. (11),(12)). We have also discussed the space-time evolution of the fluids. Since the initial conditions are not different so much, temperature and transverse flow evolution do not show significant difference between SPS and RHIC. Only the equi-temperature contour of $`T=158`$ MeV shows the qualitative difference between SPS and RHIC due to the difference in net baryon number. Steeper pressure gradient in the longitudinal direction at SPS makes the deviation from the scaling solution larger than at RHIC. Two-pion correlation functions and the HBT radii are also investigated. Our model shows good agreement with the SPS data. For RHIC data, however, the outward and longitudinal HBT radii of our result are larger than the experimental data even though the dynamical longitudinal expansion and finite size effect are incorporated.

###### Acknowledgements.

The authors would like to thank Prof. I. Ohba and Prof. H. Nakazato for their fruitful discussions and comments. They also would like to acknowledge J. Alam, R.S. Bhalerao, P. Kolb, L. McLerran, H. Nakamura, and F. Retrière for valuable discussions and comments. This work was in part supported by the Ministry of Education, Science and Culture, Japan (Grant No. 11440080 and 13135221), Waseda University Grant for Special Research Projects No.2001A-888, Waseda University Media Network Center and ERI of Tokuyama University.


text:
# Decomposition of multicomponent mass spectra using Bayesian probability theory

## I Introduction

Mass spectrometry is a standard technique for residual gas analysis and active processing control in vacuum devices. Traditional low resolution quadrupole mass spectrometers are widely used due to high sensitivity, reasonable stability, wide operational pressure range, high scan speed and low costs. To be filtered in the quadrupole field, neutral gases have to be ionized first, most commonly by electron impact. At a typical electron energy of 50-100 eV used to achieve a high ionization efficiency and stability analyte molecules can decompose into a variety of fragment ions leading to the so called cracking pattern(CP). The fragmentation is a molecular specific property which reflects the chemical structure and the atomic composition of the molecule. It can therefore be used for the decomposition of multicomponent mass spectra.

A direct decomposition of mass spectra from mixture gases by successive subtraction of component contributions from the mixture signal, the wide spread pedestrian approach, works only if there are non-interfering mass numbers. A more elaborate method, least square evaluation, is not subject to this limitation and can incorporate measurement errors into the analysisdobrozemsky72 ; dobrozemsky92 . In both cases one must assume the CP of the components to be exactly known dobrozemsky72 . The CP is, however, not an intrinsic property of molecules but depends on the particular mass spectrometer and operating parameters making the determination of CP nontrivial. The problem becomes even more severe if unstable species like radicals are among the components of a mixture, since in this case there are neither calibration measurements nor literature values for the CP available. In case of mixtures containing a single radical the CP of the radical is usually identified with the rest signal after subtraction of all component signals from the mixture signal. Because of the error propagation of measured signals the above procedure, however, provides only poor and sometimes unphysical (e.g. negative cracking coefficients) estimations.

To overcome these difficulties and to treat measured spectra consistently we introduced a novel method using Bayesian probability theory schwarz01 . This was applied to the pyrolysis of azo-methane in order to evaluate concentrations and improved cracking coefficients together with their confidence intervals for the pyrolysis products including the methyl radical. In the present work we deal with another case of decomposition, in which mixtures made up given different compositions have to be decomposed without using any calibration measurements. Because of lack of calibration measurements this method offers a practical use, in particular, in handling unknown mixture gases. But most importantly, this approach is able to treat mixtures containing unstable species, because exact CPs are not needed.

As an example we synthesize mixtures of three different hydrocarbon gases, ethane(C2H6), propane(C3H8) and n-butane(C4H10). The gas composition and the cracking coefficients of the mixture components will be estimated by independent experiments and compared with the results from the Bayesian analysis.

## II Bayesian data analysis

The procedure used in this work is similar to that used for the decomposition of the pyrolysis products of azo-methane introduced by Schwarz-Selinger et al. schwarz01 , where details of the calculation may be found. Here we present only the relevant steps for our algorithm.

Assuming linear response of the mass spectrometer the mass signal $`\stackrel{~}{𝐝}_j`$ of a mixture is the sum of contributions of all species in the mixture

$$\stackrel{~}{𝐝}_j=\stackrel{~}{𝐂}\stackrel{~}{𝐱}_j+𝜺_j,$$
(1)
where $`j=1,\mathrm{},J`$ with $`J`$ mixtures. The vector $`\stackrel{~}{𝐝}_j`$ has $`N`$ elements representing the same number of mass channels. The composition vector $`\stackrel{~}{𝐱}_j`$ has the dimension $`M`$ of the number of species in the mixture. $`\stackrel{~}{𝐂}`$ is a cracking matrix with $`M`$ column vectors corresponding to the cracking patterns of $`M`$ species. Since we are interested in compositions and cracking patterns but not in the absolute intensities, $`\stackrel{~}{𝐝}_j`$ and the column vectors of $`\stackrel{~}{𝐂}`$ are normalized with respect to the sum of intensities of all mass channels for practical reasons, which enforces in turn the normalization of $`\stackrel{~}{𝐱}_j`$

$$\underset{n=1}{\overset{N}{}}d_{nj}=\underset{n=1}{\overset{N}{}}\underset{m=1}{\overset{M}{}}\stackrel{~}{c}_{nm}\stackrel{~}{x}_{mj}=\underset{m=1}{\overset{M}{}}\stackrel{~}{x}_{mj}\underset{n=1}{\overset{N}{}}\stackrel{~}{c}_{nm}=\underset{m=1}{\overset{M}{}}\stackrel{~}{x}_{mj}=1,$$
(2)
where $`\stackrel{~}{x}_{mj}`$ is the m-th component of the j-th mixture. This sum rule has been overlooked in our previous work schwarz01 . The results in Ref. 3 approximate however closely the sum rule. Quite the same it does improve stability of the estimates on $`\stackrel{~}{𝐱}`$ and $`\stackrel{~}{𝐂}`$ if the relation is explicitly included in the analysis. The vector $`𝜺_j`$ in Eq. (1) represents the finite noise of the experiments. We assume that $`𝜺_j`$ is bias free, $`𝜺_j=0`$, and has a variance of $`\epsilon _{nj}^2=\sigma _{nj}^2`$. Note that for the absolute composition of the mixtures $`\stackrel{~}{𝐱}_j`$ in (1) has to be scaled by a sensitivity factor of the mass spectrometer for molecules considered, which can be done by independent calibration experiments. Since it is beyond the scope of this work it will not be discussed. In the following we call $`\stackrel{~}{𝐱}_j`$ concentration of a mixture.

Now we can make use of the normalization condition (2) and reduce the number of unknown concentrations by 1. We split $`\stackrel{~}{𝐂}\stackrel{~}{𝐱}_j`$ in Eq. (1) into two parts

$$\stackrel{~}{𝐂}\stackrel{~}{𝐱}_j=\underset{m=1}{\overset{M}{}}\stackrel{~}{𝐜}_m\stackrel{~}{x}_{mj}=\underset{m=1}{\overset{M1}{}}\stackrel{~}{𝐜}_m\stackrel{~}{x}_{mj}+\stackrel{~}{𝐜}_M\stackrel{~}{x}_{Mj}.$$
(3)
Using $`\stackrel{~}{x}_{Mj}=1_{m=1}^{M1}\stackrel{~}{x}_{mj}`$, $`\stackrel{~}{𝐂}\stackrel{~}{𝐱}_j`$ becomes

$$\stackrel{~}{𝐂}\stackrel{~}{𝐱}_j=\stackrel{~}{𝐜}_M+\underset{m=1}{\overset{M1}{}}(\stackrel{~}{𝐜}_m\stackrel{~}{𝐜}_M)\stackrel{~}{x}_{mj}.$$
(4)
Then Eq. (1) can be rewritten

$$\stackrel{~}{𝐝}_j\stackrel{~}{𝐜}_M=\underset{m=1}{\overset{M1}{}}(\stackrel{~}{𝐜}_m\stackrel{~}{𝐜}_M)\stackrel{~}{x}_{mj}+𝜺_j.$$
(5)
We define

$`𝐝_j`$ $`=`$ $`\stackrel{~}{𝐝}_j\stackrel{~}{𝐜}_M`$
$`𝐂`$ $`=`$ $`(\stackrel{~}{𝐜}_1\stackrel{~}{𝐜}_M,\stackrel{~}{𝐜}_2\stackrel{~}{𝐜}_M,\mathrm{},\stackrel{~}{𝐜}_{M1}\stackrel{~}{𝐜}_M)`$ (6)
$`𝐱_j^T`$ $`=`$ $`(\stackrel{~}{x}_{1j},\stackrel{~}{x}_{2j},\mathrm{},\stackrel{~}{x}_{(M1)j}),`$
and obtain the simplified model equation

$$𝐝_j=\mathrm{𝐂𝐱}_j+𝜺_j.$$
(7)
Knowing the model equation for the mass signal, which assumes $`𝐝_j\mathrm{𝐂𝐱}_j=0`$ and $`\left(𝐝_j\mathrm{𝐂𝐱}_j\right)_n^2=\sigma _{nj}^2`$, the principle of maximum entropy sivia96  leads to a Gaussian sampling distribution(likelihood) for the data set of the $`j`$-th mixture

$$p(𝐝_j|𝐱_j,𝐂,𝐒_j,I)=\frac{1}{_ns_{nj}\sqrt{2\pi }}\mathrm{exp}\left\{\frac{1}{2}\left(𝐝_j\mathrm{𝐂𝐱}_j\right)^T𝐒_j^2\left(𝐝_j\mathrm{𝐂𝐱}_j\right)\right\},$$
(8)
where $`\left(𝐒_j^2\right)_{nn}=1/s_{nj}^2`$ and $`\left(𝐒_j^2\right)_{nl}=0`$ for $`nl`$, with $`s_{nj}`$ being the n-th component of the sample variance $`𝐬_j`$ used for $`\sigma _{nj}`$. $`I`$ denotes all background information available. We assume independent measurements. This means that the outcome of an experiment is not influenced by any other data set. In such a case the likelihood for all data sets reduces to a product of the likelihoods of the individual measurements

$$p\left(𝐃|𝐗,𝐂,𝐒,I\right)=\underset{j}{}p\left(𝐝_j|𝐱_j,𝐂,𝐒_j,I\right).$$
(9)
This can easily be shown by repeated application of the product rule of the Bayesian probability theory $`p(A,B)=p(A)p(B|A)`$. For reasons of convenience we used the notation $`\left\{𝐝_j\right\}𝐃`$, $`\left\{𝐱_j\right\}𝐗`$ and $`\left\{𝐒_j\right\}𝐒`$, where $`𝐃`$ and $`𝐗`$ are matrixes with $`J`$ column vectors $`𝐝_j`$ and $`𝐱_j`$, respectively.

Our task is to find the concentration $`𝐱_j`$ ($`j=1,\mathrm{},J`$) and the cracking coefficients $`\stackrel{~}{𝐜}_m`$ ($`m=1,\mathrm{},M`$). Let us start with $`𝐱_j`$. The expectation value of $`𝐱_k`$ is derived from the probability density $`p(𝐗|𝐃,𝐒,I)`$ as

$$𝐱_k=𝑑𝐗𝐱_kp(𝐗|𝐃,𝐒,I).$$
(10)
$`p(𝐗|𝐃,𝐒,I)`$ is obtained from Bayes theorem as sivia96

$$p(𝐗|𝐃,𝐒,I)=p(𝐗|I)p(𝐃|𝐗,𝐒,I)/Norm,$$
(11)
in terms of the marginal likelihood $`p(𝐃|𝐗,𝐒,I)`$. The latter is derived by the marginalization rule

$$p(𝐃|𝐗,𝐒,I)=𝑑\stackrel{~}{𝐂}p(𝐃,\stackrel{~}{𝐂}|𝐗,𝐒,I)=𝑑\stackrel{~}{𝐂}p(\stackrel{~}{𝐂}|I)p(𝐃|\stackrel{~}{𝐂},𝐗,𝐒,I).$$
(12)
Note that $`p(\stackrel{~}{𝐂}|𝐗,𝐒,I)=p(\stackrel{~}{𝐂}|I)`$, since the composition and the sample variance do not contain any information about the cracking coefficients. We use a flat prior for $`p(𝐗|I)`$ in the range of $`0<x_{mj}<x_{max,j}`$, since there exists no prior knowledge about the composition, and obtain

$$𝐱_k𝑑\stackrel{~}{𝐂}p(\stackrel{~}{𝐂}|I)𝑑𝐗𝐱_kp(𝐃|𝐗,\stackrel{~}{𝐂},𝐒,I),$$
(13)
where the proportionality accounts for the normalization factor in Eq. (11) and the flat prior $`p(𝐗|I)`$. For the analytical integration of the inner integral in (13) we assume that the integrand is sharply peaked in the prior range for $`𝐱_j`$ so that the integral does not change if the integral limits are extended to $`(\mathrm{},\mathrm{})`$. The likelihood $`p(𝐃|𝐗,\stackrel{~}{𝐂},𝐒,I)`$ in Eq. (13) is a multivariate Gaussian function (8,9), whose exponent can be transformed into a complete square

$$\left(𝐝_j\mathrm{𝐂𝐱}_j\right)^T𝐒_j^2\left(𝐝_j\mathrm{𝐂𝐱}_j\right)=\left(𝐱_j𝐱_{j0}\right)^T𝐐_j\left(𝐱_j𝐱_{j0}\right)+𝐑_j.$$
(14)
The integral of the right hand side of Eq. (14) over $`𝐱_j`$ with range $`(\mathrm{},\mathrm{})`$ returns

$$\mathrm{exp}\left(\frac{R_j}{2}\right)|det𝐐_j|^{\frac{1}{2}},$$
(15)
and thus for the rightmost integral of (13) we obtain

$$𝑑𝐗𝐱_kp(𝐃|𝐗,\stackrel{~}{𝐂},𝐒,I)=𝐱_{k0}\underset{j}{}\mathrm{exp}(\frac{R_j}{2})|det𝐐_j|^{\frac{1}{2}},$$
(16)
where $`𝐱_{k0}`$, $`R_j`$ and $`𝐐_j`$ result from the comparison of the coefficients in Eq. (14)

$`𝐐_j`$ $`=`$ $`𝐂_j^T𝐒_j^2𝐂_j`$
$`𝐱_{k0}`$ $`=`$ $`𝐐_k^1𝐂_k^T𝐒_k^2𝐝_k`$ (17)
$`R_j`$ $`=`$ $`𝐝_j^T𝐒_j^2𝐝_j𝐱_{j0}^T𝐐_j𝐱_{j0}.`$
The calculation of the matrix equations can be considerably simplified by using singular value decomposition (SVD) of $`𝐂_j`$. For details readers are referred to Ref. 3 and 5.

To complete the integration (13) the prior probability $`p(\stackrel{~}{𝐂}|I)`$ for the cracking coefficients has to be assigned. As a prior knowledge $`I`$ we use the table of Cornu and Massot cornu79 , which provides only point estimates for CP. Then, by virtue of the maximum entropy principle sivia96 , the prior probability distribution of $`\stackrel{~}{𝐂}`$ given the tabulated value $`c_0`$ is an exponential function

$$p(\stackrel{~}{c}|I)=p(\stackrel{~}{c}|c_0)=\mathrm{exp}(\lambda \stackrel{~}{c})/Z,$$
(18)
where $`c_0`$ is the literature value for the matrix element $`\stackrel{~}{c}`$. Eq. (18) is a general form of the exponential function having the normalization constant $`1/Z`$ and the scale factor $`\lambda `$. The latter can be obtained from the requirement that the expectation value of $`\stackrel{~}{c}`$ on the support $`0<\stackrel{~}{c}<1`$ must be equal to $`c_0`$

$$\stackrel{~}{c}=_0^1𝑑\stackrel{~}{c}\stackrel{~}{c}\mathrm{exp}(\lambda \stackrel{~}{c})/Z=\frac{1}{\lambda }\frac{1(1+\lambda )\mathrm{exp}(\lambda )}{1\mathrm{exp}(\lambda )}=c_0.$$
(19)
For $`p(\stackrel{~}{𝐂}|𝐂_0)`$ the prior probabilities $`p(\stackrel{~}{c}|c_0)`$ of all elements of the matrix $`\stackrel{~}{𝐂}`$ have to be multiplied

$$p(\stackrel{~}{𝐂}|I)=p(\stackrel{~}{𝐂}|𝐂_0)=\underset{nm}{}p(\stackrel{~}{c}_{nm}|c_{0,nm}).$$
(20)
Inserting these results into Eq. (13) we have

$$𝐱_k𝑑\stackrel{~}{𝐂}𝐱_{k0}p(\stackrel{~}{𝐂}|I)\underset{j}{}\mathrm{exp}(\frac{R_j}{2})|det\text{Q}_j|^{\frac{1}{2}}.$$
(21)
The integration over $`\stackrel{~}{𝐂}`$ is performed by Markov Chain Monte Carlo(MCMC) methodgilks96  with the sampling density $`\rho (\stackrel{~}{𝐂})`$

$$\rho (\stackrel{~}{𝐂})p(\stackrel{~}{𝐂}|I)\underset{j}{}\mathrm{exp}(\frac{R_j}{2})|det\text{Q}_j|^{\frac{1}{2}}.$$
(22)
For the second moment of $`𝐱_k`$ we find

$$x_{mk}^2=𝑑\stackrel{~}{𝐂}\left[\left(𝐐_k^1\right)_{mm}+x_{mk0}^2\right]\rho (\stackrel{~}{𝐂}),$$
(23)
which allows to calculate the variance of $`𝐱_k`$

$$\mathrm{\Delta }x_{mk}^2=x_{mk}^2x_{mk}^2.$$
(24)
Now we have to assign the expectation value and the variance for $`\stackrel{~}{x}_{Mk}`$. The calculation for $`\stackrel{~}{x}_{Mk}`$ is straightforward

$$\stackrel{~}{x}_{Mk}=1\underset{m=1}{\overset{M1}{}}\stackrel{~}{x}_{mk}=1\underset{m=1}{\overset{M1}{}}\stackrel{~}{x}_{mk}.$$
(25)
Accordingly, for $`\stackrel{~}{x}_{Mk}^2`$ we have

$`\stackrel{~}{x}_{Mk}^2`$ $`=`$ $`(1{\displaystyle \underset{m=1}{\overset{M1}{}}}\stackrel{~}{x}_{mk})(1{\displaystyle \underset{l=1}{\overset{M1}{}}}\stackrel{~}{x}_{lk})`$ (26)
$`=`$ $`12{\displaystyle \underset{m=1}{\overset{M1}{}}}\stackrel{~}{x}_{mk}+{\displaystyle \underset{m,l}{\overset{M1}{}}}\stackrel{~}{x}_{mk}\stackrel{~}{x}_{lk}.`$
From Eq. (25) and (26) it follows

$`\mathrm{\Delta }\stackrel{~}{x}_{Mk}^2`$ $`=`$ $`\stackrel{~}{x}_{Mk}^2\stackrel{~}{x}_{Mk}^2`$
$`=`$ $`12{\displaystyle \underset{m=1}{\overset{M1}{}}}\stackrel{~}{x}_{mk}+{\displaystyle \underset{m,l}{\overset{M1}{}}}\stackrel{~}{x}_{mk}\stackrel{~}{x}_{lk}`$
$`\left\{12{\displaystyle \underset{m=1}{\overset{M1}{}}}\stackrel{~}{x}_{mk}+{\displaystyle \underset{m,l}{\overset{M1}{}}}\stackrel{~}{x}_{mk}\stackrel{~}{x}_{lk}\right\}`$
$`=`$ $`{\displaystyle \underset{m,l}{\overset{M1}{}}}\left\{\stackrel{~}{x}_{mk}\stackrel{~}{x}_{lk}\stackrel{~}{x}_{mk}\stackrel{~}{x}_{lk}\right\}.`$
The estimation of cracking coefficients proceeds in a similar way. First, Bayes theorem was applied to rewrite the required posterior probability $`p(\stackrel{~}{𝐂}|𝐃,𝐒,I)`$ in terms of a prior on $`\stackrel{~}{𝐂}`$ and the marginal likelihood $`p(𝐃|\stackrel{~}{𝐂},𝐒,I)`$

$$p(\stackrel{~}{𝐂}|𝐃,𝐒,I)p(\stackrel{~}{𝐂}|I)p(𝐃|\stackrel{~}{𝐂},𝐒,I),$$
(28)
and by using the marginalization rule and the product rule, the likelihood $`p(𝐃|\stackrel{~}{𝐂},𝐒,I)`$ becomes

$`p(𝐃|\stackrel{~}{𝐂},𝐒,I)`$ $`=`$ $`{\displaystyle 𝑑𝐗p(𝐃,𝐗|\stackrel{~}{𝐂},𝐒,I)}`$ (29)
$`=`$ $`{\displaystyle 𝑑𝐗p(𝐗|I)p(𝐃|𝐗,\stackrel{~}{𝐂},𝐒,I)}.`$ (30)
Finally, we arrive at

$$\stackrel{~}{𝐜}_m^\mu =𝑑\stackrel{~}{𝐂}\stackrel{~}{𝐜}_m^\mu \rho (\stackrel{~}{𝐂}),$$
(31)
where $`\mu =1,2`$ for the first and second moment for $`\stackrel{~}{𝐜}_m`$. The sampling density $`\rho (\stackrel{~}{𝐂})`$ is the same as in Eq. (22). As a consequence $`\stackrel{~}{𝐂}`$ and $`\stackrel{~}{𝐗}`$ can be determined by a same MCMC sampling.

## III Results and Discussion

All mass spectra presented in this work were taken with a quadrupole mass spectrometer (Hiden, HAL 201), which has the standard axial ion source and the secondary electron multiplier (SEM) for ion detection. The mass spectrometer chamber with a base pressure better than $`110^9`$hPa is differentially pumped and connected with the main chamber by a capillary ($`\varphi `$ 1 mm, 10 mm long). Both chambers are pumped with a separate turbo molecular pump. The gas flows into the main chamber are controlled by three individual MKS flow controllers (type 1259C), which allow a nominal error of 0.8 %. During calibration as well as mixture measurements the pressure in the mass spectrometer chamber was kept at about $`510^8`$hPa to minimize any pressure effect on the cracking patterns.

Fig. 1 depicts the mass signals of 6 mixtures of different gas compositions composed of ethane, propane and n-butane for the most intensive 14 mass channels. Note the strong similarity of the six spectra. The quoted ratio was determined from the gas flows into the main chamber and does not reflect the actual concentrations, since the gas flow through the capillary as well as the pumping power in the main and the mass spectrometer chamber strongly depend on species. However, they still may serve as a crude estimate of the mixture ratio.

The upper part in Fig. 2a shows calibration measurements of C2H6, C3H8 and C4H10 together with their literature values taken from Cornu and Massot cornu79 . The cracking coefficients are normalized with respect to the sum of the intensities, so that they all have the total intensity of one making the direct comparison of cracking patterns from different sources possible. Although the rough structures of CP from literature and calibration measurements for the three gases look very similar, nearly every peak shows a clear deviation between calibration and literature values. This is actually expected since the CP depends on the geometry of ion sources, analyzer transmission and mass dependent secondary electron multiplier gain reflecting the specification of a particular mass spectrometer. The lower plots in Fig. 2a show the results of the Bayesian analysis on the mixture data in Fig. 1. The agreement of CP between calibration and our analysis is very good. For ethane the Bayesian results coincide exactly with those from the calibration, and for propane and butane only few points show a small deviation, namely peaks at m=28 and m=29, 43, respectively. But this deviation means not necessarily an error in the measurements or analysis. In contrast, it is conceivable that CP of pure gases is changed in a mixture by an influence of gases on each other through a simultaneous presence in the ion source or by variation of partial pressures breth82 .

Fig. 2b shows the concentrations of the components in six mixtures from the Bayesian analysis and from the independent estimation by pressure measurement. For the latter the partial pressure of known gas flows of pure gases was measured by a baratron gauge, which is mounted in the main chamber, and the flux rate of a mixture was rescaled to partial pressures of each gas. The error of this estimation was assessed to 10 %, which is mainly contributed by the long-term fluctuation of the total pressure and the tolerance of the flow controller. Our results in Fig. 2b are in a satisfactory agreement with the estimation by pressure measurement and rather accurate with about 3 % confidence interval. It can also be found that there is a systematic deviation between the two sources: the concentrations from the pressure estimation for ethane are lower than those from our analysis and the concentrations for butane show exactly the opposite. This may be explained by the mass dependent conductance through the capillary reif65 , which connects the mass spectrometer and the main chamber with the gas inlet system. Since the conductance of the capillary is larger for ethane than for butane, the estimations by the baratron have to be increased for ethane and decreased for butane, which leads to a better agreement between the two determinations. In fact, the determination of the composition of a mixture is rather difficult and needs much care. A gaseous molecule in a mixture can affect the molecular kinetic of other gases and thus the pressure reading of pure gases cannot be directly combined with the mixture composition. If no calibration of pure gases is available as in the case of radicals, the pressure determination is even impossible. Our method in this case provides a unique way to yield concentrations of all components in the mixture.

The accuracy of Bayesian analysis relies on the amount of information provided, which include prior knowledge, number of measurements or the availability of calibration measurements. It is therefore interesting and important to ask how much information or how many data sets we need to obtain reliable results from the equation system (7). A single mixture delivers 13 known and 2 unknown parameters ($`141`$ for $`𝐝_j`$ and $`31`$ for $`𝐱_j`$ due to the normalization), while the number of unknown cracking coefficients (=30) is unchanged. For the latter mass numbers with no significant signal compared to the background noise are neglected. For the analysis with 3 mixture data, for example, the number of known and unknown parameters are 39 and 36, respectively, illustrating an overdetermined equation system. This consideration, however, is only valid if the input data are linearly independent. Our data by no means meet this condition, as can be seen qualitatively in Fig. 1. Quantitatively, the significant singular values of the SVD of the data matrix gives the number of sources, which effectively contribute to solving the equation system. For the data matrix with six column vectors only three singular values are significant (Table 1) indicating that at least six data sets would be needed for the reliable decomposition of the mixtures by the purely mathematical treatment. In contrast, for the Bayesian approach each data set contributes to increase the accuracy of the analysis on the basis of prior knowledge. As can be seen in Fig. 2, six mixtures were proved to contain enough information to decompose the mixtures. Let us now reduce the input data down to three mixtures. In Fig. 3 the results from three data sets are compared to those from six data sets. The point estimations of the cracking coefficients are nearly unchanged, but the error bars became somewhat broader, which seems to be reasonable because of less information. The concentrations derived from 3 and 6 mixtures also show no significant change. These very satisfactory outcomes from three mixtures would not be possible by directly resolving the equation system (7). It should be also noted that the analysis with two mixture data illustrating a superior number of unknown parameters leads to a strong uncertainty in the estimation of both cracking coefficients and concentrations.

We can also even increase the accuracy by adding calibration measurements to our analysis. Fig. 4 shows the results deduced from 7 data sets including the 6 mixture data and the calibration measurement for butane. One can see that the deviation of cracking coefficients at m=29, 43 for butane shown in Fig. 2a entirely vanishes and that the overall error significantly decreases. This demonstrates that the incorporation of calibration measurements is an efficient way to improve the accuracy of our analysis. This may be possible even for very unknown mixtures, since in most cases one can infer from the source of mixture gases one or more of its components.

## IV Summary

We introduced a method for decomposition of multicomponent mass spectra using Bayesian probability theory, which was used to derive concentrations and cracking coefficients of mixture components as well as their confidence interval without any use of calibration measurements. For synthesized mixtures from ethane, propane and butane the algorithm was able to deliver precise results, which were compared with independent estimations. This example exhibits a very challenging case, in which the mass spectra of the mixtures show only a moderate difference from each other. Since the number of available mixture data does not restrict the application of our method but merely affects the accuracy of the analysis, even cases in which limited number of mixtures is available can be processed by our Bayesian approach. For unknown mixtures, if only few components are assumed to be known from the source of the mixtures, the outcome of our analysis for the concentrations serves as a critical criterion to decide whether a certain species is among the constituents of the mixture.

Future applications include the decomposition of neutral particle fluxes from low temperature hydrocarbon plasmas including larger radicals employed in thin-film deposition and etching. The radicals are believed to play an crucial role in film formation at the plasma-surface boundary jacob98 ; sugai90 . Their quantification, which will be possible by the present method, is therefore prerequisite for understanding of the microscopic film growth process. In this case mixture data are easily available, e.g. by variation of plasma parameters jacob98 .

References

Table 1:Eigenvalues of the singular value decomposition (SVD) of the data matrix with six column vectors shown in Fig. 1.

Fig.1: Mass signals of six different gas mixtures composed of ethane(C2H6), propane(C3H8) and n-butane(C4H10). The mixture rates given for C2H6:C3H8:C4H10 were determined from gas flows through the gas-inlet device.

Fig.2: (a) upper: cracking patterns from calibration measurements and literature taken from Ref. 6. The mass signals were normalized to the sum of the intensities. lower: Calibration measurements vs. Bayesian analysis using the mixture data in Fig. 1. (b) component concentration of the 6 mixtures from pressure measurement(grey) and from the Bayesian analysis(black). E, P and B in x-axis stand for ethane, propane and butane, respectively. See text for error determination.

Fig.3: Cracking patterns (a) and concentrations (b) of mixture components from Bayesian analysis using 3(black) and 6(grey) mixtures.

Fig.4: Cracking patterns from calibration measurements and Bayesian analysis using 6 mixtures and the calibration measurement of butane.