Introduction

The circumgalactic medium (CGM), a gaseous halo enveloping galactic disks, acts as both a reservoir for star formation and a repository for galactic outflows and feedback processes. Its multiphase structure, characterized by regions of varying temperatures and densities, retains the imprints of the galaxy's formation and evolutionary history. Investigating the CGM is crucial for understanding the baryon cycle in galaxies, as it regulates the flow of gas between the galaxy and the intergalactic medium, thereby influencing star formation rates and the effectiveness of feedback processes (J. Tumlinson et al. 2017;C.-A. Faucher-Giguère & S. P. Oh 2023).

Recent observational and theoretical studies have significantly advanced our understanding of the CGM. Observations using metal absorption lines, such as those from O VI, C IV, and Mg II, have unveiled the complex multiphase nature of the CGM, highlighting substantial variations in temperature, density, and ionization states across different regions (J. X. Prochaska et al. 2017;J. Tumlinson et al. 2017). Discrepancies between observed line ratios and model predictions indicate the necessity for models beyond collisional ionization equilibrium (CIE) and photoionization equilibrium (PIE) to accurately interpret CGM properties. In particular, the observed line ratios of different ions often deviate from predictions of equilibrium ionization models, suggesting the presence of additional physical processes. B. P. Wakker et al. (2012) analyzed the Galactic CGM and found that equilibrium models alone fail to fully explain the observed O VI, C IV, N V, and Si IV column densities. Their results suggest that nonequilibrium ionization processes, such as cooling flows and conductive interfaces, may be necessary to sustain the observed transition-temperature gas. Similarly, J. K. Werk et al. (2016) studied the CGM of low-redshift starforming galaxies and found that standard photoionization and CIE models struggle to fully explain the observed distribution of O VI. Their results indicate that additional ionization sources may be required.

Efforts to address these discrepancies include several physical processes. For instance, radiative cooling flows, which incorporate gas dynamics and self-photoionization, allow for intermediate ionization states in the CGM (O. Gnat & A. Sternberg 2009;S. Ji et al. 2019). Turbulent mixing layers, where Kelvin-Helmholtz instabilities blend hot and cool gas, provide another mechanism for producing transition-temperature ions such as C IV and O VI (J. D. Slavin et al. 1993;L. Armillotta et al. 2017;D. B. Fielding et al. 2020). Conductive interfaces, where cool clouds evaporate and hot gas condenses at the surface, add thermal structure to the CGM (L. L. Cowie & C. F. McKee 1977). Finally, radiative shocks, potentially driven by galactic winds, enhance ionization and line emission in certain regions (O. Gnat & A. Sternberg 2009;B. P. Wakker et al. 2012). Each of these processes has been shown to influence the observed ionization state of the CGM, but none fully resolves the discrepancies between models and observations.

An important aspect of the CGM is the influence of cosmic rays (CRs). CRs provide additional pressure support, affect thermal balance, and may drive galactic winds, all of which shape the overall structure and evolution of the CGM. Their ability to propagate over vast distances and recoup their losses by extracting energy from ambient shocks and turbulence allows CRs to impact regions far from their origin, influencing the thermal and dynamical properties of the CGM on a larger scale (M. Salem et al. 2016; I. S. Butsky & T. R. Quinn 2018; S. Ji et al. 2020). At the same time, their interactions with the disk-halo interface can drive multiphase winds or fountain flows, with CR heating and forces dominating the dynamics in low-β halos (T. H. N. Tsung et al. 2023). However, in the absence of observational evidence for CRs in the CGM, we can only speculate upon their influence; and because direct detection of CRs in diffuse environments such as the CGM is challenging, we must resort to observational predictions based on modeling the thermal gas itself. Recent studies have shown that CRs can also play a critical role in maintaining the local thermal stability of the CGM by offsetting runaway cooling in dense regions (I. S. Butsky et al. 2020). Furthermore, simulations by T. H. N. Tsung et al. (2023) demonstrated that CRs can suppress local thermal instability via entropy mode propagation, though this effect operates within a narrow parameter regime before CR heating destabilizes the background gas. Their results also show that CR-driven dynamics at the disk-halo interface, particularly in low-β halos, can lead to the formation of single-phase winds or multiphase fountain flows with significant turbulence and mass loading.

The transfer of momentum and energy between CRs and thermal gas is mediated by magnetic field structures on a scale of order the CR gyroradius. These small-scale magnetic structures could be part of an extrinsic turbulent cascade, but can also be driven by the CRs themselves via the so-called streaming instability (R. Kulsrud & W. P. Pearce 1969). The latter case is known as CR self-confinement, and we will assume throughout this paper that it holds. An interesting consequence of the theory is that it breaks down if the projections of the CR pressure gradient and thermal plasma density gradient along the background magnetic field oppose each other. In this case-as predicted from analytical arguments by J. Skilling (1971) and shown first in a simulation by J. Wiener et al. (2017)-scattering disappears, the CRs decouple from the gas, and their pressure becomes constant. Such transport bottlenecks should be common in the CGM, with its clumpy density structure, provided the clumps are magnetically connected to the galactic disk, where most CRs presumably originate.

Bottlenecks were studied in a series of papers of increasing generality in J. Wiener et al. (2017Wiener et al. ( , 2019)), C. Bustard & E. G. Zweibel (2021), and T. H. N. Tsung et al. (2022). These works showed that the large CR pressure gradient accelerates cool clouds upward, imparting significant momentum and causing changes in gas density and temperature profiles. Additionally, CRs were found to pressurize the cold gas, reducing its density and heating the cloud's surface, leading to broader transition layers and altering ion abundances. These findings underscore the significant role of CRs in shaping the thermal and dynamical properties of the CGM. Due to the disparity in scale between the assumed cloud size and the front thickness, the transition region has only been modeled previously as a 1D structure in a steady state (J. Wiener et al. 2017).

The transition layers between the cool and hot phases, characterized by a density contrast typically around a factor of 100 as temperatures shift from 10 4 K to 10foot_0 K, are particularly affected by CRs. Models of these transition layers are crucial not only for explaining observed line ratios but, as was shown in J. Wiener et al. (2017), may also be indicators of clouds dominated by CR pressure. This dominance would signify that CRs are a significant component of the galactic structure, affecting both the local environment of the CGM and the broader galactic ecosystem; it would also be evidence for the theory of CR self-confinement.

In this study, we refine the 1D plane-parallel symmetric model presented by J. Wiener et al. (2017) by incorporating thermal conduction and gas flow in the energy equation in addition to CR heating and radiative cooling. These additions aim to broaden the range of CR-mediated front models, and provide a more comprehensive understanding of the CGM's structure and dynamics.

We acknowledge that our model omits a number of potentially important effects, chiefly photoionization, which can have a large effect on the column densities of transition-temperature ions, turbulent mixing, which brings fluid elements at very different temperatures into close proximity (M. C. Begelman & A. C. Fabian 1990; C.-A. Faucher-Giguère & S. P. Oh 2023, and references therein), and the effects of magnetic fields, both dynamical and as an enforcer of anisotropic heat and CR transport. We have kept the problem simple to minimize the number of free parameters and unknown functions, arrive at a relatively complete level of physical understanding, and, through observational tests, probe the importance of the small number of processes we do retain.

The plan of this paper is as follows. In Section 2.1, we briefly review some well-known front models. In Section 2.2, we formulate the problem, with a separate discussion of the boundary conditions in Section 2.3. Section 3 gives detailed discussions of two types of front, static (Section 3.1) and evaporative (Section 3.2). We show the effects of changing the background magnetic field and CR pressure in the cloud in Section 4. In Section 5, we compute the column densities of several transition-temperature ions across our range of front models and compare their ratios to the ratios derived by B. P. Wakker et al. (2012) from absorption-line spectroscopy. Section 6 summarizes and discusses our main conclusions. Appendix A gives a short discussion of condensation fronts. Appendix B shows some results on how the column densities of transition-temperature ions depend on the domain over which we integrate.

Physical Setup

Review of Front Models

Past work on fronts without CRs provides valuable context for understanding the distinct role of CRs in gas dynamics. L. L. Cowie & C. F. McKee (1977) investigated evaporative fronts in the interstellar medium (ISM), where thermal conduction from a surrounding hot medium drives the continuous evaporation of cooler gas clouds. Their model established a steady-state balance between thermal conduction and enthalpy flux as the heating and cooling terms, respectively, with mass being lost from the cool phase to the hot phase. They also demonstrated how the geometry of the system influences the stability and evolution of these fronts. T. Inoue et al. (2006) later examined the phase-transition layers connecting two thermal equilibrium phases, a cold neutral medium cooled by C II and a warm neutral medium cooled by Lyα, with heating provided by photoelectrons from dust. Because both studies solved for temperature assuming energy balance and heat transfer by classical conduction, they were treated as two-point boundary value problems. 6 In L. L. Cowie & C. F. McKee (1977), the temperature of the hot medium was chosen based on observations and the cloud temperature was approximated by zero. T. Inoue et al. (2006) set the boundary temperatures to their equilibrium values.

In J. Wiener et al. (2017), thermal conduction and enthalpy flux were dropped, leaving collisionless CR heating supplemented by a base heating rate per particle to balance radiative cooling. Without heat conduction, the energy equation was reduced to first order. The models were integrated outward for 100 pc, at which point T had typically increased from its chosen value of 10 4 K at the cloud boundary to a few 10 6 K. Over this temperature range it was not possible to find a reasonable high-temperature thermal equilibrium, so the model was solved as an initial value problem and the temperature was still increasing at the 100 pc cutoff. It was checked a posteriori that the fronts found from the model were too thick for heat transport by thermal conduction to play a significant role. We expand on this issue in Section 3.

It should be acknowledged explicitly that the models discussed above are the simplest ones possible from a dynamical point of view. Approximate pressure balance is maintained, and the only flows present are highly subsonic enthalpy fluxes. In particular, turbulent mixing is ignored. Magnetic fields play no role in the dynamics and due to their orientation have no effect on heat transport. While L. L. Cowie & C. F. McKee (1977) showed that global geometry (spherical versus planar) plays a role, this cannot be explored in models with uniform magnetic fields.

Gas and Cosmic-ray Equations

In this work, we assume that CRs are self-confined via the streaming instability, with their propagation governed by interactions with Alfvén waves. The streaming speed is set by the local Alfvén velocity, and the resulting heating arises from the damping of Alfvén waves by the thermal plasma. This approach neglects CR diffusion or extrinsic confinement due to magnetic turbulence, as we focus on regions of the CGM where self-confinement likely dominates. Additionally, we neglect hadronic losses, as the low gas density in the CGM (n ≲ 0.1 cm -3 ) makes these interactions infrequent and their contribution to energy losses negligible. Including diffusion or hadronic losses would alter the distribution of CR energy and heating, but these effects are unlikely to significantly affect the results under the conditions considered here.

Under these assumptions, the evolution of the gas and CRs is governed by the following set of fluid equations, which describe the conservation of mass, momentum, and energy. In 1D, they are

where ρ and v represent the gas density and velocity, respectively; v A is the Alfvén velocity / B 4 ; P g and E g denote the thermal pressure and total energy density of the gas, with

κ T stands for the thermal conduction coefficient; and L is the net cooling function. For simplicity, we assume a fully ionized hydrogen plasma and an ideal gas equation of state, and take γ g = 5/3 throughout. We use a net cooling function L n n 2 . We adopt Γ = 1 × 10 -26 erg s -1 , consistent with the value used in T. Inoue et al. (2006). Because of the sharp cooling drop-off near T = 10 4 K, this choice of Γ allows for thermal equilibrium at both T ∼ 10 4 K and T ∼ 10 6 K. We note that this value of Γ is higher than the lower limit set by Coulomb scattering of CRs in thermal plasma (e.g., T. M. Yoast-Hull et al. 2013). For the cooling function, we use

( ) ( ) ( ( )) / = × × T 1.1 10 10 erg cm s , 4 T 21 log 10 K 3 1 5 where ( ) ( ) ( ) = + + + + x x e e 0.4 3 5.2 . 5 x x 0.08 1.5 0.08

This analytical fit for the cooling function was first given in S. Imada & E. G. Zweibel (2012) and modified in J. Wiener et al. (2019). The fit is designed to match the solar-metallicity cooling rates from R. P. C. Wiersma et al. (2009), where they used the Cloudy photoionization code to calculate cooling functions for gas with varying metallicities, temperatures, and densities under a UV/X-ray background at z = 0.

v 0 400 m s -1 40-400 m s -1 Initial gas velocity. P g,0

1.38 × 10 -13 erg cm -3 Fixed Initial gas pressure. P c,0 (α -1)P g,0 2P g,0 -9P g,0 Initial cosmic-ray (CR) pressure. α 3 3-10 Initial CR-to-gas pressure ratio.

In the fluid approximation, CRs are governed by the energy equation (J. F. McKenzie & H. J. Voelk 1982;D. Breitschwerdt et al. 1991):

where Q accounts for CR sources and sinks, and F c represents the CR energy flux:

Here, κ c is the CR diffusion coefficient, and we take γ c ≡ 1 + P c /E c = 4/3 throughout the paper.

The influence of CRs is manifested through the CR pressure (P c ) gradient term in the momentum equation and the heating term in the energy equation. As shown in J. Wiener et al. (2013) and E. G. Zweibel (2017), streaming CRs heat gas at the rate of -v A ∇ P c . At a microscopic level, the heating represents the transfer of CR energy to Alfvén waves due to the streaming instability together with the transfer of wave energy to the thermal background, as must occur in a steady state. The pressure gradient force represents momentum transfer between the CRs and the magnetic field fluctuations, and occurs whether CRs are self-confined or confined by extrinsic turbulence.

We seek steady-state solutions of Equations ( 1)-( 6) under the further simplifying assumptions that CR diffusion is negligible and v/v A ≪ 1. That is, we set all time derivatives to zero, drop κ c , and wherever v + v A appears, we replace it by v A . This leads to the following three conserved quantities:

At the x = 0 boundary, gas pressure (P g,0 = n 0 k B T 0 ) and CR pressure (P c,0 ) are specified. To parameterize the relative contribution of CR pressure, we introduce the parameter α ≡ P c,0 /P g,0 + 1 to quantify the ratio of CR pressure to gas pressure at the boundary.

Substituting these relationships into the momentum and energy equations, we obtain the following relationship between number density (n) and temperature (T):

where the term αn 0 k B T 0 accounts for the total pressure (gas +CR) at the boundary. This allows us to treat T as a function of n. We then recast Equation (3) as a second-order ordinary differential equation (ODE) for n. We approach this ODE system as an initial value problem to determine the spatial profile of the variables. The boundary values are established within the warm ionized medium, as detailed in Section 2.3.

Initial Value Problem: Boundary Values

In the warm ionized medium on the left boundary (x = 0), we set the proton number density (n p,0 ) at 0.05 cm -3 , resulting in a total (electron+proton) number density (n 0 ) of 0.1 cm -3 . The initial temperature (T 0 ) is established at 10 4 K, leading to an initial gas pressure (P g,0 = n 0 k B T 0 ) of 1.38 × 10 -13 erg cm -3 . These boundary conditions are chosen to represent typical conditions in the warm ionized medium of the CGM. In this environment, the thermal conduction coefficient is given by κ T = 5.6 × 10 -7 T 5/2 erg s -1 K -1 cm -1 (L. Spitzer 1962).

In our fiducial case, we adopt a constant magnetic field strength B = 3 μG and assume the field is perpendicular to the plane of the front. This choice aligns with observations of the ∼microgauss magnetic field strength in the CGM (T.-W. Lan & J. X. Prochaska 2020; K. Böckmann et al. 2023;V. Heesen et al. 2023).

For the evaporative solutions, we set the gas velocity at the boundary (v 0 ) to 400 m s -1 , neglecting the effects of any relative motion between the cloud and the ambient medium. profiles for a static-front system, where the velocity v = 0 and only thermal conduction, CR heating, and radiative cooling interact. The profiles are differentiated by the initial gradients of number density, dn dx 0 , across a range from 10 -14 to 10 -17 cm -4 . We emphasize that except for the initial gradients, all profiles originate from the same initial conditions. The top panel shows the number density profiles, where steeper initial gradients lead to a more rapid decline in number density from the origin. In contrast, the less steep initial gradients see a substantial and sharp decrease around ∼1 pc, indicating a delayed transition. The middle panel shows the corresponding temperature profiles. The different profiles underscore the significant impact of the initial gradient on the structure of the front. The bottom panel illustrates the pressure profiles; while the initial CR pressure is twice that of the gas pressure, it gradually decreases as the density drops.

The positive value of the mass flux corresponds to the warm medium evaporating into the hot medium. We also explore a scenario with a static front, where v = 0. We show in Appendix A that a condensation front cannot form under these conditions. In our setup, we have v A ≫ v, which satisfies the condition for the validity of Equation (3). In our fiducial case, we assume P c,0 = 2P g,0 (a value consistent with the bottleneck formation studies of J. Wiener et al. 2017Wiener et al. , 2019)). Specific values of α are varied in this work to explore their effect on the structure of the transition layers. We list all relevant parameters in Table 1.

Results

We proceed by solving the energy equation (Equation (3)), reformulated as a second-order initial value problem for n, imposing two boundary conditions: n 0 = 0.1 cm -3 and dn dx 0 , the value of which is varied to explore the effect of the gradient (which determines the strength of various flux terms) on our outcomes. We defer discussions of how dn dx is selected in a natural environment to future works, except for a few remarks in Section 6. In J. Wiener et al. (2017), the front models are generated by integrating a first-order ODE, and dn dx 0 is determined through the polytropic relationship between P c and n, along with the requirement that CR heating balances radiative cooling.

Static Front

In this section, we first examine the simpler static-front case (where v = 0) to isolate and better understand the interaction between thermal conduction, CR heating, and radiative cooling. We numerically solve the energy Equation (3),

T c A using forward integration, and show the resulting number density, temperature, and pressure profiles as a function of distance in Figure 1. Following J. Wiener et al. (2017), we integrate to 100 pc. These profiles are distinguished by their initial gradients of number density, dn dx 0 , with values ranging from 10 -14 to 10 -17 cm -4 . All the models are initialized with CR pressure dominant (P c0 /P g0 = 2), and all become gas pressure dominated at large distances, typically by about an order of magnitude. This is due to the polytropic relationship between P c and ρ, whereas the decrease in ρ is compensated by an increase in T. As shown in Figure 1, steeper initial gradients (e.g., 10 -14 cm -4 , 10 -15 cm -4 ), which imply stronger influences from thermal conduction and CR heating, result in a noticeable and immediate decrease in number density from the origin. In contrast, profiles with less steep initial gradients (e.g., 10

dn dx 0 16 cm -4 ) only show a significant decrease in number density at ∼1 pc, indicating a delayed transition. Regardless of the initial gradient, the number density always decreases to of order 10 -3 cm -3 at x = 100 pc, with the steepest gradient reaching ∼1 × 10 -3 cm -3 and the least steep gradient dropping to ∼3 × 10 -3 cm -3 . The same behavior is observed in the temperature profiles in the middle panel of Figure 1. Lower ). This figure delineates the varying dynamics between different flux components under different initial number density gradients. For smaller initial number density gradients, CR heating first is balanced by conductive cooling, until it decreases to the level of radiative cooling, which is when conduction becomes the less dominant heating term. For larger initial gradients, the transition from conductive cooling to heating occurs over a more extended distance (beyond 100 pc). The interplay between these components is critical to understanding the thermal dynamics of the system.

countered by CR heating and conductive heating. This balance continues all the way to the boundary at x = 100 pc. It is not surprising that all these models require conductive cooling near the cloud, since the gradient required to balance CR heating with radiative cooling is of order 10 -21 cm -4 , much smaller than any of the models presented here.

Comparing our static-front results with previous work highlights the influence of CRs on the energy balance. In L. L. Cowie & C. F. McKee (1977), without CRs, conductive heating balances radiative cooling at the front. In contrast, when CRs are included, as in our model and that of J. Wiener et al. (2017), CR heating becomes the dominant process balancing radiative cooling. Nonetheless, thermal conduction remains significant, as strong CR heating requires conductive cooling to dissipate the excess energy. Our result highlights the crucial role of thermal conduction in maintaining thermal equilibrium within the system.

Evaporation Front

In this section, we incorporate an evaporative flow into the energy equation:

We now numerically solve Equation (11), and present the resulting number density, temperature, and pressure profiles in Figure 3. For comparison, the static-front cases are shown with dashed lines. Notably, the profiles of the evaporation fronts begin to deviate from those of the static case beyond 0.1 pc, with these deviations amplifying with increasing distance. Specifically, the number density in the evaporationfront scenario drops more significantly compared to the static case, which corresponds to a greater increase in temperature.

pressure profiles show that as we move away from the origin, the CR and gas pressure profiles evolve similarly to the static-front case. However, after the point where P c = P g , the difference between P c and P g becomes greater compared to the static-front scenario, indicating a sharper transition between the roles of CR pressure and gas pressure. The presence of flow introduces dynamical pressure, which is absent in the static case. In the evaporation-front scenario, the pressure profile shows a continuous decrease in the sum of CR and gas pressures, accompanied by an increase in dynamical pressure. However, the dynamical pressure remains relatively small, accounting for only about 3% of the total pressure at x = 100 pc, emphasizing the dominant role of CR and gas pressures in shaping the dynamics of the layer.

We also present the contributions of each flux component as a function of distance in Figure 4. For less steep initial gradients (e.g., 10 dn dx 0 16 cm -4 ), we identify three distinct regions. Near the origin, CR heating is the dominant process and is balanced by conductive cooling. Further from the origin, there is a transition region where CR heating decreases to levels comparable to radiative cooling, causing conduction to shift from cooling to heating. This transition indicates a regime where CR heating is balanced by radiative cooling. Beyond this point, as radiative cooling reaches its maximum and subsequently declines, conductive heating and enthalpy flux cooling become the dominant processes, similar to the L. L. Cowie & C. F. McKee (1977) evaporation fronts. For steeper initial gradients (e.g., 10

dn dx 0 15 cm -4 ), the dynamics change significantly. Strong CR heating prevents radiative cooling from fully balancing the thermal energy input, leading to the absence of a region where radiative cooling balances CR heating and thereby precluding the formation of regions with flat number density and temperature profiles. These results underscore the complex interplay between CR heating, radiative cooling, and conduction in regulating the thermal dynamics of the system. The middle panel shows the corresponding temperature profiles. The different profiles underscore the significant impact of the initial gradient on the structure of the front. The bottom panel illustrates the pressure profiles; while the initial cosmic-ray (CR) pressure is twice that of the gas pressure, it gradually decreases as the density drops.

Varying Physical Parameters

In this section, we examine the effects of changing physical parameters, specifically the magnetic field strength and the ratio between CR pressure and gas pressure at the x = 0 boundary, denoted as α ≡ P c,0 /P g,0 + 1. The magnetic field strength appears only as a multiplier in the CR heating term. Meanwhile, the ratio α sets the initial CR pressure and also determines the system's total pressure, with a higher α indicating increased initial CR pressure and total pressure. We investigate how variations in these parameters impact the dynamics and structure of the layer.

Cosmic-ray Pressure and Total Pressure

In this subsection, we discuss the results for α = 10. As shown in Figure 5, with this higher value, at the x = 100 pc boundary the number density does not drop as low, reaching only ∼10 -2 cm -3 for 10 dn dx 0 15 cm -4 , compared to ∼10 -3 cm -3 in the fiducial case. The impact of α on the temperature profile is smaller, however. Additionally, the role of dynamical pressure across the layer becomes much less significant. At the x = 100 pc boundary, dynamical pressure constitutes only 0.8% of the total pressure, in contrast to the α = 3 case, where dynamical pressure at the same boundary is comparable to CR pressure.

Examining the component plot (Figure 6) reveals the impact of α. For 10 dn dx 0 15 cm -4 , stronger CR heating at the origin extends the distance needed for conduction to change from cooling to heating. This difference is illustrated in the top two panels of Figure 6. We also observe that for 10 dn dx 0 16 cm -4 , the number density gradient, and hence the temperature gradient, are steeper than their α = 3 counterparts (bottom panels of Figure 6). Additionally, after the peak in radiative cooling, the α = 10 scenarios require a longer distance for radiative cooling and CR heating to fall to levels similar to enthalpy flux and conduction. This extended distance is reflected in the more gradual changes in the number density and temperature profiles in this region. The component plots demonstrate how the strength of CR heating at the origin impacts the system's thermal structure and the spatial extent of physical processes. These observations are similar to the α = 3 case, where the presence of CR heating leads to the adjustment of other terms in the energy equations to counteract CR heating.

Magnetic Field Strength

We now investigate the impact of a heightened magnetic field strength, B, set to 30 μG. Note that when only B is adjusted, the total pressure of the system remains the same; only the CR heating term increases by a factor of 10. With this stronger magnetic field, achieving number density and temperature profiles as steep as in previous cases requires steeper initial number density gradients. This change is clearly shown in Figure 7. A key difference between the B = 30 μG result and the α = 10 result is that, with constant total pressure, the temperature does not reach ∼10 6 K, as it does in the fiducial case and in the α = 10 case. Similarly, the enhanced magnetic field reduces the relative importance of dynamical pressure compared to the trend in the fiducial case (see bottom panels of Figure 7).

Turning to Figure 8, we see that a tenfold increase in the magnetic field strength causes CR heating to decrease more rapidly, aligning more quickly with radiative cooling levels. For all initial density gradients, by x = 100 pc the relatively smaller decrease in number density compared to the fiducial case means that enthalpy flux cooling and conductive heating are no longer the dominant processes. Instead, CR heating and radiative cooling become dominant at the x = 100 pc 15 cm -4 ), CR heating is so strong that radiative cooling cannot balance it, precluding the formation of regions with flat number density and temperature profiles.

boundary. This sustained dominance of CR heating and radiative cooling at larger scales is a direct outcome of the stronger magnetic field's influence on the system.

It may seem paradoxical that increasing B, and with it v A , increases the heating rate but decreases T far from the cloud. We can understand this by observing that as B increases, so does the range over which CR heating is balanced by conductive cooling. This in turn requires that T cannot increase too quickly. For example, it can readily be shown that if T has a power-law dependence on position, T ∝ x p , conductive cooling requires p < 2/7. This accounts for the relatively slow increase in T in the high-magnetic-field models.

Ion Column Density Ratios

We use the ionization fractions X i (T) from O. Gnat & A. Sternberg (2007) and the relative abundances A X from M. Asplund et al. (2009) to compute the column densities of specific ions. The column density for each ion is defined by N X = A X ∫X i (T)n dx, where n is the number density, and the integral is across the front, up to a distance of 100 pc. We discuss this choice of integration boundary in Appendix B. We then compare our calculated line ratios with measurements from high-velocity clouds (HVCs) observations in the Milky Way from B. P. Wakker et al. (2012). Discrepancies between our theoretical line ratios and the observational data could highlight potential shortcomings in the current models.

We compute the column densities for C IV, O VI, and Si IV, and derive the diagnostic line ratios Si IV/C IV, C IV/O VI, and N V/O VI. These ratios serve as key indicators for understanding the ionization processes and thermal structure within the ISM. The Si IV/C IV ratio probes cooler, photoionized gas with temperatures in the range of 10 4 -10 5 K, primarily indicating regions where UV radiation dominates the ionization processes (O. Gnat & A. Sternberg 2007;B. P. Wakker et al. 2012). In contrast, the C IV/O VI ratio spans the transition between photoionized and hotter, collisionally ionized gas, typically found at temperatures between 10 5 and 10 6 K. This ratio is particularly effective in identifying multiphase structures within the ISM, often tracing shockheated gas in more dynamic regions (J. M. Shull et al. 2012;B. P. Wakker et al. 2012). Finally, the N V/O VI ratio is sensitive to collisionally ionized gas in the warm-hot ionized for the evaporation-front scenario with a high CR-to-gas pressure ratio at the x = 0 boundary (α = 10). In the top panels, the solid lines represent α = 10, while the dashed semitransparent lines depict the profiles for α = 3. The top panels illustrate the impact of α on the number density and temperature profiles, showing that the number density reaches only ∼10 -2 cm -3 at x = 100 pc for 10 dn dx 0 15 cm -4 , compared to ∼10 -3 cm -3 in the fiducial case. The bottom panels show the pressure profiles for the highest and lowest initial density gradients in blue and red, respectively, indicating that the transition from CR pressure dominated to gas pressure dominated persists at higher α, while dynamical pressure constitutes only 0.8% of the total pressure at the x = 100 pc boundary, in contrast to the α = 3 case where it is comparable to CR pressure. medium (or WHIM), making it a useful diagnostic for gas in intermediate ionization states and the presence of shock-heated or disturbed gas (B. D. Savage & B. P. Wakker 2009; P. Richter et al. 2017).

We show in the left panel of Figure 9 our model-dependent Si IV/C IV ratio versus C IV/O VI ratio in logarithmic scale, while in the right panel we present N V/O VI ratio versus C IV/ O VI ratio, overplotted with observational data from B. P. Wakker et al. (2012). The black circular data points with error bars represent the observational measurements, and our model predictions are superimposed as colored geometric shapes. Additionally, we include a CR heating balancing radiative cooling model, as discussed in J. Wiener et al. (2017), shown in orange.

Our results reveal that most of the models lie well outside the observational range in the Si IV/C IV versus C IV/O VI parameter space. The highest magnetic field (B) value models predict Si IV/C IV ratios that are consistent with the observations; however, these models simultaneously predict C IV/O VI ratios that are too high compared to the observational data. In contrast, models with lower B values produce Si IV/C IV ratios that are significantly lower than those observed in the HVCs.

Similarly, in the right panel, the high-B models systematically overestimate the N V/O VI ratio, placing them outside the range of observational constraints. For the other models, those with the largest dn dx 0 underestimate the C IV/O VI ratio, suggesting that the structure and dynamics of the gas in these models may be insufficient to fully capture the physical processes driving the ionization conditions. However, the range of models that overlap the data plotted in (N V/O VI, C IV/O VI) space is quite large, suggesting that this is not a particularly sensitive diagnostic of conditions in the front, including the initial value of the density gradient.

A prima facie interpretation of these results suggests that if a model lies within the observational range in the N V/O VI versus C IV/O VI parameter space but fails to reproduce the correct Si IV/C IV ratio, it implies that the model is better at predicting the ionization profile for gas in the 10 5 -10 6 K temperature range than for cooler gas at 10 4 -10 5 K. This also indicates that the model more accurately captures conditions dominated by collisional ionization rather than photoionization. These discrepancies highlight potential limitations to the current model framework, particularly regarding the role of magnetic fields and CR heating in shaping the ionization state of the ISM and CGM. In fact, the line ratios are not well fit by previous models which do not include CRs (see Figure 11 of B. P. Wakker et al. 2012). Our α = 1 models, as shown in Figure 9, align with this finding.

It should not be overlooked, however, that some of our models do successfully reproduce both pairs of line ratios. As might be expected from the forgoing discussion, these are the models with intermediate field strength: B = 20 μG for α = 3. As we already noted, the fits are relatively insensitive to the initial density gradient dn dx 0 , which sets the initial CR heating rate and which we regard as the least well-determined parameter in our models.

Summary and Discussion

In this paper, we explored the evolution of the transition layer between warm (10 4 K) and hot (10 6 K) gas phases in the CGM by integrating the energy equation (Equation (3)) as an initial value problem. We showed that CR heating, thermal conduction, radiative cooling, enthalpy flux, and dynamical pressure collectively shape the temperature and density profiles of the transition layer. At large distances from the cloud, we found that conduction is the main heat source and that it can be balanced by an outward enthalpy flux, as in the models of L. L. Cowie & C. F. McKee (1977).

We also found that the temperature and density profiles across the front are quite sensitive to the value we adopt for the 15 cm -4 , stronger CR heating at the origin extends the distance needed for conduction to change from cooling to heating. The bottom two panels show that for 10 dn dx 0

16 cm -4 , the number density gradient, and hence the temperature gradient, are steeper than their α = 3 counterparts. The component plots demonstrate the impact of CR heating on the front structure and the spatial extent of physical processes. initial density gradient. When the gradient is sufficiently small, radiative cooling suffices to balance CR heating, as assumed in the front models of J. Wiener et al. (2017). For sufficiently large gradients, CR heating is balanced by conductive cooling. In view of this sensitivity, we have to ask what controls this ratio in a natural environment. While addressing this question is beyond the scope of this paper, we speculate that the answer lies either in the front formation process or in the stability of the front to small perturbations. These possibilities will be examined in future work.

Our findings emphasize the dominant role that CRs can play in regulating the structure of the transition layer between the warm and hot gas phases of the CGM. The addition of CR heating in these models dictates the thermal balance, with the other flux terms adjusting accordingly to balance CR heating. These effects are sensitive to the magnetic field strength and the ratio of CR to gas pressure (α) at the warm gas boundary.

Our analysis of ion column density ratios shows that while our model predictions align with some observational data, and some represent successful fits, discrepancies persist. Only the highest-magnetic-field models reproduce the observed Si IV/ C IV ratios, but they overestimate the C IV/O VI ratios.

Conversely, other models underestimate the Si IV/C IV ratios, although some fall within the observed range for C IV/O VI and N V/O VI. The agreement with C IV/O VI and N V/O VI, but not with Si IV/C IV, suggests that these models better represent ionization conditions at higher temperatures (10 5 -10 6 K), where collisional ionization dominates, but struggle to capture the profiles in cooler, photoionized regions. These discrepancies indicate that while CR heating significantly influences the ionization structure of the CGM, additional processes such as turbulent mixing, shocks, or self-ionization may be necessary to fully explain the observed line ratios in lower-temperature regimes. The need to include such effects would be particularly compelling if magnetic field strengths as high as 20 μG (the magnetic field strength adopted in our most successful models) could be ruled out.

In addition to the influence of CRs, our study highlights the crucial role of thermal conduction in maintaining energy balance within the transition layers. Thermal conduction dissipates excess heat from CR interactions, preventing runaway heating and stabilizing the temperature profiles.

Future work will focus on investigating the stability of these transition layers, particularly in the presence of CR pressure, magnetic fields, and thermal conduction. Studying the onset , and pressure (bottom) profiles for the evaporation-front scenario with a magnetic field strength of 30 μG and a CR-to-gas pressure ratio (α) of 3. In the top panels, the solid lines represent the current scenario, while the thin lines show the profiles with a magnetic field strength of 3 μG. The profiles indicate that with a stronger magnetic field, dynamical pressure is reduced. Additionally, achieving similarly steep profiles as seen with a weaker magnetic field requires steeper initial number density gradients.

and development of instabilities within these fronts will be crucial for understanding their behavior. However, our models do not include turbulent mixing or magnetic tension, both of which could further alter the transition layer structure. Turbulent mixing, induced by velocity shear or instabilities such as the Kelvin-Helmholtz instability, could broaden the transition layers by mixing gas from the warm and hot phases. This would enhance the presence of intermediate-temperature gas, potentially boosting ion column densities and modifying ionization ratios. Such an effect might alleviate some of the discrepancies between our predicted and observed line ratios, particularly in cases where our models underestimate the Si IV/C IV ratios. However, turbulence could also disrupt CR bottlenecks, leading to a redistribution of CR heating and broader temperature and density profiles, potentially exacerbating tensions with observational constraints.

Magnetic tension, on the other hand, acts as a stabilizing force that resists the deformation of field lines, potentially suppressing turbulent mixing and preserving sharp boundaries between the warm and hot phases. This suppression would maintain a steeper density and temperature gradient, limiting the formation of intermediate-temperature gas and thereby reducing the column densities of transition-temperature ions. 1.0 0.5 0.0 0.5 1.0 log [N(CIV)/N(OVI)] 2.0 1.5 1.0 0.5 log [N(SiIV)/N(CIV)] =1 (No CR) B=3 G, =3 B=3 G, =10 B=20 G, =3 B=30 G, =3 1.0 0.5 0.0 0.5 1.0 log [N(CIV)/N(OVI)] 1.0 0.5 0.0 0.5 log [N(NV)/N(OVI)] dn dx |0 = 10 14 cm 4 dn dx |0 = 10 15 cm 4 dn dx |0 = 10 16 cm 4 dn dx |0 = 10 17 cm 4 CR heating + radiative cooling Figure 9. Left panel: Si IV/C IV vs. C IV/O VI ratios in logarithmic scale, compared with observational data from B. P. Wakker et al. (2012, black circles with error bars). Right panel: N V/O VI vs. C IV/O VI ratios, similarly overlaid with observational data. Colored shapes denote our model predictions, with variations in magnetic field strengths and CR parameters. The orange shapes represent models where CR heating balances radiative cooling, following J. Wiener et al. (2017). Both panels reveal discrepancies between the models and observations, particularly in the Si IV/C IV vs. C IV/O VI space, where no models fit within the observed range. High-magnetic-field models match Si IV/C IV but overestimate C IV/O VI, while lower-field models miss Si IV/C IV but align with C IV/O VI and N V/O VI.

This suggests that, in general, the models are more accurate at higher temperatures but struggle with cooler, photoionized conditions. Note that models with intermediate field strength, B = 20 μG, do overlap the data on both plots, although we make no claim to their uniqueness.

Moreover, magnetic tension modifies CR transport, enforcing stronger anisotropy and concentrating CR heating along magnetic field lines. This could lead to localized enhancements in ionization, altering the ion column density ratios in a nontrivial way. These effects introduce additional complexity into the structure of the CGM transition layers, and whether they alleviate or exacerbate tensions in our model predictions depends on the dominant physical scales of mixing and field alignment. Understanding the interplay between CR pressure, magnetic fields, thermal conduction, and turbulent mixing is therefore crucial for interpreting CGM structure and ionization conditions.

The line ratios we presented are sensitive to the cutoff point chosen for the integration over the transition layer. In Figure 10, we plot the line ratios (log(C IV/O VI), log(Si IV/ C IV), and log(N V/O VI)) as functions of temperature where the integration is cut off, with the shaded regions representing the range of observational data. The vertical lines indicate the temperature at x = 100 pc for the two models: the red vertical line corresponds to the temperature for the model with = 10 dn dx 0 17 cm -4 , and the blue vertical line corresponds to the model with = 10 dn dx 0 14 cm -4 . As seen in Figure 10, the choice of temperature cutoff has a significant, nontrivial effect on the calculated line ratios. J. Wiener et al. (2017) used a cutoff at 2 × 10 6 K, which aligns more closely with the converged values in the highertemperature regime, compared to cutting off at 100 pc for the red model. However, while extending the integration to higher temperatures may yield more stable line ratios, the physical thickness of the transition layers becomes an important consideration. The models suggest that line ratios converge at higher temperatures, but the temperature at x = 100 pc remains well below this threshold for both models. This highlights a key issue: While we can theoretically integrate to very high temperatures, the question remains whether the transition layers are actually that thick in reality. Currently, we lack observational probes that can resolve the full extent of these layers. This uncertainty underscores the need for careful definition of the integration boundaries and comparison across models and observations.

ORCID iDs

Hanjue Zhuaa https:/ /orcid.org/0000-0003-0861-0922 Ellen G. Zweibelaa https://orcid.org/0000-0003-4821-713X Nickolay Y. Gnedinaa https:/ /orcid.org/0000-0001-5925-4580 17 cm -4 (red) and = 10 dn dx 0

14 cm -4 (blue), respectively. The plot illustrates the sensitivity of the calculated line ratios to the choice of integration cutoff.