Not gone with the Wind: Survival of High-Velocity Molecular Clouds in the Galactic center

High-velocity atomic clouds in the Galactic center have attracted significant attention due to their enigmatic formation process, which is potentially linked to the starburst or supermassive black hole activities in the region. Further, the discovery of high-velocity molecular clouds (HVMCs) presents a greater puzzle, because they are much denser and more massive. If the HVMCs were accelerated by the strong activities in the Galactic center, they are expected to be destroyed before they reach such a high velocity. To shed light on this phenomenon, we perform three-dimensional numerical simulations to investigate the origin and hydrodynamic evolution of HVMCs during a starburst in the Galactic center. We find that the presence of a magnetic field provides effective protection and acceleration to molecular clouds (MCs) within the galactic winds. Consequently, the MCs can attain latitudes of approximately 1 kpc with velocities around 200 km/s, consistent with the observed characteristics of HVMCs. The consistency of our findings across a wide parameter space supports the conclusion that HVMCs can indeed withstand the starburst environment in the Galactic center, providing valuable insights into their survival mechanisms.


INTRODUCTION
Galactic feedback, especially the nuclear wind, is now commonly accepted as an important process affecting the galactic evolution (e.g.Fabian 2012;Heckman & Best 2014;Heckman & Thompson 2017;Naab & Ostriker 2017;Zhang 2018, and references therein), which is, however, pretty weak in our Milky Way at present (Baganoff et al. 2003;Haywood et al. 2016).Therefore, it is expected that Milky Way had been active before, but quenched after that, which should produce some corresponding relics.Over the past tens of years, these feedback relics possibly have been discovered at radio, X-ray and -ray band, such as the Galactic Center Lobe (GCL; Sofue & Handa 1984), the microwave haze (Finkbeiner 2004;Planck Collaboration et al. 2013), the polarized lobes (Carretti et al. 2013), the Fermi bubbles (Su et al. 2010), the radio bubbles (Heywood et al. 2019), the X-ray chimneys (Ponti et al. 2019) and the eROSITA bubbles (Predehl et al. 2020).These structures have scales ranging from ∼100 pc to ∼10 kpc, indicating that they originated from a series of violent activities.In addition, in the Galactic center, many highvelocity clouds (HVCs) were detected both above and below the Galactic plane (Collins et al. 2004(Collins et al. , 2005;;Di Teodoro et al. 2018;Lockman et al. 2020;Ashley et al. 2020).Especially, two highvelocity molecular clouds (HVMCs) are also discovered insides the HVCs (Di Teodoro et al. 2020), The altitudes of the two MCs are 0.6 and 0.9 kpc, respectively.Their velocities along -axis are ∼ 180 and 150 km s−1, while the radial velocities are ∼ 240 and 300 km s −1 based on a biconical model (Di Teodoro et al. 2020).
★ E-mail: murphychang@zju.edu.cn† E-mail: miaoli@zju.edu.cn Their molecular mass are both ∼ 380 M , and their atomic mass are 220 and 800 M .The HVMCs show good coincidence with some aforementioned relics, so they possibly originate from similar process, e.g., accelerated by the Galactic nuclear wind.Although these relics and HVCs/HVMCs are commonly suggested to be produced by the feedback activity, the detailed mechanism is still the subject of intense debate.Several models have been proposed to explain their formation, some of which focus on one structure (Crocker & Aharonian 2011;Zubovas et al. 2011;Guo & Mathews 2012;Zubovas & Nayakshin 2012;Fujita et al. 2013;Mou et al. 2014;Fujita et al. 2014;Lacki 2014;Mou et al. 2015;Sarkar et al. 2015;Zhang & Guo 2020), while others attempt to simultaneously explain multiple structures (Yang et al. 2013;Crocker et al. 2015;Yang & Ruszkowski 2017;Zhang et al. 2021;Yang et al. 2022).Most of these models exhibit self-consistency, and some following simulations have provided further validation of their viability (Guo & Mathews 2012;Mou et al. 2014Mou et al. , 2015;;Sarkar et al. 2015;Zhang & Guo 2020;Yang et al. 2022, e.g.).However, simulating the acceleration of HVMCs still presents a challenge for their formation models.Compared to atomic clouds, molecular clouds are denser and cooler, making it more difficult to accelerate them to high velocity without disruption (Schneider & Robertson 2017;Cashman et al. 2021).Some simulations for extra-galactic interaction between clouds and nuclear winds show that clouds can be protected by magnetic field (Alūzas et al. 2014;McCourt et al. 2015;Banda-Barragán et al. 2016a;Zhang et al. 2017;Sparre et al. 2020;Jung et al. 2023), cooling (Gronke & Oh 2018, 2020;Kanjilal et al. 2021) and thermal conduction (Armillotta et al. 2017), which confirms that cool clouds can survive acceleration by a hot wind.Nevertheless, these simulations usually involve a constant hot wind, which is completely different from the unpredictable nuclear wind produced by starburst or AGN.Moreover, most of them focus on high-latitude atomic or even ionized clouds, so they cannot clearly explain the formation of HVMCs at ∼ 1 kpc in our Milky Way.It is therefore necessary to perform robust simulations to see whether the HVMCs observed in the Galactic center can be reproduced.
The formation of HVMCs is closely linked to the other feedback relics, and could potentially be used to distinguish among different models for their origin.While the activity of active galactic nuclei (AGN) has the capability to accelerate molecular clouds (MCs) to high velocities, it is often so powerful that the clouds usually diffuse to atomic/ionized form.Unless, there are some periodic bursts, such as those arising from accretion onto the supermassive black hole (SMBH) Sgr A* (Wang et al. 2013).These bursts should be weaker than normal AGN, but still release comparable amounts of energy, allowing the MCs to be efficiently accelerated without being quickly destroyed.Relatively speaking, a starburst is a more feasible explanation for the formation of HVMCs, as the Galactic center exhibited a higher star formation rate about 30 million years ago (Nogueras-Lara et al. 2020) and the molecular outflow is universal in active star-forming galaxies (Spilker et al. 2018;Roberts-Borsani et al. 2020;Spilker et al. 2020;Stuber et al. 2021;Butler et al. 2023).In fact, although the supernova feedback is important in galaxy formation (Kim & Ostriker 2015;Martizzi et al. 2016;Li et al. 2017;Hu 2019), its working mechanism has not been fully understood, which leads to a difficulty to understand the role of HVMCs and also limits the cosmological simulations (Li & Bryan 2020).Currently, it is known that randomly distributed SNe in the disk only drive inefficient galactic winds because most supernova remnants lose their energy radiatively before breaking out of the disc (Fielding et al. 2018), leading to a difficulty to push the HVMCs to high latitude in such a galactic wind.Nevertheless, a starburst in the Galactic center can produce much stronger galactic wind and more efficiently accelerate the HVMCs.It is expected that the starburst ended recently, but left these feedback relics and HVMCs.It is difficult to tell which model is correct, because the hydrodynamical evolution of AGN and starburst activity can be similar at large scale.Their energy input rate can be similar, as a result, the wind driven by these activities can reach a comparable velocity at high latitude.However, there should be noticeable differences at smaller scale (≤ 1 kpc), such as the acceleration process of HVMCs and the morphology of relics, because the starburst can happen more randomly in a much larger region than the AGN activity.Moreover, the metallicity of HVMCs is possibly different for AGN and starburst models, since starburst can produce more heavy elements.Although Ashley et al. (2022) indeed found different metallicitiy distribution of HVCs in the Galactic center, they explain that HVCs originate in Milky Way's disk and halo.These models can be further examined through simulations.
In this paper, we investigate whether HVMCs observed in our Milky Way can be accelerated to high latitudes by a starburst.To this end, we perform a detailed simulation of the process.We start by simulating a series of random core-collapse supernova explosions in the Galactic center with a frequency estimated based on a past star formation rate (Nogueras-Lara et al. 2020).We set a molecular cloud above the explosion region to study how the cloud is accelerated by the outflow wind and whether it can survive until it reaches 1 kpc, a position similar to the clouds detected by Di Teodoro et al. (2020).The explosion region is believed to be adjacent to the central molecular zone (CMZ), where more giant molecular clouds steadily exist.Next, we check the density, temperature, and velocity of the clouds obtained from the simulations and modify initial conditions to study the influence of different parameters.We will try to identify various HVMCs candidates obtained from the simulation by comparing with the observation and study their properties in detail.Finally, we investigate the mixture of clouds and the ejecta of supernovae to disentangle the metallicity in HVMCs.This paper will describe the simulation setup in Section 2 and show the results in Section 3. The formation of the HVMCs, their metallicity and their relation with feedback relics will be discussed in Section 4. The Section 5 is a summary.

SIMULATION
To perform the simulations, we utilized the publicly available, modular magnetohydrodynamic (MHD) code PLUTO1 (Mignone et al. 2007(Mignone et al. , 2012) ) to perform the simulations.This grid-based MHD code employs a second-order Runge-Kutta time integrator and a Harten-Lax-van Leer Riemann solver for middle contact discontinuities, making it well-suited for simulating the interaction between the SN shock and the molecular clouds.

Basic configuration
The simulation is based on a three-dimensional (3D) MHD cartesian frame with a grid of 200×200×2000, equivalent to a physical volume of 100 × 100 × 1000 pc 3 and a linear resolution of 0.5 pc pixel −1 .We set the -axis to be perpendicular to the Galactic disk (north as positive), the -axis to run along decreasing Galactic longitude, and the -axis to be parallel to the line-of-sight (the observer at the negative side).We adopted an outflow boundary condition for all directions, which means that some of the clouds' material may flow outside of the simulation box.
The simulation is governed by the ideal MHD conservation equations, where  is the mass density,  the thermal pressure, v the velocity, B the magnetic field, 1 the dyadic tensor, Φ the gravitational potential, and   the total energy density, defined as: where  is the internal energy.We use an ideal equation of state, i.e.,  = /(Γ − 1), in which the ratio of specific heats Γ = 5/3.
To accurately model the gravitational potential in the simulation volume, we assume that it is static and fully determined by the SMBH, the nuclear star cluster (NSC), and the nuclear disk (ND).A point mass of 4 × 10 6 M is taken to represent the SMBH.For the NSC and the ND, we adopt a spherical distribution following Chatzopoulos et al. (2015, Equation 5 therein).To incorporate radiative cooling in the simulation, we use a piece-wise cooling function with a lower limit of the cooling temperature set to 100 K.We assume a solar abundance (H abundance  =0.711,He abundance  =0.2741, metallicity  =0.0149) for the ISM and the initial MC.The multiphase gas in the Galactic center includes hot ionized (∼ 10 6 K) (Kataoka et al. 2013;Ponti et al. 2019), warm ionized (10 4 to 10 5 K) (Fox et al. 2015;Bordoloi et al. 2017) and cool atomic (10 3 to 10 4 K) gas (McClure-Griffiths et al. 2013;Di Teodoro et al. 2018), etc., in which the gas lower than 100 K is usually taken as molecular gas.In the simulation, temperatures below 100 K are typically not due to cooling, but rather due to adiabatic expansion.

Supernova explosion and molecular clouds
The initial conditions for our simulations are based on both observations and analytical models.Observationally, the high-velocity molecular clouds HVMCs typically exhibit densities ranging from 10 to 300 cm −3 and outflow velocities between 200 and 300 km s −1 (Di Teodoro et al. 2020).However, to account for the significant gas loss that occurs during their propagation, we assume that the initial densities of the MCs should be higher.In addition, we need to consider other parameters such as the supernova explosion frequency and the initial latitude of the MCs to ensure that they reach the observed velocities without being completely destroyed.Therefore, we perform a systematic exploration of the parameter space to identify the most plausible initial conditions for our simulations.Here, we introduce the cloud crushing time, to quantify the timescale of cloud crushing (Klein et al. 1994), in which  mc is the radius of the initial molecular cloud,  mc the density of the cloud,  sn the wind velocity produced by supernovae, and  sn the wind density.Based on general understanding, a cloud should begin to crush when the evolution time is longer than  cc , and should totally crush after a period of 2 cc .However, this estimation does not take into account the effects of the magnetic field and cooling mechanisms, which can play an important role in the cloud's evolution.
In a cylindrical region with a radius of 35 pc and a height of 10 pc, the fiducial SN birth rate is set to be 10 kyr −1 (Di Teodoro et al. 2018), which is estimated by assuming an SFR of 1 M yr −1 , a Kroupa (2001) initial mass function (IMF) and a minimum mass of 8 M for the progenitor star of a core-collapse SN.The center of the cylindrical region is set to be located at the western 100 pc of Sgr A*.Barnes et al. (2017) and Sormani et al. (2020) estimated a current SFR of 0.1 M yr −1 inside the CMZ, while Nogueras-Lara et al. (2020) found that star formation in the ND (which has a similar radial extent as the CMZ) has been relatively active in the past 30 Myr, with an SFR of 0.2 − 0.8 M yr −1 .Our assumed SFR of 1 M yr −1 is compatible with a local starburst, which may be the case if SN events have been episodic and clustered on a Myr timescale.This SFR is actually larger than the typical value in such a small region, so we also test a run with lower SN birth rate of 5 kyr −1 .We have neglected Type Ia SNe, which have a birth rate of 0.05 kyr −1 according to the enclosed stellar mass in the ND/NSC (Mannucci et al. 2005).The SNe are set to randomly explode in the cylindrical region, and we use same random seed in all runs.
The density of the MCs follows an inverse square law, n mc = n 0 /r 2 , in which n 0 is the central density, r the radius.In the fiducial simulation, n 0 = 1500 H cm −3 , the maximum radius of the initial MC is 10 pc, and the height of the MC from the Galactic plane is 50 pc.Based on these settings (r mc = 10 pc, v sn = 1000 km s −1 , n mc = 15∼50 H cm −3 , n sn = 0.01 H cm −3 ), we can estimate the  cc 1∼2 Myr, so the cloud will totally crush after 4 Myr in the classical analysis.However, in our preliminary tests, we find the cloud can survive beyond 7 Myr by including a vertical magnetic field and the cooling effect.In this scenario, after around 7 Myr, the cloud will run outside of the simulation box.Therefore, the simulation results are presented up until around 7 Myr.
In addition, once injected, the ejecta will eventually partially mix with the molecular clouds, and change their metallicity.To study the mixture, we introduce two tracer parameters,  1 and  2 , which are both evaluated at each pixel in the simulation and obey a simple conservation law: 1 has a value of 1 for pure SN ejecta and 0 for the unpolluted molecular clouds and ISM, while  2 has a value of 1 for pure molecular clouds and 0 for the unpolluted SN ejecta and ISM.The values in between indicate a mixed gas.These tracer parameters allow us to track the mixing process over time and analyze the distribution of metals in the simulated system.

The ISM and the magnetic field
We initialize our simulation with a uniform distribution of ISM density and temperature, with values of 0.01 H cm −3 and 10 6 K, respectively, over the entire simulation box.Although thermal pressure is expected to be higher at lower latitudes due to rough hydrostatic equilibrium against gravity, our preliminary tests suggest that this effect is unimportant since the shock wave from the supernovae breaks this equilibrium early on.Moreover, the stellar wind in the Galactic center is also strong and can unremittingly break this equilibrium.
The distribution of magnetic fields in the Galactic center remains a challenging problem, particularly in the central tens of parsecs (Ferrière 2009), with many different components, influencing the strength and direction of the magnetic field.There is actually a general model for the whole Milky Way (Beck 2013;Cerri et al. 2017), in which the magnetic field is parallel to the Galactic plane at lower latitude and gradually tend to be perpendicular at higher latitude, but this is only an approximation in the Galactic center.Therefore, we in this work test different runs, respectively with parallel, perpendicular and no magnetic field.
The magnetic strength range from ∼ 1 mG in the central tens of parsecs (Ferrière 2009) to few G at 1 kpc above the Galactic plane (Cerri et al. 2017).For simplicity, we adopt a homogeneous magnetic strength of 10 G over the whole simulation box.The initial parameters are summarized in Table 1.

RESULTS
In this section, we present the simulation results.We first describe in detail the evolution of the MCs in the vertical magnetic field in the fiducial run (Section 3.1) .We then examine the role of the magnetic field in the two additional runs, one with horizontal magnetic field (Section 3.2) and the other with no magnetic field (Section 3.3), to illustrate how the change affects the formation of the HVMCs.Finally, we study the influence of the cloud density and the supernovae explosion frequency (Section 3.4).
To quantitatively compare with the observation, we here parameterize the main features of the observed MCs, MW-C1 and MW-C2 (Di Teodoro et al. 2020) .The altitudes of the two MCs are 0.6 and 0.9 kpc, respectively, so we choose ∼ 1 kpc as the standard position to guarantee the simulated clouds can indeed reach the height.Their  3) ISM H density, in units of cm −3 .( 4) The MC central H density, in units of cm −3 .( 5) The MC radius, in units of pc. ( 6) The MC height, in units of pc. ( 7) The radius of the cylindrical explosion region, in units of pc. ( 8) The height of the cylindrical explosion region, in units of pc. ( 9) The direction of the magnetic field along the Galactic plane.
velocities along -axis are ∼ 180 and 150 km s−1, while the radial velocities are ∼ 240 and 300 km s −1 based on a biconical model (Di Teodoro et al. 2020).Our simulation focuses on the propagation vertical to the Galactic plane, so we take 200 +100 −50 km s −1 as the typical value.The molecular mass of MW-C1 and MW-C2 are both ∼ 380 M , but their atomic mass are 220 and 800 M .Thus we pay more attention to match the molecular mass, and the atomic mass can vary in a large range.With a diameter of ∼30 pc, their mean molecular number densities are 130 and 190 H 2 cm −3 , and the mean atomic number densities are 1 and 3 H cm −3 .In the work, we take the clouds denser than 10 H cm −3 as MCs, and the clouds with a density between 1 and 10 H cm −3 as atomic clouds.
In addition, there are also some qualitative features which are worth reproducing.Surrounding the HVMCs, there are always some atomic clouds with lower density and larger volume, which were usually taken as HVCs before the discovery of HVMCs.The number of detected HVCs is much larger than HVMCs, and most of HVCs are uniformly distributed above 250 pc (Di Teodoro et al. 2018).There are possibly more HVMCs hidden in the HVCs, so more highresolution and high-sensitivity molecular observations are necessary.

The run for the fiducial set
The column density and velocity evolution of f100n1500v is shown in Figure 1, in which the clouds will reach 1 kpc at 7 Myr.In the following text, we call the time, at which the results well match the observation, as the fiducial time.At the early stage, the supernova shock wave would blow the initial MC to be a thin filament, because the central density of the cloud was much higher than the boundary.The filamentary structures have also been investigated by Banda-Barragán et al. (2016b); Jung et al. (2023), who claim the filaments are only formed in magnetized environment and the cloud will crush to small clumps without magnetic field, consistent with our results.When the peripheral low-density material was blown to higher latitude, the central dense core was being slowly accelerated.Some pioneer high-velocity clumps broke away from the main cloud at 3 Myr, and run outside of the simulation box at 4 Myr.At this stage, the main cloud became more irregular, but kept as one cluster.After 7 Myr, the cloud would reach 1 kpc, a position consistent with the observation.During the propagation of the cloud, the supernovae shock was always being reflected by the cloud and gradually produced stronger reverse shock.This process leads to the obvious dividing line both for the density and velocity at 7 Myr.The reverse shock could roughly balance the forward shock, as a result, the cloud acceleration rate largely decreased.We also show the temperature and magnetic field evolution in Figure 2. The outflow wind interact with the MCs, heating the surrounding ISM and compressing the magnetic field, while the central cores of the clouds still contain cool gas and low magnetic field at the early stage.The shock wave from the supernovae can sweep the whole simulation box at ∼ 1 Myr and heat the ISM to high temperature.However, with the receding of a part of the outflow wind at higher latitude, the magnetic field becomes much weaker.
Figure 1 & 2 also show the starburst wind is not constant, especially at low latitude, because we adopt the random supernovae explosions in the simulations.The ever-changing wind will significantly influence the evolution of the initial cloud.However, the variation of the starburst wind is small at high latitude, where it can be taken as a constant wind.
To study whether the clouds at 7 Myr can be still taken as MCs with a velocity of ∼ 200 km −1 , we show the density-velocity and density-temperature maps in Figure 3 2017), but we replace their constant wind with the simulated starburst wind.As a result, Figure 4 contains the SNR-dominated region, i.e, the upper left band, which is absent in their work.The clouds selected based on criteria of n≥ 1 H cm −3 and T≤ 10 4 K have a total mass of ∼ 1500 M , while those selected based on criteria of n≥ 10 H cm −3 , T≤ 200 K, and  ≥800 pc are taken as molecular clouds and have a molecular mass of ∼ 850 M .However, these clouds cover a region larger than the MW-C1 and MW-C1, and we should compare parameters at same scale.If we choose the densest central clouds (diameter ∼30 pc, i.e., 60 cells) as the counterpart, the mass can better match the observation.The clustering of the clouds is also considered in the estimation, in which some cells with appropriate density and temperature will still be excluded, if there is not any cloud cell within the surrounding 0.5 pc.We also estimate the present mass-weighted mean velocity of ∼ 190 km s −1 for all clouds (n≥ 1 H cm −3 and T≤ 10 4 K), while the mean velocity over the past 7 Myr is ∼ 130 km s −1 , both a little lower than the observation.
In the vertical magnetic field, the clouds can propagate to 1 kpc without destruction, and still keep a considerable mass even larger than the observed HVMCs.However, the mean velocity is a little smaller than the typical value.To increase the velocity, a straightforward method is to increase the supernovae explosion frequency, but the frequency used in our work is already a little higher than MNRAS 000, 1-?? (0000)  the standard value.In addition, it is unexpected that including the horizontal magnetic field can also increase the velocity, which will be illustrated in the next section.Assuming a lower MCs or ISM density is also practical, so we test a case with a lower density of the initial MC in Section 3.4.In summary, the fiducial run can indeed explain the acceleration of MCs at high latitude, while some features cannot be reproduced perfectly.

The run with horizontal magnetic field
We show the column density evolution of f100n1500h in Figure 5, while the density-velocity distribution at 5 Myr is shown in Figure 6.Similar to f100n1500v, the MC was blown to be a thin filament initially, but gradually some gas was stripped.At 2 Myr, a pioneer high-velocity clump separated from the main cloud, but run outside the simulation box at 3 Myr.With the gas stripping, the MC showed a more irregular shape and finally crushed to several clumps.These clumps have lower densities and higher velocities, but can be still taken as molecular clouds.Especially, there is the second high-velocity clump separating from the main cloud after 4 Myr and reaching ∼ 1 kpc after 5 Myr, with a mass-weighted mean velocity of ∼ 340 km s −1 , a total mass of ∼ 700 M (n≥ 1 H cm −3 ) and a molecular mass of ∼ 100 M (n≥ 10 H cm −3 ).The velocity is higher, but the masses are both lower than the observation's.
By comparing with f100n1500v, we find a horizontal magnetic field can stimulate the acceleration and the crushing of the MCs, which is possibly caused by the magnetic tension force vertical to the Galactic plane, i.e., the magnetic draping, a ubiquitous mechanism already found in the launching of clouds (Cottle et al. 2020).The outflow wind can compress the MCs and the surrounding ISM, then amplify the local magnetic field, i.e., the magnetic tension.The magnetic field can help to efficiently accelerate the MCs, while some MCs material will flow along the magnetic field, even run outside of the simulation box.As a result, the MCs can be pushed to high velocity at high latitude, but lost much mass.In addition, if the magnetic field includes more horizontal components, the clouds can be further dispersed at large scale, which can produce some smaller clouds than those in f100n1500v.These clouds may be more similar to the observed MW-C1 and MW-C2.
In fact, the magnetic field in the Galactic center is complicated, while the vertical component is more important (Beck 2013;Cerri et al. 2017).At present, there is no a standard magnetic field model, so we test the two runs to study the influence of the magnetic direction on the simulation.In terms of the two runs, a mixed magnetic field would likely better explain the observed properties of the HVMCs.

The run without magnetic filed
We show the column density evolution of f100n1500n in Figure 7, and the density-velocity distribution after 5 Myr in Figure 8.There are no large clumps separation.Instead, lots of small clumps gradually diffuse from the main cloud, consistent with the simulation results of Schneider & Robertson (2017).The main cloud is slower and will be depleted after 5 Myr, roughly consistent with the crushing time estimation, as a result, the MCs will not reach 1 kpc.In other words, in comparison with f100n1500v and f100n1500h, the the magnetic field can indeed well protect the clouds.
However, Figure 7 illustrates the densest regions are significantly denser, and Figure 8 also shows there are more high-density clouds ( n ≥ 1000 H 2 cm −3 ) than the runs with magnetic field, which indicates the magnetic field stimulates the destruction of the high-density clouds.This effect may be attributed to the increased turbulence resulting from the presence of the magnetic field, facilitating a more efficient mixing of the MCs and interstellar medium (ISM).As a consequence, high-density clouds share material with low-density regions.Additionally, the clumps grow larger and exhibit prolonged survival but possess lower densities.Consequently, the local density of the clouds is diminished in the presence of a magnetic field.
In conclusion, the magnetic field plays a crucial role in the formation of HVMCs.However, it should be noted that the magnetic field does not always increase the density of MCs and can also disperse some of the densest MCs at smaller scales.

The run with lower density and lower explosion frequency
We show the results of f100n1000v in Figure 9 and Figure 10.The initial cloud was also blown to a filament, but a little wider than previous runs.A large amount of gas were stripped at 2 Myr, and dissipated at 3 Myr.Then the main cloud was divided to two clouds, in which the faster one almost reached 1 kpc at 5 Myr, but the left one gradually disappeared.After 5 Myr, the simulation box has a total mass of ∼ 1100 M (n≥ 1 H cm −3 , T≤ 10 4 K) and a molecular mass of ∼ 500 M (n≥ 10 H cm −3 , T≤ 200 K,  ≥800 pc), roughly consistent with the MW-C2.The mass-weighted mean velocity is ∼ 290 km s −1 , also similar to the observation.However, the diameter of the whole cloud is larger than the observation, so the density is lower.
By comparing with f100n1500v, we can estimate a central density between 1000 and 1500 cm −3 for the initial cloud.Meanwhile, if the magnetic field includes more horizontal components, the clouds can be further dispersed to some small clouds.In other words, if f100n1500v uses a lower central density and more horizontal magnetic field, it can better match the observation.However, the primary focus of this study is to investigate whether the MCs can be accelerated to high velocities at high latitudes, and the current findings adequately address this inquiry.Moreover, it is important to note that the parameters for MW-C1 and MW-C2 are only approximations derived from a simplified biconical wind model, and the completeness of the HVMCs sample remains uncertain.There are only two detected HVMCs in the Galactic center, so the main feature of HVMCs is actually still ambiguous.As a result, conducting an exhaustive search of the parameter space is unnecessary at this stage.
The results of f200n1000v are shown in Figure 11 and Figure 12.The gas were gradually stripped, but the cloud can still survive to reach 1 kpc at 8 Myr.At 8 Myr, the simulation box has a total clouds mass of ∼ 640 M ( n≥ 1 H cm −3 , T≤ 10 4 K), a molecular mass of ∼ 400 M (n≥ 10 H cm −3 , T≤ 200 K,  ≥800 pc) and a massweighted mean velocity is ∼ 180 km s −1 , roughly consistent with MW-C1, which indicates the cloud can also well survive, even if the explosion frequency of the supernovae is lower than 10 kyr −1 .In fact, the explosion frequency should vary with the evolution of the cloud, and the features of resultant clouds are also dependent on the variation.
We summarize all results in Table 2 and visualize it in Figure 13, which will be further discussed in Section 4. The criteria (n≥ 10 H cm −3 , T≤ 200 K,  ≥800 pc) used to choose the molecular components is not always reasonable, and some hydrogen atoms can also survive on such a criteria.We here show the results with a strict criteria (n≥ 100 H cm −3 , T≤ 150 K,  ≥800 pc) for an error estimation.It is actually difficult to accurately estimate the realistic velocity of the clouds based on the current observation.We take the velocity along the sightlines as the lower limit, and the outflow velocity as the standard velocity.The outflow velocity is estimated based on a biconical model, which is also an important reference for our simulations, so we use the outflow velocity to directly compare with the simulations.

DISCUSSION
In the preceding sections, we have presented 3D simulations that illustrate the long-term hydrodynamic evolution of MCs propelled by subsequent supernova explosions.These simulations incorporate simplified, yet sufficiently realistic physical conditions of both the MCs and the surrounding environment.The first three simulation runs, which represent the evolution with vertical, horizontal, and no magnetic field, exhibit varying degrees of success and shortcomings in replicating the primary observed characteristics of MW-C1 and MW-C2.The last two runs show the simulations work well in a wide parameter space.In this section, we analyze the outcomes of these simulations and discuss their implications for our comprehension of the enigmatic ecosystem in the Galactic center.

Formation and evolution of the HVMCs
Table 2 demonstrates that the total mass of the four runs with magnetic field aligns with the observed HVMCs, while f100n1500h exhibits a lower molecular mass and higher velocity.Figure 13 clearly indicates that f100n1000v provides the closest match to the two HVMCs, though the other two runs also show rough consistency with the observations.However, the key point we want to make is that the HVMCs can indeed be accelerated to high velocities without disruption, which is also reflected by the other three runs.We can assume the position of the central point shown in Figure 13 can be interpolated accordingly, if we change one of the parameters, such as the direction of the magnetic field, the supernovae explosion frequency and the density of the initial cloud, based on which we can roughly estimate the dependence of their positions on these parameters.For example, by comparing f100n1500v with f100n1500h and drawing a line between the two central points, we can expect a point will be located between the two MCs, when we only modify the direction of the magnetic field.Similarly, we can get a cloud with higher density and velocity than the two MCs by properly increase the supernovae frequency or decrease the initial cloud density.
All four runs with magnetic field can reproduce the HVMCs, so the HVMCs can be indeed formed by the acceleration of the starburst in the Galactic center.On the other hand, these results indicate the magnetic field is important and the MCs can well survive the shock of supernovae even at a scale of ∼ 1 kpc.Zhang et al. (2017) claim the cold gas with temperature of 10 2 ∼ 10 4 K cannot survive a hot Galactic wind, but they neglect the magnetic field and pay more attention to study the process at larger scale, which will not conflict   with our results.Of course, in our simulations, there are also some features inconsistent with the observations, so we will try to clarify them in this section.
To better study the evolution of the HVMCs, we show the mass evolution of all runs in Figure 14.The criterion for distinguishing between various components follows the description presented in Section 3.1.There are three kinds of mass, the total gas mass (atoms + molecules), the total molecular mass and the mass of molecular gas with a latitude higher than 800 pc.For simplicity, we take the last one as the molecular mass of the HVMCs.
This analysis of the mass evolution of the HVMCs shows that the total gas and molecular mass of all five runs gradually decrease over time due to ionization, stripping by the hot wind and outflows from the simulation box.However, f100n1500h and f100n1000v show a rapid declination respectively after 3.5 Myr and 5.5 Myr.For f100n1500h, this is caused by the dissipation of the pioneer high-velocity clump which quickly diffuse and run out of the left and right edges of the simulation box along the horizontal magnetic field.Similar to f100n1500h, f100n1000v also has a high-velocity clump running out of the simulation box, but from the upper edge after 5.5 Myr.As for the total molecular clouds, they are stripped and dissociated rapidly at the beginning, and maintain a steady decrement.At last, f100n1500h and f100n1000v lost most of molecular gas after 6 Myr, while a large amount still survive in f100n1500v, f100n1500n and f100n1000v.In f100n1500n, the clouds, almost totally crushed after 5 Myr, cannot approach 800 pc, so they are impossible to form the HVMCs.
In f100n1500v, the total gas mass and HVMCs mass are much higher than MC-C2 after 7 Myr, while in f200n1000v, they are comparable to MC-C1 at 8 Myr.At this stage, the total molecular mass istotally composed of the HVMCs mass, so all of the molecular components have propagate beyond 800 pc.In f100n1500v, the clouds have lower-velocity and larger volume than the observed HVMCs, but the mean density is similar.Therefore, if the clouds crushed as some higher-velocity small clouds similar to the observed HVMCs, this run can better match the observation.It happens that the clouds will diffuse to be some small clumps with a mass-weighted mean velocity of ∼ 340 km s −1 in f100n1500h, though the velocity becomes a little higher than the observations.Therefore, it is natural to expect a magnetic field including both vertical and horizontal components, will help to produce the better-matched HVMCs in the simulation.Such a configuration is actually more reasonable for the real magnetic field in the Galactic center, A general model for the whole Milky Way also shows the magnetic field is parallel to the Galactic plane at lower latitude and gradually tend to be perpendicular at higher latitude (Cerri et al. 2017), so the expectation is sensible.In addition, if a higher resolution (4 times) is applied in the simulation, the clouds will also crush to be smaller clumps (Schneider & Robertson 2017; Gronke & Oh 2020), of which velocity and total mass are similar to those in the lower resolution.Therefore, it will be more consistent with the observations, since MW-C1 and MW-C2 are both smaller than the clouds produced in f100n1500v, f100n1000v and f100n1500h.In other words, the resolution used in our work is adequate to explain the formation of HVMCs, if we do not take the volume of the HVMCs as an essential feature.Of course, using a low resolution, the simulations cannot accurately describe the instability and the mixing between the cold gas and the hot wind, which may stimulate the crushing of clouds, but the advection of hot high-enthalpy gas into the mixing layer actually can result in growth and acceleration of the cold phase (Fielding et al. 2020).
The observations show many HVCs distributed over a large latitude from ∼ 100 pc to ∼ 10 kpc (Di Teodoro et al. 2018;Lockman et al. 2020;Lehner et al. 2022), though most of the HVCs are located in the lower 2 kpc.In our simulation, we only consider the starburst happening in a small region and include only one initial cloud, which limits the number of HVCs formed in the simulation box.However, the main focus of our work is to investigate the formation mechanism of HVMCs, rather than reproducing the exact number and distribution of observed HVCs.The fact that we can reproduce the key features of HVMCs observed in the Milky Way, such as their high velocity and high density, suggests that our proposed formation mechanism is plausible and can contribute to the understanding of the origin of HVCs in general.Further studies including more initial clouds and considering the starburst happening over a larger region would be needed to fully reproduce the observed distribution of HVCs.
The MCs in the run without magnetic field will be crushed in a short term, so the magnetic field is essential for the formation of the HVMCs.The magnetic field can wrap and protect the MCs, a mechanism named as the magnetic draping, which is significant at a large scale range, from the small scale of comets to the large scale of galaxy clusters (Riedler et al. 1986;Jun & Norman 1996;Brain et al. 2006;Dursi & Pfrommer 2008).Therefore, it is possibly contributed to the survival of our HVMCs.Nevertheless, Figure 2 shows the magnetic field surrounding the cold clouds is chaotic and does not well wrap the clouds, and our zoom-in check also shows same results, which is possibly caused by the low resolution and the wrapping is only obvious at much smaller scales.Jung et al. (2023) try to study the survival of HVCs in the Galactic halo, and claim that magnetic fields suppress hydrodynamic instabilities and the growth of smallscale structures, which is also responsible for the protection of the HVMCs in our simulations.In addition, the direction of the magnetic field can also influence the evolution of the HVMCs, which can be read from Figure 1, 5 and 7.In a vertical magnetic field, the clouds can keep high density and propagate to high latitude.Sparre et al. (2020) also conclude that the vertical magnetic field can well protect a cold cloud, but the cloud they used actually has a temperature of 10 4 K and a density of 0.1 cm −3 , totally different from the parameters used in our simulations.In a horizontal magnetic field, the clouds will lose an amount of mass, but still can propagate to high latitude without crushing.If there is not magnetic field, the clouds cannot propagate to high latitude.The importance of direction is also discussed by Alūzas et al. (2014) & Cottle et al. (2020), though the properties of clouds, winds, magnetic field and ISM they used are different from ours.
Additionally, our simulations consistently demonstrate that the reverse shock generated by the interaction between the clouds and the Galactic wind effectively balances the forward shock at later stages.As the clouds propagate, the forward shock of the Galactic wind encounters resistance from the clouds, leading to the gradual formation of stronger reverse shocks.This phenomenon is clearly observed in the density-velocity distribution plots presented in Figure ?? and 11.It is expected that at this late stage, the clouds have attained their maximum velocity within the framework of our model, and further acceleration becomes inefficient.Furthermore, we note that the star formation rate (SFR) employed in our model represents an upper limit within reasonable estimations, ensuring that the supernova explosion frequency is also maximized.Among the runs, f100n1500h stands out with the highest velocity exceeding 400 km s −1 , although it should be noted that the assumption of a complete horizontal magnetic field in the Galactic center is not physically realistic.Thus, if our model accurately captures the physics, we predict that the maximum velocity attainable by the HVMCs would be approximately 400 km s −1 .
Overall, the simulation results provide a promising framework for explaining the formation of HVMCs and their connection to HVCs.The HVMCs can indeed originate from a starburst in the Galactic center, which is reasonable in a large parameter space.The magnetic field can protect the MCs and contribute to the acceleration of MCs, but the acceleration of MCs is limited at high latitdue.However, there are still many uncertainties and complexities involved in the process, such as the role of magnetic fields, the effects of different initial conditions, and the possible interactions with other structures in the Galactic center.Therefore, further investigations are needed to refine and extend the current model, and to test its validity against more detailed observations and simulations.

The metallicity of HVCs
The formation of HVMCs is tightly associated with the HVCs', but the origin of HVCs is also ambiguous.The HVCs are usually defined as the interstellar gas clouds that moving at speeds substantially different (up to several hundreds km s −1 ) to the rotation of the disk of the Milky Way, and they are mostly distributed in the whole Galactic halo.Most of them have lower metallicity than what we find in the disk, so they may come from the Galactic halo or intergalactic medium.However, some of them, especially in the Fermi bubbles, have much higher metallicity, so they may be ejected from the Galactic disk.The HVCs in the Fermi bubbles are usually called as FB HVCs, which will be primarily discussed in this section.
It has been suggested that the HVCs are composed of diffuse inflowing gas and collimated outflowing material, which are likely manifestations of a galaxy-wide gas cycle triggered by stellar feedback, known as the galactic fountain (Li & Tonnesen 2020;Marasco et al. 2022).The feedback and the interaction with surrounding galaxies both influence the material cycle in our Milky Way, in which, most of the FB HVCs should be taken as a part of the collimated outflow  (Fox et al. 2015;Bordoloi et al. 2017;Ashley et al. 2020), because the stellar activity in the Galactic center is stronger than the disk.However, Ashley et al. (2022) found the FB HVCs have a wide range of metallicities from ≤ 0.2 of solar to ∼ 3.2  , thus the gas from the halo may also mix with the local ISM and ejecta from the disk.The supersolar metallicity of ∼ 3.2  implies that the HVCs are initially metal-rich, or there is a metal-enrichment process during the acceleration of the HVCs, since the Galactic ISM metallicity is usually ∼ 1 solar (Zuo et al. 2021).Therefore, it is convenient to assume the FB HVCs with high metallicity are formed by the driven of many sequential supernovae explosions which can simultaneously accelerate the clouds and provide heavy elements, a process also possibly happening in other galaxies (Emerick et al. 2019).The SMBH   2.
For the simulation results, points depict the molecular mass and mass-weighted velocity.The total mass and the MCs mass on strict criteria (see text for details) serves as the upper limit and lower limit, while the velocity on another strict criteria represents the lower limit of the velocity.
For MW-C1 and MW-C2, the central points display the molecular mass and outflow velocity, while the total mass serves as the upper limit and the velocity along the line of sight represents the lower limit.
activity may also drive the HVCs, but a metal-enrichment process, i.e., the supernovae explosions, is always necessary.The origin of HVMCs is likely analogous to FB HVCs, but this has yet to be confirmed due to the lack of information about their metallicity.To investigate this further, we examined the ratio of ejecta mass to cloud mass in our simulation, as shown in Figure 15.The ratio generally increases over time for all runs, but there is a peak at 5 Myr for f2001000v, which may be due to the low-metallicity cloud material flowing out of the simulation box.In f100n1500n, the clouds contain more ejecta material since they are slow, resulting in a more efficient mixture.Assuming an initial cloud metallicity of 1  and a supernova ejecta metallicity of 6  , a standard ratio of 0.1 would yield a final cloud metallicity of 1.5  , still lower than the observed 3.2  in some FB HVCs (Ashley et al. 2022).This suggests that the initial clouds were possibly already metal-rich before being driven to become HVCs.While SMBH activity may also drive HVCs, a metal-enrichment process such as supernova explosions is possibly necessary to explain the high metallicity of some FB HVCs.
If the model is correct, the role of HVCs in the galaxy-wide gas cycle can be understood.The low-metallicity HVCs originating from the halo or intergalactic medium are pulled by the gravitational potential of the Milky Way and surrounding galaxies, while highmetallicity HVCs are driven by galactic fountains that are energized by supernovae explosions in our Milky Way or the SMBH in the Galactic center.The FB HVCs consist of both types of HVCs, but the HVMCs embedded in FB HVCs should be driven by the fountains, which could be further confirmed by future metallicity analysis based on new ultraviolet absorption observations.

The relation between HVMCs and feedback relics
It is interesting to ask whether the HVMCs have a causal relation with the radio bubbles (Heywood et al. 2019) and X-ray chimneys (Ponti et al. 2019) found on smaller scales, or the Fermi bubbles (Su et al. 2010) and eROSITA bubbles (Predehl et al. 2020) found on much larger scales.We note that the age of the HVMCs inferred from our simulations is a few Myr, roughly consistent with the dynamical timescale of a few Myr for both of the radio bubbles and the Fermi bubbles originally suggested by Heywood et al. (2019) and Yang et al. (2013), respectively.However, their timescales actually have not been resolved, the radio bubbles may be younger (Zhang et al. 2021, 330 kyr) and the Fermi bubbles may be much older (Crocker & Aharonian 2011, 1 Gyr).In particular, Heywood et al. ( 2019)'s estimation was based on the assumption of a constant expansion velocity of the bubbles, which is implausible, hence a shorter timescale is expected.In the context of the supernova-based model for the origin of the radio bubbles/chimneys (Zhang et al. 2021), the radio bubbles would be a dynamically younger and independent structure simply evolving in the interior of the Fermi/eROSITA bubbles, which themselves were formed by older activities in the Galactic center.However, the HVMCs should also originate from a similar activity, which implies there are three independent activities, respectively correlated with the radio bubbles/X-ray chimneys, the HVMCs and the Fermi bubbles/eROSITA bubbles.The difference is that the HVMCs will be difficult to propagate to much higher latitude in our simulations, because the acceleration rate of HVMCs at high latitude will largely decrease.If the three independent activities are not related with each other, we have to use three models to respectively explain the structures at three scales, which will lead to an inelegant physical pattern.
Alternatively, as suggested by Ponti et al. (2019), the X-ray chimney/the radio bubbles may be a channel that transports energy from the Galactic center to the high-latitude region currently occupied by the Fermi bubbles, and the HVMCs are the manifestation of the transportation process, which is a more elegant unified model.In fact, the HVCs can spread from ∼ 100 pc to ∼ 10 kpc (Di Teodoro et al. 2018;Lockman et al. 2020;Lehner et al. 2022), though most of the FB HVCs are located in the lower 2 kpc, which may be the clue connecting the feedback relics at different scale.In this case, the channel should have existed for tens of Myr, so that star formation in the Galactic center can be sufficient to supply the total energy content of the Fermi bubbles, ∼ 10 56 erg (Carretti et al. 2013).However, such a picture contradicts with the capped morphology of radio bubbles (the southern bubble is not obviously capped in X-rays; Ponti et al. 2021), which, according to our simulations, is naturally explained as the expanding shell of a newly born outflow.This picture may be reconciled if star formation in the Galactic center has been episodic on a timescale of ∼10 Myrs (Krumholz & Kruĳssen 2015), then the X-ray chimney/the radio bubbles are (re)established and the HVCs/HVMCs are (re)accelerated by consecutive generations of mini-starbursts and collapses inbetween.Of course, over such a long interval, the activity of Sgr A* can also play an important role in contributing to the formation of these relics, especially in view of the fact it was likely much more active in the recent past (Ponti et al. 2010(Ponti et al. , 2013;;Camilo et al. 2018).In a hybrid scenario, Sgr A*, with supernovae and even stellar winds, can simultaneously sustain the channel and transport energy to larger scales, implying X-ray emission beyond the edge of the radio bubbles, which is also suggested by Ponti et al. (2021).For example, a AGN activity produces the large-scale structure and triggers the surrounding starburst, then the newly-formed massive stars drive strong stellar wind and explode as supernovae to produce the small-scale structure.Possibly, the stellar winds and shock wave of supernovae can also trigger the tidal disruption event of the central SMBH, then produce a smaller-scale structure.
In conclusion, our findings suggest the existence of a potentially stable channel in the Galactic center, driven by a combination of diverse activities, which episodically accelerates gas clouds and transports energy to higher latitudes.The HVMCs/FB HVCs are also the ingredient of the channel, but the HVMCs usually exist in low latitude due to the higher possibility of crushing at higher latitude.This pattern offers a comprehensive explanation for the interrelation between various feedback remnants, without necessitating the introduction of new models.

SUMMARY
To investigate the formation of HVMCs in our Galactic center, we perform simulations utilizing a starburst model, where HVMCs originate from low-latitude molecular clouds accelerated by a subsequent supernovae explosions.Previous studies have raised concerns about the destruction of molecular clouds due to the violent activity in the Galactic center, making it challenging for them to reach higher latitudes and velocities without disruption.However, our simulation results demonstrate that this problem can be resolved within a wide parameter space, given the appropriate local environment.
The main findings are summarized as follows: • The HVMCs can indeed be formed in a starburst in the Galactic center.
• The magnetic field can protect the molecular clouds.
• The magnetic pressure, enhanced by the compression of shock wave, can contribute to accelerating the clouds.
• The acceleration rate of HVMCs will largely decrease at high latitude, because the reverse shock, generated by the interaction between the shock wave and the molecular clouds, can gradually balance the forward shock from the supernovae.Therefore, we can predict the largest velocity the HVMCs can reach is ∼ 400 km s −1 .
• The mixture between the clouds and the ejecta of the supernovae is more efficient at low latitude, and this process can significantly impact the metallicity of HVCs.
• HVMCs/FB HVCs potentially serve as ingredients in a channel sustained by diverse activities in the Galactic center, intermittently accelerating gas clouds and transporting energy to higher latitudes.
Due to the limited size of the simulation box, the subsequent evolution of HVMCs beyond 1 kpc latitude remains uncertain.Furthermore, the small box size restricts us to initializing only one cloud, resulting in inconsistent HVC number density and distribution compared to observations.Future efforts involve expanding the simulation box, simplifying supernova explosion settings, and implementing adaptive mesh refinement to provide a more comprehensive understanding of the phenomenon.
and 4. At this moment, the simulation box contains three components: the clouds, ISM-dominated and SNR-dominated region, respectively corresponding to the lower right, lower left and central part of Figure 3, and the lower right part, the central and the upper left band of Figure 4.
Figure 4 is similar to the Figure 8 of Schneider & Robertson (

Figure 1 .
Figure 1.The - column density maps of f100n1500v between 1∼7 Myr with a step of 1 Myr.The white arrows show the flow velocity in the slice through the  = 0 pc, and the scale is shown at the upper right.The main cloud can indeed survive with comparable mass with the observation until it reaches 1 kpc at 7 Myr, though it will lose a large amount mass.

Figure 2 .Figure 3 .Figure 4 .
Figure 2. The temperature and the magnetic field maps of f100n1500v in the slice through the  = 0 pc between 1∼7 Myr with a step of 1 Myr.The red arrows show the magnetic field, and the scale is shown at top right.The cold gas slowly diffuse away from the central slice, and almost dissipates at 7 Myr.The magnetic strength is amplified at the early stage, but gradually decreases after ∼ 3 Myr.

Figure 5 .Figure 6 .
Figure5.The column density maps of f100n1500h between 1∼6 Myr with a step of 1 Myr.The white arrows show the flow velocity in the slice through the  = 0 pc, and the scale is shown at the upper right.The main cloud can also survive until it reaches 1 kpc , but it will lost much more mass than f100n1500v.MNRAS 000, 1-?? (0000)

Figure 7 .
Figure 7.The - column density maps of f100n1500n between 1∼6 Myr with a step of 1 Myr.The white arrows show the flow velocity in the slice through the  = 0 pc, and the scale is shown at the upper right.The main cloud crushed quickly and cannot survive after 6 Myr.MNRAS 000, 1-?? (0000)

( 1 )Figure 8 .
Figure8.The density-velocity map for f100n1500n at 5 Myr.The little pink box shows the observed mean density and velocity range, and the larger pink box shows its zoom-in picture.The velocity is binned for every 100 km s −1 , so it can be conveniently read by counting the bins.The map shows the mass of every bin in solar mass, i.e., we can also estimate the mass of different components by counting the bins.The bins lower than 1 M are suppressed.

Figure 9 .Figure 10 .
Figure 9.The - column density maps of f100n1000v between 1∼6 Myr with a step of 1 Myr.The white arrows show the flow velocity in the slice through the  = 0 pc, and the scale is shown at the upper right.The main cloud can indeed survive until it reaches 1 kpc at 5 Myr, though it will lost a large amount mass.MNRAS 000, 1-?? (0000)

Figure 11 .
Figure 11.The - column density maps of f200n1000v between 1∼8 Myr with a step of 1 Myr.The white arrows show the flow velocity in the slice through the  = 0 pc, and the scale is shown at the upper right.The main cloud can indeed survive until it reaches 1 kpc at 8 Myr, though it will lost a large amount mass.

Figure 12 .
Figure12.The density-velocity map for f200n1000v at 8 Myr.The little pink box shows the observed mean density and velocity range, and the larger pink box shows its zoom-in picture.The velocity is binned for every 100 km s −1 , so it can be conveniently read by counting the bins.The map shows the mass of every bin in solar mass, i.e., we can also estimate the mass of different components by counting the bins.The bins lower than 1 M are suppressed.

Figure 13 .
Figure 13.Visualization of the Mass-Velocity Relation presented in Table2.For the simulation results, points depict the molecular mass and mass-weighted velocity.The total mass and the MCs mass on strict criteria (see text for details) serves as the upper limit and lower limit, while the velocity on another strict criteria represents the lower limit of the velocity.For MW-C1 and MW-C2, the central points display the molecular mass and outflow velocity, while the total mass serves as the upper limit and the velocity along the line of sight represents the lower limit.

Figure 14 .Figure 15 .
Figure14.The mass evolution of HVMCs.The solid, dashed and dotted lines respectively show the corresponding the total mass, the total molecular mass and the molecular mass with a latitude higher than 800 pc.

Table 1 .
Parameters of the simulation runs

Table 2 .
The resultant parameters at the fiducial times