text:
# An evolutionary model for pulsar-driven supernova remnants

## 1 Introduction

A supernova remnant (SNR) is the relic of a supernova explosion, which injects an energy of $`10^{51}`$ erg into the surrounding medium. The dynamics of a young SNR expanding into the interstellar medium (ISM) is determined by two shocks: a forward shock which propagates into the ISM and is being decelerated by sweeping up material from the ISM, and a reverse shock which results from the high pressure behind the forward shock and propagates back into the freely expanding ejecta of the SNR (McKee 1974; McKee & Truelove 1995).

A young SNR becomes dynamically more interesting in those cases where the collapse of the progenitor star, preceding the supernova explosion, yields a pulsar: a rapidly rotating neutron star. In those cases the dynamics of the central region of the SNR is dominated by the continuous injection of energetic particles by a relativistic pulsar wind, driven by the spin-down energy of the pulsar. The pulsar wind is terminated by a strong MHD shock (Rees & Gunn 1974), and drives a pulsar wind nebula (PWN) in the interior of the young SNR. The dynamics of the PWN is coupled to the evolution of the SNR, because the total energy release over the pulsar’s lifetime is small ($`10^{49}10^{50}`$ erg) compared with the total mechanical energy of the SNR ($`10^{51}`$ erg).

Several authors (Reynolds & Chevalier 1984; van der Swaluw et al. 2001; Blondin et al. 2001; Bucciantini et al. 2003) have considered the evolution of a centered PWN inside an evolving SNR. In these systems, the initial stage of the PWN is the supersonic expansion stage: the pulsar wind bubble is bounded by a strong PWN shock propagating through the freely expanding ejecta of the SNR. A transition to the reverse shock interaction stage takes place when the reverse shock collides with the PWN shock, and subsequently crushes the pulsar wind bubble (van der Swaluw et al. 2001, Blondin et al. 2001). This reverse shock interaction stage is characterised by an unsteady expansion of the pulsar wind bubble (van der Swaluw et al. 2001), due to the reverberations from the violent collision between the reverse shock and the PWN shock. The expansion of the PWN proceeds subsonically when these reverberations have vanished. The expansion is subsonic because the surroundings of the PWN have been reheated by the passage of the reverse shock: the PWN shock, bounding the swept-up shocked ejecta around the hot pulsar wind bubble, has disappeared.

In this paper we discuss the evolution of PWNe inside SNRs, for pulsars with a constant kick velocity: initially the PWN starts its expansion at the center of the SNR, however the pulsar motion will move the pulsar wind cavity along as it moves through the SNR interior. We present results from hydrodynamical simulations for such a system, which distinguishes all three evolutionary stages mentioned above, i.e. the supersonic expansion stage, the reverse shock interaction stage and the subsonic expansion stage. However at the end of the simulation an additional stage can be distinguished, when the head of the PWN, containing the active pulsar, deforms into a bow shock, due to the supersonic motion of the pulsar.

The supersonic expansion stage in the simulation shows a PWN which is off-centered with respect to the twofold shock structure of the SNR, due to the kick velocity of the pulsar. Therefore the timescale on which the reverse shock collides with the complete shock surface, bounding the PWN, can be a significant fraction of the total lifetime of the PWN when the reverse shock interaction stage starts. We use a semi-analytical approach to estimate the timescale of the collision process, which is shown to scale roughy with the pulsar velocity. The results from the simulation are in almost complete agreement with these semi-analytical calculations. Due to the high velocity of the pulsar, its position inside the PWN is strongly off-centered after the passage of the reverse shock. Ultimately, at the end of the simulation, when the pulsar is approaching the shell of its SNR, the head of the PWN, containing the active pulsar, is deformed into a bow shock at a time and position, in complete agreement with analytical predictions made by van der Swaluw et al. (1998).

We determine the evolutionary stage of the PWNe inside the SNRs N157B, G327.1-1.1 and W44, using our model. For the first two systems, we argue that the position of the pulsar at the head of its PWN, is a result of the passage of the reverse shock and the high velocity of the pulsar: the head of these PWNe are not bounded by a bow shock. Therefore these PWNe are either in the reverse shock interaction stage or the subsonic expansion stage. The PWN inside the SNR W44 however, is shown to be a good candidate for having a bow shock nebula around its pulsar.

## 2 The evolution of a PWN inside a SNR

### 2.1 The interaction between the reverse shock and the PWN shock

The initial stage of PWN evolution is characterised by a hot pulsar wind bubble, bounded by a strong PWN shock, propagating through the freely expanding ejecta of the SNR. A transition to the subsonic expansion stage occurs via the reverse shock interaction stage. This interaction stage starts when the reverse shock collides with the PWN shock. In the case of a centered pulsar ($`V_{\mathrm{psr}}=0`$), the reverse shock collides with the PWN shock surface instantaneously, due to the spherical symmetry of the SNR and the PWN. However, when the pulsar has a kick velocity, there will be an associated timescale on which the reverse shock collides with the complete surface of the PWN shock. This is the first stage of the reverse shock interaction stage. Next, the pulsar wind bubble oscillates back and forward due to the presence of reverberations from the passage of the reverse shock (van der Swaluw et al. 2001). The reverse shock interaction stage ends when these reverberations have vanished and the pulsar wind bubble proceeds its expansion subsonically. In this section we use a semi-analytical approach to derive a timescale on which the reverse shock collides with the complete surface of the PWN shock.

McKee & Truelove (1995) give analytical approximations for the trajectories of the forward shock and the reverse shock of a SNR in the case of a uniform ambient medium. Their equations for the trajectories of both shocks are normalised to a timescale $`t_{\mathrm{ST}}`$, which marks the age of the remnant when it has swept up roughly 1.61 the ejected mass $`M_{\mathrm{ej}}`$. Their equations descibe the expansion of an isolated SNR in the free expansion stage and the Sedov-Taylor stage. The trajectory of the forward shock converges to the Sedov-Taylor solution when the SNR age $`t>>t_{\mathrm{ST}}`$. Internal (radiative) cooling is not inlcuded in their model, therefore the pressure-driven snowplow stage is not described (see however Blondin et al. 1998).

The reverse shock hits the center of the SNR at approximately 5 times the transition time $`t_{\mathrm{ST}}`$, putting an upper limit on the age $`t_{\mathrm{col}}`$, when the reverse shock has collided with the complete shock surface bounding the PWN:

$$t_{\mathrm{col}}\mathrm{\hspace{0.33em}5}t_{\mathrm{ST}}=\mathrm{\hspace{0.33em}1\hspace{0.33em}045}E_{51}^{1/2}\left(\frac{M_{\mathrm{ej}}}{M_{}}\right)^{5/6}n_0^{1/3}\mathrm{years},$$
(1)
here $`E_{51}`$ is the total mechanical energy of the SNR in units of $`10^{51}`$ erg and $`n_0`$ is the ambient hydrogen number density assuming an interstellar composition of 10 H : 1 He. The above timescale is very close to the one given by Reynolds& Chevalier (1984).

An equation for the radius of the PWN shock, when it is interacting with the freely expanding ejecta of its SNR was given by van der Swaluw et al. (2001), where a constant pulsar wind luminosity $`L_{\mathrm{pw}}=L_0`$ was taken:

$$R_{\mathrm{pwn}}(t)\mathrm{\hspace{0.33em}0.889}\left(\frac{L_0t}{E_0}\right)^{1/5}V_0tt^{6/5},$$
(2)
here $`E_0`$ is the total mechanical energy of the SNR and $`V_0`$ is defined as:

$$V_0=\sqrt{\frac{10}{3}\frac{E_0}{M_{\mathrm{ej}}}}.$$
(3)
We have taken the case for which the adiabatic heat ratio of the pulsar wind material equals $`\gamma _{\mathrm{pwn}}=5/3`$.

We will use the equations for the trajectories of the forward and the reverse shock of the SNR from Truelove & McKee (1995) and the above equation for the radius of the PWN shock to calculate the collision time $`t_{\mathrm{col}}`$, i.e. when the reverse shock collides with the PWN shock. We make a distinction between the front of the PWN shock $`R_{\mathrm{pwn}}^\mathrm{f}=R_{\mathrm{pwn}}+V_{\mathrm{psr}}t`$, and the back of the PWN shock $`R_{\mathrm{pwn}}^\mathrm{b}=R_{\mathrm{pwn}}V_{\mathrm{psr}}t`$, which enables us to calculate $`\mathrm{\Delta }t_{\mathrm{col}}`$, the timescale on which the reverse shock collides with the complete shock surface bounding the PWN. This approach neglects the change in the PWN expansion due to the displacement of the pulsar inside the freely expanding ejecta. Figure 1 shows an example of this approach, where the position of the four different shocks are plotted as a function of time. The parameters have been taken similarly to the parameters which we will use in the simulation. The figure shows that the timescale on which the reverse shock collides with the complete PWN shock surface equals $`\mathrm{\Delta }t_{\mathrm{col}}`$ 1050 years. Figure 2 shows $`\mathrm{\Delta }t_{\mathrm{col}}`$ as a function of the pulsar velocity, where one observes that this timescale correlates almost linearly with the pulsar velocity. We conclude that the collision timescale $`\mathrm{\Delta }t_{\mathrm{col}}`$ of PWNe containing a high velocity pulsar, can be a significant fraction of the timescale associated with the supersonic expansion stage.

### 2.2 The formation of the PWN bow shock

After the passage of the reverse shock the PWN consists of two distinct parts:

* The roughly spherically symmetric PWN relic, from the initial energetic stage of the pulsar wind, which will ultimately expand subsonically, after the reverberations of the reverse shock have vanished.
* The relic PWN connects with the head of the PWN, directed towards the shell of the remnant. The head of the PWN contains the active pulsar, which propagates towards the SNR shell.

The age of the remnant when the pulsar overtakes the shell of its remnant was given by van der Swaluw et al. (2003):

$$t_{\mathrm{cr}}\mathrm{\hspace{0.33em}1.4}\times 10^4E_{51}^{1/3}V_{1000}^{5/3}n_0^{1/3}\mathrm{years}.$$
(4)
$`V_{1000}`$ denotes the velocity of the pulsar in units of 1000 km/sec. This timescale was calculated in the limit of a Sedov-Taylor SNR, i.e. it is assumed that the pulsar will break through the shell of the remnant, when radiative losses of the SNR are neglible.

We follow van der Swaluw et al. (1998) by calculating the Mach number $`_{\mathrm{psr}}`$ of the pulsar as it propagates through the SNR interior, using the Sedov-Taylor solution (Sedov 1959). Figure 3 shows the Mach number $`_{\mathrm{psr}}`$ as a function of the age of the remnant $`t`$, normalised to the crossing time $`t_{\mathrm{cr}}`$. Alternatively Figure 4 shows the Mach number $`_{\mathrm{psr}}`$ as a function of the position of the pulsar $`R_{\mathrm{psr}}`$, normalised to the position of the SNR blastwave $`R_{\mathrm{snr}}`$.

The figures show that the PWN head, containing the active pulsar, is deformed into a bow shock (i.e. when $`_{\mathrm{psr}}=1`$), at half the crossing time when the pulsar is positioned at a radius of $`0.677`$ times the radius of the SNR blastwave. The Mach number will slowly increase after the formation of the bow shock, due to the decrease of the sound speed, as the PWN is approaching the shell of the remnant. When the pulsar breaks through the shell, the Mach number equals $`_{\mathrm{psr}}=\frac{7}{\sqrt{5}}3.13`$ (van der Swaluw et al. 2003).

A lower limit for the pulsar velocity can be derived, such that the bow shock formation occurs while the SNR is in the Sedov-Taylor stage. We follow van der Swaluw et al. (2003), who already give a lower limit for the pulsar velocity in order for the pulsar to cross the SNR shell in its Sedov-Taylor stage. They use the transition time calculated by Blondin et al. (1998), as the age of the SNR when the transition from the Sedov-Taylor stage to the pressure-driven snowplow stage occurs:

$$t_{\mathrm{PDS}}=2.9\times 10^4E_{51}^{4/17}n_0^{9/17}\mathrm{yr}.$$
(5)
Our results show that the bow shock creation occurs at half the crossing time, therefore the requirement $`0.5t_{\mathrm{cr}}t_{\mathrm{PDS}}`$ yields a lower limit for the pulsar velocity to let the bow shock formation event occur when the SNR is in its Sedov-Taylor stage:

$$V_{\mathrm{psr}}325n_0^{2/17}E_{51}^{1/17}\mathrm{km}/\mathrm{s},$$
(6)
which is a reasonable fraction of the observed pulsar velocities by Hansen and Phinney (1997).

## 3 Numerical Simulations

### 3.1 Simulation method

We use a second order, properly upwinded hydrodynamics code (described in Downes & Ray,d&r99 (1999)) to simulate the dynamics of the interaction between a pulsar wind and a supernova remnant when the pulsar has a high velocity. The hydrodynamics equations are intergrated in cylindrical symmetry, and the boundary conditions are taken as gradient zero everywhere except on the $`r=0`$ boundary, where they are set to reflecting.

### 3.2 Initialising the system

The simulations are performed in the rest frame of the pulsar. The pulsar velocity equals $`V_{\mathrm{psr}}=\mathrm{1\hspace{0.33em}000}`$ km/sec. Initially, the ambient density and pressure are uniform with a uniform velocity of $`\mathrm{1\hspace{0.33em}000}`$ km s-1 in the negative $`z`$ direction. The ambient density is $`\rho _0=10^{24}`$ g/cm-3 (corresponding with $`n_00.427`$), and the pressure is chosen so that the ambient temperature is $`T=1\times 10^5`$ K. The resolution is taken as $`\mathrm{\Delta }r=\mathrm{\Delta }z=3.6\times 10^{16}`$ cm.

An exanding SNR is created by initialising a sphere of radius 0.25 pc, or 21 grid cells, with a high thermal energy and density such that the total energy contained in the sphere is $`E_0=10^{51}`$ ergs, while the ejecta mass is $`M_{\mathrm{ej}}=3M_{}`$. Note that this material also has a velocity of $`\mathrm{1\hspace{0.33em}000}`$ km s-1 in the negative $`z`$ direction.

The source of the pulsar wind is modeled as a small sphere of radius 0.175 pc (15 grid cells) into which thermal energy is deposited at a rate of $`L_{\mathrm{pw}}=10^{38}`$ erg/s. Mass is also deposited into this region at a rate chosen such that the pulsar wind terminal velocity is $`v_{\mathrm{}}=\mathrm{30\hspace{0.33em}000}`$ km s-1. The pulsar wind luminosity has been taken constant, in order to resolve the pulsar wind termination shock throughout the simulation. The results of the simulation will not change qualitatively, because the total injected energy by the pulsar during its stay in the SNR interior $`E_{\mathrm{pw}}=L_{\mathrm{pw}}t_{\mathrm{cr}}E_0`$.

### 3.3 The reverse shock interaction stage

The initial stage of the PWN evolution is the supersonic expansion stage. Figure 5 shows the density profile of the PWN/SNR system during this stage. One can clearly distinguish the four different shocks of the system: the pulsar wind termination shock, the PWN shock, the reverse shock and the forward shock. One clearly observes the off-centered position of the PWN with respect to the twofold shock structure of the SNR, caused by the motion of the pulsar. Furthermore the pulsar wind cavity is roughly spherically symmetric. This results from the high sound speed inside the pulsar wind bubble, which smooths out pressure perturbations rapidly, yielding an almost uniform pressure distribution in the PWN.

Figure 6 shows the density profile of the PWN/SNR system, when the front of the PWN shock has already collided with the reverse shock. At this stage the pulsar position starts to get off-centered with respect to the PWN structure. This is caused by the reverse shock interaction, which raises the presssure inside the disturbed parts of the bubble, which results in pulsar wind material flowing towards the back of the PWN. Furthermore the increase in the pressure inside the front part of the pulsar wind bubble, pushes the forward termination shock backwards (towards the pulsar) as well. On the other hand, the downstream pressure at the backward termination shock has not been influenced by the reverse shock interaction yet. Therefore the pulsar wind cavity structure is no longer spherically symmetric as can be seen in Figure 6. The total time on which the reverse shock collides with the complete shock structure bounding the pulsar wind bubble equals $`\mathrm{\Delta }t_{\mathrm{col}}\mathrm{1\hspace{0.33em}200}`$ years in the simulation, after which the age of the SNR equals $`t\mathrm{2\hspace{0.33em}400}`$ years. These timescales are in almost complete agreement with the analysis performed in section 2.1 of this paper, which validates the approach made in that section, to derive the collision timescale $`\mathrm{\Delta }t_{\mathrm{col}}`$.

Figure 7 shows the PWN/SNR system shortly after the passage of the reverse shock. The pulsar is positioned at the head of the PWN, due to the passage of the reverse shock and the high velocity of the pulsar. The PWN structure consists of roughly two parts, a relic PWN and the head of the PWN, as was described in section 2.2 of this paper. After the passage of the reverse shock, the reverse shock interaction stage is continued by the reflected and transmitted shock waves, which propagate through both the pulsar wind bubble and the surrounding ejecta of the SNR. This stage was described in detail by van der Swaluw et al. (2001) for a spherically symmetric PWN/SNR system. We observe a qualitatively similar process, in which the pulsar wind bubble is oscillating. The dynamics are more complicated compared with the simulations of van der Swaluw et al. (2001), due to the asymmetry of the PWN/SNR system caused by the motion of the pulsar. During this stage of the reverse shock interaction, the PWN is oscillating between expansion and compression and the simulations reveal Rayleigh-Taylor and Kelvin-Helmholtz instabilities. These instabilities lead to the mixing of ejecta and pulsar wind material, as was similarly observed in the simulations performed by Blondin et al. (2001). Therefore at the end of the reverse shock interaction stage, the relic PWN consists of a mixture of pulsar wind material and ejecta. In our simulation the relic PWN will move off the grid, before the bow shock formation occurs around the head of the PWN.

### 3.4 The formation of the bow shock

At the end of the simulation, as the pulsar approaches the shell of the SNR, the head of the PWN, containing the active pulsar, deforms into a bow shock, due to the supersonic motion of the pulsar. The simulation shows that the bow shock formation event occurs at roughly $`t0.5t_{\mathrm{cr}}`$, when the position of the pulsar $`R_{\mathrm{psr}}`$ with respect to the radius of the blastwave $`R_{\mathrm{snr}}`$ equals $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$. This clearly demonstrates the validity of the analytical approximation made in section 2.2 of this paper. Figure 8 shows the density profile of the PWN/SNR system after the bow shock formation. The curvature of the SNR shock is small compared with the bow shock structure, which validates the assumptions made by van der Swaluw et al. (2003) to model the break-through event.

## 4 Observations

After the passage of the reverse shock, the PWN consists of a relic PWN and the head of the PWN, containing the active pulsar. The figures 7 and 8 show the density profiles of the PWN/SNR before and after the formation of the bow shock. For both cases the pulsar is located inside the head of the PWN. The synchrotron maps from PWNe will therefore show a similar morphology, which makes it hard to determine the presence of a bow shock based on the observed PWN/SNR morphology. However, we have demonstrated that the following criteria can be used to determine the presence of a bow shock around the head of the PWN:

* The ratio between the position of the pulsar $`R_{\mathrm{psr}}`$ with respect to the position of the forward shock of the SNR $`R_{\mathrm{snr}}`$ has to satisfy $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$
* The age of the remnant has to be larger then half the crossing time, i.e. $`t>0.5t_{\mathrm{cr}}`$.

We will apply the above two criteria to three SNRs: N157B, G327.1-1.1 and W44. We will argue that the head of the PWNe inside the first two SNRs are not bounded by a bow shock, and therefore these PWNe are in either the reverse shock interaction stage or the subsonic expansion stage. The PWN inside SNR W44 is shown to be the best candidate of a PWN in the bow shock stage.

### 4.1 The pulsar wind nebula inside N157B

N157B is a young SNR dominated by plerionic emission from the PWN inside this remnant. The age of the remnant is approximated to be $`\mathrm{5\hspace{0.33em}000}`$ yrs (Wang & Gotthelf 1998) and contains a 16 ms pulsar (Marshall et al. 1998). The velocity of the pulsar is high ($`V_{\mathrm{psr}}\mathrm{1\hspace{0.33em}000}`$ km/sec), if one assumes that the pulsar was born at the central region of the bright radio emission (Lazendic et al. 2000).

Wang & Gotthelf (1998) argue for a bow shock interpretation of the PWN in N157B: the supersonic motion of the pulsar has deformed the PWN around the pulsar into a bow shock. In this scenario the SNR N157B is a young variant of the SNR CTB80, which is thought to have a PWN bow shock located close to the shell of the remnant (see e.g. Strom 1987; Kulkarni et al. 1988; Migliazzo et al. 2003). It is remarkable though, that the spindown luminosity of the pulsar inside N157B ($`L_04.8\times 10^{38}`$ ergs/sec) is so much higher compared with the spindown luminosity from the pulsar inside CTB80 ($`L_04.0\times 10^{36}`$ ergs/sec). If there is a bow shock around the PWN inside N157B, this implies that an upper limit for the age of the SNR, when the pulsar crosses the shell, is approximately 10 000 years (using the current age of 5 000 years and the criterion that bow shock creation occurs at half the crossing time). This is again a remarkable contrast with the current age of CTB80 of 100 000 years, which is close to the crossing time of this SNR. Furthermore from figure 2 of Wang et al. (2001) it seems like the position of the pulsar is more or less centered in the SNR. This is in contrast with what one would expect from the analysis performed in section 2.2 of this paper, which predicts $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$.

We propose a scenario for SNR N157B in which the contradictions mentioned above disappear. Following the results from sections 2 and 3, we suggest that the central bright parts of the radio and X-ray emission inside N157B correspond with the relic PWN, being blown in the initial stage of the PWN, when it was expanding supersonically. The PWN inside N157B has just collided with the reverse shock, this collision and the high velocity of the pulsar has off-centered the pulsar with respect to the PWN structure as was discussed in the sections 2 and 3.3 of this paper. Therefore the head of the PWN is not bounded by a bow shock: the PWN is in the aftermath of the reverse shock interaction stage or in the beginning of the subsonic expansion stage.

### 4.2 The pulsar wind nebula inside G327.1-1.1

Another example of a young SNR, containing a PWN with the position of the pulsar off-centered with respect to its PWN is G327.1-1.1. There has not been a pulsar detected inside this remnant so far. Sun et al. (1999) presented a radio contour map using MOST data overlaid with X-ray data from ROSAT. The X-ray emission is centered around a finger of radio emission sticking out of a central radio bright region, indicating the presence of a pulsar wind. Following Sun et al. (1999) the SNR can be modelled in X-rays by the following set of parameters: $`E_{51}=0.23`$, $`n_0=0.10`$, $`V_{\mathrm{psr}}`$=600km/sec and an age of $`t=1.1\times 10^4`$. Using equation (4) to calculate $`t_{\mathrm{cr}}`$ we get a value of $`4.3\times 10^4`$ years. The age of the system therefore is less than half the crossing time, implying the absence of a bow shock. Furthermore the displacement of the PWN finger (i.e. the head of the PWN), containing the pulsar, does not satisfy the other criterion for the presence of a bow shock, i.e. $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$. This criterion is less restrictive, since one might introduce an angle between the pulsar velocity and the observer such that the ratio $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$. Notice however that this will not influence the age of the remnant!

Therefore we propose a scenario for G327.1-1.1 in which there is no bow shock around the head of the PWN: the pulsar has been positioned at he head of the PWN due to the high velocity pulsar and the passage of the reverse shock. The central bright part of the remnant corresponds with the relic PWN. The finger of the PWN corresponds with the head of the PWN, which contains the active pulsar. Because of the larger value of the ration $`R_{\mathrm{psr}}/R_{\mathrm{snr}}`$ and the larger age with respect to the PWN inside SNR N157B we favour the scenario where the PWN is in its subsonic expansion stage.

### 4.3 The pulsar wind nebula inside W44

The SNR W44 is an older remnant then the previous two SNRs discussed. Furthermore only a small fraction of the radio emission from this remnant is characterised by plerionic emission (Frail et al. 1996). Taking the characteristic age of the pulsar, $`\mathrm{20\hspace{0.33em}000}`$ years (Wolszcan et al. 1991), as the age of this remnant yields an upper limit of $`t_{\mathrm{cr}}\mathrm{40\hspace{0.33em}000}`$ years for the crossing time of this remnant. This age is much closer to the age of the SNR CTB80. The displacement of the pulsar position only marginally violates the condition for bow shock formation, i.e. $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$. We therefore conclude that the PWN observed in the SNR W44 corresponds with the head of the PWN, which has been deformed into a bow shock.

## 5 Conclusions

We have considered the case of a PWN interacting with a SNR, for which the associated pulsar is moving at a high velocity through the interior of its SNR. The model we discussed made use of a hydrodynamics code. One could distinguish four different stages in the simulation: the supersonic expansion stage, the reverse shock interaction stage, the subsonic expansion stage and the bow shock stage. Below we summarise the most important results from our model:

* The reverse shock interaction stage starts when the reverse shock collides with the front of the PWN shock. The timescale on which the reverse shock collides with the complete PWN shock, bounding the pulsar wind bubble, scales with the pulsar velocity. This timescale is a significant fraction of the lifetime of the PWN, when the reverse shock interaction stage starts, for PWNe containing a high velocity pulsar.
* The high velocity of the pulsar results in an off-centered position of the pulsar with respect to its pulsar wind bubble, after the passage of the reverse shock.
* The morphology of a PWN, after the passage of the reverse shock, consists of a twofold structure which ultimately expands subsonically inside the relaxed Sedov-Taylor SNR:

1) a roughly spherically symmetric relic PWN

2) a head, containing the pulsar, directed towards the SNR shell
* The formation of a bow shock around the head of the PWN occurs at half the crossing time, when the position of the pulsar $`R_{\mathrm{psr}}`$ with respect to the SNR shock $`R_{\mathrm{snr}}`$ equals $`R_{\mathrm{psr}}0.677R_{\mathrm{snr}}`$. We derived a lower limit for the pulsar velocity for the bow shock formation to occur while the SNR is in its Sedov-taylor stage.

From our model it follows that a SNR containing a pulsar at the head of its PWN does not imply the presence of a bow shock, bounding the head feature of the PWN. We discussed three SNRs and interpreted their morphology, using the results from our model. For the SNRs N157B and G327.1-1.1 we argued that they contain PWNe which do not have a bow shock. The PWN inside the SNR W44 seems to be a good candidate for having a bow shock around the head of its PWN.

###### Acknowledgements.

This work was part-funded by the CosmoGrid project, funded under the Programme for Research in Third Level Institutions (PRTLI) administered by the Irish Higher Education Authority under the National Development Plan and with partial support from the European Regional Development Fund.


format: text

text:
# LCG-1 Deployment and usage experience

## 1 Brief introduction to LCG-1

LHC Computing Grid (LCG) is one of the five CERN projects at the moment. The main goal of the project is is to prepare the computing infrastructure for the simulation, processing and analysis of LHC data for all four of the LHC collaborations. This infrastructure is based on the Grid model of distributed computing . In general LCG is the middleware project, which means that it will hide much of the complexity of the environment from the end user.

LCG is based on the software developed in the EU-DataGrid project  and VDT  which trace their roots back to the Globus Toolkit .

## 2 LCG-1 Deployment scheme

### 2.1 LCG-1 Deployment overview

The deployment model for the LCG-1 includes a central CVS repository which stores the recommended configuration templates for the current LCG-1 release. The main deployment platform at the moment is RedHat Linux 7.3. LCG-1 deployment team provides binary RPM packages for all LCG-1 middleware and LHC experiments software.

The recommended installation procedure for the LCG-1 is the following. The LCFGng configuration server is installed on the site, and the hosts running LCG-1 services are configured automatically from the LCFGng server.

The configuration for the site should be based on the current LCG-1 release version, obtained from the central CVS. For each of the sites from the official deployment list a module in CVS is created, and the site configuration is then uploaded to CVS (see fig. 1).

### 2.2 LCG-1 Deployment at MOSCOW site (SINP MSU)

As a part of the deployment test LCG-1 software was installed on the SINP MSU PC farm (LCG-1 site name for the site is MOSCOW). The farm has 20 dual-CPU nodes and two 1.2 TB fileservers. The fileservers have Gigabit Ethernet uplinks and the nodes are connected with the FastEthernet network.

At the moment of the LCG-1 deployment the was an automatic software installation system running on our PC farm. It is based on Etherboot network boot package and anaconda kickstart technology. There is a configuration server running DHCP and TFTP services, which stores different network boot images with preconfigured kickstart scripts for anaconda. There is also a possibility not to use a kickstart script and to install the node manually. In general an installation of a node from scratch requires a boot ROM on the node, but for the nodes without a boot ROM there is a possibility to boot from a special boot disk which will emulate the boot ROM. Another possibility is to install a fake Linux kernel which also emulates a network boot.

For the LCG-1 deployment installation system was modified in the following way. The LCFGng server software was installed on the node running DHCP and TFTP services, and a new network boot image for LCFGng was added. This image is based on the standard LCFGng boot disk provided with the LCFGng software. During the installation process the LCFGng managed node obtains the configuration information from the server using HTTP and LCFG protocols and than installs and configures the software packages shared on the NFS server.

Only a part of the SINP MSU PC farm was configured for LCG-1 tests. The nodes selected for the LCG-1 were disconnected from the main farm. The following LCG-1 components were installed: Computing Element, Worker Node, Storage Element, Resource Broker, BDII and User Interface. The MDS and MyProxy nodes will be installed in the nearest future.

We also had some experience installing LCG-1 software on other PC farms in the MSU. First of all, SRCC MSU parallel cluster runs its own batch system developed by SRCC MSU and CMC MSU. The interface between LCG middleware and this batch system was done through the extra package to the Globus Toolkit, but this interfaces is not yet completely tested. Other general differences between a plain LCG-1 installation and installation on this cluster are:

* Manual node configuration only. The LCFGng cannot be installed on this cluster due to some technical limitations.
* Middleware for the Worker Nodes is installed on the shared file system.

In the nearest future we are going to install LCG middleware on the Physical faculty and the faculty of Computing Mathematics and Cybernetics clusters in the MSU and connect these clusters to the SINP MSU Resource Broker. However, installation on these clusters has the similar limitations as the installation on the SRCC MSU parallel cluster.

## 3 Conclusion

At present installation of LCG-1 middleware has both big advantages and disadvantages.

### 3.1 Difficulties

One of the main difficulties with the LCG-1 installation for now is the lack of available documentation. The only well-documented way to install LCG software requires using of the LCFGng configuration server. This method is not applicable in a number of cases.

Another issue is that the minimal standalone site must be running a big number of nodes just for supporting the infrastructure. The minimal set of nodes required for a standalone site is a Resource Broker, BDII & MDS, Computing Element with at least one Worker Node, at least one Storage Element, and the MyProxy server if long-time jobs support is required. This gives us at least five nodes supporting the infrastructure which cannot be used for other needs such as running jobs for example.

### 3.2 Features

The summary of the features provided by the LCG-1 for now is the following:

* Convenient way for job balancing between several sites.
* Common way of user authentication and authorization for job submission.
* Some basic accounting over different sites.
* Data replication.

## 4 Acknowledgments

The work was partially supported by CERN-INTAS grant 00-0440.


text:
# Investigating the central engine of Seyfert 2 galaxies with and without Polarized Broad Lines

## 1 Introduction

Seyfert galaxies are classified into two spectroscopic groups based on the presence (type 1) or absence (type 2) of broad permitted optical emission lines. The first discovery of Polarized Broad Lines (PBL hereafter) in NGC 1068 (Antonucci & Miller, 1985) and later in several other Seyfert 2 galaxies (Miller & Goodrich, 1990; Tran, 1995; Young et al., 1996; Heisler et al., 1997; Awaki et al., 2000; Moran et al., 2000; Alexander, 2001; Gu et al., 2001; Lumsden et al., 2001) suggests that Seyfert 2 galaxies harbor a bright Seyfert 1 nucleus hidden from our view by an optically and geometrically thick obscuring torus (“Unified models”, see the review of Antonucci (1993)). In this frame, the difference between type 1 and 2 Seyfert galaxies is only due to the viewing angle. X-ray observations support this model showing the presence of large column densities along the line of sight of Seyfert 2 galaxies (Risaliti et al., 1999; Guainazzi et al., 2001). However, other observations are inconsistent with unified models, in particular in the hard X-ray domain where Zdziarski et al. (1995) and Deluit & Courvoisier (2003) (D03 hereafter) show differences in the spectral index, energy cutoff and amount of reflection between Sy 1 and Sy 2 galaxies. Furthermore, Alexander (2001) refutes the absorption argument to explain the absence of broad lines in some Sy 2 galaxies showing that the density of the absorbing medium is not linked to the detectability of broad lines.
The polarization is presumed to be due to nuclear emission scattered by a gas of electrons located further out than the obscuring material. According to current models, the detection of polarized broad emission lines in type 2 Seyfert galaxies implies that these objects have a Hidden Broad Line Region (HBLR henceforth). However, not all Sy 2 galaxies that have been observed in polarized light show the presence of broad lines. This might imply an intrinsic difference between Seyfert 2 with PBL (PBLs hereafter) and those with No PBL (NPBLs hereafter). We propose to test this by studying their hard X-ray intrinsic emission in a spectral domain where absorption does not modify significantly the primary emission.
The existence of Sy 2 without PBL could be due to the absence of a broad line region in some Sy 2 galaxies or to the absence of a scattering medium or different physical conditions within the inner regions.
To investigate these hypotheses, we intend to describe the general emission properties of various classes of Seyfert 2 galaxies by computing their average hard X-ray spectrum.
In D03 we computed the average spectra of Sy 1, Sy 1.5 and Sy 2 galaxies. We found that the Sy 1 and Sy 2 galaxies of our sample have a different behavior, in particular concerning the presence of a cutoff which is required in the Sy 1 emission and absent in the Sy 2 average spectrum. Another point of disagreement was the role played by Compton reflection possibly influencing the Sy 2 emission, unlike Sy 1.
We pursue our investigations within the Seyfert 2 class of our sample to search for common properties between Sy 1 galaxies and the two Sy 2 subclasses. To test the argument that the absorption is the only reason for which we cannot see the HBLR, we consider Compton-thin Sy 2 galaxies. After presenting the average spectra obtained, we study physical processes responsible for the observed X$`\gamma `$ emission.

## 2 The Dataset

The initial sample of Seyfert 2 galaxies is that presented in D03. From the initial 22 Seyfert 2 galaxies we keep 16 objects for a total of 24 observations. Five objects have been observed several times. We apply selection criteria, presented in subsections 2.1. and 2.2., to the initial sample.

### 2.1 Spectropolarimetric data

This study is based on a sample composed of Seyfert 2 galaxies of two types: one presenting broad lines in polarized light, the second for which no PBL have been detected. We point out the fact that the non detection of broad lines does not prove definitively their genuine absence, in particular if we consider the limitations of current spectropolarimetric instrumentation. However, the purpose of this study is to start a comparison of the two subclasses of Seyfert 2 galaxies in the hard X-ray domain, using the best available published information on the presence of hidden broad lines.
The current classification of Seyfert galaxies is mainly determined by optical criteria, summarized in Table 1.
We collected all Seyfert 2 galaxies from the recent literature (from 1980 to 2003) for which spectropolarimetric data are available. Several objects of the D03 Sy 2 sample have been excluded: IRAS 18325-5926 because of its NED classification as Sy 2 which has been questioned by Iwasawa et al. (1995); NGC 4939 is a Sy 2 but no spectropolarimetric data is available. NGC 5674 and ESO 103-G35 have been excluded because of their Sy 1.9 classification. Indeed, current models generally claim that the detection of a broad H$`\alpha `$ line implies the presence of a hidden broad line region and explain the absence of the H$`\beta `$ line as an obscuration of this region, leading to the conclusion that Sy 1.9 galaxies harbor a hidden Sy 1 nucleus.
In addition to spectropolarimetric considerations, we use only data with an integrated signal to noise ratio higher than 2 in the PDS energy range (15-136 keV). This threshold is high enough to conclude a detection by the PDS detector for a known source. We therefore consider in our study all Sy 2 galaxies for which both spectropolarimetric information and PDS data of sufficient quality are available.

### 2.2 Sources not significantly influenced by absorption effects

As exhaustively explained in D03, we limit our study to Sy 2 galaxies having a hydrogen column density less than 7$``$1023cm-2. We have established that below this threshold, the absorption effects do not alter strongly the primary emission, allowing us a direct view of the emission produced by the central engine.

### 2.3 Sample spectral characteristics

The final sample is composed of 8 Seyfert 2 galaxies with PBL and 8 without PBL detected (Table 2).

#### 2.3.1 Absorption distribution

We present in Figure 1 the hydrogen column density distribution of PBLs and NPBLs Sy 2 galaxies.

Several contradictory results have been found during recent years. Heisler et al. (1997) and Gu et al. (2001) found a higher column density for NPBLs whereas Alexander (2001) found no link between the detectability of polarized broad lines and the density of the absorbing material.
In our study, the column density distributions of the two Sy 2 classes are not significantly different. In particular the NPBLs galaxies are not more absorbed than PBLs since they have the same number of objects with a low absorption (NH$`<`$1022cm-2) than with a higher absorption (NH$`>`$1022cm-2). The fact that we consider a small number of objects and only Compton-thin Sy 2 galaxies prevents us from generalizing this result for the overall Sy 2 class. We however show that even in Compton-thin Sy 2 galaxies the absence of polarized broad lines is not unusual.

#### 2.3.2 Luminosity distribution

The hard X-ray luminosity is a signature of the nucleus activity.
The comparison of luminosity distributions of the two Sy 2 subclasses (Figure 2) reveals that NPBLs are weaker than PBLs, indicating that the activity of their respective central engine is slightly different. Indeed, 87.5 % of PBLs have a luminosity between 1043 and 1044 erg s-1 against only 25% of NPBLs. 50% of the NPBLs population has a luminosity between 1042 and 1043 erg s-1.
The hypothesis according to which PBLs would be more luminous than NPBLs has been evoked in other wavelength domains (e.g. Kay (1994) and Tran (2003)); we discuss the implications of this result in the section 7.

## 3 Data Analysis

The Narrow Field Interments (NFI) of the BeppoSAX satellite is composed of a Low Energy Concentrator Spectrometer (LECS; (Parmar et al., 1997)), three Medium Energy Concentrator Spectrometers (MECS; Boella et al. (1997)) and a Phoswich Detector Counter (PDS; Frontera et al. (1997)). LECS and MECS have imaging capabilities and operate respectively in the 0.1-5 keV and 2.0-10 keV domains. The PDS, used in our study, is an instrument functioning in an energy band between 15-200 keV with an energy resolution of 15% at 60 keV. The detector is designed to allow a good control of background variations using rocking collimators that periodically sample source+background combinations and background alone.
The data analysis has been performed with the XSPEC (Arnaud, 1996) version 11 package using the latest PDS response matrices released by the BeppoSAX SDC111http://www.asdc.asi.it/bepposax/. To compute the average spectra, we used the task Mathpha version 5.1. combining the data files. The errors have been propagated at each step of the procedure. In this procedure, we did not consider the background since the individual spectra were already background subtracted. We computed the average count spectrum weighting the individual spectra by their respective exposure time.
We kept the same binning for the average spectrum as for the original data of each source. For a source observed several times, we first computed its average spectrum. The PDS energy band used in our study extends from 15 to 136 keV.

## 4 Spectral properties of various Seyfert 2 galaxies classes

We present in Figure 3 the average spectrum of Seyfert 2 galaxies with and without polarized broad lines.

The integrated signal to noise ratio in the 15-136 keV energy range is 98 for the PBLs and 47 for the NPBLs average spectrum. The NPBLs countrate is half as strong as the PBLs at 15 keV.
We compare the central engine signatures of the two types of Sy 2 galaxies using two methods: firstly, with a model-independent analysis comparing directly the normalized count spectra independently of any instrumental consideration; secondly, with a model-dependent method using a spectral fitting procedure.
We also intend to compare the emission of the two Sy 2 classes with the Sy 1 galaxies to investigate whether they have common properties as supposed by unified models. The Seyfert 1 average spectrum used in this study comes from D03.

### 4.1 Model-Independent Analysis: comparison of the normalized count spectra

We compare the average count spectra normalizing them at 15 keV (Figure 4).

#### 4.1.1 Comparison between Seyfert 2 galaxies with and without polarized broad lines

The average spectra of PBLs and NPBLs could be almost identical considering the error bars. Nevertheless, the PBLs emission seems to be slightly higher than that of NPBLs. To strengthen this result, we rebin the spectra to improve the bin significance (Figure 5).

The PBLs emission is higher than NPBLs along the major part of the energy domain indicating that the PBLs spectrum contains a larger hard photons population than NPBLs. The PBLs average spectrum is 45% higher than that of NPBLs around 85 keV where they differ most in the ratio of Figure 4.

#### 4.1.2 Comparison of PBLs and NPBLs Sy 2 galaxies with Sy 1 galaxies

1. Comparison of Sy 2 with PBL and Sy 1
The PBLs emission is always higher than that of Sy 1 galaxies giving an indication of the hardness of this category of Sy 2 galaxies. We note two excesses of the PBLs emission compared with Sy 1 around 55 keV and 100 keV. The deficit of hard photons in the Sy 1 class can be explained either by the presence of a cutoff in its emission or by a lower PBLs intrinsic spectral index and/or by a reflection component in PBLs. The PBLs emission is 85% greater than Sy 1 at $``$100 keV.
2. Comparison of Sy 2 without PBL and Sy 1
The two spectra are similar below 40 keV confirming that the emission of NPBLs is not significantly modified by absorption effects at these energies. The average spectra remain close up to $``$85 keV. The NPBLs emission is higher than that of Sy 1 at the two same energies (55 and 100 keV) where the PBLs emission dominates in the ratio PBLs/Sy 1. This suggests the presence of a particular spectral component in the Sy 1 average spectrum.
The maximal difference measured between the two classes is found at $``$100 keV where the NPBLs emission is 74% greater than that of Sy 1.

#### 4.1.3 Preliminary results

At first sight (Fig. 3), the spectra of Seyfert 2 galaxies with and without polarized broad lines appear very similar. However, the model-independent analysis we performed reveals small but significant differences (Fig. 5). The comparison with Sy 1 galaxies indicates that singularities exist also between Sy 1 galaxies and the two subclasses of Sy 2 galaxies.
The physical processes that can be responsible for the observed spectra and for the differences found between the various classes of Seyfert galaxies are investigated in the next subsection.

### 4.2 Model-Dependent Analysis: investigating the emission processes

To further test the above results, we compare the spectra using a fitting procedure that will allow us to investigate also the components occurring in the emission. We fit the spectra initially with simple models then with the reflection model PEXRAV. The results are shown in Tables 3 and 4.

#### 4.2.1 Fitting of the spectra with simple theoretical models

We fit the spectra successively with a simple power law and a cutoff power law (PL and CPL hereafter). To know the best model fitting the spectra, we refer to the reduced $`\chi `$2 value expected to be the nearest possible to 1.
A simple power law cannot reproduce the PBLs spectrum ($`\chi `$$`{}_{r}{}^{}{}_{}{}^{2}`$=3.58) thus at least one additional component is required. The occurrence of a cutoff in the PBLs observed spectrum is the best scheme. Indeed, the reduced $`\chi `$2 value ($`\chi `$$`{}_{r}{}^{}{}_{}{}^{2}`$=0.90) for a CPL and in particular the very low FTEST probability value (P$`_{\text{FTEST}}`$=7.9$``$10-5) indicate the high probability of having a cutoff in the PBLs spectrum.
The NPBLs average spectrum is well represented by a simple power law and a FTEST gives no justification for including a cutoff in the model.
Considering the best fit obtained for the two Sy 2 classes, NPBLs have a spectral index greater than PBLs, showing the softness of their emission ($`\mathrm{\Gamma }`$$`{}_{}{}^{\text{PL}}{}_{\text{NPBLs}}{}^{}`$=1.99$`{}_{}{}^{+0.08}{}_{0.07}{}^{}`$ against $`\mathrm{\Gamma }`$$`{}_{}{}^{\text{CPL}}{}_{\text{PBLs}}{}^{}`$=1.38$`{}_{}{}^{+0.14}{}_{0.15}{}^{}`$ with a cutoff energy of Ec=83.6$`{}_{}{}^{+35.6}{}_{20.4}{}^{}`$ keV).
PBLs have the same general behavior as Sy 1 galaxies, which also require a cutoff in their average spectrum. In particular the Sy 1 energy cutoff (Ec=70$`{}_{}{}^{+30}{}_{17}{}^{}`$ keV, see D03 for spectral parameters) is very similar to that of PBLs (Ec=83.6$`{}_{}{}^{+35.6}{}_{20.4}{}^{}`$ keV). The spectral index is possibly slightly steeper for Sy 1 ($`\mathrm{\Gamma }`$=1.49$`{}_{}{}^{+0.17}{}_{0.17}{}^{}`$) than for PBLs ($`\mathrm{\Gamma }`$=1.38$`{}_{}{}^{+0.14}{}_{0.15}{}^{}`$) but considering the errors they could be identical.
Therefore, Sy 1 and Sy 2 galaxies with polarized broad lines have spectral similarities whereas NPBLs galaxies present distinct properties with the absence of a constrained cutoff and a probable greater spectral index.

#### 4.2.2 Fitting of the spectra with the reflection model PEXRAV

To investigate if Compton reflection plays a role in the PBLs and NPBLs emission we use the PEXRAV model (Magdziarz & Zdziarski, 1995). This model calculates the expected X-ray spectrum when a point source of X-rays is incident on optically thick, mainly neutral material. The parameter R is related to the reflection component and the model is inclination/reflection dependent. The results are presented in Table 4.
We first fix the inclination from the normal of the disk at cos$`\theta `$=0.45 to compare the various spectral parameters of the two Sy 2 subclasses independently of any strong influence of the inclination. The PBLs average spectrum requires a cutoff whereas no indication of such a presence is observed in NPBLs as evoked in 4.2.1. The spectral indices of Sy 2 with and without PBL are within the same value domain. The presence of a reflection process in the NPBLs/PBLs emissions is not proven since the PEXRAV model does not improve significantly the fitting results. In particular, the reflection process could be totally absent in the PBLs spectrum (R=0.51$`{}_{}{}^{+1.48}{}_{0.51}{}^{}`$).
We now leave the inclination free to limit the dependence between the different spectral parameters. We find that PBLs would be close to the normal of the disk with an angle of cos$`\theta `$=0.81 or 36 whereas NPBLs would be seen with an angle of cos$`\theta `$=0.14 or 82.
The PBLs angle corresponds exactly to that found for Sy 1 galaxies in D03 where we obtained cos$`\theta `$=0.81. That does not prove that PBLs and Sy 1 have genuinely the same inclination but that their spectra give the same parameters when applying the PEXRAV model.

## 5 Comparison of the two different analyses

The models and parameters that best fit PBLs show that their properties are similar to those of Seyfert 1 galaxies (see 4.2.). However, the ratio of their average count spectra, shown in 4.1.2, indicates a significant difference at high energies. To understand this apparent contradiction between the two methods, we select the best fit model for each class and we calculate the ratio of these models (Figure 6).
We see in figure 6b that the PBLs model is higher than that of NPBLs over most of the energy domain by the same amount ($``$30$`\%`$) as seen in 4.1.1 in the count spectra ratio. The PBLs model is always higher than that of Sy 1 galaxies (Fig. 6c), that is explained by the Sy 1 cutoff occurring at a slightly lower energy than PBLs.
The divergence between NPBLs and Sy 1 galaxies does not exceed 20$`\%`$ below $``$85 keV whereas at high energies NPBLs are more powerful emitters (Fig. 6d). Indeed, the presence of a cutoff in the Sy 1 spectrum at $``$70 keV entails a rapid decrease of the Sy 1 emission whereas the NPBLs drop in emission is only due to the spectral index of the continuum.
The model ratios have thus the same general behavior and amplitude as those of the count spectra, confirming the results of the fits to the models.

## 6 Results

The spectra of our Seyfert 2 galaxies with and without PBL have a common general behavior at first sight. Investigating deeper, the PBLs spectrum presents a well-constrained cutoff whereas the NPBLs emission does not require this property.
Our results also indicate that Seyfert 2 galaxies with polarized broad lines have several common properties with Sy 1 galaxies (cutoff and hardness of the same order).
Therefore, the nature of Seyfert 2 galaxies without polarized broad lines remains unclear. We found significant divergences between NPBLs and the pair Sy1/PBLs. Indeed, NPBLs have no constrained cutoff, have a greater spectral index and are less luminous than PBLs, revealing intrinsic differences.

## 7 Discussion

Far from having definitively proven that NPBLs have a different central engine hence that they do not have a hidden Sy 1 nucleus, we show that their primary emission presents substantial differences compared to PBLs and Sy 1 galaxies.

### 7.1 The high-energy cutoff

Sy 1 and PBLs galaxies possess a cutoff of the same order, indicating that the temperature of the electron gas is similar for both classes (T$`_{\text{Sy1}}`$=Ec/k$`_\text{B}`$=8.27$``$108 K and T$`_{\text{PBLs}}`$=9.88$``$108 K). The presence of similar cutoffs in the PBLs and Sy 1 emission reveals that the same global physical conditions occur in the inner regions of their source. This supports unified models.
The absence of a constrained cutoff in NPBLs emphasizes differences in the primary emission of various classes of Sy 2 galaxies. However, a cutoff could arise at higher energies thus produced by very high electron gas temperatures. We therefore have to extend the energy domain of our investigation. This lack could also be explained by a wide distribution of NPBLs cutoff energies.
The genuine absence of a cutoff would suppose that different emission processes arise within the NPBLs class.
Nevertheless, the cutoff argument cannot explain the absence of polarized broad lines in the NPBLs sources. Another way to explain the existence of two Sy 2 types is by the study of the central engine activity, explored in the next subsection.

### 7.2 The central engine activity

We show in this study that the NPBLs countrate at 15 keV is two times weaker than PBLs (section 4), the NPBLs spectrum is softer than Sy 1 and PBLs (see 4.3.1 and 4.3.2) and in particular the NPBLs of our sample are less luminous (see 2.3.2). Therefore, the central engine of NPBLs appears weaker than PBLs.
We speculate that Sy 2 galaxies without polarized broad lines correspond to a state in which the accretion flow is weak and devoid of line emitting clouds. This flow may be insufficient to produce a large UV flux which in turn would explain why Sy 2 without PBL cannot cool through Compton processes between the electron population and the UV photons, thus explaining the absence of cutoff in the high energy emission.

### 7.3 The absorption/environment role in Seyfert 2 classification

The absence of polarized broad lines and of a cutoff in NPBLs galaxies cannot both be explained by very high column densities. Indeed, we consider only Compton-thin Sy 2 galaxies thus sources slightly influenced by absorption. The argument according to which the column density of our NPBLs would be underestimated is improbable here, since this would imply a large divergence between the average emission of Sy 1 and Sy 2 without PBL in the 15-40 keV domain, that is not observed.
The absence of polarized broad lines in NPBLs may be due to either the lack of a scattering medium above the absorbing gas or to the genuine absence of line emitting clouds in the nuclear regions.
However, the absence of a scattering medium alone cannot explain all the differences seen in the continuum emission of the two Sy 2 types.

### 7.4 Implications

We summarize in Table 5 the characteristics of Seyfert 2 galaxies with and without polarized broad lines.

The PBLs and Sy 1 galaxies of our sample have similar physical conditions in the inner regions of their sources. Consequently, the observed differences in the line properties and in X-ray emission can be due to a greater absorption in PBLs compared to Sy 1 galaxies, which agrees with unified models.
The NPBLs, however, are not only hidden Sy 1 nucleus since the absorption argument is unable to reproduce the different spectral properties found in this study. To confirm this assertion and to allow a meaningful analysis of the central engine signatures of the different types of Sy galaxies, we need a larger sample of Sy 2 having both spectropolarimetric observations and high quality X-ray data.
We conclude that while Seyfert 2 galaxies with polarized broad lines have continuum emission properties that match those of Seyfert 1 galaxies, Seyfert 2 galaxies without PBL do seem to present genuine differences in their nuclear properties.

###### Acknowledgements.

The author wishes to sincerely thank the anonymous referee for the fruitful report he provided, M. Gaber for his advice and P. Veron for the useful comments on the ”zoology” of Seyfert galaxies.

## Appendix A Average spectra fitted with the different models

We present the average count spectrum of each class fitted by the models applied in this study.


text:
# Using Multiwavelength Observations to Estimate the Black Hole Masses and Accretion Rates in Seyfert Galaxies

## 1. Introduction

In a series of three papers (Chiang & Blaes 2001, Chiang 2002, and Chiang & Blaes 2003), we showed that the optical/UV (OUV) and X-ray light curves obtained from multiwavelength monitoring campaigns of bright type 1 Seyfert galaxies could be readily explained by thermal reprocessing of radiation from a central X-ray source by a thin accretion disk. These observations had been puzzling since there is no apparent correlation between the OUV and X-ray light curves (Nandra et al. 1998; Edelson et al. 2000) as predicted by early models of thermal reprocessing (Krolik et al. 1991; Courvoisier & Clavel 1991; Collin-Souffrin 1991). The resolution of these difficulties were pointed to by Nandra et al. (2000) who noted that for the NGC 7469 observations the OUV variations are correlated with the X-ray spectral index as well as with an extrapolation of the measured 2–10 keV fluxes up to 100 keV. This latter result leads to the natural conclusion that the disk flux from thermal reprocessing is driven by the bolometric X-ray luminosity, which can differ substantially from that seen in narrow band observations (e.g., 2–10 keV) owing to variations in the X-ray spectral shape.

If the hard X-ray emission in type 1 Seyfert galaxies is produced by thermal Comptonization in a hot plasma, then the X-ray spectral index is determined by the plasma temperature and the optical depth to Compton scattering. Furthermore, if the thermal reprocessed emission comprises a significant fraction of the seed photons for the Comptonization process, then the feedback between the plasma and the reprocessing material acts as a thermostat, regulating the temperature of the plasma (Stern, et al. 1995). For a given optical depth, the spectral index of the emission is then set by the strength of the feedback, and that, in turn, is determined by the geometry of the hot plasma and the reprocessing material and by the amount of seed photons that are not due to reprocessing (Poutanen, Krolik, & Ryde 1997; Zdziarski, Lubiński, & Smith 1999).

## 2. The Model

There is an empirical relation between the Compton amplification factor, defined as the ratio of X-ray luminosity to seed photon luminosity entering the plasma, $`A=L_x/L_s`$, and the spectral index of the thermal Comptonization continuum:

$$\mathrm{\Gamma }=\mathrm{\Gamma }_0\left(\frac{L_x}{L_s}1\right)^{1/\delta }.$$
(1)
Here, parameter values of $`\mathrm{\Gamma }_0=2.06`$ and $`\delta =15.6`$ have been found from fits to Monte Carlo calculations of the thermal Comptonization process (Chiang & Blaes 2003; see also Malzac, Beloborodov, & Poutanen 2001, and Beloborodov 1999). Thus, the amplification factor $`A`$ is essentially determined for a given geometry and intrinsic seed photon distribution.111There is a slight complication in that the thermal Compton emission is not completely isotropic, being preferentially back-scattered for the lowest orders.

Since the shape of the OUV spectral energy distribution (SED) of type 1 Seyfert galaxies is roughly consistent with emission from a multitemperature thin disk, albeit with significant color corrections (Hubeny et al., 2001), and since the X-ray emitting region must have a fairly small size ($`<10^{14}\text{cm}𝒪(10GM/c^2)`$) from variability time scales, we model the geometry of the inner regions of these objects as consisting of a hot, centrally-located spherical plasma with a razor thin accretion disk encircling it. Some adjustable overlap exists between the disk and the central plasma that allows for variation in the amplification factor, and thus, will also allow for spectral variation. This geometry is motivated by the disk instability models of Shapiro, Lightman & Eardley (1976) as well as by more recent advection-dominated accretion flow (ADAF) models (e.g., Narayan & Yi 1995). Following these models, we will assume that the central spherical plasma region is the disk accretion flow having made a transition from a relatively cold, standard thin disk to a geometrically thick, hot inner flow.

For observations of the Seyfert galaxies NGC 5548 (Magdziarz et al. 1998), NGC 7469 (Nandra et al. 1998), and NGC 3516 (Edelson et al. 2000), the combinations of ground-based optical spectra and UV spectroscopy from IUE and/or HST provide sufficient spectral coverage to determine shape of the OUV SED, including, crucially, the roll-over at the higher frequencies. In our simple picture of the accretion disk, the emission at any given radius comprises thermal radiation from internal viscous dissipation and from thermal reprocessing. For our chosen “sphere + disk” geometry, both of these sources of emission obey a $`r^{3/4}`$ temperature distribution. Therefore, the inner truncation radius of the thin disk, $`r_{\mathrm{tr}}`$, where it makes the transition to the geometrically thick, hot central plasma region, determines the shape of the OUV roll-over, while the temperature at this radius sets the magnitude of the OUV SED. Once the disk temperature is set by the OUV observations, the source of seed photons from both thermal reprocessing and viscous dissipation is fully determined.

The temperature of the Comptonizing electrons, the optical depth and the radius of the central plasma sphere can then be found by fitting the X-ray flux and spectral index. For a single epoch of data, this leaves one parameter of this model unconstrained. If one assumes a Thomson depth of unity, then the full SED from the optical through hard X-rays can be described with reasonable parameters (Chiang & Blaes 2001). We note that this model would be completely constrained if the high energy roll-over of the thermal Comptonization continuum could be reliably measured. Such measurements would provide direct estimates of both the the temperature of the Comptonizing plasma and the overall bolometric luminosity. In analogy with Galactic low mass X-ray binaries, this roll-over is believed to occur at cut-off energies of $`100`$ keV. Unfortunately, there have not been any X-ray observing missions with the effective area required to measure the X-ray spectra at these energies on sufficiently short time scales.

Another constraint may be imposed that comes from the expectation that the hot central plasma must be less efficient at radiating accretion power than a standard optically thick disk would be in its place. This condition yields the relation

$$L_x<\eta \dot{M}c^2L_{\mathrm{disk}},$$
(2)
where $`\eta 0.1`$–0.4 is the expected efficiency for a thin disk that extends all the way to the innermost stable circular orbit, $`\dot{M}`$ is the accretion rate, and $`L_{\mathrm{disk}}=L_{\mathrm{disk}}(M\dot{M})`$ is the luminosity of the truncated disk from internal viscous dissipation. Although this condition adds a constraint it also adds yet another unknown quantity, the product $`M\dot{M}`$ (Chiang & Blaes 2003).

Multiwavelength observations of variability on time scales that are shorter than the expected viscous times for these systems ($`t_{\mathrm{visc}}>10^8`$s) can provide additional information that can be used to estimate the above model parameters uniquely. On these time scales, for typical type 1 Seyfert galaxies, we can take $`M\dot{M}`$ to be constant. Since there is one fewer constraint than there are model parameters, we can solve for one of the model parameters as a function of another for any given set of OUV and X-ray observations. In figure 2 of Chiang & Blaes (2003), we plot the plasma electron temperature, $`k_BT_e`$, versus $`M\dot{M}`$ for each of nine epochs of the OUV/X-ray data available for NGC 5548 (Magdziarz et al. 1998; Chiang & Blaes 2003). Assuming that a single value of $`M\dot{M}`$ must obtain for all nine epochs, we can infer an upper limit based on the restriction that the electron temperature cannot be smaller than $`30`$ keV, otherwise evidence of the thermal Comptonization cut-off would have been seen in the epoch 7 Ginga X-ray data. This constraint yields $`M\dot{M}<2.5\times 10^5M_{}^2`$yr-1. The dotted vertical lines indicate the black hole masses $`M=(1,2,3,4,5,6,7)\times 10^7M_{}`$ found for the epoch 1 data, which are the most constraining. Given that the black hole mass must be considered constant as well, this yields an upper limit of $`M<2\times 10^7M_{}`$. By contrast, the mass measured for NGC 5548 using reverberation mapping is $`M6\times 10^7M_{}`$ (Peterson & Wandel 2000).

We have applied similar analyses to data for NGC 3516 and NGC 7469 (Chiang 2002). Although we could not estimate the disk inner truncation radius for NGC 3516 since the OUV data do not show a high frequency decline, we can use the observations of a broad Fe K$`\alpha `$ line to identify the inner truncation radius of the thin disk with the innermost stable circular orbit. From our analyses of the multiwavelength data, we find a black hole mass of $`M<2\times 10^7M_{}`$. For NGC 7469, we find $`M\dot{M}10^6M_{}^2`$yr-1; this agrees with a value found from a completely independent method that considers interband OUV continuum time lags under the assumption of thermal reprocessing by a thin disk (Collier et al. 1998). Reverberation mapping observations place the central black hole mass of NGC 7469 in the range $`M10^6`$–10$`{}_{}{}^{7}M_{}^{}`$. Along with our estimate of $`r_{\mathrm{tr}}3\times 10^{14}`$cm for the thin disk inner truncation radius, this implies $`r_{\mathrm{tr}}>200GM/c^2`$ which is consistent with observations of a narrow, unresolved Fe K$`\alpha `$ line as seen by ASCA for this object (Guainazzi et al. 1994).

## References

Beloborodov, A. M. 1999, in ASP Conf. Ser. 161, High Energy Processes in Accreting Black Holes, ed. J. Poutanen & R. Svensson (San Francisco: ASP), 295

Chiang, J. 2002, ApJ, 572, 79

Chiang, J. & Blaes, O. 2001, ApJ, 557, L15

Chiang, J. & Blaes, O. 2003, ApJ, 586 97

Collier, S., et al. 1998, ApJ, 500, 162

Collin-Souffrin, S. 1991, A&A, 249, 344

Courvoisier, T. J.-L., & Clavel, J. 1991, A&A, 248, 389

Edelson, R. L., et al. 2000, ApJ, 534, 180

Guainazzi, M., Matsuoka, M., Piro, L., Mihara, T., & Yamauchi, M., 1994, ApJ, 436, L35

Hubeny, I., Blaes, O., Krolik, J. H., Agol, E. 2001, ApJ, 559, 680

Krolik, J. H., et al. 1991, ApJ, 371, 541

Magdziarz, P., et al. 1998, MNRAS, 301, 179

Malzac, J., Beloborodov, A. M., & Poutanen, J. 2001, MNRAS, 326, 417

Nandra, K., et al. 1998, ApJ, 505, 594

Nandra, K., et al. 2000, ApJ, 544, 734

Narayan, R., & Yi, I. 1995, ApJ, 452, 710

Peterson, B. M., & Wandel, A. 2000, ApJ, 540, L13

Poutanen, J., Krolik, J. H, & Ryde, F. 1997, MNRAS, 292, L21

Shapiro, S. L., Lightman, A. P., & Eardley, D. N. 1976, ApJ, 204, 187

Stern, B. E., et al. 1995, ApJ, 449, L13

Zdziarski, A. A., Lubiński, P., & Smith, D. A. 1999, MNRAS, 303, L11


text:
# Untitled Document

This paper has been withdrawn due to a flaw in the treatment of different spin components. It is replaced with arXiv:0810.2399 \[quant-ph\].


text:
# Photon trajectories in incoherent atomic radiation trapping as Lévy flights

## Abstract

Photon trajectories in incoherent radiation trapping for Doppler, Lorentz and Voigt line shapes under complete frequency redistribution are shown to be Lévy flights. The jump length ($`r`$) distributions display characteristic long tails. For the Lorentz line shape, the asymptotic form is a strict power-law $`r^{3/2}`$ while for Doppler the asymptotic is $`r^2\left(\mathrm{ln}r\right)^{1/2}`$. For the Voigt profile, the asymptotic form has always a Lorentz character, but the trajectory is a self-affine fractal with two characteristic Hausdorff scaling exponents.

Lévy flight, radiation trapping, radiation diffusion, truncated Lévy flight, stochastic processes.

Interest in distributions with long tails has increased over the last years, with intensive search for such laws in real physical systems. Most of the basic work on these distributions was carried out by Lévy in the 1930s Levy , but it was only relatively recently that these distributions were shown to be applicable to the description of a number of physical, biological and social phenomena Bou90 ; Shl95 ; Bar02 ; econometrics , including particle motion in turbulent media Bou90 ; Shl95 ; Shl99 , anomalous diffusion in microheterogeneous systems Ott1990 , chaotic transport in laminar fluid flow Swinney , the albatross flight Nat96  and frequency fluctuations of chromophores isolated in glassy environments Barkai . The purpose of this Letter is to show that photon trajectories in incoherent radiation trapping, a basic and common phenomenon in atomic and atmospheric physics, and in astrophysics, fit in the category of Lévy (superdiffusive) flights.

Although “the general trend nowadays is to put Lévy-type anomalous (super)diffusion on a similar footing with normal, Brownian-type, diffusion” Tsa1995 , the actual physical systems invoked to justify the practical importance of Lévy flights in physical space tend to be somewhat marginal in physical phenomena Shl95 . This contribution aims to prove that the well known case of classical incoherent radiation trapping is one of the simplest and best characterized Lévy flight found up to now. In addition, radiation trapping has a major role in fluorescent lamps and this could make radiation trapping the economically most relevant concretization of a Lévy flight. In retrospect, the 1990 Ott et al. Ott1990  contribution on the superdiffusion in microheterogeneous systems can no longer be considered the first experimental realization of a random flight with infinite moments, given that it was preceded by many studies of radiation trapping.

Radiation trapping of energy is important in areas as diverse as stellar atmospheres Mih , plasmas and atomic vapors luminescence Mol , terrestrial atmosphere and ocean optics, molecular luminescence MN , infrared radiative transfer and cold atoms Bar02 . In these optically thick media, the emitted radiation suffers several reabsorption and reemission events before eventually escaping to the exterior; the radiation is said to be imprisoned or trapped. Atomic radiation trapping is also known as imprisonment of resonance radiation, line transfer, radiation diffusion or multiple scattering of resonance radiation.

The first quantitative theory for atomic radiation trapping was presented in the 1920s by Compton and Milne, that have developed a modified diffusion equation for the (frequency) coherent spreading of excitation. It was only in 1932 that the frequency redistribution between absorption and reemission was taken into account by Kenty Ken32 . He considered a Doppler spectral distribution and arrived to the unexpected result that, for an infinite medium, the diffusion coefficient would be infinite. This was the first realization of the fundamental fact that all moments of the jump size distribution are infinite. Kenty’s result shows that a diffusion-type equation is not valid for radiation trapping with frequency redistribution effects. Nevertheless it was only in 1947 that Hosltein and Biberman independently proposed a Boltzmann-type integrodifferential equation H  which remains the starting point of the vast majority of radiation trapping models Mol .

Consider the case of inelastic scattering where, as the result of reabsorption-reemission events there is a photon frequency redistribution in the lab reference frame. The frequency distribution of the emitted photons is given by the emission spectrum $`\mathrm{\Theta }\left(x\right)`$. The absorption probability of a photon with frequency $`x`$ at a given distance from the emission point depends on the absorption spectrum $`\mathrm{\Phi }\left(x\right)`$, and is given by $`p\left(r|x\right)=\mathrm{\Phi }\left(x\right)\mathrm{e}^{\mathrm{\Phi }\left(x\right)r}`$ (Beer-Lambert law) where $`r`$ is the opacity or optical density and is a dimensionless distance. Radiation trapping can then be envisaged as a random flight in physical space with spectral shape dependent jump size distributions. The jump size distribution takes into account the absorption probability for all possible optical emission frequencies, hence

$$p\left(r\right)=_{\mathrm{}}^+\mathrm{}\mathrm{\Theta }\left(x\right)p\left(r|x\right)dx$$
(1)
Eq. (1) fully characterizes the spatial aspects of the random flight.

The moments of this distribution are

$$r^n=n!_{\mathrm{}}^+\mathrm{}\frac{\mathrm{\Theta }\left(x\right)}{\mathrm{\Phi }^n\left(x\right)}dx$$
(2)
and can be shown to be infinite for all physical reasonable atomic emission and absorption spectral distributions as concluded by Holstein H . His original analysis H  only included classical incoherent trapping, in which the emitting state is statistically unrelated with the absorbing one and there is therefore Complete Frequency Redistribution (CFR). In this case the absorption and emission spectra are identical and Eq. (1) reduces to

$$p\left(r\right)=_{\mathrm{}}^+\mathrm{}\mathrm{\Phi }^2\left(x\right)\mathrm{e}^{\mathrm{\Phi }\left(x\right)r}dx$$
(3)
Nevertheless, it can be concluded from Eq. (2) that the case of Partial Frequency Redistribution (PFR) is also characterized by an infinite moments Lévy statistics. The case of PFR is especially important in an astrophysical context. Neither case of complete coherent nor complete incoherent scattering is achieved exactly in stellar atmospheres and it is then necessary to consider the photon redistribution and to calculate redistribution functions which will give $`\mathrm{\Theta }\left(x\right)`$ Mih . In the usual laboratory conditions, vapor densities are high enough for CFR to apply Mol . This will be the case considered here. In two-level CFR atomic models both absorption and emission spectra can be described by Doppler – $`\mathrm{\Phi }_D\left(x\right)=\frac{1}{\sqrt{\pi }}\mathrm{e}^{x^2}`$ –, Lorentz – $`\mathrm{\Phi }_L\left(x\right)=\frac{1}{\pi }\frac{1}{1+x^2}`$ –, or Voigt – $`\mathrm{\Phi }_V\left(x\right)=\frac{a}{\pi ^{3/2}}_{\mathrm{}}^+\mathrm{}\frac{\mathrm{e}^{u^2}}{a^2+\left(xu\right)^2}du`$ – spectral distributions. $`x`$ is a normalized difference to the center of line frequency and $`a`$ is the Voigt characteristic width.

We now consider the asymptotic approximations valid for large jump sizes. A random flight in which the probability density of jump lengths is given by

$$p\left(r\right)\underset{r\mathrm{}}{}\frac{1}{r^{(1+\mu )}}$$
(4)
with $`\mu <2`$ is a self-similar random fractal with fractal dimension $`\mu `$. It is called a Lévy flight after Mandelbrot Man , and defines a broad distribution for which all the moments of order not smaller than $`\mu `$ are divergent.

If $`\mathrm{\Phi }\left(x\right)`$ is substituted for $`\mathrm{\Theta }\left(x\right)`$ in Eq. (2) above, it is found that $`r=\mathrm{}`$, whatever the spectral lineshape used, as long as $`\mathrm{\Phi }\left(x\right)`$ is nonzero for large $`\left|x\right|`$. In this way, $`\mu 1`$ for any $`\mathrm{\Phi }\left(x\right)`$.

In order to find the specific value of $`\mu `$ for the spectral distributions mentioned, we begin by rewriting Eq. (3) as

$$p(r)=\frac{\mathrm{d}^2J\left(r\right)}{\mathrm{d}r^2}$$
(5)
where

$$J\left(r\right)=_{\mathrm{}}^+\mathrm{}\left(1\mathrm{e}^{\mathrm{\Phi }\left(x\right)r}\right)dx$$
(6)
When the line shape is gaussian (Doppler), the integrand in Eq. (6) approaches a square wave form for large $`r`$, with inflection points at $`x_0`$ and $`x_0`$ with $`\mathrm{\Phi }\left(x_0\right)r1`$. Hence, $`x_0=\sqrt{\mathrm{ln}r}`$ for large $`r`$$`J\left(r\right)2x_0=2\sqrt{\mathrm{ln}r}`$ and therefore, from Eq. (5),

$$p\left(r\right)\underset{r\mathrm{}}{}\frac{1}{r^2\left(\mathrm{ln}r\right)^{1/2}}$$
(7)
Assuming an homogeneous scaling law one arrives at an effective $`\mu =1+1/2ln(ln(r))/ln(r)`$ which goes to $`1`$ for $`r\mathrm{}`$: although the asymptotic of Eq. (4) is only approximately valid, one can nevertheless classify Doppler trapping as a strict Lévy flight with $`\mu =1`$ with all the moments of the jump distribution being infinite.

When the line shape is Cauchy-like (Lorentz), the integrand in Eq. (6) can be simplified to $`1\mathrm{exp}\left(\mathrm{\Phi }\left(x\right)r\right)1\mathrm{exp}\left(\frac{r}{\pi x^2}\right)`$ for large $`r`$ and $`J\left(r\right)`$ becomes $`2\sqrt{r}`$. Therefore,

$$p\left(r\right)\underset{r\mathrm{}}{}\frac{1}{r^{3/2}}$$
(8)
and the asymptotic distribution is a Lévy flight with $`\mu =1/2`$.

The Voigt distribution is asymptotically coincident with the Lorentz distribution, and therefore has also $`\mu =1/2`$ for any value of its parameter $`a`$.

Eqs. (7) and (8) were already implied in Holsteins’s 1947 contribution (original equations are, up to some minor terms, integrated versions of the equations in this Letter), a fact hitherto unnoticed in the literature.

Fig. 1 show the jump probability obtained from direct numerical integration of Eq. (3). A linear fit gives for the characteristic “tail index” (Hausdorff box counting) fractal dimension ($`\mu `$) of the Doppler distribution $`\mu =1.07`$ (jump sizes $`10^3\mathrm{}10^4`$) and for the Lorentz case $`\mu =0.500`$ (jump sizes $`10^5\mathrm{}10^8`$), in agreement with the previous assertions. Fig. 1 also show that the continuous transformation from Doppler into Lorentz-like spectra as the Voigt $`a`$ width parameter changes from zero into infinite values does not manifests itself in a continuous change of the effective, $`r`$ dependent, $`\mu `$ values. There is an abrupt change instead at a jump distance which scales approximately as $`1/a`$. This is expected as the asymptotic expansion is related to the most extreme values of the wings of the spectral distribution and these change abruptly from Doppler to Lorentz-type asymptotics for $`a`$ values as low as $`10^4`$ (see Fig. 2). The data in Figs. 1 and 2 allows one to define the Voigt radiation trapping trajectories as a self-affine fractal with two different scaling exponents which manifest themselves at different lengthscales.

Consider now Fig. 3 which shows single excitation trajectories in a 3D infinite medium (compare Mandelbrot’s well known figure Man ). The Lorentz and Doppler cases display the two qualitative features characteristic of Lévy flights: (i) the longer pathlength jumps, although much less common, are of paramount importance to the overall spreading of excitation and (ii) self-similar behavior. There is an hierarchy of clusters formed at different lengthscales but with similar topology. The set of visited points constitute a fractal of characteristic dimension $`\mu `$ for the Doppler and Lorentz distributions Bou90 ; Shl95 ; Shl99 . The Voigt case shows its self-affine nature. The trajectory topology changes from a Lorentz character into a Doppler one for smaller distances. The overall topology of the trajectory as a whole is of course dictated by the higher Lorentz scaling.

Up to now we have considered frequency redistribution at each scattering event. However, for high opacity two-level systems there are presumably many elastic scattering events, before an inelastic scattering event occurs. It is therefore important to consider at least qualitatively the influence of elastic scattering. Elastic scattering events will fold up the excitation trajectories of Fig. 3 and the more the higher the ratio of elastic to inelastic scattering probabilities. From this argument alone we will expect that an increase in the elastic to inelastic probabilities ratio will lead to an increase in the fractal dimension, as it will approach the 3D brownian motion with fractal dimension $`2`$.

The discussion thus far has been restricted to the case of Lévy flights. For the vapors studied in the laboratory the time-of-flight of in-transit radiation is negligible compared to the waiting time between absorption and reemission, thus rendering the Lévy flights formalism especially adequate. From the set of visited points one can obtain spatial distribution functions of the excited density which fully characterize the spatial aspects of radiation migration. The temporal evolution can be factorized from the spatial part (Lévy flight) and then separately handled MN ; epereira . However, for interstellar gases the time between scattering events might be large compared to the absorption-reemission times. Because the speed of light is finite, a Lévy walk modification of the flights here presented will be more appropriate in an astrophysical context Shl99 .

The theoretical models for incoherent atomic radiation trapping are based either on the Holstein-Biberman multi-exponential mode expansion or on the so-called multiple scattering representation Mol . In the last case, the temporal evolution is known analytically and the spatial distributions can be separately obtained, either directly from the Eq. (1) or from its asymptotic approximation Eq. (4epereira . It is within this theoretical framework that the results reported in this Letter are most useful. On the other hand, the individual terms in the Holstein expansion have no direct relation with the one (or $`n^{th}`$) step jump probabilities. Although both approaches are ultimately equivalent MN ; Lai , we follow the multiple scattering approach, preferable for systems of low opacity MN . Also, the multiple scattering representation has a simple interpretation, since each term corresponds to a specific generation of excited atoms or molecules. In a large number of practical situations the opacity is not high enough to warrant the exclusive use of Holstein’s slowest exponential mode. On the other hand, the multi-exponential expansion remains valid, but its fundamental mode cannot be identified with Holstein’s high opacity result and should be estimated by the stationary mode associated with a non-changing spatial distribution function MN ; Lai , a point often misunderstood in the literature.

Results of Fig. 3 strictly apply to an infinite 3D medium but in a vast number of experimental situations the system is finite and the trajectory is truncated before the asymptotic Lévy expansion is able to manifest itself. Mantegna and Stanley Man94  introduced in 1994 a class of Lévy flights, the truncated Lévy flight (TLF), in which the largest steps of an ordinary Lévy flight are eliminated by a sharp cutoff in its power tail. This work allows a reassessment of Kenty’s pioneering contribution since he was the first to have used a truncation procedure: within the framework of the kinetic theory of gases he considered a truncated Maxwell distribution of speeds, with a maximum speed corresponding to a free path equal to the linear size of sample cell Ken32 . Although truncation of the jump size renders the moments of the distribution finite, the convergence to a Gaussian can be extremely slow and the random walk can exhibit anomalous behavior and multi-scaling properties in a wide range before convergence Man94 ; NakKop . Moreover, trajectory truncation is more complex to handle and does not create convergence to the Gaussian in usual situations.

In this Letter, it was shown that all photon trajectories arising from incoherent two-level complete frequency redistribution trapping are superdiffusive Lévy flights with $`\mu 1`$. In particular, for the Doppler lineshape $`\mu =1`$, whereas for Lorentz and Voigt $`\mu =1/2`$.

###### Acknowledgements.

E. Pereira would like to acknowledge fruitful discussions with François Bardou (CNRS, IPCMS). This work was supported by Fundação para a Ciência e Tecnologia (FCT, Portugal) within project POCTI/34836/FIS/2000.