GARROTXA Cosmological Simulations of Milky Way-sized Galaxies: General Properties, Hot Gas Distribution, and Missing Baryons

We introduce a new set of simulations of Milky Way-sized galaxies using the AMR code ART + hydrodynamics in a $\Lambda$CDM cosmogony. The simulation series is named GARROTXA and follow the formation of a halo/galaxy from z~$=$~60 to z~$=$~0. The final virial mass of the system is $\sim$7.4$\times$10$^{11}$M$_{\odot}$. Our results are as follows: (a) contrary to many previous studies, the circular velocity curve shows no central peak and overall agrees with recent MW observations. (b) Other quantities, such as M$\_{*}$(6$\times$10$^{10}$M$_\odot$) and R$_d$ (2.56 kpc), fall well inside the observational MW range. (c) We measure the disk-to-total ratio kinematically and find that D/T=0.42. (d) The cold gas fraction and star formation rate (SFR) at z=0, on the other hand, fall short from the values estimated for the Milky Way. As a first scientific exploitation of the simulation series, we study the spatial distribution of the hot X-ray luminous gas. We have found that most of this X-ray emitting gas is in a halo-like distribution accounting for an important fraction but not all of the missing baryons. An important amount of hot gas is also present in filaments. In all our models there is not a massive disk-like hot gas distribution dominating the column density. Our analysis of hot gas mock observations reveals that the homogeneity assumption leads to an overestimation of the total mass by factors 3 to 5 or to an underestimation by factors $0.7-0.1$, depending on the used observational method. Finally, we confirm a clear correlation between the total hot gas mass and the dark matter halo mass of galactic systems.


INTRODUCTION
One of the most challenging problems that hydrodynamical simulations has faced since the pioneering works of Evrard (1988), Hernquist & Katz (1989), Cen et al. (1990), Navarro & Benz (1991), and Navarro & White (1993) has been to produce, within the standard Λ cold dark matter (ΛCDM) hierarchical structure formation scenario, systems that look like real disk galaxies. That is, galaxies with extended disks and present-day star formation rates (SFRs), disk-to-bulge ratios, and gas and baryonic fractions that agree with observations. It has long been known that disks should form when gas cools and condenses within dark matter (DM) halos (White & Rees 1978;Fall & Efstathiou 1980;Mo et al. 1998) conserving angular momentum obtained through external torques from neighboring structures (Hoyle 1953;Peebles 1969). Yet, the first halo/galaxy hydrodynamical simulations (Navarro & Benz 1991;Navarro & Steinmetz 2000) inevitably ended up with small and compact disks and massive spheroids that dominated the mass. Two effects were found to be the responsible for this result. First, artificial losses of angular momentum caused by insufficient resolution and other numerical effects (Abadi et al. 2003;Okamoto et al. 2003;Governato et al. 2004;Kaufmann et al. 2007). Second, the dynamical friction that transfers the orbital angular momentum of merging substructures to the outer halo and causes cold baryons (stars plus gas) to sink to the center of the proto-galaxy (e.g., Hernquist & Mihos 1995). Once in the center, baryons that are still in the form of gas are quickly converted to stars. These newborn stars will join the stellar component and both will form an old spheroid, leaving no gas for the formation of new disk stars at later times (e.g., Maller & Dekel 2002). Later, the increase in resolution in the numerical simulations led to several authors (e.g., Governato et al. 2004;Robertson et al. 2004;Okamoto et al. 2005) succeeding in obtaining more realistic disk galaxies. This better resolution was achieved using the so called zoom-in technique. This technique allows simulators to obtain a high resolution by choosing a small region, almost always a sphere centered on a specific halo, from a relatively large box, and rerunning it with higher resolution (e.g., Klypin et al. 2001). To obtain a "Milky Way (MW)-like system the authors selected halos with masses similar to the one estimated for the MW (∼10 12 M e ) and halos that have not suffered a major merger since z ∼ 1.5 (e.g., Governato et al. 2004;Robertson et al. 2004;Okamoto et al. 2005;Scannapieco et al. 2009;Guedes et al. 2011;Moster et al. 2014;Vogelsberger et al. 2014).
The recent relative success in forming more realistic disk galaxies has not been limited to resolution. The implementation of new subgrid physics such as efficient supernova (SN) feedback (Stinson et al. 2006) and new star formation (SF) recipes have shown to be important (Gnedin et al. 2009). The more effective stellar feedback acts by removing low-angularmomentum material from the central part of (proto)galaxies (e.g., Brook et al. 2012) while a SF based on the molecular hydrogen abundance is more realistic. These improvements make simulated galaxies follow observed scaling relations, like the Tully-Fisher relation, for the first time (Governato et al. 2007) and thus relax the tension that existed between observations and ΛCDM predictions (Governato et al. 2010).
Until recently, simulated disk galaxies with rotation curves similar to that of the MW were scarce. Most works tended to obtain systems with rotation curves that were too peaked and declining (e.g., Scannapieco et al. 2012;Hummels & Bryan 2012). Agertz et al. (2011) showed that by fine tuning the values of the SF efficiency and the gas density SF threshold parameters it was possible to avoid the formation of galaxies with such peaked rotation curves. In more recent high resolution simulations, including ours, this is no longer a problem (Aumer et al. 2013;Marinacci et al. 2014;Agertz & Kravtsov 2015;Keller et al. 2015;Mollitor et al. 2015;Murante et al. 2015;Santos-Santos et al. 2016), in part due to an added "early-feedback" (Stinson et al. 2013) to the thermal SN one 5 , the implementation of an efficient kinetic feedback, and/or the inclusion of active galactic nucleus (AGN) feedback. The early-feedback is attributed to a feedback that is present from the very beginning, before the first SN explodes, and it may include radiation pressure (Krumholz 2015) and photoheating of the gas due to the ionizing radiation of massive stars (Trujillo-Gomez et al. 2015), stellar winds, etc. Currently the need for such effects is widely acknowledged, but their correct implementation is still a challenge (González-Samaniego et al. 2014b;Oman et al. 2015).
Obtaining systems that resemble the MW opens several possibilities for the study of galaxy formation and evolution. For example, it is interesting to study the origin of the galactic halo coronal gas, its properties, and how they can account for a fraction of the missing baryons. The lack of baryons in the universe was reported by Fukugita et al. (1998) when they found that the cosmological baryon fraction inferred from Big Bang nucleosynthesis is much higher than the one obtained by counting baryons at a redshift z = 0. This problem was confirmed when a more accurate value for the cosmological baryon fraction was obtained from the high precision data of WMAP (Dunkley et al. 2009) and Planck (Planck Collaboration et al. 2014). In fact, the cosmic ratio Ω b /Ω M is 3−10 times larger than the one observed for galaxies. To solve this missing baryons problem, following the natural prediction from structure formation models (White & Rees 1978;White & Frenk 1991), authors like Cen & Ostriker (1999), Davé et al. (2001), and Bregman (2007) proposed that galactic winds, SN feedback, or strong AGN winds ejected part of the galactic baryons to the DM halo and circumgalactic medium (CGM) as hot gas. Others authors like Mo & Mao (2004) also proposed that part of the gas never collapsed into the DM halos as it was previously heated by Population III SNe. From the point of view of large volume hydrodynamical galaxy formation simulations, some attempts have been made to investigate the presence and detectability of hot-gas halo corona. Toft et al. (2002) presented one of the first studies in this direction, however, due to inefficient feedback, the amount of hot gas in their simulated halo was too small. The global properties of hot-gas halo coronae in MW-sized simulations have also been studied (e.g., Guedes et al. 2011;Mollitor et al. 2015). More recently, models that account for SN feedback and implement new SF processes succeeded in obtaining more realistic results when compared with observations (Crain et al. 2010(Crain et al. , 2013. From an observational point of view, hot-gas halo coronae were proposed to be detected by their X-ray emission. The first detections where performed by Forman et al. (1985) in external galaxies and later several authors obtained better detections in other systems and studied them in more detail (e.g., Mathews & Brighenti 2003;Li et al. 2008;Anderson & Bregman 2010;Bogdán & Gilfanov 2011). More recently, using Chandra, XMM-Newton, FUSE, and other instruments (e.g., Hagihara et al. 2010Hagihara et al. , 2011Miller & Bregman 2013) a large reservoir of this hot gas was detected surrounding our galaxy (e.g., ∼6 × 10 10 M e ; Gupta et al. 2012) and it was proposed that this could account for a fraction of the so called missing baryons mass. These detections were achieved through the analysis of X-ray absorption lines from O VII and O VIII which only exist in environments with temperatures between 10 6 and 10 7 K; the MW halo has a temperature of about logT(K) = 6.1 −6.4 (Yao & Wang 2007;Hagihara et al. 2010). The problem of using this technique is that these X-ray absorption lines can only be observed in the direction of extragalactic luminous sources (quasi-stellar objects (QSOs), AGNs, etc.) or galactic X-ray emitters such as X-ray binaries. This observational strategy provides limited information about the hot-gas distribution and position inside the galactic halo or in the intergalactic medium. Using this limited information, several recent works (Gupta et al. 2012;Miller & Bregman 2013;Gupta et al. 2014) obtained a value for the total hot-gas mass embedded in the galactic halo that accounts for about 10%-50% of missing baryons. To obtain the total hot-gas mass from the low number of available observations, the authors had to assume a simple hot-gas density profile. It is important to be aware that this assumption can lead to biased total hot-gas mass values. Recently Faerman et al. (2016) used an analytic phenomenological model to study the warm-hot gaseous corona distribution in galaxies. They used MW hot-gas X-ray and UV observations as input for their model and found that the warm-hot halo coronae may contain a large reservoir of gas. They concluded that if metallicity is about 0.5 of solar metallicity, it can account for an important fraction of the missing baryons.
Here, we introduce a new set of MW-sized simulations with a high number of DM particles and cells (∼7 × 10 6 ). The spatial resolution, the side of the cell in the maximum level of refinement, is 109 pc. These are highly resolved simulations even by today's standards. The simulated galaxy ends up with an extended disk and an unpeaked slowly decreasing rotation curve. It is important to mention that the subgrid physics used in the present work was shown to produce realistic low-mass galaxies, at least as far as the rotation curves were concerned (e.g., Colín et al. 2013). In the GARROTXA (Galaxy High Resolution Runs in a Cosmological Context using ART) series of simulations the temperature reached by the gas in the cells, where stellar particles are born, is actually higher than the one obtained in, for example, the simulations by González-Samaniego et al. (2014b) because the SF efficiency factor is higher (0.65 versus 0.5). This, along with the fact that here the gas density threshold is lower, 1 cm −3 instead of 6 cm −3 , makes the stellar feedback in the GARROTXA simulations stronger.
Using our simulations we have also studied how the hot-gas component embedded in the DM halo can account for part of the missing baryons in galaxies. A novel aspect of the present paper is the determination, in our GARROTXA series of simulations, of the hot halo gas distribution in a full-sky view, which offers us an opportunity to detect observational biases in the determination of its mass when only a small number of lineof sight observations are used.
This paper is laid out as follows. In Section 2 we present the code and initial conditions we have used. In Section 3 we introduce the general properties of our MW-sized simulations, and in Section 4 we present the study of their hot-gas component. We summarize our conclusions in Section 5.

The Code
The numerical simulations we introduce in this work, GARROTXA, have been computed using the Eulerian hydrodynamics + N-body ART code (Kravtsov et al. 1997;Kravtsov 2003). This is the Fortran version of the code introduced in Colín et al. (2010) in the context of the formation of low-mass galaxies. This version differs from the one used by Ceverino & Klypin (2009), among other things, in the feedback recipe. Unlike them, we deposit all the energy coming from stellar winds and SNe in the gas immediately after the birth of the stellar particle. This sudden injection of energy is able to raise the temperature of the gas in the cell to several 10 7 K, high enough to make the cooling time larger than the crossing time and thus avoiding most of the overcooling. 6 Here we have used it to obtain high resolution MW-sized galaxies.
The code is based on the adaptive mesh refinement technique, which allows one to increase the resolution selectively in a specified region of interest around a selected DM halo. The physical processes included in this code are the cooling of the gas and its subsequent conversion into stars, the thermal stellar feedback, the self-consistent advection of metals, and a UV heating background source. The total cooling and heating rates used incorporate Compton heating/cooling, atomic, molecular hydrogen, and metal-line cooling, and UV heating from cosmological background radiation (Haardt & Madau 1996), and are tabulated for a temperature range of 10 2 K < T < 10 9 K and a grid of densities, metallicities (from log (Z) = −3.0 to log(Z) = 1.0, where Z is in solar units), and redshifts, using the CLOUDY code (Ferland et al. 1998, version 96b4).
The SF is modeled as taking place in the coldest and densest collapsed regions, defined by T < T SF and n g > n SF , where T and n g are the temperature and number density of the gas, respectively, and T SF and n SF are the temperature and density SF threshold, respectively. A stellar particle of mass m is placed in all grid cells where these conditions are simultaneously satisfied, and this mass is removed from the gas mass in the cell. The particle subsequently follows N-body dynamics. No other criteria are imposed. In most of the simulations presented here, the stellar particle mass, m, is calculated by assuming that a given fraction (SF local efficiency factor ò SF ) of the cell gas mass, m g , is converted into stars; that is, m = ò SF m g , where ò SF is treated as a free parameter. In the MW-sized models presented here we have used ò SF = 0.65, T SF = 9000 K, and n SF = 1 cm −3 , where the latter is the density threshold in hydrogen atoms per cubic centimeter. As shown in Colín et al. (2010), these values successfully reproduce realistic low-mass galaxies, at least as far the circular velocity is concerned. In Colín et al. (2010) the authors also found that the reduction of ò SF makes the temperature reached by the gas in the starforming cell too low and, as a consequence, the cooling became too efficient. They also showed that the strength of the stellar feedback recipe used here also depends on the value of n SF : the lower its value, the stronger the effect of the feedback. We tested this and found that most of the overcooling is avoided when n SF is around 1 cm −3 and ò SF = 0.65. Finally, the value T SF is almost irrelevant as long it is set below or close to 10 4 K because this is always achieved when the density threshold is reached. Incidentally, the Jeans length of an isothermal gas, with a density of n H = 3 cm −3 and a temperature of 3000 K, is 638 pc, which is more than four times the length of the cell of the current highest refinement level 7 (109 pc). This value has been obtained using the standard perturbative derivation of Jeans criteria (Kolb & Turner 1990;Shu 1991;Binney & Tremaine 2008), which delivers the expression shown in Equation (1) for the Jeans length, where c s is the sound speed and ρ 0 the mass density of the star-forming gas (cold and dense gas). In Appendix B we show a histogram of the Jeans length, in units of the corresponding cell, for the cold gas, In this work we have found that these values also successfully reproduce a galaxy rotation curve similar to that of the MW at z = 0 for galaxies with a mass similar to that of the MW. Since stellar particle masses are much more massive than the mass of a star, typically 10 4 -10 5 M e , once formed, each stellar particle is considered as a single stellar population, within which the individual stellar masses are distributed according to the Miller & Scalo (1979) initial mass function (IMF). Stellar particles eject metals and thermal energy through stellar winds and Type II and Ia SN explosions. Each star more massive than 8 M e is assumed to instantaneously dump 2 × 10 51 erg of thermal energy into the interstellar matter (ISM); 10 51 erg comes from the stellar wind and the other 10 51 erg from the SN explosion. Moreover, the star is assumed to eject 1.3 M e of metals. For the assumed Miller & Scalo (1979) IMF, a stellar particle of 10 5 M e produces 749 SNe II. For a more detailed discussion of the processes implemented in the code, see Kravtsov (2003) and Kravtsov et al. (2005). Stellar particles dump energy in the form of heat to the cells in which they are born. If the subgrid physics are not properly simulated, most of this thermal energy inside the cell is radiated away. Thus, to allow for outflows, it is common to adopt the strategy of turning off the cooling during a time t off in the cells where young stellar particles (age < t off ) are located (see Colín et al. 2010). This mechanism, along with a relatively high value of ò SF , allows the gas to expand and move away from the starforming region. As t off can be linked to the crossing time in the cell at the finest grid, we could see this parameter as depending on resolution in the sense that the higher the resolution, the smaller its value. In the simulations presented here the cooling is stopped for 40 Myr. Although it has been demonstrated that for spatial resolutions similar to the ones in our models this process is unnecessary when simulating dwarf galactic systems (González-Samaniego et al. 2014b), a more detailed study is necessary to ensure that the same occurs when simulating MWsized galaxies.

Simulation Technique and Halo Selection
The simulations presented here have been run in a ΛCDM cosmology with Ω 0 = 0.3, Ω Λ = 0.7, Ω b = 0.045, and h = 0.7. The CDM power spectrum was taken from Klypin & Holtzman (1997) and it was normalized to σ 8 = 0.8, where σ 8 is the rms amplitude of mass fluctuations in 8 Mpc h −1 spheres. To maximize resolution efficiency, we used the zoom-in technique. First, a low-resolution cosmological simulation with only DM particles was run, and then regions of interest (DM halos) were picked up to be re-simulated with high resolution and with the physics of the gas included. The low-resolution simulation was run with 128 3 particles inside a box of 20 Mpc h −1 per side, with the box initially covered by a mesh of 128 3 cells (zeroth level cells). At z = 0, we searched for MW-like mass halos (7.0 × 10 11 M   Ḿ 1.5 10 vir 12 M e ) that were not contained within larger halos (distinct halos). We only selected halos that had not had major mergers since z = 1.5 and that, at z = 0, do not have a similar mass companion inside a sphere of 1 Mpc h −1 . Other restrictions we have imposed are that halos need to be in a filament or a wall, not in a void or a knot. After this selection, a Lagrangian region of 3 R vir was identified at z = 60 and re-sampled with additional small-scale modes (Klypin et al. 2002). The virial radius in the lowresolution runs, R vir , is defined as the radius that encloses a mean density equal to Δ vir times the mean density of the universe, where Δ vir is a quantity that depends on Ω 0 , Ω Λ , and z. For example, for our cosmology Δ vir (z = 0) = 338 and Δ vir (z = 1) = 203. The number of DM particles in the high resolution region depends on the number of DM mass species and the mass of the halo. For models with four or five DM mass species this value varies from ∼1.5 × 10 6 (model G.242) to about 7.0 × 10 6 (model G.321). The corresponding DM first species mass per particle of modeled galaxies (m DM1sp ) is given in Table 1.
In ART, the initially uniform grid is refined recursively as the matter distribution evolves. The cell is refined when the mass in DM particles exceeds 81.92 (1-F b U , ) m p /f DM or when the mass in the gas is higher than 81.92 is the universal baryon fraction, m p = 7.75 × 10 4 M e h −1 is the mass of the lightest particle species in the DM only simulation, and f DM and f g are factors that control the refinement aggressivity in the DM and gas components, respectively. In our models F b U , = 0.15 according to the cosmology we have imposed, and the mean DM and gas densities in the box are assumed to be equal to the corresponding universal averages. For the simulations presented in this work the grid is always unconditionally refined to the fourth level, corresponding to an effective resolution of 2048 3 DM particles. A limit also exists for the maximum allowed refinement level and this is what will give us the spatial size of the finest grid cells. Here we allow the refinement to go up to the 11th refinement level, which corresponds to a spatial size of the finest grid cells of 109 comoving pc. Our main model is a re-simulation of a MW-like halo of M vir ∼ 7.4 × 10 11 M e . However, we also re-simulated other halos with masses between 13.9 × 10 11 M e and 2.0 × 10 11 M e . We have used these simulations to determine how the general properties of the galactic system change with mass. For each one of these models we have also performed simulations that involve changing initial parameters such as ò SF , n SF , resolution, or the number of DM mass species. This set of models helps us to assess how the final results depend on the initial conditions. A more detailed study of the effects of changing the simulation parameters can be found in González-Samaniego et al. (2014b). Table 1 summarizes the properties of the main MW-sized simulations presented here.
Although our MW-sized simulations are realistic, some physical processes are not well implemented in the code we have used, for example, the AGN feedback, the molecular cooling down to 200 K, the formation of molecular H 2 clouds, deep chemical treatment, the radiation pressure, and the photoionization from massive stars. However, the last two are now implemented in a new version of the code and a new set of simulations is being generated.
In this work we study MW-sized systems, however, it is not our goal to reproduce all observed properties. Here we define a MW-sized system as the one that has a total mass, v max , and disk with similar mass and scale lengths as the ones observed. Following previous works (e.g., Governato et al. 2004;Robertson et al. 2004;Okamoto et al. 2005;Scannapieco et al. 2009;Guedes et al. 2011;Aumer et al. 2013;Stinson et al. 2013;Moster et al. 2014;Vogelsberger et al. 2014) we have selected a DM halo that in the high resolution run at z = 0 has no other massive halos inside a 800 kpc sphere, and that has a quiescent recent assembling history. We note that our initial conditions do not reflect the environment in which the MW is embedded, but this is certainly not needed if we just aim to produce a system with a non-negligible disk component. This latter is almost always achieved under the requirement that the system has evolved in relative isolation for the last 10 Gyr or so (e.g., Aumer et al. 2013;Stinson et al. 2013;Moster et al. 2014;Vogelsberger et al. 2014). These specific systems will be our main MW-sized models. Details of each realization will depend on the selected subgrid physical parameters, environment, and assembling histories.
Note. r vir is assumed to be equal to r 97 . Mass related parameters: M vir is total mass inside r vir ; M * and M gas are the total stellar and gas masses inside r vir ; the hot, warm, and cold-gas masses (M hotgas , M warmgas , M coldgas ) are the gas masses inside r vir with temperatures of T > 3 × 10 5 K, 3 × 10 5 K > T > 3 × 10 4 K, and T < 3 × 10 4 K, respectively; F b U , is the baryonic fraction (M * + gas /M vir ); m DM1sp is the mass of a single particle belonging to the first DM mass species; * m min max are the lowest and highest masses of stellar particles; and m SPH is the mass of SPH particles. Structural parameters: c is the DM halo concentration parameter defined as r vir /r s ; R d is the disk scale length; h z,young old are the disk scale heights of young (0−0.5 Gyr) and old (4.0−11.0 Gyr) stellar populations; and α X is the hot-gas density power-law exponent (ρ hotgas (r) ∝ a r X ). Kinematical parameters:  V c is the circular velocity at solar position (R e = 8 kpc); R peak is the radius within which the circular velocity curve reaches its highest value (V c (R peak )). Initial condition parameters: n SF , T SF, and ò SF are the SF density and temperature thresholds, and the SF efficiency; Box is the simulated box size in Mpc h −1 ; z ini is the initial redshift; and Ω 0 , Ω Λ , Ω b , H 0 , and σ 8 define the assumed cosmological model. Comparisons with other stateof-the-art simulations are also shown (ERIS simulation of Guedes et al. 2011;and model B of Mollitor et al. 2015). Observed values for the MW have been obtained from Ferrière (2001) Table 1. In the table we see that there are parameters which are very sensitive to changes in the refinement settings as, for example, the SFR (z = 0) which changes from 0.02 to 0.32, for G.322 and G.323, respectively. Yet, most parameters do converge as can be seen in Table 1; in particular, the circular velocity profiles agree within 5% (see Appendix A for more information).
To ensure numerical convergence we also generated a lowresolution model. A comparison between models can be found in Appendix A.

Morphology
At z ∼ 0 our three G.3 models are massive spiral galaxies with several non-axisymmetric structures in the disk, such as bars or spirals (see Figures 1 and 2). Their assembling history has been quiet after z = 1.5 and the last major merger occurred at z = 3. In Figure 1 we show face-on (first column) and edge-on (last two columns) views of the stars of our model G.323 at z = 0. For the other two models the main picture is similar but with a smaller number of stellar particles. Each row correspond to a different stellar age (see the figure caption). In this figure it can be observed that the young stellar component (0−4 Gyr) is distributed in a flat disk structure. Older stellar populations (4−10 Gyr) are also distributed in a disk structure, but it is not as flat. We argue that the disk scale height and length depend on the age of the stellar population and that the former is higher and the latter is smaller when the population becomes older; this issue will be addressed in Section 3.4.2. Another interesting property that can easily be observed in the edge-on views of our models is that the stellar population has an evident flare which is more evident for older populations. The face-on view of the younger stellar population shows the presence of rings, spiral arms, and also a young bar in the disk. The spirals, rings, and young bar can better be observed in models G.322 and G.323 where the number of stellar particles is higher (see Figure 2). The young bar has grown from disk particles that have ages between 0 and 8−9 Gyr. The old stellar population (11 −13.476 Gyr) is distributed in a bar/ellipsoid component and it originated in the last major merger at z = 3 (Valenzuela, O., et al. 2016, in preparation). In Figure 3 we show the gas component as function of temperature. In this figure it is easy to see that the gas distribution is not isotropic. The gas is distributed in different regions of the system depending on its temperature: cold gas is present in the young stellar disk region and hot gas fills the out-of-plane region and is embedded in the DM halo. It is also interesting to see that the cold/intermediate gas has a warped structure that is also marginally observed in the stellar component.

Stellar, Gas, and DM Virial Masses
Here we define the virial radius (r vir ) as the one where the sphere of radius r vir encloses a mean density 97 times denser than the critical density (ρ crit ) of a spatially flat universe ρ crit = 3H 2 (z)/(8πG). We have used the value 97 as it is the value derived from the spherical top-hat collapse model for ΛCDM at z = 0 for our cosmology (Bryan & Norman 1998). This virial radius definition is borrowed from structure growth theory and thus its use is not appropriate when one wants to define a physically meaningful halo edge (Cuesta et al. 2008;Zemp 2014). To avoid this problem, we also give the properties for another commonly used definition which is r 200 , the radius that encloses a mean density equal to 200 times ρ crit . Using the first definition we have obtained that the virial radius at z = 0 is r vir ∼ 230 kpc in our three main models and that the mass enclosed in this radius is M vir = 7.20−7.61 × 10 11 M e . This value for M vir falls well inside the observational range, M vir = 0.6−2.4 × 10 12 M e , found by Xue et al. The total mass enclosed at r 200 = 175.6 kpc is M 200 = 6.71 −6.90 × 10 11 M e . The total virial mass is distributed in dark, stellar, and gaseous matter as follows: M DM = 5.86 −6.79 × 10 11 M e , M * = 6.1−6.2 × 10 10 M e , and M gas = 1.73 −2.70 × 10 10 M e . The baryonic fraction of our G.3 models is F b U , = 0.104−0.121. These values are 19%-31% smaller than the universal value for the adopted cosmology which is F b U , = 0.15. Inside r vir , and for the model G.321, we have a total number of particles of N total = 7.52 × 10 6 , where N DM1sp = 7.12 × 10 6 and N * = 3.94 × 10 5 , and a total number of gas cells of 2.0 × 10 6 . For models with a more aggressive refinement the numbers of stellar particles and gas cells grow up to 2.34 × 10 6 and 6.8 × 10 6 , respectively. All DM particles inside the virial radius belong to less massive DM species (1sp) and have a mass of 9.25 × 10 4 M e . Star particles have masses between ∼10 3 and 1.2 × 10 6 M e . All these parameters are summarized in Table 1.

Circular Velocity Curves
Circular velocity curves for our main simulations, computed using the < GM r r approach, are shown in Figure 4 (top left panel) as a solid black line. The circular velocity curves have also been computed using the real galactic potential using the Tipsy package. The results obtained using both techniques are in good agreement.  Pichardo et al. (2003). The solid lines are the total circular velocity curves while the DM, stellar, and gas components are shown as long dashed, dotted, and short dashed lines, respectively. In the top right panel we show, the G.321 velocity rotation curves of stars and gas, computed as the mean of the tangential velocity in rings centered to the galactic center (cyan solid and dashed lines), and the total circular velocity curve obtained using the V c = < GM r r approach. As expected, the gas follows the circular velocity while stars are affected by the asymmetric drift. In this panel we also compare the stellar rotation curve of our simulation (cyan solid line) with values from MW observations. We show rotation velocity curve estimations from blue horizontal-branch halo stars in the Sloan Digital Sky Survey (SDSS; Xue et al. 2008) (cyan and magenta dots), from Sofue et al. (2009) (red dots), and from the observations of López-Corredoira (2014) (green dots). In the bottom panel we compare the recent data for the circular velocity curve obtained by Reid et al. (2014), using high mass star-forming regions, with the total circular velocity curve of our G.321 model. As can be seen in the figure our G.3 models are in very good agreement with the most recent observational data. We also show the recent hypothetical gaseous rotation curve of the MW obtained by Chemin et al. (2015) after correcting for non-circular bar motions. It is important to mention that we only show results for the G.321 model because the circular velocity curves of models G.322 and G.323 do not differ significantly. The peak of the z = 0 circular velocity of our models is reached at R peak ∼ 5.69 kpc with a value of V c (R peak ) = 237.5−243.8 km s −1 , the value at a standard solar radius (R e = 8 kpc) is  V c = 233.3−239.8 km s −1 . The ratio between the circular velocity at 2.2 times disk scale radius (V 2.2 ) and the circular velocity at r 200 (V 200 ) of our simulations (V 2.2 /V2 00 1.9) is also inside the observational range for the MW, which is -+  (Xue et al. 2008;Dutton et al. 2010).
An interesting exercise has been to compare the circular velocity curves of our simulations with that of the Guedes et al. (2011) model. We can easily see that our simulations have slightly higher velocity curves and do not present a peak in the inner regions compared to Guedes et al. (2011). In simulations this internal peak is usually associated with the presence of a massive bulge in the central parts of the simulated galactic . Face-on (left) and edge-on (right) views of gas density at z = 0 in model G.323 as function of its temperature. From top to bottom: 10 6 -10 7 K, 3 × 10 5 -10 6 K, and 10 3 -3 × 10 5 K. The color scale shows the gas temperature.
system. On the other hand, a similar peak in the V c is shown in the observational work of Sofue et al. (2009). It is under discussion whether this peak detected in the Sofue et al. (2009) observations is a signature of non-circular motions inside the bar or of a massive bulge (Duval & Athanassoula 1983;Chemin et al. 2015). Non-circular motions need to be computed to solve such a question. We plan to undertake this work in the near future following that one started by Chemin et al. (2015). Finally, despite the fact that the rotation and circular velocity curves we introduce here fall well inside The total V c of our model G.321 at z = 0 is shown as a black solid line. In green we show the hypothetical gaseous rotation curve of the MW obtained by Chemin et al. (2015) after correcting for the non-circular motions in the velocity profile of the MW inferred with the tangent-point method of Kalberla et al. (2005). In gray we show the Chemin et al. (2015) curve's 1σ range. observational ranges, it is important to take into account that current observational uncertainties are still high.

DM Component
As we mentioned in the previous section, all DM particles inside the virial radius are particles of the first DM species, the less massive one. In Figure 5 we show the DM density profile as a function of radius (red circles) and the best-fit of the Navarro-Frenk-White (NFW) profile for the model G.321 (upper yellow solid line). We avoid the central region (R < 5 kpc) in the fit of the NFW profile to avoid perturbations derived from the presence of baryons (adiabatic contraction). The best NFW fit gives larges values for the halo concentration, which are c = r vir /r s = 28.5, 26.8, and 26.9 for G.321, G.322, and G.323, respectively. These values fall inside the observational range presented in Kafle et al. (2014), 14.8 . No core is observed in the central region of our G.3 models. The halo spin parameters as defined in Bullock et al. (2001) are λ' = 0.019, 0.017, and 0.022. These λ' values also fall near the range obtained by Bullock et al. (2001) after analyzing more than 500 simulated halos (see Figures 1 and 2 in Bullock et al. 2001).

Gas Component
In our models G.321, G.322, and G.323 the total gas mass inside the virial radius is M gas = 2.7, 1.73, and 1.96 × 10 10 M e . In G.321, 9.34 × 10 9 M e of the gas mass is in the cold-gas phase (T < 3 × 10 4 K). In G.322 and G.323 the cold-gas masses are 1.52 and 1.86 × 10 9 M e . The cold-gas mass in G.321 is comparable to the total mass of the molecular, atomic, and warm ionized medium inferred for the MW, ∼7.3 −9.5 × 10 9 M e (Ferrière 2001). For the rest of the models presented here the amount of cold gas is smaller than the observed value. This z = 0 cold-gas mass reduction is a consequence of an increase in SF at early ages, which in its turn is provoked by the higher aggressivity in the refinement criteria. This increase in early SF is the only significant effect we have observed when changing the refinement criteria and it is a consequence of resolving a larger number of small dense regions in which the SF criterion is accomplished. At z = 0, the SF criterion is satisfied only in the inner disk regions (R < 7 kpc). This spatially limited SF accords with the young stellar distribution shown in Figure 1 and also with the decrease of the SFR at z = 0, which can be seen in Figure 10 below.
The hot-gas masses (T > 3 × 10 5 K) are around 1.22, 0.98, and 1.32 × 10 10 M e for models G.321, G.322, and G.323, respectively. The definition of hot gas we use in this work is observationally motivated: total hot-gas mass has to be inferred from ionized oxygen observations (Gupta et al. 2012) and this is only possible when the gas has opacity in the X-ray region (T > 3 × 10 5 K). In Figure 3 we show the gas distribution as a function of temperature in model G.323. Each panel shows the distribution of gas at different temperatures (see the figure caption). As can be seen in this figure, cold gas is located in the disk region while warm-hot gas is embedded within the DM halo. Contrary to the standard assumption, hot gas does not follow the DM radial distribution. This can be seen in Figure 5, where we show the DM density distribution (red) and the hotgas density distribution (blue). The hot-gas density distribution inside r = 100 kpc can be fitted by a power law (ρ hotgas (r) ∝ a r X ). After fitting data from models G.321, G.322, and G.323, we have obtained that the hot-gas density distribution scale factors (α X ) of these models are −0.62, −0.83, and −0.67, respectively. These results are slightly higher than the values used to fit analytical models of hot gas in DM halos (−0.9 in Anderson & Bregman 2010) and the ones obtained in other similar models like ERIS (−1.1 in Guedes et al. 2011). The hot-gas density distribution power-law fit for model G.321 is shown in Figure 5 (bottom yellow straight line).
In Section 4 we present our first results on the study of the hot-gas distribution and a comparison with observational values.

Stellar Component
The total stellar mass in the virial radius of our G.3 models is M * = 6.1−6.2 × 10 10 M e , which is comparable to the values estimated for the MW (4.5− 7.2×10 10 M e in Flynn et al. 2006 andLicquia &Newman 2015). However, if we observe the relations between stellar mass and total mass (see Figure 6) derived from our models (blue dots) and also from observations of the MW (shadowed region) we see that most of them fall above the * M /M vir predicted by cosmological theories. The same result is observed when using data from other recent MW-sized simulations like the ERIS simulation (Guedes et al. 2011) or those of Mollitor et al. (2015). In this work we have not studied this mismatch in detail, we leave this to future works. Figure 5. Average DM (red dots) and hot-gas (T > 3 × 10 5 K, blue dots) density profiles at z = 0. The upper solid yellow line shows the best NFW profile fit to the DM mass distribution, in the range from 5 kpc, outside the barbulge region, to r vir . The lower solid yellow line shows the best-fit power-law profile to the hot-gas mass distribution between 5 and 20 kpc. The DM best-fit NFW profile is characterized by a large halo concentration parameter c = 28.5 as the DM halo contracts in response to the condensation of baryons in its center (adiabatic contraction). In the inner region it can clearly be observed that no DM core is present. The hot-gas best-fit power-law profile gives a slope of α X = −0.62.

Stellar Disk and Stellar Spheroid
Simulated galactic systems in our main models (G.3 series) have both stellar disks and stellar spheroidal components (halo and/or bulge). To study the properties of stellar galactic disks we need to know which stars belong to them. To undertake the selection process we have used an extension of the kinematic decomposition proposed by Scannapieco et al. (2009). In Scannapieco et al. (2009) the authors used the ratio between the real stellar particle angular momentum perpendicular to the disk plane ( j z ) and that obtained assuming stellar particles follow a circular orbit (computed from the circular velocity curve, j c ). The circular velocity curve has been computed using both the GM/rapproach and the real galactic potential; we have found that the results do not depend on the v c computation technique. Using the j z /j c ratio it is possible to distinguish particles that are rotationally supported ( j z = j c for disk stars) from ones that are not ( j z < j c for spheroid stars). Here we have decided to make additional cuts to improve this technique. We have performed complementary cuts in metallicity, in the vertical coordinate | | z , in the angle between the rotation axis of each particle and the vector perpendicular to the plane ( ( ) a cos ), and in the distance to the galactic center (R). To avoid kinematic biases we have made cuts that do not involve kinematics first. We started the selection process by making the metallicity cut; next, we made a cut in the vertical coordinate, and then in then radius, in ( ) a cos ; finally we have used the j z /j c condition. Some of these restrictions we have imposed require previous knowledge of the disk-plane position. To find this position we have used an iterative geometrical approach (Atanasijevic 1971, pp 9-11). In this approach we end up with the plane that minimizes the cumulative distance of all particles to it. In the first iteration of this process we only use stars that are inside a 20 kpc sphere centered within the center of mass of the lightest DM mass species particles. After this first iteration we only use stars that fall near the newly defined plane (i.e., | | z new < 5 kpc). We have checked that after a few iterations we obtain a plane that coincides with the disk plane defined by young stars and cold gas (see the bottom panels of Figures 1-3).
After analyzing the stellar distribution of our G.3 models in both Cartesian and phase spaces, we have found conditions that better distinguish disk from spheroid particles. In our main models these conditions are: j z /j c > 0.55, | | z < 3.5 kpc, log (ZIa/ZIa e ) > −0.5, where ZIa e = 0.00178 and is the iron solar abundance according to Asplund et al. (2009), cos(α) > 0.7, and R < 25 kpc. All particles inside r vir that do not accomplish one or more of these restrictions are considered spheroid particles (halo and/or bulge particles). To account for thick disk particles we relax the log(ZIa/ZIa e ) and j z /j c conditions; that is, stellar particles will also belong to the disk component if they fulfill all previous conditions except that they have log (ZIa/ZIa e ) < −0.5 and j z /j c deviates by less than 0.75 from the j z /j c peak. 8 The result of this selection can be seen in Figure 7 where we show the stellar mass fraction as a function of j z /j c  . Stellar mass fraction as a function of the orbital circularity parameter j z /j c , describing the degree of rotational support of a given stellar particle (black: all particles; red: spheroid; blue: disk). The stars in a centrifugally supported thin disk manifest themselves in a sharply peaked distribution about unity. We have only plotted stellar particles inside the disk region (| | z < 0.5 kpc and r < 15 kpc) of model G.321 at redshift 0. for model G.321, the result is similar for the other models. In this figure we show how well we trace the two main stellar components of the model: one centered at j z /j c = 1, i.e., rotationally supported, and another at j z /j c = 0.
The masses of the rotationally supported component (stellar disk) are around M d = 1.82, 2.17, 2.21 × 10 10 M e and those of the spheroidal one (bulge + halo) are M sp = 4.29, 3.93, 3.79 × 10 10 M e for models G.321, G.322, and G.323, respectively. Comparing these results with observations we find that the mass of the spheroidal component in our models (bulge + halo stars) is higher than the upper limit obtained by Flynn et al. (2006), M b 1.3 × 10 10 M e , and the most recent observations presented by Valenti et al. (2016), M b = 2.0 ± 0.3 × 10 10 M e . The spheroids of our galactic systems are ∼2 times heavier than that of the MW. On the other hand, the disks in our G.3 models are less massive than that estimated for the MW by Bovy & Rix (2013) (∼4.6 × 10 10 M e ). As a consequence, the kinematic disk-tototal ratio (D/T) decomposition, f disk , computed as * M disk /M * (Scannapieco et al. 2010(Scannapieco et al. , 2015, where * M disk is the mass of stellar particles with j z /j > c 0.5, results in a value (∼0.42) that is significantly lower than the one obtained in MW observations (∼0.75 in Scannapieco et al. 2011). However, as pointed out by Scannapieco et al. (2010Scannapieco et al. ( , 2015 the use of a kinematical decomposition results in a lower D/T than those obtained using a photometric decomposition, such as the one determined in MW observations. For instance, in Scannapieco et al. (2010) the authors obtained a set of MW-sized galaxies with a kinematical f disk ∼ 0.2 and they demonstrated that f disk becomes 0.4−0.7 when using a photometric decomposition.
Is the presence of this massive spheroidal component due to a large stellar concentration in the central region of our simulated galaxies, a common result in earlier hydrodynamical zoom-in simulations? As can be seen in the rotation curves shown in Figure 4 this is not the case (the curves do not present a central peak). We have analyzed this spheroidal component and we have detected that it is mainly a triaxial structure that was formed in the major merger occurring at z = 3 (see Section 3.1.1). This triaxial structure has a low density in the disk region and thus it is not a classical massive bulge.

Disk Properties
With the stellar disk/spheroid selection process we have found that the disk component has 2.0, 6.0, and 18.7 × 10 5 particles in each one of our main G.3 models. We have also computed the mean volume and surface density of the disk as function of radius. We have found that a single or a double exponential power law can be fitted to the data depending on the age of the selected population. In Figure 8 we show the total disk surface density (black) as function of the radius of the model G.321 at z = 0. We have fitted two exponential power laws (red and blue) to the surface density curve. The results from these fits show that the scale length in the inner regions (2 −6 kpc) is R d = 4.89 kpc and in the outer ones (6−12 kpc) R d = 2.21 kpc. If we fit a single exponential to the whole radial range from 2 to 12 kpc we obtain that R d = 2.56 kpc, a result that is in agreement with values obtained for the MW (e.g., 2.3 ± 0.6 kpc in Hammer et al. (2007) or 2.15 ± 0.14 kpc in Bovy & Rix (2013) using G-type dwarfs from SEGUE). It is important to note that, as can be seen in Figure 8, a small concentration of stars is present in the central regions (R < 2 kpc). This concentration is caused by the presence of a young bar.
In Figure 9 we show the volume density as function of the distance to the plane (z) of three disk stellar populations split by age. We have selected young (0−0.5 Gyr), intermediate (0.5 −4.0 Gyr), and old (4.0−11.0 Gyr) populations. The densities have been computed using vertical bins of Δz = 0.2 kpc and only stellar particles that are between R = 2 and R = 9 kpc. We have found that the scale height (h z ) of each population, when fitting a decreasing exponential to the vertical density distribution in the range 0.1 kpc <| | z < 2.0 kpc, is h z = 277 pc for the young, h z = 959 pc for the intermediate, and h z = 1356 pc for the old population. These values have been obtained from the G.321 model, results from models G.322 and G.323 can be found in Table 1. These scale height  values are compatible with those expected for the MW: observations show that the vertical distribution of stars in the MW can be fitted by two exponentials with scale heights of h z = 300 ± 60 pc (Jurić et al. 2008) and h z = 600−1100 ± 60 pc (Du et al. 2006). Our values also follow the relation proposed by Yoachim & Dalcanton (2006) who argue that the disk scale height in edge0-on galaxies increases in disks following z 0 ≈ 2h z = 610(V c /100) km s −1 .

Star Formation History (SFH)
In our simulations stars are being formed with SFR = z 0 = 0.27 M e yr −1 , 0.18 M e yr −1 , and 0.41 M e yr −1 , for models G.321, G.322, and G.323, respectively. These values are smaller than those inferred by Robitaille & Whitney (2010) using Spitzer data for the MW (SFR = 0.68−1.45 M e yr −1 ) or that from Licquia & Newman (2015), 1.65 ± 0.19 M e yr −1 . In Figure 10 (top panel) we show the SFR as function of the redshift for the total (black), spheroidal (red), and disk (blue) components. We have computed the SFH of each stellar component using the formation time of their corresponding stellar particles (saved by the code) inside r vir at z = 0. We have followed the strategy described in Section 3.4.1 to distinguish between disk and spheroid stellar particles. We also show, as the shadowed region, the predicted SFH for MW-like halos derived from semi-analytical models that combine stellar mass functions with merger histories of halos (Behroozi et al. 2013). The peak in the SFH of our G.3 models occurs at slightly earlier ages than that predicted for a MW-sized galaxy. In the bottom panel we plot the total stellar mass as function of redshift. Some important results that can be appreciated from this figure are, first, that the spheroidal component is being build at high redshifts, while the disk starts its formation at around z = 2.5 (∼11 Gyr from the present time), just after the last major merger that was at z = 3. Second, it is also important to note that the SF of the spheroidal component decreases quickly after z ∼ 2.4 and become negligible at around z = 0.5. The disk SFR also decreases fast in the last time instants except for some short periods of SF that are a consequence of the accretion of small gaseous satellites. This reduction in SF when reaching z = 0 is a consequence of the reduction of cold gas available for SF. Several processes can lead to a reduction of the cold-gas mass fraction. First, cold gas can be drastically consumed at higher redshifts due to a non-realistic implementation of the physical parameters that control SF (Liang et al. 2016). When such a problem exists, a large number of old stars is present at z = 0 in the inner disk region. In our G.3 models, as can be seen in Figure 7, an old spheroid is present in the disk region, however, this old stellar population is not large enough to account for the total reduction of the observed cold-gas component. Second, if stellar and SN feedback are too efficient, an significant amount of cold gas can be heated up and thus become unavailable for SF. In this second scenario what we would expect is that the fraction of hot gas increases with time, a behavior that is not observed in our G.3 models. Finally a decrease in the cold-gas inflow due to inhomogeneities in the CGM can also be a cause of the reduction of the SFR. In this last hypothesis we would observe only a small amount of cold gas present around our system, or falling from filaments, in the last Gyr. We argue that in our models the reduction of the amount of cold gas is a consequence of a combined effect of the first and the last hypotheses, as we observe that an old spheroid is present and also that there is not a large amount of cold gas falling from filaments at z = 0, however, we leave the more detailed study of such processes for the future.
In ART, the metallicity Z is divided according to which kind of SN produced the metals: ZII or ZIa, if they are produced by SN II or SN Ia, respectively. The metallicity Z is then = + Z Z Z Ia II. Alpha elements, such as oxygen, are mostly produced by SN II while iron is mostly produced by SN Ia; we can thus assume that the abundance of alpha elements is proportional to ZII while the iron abundance is proportional to ZIa.

Are Our G.3 Models MW-like Systems?
Although it is not the aim of this paper to obtain a realistic model of the MW galaxy, nevertheless the models we present here (i.e., G.321, G.322, and G.323) can be considered at least MW-sized as they reproduce several of the MW's observed properties. On the other hand, it is also evident that some parameters of these models are different from observations of the MW. Below, we give an overview of the parameters that match MW observations and of those that do not.

MW-like Properties
Our study is focused on a galactic system, the G.3 series, that forms spiral arms and a galactic bar. Its assembly history is quiet after z = 1.5 and there are no major mergers after z = 3, something that is also expected for our Galaxy (Forero-Romero et al. 2011). Its total mass is inside the observational range proposed for the MW by Xue et al. (2008), Boylan-Kolchin  Kafle et al. (2012Kafle et al. ( , 2014, 1.0 ± 0.4 × 10 12 M e . The stellar mass inside the virial radius (6.1-6.2 × 10 10 M e ) also falls inside the observational range for the MW, 4.5-7.2 × 10 10 M e (Flynn et al. 2006;Licquia & Newman 2015). Finally, the cold-gas mass inside the virial radius in model G.321 (9.34 × 10 9 M e ) also matches the observed values (7.3-9.5 × 10 9 M e in Ferrière 2001). Several structure parameters of our G.3 models are also within the MW observations. This is the case for the DM halo concentration parameter (c) which in our models takes values from 26.8-28.5, while observations predict it should be -+ 21.1 8.3 14.8 (Kafle et al. 2014). The hot-gas power-law profile index (α X ), defined in Section 3.3, of our G.3 models is between −0.62 and −0.83, which is also close to the value proposed by Anderson & Bregman (2010), −0.9. The disk scale lengths and heights of our galactic systems (G.3 series) also coincide or are close to observed values. The former in our models is 2.56−3.2 kpc and observed values are 2.3 ± 0.6 kpc in Hammer et al. (2007) and 2.15 ± 0.14 in Bovy & Rix (2013). The observed disk scale height of young and old stars is 300 ± 60 pc in Jurić et al. (2008) and 600−1100 pc in Du et al. (2006), respectively, and in our models we have obtained for young and old populations, 277−393 pc and 960−1356 pc, depending on the model. Finally, several kinematical parameters also agree with the MW observations; the rotation curve of our simulated galaxies, for instance, roughly match observations (see Figure 4). Some examples of such agreement are the ratio V 2.2 /V 200 , which is 1.9 in our models and 1.67 - Xue et al. (2008), and v c (R e ) which is 233.3−239.8 km s −1 in this work and 221 ± 18 in Koposov et al. (2010) and 236 ± 11 in Bovy et al. (2009).

Parameters that Do Not Match Observations
Among the parameters that do not agree with the observed values of the MW is the SFR and the bulge to disk ratio. The cold-gas fractions in the G.323 and G.323 models also do not match observations, they account for only 1.52 and 1.86 × 10 9 M e , while the observed values are 7.3-9.5 × 10 9 M e Ferrière (2001). The SFRs in our models range from 0.18 to 0.41 M e yr −1 , below the observational range 0.68-1.45 M e yr −1 from Robitaille & Whitney (2010). Finally, we have presented in Section 3.4.1 the stellar spheroiddisk decomposition and shown that a massive spheroid exists in our models. Such a massive structure has not been observed in the MW. The total masses of the spheroid (bulge+halo) and the stellar disk in our G.3 models are 3.79−4.29 ×10 10 M e and 1.82−2.21 × 310 10 M e , respectively. Moreover, observations show that the bulge and the disk masses of our Galaxy are 2.3 × 10 10 M e (Flynn et al. 2006;Valenti et al. 2016) and ∼4.6 × 10 10 M e (Bovy & Rix 2013), respectively; that is, the spheroid of our models is at least two times more massive than the bulge of the MW.

The Missing Baryons Problem
It has long been known that the cosmological baryon fraction inferred from Big Bang nucleosynthesis is much higher than the one obtained by counting baryons at redshift z = 0 (e.g., Fukugita et al. 1998). However, it was recently, with WMAP and Planck high precision data of the cosmic microwave background (Dunkley et al. 2009;Planck Collaboration et al. 2014), that the cosmic baryonic fraction was well constrained and thus the lack of baryonic mass in galaxies became evident. This data revealed that the cosmic ratio between baryonic and total matter (Ω b /Ω M ) is 3−10 times larger than the one observed in galaxies. For example, Hoekstra et al. (2005) found a baryonic fraction for isolated spiral galaxies of 0.056 and for elliptical of 0.023, while Dunkley et al. (2009) presented a cosmic baryonic fraction from WMAP five year data of about 0.171 ± 0.009.
Studies like Cappi et al. (2013) and Georgakakis et al. (2013) tried to solve this missing baryon problem by proposing that galactic winds, SN feedback, or strong AGN winds ejected baryons to the CGM. Others proposed that most of the gas never collapsed into the DM halos as it was previously heated by Population III SNe (e.g., Mo & Mao 2004), this is known as the preheating scenario. In this latter scenario, missing baryons are still in the extragalactic warm-hot intergalactic medium (WHIM), following filaments. Thus, it has become clear that studying the CGM is necessary not only to find the missing baryons, but also to understand its properties in the context of galaxy formation and evolution models. To further constrain galaxy evolution and find out which one of the proposed theories solves the missing baryons problem, it is necessary to measure the amount of warm-hot-gas phase present in the CGM and also its metallicity. The results of such measurements give information about the heating mechanism: hot gas with low metallicities is indicative of Pop III SN preheating while a more metallic gas suggests that it has been enriched in the disk via SNe and stellar winds.
Supporting the first hypothesis, the idea of very strong feedback driving galactic-scale winds causing metal-enriched hot gas to be expelled out of the star-forming disk possibly to distances comparable to or beyond the virial radius was generally accepted. This process also lies at the heart of obtaining realistic simulations of disk galaxies (Guedes et al. 2011;Marinacci et al. 2014). By this mechanism gas is temporally located far from the star-forming regions delaying SF and preventing the formation of an old bulge in the center of the system.
From the analysis of O VII and O VIII X-ray absorption lines observed toward the line-of-sight of extragalactic sources, several authors performed estimations of the hot-gas mass in the MW halo, suggesting that an important portion of the missing baryons is located in this hot-gas phase (Gupta et al. 2012(Gupta et al. , 2014. They concluded that the warm-hot phase of the CGM is extended over a large region around the MW, with a radius and mass around 100 kpc and 10 10 M e . However, these estimations are limited to a few directions toward extragalactic sources (bright QSOs) and as a consequence issues like the homogeneity or isotropy of the hot-gas distribution possibly affecting the mass estimation are difficult to address. Other studies, such as Feldmann et al. (2013), also suggest that the presence of this hot-gas phase can explain part of the isotropic gamma-ray background observed by the Fermi Gamma-ray Space Telescope.

Spatial Distribution
The amount of hot gas in our G.3 models is M hot = 0.98 −1.32 × 10 10 M e . This hot gas, which has a temperature above 3.0 × 10 5 K (having opacity in the X-rays), is embedded in the DM halo but mostly outside the stellar disk. Aiming to perform a fair comparison with observations, we estimated the hot-gas hydrogen equivalent column density and emission measure (EM) 9 as: where n H is the hydrogen particle density, dl is an element of the path in the line-of-sight, and n e is the electron density.
In Figure 11 we show a full-sky view of the hot-gas column density distribution in our model G.321. This figure has been obtained by computing the hot-gas column density from a position that is at 8 kpc from the galactic center, inside the simulated galactic disk, and resembles the Sun position, and assuming arbitrary azimuthal angles. We have obtained column density values that fall near the observational ranges, if a solar metallicity is assumed (see Table 3 and Gupta et al. 2012). However, when assuming lower metallicities, observations are consistent with values that are above the ones measured in our simulations. This change can give us information about the absorption of the local weakly interacting massive particle. We also conclude that the distribution of hot gas inferred from all estimations is far from homogeneous. In order to quantify the hot-gas anisotropy we followed two strategies. First we computed the amplitude of the first spherical harmonics (Y l m from l,m = 0 to 5 ) and later we also computed the filling factor as defined by Berkhuijsen (1998). Our results show that the dominant spherical harmonic is Y 1 0 , which is an indicator of the dipolar component, with a small contribution of several high order components. The asymmetry of the hot-gas distribution with respect to the disk plane suggest some degree of interaction with the extragalactic medium. The computation of the mean filling factor from lines-of-sight distributed all along the sky gives a value of f v ∼ 0.33 ± 0.15. This value for the filling factor also suggests that hot gas is concentrated in a few regions of the sky. It is also important to mention that in our G.3 models (see Figures 3 and 11) we do not detect a hotgas thick disk component. This is important as the contribution of the hot-gas thick disk to the observed X-ray emission/ absorption is under discussion (Savage et al. 2003). In addition to the spatial distribution we study the metallicity and velocity components of hot-gas cells. The metallicity and velocity fields provide us information about how the hot-gas component interacts with the extragalactic medium. In Figure 12 we show the metallicity distribution of the hot gas as seen from the Galactic center. It is evident that the gas metallicity distribution does not depend on the azimuthal angle (θ) (Figure 12, top panel) while the dependence on the vertical latitude (f) is clear (Figure 12, bottom panel). A low metallicity component located in the northern hemisphere is also noticeable, which is absent from the southern hemisphere, where the metallicity is on average much higher. It is important to note that in Figure 12 (bottom panel) we can distinguish a clear break in the hot-gas distribution at low latitudes. This hole coincides with the position of the stellar/cold-gas disk. Figure 11. Hot-gas (T > 3 × 10 5 K) column density (top) in a full-sky view of simulation G.321 in galactic coordinates at z = 0. All values have been computed as observed from R = 8 kpc and at an arbitrary azimuthal angle. Figure 12. Hot-gas (T > 3 × 10 5 K) view in spherical galactic coordinates of the G.321 model at z = 0. Top: metallicity as a function of the azimuthal angle. Bottom: metallicity as a function of the vertical angle. ZIa and ZII are metals from SNe Ia and II, respectively. All values have been computed as observed from the galactic center. 9 In this work we are not showing the dispersion measure values as they only differ by a factor of two from the N H ones. N H values are presented in Figure 11. Although not presented here, the dispersion measure values in our models are of the same order as those in Guedes et al. (2011).
In order to shed some light on the bimodal origin in the metal distribution, we have analyzed a set of three color maps showing the projected metallicity and velocities along each one of the principal planes (see Figures 13 and 14). Figure 13 (the X-Yplane, i.e., the disk plane) shows that in general hot gas in the disk plane rotates following the disk rotation direction. Figure 14 (X-Z plane, top, and Y-Z plane, bottom) shows a more complex scenario with several vertical motions and metallicity gradients. It becomes clear that both low and high metallicity vertical flows are present in the hot-gas halo. It is also clear from this figure that a bimodality in metallicity exists in the vertical direction and also that low-metallic hot-gas flow has a mean motion that brings it from outside to inside the system, i.e., it is falling from the IGM. High-metallic hot-gas flows observed in the southern hemisphere have a slower and more irregular motion than ones in the northern hemisphere. While it is easy to understand and interpret why low-metallic hot gas is falling inside the system, and why in some cases hot metallic gas departs from the stellar disk region (SN feedback, stellar winds, etc.), it is not so obvious why some clouds of such hot metallic gas, with sizes of tens of kpc, behave differently. We suspect that this last high-metallic component could be associated with a small gaseous satellite that is passing through the system at z = 0 (see Figure 15). To confirm all these hypotheses we have created the same velocity-metallicity maps as the ones in Figures 13 and 14 for several snapshots back to z = 0.5. From the analysis of these snapshots we have confirmed that in some instances hot gas with enhanced metallicity departs from the disk due to feedback from the stellar component. We have also observed that low metallicity gas appears to be coming from the CGM. Finally, at around z = 0, we see that a gaseous satellite approaches the main galaxy, perturbing hot-gas flows. The perturbation of the hot-gas distribution by accretion of small satellites is an interesting process and it deserves more detailed study. The case of the Large Magellanic Cloud has been used as a constraint on the hot-gas distribution (Salem et al. 2015). Finally, we have also compared the disk-plane gas mean metallicity with that of the hot gas in the halo and we have found that both components have similar mean metallicities around −0.64 dex (except for the high and low metallicity flows discussed above).
A conclusion from the metallicity and kinematic analysis is the presence of hot-gas flows from the IGM, either from the low metallicity environment or from sinking of metal-enriched satellites. We also conclude that the disk has a low hot-gas concentration as a consequence of SN explosions and stellar winds pushing hot gas up to the halo. Such results are coherent with the Dekel & Birnboim (2006) proposal of cold-gas flows feeding halos with M vir < 10 12 M e .

Halo Virial Mass and the Total Hot Gas in Galaxies
In Table 2 and Figure 16 we present results suggesting a possible correlation between halo virial mass and the total hot gas in galaxies. After analyzing MW-sized models, and also few more massive and some lighter models, we have found that M vir correlates with total hot-gas mass. The blue dots are simulated galaxies that come from a work in preparation. These simulations will be discussed in an extended study in preparation, but we show them to illustrate that the correlation persists regardless of the subgrid physics. Their corresponding host halos were drawn from a 50 Mpc h −1 boxside and chosen to be relatively free of massive companions inside a sphere of 1 Mpc h −1 . A similar result was presented in Crain et al.  (2010). In their work they found instead a correlation between X-ray luminosity and virial mass. More specifically, they reported two linear correlations, one for M vir < 5 × 10 12 and another for M vir > 5 × 10 12 . Our study is mainly focused on MW-sized galaxy systems, i.e., their low halo mass fit, and we have consistently found a simple linear correlation (see Figure 16). Although the absolute values obtained in our work cannot be directly compared with the ones in Crain et al. (2010) as we are using M hotgas and they are using L X , if finally confirmed this linear correlation will provide a new constraint to virial mass in galaxies, including the MW. As mentioned previously, it is challenging to infer the total hot-gas mass in the MW using the available observational data. Matter in the CGM is hot (T = 10 5 -10 7 K) and tenuous (n ;10 −5 −10 −4 cm −3 ). Plasma in these conditions couples with radiation mainly through electronic transitions of elements (C, N, O) in their He-like and H-like ionization stages, the strongest of these transitions falling in the soft X-rays. Due to the extreme low densities and small optical lengths of the CGM, detection of the absorption/emission features produced by this material has proved particularly challenging with the still limited sensitivity of the gratings on board Chandra and   Table 2. Blue dots: simulations outside the GARROTXA run, parameters from these simulations are not presented in Table 2. Figure 17. Spherically averaged hot-gas density (T > 3 × 10 5 K) as a function of radius for the model G.321 at z = 0. XMM-Newton. As a consequence, observations in only a few line-of-sight directions have been useful to observe X-ray absorption in quasar spectra. In this scenario it is clear that simulations can play an important role in the study of hot gas as they allow us to explore the complexities in simulated galactic systems. Including such information enables us to study how observational techniques can lead to biased results for the total hot gas in the galactic system.
With such an aim we have located our mock observer at 8 kpc from the galactic center (i.e., similar to the solar position) in the galactic plane and at a random azimuthal angle, inside our MW-sized simulations. For comparison we refer the reader to the observations of O VII column densities presented in Gupta et al. (2012) and to the total hot-gas values they obtained. In Table 3 we show the results we obtained when making observations in a similar number of randomly distributed directions into the sky as the ones in Gupta et al. (2012). As we do not have realistic chemical species like oxygen in our simulations, we have assumed that we are observing all hot gas, and present the column densities measured in each line-of-sight. As it is common in the literature (e.g., Bregman & Lloyd-Davies 2007;Gupta et al. 2012), we have assumed a spherical uniform gas distribution for the hot-gas halo. As a first attempt, we have estimated, from the column density values in each direction, the total galactic hot-gas mass, assuming a path length of L = 239 kpc = r vir . We have found that using the spherical uniform approach, the total hot-gas mass is systematically overestimated (see Table 3 last column). This is a consequence of the fact that the density distribution in our model is not uniform but decreases as a power law (see Figures 5 and 17). However, this calculation is not very useful, as it requires an a priori knowledge of the CGM path length (L) which is observationally unknown.
A popular technique among observational studies (e.g., Gupta et al. 2012;Bregman & Lloyd-Davies 2007) is to take advantage of the two observables for the hot gas in the CGM: the column density (N H ) and the EM. This technique uses these two quantities to derive the total hot-gas mass with no need to impose an a priori mean density or optical path length. We show results from using this technique in Table 4, both for the optical path length, the mean density, and the total hot-gas mass. In this case, as can be seen in the table, the results show an underestimation of the total hot-gas mass by factors of ∼0.7 −0.1, a result that is independent of the observed direction. Again this result is a consequence of the fact that the real density distribution of our model is not uniform but decreasing with radius as a power law.
As a conclusion we state that using both column density plus imposed optical depth (independent of using the value from a single observation or the mean) or a combination of column density and EM, we are not able to obtain the real hot halo gas mass, when assuming it is isotropically distributed in the galactic halo. Nevertheless, these measurements, particularly the observational method in the literature using both N H and EM, are useful to obtain an order of magnitude estimate of the total hot-gas mass of the halo. Work is in progress to determine Note. We also show the most relevant changes in the initial parameters from the ones used in the model G.321 (see Table 1 for the definition of the parameters). Note. The total hot-gas mass in the simulation, up to r vir , is M vir = 1.2 × 10 10 M e . We have defined hot gas as gas at T > 3 × 10 5 K.
We have assumed L = 230.1 kpc = r vir . the density profile that will allow us to obtain the best hot halo gas mass estimation from observations.

Total Hot-Gas Mass: Accounting for the Missing Baryons
The baryonic fraction in our simulations corresponding only to stellar and cold-gas mass falls in between 0.084 and 0.096. These numbers are far from the cosmic baryon fraction F b U , = 0.171 ± 0.009 (Dunkley et al. 2009;Planck Collaboration et al. 2014), but are comparable to the reported values for isolated spiral galaxies based on weak lensing and stellar mass of 0.056 (Hoekstra et al. 2005). Such a mismatch between the cosmic and galaxy baryon budget suggests that either gas has escaped halos, or has never reached them and is located in the large scale structure WHIM gas (Cen & Ostriker 1999), or that warm-hot-gas embedded in DM halos may account for at least for part of the missing baryons if not all. In our simulations the baryonic fraction reaches (0.107 −0.120) even after adding the hot-gas mass component inside halos. Because we do not include the feedback effect of an AGN, which may eject the gas outside the halo, we can consider our estimates as an upper limit to the halo hot-gas mass. Based on the previous discussion we argue that the missing baryons must be located both, like the hot-gas mass, inside DM halos (e.g., Gupta et al. 2012;Miller & Bregman 2013) and in the IGM, along the filaments and far from the main galactic systems (Rasheed et al. 2010;Eckert et al. 2015). This is consistent with recent studies for CGM in M31 (Lehner et al. 2015) as well as gas stripping in the Large Magellanic Cloud (Salem et al. 2015).

CONCLUSIONS
In this paper we have introduced a new set of MW-sized cosmological N-body plus hydrodynamics simulations. Their resolution and realism are high enough to allow us to study stellar kinematics (Roca-Fàbrega et al. 2013 and the evolution of the large scale structures of these kinds of galaxies. Most stellar disk parameters in our G.3 models are in agreement with the observational ranges for the MW. Disk scale length and scale height as well as their dependence on the stellar age are reproduced. The simulated disk shows a large variety of large scale structures such as spiral arms, a ring, and a bar. It also shows a flared and a warped configuration. The quality of our runs is comparable to the most recent works presented by Guedes et al. (2011) and Mollitor et al. (2015), a detailed comparison can be seen in Table 1. Our simulations have a similar or a larger number of resolution elements (stellar particles or gas cells), higher spatial resolution, and a smaller mass per particle than in the above works. In addition, as we have shown in Figure 4, our circular velocity curves are in very good agreement with recent observations. On the other hand, a spheroid (bulge+halo) that is 2−3 times more massive than the one observed in the real MW is present in our G.3 models. Since the density of this spheroidal component is relatively low in the central region of our galactic systems it does not produce centrally peaked circular velocity curves.
We have also analyzed the SFH in our simulations. We have found that our SFR values at z = 0 are lower than recent observational values. Also the SFH peak in our G.3 models occurs at higher redshifts than the one predicted by Behroozi et al. (2013). In the most recent and realistic MW-sized simulations (Guedes et al. 2011;Bird et al. 2013;Mollitor et al. 2015), the predicted SFH is also not well reproduced. We have also found that the total gas mass and its distribution in temperature slightly differ from MW observations. This mismatch in SFH and gas mass distribution is usual in simulations due to a lack of understanding of subgrid physics and to the difficulty of accounting for all physical processes involved in the formation and evolution of galaxies. Studies like ours will help to constrain such subgrid models. Recently, it has been proposed that the inclusion of kinematic feedback might solve the early SF issue and thus change the position of the SFH peak.
The work presented here is focused on the study of the distribution and the amount of potentially X-ray emitting gas (T > 3 × 10 5 K) in our simulations. The comparison of the properties of the hot gas in our simulation with the observations of the hot halo X-ray corona of the MW and external galaxies was made by producing hot-gas column density mock observations. Our conclusions are as follows: 1. Hot gas is not distributed homogeneously around the galactic halo. 2. There is no evidence of a hot-gas thick disk in our models. 3. By using the real mass of warm-hot gas in our simulations we conclude that about 53%-80% of the missing baryons can be accounted for by hot gas in the halo corona. The exact value depends on the simulated galactic system (i.e., halo mass). We argue that the missing baryons are located in filaments (IGM), not inside main galactic systems, as preheated warm-hot gas (Rasheed et al. 2010). It is important to mention that in our simulations AGN feedback from a central galactic black hole has not been implemented. It is still open to debate how a relatively low-mass black hole, such as the one in the real MW, could affect the mass and temperature distribution of the hot halo gas and thus the SFH. 4. The hot-gas density follows several streams in the halo that can be detected because of their changes in kinematics and metallicity. We have seen how low metallicity streams come from the IGM while others with higher metallicity depart from the disk. We have also detected the infall of high metallicity gas that can be interpreted as coming from a satellite accretion. Note. Using EM and N H values we have computed the path length (L) and hotgas volume density (n). Finally we show the total hot-gas mass obtained when assuming a sphere with constant density (n) and its radius equal to the path length (L). The real total hot-gas mass in the simulation is M vir = 1.2 × 10 10 M e .
5. The total hot-gas mass in the halo cannot be recovered when assuming a spherical uniform gas distribution. In our model the real hot-gas density distribution leads to an overestimation when using the mean column density plus a fixed optical depth under the spherical uniform assumption. On the other hand we find an underestimation of the total hot-gas mass when we use the column density plus the EM when using the same assumption. 6. Following previous works we have found that a clear relation exists between total hot-gas mass and the virial mass of halos. This result is important as it becomes a new method to constrain the total mass of galactic systems, such as the MW.
We thank A. Klypin and A. Kravtsov for providing the numerical codes. We thank the HPCC project, T. Quinn, and N. Katz for the implementation of the TIPSY package. This work was supported by the MINECO (Spanish Ministry of Economy)-FEDER through grants AYA2012-39551-C02-01 and ESP2013-48318-C2-1-R, and Conacyt Fronteras de la Ciencia through grant 281. The simulations were performed at the Supercomputer Miztli of the DGTIC-UNAM, and at  Table 5 for more details about the models). Black: total; red: DM; blue: stars; and green: gas. Note. All parameters shown here have been well described in Table 1.
Atocatl and Abassi2, HPC facility of the Institute of Astronomy-UNAM.

APPENDIX A RESOLUTION STUDY
In this appendix we show a test of numerical convergence using a low-resolution re-simulation of our main model. The re-simulation presented here (G.321_lr) has been obtained using the same initial conditions and initial parameters as in G.321, but with a resolution of 218 pc instead of 109, and ∼9 × 10 5 DM particles inside R v instead of ∼7 × 10 6 . In Figure 18 we show a comparison of the total DM, stellar, and gas circular velocity curves computed using ( ) GM r as a proxy for V c . We see that the resolution does not affect the total circular velocity curve profiles. However, it is interesting to note that the circular velocity curve from the G.321 gas component shows significantly higher V c values at large radii. This result is in agreement with the results of Scannapieco et al. (2012); i.e., the properties of the gaseous component seem to be the most sensitive to numerical resolution effects. We show in Table 5 the values of the main parameters of the G.321 run and compare them with the corresponding ones of the G.321_lr model. We see that, aside from the mass of cold gas and related quantities, like SFR at z = 0, most parameters do not differ from each other by more than 30%. This was also found in Scannapieco et al. (2012). From this result, we conclude that although, in general, parameters are not sensitive to changes in resolution, caution must be exercised when analyzing the coldgas component and the SFR. In addition, we notice that the cold-gas fraction is not only sensitive to numerical resolution (the maximum refinement level), but also to changes in terms of the number of cells (different aggressive refinement, see Section 3.3). On the other hand, the difference observed in the G.321 gas velocity curve at large radii is simply due to the fact that this run ends up with much more gas. Additional studies of numerical convergence of the code can be found in Avila-Reese et al. (2011) andGonzález-Samaniego et al. (2014a).

APPENDIX B JEANS LENGTH
In this section, following the condition proposed by Truelove et al. (1997), we show that in our simulations the resolution is above the minimum that ensures that the formation of stellar particles is of physical origin rather than numerical. In order to see this we have computed the local Jeans wavelength of the cold gas (λ J ) according to Equation (1). In Figure 19 we show a histogram of the Jeans length in units of the size of the corresponding cell for simulation G322. Except for a few cells, the cells satisfy the Truelove criterion.