Effect of cosmic rays and ionizing radiation on observational ultraviolet plasma diagnostics in the circumgalactic medium

The relevance of some galactic feedback mechanisms, in particular cosmic ray feedback and the hydrogen ionizing radiation field, has been challenging to definitively describe in a galactic context, especially far outside the galaxy in the circumgalactic medium (CGM). Theoretical and observational uncertainties prevent conclusive interpretations of multiphase CGM properties derived from ultraviolet (UV) diagnostics. We conduct three dimensional magnetohydrodynamic simulations of a section of a galactic disk with star formation and feedback, including radiative heating from stars, a UV background, and cosmic ray feedback. We utilize the temperature phases present in our simulations to generate Cloudy models to derive spatially and temporally varying synthetic UV diagnostics. We find that radiative effects without additional heating mechanisms are not able to produce synthetic diagnostics in the observed ranges. For low cosmic ray diffusivity $\kappa_{\rm{cr}}=10^{28} \rm{cm}^2 \rm{s}^{-1}$, cosmic ray streaming heating in the outflow helps our synthetic line ratios roughly match observed ranges by producing transitional temperature gas ($T \sim 10^{5-6}$ K). High cosmic ray diffusivity $\kappa_{\rm{cr}}=10^{29} \rm{cm}^2 \rm{s}^{-1}$, with or without cosmic ray streaming heating, produced transitional temperature gas. The key parameter controlling the production of this gas phase remains unclear, as the different star formation history and outflow evolution itself influences these diagnostics. Our work demonstrates the use of UV plasma diagnostics to differentiate between galactic/circumgalactic feedback models.

Here we focus on Milky Way-like galaxies, where stellar feedback (including cosmic rays) is likely to play a role in galactic evolution (i.e., Hopkins et al. 2012aHopkins et al. , 2014;;Ruszkowski et al. 2017;Habegger et al. 2022;Butsky et al. 2023).Radiation and CR phenomena are similar in that small-scale physics, much below the resolution scale of a galactic simulation, plays a key role in global transport.For example, two opposing limits-the free-streaming and diffusive limits-must be adequately modeled (e.g., Jiang & Oh 2018;Thomas & Pfrommer 2019;Hopkins et al. 2022).The stellar radiation field likely extends well past the galactic disk and into the inner CGM up to roughly 0.1-0.2 vir for Milky Way-like, low redshift galaxies (Holguin et al. prep).Therefore, it is important to con-sider the influence of this spatially and temporally varying field on the ionization state in the inner CGM.Three dimensional radiative transfer is incredibly computationally expensive as the radiation field acts non-locally.
CR physics is challenging to model as there are large uncertainties in the underlying assumptions about the interaction of CR and the magnetized environment.CRs are known to have a galactic confinement time much larger than the free streaming time and their arrival distribution measured at the Earth for Tev-PeV cosmic rays is measured to be isotropic to 10 −4 (Nagashima et al. 1998).CR confinement is primarily described by two ideas, extrinsic turbulence and self-confinement (see Zweibel 2013).Within both of these descriptions, many models have been developed and even implemented in galactic simulations (see Zweibel 2013Zweibel , 2017;;Thomas & Pfrommer 2019).However, distinguishing between these models based on observational evidence of CRs, such as the total column density of matter traversed by CRs ('grammage') (Strong et al. 2007), radio synchrotron (e.g., the CHANG-ES radio halo survey Irwin et al. 2012) or gamma-rays (e.g., Ackermann et al. 2012) has been difficult.Hopkins et al. (2021c) compares a wide range of CR transport parameters in dwarf and MW-type galaxy simulations and identifies a set of models that are consistent with observational evidence.These allowed models of CR transport, including more physically motivated ones, still have significant uncertainties.They require some fine-tuning, tend to produce similar effects on galaxy and inner-CGM (<10 kpc) properties, and are not well-constrained observationally farther out into the CGM (Hopkins et al. 2021b).These models are also problematic with regards to constraints from the CR spectrum (Hopkins et al. 2021a).The large uncertainty of CR transport motivates the introduction of additional constraints.
Ultraviolet (UV) ion diagnostics contain a wealth of information about the plasma state of the CGM from the cold, neutral medium to the hot, ionized galactic corona (see Tumlinson et al. 2017, and references therein).These diagnostics have highlighted the CGM as a significant reservoir of baryons existing in a multiphase medium.The exact properties of the plasma that produces some of these observed diagnostics is uncertain.Untangling the contributions from unresolved gas components requires careful statistical analysis of the spectra, including accounting for multiple, co-spatial components (Cooper et al. 2021).In particular, there is evidence based on the ionization state of metals that there exists a  ∼ 10 5 K phase, in between classical ISM and coronal stable phases (McKee & Ostriker 1977;Wolfire et al. 2003).The existence of a 10 5 K phase is perplexing as it sits near the peak of the (solar metallicity) cooling curve (Sutherland & Dopita 1993).Unfortunately, modeling of this (thermally) unstable phase is complicated by degeneracies in modeling for local ionization equilibirum; similar diagnostic signals can occur from collisionally ionized  ∼ 10 5 K and photoionized  ∼ 10 4 K plasmas.The collisional ionization explanations typically invoke additional heating physics, such as (radiative) turbulent mixing layers (e.g.Begelman & Fabian 1990;Kwak & Shelton 2010;Ji et al. 2019) or collisionless heating from cosmic ray streaming (Wiener et al. 2013).On the other hand, ionization states observed may be from photoionization (Muzahid et al. 2015).In addition to the heating, CRs can also influence the properties of the multiphase plasma (e.g.Ji et al. 2020;Butsky et al. 2020) via their pressure, potentially affecting these diagnostics even further.
Galactic computational simulations can aid in untangling these degeneracies.CR feedback models have become progressively more detailed (i.e.Werhahn et al. 2021a,b) and interest in tackling these interdependent multiphase and multiphysics questions is growing (Faucher-Giguere & Oh 2023).Several simulations now include three-dimensional radiative transfer at run time while remaining computational feasible (Rathjen et al. 2021;Kim et al. 2022;Katz 2022).At the same time, UV diagnostic data will continue to improve, such as with the James Web space telescope and 30-meter class ground-based telescopes.This data will illuminate CGM spectroscopic properties up to cosmic noon  ∼ 2. In order to take full advantage of the increased quantity and quality of observed data, it is crucial for simulations to make predictions aligned with these diagnostics (e.g.Rathjen et al. 2023).
Our goal is to provide a proof-of-concept of CGM UV absorption-line diagnostic trends derived from three-dimensional MHD simulations, focusing on radiative heating and CR feedback.In Section 2, we describe the MHD simulation domain and feedback implementations.In Section 3 we describe results from updated parallel-plane elongated slab simulations from Holguin et al. (2019) with simplified radiative heating, improved supernova (SN) feedback, and CR collisional losses.In these six simulations, we explore several limiting cases in the parameter space of CR and radiation heating feedback.In Section 4 we describe a parallel-plane post-processing approximation, including the stellar and metagalactic radiation field, that can estimate the ionization state of the plasma throughout the domain.We also examine theoretical expectations for CR and radiation feedback effects on UV diagnostics tracing gas phases between  = 10 4 K and  = 10 6 K.In Sections 5 and 6, we describe resulting synthetic diagnostic predictions and broadly compare them to typical observed values.Finally, in Section 7 we provide concluding remarks.

MHD SIMULATION METHODS
The simulations we run are similar to those from Holguin et al. (2019), with updated stellar feedback methods and additional radiative heating and CR loss models.The simulations are run with the adaptive mesh refinement MHD code FLASH 4.2 (Fryxell et al. 2000;Dubey et al. 2008) using a directionally unsplit staggered mesh solver (Lee & Deane 2009;Lee 2013), including CR transport physics (Yang et al. 2012;Ruszkowski et al. 2017;Farber et al. 2018) in an elongated box geometry with dimensions 2×2×40 kpc 3 (Hill et al. 2012;Walch et al. 2015).
We solve a coupled MHD-CR two-fluid model (Salem & Bryan 2014;Ruszkowski et al. 2017) composed of thermal gas and ultrarelativistic CR fluid with adiabatic indices  = 5/3 and  cr = 4/3 respectively.The CR fluid represents a single energy channel of typical galactic CRs with mean CR Lorentz factor  rel = 3 and momentum distribution slope  = 4.5.Both hadronic and Coulomb CR losses are included (Guo & Oh 2008).CR transport includes advection, dynamical coupling with gas, and diffusing relative to the magnetic field lines with spatially and temporallay constant diffusivities.Heating of the gas occurs due to radiative heating from both non-ionizing and ionizing UV stellar spectrum, CR streaming and hadronic losses, and supernova feedback.Cooling of gas occurs via radiative cooling.Gas dynamics is impacted by self-gravity, stellar momentum feedback, and dynamical coupling with CRs.Star formation occurs according to (Cen & Ostriker 1992), resulting in the creation of a stellar population particle and equal mass loss of gas (see 2.2 for more details).The following equations summarize the model: where  is the gas density,   is the total baryon density including both the gas and stars,  form is the density sink from stellar population particle formation,  *  feed represents the gas density source from stellar feedback (see Section 2.2),   is the gas velocity,  is the magnetic field,  is the gravitational constant,  is the gas gravitational potential,  = −∇ +  NFW is the gravitational acceleration (the sum of gas self-gravity, stellar particle, and halo dark matter contributions to the gravitational acceleration, described in Section 2.1) where  NFW is the gravity from the Navarro-Frenk-White (NFW) potential,  tot is the sum of gas ( th ), magnetic, and CR ( cr ) pressures,  SN is the momentum injection due to stellar winds and SNe.Furthermore,  =  2  +   +  cr +  2 /8 is the total energy density per volume (the sum of gas, CR, and magnetic components, respectively),  is the radiative cooling rate per unit volume,  rad is the radiative heating rate per unit volume, Γ cr is the CR energy loss rate per unit volume, Γ cr,gas is the energy rate gained by the gas from the CR losses,  cr is the CR diffusion coefficient, and  SN is the supernova heating rate per volume.

Gravity
The gravitational acceleration has three components: baryons, stellar particles, and the dark matter halo.Self-gravity from baryons is computed by solving the Poisson equation with the Barnes-Hut tree solver (Barnes & Hut 1986), which is implemend in FLASH by Wunsch et al. (2018).For the dark matter halo gravity modeled as an NFW potential (Navarro et al. 1997), we include the vertical component of gravity as the domain is thin and much smaller than the halo.The halo gravitational acceleration is where  is the gravitational constant,  200 is the halo virial mass,  is the height above the mid-plane,  = ||/ 200 ,  is the halo concentration parameter, and  200 is the virial radius.Table 1 summarizes our choices of these parameters.

Star formation and supernova feedback
Following the approach of Cen & Ostriker (1992) (see also Tasker & Bryan 2006;Bryan et al. 2014;Salem & Bryan 2014;Li et al. 2015), the creation of a stellar population particles occurs when all of the following conditions are simultaneously met: (i) gas number density exceeds 10 cm −3 (Gnedin & Kravtsov 2011;Agertz et al. 2013); (ii) the cell mass exceeds the local Jeans mass; (iii) the divergence of the gas velocity   is negative, ∇ •   < 0; (iv) gas temperature reaches the floor of the cooling function or the cooling time becomes shorter than the dynamical time  dyn = √︁ 3/(32   ).The particle is formed with a velocity equal to the surrounding gas, mass  * =  SF (/ dyn )  3 , where  SF = 0.05 is the star formation efficiency,  is the timestep, and  is the local cell size.A mass equal to the particle mass is removed from the gas.We set a minimum particle mass  * ,min = 10 5  ⊙ in order to manage the number of particles formed, but we still permit particles with  * <  * ,min to form with a probability  * / * ,min and mass  * = 0.8 3 .
In order to model feedback due to stellar winds and SNe, we inject gas thermal energy and momentum, as well as CRs energy.The gas and CR energy is injected into the particle's local cell, while the momentum is injected into the immediately surrounding cells as in Farber et al. (2022). 1 The inclusion of momentum feedback is essential when the resolution does not sufficiently resolve the cooling radius, as the energy injected can be quickly radiated away and fixes such as delayed cooling over predict the energy and momentum at late-times in SN evolution (Martizzi et al. 2015).The full details of the momentum feedback are found in Farber et al. (2022) and we briefly summarize here.We use the following model for SNe evolution in an inhomogeneous medium, found in Eq. 12 of Martizzi et al. (2015), where we assume solar metalicity: The momentum  r () deposition at a particular distance  from center of the cell containing a stellar particle is written as where Θ() is the Heaviside step function,  0 =  * ,mom √︁ 2 SN  SN 0.9 ej,0 .The parameter  * ,mom = 5 represents a boost to the momentum received by the surrounding gas (Agertz et al. 2013;Semenov et al. 2018), to compensate for advection errors, clustering of supernovae (Gentry et al. 2017(Gentry et al. , 2020) ) and backreaction of cosmic rays (Diesing & Caprioli 2018).

Radiative cooling
Radiative cooling is calcuated using the Townsend cooling method (Townsend 2009;Zhu et al. 2017), updated from Farber et al. (2018Farber et al. ( , 2022) ) to include net heating.We approximate the radiative cooling functions from Dalgarno & McCray (1972)   Using these data points, we assign a a luminosity to each stellar population particle, and the total luminosity is used to calculate the radiative heating rates in Eqs.11 and 12.
(1976) with a piecewise power law cooling function Λ() in units of cm 3 s −1 as where  is the gas temperature in K.This cooling function is appropriate for a gas of solar abundance, which is completely ionized at  = 8000 K.

Radiative heating
Radiative heating from stellar sources is calculated using a temporally-varying parallel-plane model, similar to Kim & Ostriker (2017), including from both FUV (6 eV < ℎ < 13.6 eV, far ultraviolet) and ionizing (or EUV, ℎ > 13.6 eV, extreme ultraviolet) spectrum ranges.Photoelectric heating occurs when FUV radiation interacts and frees electrons from interstellar dust grains (Weingartner & Draine 2001).Photoionization heating occurs when ionizing radiation frees electrons from an atom, particularly hydrogen, and the electron eventually recombines (Osterbrock & Ferland 2006).The radiative heating rates in those two ranges  rad,FUV and  rad,EUV , are proportional to the total stellar population particle intensities  FUV and  EUV respectively, at a given simulation time.Due to the parallel-plane approximation, the position of the particles does not matter.We calculate the individual stellar particle intensities using STARBURST99 (Leitherer et al. 1999), which gives a spectrum (erg s −1 Hz −1 M −1 ⊙ ) of a stellar population at a fixed mass  and age  p .We integrate the spectrum in the FUV and EUV ranges to derive the time-dependent luminosity-to-mass ratio Ψ( p ) shown in Figure 1.Since the hydrogen photoionization cross section drops rapidly above 13.6 eV, we consider the EUV range to be 13.6 eV < ℎ < 27.2 eV.
We calculate the total luminosity in each band by summing the luminosity of each particle, given its age and mass from Ψ. Radiative heating occurs only in neutral gas.Kim & Ostriker (2017) models the transition between neutral and ionized gas regimes with a temperature dependent mean molecular weight interpolation, following results from Sutherland & Dopita (1993).We instead use sharp transition between neutral and ionized regimes at a temperature  ion = 1.5 × 10 4 K.This transition temperature is also used to determine the regime of CR losses.We use the photoelectric heating rate Γ pe from Kim & Ostriker (2017), assuming a 4 kpc 2 cross section for our parallel-plane domain. where The photoionization heating rate for hydrogen (Osterbrock & Ferland 2006) is given by where  H is the neutral hydrogen number density, ℎ 0 ≈ 13.6 eV = 1 Ryd is the hydrogen ionization energy ,   is the specific intensity,   is the frequency-dependent ionization cross section of neutral hydrogen.We integrate in the range ℎ = [1 Ryd, 2 Ryd] and assume   is roughly constant over this interval.For a slab geometry with attenuating gas, the mean intensity is given by 4 EUV = Σ EUV (1 −  2 ( ⊥ /2))/ ⊥ = Σ  () (Ostriker et al. 2010), where A(z) is the attenuation factor accounting for absorption of ionizing radiation by neutral hydrogen as a function of height  above the galactic midplane,  2 is the second exponential integral,  ⊥ =  ion Σ()  , Σ n,gas is the neutral column density from the midplane to a given .
All together   ∼  EUV Δ ∼  EUV Δ 4 kpc 2 for our geometry.The heating rate is calculated as follows The total radiative heating is  rad ( H ,  FUV ,  EUV , ) =  rad,FUV +  rad,EUV .The overall net energy recieved by gas is ( rad − Λ), which can take a positive or negative value.
We do not include dynamical coupling between the radiation field and gas.The FUV portion of the stellar spectrum exerts pressure on the dust, which is coupled to the gas.Generally, for winds driven by radiation pressure to be efficient, the stellar flux must be near the dust Eddington limit, which is usually the case in in high redshift and/or star bursting galactic environments (Zhang 2018).In both observed and simulated Milky Way-like galaxies, the focus of our work, other parts of stellar feedback likely dominate the ac-celeration of galactic winds (Andrews & Thompson 2011;Hopkins et al. 2012b).

Cosmic ray physics
CR coupling to gas via advection and diffusive transport follows a similar model as in Farber et al. (2018).We use a constant value  cr for the CR diffusivity parallel to magnetic field lines.2As CRs travel through the ISM, they interact inelastically with particles in the ISM.We include both hadronic and Coulomb losses with a model from Guo & Oh (2008), as implemented in Chan et al. (2019), modeling Γ cr , the energy loss rate of the CRs, and Γ cr,gas , the amount of the CR loss rate that is absorbed by the gas (with the rest of the energy lost as gamma rays): where  n is the number of nucleons,   is the number of free electrons per nucleon, and  cr is the CR energy density.
In the self-confinement model of CR transport (Zweibel 2013(Zweibel , 2017)), CR streaming along the magnetic field lines is subject to the streaming instability, resulting in the generation of Alfvèn waves.The subsequent wave-particle interactions between the population of CRs and Alfvèn waves limits the bulk CR transport speed to the Alfvèn speed.This process also leads to energy transfer from the CRs to the thermal gas via collisionless heating, according to where  A is the Alfvèn speed and  cr is the CR pressure (Wiener et al. 2013).We do not include CR streaming in our simulations, only cosmic ray diffusion.Additionally, it is numerically much cheaper to handle CR transport by diffusion rather than streaming.Moreover, previous work suggest gamma-ray observations are in tension with the slower transport of CR streaming compared to more rapid CR transport by diffusion (Chan et al. 2019;Hopkins et al. 2021c with a higher than typical diffusion coefficient; cf.(Strong et al. 2007)).
Following this idea, we allow for CR streaming heating even though we do not include streaming itself.Using a diffusion model accelerates our simulations because the streaming model imposes more stringent time step constraints than those imposed by typical cosmic ray diffusivities.

Simulation domain and initial conditions
The elongated 2 × 2 × 40 kpc three-dimensional box geometry represents an approximately parallel-plane section of a galactic disk.
The maximum midplane height of ±20 kpc is chosen to ensure a realistic temperature distribution (Hill et al. 2012), while maintaining the parallel-plane approximation.The horizontal  and  coordinate sides have periodic boundary conditions, while the  coordinate 'top' and 'bottom' have diode boundary conditions (zerogradient boundary conditions which set inflowing velocities, when they arise, to zero instead), which avoids issues from gas fallback near the boundary (e.g.Sur et al. 2016).We do not include differential rotation effects as we model a patch at galactic center.We use static mesh refinement with a resolution of 31.25 pc for || < 8 and 62.5 pc for 8 kpc < || < 20 kpc.These higher resolution refinement regions are four times more extended than in Holguin et al. (2019) in order to better capture the phase structure of the plasma farther out of the galaxy, simitar to the choices made in Peeples et al. (2019).
For the gas initial conditions, we use a vertical equilibrium density solution for a stratified, isothermal self-gravitating system (Spitzer Jr 1942;Salem & Bryan 2014) projected to a stratified parallel-plane box (e.g.De Avillez & Breitschwerdt 2007;Walch et al. 2015;Rathjen et al. 2021;Holguin et al. 2019).The unperturbed density distribution is as follows where  0 is the initial midplane density and  0 is the vertical scale height.We can define the initial disk surface gas density as Σ 0 = ∫ 20kpc −20kpc ()d.The density  halo = 1.0 × 10 −28 g/cm 3 is the initial density of the halo.
In order to improve reproducibility, we reduce differences in initial star formation due to machine dependent effects by perturbing the density distribution as follows: where   represents the spatial coordinates and   is the perturbation wavelength in units of the highest resolution of 31.25 pc.We choose the density distribution parameters (see Table 1) such that Σ 0 = 100  ⊙ /pc 2 , which corresponds to the average surface density within a radius of 10 kpc in the isothermal self-gravitating solution.There are no stellar particles initially.After the bulk of star formation occurs (about 100 Myr) in the simulation, roughly 20  ⊙ /pc 2 is converted to stellar particles.The gas initial temperature is  = 10 4 K.The initial magnetic field follows the density distribution as () ∝ () 2/3 with the mid-plane value  0 ≈ 3.The field is oriented along one of the horizontal directions.
From the initial conditions the gas radiatively cools and collapses to a smaller scale height distribution, until the densities become large enough (and the other criteria satisfied as described in Section 2.2) for stellar particles to form and begin stellar feedback.We do not include any artificial pressure support, as resulting steady-state properties are insensitive to particular choices of pressure support (Kim & Ostriker 2017).By 100 Myr, stellar feedback suppresses star formation and the simulations tends towards steady state.

SIMULATION RESULTS
We run 6 simulations to 200 Myr, optimizing for computational expense and exploration of the CR and ionizing radiation parameter space.The simulations are summarized in Table 2.In the short name for the simulations 'k' is followed by either 28 or 29, referring to the log of  cr , 'Hrad' refers to the inclusion of radiative heating in both the FUV and EUV bands, and 'Hcr' refers to the inclusion of cosmic ray streaming heating.All simulations include CRs.In all simulations that include radiative heating there is a reduction in the total EUV luminosity by a factor  * ,EUV = 0.1, which is a reasonable mean escape fraction for ionizing photons from HII regions, although in practice there is a lot of scatter (e.g.Leitherer et al. 1995;Pellegrini et al. 2012;Heckman & Best 2023).Due to our increase in the size of the highest resolution region, the added computational cost of  cr = 10 29 cm 2 s −1 limits the number of simulations we can run at this higher CR diffusivity.
For illustrative purposes, in Figure 2 we highlight both of the simulations with higher CR diffusivity, 'k29_Hrad' and 'k29_Hrad_Hcr', which produce extended outflows.The figure shows slices of a  −  plane through the origin.The color in the slices denotes the value of  i , ,  cr ,  rad , and  cr along the plane.Radiative heating produces a volume filling  ∼ 10 4 K medium around the galactic disk, unlike in our previous simulations in Holguin et al. (2019), where there is a significant amount of colder phases.The scale height of radiative heating is smaller compared to that of CR streaming heating, when it is enabled.The outflow is weaker with streaming heating present, although the CR distribution does extend far outside the galaxy for the weaker wind.This case suggests that the CRs are losing significant amounts of energy to the surrounding gas.
In Figures 3, 4, 5 and 6, we compare the simulation results for all 6 of the simulations we conducted for key quantities.In Figure 4, we refer to columns and rows relative to the label orientation, in that the columns show profile quantities, and the rows are of different simulations.The first column shows the heating rate from radiative and CR streaming.The last two columns show the number density and volume filling fractions of the phases used for post-processing in Eq. 20.First, we describe simulations without CR streaming heating.Radiative heating has the effect of reducing the volume filling fraction of the cold neutral medium (phases 0, 1, and 2) while increasing the fraction of the warm neutral and ionized phases (phases 3 and 4).The hotter phases are not significantly affected, as radiative heating does not affect ionized plasma.The density profiles of the phases remain similar, as the phases are dominated by non-thermal pressure support instead of thermal pressure support, which is shown in the fourth column.The outflows for the lower  cr runs are weaker than in the higher  cr case, in agreement with previous works (e.g.Salem & Bryan 2014;Hopkins et al. 2020;Farcy et al. 2022).The outflows extend to roughly 10 kpc for lower CR diffusivity and 20 kpc for higher CR diffusivity.With higher CR diffusivity, the density profiles of the gas phases within 10 kpc are not so tightly overlapping, even though the non-thermal pressure support is dominant.The overlapping of phase densities does occur past 10 kpc.The dominant phases in the outflow are the hot phases (phases 6 and 7), unlike the lower  cr case.There is also a resurgence in star formation for higher  cr , likely due to CRs more easily leaving the galactic disk and reducing pressure support after their initial injection by the first burst of star formation.
Simulations with CR streaming heating significantly reduce the spatial extent of outflows from the galaxy for both values of  cr .The energy transfer to the gas may not have a large influence in its thermal energy for some parts of the cooling curve (see Hopkins et al. 2020), but the energy loss can be significant to the CR population.In a CGM model with a fixed CR flux Huang & Davis (2022) similarly found that CR streaming heating can be a dominant mechanism of energy transfer from the CRs.Two quantities point to the outflow weakness being caused by the reduction in CR energy density:  cr /Γ cr (the CR energy loss rate relative to hadronic and Coulomb processes in Eq. 18) and the loss timescale  cr,st ∼  cr / cr (Eq. 19).CR collisional losses are proportional to gas density, so as the density decreases farther away from the galaxy, these losses are reduced.Setting the collisional and collisionless CR loss terms equal to each other yields a critical number density, where  cr is the length scale of the global cosmic ray distribution and   is the magnetic field strength in .Below this density, collisionless losses dominate over collisional losses.For reasonable parameters in the low density medium above the galaxy, CR streaming losses should be the dominant loss mechanism below  ∼ 0.1 cm −3 .We also calculate the timescale of CR energy loss via CR streaming heating, which indicates that for reasonable parameters, the CR energy loss timescale due to collisionless processes is of order the total simulation time of 200 Myr, or shorter in lower density regions.The column for  cr / th in Figure 4 shows that in the simulations with CR streaming heating, CRs retain pressure dominance in the disk (within 3 kpc) but farther out of the galaxy, thermal pressure domi-nates, as expected by our analysis above.In the heating rate column, we note that CR streaming heating does not match the radiative heating rate until roughly 7-10 kpc in simulations with both heating mechanisms included.Although, it is likely we could be overestimating the radiative heating rate farther away from the galaxy since a non-attenuated parallel-plane radiation field does not decay with height.In global simulations, the radiative heating rate could fall faster with height above the disk (as the stellar field decays), allowing other heating mechanisms like CR streaming heating to dominate for a larger range in height.
Figure 5 provides another perspective of the radiative heating and pressure state of the galaxy.This figure highlights the varied combinations of radiative heating and pressure support states of the gas depending on the feedback mechanisms included.In the simulations without collisionless CR losses, the CR-dominated pressure remains similar across all radiative heating values.Comparing the 'k28' and 'k29' simulations, we see a loss of CR pressure support in the midplane with faster CR diffusivity, as the CRs more quickly diffuse out of the disk.The slabs approach the CR-dominated values of pressure at roughly 2 kpc in height.In the 'Hcr' cases with collisionless CR losses, the CR pressure support declines with height in between 2-5 kpc.Interestingly, the 'k28_Hrad_Hcr' case contains gas that exists at all levels of CR pressure support between the highest level and thermally dominated.Figure A3 shows the same plot, except with the  cr contribution to heating included.The main difference in this plot is that the heating rate for data points farther away from the disk, with lower  cr and  rad , is larger so that the curve decreases smoothly with the heating rate.Further discussion of predicted plasma diagnostics associated with these plots occurs in later sections.Figure 6.SFR of the 6 simulations we ran.There is an initial burst of star formation at 60 Myr.The simulations with higher  cr experience another, larger burst of star formation at 120 Myr, due to CRs leaving the dense disk more efficiently.This removal of pressure support triggers more gravitational collapse and star formation.

Division of simulation domain and Cloudy model setup
In order to produce synthetic plasma diagnostics from our simulations, we use the spectral synthesis code, Cloudy (Ferland et al. 2017).A typical Cloudy parallel-plane model takes inputs of hydrogen density, plasma temperature, metalicity, and incoming spectra (radio to x-ray), resulting in outputs of species ionization state and outgoing spectra.We use the MHD simulations to inform various Cloudy models.Cloudy is able to calculate the plasma temperature.However, this calculation does not include the various feedback mechanisms present in the three dimensional MHD simulation, so we let the model temperature be informed by the simulation.We calculate the galactic midplane spectra  , * by assigning spectra to each simulation stellar population particle based on particle age and mass using Starburst99.
Figure 7 shows a diagram that illustrates the post-processing framework.We divide the simulation domain into parallel-plane sections of two resolution cells in height.Within each section , we calculate the mean and standard deviation of  and  of each gas phase.We define eight gas phases in each section  from the temperature floor to hot coronal gas as follows, Phase 0 : Phase 3 : 8000K < T < 1.5 × 10 4 K Phase 4 : 1.5 × 10 4 K < T < 3 × 10 4 K Phase 5 : 3 × 10 4 K < T < 10 5 K Phase 6 : 10 5 K < T < 10 6 K Phase 7 : 10 6 K < T (20) which include rough matches to the traditional ISM phases (phase 3: warm, neutral medium; phase 4: warm, ionized medium) (Draine 2010).Each phase  has mean density ⟨ ,  ⟩, temperature ⟨ ,  ⟩, and associated standard deviations  , and   , respectively, as well as the filling fraction  vol;i,j of each phase in that section.A Cloudy model is assigned density  C ,  and temperature   ,  by drawing from a Gaussian distribution defined by the mean and standard deviation of density and temperature of each phase, with values limited to be within one standard deviation.We draw properties 10 times in order to produce 10 separate profiles of density and temperature in the domain.These separate profiles allow us to produce intervals in the synthetic line ratio profiles.The incoming spectra for each model is the midplane spectra  ; * if  = 0 (lowest height section) or  ;−1 (the outgoing spectra from section below the current one) otherwise.The output spectra of the Cloudy model is  ;,  .The total outgoing spectra  , for section  is the sum of the eight spectra in each phase  ;,  weighted by the phase filling fraction  vol;i,j , which becomes the input spectra to the Cloudy models in the next section  + 1.Through this framework, we obtain estimates of the ionization state of numerous species and are then able to generate plasma diagnostics that are spatially (with midplane height) and temporally dependent.
The parameter  * is a binary value that activates or suppresses the presence of radiative heating.Due to limited resolution, we do not resolve the giant star forming complexes where massive stars form.We multiply the ionizing part of the stellar spectrum by an escape fraction  * ,EUV to account for the absorption of hydrogen ionizing photons after emission from stars and out into the ISM.We also have  bkg,EUV , the fraction of the metagalactic background field present. bkg,EUV is a binary on/off option.Since the postprocessing models the outward and not inward propagation of radiation, the metagalactic background field must be included in each Cloudy model as a binary on/off option.The metagalactic field is added to a Cloudy model by specifying the desired field name.We use 'HM12' (Haardt & Madau 2012) for our work.In order to avoid adding the previous slab attenuated background field every time we transmit a spectra to a new domain section, we add the 'no isotropic' option, which automatically removes the attenuated component of the isotropic spectra sources.Overall, we have consistency between the the choices in post-processing parameters and simulation parameters (e.g., same  * ,  * ,EUV and  bkg,EUV in simulation and post-processing).

Observational absorption-line diagnostics
Absorption-line ratios are common tracers of the gas kinematic and plasma state.Their values can distinguish between various multiphase non-equilibrium models (e.g.radiative turbulent mixing layers (Ji et al. 2019), conductive interfaces (Borkowski et al. 1990;Marcolini et al. 2005;Brüggen & Scannapieco 2016).Line ratio diagnostics are convenient to use because our parallel-plane domain represents only a piece of the galaxy, preventing an accurate prediction in the total ion column densities.We focus on three relatively common diagnostic absorption-line ratios: two intermediate (ionization energy > 30 eV) ion ratios, Si-IV/C-IV and N-V/O-VI, as well as one low-to-intermediate ion ratio, C-II/C-IV (Tumlinson et al. 2017).These diagnostics span temperature ranges from the warm, neutral phases (phases 3,4) to just below the hot phase (phase 5).
In order to explore the temperature and density parameter space that the simulation post-processing will effectively draw from, we run a grid of 20 by 20 data points of Cloudy models in  i = [10 −3 , 10 1 ] cm −3 and  = [10 2 , 10 6 ] K, evenly spaced in log space, similar to common analyses done in observationallyfocused work (e.g., Fox et al. 2005).This analysis allows us to gain intuition for the potential interaction between radiative and CR feedback mechanism.The Cloudy models all assume the same domain thickness of 50 pc, roughly the size of two of the highest resolution elements in the MHD simulations.The models in this figure are are not informed by the simulation results.The density and temperature values in this analysis do not represent upper or lower limits for the Cloudy models in subsequent sections that are informed by simulation results.The input radiation field is that of the SFR = 1  ⊙ /yr stellar population spectra from Starburst99, multiplied by an escape fraction  * ,EUV = 0.1.The results are shown in Figure 8 for Si-IV/C-IV, N-V/O-VI, and C-II/C-IV from left to right.The black line contours correspond to the line ratio values without any external radiation field.The red line contours show the line ratio values with the stellar field included.The background filled color corresponds to the transmission fraction of the incoming spectra through the slab.This color is plotted to illustrate the optical thickness of the medium to the spectra.
The relatively constant value of Si-IV/C-IV ∼ 0.2 observed in the Milky Way (Savage & Wakker 2009;Werk et al. 2019) is a useful guide for examining the line ratio parameter space.In collisionally ionized gas, the possible gas phases for that line ratio are restricted to a narrow temperature range around  ∼ 10 5 K: the value of the line ratio is also quite sensitive to the temperature around  ∼ 10 5 K, with two orders of magnitude changes in value.With a SFR = 1 M ⊙ /yr stellar field added, the line ratio contours are quite different at  < 10 5.5 K. Lower temperature models in the range  = [10 3 , 10 4 ] K are able to produce Si-IV/C-IV ∼ 0.2, but only at low densities of  ∼ 10 −1.5 cm −3 .Under standard ISM thermal pressure balance (McKee & Ostriker 1977;Draine 2010), these phases of gas are required to be much denser (at lease an order of magnitude) to remain in pressure equilibrium with the hot gas phase.By providing a non-thermal pressure source, CRs can allow multiple phases to exist at the same density (see Ji et al. 2020), in particular, by allowing more tenuous cold phases.The overplotted blue circles roughly mark where (right) thermal and (left) non-thermal gas phases exist.Whereas in the thermally supported case, lower temperature phases cannot reproduce observed line ratio values even with an external radiation field, with CRs, a variety of phases from  = 10 3 K to  = 10 5 K could reproduce the observed line ratio values.
The parameter space for C-II/C-IV is similar to that of Si-IV/C-IV.The addition of a stellar radiation field allows phases with  < 10 5 K to produce lower line ratio values.With non-thermal pressure support, these tenuous, cold phases can plausibly exist and contribute to the observed signal.The N-V/O-VI line ratio is more resistant to change by the external stellar field, as expected by the higher ionization energy of the ions.Typical observed values of the line ratio are in the range [10 −1.5 , 10 0.5 ] ( Wakker et al. 2012), which are found in a small part of the parameter space in Figure 8, favoring a collisionally ionized medium.Some line values of roughly 10 1 can occur at ∼ 10 4 K with a stellar radiation field present, but only in very tenuous gas in this phase with  < 10 −2 cm −3 .We additionally tested (not shown here) the effect of only the metagalatic background field and found that it produces effects in between that of the collisional ionization and stellar field cases.The metagalactic background field itself is not sufficient to change the parameter space significantly.
Besides pressure support, CR streaming heating could have an impact on the gas thermal state.Radiative heating does not occur in a significantly ionized gas, so the radiation field does not heat this gas (although radiative effects on the line emission/absorption still can occur).Wiener et al. (2013) suggested that because the CR scale height extends well outside the galactic disk and collisionless CR streaming heating affects the ionized medium, CRs could provide a plausible heating mechanism to explain anomalously high electron densities outside the disk (Reynolds et al. 1999).Regardless of the pressure equilibrium (thermal or non-thermal), additional energy injection into the ISM/CGM plasma reduces the value of the Si-IV/C-IV and C-II/C-IV line ratios in Figure 8.  .Line ratio diagnostics from Cloudy models in density and temperature space for a 50 pc slab (∼ 2 high resolution cell widths in our MHD simulations).These models are not informed by MHD simulations.The black dotted contours show the line ratio values assuming collisional ionization, while the red contours show the values with the addition of photoionization from a constant SFR stellar population spectrum (SFR = 1 M ⊙ /yr).The background color shows the transmission fraction through the slab for reference.The addition of the stellar field significantly changes the line ratio contours at  < 10 5 K, except for N-V/O-VI.We sketch the potential effect of non-thermal pressure support on the Si-IV/C-IV plot.The right pair of blue dots are placed roughly where the traditional warm, ionized medium and cold, neutral phases exist based on thermal equilibrium.When the dominant pressure source is non-thermal pressure, multiple phases can exist at the same density (e.g.Ji et al. 2020).When these phases are shifted to lower densities, the line ratio values are closer to Milky Way values.In a realistic system the dominant pressure support will vary between thermal and non-thermal, so the analysis here provides limiting cases.

POST-PROCESSING RESULTS
Section 4.2 demonstrates influence of the gas thermal state and local radiative field on values of C-II/C-IV, Si-IV/C-IV and N-V/O-VI diagnostics.In a realistic system, many physical processes are acting in a complex, interdependent fashion, making the inverse problem of determining the gas state from spectral observations particularly challenging.We further explore this parameter space using the three-dimensional galaxy simulations from Section 3, which capture key parts of the nonlinear galactic system.The simulated CGM diagnostics can help disintagle the importance of various physical processes through testable predictions.We post-process our simulations using the method described in Section 4.1.

BPT diagram
Figure 9 shows where each simulation falls on the BPT diagram (Baldwin et al. 1981), color coded for different distances from the midplane.Each axis on the plot is a line ratio in emission.The BPT diagram is commonly used to divide observed galaxies into those whose HII regions are primarily ionized by stars and those primarily ionized by an active-galactic nuclei (AGN).In other words, the dividing line on the BPT diagram separates galaxies into those ionized by a softer or harder spectrum respectively (Veilleux & Osterbrock 1987).
The BPT diagrams are useful in that they help summarize the properties of the ionizing radiation field, as well as ISM plasma conditions and metalicity (Kewley et al. 2013b).The data points from 'k28' are below the typical star forming sequence, as expected by the lack of stellar radiation present in the run.The addition of only CR streaming heating in 'k28_Hcr' moves the data points vertically in O-III/H, overlapping star forming sequence.The addition of stellar radiation for low CR transport speed tends to move the near-disk data points horizontally and narrows values in N-II/H.For 'k28_Hrad_Hcr' that includes both radiative and CR streaming heating, we see the combined effects on the BPT diagram from adding each heating mechanism separately.The near disk data points in the higher CR transport simulation ('k29_Hrad') with only radiative heating are roughly around the star forming sequence, at a similar N-II/H and higher O-III/H compared to the lower CR transport version.The addition of CR streaming heating at higher CR transport moves near-disk data points off of the star forming sequence and closer to the dividing region between soft versus hard ionization.

Other plasma diagnostics
Figures 10, 11, and 12 show the time evolution of the resulting line ratio profiles in height for Si-IV/C-IV, N-V/O-VI, and C-II/C-IV respectively.For Si-IV/C-IV, we plot the observed Milky Way value of ∼ 0.2.For N-V/O-VI, we plot lines bounding the range of observed values.For C-II/C-IV, we plot the mean value from Fox et al. (2005).Figure 5 shows the mean radiative heating and  cr / th , with the color representing the Si-IV/C-IV value, within the inner 10 kpc slabs of the galaxy.This plot illustrates the differences in pressure balance and radiative heating between the simulations, and what diagnostic values occur in that parameter space.Similar figures for N-V/O-IV and C-II/C-IV do not yield as much information (especially N-V/O-IV due to relatively constant profiles) on the line ratio behavior, so those figures are placed in Appendix A.
Most of the simulations over predict Si-IV/C-IV, and the largest values tend to occur in  cr -dominated areas.The simulation 'k28' predicts Si-IV/C-IV ∼ 10 2 , more than two orders of magnitude above the observed values in the Milky Way (Savage & Wakker 2009;Werk et al. 2019).Given the lack of appropriate physical processes included in this simulation, its inability to reflect observational data is not unexpected.The addition of radiative effects in 'k28_Hrad' reduces the line ratio by a factor of 10 across the entire profile.The profile is also less variable in space and time than the profile from the 'k28' run.Even though the plasma pressure is entirely dominated by CRs on average at large  rad , we do not see a sufficient reduction in the line ratio towards the observed values , as Figure 8 would predict.If we instead add CR streaming heating, as in the 'k28_Hcr' run, we see that the line ratio values are even lower, including a profile at 130 Myr that is below the observed value.The outflows in this run are more suppressed, as mentioned in Section 3, so this run does not entirely match the observed data that extends up to 15 kpc.This simulation better isolates potential CR streaming heating behavior, as the radiative heating is not present.Interestingly, Figure A3 shows two regions of low line ratio values: the highest  cr ∼ 10 −25 [ergs −1 cm −3 ] data points, and then the data points with a slow decrease in heating at  cr ∼ 10 −27 [ergs −1 cm −3 ].The latter region corresponds (by comparing to the  cr / th profile) the flat line ratio profile between 2 and 8 kpc.In this region, the second most abundant phase (behind the hottest phase) is Phase 6, the intermediate gas phase.At higher  cr the impact of CR streaming heating is also a suppressed wind, favoring 'k29_Hrad' over 'k29_Hrad_Hcr'.Both runs produce Si-IV/C-IV values a factor of ∼ 5 − 10 larger than the observed value.Even though 'k29_Hrad' and 'k29_Hrad_Hcr' have different heating and pressure behavior, as Figure 5 shows, the actual line ratio profiles are both relatively flat.Furthermore, for all simulations with CR streaming heating, significant amounts of gas exist in interval 1 <  cr / th < 3, but we do not see the expected trend of higher CR pressure and lower line ratio values.
The N-V/O-VI line ratio provides another metric with which to evaluate the simulations.Values of this line ratio could favor particular underlying physical processes, such as cooling flows, shock heating, or turbulent mixing.In this work, we show as green solid lines the range of values measured observationally, as summarized by Wakker et al. (2012), in order to broadly evaluate each model.The N-V and O-VI ions have a larger ionization potential than Si-IV and C-IV ions, so they trace gas phases well-above those directly influenced by radiative heating (i.e.see Figure 8).Generally the simulations favored for N-V/O-VI are the same ones favored for Si-IV/C-IV.The simulation that least matches the observed bounds is the 'k28' run, which has N-V/O-VI several orders of magnitude above the bounds.This behavior can be traced to the lack of  ∼ 10 5−6 (phase 5) gas.O-VI sharply peaks within phase 5 temperature range, so a severe lack of phase 5 gas produces large line ratios.As expected by the higher ionization potentials, adding radiative heating as in 'k28_Hrad' does not significantly improve results.The runs with CR streaming heating do have predicted line values near observed bounds.At higher  cr , both runs with and without CR streaming heating produced reasonable N-IV/O-VI values, corresponding to phase 5 gas having a significant volume filling factor in both simulations.
C-II/C-IV traces lower temperature gas phases than the previous line ratios discussed.Similar to S-IV/C-IV, the values tend to be larger for larger  cr , as Figure A2 shows.As before, the 'k28' over predicts the line ratio value (i.e.there is too much colder, neutral material).Adding CR streaming heating, results in profiles that are much closer to the example observations.The simulation 'k28_Hcr' (without radiative heating) in particular, produces values  2).Generally, the regions closer to the galactic disk (< 5 kpc) in all of the simulations were roughly on the star forming sequence, except for 'k28', which did not include radiative or CR streaming heating.

K29_Hrad
K28_Hrad K28 K28_Hcr K29_Hrad_Hcr Figure 10.Line ratio profiles of Si-IV/C-IV for 6 simulations, colored by simulation time.The dashed, green line denotes the observed value in the Milky Way (Werk et al. 2019).For high CR diffusivity, 'k29_Hrad' (not including CR streaming heating) had the best match to Milky Way value.For low CR diffusivity, the addition of CR streaming heating reduced the line ratio values towards observed values, but the gas outflow was weak.
at 190 Myr that are quite close to the example observations, and is the only simulation at lower CR diffusivity that has an under predicting profile (120 Myr).This behavior could be due to the fact that this simulation uniquely has an outflow dominated by phase 3 gas (warm, neutral medium) as shown by the filling fraction in Figure 4. Of the simulations with higher CR diffusivity, the run 'k29_Hrad_Hcr' produces the best match to the observed line ratio values we chose.

Caveats and further work
In the following list we describe potential caveats in our work and additional research directions: • The observational values we use in Figures 10, 11, and 12 are meant as rough guides to aid simulation comparison.The synthetic line ratio data points from the simulations all most certainly trace a larger range of physical situations than would be traced by a particular observation (e.g.high velocity clouds).
• The values of synthetic UV diagnostics depend on SFR peaks that may be fluctuations.More work is needed to see how relevant these peaks are to changes in model assumptions.
• We do not focus on interstellar dust effects on the predicted line ratios.Dust preferentially reduces certain elements compared to others.Preferential depletion of Si onto dust will lower the Si-IV/C-IV ratio allowing some previously observationally inconsistent models to be viable, as Werk et al. (2019) notes.The exact depletion of elements relative to others is still debated.More detailed models of dust creation/destruction (e.g.Öberg 2016) each can influence line diagnostics based on their own resulting elemental depletion.
• CR transport is an active area of research.More detailed work in plasma turbulence and charged particle transport, especially at spatial scales below which the fluid approximation for CRs holds, are needed in order to better understand the impact of CRs on global galaxy scales (e.g., Bai 2022).
• We use a single cooling function assuming solar metalicity.We do not explicitly track ion species with the code, and we do not include non-equilibrium effects in the cooling.
• Our resolution is relatively limited compared to the scale of many important processes operating at the CGM cloud pc-scales.For example, we do not resolve turbulent mixing layers, cloud formation/destruction, or CR-cloud interaction.
• Our model of parallel-plane radiation propagation through the domain, especially with un-resolved CGM cloud structure, is a crude approximation of real ionizing photon transport through the galaxy and halo.This work is meant to be a proof-of-concept highlighting potential UV plasma diagnostics to use as metrics to evaluate models of cosmic ray and ionizing radiation within galactic simulations.
• In addition to our simulation domain being an approximation of a global galactic disk and halo, it is also not cosmological.Cosmological processes such as mergers and gas accretion are not present in our simulations.

DISCUSSION
The BPT diagrams in Figure 9 showed that the inner 5 kpc of our galaxy models (most comparable to the observed lines from the ISM of galaxies (e.g., Kewley et al. 2013a)) generally follow the star forming sequence, except 'k28' which explicitly excluded the stellar and metagalatic field.For some simulations, especially 'k28_Hrad', the regions farther out of the galaxy tend to cross the line between star forming and active galaxies.This is an indication that the spectra in these simulations are becoming harder with height, possibly due to higher photon absorption in the lower CGM.The fact that Figure 4 shows the dominant volume filling phase in 'k28_Hrad' is phase 3 (warm, neutral medium) supports this idea.
The source of the substantial reservoir of O-VI-bearing gas in the CGM, coupled with the observed N-V/O-VI and Si-IV/C-IV ratios, has been debated, especially as it related to underlying photoionzed vs. collisionally ionized gas phases.A long-lived intermediate temperature ( ∼ 10 5 − 10 6 K) could explain these observations.In addition to several other non-equilibrium processes, CR streaming heating could also help generate this gas phase.All of our simulations with CR streaming heating produce a significant intermediate temperature volume filling phase (phases 5 and 6) in outflow regions.In fact, those phases are the dominant phase besides the hot medium (phase 7).In the simulation without radiative heating and only CR streaming heating ('k28_Hcr'), the region approaching the observational values of Si-IV/C-IV corresponds to regions where phase 6 was significantly volume filling.The high CR diffusivity ('k29_Hrad') simulation did also have significant intermediate temperature gas.The exact origin of the intermediate temperature phase is unclear: is the additional thermal energy from CR collisionless heating responsible the phase, or is it that the energy loss from the CRs causes an effect (i.e.different star formation history) that eventually produces the phase?Higher CR diffusivities do lead to more stars produced, and thus more stellar feedback.
The role of photoionization is less important than expected from Figure 8.We hypothesize that the CR-dominated CGM pressure equilibrium would allow lower temperature phases to exist at low densities and naturally give rise to lower value for Si-IV/C-IV line ratio.In fact, the peak line ratio values in the simulations tend to occur within CR-dominated region.The two galaxies ('k28' and 'k28_Hcr' ) with strong density overlap of the gas phases had the worst match to observed line ratio values.Photoionization did contribute during the SFR burst in 'k29_Hrad', where contribution to the Si-IV and C-IV lines tended to come from phases 3 and 4, as well as 5 in some samples.The SFR burst itself may be due to a particular feedback choice such as star formation efficiency, but regardless, the 'k29_Hrad' model is an example of how both photoionization and collisionalizion play key roles in the plasma diagnostic values.Simulations with a fixed SFR would provide additional insight into the relative importance of direct vs. indirect effects of radiative and CR feedback.
The effect of CR streaming heating on the strength of the outflow is also important.For choices of CR diffusion, streaming heating removes enough energy from the CRs such that a weaker outflow occurs, remaining below 10 kpc.At low CR diffusivity, the galaxies with collisionless heating included produced more significant quantities of intermediate temperature phases, but with the small outflow extent.The main difference between galaxies with higher CR diffusivity and those with lower CR diffusivity is the extent of the galactic outflow.One interpretation of this result is that CR confinement occurs not via self-confinement, but instead is due to extrinsic turbulence, which is not associated with collisionless plasma heating.There are many theoretical arguments regarding self-confinement versus extrinsic turbulence models that are outside the scope of this work (see for example Yan & Lazarian 2002;Blasi et al. 2012;Zweibel 2013;Xu 2021;Kempski & Quataert 2021).
If CR streaming heating is indeed present in the real galactic environment and assuming there are not additional wind launching mechanisms, for high  cr , there must be a mechanism for streaming heating suppression that our model does not take into account.Our resolution scale is above that of expected pc-scale clouds of dense material in the CGM (Gronke & Oh 2018;Farber & Gronke 2022), so we do not model phenomena at these small scales.According to the timescale of CR losses to collisionless heating in Eq. 19, the loss timescale is proportional to  cr , the length scale of the cosmic ray distribution, which locally can change significantly at the cloud interface.Proper accounting of energy losses due to CR streaming heating requires accurate modeling of cloud-scale physics.The path of the CR distribution through a turbulent, clumpy medium has been an active area of interest (e.g., Wiener et al. 2017Wiener et al. , 2019;;Brüggen & Scannapieco 2020;Huang & Davis 2022).For example, Bustard & Zweibel (2021) focused on the propagation of CRs through a clumpy medium, accounting for effects of cloud partial ionization.Cloud ionization can change where the CR gradient forms at cloud interfaces, influencing where significant collisionless heating occurs.Also, partial ionization, with associated ion-neutral transport enhancement of CRs, allows CRs to more easily pass through clouds.Of course, in a more realistic model, the cloud ionization should be influenced by the stellar and metagalactic radiation field.Other wind launching mechanisms besides CRs may allow extended outflows with phase structures similar to 'k29_Hrad_Hcr'.
The dependence of results on CR diffusivity also demonstrates the necessity of more detailed CR transport models in order to better account for the impact of feedback effects on UV diagnostics.The improvement in the values of diagnostics at low CR diffusivity could indicate that CR streaming heating indeed improves line ratios, but only accompanied by a corresponding weaker wind.In a realistic system, the CR diffusivity could vary greatly on pc-scales as the cosmic rays traverse CGM inhomogeneities: for example, Bustard & Zweibel (2021) showed that the formation of CR bottlenecks depends on cloud ionization.The pressure anisotropy instability (Zweibel 2020) may contribute to CR confinement in these environments.Other models of CR transport will similarly cause a cascade of coupled effects.For example, in the CR transport model studied in (Holguin et al. 2019), CR transport is enhanced within roughly a kpc of the midplane due to turbulence, but is much less enhanced farther out of the disk where turbulence is weaker.This type of model causes enhanced star formation (and subsequent outflows powered by the stellar feedback) because CRs can more easily escape the dense disk.This model also reduces CR losses from dense gas, which leaves a greater reservoir of CR energy outside of the disk.Recent work has also investigated the potential significance of dust damping of Alfvén waves to CR confinement (Squire et al. 2021).The properties of the galactic radiation field may impact the subsequent coupling of dust to the CR system.

CONCLUSIONS
The scope of this work is to provide a proof-of-concept study of the interaction between two complex galactic feedback mechanisms, ionizing radiation and CRs, and the resulting effect on UV diagnostics.We sketch out how such an analysis may be conducted and highlight where complex feedback interactions may or may not occur.More detailed modeling of feedback mechanisms, galactic outflows/inflows, and methods for calculating synthetic diagnostics will reveal a wealth of information about the galactic and circumgalactic environment.
We perform MHD simulations of a section of a galactic disk, including radiative heating and CR feedback.We test the effects of omitting radiative heating, the value of CR transport parallel to the magnetic field ( cr = 10 28 or 10 29 cm 2 s −1 ), and the inclusion of collisionless CR streaming heating.We post-process the simulations using parallel-plane Cloudy models informed by the simulation results in order to produce synthetic observations of four UV diagnostics: the BPT diagram, Si-IV/C-IV, N-V/O-VI, and C-II/C-IV.We compare the gas phase structure, predicted UV diagnostics compared to observational values, and evaluate the relative ability of the simulations to reproduce observations.Our conclusions are: • We demonstrate that including radiative heating in simulations can have a substantial impact on observable Si-IV/C-IV values, while unable to improve N-V/O-VI values (due to the higher ionization potential).Although, this effect suffers some interdependence among other model parameters that is challenging to quantify.Nonetheless, variations of other model parameters, like star formation history, are likely playing a sub-dominant role.Other feedback mechanisms may be needed to better simultaneously match a variety of UV diagnostics.
• The inclusion of CR streaming heating in our simulations reduces the strength of galactic outflows due to the additional energy losses of the CR population.
• In simulations with slower CR transport, the addition of CR streaming heating results in synthetic diagnostics values closer to observed values due to the production of intermediate, transition ( ∼ 10 5−6 ) temperature gas.This finding lends credibility that this intermediate phase can exist within a CR-driven outflow.However, slower CR transport results in a less extended outflow.
• The simulations with the inclusion of faster CR transport produce the closest predictions to observed UV diagnostics due to the production of transitional temperature gas.Without CR streaming heating losses included in feedback, the outflows extend out to the edge of the domain at a height of 20 kpc above the midplane.
• The type of analysis presented in this proof-of-concept also demonstrates the potential use of UV diagnostics as additional constraints on key galactic physical models that are currently not well understood, such as cosmic ray transport and other feedback mechanisms.

Figure 1 .
Figure 1.Luminosity-to-mass ratio ( ⊙ per  ⊙ ) from Starburst99 in the FUV and EUV wavelength ranges as a function of stellar population age.Using these data points, we assign a a luminosity to each stellar population particle, and the total luminosity is used to calculate the radiative heating rates in Eqs.11 and 12.

Figure 2 .
Figure 2. Slice plots from simulations at 190 Myr, both with  cr = 10 29 cm 2 s −1 and radiative heating.The slice are (left to right) number density   , temperature , cosmic ray number density  cr , and heating rates (radiative  rad and CR streaming heating rate  cr ).The outflow in 'k29_Hrad' extends to the top of the domain at 20 kpc from the midplane.The inclusion of CR streaming heating losses reduces the CR population enough to reduce the extent of the outflow to roughly 10 kpc.

Figure 4 .]Figure 5 .
Figure 4. Profiles at 190 Myr from 6 simulations for theating rate from total radiative and CR streaming separately, , (left), gas densities per thermal phase (middle), and volume filling fraction  vol (right).The thermal phases, defined in Eq. 20, correspond to the phases used in the Cloudy modeling.

Figure 7 .
Figure 7. Diagram of MHD simulation and Cloudy post-processing framework.A slice (2 × 2 × 5 kpc) of an example simulation is shown divided up into red parallel plane sections, the properties of which are used to inform Cloudy models.Each section holds eight temperature phases, from Eq. 20, each with mean density , temperature  and volume filling factor  vol .The midplane spectra from the stellar population particles is transmitted outward through the sections.We calculate the ionization state of H, C, N, O, and Si ions in each section.

Figure 8
Figure8.Line ratio diagnostics from Cloudy models in density and temperature space for a 50 pc slab (∼ 2 high resolution cell widths in our MHD simulations).These models are not informed by MHD simulations.The black dotted contours show the line ratio values assuming collisional ionization, while the red contours show the values with the addition of photoionization from a constant SFR stellar population spectrum (SFR = 1 M ⊙ /yr).The background color shows the transmission fraction through the slab for reference.The addition of the stellar field significantly changes the line ratio contours at  < 10 5 K, except for N-V/O-VI.We sketch the potential effect of non-thermal pressure support on the Si-IV/C-IV plot.The right pair of blue dots are placed roughly where the traditional warm, ionized medium and cold, neutral phases exist based on thermal equilibrium.When the dominant pressure source is non-thermal pressure, multiple phases can exist at the same density (e.g.Ji et al. 2020).When these phases are shifted to lower densities, the line ratio values are closer to Milky Way values.In a realistic system the dominant pressure support will vary between thermal and non-thermal, so the analysis here provides limiting cases.

Figure 9 .
Figure 9. BPT diagram of 6 simulations.The solid line denotes the general line between star forming and active galaxies fitted from the Kewley et al. (2001) fit, while the dashed line denotes the rough trend for star forming galaxies.The simulations are shown at 190 Myr (same time shown in Figure2).Generally, the regions closer to the galactic disk (< 5 kpc) in all of the simulations were roughly on the star forming sequence, except for 'k28', which did not include radiative or CR streaming heating.
and Raymond et al.

Table 1 .
Simulation model parameters

Table 2 .
Summary of the 6 simulations we conducted, spanning limiting cases in CR diffusivity, as well as radiative and CR streaming heating.Sim Number Sim short name Radiative heating log  cr [cm 2 s −1 ] Profiles at 190 Myr from 6 simulations of gas and CR number densities, as well as ratio of CR to thermal pressure  cr / th , vs. midplane height.The 'k29' simulations with larger CR diffusivity, have a more extended outflow of gas compared to their 'k28' counterparts.The inclusion of CR streaming heating suppressed the outflow.The lower lower CR diffusivity simulation with CR streaming heating and radiative heating, resulted in the weakest outflows.