text:
# X-ray and optical observations of 1RXS J154814.5-452845: a new intermediate polar with soft X-ray emission Based on observations with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member states and the USA (NASA) and on observations collected at the European Southern Observatory

## 1 Introduction

The ROSAT Galactic Plane Survey project (RGPS, Motch et al. 1991) to identify X-ray sources detected at low galactic latitudes ($`|`$b$`|`$ $`<`$ 20°) led to the discovery of new cataclysmic variables (Haberl & Motch 1995; Motch et al. 1996). Cataclysmic variables (CVs) are close binary systems in which the low mass secondary star fills its Roche lobe and transfers matter onto the white dwarf primary. Depending on the magnetic field strength of the white dwarf the accretion flows via an accretion disc or directly couples onto the magnetic field lines of the white dwarf in the case of a strong magnetic field. Magnetic cataclysmic variables are usually divided into polars (AM Her systems) in which the strong magnetic field forces the two stars into co-rotation and intermediate polars with shorter white dwarf spin periods of typically a few hundred seconds (IPs, for reviews of this latter sub-class see Patterson 1994; Hellier 1995, 1999). IPs accrete from the accretion disc via an azimuthally extended accretion curtain which causes high photo-electric absorption.

The X-ray spectra of most intermediate polars show a hard, absorbed bremsstrahlung continuum plus emission lines from highly ionized iron. This emission originates from the matter which free-falls towards the white dwarf creating a strong shock near the surface with temperatures kT up to several tens of keV. The hard X-rays heat the white dwarf surface, creating a soft blackbody-like component as is generally observed in polars (Beuermann & Thomas 1993) where the accretion stream is of pencil-like geometry. Additional heating may be caused by dense blobs in the accretion flow which directly penetrate the white dwarf atmosphere. The temperature of the emission is thought to decrease with increasing size of the accretion region on the white dwarf (King & Lasota 1990). Therefore, the absence of detectable soft X-ray emission in classical IPs may be explained by a too low temperature and/or suppression of soft X-rays due to considerable absorption. The discovery of little absorbed, soft X-ray emission from a small number of ‘soft’ intermediate polars in ROSAT data (Mason et al. 1992; Haberl et al. 1994) together with relatively high magnetic field strengths (e.g. Potter et al. 1997) which collimates the accretion stream is consistent with this picture. Here we present the discovery of a new cataclysmic variable. Optical observations together with detailed X-ray follow-up using ROSAT and XMM-Newton clearly identify 1RXS J154814.5-452845 as new intermediate polar.

## 2 Discovery in the ROSAT all-sky survey

ROSAT detected 1RXS J154814.5-452845 during the all-sky survey (Trümper 1982) with 0.54$`\pm `$0.06 cts s-1 in the PSPC detector (Pfeffermann et al. 1987). The source is listed in the bright source catalogue (Voges et al. 1999) with hardness ratios (defined as HR1 = (H–S)/(S+H) and HR2 = (H2–H1)/(H1+H2) where S, H, H1 and H2 denote count rates in the 0.1–0.4 keV, 0.5–2.0 keV, 0.5–0.9 keV and 0.9–2.0 keV bands, respectively) which indicate a moderately hard X-ray spectrum in the ROSAT band. Table 1 compares the hardness ratios of 1RXS J154814.5-452845 to four other ROSAT-discovered IPs extracted from the RRA PSPC source catalogue from pointed ROSAT observations (ROSAT Consortium 2000). Classical IPs like RX J0028.8+5917 and RX J1712.6–2414 show high absorption, for which HR1 is very sensitive to with values approaching +1. The ‘soft’ IPs RE 0751+14 and RX J0558.0+5353 show little absorption (negative HR1) but also softer intrinsic spectra (negative HR2) caused by a strong soft blackbody spectral component which dominates the spectrum in the ROSAT band. The hardness ratios of 1RXS J154814.5-452845 are between those of the two groups, suggesting less absorption than seen from classical IPs and/or the existence of a soft spectral component.

## 3 Optical observations

### 3.1 Optical identification

The source was included in the RGPS bright source sample for optical identification. We observed the field of 1RXS J154814.5-452845 for the first time on 1997, June 6 with DFOSC and the LORAL 2k x 2k CCD at the ESO-Danish 1.54m telescope. The X-ray source was readily identified with a V $``$ 14.6 blue object exhibiting Balmer, He I, He II and N III-C III $`\lambda \lambda `$4640-60 lines in emission. We show in Fig. 1 the mean of two 10 mn long spectra obtained with a 2.5 ″ slit through grism #7 ($`\lambda \lambda `$ 3860 - 6810 Å; 10 Å FWHM resolution). Fig. 2 shows the position of the optical counterpart on a V band image.

### 3.2 Follow-up optical observations

We organized an optical observation campaign in May 1998 involving both photometric and spectroscopic monitoring at the ESO-Dutch 0.9m and ESO-Danish 1.54m telescopes, respectively. During three consecutive nights from May 7 to 9 UT we obtained optical V band photometry using the standard CCD camera equipped with the TK512CB # 33 CCD. In order to ensure the best compromise between time resolution and observing efficiency the CCD was windowed and the exposure time was set to 30 s. With this instrumental setting, the mean time between two consecutive exposures was $``$ 53 s. A total of 753 frames were acquired over 11.5 h. Raw images were corrected for bias and flat-field using standard MIDAS procedures. Instrumental magnitudes were computed using the SExtractor package (Bertin & Arnouts 1996) for the optical counterpart and 3 comparison stars, named A, B and C whose positions are indicated in Fig. 2. Weather conditions were not photometric during the first and the third night. However, the accuracy of the differential photometry remains good with average errors below the 1% level. From standard star measurements in June 1997 and May 1998 we derive the following magnitudes for the comparison stars; VA = 12.94 $`\pm `$ 0.04, VB = 13.51 $`\pm `$ 0.04 and VC = 15.00 $`\pm `$ 0.04. The photometric uncertainty is dominated by the colour transformation. The mean magnitude of 1RXS J154814.5-452845 recorded in May 1998 was V = 14.61, somewhat fainter than in June 1997, V = 14.46. Spectroscopic observations took place from May 15 till 17 UT. At the ESO-Danish 1.54m telescope, we used the same DFOSC setting as in 1997 but with a smaller entrance slit of 1.5 ″. With grism #7 the wavelength range was $`\lambda \lambda `$ 3860 - 6840 Å with a 5.7 Å FWHM resolution). Because of the smaller slit used in 1998 than in 1997, the spectral resolution improved but no accurate flux calibration could be established. In particular, differential refraction distorted the shape of the flux calibrated spectrum. The first spectrum had an exposure time of 20 mn and the remaining 8 were 15 mn long. Arcs from He and Ne lamps provided the wavelength scale. As for photometric data, MIDAS procedures were used to correct 2-d frames and extract and calibrate spectra in wavelength and flux.

### 3.3 Spectroscopic analysis

#### 3.3.1 The orbital period

Radial velocities of Balmer and of He II $`\lambda `$ 4686 emission lines vary relatively slowly within a single night. The total amplitude is 200 km s-1 for H$`\alpha `$. Balmer and He II velocities seem roughly in phase although there are large differences in mean velocities, higher Balmer lines being more blue shifted than H$`\alpha `$. We list in Table 2 the measured emission line velocities. In order to search for orbital motion we computed periodograms for the various lines using the method outlined in Lomb (1976) and Scargle (1982) which is suitable for unequally spaced data. Spectroscopic data are consistent with a rather long orbital period, larger than 5 h . Because of the strong one day aliasing and sparse coverage it is not possible to unambiguously derive the orbital period. We find that the two most probable periods are P = 9.37 $`\pm `$ 0.69 h and P = 6.72 $`\pm `$ 0.32 h. We show in Fig. 3 the H$`\alpha `$ radial velocity curve folded with a period of 9.65 $`\pm `$ 0.03 h which gives the best sine fit. At this period, He II $`\lambda `$ 4686 velocities have a more structured pattern of variation. However, the overall shape of the radial velocity curve changes with period and clearly more spectroscopic data are needed to determine the true orbital period and velocity structure.

A rather strong indication in favour of a long orbital period is the presence of several metallic absorption lines in the red part of the spectrum. We compare in Fig. 4 the red part of the mean optical spectrum obtained in 1998 with a template K4 III spectrum extracted from the Jacoby et al. (1984) library of stellar spectra. The resolution of the library spectra ($``$ 4.5 Å) compares well with that of our observations. In order to better highlight the similarities we diluted the late type spectrum by a factor 5 using a constant continuum and plotted together the rectified spectra normalized to the average continuum. The Mg ”b” bands and the $`\lambda `$ 5206, 5270, 6362, 6497 Å atomic line blends are visible as well as the Ca I $`\lambda `$ 6162 Å line. The Na I doublet is well marked, however, part of its strength could be due to interstellar absorption as indicated by the presence of a relatively deep $`\lambda `$ 5780 Å band. The large absorption at $`\lambda `$ 6280 Å is probably a mixture of telluric O2 and diffuse interstellar bands. It is not possible to classify the late type spectrum with great accuracy. Torres-Dodgen & Weaver (1993) have identified useful lines for classifying stars in the near-infrared. The absence of TiO bands at $`\lambda `$ 6080 - 6390 Å indicates a spectrum earlier than K7V. The marked $`\lambda `$ 6497 Å blend and Ca I line suggests an early K type. Finally, the presence of the Fe I, Ti I Cr I blend at $`\lambda `$ 6362 Å may indicate an evolved type. Unfortunately, the absorption line depths on individual 15 mn long spectra are not strong enough to measure a reliable radial velocity curve of the late type companion.

A spectral type hotter than K7 and possibly early K is compatible with the long orbital period derived from our spectroscopy. The mean density of the Roche lobe filling star can be expressed as $`\rho 115`$ P$`{}_{}{}^{2}{}_{\mathrm{hr}}{}^{}`$ g cm-3 (Frank et al. 1992). The shortest possible period of 6.72 h implies densities comparable to those of a main sequence K5 star. If the orbital period is 9.37 h a main sequence star would have a G type. Although this cannot be ruled out completely from our data, the requirement for stable mass transfer of q = M2/M1 $`5/6`$ (Frank et al. 1992) would rather suggest a slightly evolved K secondary instead. This picture is consistent with the spectral type distribution of secondary stars in cataclysmic variables compiled by Smith & Dhillon (1998) which shows that no system below 4.3 h period has a secondary earlier than K.

Standard evolutionary sequences for unevolved donor stars undergoing mass-transfer predict a $``$ K2 spectral type if the orbital period of the CV is 6.72h (Baraffe & Kolb 2000). Alternatively, a K type secondary in a 9.37h period CV would only be possible if the mass donor star is slightly evolved (Beuermann et al. 1998). Therefore, in the absence of confirmed orbital period and secondary mass estimates, the evolutionary status of the mass donor star cannot be determined.

#### 3.3.2 Spectral signature of the white dwarf ?

Another outstanding feature of 1RXS J154814.5-452845 is the presence of broad absorption lines at H$`\beta `$ and higher order Balmer lines. Troughs are seen up to H$`ϵ`$ in the blue part of the spectrum (see Fig. 5). In that respect, the spectrum of 1RXS J154814.5-452845 is strikingly similar to that of the ROSAT discovered classical intermediate polar V 709 Cas (RX J0028.8+5917 Bonnet-Bidaud et al. 2001). These authors argue that the broad absorptions are the signature of the non heated part of the accreting star and derive log g = 8 and T = 23,000 K for the DA white dwarf. Similar conclusions may be driven for 1RXS J154814.5-452845. One should however take this interpretation with some caution since similar broad Balmer lines in absorption with superposed narrow emissions are sometimes seen from accretion discs in UX UMa stars (see e.g. Warner 1995, and references therein).

### 3.4 Photometric timing analysis

We show in Fig. 6 the V light curve of 1RXS J154814.5-452845 for the three consecutive nights of observations. Significant photometric activity occurs all along the light curve on time scales of minutes. No conspicuous large long term variations indicating a possible orbital modulation is seen although the hollow at the beginning of the second night may represent some minimum. In order to search for coherent modulations we computed the Lomb-Scargle periodogram of the entire time series. Fig. 7 shows the resulting power spectrum. A clear signal is detected at $`f`$ = 124.72 $`\pm `$ 0.10 cycle/day or P = 692.75 $`\pm `$ 0.56 s. This period is consistent with the most probable X-ray period derived from the long HRI observation (P = 693.3 $`\pm `$ 1.1 s, see below). Although the two closest one day aliases of the photometric period cannot strictly be ruled out on a statistical basis, none of them appears consistent with any of the possible HRI aliases (see Fig. 8). The natural aliasing frequencies of ground based and space borne instruments are obviously different. It is therefore very likely that the optical pulsation occurs at the same frequency as in X-ray and probably represents the spin period of the accreting white dwarf unless the accretion is disc-less as supposed in RX J1712.6-2414 for instance (Buckley et al. 1997). The absence of large power at harmonic frequencies is consistent with the sinusoidal shape of the folded light curve shown in Fig. 9. The full amplitude of the folded V light curve is 0.036 mag. well in the range of that observed in other intermediate polars (see e.g. Patterson 1994; Welsh & Martell 1996).

We investigated further the presence of beat frequencies by applying the Lomb-Scargle periodogram to data with the main modulation at P = 692.75 s removed. The resulting periodogram, plotted in Fig. 10, exhibits an excess power in the interval 120 to 125 cycle/day. The highest peak is at 122.14 $`\pm `$ 0.14 d-1 (P = 707.39 $`\pm `$ 0.21 s). At this given frequency, the formal probability of the peak to be real is 1.7% assuming a white noise local distribution. The value of the peak frequency agrees very well with that expected for the synodic frequency ($`f`$ = 122.14 $`\pm `$ 0.21 d-1) computed assuming the best orbital period of P = 9.37 $`\pm `$ 0.69 h, derived from the optical spectroscopy, and a white dwarf spin period of P = 692.75 $`\pm `$ 0.56 s, derived from the optical photometry. However, an orbital period of $``$ 6.7 h cannot be strictly ruled out by the observations. At the synodic period the light curve is nearly sinusoidal and has a full amplitude of 0.023 mag. The detection of a synodic frequency confirms that the main X-ray and optical period is the true rotation period of the white dwarf.

## 4 ROSAT soft X-ray observations

We proposed two observations of 1RXS J154814.5-452845 in the guest observer program of ROSAT using the HRI detector (David et al. 1996). During the first short observation in 1996, Sept. 13 (RH400792, net exposure time of 1.97 ks) the average count rate was $`0.22\pm 0.01`$ cts s-1. It was dedicated to derive an accurate X-ray position: RA = $`15^\mathrm{h}48^\mathrm{m}14.5^\mathrm{s}`$, Dec = $`45\mathrm{°}28\mathrm{}39\mathrm{}`$ with an error of $`7\mathrm{}`$ dominated by systematic boresight uncertainty. After the optical identification as cataclysmic variable and IP candidate a second ROSAT HRI observation (RH300624) was performed in 1998, Sept. 11-14, with a net exposure of 17.89 ks. The flux with an average count rate of $`0.190\pm 0.006`$ cts s-1 was similar to the first HRI observation and also to the survey detection (the HRI to PSPC count rate conversion factor is 3.0 for a hard spectrum).

The second HRI observation, performed near the end of the ROSAT mission, suffered from unstable attitude and no accurate position could be derived. However, the data can be used for a temporal analysis. The HRI light curve (0.1–2.4 keV) shows sinusoidal intensity variations by a factor of about two. A Z2 test (Buccheri et al. 1983) yields formally a period of 693.3$`\pm `$1.1 s consistent with the value derived from the optical observations and confirming 1RXS J154814.5-452845 as IP, however, the large number of short observation intervals typically lasting less than 2000 s causes many alias peaks. Periods of 682.2 s, 627.5 s and 618.4 s are of statistically similar significance, but are inconsistent with any of the possible alias periods seen in the optical (Fig. 11). The folded soft X-ray light curve obtained from the HRI data using the 693.3 s period is presented in Fig. 12.

## 5 XMM-Newton broad band X-ray observations

XMM-Newton (Jansen et al. 2001) observed 1RXS J154814.5-452845 as part of the telescope scientist guaranteed time during satellite revolution 137. The observation (id 0105460301) on 2000, Sept. 7 lasted for about 20 ks and included the EPIC instruments. The Reflection Grating Systems were switched off due to technical problems and the Optical Monitor was not used due to a bright star in the field. Here we report on the results from the European Photon Imaging Cameras based on MOS (EPIC-MOS1 and -MOS2, Turner et al. 2001) and PN CCD X-ray detectors (EPIC-PN, Strüder et al. 2001) which are mounted behind the three telescopes (Aschenbach et al. 2000). The background in the EPIC detectors was low during the whole observation with the exception of a short interval of $``$400 s when it increased from 1-2 cts s-1 to about 10 cts s-1 per PN CCD. Source and background events were extracted from circular regions with radius 40″ around the source position and a nearby source-free area in the image.

### 5.1 Temporal behaviour

The light curves in the 0.2–8.0 keV band obtained from PN and combined MOS data is plotted in Fig. 13. Clearly seen are the periodic $``$700 s intensity variations on top of a more gradual change on a time scale of about 7000 s. A timing analysis using a Fast Fourier Transform (FFT) technique shows a single peak at the spin frequency with no detectable peaks at harmonic or sideband frequencies originating from the beat between spin and binary orbital period. This is consistent with a sinusoidal variation as is seen in the ROSAT HRI data and a long orbital period in excess of the duration of the observation of $``$20 ks as is suggested by the optical data. The FFT power spectrum is shown in Fig. 14. There is some signature for a modulation on a longer time scale of 4000–8000 s in the power spectrum as was already seen in the light curve. Such a variation might be caused by the binary orbit, but the optical data are more consistent with a longer orbital period. A $`\chi ^2`$ folding analysis was used to derive the best value of 692.3$`\pm `$12 s for the spin period. For comparison with Fig. 11 the periodogram is shown in Fig. 15. The contiguous XMM-Newton observation leads to a unique identification of the spin period, however, the shorter time base of the observation causes a larger uncertainty.

The EPIC-PN data were folded with a period of 692.3 s in different energy bands (Fig. 16). The 0.2–2.0 keV profile is consistent with that obtained from the HRI data (0.1–2.4 keV). The smaller statistical errors in the PN light curve allow to identify two narrow features. A dip before the pulse maximum (or peak on the increasing flank) and an asymmetric minimum are significant. The intensity dip is more prominent at energies 2.0–8.0 keV while the minimum is more symmetric at high energies. The hardness ratio (ratio of the count rates in the 2.0–8.0 keV to the 0.2–2.0 keV band) shows an increase during intensity minimum. The asymmetry in the soft light curve minimum reflects in a deep dip in the hardness ratio, reaching the low level at intensity maximum. The intensity dip is apparently energy independent and not visible in the hardness ratio. The behaviour of increased hardness during intensity minimum is often observed from IPs and is generally explained in the accretion curtain scenario by increased absorption during pulse minimum caused by the curtain. The sharp and short ($`<`$0.2 in phase) decline in hardness ratio near intensity minimum then suggests that the line of sight at this phase is nearly parallel to the curtain, allowing an almost free view down to the white dwarf surface. As the white dwarf rotates the curtain quickly intersects the line of sight causing increased absorption. The hardness ratio dip, the sinusoidal pulse variation with a single peak and the narrow optical emission lines suggest that we see emission from one pole only (if we indeed look down nearly parallel to the curtain the opposite pole in a dipole field should be hidden by the white dwarf at this spin phase). One should mention here that EX Hya also shows sinusoidal pulse variations while we see emission from both poles as indicated by, however, broad optical emission lines (Hellier et al. 1987).

### 5.2 The spin-phase averaged energy spectrum

To further investigate the emission properties of 1RXS J154814.5-452845 spin-averaged energy spectra from the EPIC detectors were accumulated. From the PN data single- (pattern 0) and double-pixel (pattern 1-4) events were selected into two separate spectra to utilize the better energy resolution for single-pixel events. For the two MOS spectra pattern 0-12 events were used. The four spectra were simultaneously fit with the same spectral model only allowing different normalizations between the two MOS and the PN spectra to account for possible inter-calibration problems in absolute fluxes and the different temporal coverage of the instruments (Fig. 13). Spectral response files provided by the hardware teams were used (epn\_ff20\_sY9\_medium.rmf and epn\_ff20\_dY9\_medium.rmf for PN single and double events released in March 2001; m1\_medv9q19t5r4\_all\_15.rsp and m2\_medv9q19t5r4\_all\_15.rsp for MOS released in Febr. 2001). In the past various spectral emission and absorption components were required to fit X-ray spectra of IPs. These include emission from hot thin plasma and soft blackbody radiation (dominating the spectrum of soft IPs in the ROSAT band). In the following different combinations of these models are compared with the observed EPIC spectra. In all models absorption partially covering the X-ray source is applied and an iron fluorescence line is included.

#### 5.2.1 Emission from thin plasma with constant temperature

A spectral model consisting of hot thermal bremsstrahlung suffering absorption by matter partially covering the X-ray source together with an iron line complex (a fluorescent line from nearly neutral iron and two lines from Fe XXV and Fe XXVI) as it was representing X-ray spectra from classical IPs in the past gives no acceptable fit to the EPIC spectra of 1RXS J154814.5-452845. We used the MEKAL model (Mewe et al. 1985) for thermal plasma emission from XSPEC 11.0.1 to account for the bremsstrahlung component and the two lines from highly ionized Fe (model A, see Table 3). This resulted in a reduced $`\chi ^2`$ of 3.5 with a large excess flux below 0.8 keV.

#### 5.2.2 Thin plasma emission plus blackbody component

Improved fit results were obtained by adding a blackbody component to account for the low-energy excess (model B), improving the reduced $`\chi ^2`$ to 1.30. In this model the low absorption component of the partial absorber was also applied to the blackbody emission. Around 1 keV the fit is still not acceptable, again showing excess flux. This excess indicates iron L emission as is also seen in high resolution spectra from the IP RE 0751+14 = PQ Gem (Burwitz et al. 2001).

#### 5.2.3 Two-temperature thin plasma emission plus blackbody component

For iron L emission a lower temperature of the emitting matter is required and we added another MEKAL component (soft thermal kT$`{}_{}{}^{\mathrm{s}}{}_{\mathrm{MK}}{}^{}`$ component in contrast to the hard thermal kT$`{}_{}{}^{\mathrm{h}}{}_{\mathrm{MK}}{}^{}`$ component) to the low absorbed part of the model spectrum (model C). This gives an acceptable fit over the whole energy band with reduced $`\chi ^2`$ of 1.26 for 2113 dof.

#### 5.2.4 Multi-temperature thin plasma emission plus blackbody component

The co-existence of soft (1 keV) and hard (14 keV) thermal emission indicates that we actually observe a wide range of temperatures. A multi-temperature model was used by Done & Magdziarz (1998) to fit ASCA and Ginga spectra from BY Cam. Such a model is available in the XSPEC package, based on the MEKAL model with emission measures which follow a power-law in temperature (the emission measure from temperature T is proportional to (T/Tmax)α). To fit this model we replaced the two MEKAL components by the multi-temperature component, suffering absorption again by matter partially covering the X-ray source. As first approach we fixed the parameter $`\alpha `$ at 1.0 (model D in Table 3). The fit further improved to a reduced $`\chi ^2`$ of 1.21 for 2115 dof (the $`\chi ^2`$ decreased by $``$100 although the number of free parameters was reduced by 2). Alternatively allowing $`\alpha `$ as free fit parameter the maximum temperature is not constrained and was therefore fixed at 60 keV, the value derived with fixed $`\alpha `$. This model (E in Table 3) finally yields the best fit (reduced $`\chi ^2`$ of 1.21) although statistically not different from model D. The spectra (PN singles+doubles, MOS1 and MOS2) are shown in Fig. 17 together with the best fit model E.

As derived from the normalizations of the differently absorbed multi-temperature emission components, the high column density matter covers the plasma by 47% while 53% are absorbed with a much reduced column density. The slope of the optical spectrum indicates that the column density of the low absorption component, which corresponds to E(B-V) $``$ 0.25, could entirely be of interstellar origin. The flux determined from the spectra in the 0.2–10.0 keV band was 2.28$`10^{11}`$ erg cm-2 s-1, 2.34$`10^{11}`$ erg cm-2 s-1 and 2.30$`10^{11}`$ erg cm-2 s-1 for MOS1, MOS2 and PN respectively. PN single and double spectrum were forced to have the same normalization in the fit and therefore yield the same flux. The fluxes derived from the three detectors agree remarkably well, it should be noted however that the used data does not cover identical time intervals. The spectra are not corrected for point spread function losses outside the source extraction radius. The used radius of 40″ should encompass 93% of the flux (Aschenbach et al. 2000) with little energy dependence. The intrinsic 0.2–10.0 keV flux (with absorption set to 0) of the multi-temperature emission component is 3.3$`10^{11}`$ erg cm-2 s-1 and the bolometric blackbody flux is 1.5$`10^{11}`$ erg cm-2 s-1, yielding a ratio for blackbody to plasma emission of $``$46%.

The spectra around the Fe K emission line complex are shown enlarged in Fig. 18. Iron fluorescence, Fe XXV and Fe XXVI lines are clearly resolved in the PN single-events and MOS spectra. The fluorescence line energy derived from the fits is 6.396$`\pm `$0.007 keV (models C-E). The equivalent width (Fe$`{}_{}{}^{6.4}{}_{\mathrm{eqw}}{}^{}`$) is listed in Table 3. The Fe XXV and Fe XXVI line intensities require an under-abundance of iron in models A–C. We therefore allowed to vary the metal abundance in the MEKAL components, but linked all metal abundances in all MEKAL components to one single free parameter, which is mainly determined by the iron lines. The best fit values for the metal abundance varies considerably between models B–C with $``$0.5 solar and the multi-temperature models D and E where it is consistent with the solar value (Table 3). For the equivalent hydrogen densities given in Table 3 solar abundances in the absorbing matter are assumed. At least for the high absorption column, which is intrinsic to the 1RXS J154814.5-452845 system, reduced metal abundances by a factor of two would increase the hydrogen column density by the same factor.

There is indication for a stronger Fe XXV line than the model predicts. Producing more He-like iron emission relative to H-like iron emission requires a lower temperature than the overall fit (which includes the bremsstrahlung continuum) predicts. Fits to RXTE and BeppoSAX spectra of IPs extending to energies up to nearly 100 keV yield reduced temperatures when a hard emission component as arises from reflection on the white dwarf surface is included in the spectral model (de Martino et al. 2001).

### 5.3 Phase resolved spectral analysis

Spectra were accumulated using events from spin phase 0.45–0.95 (spin maximum) and 0.95–0.45 (spin minimum, see Fig. 16). Subtracting the spin minimum from spin maximum spectrum shows an increasing count rate from 0.2 keV up to a maximum around 0.5 keV and then a decline down to 2 keV. This suggests that mainly the two soft components (blackbody and 1 keV MEKAL in model C) change their intensity during the white dwarf spin revolution. In particular the decrease of the difference spectrum towards low energies indicates that absorption changes play a minor role.

A fit to spin minimum and maximum spectra using models C and E shows that temperatures and element abundances do not change. Therefore, these parameters were fixed at the values derived from the spin averaged spectra. The normalizations and column densities derived from the fits using models C and E are compared for spin maximum and minimum spectra (for the case of PN) in Table 4. As suggested by the difference spectrum, the column density of the high absorption component shows only a marginal change between spin minimum and maximum. The parameters for absorption and temperature agree within the errors, although there is an indication for some NH decrease during spin maximum. Also the NH values of the low absorption components are consistent within errors. From the low absorption components the soft thermal component (model C) shows the strongest intensity increase during spin maximum of at least a factor of 2.0 while the blackbody and the hard thermal components vary at most by factors 1.7 and 1.3, respectively. This is illustrated in Fig. 19 where the model components are drawn for spin minimum and maximum. In model E the intensity increase of the soft thermal emission is reflected in a decrease of the power-law index $`\alpha `$ of the temperature distribution.

The equivalent width of the Fe fluorescence line varies very little between spin minimum and maximum. A small decrease during spin maximum is insignificant given the large errors. If the fluorescence line would mainly originate on the white dwarf surface a correlation of its equivalent width with other emission from the surface may be expected. However, the blackbody emission increases during spin maximum suggesting that a large part of the line intensity is probably produced in the cold pre-shock matter of the accretion curtain which strongly absorbs the hard X-ray emission. In this case little changes in the equivalent width are expected as long as the line emitting volume remains visible to the observer. In contrast to 1RXS J154814.5-452845 de Martino et al. (2001) detect variations in the equivalent width of the Fe fluorescence line from RX J0028.8+5917 during a short interval in spin phase (of duration 0.2) which we may miss due to our coarse phase resolution.

Together with the dip in hardness ratios around spin minimum the spectral variations seen between spin minimum and maximum suggest that we are viewing the white dwarf nearly pole-on. If the soft thermal component originates in the accretion curtain below the shock region, the variation in its intensity would then be caused mainly by the varying aspect angle under which the curtain is observed during the spin revolution and suggests beaming of the radiation preferentially out of the sides of the accretion curtains. This is consistent with the decrease in covering fraction of the high column density absorber from 51% at spin minimum to 46% at maximum as calculated from the norm$`{}_{}{}^{\mathrm{h}}{}_{\mathrm{MK}}{}^{}`$ values for model E given in Table 4.

## 6 Discussion

Our extensive optical and X-ray observations establish 1RXS J154814.5-452845 as an intermediate polar with rather unique properties. Optical spectroscopy and photometry indicate an orbital period of 6.72 h or 9.37 h, probably together with GK Per (47.9 h), V 1062 Tauri (9.96 h) and AE Aqr (9.88 h) among the longest known from IPs. The optical spectrum is strikingly similar to that of V 709 Cas (Bonnet-Bidaud et al. 2001) with broad absorption lines at H$`\beta `$ and higher order Balmer lines. Bonnet-Bidaud et al. (2001) argue that these are produced in the non heated part of the white dwarf atmosphere. Therefore, 1RXS J154814.5-452845 may be also a rare case where the three constituents of the system, hot white dwarf, accretion disc and mass donor star can be detected altogether at optical wavelength. This peculiarity opens prospects for an accurate determination of system parameters.

All X-ray observations of 1RXS J154814.5-452845 found the source at a similar X-ray flux level. Using the spectral model parameters derived from the EPIC spectra expected ROSAT PSPC and HRI count rates of 0.400 cts s-1 and 0.146 cts s-1, respectively, were simulated. This shows that the source was brightest during the Sept. 1996 observation, a factor of 1.5 brighter than in Sept. 2000 when it was at the lowest observed intensity level.

To fit the X-ray spectra of 1RXS J154814.5-452845 as obtained by the EPIC instruments (0.15–12 keV) a complex model is required. Thermal emission from a hot multi-temperature plasma in collisional ionization equilibrium together with absorption by matter partially covering the hard X-ray source represents the spectrum above 1 keV. An approach with a two-temperature model suggests that we see a wide range of temperatures covering at least 1 to 14 keV. High resolution spectra of magnetic CVs are consistent with a continuous temperature distribution (Mauche 2001; Burwitz et al. 2001) as is expected from the cooling matter in the post-shock accretion column. The fit of a model with a continuous temperature distribution results in a maximum plasma temperature of 60$`\pm 10`$ keV, relatively unconstrained by the limited bandpass of the EPIC instruments. The temperature of the post-shock region expected from the simple shock model kT = 3/8G$`\mu `$mHMwd/Rwd with the mean molecular weight $`\mu `$ and the hydrogen mass mH, for a 1 M white dwarf with a radius of 0.008 R (for a mass radius relation see e.g. Panei et al. 2000) is $``$60 keV. However this may be regarded only as upper limit because spectral models including other hard spectral components result in decreased shock temperatures (Beardmore et al. 2000). Fitting Ginga LAC spectra of V1223 Sgr with and without taking into account reflection of the hard bremsstrahlung from the white dwarf surface these authors found shock temperatures of 30$`{}_{4}{}^{}{}_{}{}^{+12}`$ keV and 57 keV, respectively. Reducing the shock temperature to 30 keV translates to a white dwarf mass of $``$0.7 M. A more reliable shock temperature and from this a better constrained estimate for the white dwarf mass requires the measurement of the hard X-ray spectrum at energies above 10 keV.

At energies below 1 keV a blackbody component with a temperature kT of 86 eV is needed to fit the EPIC spectra. The bolometric blackbody flux is $``$46% of the intrinsic 0.2–10 keV hot plasma emission. This fraction is lower than that observed from AM Her systems and soft IPs, where the blackbody flux exceeds the bremsstrahlung emission. A correlation of the flux ratio with the magnetic field strength at the main accreting pole suggests increasing suppression of bremsstrahlung in strong magnetic fields (Beuermann & Burwitz 1995). The low flux ratio derived for 1RXS J154814.5-452845 is consistent with a magnetic field strength below 10 MG as generally anticipated for classical IPs. The temperature of 86 eV is higher than that observed from AM Her systems and soft IPs which show typical values of 20–40 eV (e.g. Matt et al. 2000) and 50-60 eV (Duck et al. 1994; Haberl et al. 1994), respectively. Models for the accretion process onto magnetic white dwarfs predict that the effective temperature decreases with increasing size of the area covered by accretion on the white dwarf (King & Lasota 1990). The lack of soft X-ray emission from IPs supported this view and suggested temperatures too low for soft X-ray emission to be detected. The case of 1RXS J154814.5-452845 and in this sense also the soft IPs show an opposite behaviour: the temperature increases from AM Her systems via soft IPs to 1RXS J154814.5-452845, i.e. from high magnetic field strength (small accretion area) to low magnetic field strength (large accretion area). This suggests that either 1RXS J154814.5-452845 is a very peculiar case or that the temperature characterizing the soft emission is determined mainly by other accretion parameters like the mass accretion rate or mass and radius of the white dwarf. Hellier (1999) also discussed that blobby accretion im AM Her systems may lead to a larger footprint on the white dwarf, while in IPs accretion of more homogenous matter from the accretion disc probably results in a narrow radial confinement of the accretion curtain, compensating for the larger range in azimuth.

For an absorbing column density of $`10^{23}`$ cm-2 the expected equivalent width of the Fe fluorescent line should be of the order of $``$100 eV (Inoue 1985). It is however not clear if the spherical geometry of the absorbing matter around the early type star in a high mass X-ray binary as used by Inoue (1985) can be transfered to the situation we have in magnetic CVs with a more cylindrical geometry of the accretion stream. From a tall cylinder with the line of sight along the cylinder axis high absorption is expected but the cylinder contains relatively little matter to produce the fluorescence line. In a flat cylinder the absorption is low but the fluorescence line stronger. In the accretion curtain scenario of IPs the geometry may be closer to the latter case and indeed be the main reason for the high equivalent width of the Fe fluorescence line observed in 1RXS J154814.5-452845 and other IPs. To estimate the Fe line contribution from reflection by cold matter on the white dwarf surface a more detailed modeling of the Fe line production in the accreting matter is required and arguments based on the results of Inoue (1985) sometimes used in the literature should be taken with caution.

The spectral analysis of 1RXS J154814.5-452845 has shown a strong model dependence of the derived metal abundances. While one- or two-temperature plasma emission models result in abundances of 50% solar, the models with continuous temperature distribution yield values compatible with solar elemental abundances. This is likely caused by insufficient spectral resolution in particular at low energies which does not allow to properly disentangle line from continuum emission. To constrain elemental abundances of the accreted matter high resolution spectra obtained from grating spectrometers are required and values derived from broad-band, medium resolution spectra should be treated with caution.

## 7 Conclusions

The follow-up observations of 1RXS J154814.5-452845 at optical and X-ray wavelength have revealed a new intermediate polar with interesting properties. The pulsations with a period of 693 s most likely reflect the rotation period of the white dwarf and the optical data suggest a long orbital period of more than 6 hours. Optical spectra are consistent with a late type companion of spectral type K and the presence of broad absorption features underneath the Balmer emission lines could be the signature of the white dwarf photosphere, similar to the classical intermediate polar V 709 Cas. The improved spectral resolution of the XMM-Newton EPIC instruments together with the large collecting area of the X-ray mirrors allowed a detailed study of the broad band 0.1-12 keV spectrum. The spectrum can be modeled by multi-temperature emission from hot plasma originating in the post-shock accretion region with a maximum temperature of $``$60 keV, an iron fluorescence line with an equivalent width of $``$260 eV and a 86 eV blackbody component. About 47% of the hard emission is absorbed by matter partially covering the X-ray source with a column density of $`10^{23}`$ cm-2 while the remaining part and the blackbody emission are attenuated by a reduced column density of 1.5$`10^{21}`$ cm-2. Intensity variations of individual spectral components between spin minimum and maximum suggest that the spin modulation is mainly caused by viewing effects and only little changes in absorption. The X-ray spectrum of 1RXS J154814.5-452845 and in particular the detection of soft blackbody emission make this IP an interesting case between soft IPs and classical hard IPs.

###### Acknowledgements.

The ROSAT and XMM-Newton projects are supported by the German Bundesministerium für Bildung und Forschung / Deutsches Zentrum für Luft- und Raumfahrt (BMBF/DLR), the Max-Planck-Gesellschaft and the Heidenhain-Stiftung.


format: text

text:
# Monte Carlo Simulation of Ly𝛼 Scattering and Application to Damped Ly𝛼 Systems

## 1 Introduction

The high column density absorption systems (damped Ly$`\alpha `$ and Lyman limit systems) are a key part of understanding galaxy formation. As a galaxy collapses from the highly ionized intergalactic medium, the gas will inevitably go through a phase of a self-shielded cloud of atomic hydrogen before it can cool and collapse further into molecular clouds, form stars, and increase its metallicity. Many of the damped Ly$`\alpha `$ systems (DLAs) at high redshift may be gaseous halos in the process of forming galaxies rather than fully formed spiral disks, especially in view of the very low metallicities at high redshifts (Lu et al. 1996; Prochaska, Gawiser, & Wolfe 2001). If so, the internal structure of these objects should reveal to us the detailed processes by which galaxies form. Several basic questions arise in relation to this: how big are the absorption systems? Is the gas smoothly distributed, or is it in the form of clumps and an interclump medium? If so, how big are the clumps? What is the dynamical state of the gas? How does its mean rotation velocity compare to its velocity dispersion?

Some of these questions may be addressed by studying the associated metal absorption lines (e.g., Prochaska & Wolfe 1997, 1998) or by radio and optical observations of galaxies found to be associated with the absorbers (see Briggs et al. 1989; Djorgovski et al. 1996), although the majority of DLAs at high redshifts might not be associated with luminous galaxies.

An alternative observational probe of the structure of DLAs and Lyman limit systems (LLSs) may be found in their Ly$`\alpha `$ emission produced in hydrogen recombinations and the subsequent scattering of Ly$`\alpha `$ photons through the gas (Hogan & Weymann 1987; Gould & Weinberg 1996). The source of the ionization may be the external cosmic ionizing radiation, shock-heating due to the gravitational collapse of the gas in a dark matter halo, or star formation within the system. In the case of fluorescence of external radiation, the maximum surface brightness is achieved in any system with $`N_{\mathrm{HI}}10^{18}\mathrm{cm}^2`$, where all the incident external photons are absorbed.

Efforts have been made to image known DLAs and LSSs towards high-$`z`$ quasars. DLAs towards PKS0528-250 (Warren & Møller, 1996) and Q0151+048A (Fynbo, Møller, & Warren, 1999) and one LLS towards Q1205-30 (Fynbo, Thomsen, & Warren, 2000) have been successfully detected using tuned narrow band filters. There are also some spectroscopically confirmed detections of Ly$`\alpha `$ emission from DLAs (e.g. Hunstead, Pettini, & Fletcher 1990; Djorgovski 1998; Pettini et al. 1995; Leibundgut & Robertson 1998). Most of these absorber systems with Ly$`\alpha `$ emission detected are near the redshift of the background quasar, which suggests a Ly$`\alpha `$ source due to photoionization by the quasar (Fynbo, Møller, & Warren 1999; Fynbo, Thomsen, & Warren 2000).

Hydrogen Ly$`\alpha `$ line is a resonance line. The problem of radiative transfer of resonance line radiation can be approximately solved analytically under certain conditions (e.g., Harrington 1973, 1974; Neufeld 1990). These analytic solutions can be found only for a limited number of cases such as a static, extremely opaque and plane-parallel medium. In some cases, numerical methods are used to solve the transfer equation (e.g., Adams 1972; Hummer & Kunasz 1980; Urbaniak & Wolfe 1981). For a more general geometry, density distribution, and kinematics, the Monte Carlo method turns out to be very useful and thus it has been applied to many problems (e.g., Auer 1968; Avery & House 1968; Caroff, Noerdlinger, & Scargle 1972; Panagia & Ranieri 1973; Bonilha et al. 1979; Natta & Beckwith 1986; Ahn, Lee, & Lee 2000, 2001, 2002).

This paper presents the results of a numerical Monte Carlo code we have developed to compute the transfer of Ly$`\alpha `$ photons through an arbitrary distribution of hydrogen. The code predicts the two-dimensional (2-D) image and line-spectrum along any given observed direction, for a general three-dimensional distribution of gas of given geometry, Ly$`\alpha `$ emissivity, neutral hydrogen density, and bulk velocity field. In this paper, it will be applied only to spherical and axially symmetric clouds, although in the future it can be applied to results of numerical simulations of galaxy formation.

We describe the Monte Carlo code in §2. In §3, we apply the code to several simple cases of hydrogen clouds to demonstrate the effects of the neutral hydrogen column density, Ly$`\alpha `$ emissivity and bulk velocity field on the spectra of escaped Ly$`\alpha `$ photons. In §4, we model DLAs as static and rotating clouds and investigate their 2-D Ly$`\alpha `$ emission spectra. Finally, we have a brief summary and discussion in §5.

## 2 Monte Carlo Simulation of Ly$`\alpha `$ Scattering

### 2.1 Description of the Method

The method used to numerically compute the images and spectra of Ly$`\alpha `$ emission consists of generating random realizations of the trajectory of a large number of Ly$`\alpha `$ photons as they are scattered within the specified gas distribution. In a nutshell, the method proceeds through the following steps: First, we generate the point where the photon is created, with the emissivity distribution proportional to the recombination rate at every point. Second, the optical depth $`\tau `$ that the photon will travel through before it is scattered is randomly generated with a probability density $`e^\tau `$. The spatial location of the scattering at optical depth $`\tau `$ from the point of emission is then determined along a randomly chosen direction, with the knowledge of the neutral hydrogen density distribution and the scattering cross section (see eq. below). We then generate the velocity of the atom that scatters the photon, as described below, compute the new frequency of the photon after scattering, and generate the new direction of the photon. The process of propagation and scattering is repeated until the photon escapes the modeled system.

As these random realizations of photon trajectories are being carried out, a Ly$`\alpha `$ image of the system along a specified direction can be created. The image consists of a three-dimensional array of the two projected coordinates perpendicular to the direction of the image and the frequency of the escaped photons. The array contains the mean number of photons emitted at every projected position and frequency. At every scattering of a photon, the probability that the photon is re-emitted in the direction of the image and escapes is separately calculated. This probability is added on the element of the image array corresponding to the projected position and frequency of the photon. In practice, for the majority of the scatterings the optical depth is large and the probability of escape is negligible, so one can avoid the computation of the optical depth in the direction of the image on most scatterings.

We now describe how the optical depth of the Ly$`\alpha `$ photons is computed and the velocity of the atom at each scattering is generated. The scattering cross section of Ly$`\alpha `$ photons as a function of the frequency in the rest frame of the hydrogen atom is

$$\sigma _\nu =f_{12}\frac{\pi e^2}{m_ec}\frac{\mathrm{\Delta }\nu _L/2\pi }{(\nu \nu _0)^2+(\mathrm{\Delta }\nu _L/2)^2},$$
$`(1)`$
where $`f_{12}=0.4162`$ is the Ly$`\alpha `$ oscillator strength, $`\nu _0=2.466\times 10^{15}`$Hz is the line center frequency, $`\mathrm{\Delta }\nu _L=4.03\times 10^8\nu _0=9.936\times 10^7`$ Hz is the natural line width, and the other symbols have their usual meaning. For a hydrogen atom with a velocity component $`v_z`$ along the photon’s direction in a fixed “laboratory” frame, the term $`\nu \nu _0`$ in the above equation becomes $`(\nu _i\nu _0)(v_z/c)\nu _0`$, where $`\nu _i`$ is the frequency of the incident photon in the “laboratory” frame. For a Maxwellian distribution of atom velocities, the resulting average cross section is

$$\sigma (x_i)=f_{12}\frac{\sqrt{\pi }e^2}{m_ec\mathrm{\Delta }\nu _D}H(a,x_i),$$
$`(2)`$
where

$$H(a,x)=\frac{a}{\pi }_{\mathrm{}}^+\mathrm{}\frac{e^{y^2}}{(xy)^2+a^2}𝑑y$$
$`(3)`$
is the Voigt function, $`\mathrm{\Delta }\nu _D=(v_p/c)\nu _0`$ is the Doppler frequency width, $`v_p=(2kT/m_H)^{1/2}`$ is the atom velocity dispersion times $`\sqrt{2}`$, $`T`$ is the gas temperature, $`m_H`$ is the mass of the hydrogen atom, $`x_i=(\nu _i\nu _0)/\mathrm{\Delta }\nu _D`$ is the relative frequency of the incident photon in the laboratory frame, and $`a=\mathrm{\Delta }\nu _L/(2\mathrm{\Delta }\nu _D)`$ is the relative line width. The optical depth along the photon trajectory is computed by integrating the cross section in equation (2) times the atomic hydrogen density.

Once a spatial location where the photon is scattered has been determined, the velocity $`v_z`$ along the direction of the incident photon of the atom responsible for the scattering obeys the following distribution:

$$f(u_z)=\frac{a}{\pi }\frac{e^{u_z^2}}{(x_iu_z)^2+a^2}H^1(a,x_i),$$
$`(4)`$
where $`u_z=v_z/v_p`$. We describe how we generate random numbers with this distribution in the Appendix. The two velocity components perpendicular to the direction of the photon simply follow a Gaussian distribution parameterized by the local temperature.

After the total velocity of the atom is chosen according to the above distribution, we first perform a Lorentz transform of the direction and frequency of the photon to the rest frame of the atom. The direction of the scattered photon is then generated according to a dipole distribution in this frame. Its frequency differs from the incident one by the recoil effect which is taken into account in the code (it is negligibly small in the applications we will present in this paper). The direction and frequency of the scattered photon is then transformed back to the laboratory frame, a new optical depth is chosen, and the entire process is repeated until the photon escapes from the cloud. In the presence of a fluid velocity, the procedure changes only by replacing the incident frequency $`x_i`$ in equations (2) and (4) by its value in the fluid frame, $`x_{\mathrm{f}i}=x_i(v_{\mathrm{f}z}/c)(\nu _0/\mathrm{\Delta }\nu _D)`$, where $`v_{\mathrm{f}z}`$ is the fluid velocity parallel to the photon’s direction (valid in the non-relativistic regime).

As mentioned previously, the Ly$`\alpha `$ image and spectrum of the system are calculated as the photon trajectories are randomly generated. At every photon scattering, we compute first the optical depth for the photon to escape the cloud along the direction of the image, and the photon is added to the corresponding element of the image array in projected position and frequency, with the weight $`e^\tau (1+\mu ^2)\mathrm{d}\mathrm{\Omega }`$, where $`\tau `$ is the optical depth for escaping, $`\mu `$ is the cosine of the angle between the incident photon and the image direction, and $`\mathrm{d}\mathrm{\Omega }`$ is the solid angle subtended by the pixel (both $`\mu `$ and $`\mathrm{d}\mathrm{\Omega }`$ are measured in the rest frame of the atom). The factor $`1+\mu ^2`$ accounts for the dipole probability distribution of the direction of the photon after scattering. A photon also has a probability to directly escape the system along the direction of the image when it is created which is added to the corresponding element of the image array. After the contribution to the image has been added, a random direction of the new scattered photon is generated according to the same dipole distribution, and the realization of the photon trajectory is continued until the photon finally escapes from the cloud. The total number of photons in the simulation is chosen to give an acceptable signal-to-noise ratio (the typical number of photons used for simulations in this paper is on the order of $`10^310^4`$).

### 2.2 Tests of the Code

We test our numerical code for the case of a static, plane-parallel slab, for which Neufeld (1990) derived an analytic solution of the mean intensity in the limit of a large optical depth. We simulate the case of a midplane source radiating line-center Ly$`\alpha `$ photons in a plane-parallel slab of uniform density for three optical depths: $`\tau _0=10^4,10^5,`$ and $`10^6`$, where $`\tau _0=\sigma (x_i=0)N_{\mathrm{HI}}`$ is the line-center optical depth from the midplane to the boundary of the slab. We assume a temperature $`T=10K`$ ($`a1.49\times 10^2`$). As we have mentioned, the recoil effect is taken into account in the code. For the purpose of comparison, we first turn off this effect and run the simulation. We compare the spectra of the escaped Ly$`\alpha `$ photons with the analytic results (eq.\[2.24\] in Neufeld 1990; note that his definition of $`\tau _0`$ differs from ours by a factor of $`\sqrt{\pi }`$ as $`a1`$. See also Ahn, Lee, & Lee 2001). The left panel of Figure 1 shows excellent agreement, becoming better as the optical depth increases as expected since the analytic solution applies in the limit of a large optical depth ($`\sqrt{\pi }\tau _010^3/a`$, Neufeld 1990). We then turn on the recoil effect and rerun the simulation. The results are shown in the right panel of Figure 1. The double-peaked frequency distribution of the escaped photons becomes asymmetric – more photons appear at the red part. The profile from analytic solutions can be regarded as an average of the red part and the blue part. The magnitude of the recoil effect is easily understood as reflecting the thermalization of photons around frequency $`\nu _0`$ inside the cloud (Wouthuysen 1952; Field 1959), which modifies their abundance by a factor of $`e^{x/x_T}`$, where $`x\mathrm{\Delta }\nu /\mathrm{\Delta }\nu _D`$ and $`x_TkT/(h\mathrm{\Delta }\nu _D)=(m_Hc^2kT/2)^{1/2}/(h\nu _0)`$. Obviously, the recoil effect becomes more important for low temperature and large optical depth.

We also test the code for a spherically-symmetric case with a particular bulk velocity field. Loeb & Rybicki (1999) study and simulate the scattered Ly$`\alpha `$ radiation around a point source before cosmological reionization. The Ly$`\alpha `$ photons are scattered by the intergalactic medium (IGM) which is undergoing Hubble expansion. According to their approximation, the gas temperature is regarded as zero. We run a corresponding simulation with our code and find that the frequency distribution and the surface brightness profile of the escaped Ly$`\alpha `$ photons shown in Loeb & Rybicki (1999) are well reproduced.

## 3 Spherical Clouds of Uniform Density

We start applying our code to the most simple case of spherical hydrogen clouds with uniform density. We consider the following cases:

Case 1: uniform emissivity, static cloud

Case 2: uniform emissivity, expanding or contracting cloud

Case 3: central point source, static cloud

Case 4: central point source, expanding or contracting cloud

For each case, we perform two runs with different column densities ($`2\times 10^{18}\mathrm{cm}^2`$ and $`2\times 10^{20}\mathrm{cm}^2`$, corresponding to line-center optical depths $`\tau _0=8.3\times 10^4`$ and $`8.3\times 10^6`$, respectively). The temperature is assumed to be $`2\times 10^4`$ K everywhere. In the case of expansion, we set it to be Hubble-like (i.e., the velocity is proportional to the radius), with the velocity at the edge of the system fixed to be $`200\mathrm{km}\mathrm{s}^1`$. The cases of contraction are done in the same way; their line spectra can generally be obtained by simply mirror-reflecting in frequency the photons escaped from the expanding cloud about the rest frame Ly$`\alpha `$ frequency.

Figures 2 and 3 show the distribution of escaped photons for all the cases. For static cases, as expected, the escaped photons have a double-peaked frequency distribution. Except for the very small effect of atomic recoil (see Field 1959), the two peaks are exactly symmetric; the differences in the figures are due to the simulation noise.

The results obtained for the emergent spectrum can be understood by noting the typical trajectory in frequency and space followed by a photon (e.g., Osterbrock 1962; Adams 1972; Urbaniak & Wolfe 1981; Ahn, Lee, & Lee 2000). In the case of a Lyman limit system (like in Fig. 2), the photons are always most likely to be scattered by atoms that have the right velocity along the photon direction so that, in their rest frame, the photon frequency is shifted very close to the line center. This implies there is practically no correlation between the frequency of the photon before and after scattering. At any random time, a photon is most likely to be found within a frequency range $`\mathrm{\Delta }\nu _D=(v_p/c)\nu _0`$ of the line center, and will occasionally make larger excursions away from the line center when scattered by an atom having a large velocity perpendicular to the photon’s direction. After any such excursion, the photon will likely return immediately to the line center with one or a few scatterings. Therefore, the photon escapes the cloud not by diffusing in frequency, but by making a single large jump when scattered by an atom in the high-velocity tail of the Maxwellian distribution, which reduces its optical depth for escaping to a value of order unity. For a Lyman limit system with $`N_{\mathrm{HI}}=2\times 10^{18}\mathrm{cm}^2`$, the line-center optical depth is $`\tau _0=8.3\times 10^4`$, and the photons escape when their optical depth is $`\tau =\tau _0e^{x^2}1`$ (where $`x\mathrm{\Delta }\nu /\mathrm{\Delta }\nu _D`$), which implies $`|x|\sqrt{\mathrm{ln}\tau _0}=3.4`$.

In the case of a DLA system, the scattering history differs from the Lyman limit case once the photon reaches a sufficiently large excursion in frequency to make it more likely that the next scattering is caused by a random atom far from the line center, rather than an atom moving with the right velocity to shift the photon frequency to the line center in the atom frame. Then, the evolution of the photon is described by diffusion in frequency. Since the photons are now reaching their escaping frequency not by a large jump but by a series of small steps, they undergo greater spatial diffusion than in Lyman limit systems. This, plus the fact that the optical depth has now a power-law instead of Gaussian dependence at large $`x`$, broadens the width of the two emission peaks. Of course, the emission peaks also move further from the line center because a smaller scattering cross section is required for the photons to escape.

In Figure 4, we can see that the surface brightness produced by a central point source is more extended for a DLA system than for a Lyman limit system, owing to the greater spatial diffusion.

For an expanding cloud, photons should escape on average with a redshift because they are doing work on the expansion of the cloud as they are scattered. For a Lyman limit system with a central source, a photon undergoing a large positive frequency jump will be moved back to the line center relative to the fluid as it travels through the cloud, owing to the Hubble-like expansion of the cloud. On the other hand, a negative frequency jump will allow the photon to escape directly. In our case, the velocity at the cloud edge, $`200\mathrm{km}\mathrm{s}^1`$, is much larger than the atomic velocity dispersion, so a photon needs to undergo many positive frequency jumps in order to diffuse spatially through the cloud and be able to escape on the blue side of the line, which explains why the blue peak is highly suppressed (see Fig. 3 of Ahn & Lee 1998 for a similar explanation). The situation is exactly reversed for a contracting cloud. The case of uniform emissivity broadens the line, essentially because of the different velocities of the emission sites of the photons; there is also the additional effect that photons emitted near the edge of the cloud are more likely to escape on the blue (red) peak for an expanding (contracting) cloud. In DLAs, the line is also broader for a point source compared to a Lyman limit system because of the power-law dependence of the cross section on frequency and the greater degree of spatial diffusion.

To summarize, these simple models show how the frequency and spatial distributions of Ly$`\alpha `$ photons escaped from a cloud are related to the Ly$`\alpha `$ emissivity distribution, the bulk velocity field and the column density of the cloud.

## 4 Ly$`\alpha `$ Emission from Model Damped Ly$`\alpha `$ Systems

We now apply our code to a gas cloud with an isothermal density profile as a model for DLAs. The nature of DLAs is still a subject of debate. They could be protogalaxies with a rotating disk component (e.g., Prochaska & Wolfe 1997, 1998) or spherically distributed clouds of gas moving randomly in halos (e.g., Mo 1994; Haehnelt, Steinmetz, & Rauch 1998; McDonald & Miralda-Escudé 1999). The other mechanism that may give rise to DLAs is large-scale outflows due to galactic winds in dwarf galaxies (e.g., Nulsen, Barcons, & Fabian 1998; Tenorio-Tagle et al. 1999; Schaye 2001). In this paper, we focus on the first two pictures. We investigate the emergent spectra and spatial distribution of the Ly$`\alpha `$ emission, to show what it can tell us about the internal structure of the system.

### 4.1 Model Description

We model the gas cloud producing a DLA assuming that it forms in a dark matter halo of mass $`10^{11}M_{}`$, virialized at redshift $`z=3`$. In an $`\mathrm{\Omega }_M=1`$ universe with $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, the corresponding virial radius is $`r_{\mathrm{vir}}=24\mathrm{kpc}`$, and the virial velocity $`V_{\mathrm{vir}}=104\mathrm{km}\mathrm{s}^1`$ (e.g., Padmanabhan 1993). We assume the gas density profile is a singular isothermal sphere with a cutoff at the virial radius and a fraction of the halo mass in gas of 0.05.

We consider two different cases as an illustration to show how the rotation of a system could be probed by observations of Ly$`\alpha `$ emission. The first case is a spherically symmetric cloud, and the second a rotating oblate ellipsoid with an axis ratio of $`0.5`$ for the gas distribution. The rotating velocity is set to be $`V_{\mathrm{rot}}=\sqrt{2/3}V_{\mathrm{vir}}`$. In one case, we assume the gas to be static and at a temperature $`2\times 10^4`$ K (a typical temperature of gas that cools after shock-heating but stays photoionized). We also consider another case where the gas is given an additional velocity dispersion of $`V_{\mathrm{vir}}/\sqrt{3}`$ in the spherical model, and $`V_{\mathrm{vir}}/3`$ in the flattened, rotating model. This additional velocity dispersion is simply added quadratically to the thermal one, which would be valid if the gas were in optically thin clumps moving at this velocity dispersion. In practice, any relevant gas clumps in DLAs will be optically thick to Ly$`\alpha `$ photons, but their inclusion would make our model much more complex. We will discuss the effect we would expect from optically thick clumps in §5.

There are various sources to produce Ly$`\alpha `$ photons: internal dissipation, star formation, and fluorescence caused by the intergalactic UV background. The external UV background will give rise to Ly$`\alpha `$ fluorescence; at $`10^4`$ K, in an optically thick medium, recombination of the photoionized hydrogen has a 68% probability of producing a Ly$`\alpha `$ photon (Spitzer, 1978). Shock-heating of the gas will dissipate kinetic energy into heat, which will result in Ly$`\alpha `$ emission after line excitation and collisional ionization followed by recombination. In addition, internal star formation can of course also produce Ly$`\alpha `$ emission. We consider only fluorescent emission from the external background in this paper. In this case, the emissivity is proportional to the square of the ionized gas density. The other two sources of emission are more centrally concentrated, and therefore similar differences as the ones between uniform emissivity and central source in the previous section can be expected.

We take self-shielding into account in the spherical case, calculating the neutral hydrogen density of the cloud at every radius, where the total gas density profile is singular isothermal. We use an iterative algorithm similar to the one presented by Tajiri & Umemura (1998), as described in Zheng & Miralda-Escudé (2002). Here, we assume a UV background intensity of $`10^{22}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2\mathrm{s}^1\mathrm{Hz}^1\mathrm{sr}^1`$, constant at all frequencies between the H I Lyman limit, $`\nu _L`$, and the He II Lyman limit, $`4\nu _L`$. The intensity is set to zero at frequencies above $`4\nu _L`$. This is a good approximation since most higher frequency photons are likely to be absorbed by He II before reaching the hydrogen self-shielded zone. Even though in reality some very energetic photons exist, which can penetrate into the self-shielded region, their effect on the profile of neutral hydrogen is very small due to their low intensity. The profile of neutral hydrogen that is obtained is similar to that shown in Zheng & Miralda-Escudé (2002). In the case of the flattened, rotating gas distribution, we simply take the same neutral hydrogen density profile as for the spherical case and flatten it with an axis ratio of $`0.5`$ (the numerical code we have developed for computing the self-shielding correction applies only to spherically symmetric systems). This is of course not exact, although a much better approximation than neglecting self-shielding altogether.

### 4.2 Results of Simulations

The surface brightness profile for the spherically symmetric, non-rotating cloud is plotted in Figure 5. The thick solid line shows the case of including the fluid velocity dispersion, $`\sigma =V_{\mathrm{vir}}/\sqrt{3}=60\mathrm{km}\mathrm{s}^1`$, while the thin solid line is the result with only the thermal velocity dispersion, $`(kT/m_H)^{1/2}=12.8\mathrm{km}\mathrm{s}^1`$. The dashed line is what is obtained when self-shielding is not included (it has been shifted vertically). The velocity dispersion has practically no effect on the surface brightness profile, because the change in spatial diffusion of the photons is very small. The emissivity due to the recombination has a “hole” in the inner part of the cloud because of self-shielding. Therefore, we see a “core” in the surface brightness profile (with a slight decline of surface brightness toward the center in the inner part due to geometric effects). The outer part of the profile follows the optically thin case. The core of the surface brightness profile is a signature that the emission is due to fluorescence of the external background, since other sources of emission are more centrally concentrated.

Because the ionizing photons of the external background extend only up to $`4\nu _L`$ in frequency, essentially all the photoionizations and subsequent recombinations occur in the outer region of the gas cloud with a hydrogen column density of $`N_{\mathrm{HI}}10^{19}\mathrm{cm}^2`$, the inverse of the photoionization cross section at frequency $`4\nu _L`$. The damped absorption wings are not yet important at this column density, therefore the photons escape the cloud after a single large jump in frequency when scattered by an atom in the high-velocity tail of the Maxwellian distribution, as discussed in §3 for the case of Lyman limit systems.

Figure 6 shows the 2-dimensional spectra of the Ly$`\alpha `$ photons escaped from the non-rotating cloud. We assume that the slit is large enough to include the entire cloud. The double-peaked distribution of escaped Ly$`\alpha `$ photons, with the intermediate frequency range of essentially zero flux, appears at every spatial position.

Including an additional velocity dispersion has the effect of separating the two peaks. This is true in our model because we consider the case where there is no correlation between the velocity of two atoms that are spatially close, valid only when any gas clumps moving coherently are so small that they are optically thin to Ly$`\alpha `$ photons. In other words, the additional velocity dispersion is effectively thermal. In the presence of coherent motions of optically thick regions, the intermediate frequency range with a highly reduced flux disappears, as discussed in §3 in the examples of an expanding or contracting cloud.

Spectra of the flattened rotating cloud depend on the viewing angle and slit orientation. The face-on spectrum, plotted in Figure 7, is very similar to that of the spherical static case. Figure 8 shows the edge-on spectra, with the slit perpendicular (top panels) and parallel (bottom panels) to the rotational axis. The spectrum is averaged along the spatial coordinate perpendicular to the slit; in other words, we are assuming that the entire cloud is included in the slit, but that the wavelength dispersion is wide enough that the size of the cloud does not introduce any smoothing of the spectrum. When the slit is perpendicular to the rotation axis, the rotation curve pattern is clearly seen. The velocity shift of the spectrum is about the same as the rotational velocity, $`V_{\mathrm{rot}}=V_{\mathrm{vir}}\sqrt{2/3}=85\mathrm{km}\mathrm{s}^1`$, at the edge of the system \[note that the unit used in the horizontal axis is $`\mathrm{\Delta }\nu _D=(2kT/m_H)^{1/2}\nu _0/c=18.2\mathrm{km}\mathrm{s}^1(\nu _0/c)`$\]. The variation of the shape of the spectrum along the slit is easily understood by thinking of an analogy to the expanding and contracting cloud in §3. The Ly$`\alpha `$ photons are produced on the outer ionized layer of the system, near the radius $`r_{\mathrm{ss}}`$ where the hydrogen becomes self-shielding. The line-of-sight velocity at this radius is proportional to $`V_{\mathrm{rot}}r_p/r_{\mathrm{ss}}`$, where $`r_p`$ is the projected radius on the slit position. Therefore, the central trough of the spectrum shifts linearly with the projected radius. Moreover, the variation of the velocity along the line of sight causes the peak of emission that is further from the mean velocity of the system to be enhanced, for similar reasons as in the case of the expanding and contracting clouds discussed in §3.

For the case of the slit parallel to the rotational axis, the 2-dimensional spectrum displays a symmetric double-peaked pattern, which essentially results from averaging along the equator the spectrum with different projected velocity. This averaging results in a less sharp central trough and smoothing of the two peaks. In all cases, the velocity dispersion increases the distance between the two peaks in the photon distribution and broadens the width of each peak.

## 5 Summary and Discussion

We have developed a Monte Carlo code of Ly$`\alpha `$ photon scattering, to construct simulated images and spectra. The code is fully three-dimensional and adaptable to a cloud with any given Ly$`\alpha `$ emissivity, neutral hydrogen density distribution, and bulk velocity field. Several simple cases have been presented to show the images and spectra expected from a uniform gas distribution. We have applied the code to model DLA systems illuminated by the intergalactic UV background to study the spatial and frequency distribution of emitted Ly$`\alpha `$ photons. Self-shielding produces a core in the Ly$`\alpha `$ surface brightness profile. For a spherical, static cloud, Ly$`\alpha `$ photons have a double-peaked distribution which is symmetric around the line center frequency. An oblate, rotating cloud shows the same double-peaked distribution with the rotational pattern of increasing mean wavelength with the position on the slit. Observations of the fluorescent emission from damped Ly$`\alpha `$ systems with large telescopes (Hogan & Weymann 1987; Gould & Weinberg 1996) could reveal the presence of such systems and measure a velocity gradient, providing a direct measurement of the rotation rate. The Ly$`\alpha `$ image would also reveal if the emission is due to fluorescence from the external radiation (in which case we expect a large core of the surface brightness at the radius where the hydrogen becomes self-shielded) or due to internal gas dissipation or star formation.

The observations of metal lines in DLA systems show multiple absorption lines (e.g., Prochaska & Wolfe 1997, 1998), suggesting that the gas is clumpy. In this paper, we have assumed instead that the gas has a smooth distribution. The case of enhanced velocity dispersion we have considered is valid only if the gas is in very small clumps that are optically thin to Ly$`\alpha `$ photons. Clumps that can give rise to the observed metal lines are highly optically thick. The effect of a clumpy gas distribution is not difficult to imagine. If there is only one clump along a line of sight, then this clump will essentially produce the same spectral shape in its emission line, shifted to its velocity. Thus, if the emission line could be observed with sufficient angular resolution to resolve the individual clumps, the characteristic spectral shape of Figure 6 would appear at every position, although shifted to the clump velocity. If the clumps are not resolved, the emission line shape is obviously smoothed, just like in the case of an expanding or contracting cloud shown in Figure 2. Our numerical code can be applied in the future to the results of detailed hydrodynamic simulations of gaseous halos in the process of forming galaxies in cosmological realizations to predict other observable signatures that can test the mechanisms by which energy dissipation of gas leads to galaxy formation.

A possible effect that we have not considered in this paper is the absorption of Ly$`\alpha `$ photons by dust. Dust seems to be present in at least some DLAs as indicated by the reddening of the background quasars that have DLAs along the line of sight (Fall, Pei, & McMahon 1989; Pei, Fall, & Bechtold 1991) and the abundance ratios of iron-peak elements which follow a dust depletion pattern (Savage & Sembach 1996; Vladilo 1998). The resonant nature of the scattering of Ly$`\alpha `$ photons by neutral hydrogen atoms increases the probability for them to be absorbed by dust grains. This is probably relatively more important when the Ly$`\alpha `$ emission originates from a star-forming region (e.g., Meier & Terlevich 1981; Charlot & Fall 1993), where dust is likely to be abundant. If Ly$`\alpha `$ photons of DLAs are caused by fluorescent emission as cases we consider in this paper, the absorption by dust would not be very severe due to the low column densities of the ionized outer layer (Charlot & Fall, 1991). Since in the presence of dust the escape of Ly$`\alpha `$ photons depends on many factors, such as the H I column density, the distribution of Ly$`\alpha `$ source, the property of dust, and the topology of the medium (e.g., Ferland & Netzer 1979; Hummer & Kunasz 1980; Neufeld 1990, 1991; Charlot & Fall 1991), detailed modeling of DLAs is necessary for dust absorption to be taken into consideration.

The authors thank the referee Joop Schaye for helpful comments. We also thank Xuelei Chen and David Weinberg for useful discussions. This work was supported in part by NSF grant NSF-0098515.
APPENDIX

To choose the velocity along the direction of the incident photon of the atom responsible for the scattering, we need to generate a random number $`u`$ that has the following distribution (see eq.):

$$f(u)\frac{e^{u^2}}{(xu)^2+a^2},$$
$`(\mathrm{A1})`$
where $`x`$ and $`a`$ are given. We give a brief description on how we generate random numbers with this distribution. Since the distribution has the symmetry under the transformation of $`(x,u)(x,u)`$, we limit the following discussion to $`x>0`$.

We make use of the rejection method (Press et al., 1992) to generate $`u`$. The comparison distribution can be chosen to be $`g(u)[(xu)^2+a^2]^1`$, which can be integrated and inverted analytically. A value of $`u`$ with the distribution $`g(u)`$ is first generated. We keep this value only if a second random number uniformly distributed between 0 and 1 is smaller than $`e^{u^2}`$. In this way, we obtain values of $`u`$ with the distribution $`f(u)`$. In practice, when $`x1`$, the above method is inefficient because a very small fraction of values of $`u`$ are not discarded. To increase the fraction of acceptance, we modify the comparison distribution to be

$$g(u)\{\begin{array}{cc}[(xu)^2+a^2]^1,\hfill & uu_0\hfill \\ e^{u_0^2}[(xu)^2+a^2]^1,\hfill & u>u_0\hfill \end{array}$$
$`(\mathrm{A2})`$
The value of $`u_0`$ can be chosen to minimize the fraction of generated values that will be discarded. The acceptance fractions are then required to be $`e^{u^2}`$ and $`e^{u^2}/e^{u_0^2}`$ in the regions $`uu_0`$ and $`u>u_0`$, respectively.

A first random number $`R`$ uniformly distributed between 0 and 1 determines which region we use by comparing it with $`p`$, where

$$p=\frac{_{\mathrm{}}^{u_0}g(u)𝑑u}{_{\mathrm{}}^+\mathrm{}g(u)𝑑u}=(\theta _0+\frac{\pi }{2})[(1e^{u_0^2})\theta _0+(1+e^{u_0^2})\frac{\pi }{2}]^1,\theta _0=\mathrm{tan}^1\frac{u_0x}{a}.$$
$`(\mathrm{A3})`$
We generate $`u`$ through $`u=a\mathrm{tan}\theta +x`$, where $`\theta `$ is a random number uniformly distributed in $`[\frac{\pi }{2},\theta _0]`$ and $`[\theta _0,\frac{\pi }{2}]`$ for $`Rp`$ and $`R>p`$, respectively. Then another random number uniformly distributed between 0 and 1 determines whether the generated value of $`u`$ is accepted by comparing it with the corresponding fraction of acceptance.


text:
# On the coherence verification of the continuous variable state in Fock space

## Abstract

Almost all novel observable phenomena in quantum optics are related to the quantum coherence. The coherence here is determined by the relative phase inside a state. Unfortunately, so far all the relevant experimental results in quantum optics are insensitive to the phase information of the coherent state. Lack of phase information may cause serious consequences in many problems in quantum optics. For example, an ensemble of two mode squeezed states is a classical ensemble if the phase in each state is totally random; but it is a non-classical ensemble if the phase in each state is fixed. As a timly application, verification of this type of phase information in an ensemble of two mode squeezed states from the conventional laser is crucial to the validity of the continuous variavable quantum teleportation(CVQT) experiment. Here we give a simple scheme to distinguish two different ensemble of states: the Rudolph-Sanders ensemble, by which each squeezed states emitted has a uniform distribution from $`02\pi `$ on the phase value; and the van Enk-Fuchs ensemble, which emmits identical states with a fixed(but unknown) phase for every state. We believe our proposal can help to give a clear picture on whether the existing two mode squeezed states so far are indeed non-classical states which can be used as the entanglement resource.

Quantum coherence of different states plays a fundamentally important role in the whole subject of quantum optics. The squeezed states have been widely used to demonstrate various novel non-classical properties in the past. In particular, a two mode squeezed state

$`|r,\varphi >=e^{re^{i\varphi }(a_1a_2a_1^{}a_2^{})}|00>={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}\lambda ^ke^{ik\varphi }|00,`$ (1)
where $`a_i^{}`$ and $`a_i`$ are bosonic creation and annihilation operators respectively in the two mode Fock space, $`|00>`$ is the vaccum state, is hopefully to be used as the entangled resource to carry out various novel tasks in the creteria of quantum information such as the quantum teleportation. Recently, this state has been used as the enatnglement resource to teleport a quantum state between two spatially separated parties. However, as it is pointed out by Rudolph and Sanders that the the phase information $`\varphi `$ here has never been tested, therefore states as defined by the above equation from certain source could have different phase $`\varphi `$ for each different wavepackets. The lack of the phase information causes the loss of coherence of the state. Although lots of experiments have been done sucessfully to reconstruct the continuous variable states in the past, none of them is related to the phase information. That is to say, a pure state is not the only possible state that is compatible with the observed results. As it is noted that, a conventional meassurement on optical fields, such as homodyne detection using lasers, involving mixing of defferent incoherent fields and subsquent detection by energy absorption in photodetectors: all such meassurements are completely insensitive to any optical coherence. Actually, the state could be in arbitrary classical probabilistic distribution via the different phases. If we use the random uniform distribution, the observed state is actually a mixed diagonal state in Fock space, which is given by

$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{2\pi }}|r,\varphi ><r,\varphi |d\varphi ={\displaystyle \underset{0}{\overset{\mathrm{}}{}}}\lambda ^{2n}|nn><nn|.`$ (2)
Definitely, there is neither quantum coherence nor quantum entanglement in the above state. The state defined by eq(2) is a totally classical state. If the phase $`\varphi `$ of each states from certain source are indeed uniformly distributed from zero to $`2\pi `$, the so called squeezed states will have little novel properties in the practical use because in such a case all quantum coherence has been lost and all the observable phenomena are just the same as that given by classical optics. For example, we can never really take the advantage of its squeezing property for certain quadrature variable( although a certain quadrature operator ie indeed squeezed on a single wavepacket).

Due to the lack of the phase information, the observable quantum coherence property have never been verified on an ensemble of squeezed states. To answer the question whether the conventional source can produce the nontrivil quantum coherent state we must have a way to detect the phase information $`\varphi `$.

A timly application for this type detection is on the validity of the CVQT experiment. Recently, the phase information of the two mode squeezed states has drawn much attention of the physicists due to the issue of quantum teleportation of the continuous variable state in Fock space. Since the tomography result is independent of the phase information $`\varphi `$, it is also possible that the the state used as the entanglement source in the quantum teleportation experiment is a mixture of the states with different $`\varphi `$, as defined by eq.(2). If this is the case, then no quantum state can be teleported by such a separable state. Definitely, an ensemble of $`N`$ identical copies of pure states as defined by $`|r,\varphi ><r,\varphi |^N`$ with the unknown $`\varphi `$(but all states in the ensemble have the same $`\varphi `$) is totally different from the ensemble as defined by eq.(2). Unfortunetely, so far there has been no way to distinguish these two totally different cases.

Due to the unclearity of the phase $`\varphi `$, there is a very hot discussions on the validity of the entanglement resource used in the CVQT experiment. The discussions can be sumarrized as the following:

1. Rudolph and Sanders: The phase $`\varphi `$ in each of the squeezed states from a conventional laser source are uniformly distributed from zero to $`2\pi `$. The correct form of the quantum state for the ensemble is given by eq(2)

2. van Enk and Fuchs: The traditional formalism is insufficient to describe the meachanism two mode squeezed states produced by the conventional laser source. Using quantum de Finetti theorem one can see that the two mode squeezed states produced by a conventional laser is essentially an ensemble of many coppies of identical two mode squeezed states with a fixd $`\varphi `$, though the value of $`\varphi `$ is unknown.

We may see that the validity of the CVQT experiment is now reduced to which of the above statements correctly describe the property of the source which produses the two mode squeezed states in the CVQT experiment. That is to say, to know the validity of CVQT experiment, we have to distinguish the Rudolph-Sanders source and the van Enk-Fuchs source. In this letter we give a scheme to detect the phase information $`\varphi `$ in the two mode squeezed state. That is to say, the meassurement result by our scheme is $`sensitive`$ to the quantum coherence. Using our scheme, the ensemble of states defined by eq(1) or eq.(2) can be easily distinguished. Obviously, the scheme has broad potential applications in the whole subject of quantum optics. For example, given many copies of pure squeezed states, we can verify that they are indeed pure states with the same phase $`\varphi `$. For another example, given two Fuch sources, each source emmits many identical two mode squeezed states, the phase of the state from each source $`\varphi _1`$ and $`\varphi _2`$ respectively, they are fixed and unknown, our scheme can be used detect the value $`\varphi _1\varphi _2`$. As an immediate application, the scheme can be used to judge whether the states used as the entangled resource in a recent quantum teleportation is the pure state as defined in eq.(1) or a mixed state defined by eq.(2). Consequently, this detection could give a a clear judgement on whether the experiment done in ref. is essentially a quantum teleportation or a classical simulation of quantum teleportation.

Now we show how to distinguish a Roudolph-Sanders source and a van Enk-Fuchs source. Lets first see what happens to the van Enk-Fuchs source, i.e. a photon source emmits the identical pure squeezed states, all of them have the same but unknown phase $`\varphi `$ . As it is noted in ref that a two mode squeezed state can be produced with the specifically chosen polarization for each modes as the following

$`|\psi _1(\varphi )>={\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}\lambda ^le^{il\varphi }|l>_{ah}|l>_{bv},`$ (3)
where the subscripts $`a`$, $`b`$ are for the mode $`a`$ and $`b`$ respectively, $`h`$ and $`v`$ represents the horizontal and vertical polarizations respectively. Now we consider another wave packet from the same source, of which the mode $`a`$ and $`b`$ is exchanged. The quantum state for this wavepacket is $`|\psi _2(\varphi )>=_{m=0}^{\mathrm{}}\lambda ^me^{im\varphi }|m>_{av}|m>_{bh}`$. The total state is then given by

$`|\mathrm{\Psi }>=|\psi _1>|\psi _2>={\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}\lambda ^{l+m}e^{i(l+m)\varphi }|l>_{ah}|m>_{av}|l>_{bv}|m>_{bh}.`$ (4)
After taking a meassurement on the photon number of each mode, i.e., on the quantity $`n=l+m`$, the state $`|\mathrm{\Psi }>`$ is collapsed to a specific entangled state. As it has been demonstrated in ref, this type of meassurement cab be carried out by either a quantum non-demolition meassurement or by a more feasible way, the destructive photon counting with post selection. Suppose the meassurement result is $`n`$ for each mode, then the state after the meassurement is $`_{m=0}^ne^{i[(nm)\varphi +m\varphi ]}|nm>_{ah}|m>_{av}|nm>_{bv}|m>_{bh}`$. The phase information is clearly included in the state after the meassurement. For simplicity, we consider only the cases of $`n=1`$ only. In the real experiment we can set an appropriate value of $`\lambda `$ so that we have a significant probability to get the result of $`n=1`$. With such a setting, given many coppies of state $`|\mathrm{\Psi }>`$, we can always have a significant number of states with $`n=l+m=1`$ after the meassurement. In the case of $`n=1`$, the state is

$`|\mathrm{\Psi }_1>=e^{2i\varphi }(|1>_{ah}|0>_{av}|1>_{bv}|0>_{bh})+|0>_{ah}|1>_{av}|0>_{bv}|1>_{bh}.`$ (5)
There is only one photon in each mode for the satte defined above. Obviously, the state defined above can be rewritten in the following way

$`|\mathrm{\Psi }_1>=e^{2i\varphi }(|H>_a|V>_b+|V>_a|H>_b).`$ (6)
Where the states $`|H>`$, $`|V>`$ are for a horizontal or a vertical polarized photon states respectively. For this state, meassurement results for the polarization of the two modes are always $`different`$. Clearly, if the states given from the source are indeed pure states, i.e., they are many coppies of $`|\mathrm{\Psi }_1>`$, then we can rotate the polarizers by the same angle to both modes, and in principle we can find certain angle($`\pi /4`$) by which the meassurement result of the polarization in the two modes are always $`same`$ ! On the other hand, if the given states from the source are mixed states as defined in eq(2), i.e., phase $`\varphi `$ in each wave packet can be different, we will have a clearly different observation results. In such a case, by the same operation, for all cases we get the meassurement of $`l+m=1`$, the state is $`\frac{1}{2}(|HV><HV|+|VH><VH|)`$. The correlation between the two modes will linearly decreased with $`\mathrm{cos}^2\beta `$ when the polarizers are rotated, here $`\beta `$ is the angle that is rotated. Consequently, a van Enk-Fuchs source and a Rudolph-Sanders source can be distinguished in the following way:

Initially, in both cases they are totally negatively correlated, i.e., whenever a photon is detected in mode $`a`$, there most be no photon detected in mode $`b`$, and vice versa. We rotate the polarizers by $`\pi /4`$, the van Enk-Fuchs source(see eq(1)) will give a totally positive correlation for the meassurement result, while the Rudolph-Sanders source(see eq(2)) will give no correlation for the meassurement. Also, in rotating the polarizers, if the correlation does not decrease linearly with $`\mathrm{cos}^2\beta `$, the source must be not a Rudolph-Sanders one.

Obviously, the above scheme can be also used to detect the phase difference for two van Enk-Fuchs sources. Now $`|\psi _1>`$ and $`|\psi _2>`$ are collected from two different sources whose unknown fixed phase are denoted by $`\varphi _1`$ and $`\varphi _2`$ respectively. After the meassurement on the photon number basis of each mode is done, the state up to a global phase factor is

$`|\mathrm{\Psi }_1>=|H>|V>+e^{i(\varphi _2\varphi _1)}|V>|H>.`$ (7)
We rotate the polarizers continuously and test the correlations at each stage. There must be certain angle $`\beta _0`$ at which we can observe that the meassurement results of the two modes are totally positively correlated. This angle $`\beta _0`$ determines the phase difference $`\varphi _2\varphi _1`$.

Note that in a real experiment, we actually do not need to take a meassurement on the basis of of the photon numbers of each mode, i.e. $`_{n=0}^{\mathrm{}}|n><n|`$, where $`n=l+m`$. The meassure on such a basis could be difficult by our current technology. The only subtle task for us is to collect the wave packet $`|\psi _1>`$ and $`|\psi _2>`$ so that they are indistiuguishable. Once we can make $`|\psi _1>`$ $`|\psi _2>`$ indistinguishable we can simply carry out the detection scheme by the usual photon counting and obtain the correct result to a good approximation. Suppose we use port A and port B to detect the photons from mode A and B respectively. A polarizer is placed before each port. The polarizers are placed horizontally in the begining. We choose a very small $`r`$ so that $`r^2<<1`$ for the photon source. In the photon counting, we only record the data of those events where at least one port detects one photon $`and`$ no port detects more than one photon. For simplicity, we name the events satisfying this condition as ”good events”. We will analyse the correlation between the two modes only using the data of the good events. In the whole process, there is only a very small probability that a good event is caused by a state with $`l+m1`$, i.e., the wave packet has been collapsed to a state of which the photon number of each mode is not 1. It’s easy to see that this type of event can only happen with a small ptobability. First, $`l+m=0`$ is impossible, because whenever we observed a good event, we have observed one photon at least in one port, so the photon number of that mode must not be 0, consequently $`l+m0`$. Second, $`l+m=2`$ or a larger number is possible but the chance is negligible. We have already set $`r`$ to be very small. Since $`r^2<<1`$, the probability that the state $`|\mathrm{\Psi }>`$ collapsed to a state with $`l+m=1`$ is much larger than that of a state with $`l+m>1`$. For example, taking $`r=0.01`$, in average, we can have one good event from 10000 events, and the probability that this good event is caused by a stete with $`l+m=1`$ is several thousands times larger than that by a state with $`l+m1`$. Actually, to a good approximation, we even do not have to require the photon detector to distinguish 1 photon case from many photon case in the detection if $`r`$ is small enough. Because once the detector detects photons, the probability of 1 photon is much higher than other cases. For clarity, we insist on using the term ”1 photon” in the rest parts of the letter. However, whenever the term ”one photon” is used, it can be approximately understood as ”any number of photons”.

Initially, all good evnts must be the events on which 1 photon detected in certain mode and 0 photon is detected in another mode. If the source is the van Enk-Fuchs source, when the polarizers are rotated by $`\pi /4`$, a photon in each mode must be detected to all good events. If the source is the Rudolph-Sanders source, when the polarizers are rotated by $`\pi /4`$, to all good events, half of them are $`(1_A,0_B)`$ or $`(0_A,1_B)`$ and half of them are $`(1_A,1_B)`$, where the 0 or 1 represents the number of photons detected, subscripts $`A,B`$ represent the the port $`A,B`$(or the mode $`a,b`$) respectively.

The phase difference can also be detected for two Enk-Fuchs sources. In this case, we have to rotate the polarizers continuously. The phase difference $`\varphi _2\varphi _1`$ is determined by the angle $`\beta _0`$ we observed 1 photon in each mode to all good events. Thus we see, to a good approximation, our scheme can carried out by the normal photon counting technique.

Acknowledgement: We thank Prof Imai for support.