text:
# Thermodynamics of Electrolytes on Anisotropic Lattices

## I Introduction

After controversial experimental results on the nature of ionic criticality pitzer ; japas ; weingartner01  the studies of critical phenomena in ionic fluids received a fresh impetus, and over the last decade substantial progress has been made weingartner01 ; levin02 . However, a complete thermodynamic description of electrolyte systems is still unavailable.

The simplest and most frequently employed model system for both theoretical and Monte Carlo investigations of Coulombic systems is the restricted primitive model (RPM), which is a gas of equisized hard-sphere ions carrying positive and negative charges of equal magnitude. From the theoretical point of view, a reasonable description of the RPM in the critical region has been obtained at a mean-field level using integral equations approach stell95 ; yeh96 ; blum02  and Debye-Hückel theory levin02 ; levin96 . While theoretical predictions of the critical parameters levin96 ; blum02 ; kulinskii  agree reasonably well with those obtained in computer simulations panagiotopoulos99 ; caillol02 , recent Monte Carlo investigations of ionic systems indicate a non-classical Ising-like asymptotic critical behavior lujiten02 . Moreover, predictions of different mean-field theories for charge- and size-asymmetric primitive systems differ significantly from each other zuckerman01 ; banerjee02 ; artyomov03 ; netz ; sabir . In the absence of real experimental data, numerical simulations provide an important criterion for testing different theoretical methods camp99 ; panagiotopoulos02 ; yan02 ; cheong03 .

In recent years, a new theoretical approach to study the thermodynamics and criticality of Coulomb systems based on the investigation of corresponding lattice models, has been introduced dickman99 ; kobelev02-1 ; brognara ; ciach01 ; kobelev02-2 ; artyomov03 . Being obviously less realistic than their continuous-space counterparts, lattice models are much easier to handle analytically, and the information they provide is crucial for understanding the physics of criticality. Moreover, finely discretized lattices with lattice spacing several times smaller than the ion size present substantial computational advantages over continuous space simulations, while producing almost identical critical parameters and phase diagrams panagiotopoulos99 . By contrast, when the ionic diameter is equal to the spacing of the underlying lattice, the phase diagram is drastically different. The symmetry of simple cubic or body-centered cubic lattices allows for unequal charge distribution between two identical sublattices, yielding an ordered structure similar to an ionic crystal. The competition between electrostatic interactions of the two charged sublattices and the entropy of charge separation leads then to order-disorder phase transitions and a tricritical point dickman99 ; panagiotopoulos99 ; kobelev02-1 ; brognara .

A systematic study of the lattice restricted primitive model (LRPM) based on the Debye-Hückel approach has shown that, in accordance with Monte Carlo simulations, a Coulombic-governed gas-liquid phase transitions seen in continuum electrolytes is totally suppressed on simple cubic and body-centered cubic lattices due to the formation of the thermodynamically more stable charge-ordered phase kobelev02-1 . Ciach and Stell ciach01  argued that, unlike non-ionic fluids, in charged systems the most important are short-range fluctuations. Consequently, they studied a LRPM supplemented with additional short-range interactions between ions and concluded that at some strength of the short-range attraction both the normal critical and tricritical point can become thermodynamically stable ciach01 . This leads to complex phase diagrams with two distinct phase transitions, one between gas and disordered liquid phase, and another one between disordered liquid and a charge ordered phase. Recent Monte Carlo studies have confirmed qualitatively most of these theoretical predictions but further investigations are required diehl03 .

At low dimensions Coulombic interactions are stronger than in three-dimensional systems, and the possibility for gas-liquid phase transitions increases. This idea has been used recently in the investigation of anisotropic lattice models of electrolytes using Debye-Hückel method kobelev02-2 . In this study the anisotropy mimicked the lowering of the spatial dimensionality. It has been found, that for strongly anisotropic lattices gas-liquid phase coexistence is restored. However, this theory neglects an ion clustering phenomena, which are important in charged particles systems levin96 ; levin02 . Thus these theoretical predictions should be tested by computer simulations.

Current Monte Carlo studies of lattice RPM in $`d=3`$ dimensions employ the continuous space $`1/r`$ potential panagiotopoulos99 ; dickman99 ; diehl03 . However, because of the discrete lattice symmetry, the correct lattice Coulomb potential differs from $`1/r`$ at short distances and approaches it asymptotically at large distances. This raises the question of how using the correct electrostatic interactions will affect the thermodynamics and critical behavior of ionic systems. Indeed, the lattice correction to $`1/r`$ potential decays rapidly with distance and thus it can be viewed in a sense as an additional short-range interaction. In order to understand correctly the underlying mechanisms of phase transitions in Coulombic systems, lattice models should be investigated with the corresponding lattice potentials.

In this paper we study the RPM on tetragonal lattices using grand canonical Monte Carlo simulations and histogram reweighting technique ferrenberg . The correct lattice Coulomb potential is used and an Ewald-like summation is utilized to account for the long-range nature of electrostatic interactions. The paper is organized in the following way. In Sec. II the details of the model and simulation method are given. The results and discussions are presented in Sec. III, and summary and conclusions are drawn in Sec. IV. Analytical expressions and the details of lattice potential calculations are given in Appendix.

## II Anisotropic Lattice Model and Simulations Method

The system we have studied consists of $`2N`$ ions, half of them carrying charge $`q`$ and half $`q`$, positioned on a simple tetragonal lattice with lattice parameters ratio (degree of anisotropy) $`\alpha =b/a`$. The pairwise interaction potential of two ions separated by distance $`r_{ij}`$ has the form (for details see Appendix A)

$$U_{ij}=\frac{1}{4\pi ^2Db}_\pi ^\pi \frac{e^{i\mathrm{𝐤𝐫}}d^3𝐤}{2+1/\alpha ^2(\mathrm{cos}k_1+\mathrm{cos}k_2+1/\alpha ^2\mathrm{cos}k_3)},$$
(1)
where $`D`$ is the dielectric constant of the structureless solvent. Reduced quantities are obtained by scaling with the energy of the strongest ion-ion interaction in two neighboring $`xy`$ planes in the continuum limit, $`E_0=q^2/Db`$, and the unit cell volume, $`v_0=a^2b`$, namely,

$$T^{}=\frac{k_BT}{E_0},\rho ^{}=\frac{2Na^2b}{V},$$
(2)
where $`V`$ is the volume of the system.

Since the Coulomb potential is long ranged, one needs to account for interactions with particles in all images of the periodic box used in the simulations. In the case of continuous-space $`1/r`$ interactions the standard method to achieve this is the Ewald summation (see e.g. Ref.frenkel ; deLeeuw ). However, no similar method is available for the lattice Coulomb potential. Therefore, we propose the following approximate scheme, which is based on the asymptotic properties of the anisotropic lattice potential (1).

It can be shown that the potential (1) for finite $`\alpha `$ at large distances behaves as

$$U_{ij}(𝐫=(x,y,z))\frac{\alpha }{\sqrt{x^2+y^2+\alpha z^2}}$$
(3)
For this potential, an Ewald-like sum can be constructed along the lines of the derivation in Ref.deLeeuw , which yields the “continuum” part of the energy of interactions of periodic system of $`2N`$ particles

$`E_1={\displaystyle \underset{allspace}{}}{\displaystyle \frac{1}{|𝐫_{ij}^{}|}}`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{2N}{}}}q_iq_j\left({\displaystyle \underset{𝐧^{}}{}}{\displaystyle \frac{\alpha ^2\mathrm{exp}(\pi ^2|𝐧^{}|^2/\eta ^2+2\pi i𝐧^{}𝐫_{ij}^{})}{2\pi |𝐧^{}|^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐧^{}}{}}{\displaystyle \frac{\mathrm{erfc}(\eta |𝐫_{ij}+𝐧^{}|)}{|𝐫_{ij}^{}+𝐧^{}|}}\right)`$ (4)
$`+`$ $`{\displaystyle \underset{|𝐧^{}|0}{}}\left({\displaystyle \frac{\mathrm{erfc}(\eta |𝐧^{}|)}{|𝐧^{}|}}+\alpha ^2{\displaystyle \frac{\mathrm{exp}(\pi ^2|𝐧^{}|^2/\eta ^2)}{\pi ^2|𝐧^{}|^2}}\right){\displaystyle \frac{\eta }{\sqrt{\pi }}}{\displaystyle \underset{i=1}{\overset{2N}{}}}q_i,`$
where $`r^{}=(x/\alpha ,y/\alpha ,z)`$, $`n^{}=(\alpha n_x,\alpha n_y,n_z)`$, $`0n_x,n_y,n_z2\pi `$, and $`\mathrm{erfc}(\mathrm{x})`$ is the complimentary error function. Here conducting boundary conditions have been utilized, which are less sensitive to finite-size effects panagiotopoulos02 . As usual, the real-space damping parameter $`\eta `$ must be chosen in a such way that all sums in (4) converge fast. In our simulations we used $`\eta =5`$, and the Fourier-space sums were restricted to 518 wave vectors.

The difference between the correct lattice potential (1) and its asymptotic limit (3) decays rapidly with distance: see Table 1. Therefore, to obtain the total energy of electrostatic interactions, the lattice correction is added for ions in elliptic shells up to certain distance,

$$E_{total}=E_1+\underset{|𝐫^{}/L|<n_{max}}{}\left(U_{ij}(r)\frac{1}{|𝐫_{ij}^{}|}\right).$$
(5)
We performed grand canonical Monte Carlo simulations on cubic boxes of length $`L`$ stretched in $`z`$ direction by a factor of $`\alpha `$, under periodic boundary conditions. Distance-biased algorithm was employed for insertions and removals of pairs of unlike ions at each time step in order to facilitate acceptance, following Ref.orkoulas94 . To analyze the simulation data and obtain the coexistence curves, multi-histogram reweighting techniques ferrenberg  have been used.

## III Results and Discussion

To reduce possible finite-size effects, we used a system size $`L^{}L_{x,y}/a=L_z/b=12`$ for the isotropic lattice ($`\alpha =1`$) and $`\alpha =1.2`$, and $`L^{}=16`$ for stronger anisotropies. Short initial runs were used to get a rough estimate of the location of the tricritical point. After that, typical runs involved $`13\times 10^7`$ Monte Carlo steps for equilibration and $`27\times 10^8`$ for production. We did not try to locate the second order transitions Néel line and restricted ourselves with obtaining approximate location of the tricritical point by simple linear extrapolation of the coexistence lines, as expected for $`d=3`$ tricriticality.

In order to separate the effect of using the correct lattice potential in simulations, we studied first the phase coexistence in LRPM on a simple cubic lattice: see Fig.1. On a simple cubic lattice, Coulomb gas with $`1/r`$ potential is known to phase separate, with two coexisting phases being a low-density disordered phase and an ordered high-density ionic crystal-like phase panagiotopoulos99 ; dickman99 . The tricritical parameters estimates for the system size $`L^{}=12`$, obtained by Panagiotopoulos and Kumar panagiotopoulos99 , are $`T_{tri}^{}=0.15\pm 0.01,\rho _{tri}^{}=0.48\pm 0.02`$ panagiotopoulos99 . For the same system size, taking into account the lattice correction to the potential yields $`T_{tri}^{}=0.22,\rho _{tri}^{}=0.48`$ with the same accuracy. Thus the tricritical density remains the same within the error limits, but the temperature increases by the factor of $`1.5`$. Qualitatively, the increase in the tricritical temperature should be expected since lattice Coulomb interactions are stronger than $`1/r`$ at short distances, which leads to a higher stability of dense phases. Nevertheless, the tricritical point is still significantly lower (by 50%) than the predictions of the Debye-Hückel theory of lattice electrolytes $`T_{tri}0.365`$ kobelev02-1 , but is about 20% higher than the results $`T_{tri}0.202`$ from the hierarchical reference theory (HRT) of LRPM by Brognara et al.brognara , who also used the correct form of potential in their calculations. Note, that both of these theories agree on the tricritical density to be $`\rho _{tri}^{}0.38`$, which is about 20% lower than our simulation results. However, while taking into account Bjerrum pairing of ions into neutral dipoles bjerrum  in Debye-Hückel approach would supposedly produce better results by shifting the phase coexistence to lower temperatures and higher densities, the predictions of the HRT seem to be final.

The phase coexistence for different degrees of anisotropy $`\alpha `$ is shown in Fig.1. At any $`\alpha `$, the system shows only order-disorder phase separations and a tricritical point. Thus, the effect of stretching the lattice amounts merely to altering the tricritical point location, without changing qualitatively the topology of the phase diagram. This goes in contrast with the recent theory based on Debye-Hückel calculations kobelev02-2  which predicts that at the lattice parameters ratio $`b/a2.6`$ a distinct gas-liquid phase coexistence may reappear. However, as was mentioned above, the theoretical approach of Ref.kobelev02-2  disregards the Bjerrum pairing bjerrum , or, more generally, clustering of ions, which is crucial for the thermodynamics of the gas-liquid transition in electrolytes. Typically, taking into account pairing somewhat diminishes the critical temperature but significantly increases the critical density levin96 ; levin02 ; kobelev02-1 . Nevertheless, taking into consideration the ion clustering is unlikely to change the topology of the phase diagram of LRPM on the isotropic lattice. However, since the strength of the electrostatic interactions increases with anisotropy, this may lead to more profound tendency of ions to bind into neutral dipoles/clusters (mostly, in $`xy`$-planes), producing eventually a Kosterlitz-Thouless (KT) transitions kosterlitz  line at infinite stretching. Stronger pairing then could shift the gas-liquid coexistence curve much further than in the isotropic case, again putting it completely inside the order-disorder phase envelope. It should be noted that simulations of the two-dimensional lattice Coulomb gas teitel97  also reveal no normal first order gas-liquid transition, but only an order-disorder phase separations teitel97 .

Tricritical temperature and density as functions of the lattice parameters ratio $`\alpha `$ are presented in Fig.2. As anisotropy increases, both tricritical temperature and density increase. Physical account of this behavior follows from the properties of the lattice Coulomb potential. Indeed, in $`z`$-direction lattice potential (1) behaves almost as $`1/z`$ (see Table 1). However, in crystallographic directions different from $`[001]`$, the anisotropic lattice potential is much stronger than $`1/r`$ and decays slower. While the interactions in $`z`$-direction almost do not change with $`\alpha `$, stronger anisotropy brings stronger interactions inside the $`xy`$-layers. This becomes especially important at large densities, where the average distance between ions is of the order of the lattice spacing. Therefore, the average energy per particle increases with anisotropy. As a consequence, higher temperatures are required for fluctuations to achieve a critical level, and the tricritical temperature grows. At the same time, more intensive interactions between ions lead to more active clustering. The free charges density, and hence the electrostatic part of the free energy diminishes significantly. However, the entropic part of the energy remains the same. Therefore, electrostatic interactions of two charged sublattices become comparable with the entropy of charge separation at higher overall densities, and the coexistence shifts to denser phases.

The trend for tricritical temperature and density to increase with anisotropy is correctly reproduced by the Debye-Hückel treatment of anisotropic lattices kobelev02-2 : see Fig.2. However, the quantitative agreement between theoretical predictions and simulations is not very good. While the theory predicts only a few percent increase in density at large lattice stretching, our simulations showed more significant growth. For the tricritical temperature, on the contrary, simulations yield a somewhat smaller slope.

At large stretching, increasing interactions inside $`xy`$-planes result in lesser relative importance of interactions between ions with different $`z`$ coordinates. At some point, interactions with ions in the same layer start dominating, and the $`3D`$ system becomes quasi-$`2D`$. Since different layers become uncoupled, due to its probabilistic nature phase separation in each of them takes place at different values of the chemical potential, and one may have configurations with gas phase in some of the layers and high-density ordered phase in the others, which prevents reliable sampling of the system. At $`\alpha `$=10, it turned to be impossible to obtain a reliable phase separation in a periodic system of more than one layer in $`z`$-direction: at high densities the system would not equilibrate after more than $`10^9`$ Monte Carlo steps.

To this end, for $`\alpha =10`$ we carried out simulations of the system consisting of only one plane, that is, square box with $`L_x^{}=L_y^{}=32`$ periodic in $`x`$ and $`y`$ directions, with the potential (1). Using only one layer in simulations yielded a highly asymmetric phase coexistence very similar to that obtained by Teitel for Coulomb gas on two-dimensional square lattice teitel97 , with the tricritical point located at $`T_{tri}^{}=0.28,\rho _{tri}^{}=0.70`$: see Fig.1. But in two dimensions, the coefficient in the Laplace equation (6) is not $`4\pi `$ but $`2\pi `$. Since all our simulations were done with three-dimensional potential, the temperature at large anisotropies must be divided by 2 to obtain the correct $`2D`$ limit, while the density should remain intact. This yields $`T_{2D}^{}=0.14,\rho _{2D}^{}=0.70`$ which exactly coincides with the values obtained by Teitel.

It should be also noted that parts of the coexistence curve for $`\alpha =10`$ at intermediate densities lay below the corresponding curve for $`\alpha =4`$. However, for the latter case in contrast to $`\alpha =10`$, the interactions between layers are still very important and add significantly to the total free energy of the system: simulations of periodic system of only one layer produce much lower coexistence curve. Therefore, although at equal densities the energy of ions interactions inside a single layer for $`\alpha =4`$ is less than for $`\alpha =10`$, the total energy density may be higher because of the interactions between layers, and hence the temperature of phase separations would also be higher. However, this is not essential at very low or high densities where the interactions inside layers dominate anyway. Taking into account the influence of layers on each other at $`\alpha =10`$ probably would eliminate this intersection of the coexistence curves. Nevertheless, this hardly would change our estimate for the tricritical point location in the limit of infinite stretching, which, as it has been mentioned, coincides with the results by Teitel for two-dimensional square lattice.

To gain further insight on the behavior of LRPM on anisotropic lattices we also performed simulations with the continuous space Coulomb potential. The corresponding phase diagrams for box size $`L^{}=16`$ and $`\alpha 2`$ are presented in Fig.3. At any fixed lattice stretching, the phase coexistence is very similar to the case of the true lattice potential. No normal gas-liquid phase transition is found, and the phase coexistence is between a disordered low-density phase and ordered high-density phase. The tricritical temperature shows the same trend as for the correct lattice potential, increasing with the inter-layer distance. However, since for small $`\alpha `$ the interactions are weaker now, the tricritical temperatures are lower.

The scaling of the temperature $`T`$ with $`E_0=q^2/Db`$ is maintained in this paper to allow comparisons with earlier theoretical results. For the correct lattice (or continuum) two-dimensional potential, this is allowed as $`b\mathrm{}`$ because the potential becomes length scale-independent at the two-dimensional limit. However, for the $`1/r`$ potential, using a scaling with the layer-layer energy results in a divergence of the reduced critical temperature. While for the correct lattice potential the ”natural” scale of energy is set by layer-layer interactions whose asymptotic behavior does not change with $`\alpha `$ (see Table 1), for the $`1/r`$ potential this scale is given by the energy of the ion-ion interactions at contact inside a layer $`E_0^{\mathrm{}}=q^2/Da`$. Note, however, that for actual simulations this is not important, because all what matters is the Boltzmann factor, which remains the same as long as the energy and the temperature are measured in the same units.

Calculations with the temperature scaled by this energy $`T^{\mathrm{}}=kTDa/q^2`$ (where $`D`$ is the three-dimensional dielectric constant, appropriate for the $`1/r`$ potential) are shown in Fig. 4. With this ”natural” scaling, the tricritical point behavior as a function of anisotropy is drastically different: the tricritical temperature decreases with anisotropy, converging to $`T^{\mathrm{}}=0.11`$. The reason for this is that now increasing anisotropy is not accompanied by stronger interactions between ions. On the contrary, while interactions inside the $`xy`$planes remain the same, due to a weaker coupling of the layers effective energy per particle diminishes, and so does the tricritical temperature. This decay is rather fast and the temperature saturates already at the lattice parameters ratio $`\alpha 2`$. Detailed studies of ions configurations showed that at $`\alpha 2`$ different $`xy`$planes become totally independent of each other, and the system behaves as a quasi-$`2D`$ liquid. Similarly to the true lattice potential case, to obtain phase coexistence for the limit of infinite stretching ($`\alpha 2`$), simulations of periodic system of only one plane were carried out with box size $`L=32`$ and $`L=48`$.

The tricritical density gradually increases from $`\rho _{tri}^{}0.48`$ for isotropic lattice to $`\rho _{tri}^{}0.55`$ at strong anisotropy, following the trend shown by the system with correct lattice interactions: see Fig.3. This again reflects the relative importance of clustering when one approaches the two-dimensional limit of infinite stretching. However, since the potential is now weaker, the limiting density is lower than in true lattice $`2D`$ systems.

## IV Conclusions

In this study, grand canonical Monte Carlo simulations and histogram reweighting techniques have been used to investigate phase behavior of the restricted primitive model with correct lattice Coulomb potential on simple cubic and tetragonal lattices. Our results show that, at all degrees of anisotropy, only order-disorder phase separations and a tricritical point exist. For isotropic cubic lattice, the effect of using the correct lattice potential shows only in the value of the tricritical temperature. As anisotropy increases, both the tricritical temperature and density increase. At large anisotropy the system undergoes a qualitative shift from a three-dimensional to two-dimensional behavior. This shift is not accompanied by any change in the phase diagram topology, and only the tricritical parameters alter. To study in detail the transition from $`3D`$ to $`2D`$ behavior, further investigations are needed. In particular, its location may depend on the system size and boundary conditions used in Ewald summation. For continous $`1/r`$ potential, the tricritical temperature displays an opposite trend and decreases with the anisotropy. Although exact values of the tricritical temperature and density can be influenced by finite size effects, an issue totally omitted in this study, their qualitative behavior seems physically plausible. It would be also useful to develop another procedure for Ewald-like summation for lattice potential, which will not not rely on its asymptotic behavior.

The only existing theoretical treatment of Coulomb gas on anisotropic lattices based on the Debye-Hückel approach captures correctly the trend of the tricritical temperature and density behavior with anisotropy kobelev02-2 . However, the absence of the anticipated normal gas-liquid transition suggests that with increasing anisotropy, stability of both ordered and disordered dense phases increases in a similar way, and the formation of an ordered structure remains thermodynamically more favorable than possible gas-liquid phase coexistence. Taking explicitly into account the ion clustering will provide a better theoretical description of thermodynamics of ionic systems on anisotropic lattices.

## Acknowledgements

Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund (Grants No. 37867-G6 to ABK and 38165-AC9 to AZP) for support of this research. ABK also acknowledges the support of the Camille and Henry Dreyfus New Faculty Awards Program (Grant No. NF-00-056). The authors would like to thank Professor M.E. Fisher for critical comments and useful suggestions.

## Appendix A  Lattice Coulomb Potential: Exact Representation and Fast Calculation

Due to the discrete symmetry, the correct lattice Coulomb potential differs from continous $`1/r`$. The analytical expression for this potential follows from the lattice version of the Poisson equation kobelev02-1 ,

$$\mathrm{\Delta }\phi (𝐫)=\frac{4\pi }{Dv_0}\delta (𝐫),$$
(6)
where the exact form of the lattice Laplacian depends on the geometry of the lattice, $`D`$ is the dielectric constant of the media, and $`v_0`$ is the unit cell volume. For simple tetragonal lattice with lattice spacings $`a`$ in $`x`$ and $`y`$ directions and $`b`$ in $`z`$ direction one has

$$\mathrm{\Delta }\phi =\mathrm{\Delta }_x\phi +\mathrm{\Delta }_y\phi +\mathrm{\Delta }_z\phi ,$$
(7)
with

$$\mathrm{\Delta }_i\phi (𝐫)=1/a_i^2\left[\phi (𝐫a_i𝐞_i)2\phi (𝐫)+\phi (𝐫+a_i𝐞_i)\right],$$
(8)
where $`i=x,y,z`$; $`a_x=a_y=a,a_z=b`$ and $`𝐞_i`$ are the unit vectors along the corresponding lattice directions. After introducing the anisotropy parameter $`\alpha =b/a`$, the solution of the Poisson equation has the following form

$$\phi (𝐫)=\frac{1}{4\pi ^2Db}_\pi ^\pi \frac{e^{i\mathrm{𝐤𝐫}}d^3𝐤}{2+1/\alpha ^2(\mathrm{cos}k_1+\mathrm{cos}k_2+1/\alpha ^2\mathrm{cos}k_3)}.$$
(9)
However, the potential in the form of a triple integral of a periodic function is not convenient for numerical calculations. More practical expression for this integral has been obtained by Maradudin et al. maradudin

$$\phi (𝐫(l,m,n))=\frac{2\pi }{Db}_0^{\mathrm{}}\mathrm{exp}[(2+1/\alpha ^2)t]I_l(t)I_m(t)I_n(t/\alpha ^2)𝑑t,$$
(10)
where $`I_l`$ denotes a modified Bessel function of the first kind.

To enhance calculations of the lattice potential for the whole lattice, an error-free propagating algorithm, similar to that proposed by Friedberg and Martin friedberg  for isotropic lattice, has been used. First, we compute lattice potentials in planes $`x0z`$ and $`x1z`$ using Eq.(10). Now, for a given lattice cite ($`k,l,m`$), the anisotropic lattice potential satisfies the gradient equations

$`{\displaystyle \frac{\phi (l,m,n)\phi (l1,m,n)}{l}}`$ $`=`$ $`{\displaystyle \frac{\phi (l,m+1,n)\phi (l,m1,n)}{m}}`$ (11)
$`=`$ $`\alpha ^2{\displaystyle \frac{\phi (l,m,n+1)\phi (l,m,n1)}{n}}`$
which follow from Eq.(9). At fixed $`n`$, potential at any point $`(l,m,n)`$ can be obtained recursively from two lines, $`m=0`$ and $`m=1`$. For $`l=0`$ it would be indeterminate, but this is not important since due to the symmetry $`\phi (0,m,n)=\phi (m,0,n)`$. As Friedberg and Martin showed friedberg , this recursion does not lead to increase of the initial possible rounding error at $`m=0`$ and $`m=1`$ lines as long as $`|l||m|`$. For $`|l|>|m|`$ the symmetry $`\phi (l,m,n)=\phi (m,l,n)`$ can be employed.

Note that similar considerations can be done for initial $`x0z`$ and $`x1z`$ planes in attempt to generate the potentials from two lines, using the second equation in (11). In these planes the potential is not symmetric, and therefore two sectors, $`|l|\alpha ^2|n|`$ and $`|l|\alpha ^2|n|`$ need to be calculated separately, starting from the corresponding axes and then ”stitched” together at the $`|l|=\alpha ^2|n|`$ line. However, since all potentials and interactions are computed only once at the beginning of simulations and then fast look-up algorithm is used, this gives us a relatively small advantage in the overall computational time.


format: text

text:
# Detection of Low-Hard State Spectral and Timing Signatures from the Black Hole X-Ray Transient XTE J1650–500 at Low X-Ray Luminosities

## 1. Introduction

To fully understand the behavior of black hole candidate (BHC) X-ray transients, it is important to determine the source properties for the full range of observed X-ray luminosities from outburst peak at $`L10^{3839}`$ erg s-1 to quiescence at $`L10^{3033}`$ erg s-1. The brighter part of this range has been well-studied by several satellites, and, due to its broad bandpass (2-200 keV) and its excellent timing capabilities, the Rossi X-ray Timing Explorer (RXTE, Bradt, Rothschild & Swank 1993) has proved to be an extremely useful tool for studying the luminosity range above $``$$`10^{3536}`$ erg s-1 (for sources with typical distances of 3-10 kpc). For most systems, this is below the luminosity where the source makes a transition from the canonical high-soft state to the low-hard (or hard) state (McClintock & Remillard, 2003), and our group as well as others have used RXTE to study BHC transients as they undergo soft-to-hard state transitions during outburst decay (Tomsick & Kaaret, 2000; Nowak, Wilms & Dove, 2002; Kalemci, 2002; Kalemci et al., 2003a). However, as RXTE is not an imaging instrument, background and source confusion become significant at lower luminosities, making both spectral and timing studies problematic. Thus, at low luminosities, imaging observatories such as BeppoSAX, the Chandra X-ray Observatory (Weisskopf et al., 2002), and XMM-Newton are necessary, and these satellites have been used to study BHC transients in quiescence (Garcia et al., 2001). The aim of this study is to use Chandra to study the evolution of XTE J1650–500 in the luminosity range between the transition to the hard state and quiescence.

Few previous studies have been dedicated to observing BHC transients at luminosities between the hard state transition and quiescence even though the source evolution has not been well-characterized, and the physical changes that occur in this regime are unclear. The results of the studies of BHC transients that have been undertaken in this regime are intriguing. A correlation between X-ray and radio flux has been found for several BHC systems in the hard state (Corbel et al., 2003; Gallo, Fender & Pooley, 2003), and this has been interpreted as evidence for a connection between X-ray production and a compact radio jet (Markoff et al., 2003). As the sources decay, it has been observed that the X-ray energy spectra maintain a power-law shape; however, in some cases, there is evidence for spectral evolution with sources becoming softer at low luminosities (Ebisawa et al., 1994; Tomsick, Corbel & Kaaret, 2001; Corbel et al., 2003; Kalemci, 2002). X-ray timing studies immediately after the hard state transition show that, in most cases, the characteristic frequencies drop (Tomsick & Kaaret, 2000; Nowak, Wilms & Dove, 2002; Kalemci et al., 2003a), and this may indicate movement of the inner edge of the optically thick accretion disk away from the black hole. However, for the most part, the X-ray timing studies have not extended to luminosities much below the hard state transition.

The hard state accretion geometry, the emission mechanisms, and the connection between the accreting material and the compact jet outflow that is present in the hard state for at least some BHCs (Fender, 2001) are all currently active areas of debate. Several theoretical models have been suggested. One possibility is that a quasi-spherical, optically thin region forms in the inner portion of the accretion flow as in advection-dominated accretion flow (ADAF) models (Narayan, McClintock & Yi, 1996). The X-ray emission is predominantly due to thermal Comptonization in the hard state for ADAF models and for more generic “sphere+disk” models, and some success has been achieved in comparing such models to observations (Esin, McClintock & Narayan, 1997; Nowak et al., 1999). However, for a range of disk viscosities, the ADAF model is convectively unstable (Quataert & Gruzinov, 2000), and this should be considered before concluding that any ADAF or sphere+disk model provides the correct physical description. Another possibility that implies a much different site for X-ray production and a different accretion geometry is that the X-ray emission is due to magnetic flares above the disk (Galeev, Rosner & Vaiana, 1979; di Matteo, Celotti & Fabian, 1999). Although the presence of a jet is not incorporated into the basic models described above, some recent work describes how a jet may arise from a “magnetic corona” (Merloni & Fabian, 2002). Other models have focused on the emission from the compact jet, and one suggestion is that the X-ray emission is produced in the jet via a synchrotron mechanism (Markoff, Falcke & Fender, 2001). Improving the luminosity coverage for the hard state is important for providing tests of these models. In addition, it is currently unclear whether the same physical model can apply to these systems in the hard state and in quiescence. Major physical changes to the system could be seen as another state transition at very low luminosities. On the other hand, the lack of a sharp transition could indicate that quiescence is simply a low luminosity extension of the hard state as has been previously suggested for GX 339–4 (Corbel et al., 2000; Kong et al., 2000).

We chose to carry out our program by observing XTE J1650–500. This was a new BHC transient when it was discovered in outburst in 2001 (Remillard, 2001). Although the compact object mass is not well-constrained so that it has not been confirmed that this is a black hole system, its X-ray properties strongly suggest that it harbors a black hole (Tomsick et al., 2003; Homan et al., 2003; Kalemci et al., 2003b). During its outburst, the source was optically identified (Castro-Tirado et al., 2001), and a radio counterpart was found (Groot et al., 2001). In quiescence, a radial velocity study was carried out, indicating a binary orbital period of 5.1 hr and an optical mass function of $`0.64\pm 0.03`$ $`\mathrm{M}_{\mathrm{}}`$ (Sanchez-Fernandez et al., 2002). The value for the mass function is relatively low for a black hole system and could indicate that the system contains a relatively low mass black hole, has a low binary inclination ($`i<40^{}`$), or both. We previously reported on RXTE observations of XTE J1650–500, including a study of its timing properties near the soft-to-hard state transition (Kalemci et al., 2003b). We also reported on 14 day X-ray oscillations and $``$100 s X-ray flares that were detected with RXTE near the end of the outburst (Tomsick et al., 2003). The current work focuses on our studies of the source with Chandra.

## 2. Observations and Analysis

We obtained 3 Chandra observations of XTE J1650–500 near the end of its 2001-2002 outburst, and the times of the observations are marked on the long-term X-ray light curve shown in Figure 1. We compiled the 3-20 keV light curve in Figure 1 using all the observations of XTE J1650–500 that were made with the Proportional Counter Array (PCA) on RXTE. The 3 Chandra observations were made 67, 78, and 104 days after the soft-to-hard state transition (Homan et al., 2003; Kalemci et al., 2003b), during the time when the source was exhibiting large-amplitude flux oscillations with a characteristic time scale of 14 days (Tomsick et al., 2003). All 3 observations were made at flux levels well below the level of the soft-to-hard state transition.

The Chandra observations were made using the Advanced CCD Imaging Spectrometer (ACIS). In each case, the source was placed on one of the back-illuminated ACIS chips (S3). The first two observations (1 and 2), were made with the High Energy Transmission Grating (HETG) inserted, and each observation lasted approximately 10 ks (see Table 1 for exact times and other detailed information about the observations). To increase the throughput, the grating was not used for observation 3, making photon pile-up (Davis, 2001) a concern. To mitigate the effects of pile-up, we reduced the CCD exposure time from 0.4 s per frame from the nominal 3.2 s per frame by using a 1/8 CCD sub-array. In addition, the telescope was pointed $`2^{}.7`$ away from the target, slightly blurring the point spread function (PSF) for XTE J1650–500. An exposure time of 18 ks was obtained for observation 3. In addition to the Chandra observations, shorter 1.8-2.4 ks RXTE observations were co-ordinated to cover part of each Chandra observation.

### 2.1. Chandra Analysis

For the two grating observations, we used the “level 2” event list produced by the standard data processing (ASCDS v6.5.1 using CALDB v2.11). For both observations, we applied a standard correction to the ACIS-S4 chip to correct a flaw in the chip readout. For this, as well as the data processing described below, we used the CIAO (Chandra Interactive Analysis of Observations) v2.3 software and CALDB (Calibration Data Base) v2.21. Prior to extracting the grating spectra, we produced background light curves. For observation 1 only, we found a few background flares where the count rate increased by 50-100% for a time period of less than 50 s. Due to their short duration, the flares make only a small contribution to the total background. Thus, we did not remove these times in the subsequent analysis, and the full integration time was used for both observations.

For observations 1 and 2, we produced 0.3-8 keV zero order images. We used the source detection algorithm wavdetect to search for sources in a 2 by 2 arcminute region around the Chandra aim-point, and we only detect the target in each case. The observation 1 and 2 source positions are consistent with one another, and we obtain a mean source position of R.A. = 16h 50m 00s.98, Dec. = –49 57 43′′.6 (equinox 2000.0) with an uncertainty of $`0^{\prime \prime }.6`$. This position is consistent with the previously reported radio (Groot et al., 2001) and optical (Castro-Tirado et al., 2001) positions, and we note that the Chandra position is somewhat more accurate than the $`3^{\prime \prime }`$ and $`2^{\prime \prime }`$ radio and optical positions, respectively. We conclude that the detected source is indeed XTE J1650–500. Given the recent detections of X-ray jets from XTE J1550–564 (Corbel et al., 2002) and 4U 1755–338 (Angelini & White, 2003), we examined the XTE J1650–500 images for X-ray extensions in the vicinity of the source, but we do not see any features that are likely to be jets. Using the zero order ACIS response matrix and a spectral shape similar to that observed for the XTE J1550–564 jet ($`\mathrm{\Gamma }=1.6`$ with interstellar absorption), the upper limit on the 0.3-8 keV absorbed flux is $`6\times 10^{14}`$ erg cm-2 s-1 at positions $`10^{\prime \prime }`$ from XTE J1650–500.

We used the CIAO routine tgextract for the extraction of grating spectra and the CIAO routines mkgrmf and mkgarf to produce response matrices. Combining the +1 and $`1`$ MEG (Medium Energy Grating) 0.7-7 keV spectra yields 4786 counts in 10001 seconds for observation 1 and 4342 counts in 9509 seconds for observation 2. For the +1 and $`1`$ HEG (High Energy Grating) 1-9 keV spectra, there are 2706 counts and 1872 counts for observations 1 and 2, respectively. We also produced estimated background spectra; however, we find that the background is not significant with estimates of $``$20 counts for the MEG and $``$11 counts for the HEG (0.4% of the detected counts in both cases). The higher MEG and HEG orders contain between $``$30 and $``$200 counts (source plus background). For each observation, the zeroth order contains $``$2500 counts, but there is a high level of photon pile-up. The flux obtained from the grating spectra implies that the number of counts would be a factor of 3 higher without photon pile-up. Thus, due to the low number of counts in grating orders 2 and 3 and the high level of pile-up in the zeroth order, we focus on the order 1 spectra in this work. A final step in the extraction of the spectra is to apply a correction to the response matrix to account for the gradually changing low energy ACIS response. We used the script corrarf to apply this correction.

For observation 3, we used the “level 2” event list produced by the standard data processing (ASCDS v6.6.0 and CALDB v2.12). After inspecting an image including the full energy band, we produced a 0.3-8 keV image to reduce the background level. With the 1/8 sub-array and only the S3 chip active, the field of view is approximately $`1^{}`$ by $`8^{}.4`$ (128 by 1024 pixels), and we used wavdetect to search for sources. We found 4 sources detected at greater than 4-$`\sigma `$ significance. The brightest source by far is at a position consistent with the XTE J1650–500 position found above for the grating observations. Although XTE J1650–500 appears to be elongated, we used CIAO to simulate an off-axis PSF, and we found that XTE J1650–500 is consistent with a point-like source. The other 3 sources also appear to be point-like, and the closest source is more than $`1^{}`$ from XTE J1650–500. Despite using an instrument setup designed to mitigate pile-up, the measured ACIS count rate of 0.493 s-1 for XTE J1650–500 (see Table 1) indicates that the data at the PSF core is moderately affected by photon pile-up. However, we were still able to use this data for spectral and timing analysis as described below and in Appendix A.

### 2.2. RXTE Analysis

We initially performed spectral and timing analysis for the 3 RXTE observations that were made during the Chandra observations. Using the “Standard 2” PCA data, we extracted energy spectra using scripts developed at UC San Diego and the University of Tübingen that incorporate the standard software for RXTE data reduction (FTOOLS v5.2). We used the 2002 February release of the “Faint Source” model for background subtraction. Obtaining the best possible background subtraction is important for this work because of the low source count rates. To achieve this, we used Standard 2 data from the top anode layer only and excluded data from PCU 0, which has a higher background level due to the loss of its propane layer111See http://lheawww.gsfc.nasa.gov/users/craigm/pca-bkg/bkg-users.html for a detailed analysis of the background model performance..

In addition to the 3 RXTE observations, we extracted energy spectra for the final 8 observations for which the PCA count rates are shown in Figure 1. For these observations, the mean 5-8 keV PCA count rate is $`0.021\pm 0.018`$ s-1 per Proportional Counter Unit (PCU) using all three anode layers, indicating that the rate is consistent with zero at slightly more than the 1-$`\sigma `$ level. We determined this rate in order to make a comparison to the rate expected from the Galactic ridge emission at the Galactic latitude of XTE J1650–500 ($`b=3.44^{}`$). Valinia & Marshall (1998) have determined the level of Galactic ridge emission (in PCA count rate) averaged over Galactic longitude as a function of $`b`$ (see Figure 2 of Valinia & Marshall 1998). They found that the distribution could be described by narrow and broad gaussian components, and only the broad gaussian contributes at the Galactic latitude of XTE J1650–500. Using the parameters for the broad gaussian leads to a best estimate of the 5-8 keV Galactic ridge rate of 0.71 s-1 for all anode layers and 5 PCUs. Scaling the observed count rate given above (0.021 s-1) to 5 PCUs gives $`0.11\pm 0.09`$ s-1, which is significantly below the best estimate for the Galactic ridge emission. Thus, we conclude that the Galactic ridge emission can easily explain the low level of emission we detect for the final 8 observations. As this analysis shows that it is likely that XTE J1650–500 does not contribute to the flux detected during these observations, we produced an average spectrum (top anode layer only) for the final 8 observations and subtracted it from the PCA spectra for the 3 RXTE observations that were simultaneous with the Chandra observations. The effect to the continuum is small; however, the Galactic ridge emission contains a strong iron emission line at 6.7 keV, and this line is apparent in the average spectrum for the final 8 observations.

Finally, we analyzed the data from the 15 RXTE observations taken between MJD 52265 and MJD 52298, which is after the gap in RXTE coverage and before the Chandra observations began (see Figure 1). We processed the PCA data as described above, including the subtraction of the Galactic ridge emission from the energy spectrum. For the spectral analysis described in § 3.3, we also removed flares from 2 of the observations as described in Tomsick et al. (2003).

## 3. Results

### 3.1. Energy Spectra at Low Flux Levels

Using XSPEC 11.2, we performed simultaneous least-squares fits to the Chandra and PCA spectra for all three observations. For the PCA, we included 0.6% and 0.3% systematic errors for energy bins below and above 8 keV, respectively, as described in our previous work (Tomsick et al., 2003). We began by fitting the observation 1 spectrum with a power-law with interstellar absorption. We also included a constant factor to allow for the possibility that the overall normalization is different for the PCA and Chandra, but we assumed that the overall MEG and HEG normalizations agree with each other. As shown in Table 2, a power-law model with $`\mathrm{\Gamma }=1.66\pm 0.05`$ and $`N_\mathrm{H}=(6.7\pm 0.5)\times 10^{21}`$ cm-2 (errors are 90% confidence, $`\mathrm{\Delta }\chi ^2=2.7`$) provides a good fit to the spectrum ($`\chi ^2/\nu =116/155`$). The other continuum models we tried do not provide acceptable fits. Thermal models, such as a blackbody or a disk-blackbody (Makishima et al., 1986), result in reduced $`\chi ^2`$ values of 8.4 and 2.3, respectively. Figure 2 shows the observation 1 spectrum and residuals with the power-law fit. No clear features are seen in the residuals, and we note that the fit is not improved by adding a soft component, a high energy cutoff, or an iron emission line. The upper limit on the equivalent width of a narrow (where “narrow” is set by the binning shown in Figure 2) Fe K$`\alpha `$ emission line at 6.4 keV is 161 eV.

Table 1 and Figure 1 indicate that the count rates are very similar for observations 1 and 2. We fitted the observation 2 spectrum, and Table 2 shows that the spectral parameters do not change significantly between the two observations. Also, in both cases, the spectral fits indicate that the PCA/Chandra normalization factor is near 1.12. As the PCA exposure only covers part of the Chandra observation, we checked to see if this difference could be due to source variability. We extracted a grating spectrum for observation 1 that is strictly simultaneous with the PCA exposure. During this time, the MEG count rate is $`0.483\pm 0.016`$ s-1 (0.7-7 keV), which is consistent with the rate for the full observation ($`0.479\pm 0.007`$ s-1). Thus, we conclude that the overall normalizations are different for the PCA and the Chandra gratings, and we note that this is not surprising since it has been previously determined that the PCA normalization is high relative to previous instruments (Tomsick et al., 1999).

Due to the presence of photon pile-up in the PSF core, we extracted the observation 3 ACIS spectrum after removing the counts detected in the PSF core. We extracted the source counts from an elliptical annulus centered on the XTE J1650–500 position. We used an inner ellipse with a semi-major axis of $`1^{\prime \prime }.3`$ and a semi-minor axis of $`0^{\prime \prime }.9`$, and an outer ellipse with a semi-major axis of $`4^{\prime \prime }.9`$ and a semi-minor axis of $`3^{\prime \prime }.1`$. We rotated the ellipses to match the off-axis Chandra PSF. This region contains 1573 counts compared to the 9000 counts obtained when the core is included (see Table 1), and we are confident that the photons from the elliptical annulus do not suffer from pile-up. However, it is critical to realize that ACIS response matrices produced by the standard CIAO software do not account for the size or shape of the source extraction region or the fact that the Chandra PSF is energy dependent. We corrected the ACIS response matrix for this using the Chandra Ray Tracing software (ChaRT) as described in Appendix B. When producing the response matrix we also applied the low energy ACIS response correction for detector evolution as described above for the grating spectra.

We fitted the observation 3 ACIS+PCA spectrum with an absorbed power-law model as shown in Figure 3. This model provides a good fit to the data ($`\chi ^2/\nu =104/100`$) with parameter values of $`N_\mathrm{H}=(6.4\pm 0.7)\times 10^{21}`$ cm-2 and $`\mathrm{\Gamma }=1.93\pm 0.13`$ (see Table 2). The measured column density is consistent with the value obtained from the grating observations. Thus, we re-fitted the spectrum after fixing the column density to the value found using the grating spectra ($`N_\mathrm{H}=6.7\times 10^{21}`$ cm-2), and we obtain $`\mathrm{\Gamma }=1.96\pm 0.09`$. Also, for consistency with observations 1 and 2, we used the PCA normalization for observation 3 divided by the PCA/grating normalization factor of 1.12 to determine the observation 3 flux reported in Table 2. The spectral results indicate that the XTE J1650–500 keeps its power-law shape but softens somewhat at low luminosity.

### 3.2. Power Spectra at Low Flux Levels

We examined the PCA and Chandra light curves for all three observations. As shown in Figure 4, clear variability is present in the PCA light curves for observations 1 and 2 and in the ACIS light curve for observation 3. The plotted 3-20 keV PCA light curves (panels a and b) have 4 s time bins, and we note that the count rate is higher than reported in Table 1 because we used all anode layers and PCUs 0, 2, and 3 rather than just the top layer and PCUs 2 and 3. While the quality of the background estimate is slightly worse when PCU 0 and all layers are included, the higher count rate allows for tighter constraints on the timing properties. We also produced grating light curves (MEG+HEG) for observations 1 and 2. Although the statistical quality of the PCA light curves is superior to that of the grating light curves, we use the grating light curves for the analysis below as they allow us to study the energy dependence of the timing properties. Figure 4c shows the 0.5-8 keV ACIS light curve for observation 3 with $``$50 s time bins. We also examined the PCA light curve for observation 3, but it is background dominated and is not useful for timing analysis.

For observations 1 and 2, we produced power spectra using PCA data taken in a high time resolution “event” mode (labeled E\_125us\_64M\_0\_1s). These data consist of an event list with photons time tagged to within 125 $`\mu `$s and 64 energy bins. First, we used the data in the 3-20 keV band to produce Leahy normalized power spectra (Leahy et al., 1983) for observations 1 and 2 separately. In each case, we averaged power spectra from four 460 s light curve segments and used a Nyquist frequency of 64 Hz. No noise in excess of the Poisson level is present above $``$1 Hz, and there were no clearly significant differences between the power spectra for the two observations. Thus, we produced another Leahy normalized power spectrum (again with a 3-20 keV energy band) consisting of the eight 460 s light curve segments from observations 1 and 2 with a Nyquist frequency of 8 Hz. We converted to the rms normalization (Miyamoto et al., 1994), and this power spectrum is shown in Figure 5a. The power spectrum is dominated by a band-limited noise component, and such a component has typically been fitted with a broken power-law or a zero-centered Lorentzian in previous work (Belloni, Psaltis & van der Klis, 2002). A zero-centered Lorentzian with a half-width (FWHM/2) of $`0.067\pm 0.007`$ Hz, and a fractional rms amplitude of 33.1%$`\pm 1.1`$% (see Table 3) provides a reasonably good fit to the power spectrum ($`\chi ^2/\nu =93/72`$). We note that the half-width of a zero-centered Lorentzian is equivalent to the “peak” frequency that is sometimes quoted when fitting Lorentzians to power spectra (Belloni, Psaltis & van der Klis, 2002; Kalemci et al., 2003b). Positive residuals are present near 0.1 Hz, possibly indicating the presence of a quasi-periodic oscillation (QPO). The addition of a second Lorentzian at 0.11 Hz improves the quality of the fit, resulting in $`\chi ^2/\nu =83/69`$; however, this improvement corresponds to a detection at only the 96% confidence level. Although we repeated this analysis several times with different frequency binnings, we consistently found that the QPO is only detected at a level of between 2 and 3-$`\sigma `$. Thus, we take the measured rms amplitude of 12% as an upper limit on a possible QPO at 0.11 Hz.

Using the grating light curves, we produced a 0.7-9 keV power spectrum for observations 1 and 2 by averaging eight 2414 s light curve segments. The time resolution is set by the interval of time between chip read-outs, and this interval is 2.54104 s for these observations, giving a Nyquist frequency of 0.2 Hz. Although the time resolution and the statistical quality are lower for the grating light curves than for the PCA data, using the grating data allows us to extend the timing study to lower energies and frequencies. After converting to rms normalization, we fitted the grating power spectrum with a zero-centered Lorentzian, and the fit is acceptable ($`\chi ^2/\nu =25/23`$). As shown in Table 3, the half-width is $`0.079\pm 0.020`$ Hz, and the fractional rms is 30.1%$`\pm 2.5`$%. Thus, the parameters for the 0.7-9 keV power spectrum are consistent with those for the 3-20 keV PCA power spectrum.

We also produced a power spectrum using the observation 3 ACIS data. Using a 0.5-8 keV light curve with a time resolution of 0.44104 s, we produced a Leahy normalized power spectrum with a Nyquist frequency of 1.13 Hz, consisting of the average of six 3350 s light curve intervals. In the Leahy normalization, Poisson noise leads to a power level of 2.0; however, at high frequencies (0.1-1 Hz) the ACIS power spectrum is consistent with a constant level of 1.77. We performed simulations described in Appendix A showing that pile-up at the level present for observation 3 leads to a drop in the Poisson noise level from 2.0 to 1.78. This, along with a small drop in the Poisson noise level caused by deadtime during chip read-outs (see Appendix A) explains our measurement of 1.77. Thus, when converting the ACIS power spectrum from the Leahy normalization to the rms normalization, we subtracted a value of 1.77 rather than the usual value of 2.0. The rms normalized power spectrum fitted with a zero-centered Lorentzian is shown in Figure 5b, and this model provides a good fit to the data ($`\chi ^2/\nu =56/65`$). The Lorentzian half-width is $`0.0035\pm 0.0010`$ Hz (see Table 3), which is significantly lower than the value obtained for observations 1 and 2. This indicates that, as the flux drops, the band-limited noise moves to lower frequencies. This fact is illustrated in Figure 5a, where the model used to fit the observation 3 power spectrum is plotted as a dashed line on the power spectrum for observations 1 and 2. The fit to the observation 3 power spectrum indicates a fractional rms amplitude of 17.7%$`\pm 1.7`$%, but this (16%) should be taken as a lower limit on the actual value since the noise level is reduced by pile-up. From our simulations, we estimate that the actual rms amplitude is near 19%. We also produced a power spectrum using only the 1573 counts from the region of the CCD that is free of pile-up. Fitting the power spectrum with a zero-centered Lorentzian, the half-width is not well-constrained, but it is consistent with the previous fit, and the rms amplitude is 19%$`\pm 7`$%, giving an upper limit of 26%. Thus, considering the fits to the power spectra with and without pile-up, we conclude that the rms is 19%$`{}_{3\%}{}^{}{}_{}{}^{+7\%}`$.

### 3.3. Summary of XTE J1650–500 Hard State X-Ray Observations

For observations 1-3, two of our main results are that the source softens ($`\mathrm{\Gamma }`$ changes from $`1.66\pm 0.05`$ to $`1.93\pm 0.13`$) and the band-limited noise component moves to lower frequencies (the Lorentzian half-width, FWHM/2, changes from $`0.067\pm 0.007`$ Hz to $`0.0035\pm 0.0010`$ Hz) as the source flux decays. Thus, where possible, we determined the values of these parameters ($`\mathrm{\Gamma }`$ and FWHM/2) for other RXTE observations of XTE J1650–500 in the hard state during outburst decay. For the 7 observations made immediately after the soft-to-hard state transition on MJD 52231.5 (see Figure 1), we use the spectral parameters from Kalemci (2002). These energy spectra are dominated by a power-law component, but it should be noted that a soft component is also present. For these 7 observations, the timing parameters come from Kalemci et al. (2003b). In contrast to our Chandra/RXTE observations at low flux levels, several Lorentzians are required to describe the power spectra for most of the observations immediately after the soft-to-hard state transition. The lowest frequency component in those power spectra are nearly zero-centered Lorentzians, and we assume that they correspond to the band-limited noise component that we see in the low flux power spectra.

For the 15 RXTE observations after the gap in RXTE coverage and before the Chandra observations began, the 3-20 keV energy spectra are all well-described by a single power-law component with the column density fixed to the value of $`N_\mathrm{H}=6.7\times 10^{21}`$ cm-2 found for the grating observations. For these observations, the 3-20 keV flux varies from $`1.2\times 10^{11}`$ to $`2.9\times 10^{10}`$ erg cm-2 s-1, and the values of $`\mathrm{\Gamma }`$ range from $`1.69\pm 0.02`$ to $``$2.0. The data suggest a trend of softening as the flux drops consistent with that observed with Chandra, but the errors on $`\mathrm{\Gamma }`$ are large for the low flux observations. We produced power spectra for all 15 observations; however, we could not obtain reliable fit parameters for about half of the observations due to poor statistics, short observations, or flaring (Tomsick et al., 2003). In 8 cases, we obtained useful power spectra, and these are dominated by a band-limited component with a fractional rms amplitude near 30%. We fitted these power spectra with a zero-centered Lorentzian to determine the half-width for comparison to the other observations. Figure 6a shows the measured values of the Lorentzian half-width vs. 3-9 keV flux, and Figure 6b shows the half-width vs. the measured values of $`\mathrm{\Gamma }`$.

For the data points above a half-width of 0.5 Hz, there is a tight correlation between $`\mathrm{\Gamma }`$ and the half-width with a linear correlation coefficient of 0.96. The two points near 0.25 Hz also appear to be consistent with this correlation, but the correlation coefficient drops slightly to 0.91 when these points are included. It is clear that the correlation does not extend to lower frequencies as $`\mathrm{\Gamma }`$ reaches a minimum near 1.65-1.70. In fact, the data at frequencies between 0.03 Hz and 0.3 Hz suggest the possibility of an anti-correlation between $`\mathrm{\Gamma }`$ and the half-width, but the linear correlation coefficient is only –0.41. We checked to see if it is valid to compare the Chandra/RXTE point in this frequency range to the RXTE only points by fitting the RXTE spectra alone for observations 1 and 2. Fixing $`N_\mathrm{H}`$ to $`6.7\times 10^{21}`$ cm-2, we obtained values of $`\mathrm{\Gamma }`$ consistent with those found using Chandra and RXTE together. Although the anti-correlation is not statistically significant, the distribution of the values of $`\mathrm{\Gamma }`$ for the 0.03-0.3 Hz points is also not well-described by a constant at their weighted average of $`\mathrm{\Gamma }=1.70\pm 0.01`$. Thus, we conclude that there is some scatter in $`\mathrm{\Gamma }`$ for the 0.03-0.3 Hz points, but we cannot establish a definite relationship between $`\mathrm{\Gamma }`$ and the Lorentzian half-width in this frequency regime.

## 4. Discussion

One of the main motivations for this work is to measure the X-ray properties of an accreting black hole in the relatively unexplored luminosity range between the transition to the low-hard (or hard) state and the quiescent state. Such measurements allow us to address the question of whether these should be considered as two separate states or if the quiescent state is, in fact, a low luminosity version of the hard state. This question has important implications for whether there are significant physical changes in accreting black hole systems that occur at low luminosities such as changes in the geometry of the accretion flow and/or emission mechanisms.

For XTE J1650–500, we have measured the spectral and timing properties of the source at 1-9 keV (absorbed) luminosities of $`9\times 10^{34}`$ erg s-1 (observations 1 and 2) and $`1.5\times 10^{34}`$ erg s-1 (observation 3). The quoted luminosities assume a fiducial source distance of 4 kpc, and we note that the unabsorbed luminosities are about 25% higher. As discussed in Tomsick et al. (2003), the likely distance for XTE J1650–500 is between 2 kpc and 6 kpc. For the RXTE observation immediately after the soft-to-hard state transition, we used the spectral model described in Kalemci (2002) to extrapolate the PCA-measured flux into the 1-9 keV energy band, and we infer a transition luminosity of $`7\times 10^{36}`$ erg s-1, which is in-line with typical transition luminosities (Tanaka & Shibazaki, 1996). The observation 3 luminosity, which is a factor of 470 lower than the transition flux, is only about an order of magnitude above the quiescent luminosity of the black hole system that is brightest in quiescence (V404 Cyg, Garcia et al. 2001). Thus, our XTE J1650–500 observations cover a significant portion of the luminosity range between the soft-to-hard state transition and quiescence, but we note that the quiescent luminosity for XTE J1650–500 could be considerably lower than for V404 Cyg.

The trend observed for XTE J1650–500 that the band-limited noise gradually shifts to lower frequency immediately after the transition to the hard state (Kalemci et al., 2003b) has also been observed for several other black hole systems (Kalemci et al., 2003a; Nowak, Wilms & Dove, 2002). The drop in frequency can be explained by accretion models consisting of an inner optically thin corona and an outer optically thick disk. As the mass accretion rate drops, the inner radius of the optically thick disk ($`R_{in}`$) moves away from the black hole, causing a drop in the dynamical time scales and, thus, the characteristic frequencies for the system (di Matteo & Psaltis, 1999). After the XTE J1650–500 transition to the hard state, the Lorentzian half-width is correlated with the value of the power-law index down to at least a half-width of 1 Hz and possibly down to 0.25 Hz (see §3.3). Correlations between frequency (either half-width or QPO frequency) and $`\mathrm{\Gamma }`$ have been seen before during RXTE observations of Cyg X-1, GX 339–4, and other black hole systems in their hard states (Gilfanov, Churazov & Revnivtsev, 1999; Revnivtsev, Gilfanov & Churazov, 2001; Kalemci, 2002). As these authors argue, this correlation is expected for the disk and inner corona picture because the soft photons from the optically thick disk cool the corona, leading to a softer power-law index. Thus, as $`R_{in}`$ increases, and the characteristic frequencies drop, the corona is subject to a lower flux of soft photons from the disk, and the coronal spectrum hardens. This effect is predicted in both ADAF (Esin, McClintock & Narayan, 1997) models and in general disk and inner corona models such as the sphere+disk model (Nowak, Wilms & Dove, 2002) where the X-ray emission from the corona is predominantly due to thermal Comptonization.

While the half-width and $`\mathrm{\Gamma }`$ are correlated above 0.25-1 Hz, $`\mathrm{\Gamma }`$ does not get any harder than 1.65-1.70 for XTE J1650–500, and the spectrum softens to $`1.93\pm 0.13`$ as the half-width and flux continue to drop. This behavior can be interpreted within the ADAF or sphere+disk (henceforth ADAF/SD) models. The decrease in the half-width frequency indicates that the disk inner radius continues to recede; however, the disk flux is at such a low level that it no longer significantly cools the corona, and the spectrum does not continue to get harder. Instead, the hardness of the spectrum depends on the “Compton y-parameter,” which is a function of the optical depth and the temperature of the corona (Rybicki & Lightman, 1979). As the mass accretion rate drops, the coronal optical depth decreases, and a softer spectrum is predicted. Softer spectra have been observed at low luminosities for several black hole sources besides XTE J1650–500. XTE J1118+480 shows significant softening that starts at an X-ray luminosity of $`2\times 10^{35}`$ erg s-1 ($`d`$/1.8 kpc)2, and, at the end of its 2001 outburst, XTE J1550–564 begins gradual softening at a luminosity of $`2\times 10^{36}`$ erg s-1 ($`d`$/5.3 kpc)2 (Kalemci, 2002). Other examples of spectral softening at low luminosities include GS 1124–68 (Ebisawa et al., 1994), as well as some black hole sources in or close to quiescence (McClintock et al., 2003; Corbel et al., 2003; Kong et al., 2002; Tomsick, Corbel & Kaaret, 2001).

Although spectral fits with the ADAF model imply relatively large changes in $`R_{in}`$ between the soft and hard states, some previous black hole timing studies conclude that such changes, which can be as large as factors of 30-3000 (Esin, McClintock & Narayan, 1997), do not occur. For several black hole systems, di Matteo & Psaltis (1999) find that the characteristic frequencies do not vary by more than a factor of 50 between states. For Keplerian orbits, $`R_{in}\nu ^{2/3}`$, where $`\nu `$ is the characteristic frequency, so that a factor of 50 in frequency implies a change in $`R_{in}`$ by a factor of 14, which is significantly less than indicated by ADAF fits. However, for XTE J1650–500, we find that the Lorentzian half-width changes by a factor of nearly 1200 while the source remains in the hard state, implying that $`R_{in}`$ can change by as much as a factor of 113. Results of fast optical photometry for quiescent black hole systems suggest that $`R_{in}`$ may change by an even larger factor. For A 0620–00, Hynes et al. (2003) find that the optical power spectrum is dominated by band limited noise with a break frequency of $`9.5\times 10^4`$ Hz, which is a factor of 3.7 lower than the lowest frequency we find for XTE J1650–500. If the frequency is, in fact, related to $`R_{in}`$, then our result for XTE J1650–500 and the results of Hynes et al. (2003) imply that very large values of $`R_{in}`$ occur at low luminosities, and that the $`R_{in}`$ evolution occurs gradually and continues to change all the way down to quiescent luminosities. These results also imply that an extremely large value of $`R_{in}`$ is not required for a source to show hard state properties.

While the ADAF/SD models described above can explain many of the X-ray properties we see for XTE J1650–500, radio observations of XTE J1650–500 made after the transition to the hard state indicate the likely presence of a compact jet (Corbel et al., in prep.), which is not accounted for in these models. However, the “magnetic corona” model of Merloni & Fabian (2002) also explains several of the XTE J1650–500 X-ray properties while predicting that energy stored in the corona should be available to drive strong outflows. In this model, the X-ray emission mechanism is also Comptonization, so that a high flux of disk photons leads to a softer energy spectrum. Thus, Merloni & Fabian (2002) predict that the spectrum will harden as the mass accretion rate drops, and they also predict that the spectrum will soften at very low mass accretion rates as more of the coronal energy is carried away by the outflow. Recent work by Fender, Gallo & Jonker (2003) suggests that we are observing XTE J1650–500 in a luminosity regime where this could occur. While the XTE J1650–500 spectral evolution is consistent with the Merloni & Fabian (2002) model, and we also previously argued that this model may be able to account for the X-ray flares reported in Tomsick et al. (2003), it is unclear if the magnetic corona model can account for the fact that the observed characteristic frequency for XTE J1650–500 changes by a factor of 1200 in the hard state. Perhaps this could indicate large variations in the size of the magnetic corona or in the sizes of the active coronal regions, but since $`R_{in}`$ is not expected to change by a large factor in this model, we do not see an obvious explanation for the large frequency range.

A jet-dominated model has also been suggested for black holes in their hard state. In this model, the compact jet contributes a large fraction of the emission from the radio band up to X-ray energies, and a synchrotron emission mechanism is invoked (Markoff, Falcke & Fender, 2001). Evidence in favor of this model includes observations of correlated radio and X-ray flux for several black hole systems (Corbel et al., 2003; Gallo, Fender & Pooley, 2003), and good agreement between the predictions of the model and the broadband GX 339–4 energy spectra over a large range of mass accretion rates (Markoff et al., 2003). In the model, X-rays are produced when electrons are accelerated at a shock that is $`z_{sh}10^3R_g`$ from the black hole, where $`R_g=GM/c^2`$ and $`M`$ is the black hole mass. The X-ray spectrum is predicted to have a power-law shape out to a cutoff energy, and the power-law index is set by the index of the electron energy distribution, $`p`$, where $`dN/dEE^p`$. Markoff et al. (2003) find that they can explain the radio/X-ray flux correlation for GX 339–4 by varying only the jet power, keeping $`p`$ and $`z_{sh}`$ constant. The fact that $`\mathrm{\Gamma }`$ changes for XTE J1650–500 and also for GX 339–4 (Corbel et al., 2003), indicates that if this model is correct, $`p`$ likely does change with source luminosity, but we suspect that changes in $`\mathrm{\Gamma }`$ are probably not a major problem for this model as one might expect $`p`$ to vary somewhat with jet power. On the other hand, one would expect that the characteristic frequencies of the system would be related to the size of the source. Assuming that the size of the X-ray emitting region is related to the location of the shock, $`z_{sh}`$, it may be difficult for the jet-dominated model to explain the large range of characteristic frequencies we observe for XTE J1650–500 without varying $`z_{sh}`$.

## 5. Conclusions

We have compared the emission properties of XTE J1650–500 in the hard state over a large range of luminosities to ADAF/SD, magnetic corona, and jet-dominated models. Although the ADAF/SD model provides natural explanations for the X-ray spectral and timing properties, these models do not include outflows and thus cannot provide a complete description of the system. In addition, ADAFs may be convectively unstable as mentioned above. The magnetic corona model, which does incorporate a jet, predicts a similar spectral evolution to the ADAF/SD model, but it it unclear if it can explain the large change in the Lorentzian half-width that we see for XTE J1650–500. The jet-dominated model can explain X-ray/radio correlations that are seen for several black hole sources without varying $`z_{sh}`$; however, it is doubtful that the model can reproduce the observed X-ray timing properties with a constant value of $`z_{sh}`$.

Regardless of which model is correct, the fact that the general shape of the energy and power spectra remain the same while the parameters evolve indicates that the X-ray emitting region of the system likely keeps the same overall structure even though the parameter changes indicate that some physical properties change, such as the inner radius of the accretion disk, the size of the corona, or the shock location, considering the models described above. It is clear that the source retains its hard state X-ray properties down to a luminosity of $`1.5\times 10^{34}`$ erg s-1, and we do not see any evidence for a sharp transition to the quiescent state. Furthermore, our results for XTE J1650–500 combined with the determination that the optical power spectrum for A 0620–00 in quiescence exhibits band limited noise with an even lower characteristic frequency (Hynes et al., 2003) suggests that black hole transients keep the same overall structure into quiescence. Thus, these results favor the hypothesis that the quiescent state is simply a low luminosity version of the hard state.

JAT acknowledges useful discussions with M. A. Nowak, J. Wilms, and W. A. Heindl. We thank M. Garcia and H. Tananbaum for assistance in scheduling the Target of Opportunity observations with Chandra. We would like to thank all scientists who contributed to the Tübingen Timing Tools. JAT acknowledges partial support from NASA grant NAG5-13055 and Chandra award number GO2-3056X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073. EK acknowledges partial support of TÜBṪAK. PK acknowledges partial support from NASA grant NAG5-7405.

APPENDIX A

We performed simulations to study the effects of pile-up on the observation 3 ACIS power spectrum. The simulations have 3 main goals: 1. To understand why the Poisson noise level for the actual Leahy normalized ACIS power spectrum is 1.77; 2. To determine how much the fractional rms amplitude is reduced by pile-up; and 3. To check that pile-up does not alter the overall shape of the power spectrum. Our goals were accomplished via 3 simulations.

Simulation #1: We produced a simulated Poisson noise light curve with the same time resolution as the actual data (but with much longer duration to improve the statistics). According to the binomial probabilities, we produced a piled-up light curve. The 2 adjustable parameters are the input count rate (before pile-up) and the probability that any pair of events occurring in the same time bin will pile-up. We adjusted these parameters to match the actual ACIS observation as shown in Figure 7. We produced a Leahy normalized power spectrum using the simulated light curve, and the measured Poisson level is 1.78.

Simulation #2: In addition to the deadtime caused by pile-up, there is a fixed, instrumental deadtime that is independent of the event rate and is set by the CCD readout time of 41.04 milliseconds. To determine the effect of the readout time, we simulated a light curve where there is a probability (equal to the readout time divided by the 0.4 s exposure time per frame) that a photon that hits the detector will not be counted. We produced a Leahy normalized power spectrum using the simulated light curve, and the Poisson level is 1.97.

Simulation #3: Here, we made simulated light curves similar to those used for simulation #1 except that we included white noise in excess of the Poisson noise. We produced a light curve free of pile-up and also a light curve with pile-up and made rms normalized power spectra for each. The fractional rms amplitude changes from 19% to 18% when pile-up effects are included. The shapes of the power spectra with and without pile-up are not significantly different.

In summary, these simulations demonstrate that the low Poisson noise level measured for the actual data is mainly due to pile-up effects, and there is also a small drop due to the CCD readout time. For noise levels similar to those we measure for observation 3 ($``$18%), pile-up only causes the fractional rms amplitude to drop by about 1%, and the pile-up does not change the shape of the power spectrum.

APPENDIX B

Here, we describe how we corrected the observation 3 ACIS response matrix for the energy dependence of the Chandra PSF. We generated this correction using the web-based Chandra calibration tool “ChaRT” (Chandra Ray Tracing). For 13 different energies from 0.2 to 10 keV, we used ChaRT to simulate a monochromatic source at the position of XTE J1650–500, $`2^{}.7`$ off-axis. For each energy, we determined the fraction of counts in the same elliptical annulus that we used to extract the XTE J1650–500 spectrum (see text), and the fractions are shown in Figure 8. We then used a spline interpolation to determine the fraction for each energy bin in the response matrix, and we applied the interpolated values to the response matrix. While this correction is necessary, we note that, in this case, it only changes the measurement of the power-law index from $`\mathrm{\Gamma }=1.79\pm 0.13`$ to $`\mathrm{\Gamma }=1.93\pm 0.13`$. Although we do not have a method of accurately determining the systematic error on this correction, even if the correction is only good to 10-20%, the systematic error is much smaller than the statistical error.


text:
# Transverse momentum spectra of hadrons in pp, pA and AA collisions

## I INTRODUCTION

In high energy heavy-ion collisions it is observed that the average transverse momenta of the produced hadrons depends strongly on the mass of hadrons tflow ; bdm . Further it is found to be substantially larger than in nucleon-nucleon collisions at a given energy. One of the possible physical effect responsible for this is transverse flow heinz . The heavier the mass of the particle more is the gain in momentum and hence more transverse flow. This has been shown by extracting the effective temperature (effective because it includes the true freeze-out temperature and the transverse flow) from transverse momentum ($`p_T`$) distribution of various hadrons. It is observed that protons have a larger effective temperature than kaons, which in turn has a larger effective temperature than pions. The success of hydrodynamic models in explaining the transverse momentum distributions of various hadrons at least upto lower momentum, is an indication of presence of radial flow radial . Fig. 1 depicts this effect for AA collisions at SPS energies. However similar effect is also observed in nucleon-nucleus (pA) collisions (also shown in Fig. 1). Hence broadening of $`p_T`$ even in pA collisions, where one does not expect any collective effect like radial flow, indicates there exists a “normal” $`p_T`$ broadening observed in all heavy-ion reactions. This needs to be understood properly before making any quantitative conclusions regarding the radial flow. Such broadening of $`p_T`$ is thought to arise due to successive collisions in nuclear reactions. It is basically an initial state collision broadening of $`p_T`$. This effect in high $`p_T`$ (typically above 2 GeV) is referred to as Cronin effect cronin . Apart from the commonly produced hadrons, the Drell-Yan dileptons and quarkonium resonances also show in pA collisions a broadening of $`p_T`$ compared to that observed in pp collisions. The virtual photons in Drell-Yan production do not interact in the nuclear medium and the quarkonia leave the medium before the collective effects develop. This makes the case for the $`p_T`$ broadening as an initial state effect stronger leonidov ; satz .

In leonidov ; jane , it has been shown that such an effect is able to explain the transverse mass ($`m_T`$) distribution of produced hadrons at SPS energies. Thereby questioning the presence of true collective effects (such as radial flow) in relativistic heavy ion collisions. The aim of the present work is to have a detail study of the effect of $`p_T`$ gain through successive collisions in nuclear reaction within the framework of a random walk model leonidov ; jane  (discussed in the next section). In this paper we systematically analyze the available pp, pA and AA data and come to the conclusion that collective effects does exists in relativistic heavy-ion collisions.

The paper is organized as follows. In the next section we discuss the random walk model. In section III we fix the freeze-out temperature for AA collisions from the inverse slope of the $`p_T`$ spectra from the pp collisions. Section IV deals with the estimation of the gain in average transverse momentum per collision from pA collision data. We compare the results of the random walk model calculation with the SPS AA collision data in section V. In section VI we give a discussion of the results. Finally a summary of the work is given in section VII.

## II RANDOM WALK MODEL

Consider the nucleus-nucleus collision to be a superposition of nucleon-nucleon collision. In this picture, let us assume that in each successive interactions of a nuclear collision one creates a fire ball just like that formed in a nucleon-nucleon collision. If a nucleon starts with zero $`p_T`$, after the first collision the next one will generally occur at some non-vanishing transverse velocity. Thus there is a gain in transverse momentum through successive collisions. The propagation of nucleon through successive collisions is assumed to follow a $`randomwalk`$ pattern. It is of interest to find if such a model can explain the $`p_T`$ spectra of produced hadrons in AA collisions and more specifically if it can account for the observed $`p_T`$ broadening with increase in mass of the hadrons. The random walk pattern is modeled by a Gaussian leonidov  in the present calculation as,

$$f_{pA}(\rho )=\left[\frac{4}{\pi \delta _{pA}^2}\right]^{1/2}\mathrm{exp}(\rho ^2/\delta _{pA}^2),$$
(1)
where $`\rho `$ is the transverse rapidity and

$$\delta _{pA}^2=(N_A1)\delta ^2,$$
(2)
denotes the kick per collision $`\delta `$ as determined from pA interactions. $`N_A`$ is the number of nucleons which the incident proton encounters on its path through the target nucleus. It is given by leonidov

$$N_A(3/4)(2\pi r_0^2R_A)n_0,$$
(3)
where $`r_00.8`$ fm is the nuclear radius, $`R_A`$ (= 1.12 $`A^{1/3}`$ fm) is the radius of the nucleus and $`n_0`$ = 0.17 $`fm^3`$ is the nuclear density.

The corresponding distribution for an $`AB`$ collision can be shown to of the form,

$$f_{AB}(\rho )=\left[\frac{4}{\pi \delta _{AB}^2}\right]^{1/2}\mathrm{exp}(\rho ^2/\delta _{AB}^2),$$
(4)
with

$$\delta _{AB}^2=(N_A+N_B2)\delta ^2.$$
(5)
The final expression for the transverse mass distribution is given as jane

$$\left(\frac{dN}{dym_Tdm_T}\right)_{y=0}=\frac{gV}{2\pi ^2}\left[\frac{4}{\pi \delta _{AB}^2}\right]^{1/2}𝑑\rho \mathrm{exp}(\rho ^2/\delta _{AB}^2)m_TI_0\left(\frac{p_T\mathrm{sinh}\rho }{T}\right)K_1\left(\frac{m_T\mathrm{cosh}\rho }{T}\right).$$
(6)
Where $`I_0`$ and $`K_1`$ are modified Bessel functions, T is the freeze-out temperature and $`m_T=\sqrt{p_T^2+m^2}`$ is the transverse mass.

It should be noted that the volume in the Eqn. 6 refers to the volume of the system as observed in a $`pp`$ collision, since each collision in the random walk produces a $`pp`$ type of fireball. If we now introduce a boost-invariant distribution of fireballs along the longitudinal rapidity axis, we finally obtain by integrating over the fireball distributions

$`\left({\displaystyle \frac{dN}{dym_Tdm_T}}\right)_{y=0}`$ $`=`$ $`{\displaystyle \frac{gV}{2\pi ^2}}\left[{\displaystyle \frac{4}{\pi \delta _{AB}^2}}\right]^{1/2}{\displaystyle 𝑑\rho \mathrm{exp}(\rho ^2/\delta _{AB}^2)}`$ (7)
$`{\displaystyle _{Y_L}^{Y_L}}𝑑Ym_T\mathrm{cosh}YI_0\left({\displaystyle \frac{p_T\mathrm{sinh}\rho }{T}}\right)K_1\left({\displaystyle \frac{m_T\mathrm{cosh}Y\mathrm{cosh}\rho }{T}}\right).`$
The above expression has two parameters, T and $`\delta `$. The temperature is assumed to be same for the whole fireball. The temperature is fixed from the pp collisions and $`\delta `$ is calculated from pA collisions as mentioned before. Fixing these two parameters from pp and pA collisions we will attempt to describe the data for AA collisions.

## III TRANSVERSE MOMENTUM SPECTRA IN pp COLLISIONS

As discussed in the previous section, we would fix the freeze-out temperature for the AA collisions from the inverse slope of the hadronic transverse momentum spectra originating from pp collisions. At high energies the transverse momentum spectra of produced hadrons can be described by the following expression,

$$(\frac{d^2N}{dydp_T^2})^{pp}const.m_TK_1(m_T/T).$$
(8)
Using the expression given in Eq. 9, we fitted the experimental data of pp collisions of pions, kaons and protons at different centre-of-mass (CMS) energies varying from 23 GeV to 63 GeV isr . The results are shown in Fig. 2. The extracted temperature for the various hadron species are shown in Fig. 3. The error bars shown corresponds to the error on the slope parameter due to the fitting procedure adopted here. One observes that there is very little variation in T for all the produced hadrons in the range of CMS energy for which data is available and analyzed. We will consider two values of temperature for our study of AA data, one corresponding to 140 MeV and the other corresponding to 150 MeV.

## IV TRANSVERSE MOMENTUM SPECTRA IN pA COLLISIONS

Having fixed the the freeze-out temperature for AA collisions from the inverse slope of $`p_T`$ spectra of produced hadrons in pp collisions, we now try to get the average transverse rapidity shift per collision, $`\delta `$. This we will determine empirically from transverse mass spectra of the produced hadrons in pA collisions. To get a proper understanding of this parameter, we will study pA data for various available beam energies, different species of produced hadrons and several targets. Once this parameter is fixed, along with the temperature fixed from pp collisions, we will be able to predict the initial state transverse momentum broadening in AA collisions as has been discussed in the section II. The results will be presented in the next section.

The transverse momentum spectra for different species of produced hadrons for pA collisions can be obtained using the Eqn. 7. Where $`\delta _{pA}`$ is the only parameter, with T being fixed from pp collisions. The results are shown in Figs. 4 and 5. Fig. 4 shows the comparison of transverse mass distribution for pion, kaon and proton at 14.6 GeV for various targets e802  with those obtained from model calculations. One finds that the data is reasonably well explained by the model calculations. The values of $`\delta _{pA}`$ obtained for pion, kaon and proton are 0.05, 0.25 and 0.25 respectively. The corresponding $`\delta `$ values for the p+Cu collisions are obtained using Eqn. 2 to be 0.04, 0.22, 0.22 for pions, kaons and protons respectively. Also shown in Fig. 4 is the transverse mass spectra for pions at 14.6 GeV for various targets. We find that there is very little atomic mass dependence of $`\delta _{pA}`$. It was found to be 0.05 for the various targets shown in the figure. However one observes from Eqn. 5, the average transverse rapidity shift per collision, $`\delta `$ decreases with increase in atomic mass as we go from beryllium to gold target. Fig. 5 shows the transverse mass spectra for pion, kaon and proton for p+Pb collisions at 450 AGeV na44 . The solid lines are the result of model calculation taking T = 150 MeV. The resulting $`\delta _{pA}`$, which satisfactorily fits the data as shown, turns out to be of the order of 0.45 for three hadronic species. This is a clear increase from the lower energy value, indicating that $`\delta _{pA}`$ has an energy dependence. The corresponding value $`\delta _{pA}`$ for a T = 140 MeV (dotted lines) lies between 0.35-0.4.

For the description of the AA data at SPS, we will consider two values of temperature, of 140 MeV and 150 MeV fixed from the pp data. Since the collisions at SPS are lead on lead target, we will use $`\delta _{pPb}`$ values as discussed above for pPb collisions at 450 AGeV. We will also consider one case with $`\delta _{pBe}`$ at 14.6 AGeV in order to get an idea of the effect at lower energies.

## V TRANSVERSE MOMENTUM SPECTRA IN AA COLLISIONS

Having fixed the temperature from the pp spectra and knowing the value of $`\delta `$ from the available pA data, we now try to see if the random walk model is able to explain the transverse mass spectra of pions, kaons and protons at SPS energies. The results of the model calculation using Eqn. 7 along with the available data na49  are shown in Fig. 6.

The results for pions are shown in top panel of Fig. 6. The lines corresponds to the random walk model calculations for two different temperatures of 140 MeV and 150 MeV. One observes that the model is able to explain the pion spectra satisfactorily over the energy range of 40 AGeV to 158 AGeV for both 140 and 150 MeV temperatures. The $`\chi ^2/ndf`$ being less than 1.0 for all the cases. Taking a $`\delta `$ value corresponding to pA results from AGS energies (14.3 AGeV) fails to explain the data. This is shown as dotted line for 40 AGeV beam energy. Similar results are obtained for higher beam energies and not shown in the figure for clarity of presentation.It may be mentioned that the absolute normalisation is adjusted to get the best possible fit.

For the kaons (middle panel of Fig. 6), one finds that the model fails to explain the lower transverse momentum region. Further one notices that the $`\chi ^2/ndf`$ worsens as we go to higher beam energies. It varies from 3.5 at 40AGeV to 10 at 158 AGeV. There is not much difference between the results for the model calculation for temperature of 140 and 150 MeV, hence the spectra corresponding to 140 MeV is not shown in the figure. The model also fails to explain the observed transverse mass spectra of protons at both lower and higher transverse momentums. The results for two different temperatures of 140 MeV and 150 MeV are shown in the bottom panel of Fig. 6. The $`\chi ^2/ndf`$ for the three cases lies between 3.5 to 7.5.

## VI DISCUSSION

The results indicate the following : (i) the random walk model fails to explain the transverse mass spectra as the mass of the hadron increases where the effect of possible transverse flow will be more and (ii) the $`\chi ^2/ndf`$ indicates that the disagreement of the model calculation with data for kaons and protons increases with increase in beam energy. In this section we discuss the limitations of the model and try to estimate the relative contribution of initial state $`p_T`$ broadening and transverse flow.

### VI.1 Random walk pattern

The random walk pattern was chosen to be a Gaussian (Eqn. 1). In principle this can be any other type of statistical distribution, such as a Lorentzian. We have studied the sensitivity of the result to such a distribution. We find that the p+Pb data for pion, kaon and proton is well explained by a $`\delta _{pPb}`$ of 0.2, 0.15 and 0.35 respectively. The values are lower compared to those obtained considering a Gaussian random walk pattern. Extrapolating these values to Pb+Pb collisions and taking the value of $`T`$ to be 150 MeV, we obtain the $`m_T`$ spectra for the available SPS data at the highest energy. The results are shown in Fig. 7. We observe that the model fails to explain the data satisfactorily, specially at low $`m_T`$ for kaons and protons. However it agrees fairly well with the observed pion spectra. It seems that the results are insensitive to the choice of random walk pattern.

### VI.2 Parameters of random walk model

One of the disadvantage of the present work is the absence of the pp and pA data in literature corresponding to the exact SPS energies for which AA collision data is available. However one must mention that the near constant temperature exhibited by the available pp data for a wide range of CMS energy more or less fixes the temperature parameter of the model. Contrary to this one observes a considerable variation of $`\delta `$ with beam energy. In view of these, it is obvious to ask, if the SPS kaon and proton transverse mass spectra can be explained within the framework of random walk model for arbitrary values of temperature and $`\delta `$.

The results are shown in Fig. 8. One observes that the kaon spectra for the three different beam energies are well explained by the model taking a common temperature of 235 MeV and $`\delta _{PbPb}`$ of 0.25. The $`\chi ^2/ndf`$ is less than 1.0 for all the three cases. It may be mentioned that this temperature although very high compared to that obtained from pp (section III), is very close to the temperature obtained taking radial flow into account na49 . For the protons the spectra is satisfactorily reproduced over a large range of transverse momentum. For this we needed a common $`\delta _{PbPb}`$ values of 0.15 but the temperature for 40 and 80 AGeV is 260 MeV, while that for 158 AGeV energy it is 360 MeV.

### VI.3 Transverse flow

The failure of random walk model to fully explain the transverse momentum spectra of produced hadrons indicate that initial state $`p_T`$ broadening alone cannot account for experimentally measured transverse mass spectra. One observes that the model fails to explain the lower momentum part for kaons and protons. This may be attributed to the lack of chemical equilibrium, which can be accounted for by introducing a chemical potential for the hadrons But the formulation of the model is such that it affects only the normalization and not the slope. The other possibility is that, there may be both initial state $`p_T`$ broadening and collective effect like transverse flow. In other words the resultant spectra of produced hadrons is the effect of both the above processes. In order to look at this possibility, we express the transverse mass spectra of produced hadrons as;

$$\left(\frac{dN}{dym_Tdm_T}\right)_{total}=C_1\left(\frac{dN}{dym_Tdm_T}\right)_{random}+C_2\left(\frac{dN}{dym_Tdm_T}\right)_{flow}$$
(9)
where, the contribution from flow flow ; na49  is defined as, $`\left(\frac{dN}{dym_Tdm_T}\right)_{flow}m_TI_0\left(\frac{p_T\mathrm{sinh}\rho }{T}\right)K_1\left(\frac{m_T\mathrm{cosh}\rho }{T}\right)`$, mean transverse velocity $`\beta _T`$ is defined through the relation $`\rho =atanh\beta _T`$ and T is the freeze-out temperature here taken to be 150 MeV as seen from the slope of pp spectra.

We tried to see if such a parametrization can explain the SPS data at 158 AGeV beam energy and get the values of $`C_1`$ and $`C_2`$ which will tells the relative contribution of the two effects. The results are shown in the Fig. 9. Since pion data is very well explained by the random walk model, we tried not to fit it with the above parametrization. The aim here is to see how much contribution of transverse flow in addition to initial state $`p_T`$ broadening through the random walk model can explain the data. From the figure we find that the results for the kaons and protons shows that such a parametrization works. We have taken a value of $`\beta _T`$ of 0.45 for kaons and 0.48 for protons. The values of $`C_1`$ and $`C_2`$ are 0.5 for kaons and protons, $`C_1`$ = 0.2 and $`C_2`$ = 0.8. The results indicate the flow effects dominates with increase in hadron mass.

## VII SUMMARY

In summary, a systematic study has been carried out to understand the transverse mass spectra of the produced hadrons in nucleus-nucleus collisions within the framework of a simple random walk model. The model is based on the assumption that the nucleus-nucleus collision is a superposition of nucleon-nucleon collision where for each successive interactions of a nuclear collision one creates a fire ball. The temperature of the fire ball is fixed from the available nucleon-nucleon collision data. There is a gain in transverse momentum through successive collisions, the propagation of which is assumed to follow a random walk pattern. The average transverse rapidity shift per collision or the gain in transverse momentum per collision is obtained from the nucleon-nucleus collision data. Having fixed the two parameters of the model, we then apply it to nucleus-nucleus collisions. It is observed that, although the model fairly well explains the transverse mass spectra of pions at SPS energies, it fails to do so for the higher mass hadrons like kaons and protons. We find the choice of different distribution of random walk pattern yields the same result. This indicates the presence of true collective effects in the nucleus-nucleus collisions. However it seems a considerable portion of the gain in transverse momentum comes from initial state effects. Considering both the initial state $`p_T`$ broadening and the transverse flow effect, we are able to explain the SPS data. We find that the contribution of transverse flow increases with increase in hadron mass. The avability of systematic data from pp to pA to AA collision at a fixed CMS energy would help in quantifying the effect of initial state broadening of transverse momentum and hence the true amount of collective effect.

###### Acknowledgements.

I am grateful to the Board of Research on Nuclear Science and Department of Atomic Energy, Government of India for financial support. I would like to thank Micheal Murray and Subrata Bhattacharyya for providing me the experimental data. I would like to thank Jan-e Alam for critical reading of the manuscript and for many helpful discussions.