Low- and High-velocity O vi in Milky Way-like Galaxies: The Role of Stellar Feedback

Milky Way-type galaxies are surrounded by a warm-hot gaseous halo containing a considerable amount of baryons and metals. The kinematics and spatial distribution of highly ionized ion species such as O vi can be significantly affected by supernova (SN) explosions and early (pre-SN) stellar feedback (e.g., stellar winds, radiation pressure). Here we investigate effects of stellar feedback on O vi absorptions in Milky Way−like galaxies by analyzing the suites of high-resolution hydrodynamical simulations under the framework of SMUGGLE, a physically motivated subgrid interstellar medium and stellar feedback model for the moving-mesh code Arepo. We find that the fiducial run with the full suite of stellar feedback and moderate star formation activities can reasonably reproduce Galactic O vi absorptions observed by space telescopes such as the Far-Ultraviolet Spectroscopic Explorer, including the scale height of low-velocity (∣v LSR∣ < 100 km s−1) O vi, the column density–line width relation for high-velocity (100 km s−1 ≤ ∣v LSR∣ < 400 km s−1) O vi, and the cumulative O vi column densities. In contrast, model variations with more intense star formation activities deviate from observations further. Additionally, we find that the run considering only SN feedback is in broad agreement with the observations, whereas in runs without SN feedback this agreement is absent, which indicates a dominant role of SN feedback in heating and accelerating interstellar O vi. This is consistent with the current picture that interstellar O vi is predominantly produced by collisional ionization where mechanical feedback can play a central role. In contrast, photoionization is negligible for O vi production owing to the lack of high-energy (≳114 eV) photons required.


INTRODUCTION
The multiphase gas within and surrounding galaxies including the Milky Way (MW) is an essential ingredient of galactic ecosystems that govern the galaxy evolu-tion, and may contain a significant amount of baryons and metals in the form of the cold (T ≲ 10 4 K), warm (T ∼ 10 5 − 10 6 K), and hot (T ≳ 10 6 K) gaseous phases (e.g., Putman et al. 2012;Tumlinson et al. 2017, and references therein).The existence of a warm-hot Galactic corona was originally proposed by Spitzer (1956) to provide pressure confinement to the neutral clouds that are ∼ 1 kpc above the Galactic plane, and was later confirmed by detections of the soft X-ray background (e.g., Bowyer et al. 1968) and interstellar O VI absorptions (e.g., Jenkins & Meloy 1974;York 1974).Spitzer (1956) also pointed out that such diffuse gas can be studied via the resonance doublet absorption lines of lithiumlike ions, e.g., O VI, N V, and C IV. Plasmas in the temperature range of about (1 − 5) × 10 5 K traced by these species can be produced via moderate shocks or rapid cooling of hotter coronal gas probed in X-rays.O VI λλ1032, 1038 doublet is of special significance owing to the large oscillator strengths (Morton 2003) and high cosmic abundance of oxygen.Under the condition of collisionally ionized equilibrium, O VI peaks in abundance at the temperature of ∼ 3 × 10 5 K (Sutherland & Dopita 1993).
The first large-scale surveys of O VI absorption in the Milky Way were made by the Far Ultraviolet Spectroscopic Explorer (FUSE; Moos et al. 2000;Sahnow et al. 2000).FUSE detections of O VI absorption lines toward extragalactic objects (e.g., active galactic nuclei (AGNs)/quasars) and stars in the Galactic disk, Galactic halo, and Magellanic clouds, reveal a widespread but irregular distribution of interstellar O VI with a column density of log(N/cm −2 ) ∼ 13.0 − 14.8 (e.g., Savage et al. 2000;Howk et al. 2002a;Wakker et al. 2003;Savage et al. 2003;Sembach et al. 2003;Oegerle et al. 2005;Bowen et al. 2008;Sarma et al. 2017).The O VI absorbers detected by FUSE and Hubble Space Telescope (HST) along the lines of sight (LOS) of quasars/stars move at various velocities with respect to the local standard of rest (LSR), i.e., |v LSR | ranges from < 100 to ≳ 400 km s −1 (e.g., Sembach et al. 2000;Murphy et al. 2000;Sembach et al. 2003;Fox et al. 2006;Collins et al. 2007;Shull et al. 2011).
Low-velocity (e.g., |v LSR | < 100 km s −1 ) O VI is believed to be an extension of the Galactic disk, inflated due to its relatively high temperature, and can be well approximated by an exponentially declined layer with a midplane density of ≲ 2 × 10 8 cm −3 and a scale height of ∼ 2.3 − 4 kpc (e.g., Savage et al. 2000Savage et al. , 2003;;Zsargó et al. 2003;Indebetouw & Shull 2004;Savage & Wakker 2009).In contrast, the nature of high-velocity (e.g., 100 ≤ |v LSR | < 400 km s −1 ) O VI as well as intermediate-and low-ions (e.g., O I, C II, Si II, Mg II, Si III, C IV, Si IV) and atoms, the so called high-velocity clouds (HVCs), is still debated, largely owing to the highly uncertain distances for most cases.While some high-velocity O VI features are spatially and kinematically associated with known H I structures (e.g., Complex C and the Magellanic Steam), some have no neutral counterparts detected (e.g., Sembach et al. 2003;Nicastro et al. 2003;Fox et al. 2004;Collins et al. 2004Collins et al. , 2005)).In addition, the covering fraction of high-velocity O VI (≳ 60%; e.g., Sembach et al. 2003;Fox et al. 2006) is found to be higher than that of neutral and moderately ionized HVCs (∼ 20% − 40% for H I, C IV, and Si IV; e.g., Lockman 2002;Herenz et al. 2013), indicating a spatially more extended distribution for highly ionized HVCs.Despite the multiple origins proposed for HVCs, for example, the Galactic fountain (e.g., Shapiro & Field 1976;Bregman 1980;Fraternali & Binney 2006), materials stripped or ejected from satellite galaxies (e.g., Putman 2004;Herenz et al. 2013), and accretion from the intergalactic medium (IGM; e.g., Fraternali et al. 2015;Kereš & Hernquist 2009), the spatial distribution and kinematics of high-velocity O VI are probably dominantly governed by the fountain model, which proposes that gas circulation in the halo is powered by stellar feedback, e.g., stellar winds and supernova (SN) explosions.Such a scenario is also supported by recently observed rain-like inflows and collimated outflows (e.g., Lehner et al. 2022;Marasco et al. 2022).
O VI absorptions for low-redshift galaxies have been extensively studied by HST and FUSE (e.g., Tripp et al. 2000;Tripp & Savage 2000;Danforth & Shull 2005;Lehner et al. 2006;Tumlinson et al. 2011;Prochaska et al. 2011;Savage et al. 2011;Stocke et al. 2013;Fox et al. 2013;Peeples et al. 2014;Mathes et al. 2014;Johnson et al. 2015;Kacprzak et al. 2015;Prochaska et al. 2019;Tchernyshyov et al. 2022).Strong O VI absorptions have been preferentially detected around starforming galaxies, with an average O VI column density of log(N/cm −2 ) ∼ 14.5 (e.g., Tumlinson et al. 2011), indicating a strong impact of star formation activities on the global properties of warm gaseous halo traced by O VI .Additionally, the covering fraction of O VI was found to depend on the inclination angle of galaxies and to follow a bimodal distribution that peaks within ∼ 30 • of the galaxy minor axis and ∼ 10 • − 20 • of the major axis (e.g., Kacprzak et al. 2015), consistent with a circumgalactic medium (CGM) originating from majoraxis-fed inflows/recycled gas and from minor-axis-driven outflows, i.e., a scenario also revealed by cooler gaseous phases traced by Mg II absorptions (e.g., Bouché et al. 2012;Kacprzak et al. 2012).Those observational evidences highlight the influence of star formation activities and stellar/AGN feedback in shaping the spatial distribution of O VI-bearing gas of external galaxies.
Stellar feedback, i.e., injection of substantial amounts of energy and angular momentum into the interstellar medium (ISM) via early (pre-SN) feedback and SN explosions, is likely to leave imprints on the properties of gaseous halos (e.g., Appleby et al. 2021;Mina et al. 2021).O VI ions probably trace matter in the interfaces between the cooler ionized/neutral clouds and hotter gas, and can thus serve as indirect probes of stellar feed-back (e.g., Lehner et al. 2011).Hydrodynamical simulations are powerful tools for studying the ISM/CGM, and given that multiscale physical processes are involved in galaxy formation, subgrid models are often used to implement small-scale processes such as star formation, metal mixing and transport, and stellar feedback (e.g., Cen & Ostriker 1992;Agertz et al. 2013;Hopkins et al. 2014;Li et al. 2017Li et al. , 2018;;Smith et al. 2018;Hopkins et al. 2018;Marinacci et al. 2019).Those subgruid models are parameterized and tuned to reproduce the observations, which means feedback energy is treated as a free parameter despite its well-known importance (e.g., Barbani et al. 2023).A variety of simulations have shown that stellar feedback can have a significant impact on the physical properties such as kinematics, column densities, and total content of O VI (e.g., Hummels et al. 2013;Marasco et al. 2015;Liang et al. 2016;Fielding et al. 2017;Li & Tonnesen 2020).
Stars and MUltiphase Gas in GaLaxiEs (SMUGGLE; Marinacci et al. 2019) is a physically motivated subgrid ISM and stellar feedback model for the movingmesh code Arepo (Springel 2010) and has been widely used since its development (e.g., Kannan et al. 2020;Burger et al. 2022;Sivasankaran et al. 2022;Barbani et al. 2023).It has successfully reproduced hydrogen emission line profile (Smith et al. 2022), constant density cores in dwarf galaxies (Jahn et al. 2023), and in particular, a realistic cold ISM and star cluster properties in isolated and merging galaxies (Li et al. 2020(Li et al. , 2022)).In this paper, we test whether the SMUGGLE model can reproduce observations of warm O VI gas in and around the Milky Way, and investigate how the properties of O VI gas are affected by stellar feedback (e.g., early feedback and SN explosions) by analyzing the suites of simulations presented in Li et al. (2020, hereafter L20).
The paper is structured as follows.Section 2 briefly introduces the SMUGGLE model and L20's simulation, and generates synthetic observations of O VI absorptions.Section 3 presents results and discussion on O VI properties for different feedback model variations, as well as comparison with the observations and some caveats.Section 4 summarizes the main conclusions.

METHODOLOGY
We analyze a suite of hydrodynamic simulations of isolated MW-sized galaxies presented in L20 under the SMUGGLE framework (Marinacci et al. 2019).We refer the reader to the original papers for details of the model and the simulations.Below we give a brief overview of the SMUGGLE model and L20's simulations, and describe the methodology we use to create mock observations of O VI absorptions toward background sources, following Fang et al. (2002).The basic idea is to generate random LOS across the simulated region and obtain the temperature, baryon density, and velocity distributions along the LOS.Then O VI ion density can be derived from the metallicity and ionization fraction, from which the optical depth along the LOS and thus the synthetic spectrum can be obtained.Finally, the column densities O VI and Doppler b-parameters for high-and low-velocity clouds can be calculated from the profile of O VI absorption line.

The SMUGGLE galaxy formation model
The SMUGGLE model incorporates physical processes such as gravity, hydrodynamics, gas cooling and heating, star formation, and stellar feedback, and is able to resolve the multiphase gas structure of the ISM.Star particles are formed in cold, dense, and self-gravitating molecular gas reaching a density threshold of n th = 100 cm −3 .The local star formation rate for star-forming gas cells is controlled by the star formation efficiency per free-fall time ϵ ff , i.e., Ṁ * = ϵ ff M gas /τ ff , with M gas the gas mass and τ ff the free-fall time of the gas cell.
The model implements various channels of stellar feedback including photoionization, radiation pressure, energy and momentum injection from stellar winds and from supernovae, which are categorized into two types: (i) SN feedback − injects large amounts of energy and momentum to the ISM.The event number of type II SNe at each time-step is obtained by integrating the Chabrier (2003) initial mass function, and the event rate of type Ia SNe is calculated using a delay time distribution (Vogelsberger et al. 2013).
(ii) Early (pre-SN) feedback − includes radiative feedback and stellar winds.Photoionization and radiation pressure from young massive stars, namely radiative feedback, can impact the ionization state and offer pressure on surrounding gas and thus represent a source of momentum.The energy and momentum injection via stellar winds from young massive OB stars and older populations − asymptotic giant branch (AGB) stars are calculated from the mass loss of the two types of stars, and the former provides another channel of early feedback.
L20 performed a suite of high-resolution, isolated galactic disk simulations using the SMUGGLE model.Incorporated with explicit gas cooling and heating over a wide range of temperatures (10 − 10 8 K), the thermodynamical properties of the multiphase ISM is well studied.The simulations encompass a cubic region of 600 kpc on each side and cover the entire galaxy with the z-axis perpendicular to the galactic disk plane.The initial conditions of the simulation are the same as those of Marinacci et al. (2019).It contains a MW-sized galaxy of 1.6×10 12 M ⊙ , which is composed of a stellar bulge and disk, a gaseous disk, and a dark matter halo, all with masses similar to those of the Milky Way (see Bland-Hawthorn & Gerhard 2016, and references therein).The gaseous disk has an initial mass of ∼ 9 × 10 9 M ⊙ and the density decreases exponentially with a scale length of 6 kpc.The initial setup leads to a gas fraction of roughly 10 percent within a radius of R = 8.5 kpc.The mass resolution reaches 1.4 × 10 3 M ⊙ per gas cell, corresponding to the highest resolution run in Marinacci et al. (2019).Gravitational softening is adaptive for gas cells, with a minimum value of ∼ 3.6 pc.Table 1 lists the main parameters that characterize the initial condition of the simulations.
In L20, we performed six simulations with different subgrid models (feedback channels) and parameters (ϵ ff ).The model variations are summarized in Table 2 and detailed below.
(i) SFE1 − fiducial run in M19 with star formation efficiency of ϵ ff = 0.01 and all stellar feedback channels.
(iv) Nofeed − the same as SFE1 except with no stellar feedback.
(v) Rad − the same as SFE1 except with only early feedback via stellar winds and radiation.
(vi) SN − the same as SFE1 except with only SN feedback.

Mock observations
To generate synthetic observational data, we build a mock galactic coordinate system consistent with that of the Milky Way.Specifically, we place the observer at the location of the Sun, i.e., 8.2 kpc away from the center of the simulated galaxy (Bland-Hawthorn & Gerhard 2016).To avoid selection effects due to a single observer in a specific location, four observers at different off-center locations are situated on the galactic disk, each 8.2 kpc away from the galactic center and 90 • apart from each other (similar to the eight off-center locations in Zheng et al. 2020).We define galactic longitude l and latitude b similar to those of the Galaxy.
For each given set of (l, b) and distance D of the star/quasar to the observer, we trace the LOS across the simulated region using the yt analysis toolkit (http://ytproject.org;Turk et al. 2011) that enables us to obtain gas properties such as temperatures, velocities in the LSR reference frame, and baryon densities along LOS.To directly compare the properties of the warm gas in the simulated galaxy with observations, we convert hydrogen density to O VI density, and then to O VI column density.For a grid of gas temperatures (T ∼ 10 3 − 10 7 K) and hydrogen densities (n H ∼ 10 −8 − 10 6 cm −3 ), O VI in Milky Way-like galaxies v we adopt the Cloudy code (version C17.02 ;Ferland et al. 1998Ferland et al. , 2017) ) to calculate the ionization fraction f OVI (T, n H ) of O VI, taking into account the ultraviolet (UV) background radiation from quasars and galaxies (Haardt & Madau 1996).The number density of O VI can be derived via where Z is the gas metallicity that is set to be the solar, i.e., Z = Z ⊙ , and A O = 4.9 × 10 −4 is the abundance of oxygen (Asplund et al. 2009).Take a random sightline at (l, b) = (0 • , 30 • ) as an example.Figure 1 shows the gas temperature, LSR velocity, baryon number density, and O VI number density along the sightline for the fiducial SFE1 run.The temperature of the gas along the LOS spans a wide range of ∼ 10 3 −10 7 K, and the LSR velocity ranges from roughly −200 to 400 km s −1 .The baryon density exhibits a downward trend as the distance increases, reaching the cosmic mean value of ∼ 2.1 × 10 −7 cm −3 at D ∼ 100 kpc (the blue dashed line).O VI density generally traces the variation of gas temperature for distances ≲ 50 kpc.The reason is that for temperature of ≲ 5 × 10 5 K, the ionization fraction of O VI is a monotonic function of the gas temperature.
The number density of O VI ions along the LOS provides a direct measure to the optical depth (τ ) around O VI 1032 Å line, from which the mock spectrum of a background source (e.g., a star or quasar) can be obtained (e.g., Spitzer 1978;Zhang et al. 1997;Fang et al. 2002).We consider the effects of line broadening and line-center shift caused by both the Hubble velocity and the peculiar velocity of the gas along the LOS (e.g., Fang et al. 2002) We adopt the apparent optical depth (AOD) method to calculate the column density (N ), centroid velocity, and Doppler b-parameter for low-and high-velocity O VI by assuming that the absorption profile is not saturated (e.g., Savage & Sembach 1991;Sembach et al. 2003).For low-velocity O VI, the column density, centroid velocity, and line width are calculated with fixed integration limits of (v − , v + ) = (−100, 100) km s −1 (the magenta band in Figure 2).For high-velocity O VI, the integration limits rely on O VI velocity structures (e.g., Sembach et al. 2003) and are extracted from the cyan bands in Figure 2, which include regions with exp(−τ ) < 0.95 and 100 ≤ |v LSR | < 400 km s −1 .Here 0.95 is chosen somehow arbitrarily to exclude false "absorption features" caused by noise.The two cyan bands indicate two highvelocity components with integration limits of the velocity set by the boundaries of each cyan band.We define |v LSR | < 100 km s −1 as low-velocity clouds (LVCs) and 100 ≤ |v LSR | < 400 km s −1 as HVCs.
Figure 3 displays all-sky map of O VI column density derived from the fiducial SFE1 run, for low-velocity (top), high-velocity (middle), and total gas (bottom), respectively, observed at the Sun's location.As can be seen, low-velocity O VI are widespread in the sky with column densities of log(N/cm −2 ) ≳ 14, and stretch to high galactic latitude of |b| ≳ 60 • .In comparison, highvelocity O VI is generally located near the galactic disk with log(N/cm −2 ) ≳ 15, and gradually decline toward higher galactic latitudes.Such a "disk-like" structure for high-velocity O VI does not appear in the other five runs, which suggests that the spatial distribution of O VI strongly depends on the sub-grid models of stellar feedback.However, it is challenging to detect O VI absorption at low latitudes (e.g., |b| ≲ 25 • ) due to severe ultraviolet extinction for extragalactic objects (e.g., Wakker et al. 2003;Sembach et al. 2003).Therefore, currently it is unavailable to distinguish those models via the observed spatial distribution of O VI.

RESULTS AND DISCUSSION
We derive O VI properties of the simulated galaxy viewed by four internal off-center observers for the fiducial SFE1 run in Section 3.1 − 3.3.The fiducial results are then compared with the other five model variations in Section 3.4 where the impact of sub-grid model/parameter variations are illustrated.Section 3.5 presents results from an external view to compare with observations of external galaxies, which is followed by some caveats in Section 3.6.

The scale height of LVCs
The mock observation of the simulated galaxy in Section 2 provides measurement to the O VI column density along arbitrary LOS across the simulated region.A certain number of sightlines allow us to explore the spatial distribution of the O VI-bearing gas, which can be compared with the observations.Savage & Wakker (2009) collected column densities of O VI as well as other species along 139 LOS toward stars and quasars, and found that low-velocity O VI in the Milky Way is well fitted by an exponentially declined disk model with a scale height of h ∼ 2.6 ± 0.6 kpc.To compare our results with those observations, we generate random LOS according to the following settings.
For each of the four off-center observers, we randomly generate 4×139 LOS across the simulated galaxy toward quasars or stars, and thus the total sightline number is 4×4×139 = 2224.Here 139 is the LOS number collected by Savage & Wakker (2009), among which 109 (30) are toward stars (quasars).We assign the same ratio of numbers for sightlines toward quasars to that toward stars, i.e., 480 (1744) out of the 2224 LOS are toward quasars (stars).The quasars are situated at a distance of 260 kpc (e.g., the virial radius of the Galaxy) and the galactic latitude is randomly drawn at |b| > 20 • since detectable sightlines toward quasars are usually observationally unavailable at |b| ≲ 20 • .The stars are placed in random directions with a distance randomly drawn from 1 − 10 kpc in logarithmic space.This distance range is consistent with the observational data collected by Sav-age & Wakker (2009).O VI column density for LVCs are derived according to the AOD method presented in Section 2.2.We further set a detection limit of O VI column density log(N/cm −2 ) ≥ 13.23 for sightlines toward both quasars and stars (e.g., Savage & Wakker 2009), and there are 1601 O VI absorbers detected along the 2224 LOS, as the gray dots in Figure 4 show.Our results generally agree with the observations of low-velocity O VI (the blue symbols; Savage & Wakker 2009).
To quantify the distribution of low-velocity O VI, we adopt a simple disk model (e.g., Savage et al. 1990Savage et al. , 2000;;Yao & Wang 2005;Savage & Wakker 2009), i.e., the number density declines exponentially away from the mid-plane (or the galactic disk), and the density at a height z below/above the mid-plane can be expressed as where n 0 is the mean density in the mid-plane, h is the scale height of O VI disk.Then the column density (N ) along the LOS can be simply derived from the integration of Equation ( 2), and its projection along z-axis is The formula reveals a monotonic relation between N sin |b| and |z|, for |z| ≲ h; and for |z| ≫ h, N sin |b| eventually approaches a stable value of n 0 h.To perform a reasonable minimum−χ 2 fitting to our mock data (the gray dots in Figure 4), we divide the x-axis into bins.Red circles with error bars show the mean values and standard deviations in each bin.These binned data are then fitted to the disk model, with the best-fitting profile shown as the red solid line.Our bestfitting model has a scale height of h = 2.9 +1.9 −1.2 kpc, well consistent with the MW observations of 2.6 ± 0.5 kpc (Savage & Wakker 2009) considering the 1σ errors.The best-fit parameters and a comparison with observations are listed in Table 3.

Column density − line width relation
A correlation between the column density and the line width of O VI absorbers was first reported by Heckman et al. (2002) and has been found in various environments including the Galactic disk and halo (Jenkins 1978a,b;Savage et al. 2003;Bowen et al. 2008;Lehner et al. 2011;Sarma et al. 2017), HVCs (Sembach et al. 2003), and Magellanic clouds (Howk et al. 2002a,b;Hoopes et al. 2002;Pathak et al. 2011).Collisional processes should be responsible for the linear proportionality between the column density and b-parameter since the column den-sity linearly scales with the gas flow velocity in collisional ionization scenario (see discussions in Heckman et al. 2002;Sembach et al. 2003).Below we investigate the column density − line width relation for low-and high-velocity O VI, which are then compared with the observations.
Figure 5 depicts the column density vs. Doppler parameter distribution for low-velocity O VI.Each data point is obtained from a randomly drawn sightline at |b| > 20 • , from four off-center observers toward quasars/stars, as described in Section 3.1.O VI column densities span log(N/cm −2 ) ∼ 13.2−15.2with a median value of 13.8, and the line width follows a Gaussianlike distribution within b ∼ 13 − 106 km s −1 and peaks at ∼ 47 km s −1 .The median values are listed in the fourth and fifth columns in Table 3.The distribution and median value of the column density are well consistent with the observations (the gray histogram and dotted line; Savage & Wakker 2009).For most of the sightlines, the line width is broader than that caused by thermal motion of O VI ions, which corresponds to b therm ∼ 17.7 km s −1 for gas temperature of 3 × 10 5 K, implying significant non-thermal motions, e.g., inflows, outflows and turbulence.This could be responsible for the distorted or no relation between the column density and line width.Although no correlation between N and b is also expected for photoionized gas, given the high energy (∼ 114 eV) required for ionizing photons, most of O VI ions are implausible to be produced by photoionization except for extreme conditions with a hard radiation field and a very low gas density (see e.g., Sembach et al. 2003).
For HVCs, O VI 1032 Å absorptions have been detected by FUSE at ≥ 3σ confidence levels along 59 out of the 102 sightlines, among which 100 are toward extragalactic objects and two toward halo stars (Sembach et al. 2003).To make a direct comparison with their results, we randomly generate a total of 16 × 59 = 944 sightlines from four off-center observers, where 59 is the number of sightlines with detected O VI absorption reported by Sembach et al. (2003).The background quasars are placed at a distance of 260 kpc, and the galactic latitude is limited to |b| > 20 • .Accounting for the detected high-velocity O VI properties (see Table 1 in Sembach et al. 2003), our mock detections need to satisfy the following conditions: (i) the integration interval v + − v − ≥ 50 km s −1 ; (ii) O VI column density log(N HVC /cm −2 ) ≥ 13.06; and (iii) O VI line width b HVC ≥ 16 km s −1 .This results in 339 detections (the blue filled circles in Figure 6) of high-velocity O VI out of the 944 sightlines, with a detection rate (339/944) lower than that (84/102) given by observations.Our best fit Our results   The gray triangles denotes FUSE observations of HVCs (Sembach et al. 2003).

Our binned results Savage+09
Our column densities and line widths of high-velocity O VI occupy similar parameter space as the observations (Sembach et al. 2003), with b ∼ 16 − 107 km s −1 and log(N/cm −2 ) ∼ 13.1 − 14.8.The median values are also consistent with the real data considering the 1σ uncertainties, i.e., b ∼ 33.0 ± 18.8 vs. 40.0± 13.1 km s −1 , and log(N/cm −2 ) ∼ 13.8 ± 0.4 vs. 14.0 ± 0.3 (see the seventh and eighth columns in Table 3).Unlike the symmetric distribution for LVCs, the line widths for highvelocity O VI peaks at b ≲ 20 km s −1 , suggesting nonthermal motions for high-velocity O VI might be less significant than its low-velocity counterparts.In addition, unlike the random distribution for LVCs, there is a significant positive correlation between the column density and Doppler b-value for high-velocity O VI, still in line with the FUSE observations of the MW (Sembach et al. 2003).Such a correlation may support collisional ionization instead of photoionization as the dominant mechanism for the production of high-velocity O VI (see e.g., Heckman et al. 2002).Moreover, photoionization models underproduce observed OVI column densities by order of magnitude (e.g., Sembach et al. 2003), also landing support to the collisional ionization origin.

Cumulative column density
We note that Zheng et al. (2020) investigated cumulative O VI column densities from an inside-out view of MW analogs selected from the Figuring Out Gas & Galaxies In Enzo (FOGGIE) simulation (Peeples et al. 2019).To make a direct comparison with their results, we adopt the same method as that of Zheng et al.  (Savage et al. 2003), and high-velocity O VI toward quasars (Sembach et al. 2003).The gray dashed line and gray band are the simulation results given by Zheng et al. (2020).
(2020) and randomly generate a total of 1000 LOS with |b| > 20 • for the four off-center observers.For each of the sightline, we calculate the column density as a function of the distance r to the observer by integrating Equation (1) over r.
The median profile as well as the 16th and 84th percentiles are displayed as the blue solid line and band in Figure 7.Despite a systematic offset between our results (blue solid line) and observations of LVCs toward quasars/stars (green crosses and magenta circles; Savage et al. 2003;Savage & Wakker 2009), about half of the observational data points are consistent with our 1sigma uncertainties (blue band).The extrapolation of our results to larger distances, i.e., log(N/cm −2 ) ∼ 14, also agrees with HVC observations toward quasars/stars (Sembach et al. 2003) at r > 100 kpc.Given that those observations only include low-or high-velocity O VI, each set of the observations may represent a lower limit when compared to our results.The large discrepancy at smaller distances (r ≲ 0.3 kpc) could arise from smallscale clumps and cavities in the ISM induced by SN explosions and other feedback processes (Li et al. 2020), despite the small number statistics.
In contrast, Zheng et al. (2020) underproduced O VI in the halos by 1 − 2 orders of magnitude in the column density (the gray dashed line in Figure 7), comparing to our results and to the observations.The reason, as they have pointed out, could be that their simulated dark matter halos are smaller than the real case (Bland-Hawthorn & Gerhard 2016), and/or that they only consider the thermal feedback from SNe, which is unable to expel enough metals into the ISM/CGM.The consideration of the full suite of feedback processes (e.g., stellar winds, radiative feedback, and SN explosions) by the SMUGGLE model and by the fiducial run of L20's simulation could be responsible for our agreement with the real data.

Other simulation models
Results presented in Section 3.1 to 3.3 are derived from the fiducial SFE1 run of L20's simulation.To explore how O VI absorption features are affected by different subgrid models, we consider the other five variations listed in Table 2 for comparison.
Similar to Figure 4 for the SFE1 run, Figure 8 shows log(N LVC sin |b|) versus log |z| for the other five models.While the scale height for low-velocity O VI derived from the fiducial SFE1 run is comparable to the observations (Savage & Wakker 2009), the runs with higher star formation efficiency, e.g., SFE10 and SFE100 runs, result in smaller scale height for O VI-bearing gas.The scale height does not always decrease with increasing star formation efficiency, which is attributed to the degeneracy between the scale height (h) and mid-plane density (n 0 ).Meanwhile, the projected column density log(n 0 h) at |z| ≫ h decreases slightly as the star formation activity weakens.The reason is that early feedback (e.g., stellar winds, radiation pressure) that is enhanced by intense star formation blows gas and metals away.Higher star formation efficiency also leads to more SN events at a given time-step, and SN feedback could also play a role.The Rad run considering only radiative feedback and stellar winds results in much lower scale height and projected column density, compared to the SFE1 run with the full suite of stellar feedback, revealing that SN feedback plays an important role in reproducing the observed spatial distribution of low-velocity O VI, i.e., SN energy and momentum injections collisionally ionize more O VI and push the warm gas further out of the galactic disk.Indeed, the fitting result for the SN run is in nice agreement with the observations.In contrast, the run without feedback − "Nofeed", gives an overall lower density and a low scale height for low-velocity O VI.
Table 3 (the second to fourth columns) lists the bestfit parameters for the six runs and the values derived from observations (Savage & Wakker 2009), which, for a more clear view, are compared in Figure 9.The gray band and dashed line are the observational constraints provided by Savage & Wakker (2009).See also the red lines in Figure 4 (the SFE1 run) and Figure 8 (the other five runs) for the fitting results and Table 3 for the best-fit parameters.
none of the five runs exhibit obvious correlations.The distribution and median value of O VI column density for the SN run agree excellent well with the observations.In contrast, the runs lack of SN feedback (e.g., Nofeed and Rad runs) underproduce O VI, with median column densities ∼ 0.6 dex lower.The critical impact of SN feedback is once again highlighted.
For high-velocity O VI, the column density − line width relation are displayed in Figure 11 for the other four runs, as compared to that of the SFE1 run shown in Figure 6.The median values of the column densities and line widths for different runs are listed in the seventh and eighth columns of Table 3.The Nofeed run is not displayed because no high-velocity O VI components are detected, indicating the necessity of feedback processes to accelerate O VI particles.While the column densities of high-velocity O VI derived from different runs generally agree with the observations (Sembach et al. 2003) accounting for the uncertainties, the median values and distributions of the Doppler parameter support the SFE1 and SN runs, both including SN feedback.
The log(n 0 h) values for different model variations listed in Table 3 represent the simulated galaxy at the "present" time when the simulation is terminated.In fact, for each of the snapshot of the simulation, we can similarly obtain its log(n 0 h) value.In Figure 12, we present the evolution of log(n 0 h) across the simulation time for the six runs.As can be seen by comparing the SFE1, SFE10, and SFE100 runs, a larger star formation efficiency results in a downward tendency of log(n 0 h) over time.This could be attributed partly to the fast conversion of cold gas to stars and thus less oxygen is left for O VI production via heating.Meanwhile, the stellar winds from young massive stars have an important impact on the ISM gas (Lamers & Cassinelli 1999;Muijres et al. 2012), e.g, dispersing the gas and impeding the generation of O VI ions via SN feedback heating.Consequently, the "present" value of log(n 0 h) and its 2σ confidence region for the SFE1 run marginally agrees with the observations (Savage et al. 2003;Bowen et al. 2008;Savage & Wakker 2009), yet the SFE100 run deviates further.For the SFE1 and SN runs, the simulation data are available only for runtime within 0.8 and 0.5 Gyr, respectively.Based on the currently available data, the "present" log(n 0 h) value for these two runs are in better agreement of with the observations than the other four model variations.
To summarize, comparison of different runs in L20's simulation with the observations of low-and highvelocity O VI favors the SFE1 and SN runs, suggesting that SN feedback is required to reproduce the O VI observations, and meanwhile early feedback associated with star formation activities should be moderate (not too strong), e.g., with star formation efficiency of ϵ ff ∼ 0.01.

Comparison with external galaxies
Results in Section 3.1 -3.4 are viewed from off-center observers inside the simulated galaxy.Here we present results for the SFE1 run viewed from an external observer and compare with observations of external galaxies.
The left panel of Figure 13 shows the face-on view of O VI column density map (on xy−plane).For each grid of coordinates (x, y), the column density is obtained by integrating O VI number density in Equation ( 1) along the z-axis with path length of 600 kpc, i.e., the size of the simulation box.The white dashed line denotes the virial radius of 260 kpc.The column density peaks at the center with log(N/cm −2 ) ∼ 15.3 and gradually declines toward outer region, approaching a background value of log(N/cm −2 ) ∼ 4.4.Besides that, there is tentative evidence for structures spanning tens of kpc.
To make a direct comparison with observations of external galaxies, we plot the column density versus the impact parameter in the right panel of Figure 13.To achieve that, we generate random sightlines for both face-on and edge-on views of the simulated galaxy.The blue solid line shows the median column density, and the blue band shows the range of 5th to 95th percentiles.Our results are consistent with the observations of suband super-L * galaxies (Prochaska et al. 2011) for impact parameter ≲ 50 kpc.Beyond that, the column density sharply declines and drops below the observational values.This happens as expected because the simulation performed by L20 as well as the SMUGGLE galaxy formation model is for an isolated galaxy without gas supply from the IGM, which is also the shortcoming of this study.Indeed, the Galactic halo density (∼ 10 −4 cm −3 ) suggested by observations of the Magellanic Stream (Weiner & Williams 1996) is more than one order of magnitude higher than our results (≲ 10 −5 cm −3 ; the panel (c) of Figure 1) at a radius of ∼ 50 kpc . 3.6.Caveats

Isolated galaxy simulation
Our results are based on L20's simulations for an isolated galaxy without gas fueling from the IGM and interactions with companion galaxies.This could lead to an underestimation of O VI column density at outer regions, e.g., r ≳ 100 kpc (see Figure 13).In addition, our high-velocity O VI clouds can only be produced via the galactic fountain, i.e., triggered by stellar feedback (e.g.,    Prochaska et al. (2011).Shapiro & Field 1976;Bregman 1980;Fraternali & Binney 2006).If other mechanisms such as accretion from the IGM (e.g., Fraternali et al. 2015;Kereš & Hernquist 2009) and materials stripped or ejected from satellites (e.g., Putman 2004;Herenz et al. 2013), are also responsible for the formation of high-velocity O VI, our simulation (Table 3) may underproduce high-velocity O VI content and distort its spatial distribution.

The metallicity
Our results in this work are obtained under the assumption of solar metallicity for the gas when converting the number density of hydrogen to that of O VI in Equation ( 1).Constant metallicity is often assumed for simplicity despite the fact that the metallicity could differ by orders of magnitudes for different regions of the galaxy (e.g., Gutcke et al. 2017;De Cia et al. 2021).Alternatively, we quantify the effect of metallicity on the scale height evolution of low-velocity O VI for the SFE10 run in Figure 14, by setting three constant metallicities of 1 Z ⊙ , 3 Z ⊙ , and 5 Z ⊙ .As expected, a higher metallicity results in a larger scale height of O VI, which differs by a factor of ≲ 2 for 1 Z ⊙ and 5 Z ⊙ cases, comparable to the variations of the scale height across the simulation time of ∼ 1 Gyr.While the SFE10 run is ruled out under the assumption of solar metallicity when compared to the observations (Savage et al. 2003;Bowen et al. 2008;Savage & Wakker 2009), higher metallicity of 5 Z ⊙ makes the SFE10 run's results (h = 2.3 +2.2 −1.2 kpc) well consistent with the observations considering the errors.This indicates that to some extent, a higher metallicity can compensate for lower O VI content caused by strong early feedback (e.g, stellar winds, radiation pressure) launched by short-lived massive stars.

The UV background and other ionizing sources
Our results on O VI properties of the simulated galaxies are based on Equation (1), where the ionization fraction of O VI is derived via Cloudy (Ferland et al. 2017) modeling by taking into account extragalactic UV background radiation (Haardt & Madau 1996).While such UV background is typically applied to intergalactic medium regions (e.g., Fang & Bryan 2001), there are alternative versions of the UV background in the literature and other potential contribution of ionizing sources, e.g., stellar radiation within the galaxy and cosmic ray heating (Werk et al. 2014).
The spectral shape of UV background has been shown to affect oxygen abundance (Aguirre et al. 2008) and statistics of O VI absorbers in the IGM (Oppenheimer & Davé 2009).A comparison of various UV background have been presented in the Figure 1  adopted.Energy of ∼ 114 eV required for photoionizing O V corresponds to the high-energy tail of the spectral energy distribution (SED) of background radiation field.Consequently, most O VI could be produced from collisional ionization at temperatures of ∼ 3 × 10 5 K rather than from photoionization at lower temperatures (e.g., Cox 2005).Moreover, the flux difference at the ionizing energy is at most ∼ 0.5 dex for various frequently used UV background and should not make much difference.
Photoionization is considered as a channel of radiative stellar feedback in the framework of the SMUG-GLE model, and O VI distribution that is controlled by the ionization fraction f O VI (T, n H ) in Equation ( 1) is affected by feedback processes in terms of heating (increasing the temperature T ) and/or blowing gas away (decreasing hydrogen density n H ).However, the stellar radiation inside the galaxy is not directly considered in the Cloudy modeling.The contribution of a starburst galaxy to the total ionizing photons is evaluated in the Figure 13 of Werk et al. (2014).The SED of the radiation field with the contribution of the starburst galaxy with a star formation rate (SFR) of 1 M ⊙ yr −1 at a dis-O VI in Milky Way-like galaxies xvii tance of d = 72 kpc has similarly flat slope toward higher energies (E ≳ 70 eV), and deviates ∼ 0.1 dex from the Haardt & Madau (2001) UV background.While the assumed SFR is typically true for our simulated MW-like galaxy (Li et al. 2020), the distance of O VI clouds to the star forming region spans a wide range across the halo (0 − 260 kpc).Figure 8 of Fox et al. (2005) compares the UV background radiation with the radiation from the MW at different locations within the Galaxy, and reveals that the radiation field from the Galaxy is similar to the Haardt & Madau (1996) UV background at E ∼ 114 eV for a distance of d ∼ 20 − 30 kpc.The sharp decrease of the radiation flux at 54 eV arising from He II edge in hot stars indicates that high-energy photons of E > 114 eV are very limited.Detailed Cloudy modeling suggests that photoionization is negligible for the production of O VI despite its dependence on the radiation field adopted.
Cosmic ray heating could serve as a crucial supplementary source of ionization and heating within the Galactic virial radius (Wiener et al. 2013).For gas densities ≳ 10 −2 cm −3 , the cosmic ray background (CRB) can dominate over photoelectric heating for gas accounting for a weaker dependence on the gas density for the CRB heating.Therefore, CRB could significantly enhance the density of O VI in low density regions, although the precise number is challenging to determine because of the poorly constrained local CRB (see discussions in Werk et al. 2014).

CONCLUSIONS
We study O VI properties in MW-like galaxies by analyzing the suites of simulation performed by L20 in the framework of SMUGGLE galaxy formation model.We find that the SMUGGLE model is capable of producing consistent global properties of Galactic warm gas traced by O VI.In addition, mechanical stellar feedback is shown to have a crucial impact on the spatial distribution and kinematics of O VI absorbers.Particularly, SN feedback is necessarily required, and early feedback associated with star formation activities needs to be moderate to reproduce O VI observations.Our main findings are detailed as follows.
(i) Low-velocity O VI distribution is well described by an exponentially declining disk with a scale height of 2.9 +1.9 −1.2 kpc and log(n 0 h/cm −2 ) = 13.79 ± 0.16 for the fiducial SFE1 run (with full suites of feedback processes), generally consistent with the observations.The SN run (with SN feedback only) results in a scale height well consistent with observations as well.Other runs turning off SN feed-back or with higher star formation efficiencies lead to smaller values for the scale height.(v) We cannot reproduce observations of column density profile for external galaxies due to the lack of accretion in our simulations, suggesting that accretion is an important part of galaxy evolution modeling.
Overall, the observed Galactic O VI properties can be reasonably reproduced with simulations of isolated MWlike discs based on the SMUGGLE model with novel treatment of ISM and stellar feedback, in complement to L20's findings of its success in producing realistic cold ISM.A test of its ability in reproducing hotter Galactic gas traced by highly ionized metal species such as O VII and O VIII is deferred to a future work.One shortcoming of the SMUGGLE model could be the lack of cosmological gas accretion.The next generation of the SMUGGLE model intends to involve cosmological simulations with zoom-in of individual objects, and will serve as a powerful tool for predicting galactic structure, outflows, and CGM properties.
. The original spectrum is convolved with the FUSE line spread function to account for the instrumental broadening (b inst ), i.e., b = b 2 therm + b 2 inst , with b therm the thermal broadening and b inst ∼ 12 − 15 km s −1 (Moos et al. 2000; Sembach et al. 2003).Gaussian noise is further considered with the mean of 0 and the standard deviation of 0.01.The final synthetic spectrum, or the transmission exp(−τ ), is shown as the black solid line in Figure 2, for a background source placed at a distance of 260 kpc and in the direction of (l, b) = (0 • , 30 • ), the same LOS as in Figure 1.Multiple absorption components can be seen with |v LSR | ∼ 30 − 300 km s −1 , consistent with a wide spread of velocity for gas along the LOS shown in Figure 1(b).

Figure 1 .
Figure 1.Panels from top to bottom respectively represent the gas temperature (a), velocity in the LSR reference frame (b), baryon number density (c), and O VI number density (d) along the sightline (l, b) = (0 • , 30 • ) in the galactic coordinate to an observation at the Sun's location.The blue dashed line in panel (c) denotes the mean baryon density of the universe (see the text for details).

Figure 2 .
Figure 2.An example of synthetic O VI absorption spectrum (the black solid line) and identification of low-and high-velocity O VI components.Velocities between −100 and 100 km s −1 are identified as low-velocity component by the magenta band marks.High-velocity components are shown as the cyan bands, with velocities exceeding 100 km s −1 and exp(−τ ) < 0.95.The gray dashed line marks exp(−τ ) = 0.95.
.00 ± 13.14 13.97 ± 0.33 Note-Columns from left to right represent: (1) the six model variations of the simulation presented in L20 (the second to seventh rows), or the observations of LVCs and HVCs by Savage & Wakker (2009) and Sembach et al. (2003), respectively (the last two rows); (2)-(4) the best-fit parameters of the disk model for low-velocity O VI; (5)-(6) the median column density and line width for low-velocity O VI; (7)-(8) the median column density and line width for high-velocity O VI.

Figure 4 .
Figure 4. Projected O VI column density along z-axis versus height to the galactic plane for LVCs in the SFE1 run.The gray dots represent our mock results, and the red diamonds are our binned results with 1σ error bars.The red solid line is the best-fit model to the red diamonds, and the dashed lines enclose 1σ confidence region accounting for the errors of the scale height.The blue circles are results revealed by observations (Savage & Wakker 2009).

FractionFigure 5 .
Figure 5. Column density and line width distribution of O VI-bearing LVCs in the fiducial SFE1 run.The blue histograms give the probability distribution of the column density (right) and line width (top), with the median values denoted by the blue dashed lines, as compared with the observations shown as the gray histogram and gray dotted line (Savage & Wakker 2009).

FractionFigure 6 .
Figure 6.Similar to Figure 5 but for high-velocity O VI.The gray triangles denotes FUSE observations of HVCs(Sembach et al. 2003).

xFigure 7 .
Figure 7. Profile of the cumulative O VI column density for the SFE1 run, as represented by the blue solid line and the blue band.Purple open circles, green pluses, and gray triangles are observations of low-velocity O VI toward stars (Savage & Wakker 2009), toward quasars(Savage et al. 2003), and high-velocity O VI toward quasars(Sembach et al. 2003).The gray dashed line and gray band are the simulation results given byZheng et al. (2020).

Figure 8 .Figure 9 .
Figure 8. Distribution of the projected column densities along z-axis vs. heights above the galactic disk for low-velocity O VI, for the other five runs as labeled in each panel.Legends are similar to those in Figure 4.

FractionFigure 10 .
Figure 10.The column density and line width distribution for low-velocity O VI, for the other five runs as labeled on the top left of each panel.Legends are similar to Figure 5.

FractionFigure 11 .Figure 12 .
Figure 11.Column densities vs. line width distribution for high-velocity O VI for the other four model variations.Legends are similar to Figure 6.The Nofeed run is not displayed because no high-velocity O VI absorptions are detected.

Figure 13 .
Figure 13.Left: Face-on projection of O VI column density (along z-axis) for the SFE1 run, as the color bar denotes.The white dashed line marks the galaxy's virial radius of 260 kpc.Right: O VI column density as a function of the impact parameter derived from random face-on and edge-on sightlines.The blue solid line is our median result, and the blue area denotes the 5th − 95th percentiles range.The black filled circles and open triangles are the observational results for sub-and super-L * galaxies given by Prochaska et al. (2011).

Figure 14 .
Figure 14.The evolution of the exponential scale height h of low-velocity O VI as a function of time for the SFE10 run with metallicities of 1 Z⊙ (black solid line), 3 Z⊙ (blue dashed line), and 5 Z⊙ (green dotted line), respectively.The lines represent the median values for random sightlines described in Section 3.1, and the regions with corresponding colors represent 1σ uncertainties.The symbols with error bars are observational results reported by Savage et al. (2003), Bowen et al. (2008), and Savage & Wakker (2009), as labeled.In particular, Bowen et al. (2008) provides the scale heights of low-velocity O VI for the northern (b > 20 • ) and southern (b < −20 • ) hemisphere of the MW, respectively.

(
ii) For the SFE1 run, the column density of lowvelocity O VI is distributed in the range of log(N/cm −2 ) ∼ 13.2 − 15.2 with a median value of ∼ 13.8, consistent with observations within 1σ uncertainties.The line width of low-velocity O VI follows a Gaussian-like distribution over b ∼ 13−106 km s −1 with a median value of 47.4 km s −1 .No correlations are found between the column density and line width of low-velocity O VI for all of the model variations.(iii) For high-velocity O VI in the SFE1 run, the column density spans log(N/cm −2 ) ∼ 13.1 − 14.8 with a median of ∼ 13.8, and line width covers b ∼ 16 − 107 km s −1 with a median of ∼ 33 km s −1 .A positive correlation are found between the column density and line width of high-velocity O VI, supporting collisional ionization as the dominant mechanism for the production of high-velocity O VI.No high-velocity O VI clouds are found in the run turning off all channels of stellar feedback.(iv)The profile of cumulative O VI column density generally agrees with observations for the SFE1 run.The evolution of log(n 0 h) as a function of simulation time also supports the SFE1 and SN runs when comparing to observations.

Table 1 .
Initial setup of the simulation performed in L20.

Table 3 .
Properties low-and high-velocity O VI clouds for the six runs in L20's simulation and comparison with the observations.