text:
# XMM-NEWTON High Resolution Spectroscopy of NGC 5548

## 1 Introduction

The immediate environment of active galactic nuclei (AGN) is poorly understood. The black hole is fed from an accretion disk, which becomes visible through high energy radiation. Possibly powered by radiation pressure, an outflowing wind is formed from the accretion disk. In some AGN, the wind is visible through its continuum and line absorption of the radiation from the nucleus. These absorption lines were already detected in the UV part of the spectrum (Mathur, Elvis & Wilkes mathur (1995)), and have now been confirmed in the X-ray spectra (Kaastra et al. kaastra00 (2000)).

NGC 5548, a Seyfert 1 galaxy, has been well studied in the X-rays due to its relative X-ray brightness, and the fact that it is a nearby active galactic nucleus with a redshift of z=0.01676 (Crenshaw & Kraemer crenshaw (1999)). This Seyfert galaxy was studied in detail with ROSAT (Done et al. done (1995)). They detected an O vii and O viii K-shell absorption edge, and concluded that there was a warm absorber present. Therefore, NGC 5548 is an ideal Seyfert 1 galaxy to study with the high resolution spectroscopy instruments on-board of XMM-Newton and Chandra. Table 1 list the observations made with these instruments. In this paper we discuss a 137 ks exposure time Reflection Grating Spectrometer (RGS) spectrum from the guest observation program of XMM-Newton.

## 2 Observations and data reduction

Here we present data taken with the RGS instrument (den Herder et al. denherder (2001)). The EPIC data (the continuum spectrum, and the Fe K$`\alpha `$ line) are presented separately (Pounds et al. pounds (2002)). The time variability of the combined XMM-Newton and RXTE data sets will be presented by Markowitz et al. (mark (2002)).

Shortly after launch two of the eighteen CCD chains failed. For RGS 1 this causes a data gap in the $`1014`$ Å range (where most Fe-L and Ne absorption lines are), for RGS 2 the gap is between $`2024`$ Å (where the O vii triplet is).

The data were reduced using the developer’s version of the XMM-Newton data analysis package SAS version 5.3, which comprises a correct calibration of the effective area and the instrumental O-edge. The latter is important as the 1s-2p absorption lines of O ii and O iii are cosmologically redshifted toward the instrumental O i edge.

In this paper we analyze the first and second order spectra obtained from both RGS 1 and 2. Although the second order spectrum has a lower count rate and only a wavelength coverage from 5 to 19 Å, it does have a higher wavelength resolution than the first order. Throughout the paper we have fitted RGS 1 and 2 simultaneously, but for clarity of representation we have added the RGS 1 and RGS 2 spectra in the figures.

All the spectral analysis was done with the SPEX package (Kaastra et al. 2002b ). Fig. 1 shows the fluxed spectrum of NGC 5548, indicating the Unresolved Transition Array (UTA) of inner-shell iron absorption lines and some other prominent features.

### 2.1 Wavelength scale accuracy

In general the RGS wavelength scale can be reconstructed to an accuracy of 8 mÅ (den Herder et al. denherder (2001)). However, in our case larger uncertainties exist, due to errors in the recorded satellite attitude and right ascension in the data files available. From the analysis of the first order data we find that the O vii forbidden line is blueshifted by 285 km s-1 ($``$0.021 Å). The O vii forbidden line is not blueshifted in the earlier LETGS data of NGC 5548 (Kaastra et al. 2002a ). No blueshift for the forbidden line can be determined in second order. However, comparing first and second order results for the warm absorber, we find a smaller outflow velocity in the second order data set. The difference in wavelength between first and second order, for the strongest line in the second order, O viii Ly$`\beta `$ (both orders were fitted with the xabs model, see Sect. 3) is $``$0.019 $`\pm `$ 0.005 Å. Finally, in the first order data set there is a weak Galactic O vii recombination line, which has a rest wavelength of 21.6 Å. This absorption line also shows an apparent velocity shift of $``$230 km s-1 ($``$0.017 Å), although no blueshift should be expected. To correct for this wavelength calibration problem we applied a constant wavelength shift of $``$0.017 Å to the model for the first order spectra. If not corrected, this wavelength shift would result in a velocity error between 130 km s-1 for the longest wavelengths to 620 km s-1 for the shortest wavelengths. The residual error in the wavelength scale should be within the 8 mÅ uncertainty of the calibration. All velocities quoted in this paper are corrected for this wavelength shift.

## 3 Modeling the warm absorber

In NGC 5548 there is an intricate pattern of densely spaced absorption lines as well as the possible presence of weak relativistic emission lines (Kaastra et al. 2002a ). Blending, mostly of iron lines, is important in many wavelength ranges, making equivalent width (EW) estimations rather complicated. On the other hand, in the X-ray band many of the observed ions have more than one detected absorption line, significantly improving the column density and velocity outflow measurements. This is also important in order to disentangle possible partial covering factors and column densities of the ions. The possibility of relativistic emission lines complicates the analysis of the continuum spectral shape, which is already complicated to derive due to the presence of many edges and absorption lines. We consider here two approaches for fitting the warm absorber. In the first approach all absorption lines and edges of each individual ion are calculated, using a separate value for the column density for each ion. The transmission of all ions is then combined to calculate the total transmission of the warm absorber. Outflow velocity and velocity broadening are coupled for the different ions. This is the slab model in SPEX (Kaastra et al. 2002b ).

In the second approach the column densities of the ions are coupled through a photoionization model. In this case only the total column density, ionization parameter, elemental abundances and the outflow and broadening are free parameters. This is the xabs model. For a full description of both methods see Kaastra et al. (2002b ).

The results obtained with the xabs model are more dependent on theory, because the ionization parameter determines the column density ratios of the ions. An advantage of the xabs model is that different ionization components, each with their own velocity structure can be combined. This is harder to accomplish with slab components, in particular if the velocity components are barely separable with the current spectral resolution. This is the case in NGC 5548, and would lead to strongly correlated parameters in the spectral fitting procedure. One further difference between the models is that ions which have a rather small cross-section are automatically included in xabs, but since these are more difficult to discern individually, they are not usually included in the slab modeling.

## 4 Data analysis

### 4.1 The continuum

The continuum model consists of a power-law spectrum with Galactic absorption, and a modified black body (see Kaastra & Barr kaastrab (1989) for more details). A second modified black body with a high temperature is required to describe the EPIC data, but is not detected in the limited RGS band (Pounds et al. pounds (2002)). This indicates that RGS data alone cannot constrain the continuum accurately. Table 2 gives the continuum parameters used for the analysis of the RGS data.

The temperature ($`T`$) of the modified black body was allowed to vary as well as the normalizations of the power-law and the modified black body. However, the photon index ($`\mathrm{\Gamma }`$), which is very sensitive to small changes in the model, as well as the Galactic H i column density (1.65 $`\times `$ 1024 m-2, Nandra et al. nandra (1993)) were kept fixed throughout the fitting process. The obtained parameters are within calibration uncertainties consistent with those obtained from simultaneous EPIC and BeppoSAX observations.

### 4.2 Spectral fit using the slab model

As a first approach to the analysis of the warm absorber, we modeled it by adding a slab component to the continuum model. We refitted the data, leaving the same parameters free as above for the continuum, as well as the overall outflow velocity, and all the relevant ion column densities. Finally, after obtaining this overall fit we also let the photon index vary, and now obtain a reasonable result. The results of this fit were then used to fine-tune the column density and outflow velocity of each ion which has a detectable absorption line or edge in the RGS wavelength range. To do so we decoupled subsequently for each ion its outflow velocity from the outflow velocity of all the other ions in the model. This method was followed to minimize the effects from blending. The significance of the different ionic columns can be derived from the best fit column densities and their error bars. The result of this fitting is summarized in Table 3. This Table also lists the ionization state for the maximum column density of the particular ion. Note however, that each ion exists for a range of ionization states, and thus the quoted ionization is only an indication. Adding the slab component significantly improves the $`\chi _{\nu }^{}{}_{}{}^{2}`$ from 1.96 to 1.26 (with 2404 degrees of freedom).

For the ions with only an upper limit for their column densities, absorption is only detected through continuum absorption (Ar xiv, Fe v and Fe vii) or a single/few weak lines (O iv, Si x, Si xi, Fe i, Fe xii and Fe xiv). As a result the column density as well as the velocity of these ions are poorly determined. The other ions with poorly determined velocity shifts, but with solid column densities, are those for which the column density is determined from an edge (C iv, N v and Si ix), for which the absorption lines are strongly blended (Fe xii), or because only a few weak absorption lines are detected.

We find a very well constrained column density for C iv and N v, although both these ions have in the RGS band only an absorption edge, at 36 Å and 25 Å (rest wavelength) respectively. As the depth of the edge correlates with the assumed underlying continuum, these column densities are rather more uncertain than the fit results indicate. To ascertain how these column densities change the overall continuum shape we refitted the data, but now forcing a column density for C iv and N v of zero. For the new fit, which gives a similarly good $`\chi _{\nu }^{}{}_{}{}^{2}`$, the normalization of the modified black body increased by an order of magnitude, while the temperature of the same component decreased by an order of magnitude. As a result of this change in the continuum model, the column densities found for O iii and Ne ix also changed significantly. In this respect the xabs modeling is more reliable, as the column densities determined for ions that are only detected through an edge, are based upon the determined ionization state.

The same method for modeling the warm absorber was used for fitting the second order data. The second order data for RGS 1 and 2 were not fitted together with the first order data, as the systematic shift in the first order spectra is corrected for in the applied model and not in the data itself. However, both orders are in excellent agreement, see Fig. 2.

The average outflow velocity found, including only those ions with well determined velocities in both orders, is $``$350 $`\pm `$ 250 km s-1 for the first order data set, in comparison with $``$420 $`\pm `$ 340 km s-1 for the second order data set. For further analysis we focus on the first order data, because the better statistics outweigh the better spectral resolution of the second order data.

### 4.3 Spectral fit using the xabs model

Using the same continuum model as in the previous section, we replaced the slab model by the xabs model for the warm absorber. First we fitted only the ionization parameter ($`\xi `$), the hydrogen column density ($`N`$H) and the outflow velocity ($`v`$), in addition to the continuum parameters that we leave free throughout the fitting procedure. The abundances were kept to the solar values (Anders & Grevesse anders (1989)). Later, the abundances were left as free parameters. However, as hydrogen produces no X-ray lines, we normalized the abundance relative to oxygen, which was frozen to solar abundance. Finally we added two more xabs components to the model, both with lower ionization parameters, significantly improving the spectral fit.

Adding the first xabs component improved $`\chi _{\nu }^{}{}_{}{}^{2}`$ from 1.97 to 1.66. The addition of the second xabs component, which was the one with the lowest ionization parameter, further improved $`\chi _{\nu }^{}{}_{}{}^{2}`$ to a value of 1.46. However, for this xabs component we found two abundances which have unphysical values. Namely, Ne and S are overabundant by a factor of 43 and 106, respectively. From detailed fitting, it is clear that both these abundances are mainly determined from the absorption edges in the spectrum, and are thus rather uncertain. Fitting the model again, but fixing the abundance of Ne and S in the lowest ionization component to solar, we find a similar good fit, namely with $`\chi _{\nu }^{}{}_{}{}^{2}`$ of 1.48. The addition of the last xabs component gives a final $`\chi _{\nu }^{}{}_{}{}^{2}`$ = 1.22. Table 4 lists the best fit parameters, with the xabs components labeled to be consistent with the components found by Kaastra et al. (2002a ). The significance, according to F-tests, is at least 100, 99.8 and 99.9 percent for components A, B and C, respectively. Only in the last step we left the photon index free, but this did not improve the $`\chi _{\nu }^{}{}_{}{}^{2}`$. The continuum parameters for this model correspond within 1 $`\sigma `$ with the continuum parameters as found from the slab fit.

### 4.4 Emission and absorption by oxygen

Fig. 3 shows the spectrum around the O vii triplet. In Table 5 we list the parameters of the O vii forbidden line in emission. For the O vii forbidden line Table 6 gives the fluxes, outflow velocity as well as an upper limit to the width of the line, as measured from the HETGS and LETGS spectra (data taken from Kaastra et al. 2002a ) and our data set. Note that this line has the same measured flux over a period of 1.5 years, despite large variations in the continuum. In all three observations the line is unresolved. The above indicate that the O vii forbidden line is formed at least 0.46 pc from the ionizing source.

Interestingly, the inter-combination line is not present in the data. An excess is observed at 21.74 Å rest wavelength, shifted by 0.06 Å from the rest wavelength of the inter-combination line. The relative wavelength accuracy of RGS is 2 mÅ. Assuming that the inter-combination and forbidden lines are formed at the same distance from the ionizing source and under the same conditions, a blueshift of 620 km s-1 is inconsistent with the lack of blueshift for the forbidden line. Because the forbidden line is strong, models predict that even tenuous gasses should produce a weak but detectable inter-combination line. One possible explanation is that the inter-combination line blends with an O vi absorption line at 21.79 Å.

Another emission feature for oxygen is the radiative recombination continuum (RRC) of O vii between 16.9 and 17 Å (see Fig. 4). As there is severe blending with iron absorption lines, we approximated the RRC simply with a delta function. We find a 2 $`\sigma `$ detection for this RRC and the narrowness of the feature (unresolved by the instrument) indicates a low temperature. The flux is 0.06 $`\pm `$ 0.03 ph m-2 s-1 which is consistent with the 0.06 $`\pm `$ 0.10 ph m-2 s-1 value as determined from the LETGS spectrum (Kaastra et al. 2002a ).

The presence of inner-shell absorption lines of oxygen together with the higher ionized O vii and O viii absorption lines allow for an accurate ionization determination, independent of elemental abundances. The ionization parameter determined for oxygen spans three orders of magnitude from log $`\xi =1.8`$ for O iii to log $`\xi =+1.8`$ for O viii.

The deepest absorption line in the spectrum is the O viii Ly$`\alpha `$ line. For O viii we also detect the Ly$`\beta `$, Ly$`\gamma `$ and Ly$`\delta `$ lines and some hints for the Ly$`ϵ`$ line in absorption (see Fig. 2). The multitude and depth of these absorption lines allow thus for an accurate column density and outflow velocity determination (see Table 3). These higher order lines are important for the detection of variability in the warm absorber, as they have an optical depth close to unity. However, we did not detect significant variability in the measured column densities as compared to the previous LETGS observation, which has rather large error bars. Also within this observation no variability in the warm absorber was detected.

Fig. 3 shows that our model overestimates the O vii resonance line. A probable explanation is that the emission component of this line partly refills the absorption line due to the imperfect spectral resolution. In the LETGS data (Kaastra et al. 2002a ) this line has a clear P Cygni profile, and also in the RGS there is excess emission at wavelengths just longer than the resonance line. Both models for the warm absorber don’t include possible P Cygni profiles. The excess emission just long-ward of 22.25 Å is noise, as the width of the feature is narrower than the point spread function of the instrument. The Galactic absorption by O vii, detected with 1.5 $`\sigma `$ significance is also indicated in Fig. 3.

Shortward of the forbidden of O vii there is a deep O vi absorption line, and 0.25 Å longward of the forbidden line there is a deep O v absorption line. Both are inner-shell absorption lines, indicating the importance of these transitions in AGN. For O v two absorption lines are detected. However, the line at 22.33 Å is redshifted by about 700 km s-1 relative to the weaker line at 19.92 Å, and the other oxygen absorption lines. To check the possible uncertainty on the wavelength for the strongest O v line, the O v wavelength and oscillator strength were also calculated using the Cowan code (cowan (1981)) (Raassen, private communication). This gives a rest wavelength of 22.38 Å, closer to the observed rest wavelength of 22.39 Å, and indicates that uncertainties in the rest wavelength can be of order 0.05 Å; i.e. equivalent to the instrumental FWHM.

### 4.5 Absorption by iron

Iron is the best indicator for the ionization structure of the warm absorber, as it has absorption lines in the RGS band from neutral Fe i to the highly ionized Fe xxiv. For NGC 5548 we detect absorption from Fe vi to Fe xxiv (see Table 3). The lowest ionization component is represented by the UTA (Sako et al. sako (2001)) formed from inner-shell transitions of Fe vi through Fe xi. These ions represent the lower ionized iron ions. Fig. 4 gives a detailed view of the spectrum between 16 and 17.5 Å, where the UTA is located. Note that although there are some individual absorption features visible, the dominant effect of the UTA is the depression of the continuum due to an unresolved blend of transitions with small cross-sections. The UTA clearly shows the importance of these inner-shell transitions in determining the continuum correctly.

The moderately ionized component is represented by Fe xiii through Fe xix, where we only find an upper limit for Fe xii and Fe xiv. The highest ionized component is represented by absorption by Fe xix through Fe xxiv. For Fe xix through Fe xxii we detect a multitude of absorption lines between 8 (lowest wavelength that we included in the fit) and 15 Å. For Fe xxiii and Fe xxiv only two absorption lines are detected, resulting in uncertain velocity determinations.

It is thus clear that all three ionization components as fitted with xabs are amply confirmed by iron absorption. The question arises whether a continuous ionization structure would be a better fit to the data. From Fe vi through Fe xxiv we detect all iron ions with the exception of Fe vii, Fe xii and Fe xiv. Between Fe xii and Fe xv all column densities are low. Forcing Fe xii to Fe xvi to have a similar column densities as Fe xi or Fe xvii worsens the fit, increasing $`\chi _{\nu }^{}{}_{}{}^{2}`$ by 0.24. This indicates that for log $`\xi `$ between 1 and 1.4 ($`\xi `$ between 10 and 25) (in 10-9 W m) the column densities are lower than for the other ionization states. However, no clear depression in column densities is seen for the transition between the middle and highly ionized component. A more continuous ionization structure than presented by the xabs model for the warm absorber is not excluded by the data.

### 4.6 Absorption by carbon and nitrogen

For carbon we detect a strong C vi Ly$`\alpha `$ line (see Fig. 5), the Ly$`\beta `$ line (see Fig. 7), Ly$`\gamma `$ and hints of Ly$`\delta `$ lines (see Fig. 6). These higher order lines are important in the correct determination of the column density as well as any possible time variability, as their optical depth is around unity. C vi and N vii are produced mainly in the middle ionization component, while C v and N vi are produced in the middle and lowest ionization component. Finally, C iv and N v are produced only in the lowest ionization component (see Table 3 for an indication of the ionization state). The C vi Ly$`\alpha `$ line is rather deep and at a relatively long rest wavelength of 33.736 Å, resulting in a very good velocity shift determination. The velocity component found from the X-rays is consistent with component 3 from the UV, which has an outflow velocity of $`540`$ km s-1 instead of $`510`$ km s-1. There is a hint in our data that the line has substructure, consistent with the earlier LETGS results.

### 4.7 Relativistically broadened emission lines

After finding a model for the warm absorber that is as complete as possible, we added relativistically broadened emission lines of O viii, N vii, and C vi Ly$`\alpha `$ to the model. Such lines were observed in MCG$``$6-30-15 and Mrk 766 (Branduardi-Raymont et al. brand (2001), Sako et al. sako02 (2002)).

In NGC 5548 Kaastra et al. (2002a ) found evidence for weak broadened emission lines of O viii Ly$`\alpha `$ and N vii Ly$`\alpha `$. However, no indication for a relativistically broadened C vi Ly$`\alpha `$ line was detected. We applied the best fit parameters found from the LETGS model to the RGS data but found no significant improvement for any of the relativistically broadened emission lines. Leaving the parameters for the Laor line profile (Laor laor (1991)) free, we find that a broad Gaussian would better describe the line/excess. However, the improvement in fit is still small for all three emission lines. The results are given in Table 7, together with the best fit values for the LETGS data of NGC 5548 (Kaastra et al. 2002a ), and the RGS data of MCG$``$6-30-15 and Mrk 766 (Branduardi-Raymont et al. brand (2001), Sako et al. sako02 (2002)). Note, that the values found for the EW between our RGS results and the LETGS results are consistent within 3 $`\sigma `$, except for the inner radius. These relativistic lines should be time variable, considering that they are formed in the accretion disk close to the black hole. Possibly, time variability explains the differences between the previous LETGS results and ours. Given their weakness, in the further analysis and in all the plots we did not include these broadened emission lines. Fig. 8 gives the wavelength band where the relativistic O viii Ly$`\alpha `$ line should be situated, and a small excess is seen between 18 and 18.4 Å.

## 5 Discussion

### 5.1 Outflow velocity versus ionization

No evidence is found for a correlation between the outflow velocity and the ionization parameter (see Fig. 9). Rather the outflow velocity, for first order, is consistent with $``$380 km s-1 and $``$350 km s-1 for the slab and xabs model respectively over the three orders of magnitude sampled in ionization scale. This is consistent with the UV data (Crenshaw & Kraemer crenshaw (1999)), and the results for the C vi Ly$`\alpha `$ line in the LETGS data (Kaastra et al. 2002a ). Kaastra et al. (2002a ) found a trend of decreasing outflow velocity for higher ionized iron and oxygen ions, although all outflow velocities were consistent with an outflow velocity of $``$340 km s-1. This result is not reproduced here. The outflow velocities observed in both data sets are consistent with the five UV velocity components, suggesting a similar origin of both absorbers.

### 5.2 Hydrogen column density versus ionization

The derived absorbing column density per ion, as determined from the slab modeling, versus the ionization parameter is plotted in Fig. 10. In the Figure we use the ionization parameter corresponding to the state with maximum column density. In reality, however, the ion is produced over a range of ionization states. Fig. 10 thus assumes a continuous ionization distribution, and is dependent on the spectral energy distribution (SED) in our models with XSTAR. We tested a model with and without a UV bump, and found a maximum displacement in log $`\xi `$ of 0.3 for very highly ionized iron, and a maximum change in the derived absorbing column density of 0.1 on the log scale used. However, as most ions at a particular ionization state have a similar shift, the derived abundances and the derived absorbing column density are rather insensitive to the input SED.

In general, the derived absorbing column density increases for higher ionization parameters. The correlation coefficient, for all data points, excluding upper limits, is 0.71, and the related significance is 0.999997; for iron these values are 0.82 and 0.9982, respectively. The derived absorbing column density from the three component xabs model is very similar to the absorbing column density as derived from the slab model, and a similar increase for higher ionization is observed.

In general the abundances found are consistent with solar. An exception is iron, for which both models have lowly ionized iron overabundant compared to lowly ionized oxygen, while a reverse trend occurs for highly ionized iron and oxygen. As carbon and nitrogen are consistent with oxygen, the effect is probably due to uncertainties for the lowly ionized iron ions. A similar overabundance for iron in the lowest ionization state was noted by Blustin et al. (blustin (2002)) for NGC 3783. In the LETGS observation of NGC 5548 (Kaastra et al. 2002a ) found a similar abundance pattern for iron, throughout the ionization range. The explanations listed by Kaastra et al. (2002a ) for this discrepancy have been tested using this data set. First, a calibration effect is unlikely as RGS and LETGS give consistent results. Secondly, we determined the column density for lowly ionized oxygen ions from our own X-ray spectrum; and thus we do not need to rely on non-simultaneous UV measurements. Spectral variability is thus not a possible explanation for the apparent overabundance of lowly ionized iron. Thirdly, we tested the sensitivity of the model to different values for the line broadening, ensuring that saturation due to narrow line widths cannot be responsible for the overabundance of lowly ionized iron. Reducing the width of the lowly ionized absorption lines, has the effect that iron saturates first, worsening the overabundance already found. The maximum width of the lowly ionized lines is set by the width of the inner-shell oxygen lines, which are unresolved.

Possible explanations for the iron overabundance are an uncertainty in the wavelength of the M-shell iron absorption lines; or if the optical depth, $`\tau `$ 1, an uncertainty in the radiative transfer models used. Another possibility is the omission of certain processes like re-emission and Bowen fluorescence in the radiative transfer models used. Finally, the iron abundance could be non-solar, although the physical process causing such an overabundance of iron in Seyfert galaxies is unknown. Non-solar abundances have been detected in IRAS 13349+2438 (Sako et al. sako (2001)) and NGC 1068 (Kinkhabwala et al. kink (2002)).

### 5.3 Outflow geometry

The geometry and physics of the warm absorber as observed in some AGN is poorly understood. One of the outstanding questions is whether this outflow is spherical, in which case the ionization of the gas is mainly dependent on the distance. Alternatively, the outflow could be in localized streams formed due to instabilities in the accretion disk. The ionization parameter is then mainly dependent on the density. As distance and density cause a degeneracy in the determination of the ionization parameter, the location of the outflowing X-ray and UV warm absorber is uncertain. Constraining the geometry of this outflow will help solve these outstanding problems.

From mass conservation ($`\dot{M}_{\mathrm{loss}}`$ = 4 $`\pi r^2nm_\mathrm{p}v\mathrm{\Omega }`$), the ionization parameter ($`\xi `$ = $`L/nr^2`$), the outflow velocity ($`v`$), and the bolometric luminosity ($`L`$), one can determine the mass loss rate, $`\dot{M}_{\mathrm{loss}}`$, as a function of the opening angle ($`\mathrm{\Omega }`$) of the outflowing wind. This is summarized in eq. (1), where $`m_\mathrm{p}`$ is the proton mass. Assuming that the mass loss rate of the wind cannot be higher than the mass accretion rate onto a black hole and that the system is stationary, one can equate both and obtain an upper limit for the opening angle. Eq. (2) gives the mass accretion rate, $`\dot{M}_{\mathrm{acc}}`$, with a mass conversion efficiency of $`\eta `$.

$$\dot{M}_{\mathrm{loss}}=\frac{4\pi Lm_\mathrm{p}v\mathrm{\Omega }}{\xi }$$
(1)
$$\dot{M}_{\mathrm{acc}}=\frac{L}{c^2\eta }$$
(2)
In the case of a Schwarzschild black hole, as considered here, $`\eta `$ = 0.057. For our average outflow velocity of 380 km s-1 the upper limits of the opening angle vary between 3.8 $`\times `$ 10-7 sr for log $`\xi `$ = $``$1.8 and 0.024 sr for log $`\xi `$ = 3, with log $`\xi `$ in 10-9 W m. For an average ionization state of log $`\xi `$ = 1.7 in 10-9 W m, we find an upper limit to the opening angle of 0.0019 sr. The upper limits to the opening angle indicate that the outflow is mostly in a narrow stream, where the densest part of the stream is the narrowest and lowest ionized. For higher outflow velocities or a higher mass conversion efficiency, i.e. in the case of a Kerr black hole, the stream is even narrower.

From Fig. 10 we determine that the derived absorbing column density scales with the ionization parameter, $`N_\mathrm{H}`$ $``$ $`\xi ^\alpha `$. $`\alpha `$ ranges from 0.25 to 0.5 depending on whether one uses the $`N_\mathrm{H}`$ obtained for the iron or oxygen ions at the lower ionization states. We define $`N_\mathrm{H}nd`$, where $`n`$ is the density and $`d`$ is the thickness of the outflow in the line of sight. Equating $`N_\mathrm{H}`$ to $`\xi ^\alpha `$, where $`\xi `$ = $`L/nr^2`$, one can determine a relation between the density, thickness and distance. If one further assumes that for a fixed distance ($`r`$) from the ionizing source, there is a range in densities, resulting in a range in ionization states, one finds that the density is related to the thickness $`nd^{5/4}`$ or $`nd^{3/2}`$, respectively. From both simple models we conclude that the outflow occurs in streams with a very small opening angle for the dense and lowly ionized gas, and with a larger opening angle for the less dense and higher ionized gas.

### 5.4 Comparison of UV and X-ray detection of O vi

O vi has absorption lines in both the UV ($`\lambda \lambda `$1032, 1038 Å) and X-ray band (22.01 Å, 19.34 Å and 21.79 Å). Our X-ray column density of O vi of 1020.5 m-2 is nearly an order of magnitude larger than the UV column density (1019.62 m-2) of Brotherton et al. (brotherton (2002)). However, the outflow velocities in the UV and the X-rays are comparable. The difference in derived absorbing column density indicates that the UV data are severely saturated, leading to an underestimate in the column density. Arav et al. (arav (2002)) compare the column density of O vi as observed in this data set and an earlier FUSE UV measurement, and conclude that the column densities could be equal if velocity dependent covering is considered. However, simultaneous observations are necessary to exclude the possibility of column density variability (see Arav et al. arav (2002) for more details). Generally, in the X-rays, we find significantly more low ionization gas than previously deduced from UV data. The discovery of O iii to O v absorption lines, in addition to O vi, gives us for the first time an estimate for the column densities of these ions, which are inaccessible in the UV.

## 6 Summary

We have presented here the highest signal to noise high resolution X-ray spectrum of NGC 5548 obtained yet. The spectrum shows a rich structure in narrow and broad spectral features. We detected a very weak O vii RRC, consistent with a low temperature. Inner-shell oxygen lines, together with the higher order absorption lines from O viii Ly$`\alpha `$ and C vi Ly$`\alpha `$ have an optical depth near unity, important in the detection of spectral variability of the warm absorber, although no spectral variability was detected in this data set or from comparison with the earlier LETGS spectrum. Uncertainties in rest wavelengths for the inner-shell ions complicate the study of these absorption lines, certainly for the lowly ionized ions. Possibly this or uncertainties and the omission of certain processes in the radiation transfer models lead to an deviant iron to oxygen abundance ratio.

The detected warm absorber spans three orders of magnitude in ionization, from ions as lowly ionized as those measured in the UV band, to nearly fully ionized iron. Only a few percent of the X-ray warm absorber column density is lowly ionized, represented by oxygen and iron inner-shell transition lines. The outflow velocities determined from the X-rays are all consistent with those measured in higher resolution UV spectra, and no correlation between the ionization state and velocity was detected. The similarity in outflow velocity and the detection of O vi give strong evidence that the UV and X-ray warm absorber are different manifestations of the same outflowing wind phenomenon. However, we detect an order of magnitude more O vi in the X-rays than previous UV measurements. In general, we detect substantially more low ionized gas, than is deduced to be present from UV observations. Probably the UV band underestimates the true column densities due to saturation and velocity dependent covering factor. From simple arguments, we found that the outflow detected is consistent with a narrow stream, with smallest opening angles for the densest and lowest ionized gas.

## ACKNOWLEDGMENTS

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). SRON National Institute for Space Research is supported financially by NWO, the Netherlands Organization for Scientific Research. RE acknowledges support from the NASA XMM-Newton grant NAG5-10032. We thank Ehud Behar (Columbia University), Ton Raassen (SRON) for helpful discussions on the uncertainties of inner-shell absorption line wavelengths.


format: text

text:
# WHY STUDY INTERSTELLAR MATTER VERY CLOSE TO THE SUN?

## INTRODUCTION

Interstellar gas and dust within $``$30 pc of the Sun provide a unique opportunity to evaluate the chemical composition and physical properties of diffuse interstellar clouds in a small region of space. Several motivations for studying local interstellar matter (ISM, LISM) stand out.

1. The ISM sets the boundary conditions of the heliosphere. By analogy, the ISM also sets the boundary conditions for the astrospheres of extrasolar planetary systems. Through the heliosphere-cosmic ray flux (and by analogy the astrosphere-cosmic ray flux) connection, the distribution of stable planetary climates may, in part, depend on the present and historical physical properties of the ISM surrounding each system.
2. Although many details of the chemical composition of the ISM are known, the overall picture is poorly understood because elements are observed only in the gas phase of the ISM, line blending in distant stars confuses the interpretation of data, H+ is unobserved in low column density sightlines, and the interstellar grain composition, gas-dust mixing and chemical homogeneity are unknown.
3. The configuration of the heliosphere depends on the physical properties of the surrounding ISM. The Sun moves through space, and currently the Sun is embedded in an outflow of material from the Scorpius constellation. As the Sun travels through space, the heliosphere is responsive to “galactic weather” patterns. Also, the plane of the ecliptic is tilted by $``$60o with respect to the galactic plane, introducing a north/south asymmetry in the heliosphere because the local interstellar magnetic field appears to be parallel to the galactic plane.
4. Turbulence in the ISM, whether macroscopic or microscopic, is poorly understood. Interstellar absorption lines at the velocity of the cloud surrounding the solar system (also known as the Local Interstellar Cloud, LIC) are not observed towards the nearest star, $`\alpha `$ Cen ($``$45o from the nose of the heliosphere), suggesting either a turbulence-driven velocity discontinuity, or a physically distinct interstellar cloudlet in this direction. The heliosphere responds differently to the case where velocity discontinuities, versus “empty” space, distinguish cloudlets.
5. Lastly, dust and neutral gas from the LIC enter the solar system, where observations of the interaction products directly sample the chemical composition and ionization levels of the parent ISM. The combination of ISM observations towards nearby stars and $`insitu`$ observations of He, H, pickup ions, and anomalous cosmic rays in the heliosphere provides strong constraints on the physical conditions of the surrounding cloud.

These unique qualities of the ISM close to the Sun offer a space laboratory for studying cosmic processes that is not duplicated elsewhere in the Universe, and are elaborated on below.

## BOUNDARY CONDITIONS OF HELIOSPHERE AND ASTROSPHERES

The first spectral observations of the interplanetary Lyman-$`\alpha `$ glow established that neutral gas at interstellar velocities has penetrated the heliosphere (Adams & Frisch 1977). The heliosphere responds to the relative Sun-cloud velocity and the partial pressure of the constituents of the LIC. As the Sun journeys through the Milky Way Galaxy, the heliosphere configuration varies in response to the properties of the surrounding ISM (Frisch 1993, 1997; Zank & Frisch 1999; Mueller et al. 2001, 2002). Using recent Hipparcos results for the solar apex motion (Dehnen & Binney 1998), the Sun moves through the local standard of rest at $``$14 parsecs per million years, and has a relative velocity with respect to the LIC of $``$26 parsecs per million years.

The Sun is currently immersed in an outflow of ISM from the region of the Scorpius-Centaurus Association (see below), and the distance to the outflow surface in the heliosphere tail direction is $``$0.9 pc. Hence the Sun is likely to have entered this outflow within the past 105 years, and possibly within the past $`10^4`$ years depending on cloud geometry (Frisch 1994). Prior to that the Sun was immersed in the Local Bubble interior, a region of low density hot plasma ($``$0.005 cm-3, 106 K) where the heliosphere was probably large (heliopause radius in nose direction, $`R_{\mathrm{HP}}`$$``$300 AU, Frisch et al. 2003, FMZL03). An encounter with a diffuse cloud of somewhat higher density, $`n`$(Ho)$``$10 cm-3, would contract the heliosphere ($`R_{\mathrm{HP}}`$$``$15 AU) with dramatic consequences for galactic cosmic ray modulation, and cosmic ray fluxes at the Earth’s surface. In contrast, todays heliosphere has $`R_{\mathrm{HP}}`$$``$100 AU (see other papers in this issue). Fig. 1 illustrates the motion of the Sun through the Local Bubble and with respect to the Cluster of Local Interstellar Clouds (CLIC).

The influence of the galactic cosmic ray flux on the global electrical circuit (Roble 1991,Tinsley 2000) provides a tentative link between heliosphere dimensions and the 1 AU cosmic ray flux. A positive correlation between low altitude tropospheric cloud cover ($``$3 km altitude) and neutron monitor counts (a direct tracer of $``$0.5 GeV cosmic ray fluxes at Earth, FMZL03) has been established over one solar cycle ($``$1985–1995, Marsh and Svensmark 2000). Evidently cosmic rays interacting with the Earth’s atmosphere produce aerosols which seed water vapor condensation and cloud formation. During the period of time the Sun was immersed in the Local Bubble interior (probably for many millions of years), the heliosphere should have been large with maximum galactic cosmic ray modulation, and in the absence of surrounding ISM neutrals anomalous cosmic rays would have been minimal. If further investigation confirms the correlation between low altitude cloud cover and the 1 AU cosmic ray flux, then it would appear that heliospheric modulation of the galactic and anomalous cosmic ray fluxes has a direct influence on terrestrial climate, and a large heliosphere may be a factor in climate stability.

By analogy, the ISM also sets the boundary conditions for the astrospheres of extrasolar planetary systems (Frisch 1993), although variations in stellar winds must also be considered (Wood et al. 2001). Through the heliosphere-cosmic ray (and by analogy astrosphere-cosmic ray) connection, the physical properties of the present and historical ISM surrounding each system may be a factor in the distribution of stable planetary climates. Although Earth-mass planets have not yet been discovered, over 10% of the more than 100 planets now discovered have semimajor axes near $``$1 AU (for a list of discovered planets, see URL http://exoplanets.org). Spectral types for the central star are typically F, G, or K, and include both giant and main sequence objects. Mass-loss characteristics of cool stars vary with stellar activity levels, which are traced by X-ray surface flux, and decrease with stellar age (Wood et al. 2002). The astrospheres of extrasolar planetary systems may thus vary substantially from the heliosphere, both because of the surrounding ISM and stellar activity levels. These considerations will be important in evaluating the yet-to-be-discovered climates of extrasolar planets, as well as for understanding the climatic history of the Earth.

## CHEMICAL COMPOSITION, PHYSICAL PROPERTIES, AND GAS-TO-DUST MASS RATIO OF SURROUNDING CLOUD

The low opacity of the CLIC to the cloud surface in the heliosphere tail direction yields relatively high average ionizations of Ho and Heo at the solar location. Observations of the ISM towards more distant stars yields sight-line averaged abundances which may be misleading in the presence of ionization gradients. H+ is generally unobserved in low column-density sightlines. Both Ho and H+ must be included for determining the chemical composition of the LIC, since most of the observed ions are present in both neutral and ionized regions (i.e. He, Ne, Ar, and elements with first ionization potentials less than 13.6 eV). The correct gas-phase abundances are required to evaluate gas-dust mixing and the chemical homogeneity of the ISM.

The LIC is unique in that radiative transfer models are constrained by both observations of the ISM towards nearby stars, and by $`insitu`$ observations of ISM interaction products inside of the solar system (Slavin & Frisch 2002, SF02). The ISM composition and ionization at the solar location has been studied using the radiative characteristics of the material, with the models constrained by observations of ISM towards $`ϵ`$ CMa (using data of Gry & Jenkins 2001) and by pickup ion data for an assumed $`n`$(Heo)=0.016$`\pm `$0.002 cm-3 in the LIC. Within the past year new data offers an improved basis for selecting among the 25 models generated by the SF02 study. Studies of Oo/Ho in the generic ISM now shows that the ratio is relatively constant and $``$410 (Andre et al. 2002). Voyager observations of anomalous cosmic rays yield data on interstellar argon (Cummings & Stone 2002) which is partially ionized in local ISM (Ar0/Ar$``$0.20–0.24). More extensive data on pickup ions are also available (Gloeckler & Geiss 2002), as well as better $`insitu`$ observations of Heo (Witte et al. 1996, and private communication). Two models provide the best match to the range of the available data for an assumed $`n`$(Heo)=0.016$`\pm `$0.002 cm-3. These two best models are Models 2 and 18, and the parameters for these models are shown in Table 1. These models, based on the assumption that O/H$``$400 PPM in d$`<`$3 pc ISM (applying results of Andre et Andre et al. 2002 to the LISM), show that N and S are present at $``$60% of solar abundance while C is mildly subsolar (SF02). For this comparison, solar abundances from Holweger (2001) are assumed. Alternatively, assuming O and C solar abundances from Allende Prieto et al. (2001) would imply unrealistic supersolar C abundances.

The filtration of neutrals in the heliosheath regions is a key unknown quantity. For example, filtration factors for Ar range from $``$0.64 (Cummings & Stone 2002) to $``$0.99 (Mueller and Zank 2002). Hydrogen filtration factors are predicted on a purely theoretical basis (e.g. Ripken & Fahr 1983, Izmodenov et al. 1999), however the radiative transfer models also yield an estimate of hydrogen filtration in the heliosheath regions. Gloeckler & Geiss (2002) derive the filtration factor for Ho by comparing H/He measured in the pickup ion population with Ho/Heo for the nearby ISM. Using their approach, Frisch and Slavin (2002) have shown that Models 2 and 18 predict H filtration factors ($``$0.43) which yield agreement between pickup ion $`n`$(Ho) values and the model predictions for the solar location. This is illustrated in Fig. 2.

The gas-to-dust mass ratio ($`R_{\mathrm{g}/\mathrm{d}}`$) may be calculated from the radiative transfer models, and the results suggest that interstellar gas and dust are not well mixed in the LIC (Frisch & Slavin 2002). When missing mass arguments and an assumed solar reference abundance are applied to Model 18, a gas-to-dust mass ratio of $``$180 is found; subsolar ISM reference abundances yield values a factor of $``$3 larger. However, $`insitu`$ dust observations of interstellar dust by Ulysses, Galileo, and Cassini yield a gas-to-dust mass ratio of $``$125$`\pm `$18 (Frisch et al. 1999; Landgraf et al. 2000; Frisch & Slavin 2002). One explanation for the difference between these two values for $`R_{\mathrm{g}/\mathrm{d}}`$ would be poor gas-dust mixing as a superbubble shell fragment sweeps up both gas and dust in the ISM. Most of the mass of the interstellar dust grains detected within the solar system is carried by grains with masses $`>10^{13}`$ g (or radius $`>`$0.2 $`\mu `$m). The gyro-radius for a grain with mass $`10^{13}`$ g (radius $``$0.2 $`\mu `$m) in a weak magnetic field ($``$3 $`\mu `$G) is $``$0.3 pc (Gruen & Landgraf 2000), which is greater than the distance to the LIC edge in the upstream direction. Thus massive dust grains decouple from gas in the upstream LIC, and the $`R_{\mathrm{g}/\mathrm{d}}`$ over the small scales sampled by the in situ data differs from the value averaged over a typical long interstellar sightlines. The result is that $`R_{\mathrm{g}/\mathrm{d}}`$ may vary over sub-parsec length scales. Other explanations, such as poorly understood ISM abundances, are also possible.

## SOLAR MOTION AND ORIENTATION WITH RESPECT TO THE INTERSTELLAR MAGNETIC FIELD

The Sun moves through space and has recently emerged from the interior of the Local Bubble and entered an outflow of ISM from the Scorpius-Centaurus Association (Fig. 1). Hipparcos results yield a solar motion in the local standard of rest of $``$13.4 km s-1 towards the galactic coordinates $`l`$=28o, $`b`$=+32o (Dehnen & Binney 1998).

The north ecliptic pole is directed towards $`l`$=96o, $`b`$=30o, giving an ecliptic plane tilted by $``$60o with respect to the galactic plane. Observations of the interstellar magnetic field directly outside of the heliosphere have been elusive, but evidence for a weak magnetic field ($`<`$3 $`\mu `$G) is given by starlight polarization caused by interstellar dust grains aligned by nearby magnetic fields (Tinbergen 1982, Frisch 1990). The patch of ISM causing this polarization is located towards $`l`$=0$`\pm `$20o, $`b`$=–20$`\pm `$25o, and is close to the ecliptic plane towards the nose direction of the heliosphere (see Fig. 3). The polarization data show a nearby magnetic field orientation which is approximately parallel to the galactic plane and directed towards $`l`$70o. The classical interstellar dust grains which cause the polarization are charged, and pile up in the heliosheath regions as they are deflected around the heliosphere (Frisch et al. 1999). Although the region of maximum polarization closely follows the ecliptic plane, it also coincides with the CLIC LSR upstream direction (Fig. 3). Observations by Voyager 1 and 2 observations of $``$3 kHz emission in the outer heliosphere also indicate the interstellar magnetic field direction is parallel to the galactic plane (Kurth & Gurnett 2003). The combined ecliptic tilt and interstellar magnetic field direction produce a north/south asymmetry in the heliosphere configuration which is visualized in Fig. 4.

## KINEMATICS AND TURBULENCE OF NEARBY ISM

The LIC is one of a group of cloudlets which have a bulk flow motion in the Local Standard of Rest (LSR) of –17$`\pm `$5 km s-1, from the upstream direction $`l`$2.3o, $`b`$–5.2o (Frisch et al. 2002). This group of cloudlets has been referred to as the Cluster of Local Interstellar Clouds, or CLIC. The dispersion of the cloudlets about the bulk flow represents turbulence in the nearby ISM. Models 2 and 18 indicate that $`n`$(Ho)=0.208–0.242 cm-3 in the LIC. The hydrogen column density towards the nearest star, $`\alpha `$ Cen (1.3 pc), has been reported as $`N`$(Ho)=0.41–1.04 1018 cm-2 (Linsky & Wood 1996), yielding a filling factor of $`f`$=0.5–1 for ISM in this sightline. Two cloudlets near the Sun have velocities consistent with the absorption component towards $`\alpha `$ Cen, the “G” cloud (Lallement & Bertin 1992) and the Apex Cloud seen towards Aquila and Ophiuchus (Frisch 2003). The LSR upstream direction velocity vectors for the G and Apex Clouds are, respectively, ($`V,l,b`$)=(–17.7 km s-1,350.5o,9.2o) and (–23.4 km s-1,4.8o,2.1o). For comparison the LSR upstream LIC motion is (–15.8 km s-1,346.1o,0.2o). The Sun is located within a few parsecs of the downstream edge of the bulk flow. The smaller velocity of the LIC, located at the leading edge of the CLIC bulk flow, suggests flow deceleration, while the dispersion of the CLIC velocities reflects turbulence. These flow parameters are derived from observations of 96 interstellar absorption line components (mainly Ca+ data) towards 60 stars which sample the ISM within 30 parsecs of the Sun (Frisch, Grodnicki & Welty 2002). The LSR upstream direction of the CLIC bulk flow is within $``$20o of the observed tangent point of Loop I, which has been formed by successive epochs of star formation in the Scorpius-Centaurus Association over the past 107 years. The CLIC gas has been attributed to a superbubble shell fragment associated with this star formation, with possible implications for the chemical abundance pattern for this nearest gas. For such an origin, the CLIC would be composed of swept up material consisting of ejecta from previous epochs of star formation in the SCA, the ambient ISM, and recently ejected nucleosynthetic products. This scenario provides a hybrid method for the dispersal of nucleosynthetic products – e.g. cloud collisions and turbulence within the outflow.

## CONCLUSIONS

Why study nearby ISM? Taking an anthropomorphic viewpoint, the ISM dominates the interplanetary medium throughout most of the heliosphere. Heliosphere models show that heliosphere properties, and by analogy exoplanetary system astrospheres, are dramatically affected by the physical properties of the surrounding interstellar cloud. From an astronomical viewpoint, the d$`<`$30 pc ISM is the best region for studying the physical processes which operate in the ISM because of the unusual range of data available for constraining radiative transfer models. These data include the ISM properties within the solar system, which trace the neutral component. Accurate ISM models permit basic axioms underlying ISM studies to be tested, such as the assumption that both gas and dust must be included when evaluating the chemical composition of the ISM.

## ACKNOWLEDGEMENTS

The author would like to gratefully acknowledge research support from NASA grants NAG5-6405, and NAG5-11005, and NAG5-8163.

Manuscript received $`3December2002`$, revised $`13February2003`$, accepted $`13February2003`$.


text:
# Influence of Spatial Correlations on the Lasing Threshold of Random Lasers

## Abstract

The lasing threshold of a random laser is computed numerically from a generic model. It is shown that spatial correlations of the disorder in the medium (i. e., dielectric constant) lead to an increase of the decay rates of the eigenmodes and of the lasing threshold. This is in conflict with predictions that such correlations should lower the threshold. While all results are derived for photonic systems, the computed decay rate distributions also apply to electronic systems.

In the theory of disordered media two important regimes, diffusive and localised, are distinguished Beenakker (1997). In the diffusive regime (for weak or moderate disorder), eigenstates are extended and efficient transport is possible. In the localised regime (for strong disorder), eigenstates become localised and transport is strongly inhibited. Many experimental findings for random lasers are more consistent with the assumption of a lasing mode that is localised while a direct experimental analysis of the sample shows that it is in the diffusive regime.

To determine whether a sample is in the localised or in the diffusive regime, a transport property is measured. The most efficient way to achieve this is to check for the rounding of the backscattering cone Schuurmans et al. (1999). Such a rounding is not reported from experiments Cao et al. (1999a); van Soest et al. (2001). Transport is, however, dominated by extend eigenstates, and the simultaneous existence of a few localised eigenstates in a sample in the diffusive regime, i. e., that is *on average* diffusive, would not be noticed Burin et al. (2001). Such localised modes have recently been detected experimentally in a diffusive sample Cao et al. (2002).

The important question is to explain under which conditions such localised eigenstates can exist in a diffusive sample. \[These states have been termed anomalously localised states (ALS) or prelocalised states. For a recent review see Ref. Mirlin, 2000.\] One-dimensional disordered systems are always in the localised regime, i. e., they can never show diffusive behaviour. Theoretical studies on such systems thus cannot give information on the interplay between extended and localised modes. The situation is different in two- and three-dimensional samples. Two-dimensional samples shorter than the localisation length behave similarly to three-dimensional samples, and one is allowed to replace three-dimensional systems by their computationally cheaper two-dimensional counterparts.

The computational cost of treating a two-dimensional sample is significantly higher than for a one-dimensional sample, and only few studies have been published footnote1 . Reference Burin et al., 2001 models the scatterers in the disordered media as dipoles, Ref. Sebbah and Vanneste, 2002 studies circular particles using the finite-difference time-domain (FDTD) method. Both publications do not state explicitly whether their samples are in the diffusive or in the localised regime but the parameters given strongly suggest that the samples are in the localised regime.

The only publication so far on the interplay between diffusive and localised eigenstates inside a diffusive sample seems to be Ref. Apalkov et al., 2002. There it was estimated that localised states become exponentially more frequent when the disorder inside the sample is spatially correlated. Motivated by the picture of photons travelling in closed loop inside a ring-shape structure Cao et al. (1999a), they study a ring-shaped area of higher dielectric constant. This is a very special situation, and it is not obvious how characteristic such a special situation is for the entire behaviour. (It should be noted that the opposite effect, namely that in a localised sample a few modes become extended when spatial correlations in the disorder are introduced, is well-understood, see e. g. Ref. Izrailev and Krokhin, 1999.)

In this paper we will study this problem from a more generic approach. The lasing threshold of a sample is determined by the decay rates of the eigenstates of the system since the loss (=decay) of photons in the mode has to be compensated by pumping if the sample is to start lasing action. Following the approach of Ref. Patra, 2003 we numerically compute the decay rate distribution of a two-dimensional sample on a suitable grid. (Earlier work on the lasing threshold of chaotic cavities Patra et al. (2000) cannot be applied since by construction all eigenstates are extended Beenakker (1997).) We improve on previous work by including spatial correlations.

We use the Anderson Hamiltonian which describes the motion of an uncharged particle in a spatially varying potential. The results can directly be applied also to photonic systems since the Helmholtz equation with a spatially varying dielectric constant has the same form. The sample is discretised with lattice spacing $`\mathrm{\Delta }`$, where for electronic systems $`\mathrm{\Delta }=\pi /k_\mathrm{F}`$ ($`k_\mathrm{F}`$ is the wave vector at the Fermi level) and for photonic systems $`\mathrm{\Delta }=2\lambda /\pi `$ ($`\lambda `$ is the wave length of the light). This is a natural choice in that there is then one propagating mode per transversal lattice point, and the width of the sample is best measured in terms of the number $`N`$ of propagating modes.

Transport is modelled by nearest-neighbour hopping with rate $`1`$. (The results are easily adapted to arbitrary speed $`c`$ of transport.) With a spatially varying potential $`P(x,y)`$ the Hamiltonian for a sample of length $`L=\stackrel{~}{L}\mathrm{\Delta }`$ becomes Patra (2003)

$$\begin{array}{c}_{(x,y),(x^{},y^{})}=\delta _{xx^{}}\delta _{yy^{}}\left[P(x,y)i(\delta _{1y}+\delta _{Ly})\right]\hfill \\ \hfill +\delta _{yy^{}}\left(\delta _{x+1,x^{}}+\delta _{x1,x^{}}\right)+\delta _{x,x^{}}\left(\delta _{y+1,y^{}}+\delta _{y1,y^{}}\right),\end{array}$$
(1)
with $`x=1,\mathrm{},\stackrel{~}{L}`$ and $`y=1,\mathrm{},N`$. The imaginary part of $``$ models coupling of the sample to the outside where we assume that we operate at the centre of the conduction band.

The spatial correlations are assumed to fall of exponentially such that $`P(x,y)`$ takes on random values, normal-distributed with zero mean and correlator

$$P(\stackrel{}{r})P(\stackrel{}{r}^{})=w^2\mathrm{exp}\left(\frac{|\stackrel{}{r}\stackrel{}{r}^{}|}{R_\mathrm{c}}\right).$$
(2)
$`w`$ measures the strength of the disorder, and $`R_\mathrm{c}`$ is the correlation radius. Since we need to generate a large number of mutually correlated random numbers, a Fourier based method has to be employed Kozintsev (1999).

The eigenvalues of the matrix $``$ correspond to the (quasi-)eigenmodes of the system. Their real part $`\omega `$ gives the energy (or, for photonic systems, the frequency) of the mode, and their imaginary part $`\gamma `$ the decay rate footnote2 . We thus have an eigenvalue problem of a non-Hermitian complex symmetric matrix but an eigensolver specifically adopted to this structure exists Patra (2003). Even with this efficient eigensolver, this is still a numerically expensive task, and it is impossible to analyse so many samples that there would be no more noise in the results.

While the model is described in terms of the disorder strength $`w`$, contact with experiments or analytical theories is best made by introduction of the mean-free path $`l`$. It can be computed from the length-dependence of the transmission probability $`T`$ through the sample. In the diffusive regime, $`lLNl`$, it is given by Beenakker (1997)

$$\frac{1}{T}=1+\frac{L}{l}.$$
(3)
The transmission probability has been computed using the method of recursive Green’s functions Baranger et al. (1991) for variable disorder strength $`w`$ and correlation length $`R_\mathrm{c}`$. We determined the mean-free path by fitting the numerically computed $`T(L)`$ to this functional form self-consistently in the interval $`[l;10l]`$. (Picking some other interval, e. g. $`[0;l]`$, changed the result only by about $`1\%`$.) The rescaled results are depicted in Fig. 1 for both $`N=51`$ and $`N=81`$, i. e., for samples of different width. As the figure shows, both sets of curves are almost identical, thereby demonstrating that we are operating in the wide-sample regime. The mean-free path increases significantly as $`R_\mathrm{c}`$ increases. This is immediately obvious since with increasing $`R_\mathrm{c}`$ the potential changes less within a given distance; hence, there is less scattering.

We would like to point out two “curiosities”. The numerical data suggest that the mean-free path $`l`$ factorises as $`l(w,R_\mathrm{c})=f_1(w)f_2(R_\mathrm{c})`$. We did not manage to find an explanation for this observation. Furthermore, the mean-free path seems to scale as $`lw^{1.71}`$ where $`1.71`$ is a numerical parameter. For uncorrelated random order that is uniformly distributed in the interval $`[w;w]`$, a scaling $`lw^{1.5}`$ was found numerically Patra (2003). An analytical theory is available only for one-dimensional systems in the limit $`w0`$ where $`l1/w^2`$ is found Kramer and MacKinnon (1993), so that a universal scaling for finite $`w`$ might not exist at all.

The increase of $`l`$ with increasing $`R_\mathrm{c}`$ poses a problem for a systematic study of the effects of correlations: One has to decide whether to compare samples with identical $`l`$ (and thus variable $`w`$) or samples with identical $`w`$ (and thus variable $`l`$). The final results must depend (apart from trivial prefactors) only on the ratios $`L/l`$ and $`R_\mathrm{c}/l`$ — not on any of those quantities separately. This decision is thus “only” one of numerical efficiency and minimisation of finite-size effects.

For most of our simulations, we have decided to keep $`l`$ constant at $`l=12.5\mathrm{\Delta }`$. For each value of $`R_\mathrm{c}`$, the needed value for $`w`$ was determined by interpolation of the numerical data presented in Fig. 1. The choice of constant $`l`$ offers the advantage that, even if $`R_\mathrm{c}`$ is changed, samples with identical “physical” length $`L/l`$ occupy the same number of lattice points, and thus need the same amount of numerical work. (For constant “physical” length $`L/l`$, the needed computing time scales as $`𝒪(l^2)`$. With constant $`w`$, this would impose severe restrictions on the range of $`R_\mathrm{c}`$ that could be treated.)

We have computed the decay rates for samples of width $`N=50`$ for length $`L/l=1,2,3,4,5,6,9,12,15,18`$, and correlation radius $`R_\mathrm{c}/\mathrm{\Delta }=0.0,0.2,0.5,1.0,1.5,\mathrm{},7.5`$. For each set of parameters, approximately 2000 samples were generated. The maximum value of $`L`$ is limited because we are interested in the diffusive regime, hence $`L`$ has to be sufficiently smaller than the localisation length $`L_{\mathrm{loc}}=(N+1)l/2`$. We did not consider larger values of $`R_\mathrm{c}`$ than $`7.5\mathrm{\Delta }`$ since the sample should be much wider than the characteristic length scale of the disorder. Otherwise the sample would effectively become one-dimensional.

To check the results, we have computed the decay rate distribution also for $`N=80`$ for a few selected values of $`L/l`$ and $`R_\mathrm{c}`$. To complement the other simulations, we have kept $`w`$ constant. As explained above, this implies that we could only include $`R_\mathrm{c}2\mathrm{\Delta }`$.

Following the approach introduced in Ref. Patra, 2003 for samples in the diffusive regime, we fit the numerically computed decay rate distribution to the functional form

$$𝒫(\gamma )=\frac{\gamma _0^2}{\gamma ^2}\left[1Q(M+1,\frac{M\gamma }{\gamma _0})\right],$$
(4)
where the fitting parameters $`M`$ and $`\gamma _0`$ depend on $`N`$, $`L`$ and $`R_\mathrm{c}`$, and $`Q(a,x)\mathrm{\Gamma }(a,x)/\mathrm{\Gamma }(a)`$ is the regularised Gamma function .

All numerically computed histograms fit well to the form (4). The dependence of $`P(\gamma )`$ onto $`M`$ is only weak for $`M1`$, making a precise determination of $`M`$ difficult. Within this error limit, we did not find a significant dependence of $`M`$ on $`R_\mathrm{c}`$, and M is approximately given by the $`R_\mathrm{c}=0`$ result $`M=N/[1+L/(6l)]`$ Patra (2003).

The fitting parameter $`\gamma _0`$, marking the typical value of the decay rates, can be determined to much better precision. $`\gamma _0`$ is much more important for the lasing threshold than $`M`$, so the limited precision of $`M`$ does not pose a problem. The determined $`\gamma _0`$ is shown in Fig. 2.

For $`R_\mathrm{c}=0`$, $`\gamma _0L^2`$ seems to approach a constant value as $`L`$ is increased. This value is about $`20\%`$ larger than the value $`1/(2cl)`$ found numerically for equi-distributed disorder in the interval $`[w;w]`$ Patra (2003). Trying to approach the limit $`L\mathrm{}`$ numerically is not possible since then the sample would become localised.

The important conclusion from Fig. 2 is that for samples of arbitrary length, the introduction of correlations in the disorder leads to an increase of the decay rates. This increase is quick as $`R_\mathrm{c}`$ is increased starting from $`0`$, and becomes slower for large $`R_\mathrm{c}`$. The same behaviour is seen in the control simulations with $`N=80`$ and fixed $`w`$ (and thus variable $`l`$).

Until now, all results are valid for both electronic and photonic systems. Now we will specialise to random lasers. The light inside a random laser is amplified by a laser dye. This dye is able to amplify light within a certain range of frequencies, so only $`K1`$ eigenmodes out of all eigenmodes of the system are amplified. This number varies only slightly between different realisations of the same ensemble due to an effect known as spectral rigidity Beenakker (1997). The lasing threshold is given by the smallest decay rate out of the $`K`$ modes within the amplification window Misirpashaev and Beenakker (1998). This is immediately obvious since the lasing threshold is passed when photons are created faster than they can decay (=escape from the sample).

There are two different approaches to computing the lasing threshold of a random laser. The direct approach is to compute the eigenmodes of a certain number of realisations of the disordered systems, then for each realisation to determine the smallest decay rate inside the amplification window, and finally collect statistics for those values. Since this process yields only a single datum per sample, a very large number of realisations needs to be computed to arrive at data of sufficient quality. The average lasing threshold determined in this way is depicted in Fig. 3 as dashed line.

Frequently more efficient is the second approach where one starts with the computation of the distribution $`P(\gamma )`$ of the individual decay rates. The intermediary result is either a numerical histogram, or, by fitting the histogram to an analytical form, a distribution function that can be evaluated directly for arbitrary argument. We adopt the latter and use Eq. (4)) together with the values of $`M`$ and $`\gamma _0`$ computed by fitting.

The distribution $`P_\mathrm{l}(\gamma _\mathrm{l})`$ of the lasing threshold is the distribution of the smallest value out of $`K`$ values, each distributed according to $`P(\gamma )`$. This assumes that the decay rates of different modes are uncorrelated. For $`K1`$ this seems logical but to our best knowledge no explicit check of this assumption has been published so far. As a side-effect of our computations, we will fill this gap.

$`P_\mathrm{l}(\gamma _\mathrm{l})`$ is difficult to evaluate numerically for $`K1`$ since it is sharply peaked. The position $`\gamma _\mathrm{m}`$ of the maximum of $`P_\mathrm{l}`$ is immediately seen to be given by

$$0=\frac{dP(\gamma _\mathrm{m})}{d\gamma _\mathrm{m}}\left[1_0^{\gamma _\mathrm{m}}P(\gamma ^{})𝑑\gamma ^{}\right](K1)[P(\gamma _\mathrm{m})]^2.$$
(5)
Since $`P_\mathrm{l}`$ is that sharply peaked, $`\gamma _\mathrm{m}`$ already contains all the relevant information, and nothing relevant is lost by not computing the entire distribution. The lasing threshold computed from Eq. (5), after inserting the fitting parameters $`M`$ and $`\gamma _0`$ computed from the numerical histograms into Eq. (4), is shown as solid line in Fig. 3.

From the figure, two important conclusions can be drawn. First, the lasing threshold computed via the two separate methods agrees well. (The “noise” of the two sets of curves is correlated since the same raw data was used as input for both methods.) This means that the decay rates of different modes indeed are uncorrelated. Furthermore, also fitting the numerical histogram to the form (4) is a valid procedure.

The second conclusion — the heart of this paper — is that introducing spatial correlations into the disorder of a random laser increases the lasing threshold, in contradiction to predictions Apalkov et al. (2002).

In this paper, we have thus arrived at two related — but not identical — results. We have shown that the decay rates increase if spatial correlations of the disorder are introduced. The computed decay rate distribution possesses the same form, just with different parameter, as earlier observed for diffusive samples with uncorrelated disorder Patra (2003). This first result means that the “typical” eigenstates become more lossy.

Our second result is that also the lasing threshold increases. This means that also the “special” eigenstates with lower-than-average loss, which are selected by mode competition to become the lasing modes, become more lossy. Even though we did not directly compute the spatial extend of the eigenstates, this still clearly demonstrates that no (pre-)localised eigenstates are formed by the introduction of spatial disorder. We thus fail to observe the prediction that such states should be created Apalkov et al. (2002).

There are several explanations for the difference between our results and Ref. Apalkov et al., 2002. One explanation is that a single ring-shaped area of increased dielectric constant does lead to the formation of a localised state, as suggested by the authors, but the influence of the disorder around that ring-shaped area significantly reduces this effect. Another, equally likely, explanation is that in our simulations we are only able to treat samples of finite size, with a finite number of eigenstates. The creation of a localised state may be an event that is so rare the we fail to see such an event occur in our finite-size simulations. On the other hand, the typical length scale is given by the area per lasing mode, measured experimentally to be a few $`10\mu \mathrm{m}^2`$ in two-dimensional ZnO films Ling et al. (2001), and our samples are larger than this.

To give a definite answer on whether spatial correlations can explain the formation of localised states, more numerical studies are needed, preferably using different methods. Specialised but numerically efficient models Burin et al. (2001) cannot incorporate spatial correlations of the dielectric constant. Two-dimensional FDTD simulations have already been used to describe random lasers Sebbah and Vanneste (2002). They need to make only minimal assumptions and they can be extended to include arbitrary spatial correlations. FDTD simulations thus might be a good candidate but diffusive samples need to be larger than the localised samples studied so far. Given that FDTD is computationally very expensive, it is not obvious to us whether this would still be numerically feasible.

This work was supported by the European Union Marie Curie fellowship HPMF-CT-2002-01794.