text:
# Measuring Neutron Star Radii with Gravitational Wave Detectors

## Abstract

Coalescing binary neutron stars (NS) are expected to be an important source of gravitational waves (GW) detectable by laser interferometers. We present here a simple method for determining the compactness ratio $`M/R`$ of NS based on the observed deviation of the GW energy spectrum from point-mass behavior at the end of an inspiral event. Our method is based on the properties of quasi-equilibrium binary NS sequences and does not require the computation of the full GW signal $`h(t)`$. Combined with the measurement of the NS masses from the GW signal during inspiral, the determination of $`M/R`$ will allow very strong constraints to be placed on the equation of state of nuclear matter at high densities.

Coalescing compact binaries containing two neutron stars (NS) are among the most important sources of gravitational waves (GW) for LIGO LIGO , VIRGO VIRGO , and other laser interferometers. Should the inspiral of such a binary be detected, the frequency evolution of the GW signal will immediately yield the system’s “chirp mass” $`M_{\mathrm{ch}}\mu ^{3/5}^{2/5}`$, where $`\mu `$ and $``$ are the reduced and total mass of the binary, respectively. Higher-order post-Newtonian effects on the phase evolution of the signal also allow for the determination of the reduced mass $`\mu `$, and thus the individual masses $`M_1`$ and $`M_2`$ of the two NS CF . The determination of the NS radii in addition to their masses would yield important information about the equation of state (EOS) at nuclear densities, and could even indicate the presence of more exotic phases, such as strange quark matter instead of ordinary nuclear matter Strange . The GW signal of a coalescing binary could yield such information but this is limited in two different ways. During the slow inspiral phase at large separations, i.e., low frequencies ($`f1`$kHz), the stars behave like point masses, and finite-size effects are not expected to leave any signature in the GW signal CF ; Koch ; Lai . During the final hydrodynamic merger, characteristic frequencies of GW emission could yield important information about the fluid EOS Zhu ; FR ; Oech2 , but these frequencies are expected to lie well above the photon shot noise limit of current interferometers ($`f1`$kHz). Thus, it is only during the last few orbits of the inspiral, just prior to merger, that we can hope to see the imprint of the NS radii on a measurable GW signal (with $`f1`$kHz).

Several groups have studied this terminal phase of inspiral by constructing quasi-equilibrium sequences of close NS binaries in the conformal flatness approximation of general relativity (GR) Baumg ; Marro ; Uryu00 ; Gourg01 . In this approximation, it is assumed that the binary system evolves along a sequence of appropriately constructed equilibrium states with decreasing binary separation as energy is radiated away in GW. From the binary equilibrium energy curve $`E_{\mathrm{eq}}(r)`$, which gives the total system energy as a function of binary separation $`r`$, and the GW luminosity $`L_{\mathrm{GW}}`$, one can derive the radial infall rate as $`v_r=L_{\mathrm{GW}}(dE_{\mathrm{eq}}/dr)^1`$. With $`v_r=dr/dt`$ known, this provides the time evolution along the equilibrium sequence and the GW signal $`h(t)`$. This approach should remain accurate as long as the radial infall timescale $`r/v_r`$ is longer than the dynamical timescale of the system, i.e., until the point where dynamical instability sets in, and the two stars plunge inward rapidly and merge.

Unfortunately, calculating the correct GW luminosity for a given matter configuration in GR is an extremely difficult task. Different approaches have required either time integration of the full non-linear equations of GR Duez , or the solution of a complicated wave equation for terms representing the spherical harmonic expansion of the GW metric perturbation SU . The great complexity of these approaches is in stark contrast with the simplicity of the quasi-equilibrium approximation. However, we point out here that the GW energy spectrum $`dE_{\mathrm{GW}}/df`$ can be calculated directly and very simply from the equilibrium energy curve, independent of any knowledge about the GW luminosity. Indeed, by definition of the quasi-equilibrium approximation, the energy decrease $`dE_{\mathrm{eq}}`$ between two neighboring binary configurations along the sequence is equal to the energy $`dE_{\mathrm{GW}}`$ radiated away as the wave frequency sweeps up by $`df`$, where the GW frequency is twice the orbital frequency, $`f=2f_{\mathrm{orb}}`$. Thus, one should simply compute the total energy $`E_{\mathrm{eq}}`$ as a function of frequency $`f`$ along the equilibrium sequence, and the GW energy spectrum is then given by the derivative $`dE_{\mathrm{GW}}/df=dE_{\mathrm{eq}}/df`$. As a trivial example, consider the inspiral of two point masses in the Newtonian limit, where we have $`E_{\mathrm{eq}}r^1`$ and $`fr^{3/2}`$. It follows that $`E_{\mathrm{eq}}f^{2/3}`$ and thus $`dE_{\mathrm{GW}}/dff^{1/3}`$, a well-known result. In addition to the assumptions underlying the quasi-equilibrium approximation, the validity of this simple approach relies on the additional assumption that the GW emission during the later merger phase has no effect on the energy spectrum at lower frequencies. Indeed, this has been demonstrated in numerical hydrodynamic calculations of binary mergers Zhu ; FR , which show a clear separation between the inspiral and merger components of the emission in frequency space.

We have investigated the properties of the GW emission during the final phase of binary NS inspiral using new, highly accurate equilibrium sequences calculated with the Meudon code Lorene . The formalism is based on the conformally flat approximation, which is expected to yield accurate equilibrium matter configurations for this phase Baumg ; Marro ; Uryu00 ; Gourg01 . The resulting five non-linear, coupled elliptic equations are solved using a multi-domain spectral method Bonaz ; Grand01 . This approach has already been used successfully in various astrophysical applications Gourg01 ; Gourg99 ; Grand02 . Typically, the computed fields satisfy the constraints of full GR to within $`1\%`$ MMW . The code has been improved recently, especially with regard to the treatment of the external compactified domain, and the numerical error in computing equilibrium configurations, measured in terms of how well the virial theorem is satisfied, is of order one part in $`10^5`$. The equilibrium sequences presented here are the natural extension of the work already published in Gourg01 ; TaniII . They will be discussed in more detail in TaniIV . Here we only show the variation of the ADM mass of the system (total binary mass-energy $`c^2E_{\mathrm{eq}}`$) and the GW frequency (twice the orbital frequency), which are sufficient for our purposes to determine the GW energy spectrum.

The equilibrium configurations have been calculated for irrotational binaries, i.e., assuming that the fluid has zero vorticity in the inertial frame. Indeed, the NS should be spinning slowly at large separations and the viscosity of NS matter is too small for tidal spin-up to become significant on the coalescence timescale Koch ; BC . Based on the current set of well-measured NS masses in relativistic binary radio pulsars, it is expected that all NS in coalescing binaries will have masses $`M1.35M_{}`$ TC . Hence, for simplicity, we consider only equal-mass binaries where $`M_1=M_2=1.35M_{}`$. Also for simplicity, we model the NS EOS with a simple polytropic form $`P=K\rho ^\mathrm{\Gamma }`$, where $`P`$ is the pressure and $`\rho `$ the rest-mass density. The constant $`K`$ represents the overall compressibility of the matter and largely sets the value of the stellar radius for a given mass, while the adiabatic exponent $`\mathrm{\Gamma }`$ measures the stiffness of the EOS and affects the degree of central concentration of the NS interior. Based on our experience with hydrodynamic calculations FR , we expect that the GW energy spectrum just prior to merger is determined primarily by the NS radius through $`K`$, with relatively little sensitivity to $`\mathrm{\Gamma }`$. For this reason, in this initial study, we allow $`K`$, and therefore also the stellar radius $`R`$, to vary for different NS models, but we set $`\mathrm{\Gamma }=2`$ for all models, as this value fits well most published NS EOS (see, e.g., Lattimer  and references therein). Specifically, we consider NS models with compactness ratios $`M/R=`$0.12, 0.14, 0.16, and 0.18 (setting $`G=c=1`$), where $`M`$ is the ADM (gravitational) mass measured by an observer at infinity for a single isolated NS, and $`R`$ is the circumferential radius of the NS. For $`M=1.35M_{}`$, the corresponding radii are $`R=`$16.6, 14.2, 12.4, and 11.1 km, respectively, spanning the range of values for NS radii calculated from various physical EOS. Note, however, that our results are unchanged under the rescaling given by $`R^{}=\kappa R`$, $`M^{}=\kappa M`$, $`f^{}=f/\kappa `$, for any constant $`\kappa `$.

For each NS model, about 12 binary equilibrium configurations are computed with decreasing separations, until a cusp develops on the NS surface. Equilibrium configurations for smaller separations do not exist. To each sequence we fit a curve of the form

$$(f)=2.7M_{}k_\mathrm{N}f^{2/3}+k_1f+k_2f^2$$
(1)
to represent the variation of total mass-energy as a function of GW frequency. The first term gives the total gravitational mass of the system at infinite separation, while the second term represents the Newtonian point-mass behavior, with $`k_\mathrm{N}=2^{4/3}\pi ^{2/3}G^{2/3}M^{5/3}=4.056\times 10^4M_{}\mathrm{Hz}^{2/3}`$. The term $`f`$ was introduced heuristically to represent the lowest-order post-Newtonian and finite-size corrections to the point-mass behavior at intermediate frequencies. The term $`f^2`$ represents the tidal interaction energy, which causes the equilibrium energy curve to flatten at high frequencies. Our best (least-squares) fit values of $`k_1`$ and $`k_2`$ for each sequence are listed in Table 1 and the results are illustrated in Fig. 1. The asterisks show the data points along each sequence, with a typical error between the data points and the fit of $`\delta 10^4M_{}`$. We find in all cases that $`k_2`$ is positive, as we would expect: tidal deformations and relativistic gravitational effects increase the equilibrium energy LRS3 ; LomRS . We note that none of the equilibrium curves shows evidence of an energy minimum, which would have implied the onset of dynamical instability Lai ; TaniII ; LRS3 ; LRS2 .

Computing the GW energy spectrum for each model now only requires differentiating the fitted curves with respect to frequency. The results are shown in Fig. 2. In each case, we see a characteristic frequency range where the spectrum plunges rapidly below the extrapolation of the low-frequency result. This corresponds to the flattening of the energy curves in response to the growing tidal interaction and PN effects. Also shown is the energy spectrum of a 3PN, irrotational, point-mass binary, computed according to the results found in Ref. Damour , which closely tracks our most compact model, indicating that the differences we see in the energy spectra result from finite-size effects associated with the NS radius. To quantify the deviations from the Newtonian case, we define a set of break frequencies, at which the energy spectrum has dropped by some factor below the point-mass result. The values we find for $`f_{10},f_{25}`$, and $`f_{50}`$, where $`dE_{\mathrm{GW}}/df`$ has dropped by $`10\%,25\%`$, and $`50\%`$, respectively, are listed in the last three columns of Table 1. We see that all these characteristic frequencies lie within the frequency range accessible by LIGO-type detectors, with perhaps the exception of $`f_{50}`$ for the more compact sequences. Note that the calculation of $`f_{50}`$ values requires extrapolating the equilibrium energy curves slightly beyond the last equilibrium model (where a cusp develops), and may therefore be less reliable. However, we find that $`f_{50}`$ has a particularly steep, quasi-linear dependence on the NS compactness, given by $`f_{50}[10^4(M/R)460]\mathrm{Hz}`$ within the range of NS radii we considered. For comparison, $`f_{25}`$ values can be determined safely within the quasi-equilibrium approximation, and the sensitivity on compactness is only slightly reduced, with $`f_{25}[5000(M/R)85]\mathrm{Hz}`$.

The best definition of the break frequency will be a tradeoff between higher signal-to-noise ratio at lower frequencies LIGOnote , and greater ease of discrimination between different EOS at the higher frequencies. In addition, the quasi-equilibrium approximation is expected to be most accurate at lower frequencies, where the inspiral rate is slower. However, we doubt that this could become a major issue: if we adopt, for simplicity, the point-mass formula for the GW luminosity, and compute the corresponding radial infall velocity along our equilibrium sequences, we find that $`v_r`$ never exceeds $`5\%`$ of the orbital velocity, even at the point where we define $`f_{50}`$ (the corresponding fraction at the point where we define $`f_{25}`$ is about 2%) NOTEvr . Ultimately, the optimum choice should be determined by data analysts, taking into account the accuracy with which the break frequencies can be extracted using matched filtering or other techniques CF ; Match . Preliminary studies of this problem have already been performed for both broad-band and narrow-band interferometer configurations Vallis ; Hughes . Defining a precise break frequency may not even be necessary. Instead, the GW inspiral templates could be terminated at high frequency in a manner that reproduces the energy spectrum given by a simple analytic form, such as our Eq. 1. The free parameters $`k_1`$ and $`k_2`$ could then be measured experimentally and compared directly to the predictions of various realistic nuclear EOS used in computing binary equilibrium sequences. In future work, we plan to compute such sequences, and the corresponding energy curves, for a wide variety of published, realistic NS EOS.

We are grateful to Kip Thorne for originally pointing out the importance of the break frequencies in GW spectra. This work was supported by NSF Grants PHY-0070918, PHY-0121420, and PHY-0133425.


format: text

text:
# Does parton saturation at high density explain hadron multiplicities at RHIC ?

## Abstract

We discuss the recent claim that hadron multiplicities measured at RHIC energies are directly described in terms of gluon degrees of freedom fixed from the initial conditions of central heavy ion collisions. The argument is based on the parton saturation scenario expected to be valid at high parton densities and on the assumption of conserved gluon number. Alternatively we conjecture that ”bottom-up” equilibration before hadronization modifies this picture, due to nonconservation of the number of gluons.

At RHIC and LHC energies in the central rapidity region of heavy ion collisions a high density of energy is deposited mainly in the form of gluons. In recent papers by Kharzeev et al.  (hereafter referred to as KLN) it is assumed that the initial gluon density determines the number of produced (charged) particles, i.e. that there is a direct correspondence between the number of partons in the initial state and the number of particles in the final state. The argument is formulated in the framework of the saturation scenario , which determines the initial gluon distribution inside the colliding nuclei, and where the characteristic momentum scale is the hard scale $`Q_s>>\mathrm{\Lambda }_{QCD}`$.

Based on this scenario the initial gluon multiplicity immediately after the high energy nuclear collisions can be calculated in the McLerran-Venugopalan model . At $`\tau _01/Q_s`$ this initial gluon density is ,

$$n_{hard}(\tau _0)=c\frac{(N_c^21)Q_s^3}{4\pi ^2N_c\alpha _s(Q_s\tau _0)}.$$
(1)
In  this relation is used to obtain the hadron multiplicities in A-A collisions, which are then compared with the corresponding RHIC data .

In this note we propose a possibly more realistic description, which departs from taking the initial condition as the only one determining ingredient for the final hadronic state. This description is based on the “bottom-up” scenario , which leads to thermalization of the gluons produced after the collision. In contrast to KLN, this scenario stresses the importance of branching processes of gluons which may allow for a short enough equilibration time, at least at high energies. As a consequence the number of gluons is increasing between initial and equilibration times. In this context let us note that already in  a phenomenological factor $`\kappa _{inel}>1`$ is introduced accounting for gluon number changing processes which may occur at late times beyond when the classical approach is applicable.

In the following we investigate in the “bottom-up” picture the resulting hadron multiplicities and compare with the predictions of the KLN model. We discuss in some detail the conceptual differences between these two interesting possibilities.

The arguments of KLN  go through a number of steps, which we first critically review:

i) The number of gluons in the initial state, and at the time when gluons transform into hadrons, is assumed to be equal, which can be true when only $`22`$ processes are taken into account. Only in this circumstance the measured hadron multiplicities are reflecting directly the initial conditions.

ii) In (1) the parameter $`c`$ is a constant linking the number of gluons in the nucleus wave function to the number of gluons which are freed during the collision. It is expected to be of $`O(1)`$ . An (approximate) analytical calculation gives $`c=2\mathrm{ln}21.39`$ . However, extracting $`c`$ from the numerical simulation in  leads to the estimate ***This estimate is based on the relation $`c=4\pi ^2f_N/[(N_c^21)\mathrm{ln}Q_s^2/\mathrm{\Lambda }_{QCD}^2]`$, where $`f_N=0.3`$ .  $`c0.5`$ for $`Q_s=1GeV`$.

iii) The hard saturation scale $`Q_s^2`$ in the case of one nucleus has been determined to be ,

$$Q_s^2(\stackrel{}{s})=\frac{4\pi ^2N_c}{N_c^21}\alpha _s(Q_s^2)xG(x,Q_s^2)\rho _{nucleon}(\stackrel{}{s}),$$
(2)
where $`\stackrel{}{s}`$ is the impact parameter, $`xG(x,Q^2)`$ is the gluon structure function in the nucleon, and $`\rho _{nucleon}`$ is the transverse density of nucleons in the nucleus. Notice that in the above equation an additional numerical multiplicative factor may come as a consequence of the inherent uncertainty in the precise determination of the saturation momentum outside of the McLerran-Venugopalan model . We shall call this factor $`K`$ and introduce it a little later on. KLN generalize (2) to the case of two colliding nuclei. Let us first write Notice in (3), consistently with (2), the factor $`\frac{1}{2}`$ compared to eq. (14) in . This has the consequence of lower values of $`Q_s^2`$ appearing in our discussion of the KLN approach.

$$Q_s^2(\stackrel{}{s},\stackrel{}{b})=\frac{4\pi ^2N_c}{N_c^21}\alpha _s(Q_s^2)xG(x,Q_s^2)\frac{\rho _{part}(\stackrel{}{s},\stackrel{}{b})}{2}.$$
(3)
where $`\rho _{part}(\stackrel{}{s},\stackrel{}{b})`$ is the density of participating nucleons in the transverse plane as a function of the impact parameter $`\stackrel{}{b}`$ of the A-A collision, and of $`\stackrel{}{s}`$ the transverse coordinate of the produced gluon. The factor $`\frac{1}{2}`$ is required by the proper definition of $`Q_s^2`$, relative to one nucleus . Integrating the density with respect to $`\stackrel{}{s}`$ leads to

$$d^2s\rho _{part}(\stackrel{}{s},\stackrel{}{b})=N_{part}(\stackrel{}{b}),$$
(4)
with $`N_{part}(\stackrel{}{b})`$ the number of participants in the A-A collision at fixed $`\stackrel{}{b}`$.

iv) Focusing on the distribution of freed gluons $`dN/d\eta `$ for $`\eta =0`$ at given $`\stackrel{}{b}`$, we can write

$`{\displaystyle \frac{dN}{d\eta }}(\stackrel{}{b})`$ $`=c{\displaystyle \frac{N_c^21}{4\pi ^2N_c}}{\displaystyle d^2s\frac{1}{\alpha _s}Q_s^2(\stackrel{}{s},\stackrel{}{b})}`$ (6)
$`cxG(x,\overline{Q}_s^2){\displaystyle \frac{N_{part}(\stackrel{}{b})}{2}}.`$
The corresponding charged hadron multiplicity for the most central collisions is then obtained as

$$\frac{2}{N_{part}}\frac{dN_{ch}}{d\eta }\frac{1}{3}c\left[\mathrm{ln}\frac{\overline{Q_s^2}}{\mathrm{\Lambda }_{QCD}^2}\right],$$
(7)
using the gluon structure function (9) given below. As reference energy we take $`\sqrt{s}=130GeV`$.

One should remark that in (6) and (7), $`\overline{Q}_s^2`$ shows up as an effective average over the variable $`\stackrel{}{s}`$ in (3). As a first approximation the following is used:

$$\overline{Q_s^2}(\stackrel{}{b})=K\frac{4\pi ^2N_c}{N_c^21}\alpha _s(\overline{Q_s^2})xG(x,\overline{Q_s^2})\frac{\rho _{part}(\stackrel{}{b})}{2}.$$
(8)
v) Information on the gluon structure function is necessary. In the small $`x`$ regime, the main feature is that it increases with $`Q^2`$ at fixed $`x`$. Following KLN , it is reasonable to take

$$xG(x,Q^2)=0.5\mathrm{ln}(\frac{Q^2}{\mathrm{\Lambda }_{QCD}^2}),$$
(9)
with $`\mathrm{\Lambda }_{QCD}=200MeV`$, such that $`xG(x,Q^2)2`$ at $`Q^2=2GeV^2`$ (at $`x=0.02`$).

vi) The strong coupling constant is

$$\alpha _s(Q^2)\frac{1}{\beta _0\mathrm{ln}(\frac{Q^2}{\mathrm{\Lambda }_{QCD}^2})},$$
(10)
with $`\beta _0=(112n_f/3)/4\pi `$ for $`N_c=3`$; we take $`n_f=3`$.

vii) Using (8) together with the above choices for $`\alpha _s`$ and $`xG(x,Q^2)`$, one finds for central $`AuAu`$ collisions at RHIC (at $`\sqrt{s}=130GeV`$), taking for the moment $`K=1`$:

$`\overline{Q_s^2}(\stackrel{}{b}=0)0.63GeV^2`$, using $`\rho _{part}(\stackrel{}{b}=0)=3.06fm^2`$,

as quoted in Table 2 by . For larger values of $`\stackrel{}{b}`$, $`\overline{Q_s^2}`$ becomes even smaller, e.g. for $`b=10fm`$, $`\overline{Q_s^2}0.32GeV^2`$ ! This casts some doubts on the applicability of this model, even for central collisions at RHIC energies.

As already noticed, the missing factor $`\frac{1}{2}`$ in KLN and a larger value used for $`\alpha _s`$, namely $`\alpha _s(Q^2)=0.6`$ at $`Q^2=2GeV^2`$, allows them to get larger values of $`Q_s^2`$, i.e. $`Q_s^22GeV^2`$ for the central collisions.

Taking in the following an optimistic point of view, we use (8) with a multiplicative $`K`$ factor of order $`O(1)`$, explicitely $`K1.6`$, such that for $`\stackrel{}{b}=0`$, $`\overline{Q_s^2}(\stackrel{}{b}=0)=1GeV^2`$. Our use of $`K`$ is to a large extent cosmetic. This factor only appears in calculating $`Q_s^2`$, but does not change (7) or (21). Such a factor will affect the average transverse momentum per produced gluon but not the total number of produced gluons. The difference between $`\overline{Q_s}(\stackrel{}{b}=0)=1GeV`$ and $`\overline{Q_s}(\stackrel{}{b}=0)=0.8GeV`$, corresponding to $`K=1.6`$ and $`K=1`$, respectively, has little effect on the equilibration temperature $`T_{eq}`$ and time $`\tau _{eq}`$ quoted later in this note.

viii) Finally, there is an equality assumed between the number of partons in the final state and the number of observed hadrons (”parton-hadron duality” ).

Let us now turn to the “bottom-up” scenario , which in contrast to the KLN  prescription, does not relate the multiplicities to the initial condition only, but also to the way gluons are thermalized.

In the framework of perturbative QCD the time evolution of the gluonic system, when described by a non-linear Boltzmann equation based on $`22`$ processes, a relatively long time ($`Q_s^1\mathrm{exp}(\alpha _s^{1/2})`$) is required for the approach to kinetic equilibration . Fast and efficient thermalization occurs when inelastic processes, namely gluon splittings $`23`$ are taken into account, and kinetic equilibration occurs much faster at times $`\alpha _s^{13/5}Q_s^1`$ . In this ”bottom-up” scenario the time evolution of the system proceeds through several regimes (Fig. 1), with $`Q_s\tau \alpha _s^{3/2},\alpha _s^{5/2}`$ and $`\alpha _s^{13/5}`$.

The difference between the ”bottom-up” and the KLN picture is schematically illustrated in this Fig. 1. Under the KLN assumption the hard gluons (on the momentum scale $`Q_s`$) are conserved in number, i.e. $`n_{hard}\tau =const`$, and they finally hadronize, after passing through a hydrodynamical stage. The ”bottom-up” scenario is characterized by the fact that hard gluons are degrading, soft ones are formed and start to dominate the system. As a result the interactions of gluons in this kinetic scenario modify strongly the initial gluon spectrum. The gluons are redistributed and thermalizing, such that a quark gluon plasma is formed. The number of gluons, together with the entropy , is increasing with proper time $`\tau `$, such that the ratio is,

$$R=[n_{soft}(\tau )(Q_s\tau )]|_{\tau _{eq}}/[n_{hard}(\tau )(Q_s\tau )]|_{\tau _0}\alpha _s^{2/5}1.$$
(11)
The following processes are expected to be present (at RHIC and higher energies in the central region of pseudorapidity $`\eta 1`$): at $`Q_s\tau 1`$ saturated hard gluons $``$ elastic scatterings and branching/production of soft gluons $``$ at $`Q_s\tau |_{eq}\alpha _s^{13/5}`$ thermalization of these soft gluons with temperature $`T_{eq}\alpha _s^{2/5}Q_s`$ hydrodynamic expansion $``$ hadronization (Fig. 1).

In more detail the parametrically estimated time scales are determined as follows .

Early times $`1Q_s\tau \alpha _s^{3/2}`$:

At the earliest time, $`\tau Q_s^1`$, gluons, i.e. hard gluons, have typical large transverse momentum of order $`Q_s`$ and occupation number of order $`1/\alpha _s`$. Later on gluons with smaller momenta, but still larger than $`\mathrm{\Lambda }_{QCD}`$ will be produced (nevertheless, they are denoted as soft gluons). The density of hard gluons $`n_{hard}`$ decreases with time due to the one-dimensional expansion. Gluons interact by elastic scatterings at small angle, with exchange momentum $`Q_s`$. The typical occupation number is large until $`Q_s\tau \alpha _s^{3/2}`$, when it becomes of $`O(1)`$.

This regime is the transition region from the non-linear classical gluon field to the one where the transport description by Boltzmann equations should become applicable.

Times $`\alpha _s^{3/2}Q_s\tau \alpha _s^{5/2}`$:

Inelastic scatterings produce (soft) gluons $`n_{soft}`$ with characteristic momentum estimated to be of order $`\alpha _s^{1/2}Q_s>\mathrm{\Lambda }_{QCD}`$ (Fig. 1), namely via $`hard+hardhard+hard+soft`$. The number of these soft gluons becomes comparable to that of hard ones at $`Q_s\tau \alpha _s^{5/2}`$, namely $`n_{soft}n_{hard}`$.

Times $`Q_s\tau \alpha _s^{5/2}`$:

After $`Q_s\tau \alpha _s^{5/2}`$ most gluons are soft, $`n_{soft}n_{hard}`$; they achieve thermal equilibration amongst themselves. Although the whole system is still not in thermal equilibrium, the soft part is characterized by the temperature $`T`$, $`n_{soft}T^3`$. The few hard gluons collide with the soft ones of the thermal bath and constantly loose energy to the latter. A hard gluon emits one with a softer energy, which splits into gluons with comparable momenta. The products of this branching quickly cascade further, giving all their energy to the thermal bath. This increases the temperature in the thermal bath found to be,

$$T=c_T\alpha _s^3Q_s^2\tau ,$$
(12)
i.e. it increases linearly with time, even when the system is expanding, due to the hard gluons which serve as an energy source. $`c_T`$ is a numerical constant to be discussed later. One can verify a posteriori that soft gluons are indeed equilibrated due to many interactions: the size of the expanding system of $`O(\tau )`$ is indeed much larger than the mean free path of the soft gluons, i.e.

$$\tau /\lambda _{soft}\tau n_{soft}(\tau )\sigma c_T\alpha _s^5(Q_s\tau )^21,$$
(13)
where the cross section is estimated as $`\sigma \alpha _s^2/T^2`$.

In the following, in order to provide predictions for particle multiplicities in this scenario we go as much as possible beyond the parametric estimates, ending up with consistency arguments which specify the allowed range of parameters.

First we start by considering the constant $`c_T`$ which may be written in terms of the parameter $`c`$ as ,

$$c_T\frac{15}{8\pi ^5}cN_c^30.16c.$$
(14)
The linear growth of $`T`$ as shown in Fig. 1 terminates, when the hard gluons loose all of their energy, i.e. when

$$\tau =\tau _{eq}=c_{eq}\alpha _s^{13/5}Q_s^1,$$
(15)
where the parameter $`c_{eq}`$ is unknown for the moment, although in principle it could be calculated in the “bottom-up” framework. The temperature achieves a maximal value, i.e. the equilibration temperature, which is expressed by

$$T_{max}=T_{eq}=0.16cc_{eq}\alpha _s^{2/5}(Q_s^2)Q_s.$$
(16)
Subsequently the temperature decreases as $`\tau ^{1/3}`$ , such that $`\tau n_{soft}(\tau )=const`$.

We now derive the charged hadron multiplicity from the number density of the equilibrated soft gluons,

$$n_{soft}(\tau _{eq})=2(N_c^21)\frac{\zeta (3)}{\pi ^2}T_{eq}^3,$$
(17)
which is much larger than $`n_{hard}`$ in (1) used by KLN. One finds for the ratio $`R`$ defined in (11),

$`R=`$ $`8\zeta (3)N_c{\displaystyle \frac{c_T^3}{c}}\alpha _s^{10}(Q_s^2)(Q_s\tau )^4|_{\tau _{eq}}`$ (18)
$``$ $`0.13c^2c_{eq}^4\alpha _s^{2/5}(Q_s^2).`$ (19)
For the charged hadron multiplicity the result for $`\sqrt{s}=130GeV`$ and the most central collisions, is

$`{\displaystyle \frac{2}{N_{part}}}{\displaystyle \frac{dN_{ch}}{d\eta }}`$ $`{\displaystyle \frac{Rc}{3}}\mathrm{ln}{\displaystyle \frac{\overline{Q_s^2}}{\mathrm{\Lambda }_{QCD}^2}}`$ (21)
$`0.04c^3c_{eq}^4\left[\mathrm{ln}{\displaystyle \frac{\overline{Q_s^2}}{\mathrm{\Lambda }_{QCD}^2}}\right]^{7/5},`$
which replaces (7).

In order to compare the predicted charged hadron multiplicity, (7) and (21), with RHIC data we use as a reference the result by the PHOBOS Collaboration , namely $`3.24\pm 0.1(stat)\pm 0.25(syst)`$, which is in good agreement within errors with the experimental measurements of the other collaborations at RHIC. This is best seen by comparing the values for $`dN_{ch}/d\eta `$ at midrapidity at $`\sqrt{s}=130GeV`$: $`dN_{ch}/d\eta =555\pm 12(stat)\pm 35(syst)`$ , $`609\pm 1(stat)\pm 37(syst)`$ , $`549\pm 1(stat)\pm 35(syst)`$ , $`567\pm 1(stat)\pm 38(syst)`$ , respectively.

As reference value we take

$$\frac{2}{N_{part}}\frac{dN_{ch}}{d\eta }=3.24.$$
(22)
First we note that in the KLN approach agreement with experimental data can only be achieved with the value $`c3`$, which is different from the current numerical estimate quoted before, $`c0.5`$.

Turning to the “bottom-up” approach, the experimental value (22) meets the theoretical expectation (21) for $`Rc3`$, or equivalently

$$c_{eq}\frac{2.0}{c^{3/4}}.$$
(23)
This relation is to be confronted with the consistency requirement for the ”bottom-up” scenario , that the ratio $`R`$ (19) be larger than $`2`$, implying

$$c^2c_{eq}^410,i.e.c_{eq}\frac{1.8}{\sqrt{c}}.$$
(24)
It is not too difficult to see that (23, 24) constrain the two parameters:

$$c1.5andc_{eq}1.5.$$
(25)
This in turn allows us to infer the actual properties of the medium, especially to answer tentatively the question of the formation of the equilibrated plasma. We do this via discussing the temperature $`T_{eq}`$ and the equilibration time $`\tau _{eq}`$. In order that the quark gluon plasma is produced, $`T_{eq}`$ should be bigger than the phase transition temperature $`T_{deconf}`$, which is of order $`T_{deconf}=173\pm 8MeV`$ and $`154\pm 8MeV`$ for 2 and 3 flavour QCD, respectively , i.e.

$$T_{eq}T_{deconf}.$$
(26)
This constrains

$$cc_{eq}>1.3,orc>0.2.$$
(27)
Finally, we may correlatively discuss $`\tau _{eq}`$. Under the condition (25) one finds (for central collisions)

$$\tau _{eq}2.6fm1/2R_{Au},$$
(28)
which is much bigger than the current estimate of $`0.7fm`$  This short equilibration time is required in order to describe the hadron spectra for different particles as measured at RHIC by hydrodynamic calculations, assuming the existence of the quark gluon plasma  .. Within the uncertainties inherent to these estimates the formation of the equilibrated plasma may indeed be realized.

One may consider the energy dependence at RHIC energies of the charged multiplicities. Following  it is controlled by the energy dependence of the saturation scale, i.e. $`Q_s^2(s)/Q_s^2(s_0)=(s/s_0)^{\lambda /2}`$ with $`\lambda =0.25`$ . As discussed we choose as the reference value $`\overline{Q_s^2}(s_0)=1GeV^2`$ at $`\sqrt{s_0}=130GeV`$. For the most central collisions the expression for the energy dependence of the pseudorapidity density of charged particles at midrapidity reads in the ”bottom-up” scenario,

$$\frac{2}{N_{part}}\frac{dN_{ch}}{d\eta }0.64\left(\frac{\sqrt{s}}{\sqrt{s}_0}\right)^\lambda \left[\mathrm{ln}\frac{\overline{Q_s^2}(s)}{\mathrm{\Lambda }_{QCD}^2}\right]^{7/5}.$$
(29)
Based on (29) the multiplicity at $`\sqrt{s}=200GeV`$ is obtained:

$$\frac{2}{N_{part}}\frac{dN_{ch}}{d\eta }=3.84,$$
(30)
compared to the value given by the PHOBOS Collaboration ,

$$\frac{2}{N_{part}}\frac{dN_{ch}}{d\eta }=3.78\pm 0.25(syst).$$
(31)
Finally we consider the centrality dependence, i.e. the dependence of $`dN_{ch}/d\eta `$ as a function of $`N_{part}`$. In Fig. 2 the comparison of the ”bottom-up” expectation as derived from (29) with data from the PHOBOS Collaboration  at $`\sqrt{s}=130GeV`$ and at $`\sqrt{s}=200GeV`$ at RHIC is shown. For this comparison $`\overline{Q_s^2}(\stackrel{}{b})`$ is calculated from (8) as a function of $`N_{part}`$, using $`\rho _{part}(\stackrel{}{b})`$ given in Table 2 of , together with the scaling relation for the energy dependence of $`Q_s`$ quoted above . One has to note, that for $`N_{part}<100`$ the values of $`Q_s^2`$ are becoming smaller than $`0.6GeV^2`$. The structure of the shape of $`dN_{ch}/d\eta `$ as a function of $`N_{part}`$ seen in the data (Fig. 2) could be attributed, cf. with (6), to details of $`xG(x,Q_s^2)`$ at small $`x`$ as a function of $`Q_s^2`$.

In summary, the description provided by the “bottom-up” scenario is in agreement with RHIC data, provided the parameters $`c`$ and $`c_{eq}`$, which are not determined in the picture, lie in a given, limited range. In particular for $`c1`$ and $`c_{eq}2`$ the picture looks both reasonable and attractive. For these values the results (for the most central collisions) for $`R,T_{eq}`$ and $`\tau _{eq}`$ are $`3,230MeV`$ and $`3.6fm`$, respectively, for $`Q_s=1GeV`$; for $`Q_s=0.8GeV`$, i.e. for $`K=1`$, the values are $`R=3,T_{eq}=210MeV`$ and $`\tau _{eq}=3.2fm`$. However, the picture probably does not make much sense for a small value of $`c1/2`$, which is currently favoured by numerical calculations of classical field equations, because the ratio $`R`$ turns out to become $`R6`$, which is too much inelasticity when $`Q_s=1GeV`$. This discrepancy deserves further examination.

On the other hand, in order to accomodate the present RHIC data, the KLN description requires $`c3`$. For the moment, due to a number of ambiguities, it is not yet clear whether this large value is compatible with the various constraints of the saturation picture.

The “bottom-up” scenario is likely to provide a more convincing agreement with data. However, this analysis will have to be supplemented by further ingredients, e.g. the calculation of $`c_{eq}`$, a more precise estimate of $`c`$, before one can finally claim a significant agreement with data.

Discussions with D. Kharzeev, M. Nardi, K. Redlich and R. Venugopalan are kindly acknowledged. RB acknowledges support, in part, by DFG, project FOR 339/2-1. The work of AHM is supported, in part, by a DOE Grant. The work of DTS is supported, in part, by a DOE Grant No. DOE-ER-41132 and by the Alfred P. Sloan Foundation. This research was supported in part by the National Science Foundation under Grant No. PHY99-07949; RB and DS thank ITP, UCSB, for kind hospitality, and for the facilities during the completion of this work.


text:
# 𝐽𝐻𝐾' Imaging Photometry of Seyfert 1 AGNs and Quasars I: Multi-Aperture Photometry

## 1 Introduction

Near-infrared (NIR) observations are useful to study the dust distribution in Active Galactic Nuclei (AGNs). The unified model assumes a dust torus around the central engine in which the torus viewed at different angles from the line of sight explains the difference between Seyfert 1 AGNs and Seyfert 2 AGNs (Antonucci & Miller 1985).

NIR observations are also useful to derive general features in the spectral energy distributions (SEDs) of AGNs, and to examin how much these SEDs are contaminated by the light from host galaxies. Previous authors (Sanders et al. 1989; Neugebuer et al. 1979) presented the SEDs of 109 AGNs which show that the general features of a 1$`\mu m`$ minimum and a 3$`\mu m`$ bump were present in the SEDs of many AGNs. Kobayashi et al. (1993) obtained the SEDs of 14 quasars by 16 channel prism spectrophotometry between $`0.95\mu m`$ and $`2.5\mu m`$, and found that all the SEDs are characterized as having a black body SED with a typical temperature of 1500K correspondng to that of dust sublimation, seperated by a power-law component with a variety of power indicies. Elvis et al. (1994) presented the SEDs of 47 AGNs (29 radio-quiet and 18 radio-loud AGNs) over a wide range of wavelength from the X-ray to the radio region.

Hunt et al. (1997) obtained $`JHK`$ imags of 26 nearby AGNs. Alonso-Herrero, Ward & Kotilainen (1996) obtained the $`JHKL`$ imags of 13 Seyfert 2 AGNs and decomposed their SEDs into the stellar and non-stellar components. Because NIR imaging photometry has been limited to nearby AGNs, the sample size has also been limited. Motivated by the need to expand the sample size at least by an order of magnitude, we undertook a program of imaging a few hundred AGNs in the NIR region, carried out a statistical study AGN SEDs with the largest sample ever made. Many of them were observed on two different nights seperated by a year or more, for the purpose of detecting the NIR variability of the AGNs in our sample.

In a series of three papers, we give our results from three-years of observations. In this Paper I, we present NIR magnitudes of more than 300 AGNs derived by multi-color, multi-aperture imaging photometory. Analysis and discussion of AGN variability will be presented in the forthcoming Papers II and III. We are currently conducting a project called MAGNUM (an acronym of Multicolor Active Galactic NUclei Monitoring, Kobayashi et al. 1998a, 1998b) which monitors AGN in the 11 passbands of $`UBVRIZJHKL`$, and aims to determine the distances to many AGNs by measurements of the delay time between optical and NIR variabilities. Therefore, this paper is not only an individual study, but also a preliminary study to select target AGNs for the MAGNUM Project.

## 2 Sample Selection and Observations

### 2.1 sample selection

In the beginning of this study, all AGNs were selected from the 6th version of the Quasars and Active Galactic Nuclei catalog (hereafter referred to the VV catalog; Veron-Cetty and Veron 1993). Additional AGNs were selected from new versions that were released during this sudy, the 7th version (Veron-Cetty and Veron 1996) and the 8th version (Veron-Cetty and Veron 1998). We summarize the selection criteria such as (1) coordinates, (2) AGN types, (3) absolute $`B`$-magnitudes, and (4) redshifts:

(1)Declinations were selected from $`\delta =10`$ degrees to $`+50`$ degrees, allowing for a wide coverage of right ascension, $`\alpha `$. This is necessary in order for the MAGNUM Project to observe many AGNs under good conditions during the entire year from Haleakala, on the Hawaiian Island of Maui at a latitude of $`+20`$ degrees where the MAGNUM Observatory is situated.

(2)Seyfert 1 AGNs and quasars were selected, excluding Seyfert 2 AGNs and BL Lac objects. This is necessary in order for the MAGNUM Project to observe thermal radiation from the innermost region of the dust torus which surrounds the central engine of AGNs. We excluded Seyfert 2 AGNs, because they are aligned with the dust torus edge-on so that the inner region is obscured. We also excluded BL Lac objects, because their SEDs are known to exhibit only weak thermal radiation from hot dust.

(3)Absolute $`B`$-magnitudes were selected to span from $`M_B=29`$ to $`18`$, enabling a statistical study of the $`M_B`$-dependence of various other spectral features of AGNs. This is necessary in order to discuss a Malmquist-type bias which affects an interpretation of any statistical study from bright, distant AGNs in a sample.

(4)Redshifts were selected to span from $`z=0`$ to 1, enabling a statistical study of $`z`$-dependence of various other spectral features of AGNs. By considering that the maximum wavelength covered by the MAGNUM camera is the $`L`$ band, the redshifts were limited to below unity. Otherwise the thermal radiation peaked at $`2\mathrm{\mu m}`$ corresponding to the 1500K temperature of dust sublimation shifts to much longer wavelengths, beyond the $`L`$-band filter.

Table 1 tabulates the basic quantities of 331 AGNs selected in this study. For the purpose of illustration, Fig. 1 shows the distribution of $`\alpha `$ and $`\delta `$ for all AGNs in the sample, and Fig. 2 shows their distribution of $`M_B`$ and $`z`$.

### 2.2 observation

All AGNs and quasars in the sample were observed with the 1.3m infrared telescope at the Institute of Space and Astronautical Science (ISAS), Japan. We used the PICNIC camera developed at the National Astronomical Observatory, Japan (NAOJ) for multi-color broard-band imaging (Kobayashi et al. 1994). The PICNIC camera has a NICMOS3 array of $`256\times 256`$ pixels, corresponding to a field of view of $`4.57^{}\times 4.57^{}`$ and a pixel scale of $`1.07^{\prime \prime }`$ pixel<sup>-1</sup>. Imaging photometry was done in the NIR broad bands with the $`J`$, $`H`$, and $`K^{}`$ filters. In order to decrease the thermal sky background, we used $`K^{}`$ filter which has the same transmission curve as the 2MASS $`K_s`$ filter (McLeod et al. 1995).

Our observational runs consist of three periods (January 1996$``$April 1996, November 1996$``$February 1997, and December 1997$``$April 1998). More than 300 AGNs were observed in the first and second periods. In the third period, however, more than 200 AGNs that had been observed in previous two periods were again observed in order to determine their variability. Analysis and discussion of the variability of AGNs will be presented elsewhere (Papers II and III).

The AGNs and quasars were imaged in the $`J`$, $`H`$, and $`K^{}`$ bands by stepping the telescope in a raster pattern. The tipical exposure time of each single frame was 35 sec ($`J`$), 17 sec ($`H`$), and 8 sec ($`K^{}`$), while a shorter exposure time was used if saturation might occur because of either a high thermal background, or a large flux from an object in the frame. The minimum number of frames for one object was 4 frames with $`2\times 2`$ positions ($`J`$), 4 frames with $`2\times 2`$ positions ($`H`$), and 9 frames with $`3\times 3`$ positions ($`K^{}`$). More frames were obtained for fainter objects. The maximum number of frames for one object was 50 frames ($`J`$), 200 frames ($`H`$), and 200 frames ($`K^{}`$), with $`5\times 5`$ positions in common.

Two photometric standard stars with different elevations were observed three times in one night, that is, at the beginning of observations, midnight and at the end of observations. These standard stars were imaged with a $`3\times 3`$ raster, and two frames were obtained at each position. In this way a total of 18 frames were obtained for one star, and the acquisition of all the frames in the $`J`$, $`H`$, and $`K^{}`$ bands took about 15 minutes per star. If such schedule was impossible because of bad weather or other reasons, the standard stars were observed before and after the observations of the AGNs.

Each night after the AGNs and standard stars were observed, dome flat images were obtained by turning a calibration lamp on and off in front of white board, then the dark current was measured with a cold opaque shutter blocking all external radiation.

## 3 Reduction

### 3.1 image reduction

All the frames of AGNs and standard stars were obtained with short integration times and with dithering, which produced a large number of the frames in the end. A short integration time was used to avoid saturation during measurements with high sky background. The dithering was used to minimize the effect of bad pixels and to make a sky flat of high quality.

The software system specialized to analyze the data obtained by the PICNIC camera (hereafter PICRED) was used for the reduction. PICRED is a semi-automated system, requiring manual operations and human decisions in order to deal with various types of data, from star formong regins to quasars. For our case of reducing an enormous amount of the data by repeating much the same procedure, the software system was made fully automated.

### 3.2 photometric calibration

Each time a pair of standard stars was observed, 18 frames were taken first for one star at high elevation for each passband, in order, from $`J`$ to $`H`$ and then to $`K^{}`$, and then similarly for another star at low elevation. This procedure was repeated three times in one night, mostly for different pairs. Abnormal data, deviating remarkably from others, would possibly occur due to major three factors, such as obscuration by thin clouds that were not detected during the observation, the effect of bad pixels, or misidentification of the target.

In the beginning of reduction process, bad frames, if any, were excluded, and only the remaining frames were used to determine the instrumental magnitudes from which the median $`m(\mathrm{inst})`$ and dispersion $`\sigma _d`$ were obtained for each star in each passband. The median $`m(\mathrm{inst})`$ was then transformed to the already calibrated(known) magnitude of each star, by taking a linear fit in a plot of $`\mathrm{\Delta }m=m(\mathrm{inst})m(\mathrm{calib})`$ against airmass. The error $`\sigma _a`$ in transformation was also estimated. In this way, we determined the aperture $`J`$, $`H`$, and $`K^{}`$ magnitudes of AGNs using four different apertures of 7, 10, 12, and 15 pixels in radius. The total error in magnitude is given by $`\sigma _m^2=\sigma _d^2+\sigma _a^2`$. These aperture magnitudes and errors are tabulated in Tables 2a ($`J`$ band), 2b ($`H`$ band) and 2c($`K^{}`$ band).

Figure 3 shows the frequency distribution of $`\sigma _a`$ for the aperture of 15 pixels in the $`J`$, $`H`$, and $`K^{}`$ bands. For this large aperture of 15 pixels, the error $`\sigma _a`$ would originate from the variation of atmospheric transmissivity rather than the variation in the seeing during the night. The passband of longer wavelength has the smaller $`\sigma _a`$ distribution. This feature indicates that detected photon counts in the $`J`$ band as compared to the $`H`$ and $`K^{}`$ bands is more sensitive to the variation of atmosphere transmissivity.

Figure 4 shows the frequency distribution of $`\sigma _d`$ for the aperture of 15 pixels in the $`J`$, $`H`$, and $`K^{}`$ bands. In each panel, the distributions shown are based on our observations in the first and third periods (solid line; January 1996$``$April 1996, December 1997$``$April 1998) and in the second period (dashed line; November 1996$``$February 1997). These distributions are similar to each other, except that the peak for the first and third periods occurs at larger $`\sigma _d`$ than that for the second period.

Irrespective of the observational period, however, the dispersion $`\sigma _d`$ is larger than that expected from the high S/N ratio ($`0.01`$ mag), or from small changes of atmospheric transmissivity and airmass during the short exposure time of about 15 minutes. Therefore, the dispersion $`\sigma _d`$ may originate from systematic errors in flat fielding.

Our observations in the second period were made by avoiding use of the fourth quadrant of the detector because it was out of order. The dithering shifts in the second period were at most about 30 arcsec, only half of those for normal observations in the first and third periods. The resulting difference in dithering patterns gives a measure of the flat fielding error over some tens of arcsec that is comprable to $`\sigma _d`$.

The average and median of $`\sigma _d`$ for the first and third periods are 0.033 mag and 0.027 mag, respectively. Corrresponding values for the second period are smaller. This result is almost independent of passband, in contrast to that for $`\sigma _a`$. Thus, $`\sigma _a`$ is larger than $`\sigma _d`$ in the $`J`$ band, while the converse is true in the $`K^{}`$ band.

## 4 Results

The ratio of radio flux $`f_\nu (6\mathrm{c}\mathrm{m})`$ relative to optical $`V`$-band flux $`f_\nu (V)`$ is used as a measure of radio strength of the AGNs in our sample. Figure 5 shows the distribution of this ratio based on the data taken from the VV catalog. The values of this ratio range over several orders, but are localized around 1 and 1000. Here, in this paper, the AGNs with $`f_\nu (6\mathrm{c}\mathrm{m})/f_\nu (V)<10`$ and no radio detection are classified as radio quiet, and those with $`f_\nu (6\mathrm{c}\mathrm{m})/f_\nu (V)>100`$ as radio loud.

It is known that there are two typical locations in the two-color $`JH`$ to $`HK^{}`$ diagram, such as ($`JH`$, $`HK^{}`$)=(0.8, 1) for the “QSOs” (Hyland & Allen 1982), and (0.7, 0.3) for “galaxies” (Willer et al. 1984). AGNs in our sample are found to be distributed in a wide region from “QSOs” to “galaxies”. In fact, brighter AGNs tend to populate the diagram near the “QSOs”, while fainter AGNs or Syfert $`1.82`$ AGNs near the “galaxies”. No such localization, however, occurs if the sample is divided into the radio-quiet and radio-loud AGNs.

All these features are more clearly seen in Fig. 6, where averages of $`JH`$ and $`HK^{}`$ colors and their errors for low-$`z`$ AGNs with $`0<z<0.3`$ are plotted as a function of $`M_B`$, Seyfert type, and $`f_\nu (6\mathrm{c}\mathrm{m})/f_\nu (V)`$. Shown are the results for the four different apertures of 7, 10, 12, and 15 pixels. Here, averages are taken with no weights, and errorbars are the standard deviations of the colors for an aperture of 7 pixels.

The $`HK^{}`$ color becomes monotonically bluer from $`M_B=27`$ to $`21`$, irrespective of aperture. AGNs with $`M_B<25`$ have $`HK^{}1`$ and are QSO-like, while those with $`M_B>22`$ have $`HK^{}0.5`$ and are galaxy-like. On the other hand, the $`JH`$ color stays at about 0.8, and its $`M_B`$-dependence is much weaker than that for $`HK^{}`$.

Furthermore, for AGNs with $`M_B>24`$, the $`HK^{}`$ color determined with larger aperture is more galaxy-like. This feature originats from the color gradient in the central region where the AGN dominates to the outer region where the host galaxy dominates. It should be noticed that the blueward color shift with the use of larger aperture is more significant for fainter AGNs with $`M_B>24`$ and for the $`HK^{}`$ color rather than for $`JH`$. This feature can be explained naturally by considering that the host galaxy becomes more visible within larger aperture and its SED is enhanced over the AGN at NIR wavelengths in rest frame.

The monotonical trend in the average colors is also seen by changing Seyfert type from 1 to 2. Although Seyfert $`11.5`$ AGNs are in between “QSOs” and “galaxies”, Seyfert $`1.82`$ AGNs are galaxy-like. This partly reflects the tendency that Seyfert $`1.82`$ AGNs are, on the whole, faint in $`M_B`$, because of obscuration of the central AGN component by the dust torus. On the contrary, no such trend is seen with changing radio strength. The $`HK^{}`$ color for radio-quiet AGNs is much the same as that for radio-loud AGNs and is in between “QSOs” and “galaxies”. This indicates that the radio strength has no significant correlation with $`M_B`$ in its range considered here.

In order to see the $`z`$-dependence, we show the $`JH`$ and $`HK^{}`$ colors of intermediate-$`z`$ AGNs with $`0.3<z<0.6`$ by dashed lines in Fig. 6 only for the case of an aperature of 10 pixels. It is apparent that the QSO-like colors of brighter AGNs with $`M_B<25`$ shift bluewards for larger $`z`$, while galaxies-like colors of Seyfert 1.8 shift redwards. This opposite color trend is consistent with the opposite of $`k`$-corrections between AGNs and galaxies, and calculations based on the typical SEDs have confirmed that the values of the $`k`$-corrections for AGNs and galaxies indeed agree with their respective color shifts as seen in Fig. 6. Thus, a criterion of $`M_B<27`$ is regarded as a discriminator of QSOs, and a criterion of Seyfert type later than 1.8 as a discriminator of AGNs dominated the galaxy SED component. Otherwise, the intermediate colors, as a result of being contributed equally from QSO and galaxy components, do not show any significant $`z`$-dependence, because the $`k`$-corrections of the different conponents cancel any z-dependence.

We are grateful to H. Okuda, M. Narita, and other staff of the infrared astronomy group of the Institute of Space and Astronautical Science (ISAS) for their support in using their 1.3m telescope. We thank the staff of the Advanced Technology Center of the National Astronomical Observatory of Japan (NAOJ) for their new caoting on the mirror of the 1.3m telescope at the ISAS. Gratitude is also extended to the Computer Data Analysis Center of the NAOJ. This work has made use of the NASA/IPAC Extra Galactic Database (NED), and has been supported partly by the Grand-in-Aid (07CE2002, 10304014) of the Ministry of Education, Science, Culture, and Sports of Japan and by the Torey Science Foundation.


text:
# The Scintillation Velocity of the Relativistic Binary Pulsar PSR J1141–6545

## 1 Introduction

The scintillation observed in the radio wavelength emission from pulsars is considered to be a result of diffraction of the incident wave front around small scale irregularities in the interstellar medium (ISM) (Rickett 1969). The characteristic timescale for these scintillation events can be hours, minutes, or even seconds and is highly dependent upon a number of factors: the distance to the source; the relative velocities of the source, observer and scattering material; the observing frequency; and the structure and distribution of the ISM (Cordes, Pidwerbetsky, & Lovelace 1986; Cordes & Rickett 1998). It has long been considered possible that observations of binary pulsars would reveal a modulation in scintillation velocity due to their orbital motion relative to the scattering medium (Lyne & Smith 1982). Such orbital modulation in the scintillation velocity of a binary pulsar has been significantly detected only in the circular binary PSR B0655+64 (Lyne 1984). Dewey et al. (1988) have reported a weak orbital modulation of scintillation velocity in observations of PSR B1855+09 and PSR B1913+16, but were unable to constrain any system parameters.

The ideal pulsar for such measurements would scintillate on timescales much smaller than the orbital period and have enough flux to obtain high signal-to-noise ratios within the scintillation timescale. In addition, the orbital period should be much less than the timescale for changes in the scintillation parameters due to macroscopic variations in the structure of the ISM. Finally, the pulsar should display large variations in its transverse orbital velocity.

The relativistic binary pulsar, PSR J1141–6545, was recently discovered in the Parkes multibeam pulsar survey (Kaspi et al. 2000). It has a spin period ($`P`$) of 394 ms, a very narrow ($`0.01P`$) pulse, a dispersion measure of 116 pc cm-3 and is in an eccentric, 4.7 hr orbit. It is thought to be a member of a new class of object that has a young neutron star in an eccentric orbit with a massive $``$1 M white dwarf companion. Kaspi et al. (2000) found no evidence for scintillation. We have searched specifically for short–timescale scintillation and found that the pulsar scintillates on timescales of minutes with a characteristic bandwidth of just $``$1 MHz at a central observing frequency of 1390 MHz. At the Parkes 64 m observatory it completes two orbits during its transit time, changes velocity by over 200 km s-1, and has a short orbital period. In every respect, this pulsar is an ideal target for scintillation work. Recent neutral Hydrogen observations (Ord, Bailes & van Straten 2002) place the binary at least as distant as the Carina–Sagittarius spiral arm, which is 3.7 kpc away in the direction of PSR J1141–6545.

Unlike their progenitors, the massive O and B stars, radio pulsars have very large space velocities of up to 1000 km s-1 (Lyne, Anderson, & Salter 1982; Bailes et al. 1989; Harrison, Lyne, & Anderson 1992). The origin of pulsar velocities is difficult to ascertain as individual pulsars give little clues regarding the responsible physical mechanism (Radhakrishnan & Shukre 1986; Dewey & Cordes 1987; Bailes 1989). Eccentric binary pulsars, on the other hand, are fossil records of the state of a binary before detonation of the pulsar progenitor. In the absence of a natal kick, a binary pulsar has a well-defined relationship between its eccentricity, the amount of mass ejected during the formation of the neutron star, and the runaway velocity (Radhakrishnan & Shukre 1986). However, eccentric binary pulsars are rare and it has not been possible to make a definitive statement about kicks from an analysis of their space velocities and orbital configurations (Cordes & Wasserman 1984; Hughes & Bailes 1999; Wex, Kalogera, & Kramer 2000). Kaspi et al. (1996) have shown that the massive binary PSR J0045–7319 has a precessing orbit consistent with the pulsar receiving an impulsive kick at birth. In the case of PSR J1141–6545, the space velocity, had the pulsar received no natal kick, should be small as the eccentricity is low. But evolutionary arguments concerning the minimum diameter of the progenitor (Tauris & Sennels 2000) suggest that its runaway velocity should be at least 150 km s-1.

In this paper, we present the detection of an orbital modulation of the scintillation velocity in observations of PSR J1141–6545. The effect is very significant, permitting the determination of a number of system parameters. Measurements have been obtained of the most probable space velocity and inclination of this system, as well as the longitude of periastron. As with the analysis of PSR B0655+64 by Lyne (1984), there is a clear degeneracy between two distinct solutions of equal significance. In this case, however, the degeneracy can be broken: one of the solutions indicates an unreasonably low neutron star mass. The structure of this paper is as follows: in section 2 we describe our observations and methodology for determining the scintillation parameters; section 3 describes the model and our fitted results; finally, in section 4 we discuss the results and their implications for future relativistic observables and the origin of pulsar velocities.

## 2 The Observations and Data Reduction

PSR J1141–6545 was observed on the 27th of January 2002 for approximately ten hours, using the Parkes 64 m radio telescope. The $`512\times 0.5`$ MHz filterbank centred at 1390 MHz was used with the centre element of the Parkes multi-beam receiver (Staveley-Smith et al. 1996). In each filterbank channel, the power in both linear polarisations was detected and summed before 1 bit sampling the total intensity every 250 $`\mu `$s. The results were written to magnetic tape for offline processing.

Average pulse profiles were produced by folding the raw data in each 500 kHz frequency channel modulo the topocentric pulse period. An integration time of 29 seconds provided sufficient resolution to resolve the scintillation structure throughout the orbit. For each profile, the mean level of the off-pulse region was subtracted from that of the on-pulse region to produce a dynamic spectrum. Figure 1 shows the clear orbital modulation of the dynamic spectrum. In each orbital period, there is an obvious apparent vertical “blurring” as the scintillation timescale varies.

In order to model the physical parameters of the binary system, we require estimates of the pulsar’s velocity as a function of orbital phase. Individual measurements of the transverse velocity were obtained at different epochs throughout the observation by determining the scintillation timescale and bandwidth. These were derived from the two-dimensional auto-correlation function (ACF).

Individual ACFs were calculated over a small region of the dataset that spanned approximately 12 minutes in the time domain and 24 MHz in the frequency domain. A two-dimensional fast Fourier transform (FFT) was calculated and multiplied by its complex conjugate to form the ACF. To improve the signal-to-noise ratio, a single ACF was formed at each timestep by combining all the ACFs produced across the bandpass. The scintillation timescale was defined to be the $`1/e`$ full width of a Gaussian fit to the central peak of the ACF in the direction of increasing time-lag. Following convention (Cordes 1986), the scintillation bandwidth was defined to be half the width at half-height of a Gaussian fit across the frequency lags.

By moving the ACF window through the dynamic spectrum with a timestep of one time lag, or 29 seconds, scintillation parameters were assigned a time corresponding to the central lag of the ACF. These values were integrated into 32 binary phase bins using the ephemeris presented by Kaspi et al. (2000). The value of the scintillation velocity was determined using the following relationship.

$$V_{\mathrm{ISS}}=2.53\times 10^4\frac{\left(D\mathrm{\Delta }\nu _\mathrm{d}\right)^{1/2}}{f\tau _\mathrm{d}}$$
(1)
Equation 1 is after Cordes and Rickett (1998). $`V_{\mathrm{ISS}}`$ is the scintillation velocity, $`D`$ is the earth-pulsar distance in kpc, $`\mathrm{\Delta }\nu _\mathrm{d}`$ is the scintillation bandwidth in MHz and $`\tau _\mathrm{d}`$ is the scintillation timescale in seconds. The constant multiplicative factor is the value determined by Cordes and Rickett (1998) for a uniform Kolmogorov scattering medium. The distance used was 3.7 kpc (Ord, Bailes & van Straten 2002). Any uncertainty is incorporated into a scaling parameter within the model.

This process produced the scintillation velocity as a function of mean orbital anomaly. The system is significantly eccentric ($`e0.17`$) and the model described in §3 constructs the velocity as a function of true, not mean, orbital anomaly. We therefore transform the observed scintillation velocity into a function of eccentric anomaly, $`\eta `$, by an iterative solution of Kepler’s equation, $`ϵ=\eta +e\mathrm{sin}\eta `$, where $`ϵ`$ is the mean anomaly. Eccentric anomaly is converted into true anomaly, $`\theta `$, via a simple trigonometric relationship:

$$\mathrm{tan}\frac{\theta }{2}=\sqrt{\frac{1+e}{1e}}\mathrm{tan}\frac{\eta }{2}.$$
(2)
## 3 The Model

To determine the most probable runaway velocity and inclination angle of this system, we constructed a model that calculated the transverse velocity as a function of five free parameters. These were: the components of the transverse velocity along and perpendicular to the line of nodes, $`v_{\mathrm{plane}}`$, and $`v_{\mathrm{per}}`$, respectively; the orbital inclination angle, $`i`$; a scaling factor, $`\kappa `$; and the longitude of periastron, $`\omega `$. The latter is accurately determined by pulse timing, and allows us to confirm the validity of our results.

The orbital velocity was calculated as a function of true orbital anomaly $`(\theta )`$ and broken into a radial component directed toward the focus of the orbital ellipse in the line of nodes ($`v_\mathrm{r}`$), the second component was perpendicular to $`v_\mathrm{r}`$ and in the direction of the orbital motion, $`v_\theta `$.

In this framework, the space velocity as a function of true orbital anomaly is given by:

$`v_\mathrm{r}`$ $`=`$ $`{\displaystyle \frac{2\pi xc}{\mathrm{sin}i(1e^2)^{1/2}P_\mathrm{b}}}e\mathrm{sin}\theta `$ (3)
$`v_\theta `$ $`=`$ $`{\displaystyle \frac{2\pi xc}{\mathrm{sin}i(1e^2)^{1/2}P_\mathrm{b}}}(1+e\mathrm{cos}\theta )`$ (4)
where the observable $`x=(a/c)\mathrm{sin}i`$ is the projected semi-major axis in seconds, $`a`$ is the semi-major axis of orbit, $`i`$ is the inclination of the system, $`e`$ is the orbital eccentricity and $`P_\mathrm{b}`$ the binary period.

After the construction of the binary model, the next stage was to transform this orbital velocity into a predicted scintillation velocity.

Assuming that the Earth had a near constant velocity throughout the observation and that the interstellar medium velocity was small compared to that of the pulsar, the model of the transverse velocity ($`V_{\mathrm{model}}`$) was described by the following equations:

$`v^{}`$ $`=`$ $`\left(\left(v_r\mathrm{cos}\varphi v_\theta \mathrm{sin}\varphi \right)+v_{\mathrm{plane}}\right),`$ (5)
$`v^{\prime \prime }`$ $`=`$ $`\left(\left(v_\theta \mathrm{cos}\varphi +v_r\mathrm{sin}\varphi \right)\mathrm{cos}i+v_{\mathrm{per}}\right),`$ (6)
$`V_{\mathrm{model}}`$ $`=`$ $`\kappa \sqrt{v^2+v^{\prime \prime 2}}.`$ (7)
The angle $`\varphi =\omega +\theta `$ is the true anomaly measured with respect to the line of nodes. The model parameter $`\kappa `$ was included to incorporate any errors introduced in the conversion of scintillation timescale to velocity into the model. This accounts for errors in the distance to the pulsar and the fact that the scattering medium is not uniform. The model was then evaluated at the same 32 values of true anomaly presented by the dataset. The best fit to the model was found by evaluating the chi-squared statistic throughout a wide range of each of the model parameters. This produced a chi-squared volume of 5 dimensions with a global minimum representing the most probable values of the 5 parameters.

A box-car smoothing operation was also applied to the model in order to compensate for the averaging introduced by the auto-correlation process. This had little effect on the values of the fitted parameters but did appreciably lower the chi-squared minimum.

The relative errors in each binary phase bin were calculated by first determining the scatter in each phase bin. The best fit to the data was found by location of the global minimum in the chi-squared space, and the error bars were normalised to ensure the reduced chi-squared statistic was unity. This enabled errors to be placed upon the estimates of the fitted parameters by examining the chi-squared distribution. It should be noted that this method assumes a priori that the model was a valid description of the data. A Monte Carlo analysis was then performed in order to ascertain the precision of the fitted parameters. This analysis demonstrated that the model described the observations particularly well and the range of allowable parameter values was very narrow.

Both the measured values of $`V_{\mathrm{ISS}}`$ and the best-fit model velocity, $`V_{\mathrm{model}}`$, are plotted as a function of true anomaly in Fig 2.

### 3.1 The effects of Refractive Interstellar Scintillation

The analysis described above would not have revealed any underlying systematic error introduced by either random or long-term variations in the scintillation timescale that are not associated with changes in the velocity of the pulsar. For example, the interstellar scintillation effects detailed here are a result of diffractive interstellar scintillation (DISS), but another regime of scintillation is refractive interstellar scintillation (RISS). RISS, which has a characteristic timescale generally much longer than that of DISS, was considered by Romani, Narayan and Blandford (1986) to be the explanation for the long term flux variability in pulsars observed by Sieber (1982). Although the effect of RISS is best characterised via multiple flux measurements, it is possible to estimate the RISS timescale, $`T_{\mathrm{RISS}}`$, using the measured DISS parameters:

$`T_{\mathrm{RISS}}`$ $`=150.6`$ $`\left({\displaystyle \frac{D(\mathrm{kpc})}{\mathrm{\Delta }\nu (\mathrm{kHz})V_7^2}}\right)^{1/2}\mathrm{days}.`$ (8)
Equation 8 is after Blandford and Narayan (1985). It assumes a uniform Kolmogorov medium and that the pulsar is travelling at $`100V_7`$ km s-1. In Equation 8 $`D`$ represents distance and $`\mathrm{\Delta }\nu `$ is the scintillation bandwidth. If the velocity of the source can be approximated by the measured scintillation velocity, then $`T_{\mathrm{RISS}}`$ is approximately 7 days. It is therefore unlikely that a modulation of this nature was present in the observation. Even if it was present, the modulation would be weak as RISS effects are considerably weaker than DISS events (Rickett 1990). For these reasons, the effect of RISS was discounted in this analysis. Future observations, if they span multiple days, may be prone to RISS effects that limit the precision of derived parameters.

## 4 Results

As discussed by Lyne (1984), analyses of this nature can produce two indistinguishable solutions. The degeneracy arises as the same apparent transverse velocity can be displayed by the pulsar in either a comparatively face-on binary system with a low runaway velocity, or a more edge-on system with a higher runaway velocity. An example of this can be seen in the chi-squared projection presented in Figure 3. The values of the fit parameters for both solutions after the Monte Carlo error analysis are presented in Table 1. The degeneracy between solutions can be broken by a reasonable consideration of the mass function and other timing parameters. Kaspi et al. (2000) derive a value for the mass function, $`f(M)`$ = 0.176 M, and an estimation of the sum of the system component masses, $`M_c=`$ 2.3 M, from timing measurements. Using these measurements, together with the system inclination presented here, the component masses may be obtained. The more edge-on solution indicates a pulsar and companion mass of 1.29 $`\pm `$ 0.02 M and 1.01 $`\pm `$ 0.02 M respectively; whereas the more face-on solution, suggests a pulsar mass of 1.17 $`\pm `$ 0.02 M and a companion mass of 1.13 $`\pm `$ 0.02 M. The edge-on solution is more likely as the indicated pulsar mass is more consistent with the known neutron star mass distribution, $`1.35\pm 0.04`$ M (Thorsett & Chakrabarty 1999). The most probable transverse velocity is therefore found to be of modulus 115 $`\pm `$ 15 km s-1, and the most likely inclination angle is 76 $`\pm `$ 2.5 .

In order to test the validity of the experimental method employed and the accuracy of the binary model, we allowed the angle of periastron to vary as a free parameter in the fit. The value of $`\omega `$ is already precisely determined by timing measurements and its value at the epoch of these observations is expected to be approximately 55.13. The consistency of this with our measurement of 58 $`\pm `$3.5 is very encouraging and increases our confidence in the solution which is extremely statistically significant.

## 5 Discussion

PSR J1141–6545 has been shown to have a transverse velocity of modulus 115 $`\pm `$ 15 km s-1 and an inclination of 76 $`\pm `$ 2.5. The component masses inferred from these values indicate, as first presented by Kaspi et al. (2000), that this system is most likely a near edge-on neutron star–CO white dwarf binary, with the observed neutron star having been formed most recently.

Tauris and Sennels (2000) have predicted, from evolutionary arguments, that PSR J1141–6545 would display a systemic velocity in excess of 150 km s-1. This velocity comes from both the recoil of the binary due to the rapid ejection of the exploding star’s envelope, and an impulsive kick imparted to the neutron star to leave it in an orbital configuration resembling the current state of PSR J1141–6545. The “no-kick” solution, suggests a velocity nearer 40 km s-1. The measured transverse velocity is only consistent with that predicted by Tauris & Sennels (2000), if the undetected radial velocity exceeds 96 km s-1. The simulations undertaken by Tauris & Sennels (2000) suggested that 150 km s-1 was the lower limit, and that PSR J1141–6545 would probably have a velocity much greater than this. The low tranverse velocity indicates that, as with the other eccentric binary pulsars, the kick was modest (Hughes & Bailes 1999).

It is also interesting to note that Ord, Bailes & van Straten (2002) give a lower distance limit to this system as the tangent point distance, 3.7 kpc. At this distance, and with a transverse velocity of 115 $`\pm `$ 15 km s-1, the expected contribution to any measured orbital period derivative due to proper motion and Galactic kinematic effects is expected to be at the 1 percent level. Furthermore, measurements of the range and shape of Shapiro delay, advance of periastron, gravitational redshift parameter and gravitational period derivative will over-determine the system parameters. Together with the determination of the system inclination presented here, estimates of these parameters will allow precise and unique tests of the validity of General Relativity.

This result suggests that long-term monitoring of its scintillation properties throughout a year may ultimately provide the orientation of the proper motion vector on the sky. An independent measure of the orbital period, eccentricity, longitude of periastron and inclination angle of the binary should also be obtained.