text:
# Non-linear shock acceleration in the presence of seed particles

## 1 Introduction

Diffusive acceleration at newtonian shock fronts is an extensively studied phenomenon. Detailed discussions of the current status of the investigations can be found in some recent excellent reviews . While much is by now well understood, some issues are still subjects of much debate, for the theoretical and phenomenological implications that they may have. The most important of these is the backreaction of the accelerated particles on the shock: the violation of the test particle approximation occurs when the acceleration process becomes sufficiently efficient to generate pressures of the accelerated particles which are comparable with the incoming gas kinetic pressure. Both the spectrum of the particles and the structure of the shock are changed by this phenomenon, which is therefore intrinsically nonlinear.

At present there are three viable approaches to account for the backreaction of the particles upon the shock: one is based on the ever-improving numerical simulations  that allow a self-consistent treatment of several effects. The second approach is a fluid approach, and treats cosmic rays as a relativistic second fluid. This class of models was proposed and discussed in . These models allow one to obtain the thermodynamics of the modified shocks, but do not provide information about the spectrum of accelerated particles.

The third approach is analytical and may be very helpful to understand the physics of the nonlinear effects in a way that sometimes is difficult to achieve through simulations, due to their intrinsic complexity.

In Ref.  a perturbative approach was adopted, in which the pressure of accelerated particles was treated as a small perturbation. By construction this method provides an answer only for weakly modified shocks.

An alternative approach was proposed in , based on the assumption that the diffusion of the particles is sufficiently energy dependent that different parts of the fluid are affected by particles with different average energies. The way the calculations are carried out implies a sort of separate solution of the transport equation for subrelativistic and relativistic particles, so that the two spectra must be somehow connected at $`pmc`$ a posteriori.

Recently, in , the effects of the non-linear backreaction of accelerated particles on the maximum achievable energy were investigated, together with the effects of geometry. The maximum energy of the particles accelerated in supernova remnants in the presence of large acceleration efficiencies was also studied in .

The need for a practical solution of the acceleration problem in the non-linear regime was recognized in , where a simple analytical broken-power-law approximation of the non-linear spectra was presented.

Recently, some promising analytical solutions of the problem of non-linear shock acceleration have appeared in the literature . These solutions seem to avoid many of the limitations of previous approaches.

Numerical simulations have been instrumental to identify the dramatic effects of the particles backreaction: they showed that even when the fraction of particles injected from the thermal gas is relatively small, the energy channelled into these few particles can be an appreciable part of the kinetic energy of the unshocked fluid, making the test particle approach unsuitable. The most visible effects of the backreaction of the accelerated particles on the shock appear in the spectrum of the accelerated particles, which shows a peculiar flattening at the highest energies. The analytical approaches reproduce well the basic features arising from nonlinear effects in shock acceleration.

While several calculations exist of the nonlinear effects in the shock acceleration of quasi-monochromatic particles injected at a shock surface, there is no description at present of how these effects appear, if they do, when the shock propagates in a medium where (pre)accelerated particles already exist. The linear theory of this phenomenon was developed by Bell , but has never been generalized to its nonlinear extension. In fact, the backreaction can severely affect the process of re-energization of preaccelerated particles: Bell already showed that for strong shocks the energy content of a region where cosmic rays were present could be easily enhanced by a factor $`10`$ at each shock passage, so that equipartition could be readily reached. In these conditions the backreaction of the accelerated particles should be expected.

We report here on the first analytical treatment of the shock acceleration in the presence of seed nonthermal particles, with the inclusion of their nonlinear backreaction on the shock. Our approach is a generalization of the analytical method introduced in  to describe the nonlinear shock acceleration with monochromatic injection of quasi-thermal particles. In fact, we present here also a general calculation that accounts for both thermal particles and seed nonthermal particles. This case may be of interest for the study of supernova shocks propagating through the interstellar medium (ISM) where pressure balance exists between gas and cosmic rays.

Nonlinear effects in shock acceleration of thermal particles result in the appearance of multiple solutions in certain regions of the parameter space. This behaviour resembles that of critical systems, with a bifurcation occurring when some threshold is reached. In the case of shock acceleration, it is not easy to find a way of discriminating among the multiple solutions when they appear. Neverthless, in , a two fluid approach has been used to demonstrate that when three solutions appear, the one with intermediate efficiency for particle acceleration is unstable to corrugations in the shock structure and emission of acustic waves. Plausibility arguments may be put forward to justify that the system made of the shock plus the accelerated particles may sit at the critical point, but the author is not aware of any real proof that this is what happens. The physical parameters that play a role in this approach to criticality are the maximum momentum achievable by the particles in the acceleration process, the Mach number of the shock, and the injection efficiency, namely the fraction of thermal particles crossing the shock that are accelerated to nonthermal energies. The last of them, the injection efficiency, hides a crucial physics problem by itself, and may play an important role in establishing the level of shock modification. This efficiency parameter in reality is defined by the microphysics of the shock and should not be a free parameter of the problem. Unfortunately, our poor knowledge of such microphysics, in particular for collisionless shocks, does not allow us to establish a clear and universal connection between the injection efficiency and the macroscopic shock properties.

The paper is structured as follows: in §2 we describe the effect of non linearity on shock acceleration and our mathematical approach to describe it. In particular we generalize previous calculations to the case in which seed particles exist in the region where the shock is propagating; in §3 we describe the gas dynamics in the presence of a non-negligible pressure of accelerated particles. In §4 we describe our results, with particular attention for the onset of the particle backreaction and for the appearance of multiple solutions. We conclude in §5.

## 2 Non linear shock acceleration

In this section we solve the diffusion-convection equation for the cosmic rays in the most general case in which particles are injected according to some function $`Q(x,p)`$ and the shock propagates within a region where some cosmic ray distribution exists, that we denote as $`f_{\mathrm{}}(p)`$, spatially homogeneous before the shock crossing (upstream).

For simplicity we limit ourselves to the case of one-dimensional shocks, but the introduction of different geometrical effects is relatively simple, and in fact many of our conclusions should not be affected by geometry.

The equation that describes the diffusive transport of particles in one dimension is

$$\frac{}{x}\left[D\frac{}{x}f(x,p)\right]u\frac{f(x,p)}{x}+\frac{1}{3}\frac{du}{dx}p\frac{f(x,p)}{p}+Q(x,p)=0,$$
(1)
where we assumed stationarity ($`f/t=0`$). The $`x`$ axis is oriented from upstream to downstream, as in fig. 1. The presence of pre-existing cosmic rays is introduced as a boundary condition at upstream infinity, by requiring that $`f(x=\mathrm{},p)=f_{\mathrm{}}(p)`$. One should keep in mind that the common picture of a fluid whose speed is constant ($`u_1`$) until it hits the shock surface is not appropriate for modified shocks. In fact, in this case, the pressure of the accelerated particles may become large enough to slow down the fluid before it crosses the shock surface. Therefore in general at upstream infinity the gas flows at speed $`u_0`$, different from $`u_1`$ (fluid speed immediately upstream of the shock). The two quantities are approximately equal only when the accelerated particles do not have dynamical relevance. The injection term is taken in the form $`Q(x,p)=Q_0(p)\delta (x)`$.

As a first step, we integrate eq. 1 around $`x=0`$, from $`x=0^{}`$ to $`x=0^+`$, denoted in fig. 1 as points “1” and “2” respectively, so that the following equation can be written:

$$\left[D\frac{f}{x}\right]_2\left[D\frac{f}{x}\right]_1+\frac{1}{3}p\frac{df_0}{dp}(u_2u_1)+Q_0(p)=0,$$
(2)
where $`u_1`$ ($`u_2`$) is the fluid speed immediately upstream (downstream) of the shock and $`f_0`$ is the particle distribution function at the shock location. By requiring that the distribution function downstream is independent of the spatial coordinate (homogeneity), we obtain $`\left[D\frac{f}{x}\right]_2=0`$, so that the boundary condition at the shock can be rewritten as

$$\left[D\frac{f}{x}\right]_1=\frac{1}{3}p\frac{df_0}{dp}(u_2u_1)+Q_0(p).$$
(3)
We can now perform the integration of eq. (1) from $`x=\mathrm{}`$ to $`x=0^{}`$ (point “1”), in order to take into account the boundary condition at upstream infinity. Using eq. (3) we obtain

$$\frac{1}{3}p\frac{df_0}{dp}(u_2u_1)u_1f_0+u_0f_{\mathrm{}}+Q_0(p)+_{\mathrm{}}^0^{}𝑑xf\frac{du}{dx}+\frac{1}{3}_{\mathrm{}}^0^{}𝑑x\frac{du}{dx}p\frac{f}{p}=0.$$
(4)
We can now introduce a quantity $`u_p`$ defined as

$$u_p=u_1\frac{1}{f_0}_{\mathrm{}}^0^{}𝑑x\frac{du}{dx}f(x,p),$$
(5)
whose physical meaning is instrumental to understand the nonlinear backreaction of particles. The function $`u_p`$ is the average fluid velocity experienced by particles with momentum $`p`$ while diffusing upstream away from the shock surface. In other words, the effect of the average is that, instead of a constant speed $`u_1`$ upstream, a particle with momentum $`p`$ experiences a spatially variable speed, due to the pressure of the accelerated particles that is slowing down the fluid. Since the diffusion coefficient is in general $`p`$-dependent, particles with different energies feel a different compression coefficient, higher at higher energies if, as expected, the diffusion coefficient is an increasing function of momentum.

The role of $`u_p`$ can also be explained as follows: the distribution function $`f(x,p)`$ at a distance $`x`$ from the shock surface can be written as

$$f(x,p)=f_0(p)\mathrm{exp}\left[\frac{q(p)}{3D(p)}_0^x𝑑x^{}u(x^{})\right],$$

where $`q(p)=d\mathrm{ln}f_0(p)/d\mathrm{ln}p`$ is the local slope of $`f_0(p)`$ and the diffusion coefficient $`D(p)`$ has been assumed independent of the location $`x`$. In first approximation, we can assume that the exponential factor remains important when it is of order unity, namely when its argument is much less than unity. We can therefore introduce a distance $`x_p`$, which is the distance at which the exponential equals one. This means that

$$u_pu_1_{x_p}^0𝑑x\frac{du}{dx}u_1\left[u_1u(x_p)\right]u(x_p),$$

so that the speed $`u_p`$ can be interpreted as the fluid speed at the point $`x_p`$ where the particles with momentum $`p`$ reverse their motion in the upstream fluid and return to the shock. Note that the diffusion coefficient enters the calculation of the distance $`x_p`$ but does not enter directly the calculation of $`u_p`$. In other words, different diffusion coefficients may move the point where the fluid speed is $`u_p`$ closer to or farther from the shock surface, but do not affect the value of $`u_p`$. This approach is similar to that introduced in .

With the introduction of $`u_p`$, eq. (4) becomes

$$\frac{1}{3}p\frac{df_0}{dp}(u_2u_p)f_0\left[u_p+\frac{1}{3}p\frac{du_p}{dp}\right]+u_0f_{\mathrm{}}+Q_0(p)=0,$$
(6)
where we used the fact that

$$p\frac{d}{dp}_{\mathrm{}}^0^{}𝑑x\frac{du}{dx}f=p\left[\frac{df_0}{dp}(u_1u_p)f_0\frac{du_p}{dp}\right].$$

Eq. (6) can be written in a way that resembles more the equation for shocks with no particles backreaction but in the presence of seed particles with distribution $`f_{\mathrm{}}(p)`$ and injection function $`Q_0(p)`$:

$$p\frac{df_0}{dp}=\frac{3}{u_pu_2}\left\{f_0\left[u_p+\frac{1}{3}p\frac{du_p}{dp}\right]u_0f_{\mathrm{}}Q_0(p)\right\}.$$
(7)
The solution of this equation can be written in the following implicit form:

$$f_0(p)=f_0^{reacc}(p)+f_0^{inj}(p)=$$

$$_{p_0}^p\frac{d\overline{p}}{\overline{p}}\frac{3\left[u_0f_{\mathrm{}}(\overline{p})+Q_0(\overline{p})\right]}{u_{\overline{p}}u_2}\mathrm{exp}\left\{_{\overline{p}}^p\frac{dp^{}}{p^{}}\frac{3}{u_p^{}u_2}\left[u_p^{}+\frac{1}{3}p^{}\frac{du_p^{}}{dp^{}}\right]\right\}.$$
(8)
In the case of monochromatic injection with momentum $`p_{inj}`$ at the shock surface, we can write

$$Q_0(p)=\frac{\eta n_{gas,1}u_1}{4\pi p_{inj}^2}\delta (pp_{inj}),$$
(9)
where $`n_{gas,1}`$ is the gas density immediately upstream ($`x=0^{}`$) and $`\eta `$ parametrizes the fraction of the particles crossing the shock which are going to take part in the acceleration process. The injection term in eq. (8) becomes

$$f_0^{inj}(p)=_{p_0}^p\frac{d\overline{p}}{\overline{p}}\frac{3Q_0(\overline{p})}{u_{\overline{p}}u_2}\mathrm{exp}\left\{_{\overline{p}}^p\frac{dp^{}}{p^{}}\frac{3}{u_p^{}u_2}\left[u_p^{}+\frac{1}{3}p^{}\frac{du_p^{}}{dp^{}}\right]\right\}=$$

$$\frac{3R_{sub}}{R_{sub}1}\frac{\eta n_{gas,1}}{4\pi p_{inj}^3}\mathrm{exp}\left\{_{p_{inj}}^p\frac{dp^{}}{p^{}}\frac{3}{u_p^{}u_2}\left[u_p^{}+\frac{1}{3}p^{}\frac{du_p^{}}{dp^{}}\right]\right\}.$$
(10)
Here we introduced the two quantities $`R_{sub}=u_1/u_2`$ and $`R_{tot}=u_0/u_2`$, which are respectively the compression factor at the gas subshock and the total compression factor between upstream infinity and downstream. The two compression factors would be equal in the test particle approximation. For a modified shock, $`R_{tot}`$ can attain values much larger than $`R_{sub}`$ and more in general, much larger than $`4`$, which is the maximum value achievable for an ordinary strong non-relativistic shock. The increase of the total compression factor compared with the prediction for an ordinary shock is responsible for the peculiar flattening of the spectra of accelerated particles that represents a feature of nonlinear effects in shock acceleration. In terms of $`R_{sub}`$ and $`R_{tot}`$ the density immediately upstream is $`n_{gas,1}=(\rho _0/m_p)R_{tot}/R_{sub}`$.

In eq. (10) we can introduce a dimensionless quantity $`U(p)=u_p/u_0`$ so that

$$f_0^{inj}(p)=\left(\frac{3R_{sub}}{R_{tot}U(p)1}\right)\frac{\eta n_{gas,1}}{4\pi p_{inj}^3}\mathrm{exp}\left\{_{p_{inj}}^p\frac{dp^{}}{p^{}}\frac{3R_{tot}U(p^{})}{R_{tot}U(p^{})1}\right\}.$$
(11)
Introducing the same formalism also for the reacceleration term in eq. (8), we obtain the general expression

$$f_0(p)=\frac{3R_{tot}}{R_{tot}U(p)1}_{p_0}^p\frac{d\overline{p}}{\overline{p}}f_{\mathrm{}}(\overline{p})\mathrm{exp}\left\{_{\overline{p}}^p\frac{dp^{}}{p^{}}\frac{3R_{tot}U(p^{})}{R_{tot}U(p^{})1}\right\}+$$

$$\left(\frac{3R_{sub}}{R_{tot}U(p)1}\right)\frac{\eta n_{gas,1}}{4\pi p_{inj}^3}\mathrm{exp}\left\{_{p_{inj}}^p\frac{dp^{}}{p^{}}\frac{3R_{tot}U(p^{})}{R_{tot}U(p^{})1}\right\}.$$
(12)
The solution of the problem is known if the velocity field $`U(p)=u_p/u_0`$ is known. The nonlinearity of the problem reflects in the fact that $`U(p)`$ is in turn a function of $`f_0`$ as it is clear from the definition of $`u_p`$. In order to solve the problem we need to write the equations for the thermodynamics of the system including the gas, the reaccelerated cosmic rays, the cosmic rays accelerated from the thermal pool and the shock itself. We write and solve these equations in the next section.

## 3 The gas dynamics of modified shocks with seed particles

The velocity, density and thermodynamic properties of the fluid can be determined by the mass and momentum conservation equations, with the inclusion of the pressure of the accelerated particles and of the preexisting cosmic rays. We write these equations between a point far upstream ($`x=\mathrm{}`$), where the fluid velocity is $`u_0`$ and the density is $`\rho _0=mn_{gas,0}`$, and the point where the fluid upstream velocity is $`u_p`$ (density $`\rho _p`$). The index $`p`$ denotes quantities measured at the point where the fluid velocity is $`u_p`$, namely at the point $`x_p`$ that can be reached only by particles with momentum $`p`$.

The mass conservation implies:

$$\rho _0u_0=\rho _pu_p.$$
(13)
Conservation of momentum reads:

$$\rho _0u_0^2+P_{g,0}+P_{CR,0}=\rho _pu_p^2+P_{g,p}+P_{CR,p},$$
(14)
where $`P_{g,0}`$ and $`P_{g,p}`$ are the gas pressures at the points $`x=\mathrm{}`$ and $`x=x_p`$ respectively, and $`P_{CR,p}`$ is the pressure in accelerated particles at the point $`x_p`$ (we used the symbol $`CR`$ to mean cosmic rays, in the sense of accelerated particles). The mass flow in accelerated particles has reasonably been neglected.

Our basic assumption, similar to that used in , is that the diffusion is $`p`$-dependent and more specifically that the diffusion coefficient $`D(p)`$ is an increasing function of $`p`$. Therefore the typical distance that a particle with momentum $`p`$ moves away from the shock is approximately $`\mathrm{\Delta }xD(p)/u_p`$, larger for high energy particles than for lower energy particles111For the cases of interest, $`D(p)`$ increases with $`p`$ faster than $`u_p`$ does, therefore $`\mathrm{\Delta }x`$ is a monotonically increasing function of $`p`$.. As a consequence, at each given point $`x_p`$ only particles with momentum larger than $`p`$ are able to affect appreciably the fluid. Strictly speaking the validity of the assumption depends on how strongly the diffusion coefficient depends on the momentum $`p`$.

The cosmic ray pressure at $`x_p`$ is the sum of two terms: one is the pressure contributed by the adiabatic compression of the cosmic rays at upstream infinity, and the second is the pressure of the particles accelerated or reaccelerated at the shock ($`\stackrel{~}{P}_{CR}(p)`$) and able to reach the position $`x_p`$. Since only particles with momentum $`\stackrel{>}{}p`$ can reach the point $`x=x_p`$, we can write

$$P_{CR,p}=P_{CR,0}\left(\frac{u_0}{u_p}\right)^{\gamma _{CR}}+\stackrel{~}{P}_{CR}(p)$$

$$P_{CR,0}(p)\left(\frac{u_0}{u_p}\right)^{\gamma _{CR}}+\frac{4\pi }{3}_p^{p_{max}}𝑑pp^3v(p)f_0(p),$$
(15)
where $`v(p)`$ is the velocity of particles with momentum $`p`$, $`p_{max}`$ is the maximum momentum achievable in the specific situation under investigation, and $`\gamma _{CR}`$ is the adiabatic index for the accelerated particles.

Let us consider separately the case of a strongly modified and weakly modified shock, in order to determine the best choice for $`\gamma _{CR}`$. In the case of strongly modified shocks, we will show that most energy is piled up in the region $`pp_{max}m`$, therefore in this case we can safely adopt $`\gamma _{CR}=4/3`$, appropriate for a relativistic gas. For weakly modified shocks, the accelerated particles have an approximately power law spectrum with a slope $`\alpha `$. It can be shown that in this case $`\gamma _{CR}=\alpha /3`$, so that the relativistic result $`\gamma _{CR}=4/3`$ still applies for $`\alpha =4`$ (strong shocks). For steeper spectra ($`\alpha >4`$) a larger adiabatic index should be adopted, but in those cases the solution is basically independent of the choice of $`\gamma _{CR}`$ because of the weak backreaction of the particles. For the purpose of carrying out our numerical calculation we will therefore always take $`\gamma _{CR}=4/3`$.

The pressure of cosmic rays at upstream infinity is simply

$$P_{CR,0}=\frac{4\pi }{3}_{p_{min}}^{p_{max}}𝑑pp^3v(p)f_{\mathrm{}}(p),$$
(16)
where $`p_{min}`$ is some minimum momentum in the spectrum of seed particles. For simplicity, we assume that $`p_{min}=p_{inj}`$, namely the minimum momentum of the seed particles coincides with the momentum at which particles are injected in the shock and are accelerated.

From eq. (14) we can see that there is a maximum distance, corresponding to the propagation of particles with momentum $`p_{max}`$ such that at larger distances the fluid is unaffected by the accelerated particles and $`u_p=u_0`$.

We will show later that for strongly modified shocks the integral in eq. (15) is dominated by the region $`pp_{max}`$. This improves even more the validity of our approximation $`P_{CR,p}=P_{CR}(>p)`$. This also suggests that different choices for the diffusion coefficient $`D(p)`$ may affect the value of $`p_{max}`$, but at fixed $`p_{max}`$ the spectra of the accelerated particles should not change in a significant way.

Assuming an adiabatic compression of the gas in the upstream region, we can write

$$P_{g,p}=P_{g,0}\left(\frac{\rho _p}{\rho _0}\right)^{\gamma _g}=P_{g,0}\left(\frac{u_0}{u_p}\right)^{\gamma _g},$$
(17)
where we used mass conservation, eq. (13). The gas pressure far upstream is $`P_{g,0}=\rho _0u_0^2/(\gamma _gM_0^2)`$, where $`\gamma _g`$ is the ratio of specific heats for the gas ($`\gamma _g=5/3`$ for an ideal gas) and $`M_0`$ is the Mach number of the fluid far upstream.

We introduce now a parameter $`\xi _{CR}`$ that quantifies the relative weight of the cosmic ray pressure at upstream infinity compared with the pressure of the gas at the same location, $`\xi _{CR}=P_{CR,0}/P_{g,0}`$. Using this parameter and the definition of the function $`U(p)`$, the equation for momentum conservation becomes

$$\frac{dU}{dp}\left[1\frac{\gamma _{CR}}{\gamma _g}\frac{\xi _{CR}}{M_0^2}U^{(\gamma _{CR}+1)}\frac{1}{M_0^2}U^{(\gamma _g+1)}\right]+\frac{1}{\rho _0u_0^2}\frac{d\stackrel{~}{P}_{CR}}{dp}=0.$$
(18)
Using the definition of $`\stackrel{~}{P}_{CR}`$ and multiplying by $`p`$, this equation becomes

$$p\frac{dU}{dp}\left[1\frac{\gamma _{CR}}{\gamma _g}\frac{\xi _{CR}}{M_0^2}U^{(\gamma _{CR}+1)}\frac{1}{M_0^2}U^{(\gamma _g+1)}\right]=\frac{4\pi }{3\rho _0u_0^2}p^4v(p)f_0(p),$$
(19)
where $`f_0`$ depends on $`U(p)`$ as written in eq. (12). Eq. (19) is therefore an integral-differential nonlinear equation for $`U(p)`$. The solution of this equation also provides the spectrum of the accelerated particles.

The last missing piece is the connection between $`R_{sub}`$ and $`R_{tot}`$, the two compression factors appearing in eq. (8). The compression factor at the gas shock around $`x=0`$ can be written in terms of the Mach number $`M_1`$ of the gas immediately upstream through the well known expression

$$R_{sub}=\frac{(\gamma _g+1)M_1^2}{(\gamma _g1)M_1^2+2}.$$
(20)
On the other hand, if the upstream gas evolution is adiabatic, then the Mach number at $`x=0^{}`$ can be written in terms of the Mach number of the fluid at upstream infinity $`M_0`$ as

$$M_1^2=M_0^2\left(\frac{u_1}{u_0}\right)^{\gamma _g+1}=M_0^2\left(\frac{R_{sub}}{R_{tot}}\right)^{\gamma _g+1},$$

so that from the expression for $`R_{sub}`$ we obtain

$$R_{tot}=M_0^{\frac{2}{\gamma _g+1}}\left[\frac{(\gamma _g+1)R_{sub}^{\gamma _g}(\gamma _g1)R_{sub}^{\gamma _g+1}}{2}\right]^{\frac{1}{\gamma _g+1}}.$$
(21)
Now that an expression between $`R_{sub}`$ and $`R_{tot}`$ has been found, eq. (19) basically is an equation for $`R_{sub}`$, with the boundary condition that $`U(p_{max})=1`$. Finding the value of $`R_{sub}`$ (and the corresponding value for $`R_{tot}`$) such that $`U(p_{max})=1`$ also provides the whole function $`U(p)`$ and, through eq. (8), the distribution function $`f_0(p)`$ for the particles resulting from acceleration and reacceleration in the nonlinear regime. If the backreaction of the accelerated particles is small, the test particle solution must be recovered.

## 4 The nonlinearity of acceleration and reacceleration

In this section we investigate the shock modification due to the backreaction of the accelerated particles. We split the discussion in two parts: in §4.1 we consider the case of pre-existing seed particles populating the region where the shock is propagating, and re-energized by the shock. We find the conditions for the onset of the nonlinear regime in which the shock gets modified by the re-accelerated particles. In particular we show that the critical behaviour found for shock acceleration and manifesting itself through the appearance of three solutions, takes place also in this case. In §4.2 we explore the more complicated situation in which a shock accelerates a new population of particles while possibly reaccelerating a pre-existing population of seed particles.

### 4.1 Shocks modified by re-energized seed particles

Here we assume that a fluid is moving with speed $`u_0`$ in a region where the temperature is $`T_0`$. The Mach number is some pre-defined value $`M_0`$, which can be easily related to $`u_0`$ and $`T_0`$. In the region of interest we assume that a population of seed particles is present with energy per particle already higher than some injection energy necessary for the particles to feel the shock as a discontinuity. For simplicity we assume that these seed particles have a spectrum which is a power law in momentum in the form:

$$f_{\mathrm{}}(p)=f_{\mathrm{},0}p^\alpha .$$
(22)
In this section we explore the situation in which the shock reacceleration of a pre-existing population of cosmic rays can modify the shock structure, but there is no acceleration of particles different from those that are already present as seed particles. This is equivalent to turn off the injection term ($`\eta =0`$ in Eq. (12)). In the next section we discuss what happens when both components are present. The crucial difference between the two components is that the seed particles are by definition already above the threshold for injection, so that there is no injection efficiency that instead represents such a crucial parameter for the case of acceleration of particles extracted from the thermal pool.

Having in mind the case of shocks propagating in the interstellar medium of our Galaxy, we consider here the case in which the inflowing gas and the seed particles are in pressure equilibrium, namely $`P_{g,0}=P_{CR}`$. In terms of the parameter $`\xi _{CR}=P_{CR}/P_{g,0}`$ this implies $`\xi _{CR}=1`$. The results can then be easily repeated for a generic value of $`\xi _{CR}`$.

Eq. (19) gives the quantity $`U(p)`$ as a function of the momentum $`p`$ for any choice of $`R_{sub}`$ and $`R_{tot}`$, while the acceptable solutions are those with the right matching conditions at $`p_{max}`$. For simplicity, let us assume that the $`p_{max}`$ in eq. (22) is the same maximum momentum that particles injected at the shock would achieve: this value only depends on the environmental conditions (energy losses of the particles) and/or on the geometry of the shock, which may allow the escape of the particles. The momentum $`p_{max}`$ is the same as in the distribution of the pre-exisiting seed particles if, for instance, the seed particles have been accelerated by a shock identical to the one we are considering. In any case, this assumption is not needed for the validity of our conclusions and may be easily relaxed, it simply serves to avoid parameter proliferation. Within this assumption, particles re-energized at the shock are simply redistributed in the momentum range between $`p_{inj}`$ and $`p_{max}`$.

The solution, namely the right pair of values for the compression parameters $`R_{sub}`$ and $`R_{tot}`$ \[related through eq. (21)\], is obtained when the solution corresponding to $`U(p_{max})=1`$ is selected.

In Fig. 2, we plot $`U(p_{max})`$ as a function of $`R_{tot}`$ for a shock having Mach number at infinity $`M_0=150`$. The different lines are labelled by a number representing the $`\mathrm{log}_{10}`$ of $`p_{max}`$ in units of $`mc`$. In other words the parameter $`p_{max}`$ changes between $`10^2mc`$ and $`10^9mc`$. The physical solutions are found by determining the intersections of each curve with the horizontal line corresponding to $`U(p_{max})=1`$. The minimum momentum of the seed particles is taken as $`10^3mc`$.

The intersection at $`R_{tot}4`$ is close to the well known linear solution, namely the solution that one would obtain in the test particle approximation. The energy channelled into the nonthermal particles for this solution is very small, and the shock remains approximately unmodified. Increasing the value of $`p_{max}`$, the solution moves toward slightly larger values of $`R_{tot}`$ (namely the shock becomes more modified). For some values of $`p_{max}`$, multiple solutions appear. In particular three regions of $`p_{max}`$ can be identified:

1) $`p_{max}p_{cr}^{(1)}`$: in this region (low values of $`p_{max}`$) the only solution is very close to the one obtained in the test particle approximation. For the values of the parameters in Fig. 2, $`p_{cr}^{(1)}10^3mc`$.

2) $`p_{cr}^{(1)}p_{max}p_{cr}^{(2)}`$: in this region (intermediate values of $`p_{max}`$ three solutions appear, two of which imply a strong modification of the shock, namely an appreciable part of the energy flowing through the shock is converted into energy of the accelerated particles. For the values of the parameters in Fig. 2, $`p_{cr}^{(2)}10^9mc`$.

3) $`p_{max}p_{cr}^{(2)}`$: in this region (high values of $`p_{max}`$) the solution becomes one again, and the shock is always strongly modified ($`R_{tot}4`$).

This critical behaviour appears also when one fixes the maximum momentum and uses the Mach number of the shock as the order parameter. The results are shown in Fig. 3, where Mach numbers between 10 and 500 have been considered, at fixed $`p_{max}=10^5mc`$. One can see that for Mach numbers below $`100`$ there is only one solution. For Mach numbers between $`100`$ and $`500`$ three solutions appear, one of which roughly corresponds to an unmodified shock. For larger Mach numbers, only this linear solution remains, and the shock is always only weakly modified.

The same situation is plotted in Fig. 4, where instead of $`R_{tot}`$ on the x-axis there is $`4R_{sub}`$, and $`R_{sub}`$ is the compression coefficient at the gas subshock. For unmodified shocks one expects $`4R_{sub}0`$.

Again three regions can be identified, in the parameter $`M_0`$.

1) For $`M_0M_{cr}^{(1)}`$ only one solution exists and does not necessarily correspond to unmodified shocks. In fact one can see that for $`M_0=10<M_{cr}^{(1)}100`$, the compression coefficient at the subshock is $`2`$, and $`R_{tot}8`$, for the situation plotted in Figs. 3 and 4. Even smaller values of $`M_0`$ do not give a perfectly unmodified shock.

2) For $`M_{cr}^{(1)}M_0M_{cr}^{(2)}`$ three solutions appear, two of which imply a strong modification of the shock. For the parameters used in Figs. 3 and 4, $`M_{cr}^{(2)}500`$.

3) For $`M_0M_{cr}^{(2)}`$ only the nearly unmodified solution exists.

In order to emphasize the fact that the multiple solutions actually correspond to physical solutions with very different spectral characteristics for the accelerated particles we plot in Fig. 5 the spectra for $`M_0=50`$ (for which there is only one solution) and $`M_0=150`$ (for which there are three solutions). In both cases we used $`P_{max}=10^5mc`$.

The solution corresponding to Mach number $`M_0=50`$ has $`R_{sub}=2.683`$ and $`R_{tot}=25.588`$, resulting in a total pressure of the accelerated particles $`P_{CR}/\rho _0u_0^2=0.8946`$ (solid line in Fig. 5). The three dashed lines in Fig. 5 are for $`M_0=150`$ and represent the spectra for the three solutions. In numbers, these solutions (from top to bottom in Fig. 5) are summarized as follows:

first solution: $`R_{sub}=3.0775`$, $`R_{tot}=55.607`$, $`P_{CR}/\rho _0u_0^2=0.9611`$

second solution: $`R_{sub}=3.9839`$, $`R_{tot}=14.318`$, $`P_{CR}/\rho _0u_0^2=0.734`$

third solution: $`R_{sub}=3.999371`$, $`R_{tot}=4.2548`$, $`P_{CR}/\rho _0u_0^2=0.06`$

In the cases of strong shock modification, the asymptotic spectra for $`(p/mc)1`$ are $`p^2f(p)p^{3/2}`$, while the linear theory, in case of strong shocks ($`M_01`$), would predict $`p^2f(p)p^2`$. In the region of very low energies, the spectra of the reaccelerated particles tend to zero, as known from the linear theory as well.

### 4.2 Modified shocks: the general case

Nonlinear effects in shock acceleration were first investigated in the case when a fraction of the thermal gas crossing the shock surface is energized to nonthermal energies. In this section we wish to apply the analytical method discussed in §2 and §3 to the most general case in which thermal particles are accelerated but seed particles may already be present in the environment.

We start our discussion with the case of acceleration of particles from the thermal distribution, in order to show that the multiple solutions already found in other analytical approaches  are also obtained by adopting the approach illustrated here. In Fig. 6, we plot the $`U(p_{max})`$ as a function of $`R_{tot}`$ for the case in which a fraction $`\eta `$ of the particles (as indicated in the plots) crossing the shock is actually accelerated to suprathermal energies. The Mach number is chosen to be $`M_0=150`$ and the maximum momentum is taken as $`P_{max}=10^5mc`$. One can easily see that $`R_{tot}4`$, the value for unmodified shocks, for small values of $`\eta `$ (low efficiencies), while $`R_{tot}`$ increases above 4 for $`\eta >10^4`$. Three solutions appear for intermediate efficiencies, while the solution for the accelerated particles always predicts a strongly modified shock for high efficiencies, $`\eta >3\times 10^3`$. An asymptotic value of $`R_{tot}60`$ is achieved for the parameters used here.

This example demonstrates that the critical behaviour shown to appear for high values of $`\eta `$ in several analytical or semi-analytical calculations is in fact also predicted by the approach proposed here. Next question to ask is however whether the presence of seed particles can change the critical behavior of shocks. In order to answer this question we consider a situation in which particles are accelerated from the thermal pool with efficiency $`\eta `$, and at the same time the shock propagates in a medium where the preshock pressure in seed particles equals the gas pressure. This situation is considered to resemble that of a supernova explosion in our Galaxy, where cosmic rays fill the volume remaining in quasi-equipartition with the gas. In order to test the critical structure of the shock we calculate $`U(p_{max})`$ as a function of the compression factor $`R_{tot}`$ between upstream infinity and downstream. Our results are plotted in Fig. 7. The lines refer to the cases $`\eta =10^5,10^4,5\times 10^4,10^3,5\times 10^3`$ from bottom to top, as in Fig. 6.

It is clear from this plot that the shock may be modified by the backreaction of the accelerated particles even for very low values of $`\eta `$, because of the presence of seed particles. In other words, the nonlinearity of the shock can well be dominated by the presence of preaccelerated particles rather than by the acceleration of particles from the thermal pool. Clearly this depends however on the values of the parameters (injection momentum, $`\eta `$, maximum momentum achievable, Mach number). The injection momentum, in principle, could be related to $`\eta `$ in order to reduce the number of free parameters of the problem. This is possible in the assumption that the particles downstream keep a thermal distribution. Simulations suggest that the particles injected in the accelerator are the ones with momentum a few times the thermal momentum of the particles downstream, so that $`\eta `$ can be simply calculated by integration of the Maxwell-Boltzmann (MB) distribution above this minimum momentum. Unfortunately, we do not know whether the particle distribution downstream is in fact Maxwell-Boltzmann-like. Moreover, we do not know exactly what is the threshold to impose on the injection momentum, which is a delicate issue because the number of particles taking part to the acceleration process would result from the integral of the MB distribution over its exponentially decreasing part. Therefore we preferred here to keep $`\eta `$ and the injection momentum as separate free parameters.

## 5 Discussion and Conclusions

We proposed a semi-analytical approach to show that the backreaction of particles accelerated at a shock is able to affect the shock itself, in such a way that the shock and the accelerated particles become parts of a nonlinear system. In particular, for the first time we included in this kind of calculations the seed particles that may be present in the region where the shock is propagating and that can be re-energized by the shock.

While a test-particle approach to this problem was first presented in , a nonlinear treatment was never investigated. In the pioneering work of Ref. , it was recognized that the energy of the seed particles could be enhanced by about one order of magnitude at each shock passage, and that after an infinite number of strong shocks passing through the region, the spectrum of the particles would tend to the asymptotic spectrum $`E^{3/2}`$. Two comments are in order. First, the continuous increase of the cosmic ray energy due to the re-energization of seed particles leads unavoidably the shock to be modified by the nonthermal pressure, unless one starts from an uninterestingly small pressure of seed particles at the beginning. Second, the fact that the spectrum becomes flatter than $`E^2`$, which is the result for a strong non-relativistic shock, implies that most of the nonthermal energy is pushed to the highest energies, therefore the shock again can be more easily modified. Both these points suggest that a nonlinear treatment of shock re-acceleration is required.

In our Galaxy, cosmic rays are observed to be in rough equipartition with the gas pressure and with magnetic fields, therefore supernova shocks or shocks generated in other environments propagate in a medium in which the seed particles (cosmic rays) are non-negligible. In these circumstances the non-linear effects may be very important. We showed here that in fact for some regions of the parameter space, the shock is modified mainly by the backreaction of the seed particles rather than by the cosmic rays accelerated at the shock from the thermal pool. The spectra of the re-accelerated particles have also been calculated.

The interesting phenomenon of the appearance of multiple shock solutions, already known for the case of shock acceleration, appears also for the case of reacceleration. This puzzling phenomenon may suggest that the shock behaves as a self-regulating system settling on the critical point, as proposed in . On the other hand, it is possible that the multiple solutions may be the artifact of some of the assumptions used in the analytical approach, in particular the request for time independent (stationary) solutions and the fact that the role of the self-generated waves on the diffusion coefficient is not taken into account. Further investigation, in particular in the direction of a detailed comparison of our results with numerical simulations of shock acceleration is required in order to unveil the physical meaning of the multiple solutions for modified shocks.

From the phenomenological point of view, it would be certainly worth to study the implications of nonlinear shock reacceleration on the nonthermal activity in astrophysical environments where the effect is expected to play an important role, in particular in the case of supernova remnants. In particular, as pointed out by the referee, the suggested dominance of reaccelerated ambient seed particles over freshly injected particles may have serious implications for the spectra of secondary nuclei resulting from spallation processes.

## Acknowledgments

The author gratefully acknowledges useful discussions with M. Baring, D. Ellison and M. Vietri. The author is also grateful to the referee, L.O’C. Drury for his very useful remarks. This work was partially supported through grant COFIN 2002 at the Arcetri Astrophysical Observatory.


format: text

text:
# Variability Tests for Intrinsic Absorption Lines in Quasar Spectra

## 1 Introduction

The absorption lines in quasar spectra can place important constraints on the basic properties of quasar environments, such as the outflow velocities, column densities, mass loss rates, and elemental abundances. Additionally, if quasars reside in the nuclei of massive galaxies, then we can use the chemical abundances to probe indirectly the extent and epoch of star formation in young galaxies (Schneider 1998, Hamann & Ferland 1999).

Quasar absorption lines can be divided into three categories based on the line widths: broad absorption lines (BALs), narrow absorption lines (NALs), and intermediate mini-broad absorption lines (mini-BALs). The divisions between these categories are arbitrary. Classic BALs typically have full widths at half minimum (FWHMs) of order 10,000 km s-1, but lines as narrow as 2000 to 3000 km s-1 are often still considered BALs (Weymann et al. 1991). A useful working definition of the NALs is that they have FWHMs smaller than the velocity separation of major absorption doublets (e.g., 500 km s-1 for C IV $`\lambda \lambda `$1548,1551, or 960 km s-1 for N V $`\lambda \lambda `$1239,1243), but often these lines are much narrower. Mini-BALs have widths intermediate between the NALs and BALs.

BALs appear in about 10% to 15% of optically selected quasars (Weymann et al. 1991). They appear at blueshifted velocities (relative to the emission redshift) ranging from near 0 km s-1 to $`>`$30,000 km s-1, and they clearly form in high velocity outflows from the quasar engines. Mini-BALs appear to be less common, although to our knowledge no one has yet done a quantitative inventory. Nonetheless, mini-BALs appear at the same range of blueshifted velocities as the BALs and it is thought that they too form in quasar winds (Hamann et al. 1997a, Jannuzi et al. 2003). Even though they may have complex profiles, high resolution spectra show that both BAL and mini-BAL profiles are “smooth” compared to thermal line widths, and therefore these lines are not simply blends of many NALs (Barlow & Junkkarinen 1994, Hamann et al. 1997a, Junkkarinen et al. 2001, Hamann et al. 2003).

NALs appear at a wide range of velocity shifts. They are further classified as associated absorption lines (AALs) if the absorption redshift, $`z_{\mathrm{abs}}`$, is within $`\pm `$5000 km s-1 of the emission line redshift ($`z_{\mathrm{abs}}z_{\mathrm{em}}`$; Weymann et al. 1979, Foltz et al. 1986, Anderson et al. 1987, Foltz et al. 1988). A significant fraction of AALs are believed to have a physical relationship with the quasars, based on statistical correlations between the occurrence of AALs and quasar properties (see also Aldcroft et al. 1994, Richards et al. 1999). However, in general, AALs and other NALs can form in a variety of locations, such as cosmologically intervening clouds, galaxies that are unrelated to the quasar, and clouds that are physically associated with, or perhaps ejected from, the quasar (Weymann et al. 1979). NALs must therefore be examined individually to determine if they are intrinsic to the quasar. Several diagnostics have been proposed for this purpose, including i) line strength variations over time, ii) profiles that are smooth and broad compared to thermal line widths, iii) partial line of sight coverage of the background emission source, and iv) high densities based on excited-state absorption lines (e.g., Barlow & Sargent 1997, Barlow, Hamann & Sargent 1997, Hamann et al. 1997a and 1997b).

Line strength variations can be caused by bulk motions across the line of sight or by changes in the ionization state of the gas. A change in the radial velocity of the absorbing material could also produce a shift in the wavelength (redshift) of the measured lines (although this type of shift appears to be extremely rare, Gabel et al. 2003, this paper). In either case, variations over short time scales are incompatible with absorption in large intergalactic clouds. A firm lower limit on the size of intergalactic HI-absorbing clouds is approximately a few kpc (Foltz et al. 1984, McGill 1990, Bechtold & Yee 1995, Rauch 1998), which is similar to recent direct estimates of the sizes of intergalactic C IV-absorbing clouds (Tzanavaris et al. 2003). Assuming the clouds are uniform and do not have sharp edges, the time scale for absorption line strength variations occurring via bulk motions will be roughly the time needed for this minimum cloud to cross our line of sight. If the clouds have maximum transverse velocities of $``$1000 km s-1 (comparable to the line-of-sight velocity dispersion in massive galaxy clusters), then in the $``$30 years spanned by modern observations the cloud would travel just $``$$`3\times 10^5`$ kpc, roughly five orders of magnitude less than its diameter. Variability due to changes in the ionization state is also problematic for intergalactic absorbing clouds. Changes in the ionization state are nominally limited by the time scale for recombination, e.g., if the cloud is in ionization equilibrium. Recombination on time scales of $``$30 years requires densities of $`n_e100`$ cm-3 (see §5.1 below), which is exceedingly high compared to the densities $`n_e10^3`$ cm-3 expected in intergalactic H I clouds (Miralda-Escude et al. 1996, Rauch 1998). The variation time scales might be shorter for lines forming in the (relatively) dense interstellar medium of intervening galaxies, but these absorption systems should be recognizable, e.g., via damping wings in Ly$`\alpha `$. Therefore, line strength variations in NAL systems that do not include damped Ly$`\alpha `$ strongly suggest that the absorption is intrinsic to the quasar.

We are involved in a multi-faceted program to identify intrinsic NALs and use them to place constraints on the physical properties of quasar environments. In this paper, we examine multi-epoch rest frame UV spectra of eight redshift $``$2 quasars known to have AALs and mini-BALs. We find variability in 2 out of 8 AALs and 2 out of 2 mini-BALs. We use these results to constrain the densities and locations of the absorbing gas, and briefly discuss the implications for quasar wind models.

## 2 Observations and Data Reductions

We obtained spectra of eight quasars at the Shane 3m telescope at the University of California Observatories (UCO) Lick Observatory as well as at the Smithsonian Institution and University of Arizona Multiple Mirror Telescope (MMT) Observatory. Additional spectra obtained at the Palomar Observatory were generously provided to us in digital form by C. C. Steidel, first published in Sargent et al. (1988, hereafter SBS88). Table 1 shows a log of the different observing runs, including the dates, approximate wavelength ranges, and resolutions. The data span $``$17 years ($``$6 years in the quasar rest frames) with two to four observations per quasar.

In Table 2 we list the objects observed, their emission redshifts ($`z_{\mathrm{em}}`$, from Hewitt & Burbidge 1993), the absorption redshifts ($`z_{\mathrm{abs}}`$) and velocity shifts (relative to $`z_{\mathrm{em}}`$) of their AALs, mini-BALs, or, in the case of Q0151$`+`$048 only, a high-velocity NAL system that appeared to vary. All of the AAL and NAL redshifts are from SBS88, except for PG0935+417 (from Hamann et al. 1997a) and Q0848+163 (see §4 below). Note that the AAL redshift given for PG0935+417 is an approximate average of a complex of lines that is not resolved in our Lick spectra. Please see SBS88 for a more complete list of absorption lines in the other quasars. Table 2 also lists the observation numbers (from Table 1), and a note on the variability of each system. The variabilities are discussed further in §4 below.

We chose the quasars for this study because they are bright, they are known to have AALs, and their redshifts allow ground-based observations of strong UV absorption lines, such as Ly$`\alpha `$ $`\lambda `$1216, N V $`\lambda \lambda `$1239,1243, Si IV $`\lambda \lambda `$1394,1403, C IV $`\lambda \lambda `$1548,1551, and potentially others in this wavelength range (see SBS88 and Junkkarinen et al. 1991).

We reduced the data using standard techniques with the IRAF111The Image Reduction and Analysis Facility (IRAF) is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc. under cooperative agreement with the National Science Foundation. software. We bias subtracted the data and then divided the raw two dimensional images by flat fields. Wavelength calibration was achieved using internal lamp spectra. We flux calibrated the spectra using standard stars measured the same night. Because the weather was sporadically cloudy, we used the flux calibrations primarily to recover an accurate spectral shape, rather than absolute fluxes. Finally, we added together spectra obtained at the same wavelengths on the same dates to increase the signal to noise ratio. The SBS88 spectra have a higher resolution than our Lick spectra; we therefore smoothed the SBS88 data using a boxcar routine to match the resolution of the Lick spectra and facilitate comparisons.

## 3 Results

We compared the spectra from the different observations of each quasar to look for variations in the absorption line strengths. Every absorption line within our wavelength coverage was examined, except those in the Ly$`\alpha `$ forest. Table 2 presents the main results. Two out of the eight quasars showed AAL variability. Two out of two mini-BALs varied, and one quasar showed probable variations in a high-velocity C IV NAL. Notes on the variable systems are provided in §4 below. Spectra showing the variations are also plotted in Figures 1–4. The wavelength ranges in these plots were chosen to show only the specific NALs or mini-BALs that varied. The peak flux in the plotted wavelength range of the reference spectrum (solid curve) is normalized to unity in all cases. Other spectra are over-plotted and scaled to match approximately the reference spectrum in the continuum near the lines of interest.

Our ability to detect line variability depends on the line strength and the signal to noise ratio of the spectra. For the six AALs that did not vary, we estimate upper limits on the changes, by inspection, of the strongest measured lines (generally C IV and N V). In particular, we estimate that the AAL equivalent widths varied by $``$20 % in Q0207$``$003, Q0348$`+`$061, PG0935$`+`$417, and Q0958$`+`$551, and $``$15 % in Q1159$`+`$123.

## 4 Notes on Variable Sources

### 4.1 Q0151$`+`$048

Figure 1 shows four epochs of observations for Q0151+048. We measured clear variations in a C IV mini-BAL at $`z_{\mathrm{abs}}=1.6581`$ and tentative changes in a high-velocity C IV NAL at $`z_{\mathrm{abs}}=1.6189`$. No other lines are detected at these redshifts. SBS88 report that both components of the narrow C IV doublet at $`z_{\mathrm{abs}}=1.6189`$ are blended with Si II $`\lambda `$1526 absorption (see Fig. 1). However, much higher resolution Keck spectra (Hamann et al. 2003) obtained in December 1996 provide no evidence for these Si II lines. In particular, we do not detect other low ionization lines such as FeIV $`\lambda `$1608 or CII $`\lambda `$1334. Also, the redshifts and profiles of the individual C IV doublet lines at $`z_{\mathrm{abs}}=1.6189`$ match each other very well, indicating that all of the absorption is due to C IV (at least for the epoch December 1996). The data shown in Figure 1 suggest that these lines varied. In particular, the equivalent width of the short-wavelength component of the C IV doublet is roughly twice as strong in 1981 and 1996 compared to 1997. However, the two doublet components seem to have varied by different amounts, which is unphysical. Therefore, the variability in this system is questionable. If this system is intrinsic, its displacement from the emission lines implies an outflow velocity of $``$32,900 km s-1.

SBS88 describes the C IV mini-BAL in Q0151$`+`$048 as a blend of four narrow C IV absorption systems with redshifts ranging from $`z_{\mathrm{abs}}=1.6533`$ to $`z_{\mathrm{abs}}=1.6600`$. We consider this blend rather to be a mini-BAL, based on its appearance in the SBS88 spectrum and on the absence of discrete narrow components in our high resolution Keck spectrum (Hamann et al. 2003). We measure the observed frame equivalent widths in this feature over the epochs monitored from 0.38 Å to 3.8 Å, signifying a factor of $``$10 variation. The FWHM of the mini-BAL (as observed in Oct. 1981) is $``$1020 km s-1, and its displacement from the emission redshift is $``$28,400 km s-1.

The AALs in this quasar did not appear to vary.

### 4.2 PKS0424$``$131

Figure 2 shows the variable N V and C IV AALs at $`z_{\mathrm{abs}}=2.1330`$. The change in the N V doublet is obvious. The short wavelength component of the C IV doublet is blended with an unrelated line of FeII at $`z_{\mathrm{abs}}=1.0352`$. Ly$`\alpha `$ in this AAL system is not clearly detected because of blending with another unrelated line. We also do not detect Si IV absorption at the AAL redshift. The equivalent widths of the N V and C IV AALs varied by 55% and 17%, respectively, between the two observations.

Petitjean et al. (1994) published high resolution (0.3 Å) spectra of PKS0424$``$131 obtained in 1992. They report a total observed frame equivalent width in the N V doublet of 1.62 Å, compared to 3.57 Å and 1.72 Å in our data from 1981 and 1997, respectively. Petitjean et al. (1994) also measured an equivalent width of 0.31 Å in the unblended long-wavelength component of the C IV doublet, compared to our equivalent width measurements of 3.9 Å and 0.36 Å in our 1981 and 1997 spectra. These comparisons suggest that the largest variations in these AALs occurred between the 1981 SBS88 and 1992 Petitjean et al. observations.

### 4.3 Q0848$`+`$163

Higher resolution spectra (Hamann et al. 2003) confirm that there are several AAL systems in this quasar. The two systems that SBS88 identified at $`z_{\mathrm{abs}}1.916`$ and $`z_{\mathrm{abs}}1.917`$ are blended together in our Lick spectra. We will therefore refer to them hereafter as one system at $`z_{\mathrm{abs}}`$ = 1.9165. There is another AAL system at $`z_{\mathrm{abs}}`$ = 1.9105 that SBS88 did not identify, but that is unmistakable in our high resolution Keck spectrum (Hamann et al. 2003).

Figure 3 provides some evidence for variability in these AAL systems222The spectral region covering the CIV AALs in 0848+163 was measured and reported by SBS88. However, these data were not available for the present study, and are therefore excluded from Figure 3.. Comparisons involving the February 1983 spectrum are ambiguous because the underlying emission lines varied substantially between that epoch and the more recent Lick observations. However, the two 1997 observations suggest that similar variations occurred in the N V and C IV doublets. The evidence for variability is perhaps stronger for the system at $`z_{\mathrm{abs}}`$ = 1.9105, where the equivalent width in the unblended short-wavelength component of the N V doublet declined by approximately 50%. The AALs at $`z_{\mathrm{abs}}=1.9165`$ might also have varied in 1997, but the changes coinciding with the short-wavelength doublet components in this system might be due at least in part to the variability at $`z_{\mathrm{abs}}=1.9105`$. The variability in the accompanying Ly$`\alpha `$ lines is masked by blending with Ly$`\alpha `$ forest lines. We do not detect the Si IV (or other) lines in these systems.

### 4.4 PG0935$`+`$417

Figure 4 shows distinct variation in a C IV mini-BAL that is blueshifted by $``$52,000 km s-1 with respect to the emission lines (see also Hamann et al. 1997a). The identification of this feature with highly blueshifted C IV is confirmed by the appearance of N V and O VI absorption at the same velocity (Rodriguez Hidalgo & Hamann 2003). This is the second highest observed velocity for a mini-BAL, at nearly $``$0.2$`c`$, after PG2302$`+`$029 (Jannuzi et al. 2003). High-resolution Keck observations of PG0935$`+`$417 (Hamann et al. 1997a) show that the mini-BAL remains smooth down to a resolution of 7 km s-1. This mini-BAL is therefore truly a broad line with a continuous range of absorption velocities, and not a blend of many NALs. The most dramatic variation occurred between January 1993 and March 1996 (Figure 4), where the observed frame equivalent width in the mini-BAL increased from $``$5.12 Å in 1993 to $``$6.14 Å in 1996 ($``$19%). However, the change is most obvious in the shape and centroid of the mini-BAL profile. The centroid of the line moved from $``$3917 Å to $``$3865 Å (corresponding to a shift of $``$4000 km s-1) between 1993 and 1996. The FWHM of the mini-BAL was $``$1885 km s-1 in 1993 and $``$1320 km s-1 on the other dates.

The complex of AALs in this quasar did not vary (see also Hamann et al. 2003).

## 5 Analysis

The line variability time scales place constraints on the physical properties of the absorbing clouds. Here we assume the variations were caused by changes in the ionization state, which leads to estimates of the minimum electron density and maximum distance from the continuum source (for a photoionized plasma). The results are listed in Table 3, where $`t_{\mathrm{var}}`$ is the smallest observed variability time in the quasar rest frame, $`n_\mathrm{e}`$ is the electron density, and $`R`$ is the distance from the continuum source. Note that the values of $`n_\mathrm{e}`$ and $`R`$ are limits because we measure only upper bounds on the variability times.

### 5.1 Minimum Electron Densities

We assume for simplicity that the gas is in ionization equilibrium and the ions we measure are the dominant ionization stages. In that case, the variability time is limited by the recombination time given by,

$$n_e\frac{1}{\alpha _{i1}t_{i1}^{\mathrm{recom}}}$$
(1)
where $`\alpha _{i1}`$ is the rate coefficient for recombination from the observed ion stage $`i`$ to the next lower stage $`i1`$, and $`t_{i1}^{\mathrm{recom}}`$ is the corresponding recombination time (see Hamann et al. 1997b). Arnaud and Rothenflugh (1986) give recombination rates for C IV $``$ CIII and N V $``$ NIV as 2.8 $`\times `$ $`10^{12}`$ cm3 s-1 and 5.5 $`\times `$ $`10^{12}`$ cm3 s-1, respectively. These values, together with a nominal temperature of 20,000 K (Hamann et al. 1995) and the maximum recombination times set by our observations, give the minimum electron densities presented in Table 3. Because of the particular values of the recombination rates, the calculated minimum electron densities for N V are consistently a factor of $``$2 less than those for C IV. Thus, we list constraints on the minimum electron densities in Table 3 based only on the C IV recombination rates.

### 5.2 Maximum Distances from Continuum Source

We derive the maximum distance between the absorbing clouds and the continuum source by assuming the gas is in photoionization equilibrium with an ionization parameter given by,

$$U=\frac{1}{4\pi R^2n_H\mathrm{c}}_0^{\lambda _{\mathrm{LL}}}\frac{\lambda L_\lambda }{hc}d\lambda $$
(2)
where $`n_Hn_e`$ is the hydrogen density, $`\lambda _{\mathrm{LL}}=912`$ Å is the wavelength at the Lyman limit, and $`L_\lambda `$ is the quasar luminosity distribution. We assume all of the quasars in our sample have the same continuum shape characterized by a segmented power law, $`L_\lambda \lambda ^\alpha `$, where $`\alpha =1.6`$ from 1000 Å to 100,000 Å , $``$0.4 from 10 Å to 1,000 Å , and $``$1.1 from 0.1 Å to 10 Å (Zheng et al. 1997, Telfer et al. 2002). We also assume that the absorbing clouds are not “shielded” from the continuum emission source. Thus the continuum spectrum incident on the absorbing clouds is diminished only by the geometric 1/$`R^2`$ dilution. Integrating over our adopted continuum shape in Equation 2 yields a convenient expression:

$$U0.09\frac{L_{47}}{n_{10}(R_1)^2}$$
(3)
where L47 is the bolometric luminosity in units of 1047 erg s-1, $`n_{10}`$ is the density in units of 1010 cm-3, and $`R_1`$ is the distance from the continuum source in parsecs. Note that for this continuum shape, the quasar bolometric luminosity is $`L_{\mathrm{bol}}4.4\lambda L_\lambda `$ at $`\lambda =1450`$ Å.

We derive bolometric luminosities for the quasars in our sample based on the continuum shape given above, observed $`B`$ magnitudes from Junkkarinen et al. (1991) and Hewitt & Burbidge (1993), and a cosmology with $`H_0=72`$ km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. For PKS0424$``$131, we use a measured $`V`$ magnitude and a representative value of $`BV=0.12`$ (based on other luminous quasars at this redshift, Hewitt & Burbidge 1993). The resulting values of $`L_{47}`$ are listed in Table 3. Finally, we adopt an ionization parameter of U $``$ 0.02, which is approximately optimal for the C IV and N V ions (e.g., Hamann et al. 1995, Hamann 1997a). Plugging into Equation 3 with the previously derived density limits yields the upper limits on $`R`$ listed in Table 3.

## 6 Discussion

Two of the 8 quasars in our study showed variability in their AALs. This fraction is similar to the 3 out of 15 AAL systems showing variability in a sample of lower redshift quasars studied by Wise et al. (2003). Note that these results provide only minimum variability fractions because of the limited time sampling. The fraction of AALs that are intrinsic based on variability is therefore $``$25%. We are now analyzing much higher resolution spectra of all of the quasars in our sample (Hamann et al. 2003) to look for other evidence of intrinsic absorption (§1).

The mini-BALs appear to vary more often than the AALs. Our study showed variations in 2 out of 2 systems. The only other mini-BALs tested for variability, that we know of, are PG 2302+029 (Jannuzi et al. 2003) and Q2343+125 (which is marginally a mini-BAL with FWHM $``$ 400 km s-1, Hamann et al. 1997c). Both of those systems also varied between just two observations.

The marginal detection of variability in the high-velocity NAL system of Q0151+048 is surprising but not unprecedented. Richards et al. (1999) argued that a significant fraction of NALs at blueshifted velocities from 5000 to 75,000 km s-1 are intrinsic to quasars, based on statistical correlations between the appearance of these lines and the quasar radio properties. The high-velocity NAL system in Q0151+048 might be one specific example at a velocity shift of $``$32,900 km s-1.

It is interesting to note that, with one exception, the variable absorption lines remained fixed in velocity. The exception is the mini-BAL in PG0935+417 whose changing centroid and FWHM was noted previously by Hamann et al (1997a). However, it is not clear that these velocity changes represent real changes in the absorber kinematics. It could be that the apparent velocity changes are caused by variations in the line strength (optical depth) across a fixed range of absorption velocities. A more likely example of a real velocity shift involves an AAL system in the Seyfert I galaxy, NGC 3783 (Gabel et al. 2003).

All of the variable absorption lines in our sample appear at blueshifted velocities, clearly indicating outflows from the quasars. However, the wind geometry implied by these AALs and mini-BALs must be different from the classic BALs. In particular, the line widths are small compared to the wind speed. BALs typically have $`V`$/FWHM of order unity (where $`V`$ is the outflow speed). In contrast, the mini-BAL in Q0151+048 has $`V`$/FWHM $``$ 28, and the high-velocity NAL in that quasar has $`V`$/FWHM $``$ 55. Another intrinsic absorber in Q2343+125 (Hamann et al. 1997c) has $`V`$/FWHM $``$ 60. Absorption lines like these with $`V`$/FWHM $``$ 1 probably do not form in continuous outflows that accelerate from rest along our line of sight to the continuum source. If that were the case, we should see absorption at all velocities from $`V`$ to 0, unless, perhaps, there is a peculiar ionization structure that somehow favors C IV and similar ions only at the very narrow range of observed line velocities. A more likely possibility is that the observed lines form in discrete “blobs” of gas that were ejected/accelerated sometime prior to the observations. However, the nature of these blobs as coherent entities is unclear. The large velocity dispersions implied by the line widths (FWHM $`>`$ 1000 km s-1 for the mini-BALs) should quickly lead to a spatially extended structure. Another possibility is that the flows intersect our line of sight to the continuum source (and become observable in absorption) only after they reach the observed high speeds (see also Hamann et al. 1997c, Elvis 2000). Thus the acceleration occurs somewhere outside of our sight line to the continuum source. In this case the flow could be continuous – simultaneously spanning a wide range of velocities and radial distances, $`R`$, but the absorption lines sample only a limited range of these parameters, depending on the specific geometry and orientation.

One plausible geometry involves a wind emanating from the quasar accretion disk (Murray et al. 1995, Proga, Stone & Kallman 2000). If the lines we measure form anywhere near the point of origin of these winds, within several parsecs of the central black hole and possibly much closer, then simple considerations based on Equation 3 suggest that the actual gas densities are $`n_e10^9`$ cm-3. Unfortunately, the lower limits on $`n_\mathrm{e}`$ in Table 3 are not very constraining. They are inconsistent with only the most extreme lower densities and large radial distances ($``$10 kpc) derived for some AAL systems (based on excited-state absorption lines, e.g., Morris et al. 1986, Tripp, Lu & Savage 1996, Hamann et al. 2001).

## 7 Conclusions

We observed eight AAL quasars to test for variability in their absorption line strengths as an indicator of intrinsic absorption. Two of these quasars were also known previously to have high-velocity mini-BALs. In our limited time sampling of 2–4 observations per quasar, we found variability in both of the mini-BALs, 2 out of 8 AAL systems, and possibly one additional high-velocity NAL system. These results agree with previous reports of frequent variability in mini-BALs (Hamann et al. 1997a, Jannuzi et al. 2003), and with the recent finding of AAL variability in 3 out of 15 quasars studied by Wise et al. (2003). The short time scales over which the lines varied in our study (sometimes $`<`$0.28 years in the quasar rest frame) implies that the absorbers are dense, compact, and physically associated with the quasars. We estimate minimum electron densities (from the recombination time) ranging from $``$1900 cm-3 to $``$40,000 cm-3. We conclude that the fraction of AALs that are intrinsic based on variability is at least $``$25%. The fraction of mini-BALs that are intrinsic may be $``$100%. The outflow velocities implied by the intrinsic systems range from $``$1500 km s-1 to $``$52,000 km s-1.

We are grateful to C.C. Steidel for generously providing the Palomar data. This work was supported by NSF grant AST 99-84040. D.N. thanks Craig Warner for helpful discussions, and acknowledges the University Scholars Research Program at the University of Florida for financial support.


text:
# Towards Noised-Induced Phase Transitions in Systems of Elements with Motivated Behavior

## 1 Introduction

Systems of elements with motivated behavior (systems with motivation), e.g., fish and bird swarms, car ensembles on highways, stock markets, etc. often display noise-induced phase transitions (for a review see Ref. ). The ability of noise to induce phase transitions is now well established (see, e.g., Refs ). However, in systems with motivation there is a special mechanism endowing the corresponding noise-induced phase transitions with distinctive properties.

For example, people as elements of a certain system cannot individually control all the governing parameters. Therefore one chooses a few crucial parameters and focuses on them the main attention. When the equilibrium with respect to these crucial parameters is attained the human activity slows down retarding, in turn, the system dynamics as a whole. For example, in driving a car the control over the relative velocity $`v`$ is of prime importance in comparison with the correction of the headway distance $`x`$. So, under normal conditions a driver should eliminate the relative velocity between her car and a car ahead first and only then correct the headway.

These speculations lead us to the concept of dynamical traps, a certain “low” dimensional region in the phase space where the main kinetic coefficients specifying the characteristic time scales of the system dynamics become sufficiently small in comparison with their values outside the trap region . The present paper analyzes the effect of noise on such a system and demonstrates that additive noise in a system with dynamical traps is able to give rise to new phases.

## 2 Noised oscillatory system with dynamical traps

By way of example, the following dimensionless system typically used to describe the oscillatory dynamics is considered:

$$\frac{dx}{dt}=v,\frac{dv}{dt}=\mathrm{\Omega }(v)\left[x+\sigma v\right]+ϵ\xi (t),$$
(1)
where $`\sigma `$ is the damping decrement and the term $`ϵ\xi (t)`$ is a random Langevin “force” of intensity $`ϵ`$ proportional to the white noise $`\xi (t)`$ with unit amplitude. The function $`\mathrm{\Omega }(v)`$ describes the dynamical trap effect arising in the vicinity of $`v=0`$. For this function, the following simple Ansatz

$$\mathrm{\Omega }(v)=\frac{v^2+\mathrm{}^2}{v^2+1}.$$
(2)
is used. In the chosen scales the thickness of the trap region is equal to unity and the parameter $`\mathrm{}1`$ measures the trapping efficacy. When $`\mathrm{}=1`$ the dynamical trap effect is ignorable, for $`\mathrm{}=0`$ it is most effective.

System (1), (2) was analyzed numerically using the algorithm described in . Figure 1 shows the distribution function $`𝒫(x,v)`$ of the system on the phase plane $`\{x,v\}`$, depending on the parameter $`\mathrm{}`$. As can be seen this system undergoes a second order phase transition manifesting itself in the change of the shape of the phase space density $`𝒫(x,v)`$ from unimodal to bimodal as the trap parameter $`\mathrm{}`$ decreases.

### Mechanism of the phase transition

Figure 2 illustrates the mechanism of this phase transition. It depicts a typical fragment of the system motion through the trap region for $`\mathrm{}1`$. When the path goes into the trap region $`𝒬_t`$ ($`|v|1`$) the regular “force” $`\mathrm{\Omega }(v)(x+\sigma v)`$ is depressed. So inside this region the system dynamics is mainly random due to the remaining weak Langevin “force” $`ϵ\xi (t)`$. Crucial is the fact that the boundaries $`_+𝒬_t`$ (where $`v1`$) and $`_{}𝒬_t`$ (where $`v1`$) are not identical in properties with respect to the system motion. At the boundary $`_+𝒬_t`$ the regular “force” leads the system inwards the trap region $`𝒬_t`$, whereas at the boundary $`_{}𝒬_t`$ it causes the system to leave the region $`𝒬_t`$. Outside the trap region $`𝒬_t`$ the regular “force” is dominant. Thereby, from the standpoint of the system motion inside the region $`𝒬_t`$, the boundary $`_+𝒬_t`$ is “reflecting” whereas the boundary $`_{}𝒬_t`$ is “absorbing”.

As a result the distribution of the residence time at different points of the region $`𝒬_t`$ should be asymmetric as shown in Fig. 1 (it is most clear in the lower right window). Therefore, during location inside the trap region the mean velocity must be positive, causing the system to go away from the origin. Outside the trap region the system motion is damping. So, when the former effect becomes sufficiently strong the distribution function $`𝒫(\eta ,u)`$ becomes bimodal.

## 3 Chain of oscillators

A similar noise-induced phase transition for an ensemble of oscillators with dynamical traps is analyzed. The following one-dimensional model is considered, $`N`$ balls can move along $`x`$-axis interacting with the nearest neighbors, which is described by the equations for $`i=1,2,\mathrm{}N1`$

$$\frac{dx_i}{dt}=v_i,\frac{dv_i}{dt}=\mathrm{\Omega }_s(v_{i1},v_i,v_{i+1})\left[\eta _i+\sigma \vartheta _i\right]+ϵ\xi _i(t).$$
(3)
Here $`x_i`$ is the coordinate of ball $`i`$, $`v_i`$ is its velocity, the variables $`\eta _i`$ and $`\vartheta _i`$ are given by the expression:

$$\eta _i=x_i\frac{1}{2}(x_{i1}+x_{i+1}),\vartheta _i=v_i\frac{1}{2}(v_{i1}+v_{i+1}),$$
(4)
and $`\{\xi _i(t)\}`$ is the collection of white noise sources of unit amplitude and being mutually independent. The damping decrement $`\sigma `$ and the noise intensity $`ϵ`$ are assumed to be the same for all the oscillators. The boundary balls ($`i=0`$ and $`i=N`$) are let to be fixed to prevent the ball system from moving as a whole. The function $`\mathrm{\Omega }_s(v_{i1},v_i,v_{i+1})`$ measuring the trap effect due to the nearest neighbor interaction is given by the Ansatz

$$\mathrm{\Omega }_s(v_{i1},v_i,v_{i+1})=\frac{(v_{i1}v_i)^2+(v_iv_{i+1})^2+\mathrm{}^2}{(v_{i1}v_i)^2+(v_iv_{i+1})^2+1},$$
(5)
where the parameter $`\mathrm{}`$ has the same meaning as previously, it measures the intensity of trapping. The balls are assumed to be either mutually permeable or impermeable. In the latter case the absolutely elastic collision approximation is used.

The system of equations (3)–(5) was analyzed numerically. Below, the results are presented for the impermeable balls. Again integration of the stochastic differential equations was performed with the algorithms described in . We analyzed an ensemble of 500 balls initially spaced 5 units apart and counted the number of balls falling into 100 identical intervals covering the system location region in the corresponding space.

Figure 3 actually depicts the ball distribution along the spatial axis and the velocity space at a fixed moment of time. As seen, the trap effect leads to the formation of a sufficiently inhomogeneous spatial distribution of balls whereas for the same system but without trapping, $`\mathrm{}=1`$, the balls are actually uniformly distributed over the space region. In the velocity space the trap effect causes the balls to spread over a much wider domain.

To make the trap effect in the velocity space more clear we averaged the velocity distribution over 500 time units. The result is presented in Fig. 4. As seen, the strong trap effect causes a nonanalytic behavior of the velocity distribution at zero value, it takes the form of a cusp.

Finally, Fig. 5 illustrates the time pattern of the system dynamics. It is clear that the shown spatial structure is characterized by a long life time, because, we recall, the period of individual ball oscillations without traps is about $`2\pi `$.

## 4 Conclusion

A new type of noise-induced phase transitions in systems of elements with motivated behavior is considered. Dynamics of such systems, we think, should exhibit a number of anomalies due to dynamical traps. The dynamical traps form a low dimensional region in the phase space where the kinetic coefficients become sufficiently small and, as a result, the system spends a long time in it. The cause of the dynamical traps is, e.g., the inability of people or animals to control all the system governing parameters. So they have to focus the main attention on a few crucial ones and the intensity of their activity decreases when the system attains a local quasi-equilibrium with respect to these parameters.

To illustrate this effect a simple oscillatory system $`\{x,v=\dot{x}\}`$ is studied when the trap region is located in the vicinity of the $`x`$-axis and without noise the stationary point $`\{x=0,v=0\}`$ is absolutely stable. For this system as shown numerically an additive white noise can cause the phase-space density to take a bimodal shape. For the chain of such oscillators the trap effect gives rise to a substantially nonuniform spatial distribution and leads to nonanalytic behavior of the velocity distribution near zero value.

It should be underlined that a possible phase state that could be ascribed in this case to a maximum of the distribution function in the phase space does not match any stationary point of the “regular” or “random” forces.

## Acknowledgements

These investigations were supported in part by RFBR Grants 01-01-00389, 02-02-16537, and Russian Program “Integration”, Project B0056.


text:
# High Resolution Optics & Detector Systems for Hard X-rays

## References

Gaskin, J., Sharma, D., Ramsey, B.D. & Seller, P. 2003, SPIE 5164, in press.

Gubarev, M.V., Ramsey, B.D., & Apple, J. J. SPIE 5165, 17, 2003, In Press.

Sharma, D., Gaskin, J.A., Ramsey, B.D. & Seller, P. 2002, SPIE 4851, 1019.

Ramsey, B. et al. 2002 ApJ, 568, 432.

Ramsey, B.D. et al. 2003, SPIE 5168, 432.

Weisskopf, M., et al. 2002, pasp, 114, 1.