Lack of thermal energy in superbubbles: hint of cosmic rays?

TL;DR

This study investigates how cosmic rays influence the thermal and X-ray properties of interstellar superbubbles, revealing that even small CR energy fractions can significantly lower bubble temperatures and potentially explain observed temperature discrepancies.

Contribution

The paper introduces a combined analytic and simulation approach to demonstrate the impact of cosmic rays on superbubble dynamics and thermal profiles, highlighting conditions for efficient CR energy transfer.

Findings

CRs can significantly lower superbubble temperatures by a factor of 2-10.

CR energy transfer is efficient when Mach number exceeds 12 and acceleration time is short.

CR-influenced bubbles exhibit lower X-ray temperatures, addressing the temperature discrepancy problem.

Abstract

Using analytic methods and $1$-D two-fluid simulations, we study the effect of cosmic rays (CRs) on the dynamics of interstellar superbubbles (ISBs) driven by multiple supernovae (SNe)/stellar winds in OB associations. In addition to CR advection and diffusion, our models include thermal conduction and radiative cooling. We find that CR injection at the reverse shock or within a central wind-driving region can affect the thermal profiles of ISBs and hence their X-ray properties. Even if a small fraction ($10-20\%$) of the total mechanical power is injected into CRs, a significant fraction of the ram pressure at the reverse shock can be transferred to CRs. The energy transfer becomes efficient if (1) the reverse shock gas Mach number exceeds a critical value ($M_{\rm th}\gtrsim 12$) and (2) the CR acceleration time scale $\tau_{\rm acc}\sim \kappa_{\rm cr}/v^2$ is shorter than the dynamical time, where $\kappa_{\rm cr}$ is CR diffusion constant and $v$ is the upstream velocity. We show that CR affected bubbles can exhibit a volume averaged hot gas temperature $1-5\times10^{6}$ K, lower by a factor of $2-10$ than without CRs. Thus CRs can potentially solve the long-standing problem of the observed low ISB temperatures.

Tables

Table 1. Table 1: The constant factors presented in equations ( 2.1.1 ) and ( 7 ).

$\gamma$	$s$	$C_1$	$C_2$	$\xi$
5/3	$0.0$	$0.763$	$0.163$	$0.78$
	$1.0$	$0.607$	$0.285$	$0.83$
	$2.0$	$0.376$	$1.000$	$1.00$
13/9	$0.0$	$0.732$	$0.150$	$0.78$
	$1.0$	$0.580$	$0.272$	$0.83$
	$2.0$	$0.357$	$1.000$	$1.00$
4/3	$0.0$	$0.708$	$0.140$	$0.78$
	$1.0$	$0.558$	$0.262$	$0.83$
	$2.0$	$0.341$	$1.000$	$1.00$

Table 2. Table 2: Details of the runs.

$\left[1\right]$	$\left[2\right]$	$\left[3\right]$	$\left[4\right]$	$\left[5\right]$ Micro-physics			$\left[6\right]$ Simulation box details
Model	No. of OB stars	${\rho }_c$	Two-fluid	CR diffusion	Therm Conduction	Cooling	Box size	Grids
	(N)	($m{}_H{}^{}\mathrm{cm}_{}^{-3}$)		(crd)	(t)	(c)	[$r_{\min },r_{\max }$]	[$n$]
$\mathrm{N3}\_\mathrm{d2}$	${10}^3$	${10}^2$	N	N	N	N	[0.1,400.1]	[8000]
$\mathrm{N3}\_\mathrm{d2}\_t$	${10}^3$	${10}^2$	N	N	Y	N	[0.1,300.1]	[6000]
$\mathrm{N3}\_\mathrm{d2}\_c$	${10}^3$	${10}^2$	N	N	N	Y	[0.1,300.1]	[6000]
$\mathrm{N3}\_\mathrm{d2}\_c\_t$	${10}^3$	${10}^2$	N	N	Y	Y	[0.1,300.1]	[6000]
$\mathrm{N3}\_\mathrm{d2}\_\mathrm{cr}$	${10}^3$	${10}^2$	Y	N	N	N	[0.1,400.1]	[8000]
$\mathrm{N3}\_\mathrm{d2}\_\mathrm{crd}$	${10}^3$	${10}^2$	Y	Y	N	N	[0.1,1000.1]	[20000]
$\mathrm{N3}\_\mathrm{d2}\_\mathrm{crd}\_t$	${10}^3$	${10}^2$	Y	Y	Y	N	[0.1,300.1]	[6000]
$\mathrm{N3}\_\mathrm{d2}\_\mathrm{crd}\_c$	${10}^3$	${10}^2$	Y	Y	N	Y	[0.1,1000.1]	[20000]
$\mathrm{N3}\_\mathrm{d2}\_\mathrm{crd}\_c\_t$	${10}^3$	${10}^2$	Y	Y	Y	Y	[0.1,500.1]	[10000]
$\mathrm{N3}\_\mathrm{d0}\_\mathrm{crd}\_c$	${10}^3$	$1$	Y	Y	N	Y	[0.1,1400.1]	[28000]
$\mathrm{N3}\_\mathrm{d3}\_\mathrm{crd}\_c$	${10}^3$	${10}^3$	Y	Y	N	Y	[0.1,1000.1]	[20000]
$\mathrm{N2}\_\mathrm{d2}\_\mathrm{crd}\_c$	${10}^2$	${10}^2$	Y	Y	N	Y	[0.1,1000.1]	[20000]
$\mathrm{N4}\_\mathrm{d2}\_\mathrm{crd}\_c$	${10}^4$	${10}^2$	Y	Y	N	Y	[0.1,1000.1]	[20000]

Table 3. Table 3: The observed ISB parameters.

Equations

$$\frac{d}{dt}(M\dot{R})=4\pi R^{2}P$$

$$\frac{d}{dt}\left(\frac{4}{3}\pi R^3\frac{P}{\gamma -1}\right)=L_{\mathrm{w}}-4\pi R^2\dot{R}P-L_{\mathrm{loss}},$$

$$\rho (r)={\rho }_{\mathrm{c}}{\left(\frac{r_c}{r}\right)}^s,$$

$$R=C_1{L}_{\mathrm{w}}^{1/(5-s)}({\rho }_{\mathrm{c}}{r}_{\mathrm{c}}^s{)}^{-1/(5-s)}{t}_{\mathrm{dyn}}^{3/(5-s)}$$

$$P=C_2{L}_{\mathrm{w}}^{(2-s)/(5-s)}({\rho }_{\mathrm{c}}{r}_{\mathrm{c}}^s{)}^{3/(5-s)}{t}_{\mathrm{dyn}}^{-(4+s)/(5-s)}$$

$$C_1$$

$$C_2$$

$$P=\xi \rho v_{\mathrm{sh}}^2,$$

$$R_{\mathrm{rs}}={\left(\frac{1}{4\pi \xi }\right)}^{1/2}{\rho }^{-1/2}{\dot{M}}^{1/2}{v}_{\mathrm{w}}^{1/2}{v}_{\mathrm{sh}}^{-1},$$

$$R_{\mathrm{rs}}=2.45{\rho }_2^{-3/10}{\dot{M}}_{-4}^{1/2}{L}_{39}^{-1/5}{v}_3^{1/2}{t}_6^{2/5}\mathrm{pc}.$$

$${\rho }_{\mathrm{sw}}=\left(\frac{\gamma +1}{\gamma -1}\right)\xi \rho v_{\mathrm{sh}}^2{v}_{\mathrm{w}}^{-2}$$

$$T_{\mathrm{sw}}=\frac{2(\gamma -1)}{(\gamma +1{)}^2}\frac{\mu m_{\mathrm{H}}}{k_{\mathrm{B}}}{v}_{\mathrm{w}}^2$$

$${\rho }_{\mathrm{sw},\mathrm{tc}}\simeq 0.27L_{39}^{6/35}{\rho }_2^{19/35}{t}_6^{-22/35}{\left[1-R_{\mathrm{rs}}/R_{\mathrm{cd}}\right]}^{2/5}m{}_{\mathrm{H}}\mathrm{c}{\mathrm{m}}^{-3}$$

$$T_{\mathrm{sw},\mathrm{tc}}\simeq 1.26\times 10^7{L}_{39}^{8/35}{\rho }_2^{2/35}{t}_6^{-6/35}{\left[1-R_{\mathrm{rs}}/R_{\mathrm{cd}}\right]}^{-2/5}\mathrm{K}$$

$$L_{\mathrm{loss}}={\int }_0^{R}dr4\pi R^2{n}_{\mathrm{e}}{n}_{\mathrm{i}}{\Lambda }_{\mathrm{N}},$$

$$w=\frac{p_{\mathrm{cr}}}{p_{\mathrm{th}}+p_{\mathrm{cr}}}$$

$${ϵ}_{\mathrm{cr}}=\frac{E_{\mathrm{cr}}}{E_{\mathrm{IN}}},$$

$$\frac{p_{\mathrm{cr}}}{{\dot{M}}^{1/2}{L}_{\mathrm{w}}^{1/2}{R}_{\mathrm{src}}^{-2}}\simeq 0.011{ϵ}_{\mathrm{cr}}{\left(\frac{r}{R_{\mathrm{src}}}\right)}^{-8/3}$$

$$M_{\mathrm{th},\mathrm{rs}}=\frac{v_{\mathrm{w}}}{a_{\mathrm{th},\mathrm{rs}}}\simeq 8.15{\eta }^{-2/15}{R}_{\mathrm{src},\mathrm{pc}}^{-2/3}{\rho }_2^{-1/5}{\dot{M}}_{-4}^{1/6}{L}_{39}^{1/30}{t}_6^{4/15}$$

$${\tau }_{\mathrm{cri}}\simeq 4.65{\eta }^{1/2}{R}_{\mathrm{src},\mathrm{pc}}^{5/2}{\rho }_2^{3/4}{\dot{M}}_{-4}^{-5/8}{L}_{39}^{-1/8}\mathrm{Myr}$$

$$M_{\mathrm{th},\mathrm{fs}}=\frac{v_{\mathrm{sh}}}{a_{\mathrm{th},\mathrm{fs}}}\simeq 10v_{\mathrm{sh},2}{a}_{\mathrm{th},1}^{-1}$$

$$\frac{\partial \rho }{\partial t}+\vec{\nabla }.(\rho \vec{v})=S_{\rho }$$

$$\frac{\partial }{\partial t}(\rho \vec{v})+\vec{\nabla }.(\rho \vec{v}\otimes \vec{v})+\vec{\nabla }(p_{\mathrm{th}}+p_{\mathrm{cr}})=0$$

$$\frac{\partial }{\partial t}\left(e_{\mathrm{th}}+e_{\mathrm{k}}\right)+\vec{\nabla }.\left[\left(e_{\mathrm{th}}+e_{\mathrm{k}}\right)\vec{v}\right]+\vec{\nabla }.\left[\vec{v}(p_{\mathrm{th}}+p_{\mathrm{cr}})\right]$$

$$=p_{\mathrm{cr}}\vec{\nabla }.\vec{v}-\vec{\nabla }.{\vec{F}}_{\mathrm{tc}}+q_{\mathrm{eff}}^{-}+S_{\mathrm{e}}$$

$$\frac{\partial e_{\mathrm{cr}}}{\partial t}+\vec{\nabla }.\left[e_{\mathrm{cr}}\vec{v}\right]=-p_{\mathrm{cr}}\vec{\nabla }.\vec{v}-\vec{\nabla }.{\vec{F}}_{\mathrm{crdiff}}+S_{\mathrm{cr}}$$

$$\Delta t_{\mathrm{CFL}}=(\mathrm{CFL}\mathrm{number})\times \mathrm{MIN}\left(\frac{\Delta x_{\mathrm{i}}}{u_{\max }^{\mathrm{i}}}\right)$$

$$\dot{M}\simeq 5\times 10^{-4}\left(\frac{N_{\mathrm{OB}}}{10^3}\right){\left(\frac{{\tau }_{\mathrm{OB}}}{30\mathrm{Myr}}\right)}^{-1}M{}_{\odot }\mathrm{y}{\mathrm{r}}^{-1}$$

$$L_{\mathrm{w}}\simeq 10^{39}\left(\frac{N_{\mathrm{OB}}}{10^3}\right)\left(\frac{E_{\mathrm{SN}}}{10^{51}\mathrm{erg}}\right){\left(\frac{{\tau }_{\mathrm{OB}}}{30\mathrm{Myr}}\right)}^{-1}\mathrm{erg}{\mathrm{s}}^{-1}$$

$$q_{\mathrm{eff}}^{-}=-n_{\mathrm{i}}{n}_{\mathrm{e}}{\Lambda }_{\mathrm{N}}+\mathrm{Heating}$$

$${\vec{F}}_{\mathrm{tc}}=\chi {\vec{F}}_{\mathrm{classical}}$$

Full text

1 Introduction

Superbubbles driven by stellar winds and supernovae from OB associations are the instruments of stellar feedback that regulate galaxy evolution. These expanding shells form the crucial link between stars and the interstellar medium (ISM) by depositing thermal and kinetic energy, and thereby influencing the star formation process. On a larger scale, they can launch galactic scale outflows if certain conditions are fulfilled (Nath & Shchekinov 2013; Sharma et al. 2014).

The classical model of (1977) provided the basic framework for wind driven bubbles. They described the shock structure expected in ISBs, and worked out the dynamics in the self-similar phase of evolution. Mac Low & McCray (1988) proposed an ISB model by considering correlated SNe and discussed its effects on the galactic scale. All these ideas have been studied in detail with numerical simulations (Keller et al. 2014; Yadav et al. 2017; Vasiliev, Shchekinov & Nath 2017).

In recent years, X-ray observations of various bubbles have led to a closer look into their dynamics ( 2006 and references therein). The presence of a dominant soft X-ray component at $\sim 2\times 10^{6}$ K has also been highlighted in some of these studies (e.g., Chu et al. 2003; 2009). Most of these studies provide the temperature and density of the X-ray emitting plasma which are often used to understand the effective driving force acting on the dense shell (Pellegrini et al. 2011; Lopez et al. 2011; Lopez et al. 2014). However, these analyses show that the best-fit X-ray luminosity and temperature are often lower than that expected from the classical bubble model (Chu, Gruendl & Guerrero 2003; 2009). Even with thermal conduction, which makes the bubble denser and cooler, the discrepancy is not fully resolved. Therefore, we investigate the role of CRs in modifying the ISB properties.

Supernova remnants also show results similar to ISBs. Chevalier (1983) re-examined the blast wave solution to model the effect of relativistic particles (see also Vink et al. 2010; Bell 2014). He showed that the injection of CRs can reduce the thermal energy inside the blast wave. This idea has been confirmed by analyzing the post shock temperature in RCW $86$ supernova remnant (Helder et al. 2009).

Recently, It has been suggested that ISBs are the preferred sites of CR acceleration instead of isolated supernova remnants (SNRs). Since massive stars usually form in clusters, most of the power injected by supernovae into the ISM is mediated through superbubbles and not isolated SNRs (Higdon & Lingenfelter 2005). Binns et al. (2005) have suggested that the isotopic anomalies in the composition of Galactic cosmic rays, in particular, the enhanced ${}^{22}Ne/^{20}Ne$ ratio, are suggestive of CRs being accelerated out of the matter inside ISBs. Parizot & Drury (1999) argued that a superbubble model for the origin of Galactic cosmic rays can explain the evolution of light elements Li, Be and B. Gamma-ray observations of the Cygnus superbubble have also shown that CRs are accelerated in ISBs (Ackermann et al. 2011). Recently Eichler (2017) has suggested that cosmic ray grammage traversed is correlated with the properties of the source, meaning that the escape occurs near the production site, suggestive of ISBs. Observations of extragalactic superbubbles in IC 10 (Heesen et al. 2014) and the Large Magellanic Cloud (Butt & Bykov 2008) have also suggested the ISB origin of CRs. Therefore the dynamics of ISBs, and the possible effects of CR on them, deserve to be studied in detail.

The role of CR feedback in galaxy evolution has also been discussed by several authors (e.g., Salem & Bryan 2013; Booth et al. 2013; Wiener, Pfrommer & Oh 2017). Although these studies included various physical processes e.g. self-gravity, cooling, star formation, the effects of the individual processes are difficult to disentangle.

In this paper, we present a two-fluid model of ISB. We start with a standard ISB model and include thermal conduction, radiative cooling and CR diffusion one by one to understand the role of each processes. We show that in the absence of CRs the thermal pressure of hot gas depends on the ambient density (almost independent of density profile) and the shell speed. Then we explore the effect of CRs. We have found that CRs can affect ISB via shock interactions and their effects mainly depend on the shock Mach number (Drury & Volk 1981; Becker & Kazanas 2001) and CR diffusion coefficient.

The contents of this paper are organized as follows. We present a broad analytical frame work in section 2. The details of the simulation set-up are given in section 3. The results from various runs are presented in section 4. In section 5 we discuss the dependence of the results on various parameters. In section 6 we comment on the astrophysical implication of this study. Finally, we conclude in section 7 by highlighting the main results of this paper.

2 Analytical prelude

Consider the powerful wind from OB association that drives an ISB. The dynamics and the structure of an ISB is usually understood by the momentum and energy conservation equations:

[TABLE]

[TABLE]

where $P$ is the pressure inside the bubble (assumed to be $\gg$ ambient pressure), $R$ is the position of the swept-up ISM (hereafter ‘shell’) w.r.t. the central source111By the source region, we mean a region within which most of the stars are located, i.e., the region which is driving the ISB., $M=\int_{0}^{R}dr\,4\,\pi\,r^{2}\,\rho(r)$ is the swept-up ambient mass, $\rho(r)$ is the ambient density (c.f. equation (3)), $L_{\rm w}$ is the wind power, and $L_{\rm loss}$ is the loss of energy due to radiative cooling.

In order to find a general solution of equations (1) - (2), we choose the ambient density profile as

[TABLE]

where the choice of the parameter ‘$s$’ determines the ambient density profile. In next sections, we discuss the solutions in different cases.

2.1 One-fluid standard ISB

2.1.1 Adiabatic evolution

Consider the case when the dynamical time ($t_{\rm dyn}$) of an ISB is shorter than the cooling time-scales of its different regions (for details see 1975, sections $3$ and $5.2$ in Gupta et al. 2016). Under this assumption, the contribution of $L_{\rm loss}$ in the r.h.s of equation (2) is negligible. Assuming that at any given instant, $R\propto t^{\alpha}_{\rm dyn}$ (where $\alpha>0$) and substituting $P$ from equation (1) to equation (2), we obtain

[TABLE]

[TABLE]

where

[TABLE]

Note that, the choice $s=0$ corresponds to an ISB expanding in a uniform medium ( 1977).

The pressure $P$ (the shocked wind (SW) pressure) given in equation (5) can be re-written in following form,

[TABLE]

where $v_{\rm sh}$ is the shell velocity ($=dR/dt$; see equation (4)), $\rho$ is the unshocked ISM density (i.e., the upstream density; equation (3)) and $\xi=C_{\rm 1}^{(s-2)}\,C_{\rm 2}\left[(5-s)/3\right]^{2}$. The values of $C_{\rm 1}$, $C_{\rm 2}$ and $\xi$ for different adiabatic constants ($\gamma$) and $s$ are given in Table 1. Note that, $\xi$ depends weakly on the density profile ($s$), and therefore, equation (7) is a robust estimate of the interior pressure of an ISB.

Using equation (7), the position of the reverse shock can be readily obtained by equating the ram pressure and the gas pressure of the SW region. This gives,

[TABLE]

where $\dot{M}$ is the mass loss rate of the source, and $v_{\rm w}\approx(2L_{\rm w}/\dot{M})^{1/2}$ (Chevalier & Clegg 1985, hereafter CC85 ) is the wind velocity. For a uniform ambient medium, the position of the reverse shock is given by,

[TABLE]

Here $\rho_{\rm 2}=\rho/(10^{2}\,m{\rm{}_{H}\,cm^{-3}})$, $\dot{M}_{\rm-4}=\dot{M}/(10^{-4}\,\rm{M_{\odot}\,yr^{-1}})$, $L_{\rm 39}=L_{\rm w}/(10^{39}\,\rm{erg\,s^{-1}})$, $v_{\rm 3}=v_{\rm w}/(10^{3}\,\rm{km\,s^{-1}})$ and $t_{\rm 6}=t_{\rm dyn}/(1\,\rm{Myr})$. Neglecting thermal conduction, the density of the shocked wind region ($\rho_{\rm sw}$) can be estimated by applying the shock jump condition at the reverse shock. This gives

[TABLE]

[TABLE]

If one uses a typical wind velocity $v_{\rm w}\sim 2000\ {\rm km\,s^{-1}}$ (Leitherer et al. 1999), then $T_{\rm sw}\approx 5\times 10^{7}$ K, which is higher than the observed ISB temperatures. At these temperatures and densities, the electron and ion temperatures are equalized by Coulomb collisions.

2.1.2 Thermal conduction

Thermal conduction transfers heat from the hot bubble interior to the cooled shell, which causes the evaporation of the swept-up shell mass into the interior of the bubble (Cowie & McKee 1977). The detailed self-similar analysis of (1977) showed that for a uniform ambient medium, thermal conduction decreases the temperature of the hot gas and enhances the density of the interior. Choosing the classical isotropic thermal conductivity $\kappa_{\rm th}=6\times 10^{-7}T^{5/2}$ cgs (Mac Low & McCray 1988), we have,

[TABLE]

[TABLE]

The use of classical thermal conduction, however, can be questioned because the assumption of isotropy is not valid in the presence of magnetic field. Furthermore, at the early times, the mean free path of electrons ($\lambda_{\rm m}$) $\gtrsim$ the temperature gradient scale ($l_{\rm T}$), which can cause conduction to be saturated. Therefore, equations (12) and (13) represent the upper limiting case of thermal conduction.

2.1.3 Radiative cooling

For a radiative bubble, cooling can delay the formation of the reverse shock (see section 5.2 in Gupta et al. 2016). However, once it forms, the qualitative description remains the same as above but their quantitative results change because of the term $L_{\rm loss}$ (see equation (2)). The cooling loss rate ($L_{\rm loss}$) is defined as

[TABLE]

where $n_{\rm i}/n_{\rm e}$ is the electron/ion number density, and $\Lambda_{\rm N}$ is the normalized (w.r.t. ($n_{\rm i}n_{\rm e})$) cooling rate.

Note that, although the contact discontinuity (hereafter ‘CD’) occupies a much smaller volume compared to the size of the bubble, the contribution of the term $L_{\rm loss}$ mainly comes from the CD, since this is the region where the temperature passes though the peak of the cooling curve. It can be shown that, if $t_{\rm dyn}$ is longer than the shell cooling time scale ($\tau_{\rm shell}$, see equation ($4$) in Gupta et al. 2016) then the term $(L_{\rm w}-L_{\rm loss})$ in equation (2) can be approximated as $\eta L_{\rm w}$ where $\eta$ ($\leq 1$) can be thought of as an energy efficiency parameter. This parameter ($\eta$) depends mainly on the ambient density (see section 6.3 in Gupta et al. 2016). The expressions in equations (4) and (5) are accordingly changed by replacing $L_{\rm w}$ by $\eta L_{\rm w}$.

2.2 Two-fluid ISBs

The discussions so far have been confined to the case of one-fluid ISB. Here we discuss the modifications due to CRs.

The effects of CRs depend on three main parameters : (1) the regions where CRs are accelerated/injected, (2) the fraction of energy that goes into CRs and (3) the CR diffusion coefficient. Note that all these parameters are uncertain and therefore we cover a range of models and parameters.

For CR affected ISBs, the interior gas can be considered as a mixture of thermal and non-thermal particles. It is useful to take an effective adiabatic index $\gamma_{\rm cr}(=4/3)<\gamma<\gamma_{\rm th}(=5/3)$ to infer the dynamics. The adiabatic index of the gas mixture depends on the fraction of the thermal and non-thermal (CR) particles in the gas, which may vary between different regions in the ISB. Moreover, the description becomes more complicated when CR diffusion is considered. For a preliminary understanding, we divide the discussion based on the CR injection region.

2.2.1 Injection at the shocks

Consider a scenario in which the CRs are accelerated at the forward shock (hereafter FS; equivalent to CR injection at FS). In this case, we do not expect to see any change in the interior because CRs do not penetrate CD. However, if CR acceleration happens at the reverse shock (hereafter RS) then CRs can diminish the thermal pressure in the SW region. It results in a decrease of SW temperature which reduces the effect of thermal conduction.

For injecting CR at the shock, a parameter is commonly used (Chevalier 1983; Bell 2014), which is defined as

[TABLE]

where $p_{\rm cr}$ and $p_{\rm th}$ are the CR and thermal pressure. The implications for CR injection at shocks (with and without CR diffusion) are discussed in detail in section 4.2.1.

2.2.2 Injection at the source region

CR injection in the source region presents a different scenario. In this case, it is assumed that some fraction of total deposited energy by stellar winds and/or SNe goes into CRs (Salem & Bryan 2013). Here, the injection parameter is defined as

[TABLE]

where $E_{\rm cr}$ is the energy deposited in CRs and $E_{\rm IN}$ is the total deposited energy.

The main difference between equations (15) and (16) is that the former parameter is defined at a shock whereas the latter one is defined at the source. Unlike the previous case, CR injection in the source region causes a free wind CR pressure profile ($p_{\rm cr}$) which makes the bubble structure different from the standard one-fluid ISB.

To obtain CR profile in the free wind region, one has to solve a two-fluid steady state model along the lines of CC85 . Inside the free wind, we expect that $p_{\rm cr}\propto\rho^{4/3}$ (adiabatic expansion of the wind) i.e. $p_{\rm cr}\propto r^{-8/3}$. The normalization constant depends on the fraction of energy injected as CRs and also on the input source parameters (see Table 1 in CC85 ). Our numerical simulation shows that the power law profile, which is valid only in the absence of CR diffusion, and can be written as

[TABLE]

where $R_{\rm src}$ is the radius of the source region.

As the free wind reaches the reverse shock, one does not have a one-fluid shock. Two-fluid shocks have been previously studied by several authors (Drury & Volk 1981; Drury & Falle 1986; Wagner, Falle & Hartquist 2007). CR diffusion plays a key role in two-fluid shock. For some upstream parameters, CRs can diffuse across the shock and the downstream CR pressure can increase. The downstream CRs then diffuse further, which changes the upstream CR pressure. As a result, the CR particles cross the shock multiple times (which is also known as diffusive shock acceleration, Becker & Kazanas 2001) and CRs can modify the shock.

Drury & Volk (1981) showed that, for a two-fluid shock, parameters of the downstream fluid can have three possible solutions: (1) Globally smooth, (2) Discontinuous solution and (3) Gas mediated sub-shock. Becker & Kazanas (2001) classified the solution parameter space based on the upstream gas and CR Mach numbers (see their equation ($25$)), which are defined as $M_{\rm th}=v/a_{\rm th}$ and $M_{\rm cr}=v/a_{\rm cr}$ respectively where $a_{\rm th}=\sqrt{\gamma_{\rm th}p_{\rm th}/\rho}$ and $a_{\rm cr}=\sqrt{\gamma_{\rm cr}p_{\rm cr}/\rho}$ are the adiabatic sound speed for the respective fluids. They showed that if the gas Mach number exceeds a critical value $12.3$, then the downstream (post shock) CR pressure dominates over thermal pressure and the shock structure is globally smooth (see their Figures $8$ and $16$). We, therefore, estimate the reverse shock and forward shock Mach numbers to predict the shock structure.

In case of the reverse shock (RS), the upstream wind velocity can be assumed to be the same as the wind velocity ($v_{\rm w}$) because the RS velocity $\ll v_{\rm w}$ (see equation (8)). The upstream sound speed is estimated using the steady state free wind pressure and density profiles (CC85 ), which gives $a_{\rm th,rs}=0.56(R_{\rm rs}/R_{\rm src})^{-2/3}\dot{M}^{-1/2}\,L_{\rm w}^{1/2}$. For the uniform ambient, this yields

[TABLE]

The most interesting thing about the RS is that the Mach number increases with time. Therefore, the RS is expected to show a globally smooth profile after a timescale $\tau_{\rm cri}$. For a uniform ambient medium, we obtain

[TABLE]

Note that, in equation (18) we have included the energy efficiency term $\eta$ (discussed in section 2.1.3). This shows that a radiative ISB can satisfy the critical Mach number criterion ( $M_{\rm th,rs}>12.3$) at an earlier time than the adiabatic ISB.

The forward shock (FS) Mach number can easily be found by assuming the ambient medium is at rest. This gives

[TABLE]

It is worth mentioning that depending on the ambient temperature, $M_{\rm th,fs}$ can exceed $12.3$. Therefore, the FS can also show a globally smooth shock structure. However, note that the Mach number of the forward shock decreases with time.

The important point is that, for a globally smooth solution, the upstream kinetic energy is mostly transferred to the downstream CR pressure. This diminishes the thermal pressure and can change the density and temperature profiles of ISBs.

3 Simulation set-up

In order to study the detailed effects of CRs in ISBs, we have developed a 1-D, two-fluid code that we call TFH (standing for Two Fluid Hydrodynamics) code. TFH solves the following set of equations :

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $\rho$ is the mass density, $\vec{v}$ is the fluid velocity, $p_{\rm th}/p_{\rm cr}$ is the thermal/CR pressure, $e_{\rm k}=\rho\,v^{2}/2$ is the kinetic energy density, $e_{\rm th}=p_{\rm th}/(\gamma_{\rm th}-1)$ and $e_{\rm cr}=p_{\rm cr}/(\gamma_{\rm cr}-1)$ are the thermal (gas) and CR energy densities respectively (Drury & Volk 1981; Wagner, Falle & Hartquist 2007; Pfrommer et al. 2017). The terms $S_{\rm\rho}$, $S_{\rm th}$ and $S_{\rm cr}$ in equations (21), (3) and (24) represent mass and energy terms deposited by the driving source. $F_{\rm tc}$ and $F_{\rm crdiff}$ represent thermal conduction flux and CR diffusion flux respectively. The term $q^{-}_{\rm eff}$ accounts for the radiative energy loss of the thermal gas. For simplicity, we exclude the gas heating due to CR streaming (Guo & Oh 2008).

Currently TFH has two solvers: (1) a ZEUS-like solver ( 1992) which performs transport and source terms separately (hereafter, TS; for details, see chapter $5$ in Dullemond 2009222See http://www.mpia.de/homes/dullemon/lectures/fluiddynamics) and (2) HLL solver ( 1994; for details, see chapter $10$ in 2009). Both solvers use the finite volume method and are first order accurate. Because of the first order scheme, all runs are done with very high resolution (typical resolution is $0.05$ pc). Since TFH performs two-fluid simulation, we have defined an effective sound speed as $c_{\rm s,eff}=\sqrt{(\gamma_{\rm cr}p_{\rm cr}+\gamma_{\rm th}p_{\rm th})/\rho}$ to estimate the CFL time step (Courant, Friedrichs & Lewy 1928). The CFL time step is defined as

[TABLE]

where $\Delta x_{\rm i}$ is the separation between $i^{\rm th}$ and $(i+1)^{\rm th}$ grids, $u^{\rm i}_{\rm max}={\rm MAX}(|S_{\rm L}|,|S_{\rm R}|)$, $S_{\rm L}={\rm MIN}[(v^{\rm i-1}-c_{\rm s,eff}^{\rm i-1}),(v^{\rm i}-c_{\rm s,eff}^{\rm i})]$ and $S_{\rm R}={\rm MAX}[(v^{\rm i-1}+c_{\rm s,eff}^{\rm i-1}),(v^{\rm i}+c_{\rm s,eff}^{\rm i})]$ is the maximum wave speed between the left and right moving waves from the interface. TFH has gone through various test problems, and the comparisons with the publicly available code PLUTO ( 2007) are shown in Appendix A.

For the runs in this paper, we set the solver to TS (see Appendix B) and the CFL number to $0.2$. The details of important runs are given in Table 2. In the following sections, we discuss the terms on the RHS of equations (21) - (24).

3.1 Ambient medium

We consider a uniform ambient density and temperature. The ambient temperature is chosen such that the thermal pressure is $1$ eV cm*-3*. For the runs with CRs, we set the CR pressure to $1$ eV cm*-3* (i.e., assuming equipartition between CRs and thermal gas). For all cases, the initial velocity of gas is set to zero. We assume the metallicity of the injected materials to be solar, $Z=\,Z_{\rm\odot}$, the same as in the ambient gas, where $Z_{\rm\odot}$ is the solar metallicity.

3.2 Mass and energy deposition

The source terms $S_{\rho}$, $S_{e}/S_{\rm cr}$ represent the deposited mass and thermal/CR energy per unit time per unit volume. We add mass and energy in a region of radius $R_{\rm src}=1$ pc. The radius ($R_{\rm src}$) is chosen such that the energy loss rate is much less than the energy injection rate. This condition ensures that cooling at the initial stage will not suppress the production of strong shocks (for details, see section 4.2.1 in Sharma et al. 2014).

The mass loss rate ($\dot{M}$) and the mechanical (wind) power of the source ($L_{\rm w}$) are related to the star formation rate (SFR) in the cloud and life-time of the individual stars. If $M_{\rm st}$ is the stellar mass in the cloud and $N_{\rm OB}$ is the number of OB stars with the mean-sequence lifetime $\tau_{\rm OB}$, then the ${\rm SFR}$ can be assumed to be $M_{\rm st}/\tau_{\rm OB}$ (Mac Low & McCray 1988). Assuming $M_{\rm st}\approx 50\,N_{\rm OB}\,M{\rm{}_{\odot}}$ and $\tau_{\rm OB}\approx 30\,{\rm Myr}$, we get ${\rm SFR}\simeq 1.67\times 10^{-6}\,N_{\rm OB}\,M{\rm{}_{\odot}\,yr^{-1}}$. We have also assumed that $\dot{M}=0.3\,{\rm SFR}$ and $10^{51}$ erg is the energy released in each supernova which yield

[TABLE]

[TABLE]

If one assumes that the stellar winds and/or SNe deposit energy and mass continuously then the equations (26) and (27) can be considered in the terms of SNe. For such case, the time interval between consecutive SNe is equivalent to $\delta t_{\rm SN}\simeq\tau_{\rm OB}/N_{\rm OB}=0.03\,(N_{\rm OB}/10^{3})^{-1}$ Myr, and the total energy and mass deposited during that duration are $10^{51}$ erg and $15$ $M{\rm{}_{\odot}}$ respectively. Our choice gives the wind velocity $v\simeq\sqrt{(}2L_{\rm w}/\dot{M})\approx 2500$ km s*-1* which is consistent with Starburst99 (Leitherer et al. 1999) for a constant SFR or for a coeval star cluster.

To inject cosmic rays in the source region (c.f., section 4.2), we use a injection parameter $\epsilon_{\rm cr}$ (see equation (16)) to specify the fraction of the total input energy that given to CR (Salem & Bryan 2013; Booth et al. 2013). Therefore, $S_{\rm\rho}=\dot{M}/V_{\rm src}$, $S_{\rm th}=(1-\epsilon_{\rm cr})L_{\rm w}/V_{\rm src}$ and $S_{\rm cr}=\epsilon_{\rm cr}L_{\rm w}/V_{\rm src}$. For fiducial runs, whenever CR injection at the source is mentioned, we set $\epsilon_{\rm cr}=0.2$ otherwise it is set to zero. We explore the dependence of the results on various parameters in section 5.

3.3 Cooling & Heating

TFH uses the operator splitting method to include radiative cooling and heating. We use standard definition of cooling function as given below

[TABLE]

where $\Lambda_{\rm N}$ is the normalized cooling function (CLOUDY, 1998) and $n_{\rm i}/n_{\rm e}$ are the electron/ion number density. The region where the cooling energy loss becomes comparable to its thermal energy, TFH subcycles cooling. The number of sub-steps depends on the stiffness of cooling. We artificially stop cooling when the gas temperature goes below $10^{4}$ K. This corresponds to the shell temperature which is maintained due to the photo-heating by the radiation field of the driving source (for details see Figure 4 in Gupta et al. 2016).

3.4 Diffusion terms

Thermal conduction and CR diffusion are also treated using operator splitting. Note that, for stability, the diffusion terms can have a much smaller time step compared to the CFL time step. TFH handles this by performing sub-cycling at each CFL time step.

3.4.1 Thermal conduction

The thermal conduction flux $F_{\rm tc}$ is defined as:

[TABLE]

where $\vec{F}_{\rm classical}=-\kappa_{\rm th}\,\vec{\nabla}T$ and $\kappa_{\rm th}=C\,T^{5/2}$ is the coefficient of thermal conduction ( 1962) and $C$ is chosen to be $6\times 10^{-7}$ in cgs unit. The factor $\chi$ in equation (31) limits the conduction flux when it approaches $\vec{F}_{\rm staturated}\approx(2\,k_{\rm B}\,T/\pi m_{\rm e})^{1/2}n_{\rm e}k_{\rm B}T_{\rm e}$ (see Cowie & McKee 1977) which is defined as $\chi=F_{\rm staturated}/(F_{\rm classical}+F_{\rm staturated})$.

The thermal conduction time step ($\Delta t_{\rm tc}$) is chosen to be

[TABLE]

where $\kappa_{\rm th,avg}=(\kappa_{\rm th,i-1}+\kappa_{\rm th,i})/2$ is the cell averaged thermal conductivity.

3.4.2 CR diffusion

CR diffusion follows a similar method as thermal conduction (section 3.4.1) except that the cosmic ray diffusion flux $F_{\rm crdiff}$ is defined as:

[TABLE]

where $\kappa_{\rm cr}$ is the cosmic ray diffusion constant. Here we are only concerned with the hydrodynamical effects of CRs, but one should remember that $\kappa_{\rm cr}$ is an integral over the CR energy distribution function (equation (7) in Drury & Volk 1981). We consider $\kappa_{\rm cr}=5\times 10^{26}\,{\rm cm^{2}\,s^{-1}}$ as a fiducial value which is consistent with the recent findings in the star forming/SNe region (e.g., Ormes, Ozel & Morris 1988; Gabici et al. 2010; Li & Chen 2010; Giuliani et al. 2010) but smaller than the value usually adopted in the global ISM. We discuss the dependence of simulation results on this choice in section 5.

4 Results

4.1 Ideal one-fluid ISB

To begin with, we recall the structure of a standard ISB ( 1977). We run an ideal one-fluid model by turning-off all microphysics (Table 2). Figure 1 displays the density profile at $t_{\rm dyn}=2.5$ Myr from this run.

Starting from the left, Figure 1 shows that the ISB consists of four distinct regions - (1) the source region where mass and energy are deposited, (2) the free-wind region where the wind expands adiabatically, (3) the shocked-wind region and (4) the swept-up ambient medium shell. Between the regions (3) and (4), there is a contact discontinuity (CD) which separates the ambient material from the ejecta material.

In following sections, we include various microphysical processes one by one, for a better understanding of each process separately.

4.2 CR injection in different regions

4.2.1 At the forward/reverse shock

In supernova remnant studies, the downstream CR pressure fraction is taken as $p_{\rm cr}=w\,p_{\rm tot}$ (see equation (15)) where the physical origin of $w$ is poorly understood (Bell 2014). For the purpose of illustration, we choose $w=1/2$. For implementing CR pressure fraction, we have written a shock detection module which identifies the shock location and fixes the CR pressure accordingly. We have tested this module by comparing it with the blast wave self-similar solution of Chevalier (1983) (see Appendix A.2).

The results with CR injection at forward (top panel) and reverse (middle panel) shocks are displayed in Figure 2. Black line stands for a standard one-fluid ISB (section 4.1). The red curve represents the case where only CR advection is taken into account. For the blue curve, in addition to advection, CR diffusion has been turned on.

With injection at forward shock, in the absence of CR diffusion (red curve), a fraction of the post-shock thermal pressure appears as CR pressure. In this case, the contact discontinuity does not allow CR to be advected into the bubble i.e. the total pressure in the shell remains the same. Therefore, it changes the effective adiabatic index of the swept-up shell ($\gamma_{\rm eff}=[5+3\,w]/[3(1+w)]$), resulting in a higher density jump ($\simeq 5.5$) compared to the one-fluid case (for which the compression ratio is $4$). A more tangible difference appears when CR diffusion is turned on (blue curves). In this case, CRs enter inside the bubble and also diffuse out of the shell. Therefore, the total energy in the shell is now reduced. This results in a much higher density jump, similar to the case with radiative cooling.

In the middle panel of Figure 2, CRs have been injected at the reverse shock (RS). In the absence of CR diffusion (red curves), the CR pressure in the SW region remains almost the same as at the RS. As a result, the bubble size is slightly smaller than the standard bubble because of $\gamma_{\rm eff}=13/9$. In this case, one can use an effective $\gamma$ to determine the size of the bubble (see Table 1). However, this conclusion is not valid if one considers CR diffusion. With CR diffusion the density jump at the RS is not sharp and the size of the bubble is smaller (compare the blue and black curves).

The conclusions from this section are: (1) the CR diffusion plays an important role and (2) injection of CRs at the reverse shock can change the ISB structure.

4.2.2 At the driving source region

The numerical set-up for CR injection in driving source region is discussed in section 3.2 and the results are shown in bottom panel of Figure 2.

In absence of CR diffusion (red curves in panel (c) of Figure 2), the free wind CR pressure profile follows $p_{\rm cr}\propto r^{-8/3}$ (see equation (17)). The grey line in the bottom right most panel of Figure 2 displays this relation. The injected CR particles are advected up to the CD and the structure does not show any significant difference from that of the single fluid case.

When diffusion is turned on (blue curves in panel (c)), we find that the CR pressure at the source region decreases which is expected because CRs diffuse from the high to low pressure region. The striking feature is that, at the RS, the CR pressure is quite large compared to the case of no diffusion. This jump in CR pressure is a property of two-fluid shocks as discussed in section 2.2.2.

Figure 3 displays the time ($t_{\rm dyn}$) evolution of the reverse shock Mach number for three different values of $N_{\rm OB}$, where the blue colour stands for our fiducial choice. The solid lines in this figure display the analytical estimate of the RS gas Mach number (equation (18)) and the blue circles represent the numerical results obtained from our fiducial run. This model shows that the RS satisfies the globally smooth condition after $t_{\rm dyn}\approx 1.7$ Myr. The Mach number analysis for the FS shows that at early times, the FS also satisfies $M_{\rm th}>12.3$ (not displayed in this figure). However in our simulations, we do not see globally smooth shock at FS. The reason is discussed as follows.

For a large post shock CR pressure i.e. for efficient CR acceleration, the CR particles should cross the shock multiple times. Becker & Kazanas (2001) found the critical Mach number for CR dominated post shock in steady state. Steady state assumes that CR particles have sufficient time to cross the shock multiple times. For time dependent calculation, we must consider CR acceleration time scale ($\tau_{\rm acc}$). Drury (1983) discussed that if $\Delta t$ is the average time taken by the upstream CR particles to cross the shock and to return back from downstream (thereby, CRs complete one cycle), and $\Delta\mathcal{P}$ is the average momentum gain in one complete cycle, then

[TABLE]

where subscripts $1,2$ stands for the upstream and down stream flow respectively (also see Blasi, Amato & Caprioli 2007). Therefore, we can see the globally smooth solution only if $\tau_{\rm acc}\ll\tau_{\rm cri}<t_{\rm dyn}$.

Equation (32) can be written as $\tau_{\rm acc}\sim 6\,\kappa_{\rm cr}/v_{\rm 1}^{2}$ ($\equiv$ CR diffusion time scale). For the FS and RS, the condition $t_{\rm dyn}\gg\tau_{\rm acc}$ yields

[TABLE]

[TABLE]

where we have taken $v_{\rm 1}$ as the FS velocity in a uniform ambient medium ($\approx 0.6\,L_{\rm w}^{1/5}\,\rho^{-1/5}\,t_{\rm dyn}^{-2/5}$) and the wind velocity ($v_{\rm w}$) respectively. To illustrate this, we have estimated $\tau_{\rm acc}$ for the FS and RS (Model:${\rm N3\_d2}$; Table 2), and the results are shown in Figure 4.

Figure 4 plots $\tau_{\rm acc}$ as a function of dynamical time for FS and RS where CR diffusion constant is taken as $\kappa_{\rm cr}=5\times 10^{26}$ cm2 s*-1*. The dotted and dashed-dotted lines represent the analytic estimate of acceleration time (i.e., equation (32)), and superposed on that, the blue solid and empty circles represent the simulation result for the respective shocks. This figure shows that $\tau_{\rm acc}$ is longer than $t_{\rm dyn}$ for the FS. This is so because the velocity of the upstream and down stream flow is very small (at least a factor of ten compared to RS) and decreases with time. Therefore a smooth solution is not achieved although it satisfies the upstream condition of the steady state model by Becker & Kazanas (2001). We conclude that for globally smooth solution, the condition is $\tau_{\rm acc}\ll\tau_{\rm cri}<t_{\rm dyn}$. This criterion limits the choice of $\kappa_{\rm cr}$ which is discussed in section 5.1. Our two-fluid simulation results are consistent with this conclusion and the results are shown in Figure 5.

Figure 5 shows the snapshots of density and thermal/CR pressure profiles at two different dynamical times. During the early evolution, the RS Mach number is less then $12.3$ (see Figure 3) and therefore the upstream CR does not dominate over the thermal pressure. As time evolves, after $\tau_{\rm cri}$, the gas Mach number exceeds the critical value and the shock becomes dominated by the CRs. This CR dominated shock is representative of globally smooth solution as first predicted by Drury & Volk (1981). In this case, the maximum CR pressure ($p_{\rm cr,max}$) in the SW region depends on the input source parameters, the ambient density and the dynamical time. We have found that, $p_{\rm cr,max}$ does not exceed the thermal pressure for the one-fluid ISB (i.e., $p_{\rm cr,max}\leq P$; see equation (5)), consistent with the total energy conservation.

The results discussed in this section do not include radiative cooling and thermal conduction i.e., all the changes are only due to CRs. We will discuss more realistic cases below.

4.3 Toward a realistic model

In the previous section, we have seen that the energy exchange between thermal and non-thermal particles becomes significant at the reverse shock. We have also noticed that when $p_{\rm cr}\gtrsim p_{\rm th}$, the cosmic ray fluid starts affecting the inner structures (see Figure 5). For a radiative bubble, energy loss from the dense shell reduces its thermal energy. Therefore, we expect to see the impact of CRs at an earlier time than that in the adiabatic case. In order to get a more realistic picture, we discuss the effect of thermal conduction and then we turn on radiative cooling.

4.3.1 Effect of thermal conduction

To discuss the role of the thermal conduction in CR affected ISBs, we present a comparison of one-fluid and two-fluid models. Figure 6 displays the density (blue curves) and temperature (red curves) profiles at $t_{\rm dyn}=2$ Myr. For a one-fluid ISB without thermal conduction (solid lines), the bubble temperature is high ($\sim 10^{8}$ K). With thermal conduction (dashed lines), the temperature drops to $\sim 10^{7}$ K, and the SW density increases. For a two-fluid ISB, the temperature is noticeably smaller than a one-fluid ISB even without thermal conduction (dotted lines). With thermal conduction it does not show a significant difference, except that, it smoothens the temperature near the CD (dash-dotted lines). Therefore, we conclude that CRs reduce the effect of thermal conduction. We have reached the same conclusion for ISBs with radiative cooling. In the following sections, we continue our discussion without thermal conduction.

4.3.2 Volume averaged quantities

Figure 7 displays the volume averaged hot gas ($>10^{5}$ K) pressure and the CR pressure as functions of time for various models. The solid and empty circles display the thermal and CR pressure respectively. The black solid lines (top curve : adiabatic and bottom curve : radiative) display equation (7) ($\gamma=5/3$, $s=0$) where the shell velocity is estimated from the simulation. In this case, the analytical result agrees with the simulation (as shown by the concurrence of red points with the top black line, and green points with the bottom black line). For two-fluid ISBs (yellow and blue curves), as time evolves, the thermal pressure deviates from this relation. This deviation is mainly because in a two-fluid ISB, a major fraction of the free wind kinetic energy goes to CR.

4.4 Evolution of different energy components

Now we come to the evolution of the kinetic/thermal/CR energy and the radiative loss for four different models. For all cases, we estimate the change of the entity ${\rm x}$ in the simulation box and normalized it w.r.t to the total deposited energy by that epoch. These entities have been estimated using

[TABLE]

[TABLE]

where ${\rm x}$ refers to the kinetic/thermal/CRs energy, $E_{\rm loss}$ is the radiative loss and $E_{\rm IN}^{\rm t}$ is the total deposited energy at the source. Our results are displayed in Figure 8.

The top panel of this figure displays the dynamical evolution of the kinetic energy (KE, red curve) and thermal energy (TE, green curve) for a one-fluid adiabatic ISB. The cooling is turned on in the second panel. In this case, the magenta curve represent the cooling losses which shows that almost $85-90\%$ energy is radiated from the ISB. Therefore, the total energy retained in the ISB is $10-15\%$. One should note that this fraction may change depending on the density and metallicity of ISM. The lower two panels show the results for a two-fluid ISB. The third panel shows the adiabatic bubble with CR diffusion and fourth panel shows radiative two-fluid bubble with CR diffusion. A comparison between second and bottom (magenta colour) panels demonstrate that CRs suppress cooling losses.

5 Discussion

In this section we go beyond the fiducial models and explore the parameter space (sections 5.1, 5.2, 5.3). We show a comparison of one-fluid and two-fluid runs, in terms of the distribution of hot gas at different dynamical times (section 5.4). We also show the total energy gain by the CRs and discuss its dependence on various parameters (section 5.5).

5.1 Choice of CR parameters

To begin with, we explore the parameter space of CR injection parameter $\epsilon_{\rm cr}/w$ and diffusion constant ($\kappa_{\rm cr}$). The CR injection parameter $w$ (at shock, see equation (15)) or $\epsilon_{\rm cr}$ (at driving source, equation (16)) are crucial parameters of two-fluid ISBs, although their origin and values are not known (Bell 2014). We treat them as a free parameters. As the dependence of the result on $w$ is easily predictable in the one-fluid model (see section 4.2), here, we present the result for different values of $\epsilon_{\rm cr}$ and $\kappa_{\rm cr}$.

To visualize the dependence, we have estimated the ratio of volume averaged CR pressure to thermal pressure of the hot gas (temperature $>10^{5}$ K) and displayed them as a function of dynamical time in Figure 9. Any point above the black horizontal line represents a CR dominated ISB. Panel (a) shows that, for a larger value of $\epsilon_{\rm cr}$, the bubble becomes CR dominated at an early time. This is consistent because a larger $\epsilon_{\rm cr}$ increases the upstream CR pressure at the reverse shock. Panel (b) shows that, ISB becomes CR dominated if $10^{26}\lesssim\kappa_{\rm cr}/({\rm cm^{2}\,s^{-1}})\lesssim 3\times 10^{27}$. If $\kappa_{\rm cr}$ is below $<10^{26}\,{\rm cm^{2}\,s^{-1}}$, then CR diffusion is almost negligible resulting in an unaffected ISB. Whereas, as $\kappa_{\rm cr}$ increases, CRs diffuse out of the ISB, and therefore, a one-fluid ISB model (discussed in section 2.1) is good enough to describe its structure.

5.2 Dependence on the ambient density

The choice of the ambient density is an important parameter while comparing the theoretical result with observations. Most of the observations provide the density of the photo-ionized shell, but, beyond that, it is decipher to predict it from observations. Here we discuss the role of the ambient density for the two-fluid ISBs.

To scan the ambient density parameter space, we select two densities: $\rho=1\,m{\rm{}_{H}\,cm^{-3}}$ and $\rho=10^{3}\,m{\rm{}_{H}\,cm^{-3}}$ (recall that our fiducial density is $\rho=10^{2}\,m{\rm{}_{H}\,cm^{-3}}$). The nature of the ISB can be inferred from equation (19), which states that for a high ambient density, the ISB takes a longer time to attain a globally smooth solution. Simulations agree with analytic estimates as shown in panel (c) of Figure 9. Figure shows that for a lower ambient density affects the interior at an earlier time.

5.3 Dependence on $N_{\rm OB}$

The dependence on $N_{\rm OB}$ is quite clear from equation (19). A larger number of $N_{\rm OB}$ corresponds to a higher wind luminosity and mass-loss rate ($\tau_{\rm cri}\propto\dot{M}^{-5/8}L_{\rm w}^{-1/8}\sim N_{\rm OB}^{-3/4}$). This means that the reverse shock can become smooth at early times. The simulation results shown in panel (d) of Figure 9 are consistent with this. Therefore, superbubble reverse shocks, with high Mach numbers, are a promising site for CR acceleration.

5.4 Density-temperature of hot gas

Figure 10 shows the density-temperature distribution of hot gas. To estimate this, first, we divide the data into different density channels from $\rho=10^{-3}$ to $2\,m{\rm{}_{H}\,cm^{-3}}$ $(\delta\rho=10^{0.1}\,m{\rm{}_{H}\,cm^{-3}})$. For each density channel, we create temperature bins ($\delta T=10^{0.1}\,{\rm K}$ and $10^{5}\leq T/{\rm K}\leq 2\times 10^{8}$), and then, we calculate the hot gas volume within a given density-temperature bin. The normalization is such that the same colour corresponds to an identical volume fraction in all panels.

For all panels of Figure 10, cooling is on. The left most panel displays a one-fluid ISB without thermal conduction. The temperature of the gas is $\sim 10^{8}$ K and density is low (as expected from equations (11) and (10)). When thermal conduction is turned on (second panel from the left), the temperature drops to $10^{7}$ K. In this case, the mass density is high because of the mass evaporation from the dense shell (see section 3.4.1). Most interesting processes take place when we switch to a two-fluid ISB (third and fourth panels; $\kappa_{\rm cr}=5\times 10^{26}\,{\rm cm^{2}\,s^{-1}}$). In this case, even without thermal conduction (third panel), the temperature drops to $(1-5)\times 10^{6}$ K. The thermal conduction (right most panel) does not show any significant difference because the lowered temperature diminishes the effect of thermal conduction. Figure shows that the CR affected ISBs can have temperature much lower than that of a one-fluid ISB.

5.5 Total energy gain by CRs

In previous sections, we have seen that for two-fluid ISBs, CR can gain energy from thermal particles. Here we show the energy gain by CRs and plot it as a function of dynamical time. To the best of our knowledge, this is the first time that the net gain of CR energy ($\Delta E^{\rm t}_{\rm cr}$) has been presented.

To obtain CR energy gain, we estimate $\Delta E^{\rm t}_{\rm cr}$ by following the same method as discussed in section 4.4. The fraction of energy gain by the CRs $\Delta E^{\rm t}_{\rm cr}/E^{\rm t}_{\rm IN,cr}$ (where, $E^{\rm t}_{\rm IN,cr}=\epsilon_{\rm cr}E^{\rm t}_{\rm IN}$, see equation (35)) and its dependence on various parameters are displayed in Figure 11. To prevent CR overflow from the computational box (especially for the low density ambient, large $\kappa_{\rm cr}$ and $\epsilon_{\rm cr}$), we set $\rm r_{max}=1400.1$ and choose number of grid points $n=28000$ (Table 2). We have tested the total energy efficiency of the simulation box which is defined as

[TABLE]

where ${\rm x}$ is the kinetic/thermal/CRs energy and $E_{\rm loss}$ is the radiative energy loss (also see equation (35)). All the runs fulfill the energy conservation with an accuracy $>97\%$, shown by the horizontal lines close to unity in Figure 11.

The panel (a) of Figure 11 shows that the CR energy gain fraction increases as $\epsilon_{\rm cr}$ decreases. This is reasonable because at the shock the upstream kinetic energy gets converted to the downstream CR energy via CR diffusion (although the energy transfer is insufficient for it to become a CR dominated ISB, see panel (a) in Figure 9). Panel (b) shows that CRs gain energy if $10^{26}\lesssim\kappa_{\rm cr}/({\rm cm^{2}\,s^{-1}})\lesssim 3\times 10^{27}$, consistent with the conclusion of section 5.1. Panel (c) shows that, an ISB expanding in a low density ISM achieves maximum energy at an early time. Panel (d) confirms that a large $N_{\rm OB}$ can be a promising source of CRs.

6 Astrophysical implications

X-ray diagnostic is an important tool to study the interior of ISBs. Since the bremsstrahlung emissivity decreases exponentially above the gas temperature, one can find the temperature of the plasma from the X-ray spectrum. This method is very useful to determine the best-fit temperature of the X-ray emitting gas (Dunne et al. 2003; 2006; Lopez et al. 2011). The shape of spectra also tells whether emissions are coming from the thermal particles or from non-thermal particles (Dunne et al. 2003). The non-thermal emission can be used for modeling CR energy (Helder et al. 2009).

For most of the observed ISBs, the swept-up shell is quite clumpy which allows the X-ray radiation to come out, and helps one to estimate the hot gas density. The density of the hot gas is mostly obtained from the emission measure. This technique requires assuming a hot gas volume filling factor, which is taken to be $1/2-1$. However, we should note that this may introduce an error in the hot gas density.

The best fit temperature ($T_{\rm x}$) and the density ($n_{\rm x}$) of the X-ray emitting hot gas in the observed ISBs are given in Table 3. All these data points are taken from published results (references given in the caption of Table 3). Figure 12 visualizes the content of this table. For all these ISBs, the dynamical age is less then $10$ Myr (Lopez et al. 2014). This figure shows that the number density of X-ray gas is $n_{\rm x}=(0.01\hbox{--}0.3)\,{\rm cm^{-3}}$ and the temperature is $\approx(1\hbox{--}9)\times 10^{6}$ K.

If one considers that thermal conduction is the only physical process controlling the interior temperature, then the one-fluid model can be used to explain some of the observed data points. However, as we have mentioned in section 2, thermal conduction can be affected by magnetic fields. Alternatively, CRs can be responsible for lowering the temperature of the hot gas which diminishes the mixing between the ambient and ejecta material via evaporation. Therefore, if an observation finds that the temperature of the hot gas does not agree with the wind velocity (see equation (11)) and the hot gas material is different from the ambient material (usually done by determining metallicity), then that would be a promising evidence for a CR affected ISB.

A comparison of simulation results with observations shows that the temperature matches well if CR is injected at the driving source region with a CR energy fraction $\epsilon_{\rm cr}\sim 10-20\%$ of total input energy and a CR diffusion constant $\kappa_{\rm cr}\approx(0.1-3)\times 10^{27}\,{\rm cm^{2}\,s^{-1}}$. A higher value of CR diffusion constant decouples the CR energy and the thermal energy resulting in an unmodified superbubble. Therefore, if the diffusion constant is very high $(\gg 10^{27}\,{\rm cm^{2}\,s^{-1}})$, then, regardless of whether or not ISBs are sites of CR acceleration, the thermal X-ray temperature will not be a good diagnostic of the presence of CR.

It is worth mentioning that a one-dimensional hydrodynamic analysis can not capture magnetized thermal conduction and cosmic ray diffusion. These need further investigation have been left for future work.

7 Summary

In this paper we have presented a two-fluid model of the interstellar bubbles by considering CR as a second fluid. Our work can be seen as a generalization of two standard theories of outflows (1) CR affected blast wave (first modeled by Chevalier 1983) and (2) one-fluid interstellar bubble ( 1977). The main results from this work are given below:

Dynamics without CRs: We have found that the thermal pressure inside the bubble follows a robust relation which holds for an arbitrary density profile (equation (3)) even with radiative cooling and thermal conduction. According to this, the volume averaged pressure inside an ISB is $\approx\rho\,v_{\rm sh}^{2}$ (see equation (7)), where $v_{\rm sh}$ is the velocity of the expanding shell and $\rho$ is the ambient density. Therefore, the deviation from this relation can be considered as an indication of the presence of CRs (see Figure 7). 2. 2.

CR injection and its effect: The effect of CRs depends on (1) where the seed relativistic particles are injected and (2) what fraction of the total injected energy/post shock pressure goes into it. Since an ISB consists two shocks, one can inject CRs at the forward and/or at the reverse shock. The injection at the forward shock does not change the interior structure (panel (a) in Figure 2). This makes ISBs different from a two-fluid SNe shock because the center of the blast wave is dominated by CRs if the injection is done at the shock (Chevalier 1983; also see Figure $1$ in Bell 2014). The injection at the reverse shock can reduce the thermal pressure inside the ISB. The CR injection at shocks is described by an ad-hoc parameter (denoted by $w$, definition is given in section 4.2) which does not capture the actual physics of two-fluid shock. However, a self-consistent evolution is obtained when CRs are injected in the driving source region (panel (c) in Figure 2). In this case, depending on the CR diffusion constant, the reverse shock can show all possible solutions of a two-fluid shock predicted by Drury & Volk 1981 (Figure 5). 3. 3.

The importance of CR diffusion constant : The key parameter in two-fluid model is the CR diffusion constant ($\kappa_{\rm cr}$). One can see a significant difference between the ISBs with and without CRs only if $\tau_{\rm acc}\ll\tau_{\rm cri}<t_{\rm dyn}$ (section 4.2.2) where $\tau_{\rm acc}$ is the CR acceleration time scale (equation (32); analogous to diffusion timescale), $\tau_{\rm cri}$ (equation (19)) is the time taken by the reverse shock to exceed the the critical Mach number (Becker & Kazanas 2001) for a globally smooth reverse shock profile and $t_{\rm dyn}$ is the dynamical time. We have found that this condition is fulfilled if $10^{26}\lesssim\kappa_{\rm cr}/({\rm cm^{2}\,s^{-1}})\lesssim 3\times 10^{27}$ and in the case of ISBs with a large number OB stars ($N_{\rm OB}\geq 10^{2}$) (panels (b) and (d) in Figures 9 and 11). This supports the suggestion in the literature from phenomenogical studies, that massive compact stellar associations can be promising source of CRs (Higdon & Lingenfelter 2006; Ferrand & Marcowith 2010; Lingenfelter 2012). 4. 4.

Observational signatures: An indirect evidence for the presence of CRs in ISBs can be inferred from the temperature of the X-ray emitting plasma. We have found that the CR affected bubble can have temperature $1-5\times 10^{6}$ K even in the absence of thermal conduction (Figure 10) which can explain the X-ray temperature in the observed ISBs.

The model presented in this paper is admittedly idealized, which can be extended to a more realistic scenario. $3$D MHD simulations of two-fluid ISBs will be important to shed further light on the question of CR origin.

Acknowledgements

DE acknowledges support from an ISF-UGC grant, the Israel-US Binational Science Foundation, and the Joan and Robert Arnow Chair of Theoretical Astrophysics. PS acknowledges an India-Israel joint research grant (6-10/2014[IC]). SG thanks SPM fellowship of CSIR, India for financial support.

Appendix A Code check

We have performed several standard test problems to check our code TFH. Here we present two of them: (1) the shock tube problem (Sod 1978) in cartesian geometry (section A.1) and (2) the blast wave problem (Sedov 1946; 1950) in spherical geometry (section A.2). For both cases, we present one-fluid and two-fluid solutions. We also present a test problem for the diffusion module (thermal conduction/CR diffusion) in section A.3. We have compare our results with analytical solutions and also with the publicly available one fluid code PLUTO ( 2007).

A.1 Shock-tube

This test problem is identical to the problem described in Sharma et al. (2009) (also see Sod 1978; Pfrommer et al. 2006). Problem set-up : We set geometry to cartesian coordinate and choose total $2000$ grid points in a domain $[1,2]$. As the initial condition, the left state ($x<1.5$) is defined as $(\rho,v,p_{\rm th},p_{\rm cr})_{\rm L}=(1,0,6.7\times 10^{4},1.3\times 10^{5})$. The right state ($x\geq 1.5$) is defined as $(\rho,v,p_{\rm th},p_{\rm cr})_{\rm R}=(0.2,0,2.4\times 10^{2},2.4\times 10^{2})$. The CFL number is set to $0.4$. The snap shot of various profiles are shown in Figure 13. The left panels display the profiles for one fluid shock at $t=1\times 10^{-3}$ and the right panels display the profiles for two-fluid shock at $t=5\times 10^{-4}$.

A.2 Blast wave

Problem set-up : We set geometry to spherical coordinate and choose total $1000$ grid points in a domain $[0.1,4]$ pc. The initial profiles are: $\rho=1\,m{\rm{}_{H}\,cm^{-3}}$ (uniform), $v=0.0$, $p_{\rm th}=(\rho/\mu m_{\rm H})k_{\rm B}T$ where $\mu=0.6$ and $T=10^{4}$ K. At $t=0$, at the first computational zone (say, volume $\delta V_{\rm src}$), we set energy density $10^{51}/\delta V_{\rm src}$. The CFL number is set to $0.2$. The snap shot of various profiles are shown in Figure 14. Left and right panels display the profiles of one fluid and two-fluid blast wave where all variables are scaled w.r.t to the self-similar variables (Chevalier 1983). For the two-fluid run, first, we identify the shock and set the CR pressure fraction to $w=1/2$ (see equation (15)). The profiles match quite well with the ODEs results (shown by black curves).

A.3 Diffusion module

One useful test problem to check diffusion module was proposed by Reale et al. (1995).

Problem set-up : Recall their equations ($15$) and ($16$). We set $n=5/2$ , $a=4.412$ and $Q=1.2\times 10^{15}$ K cm. The length and temperature units are defined as $10^{8}$ cm and $\simeq 6.057\times 10^{7}$ K. Total $400$ grid points are uniformly set in a domain [$0,5$]. The initial time is chosen as $t=0.1$ and the simulation is run upto $t=3.1$. For a detailed setup, see $\rm PLUTO/Testproblem/MHD/Thermal\_Conduction/TCfont$ ( 2007). We turn on only the diffusion module and our results are shown in Figure 15. For cartesian coordinate, we have compared the results with the analytical solution (black curves in left panel) given in equation ($16$) of Reale et al. (1995). The comparison between PLUTO and TFH is also shown in Figure 15.

Appendix B Solver selection

TFH has two solvers (1) TS and (2) HLL (see section 3). Here we show a comparison of density profile of an adiabatic one-fluid ISB at $t=0.5$ Myr between HLL/TS solver of TFH and HLL solver of PLUTO in Figure 16. The set-up is identical for all cases. This figure shows that 1st order HLL solver is more diffusive than 1st order TS solver explaining the reason for selecting ‘TS’ solver.

Appendix C Resolution test

Here we present the convergent tests of a two-fluid ISB by considering four different spatial resolutions: $\Delta r=0.1,0.05,0.025,0.0125$ pc. Figure 17 displays the thermal and CR pressure profiles of an adiabatic two-fluid ISB (with CR diffusion) at $t_{\rm dyn}=0.9$ Myr and $t_{\rm dyn}=2.1$ Myr. The blue colours in both panels denote our fiducial resolution. This figure displays that the results are converged only for high resolutions ($\Delta r<0.1$ pc). We have found that, the time when $P_{\rm cr}>P_{\rm th}$ increases with the decrease in spatial resolution (drastically when $\Delta r\gtrsim 0.2$ pc). Resolution test for a more realistic bubble (i.e.; in addition to CR diffusion, cooling is on) is shown in Figure 18. For our fiducial choice (resolution $\Delta r=0.05$ pc), the results are well converged and the conclusions remain same.