Constraints on cosmic rays in the Milky Way circumgalactic medium from OVIII observations

Galaxies have two components: galactic disk surrounded by a gaseous and a dark matter halo. The gaseous halo, which extends up to the virial radius (sometimes even beyond) is known as the circumgalactic medium (CGM). It is a reservoir for most of the baryons and plays a crucial role in galaxy formation and evolution by various feedback processes such as outflowing, recycling of gas (Tumlinson et al., 2017). Soft $X$-ray observation of OVII and OVIII absorption (Gupta et al., 2012; Fang et al., 2015) and emission lines (Henley et al., 2010; Henley & Shelton, 2010), and SZ effect (Planck Collaboration et al., 2013; Anderson et al., 2015) indicate the presence of a hot phase (T $\geq 10^{6}$ K) of CGM. Recently, a warm ($10^{5}\rm{K}<T<10^{6}$K) and a cool phase ($10^{4}\rm{K}<T<10^{5}$K) of CGM have also been discovered through absorption lines of low and intermediate ions ¹¹1Low Ions: Ionization Potential (IP) $<40$ eV, T $=10^{4-4.5}$ K; Intermediate Ions: $40\gtrsim$ IP (eV) $\lesssim 100$, T $=10^{4.5-5.5}$ K; High Ions: IP$\gtrsim 100$, T$>10^{5.5}$K. (Tumlinson et al., 2017) at low redshifts (Stocke et al., 2013; Werk et al., 2014, 2016; Prochaska et al., 2017a) and Ly$\alpha$ emission at high redshifts (Hennawi et al., 2015; Cai et al., 2017). It is not yet clear how this cool phase coexists with the hot phase and survives the destructive effects of various instabilities (McCourt et al., 2015; Ji et al., 2018), but this discovery has led to a picture of multiphase temperature and density structure of the CGM (Tumlinson et al., 2017; Zhang, 2018). It also partially solves the problem of missing baryon in the galaxies (Tumlinson et al., 2017).

If the CGM is considered to be in hydrostatic equilibrium by means of only thermal pressure, then in order to maintain pressure equilibrium the cold component ($\sim 10^{4}$ K) of CGM is expected to have a higher density than the hot component ($\sim 10^{6}$ K). However, Werk et al. (2014) found the density of cold gas to follow the hot gas density distribution. This has given rise to the idea of a non-thermal pressure component, consisting of cosmic-ray (CR) and magnetic pressure which can add to the thermal pressure. Without a non-thermal component, the abundances of low and mid-ions (seen in cool and warm phases, respectively) are under-estimated even in high-resolution simulations (van de Voort et al., 2019; Hummels et al., 2019; Peeples et al., 2019). Ji et al. (2020) have suggested that CR can explain these two problems. It has been suggested that CR can provide pressure support to cool diffuse gas (Salem et al., 2015; Butsky & Quinn, 2018), help in driving galactic outflows (Ruszkowski et al., 2017; Wiener et al., 2017), as well as excite Alfvén waves that can heat the CGM gas (Wiener et al., 2013).

The magnitude of the CR population in the CGM is, however, highly debated. Some simulations (Butsky & Quinn, 2018; Ji et al., 2020; Dashyan & Dubois, 2020) claimed a CR dominated halo, which increases feedback efficiency of the outflowing gas by increasing the mass loading factor and suppressing the star formation rate. The ratio of CR pressure and thermal pressure ($\rm{P}_{\rm{CR}}/\rm{P}_{\rm{th}}=\eta$) controls this effect and can reach a large value exceeding 10 over a huge portion of Milky-way sized halo (Butsky & Quinn, 2018). Another simulation claimed that MW-sized halo at low redshift (z$<1$) can be CR populated with $\eta$ nearly $10$ in outflow regions, although in warm regions (T=$10^{5}$ K) of halo, $\eta\leq 1$ (Ji et al., 2020). Dashyan & Dubois (2020) found in their simulation that dwarf galaxies, with virial mass $10^{10}-10^{11}$ M_⊙, can have a high value of $\eta$ ($\sim 100$) at $1-5$ kpc from the mid-plane when isotropic diffusion is incorporated (their figure 1 and 2 ).

It is therefore an important question as to how much CR can be accommodated in the CGM in light of different observations. One of the most obvious effects of CR population in CGM is to suppress the thermal pressure of the hot phase which is also seen in the simulation of Ji et al. (2020). The recent work by Kim et al. (2021) pointed out that this reduction in thermal pressure would lower the thermal SZ (tSZ) signal, as tSZ probes the integrated thermal pressure in the halo. This diminished value of thermal pressure would also significantly change the abundance of high ionization species, and thereby jeopardise the interpretation of their column densities. In a recent work by Faerman et al. (2021), they modified their previous isentropic model Faerman et al. (2019) with significant non-thermal pressure ($\alpha=(P_{nth}/P_{th})+1.0=2.9$) and found that the gas temperature in the central region for the non-thermal model becomes lower than the value at which OVIII CIE temperature peaks, which in turn results in a lower OVIII column density. At low CGM masses, photoionization compensated for this effect by the formation of OVIII at larger radii, but for large gas masses and large mean densities (as in the Milky Way), photoionization effect is negligible and the total OVIII column density is low. In a previous work (Jana et al., 2020), we constrained the CR pressure ($\eta=P_{CR}/P_{th}$) in the CGM in light of the isotropic $\gamma$-ray background (IGRB) and radio background, using different analytical models (Isothermal Model: IT, Precipitation Model: PP) of CGM. In the case of IT model, the value of IGRB flux puts an upper limit of 3 on $\eta$, whereas all values of $\eta$ are ruled out if one considers the anisotropy of the flux due to off-center position of solar system in Milky Way. However, in the case of PP model, IGRB flux value and its anisotropy allow a range of $\eta$ from 100 to 230. In comparison, the constraints from radio background are not quite robust. In this paper, we use the same analytical models (IT and PP) and put limits on CR population in the CGM of Milky Way by comparing the predicted OVIII absorption column densities (N_OVIII) with the available observations.

It might still be asked, why use OVIII as a probe? Previous works (Faerman et al., 2017; Roy et al., 2021) explained the observed OVII, OVIII absorption column densities in the Milky Way without the incorporation of CR component using the analytical models of CGM used here. It is, therefore, important to study the effect of CR on these column densities with the inclusion of CR component in these models. Whereas the OVII ion has a plateau of favourable temperatures ($5\times 10^{5}$ K - $1\times 10^{6}$ K), the suitable temperature for OVIII production ($2\times 10^{6}$K) peaks near the temperature of the hot halo gas of Milky Way CGM. This particular fact makes OVIII abundance a sensitive probe of the CGM hot gas and motivates the choice of OVIII column density for comparison with the observations in order to constrain CR population in Milky Way CGM.

We study two widely discussed analytical models: 1) Isothermal Model (IT), 2) Precipitation Model (PP). These models have previously been studied in the context of absorption and emission lines from different ions (Faerman et al., 2017; Voit, 2019; Roy et al., 2021) We assume the CGM to be in hydrostatic equilibrium within the potential of dark matter halo ($\phi_{\rm DM}$) , with a metallicity of $0.3$ Z_⊙ (Prochaska et al., 2017b), so that,

where $P_{\rm total}$ is the total pressure and $\rho$ and r are the density and radius, respectively. We consider the Navarro-Frenk-White (NFW) (Navarro et al., 1996) profile as the underlying dark matter potential in the IT model. However we slightly modify this potential in the case of the PP model as suggested in Voit (2019) by considering the circular velocity (v_c) to be constant (v${}_{\rm c,max}=220$ km s^-1) up to a radius of $2.163$ $\rm r_{\rm vir}/\rm c$ for a halo with M${}_{\rm vir}=$ $2\times 10^{12}$ M_⊙ and concentration parameter $c=10$. We include non-thermal pressure support with CR population and magnetic field along with thermal pressure in these models. We consider the magnetic energy to be in equipartition with thermal energy (P_mag=$0.5$P_th) so that the total pressure $\rm{P}_{\rm total}=\rm{P}_{\rm{th}}(1.5+\eta)$, which leads to,

Note that the magnitude of magnetic field in CGM is rather uncertain. While Bernet et al. (2008) claimed the magnetic field in CGM of galaxies at $z=1.2$ to be larger than that in present day galaxies, the observations of Prochaska et al. (2019) suggested a value less than that given by equipartition. Yet another study claimed near-equipartition value after correlating the rotation measures of high-$z$ radio sources with those in CGM of foreground galaxies (Lan & Prochaska, 2020).

The inclusion of non-thermal pressure suppresses the thermal pressure, and in turn the temperature or density or both of the CGM gas (also seen in recent simulation of Ji et al. (2020)). Our goal is to determine the effect of the inclusion of CR population on N_OVIII therefore constrain the CR population in Milky Way CGM. In addition to models with uniform $\eta$, we have studied the effect of varying $\eta$ as a function of radius. For this, we have explored two contrasting variations, $\eta=A\times x$ (rising profile) and $\eta=A/x$ (declining profile) where x=r/r_s, where A is normalization of the profiles. We primarily studied the case with A=1, where the maximum value (for $\eta=x$) is r${}_{vir}/r_{s}$=c=10. For $\eta=1/x$, to avoid divergence at r=0 kpc we have started our calculation from r=1 kpc, therefore the maximum value in this case is r${}_{s}/1$ kpc$=26$. We have shown the profiles of $\eta=x,1/x$ in the figure 1. We have explored the effects of scaling up and down the values with these profiles as well (e.g., $\eta=Ax,B/x$, with $A,B\geq 1$). Although $\eta$ can vary in a more complicated manner, we study these two particular cases (increasing and decreasing linearly) as the simplest representatives of a class of models in which $\eta$ is allowed to vary with radius. It should be noted that the $\eta$ here denotes the CR pressure in the diffuse hot gas with respect to thermal pressure of the hot gas.

There are several ways to proceed, starting from equation 2. One way is to simply consider the temperature all over the CGM to be constant, which implies the LHS of the equation 2 to be zero. Another way is to use specific entropy to relate the two unknown quantities in the equation 2, temperature and density. We discuss these two ways in details in the following subsections.

We calculate the ionization fraction of OVIII using CLOUDY (Ferland et al., 2017) with the input of density and temperature profiles derived from the CGM models. In order to incorporate temperature fluctuation, we calculate the ionization fraction for all the temperatures in a log-normal distribution using CLOUDY and integrate them over the corresponding log-normal distribution to get a mean ionization fraction. That means, for IT model we get a single mean ionization fraction corresponding to the log-normal distribution around the single temperature of the halo. But for PP model, we will get mean ionization fractions at each radius corresponding to the log-normal distributions with constant value of $\sigma_{\rm lnT}$ around the mean temperature $T(r)$ at each radius. We consider two cases: collisional ionization and photoionization along with collisional ionization. We use the extragalactic UV background (Haardt & Madau, 2012) at redshift $z=0$ for photoionization. We take into account our vantage point of observation i.e the solar position and calculate the variation of column density with the Galactic latitude and longitude. We consider the median of the column densities for each value of $\eta$ in order to compare with the observations.

In Figures 3 (for IT model) and 4 (for PP model), we show the variation of N_OVIII with $\eta$, considering collisional ionization (CI) and photoionization (PI) for different models. One can clearly see from these figures that with an increase in $\eta$, the value of N_OVIII decreases due to the decrease in the temperature. In Figure 3, the brown and blue lines denote the cases with CI and PI respectively in the IT model. In this figure, the orange colours denote the effect of log-normal fluctuations ($\sigma_{\rm lnT}=1.0$). The cases of CI and PI do not differ for $\eta\leq 1$, since at $T\sim 10^{6}$K, CI dominates over PI. However, for $\eta\geq 1$, PI leads to slightly larger values of OVIII column density than the CI case, as further decrease in temperature leads to more production of OVIII by PI. It should be noted that if one considers temperature fluctuation around a single favourable temperature of OVIII, there is less OVIII production with the increase in $\sigma_{\rm lnT}$.

In Figure 4 the orange shaded region indicates the boundary temperature near to virial temperature of halo $1.8\times 10^{6}$ K (T_bc1) in PP model, whereas the blue shaded region denotes the boundary condition of $8.8\times 10^{5}$ K (T_bc2) at virial radius. For the PP model, we show the cases with CI without ($\sigma_{\rm lnT}=0.0$) and with fluctuation ($\sigma_{\rm lnT}=1.0$) and PI in left, middle and right panel of Figure 4. The shaded regions in this figure refer to the CGM mass range mentioned earlier in the case of PP model. A boundary temperature that is close to the favourable temperature for OVIII yields more OVIII than one gets with lower boundary temperatures. However, with the inclusion of temperature fluctuations, one gets almost similar OVIII for both cases. With lower boundary temperature, the production of OVIII decreases for $\eta\geq 2$, because the inclusion of $\eta$ shifts the temperatures from the CI peak of OVIII. However, considering PI in this case can increase the amount of OVIII.

In Figure 5, we show N_OVIII with different $\eta$ profiles for both of the models. For IT model, we use brown and orange colours to denote $\sigma_{lnT}=0.0$ and $1.0$ respectively whereas circles denote $\eta=$ x and squares denote $\eta=$ 1/x. For PP model, different temperature boundary conditions are shown by orange (T${}_{bc}=1.8\times 10^{6}$ K) and red (T${}_{bc}=8.8\times 10^{5}$ K) shaded regions where shaded region (above the lines for $\eta$=x and below the lines for $\eta$=1/x) denote the mass range for each case in PP model. One important point to notice here is that declining $\eta$ profiles produce more OVIII than rising profiles.

Our constraints allow larger values of $\eta$ for most of the cases in comparison to the isentropic model by Faerman et al. (2019, 2021). In their recent work (Faerman et al., 2021), they have considered three cases for their model : 1) with only thermal pressure i.e, $\alpha=(P_{nth}/P_{th})+1.0=((P_{CR}+P_{B})/P_{th})+1.0=1.1$, 2) with standard case as considered in Faerman et al. (2019) i.e, $\alpha=2.1$ and 3) with significant non-thermal pressure where $\alpha=2.9$. Note that the the above mentioned values of $\alpha$ are the values at the outer boundary and the profiles of $\alpha$ will follow a declining pattern as radius decreases (see blue curve in the right panel of Figure 2 by Faerman et al. (2019)). With our definition of $\eta$, along with the magnetic field value used by us, we can convert these $\alpha$ values to $\eta\sim 0,0.6,1.4$, respectively, for their three cases. The thermal case matches N_OVIII observations up to CGM mass of $10^{11}\rm M_{\rm solar}$, and for the standard case, the value of N_OVIII is comparable to the observed value. However, for the significant non-thermal case, the OVIII value is lower than the observations by a factor of $\sim 10$ for a CGM mass of $10^{11}\rm M_{\rm solar}$. Therefore, one can say that $\eta<0.6$ is allowed for the isentropic model. For PP model, we get upper limit of $\eta\sim 0.5$ for two cases which are therefore in agreement with the constraint from the isentropic model. However, for most of our models, the upper limits derived in the present work are larger than the isentropic model.

It should be noted that the limits of $\eta$ derived here translate to constraints on cosmic ray pressure in the hot diffuse CGM, and one may wonder about the CR pressure in the cold gas. CR pressure scales with density ($\rm P_{\rm CR}\propto\rho^{\gamma_{c,eff}}$, where $\gamma_{c,eff}$ is the effective adiabatic index of the CR) of the gas, where the adiabatic index depends on the transport mechanism (Butsky et al., 2020), the limits of CR pressure in the cold gas would be different from what we have found for the hot gas. If CRs are strongly coupled to the gas, then in the limit of slow CR transport, i.e, if the only CR transport mechanism is advection, $\gamma_{c,eff}=\gamma_{c}=4/3$. In this case, CR pressure can be higher in cooler, denser gas than in the hot gas. However, in the limit of efficient CR transport i.e, if there are CR diffusion, streaming along with advection, the CR pressure will be redistributed from high density, cold region to diffuse, hot region of CGM which will lower the $\gamma_{c,eff}$. In this limit, $\gamma_{c,eff}\rightarrow 0$, and the limits of CR pressure in the cold phase are nearly equal to the CR pressure in hot phase. In the context of a single phase temperature with temperature gradient, if one considers CR pressure in the hot phase to scale as $\rho_{g}^{\gamma_{c,eff}}$ then $\eta$ is proportional to $\rho_{g}^{\gamma_{c,eff}-(5/3)}$, as adiabatic index for gas is $5/3$. Then $\eta$ will be proportional to $\rho_{g}^{-1/3}$ and $\rho_{g}^{-5/3}$ respectively depending on slow and effective transport mechanism. This implies that $\eta$ is smaller in denser gas and hence mimic a rising profile of $\eta$ in case of the single hot phase gas, which we have shown in the present paper.

However, these scaling relations are for adiabatic situations, which is not the case for PP model which involves energy loss. In the context of cold gas that may have condensed out of the hot gas, it is possible for the cold gas to contain a significant amount of CR pressure, with a large value of $\eta$ for cold gas. Simulations by Butsky et al. (2020) (their figure 10) showed that in the limit of slow cosmic ray transport, the value of $\eta$ is quite different for hot and cold gas. In the limit of efficient transport, cosmic ray pressure is decoupled from the gas, and $\eta$ has similar values in hot and cold gas. If we recall the observation by Werk et al. (2014), that the cold phase of CGM nearly follows the hot gas density profile then for small density contrast, the CR pressure in the cold phase and hot phase may not differ much even with different transport mechanisms. But even with small density contrast, thermal pressure of cold gas will be smaller than hot phase by two order of magnitude due to the temperature difference, which allows for a larger upper limit on $\eta$. However, packing more CR in cold gas will not have any effect on the constraints derived here. In fact, this may be one way to circumvent the constraints on $\eta$ described here.

We have not included any disc component in the present work. A recent work by Kaaret et al. (2020) took into account an empirical disc density profile motivated by the measured molecular profile of MW along with the halo component. They pointed out that the high density disc model can well fit the MW’s soft X-ray emission, whereas it under predicts the absorption columns. Their conclusion is that the dominant contribution of X-ray absorption comes from the halo component due to large path length. Also this difference in contribution comes from the fact that the absorption is proportional to density whereas emission measure is proportional to density square. Therefore, considering only the halo component and making the comparison with the absorption column is justified.

Note that the limits derived here are sensitive to some of the parameters such as the magnetic field and metallicity. Increase (decrease) in the magnetic pressure would increase (decrease) the total non-thermal pressure, which would decrease (increase) the limit on $\eta$. However, given the uncertainty in the magnetic field strength in the CGM, equipartition of magnetic energy with thermal energy is a reasonable assumption. For a collisionally ionized plasma, which we assume for the CGM gas, the column density is proportional to the metallicity ($Z$) in the case of IT model, whereas it is proportional to $Z^{0.3}$ for the PP model (see eq 17 of Voit (2019)). In the photoionized case, the dependence on metallicity will be stronger since the ion abundance would also depend on the metallicity. Also note that the limits derived here refer to CR population spread extensively in the CGM. Localized but significant CR population may elude the above limits. However, we also note that most of the contribution of the limits come from the inner region, where OVIII exists otherwise and is suppressed due to CR. Therefore, our limits also pertain to localized CR populations in the inner regions.

We have studied two analytical models of CGM (Isothermal and Precipitation model) in hydrostatic equillibrium in light of the observations of N_OVIII in order to constrain the CR population in the Milky Way CGM. We found that $\eta\leq 10$ in the case of no photoionization and including temperature fluctuations in the PP model whereas IT model allows $\eta\leq 6$ with inclusion of PI. However, it should be noted that IT is an extremely simplistic model and the inclusion of non-thermal component leads to a very low density. The limits from IT may therefore be only of academic interest.

Combined with the limits derived from $\gamma$-ray background Jana et al. (2020), our results make it difficult for a significant CR population in our Milky Way Galaxy. However, models in which the ratio of CR to thermal pressure varies with radius (preferably declining with radius) can be accommodated within observational constraints, if they include suitable temperature fluctuation (larger $\sigma_{\rm lnT}$ for PP and lower for IT), although making allowances for large value of $\eta$ only in the central region.