Cosmic ray driven Galactic winds

Recchia, S.; Blasi, P.; Morlino, G.

doi:10.1093/mnras/stw1966

Abstract

The escape of cosmic rays from the Galaxy leads to a gradient in the cosmic ray pressure that acts as a force on the background plasma, in the direction opposite to the gravitational pull. If this force is large enough to win against gravity, a wind can be launched that removes gas from the Galaxy, thereby regulating several physical processes, including star formation. The dynamics of these cosmic ray driven winds is intrinsically non-linear in that the spectrum of cosmic rays determines the characteristics of the wind (velocity, pressure, magnetic field) and in turn the wind dynamics affects the cosmic ray spectrum. Moreover, the gradient of the cosmic ray distribution function causes excitation of Alfvén waves, that in turn determines the scattering properties of cosmic rays, namely their diffusive transport. These effects all feed into each other so that what we see at the Earth is the result of these non-linear effects. Here, we investigate the launch and evolution of such winds, and we determine the implications for the spectrum of cosmic rays by solving together the hydrodynamical equations for the wind and the transport equation for cosmic rays under the action of self-generated diffusion and advection with the wind and the self-excited Alfvén waves.

diffusion, hydrodynamics, cosmic rays, galaxies: haloes

INTRODUCTION

The possibility of a galaxy to launch winds has attracted attention for many different reasons. For instance, star formation is regulated by the amount of gas available, and winds affect the availability of such gas. In fact, galactic models that do not include feedback processes suffer from overpredicting the amount of baryons and star formation rates (Crain et al. 2007; Stinson et al. 2013). Winds also pollute Galactic haloes with hot dilute plasma that may provide an important contribution to the number of baryons in the Universe (Kalberla et al. 1998; Kalberla & Dedes 2008; Miller & Bregman 2013). Such gas might in fact have already been detected (Snowden et al. 1995) in the form of an X-ray emitting plasma with temperature of several million degrees, and possibly associated with a Galactic wind by Breitschwerdt & Schmutzler (1994, 1999) (see also Everett et al. 2008). Finally winds can affect the transport of cosmic rays (CRs) in a galaxy, by advecting them away from their sources.

Galactic winds may be thermally driven, namely powered by core-collapse SNe (see e.g. Chevalier & Clegg 1985) or momentum-driven, powered by starburst radiation (Scoville 2003; Murray, Quataert & Thompson 2005). These two mechanisms of wind launching are thought to be at work in starburst galaxies and galaxies with active nuclei (see e.g. Veilleux, Cecil & Bland-Hawthorn 2005, for a comprehensive review). On the other hand in a galaxy like the Milky Way, winds are unlikely to be due to such processes because thermal and radiation pressure gradients are expected to be too small. A possible exception is the innermost part of the Galactic Centre region where the recently discovered Fermi Bubbles may have their origin through direct bursts from Sgr A* (Cheng et al. 2011; Zubovas, King & Nayakshin 2011) or because of past starburst activities (Lacki 2014). On the other hand, CRs can play an important role in launching winds because of the gradient that their pressure develops as a consequence of the gradual escape of CRs from the Galaxy. The force −∇P_CR associated with such gradient is directed opposite to the gravitational force, and in certain conditions the plasma above and below the disc can be lifted off to form a CR-driven wind. Notice that the gravitational force may be dominated by the dark matter component or the baryonic (gas and stars) components depending on the location. The force exerted by CRs depends in a complicated manner not only on the density of sources of CRs but also on non-linear processes of excitation of Alfvén waves through streaming instability. Both the force induced by CRs on the background plasma and the streaming instability induced by CRs depend on the gradient in CR density. In turn the distribution function of CRs is affected by their transport: diffusion is self-regulated through the production of Alfvén waves, and advection is determined by the velocity of the wind, if any is launched, and by the Alfvén waves’ velocity, directed away from the sources of CRs. This complex interplay makes the problem non-linear and in fact the full problem of calculating at the same time the properties of the wind and the spectrum of CRs as a function of the location in the Galaxy has never been solved before. The first pioneering attempt to describe the hydrodynamics of a CR-driven wind was described in a paper by Ipavich (1975), where the author used a spherically symmetric model of the Galaxy and considered only baryons and stars for the calculation of the gravitational potential. Later Breitschwerdt, McKenzie & Voelk (1991) presented an extensive discussion of the hydrodynamics of CR-driven winds: dark matter was included and a realistic geometry of the wind was considered, in which the launch takes place at some distance from the Galactic disc and proceeds in a roughly cylindrical symmetry out to a distance of about ∼15 kpc, where the flow opens up into a spherical shape. Such a geometry became a milestone for future calculations of CR-driven winds and is adopted also in our paper, as discussed below. The calculations of Breitschwerdt et al. (1991) treated CRs as a fluid, hence no information on the spectrum of CRs was retained. The important role of wave damping in the wind region was also discussed by Breitschwerdt et al. (1991), although only in the simplified case of a spherical outflow. As mentioned above, Breitschwerdt et al. (1991) assumed that the wind is launched some distance away from the disc of the Galaxy and this assumption raises the issue of what happens in the region between the disc and the base of the wind, a problem of both mathematical and physical importance discussed by Breitschwerdt, McKenzie & Voelk (1993) and that will be central to our paper as well.

The paper by Breitschwerdt et al. (1991) represented a milestone in the investigation of CR-driven winds and their calculation was adopted in much of the future literature in the field, including the recent work by Everett et al. (2008) that used the model of a CR-driven wind (including damping) to explain the observed Galactic soft X-ray emission.

The dynamical role of CRs in launching the winds was also studied via purely hydrodynamical simulations (Uhlig et al. 2012; Booth et al. 2013; Salem & Bryan 2014) and through MHD simulations (Peters et al. 2015; Girichidis et al. 2016; Ruszkowski, Yang & Zweibel 2016). These simulations, with their progressive level of sophistication, demonstrated that CRs play an important role in wind launching. Nevertheless, all this bulk of work treated CRs as a fluid, thereby not providing any information on the CR spectrum.

The first calculation of the spectrum of CRs in a Galaxy with a CR-driven wind was made by Ptuskin et al. (1997), where the rotation of the Galaxy over cosmological time-scales was also taken into account Zirakashvili et al. (1996). The authors used a simplified approach to the characteristics of the wind (for instance the assumption that the advection velocity for small heights above the disc scales linearly with distance from the disc) to infer some general implications for the spectrum of CRs. To our knowledge, the paper by Ptuskin et al. (1997) represents the only attempt to account for both a CR-driven wind and its effects of the transport of CRs. Much work has been done in later years on treating the wind in a more realistic way, including the recent time-dependent approach to the wind evolution (Dorfi & Breitschwerdt 2012) and numerical simulations of the wind dynamics (Uhlig et al. 2012) and chemical evolution (Girichidis et al. 2016), but none of these efforts included the kinetic description of CRs.

In this paper we present the first semi-analytical calculation of the hydrodynamics of a CR-driven wind, including the self-generation of Alfvén waves through streaming instability and their damping through non-linear Landau damping (NLLD), and the calculation of the diffusive-convective transport of CRs in such wind. Throughout the paper, we assume that the general topology of the magnetic field lines is such that the transport of the wind and of CRs occurs along the field lines (parallel geometry).

At present, for the sake of simplicity, we decided to leave rotation out of this calculation, but its introduction is straightforward from the technical point of view. The spectrum of CRs at the position of the Earth is calculated for several cases of interest, so as to emphasize that several of the cases that are considered feasible from the point of view of hydrodynamics actually lead to spectra that are not consistent with observations. Special attention is devoted to the discussion of the role of the near disc region, where waves cannot be produced because of ion-neutral damping, and CR scattering, if any is present, is most likely guaranteed by some sort of turbulence directly injected by the same sources of CRs, for instance supernova explosions.

We find that the spectrum of CRs observed at the Earth in the low-energy regime (E ≲ 1 TeV) is mainly affected by the wind launching and evolution, because of the interplay between advection with the wind (and Alfvén waves) and self-generated diffusion. In general, at E ≲ 10 GeV advection is dominant, while a relatively steep spectrum is found at higher energies. For E ≳ 1 TeV the CR transport is dominated by what happens in the near-disc region. In general, this situation may lead to a spectral change in the TeV range, which under some conditions may resemble the hardening observed by PAMELA (Adriani et al. 2011) and AMS-02 (Aguilar et al. 2015), at least qualitatively.

The paper is organized as follows: in Section 2 we describe the basic formalism and the mathematical methods used to compute the structure of the wind in the presence of self-generated waves. The details of how to solve the hydrodynamical equations and the CR transport equation are provided in Sections 3 and 4, respectively. Our results are discussed in Section 5 for our reference model (Section 5.1) and for an alternative scenario devised to illustrate the importance of the near-disc region (Section 5.2). We draw our conclusions in Section 6.

A SEMI-ANALYTICAL APPROACH TO CR-DRIVEN WINDS

The dynamics of winds in the presence of CRs is described by the equations of conservation of mass, momentum and energy for the wind itself and by the transport equation for CRs. The two are coupled to each other in three ways: (1) the gradient of the CR pressure acts as a force on the plasma; (2) Alfvén waves excited by CR streaming are quickly damped through NLLD, thereby resulting in heating of the plasma in the wind; (3) the diffusion term in the CR transport equation leads to an effective contribution to the enthalpy. In the assumption that damping occurs on time-scales much shorter than any other process, definitely justified for NLLD, one can write the equations of hydrodynamics of the wind (see Appendix A for their derivation) as:

\begin{eqnarray} {\nabla }\cdot (\rho {\boldsymbol u})=0, \end{eqnarray}

(1)

\begin{eqnarray} \rho ({\boldsymbol u}\cdot {\nabla }){\boldsymbol u}=-{\nabla }(P_{\rm g}+P_{\rm c})-\rho {\nabla }\Phi , \end{eqnarray}

(2)

\begin{eqnarray} {\boldsymbol u}\cdot {\nabla }P_{\rm g} = \frac{\gamma _{\rm g} P_{\rm g}}{\rho }{\boldsymbol u}\cdot {\nabla }\rho - (\gamma _{\rm g}-1){{\boldsymbol v}_{\rm A}}\cdot {\nabla }P_{\rm c}, \end{eqnarray}

(3)

\begin{eqnarray} {\nabla }\cdot \left[\rho {\boldsymbol u}\left(\frac{u^2}{2}+\frac{\gamma _{\rm g}}{\gamma _{\rm g}-1}\frac{P_{\rm g}}{\rho } + \Phi \right)\right]= -({\boldsymbol u}+{\boldsymbol v}_{\rm A})\cdot {\nabla }P_{\rm c}, \end{eqnarray}

(4)

\begin{eqnarray} {\nabla }\cdot \left[({\boldsymbol u}+{\boldsymbol v}_{\rm A})\frac{\gamma _{\rm c} P_{\rm c}}{\gamma _{\rm c} -1} -\frac{\overline{D}{\nabla }P_{\rm c}}{\gamma _{\rm c}-1}\right]= ({\boldsymbol u}+{\boldsymbol v}_{\rm A})\cdot {\nabla }P_{\rm c}, \end{eqnarray}

(5)

\begin{eqnarray} {\nabla }\cdot {\boldsymbol B}=0 , \end{eqnarray}

(6)

where ρ(z), |${\boldsymbol u}(z)$| and P_g(z) are the gas density, velocity and pressure and |${\boldsymbol B}$| is the magnetic field in the wind, while P_c(z) is the CR pressure and γ_c is the adiabatic index of the CR gas. Notice that γ_c(z) is actually calculated locally (as a function of z) from the distribution function f(p, z) that solves the CR transport equation.

We introduced the Alfvén velocity |${\boldsymbol v}_{\rm A}(z)$|⁠, while Φ(R₀, z) is the gravitational potential of the Galaxy (see Section 2.1). The fact that the wave pressure does not appear in the equations reflects the assumption of fast damping, which results in small wave pressure compared to the gas and CR pressures.

The transport of CRs is described by the advection–diffusion equation:

\begin{eqnarray} {\nabla }\cdot \left[D{\nabla }f\right]-({\boldsymbol u}+{\boldsymbol v}_{\rm A})\cdot {\nabla }f+{\nabla } \cdot ({\boldsymbol u}+{\boldsymbol v}_{\rm A})\frac{1}{3}\frac{\mathrm{\partial} f}{\mathrm{\partial} \ln p}+Q=0,\nonumber\\ \end{eqnarray}

(7)

where |$f({\boldsymbol r},p)$| and |$D({\boldsymbol r},p)$| are the CR distribution function and diffusion coefficient as functions of position |${\boldsymbol r}$| and momentum p. The term Q represents injection of CRs in the Galaxy, which we assume to be limited to the Galactic disc. The average diffusion coefficient that appears in equation (5) is defined as:

\begin{equation} \overline{D} ({\boldsymbol r})= \frac{\int _{0}^{\infty } \text{d}p\ p^{2} T(p) D({\boldsymbol r},p) \nabla f}{\int _{0}^{\infty } \text{d}p \ p^{2} T(p) \nabla f}, \end{equation}

(8)

and T(p) is the kinetic energy of particles with momentum p.

In the following we adopt the same geometry of the wind flow as initially introduced by Breitschwerdt et al. (1991), and used by many authors afterwards. The assumption is that the wind is launched from a surface at distance z₀ above (and below) the disc, and that it retains a roughly cylindrical geometry out to a distance Z_s of the order of the radius of the disc. At larger distances the flow opens into a spherical shape. The surface of the wind is then assumed to be in the form (see Fig. 1):

\begin{equation} A(z)=A_{0} \left[ 1+ \left( \frac{z}{Z_s}\right)^{2}\right], \end{equation}

(9)

which is only function of one spatial coordinate, z. This makes the problem effectively one-dimensional as long as all quantities are assumed to depend only on z. This simplifies the conservation equations that can be easily shown to lead to the following expressions (see Appendix B):

\begin{eqnarray} \rho u A = {const}, \end{eqnarray}

(10)

\begin{eqnarray} AB= {const}, \end{eqnarray}

(11)

\begin{eqnarray} \frac{{\rm d}u}{{\rm d}z}=u\frac{c_{{\ast }}^2\frac{1}{A} \frac{{\rm d}A}{{\rm d}z} -\frac{{\rm d}\Phi }{{\rm d}z}}{u^2-c_{{\ast }}^2}, \end{eqnarray}

(12)

\begin{eqnarray} \frac{{\rm d}P_{\rm g}}{{\rm d}z}=\gamma _{\rm g}\frac{P_{\rm g}}{\rho }\frac{{\rm d}\rho }{{\rm d}z}-(\gamma _{\rm g} -1)\frac{v_{\rm A}}{u}\frac{{\rm d}P_{\rm c}}{{\rm d}z} \end{eqnarray}

(13)

\begin{eqnarray} \frac{{\rm d}P_{\rm c}}{{\rm d}z}=\frac{\gamma _{\rm c} P_{\rm c}}{\rho }\frac{2u+v_{\rm A}}{2(u+v_{\rm A}-\mathcal {D})}\frac{{\rm d}\rho }{{\rm d}z}, \end{eqnarray}

(14)

\begin{eqnarray} c_{{\ast }}^2=\gamma _{\rm g} \frac{P_{\rm g}}{\rho } + \gamma _{\rm c} \frac{P_{\rm c}}{\rho } \left[ 1-(\gamma _{\rm g}-1)\frac{v_{\rm A}}{u} \right]\frac{2u+v_{\rm A}}{2(u+v_{\rm A} -\mathcal {D})}, \end{eqnarray}

(15)

\begin{eqnarray} \mathcal {D}=\frac{\frac{{\rm d}}{{\rm d}z}\left(A\overline{D}\frac{{\rm d}P_{\rm c}}{{\rm d}z}\right)}{A\frac{{\rm d}P_{\rm c}}{{\rm d}z}}. \end{eqnarray}

(16)

One can easily recognize the generalized sound speed c_* introduced by Breitschwerdt et al. (1991), where however two important differences appear: first, the non-adiabatic heating induced by wave damping, which has important implications for the wind launching. Secondly, the effective velocity term |$\mathcal {D}$| which accounts for the finite diffusivity of CRs. In the calculations of Breitschwerdt et al. (1991), it was assumed that |$\overline{D}=0$|⁠, although one can check a posteriori that in general this assumption is not justified, in that there are locations where such term is at least comparable with all others, so that neglecting it leads to an incorrect solution of the problem.

The transport equation of CRs also gets simplified with the flux geometry assumed above:

\begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} z}&& \left[ A(z)D(z,p)\frac{\mathrm{\partial} f(z,p)}{\mathrm{\partial} z}\right] -\left[A(z)U(z)\frac{\mathrm{\partial} f(z,p)}{\mathrm{\partial} z}\right] \nonumber \\ && +\;\frac{{\rm d}[A(z)U(z)]}{{\rm d}z}\frac{1}{3}\frac{\mathrm{\partial} f(z,p)}{\mathrm{\partial} \ln p}+A(z)Q(z,p)=0, \end{eqnarray}

(17)

where we introduced the advection velocity U(z) = u(z) + v_A(z). In general, the diffusion coefficient D(z, p) is assumed to be determined by the process of self-generation of Alfvén waves excited by CRs through streaming instability, although we will consider situations in which we relax this assumption. Since the self-generated perturbations in the magnetic field are relatively weak, one can still use quasi-linear theory to write the diffusion coefficient as:

\begin{equation} D(z,p)=\left.\frac{1}{3}\frac{v(p)r_L(z,p)}{\mathcal {F}}\right|_{k_{\rm res}=1/r_L}, \end{equation}

(18)

where |$\mathcal {F}$| is the normalized energy density per unit logarithmic wavenumber k, calculated at the resonant wavenumber k_res = 1/r_L(p). The local value of |$\mathcal {F}$| is determined by the balance between the CR-driven growth of Alfvén waves and their damping. In the region where the background gas is totally ionized, waves are damped through NLLD at a rate Ptuskin & Zirakashvili (2003):

\begin{equation} \Gamma _{\rm D} = (2c_k)^{-3/2}kv_{\rm A} \mathcal {F}^{1/2}. \end{equation}

(19)

On the other hand the growth occurs at a rate that is given by Skilling (1971):

\begin{equation} \Gamma _{\rm CR}=\frac{16\pi ^2}{3}\frac{v_{\rm A}}{\mathcal {F} B^2}\left[p^4v(p)\left| \frac{\mathrm{\partial} f}{\mathrm{\partial} z} \right| \right]_{p=p_{\rm res}}. \end{equation}

(20)

Equating Γ_D = Γ_CR and using equation (18) one obtains:

\begin{equation} \mathcal {F}(z,p)=2c_k\left[\frac{p^4v(p) \left| \frac{\mathrm{\partial} f}{\mathrm{\partial} z} \right| \frac{16\pi ^2}{3}r_L(z,p)}{B^2(z)} \right]^{2/3}. \end{equation}

(21)

In Sections 3 and 4 we will describe the method adopted to solve the hydrodynamical equations and the transport equation, respectively. Here we describe the way the two pieces are linked together in an iterative procedure. We initiate our computational procedure by solving the equation for the wind structure, by assuming that all quantities are known at the base of the wind, and by requiring that |$\overline{D}=0$|⁠. This first step allows us to derive a guess for the launching velocity of the wind u₀ and the wind properties (velocity, pressure of the gas, CRs and waves). At this point this guess is used as an input for the quantities entering the transport equation (17): most important, the wind velocity u(z) and the Alfvén speed |$v_{\rm A}(z)=B(z)/\sqrt{4\pi \rho (z)}$| are calculated from the output of the hydrodynamical part of the calculation. The solution of equation (17) leads to the spectrum of CRs as a function of location z and momentum p, f(z, p), and to knowing the diffusion coefficient D(z, p) as due to self-generated Alfvén waves. Notice that in general the CR pressure obtained by integrating the distribution function f(z, p) (and especially the CR pressure at the base of the wind) is different from the one assumed in order to determine the wind structure. This leads to the second iteration: the updated CR pressure at the base is used to compute the structure of the wind and the whole calculation is repeated. Moreover, the diffusion coefficient and the distribution function are used to compute γ(z) and |$\overline{D}$|⁠. Notice that the latter quantity is no longer bound to vanish, as assumed in previous calculations by Breitschwerdt et al. (1991). This iterative scheme continues until the spectrum of CRs reaches convergence at all distances and at all momenta. As a consequence the CR pressure at the base converges to the actual solution, so it does the whole structure of the CR induced wind. Clearly, if the goal is to apply this calculation to the case of our own Galaxy at the location of the Sun, R_⊙, then the CR injection is modified so as to reach a solution where the CR pressure is the same as observed. It is very important to stress that this condition alone is not sufficient to guarantee that the CR spectrum observed at the Earth is also consistent with observations: in fact, a generic wind solution leads to spectra of CRs at the Earth location that are quite different from the observed ones, mainly as a result of the strong advection with the wind, which corresponds to exceedingly hard spectra at low energy. We will discuss the qualitative comparison between the wind solutions and the observed CR spectra in Section 5.

Since the properties of the wind depend in a rather sensible way upon the properties of the gravitational potential, in Section 2.1 we describe in detail our assumptions on the contribution of disc and bulge content and dark matter to the gravitational potential.

Galactic gravitational potential

The gravitational potential of the bulge–disc component is modelled following Miyamoto & Nagai (1975):

\begin{equation} \Phi _{\rm B,D}(R_0, z)= - \sum _{i=1}^2 \frac{G M_i}{\sqrt{R_0^2 + \left(a_i + \sqrt{z^2 + b_i^2}\right)^2}}, \end{equation}

(22)

where z and R₀ are the distance from the galactic disc and the galactocentric distance, respectively. Following Breitschwerdt et al. (1991), we assumed for the bulge (M₁, a₁, b₁) = (2.05 × 10¹⁰M_⊙, 0.0 kpc, 0.495 kpc) and (M₂, a₂, b₂) = (2.547 × 10¹¹M_⊙, 7.258 kpc, 0.520 kpc) for the disc.

Dark matter distribution is assumed to be well described by a Navarro–Frenk–White (NFW) profile Navarro, Frenk & White (1996):

\begin{equation} \rho _{\rm DM} = \frac{\rho _{0}}{x(1+x)^{2}}, \end{equation}

(23)

where x = r/r_c and r_c is the radius of the core of the distribution. The two quantities ρ₀ and r_c are calculated by requiring that the total mass of the halo is M_vir = 10¹²M_⊙ Nesti & Salucci (2013) and the energy density of dark matter at the position of the Sun is 0.3 GeV cm⁻³. The dark matter halo is assumed to extend out to a maximum distance that equals the viral radius r_vir ≈ 300 kpc.

The Dark Matter halo potential corresponding to this spatial distribution reads:

\begin{eqnarray} \Phi (x)&= & -4\pi G\rho _s r_s^2\left[ \frac{\ln (1+x)}{x} -\frac{\ln (1+x_{\rm vir})}{x_{\rm vir}}\right]- \frac{GM_{\rm vir}}{r_s}\frac{1}{x_{\rm vir}} \nonumber \\ && \qquad\qquad\qquad\qquad \qquad \qquad \qquad\qquad\quad\!\! (r \le r_{\rm vir}) \nonumber \\ \Phi (x)&= & - \frac{GM_{\rm vir}}{r_s}\frac{1}{x}. \qquad\qquad\qquad \qquad \qquad \qquad (r > r_{\rm vir}). \end{eqnarray}

(24)

A plot of the acceleration associated with the different contributions to the Galactic gravitational potential at the position of the Sun R_⊙ = 8.5 kpc is shown in Fig. 2 as a function of the height z from the disc. One can see that at R_⊙ the potential is dominated by the disc component out to a distance z ∼ 10 kpc, while at larger distances the dark matter contribution becomes more important.

IDENTIFICATION OF THE WIND SOLUTION

Here we discuss the procedure adopted for the determination of the wind solution for a given set of conditions at the base of the wind. We recall that, as discussed above, this calculation is embedded in an iterative computation that is repeated until convergence of the overall structure of the wind and of the CR spectrum is achieved.

We look for a wind solution of equations (10)–(16), namely a solution for the flow velocity u(z) that shows a transition from subsonic (u < c_*) to supersonic (u > c_*) motion, where c_*(z) is the compound sound speed as a function of z.

The topology of the solutions (see also Breitschwerdt et al. 1991) of equation (10)–(16) depends on the nature of the critical points of the wind equation (12), namely the points in which the velocity derivative has zero numerator |$(c_{{\ast }}^2= \frac{{\rm d}\Phi }{{\rm d}z}/\frac{1}{A} \frac{{\rm d}A}{{\rm d}z})$| and/or zero denominator (⁠|$u^2=c_{{\ast }}^2$|⁠), and is schematically represented in Fig. 3. The point for which both numerator and denominator are zero is the critical (sonic) point and it corresponds to the location where the flow velocity equals the compound sound speed, u = c_*. As shown in Fig. 3, there are two curves passing through the critical point. The one corresponding to a subsonic flow at z₀ is the one relevant for our problem. The other one corresponds to an accretion (infall) solution. The curves in the lower branch of solutions shown in Fig. 3 correspond to flows that remain subsonic, the so-called ‘breezes’. The upper branch corresponds to supersonic flow and is physically irrelevant. For both families of curves, there is a point where the numerator of equation (12) vanishes. The other two branches are unphysical and for both there is a point where the denominator of equation (12) vanishes.

For given values of the magnetic field B₀, gas density ρ₀, gas pressure P_g0 and CR pressure P_c0 at z = z₀, we compute the velocity at the base of the wind u₀ that corresponds to the wind solution, namely the value of u₀ for which the flow starts as subsonic and then smoothly becomes supersonic at the critical point z_c (see Fig. 3). Both the location of the sonic point and the launching velocity u₀ are outputs of the calculation. For each iteration, P_c0 and |$\mathcal {D}(z)$| in equation (16) are computed based on the solution of the transport equation in the previous iteration.

From the technical point of view, we start the search for the wind solution by bracketing the range of values of u₀ that may potentially correspond to a wind solution. In order to do so, we impose three conditions: (1) the flow starts as subsonic (c_*0 > u₀); (2) u(z) is a growing function of z at z₀, namely |$\frac{{\rm d}u}{{\rm d}z}\vert _{z_0} >0$|⁠; (3) the solution leads to a final wind velocity at infinity that is physical, namely |$u_{f}^{2}>0$|⁠.

This last point corresponds to imposing energy conservation between z = z₀ and infinity and requiring that at z = ∞ all fluxes vanish with the exception of the kinetic flux associated with the wind bulk motion. This condition selects the solutions that correspond to |$u_{f}^{2}>0$|⁠.

Once a closed interval for u₀, say |$[u_{0,\rm min}, u_{0,\rm max}]$|⁠, has been determined, we numerically integrate the hydrodynamic equations starting at z₀ and for values of u₀ in this range, stepped by Δu₀. In this way we sample the topology of the solutions and can identify the transition between the last breeze and the first unphysical solution by gradually increasing the value of u₀. This procedure can be repeated to narrow down the region where the actual wind solution is located (see Fig. 4).

The procedure of identification of the wind solution also returns the location of the critical point z = z_c. We can then compute both the flow speed and its derivative at the sonic point:

\begin{eqnarray} u^2_{\rm c}= c^2_{{\ast }c}=\frac{{\rm d}\Phi }{{\rm d}z}{\big/} \frac{1}{A} \frac{{\rm d}A}{{\rm d}z}, \end{eqnarray}

(25)

\begin{eqnarray} \left(\frac{{\rm d}u}{{\rm d}z}\right)_{z_{\rm c}}&=& \frac{u_{\rm c}}{2}\left(\frac{1}{A}\frac{{\rm d}A}{{\rm d}z}\right)_{z_{\rm c}} \nonumber \\ &&+\, \sqrt{\left(-\frac{u^2_{\rm c}}{4}\left(\frac{1}{A}\frac{{\rm d}A}{{\rm d}z}\right)^2 + \frac{u^2_{\rm c}}{2}\frac{1}{A}\frac{{\rm d}^2A}{{\rm d}z^2} -\frac{1}{2}\frac{{\rm d}^2\Phi }{{\rm d}z^2} \right)_{z_{\rm c}}}.\nonumber\\ \end{eqnarray}

(26)

The last expression has been obtained by Taylor expanding equation (12) around the critical point. In addition, using the conservation relations (10), equations (11) and (B12) and the fact that at the sonic point u_c = c_*(z_c), we can compute the values of all physical quantities (ρ, B, P_g, P_c) at the sonic point. The wind solution is finally found by integrating the hydrodynamic equations starting at the sonic point towards z = z₀ and towards z → ∞. The hydrodynamic part of our calculation has been checked versus the results of Breitschwerdt et al. (1991) for the cases considered there and for the same set of parameters.

KINETIC CALCULATION

The stationary CR transport equation is as reported in equation (17). Injection of CRs is assumed to take place only in an infinitely thin Galactic disc of radius R_d, so that the injection term is written as Q(z, p) = Q₀(p)δ(z), where

\begin{equation} Q_0(p)=\frac{\mathcal {N}_{\rm SN}(p)\mathcal {R}_{\rm SN}}{\pi R_{\rm d}^2}, \end{equation}

(27)

and |$\mathcal {N}_{\rm SN}(p)$| is the spectrum contributed by an individual source occurring with a rate |$\mathcal {R}_{\rm SN}$|⁠. We have in mind supernova remnants (SNRs) as the sources of Galactic CRs, although the calculations presented here do not depend crucially on such an assumption. The spectrum of each SNR can be written as:

\begin{equation} \mathcal {N}_{\rm SN}(p)=\frac{\xi _{\rm CR}E_{\rm SN}}{I(\gamma )c(mc)^4}\left(\frac{p}{mc} \right)^{-\gamma }, \end{equation}

(28)

where ξ_CR is the CR injection efficiency (typically ∼10 per cent), E_SN is the energy released by a supernova explosion (∼10⁵¹erg), |$\mathcal {R}_{\rm SN}$| is the rate of SN explosions (∼1/30 yr⁻¹), and I(γ) is a normalization factor chosen in such a way that

\begin{equation} \int _{0}^{\infty }\mathcal {N}_{\rm SN}T(p) \, d^3p=\xi _{\rm CR}E_{\rm SN}, \end{equation}

(29)

where T(p) is the kinetic energy of a particle of momentum p.

Equation (17) can be recast in a simpler form by noting that

\begin{equation} UA\frac{\mathrm{\partial} f}{\mathrm{\partial} z}=\frac{\mathrm{\partial} (AUf)}{\mathrm{\partial} z}-f\frac{{\rm d}(AU)}{{\rm d}z}, \end{equation}

(30)

so that equation (17) becomes

\begin{eqnarray} \frac{\mathrm{\partial} }{\mathrm{\partial} z} \left[ AD\frac{\mathrm{\partial} f}{\mathrm{\partial} z} \!-\! UAf\right] \!+\!\frac{{\rm d}[AU]}{{\rm d}z}\left[f\!+\! \frac{p}{3}\frac{\mathrm{\partial} f}{\mathrm{\partial} p}\right] \!+\! AQ_0(p)\delta (z)=0.\nonumber\\ \end{eqnarray}

(31)

In addition,

\begin{equation} \frac{1}{3}\frac{1}{p^3} \frac{\mathrm{\partial} (fp^3) }{\mathrm{\partial} \ln p} = \frac{f}{3}\frac{\mathrm{\partial} \ln (fp^3)}{\mathrm{\partial} \ln p} = \frac{1}{3} \frac{\mathrm{\partial} f}{\mathrm{\partial} \ln p} + f \end{equation}

(32)

so that we can finally rewrite equation (17) as

\begin{equation} \frac{\mathrm{\partial} }{\mathrm{\partial} z} \left[ AD\frac{\mathrm{\partial} f}{\mathrm{\partial} z} -UAf\right] -\frac{{\rm d}[AU]}{{\rm d}z}\frac{f}{3}q +AQ_0(p)\delta (z)=0, \end{equation}

(33)

where we introduced the slope:

\begin{equation} q(z,p)=-\frac{\mathrm{\partial} \ln (fp^3)}{\mathrm{\partial} \ln p}. \end{equation}

(34)

We look for an implicit solution of the transport equation (33) that satisfies the boundary conditions f(0, p) = f₀(p) at the disc and f(p, H) = 0 at some outer boundary H. Note that here H is set just for numerical purposes, and can be also located at spatial infinity. In other words the transport equation in the presence of a wind does not require the existence of a predefined physical halo of given size, in contrast with the standard calculations of CR transport, where such quantity plays a crucial role.

The equation for f₀(p) can be obtained by integrating equation (33) around the disc, namely between z = 0⁻ and z = 0⁺, and using the fact that the problem is symmetric around the disc so that:

\begin{eqnarray*} U(0^+)&=& -U(0^-)=U_0 \nonumber\\ \left. \frac{{\rm d}U}{{\rm d}z}\right|_0&=&2U_0\delta (z)\nonumber\\ f(0^+)&=&f(0^-)=f_0(p)\nonumber\\ \left. \frac{{\rm d}U}{{\rm d}z}\right|_{0^+}&=& -\left. \frac{{\rm d}U}{{\rm d}z}\right|_{0^-}. \end{eqnarray*}

Integration around the disc leads to

\begin{equation} 2A_0D_0 \left. \frac{{\rm d}U}{{\rm d}z}\right|_{0^+} - 2 A_0U_0\frac{f_0}{3}\tilde{q}_0 + A_0Q_0(p)=0, \end{equation}

(35)

where we used equations (34) and (32) and introduced

\begin{equation} \tilde{q}_0(p)=-\frac{\mathrm{\partial} \ln f_0(p)}{\mathrm{\partial} \ln p}. \end{equation}

(36)

Integrating equation (33) between z = 0⁺ and a generic height z one gets the equation:

\begin{eqnarray} \frac{\mathrm{\partial} f}{\mathrm{\partial} z}(z,p)&=& \frac{U(z)}{D(z,p)}f(z,p)+\,\frac{\mathcal {G}(z,p)}{A(z)D(z,p)}\nonumber \\ &&-\;\frac{A_0 Q_0(p)}{2A(z)D(z,p)} -\frac{U_0A_0f_0(p)\overline{q}_0(p)}{A(z)D(z,p)}, \end{eqnarray}

(37)

where

\begin{eqnarray} \mathcal {G}(z,p)=\int _0^z \text{d}z\; \frac{{\rm d}(\text{AU})}{{\rm d}z}\frac{f}{3}q \end{eqnarray}

(38)

\begin{eqnarray} \overline{q}_0(p)=1-\frac{\tilde{q}_0}{3}. \end{eqnarray}

(39)

Integration of this equation between z and H with the boundary condition f(H, p) = 0 leads to the following implicit solution for f(z, p):

\begin{equation} f(z,p)={\int _z}^H \frac{\text{d}z^{\prime }}{A(z^{\prime })D(z^{\prime },p)} G(z^{\prime },p) \exp ^{-{\int _z}^{z^{\prime }}\text{d}z^{\prime \prime }\frac{U(z^{\prime \prime })}{D(z^{\prime \prime }, p)}}, \end{equation}

(40)

where

\begin{equation} G(z,p)=\frac{A_0 Q_0(p)}{2} + A_0U_0f_0(p)\bar{q}_0 - \mathcal {G}(z,p). \end{equation}

(41)

This rewriting of the transport equations hides the non-linearity of the problem in the function G(z, p), which depends on f(z, p) and on the CR diffusion coefficient D(z, p) (see equation 18). The solution of equation (40) is computed with an iterative procedure.

The procedure reaches convergence when, for a given iteration j, f ^j(z, p) and f ^{j − 1}(z, p) are close to each other within a desired precision. Note that the advection velocity U(z) is computed from the hydrodynamic equations and is fixed while iterating upon the distribution function f ^j.

RESULTS

In this section we illustrate some selected cases of CR induced winds, aimed at addressing different issues that arise in this type of problem. The first case is what we will refer to as our reference case, namely a case that illustrates the most basic characteristics of a CR-induced wind, with a minimal number of physical parameters introduced. The wind is launched very close to the disc of the Galaxy, and we consider the specific situation expected at the solar radius, namely at the Galactocentric distance R = R_⊙ = 8.5 kpc, where the Sun is located. This information is crucial in that it determines the gravitational potential that the wind has to fight against. The second model considered below is that in which a wind is launched at some distance from the Galactic disc, while particle transport in the near-disc region is assigned. This latter situation is physically motivated by the fact that ion-neutral damping is expected to damp self-generated Alfvénic turbulence within ∼ 0.5–1 kpc from the disc, because of the presence of a substantial amount of neutral hydrogen. We will show that the consequences of this setup for the spectrum of CRs observed at the Earth are very prominent.

Reference case

Our reference case corresponds to launching a CR-induced wind at a Galactocentric distance R = R_⊙ = 8.5 kpc. The base of the wind is assumed to be at the edge of the Galactic disc, z₀ = 100 pc, where we assumed that the ionized gas has a density n₀ = 0.003 cm^{− 3} and the magnetic field is |$B_{0}=1.5\, \mu \text{G}$| (to be interpreted as the component of the field along the z direction). The CR pressure at z₀ is taken to equal the observed CR pressure, P_c0 = 6 × 10^{− 13} erg cm^{− 3}. We solve simultaneously the hydrodynamical equations for the wind and the transport equation for CRs, with self-generated diffusion and advection taken into account. In order to get the desired CR pressure at the Sun's location (see above) we are bound to take

\begin{equation*} \frac{\xi _{\rm CR}}{0.1} \frac{\mathcal {R}_{\rm SN}}{1/30 \rm \, yr^{-1}} \approx 1.8, \end{equation*}

for an injection spectrum with slope γ = 4.3. The density of plasma in the wind and the temperature of the wind are shown in Fig. 5. The temperature of the wind is maintained by the continuous damping of wave energy into thermal energy of the gas. In Fig. 6 we show the wind velocity u(z) (green dashed line), the Alfvén speed v_A(z) (blue dotted line) and the sound speed c_*(z) as functions of the distance from the Galactic disc. The wind is launched with a speed of 41 km s⁻¹ and asymptotically reaches a speed of 353 km s⁻¹, while it becomes supersonic at ∼15 kpc.

Figure 1.

Flux-tube geometry for the magnetic field.

Open in new tab Download slide

Figure 2.

Gravitational acceleration as due to the three components: the bulge, the disc and the dark matter halo. All curves are computed at the Sun's position (R_⊙ = 8.5 kpc).

Open in new tab Download slide

Topology of the solutions of the hydrodynamics equations: z0 indicates the base of the wind, while zc is the critical point, where u(zc) = c*(zc). The (red) solid line represents a wind solution that starts as subsonic and becomes supersonic at the critical point. The branches of the solutions that are unphysical are labelled as such.

Figure 3.

Topology of the solutions of the hydrodynamics equations: z₀ indicates the base of the wind, while z_c is the critical point, where u(z_c) = c_*(z_c). The (red) solid line represents a wind solution that starts as subsonic and becomes supersonic at the critical point. The branches of the solutions that are unphysical are labelled as such.

Open in new tab Download slide

Schematic picture of the solution method for the hydrodynamic equations: the wind curve lies between a breeze, obtained by integrating with u0 = u0B, and an unphysical solution, obtained by integrating with u0 = u0U, where u0B < u0U. The location of the sonic point zc can be estimated as the location of the maximum of the breeze. The small box shows a zoomed region around the sonic point.

Figure 4.

Schematic picture of the solution method for the hydrodynamic equations: the wind curve lies between a breeze, obtained by integrating with u₀ = u_0B, and an unphysical solution, obtained by integrating with u₀ = u_0U, where u_0B < u_0U. The location of the sonic point z_c can be estimated as the location of the maximum of the breeze. The small box shows a zoomed region around the sonic point.

Open in new tab Download slide

Density (red solid line) in units of 10− 3 cm− 3 and temperature (greed dashed line) in units of 106 K for the wind solution in the reference case.

Figure 5.

Density (red solid line) in units of 10^{− 3} cm^{− 3} and temperature (greed dashed line) in units of 10⁶ K for the wind solution in the reference case.

Open in new tab Download slide

Figure 6.

Wind velocity (red solid line), Alfvén velocity (green dashed line) and sound speed (blue dotted line) as a function of the height z above the disc, for the reference case.

Open in new tab Download slide

For the sake of future discussion it is important to notice here that the CR advection velocity at the base of the wind is dominated by the Alfvén speed and that the latter is non-zero at the base of the wind (because of a finite density and magnetic field).

The pressure of CRs as a function of the distance z from the disc is shown in Fig. 7, together with the gas pressure and the CR pressure as derived from the kinetic calculation: the fact that the latter is basically overlapped to P_c(z) as derived from hydrodynamics shows that the system of equations (hydro plus kinetic) reached convergence (with accuracy of ∼0.1 per cent at heights z close to the base of the wind and ∼10–20 per cent at z close to the outer edge of the box used for numerical computation; here all quantities, with the exception of the wind velocity, are bound to vanish).

The gas temperature obtained for z ≲ 10 kpc (see Fig. 5) suggests that cooling may be important. At such temperature the gas is cooled by the emission of forbidden lines and soft X-rays (Dalgarno & McCray 1972) but, on the other hand, it is also heated by SN explosions through the injection of hot gas and magnetic turbulence that will eventually decay into thermal energy. Recently, it was suggested that also Coulomb losses of CRs themselves might be a substantial source of heating (Walker 2016). The relevance of cooling was already recognized by Breitschwerdt et al. (1991) which, nevertheless, assumed that energy balance exists between heating and cooling processes. Which process dominates the heating is at the moment unclear but observations reveal that the temperature in the halo reaches ∼10⁵–10⁶ K (see e.g. Miller & Bregman 2013), hence supporting the idea that either the cooling is negligible or it is balanced by heating processes. In this paper we decided to neglect the role of cooling, deferring a more complete study to a future work.

Figure 7.

CR pressure (red solid line) and gas pressure (blue dotted line) in the reference case. The green dashed line shows the CR pressure as obtained from the solution of the transport equation. The same calculation also returns the wave pressure P_w(z) (pink dotted line) that is assumed to vanish in the hydrodynamical equation.

Open in new tab Download slide

The spectrum of CRs as a function of particle momentum is reported in Fig. 8 for several distances from the Galactic plane. The most striking feature of the spectrum at the disc is the pronounced spectral hardness: the spectrum that should be measured at the Earth location at momenta below ∼1 TeV/c is roughly ∝ p^−4.4, only slightly steeper than the injection spectrum (Q(p) ∝ p^−4.3). This finding reflects the fact that the Alfvén speed at the base of the wind is very large and dominates CR advection up to high energies. One could argue that a larger value of the density at the base of the wind would make the Alfvén speed smaller, but in that case two problems appear: (1) it may become harder if not impossible to launch the wind because of excessive baryonic load. In other words, for a given CR pressure at the base, there may be cases in which the wind is not launched. (2) When the wind is in fact launched, its initial velocity may easily be super-Alfvénic, so that again advection dominates up to relatively high energies.

Figure 8.

CR spectrum, f(p) × p^4.7, in the reference case, at different locations in the wind: base of the wind (red solid line), z = 10 kpc (green dashed line), z = 50 kpc (blue short-dashed line) and z = 100 kpc (pink dotted line).

Open in new tab Download slide

This result on the spectrum of Galactic CRs appears to be at odds with previous results by Ptuskin et al. (1997), where a toy model for the velocity scaling with z provided different results, of apparently easy interpretation. Since the argument of Ptuskin et al. (1997) is rather simple, we report it here and we explain why this simple approach does not apply to realistic cases of CR-driven winds. The basic assumption of Ptuskin et al. (1997) is that the advection velocity (dominated by the Alfvén speed) scales approximately linear with z, v_A ∼ ηz. Now, it is easy to imagine that while the advection velocity increases with z, it reaches a critical distance, s_*, for which advection dominates upon diffusion. This happens when

\begin{equation} \frac{s_{{\ast }}^{2}}{D(p)}\approx \frac{s_{{\ast }}}{v_{\rm A}(s_{{\ast }})} \, \Rightarrow \, s_{{\ast }}(p) \propto D(p)^{1/2}, \end{equation}

(42)

where we used the assumption of linear relation v_A ∼ ηz. Now, when diffusion dominates, namely when z ≲ s_*(p), one can neglect the advection terms and make the approximate statement (as in the standard diffusion model) that |$D(p)\frac{\mathrm{\partial} f}{\mathrm{\partial} z}|_{z=0}\approx -Q_{0}(p)/2 \propto p^{-\gamma }$|⁠. Now, using equations (21) and (18) one can easily show that D(p) ∝ p^2γ−7 (for a p⁻⁴ injection one gets a linear scaling of the diffusion coefficient with momentum). The quantity s_*(p) plays the role of the size of the diffusion volume and one can show that, similar to a leaky box-like model, the equilibrium spectrum in the disc is

\begin{equation} f(p)\sim \frac{Q_{0}(p)}{s_{{\ast }}(p)}\frac{s_{{\ast }}^{2}}{D(p)} \sim Q_{0}(p)\, D(p)^{-1/2} \sim p^{-2\gamma + 7/2}. \end{equation}

(43)

For injection Q₀(p) ∼ p^−4.3 one infers an equilibrium spectrum f₀(p) ∼ p^−5.1 (notice the contrast with the standard diffusion model within a given halo of size H, which predicts f₀(p) ∼ Q₀(p)H/D(p)).

The problem with this argument, put forward by Ptuskin et al. (1997), is that it is strictly valid only when the advection velocity vanishes while approaching the base of the wind. One can see from equation (35) that the assumption that diffusion is the dominant process at z → 0 holds true only if U₀ = 0. The solution of the combined hydrodynamical equations and transport equation of CRs as presented above clearly shows that this is not the case. At low energies the slope of the CR spectrum at the base of the disc, as shown in Fig. 8, is ∼4.3 (see also Fig. 14, lower panel), basically the same as the injection spectrum. The same point can be also made by plotting s_*(p) (Fig. 9) and the diffusion coefficient D(p) (Fig. 10): the simple scaling s_*(p)∝D(p)^1/2 can be easily seen to be not satisfied by the actual solution of the problem. Notice that s_*(p) becomes larger than Z_s, the location where there is a transition to a spherical-like flow, at p ∼ 1 TeV/c. At energies much larger than this one can assume that the advection velocity tends to a constant, u_f. The equilibrium spectrum observed at the Earth can then be written as |$f_{0}(p)\sim Q_{0}(p)\pi R_{\rm d}^{2}/(u_{f} s_{{\ast }}^{2})$|⁠, while s_*∝D(p), so that the equilibrium spectrum has a slope −3γ + 7. This effect corresponds to a steepening of the equilibrium spectrum at the transition energy, which for the values of the parameters used in Fig. 8 corresponds to about ∼1 TeV.

Figure 9.

Effective boundary, s_*(p), between the diffusion-dominated and the advection-dominated region of the wind in the reference case.

Open in new tab Download slide

Figure 10.

Self-generated diffusion coefficient, D(p), in the reference case at different locations in the wind: base of the wind (red solid line), z = 10 kpc (green dashed line), z = 50 kpc (blue short-dashed line) and z = 100 kpc (pink dotted line). The transition from the cylindrical to the spherical geometry of the wind flow is clearly visible in the momentum dependence of D(p).

Open in new tab Download slide

Even qualitatively one can see that the CR spectrum in Fig. 8 is quite different from the one actually observed at the Earth: it is harder than observed at low energies, and it is softer than observed at high energies, even though the pressure carried by these CRs is as observed.

This example clearly shows that it is possible to construct solutions of the hydrodynamical equations that correspond to CR-driven winds, with pressures at the base of the wind that are compatible with observations and yet leading to CR spectra that are not compatible with the CR spectra observed at the Earth. In particular, the basic wind model discussed in this section does not lead to any hardening of the spectrum at high energy, hence it is not possible to fit spectral hardening such as the ones observed by PAMELA (Adriani et al. 2011) and AMS-02 (Aguilar et al. 2015). In Section 5.2, we discuss a situation in which this conclusion may not apply.

The importance of the near-disc region

As already pointed out in the original work on CR-driven winds Breitschwerdt et al. (1991) the region close to the disc may be plagued by severe ion-neutral damping which suppresses the generation of Alfvén waves Kulsrud & Cesarsky (1971). Since waves provide the coupling between CRs and the ionized plasma, their severe damping leads to quenching of the wind. In fact, ion-neutral damping was recognized as a hindering factor for diffusion even in the absence of winds Skilling (1971); Holmes (1975): in these pioneering papers, the near disc region was assumed to be wave-free and the propagation of CRs in that region was taken to be ballistic. Diffusion in the outer halo, where the density of neutral hydrogen is expected to drop and the role of ion-neutral damping to become negligible, was considered as the actual diffusion region. One could, however, speculate that some type of turbulence may be maintained in the near-disc region, perhaps due to SN explosions themselves, though the waves may be considered to be isotropic, so that the effective Alfvén speed vanishes.

In this section we discuss a scenario for the wind launching constructed in the following manner: the wind is assumed to be launched at a distance z₀ = 1 kpc from the disc and the near-disc region (|z| < z₀) is assumed to be characterized by a given diffusion coefficient, in the following form:

\begin{equation} D(p) = 3\times 10^{28} \left(\frac{p}{3 m_{p} c}\right)^{1/3} \ \rm cm^{2} \, s^{-1}, \end{equation}

(44)

and by an Alfvén velocity v_A = 0. At z ≥ z₀, namely in the wind region, the diffusion coefficient is calculated as due to self-generation through streaming instability, saturated by NLLD, as in the reference case (Section 5.1).

It is important to emphasize that the near-disc region is crucial to establish a connection between the sources and the wind region. From the mathematical point of view this is evident from equation (35), where the solution of the transport equation in the disc, f₀(p), is related to the injection rate through the value of the advection velocity at z = 0. When such advection velocity is non-zero, there is always a range of values of the particle momentum (at low momenta) where advection is more important than diffusion and the equilibrium spectrum turns out to have roughly the same slope as the injection spectrum. In the case that we consider in this section, we are assuming that the near-disc region is characterized by a vanishing advection velocity and finite assigned diffusion coefficient. From the technical point of view, the computation is the same as described above. The only minor difference is that, in order to avoid discontinuities in the advection velocity, that is zero in the near disc region but is not zero at z = z₀, we assume that both the wind velocity u(z) and the Alfvén speed v_A(z) have a low z cutoff at z < z₀, so that both velocities drop to zero in a continuous manner in the near disc region. We checked that the details of such assumption do not have serious implications for the solution of the problem, provided that the velocity drops to zero fast enough inside the near-disc region. It is also worth stressing that, contrary to the reference case illustrated in Section 5.1, the CR pressure at the base of the wind (z = z₀) does not correspond to the CR pressure measured at the Earth (z = 0). The criterion for convergence is still that the pressure at z = 0 equals the observed CR pressure at the Earth location (clearly this would be different if we were interested at a different Galactocentric distance). In order to recover the observed CR pressure in the Galactic disc at the Sun's location, in the model discussed here we need to require that:

\begin{equation*} \frac{\xi _{\rm CR}}{0.1} \frac{\mathcal {R}_{\rm SN}}{1/30 \, \rm yr^{-1}} \approx 0.8. \end{equation*}

The density and temperature of the wind, limited to the region z ≥ z₀ = 1 kpc where the wind can be launched, are plotted in Fig. 11. The corresponding wind velocity, Alfvén speed and sound speed in the wind region are plotted in Fig. 12: the advection velocity is again dominated by the Alfvén velocity at the base of the wind. The wind becomes supersonic at z_c ∼ 20 kpc, and eventually reaches u_f ∼ 400 km/s at large distances from the disc. The pressures of the gas and CRs are shown in Fig. 13, where the gas pressure is limited to the wind region while the CR pressure extends to the whole diffusion region. We also plotted there the CR pressure that is returned by the kinetic calculation to show that with good accuracy it coincides with the one derived from the hydrodynamical calculation. The wave pressure as derived from the kinetic calculation is also shown: the plateau at z < z₀ has been estimated from the assigned diffusion coefficient. The fact that P_w(z) is always much smaller than all other terms justifies the assumption of instantaneous damping (P_w = 0) in the hydrodynamical computation.

Density (red solid line) in units of 10− 3 cm− 3 and temperature (greed dashed line) in units of 106 K limited to the region where the wind is launched, z ≥ z0 = 1 kpc.

Figure 11.

Density (red solid line) in units of 10^{− 3} cm^{− 3} and temperature (greed dashed line) in units of 10⁶ K limited to the region where the wind is launched, z ≥ z₀ = 1 kpc.

Open in new tab Download slide

Figure 12.

Wind velocity (red solid line), Alfvén velocity (green dashed line) and sound speed (blue dotted line) as a function of the height z above the disc, limited to the region where the wind is launched, z ≥ z₀ = 1 kpc.

Open in new tab Download slide

Figure 13.

CR pressure and gas pressure in the case of Kolmogorov turbulence in the near-disc region. The green dashed line shows the CR pressure as obtained from the solution of the transport equation. The same calculation also returns the wave pressure P_w(z) that is assumed to vanish in the hydrodynamical equation.

Open in new tab Download slide

Using the approach described in Sections 2, 3 and 4, we also calculated the spectrum of Galactic CRs: in the upper panel of Fig. 14 we show the CR distribution in the Galactic disc for the reference case compared to the case where the wind is launched at a finite distance z₀ from the disc. One can immediately see the dramatic difference that the near-disc region makes in terms of CR spectrum observed at the Earth (or, for that matter, anywhere). In the second case there are at least two spectral breaks that can be identified and explained in a reasonably easy way: the low-energy part of the spectrum is affected by advection (below ∼10 GeV) as well as from self-generation of waves at z ≥ z₀ (at 10 GeV ≲ E ≲ 1000 GeV). At energies higher than ∼1 TeV the spectrum hardens. All these changes of slope are illustrated more clearly in the lower panel of Fig. 14, where we show the slope of the CR distribution in the disc for the reference case (red solid line) and for the case with a wind launched at a finite distance z₀ from the disc (green dashed line).

Figure 14.

Upper panel: spectrum of CRs in the disc of the Galaxy at the position of the Sun for the reference case (red solid line) and for the case in which the wind is launched at z₀ = 1 kpc from the disc and there is a near-disc region where the diffusion coefficient is assigned and advection is absent (green dashed line). The lower panel shows the corresponding slopes.

Open in new tab Download slide

The hard spectrum (slope ∼4.5 below ∼10 GeV reflects advection and is similar to the effect already seen in non-linear models of galactic CR propagation (Aloisio, Blasi & Serpico 2015, see for instance). On the other hand, the spectrum has a continuous steepening at higher energies up to ≲1 TeV: this effects is also similar to the one claimed by Aloisio et al. (2015) and reflects the very rapid energy dependence of the self-generated diffusion coefficient. Finally, the slope of the spectrum becomes ∼4.63 at E ≳ 1 TeV, because of the assumption of Kolmogorov turbulence in the near disc region (z ≤ z₀). This hardening is qualitatively similar to the one found by Tomassetti (2012) in a simple two-zone model of the Galactic halo: the low-energy part is sensitive to the far zone, while the high-energy part is sensitive to the near zone, where the assumption of Kolmogorov turbulence leads to an expected slope 4.3(injection) + 1/3(Kolmogorov) ≃ 4.63. The important point we wish to make here is that the diffusion coefficient that determines the low-energy behaviour of the bulk of CRs is not assumed here, but rather derived from the condition of self-generation of waves and advection with the CR-induced Galactic wind, as calculated above.

Although here we did not make any attempt to fit our curves to real data, it is interesting to notice that qualitatively the scenario presented in this section (contrary to the reference model in Section 5.1) has the same feature as the spectrum of CR protons as recently measured by PAMELA (Adriani et al. 2011) and AMS-02 (Aguilar et al. 2015).

CONCLUSIONS

We presented a semi-analytical calculation of the structure of a CR-driven Galactic wind that returns the hydrodynamical structure of the wind (velocity and pressure of the plasma), and at the same time the spectrum of CRs at any location in the disc or in the wind region. The transport of CRs is described as advective and diffusive: advection occurs at a local speed that is the sum of the Alfvén speed and the wind speed at that point. Diffusion is due to scattering off self-generated turbulence, as generated by the same CRs due to streaming instability. The non-linear nature of the problem is manifold: the wind is driven by the CR gradient in the z direction, perpendicular to the Galactic disc, which results from the gradual escape of CRs from the Galaxy. The force −∇P_CR acts in the direction opposite to gravity, so that, under certain conditions, it can lift ionized gas and launch a wind. On the other hand, the amount of waves generated by CRs depends on their transport, which in turn depends on the wind speed. The diffusion coefficient is calculated in a self-consistent manner from the local spectrum of waves excited by streaming instability and damped through NLLD.

Qualitatively, we confirm earlier findings that CRs can drive winds with asymptotic velocities of the order of several hundred km s⁻¹. Contrary to previous calculations, we also derived the spectrum of CRs as observed at the Earth (or any other location): for most of the parameter space for which a CR-driven wind is launched, the spectrum is in disagreement with observations: the main reason for such a finding is that the CR advection with the wind is very strong and leads to hard spectra that are unlike the one observed at the Earth. However, this conclusion is found to depend in a critical way on the conditions in the near-disc region (z ≲ 1 kpc). The conditions in such region affect the high-energy behaviour of the spectrum, and change the launching speed of the wind at z ∼ z₀ ∼ 1 kpc, so that the low-energy region of the spectrum is also changed. For certain assumptions on the diffusion coefficient in the near-disc region, the predicted CR spectrum may resemble the observed one: as an instance we showed one case in which the spectrum is relatively hard at E ≲ 10 GeV (as previously found by Aloisio et al. 2015), steeper for 10 ≲ E ≲ 1000 GeV, and show a hardening at higher energies as recently claimed by Adriani et al. (2011) and Aguilar et al. (2015). From the physical point of view, we expect that the region at z ≲ 1 kpc is rich in neutral gas, which may damp Alfvén waves through ion-neutral damping. This implies that the coupling between CRs and gas may be exceedingly weak in this region so as to inhibit the launch of the wind. In the near-disc region, diffusion, if any is present, should be guaranteed by some source other than CR self-generation.

The calculation presented here is very general and can be applied to a variety of conditions. This will allow us to explore the conditions in which a CR-driven wind could develop with a CR spectrum at the Earth resembling the observed one. This type of calculations is also useful to develop possible diagnostics of the existence of CR-driven winds in our own Galaxy as well as in others.

The authors are glad to acknowledge a continuous ongoing collaboration with the rest of the Arcetri high-energy Group. This work was partially funded through Grant PRIN-INAF 2012.

REFERENCES

Adriani

O.

et al.

2011

Science

332

69

Crossref

Search ADS

PubMed

Aguilar

M.

et al.

2015

Phys. Rev. Lett.

114

171103

Crossref

Search ADS

PubMed

Aloisio

R.

Blasi

P.

Serpico

P. D.

2015

A&A

583

A95

Crossref

Search ADS

Booth

C. M.

Agertz

O.

Kravtsov

A. V.

Gnedin

N. Y.

2013

ApJ

777

L16

Crossref

Search ADS

Breitschwerdt

D.

Schmutzler

T.

1994

Nature

371

774

Crossref

Search ADS

Breitschwerdt

D.

Schmutzler

T.

1999

A&A

347

650

Breitschwerdt

D.

McKenzie

J. F.

Voelk

H. J.

1991

A&A

245

79

Breitschwerdt

D.

McKenzie

J. F.

Voelk

H. J.

1993

A&A

269

54

Cheng

K.-S.

Chernyshov

D. O.

Dogiel

V. A.

Ko

C.-M.

Ip

W.-H.

2011

ApJ

731

L17

Crossref

Search ADS

Chevalier

R. A.

Clegg

A. W.

1985

Nature

317

44

Crossref

Search ADS

Crain

R. A.

Eke

V. R.

Frenk

C. S.

Jenkins

A.

McCarthy

I. G.

Navarro

J. F.

Pearce

F. R.

2007

MNRAS

377

41

Crossref

Search ADS

Dalgarno

A.

McCray

R. A.

1972

ARA&A

10

375

Crossref

Search ADS

Dorfi

E. A.

Breitschwerdt

D.

2012

A&A

540

A77

Crossref

Search ADS

Everett

J. E.

Zweibel

E. G.

Benjamin

R. A.

McCammon

D.

Rocks

L.

Gallagher

J. S.

2008

ApJ

674

258

Crossref

Search ADS

Girichidis

P.

et al.

2016

ApJ

816

L19

Crossref

Search ADS

Holmes

J. A.

1975

R. Astron. Soc.

170

251

Crossref

Search ADS

Ipavich

F. M.

1975

ApJ

196

107

Crossref

Search ADS

Jones

T. W.

1993

ApJ

413

619

Crossref

Search ADS

Kalberla

P. M. W.

Dedes

L.

2008

A&A

487

951

Crossref

Search ADS

Kalberla

P. M. W.

Westphalen

G.

Mebold

U.

Hartmann

D.

Burton

W. B.

1998

A&A

332

L61

Kulsrud

R. M.

Cesarsky

C. J.

1971

Astrophys. Lett.

8

189

Lacki

B. C.

2014

MNRAS

444

L39

Crossref

Search ADS

Miller

M. J.

Bregman

J. N.

2013

ApJ

770

118

Crossref

Search ADS

Miyamoto

M.

Nagai

R.

1975

PASJ

27

533

Murray

N.

Quataert

E.

Thompson

T. A.

2005

ApJ

618

569

Crossref

Search ADS

Navarro

J. F.

Frenk

C. S.

White

S. D. M.

1996

ApJ

462

563

Crossref

Search ADS

Nesti

F.

Salucci

P.

2013

J. Cosmol. Astropart. Phys.

7

016

Crossref

Search ADS

Peters

T.

et al.

2015

ApJ

813

L27

Crossref

Search ADS

Ptuskin

V. S.

Zirakashvili

V. N.

2003

A&A

403

1

Crossref

Search ADS

Ptuskin

V. S.

Voelk

H. J.

Zirakashvili

V. N.

Breitschwerdt

D.

1997

A&A

321

434

Ruszkowski

M.

Yang

H.-Y. K.

Zweibel

E.

2016

ApJ

preprint (arXiv:1602.04856)

Salem

M.

Bryan

G. L.

2014

MNRAS

437

3312

Crossref

Search ADS

Scoville

N.

2003

J. Korean Astron. Soc.

36

167

Crossref

Search ADS

Skilling

J.

1971

ApJ

170

265

Crossref

Search ADS

Snowden

S. L.

et al.

1995

ApJ

454

643

Crossref

Search ADS

Stinson

G. S.

Brook

C.

Macciò

A. V.

Wadsley

J.

Quinn

T. R.

Couchman

H. M. P.

2013

MNRAS

428

129

Crossref

Search ADS

Tomassetti

N.

2012

ApJ

752

L13

Crossref

Search ADS

Uhlig

M.

Pfrommer

C.

Sharma

M.

Nath

B. B.

Enßlin

T. A.

Springel

V.

2012

MNRAS

423

2374

Crossref

Search ADS

Veilleux

S.

Cecil

G.

Bland-Hawthorn

J.

2005

ARA&A

43

769

Crossref

Search ADS

Walker

M. A.

2016

ApJ

818

23

Crossref

Search ADS

Zirakashvili

V. N.

Breitschwerdt

D.

Ptuskin

V. S.

Voelk

H. J.

1996

A&A

311

113

Zubovas

K.

King

A. R.

Nayakshin

S.

2011

MNRAS

415

L21

Crossref

Search ADS

APPENDIX A: DERIVATION OF THE HYDRODYNAMIC EQUATIONS

In this section we derive the hydrodynamic equations used throughout this paper. They describe the macroscopic behaviour of the thermal plasma and CRs and their mutual interactions. As we already pointed out in Section 2, the wave component is not present in these equations because we are assuming that all Alfvén waves (generated through CR streaming instability) are locally dumped on short time-scales, so that the energy fed by CRs into waves just results in the heating of the gas. Here, we are also neglecting any external loss of mass, momentum and energy. A general form of those equations, with and without damping of waves, is reported in Breitschwerdt et al. (1991).

Gas mass conservation and equation of motion

In absence of any mass and momentum external loss, the time-dependent gas mass conservation and equation of motion read

\begin{equation} \frac{D\rho }{Dt}=-\rho {\nabla }\cdot {\boldsymbol u}, \end{equation}

(A1)

\begin{equation} \rho \frac{D{\boldsymbol u}}{Dt}=-{\nabla }(P_{\rm g} +P_{\rm c}) - \rho {\nabla }\Phi , \end{equation}

(A2)

where ρ, |${\boldsymbol u}$|⁠, P_g, P_c and Φ are the gas density, velocity and pressure, the CR pressure and the galactic gravitational potential respectively. |$D/Dt = \mathrm{\partial} /\mathrm{\partial} t + {\boldsymbol u}\cdot {\nabla }$| is the convective derivative.

Gas internal energy and total energy

The internal energy of the gas in a volume V is given by

\begin{equation} \epsilon _{\rm g}=\frac{Pg V}{\gamma _{\rm g} -1} \end{equation}

(A3)

and its time evolution is given by

\begin{equation} \frac{D\epsilon _{\rm g}}{Dt}= -Pg\frac{DV}{Dt}- {\boldsymbol v}_{\rm A}\cdot {\nabla }P_{\rm c} V \end{equation}

(A4)

where on the right hand side we find the work done on the gas when V changes and the energy input due to damping of the self-generated Alfvén waves. In fact, as it will be shown in Section A4, the hydrodynamic counterpart of the CR streaming instability is given by the term |${\boldsymbol v}_{\rm A}\cdot {\nabla }P_{\rm c}$|⁠, where |${\boldsymbol v}_{\rm A}$| is the Alfvén speed and P_c the CR pressure. In terms of ρ, and using equation (A1), we get

\begin{equation} \frac{D}{Dt}\left(\frac{Pg}{\gamma _{\rm g} -1}\frac{1}{\rho } \right)= -\frac{P_{\rm g}}{\rho }{\nabla }\cdot {\boldsymbol u} - \frac{1}{\rho } {\boldsymbol v}_{\rm A}\cdot {\nabla }P_{\rm c}. \end{equation}

(A5)

The gas total energy per unit volume is the sum of the kinetic and internal energy per unit volume

\begin{equation} \epsilon = \frac{1}{2}\rho u^2 +\frac{Pg}{\gamma _{\rm g} -1}. \end{equation}

(A6)

Making use of equations (A1), (A2) and (A5), we get the equation for the gas total energy

\begin{eqnarray} &&{\frac{D}{Dt}\left(\frac{1}{2}\rho u^2 +\frac{Pg}{\gamma _{\rm g} -1}\right) = -{\nabla }\cdot {\boldsymbol u}\left(\frac{1}{2}\rho u^2 + \frac{\gamma _{\rm g} Pg}{\gamma _{\rm g} -1}\right)} \nonumber \\ &&{\quad\quad\quad -\;{\boldsymbol u}\cdot {\nabla }(P_{\rm g} + P_{\rm c}) - {\boldsymbol v}_{\rm A}\cdot \nabla P_{\rm c} - \rho {\boldsymbol u}\cdot {\nabla }\Phi .} \end{eqnarray}

(A7)

CR energy

The time-dependent CR transport equation reads

\begin{equation} \frac{\mathrm{\partial} f}{\mathrm{\partial} t} +({\boldsymbol u}+{\boldsymbol v}_{\rm A})\cdot {\nabla }f= {\nabla }\cdot \left[D{\nabla }f\right] + \left[ {\nabla } \cdot ({\boldsymbol u}+{\boldsymbol v}_{\rm A}) \right] \frac{p}{3} \frac{\mathrm{\partial} f}{\mathrm{\partial} p} \end{equation}

(A8)

where f is the CR distribution function and D the CR diffusion coefficient. We also introduce the CR energy density and pressure

\begin{eqnarray} \epsilon _{\rm c} = \int _0^\infty \text{d}p\; 4\pi p^2 T(p)f(z,p), \end{eqnarray}

(A9)

\begin{eqnarray} P_{\rm c} = \int _0^\infty \text{d}p\; \frac{4\pi }{3} p^2 pvf(z,p), \end{eqnarray}

(A10)

where |$T(p)=\sqrt{p^2c^2 + (mc^2)^2}-mc^2$| is the kinetic energy. Multiplying the transport equation by 4π p²T(p) and integrating in p we get

\begin{eqnarray} &&{\frac{\mathrm{\partial} \epsilon _{\rm c}}{\mathrm{\partial} t} +({\boldsymbol u}+{\boldsymbol v}_{\rm A})\cdot {\nabla }\epsilon _{\rm c} = {\nabla } \cdot \left[\overline{D}{\nabla }\epsilon _{\rm c}\right] }\nonumber \\ &&{\qquad+\; {\nabla } \cdot ({\boldsymbol u}+{\boldsymbol v}_{\rm A}) \int _0^\infty \frac{4}{3}\pi \text{d}p p^3 T(p)\frac{\mathrm{\partial} f}{\mathrm{\partial} p}} \end{eqnarray}

(A11)

where we introduced the average diffusion coefficient

\begin{equation} \overline{D}(z)=\frac{\int _0^\infty dp\; 4\pi p^2 T(p)D(z,p){\nabla }f}{\int _0^\infty dp\; 4\pi p^2 T(p){\nabla }f}. \end{equation}

(A12)

The last integral can be performed by parts, using the fact that dT(p)/dp = v. After rearranging terms and introducing the effective CR adiabatic index ϵ_c = P_c/(γ_c − 1) we obtain:

\begin{eqnarray} \frac{\mathrm{\partial} \epsilon _{\rm c}}{\mathrm{\partial} t} \!+\! {\nabla }\cdot \left[({\boldsymbol u}+{\boldsymbol v}_{\rm A})\frac{\gamma _{\rm c} P_{\rm c}}{\gamma _{\rm c} -1} \right]\!=\! {\nabla }\cdot \left[\overline{D}\frac{{\nabla } P_{\rm c}}{\gamma _{\rm c} -1} \right] \!+\! ({\boldsymbol u}+{\boldsymbol v}_{\rm A})\cdot {\nabla }P_{\rm c}.\nonumber\\ \end{eqnarray}

(A13)

Wave energy

The time evolution of the Alfvén wave spectral energy density |$W({\boldsymbol z},p)$|⁠, taking into account wave advection by the thermal plasma and wave generation by CR streaming instability, is given in Jones (1993)

\begin{equation} \frac{\mathrm{\partial} W}{\mathrm{\partial} t} + {\nabla }\cdot \left[\left(\frac{3}{2}{\boldsymbol u}+{\boldsymbol v}_{\rm A}\right)W\right] =\frac{1}{2}{\boldsymbol u}\cdot {\nabla }W -\Gamma _{\rm CR}W \end{equation}

(A14)

where Γ_CR (see also equation 20) is the wave growth rate due to CR streaming instability (k is the wavenumber and B₀ the average magnetic field)

\begin{equation} \Gamma _{\rm CR}=\frac{16\pi ^2}{3}\frac{{\boldsymbol v}_{\rm A}}{kWB_0^2}\cdot \int _0^\infty \text{d}p p^4 v(p)\delta \left(p- \frac{qB_0}{kc}\right){\nabla }f. \end{equation}

(A15)

We also define the wave energy density and pressure

\begin{eqnarray} \epsilon _{\rm w} = \int _{k_0}^\infty \frac{{B_0}^2}{4\pi } W(z,k)\text{d}k \end{eqnarray}

(A16)

\begin{eqnarray} P_{\rm w}=\frac{\epsilon _{\rm w}}{2}. \end{eqnarray}

(A17)

The hydrodynamic version of the wave equation can be obtained by integrating in |$\int _{k_0}^\infty \frac{{B_0}^2}{4\pi }\text{d}k$|⁠. The integration of the streaming instability term gives

\begin{eqnarray} &&{\int _{k_0}^\infty \frac{{B_0}^2}{4\pi }\Gamma _{\rm CR}W(z,k)\text{d}k} \nonumber \\ &&{\quad=\int _{0}^\infty \int _{k_0}^\infty \frac{4\pi }{3} \frac{{\boldsymbol v}_{\rm A}}{k}\cdot {\nabla }f p^4 v(p) \, \delta \left(p- \frac{qB_0}{kc}\right)\text{d}p\text{d}k} \nonumber \\ &&{\quad= \int _{0}^\infty \int _{k_0}^\infty \frac{4\pi }{3} {\boldsymbol v}_{\rm A}\cdot {\nabla }f p^4 v(p)\frac{kc}{qB_0} \, \delta \left(k- \frac{qB_0}{pc}\right)\text{d}p\text{d}k }\nonumber \\ &&{\quad={\boldsymbol v}_{\rm A}\cdot {\nabla }\int _{0}^\infty \frac{4\pi }{3} f v(p)p^3\text{d}p}\nonumber \\ &&{\quad={\boldsymbol v}_{\rm A}\cdot {\nabla }P_{\rm c}} \end{eqnarray}

(A18)

where we used the properties of the Dirac's delta and the definition of P_c given in equation (A10). Finally, the equation for the wave energy is

\begin{equation} \frac{\mathrm{\partial} \epsilon _{\rm w}}{\mathrm{\partial} t} + {\nabla }\cdot \left[\left(3{\boldsymbol u}+2{\boldsymbol v}_{\rm A}\right)P_{\rm w}\right] ={\boldsymbol u}\cdot {\nabla }P_{\rm w} -{\boldsymbol v}_{\rm A}\cdot {\nabla }P_{\rm c}. \end{equation}

(A19)

In our approach we assume that all generated waves are locally and immediately dumped, so that the hydrodynamic equation for waves is not necessary. On the other hand, the term of wave generation by CR streaming instability, |${\boldsymbol v}_{\rm A}\cdot {\nabla }P_{\rm c}$|⁠, enters the right hand side of the equation for the gas pressure (A5) because the damped waves result in gas heating.

APPENDIX B: STATIONARY HYDRODYNAMIC EQUATIONS IN THE FLUX TUBE

We are dealing with a steady one-dimensional flow, in which the flux-tube geometry is preassigned as shown in Fig. 1. The stationary regime is obtained by setting ∂/∂t = 0 in the hydrodynamic equations. The stationary hydrodynamic equations are listed in Section 2. As for the flux-tube geometry, referring to Fig. 1, if z is the vertical coordinate and A(z) the flux-tube area, the divergence and gradient operators become

\begin{eqnarray} {\nabla }S=\frac{{\rm d}S}{{\rm d}z} \hat{z} \end{eqnarray}

(B1)

\begin{eqnarray} {\nabla }\cdot {\boldsymbol V}=\frac{1}{A(z)}\frac{{\rm d}}{{\rm d}z} \left[ {\boldsymbol A}(z)\cdot {\boldsymbol V} \right]. \end{eqnarray}

(B2)

By applying this prescription for the gradient and divergence operators to the stationary equations we get

\begin{eqnarray} \rho u A = {const}, \end{eqnarray}

(B3)

\begin{eqnarray} AB = {const}, \end{eqnarray}

(B4)

\begin{eqnarray} \rho u \frac{{\rm d}}{{\rm d}z}\left(\frac{u^2}{2} + \Phi \right) = -u\frac{{\rm d}}{{\rm d}z}\left(P_{\rm g} + P_{\rm c}\right), \end{eqnarray}

(B5)

\begin{eqnarray} \frac{{\rm d}P_{\rm g}}{{\rm d}z}=\gamma _{\rm g}\frac{P_{\rm g}}{\rho }\frac{{\rm d}\rho }{{\rm d}z}-(\gamma _{\rm g} -1)\frac{v_{\rm A}}{u}\frac{{\rm d}P_{\rm c}}{{\rm d}z} \end{eqnarray}

(B6)

\begin{eqnarray} \rho uA \frac{{\rm d}}{{\rm d}z}\left(\frac{u^2}{2} +\frac{\gamma _{\rm g}}{\gamma _{\rm g} -1} \frac{P_{\rm g}}{\rho } + \Phi \right) = -A(u+v_{\rm A})\frac{{\rm d}P_{\rm c}}{{\rm d}z}, \end{eqnarray}

(B7)

\begin{eqnarray} \frac{{\rm d}P_{\rm c}}{{\rm d}z}=\frac{\gamma _{\rm c} P_{\rm c}}{\rho }\frac{2u+v_{\rm A}}{2(u+v_{\rm A}-\mathcal {D})}\frac{{\rm d}\rho }{{\rm d}z}. \end{eqnarray}

(B8)

In the derivation of equation (B8) we defined

\begin{equation} \mathcal {D}=\frac{\frac{{\rm d}}{{\rm d}z}\left(A\overline{D}\frac{{\rm d}P_{\rm c}}{{\rm d}z}\right)}{A\frac{{\rm d}P_{\rm c}}{{\rm d}z}} \end{equation}

(B9)

and we used

\begin{eqnarray} \frac{{\rm d}(uA)}{{\rm d}z} = - \frac{uA}{\rho }\frac{{\rm d}\rho }{{\rm d}z} \end{eqnarray}

(B10)

\begin{eqnarray} \frac{{\rm d}(v_{\rm A}A)}{{\rm d}z} = - \frac{v_{\rm A}A}{2\rho }\frac{{\rm d}\rho }{{\rm d}z}. \end{eqnarray}

(B11)

First integrals

It is possible to derive from the hydrodynamic equations two first integrals which come in handy when searching for the wind solution. The first one derives directly from the total energy conservation: summing equations (4) and (5) and moving to the flux-tube geometry, we get

\begin{eqnarray} \frac{u^2}{2} &+&\frac{\gamma _{\rm g}}{\gamma _{\rm g} -1} \frac{P_{\rm g}}{\rho } + \Phi +\frac{\gamma _{\rm c}}{\gamma _{\rm c} -1} \frac{P_{\rm c}}{\rho }\frac{u+v_{\rm A}}{u} \nonumber\\ &&-\frac{\overline{D}\frac{{\rm d}P_{\rm c}}{{\rm d}z}}{(\gamma _{\rm c}-1)\rho u} = {const}. \end{eqnarray}

(B12)

The other first integral is obtained from equation (B8) when |$\overline{D}=0$|⁠, i.e. when |$\mathcal {D}=0$|⁠. Due to this limitation, this first integral is only used in the very first iteration of our calculation, where we assume |$\overline{D}=0$|⁠.

In fact, using equations (B10) and (B11)

\begin{equation} \frac{{\rm d}}{{\rm d}z}\left(A(u+v_{\rm A})\right)=-\frac{A(2u+v_{\rm A})}{2\rho }\frac{{\rm d}\rho }{{\rm d}z}, \end{equation}

(B13)

thus dP_c/dz in equation (B8) becomes

\begin{equation} \frac{{\rm d}P_{\rm c}}{{\rm d}z} = - \frac{\gamma _{\rm c} P_{\rm c}}{A(u+v_{\rm A})} \, \frac{{\rm d}}{{\rm d}z} \left[ A(u+v_{\rm A}) \right] \end{equation}

(B14)

from which we finally get

\begin{equation} P_{\rm c}\left[A(u+v_{\rm A})\right]^{\gamma _{\rm c}} = {const}. \end{equation}

(B15)

The wind equation

Starting from equation (B5) and using equations (B6) and (B8) we get the so-called wind equation, i.e. an equation for the flow speed,

\begin{equation} \frac{{\rm d}u}{{\rm d}z}=u\frac{c_{{\ast }}^2\frac{1}{A} \frac{{\rm d}A}{{\rm d}z} -\frac{{\rm d}\Phi }{{\rm d}z}}{u^2-c_{{\ast }}^2} \end{equation}

(B16)

where we defined the ‘compound sound speed’

\begin{equation} c_{{\ast }}^2=\gamma _{\rm g} \frac{P_{\rm g}}{\rho } + \gamma _{\rm c} \frac{P_{\rm c}}{\rho } \left[ 1-(\gamma _{\rm g}-1)\frac{v_{\rm A}}{u} \right]\frac{2u+v_{\rm A}}{2(u+v_{\rm A} -\mathcal {D})}. \end{equation}

(B17)

Download all slides

Month:	Total Views:
February 2017	3
March 2017	19
April 2017	10
May 2017	12
June 2017	12
July 2017	17
August 2017	4
September 2017	6
October 2017	20
November 2017	10
December 2017	6
January 2018	11
February 2018	4
March 2018	4
April 2018	3
May 2018	5
June 2018	6
July 2018	10
August 2018	4
September 2018	2
October 2018	2
November 2018	9
December 2018	4
January 2019	4
February 2019	4
March 2019	2
April 2019	7
May 2019	3
June 2019	8
July 2019	8
August 2019	7
September 2019	7
October 2019	7
November 2019	8
December 2019	8
January 2020	7
February 2020	12
March 2020	3
April 2020	13
May 2020	3
June 2020	5
July 2020	7
August 2020	3
September 2020	8
October 2020	7
November 2020	11
December 2020	7
January 2021	11
February 2021	14
March 2021	14
April 2021	14
May 2021	7
June 2021	14
July 2021	12
August 2021	8
September 2021	21
October 2021	5
November 2021	11
December 2021	10
January 2022	6
February 2022	6
March 2022	8
April 2022	13
May 2022	6
June 2022	2
July 2022	20
August 2022	16
September 2022	12
October 2022	3
November 2022	2
December 2022	6
January 2023	10
February 2023	6
March 2023	8
April 2023	13
May 2023	7
June 2023	5
July 2023	10
August 2023	9
September 2023	8
October 2023	2
November 2023	5
December 2023	19
January 2024	15
February 2024	3
March 2024	4
April 2024	9
May 2024	12
June 2024	21
July 2024	20
August 2024	9
September 2024	11
October 2024	12
November 2024	8
December 2024	11
January 2025	13
February 2025	3
March 2025	20
April 2025	6
May 2025	8

Article Contents

Cosmic ray driven Galactic winds

Abstract

INTRODUCTION

A SEMI-ANALYTICAL APPROACH TO CR-DRIVEN WINDS

Galactic gravitational potential

IDENTIFICATION OF THE WIND SOLUTION

KINETIC CALCULATION

RESULTS

Reference case

The importance of the near-disc region

CONCLUSIONS

REFERENCES

APPENDIX A: DERIVATION OF THE HYDRODYNAMIC EQUATIONS

Gas mass conservation and equation of motion

Gas internal energy and total energy

CR energy

Wave energy

APPENDIX B: STATIONARY HYDRODYNAMIC EQUATIONS IN THE FLUX TUBE

First integrals

The wind equation

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Cosmic ray driven Galactic winds Free

Abstract

INTRODUCTION

A SEMI-ANALYTICAL APPROACH TO CR-DRIVEN WINDS

Galactic gravitational potential

IDENTIFICATION OF THE WIND SOLUTION

KINETIC CALCULATION

RESULTS

Reference case

The importance of the near-disc region

CONCLUSIONS

REFERENCES

APPENDIX A: DERIVATION OF THE HYDRODYNAMIC EQUATIONS

Gas mass conservation and equation of motion

Gas internal energy and total energy

CR energy

Wave energy

APPENDIX B: STATIONARY HYDRODYNAMIC EQUATIONS IN THE FLUX TUBE

First integrals

The wind equation

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

Cosmic ray driven Galactic winds