text:
# The line-of-sight warp of the spiral galaxy ESO 123-G23

## 1 Introduction

The outer parts of the disks of spiral galaxies are frequently warped (Bosma bosma (1991), García-Ruiz et al. garcia-ruiz (2002)). Warps are best seen with HI observations, since the neutral hydrogen extends out to much larger radii than the optical disk. HI warps obey some general rules, one of them being that they usually start around R25 (Briggs briggs (1990)). The measured angle between the inner plane and the outermost observed HI ring spans a wide range of values, from a few degrees to a few tens in some exceptional cases. Of the various possible orientations that warps can have, one is considered as “unfavourable”, as the warp becomes less visible: this is when the largest deviation from the central plane occurs along the line-of-sight.

The origin of warps is still poorly understood. Different models have been proposed, for instance: the precession of a disk embedded in a flattened dark halo, misaligned with the disk itself (e.g. Toomre toomre (1983)), or the accretion of gas with angular momentum different from that of the disk (e.g. Ostriker & Binney ostriker (1989)), but none is completely satisfactory.

In an ongoing project to study the dark matter distribution of spiral galaxies a sample was culled from the list of Persic and Salucci (persic (1995)) selecting the best determined optical rotation curves; the selection criteria were: symmetry of the H$`\alpha `$ rotation curve, high angular extent, high HI flux and low luminosity (Gentile et al., in prep.) The aim was to combine these optical data with HI measurements tracing the gravitational potentials much further out. One of these galaxies, the edge-on spiral ESO 123-G23, turned out to look remarkable in its HI distribution. Indeed the HI emission extends far away (up to about 15 kpc) from either side of the plane. Such a pattern can be produced by different phenomena including a flaring of the outer disk, a thick layer of gas or a warp along the line-of-sight. We have therefore modelled the HI layer of ESO 123-G23 considering these possibilities. Our analysis shows that the correct explanation is a warping of the disk along the line-of-sight. The other possibilities are clearly ruled out.

## 2 Observations

ESO 123-G23 is an edge-on spiral galaxy located at a distance of 38 Mpc ($`H_0`$=75 km s-1 Mpc -1). Not much was known about this galaxy except that it exhibits a slight warp in the outer optical disk (Reshetnikov & Combes reshetnikov (1999)). We have observed ESO 123-G23 in the HI line for 2 $`\times `$ 12 hours using the Australia Telescope Compact Array (ATCA)111The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. and combining the data from two array configurations (750D and 1.5A), providing a beam of 30$`\stackrel{}{.}`$6 $`\times `$ 25$`\stackrel{}{.}`$0 and a velocity resolution (after Hanning smoothing) of 6.6 km s-1. In Fig. 1 (upper panel) we display an optical - HI overlay showing the peculiar morphology of this galaxy: while the optical image resembles an almost perfectly edge-on disk, the HI is much more extended perpendicularly to the major axis, displaying what looks like a very inclined central disk and diffuse emission at large (up to 15 kpc) projected distances from the optical disk; this diffuse emission accounts for about 40 % of the total HI mass. The left column of Fig. 2 shows some channel maps of the HI data cube. Every fourth channel is shown, so the displayed maps are spaced by about 26 km s-1. The central channel map, at a velocity of 2855.4 km s-1, is the one closer to the systemic velocity. The channel maps of ESO 123-G23 show the emission from a central edge-on disk, plus fainter emission from regions above and below the plane. Such an emission is more visible in the central channels, close to systemic velocity. This may indicate either a lower rotation velocity of the gas above the plane or a lower inclination angle (i.e. towards face-on) of the outer parts of the disk. A warp in position angle is visible, especially at v = 2723.5 and 2987.3 km s-1, and it is quite symmetric. ESO 123-G23 is reported to have no optical companion of size between Ropt/2 and 2 Ropt within a radius of 10 Ropt (Karachentsev karachentsev (1993)). An inspection of the HI data cube did not disclose any companion in HI either. The total size of the data cube is $`33\mathrm{}\times 33\mathrm{}`$ (dimensions of the primary beam) $`\times 442`$ km s-1. The position-velocity diagram along the major axis is shown in Fig. 3 (upper panel), where the east side of the galaxy is slightly more extended but no clear kinematical lopsidedness is visible. On the p-v diagram is overlapped (dots) the rotation curve of the galaxy (see Section 3).

We considered three phenomena that might mimic the observed HI structure and kinematics: a) a line-of-sight warp; b) a large flare in the outer part of the disk; and c) a two-component structure, thin + thick (slowly rotating) disk, similar to that proposed for NGC 891 (Swaters et al. swaters (1997)) and NGC 2403 (Fraternali et al. fraternali (2001)). These phenomena have in common the fact of displaying emission away from the thin disk and are hard to disentangle except via detailed 3-D modelling of the HI layer.

## 3 Models

As mentioned above, we modelled the HI data cube in three different ways. All these models consist of a set of concentric circular rings (Rogstad et al. rogstad (1974)), each described by a set of geometrical and physical parameters, chosen such as to provide the best agreement with the observed HI data cube. Some parameters were kept fixed for the three models.

The parameters that the models had in common were: the radial surface density distribution, derived using the Lucy method (Warmels warmels (1988)); the HI velocity dispersion, chosen to go from 12 km s-1 in the inner parts to 7 km s-1 in the outer parts (Kamphuis kamphuis (1993)), in order to have values as realistic as possible; a constant value of 10 km s-1 would not have affected much the results, while values different from 10 km s-1 by more than 5 km s-1 would not give a good account of the observed data cube; the shape of the vertical distribution of the HI layer, chosen to be Gaussian with different scale-heights for the various models (the choice of other distribution shapes does not significantly change the results); the position angle as a function of radius, taken to slightly change counter-clockwise after $`R_{25}`$ = 80$`\stackrel{}{.}`$7 ($``$15 kpc), and then to significantly (10) change clockwise further out (see third panel of Fig. 5). By trying models with different values of the position angle we estimate an accuracy of about 2 (errorbars in Fig. 5), beyond which the disagreement with the observations is flagrant.

The rotation curve was derived with the WAMET (WArped Modified Envelope Tracing) method (Gentile et al. gentile (2002), Vergani et al. vergani (2002)), a modification of the Envelope Tracing method (Sofue & Rubin sofue (2001)); a detailed description of this method will be given in a forthcoming paper (Gentile et al., in prep.). Only small corrections ($`<`$ 5 km s-1) were necessary to refine the rotation curve derived with our WAMET method. These corrections were made in order to give a better representation of the observed HI data cube. This method, however, like the methods that derive the rotation curve from the position-velocity diagram, provides the product $`v_{\mathrm{rot}}sin(i)`$. Therefore we have to account for the inclination, that can vary with radius within the same model. For the inner parts of the galaxy we used, in all the models, an inclination angle of 85, a value similar to optical estimates (e.g. the Lyon-Meudon Extragalactic Database, LEDA, gives an inclination of 82.3).

In the warp model the inclination angle was lowered by up to $``$ 30 beyond $`R_{25}`$ in order to fit the observations (see second panel in Fig.5). With the same procedure we used for the position angle, we estimated an error for the inclination angle of 5. The rotation curve was obtained with the WAMET method with the assumed values for the inclination angle. The FWHM of the Gaussian vertical distribution of the HI was chosen to be 300 pc.

The flare model differs from the warp model in having a constant inclination of 85 and a rotation curve obtained with such a constant value. The FWHM of the HI vertical distribution is 300 pc inside $`R_{25}`$ increasing to a value of 25 kpc at $`R100`$″. Such an unrealistic value is necessary to reproduce the HI distribution (lower contours) visible in the total HI map (Fig. 1, upper panel).

Finally the two-component model is the sum of a thin and a thick disk in analogy to the previously cited models proposed for NGC 891 and NGC 2403. In this model the thin disk has a rotation velocity taken from the WAMET method with a constant inclination and a thickness of 300 pc (FWHM). For the thick disk (40% of the total HI mass) the rotation velocity was lowered by 30 km s-1 with respect to the thin disk and the thickness of the HI layer was taken to be 20 kpc (FWHM). Similarly to the flare model, this last unrealistic value is necessary to reproduce the total HI distribution. For both the flare and the two-component model a different choice of the inclination does not improve the agreement with the observed data cube.

## 4 Discussion and conclusions

In Fig. 2 we show some channel maps of the data cube compared with the channel maps of the line-of-sight warp model, the flare model, and the two-component model (from left to right). This figure illustrates why the flare model and the two-component model have to be rejected: apart from their somewhat unphysical choice of parameters (although necessary to produce the vertical HI distribution), they cannot reproduce the observed channel maps. For instance, the flare model cannot reproduce at the same time the large north-south extension in the central channels and the relatively small north-south extension visible in the extreme channels (v = 2723.5 and 2987.3 km s-1). The two-component model is not able to reproduce the ’elbow’ present at v = 2749.9 and 2802.6 km s-1, which is also indicated on the receding side, at v = 2908.2 and 2960.9 km s-1.

The data are indeed fitted best by the model of a warp that shows maximum deviation from the plane close to the line-of-sight; Fig. 2 shows that such a model reproduces in detail all the key features of the observed HI data cube. See the upper panel of Fig. 5 for the adopted rotation curve for this model; the errors are the difference between the approaching and the receding side, considering also a minimum error of $`2/sin(i)`$ km s-1. In the upper panel of Fig. 3 we show the derived rotation curve of ESO 123-G23 (separately for the two sides) projected onto the observed position-velocity diagram, and in the lower panel the adopted rotation curve (the average between the two sides) projected onto the position-velocity diagram of the line-of-sight warp model. It is remarkable that the de-warped rotation curve of this galaxy is in good agreement with the Universal Rotation Curve of Persic et al (persic2 (1996)) while the raw curve, instead, had a very unusual decrease of 30 km s-1 from 7 kpc to 15 kpc, never found in galaxies with the maximum circular velocity of ESO 123-G23. Fig. 4 shows the HI rings of the line-of-sight warp model as seen by the observer and Fig. 1 (lower panel) shows the total HI map obtained with this model. The overall HI distribution is well reproduced. We have also tried combinations of the various models, but the emission at projected distances of up to 15 kpc from the central disk can be explained only with the warp model, and combining the warp model with realistic values of an outer flare and a thick disk gives irrelevant changes to the model.

In order to describe the geometry of the warp, as previously done by Schwarz (schwarz (1985)), we define the angle $`\theta `$ between the plane of the inner disk and the plane of the warped ring, and the angle $`\beta `$ between the line of nodes (at the intersection of the plane of the ring and the plane of the inner disk) and the line of intersection between the plane of the sky and the plane of the inner disk. The uniqueness of this warp now shows up in Fig. 5, bottom panel: $`\beta `$ is always relatively close to zero, meaning that the maximum tilt of rings is almost exactly along the line-of-sight. We point out that a negative value of $`\beta `$ means a position angle varying clockwise and that in unwarped sections of a galaxy $`\beta `$ is of course not defined. The fourth panel from top of Fig. 5 shows the monotonically increasing warp angle $`\theta `$ as a function of radius.

There is, however, an ambiguity due to our assumption of a central inclination of 85 and the unknown true orientation of the disk: the inclination varies from 85 to about 55, which could mean a variation of 30 or 40, depending on the viewing angle. If we look at both the central disk and the outer rings “from above” or “from below”, the warp angle is about 30; if we look at one “from above” and the other “from below”, the warp angle is about 40. We plot the warp angles with the assumption of viewing the galaxy in such a way that the warp is as weak as possible (30); considering the uncertainties in position angle and inclination, we estimate that the minimum warp angle necessary to reproduce the data cube is about 25. Even with this assumption, this warp remains exceptional: the warp angle, measured at the last point, is very large. In fact it is larger for example than all the warp angles found in the systematic study of galaxy disks in HI by García-Ruiz et al. (garcia-ruiz (2002)). This is probably the reason why this warp is rendered visible: for smaller angles a warp oriented along the line-of-sight is almost invisible.

We have thus found very strong evidence for ESO 123-G23 to be an almost perfectly edge-on disk galaxy with a strong warp along the line-of-sight. This is demonstrated by the comparison of a series of models with the observations. We therefore report the first clear case of a galaxy with such a warp. We note that Becquaert & Combes (1997) interpreted the morphology and kinematics of the HI disk of NGC 891 in terms of a warp close to the line-of-sight. However, the fact that Swaters et al. (1997) interpreted their HI observations of NGC 891 predominantly in terms of a thin + (slowly rotating) thick disk shows that the case was not as evident as the present one. Statistically, the frequency of occurrence of line-of-sight warps should be much higher, while most of them remain undisclosed owing to their small warp angles. The only way to discriminate between different phenomena (warp, flare, etc.) is detailed modelling in the way described here.

###### Acknowledgements.

We thank the referee, Rob Swaters, for valuable comments that improved the quality of this paper. We are grateful to Renzo Sancisi for giving us the first hints for this work. GG wants to thank Christian Brüns and Daniela Vergani for performing the first run of the observations. GG is also grateful for financial support of the Deutsche Forschungsgemeinschaft under number GRK 118 “The Magellanic System, Galaxy Interaction and the Evolution of Dwarf Galaxies”.


format: text

text:
# Transverse momentum and collision energy dependence of high 𝑝_𝑇 hadron suppression in Au+Au collisions at ultrarelativistic energies

STAR Collaboration

## Abstract

We report high statistics measurements of inclusive charged hadron production in Au+Au and p+p collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$=200 GeV. A large, approximately constant hadron suppression is observed in central Au+Au collisions for $`5<p_T<12`$ GeV/c. The collision energy dependence of the yields and the centrality and $`p_T`$ dependence of the suppression provide stringent constraints on theoretical models of suppression. Models incorporating initial-state gluon saturation or partonic energy loss in dense matter are largely consistent with observations. We observe no evidence of $`p_T`$-dependent suppression, which may be expected from models incorporating jet attentuation in cold nuclear matter or scattering of fragmentation hadrons.

High energy partons propagating through matter are predicted to lose energy via induced gluon radiation, with the total energy loss strongly dependent on the color charge density of the medium EnergyLoss . This process can provide a sensitive probe of the hot and dense matter generated early in ultrarelativistic nuclear collisions, when a plasma of deconfined quarks and gluons may form. The hard scattering and subsequent fragmentation of partons generates jets of correlated hadrons. In nuclear collisions, jets may be studied via observables such as high transverse momentum (high $`p_T`$) hadronic inclusive spectra WangGyulassyLeadinHadron  and correlations. Several striking high $`p_T`$ phenomena have been observed at the Relativistic Heavy Ion Collider (RHIC) STARHighptSpectra ; PHENIXhighpt ; PHOBOSNpartScaling ; STARHighptFlow ; STARHighptCorrelations , including strong suppression of inclusive hadron production STARHighptSpectra ; PHENIXhighpt ; PHOBOSNpartScaling . These phenomena are consistent with large partonic energy loss in high energy density matter EnergyLoss ; Suppression ; VitevGyulassy ; SalgadoWiedemann ; HiranoNara  though other mechanisms have been proposed, including gluon saturation in the initial nuclear wavefunction KharzeevLevinMcLerran , attenuation of jet formation in cold nuclear matter TomasikPartonFormationTime , and scattering of fragmentation hadrons GallmeisterEtAl . Additional measurements are required to discriminate among these pictures and to isolate effects due to final state partonic energy loss.

We report high statistics measurements by the STAR Collaboration of the inclusive charged hadron yield $`(\mathrm{h}^++\mathrm{h}^{})/2`$ (defined as the summed yields of primary $`\pi ^\pm `$, K<sup>±</sup>, p and $`\overline{\mathrm{p}}`$) in Au+Au collisions and non-singly diffractive (NSD) p+p collisions at nucleon-nucleon center of mass energy $`\sqrt{s_{_{\mathrm{NN}}}}`$=200 GeV. The Au+Au data extend considerably the $`p_T`$ range of earlier charged hadron suppression studies, and the p+p data are the first such measurement at this energy. Comparisons are made to several theoretical models. The high precision and broad kinematic coverage of the data significantly constrain the possible mechanisms of hadron suppression. In addition, the energy dependence of the yields may be sensitive to gluon shadowing at low Bjorken $`x`$ in heavy nuclei.

We compare the data to two calculations based on hard parton scattering evaluated via perturbative QCD (pQCD-I XNWangPrivateComm  and pQCD-II VitevGyulassy ) and to a calculation extending the saturation model to high momentum transfer KharzeevLevinMcLerran . Both pQCD models for Au+Au collisions incorporate nuclear shadowing, the Cronin effect CroninEffect , and partonic energy loss, but with different formulations. pQCD-I results excluding one or more nuclear effects are also shown. Neither pQCD calculation includes non-perturbative effects that generate particle species-dependent differences for $`p_T`$$`<`$5 GeV/cVitevGyulassyBaryonEnhancement ; XNWangPrivateComm .

Charged particle trajectories were measured in the Time Projection Chamber (TPC) TPC . The magnetic field was 0.5 T, resulting in a factor of three improvement in $`p_T`$ resolution at high $`p_T`$ relative to STARHighptSpectra ; STARHighptFlow . After event selection cuts, the Au+Au dataset comprised 1.7 million minimum bias events ($`97\pm 3\%`$ of the geometric cross section $`\sigma _{geom}^{AuAu}`$) and 1.5 million central events (10% of $`\sigma _{geom}^{AuAu}`$). Centrality selection and analysis of spectra follow Ref. STARHighptSpectra . Background at high $`p_T`$ is dominated by weak decay products, with correction factors calculated using STAR measurements of $`\mathrm{\Lambda }`$(+$`\mathrm{\Sigma }^0`$) and $`K_s^0`$ for $`p_T`$$`<`$6 GeV/c STARLambdaK0sToBePublished  and assuming constant yield ratios $`\mathrm{\Lambda }`$(+$`\mathrm{\Sigma }^0`$)/($`\mathrm{h}^+`$+$`\mathrm{h}^{}`$) and $`K_s^0`$/($`\mathrm{h}^+`$+$`\mathrm{h}^{}`$) for $`p_T`$$`>`$6 GeV/c. The $`\mathrm{\Lambda }`$(+$`\mathrm{\Sigma }^0`$) yield was scaled by a factor 1.4 to account for $`\mathrm{\Sigma }^+`$ decays. Table 1 summarizes the correction factors at high $`p_T`$.

After event selection cuts, the p+p dataset comprised 5 million mainly NSD events, triggered on the coincidence of two Beam-Beam counters (BBCs). The BBCs are annular scintillator detectors situated $`\pm `$3.5 m from the interaction region, covering pseudorapidity 3.3$`<`$$`|\eta |`$$`<`$5.0. A van der Meer scan VernierScan  measured the BBC trigger cross section to be 26.1$`\pm `$0.2(stat)$`\pm `$1.8(sys) mb. The BBC trigger was simulated using PYTHIA PYTHIA  and HERWIG HERWIG  events passed though a GEANT detector model. The PYTHIA trigger cross section is 27 mb, consistent with measurement, of which 0.7 mb results from singly diffractive events. The PYTHIA and HERWIG simulations show that the trigger accepts 87$`\pm `$8% of all NSD events containing a TPC track, with negligible track $`p_T`$-dependence. Non-interaction backgrounds contributed 3$`\pm `$2% of the trigger rate. The high p+p interaction rate generated significant pileup in the TPC. Valid tracks matched an in-time hit in the Central Trigger Barrel (CTB TPC ) surrounding the TPC and projected to a distance of closest approach $`DCA`$$`<`$1 cm to the average beam trajectory. To avoid event selection multiplicity bias, an approximate event vertex position along the beam ($`z_{vert}`$) was calculated by averaging $`z_{DCA}`$ over all valid tracks. Accepted events were required to have $`|z_{vert}|`$$`<`$75 cm, corresponding to 69$`\pm `$4% of all events. The track $`p_T`$ fit did not include the event vertex. The CTB track-matching efficiency is 94$`\pm `$2% and combinatorial background is 2$`\pm `$2%. Other significant p+p tracking backgrounds result from weak decays and antinucleon annihilation in detector material, with corrections calculated using HIJING Hijing  and preliminary STAR measurements. Correction factors at high $`p_T`$ are given in Table 1. For p+p collisions relative to peripheral Au+Au, exclusion of the event vertex from the $`p_T`$ fit results in poorer $`p_T`$ resolution, while the CTB matching requirement results in lower tracking efficiency. The p+p inclusive spectrum was also analysed for $`p_T`$$`<`$3.5 GeV/c by an independent method in which a primary vertex is found and incorporated into the track fit, with consistent results.

The p+p NSD differential cross section is the product of the measured per-event yield and the BBC NSD trigger cross section, and has a normalization uncertainty of $`\pm `$14%. The charged hadron invariant cross section has been measured in $`\overline{\mathrm{p}}+\mathrm{p}`$ collisions at $`\sqrt{s}`$=200 GeV UA1 . The p+p cross section reported here is smaller by a factor of $`0.79\pm 0.18`$, approximately independent of $`p_T`$, where the uncertainty includes the two spectrum normalizations and the correction for different acceptances STARHighptSpectra . The difference is due in large part to differing NSD cross section, which is $`35\pm 1`$ mb in UA1  but is measured here to be $`30.0\pm 3.5`$ mb.

Figure 1 shows inclusive invariant $`p_T`$ distributions of $`(\mathrm{h}^++\mathrm{h}^{})/2`$ within $`|\eta |`$$`<`$0.5 for Au+Au and p+p collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$=200 GeV. The Au+Au spectra are shown for percentiles of $`\sigma _{geom}^{AuAu}`$, with 0-5% indicating the most central (head-on) collisions. Error bars are the quadrature sum of the statistical and systematic uncertainties and are dominated by the latter except at the highest $`p_T`$.

Figure 2 shows $`R_{200/130}(p_T)`$, the ratio of charged hadron yields at $`\sqrt{s_{_{\mathrm{NN}}}}`$=200 and 130 GeVSTARHighptSpectra , for centrality selected Au+Au collisions. Error bars are the quadrature sum of the statistical and systematic uncertainties, dominated for $`p_T`$$`>`$4 GeV/c by statistics at 130 GeV. In the absence of nuclear effects, the hard process inclusive yield in nuclear collisions is expected to scale as $`\mathrm{N}_{\mathrm{bin}}`$, the average number of binary collisions for the respective centrality selection. $`R_{200/130}(p_T)`$ has not been scaled by the ratio $`\mathrm{N}_{\mathrm{bin}}`$(200)/$`\mathrm{N}_{\mathrm{bin}}`$(130), which Glauber model calculations STARHighptSpectra ; Glauber  give as $`1.02`$ for all centralities. Figure 2 also shows the saturation model calculation and pQCD-I calculations for p+p and centrality-selected Au+Au collisions (shadowing-only and full). Both models approximately reproduce the $`p_T`$-dependence of the ratio for Au+Au for $`p_T`$$`>`$2 GeV/c, with pQCD-I slightly better for more peripheral collisions. The various pQCD-I calculations shown illustrate that in this model the reduction in $`R_{200/130}(p_T)`$ for Au+Au relative to p+p is predominantly due to nuclear shadowingXNWangPrivateComm . This sensitivity arises because the shadowing is $`x`$-dependent and at fixed $`p_T`$, different $`\sqrt{s}`$ corresponds to different $`x_T=2p_T/\sqrt{s}`$. The quantitative agreement of pQCD-I with the data improves for more peripheral collisions, suggesting that the prescription for the centrality dependence of shadowing in XNWangPrivateComm  may not be optimal. Alternatively, introduction of $`\sqrt{s}`$-dependent energy loss to the model in XNWangPrivateComm  may also improve the agreement.

Nuclear effects on the inclusive spectra are measured by comparison to a nucleon-nucleon (NN) reference via the nuclear modification factor:

$$R_{AA}(p_T)=\frac{d^2N^{AA}/dp_Td\eta }{T_{AA}d^2\sigma ^{NN}/dp_Td\eta },$$
(1)
where $`T_{AA}`$=$`\mathrm{N}_{\mathrm{bin}}`$/$`\sigma _{inel}^{NN}`$ from a Glauber calculation accounts for the nuclear collision geometry STARHighptSpectra ; Glauber  and we adopt $`\sigma _{inel}^{NN}`$ = 42 mb. $`d^2\sigma ^{NN}/dp_Td\eta `$ refers to inelastic collisions, whereas we have measured the p+p NSD differential cross section. However, singly diffractive interactions contribute predominantly to low $`p_T`$ ZeusDiffraction . A multiplicative correction based on PYTHIA, applied to $`d^2\sigma ^{NN}/dp_Td\eta `$ in Eq. 1, is 1.05 at $`p_T`$=0.4 GeV/c and unity above 1.2 GeV/c.

Figure 3 shows $`R_{AA}(p_T)`$ at $`\sqrt{s_{_{\mathrm{NN}}}}`$=200 GeV for centrality-selected Au+Au spectra relative to the measured p+p spectrum. Horizontal dashed lines show Glauber model expectations STARHighptSpectra ; Glauber  for scaling of the yield with $`\mathrm{N}_{\mathrm{bin}}`$ or mean number of participants $`\mathrm{N}_{\mathrm{part}}`$, with grey bands showing their respective uncertainties summed in quadrature with the p+p normalization uncertainty. The error bars represent the quadrature sum of the Au+Au and remaining p+p spectrum uncertainties. For $`p_T`$$`<`$6 GeV/c, $`R_{AA}(p_T)`$ is similar to that observed at $`\sqrt{s_{_{\mathrm{NN}}}}`$=130 GeV STARHighptSpectra . Hadron production for 6$`<`$$`p_T`$$`<`$10 GeV/c is suppressed by a factor of 4-5 in central Au+Au relative to p+p collisions.

Figure 3 also shows the full pQCD-I calculation and the influence of each nuclear effect. The energy loss for central collisions is a fit parameter, with the $`p_T`$ and centrality dependence of the suppression constrained by theory. The Cronin enhancement and shadowing alone cannot account for the suppression, which is reproduced only if partonic energy loss in dense matter is included. The full calculation generally agrees with data for $`p_T`$$`>`$5 GeV/c if the initial parton density in central collisions is adjusted to be $``$15 times that of cold nuclear matter dEdxColdMatter . pQCD-II exhibits similar agreement for central collisions. In Ref. VitevGyulassy , the pQCD-II calculation was used to predict a $`p_T`$-independent suppression factor in this $`p_T`$ range from the interplay between shadowing, the Cronin effect, and partonic energy loss.

Figure 4 shows $`R_{CP}(p_T)`$, the $`\mathrm{N}_{\mathrm{bin}}`$-normalized ratio of central and peripheral Au+Au spectra. $`R_{CP}(p_T)`$ extends to higher $`p_T`$ than $`R_{AA}(p_T)`$, with smaller point-to-point uncertainties. The error bars show the quadrature sum of statistical and uncorrelated systematic uncertainties. Statistical errors dominate the uncertainties for $`p_T`$$`>`$8 GeV/c. Lines and grey bands are as in Fig. 3. $`R_{CP}(p_T)`$ for $`p_T`$$`<`$6 GeV/c is similar to measurements at $`\sqrt{s_{_{\mathrm{NN}}}}`$=130 GeV STARHighptSpectra , but is now seen to be approximately constant for 5$`<`$$`p_T`$$`<`$12 GeV/c. It is consistent with $`\mathrm{N}_{\mathrm{part}}`$ scaling at $`p_T4`$ GeV/c as reported in PHOBOSNpartScaling , but is below $`\mathrm{N}_{\mathrm{part}}`$ scaling at higher $`p_T`$.

The $`p_T`$-dependence of the suppression in Figure 4 is well reproduced for $`p_T`$$`>`$5 GeV/c by the full pQCD-I and pQCD-II calculations in both panels and the saturation calculation in the upper but not the lower panel. The magnitude of suppression is fitted to the central collision data in the pQCD models but is predicted in the saturation calculation. Attenuation of initial jet formation due to multiple nucleon interactions TomasikPartonFormationTime  generates an increase in partonic $`R_{AA}(p_T)`$ for central collisions of a factor $`2`$ in $`5<E_T<12`$ GeV. A similar $`p_T`$-dependence would be expected for high $`p_T`$ hadrons, in contrast to observations. Suppression in the final state due to in-medium scattering of fragmentation hadrons should also result in a rising $`R_{AA}(p_T)`$ with increasing $`p_T`$ due to the dependence of hadron formation time on the total jet energy GallmeisterEtAl , though detailed comparison of this model to data requires further theoretical development.

In summary, STAR has measured inclusive charged hadron yields from Au+Au and p+p collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$=200 GeV, at higher precision and over a much broader $`p_T`$ range than previous measurements. Large, constant hadron suppression is observed in central nuclear collisions at high $`p_T`$. The systematic behaviour of the suppression at high $`p_T`$ is well described both by pQCD calculations incorporating final-state partonic energy loss in dense matter and a model of initial-state gluon saturation, though the latter model provides a poorer description of peripheral collision data. The isolation of initial state effects on high $`p_T`$ hadron production may be achieved through the study of d+Au collisions at RHIC, allowing a quantitative measurement of final state effects from the data presented here.

###### Acknowledgements.

We thank C. Greiner, M. Gyulassy, D. Kharzeev, C. Salgado, B. Tomasik, I. Vitev and X.N. Wang for valuable discussions. We thank the RHIC Operations Group and RCF at BNL, and the NERSC Center at LBNL for their support. This work was supported in part by the HENP Divisions of the Office of Science of the U.S. DOE; the U.S. NSF; the BMBF of Germany; IN2P3, RA, RPL, and EMN of France; EPSRC of the United Kingdom; FAPESP of Brazil; the Russian Ministry of Science and Technology; the Ministry of Education and the NNSFC of China; SFOM of the Czech Republic, DAE, DST, and CSIR of the Government of India; the Swiss NSF.


text:
# An absorption event in the X-ray lightcurve of NGC 3227

## 1 Introduction

The presence of broad, permitted lines in the optical spectra of type 1 AGN, indicates the presence of dense ($`n_\mathrm{H}>10^8`$ cm-3), relatively low-ionisation gas close to the central engine (e.g. Peterson 1997, Krolik 1999). The covering factor of the ‘broad line region’ (or BLR) is estimated from line equivalent widths to be around 10%, while the strengths of low-ionisation and Balmer lines imply column densities $`N_\mathrm{H}>10^{22}`$ cm-2 (see Peterson 1997 and Krolik 1999 for reviews). ‘Reverberation mapping’ measurements of time delays in the response of optical lines to changes in the optical continuum indicate that the BLR lies within a few light-weeks of the continuum source (e.g. Peterson et al. 1998, Kaspi et al. 2000).

It is perhaps surprising that clearcut evidence for the low-ionisation BLR clouds in type 1 AGN has not previously emerged in the X-ray band. The soft X-ray band ($`<2`$ keV) is particularly sensitive to absorption by gas along the line of sight to the central X-ray source. Numerous studies using missions sensitive to the soft X-ray band (e.g. ROSAT, ASCA, BeppoSAX and most recently the grating instruments on Chandra and XMM-Newton) have shown the presence of absorption due to more highly ionised gas (the ‘warm absorber’) along the line of sight in many Seyfert 1s, but the precise location of this gas is still fairly uncertain (e.g. McHardy et al. 1995, George et al. 1998a, Netzer et al. 2002, Schurch & Warwick 2002). Column density variations of cold X-ray absorbing gas have been reported on time-scales of months-years in both type 1 and type 2 Seyferts on time-scales of months to years (Malizia et al. 1997, Risaliti, Elvis & Nicastro 2002). If such variations are due to clouds crossing the line of sight with Keplerian velocities, the obscuring clouds may lie within or not much beyond the BLR. Unfortunately, due to the sparse temporal sampling of X-ray spectra of AGN available with most X-ray missions, it is difficult to put these apparent absorption changes into context, and assess if they really are caused by the motion of dense absorbing clouds across the line of sight to the X-ray source.

In this letter, we present a study of the long-term X-ray spectral variability of the Seyfert 1.5 galaxy NGC 3227, using data from a monitoring campaign carried out by the Rossi X-ray Timing Explorer (RXTE). A previous study of spectral variability in NGC 3227 showed evidence of an order of magnitude increase in the cold/low-ionisation column (from $`10^{21}`$ cm-2 to $`10^{22}`$ cm-2) between ASCA observations obtained in 1993 and 1995 (George et al. 1998b). Using data from our monitoring campaign, we show that in 2000/2001 NGC 3227 underwent an unusual event lasting $`3`$ months, during which the spectrum became exceptionally hard. The symmetry of this event suggests it is due to absorption by a high column density cloud ($`N_\mathrm{H}3\times 10^{23}`$ cm-2) moving across the line of sight to the X-ray source. This interpretation is confirmed by an XMM-Newton observation obtained serendipitously, early in the event, which shows that the cloud is not substantially ionised. We use dynamical arguments, along with constraints on the ionisation state of the cloud, to infer that the cloud most likely originates in the BLR.

## 2 Observations and data reduction

We have been monitoring NGC 3227 with RXTE since January 1999. The observations were carried out as pointings of duration about 1 ksec. The frequency of the monitoring varied between weekly observations to four observations per day. Here we present data from the Proportional Counter Array (PCA), which consists of 5 Xenon filled proportional counter units (PCUs) and is sensitive to X-rays in the 2-60 keV band. For each pointing we extracted the Xenon layer 1 data from all PCUs that were switched on during the observation, and extracted spectra (and model background spectra) according to the standard method outlined by Lamer et al. (2002). PCA response matrices were calculated invididually for each observation using pcarsp v2.37, taking into account temporal variation of the detector gain and the changing numbers of detectors used.

In order to increase the signal to noise ratio of the PCA spectra, we added all spectra within weekly intervals. Each of these weekly spectra was fitted in xspec 11.0.1 with a model consisting of a power law, a Gaussian emission line, and the Wisconsin absorption model. The energy range used for spectral fitting was restricted to 2-12 keV. Since the data were extracted from varying combinations of PCU’s and the detector gain was changed on several occasions during the time span of the monitoring campaign, simple count rate lightcurves cannot be used to measure the variability of the source. We therefore calculated model fluxes from the spectral fits in the bands 2-10 keV and 8-10 keV in order to obtain detector independent light curves over the full time range of the monitoring campaign.

NGC 3227 was observed by XMM-Newton on 28-29 November 2000. We retrieved the EPIC and RGS pipeline products from the XMM science data archive (XSA). Using the XMM Science Analysis System (SAS v5.3) we produced EPIC MOS and PN spectra for NGC 3227. The extraction radius was 30” and only events in single and double patterns with $`FLAG=0`$ were selected.

## 3 A transient absorption event

### 3.1 RXTE monitoring

During the first quarter of the year 2001, transient hardening of the 2-12 keV spectrum was observed. The event lasted for about 100 days and peaked around MJD 51950. Due to the limited spectral resolution and signal to noise ratio of the spectra the spectral fits alone cannot distinguish whether the hardening is caused by intrinsic slope variation or by additional absorption. When we fit the weekly binned spectra with the power law + Gaussian line model and the absorbing column density set to the galactic value of $`N_\mathrm{H}=2.210^{20}\mathrm{cm}^2`$, we find a flattening of the photon index to $`\mathrm{\Gamma }0`$, harder than so far observed in any AGN X-ray spectrum (see Fig. 1).

We now consider the two simplest possibilities: that the transient hardening is caused by a change in the primary continuum shape, or that it is caused by obscuration by a cloud of large column density. Since the hard-spectrum event occurs during a period when the source long-term X-ray flux is quite low, and the X-ray spectra of a number of Seyfert 1 are known to become harder at low fluxes (e.g. Uttley et al. 1999, Vaughan & Edelson 2001, Lamer et al 2003), we first test whether the spectral shape is simply correlated with flux. If this hypothesis is true, we would expect the relation between photon index and flux to be the same both within the 100 d hard-spectrum event and at other times. Specifically, we need to measure photon index versus 8-10 keV flux (i.e. where absorption effects are minimal), in order to remove any spurious correlation produced by varying absorption, which would naturally cause the flux in softer bands (which are absorbed) to correlate with the fitted photon index. Fig. 2 shows the photon indices of the weekly binned spectra plotted against the 8-10 keV flux. We find that the ’normal’ spectral variability of NGC 3227 is very similar to that found in other Seyfert galaxies. At low flux levels the spectral index steepens rapidly with increasing flux; at higher flux levels the increase of the spectral slope levels off and saturates at $`\mathrm{\Gamma }2.0`$. Similar flux - slope relations have been found in NGC 4051 (Lamer et al. 2002) and MCG-6-30-15 (Shih, Iwasawa & Fabian 2002). On the other hand, the data points from the 2001 event clearly fall outside the normal range in the flux - index plane. The indices are much harder and do not strongly correlate with changes in the 8-10 keV flux.

Having ruled out simple, flux-correlated variability as the source of the hard-spectrum event, we next consider the possibility that it is caused by transient obscuration by a large column of gas. The apparent symmetry of the event supports this possibility, suggesting that a symmetric (i.e. spherical) cloud passed in front of the X-ray source. From the flux of hard X-rays and the normal flux - index relation in NGC 3227 (Fig. 2) we estimate an intrinsic photon index of $`\mathrm{\Gamma }=1.6`$ during the absorption event. We repeated the model fitting with the power law + Gaussian line model with the photon index fixed at $`\mathrm{\Gamma }=1.6`$ and the absorption column density left free to vary. Fig. 3 shows the resulting column density for the assumption of a neutral absorber. We fitted the profile of column density changes with two simple models, a uniform sphere and a $`\beta `$-model of the form $`N_\mathrm{H}=N_{\mathrm{H},\mathrm{max}}/(1+(r/R)^2)^{0.5}`$ (see Fig. 3). The $`\beta `$-model gave a much better fit (reduced $`\chi ^20.98`$ versus $`\chi ^21.85`$ for the uniform sphere model), for a $`\beta `$ index of $`0.5`$, suggesting that the cloud density is not uniform but increases towards the centre. The column density peaks at $`N_\mathrm{H}2.610^{23}`$ cm-2. However, based on the RXTE data alone, we cannot absolutely discount the possibility that the hard-spectrum event is caused by very unusual variability of the primary continuum, which is not simply correlated with flux. Fortunately, NGC 3227 was observed by XMM-Newton early in the event, and it is to this data that we now turn.

### 3.2 The XMM-Newton spectrum

NGC 3227 was observed by XMM-Newton around MJD 51877, during an early phase of the hard spectrum event (see Fig 3). We performed joint model fitting on the EPIC MOS1, MOS2 and PN spectra. The hard bands ($`>`$ 2 keV) can be modeled by a strongly absorbed ($`N_\mathrm{H}510^{22}\mathrm{cm}^2`$) power law spectrum and a narrow Fe $`\mathrm{K}_\alpha `$ fluorescence line. The absorbing column density is consistent with the value measured by RXTE at the same time (see Fig. 3). In the soft range below 1 keV an unabsorbed power law spectrum dominates (see below and Fig. 4).

Using xspec, we fitted 4 different partial covering models to the XMM EPIC spectra:

1. Single power law model with partial absorption: wabs(powl + zwabs(powl + zgauss)), with the two power law components having identical slopes.

2. Double power law model: wabs(powl1 + zwabs(powl2 + zgauss)), both power law indices are free to vary.

3. Double power law model with warm absorber as partial coverer: wabs(powl1 + absori(powl2 + zgauss))

4. Single power law with Compton reflection and warm absorber: wabs(powl + absori(hrefl(powl) + zgauss))

The results of the spectral fitting are summarized in table 1, and the double power law + warm absorber model is shown in Figure 4. In all of the three models a narrow iron $`\mathrm{K}_\alpha `$ fluorescence line with rest frame energy 6.4 keV and equivalent width $``$ 240 eV is required. The neutral, unredshifted absorption in all three models is consistent with the Galactic value ($`2.2\times 10^{20}`$ cm-2). From model 3, we find that constraints on $`N_\mathrm{H}`$ and the ionisation parameter, $`\xi `$ of the absorbing gas suggest that the obscuring cloud is not neutral, but ionised, albeit with a low ionisation parameter (see Fig. 5).

Note that the absorbed power-law component is rather hard ($`\mathrm{\Gamma }1.3`$), as expected due to the low 8-10 keV flux ($`2.210^{12}`$ erg cm-2 s-1) measured during the XMM-Newton observation (see Fig. 2 for comparison with the ‘normal’ photon index-flux relation). The unabsorbed power law component in the soft range has a normalization of about 10% of the absorbed power law component. The best fit parameters of models 2 and 3 (also see Fig. 4) also show that the unabsorbed component is intrinsically softer than the absorbed component. There are three possible reasons for this finding: 1. We see two components emitted from different regions of the AGN, e.g. the softer, unabsorbed power law arises from a more extended emitting region, which is not completely obscured by the cloud while the inner, hard X-ray emitting region is obscured. 2. The absorbed and unabsorbed components are intrinsically identical, but have a steeper slope in the soft range below $`1`$ keV. This could be due to a soft excess or a hard reflection spectrum at higher energies. The reflection spectrum (model 4) fits the data well with a rather high ratio ($`3.2\pm 0.6`$) of reflection covered to escaped component. The large reflected fraction may be due to the low flux state during the time of the XMM observations. 3. The slope of the underlying power law spectrum is highly variable, and therefore the slope of the scattered and delayed component differs from the slope of the aborbed, nuclear component.

We produced XMM EPIC light curves for both the unabsorbed component in the 0.2-1 keV range and for the absorbed component at 2-10 keV (Fig. 6). While the absorbed emission is clearly variable, no significant variability is found in the soft component (a constant model yields $`\chi ^2=36.7`$ for 33 d.o.f.). We therefore conclude that the two components are emitted in different regions of the AGN. Since the soft power-law component does not show any intra-day variability, it could arise from a relatively large volume, and might represent radiation scattered from an extended, ionised medium.

We note that, since submission of this paper the Nov 2000 XMM-Newton observation of NGC 3227 has also been published by Gondoin et al. (2003). They fit a slightly different model to us, we refer readers to Gondoin et al. for a full discussion of their model.

## 4 Physical constraints on the absorbing matter

Assuming that the hard-spectrum event is caused by the passage of an absorbing cloud along the line of sight to the X-ray source, we can place lower limits on the distance of the cloud from the source and its density by making a few basic assumptions.

For simplicity, we shall assume a cloud of uniform density and ionisation parameter, following a circular orbit around the X-ray source. Although our assumption that the cloud has a uniform density is almost certainly not true (see Figure 3), the discrepancy in $`N_\mathrm{H}`$ profiles between the $`\beta `$-model and uniform sphere model is typically only a factor of a few. Therefore we do not expect any significant (i.e. order of magnitude) error in our estimates of cloud parameters to result from this simplifying assumption. The ionisation parameter $`\xi `$ is given by

$$\xi =\frac{L_{\mathrm{ion}}}{n_\mathrm{H}R^2},$$
(1)
where $`L_{\mathrm{ion}}`$ is the ionising luminosity in the range 13.6 eV–13.6 keV, $`n_\mathrm{H}`$ is the hydrogen number density of the gas and $`R`$ is its distance from the ionising source. The diameter of the cloud is given by $`N_\mathrm{H}/n_\mathrm{H}`$, so that the velocity of the cloud across the line of sight is $`v_{\mathrm{cloud}}=N_\mathrm{H}/(n_\mathrm{H}t_{\mathrm{cross}})`$ (where $`t_{\mathrm{cross}}`$ is the crossing time of the cloud across the line of sight). Assuming a circular orbit of radius $`R`$ around the central black hole of mass $`M_{\mathrm{BH}}`$, we can obtain an expression for the hydrogen number density:

$$n_\mathrm{H}=\frac{N_\mathrm{H}}{t_{\mathrm{cross}}}\sqrt{\frac{R}{GM_{\mathrm{BH}}}}.$$
(2)
where $`G`$ is the gravitational constant. Finally, we can combine equations 1 and 2 and rearrange to obtain $`R`$:

$$R410^{16}M_7^{\frac{1}{5}}\left(\frac{L_{42}t_{\mathrm{days}}}{N_{22}\xi }\right)^{\frac{2}{5}}\mathrm{cm},$$
(3)
where $`L_{42}=L_{\mathrm{ion}}/10^{42}`$ erg s-1, $`t_{\mathrm{days}}`$ is the crossing time in days, $`M_7=M_{\mathrm{BH}}/10^7`$ M and $`N_{22}=N_\mathrm{H}/10^{22}`$ cm-2. We next take the conservative 99% confidence limits imposed on ionisation parameter by the XMM-Newton spectrum of $`\xi =`$0.1-30 (see Fig. 5), a cloud maximum column density of $`N_\mathrm{H}=1.8710^{23}cm^2`$ and crossing time of 62 days (the best fitting uniform sphere parameters). The best estimate of black hole mass from reverberation mapping is $`M_{\mathrm{BH}}=510^7`$ M  (Wandel, Peterson & Malkan 1999), and we estimate $`L_{\mathrm{ion}}=10^{42}`$ erg s-1 ($`H_0=70`$ km s-1 Mpc-1, using the observed flux during the event and assuming a typical continuum photon index $`\mathrm{\Gamma }=1.6`$). Combining these values in equation 3, we calculate that $`R`$ lies in the range 9-86 light-days. Using these limits in equation 2, we find a corresponding range of possible densities, $`n_\mathrm{H}0.72.0\times 10^8`$ cm-3, so that the cloud must be of order a light-day in diameter or less. Note that these estimates are not strongly dependent on the central black hole mass, to which no strong lower limit can be set in NGC 3227 (Wandel et al. 1999). A factor 10 reduction in black hole mass will only reduce the lower and upper limits on $`R`$ by 37% and increase $`n_\mathrm{H}`$ by 250%.

## 5 Discussion

We now consider how the physical constraints we can place on the properties of the eclipsing cloud in NGC 3227 compare with the constraints expected for gas in various regions of the AGN. Most obviously, the location of the cloud (roughly 10-100 light-days from the central X-ray source) is consistent with an origin in the BLR of NGC 3227, the size of which is estimated from optical reverberation studies to be between 10-20 light-days (Wandel et al. 1999, Winge et al. 1995, Salamanca et al. 1994), although note that since these estimates refer to the region of the BLR which emits most of the Balmer lines, the full extent of the BLR may be significantly larger. Similarly, the density we estimate for the cloud is also consistent with the lower limit on densities of BLR clouds estimated from optical line ratios (e.g. see Krolik 1999 for discussion) and rules out an origin in the lower density, more distant, narrow line region (NLR) or large-scale molecular torus envisioned by the standard unification model. This last constraint is particularly strong, since even in the absence of any lower bound on ionisation parameter, the density of the cloud can only increase with assumed orbital radius, due to the well constrained total column density and the strong dynamical constraint on cloud size imposed by the cloud crossing time.

By the same token, the cloud cannot be very dense, as according to our equation 2, the density scales only with the square root of distance from the central X-ray source, so if the assumed density were significantly larger than $`10^8`$ cm-3 the cloud would be located much further from the X-ray source than is consistent with any known component of the AGN (the same reasoning can be used to strongly constrain cloud size). Our strong constraint on cloud density is in conflict with the estimate of much denser BLR clouds ($`n_\mathrm{H}10^{11}`$ cm-3) in the luminous Seyfert 1 NGC 5548, by Ferland et al. (1992), who required a large hydrogen number density to prevent the clouds becoming too heavily ionised in the high radiation density environment implied by the relatively small BLR size estimated from reverberation mapping. However, we note that although NGC 3227 has a similar sized BLR to NGC 5548 (e.g. see Peterson et al. 2002), it has a much lower X-ray luminosity ($`10^{42}`$ erg s-1 versus $`5\times 10^{43}`$ erg s-1), so that the radiation density in the BLR of NGC 3227 is likely to be much lower than in NGC 5548. In fact, if for comparison, we convert from our ionisation parameter $`\xi `$ to the equivalent formulation $`U`$ used by Ferland et al. (1992) (e.g. see George et al. 1998a for conversion factors), we find $`U<0.1`$ in NGC 3227, consistent with the value estimated for NGC 5548 and the BLR in other AGN (see Ferland et al. 1992, Peterson 1997, and references therein).

The size inferred for the eclipsing cloud is quite large, implying that only of order $`100`$ clouds are required to produce the $`10`$% covering fraction estimated for the BLR (Peterson 1997, Krolik 1999). This small number of clouds conflicts with the large number ($`>10^4`$) required to produce the smooth broad line profiles observed in type 1 AGN (e.g. Atwood, Baldwin & Carswell 1982), assuming that the cloud intrinsic line widths are small, thermal widths only. This latter assumption might be relaxed if the cloud line widths are intrinsically broadened, e.g. if the clouds are dynamically unstable, which might be expected given the large cloud size. Note also that our cloud size estimates can be used to constrain the size of the primary X-ray emitting region to be less than $`1`$ light-day, or around 170 Schwartzschild radii for a $`510^7`$ M black hole (assuming the primary X-ray emitting region is totally covered, as suggested by the lack of variability in the unabsorbed spectral component).

When considering the implications of these results for models of the BLR in AGN in general, one must bear in mind the caveat that strong observational selection effects almost certainly exist in any detection of BLR clouds by X-ray absorption. For example, eclipses of the X-ray emitting region by much smaller clouds, of the higher densities envisaged by Ferland et al. (1992), would only be noticable if the X-ray emitting region was itself very compact (smaller than hundreds of light-seconds or less than one Schwartzschild radius for a $`510^7`$ M black hole). Furthermore, the duration of such events would be much shorter than the observed $`100`$ d (perhaps only a few thousand seconds), so would only be detectable with continuous monitoring, which we do not have. Therefore we cannot rule out the existence of such small clouds in NGC 3227. Other selection effects may also be at work. For example, the relatively low radiation density in the BLR of NGC 3227 may permit the existence of large, relatively low-density clouds which might not exist in other AGN. The presence of such clouds may also account for the smaller amplitude column-density variations reported in NGC 3227 by George et al. (1998b). Alternatively, the fact that column variations are apparently common in NGC 3227 might be explained if the BLR has a planar geometry and we are viewing the X-ray source through the BLR plane. Even taking into account the above caveats, our results imply that a population of large BLR clouds exists in at least some AGN. Taking the estimate of much smaller, higher density clouds in NGC 5548 (Ferland et al. 1992) at face value, it seems that BLR cloud densities and sizes in different AGN must span at least a three order-of-magnitude range.

## Acknowledgments

We would like to thank Hagai Netzer and Brad Peterson for helpful discussions.


text:
# Detailed comparison of the structures and kinematics of simulated and observed barred galaxies

## 1 Introduction

The dichotomy between the regular and the barred spirals in Hubble’s original galaxy classification scheme is a long standing problem in galaxy evolution and dynamics. It is unclear what causes a spiral galaxy to become barred or not, or more quantitatively, why barred galaxies represent about 30 to 70 per cent (deVaucouleurs 1963; including weak bars, Sellwood & Wilkinson 1993; viewed in IR, Eskridge et al. 2000) of nearby galaxies. This picture becomes even more confusing at high redshifts, where Abraham et al. (1999) have found the fraction to drop beyond a redshift of 0.5. Bars undoubtedly form from dynamical instabilities inherent in self-gravitating axisymmetric discs. Early simulations (e.g. Ostriker & Peebles 1973) have shown that purely self-gravitating discs are instantly susceptible to a bar instability, while the addition of a surrounding spheroidal gravitational potential can prevent bar formation; a more recent study suggests this is merely a slow down of formation, and the bar is actually enhanced at later stages (Athanassoula 2002a). The surrounding stabilizing spheroids can be identified with the dark haloes believed to surround spirals and are responsible for the nearly flat rotation curves of all galaxies. Cold dark matter models of cosmology predict that most spirals are embedded within nearly isothermal haloes (Dubinski & Carlberg 1991; Navarro, Frenk & White 1996, 1997). It is therefore useful to examine the consequences of dark haloes on the formation and evolution of barred galaxies to seek consistency with the prevailing world model and the observable universe.

The pattern speed of the bar is an important although difficult to measure indicator of disc dynamics and dark halo structure. Observations of early-type barred galaxies (Merrifield & Kuijken 1995 (M&K); Gerssen et al. 1999) through the application of the Tremaine-Weinberg method (Tremaine & Weinberg 1984b) (T&W) have shown the pattern speed to be ‘fast’, since the bars end near corotation. N-body models of disc dynamics with static background haloes (summarized in Sellwood 1981) have shown the pattern speed also to be fast. However, soon afterward it was recognized (Tremaine & Weinberg 1984a; Weinberg 1985) that rotating bars would spin down significantly due to dynamical friction from the dark halo. More recent simulations with self-consistent disc and halo distributions clearly showed this effect: angular momentum is transferred from the bar to the halo through torques due to gravitational wakes in the halo, which result in surprisingly low pattern speeds for simulated bars (Debattista & Sellwood 2000 (D&S); Misiriotis & Athanassoula 2000). D&S conducted parametric studies by varying the halo to disc mass ratio of their models and found that a maximal disc yields the least pattern speed slow down. They argued that the observational evidence of fast bars implies a maximum disc in all barred galaxies.

Bar structure and dynamics can also be examined in edge-on systems. The thin bars that form in discs are subject to a buckling instability which causes a bar to bend vertically and thicken into a bulge-like object within a few dynamical times, as was shown numerically by Raha et al. (1991), Pfenniger & Friedli (1991), and Combes & Sanders (1981). Kuijken & Merrifield (1995) (K&M) have shown that the kinematics of peanut-shaped bulges of NGC 5746 and NGC 5965 are consistent with orbits in a bar potential, the theory of which was confirmed by Bureau & Athanassoula (1999) (B&A). The observational test described by K&M can also be applied to simulations, which can be viewed at any angle for more complete results.

Through B and I band observations of 15 galaxies, Elmegreen & Elmegreen (1985) have shown that early- (SB0, SBa) and late-type bars have significantly different properties. Although resonance positions and bar lengths are hard to quantify observationally, there is evidence to suggest that early-type bars are longer, out to the corotation radius, have a flat surface brightness profile along their length, and end beyond the turnover radius of the rotation curve. The late-type bars on the other hand are much shorter, with exponential surface brightness profiles.

In this paper, we revisit the problem of the bar and buckling instabilities in self-consistent disc galaxy models with live haloes. Our goal is to study the pattern speed evolution of numerical bars as well as quantify the structural and kinematic properties for comparison with observed face-on and edge-on barred galaxies. We have taken care to re-scale our models to observed systems for comparison. We extend previous work by introducing a new, more realistic mass model and going to much higher resolution with a simulation containing N=20M particles, a factor of 20 times larger than most work on the subject. Our motivation for using high resolution is to examine the numerical convergence of results, since subtleties of dynamical interactions between discs and haloes may not be captured correctly even with haloes with more than $`1`$M particles (e.g. Weinberg 2001). The disc of our large simulation contains 10M particles, taking us well out of the regime where disc self-heating contaminates results. We are confident that the dynamics of these models reliably represents the gravitational behaviour of real galaxies. There are also disagreements between various simulations of bars in the literature which use different N-body codes and initial conditions. We attempt here to achieve convergence by comparing three of our simulations with those of other groups, as well as to observations, by examining the bar structure, pattern speed, and edge-on kinematics.

## 2 Mass Models

We simulate bars by setting up an initially axisymmetric system including a disc and dark halo which is formally in equilibrium but unstable. Models with rotation curves which are disc-dominated, approaching the maximal disc as defined by van Albada & Sancisi (1986), usually form a bar within a few dynamical times, and the Ostriker & Peebles (1973) empirical criterion is still a useful indicator of the instability.

We use the method of Kuijken & Dubinski (KD) (1995) to generate self-consistent disc+halo models with a nearly flat rotation curve. We consider models with disc-dominated rotation curves with and without a central bulge. The KD models are generated from a distribution function (DF) that is the sum of up to 3 functions: a three-integral disk DF, a bulge DF modelled as King model with an energy cut-off $`E<0`$ and a halo DF that is a flattened King model DF (or lowered Evans model) with the usual truncation at a tidal radius at $`E=0`$. We examine models with and without a bulge as described below. The disk DF is a function of E, $`z`$-angular momentum, $`L_z`$ and a third approximate integral, the $`z`$ energy, $`E_z=\dot{z}^2/2+\mathrm{\Phi }(z)`$. The disk DF is constructed assuming a exponential radial surface density profile and an approximately $`\mathrm{sech}^2\mathrm{z}`$ edge-on profile with fixed vertical scale length $`z_d`$. The disk squared radial velocity dispersion, $`\sigma _R^2`$ is assumed to decline exponentially with the same scale-length as the disk like real galaxies. We can generate N-body realizations of these models and they are formally in equilibrium with an initial virial ratio $`2T/W=1.0`$.

The bulgeless model is nearly a formal maximal disc model, with the rotation curve rising to a roughly flat profile within two scale lengths (Figure 1). Although the halo DF is a King model, the mass profile does not have the usual core since the disk mass dominates the centre. The halo density profile has a mild cusp with $`\rho r^{0.7}`$ to within 0.1 disk exponential scale lengths similar to an NFW profile. Whether or not the haloes of spiral galaxies have a central core is still controversial but our model halo here is consistent with rotation curve decompositions (e.g. Kent 1986). The disc is in formal equilibrium with a Toomre $`Q`$1.1 measured at $`R`$=2.0$`R_\mathrm{d}`$ and is approximately constant in the range $`1.0<R/R_d<4.0`$ rising to higher values beyond these limits. While the model is stable against axisymmetric perturbations it is inherently unstable to bar formation because of the dominance of the disk in the rotation curve. The model is scaled such that $`G`$=$`R_d`$=$`v_{max}`$=1, and the total disc+halo mass is about 12 (disc mass of 1.98 and halo mass of 10.57). We realize these models as N-body systems at high resolution with N=20M particles (10M each in the disc and the halo) and low resolution with N=500K particles (250K each in the disc and halo). Thus, all quoted simulation lengths are in terms of the exponential disc length, and all velocities are fractions of the maximum. We also consider a N=500K model with a small compact bulge of mass $`M=0.32`$ with halo and disk mass profiles essentially unchanged. The Toomre $`Q`$ also remains comparable at $`1.06`$. The code and parameters for generating these models are available by request. Our models differ from D&S who use a Kuz’min-Toomre disc with a vertical Gaussian profile and a lowered polytrope distribution function (DF) for the halo. The recent publication by Athanassoula & Misiriotis (2002) (A&M) contains models based on Hernquist’s (1993) method.

As well as specifying the following results in the original simulation units, the appropriate lengths and times are scaled to match the SBa galaxy NGC 4596 for easier comparison with real galaxies. We use the scale length of $`3.2\mathrm{kpc}`$ and rotational velocity of $`120\mathrm{km}\mathrm{s}^1`$ given by Gerssen et al. (1999). The $`\mathrm{sech}^2`$ scale height for the disc is 0.1 simulation units or 320 pc. For these scales the unit of time is 26.3 Myr. For the runs, we use a time step of $`\mathrm{\Delta }t=0.05`$ (1.32 Myr) for 10,000 time steps, for a total $`t`$=500 (13.2 Gyr) or approximately the Hubble time. We soften gravitational forces using a Plummer law with softening length $`ϵ=0.0125`$ (40 pc) for the halo and $`ϵ=0.005`$ (16 pc) for the disc. All runs are executed with a parallel tree code (Dubinski 1996) running on a PC cluster or a 32 processor Compaq GS320. The 20M particles simulations needed 2 minutes per step on the GS320. The total energy drifts by no more than 1 per cent, and total angular momentum is conserved to typically within 0.1 per cent. The results discussed below are for the high resolution simulation unless otherwise noted.

## 3 Bar Formation and Structural Evolution

Although these models are formally in equilibrium they are strongly unstable. A bar develops by $`t`$=26 (686 Myr), and is sustained for the duration of the run. Visually, the bar grows to its greatest extent at around $`t`$=50 (1.32 Gyr). The bar then bends and buckles between $`t`$=70 and $`t`$=90 (a period of 526 Myr) to leave a slightly slower bar of comparable length. There are isophotal twists seen during the bar buckling phase, but as noted in Shaw et al. (1993) these disappear with the increased heating of the disc seen during the bar evolution. After this short phase, the bar shows elliptical isophotes with a butterfly- or dumbbell-shape when viewed face-on. Fig. 3 shows the face-on views for $`t`$=50 to 500.

The evolution of the bar viewed edge-on clearly reveals the buckling instability. Initially the disc is thin; the formation and motion of the bar heat the disc vertically and azimuthally. At the onset of buckling there is a very short lived phase of about 260 Myr in which the bar bends into an arc before buckling into a boxy shape object (edge-on views are shown in Figs 4 and 5). The boxy-shaped bulge remains after buckling, and slowly evolves into the familiar peanut shape, which is most prominent when the bar is viewed edge-on, with its long side in the plane of the sky. When viewed end-on, the galaxy resembles a disc with a spherical bulge.

The bar remains an obvious component of the disc and shows no sign of collapsing into a spherical bulge. There is an animation of both face-on and edge-on views available at www.astro.utoronto.ca/ oneill.

There is an increase in central density following the initiation of the bar, which influences the initial slope of the rotation curve. Fig. 6 shows the rotation curve at $`t`$=0, 100, 200, and 400, and it shows the increased slope and peak during the simulation. The surface density along the bar is flat at small $`R`$ for $`t`$=50, but rapidly turns exponential following buckling. These data will be discussed in more detail in the following sections.

Even without the gaseous component, this simulated galaxy can be compared with observations, and some obvious features stand out. The face-on views look very much like observed bars, and the edge-on views are similar to several studies comparing peanut-shaped bulges and bars (K&M; Merrifield & Kuijken 1999). More detailed comparison with the latter follows in Section 5.

Recently, A&M conducted a series of tests to compare N-body simulations of different mass models, which are ideally suited to comparison with our results here. They simulated three different models: one disc-dominated (MD), which is similar to ours here, another disc-dominated model with an added bulge (MDB), and one halo-dominated (MH): all created using the Hernquist (1993) models.

Our work agrees quite well with A&M on isophote shape and rotation curve progression for the first half of the simulations, which span approximately the same time. Our evolution here disagrees to some extent with their disc dominated model for the latter half of the run: the mass build up at the centre of the disc causes the inner slope of the rotation curve to become more steep (albeit slightly) and the maximum velocity to increase even during the final stages; whereas, A&M found little difference between the rotation curves for the last quarter of their run. This also carries over to the bar isophote structure, as our increasing central density causes the bar shape to become similar to their MDB model. Therefore, the latter half of our simulation is more readily compared to dumbbell-shaped, face-on and peanut-shaped, edge-on isophotes, which they produce with the MDB model. Our low-resolution model, with a similar number of disc particles as were used by A&M, shows only boxy isophotes, but was run for half as long as the high-resolution model. Since the buckling stage was later in this simulation, the evolution into a peanut-shaped bulge would be expected after the end of the simulation, so it is uncertain whether we would agree with A&M in this case or not. The rotation curve for our low-resolution model does not change slope for the last 100 simulation time units, which may suggest that the mass redistribution does not continue as long as the high-resolution model run.

We also investigated the evolution of the halo. Unlike Weinberg & Katz (2002), the central density distribution of our halo remained virtually unchanged. There were several differences between our models which may account for this. (1) Their bars were much stronger than ours (30 per cent of the disk mass), appearing earlier in the formation of the galaxy out of a cold, gaseous disc (ours formed in the stellar disc of a galaxy initially in equilibrium), which will have a greater dynamical friction effect on the halo. (2) Our halo is based on the King model, which according to Weinberg & Katz (2002) is missing the key resonance for redistributing the halo mass. However, the KD models are composite DFs and the resulting halo has a mild cusp with an effective power law slope of -0.7 at the centre, which is similar to an NFW profile. (3) They stress the need for high resolution, using 4M particles in the halo; they are concerned that tree codes induce more small scale noise and require higher resolution – they did not outline how much more, but the 10M particles in our halo should be enough to see some indication of halo density redistribution.

## 4 Bar Pattern Speed

T&W describe the only method available to measure pattern speeds of real galaxies from observable quantities alone. This method uses luminosity, line-of-sight velocity, and tilt angle to determine the rotation rate of the bar. To date there have been a few early-type galaxy bars measured this way, and every one of them has been found to be fast, with the corotation distance to bar length ratio between 1 and 1.5. This result has been difficult to reproduce in simulations, as all the bars are greatly slowed due to dynamical friction between the bar and the halo. The amount of mass in the halo compared with the disc is still a point of contention, because the existence of a flat rotation curve constrains, but does not pinpoint, this ratio. Studies on this relation refer to ‘maximal’ or ‘sub-maximal’ discs, where the former has the largest disc to halo mass ratio in the inner regions of the disc that is consistent with observed rotation curves (van Albada & Sancisi 1986). D&S used the dynamical friction between the bar and halo to show that only a maximal disc is likely in barred galaxies, since these bars slow down the least. The final state of their maximal disc simulation was nearly a fast bar, with a corotation to bar length ratio of $`1.6\pm 0.3`$.

We used two different methods to determine the pattern speed of the bar ($`\mathrm{\Omega }_\mathrm{p}`$). The first was a direct measurement of bar rotation by determining the orientation angle of the bar as a function of time. Fig. 7 is a plot of pattern speed versus time in the original simulation units. Early on, the pattern speed is 0.44 at $`t`$=26 but there is a rapid reduction during the initial bar growth, and by $`t`$=50 it has slowed to 0.35, then there is another drop during the buckling phase from $`t`$=65 to 95, with the pattern speed settling to 0.33 (the following section compares these with observations). After that initial phase of rapid evolution, the pattern speed is influenced by dynamical friction, as it gradually slows to less than 0.21 by the end of the simulation. Even after approximately a Hubble time, it does not reach a steady state; it is still transferring angular momentum to the halo (see Fig. 8). We note in Fig. 7 that the pattern speed shows significant scatter early on, during bar formation and buckling phases. After formation, the pattern speed smoothly, but quickly, decreases. During buckling, the scatter is probably due to difficulty in precisely measuring bar position in the midst of twisted isophotes. There is also a series of oscillations which continue after the buckling phase which we believe to be real, whose period corresponds to slightly more than the circular period of the disc at corotation. We are uncertain at this time of the cause of these oscillations.

In order to compare these results with both observations and other simulations, they must be scaled to comparable units. We have taken the R (scale length) divided by v (circular velocity) to determine a scaling factor (DS subscripts for Debattista & Sellwood, OD subscripts for O’Neill & Dubinski):

$$\frac{\mathrm{\Omega }_{pDS}}{\mathrm{\Omega }_{pOD}}=\frac{v_{DS}}{v_{OD}}\times \frac{R_{OD}}{R_{DS}}.$$
(1)
In comparison with the maximal disc simulation (run 68) from D&S our scaled pattern speed is slightly lower than theirs for the duration of the simulation, $`0.086\mathrm{rad}\mathrm{t}^1`$ instead of about $`0.1`$ after the bar is fully formed, and $`0.054\mathrm{rad}\mathrm{t}^1`$ instead of around $`0.057`$ at the end of our simulation. Although certainly disc dominated, our mass model would be more correctly compared with D&S control run, which has a similar disc to halo velocity ratio (their $`\eta `$ parameter). Here, our bar pattern speed is slower initially ($`0.093\mathrm{rad}\mathrm{t}^1`$ instead of around $`0.15`$), but does not slow down as quickly or as much as their simulation (final speed $`0.058\mathrm{rad}\mathrm{t}^1`$ instead of around $`0.041`$). Because we are using different simulation techniques and mass models, the agreement here is promising, although it should be noted that our pattern speed degradation is consistently lower than D&S. This agrees with Athanassoula’s (2002b) assertion that bars slow down at a faster rate in colder disks.

We also measured the bar pattern speeds using the observational technique described by T&W. Taking observations from multiple slits laid parallel to the apparent major axis of the tilted disc and finding the luminosity-weighted position and average line of sight velocities determines the pattern speed:

$$\mathrm{\Omega }_p\mathrm{sin}(i)=\frac{<\overline{v_{los}}>v_o}{<x>x_o}$$
(2)
(M&K). This method, although direct, is not simple to execute observationally. Since the basis of this relation is the continuity equation, common observations using gaseous emissions are disqualified, as the gas cycle through stars is inherently discontinuous. The method is therefore most accurate for those early-type barred spiral galaxies with the least amount of gas and star formation. There are also preferred tilt angles to the galaxy to minimize error and difficulties surrounding determining the line-of-sight velocities. It does, however, remove many of the uncertainties inherent in the indirect methods.

Conversely, this technique is relatively simple to perform on simulated data. At various points during the simulation, snapshots were rotated so the bar was at a $`45^{}`$ angle to the x-axis, then tilted about this axis by $`30^{}`$, similar to the observed galaxies listed below. Slits were laid through the centre of the nucleus and at $`\pm `$ 1 scale length, parallel to the major axis. Fig. 9 plots the luminosity weighted velocity versus position, with the assumption of a constant $`M/L`$ ratio.

This method works very well with the simulated galaxy, with speeds agreeing to within around 2 per cent of actual for the duration of the simulation when using the known $`30^{}`$ inclination of the disc. However, in observational analysis the tilt of the galaxy is not known, and these simulations show that the outer isophotes of face-on systems are not exactly circular. Therefore, determining the tilt of the galaxy by assuming this regular shape will lead to errors. Based on the distance to NGC 4596, the outermost detectable isophote is at 5.8 scale lengths. At t=200 (5.4 Gyr), the axis ratio $`b/a`$ at this radius is 0.72. If assumed circular, the computed tilt angle would be $`44.2^{}`$ instead of $`30^{}`$, which would result in a 26 per cent error on the pattern speed as opposed to the 2.6 per cent error found with the correct inclination. These errors systematically underestimate the pattern speed and are reduced as time progresses and the outer isophotes become more regular with the disappearance of any spiral pattern or time-dependent instabilities.

The definition of a fast or slow bar depends on the length of the bar and the positions of the resonances. The former is a point of contention. A&M have shown that any one method of determining bar length will not work equally well for the duration of the galaxy’s lifetime. In this paper we use several methods to determine the length of the bar: the radius where the phase angle of the elliptical isophotes twists, the m=2 component of the Fourier amplitude and the m=2 phase component. The latter two are used by D&S, and all three are used observationally by Aguerri et al. (2003). Originally we had attempted to use Abraham et al.’s (1999) description for finding the bar length by choosing the elliptical isophote with the smallest $`b/a`$ axis ratio (option (i) or $`L_{b/a}`$ in A&M’s comparison of bar length determination methods), but we found this method consistently underrepresented the bar in comparison to the other methods. A large drop in this measure during the last third of the simulation causes the difference to jump from around 25 per cent to 60 per cent shorter, although we find similar ellipticity plots as A&M for their MD model, where the $`1b/a`$ measure reaches a peak before the end of the bar, then gradually declines. The peak occurs at a smaller radius for the later time steps. The first three methods and the resulting bar length determination from each are plotted in Fig. 10. Fig. 11 shows the average bar length using the first three methods as a function of time, with error bars corresponding to the maximum and minimum lengths calculated. The bar length appears to be constant or slightly increasing with time. Also plotted are the corotation to bar length ratios, which are all under 1.5, indicating the bar is fast throughout the galaxy’s evolution. D&S quote a final corotation to bar length ratio ($`D_L:a_b`$) of $`1.57\pm 0.27`$ which they call ‘acceptably fast … only barely so’. There is some agreement between our results and theirs, especially if we only take the Fourier component methods into account, where we find a final ratio of $`1.44\pm 0.1`$.

### 4.1 Comparison with observations

The actual pattern speed of the bar is scaled to the dimensions of real galaxies and compared with observed pattern speeds listed in Table 1.

The T&W method directly measures pattern speed from observables, and is thus more accurate than the other, indirect, methods. It has thus far only been applied to a handful of galaxies, most of which are listed in Table 1. The T&W method depends on continuity, which the life-cycle of gas inherently contradicts, so reliable use of this measurement is restricted to early-type galaxies.

The scaled pattern speeds found in this simulation are initially too slow for all but one galaxy (although two more are just barely within the lower bound of the observed error bars). However, even though the pattern speed decreases with time, the redistribution of the disk particles by the bar causes the scale length of the disk to increase with time as well. Measured with a double exponential fit to the surface density profile, this increase in the disk scale length is fairly rapid initially, jumping to 2.3 from 1.0 between $`t=50`$ and $`t=110`$, and slowly increasing to 2.8 by the end of the simulation. This increase more than offsets the pattern speed slow-down, with only one comparison still too low at $`t=110`$ while two are higher than observed, and the rest are within the accepted ranges. By the end of the simulation three are again too low, only one is too high, and the rest are acceptable (with two at the lower limit of the error bars, and two well within the measured range). Therefore, for a considerable duration, these long-lived bar is rotating at speeds comparable to those observed. This reiterates the corotation to bar length ration finding above: both the simulated and observed bars are ”fast” (ratio less than 1.5).

## 5 Edge-on Kinematics

Some attention has recently been placed on the identification of edge-on barred galaxies. K&M, Athanassoula & Bureau (1999) (A&B), and B&A have studied the line-of-sight velocity profile of peanut-shaped, edge-on galaxies and have shown that their kinematics are explained by the presence of a bar. The characteristic plot of position versus velocity gives a unique figure-eight pattern for barred galaxies, since some orbits are depleted in these systems. Our simulation results have been similarly plotted. After initial bar formation, the plot appears similar to the unbarred plots; as the bar develops, the central positions darken while the density increases, and the areas immediately to the outside of this central pole suffer some depopulation, similar to the observational plots in K&M and A&B. With bar buckling, though, comes increased scatter in the plot, removing any forming pattern, and once again leaving a plot which would be considered unbarred if viewed observationally. Plots at $`t`$=50 and $`t`$=76 are shown in Fig. 12.

The original simulation is composed only of a disc and halo and was constructed to be in equilibrium with a Toomre Q greater than 1. This generates a higher velocity dispersion than would be observed: 0.23 (times maximum velocity) tangential and 0.215 radial dispersion at 2.4 scale lengths. The radial dispersion increases to 0.258 by $`t`$=76. This blurs the figure-eight pattern in the K&M plot. Since all the observational results used cool gas to determine velocities, our results would be much closer to the analytically predicted orbits, and the distinctive figure-eight plot.

## 6 Bulge Model

To counter part of the possible velocity dispersion problem stated above, a model containing a small, compact bulge was constructed in the lower resolution range (500K particles), keeping a similar halo and rotation curve (see Section 2 for a more complete description). The evolution of the disc is similar to the earlier model in that the bar is initiated a little before $`t`$=20 (540 Myr); it differs in that the initial bar is about 20 per cent smaller, with dumbbell-shaped inner isophotes when viewed face-on. The bar remains smaller than the the previous model, by about the same amount for the duration of the simulation. This is in direct contrast with the findings of A&M who found their bulge model to have a longer bar, although the shape is similar. However, their bulge is more massive in comparison to the disk, and more extended than ours. Fig. 13 shows $`t`$=26 to 250. The initial velocity dispersion (0.146 versus 0.23 above for tangential and 0.18 versus 0.215 for radial) was lower than our original simulation; however, the line-of-sight velocity distribution (LOSVD) plots remained similar. The high resolution simulation showed much more detail and followed trends toward the figure-eight shape not seen in the lower resolution plots for the same initial conditions. Creating a simulation of the bulge mass model with higher resolution, with its lowered velocity dispersion, should yield clearer LOSVD plots than the original high resolution run, therefore we cannot make any conclusions as to whether or not the inflated velocity dispersion in the original simulation is responsible for masking the figure-eight pattern.

The pattern speed of the bulge model is higher than the original for the duration of the simulation, is reasonably steady after $`t`$=100 and can be compared with observations, as earlier. The result is that the pattern speed scales up to an acceptable rate for NGC 1023, ISO 139-G009, and NGC 1308, is too low for NGC 4596, and too high for NGC 936, IC 874, NGC 1440, and NGC 3412 for the whole run. The final measured pattern speed for the bulge model (at $`t`$=250) is the same as that for $`t`$=50 in the original disc+halo model. Some speed increase is to be expected with the higher central density, as was mentioned in D&S; however, too much central density means the disc is no longer unstable to bar formation. The bar was even closer to corotation, with a $`D_L/a_B`$ of near one for most of the run, ending with $`1.3\pm _{0.17}^{0.33}`$

## 7 Halo Bar

A slightly elongated shape was excited in our halo, with an axis ratio of around 0.88. If the alignment of this halo bar were to coincide with the disc bar, the dynamical friction should be reduced. Although the $`\mathrm{\Omega }_p`$ deceleration slows after about $`t`$=150, suggesting a decrease in torque on the bar, the halo bar takes much longer to line up with the disc bar. The bar orientations are mostly within $`10^{}`$ of one another after t=450, although not locked in at the same relative angle. D&S also found an m=2 response in the halo which slowly aligns with the disc bar.

The pattern speed degradation during the later time steps may not be solely due to the halo friction. Misiriotis & Athanassoula (2000) have shown that the $`\mathrm{\Omega }_p`$ is slower in thicker discs. As the bar buckles and the simulation progresses, the scale height of the disc increases. This should then also slow the bar, even if the halo bar is locked into synchronous rotation with the disc bar, and the dynamical friction between the disc and halo is at a minimum.

## 8 Comparison of Early and Late Type Bars

Elmegreen & Elmegreen (1985) studied the kinematic and dynamic properties of different Hubble-types. They found distinct photometric differences between early (SB0-SBbc) and late-type galaxies. Early-type galaxy bars appear to end at or near corotation and are twice the length of the turnover radius of the rotation curve. Their luminosity profile is flat compared with the exponential decrease of the rest of the disc. Late-type bars are shorter, seem to end near the inner Lindblad resonance, and are shorter than the turnover radius. They have an exponential luminosity profile with a different slope than that of the disc.

The properties of the bar evolves during our simulation. By assuming a constant $`M/L`$, the density profile of the bar can be taken as representative of a luminosity profile (Fig. 14). Soon after the bar forms, a small section of the bar near the core has a flat density profile, with a near exponential decrease for the remainder of the bar. Later time steps show only exponential profiles. The bar length is near the corotation resonance for the duration of the run, and the bars are always longer than the turnover radius and the inner Lindblad resonance. Fig. 15 contains plots of the Lindblad resonances for the bar: $`\mathrm{\Omega }`$, $`\mathrm{\Omega }\kappa /2`$, and $`\mathrm{\Omega }+\kappa /2`$, which are determined based on an azimuthally averaged redistribution of the particles to make the disk axisymmetric. The bar length criteria agree with an early-type galaxy, but the density profile is late-type. Another point to note is that early-type galaxies are also characterized by a large bulge to disk ratio, whereas our models, at most, have a very small bulge.

The main results we are comparing in this study (evolution of pattern speed, corotation to bar length ratio) may actually be similar for all barred galaxies. Kinematical studies have shown that bars should be fast even in late-type galaxies (e.g. due to dust lane placement Athanassoula 1992, and hydrodynamical simulations Weiner et al. 2001). Therefore, our small or no bulge galaxy with fast bar with an exponential density profile is fully consistent with a real late-type galaxy.

## 9 Discussion

### 9.1 Resolution

Along with the 20M particle simulation, a 500K particle run with the same initial conditions was conducted for comparison. The pattern speed is very similar, although after t=350 the lower resolution model slows a little more than quoted above ($`\mathrm{\Omega }_\mathrm{p}=0.18`$ vs. $`0.20`$). The bar lengths, based only on the $`m=2`$ amplitude and phase change, are in general shorter than in the 20M simulation before t=350, and longer thereafter, showing a distinct trend to increasing the bar length over time. Most of the values are within the error bars of Fig. 11, and there is therefore at least in nominal agreement between the two simulations.

The high resolution simulation showed much more detail in the bar. For example, isophote twists during the formation and buckling of the bar were clearly seen, whereas they were absent in the lower resolution model. The T&W method of finding pattern speed was also more accurate with the larger number of particles. The 500K model results were within 5-10 per cent of the actual pattern speed as opposed to around 2 per cent for the high resolution results.

By the midpoint of the runs, scale height calculations were comparable with the original simulation, but at t=100 the average scale height for the high resolution simulation seems to be higher than the 500K model. This is most probably due to the fact that the bar in the high resolution simulation has already buckled at this point, whereas the low resolution model has not; buckling heats the disc. The scale height is of some interest, as it shows the interaction of the bar with the disc, and an increased scale height has been shown to decrease the pattern speed (Misiriotis & Athanassoula 2000).

Therefore, it can be concluded that the lower resolution is capable of outlining many dynamical properties, but higher resolution is needed for some detailed comparisons with real galaxies.

### 9.2 Code Comparison

Victor Debattista and Jerry Sellwood generously offered the use of their mass model initial conditions to run in our simulation code for comparison. Their simulation was run using a 3-D Cartesian particle-mesh code, with a grid spacing of $`0.2R_\mathrm{d}`$, and a total particle count of 600K (102K in the disc). Here we run the same particles (run 68 in D&S, maximum disc model) with a tree code (Dubinski 1996). We initially had problems with our small softening length, until we changed this to correspond with the grid resolution in their original simulation (equivalent to a Plummer softening of 0.4 scale lengths (1.3 kpc) instead of 0.0125 for the halo and 0.005 for the disc), as their particles are set up in equilibrium with respect to the mesh. They ran many different initial conditions during their research, and we found the $`Q`$=1.5 model (unpublished) to evolve similarly to our code. We found a pattern speed of $`0.095\mathrm{rad}\mathrm{t}^1`$ at $`t`$= 1000 compared to around 0.07 plotted in their paper (their units), and a shorter bar than the original, although when Debattista and Sellwood (private communication) ran the same model in a higher resolution mesh code this result was corroborated. The equilibrium of the particles with the original mesh is most likely the reason for the disparity. The higher resolution leads to stronger forces in the central region of the galaxy, which the initial velocity dispersion does not support, causing a slight collapse of the disc particles to the centre at the beginning of the simulation. This causes the rotation curve to rise more steeply, which leads to a shorter bar (Combes & Elmegreen 1993). The shorter bar feels less friction from the halo, which results in less pattern speed deceleration and the higher pattern speeds found here.

Run 68 in D&S was actually conducted with a $`Q`$=0.05 model. Due to earlier attempts on our models, we were concerned about the low $`Q`$ value and the stability of the disc. We were unsuccessful in reproducing the output from this model, as we were in our own discs with a $`Q`$ much below 1.0. Their results may be due to the way the particles were set up in relation to their grid, or the grid method may be supplying some artificial stability to the model.

A copy of our low resolution initial conditions were also forwarded back to Debattista and Sellwood. Using a new hybrid grid code with softening of $`0.025R_d`$ and a time step of 0.025 they found bar lengths between 2.5 and 3, well within the error bars of Fig. 11. They measured the strength of the bar by calculating the $`m=2`$ amplitude of the whole disc as a function of time. The resulting strength is very similar to ours for the same resolution, however we find a greater bar strength for the 20M run between $`t=150`$ and $`t=450`$. The increase may be due to the larger number of particles available for angular momentum transfer in that model. Although we agreed well until $`t=150`$, they found a significantly higher pattern speed thereafter, as their simulation did not slow down as much as ours. The reason for this is unclear at the moment.

While initial differences were seen, the comparison between these codes has shown most of the results to be reproducible. That said, direct code comparison is very much needed in this field to determine the influence of method versus initial conditions.

## 10 Conclusions and Future Work

The three N-body galaxy models created here generally concur with previous research in structure (face-on and edge-on shapes of bars).

The bars measured here are all fast, with pattern speeds in the range of those observed directly in early-type galaxies. Although the pattern speeds of late-type galaxies, which would be a more appropriate comparison to these simulations, cannot be directly measured, kinematical studies have agreed that their bars are also fast (e.g. van Albada & Sanders 1982; Hunter et al. 1988; Athanassoula 1992; Lindblad et al. 1996; Weiner et al. 2001).

The corotation to bar axis ratio in the high resolution simulation is less than 1.5, which is designated a fast bar. This is marginally consistent with the D&S simulations, although our models did not require a fully maximal disk to remain fast.

A recent preprint by Valenzuela & Klypin (astro-ph/0204028) outlined a simulation which used 780K particles (model B), Hernquist initial conditions, and an adaptive tree-mesh code that has a higher resolution at the centre than the less active outer regions. Their spatial resolution in the central regions is slightly lower than ours, and they use many fewer particles than our high resolution model. Many of their results are in agreement with those presented here: namely, comparable angular momentum transfer, the central density build up created by the bar and its affect on the rotation curve, the lack of alteration of the dark matter density distribution, and the double exponential disk density profile. We disagree in bar length (theirs are much shorter), actual pattern speed (they find a more rapid rotation), and pattern speed degradation (their pattern speed is constant). Our disc to halo mass ratio and Q parameter are different, though, which may account for some of the discrepancies.

By comparing our 20M particle simulation with a 500K particle model, it was determined that the basic dynamics of the system can be reproduced with lower resolution (pattern speed, overall evolution), but if observational techniques are to be used on the results, higher resolution is needed. The techniques used to varying degrees of success here were the T&W method to find pattern speed and the K&M edge-on line-of-sight velocity distribution plots to find bar orbits.

The results presented here agree with much of the previous work done on isolated, bar unstable galaxy simulations. We have shown resolution is no longer a large source of error. We suggest further investigations should concentrate on generating early-type barred spiral galaxies with high resolution low softening codes. Also, more observational data of real pattern speeds, especially for late-type galaxies, would be a great asset to the field.

## 11 Acknowledgments

We would like to acknowledge the assistance of Victor Debattista, Jerry Sellwood, and the anonymous referee for detailed and thought provoking comments, as well as the NSERC and the Canadian Foundation for Innovation.