Arkenstone is a new scheme that allows multiphase stellar feedbackdriven winds to be included in coarse resolution cosmological simulations. The evolution of galactic winds and their subsequent impact on the circumgalactic medium are altered by exchanges of mass energy momentum and metals between their component phases. These exchanges are governed by complex smallscale physical processes that cannot be resolved in cosmological simulations. In this second presentation paper we describe Arkenstones novel cloud particle approach for modelling unresolvable cool clouds entrained in hot fast winds. This general framework allows models of the cloudwind interaction derived from stateoftheart highresolution simulations to be applied in a largescale context. In this work we adopt a cloud evolution model that captures simultaneous cloud mass loss to and gain from the ambient hot phase via turbulent mixing and radiative cooling respectively. We demonstrate the scheme using noncosmological idealized simulations of a galaxy with a realistic circumgalactic medium component using the Arepo code. We show that the ability of a highspecific energy wind component to perform preventative feedback may be limited by heavy loading of cool clouds coupled into it. We demonstrate that the diverging evolution of clouds of initially differing masses leads to a complex velocity field for the cool phase and a cloud mass function that varies both spatially and temporally in a nontrivial manner. These latter two phenomena can manifest in the simulation because of our choice of a Lagrangian discretisation of the cloud population in contrast to other proposed schemes.
['Astrophysics - Astrophysics of Galaxies']
Arkenstone  II. A model for unresolved cool clouds entrained in galactic winds in cosmological simulations
{'A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc - II. A model for unresolved cool clouds entrained in galactic winds in cosmological simulations': 'Matthew C. Smith , 1 ★ Drummond B. Fielding , 2 , 3 Greg L. Bryan , 4 Jake S. Bennett , 5 Chang-Goo Kim , 6 Eve C. Ostriker 6 and Rachel S. Somerville 2 1 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748, Garching, Germany 2 Center for Computational Astrophysics, Flatiron Institute, 162 5 th Avenue, New York, NY 10010, USA 3 Department of Astronomy, Cornell University, Ithaca, NY 14853, USA 4 Department of Astronomy, Columbia University, 550 West 120 th Street, New York, NY 10027, USA 5 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  \n6 Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': "A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is a new scheme that allows multiphase, stellar feedback-driven winds to be included in coarse resolution cosmological simulations. The evolution of galactic winds and their subsequent impact on the circumgalactic medium are altered by exchanges of mass, energy, momentum, and metals between their component phases. These exchanges are governed by complex, small-scale physical processes that cannot be resolved in cosmological simulations. In this second presentation paper, we describe A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc's novel cloud particle approach for modelling unresolvable cool clouds entrained in hot, fast winds. This general framework allows models of the cloud-wind interaction, derived from state-of-the-art high-resolution simulations, to be applied in a large-scale context. In this work, we adopt a cloud evolution model that captures simultaneous cloud mass loss to and gain from the ambient hot phase via turbulent mixing and radiative cooling, respectively. We demonstrate the scheme using non-cosmological idealized simulations of a galaxy with a realistic circumgalactic medium component, using the A/r.pc/e.pc/p.pc/o.pc code. We show that the ability of a high-specific energy wind component to perform preventative feedback may be limited by heavy loading of cool clouds coupled into it. We demonstrate that the diverging evolution of clouds of initially differing masses leads to a complex velocity field for the cool phase and a cloud mass function that varies both spatially and temporally in a non-trivial manner. These latter two phenomena can manifest in the simulation because of our choice of a Lagrangian discretisation of the cloud population, in contrast to other proposed schemes. This is a Learning the Universe publication.  \nKey words: galaxies: evolution - methods: numerical - hydrodynamics", '1 INTRODUCTION': "Galactic winds are observed to be a fundamental component of the cosmic baryon cycle (Tumlinson et al. 2017) and therefore play an essential role in galaxy formation theory (Somerville & Davé 2015; Naab & Ostriker 2017). As gas inflows transport mass from the intergalactic medium (IGM), through the circumgalactic medium (CGM) to the interstellar medium (ISM) of galaxies (Kereš et al. 2005; Dekel & Birnboim 2006), large-scale winds driven by feedback from stars or active galactic nuclei (AGN) provide a counterbalance, transporting mass and energy outwards. They can therefore regulate galactic star formation as an ejective feedback process, removing gas from the ISM before it can be converted into stars (e.g. Mathews & Baker 1971; Larson 1974; Chevalier & Clegg 1985; Dekel & Silk 1986; White & Frenk 1991), or preventative, reducing inflows from the CGM/IGM (e.g. Oppenheimer et al. 2010; van de Voort et al. 2011; Davé et al. 2012; Lu et al. 2015, 2017; Pandya et al.  \n2020; Carr et al. 2023). They also shape the properties of the CGM, influencing its baryon fraction, phase structure, and metallicity.  \nGalactic winds are often parameterised in terms of loading factors. The mass loading factor of a stellar feedback driven wind is the ratio of the emergent mass outflow rate to the star formation rate (SFR) of the galaxy driving it. Likewise, the energy loading is the emergent power of the wind relative to the power generated by stellar feedback (i.e. an efficiency factor, generally relative to the supernova input). The specific energy of the wind (i.e. its temperature and velocity) depends on the ratio of the energy and mass loading factors. The impact of a wind on its host CGM and the manner in which it regulates star formation depends not only on the mass outflow rate but on its specific energy (Carr et al. 2023; Voit et al. 2024a,b). Two winds can have the same energy loading but different combinations of mass loading and specific energy. A highly mass loaded wind is a source of ejective feedback, removing mass from the ISM. However, for a fixed energy loading, increasing the mass loading reduces the specific energy of the wind; the available energy is diluted by being coupled to increasing amounts of mass. Such a wind has limited  \nability to impact the CGM, tending to form fountain flows that return to the galaxy. On the other hand, a high specific energy wind (with a low mass loading factor) can transmit energy to large volumes of the CGM and beyond, heating it and suppressing inflows (i.e. preventative feedback), but has little impact on the gas which is already in the ISM.  \nThis picture is complicated by the fact that characterising a wind with a single mass outflow rate and specific energy is a substantial oversimplification. In fact, they are observed to be highly multiphase with detections of components with temperatures of ≲ 100K (e.g., Rupke et al. 2005; Bolatto et al. 2013; Martini et al. 2018), ∼ 10 4 K (e.g., Martin & Bouché 2009; Westmoquette et al. 2009; Nielsen et al. 2015), ∼ 10 5 . 5 K(e.g., Steidel et al. 2010; Kacprzak et al. 2015; Chisholm et al. 2018), and ≳ 10 7 K (e.g., Strickland & Heckman 2009; Lopez et al. 2020; Hodges-Kluck et al. 2020). High resolution simulations that can resolve the generation of multiphase outflows by stellar feedback show that different components of the wind carry different fractions of the total mass and energy loadings (e.g. Kim & Ostriker 2018; Schneider et al. 2020; Rathjen et al. 2023; Steinwandel et al. 2024). For example, Fielding et al. (2018) and Kim et al. (2020a) find that as the wind leaves the ISM a hot ( ≳ 10 6 K), fast component carries the majority of the energy while a cool ( ≲ 10 4 K), slow component carries most of the mass. With idealised simulations, Tan &Fielding (2024) demonstrate the seeding of individual cool clouds into the hot wind as the swept-up ISM fragments while the wind is generated by clustered supernovae (SNe).  \nA further complication is that as the wind flows away from the galaxy, the various phases will interact, redistributing mass, momentum and energy between them and impacting the large-scale evolution of the wind (Fielding & Bryan 2022, hereafter FB22). As demonstrated in 'cloud-crushing' simulations (Klein et al. 1994; Xu & Stone 1995) a cool cloud with radius 𝑟 cl located within a hot, fast moving ambient wind is accelerated by drag/ram pressure on a time-scale 𝑡 drag = 𝜒𝑟 cl / 𝑣 rel but tends to be destroyed by Kelvin-Helmholtz and Rayleigh-Taylor instabilities on a cloudcrushing time-scale 𝑡 cc = 𝜒 1 / 2 𝑟 cl / 𝑣 rel , where 𝜒 and 𝑣 rel are the density contrast and relative velocity between the cloud and wind, respectively. For clouds in hot galactic winds, 𝜒 ≫ 1, so the timescale for cloud destruction is shorter than that for cloud acceleration (Zhang et al. 2017). However, the picture is modified by the inclusion of various additional physical processes such as radiative cooling (e.g. Mellema et al. 2003; Cooper et al. 2009; Scannapieco & Brüggen 2015), thermal conduction (e.g. Marcolini et al. 2005; Orlando et al. 2005; Brüggen & Scannapieco 2016), magnetic fields (e.g. Mac Low et al. 1994; McCourt et al. 2015; Schneider & Robertson 2017; Cottle et al. 2020) and cosmic rays (e.g. Wiener et al. 2019; Brüggen & Scannapieco 2020). Under particular circumstances, radiative cooling can cause the cloud to accrete mass from the hot phase (e.g. Marinacci et al. 2010; Armillotta et al. 2016; Gronke & Oh 2018). Along the cloud-wind boundary, shearing motions drive turbulence within a mixing layer between the two phases. The mixing of the hot wind and cool cloud material generates gas of intermediate temperature ( ∼ 10 5 K) which can efficiently cool into the cloud. Recent high resolution studies of these turbulent radiative mixing layers (TRMLs) have begun to constrain the complex physics which governs the rate at which hot wind material is accreted into the cloud (e.g. Ji et al. 2019; Mandelker et al. 2020; Fielding et al. 2020; Tan et al. 2021; Abruzzo et al. 2024). Importantly, this material brings momentum along with its mass, providing a potentially significant form of acceleration for the clouds. This was demonstrated by Melso et al. (2019) in cloud inflow simulations, by Vĳayan et al. (2020) from analysis of outflows in simulations of star-forming disks, and  \nin Schneider et al. (2020) on kpc scales with 5 pc uniform spatial resolution idealised simulations of a starbursting galaxy.  \nUnfortunately, the resolution requirements for capturing the generation and subsequent evolution of a multiphase wind in an a priori manner in a simulation are very stringent. When injecting energy from individual SNe directly into the ISM, a mass resolution ≲ 100 M ⊙ is required just to get convergent bulk mass and energy loadings (Smith et al. 2018; Hu 2019), let alone obtain the correct multiphase structure. In order to approximately capture the mass balance between the phases, involving the interaction between cool clouds and hot winds, a mass resolution ≲ 1M ⊙ is required (e.g. Gronke & Oh 2020; Abruzzo et al. 2022; Gronke et al. 2022). This resolution is barely achievable in cosmological 'zoom-in' simulations of the lowest mass dwarf galaxies. Meanwhile, cosmological volume simulations, necessary to build statistically significant samples of galaxy evolution, typically have mass resolutions of ∼ 10 5 -10 9 M ⊙ . The highest resolution simulations in this class do see multiphase structure in galactic outflows and the CGM in general (e.g. Nelson et al. 2019; Mitchell et al. 2020), but these cannot reliably resolve the interactions between the phases. The situation improves significantly if computing resources are concentrated on resolving the CGM of a single galaxy, producing a rich multiphase structure of cool clouds suspended in a volume filling medium (van de Voort et al. 2019; Peeples et al. 2019; Suresh et al. 2019; Hummels et al. 2019; Bennett & Sĳacki 2020; Ramesh & Nelson 2024), but the total mass of cool material in the CGM typically continues to increase as the spatial resolution is improved. These simulations do begin to allow individual, relatively massive cloud-like structures to be identified and processes affecting their evolution to be studied (e.g. Ramesh et al. 2024), but the spatial resolution is still several orders of magnitude coarser than what is required to properly capture the complex mixing processes that should occur on their boundaries. Regardless, the CGM zoom-in technique is computationally intractable to apply to cosmological volumes.  \nEven assuming the most optimistic estimates of increases in computing power and code efficiency, we will not gain the ability to fully resolve multiphase galactic wind and CGM material in cosmological volume simulations for many years. An alternative approach is to avoid attempting to resolve the multiphase material natively and instead develop effective models that represent the large-scale impact of the unresolvable multiphase gas in a 'subgrid' manner. This is the approach we take in this work. While the small-scale physics must necessarily be modelled in a simplified manner, this permits them to be included in a clean, interpretable fashion. For many science questions, this is preferable to poorly resolving or entirely omitting the relevant processes. Huang et al. (2020) presents a subgrid scheme for cool clouds in galactic winds with a particle based approach that models their disintegration, modulated by the effects of thermal conduction. Weinberger & Hernquist (2023) and Butsky et al. (2024) present Eulerian multifluid approaches to treat two phases of gas within the same resolution element. We will discuss these three schemes in detail in Section 4.2.  \nA/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is a new subgrid model for stellar feedback-driven galactic winds that emphasises their multiphase nature. The scheme is intended for large volume cosmological simulations where neither the multiphase structure of the ISM, galactic winds or CGM can be resolved. A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc has been specifically designed to provide a flexible framework for coarse-graining results from analytical models and very high resolution numerical studies of relevant small-scale physics. A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is developed in concert with the Learning  \nthe Universe (LtU) Collaboration. 1 LtU will deploy A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc in the next generation of cosmological volume simulations, alongside improved models for the ISM and star formation. While distinct from A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc, these new ISM models (Hassan et al. in press will present a prototype) will be calibrated to ISM patch simulations (Kim et al. 2024) that can also be used to inform the input parameters for A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc (see e.g. Kim et al. 2020b).  \nA/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc uses a wind particle propagation scheme to inject winds immediately outside of the ISM, granting very fine control over the properties of the wind at launch irrespective of the coarse resolution that must be adopted in this type of simulation. The model has three novel features:  \n- (i) Winds are launched with hot and cool components with separate mass and energy loadings, inspired by the results from high resolution simulations, as mentioned above.\n- (ii) The hot, fast phase of the wind is injected and evolved with a new 'displacement recoupling' and refinement scheme that properly treats high-specific energy, low density flows.\n- (iii) The cool phase is modelled using 'cloud particles' to represent clouds embedded in the hot flow. These particles exchange mass, energy, momentum, and metals bidirectionally with the ambient hot wind.  \nWe first presented A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc in Smith et al. (2024) (hereafter Paper I), detailing the first two aspects of the scheme. In particular, we demonstrated that without the techniques employed by our scheme it is impossible to resolve the high-specific energy winds expected to drive preventative feedback processes. This is because the low densities inherent to high-specific energy winds cause poor spatial resolution when quasi-Lagrangian refinement strategies (ubiquitous in cosmological simulations carried out with both Lagrangian and Eulerian codes) are used. In particular, failing to resolve the sonic point of a wind results in an incorrect evolution of the balance of kinetic and thermal energy in the wind as it flows outwards. A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc makes it possible to properly consider preventative feedback that originates from SNe in cosmological simulations. We refer the reader to Paper I for full details of the scheme, as well as a detailed discussion of other stellar feedback and galactic wind schemes in the literature.  \nIn this work, we present the remaining aspect of A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc: the cloud particle scheme. This is a framework with which results drawn from state-of-the-art cloud-crushing and TRML simulations, that capture the fine details of the relevant small-scale physics, can be included in large-scale cosmological simulations. The implementation is agnostic to the choice of cloud-wind interaction model, so long as fluxes of mass, momentum, energy and metals between the phases can be predicted. For this first demonstration, we use the model of FB22. In Section 2 we present our numerical methodology. In Section 3 we apply our scheme to a series of idealised simulations of isolated galaxies (carried out at a resolution achievable in a cosmological volume simulation), demonstrating the general behaviour of the model and highlighting various interesting regimes. In Section 4 we discuss future applications of the model, the interpretation of its predictions and possible extensions to included physical processes. We also compare our scheme to other relevant approaches in the literature. We summarise our findings in Section 5.", '1.1 Nomenclature': "Weare concerned with multiphase winds comprised of various components, as described above. However, much of the cloud-crushing literature refers to cloud-wind interactions. Of course, once entrained, one can consider cool material as being part of the multiphase wind. For clarity, we adopt the following conventions. 'Hot wind' is used as a short hand for high-specific energy galactic winds (which are hot, fast and low density). Aspects of the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc scheme relating to the modelling of hot winds were described in detail in Paper I and are referred to as A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc-H/o.pc/t.pc. We generate the hot wind component with the use of 'wind particles' (a summary of the methodology can be found in Section 2.3). Cool clouds are modelled with 'cloud particles'. Here, we use 'cool' to refer to any gas with temperature ≤ 10 4 K; while 'cold' is often used for this thermal range in the field of galaxy formation, in the ISM literature the term 'cold' is reserved for atomic and molecular gas at ≲ 10 2 . 5 K. The term 'cloud-wind interaction' and similar are used as a shorthand for interactions between the clouds and hot wind components in the wind. We stress, however, that the cloud particles interact with all ambient gas that they encounter, not just wind material. The label 'hydro' is applied to quantities and measurements related to the resolved gas treated by A/r.pc/e.pc/p.pc/o.pc's finite volume scheme. This is often contrasted in figure legends with 'Clouds' which refers to gas modelled with cloud particles. We adopt lowercase 𝑟 for a spherical radius, using 𝑅 for a cylindrical radius in the plane of a galactic disc.", '2 NUMERICAL METHODS': 'A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is implemented in the A/r.pc/e.pc/p.pc/o.pc code (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020). In Section 2.1 we describe relevant features of A/r.pc/e.pc/p.pc/o.pc and other details of the code setup used in this work that are not specific to A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc. In Section 2.2 we describe the new cloud particle scheme. In Section 2.3 wedescribe the launching and recoupling of wind and cloud particles, briefly summarising relevant details of A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc-H/o.pc/t.pc that were presented in Paper I.', '2.1 Hydrodynamics, gravity, cooling and the ISM': 'A/r.pc/e.pc/p.pc/o.pc uses a finite volume scheme, solving hydrodynamics on an unstructured, moving mesh. The mesh is defined by the Voronoi tessellation of mesh-generating points which move with the local fluid velocity with small corrections to maintain cell regularity. This means that cells tend to maintain constant mass over time, giving the scheme quasi-Lagrangian properties. However, while minimised by the mesh motion, mass fluxes between cells are non-zero, so a (de)refinement scheme is typically used to (merge) split cells to keep them within a factor of two of a desired mass resolution. In addition to the constant mass (de)refinement scheme, other criteria can be used to enforce varying mass or spatial resolution within the simulation domain. A/r.pc/e.pc/p.pc/o.pc can solve magnetohydrodynamics (MHD) (Pakmor et al. 2011) but we do not include magnetic fields in this work. None the less, we remark that A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is fully compatible with the MHDscheme. Gravity is included with a tree-based algorithm. 2  \nWe include radiative cooling as described in Vogelsberger et al. (2013). This includes cooling from both primordial species (Cen 1992; Katz et al. 1996) and metal lines (in pre-calculated lookup  \ntables) in the presence of a 𝑧 = 0 UV background (Faucher-Giguère et al. 2009), with corrections for self-shielding in dense gas (Rahmati et al. 2013). While we do not impose a formal temperature floor, we donot radiatively cool below 10 4 K. At the typical resolution at which A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is intended to operate, the multiphase ISM cannot be well resolved. Therefore, in this work, we use the model of Springel &Hernquist (2003) (hereafter SH03), adopting an effective equation of state (eEoS) to represent the large-scale impacts of small-scale ISM physics (e.g. local stellar feedback, turbulence, molecular cloud formation and destruction etc.) in an abstract manner. The eEoS and star formation are switched on for gas denser than a threshold value of 𝜌 SF / 𝑚 p = 0 . 2 cm -3 . For the fine details of our parameter choices for this model (which are essentially the same as the TNG suite), see Paper I. Relevant for this work, however, is that this ISM model predicts a star formation rate (SFR), / 𝑚 ★ , for each cell denser than 𝜌 SF . This rate is stochastically sampled to create star and wind particles (as described in Paper I), as well as cloud particles (as described in Section 2.3). We emphasise that A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is agnostic as to the choice of subgrid ISM model or the method used to assign an SFR to each cell, except that it assumes that internal ISM structure is unresolved. In future, we will explore alternative models to SH03 (such as that proposed by Ostriker & Kim 2022; Kim et al. 2024, Hassan et al. in press, as calibrated from resolved star-forming, multiphase ISM simulations).', '2.2 Cloud particle evolution': 'We now describe the implementation of the cloud particle scheme. Wewill first outline how A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc implements cloud particle and background gas interactions in general terms - i.e. the parts of the scheme that are independent of any particular theory or assumptions of cloud-wind interactions. We then give specific details on how we have implemented the FB22 cloud evolution model.', '2.2.1 General model': "A cloud particle is interpreted as carrying some number of identical clouds, 𝑁 cl . The constraint that the clouds must be identical (meaning that they have the same properties as all other clouds carried by the particle in which they are hosted at all times) arises because the trajectories of clouds with different properties will diverge. The cloud number relates the masses of individual subgrid clouds, 𝑚 cl , to the mass of the cloud particle, 𝑚 part :  \n𝑁 cl = 𝑚 part 𝑚 cl . (1)  \n𝑁 cl is assigned to a cloud particle at its creation and remains constant throughout its evolution. 3 Note that we do not constrain 𝑁 cl to be an integer nor to be greater than unity. Thus, an ensemble of cloud particles represents the population of individual clouds in a statistical sense.  \nWith the constraint that all clouds hosted by a given cloud particle are identical to each other, we can use 𝑁 cl to relate all changes of  \n3 While not currently implemented, any form of particle splitting or merging would naturally also involve changing 𝑁 cl in order to conserve cloud number. Likewise, while not examined in this work, cloud evolution models could be easily implemented that involve splitting or merging clouds within the particle, modifying 𝑁 cl . However, the constraint that all clouds within the particle are identical persists.  \nconserved quantities of clouds to those of the particle. For example, the rate of change of mass of the cloud particle is  \n/ 𝑚 = 𝑁 cl / 𝑚 cl , (2)  \nwhere / 𝑚 cl is the rate of change of mass of an individual cloud in the particle.  \nThe mass of a particle changes as a result of exchanges with the gas cell in which it is located. These exchanges are bidirectional in the sense that the particle can be accreting mass from the background gas at a growth rate / 𝑚 grow while it is simultaneously losing mass to the background gas at a loss rate / 𝑚 loss . The determination of / 𝑚 grow and / 𝑚 loss depends on the adopted cloud-wind interaction model. The net transfer rate from cell to particle is therefore  \n/ 𝑚 = / 𝑚 grow - / 𝑚 loss . (3)  \nTherefore, at any given moment, a particle may be experiencing a net inflow or outflow of mass from the background gas or no net mass transfer. However, it is important to note that / 𝑚 = 0 does not necessarily mean there are no net transfers of other quantities taking place. For example, the net rate of metal mass transfer from cell to particle is  \n/ 𝑚 𝑍 = 𝑍 cell / 𝑚 grow -𝑍 part / 𝑚 loss , (4)  \nwhere 𝑍 cell is the metallicity of the cell and 𝑍 part is the metallicity of the particle (which is equal to the metallicity of the subgrid clouds). This means that even if the net mass transfer is zero, if the particle and cell have different metallicities there will still be a net transfer of metals (unless both / 𝑚 grow and / 𝑚 loss are zero). In general, all transfers of passive scalars (e.g. individual metal species, tracer dyes etc.) are handled using an equivalent version of equation 4.  \nMomentum transfer rates from cell to particle are  \n/ p = / p trans + / p drag , (5)  \nwhere the first term captures momentum transfer associated with the mass transfers between the cell and the particle and the second term allows for the inclusion of an additional drag/ram pressure acceleration force (the specifics of which depend on the adopted cloud-wind interaction model). By conservation of momentum, we can derive  \n/ p trans = / 𝑚 grow v cell - / 𝑚 loss v part , (6)  \nwhere the cell and particle velocities are v cell and v part , respectively. The resulting acceleration of the cloud particle due to this component of the momentum transfer is  \n/ v trans = / p trans - / 𝑚 v part 𝑚 part = -GLYPH<0> v part -v cell GLYPH<1> / 𝑚 grow 𝑚 part . (7)  \nNote that this acceleration depends solely on / 𝑚 grow , not on the net mass transfer. Note also that the direction of this acceleration is always anti-parallel to the relative velocity between the particle and the cell, v rel = v part -v cell . Since this is also usually true for definitions of the drag force, / p drag , the momentum transfer between cell and particle always acts to reduce their relative velocity, as might be intuitively expected.  \nWe define the rate / 𝐸 th for the net transfer rate of thermal energy from the cell to the particle. Kinetic energy transfer is implicitly mediated via the momentum transfer described above.  \nThe source terms arising from the exchange of mass, momentum, energy and passive scalars (e.g. metals) are integrated in an explicit fashion by half timesteps immediately before and after the cloud  \nparticle receives its first and second gravity kicks, respectively. 4 We limit the timestep of the cloud particle, Δ 𝑡 , such that:  \nΔ 𝑡 = MIN GLYPH<0> 𝑡 grav , 𝑡 cell , 𝑡 cellmod , 𝑡 grow , 𝑡 loss , 𝑡 stop , 𝑡 cross GLYPH<1> , (8)  \nwhere each of these time-scales is defined below. The first timestep limit is the gravitational timestep limiter, 𝑡 grav , which is the standard limiter for all particles (e.g. dark matter, stars, black holes) and gas cells in the simulation. The default choice in A/r.pc/e.pc/p.pc/o.pc is  \n𝑡 grav = √︄ 2 𝐶 grav 𝜖 soft | a | , (9)  \nwhere 𝐶 grav = 0 . 012 (the typical choice), 𝜖 soft is the gravitational softening length and | a | is the magnitude of the gravitational acceleration.  \nThe timestep limit 𝑡 cell corresponds to the current timestep of the gas cell in which the cloud particle is located, which is itself limited by the gravitational timestep limit, the hydrodynamical timestep limit (the Courant criterion) and any other limiters that may be applied by additional physics.  \nWe ensure that we resolve the timescale on which the host cell's properties are modified by enforcing the timestep limit  \n𝑡 cellmod = 𝑓 cellmod 𝑚 cell | / 𝑚 | , (10)  \nwhere 𝑓 cellmod is a free parameter. For the case of net cloud mass growth, this limiter ensures that the cell consumption time is resolved (in particular that the particle will not attempt to accrete more mass than is available). In the case of net cloud mass loss, this assists in allowing the cell's mass to grow smoothly, rather than in a sudden injection of material. However, this latter case depends more strongly on the timestep with which the cell is being integrated, which we discuss later in this section.  \nWe also define timestep limits associated with the time-scales for mass growth and loss relative to the cloud particle mass as follows:  \n𝑡 grow = 𝑓 grow 𝑚 part / 𝑚 grow , (11)  \n𝑡 loss = 𝑓 loss 𝑚 part / 𝑚 loss , (12)  \nwhere 𝑓 grow and 𝑓 loss are free parameters. Note that we resolve the growth and loss time-scales independently, rather than purely considering the net growth/loss time-scale (which necessarily cannot be shorter than either the independent growth or loss time-scales). This allows us to resolve changes in properties of the cloud particle that are correlated with the fluxes in and out of the particle (e.g. metal transfer) even when the net mass transfer rate is small (or zero). It can be seen that 𝑡 loss can tend to zero when a cloud is being destroyed, which would obviously be intractable to resolve. However, this is avoided by simply fully recoupling a cloud particle once its mass drops below some threshold.  \nAs noted above, all cloud - wind interactions reduce the relative velocity between the cloud and the wind. We wish to resolve the stopping time of the cloud particle i.e. the time for the relative velocity between the cloud particle and the cell to reach zero. Failing to do so sufficiently results in spurious oscillations of the velocity of the  \ncloud particle as it repeatedly overshoots the velocity of the cell. We therefore define the timestep limit  \n𝑡 stop = 𝑓 stop 𝑣 rel / 𝑣 rel = 𝑓 stop 𝑚 part 𝑣 rel | / p | -𝑣 rel / 𝑚 , (13)  \nwhere 𝑓 stop is a free parameter and 𝑣 rel is the magnitude of the relative velocity. Note that this estimate of the stopping time only accounts for changes in the relative velocity arising from equation 5. However, the gravitational acceleration will be approximately the same for the particle and the cell due to their proximity, so this will be a minor contributor to changing 𝑣 rel . The cell may feel additional forces that the particle does not (e.g. hydrodynamics) but if these are dominant then this will already be captured via the 𝑡 cell timestep limiter. It is possible for 𝑡 stop to tend to zero as 𝑣 rel tends to zero, depending on the dependence of / p on 𝑣 rel . However, in our applications of the model so far, we have not encountered this problem.  \nFinally, we wish to resolve the crossing time of the cell by the particle to minimise the skipping of cells in its path, which would otherwise happen for very high relative velocity between the background wind and the cloud particles. We therefore define the timestep limit:  \n𝑡 cross = 𝑓 cross GLYPH<18> 3 𝑉 cell 4 𝜋 GLYPH<19> 1 3 𝑣 -1 rel , (14)  \nwhere 𝑓 cross is a free parameter and 𝑉 cell is the cell volume.  \nThe timestep limits defined above are applied to the cloud particle. However, if it happens that the host gas cell will be on a timestep longer than 𝑡 cellmod , we move it to a shorter timestep at the next available opportunity. 5 This is primarily so that the radiative cooling of the cell can respond to changes in its properties. This approach works well and required no modifications to A/r.pc/e.pc/p.pc/o.pc's cooling routines, although in future we could instead choose to sub-cycle the cooling. We conservatively adopt 𝑓 cellmod = 𝑓 grow = 𝑓 loss = 𝑓 stop = 0 . 1 and 𝑓 cross = 0 . 3 in this work, finding this to provide accurate integration (we demonstrate this in Section 2.2.3) while still not contributing noticeably to the cost of the simulation. Except for the early stages of the cloud particle's evolution, when the relative velocity with the ambient wind is at its highest, cloud particles are typically limited by 𝑡 grav .", '2.2.2 Implementation of the FB22 model': "FB22details a model for cloud-wind interactions that can be adopted by our cloud particle scheme. Clouds grow by accreting hot wind material via a turbulent radiative mixing layer (TRML) and lose mass by turbulent shredding. For the specifics of the model, its derivation and for physical interpretation of its predictions, we refer the reader to FB22 itself. Here, we give the essential outline of the model and describe its implementation within our cloud particle scheme.  \nFB22 assumes that the cloud is in pressure equilibrium with the ambient medium. We can therefore obtain the density contrast between a cloud and the ambient medium in terms of a specific internal  \nenergy contrast,  \n𝜒 ≡ 𝜌 cl 𝜌 cell = 𝑢 cell 𝑢 cl . (15)  \nHere, the density, 𝜌 cell , and specific internal energy, 𝑢 cell , of the host cell are known. The specific internal energy of the clouds within the cloud particle, 𝑢 cl , is set by the assumption that the cloud is in thermal equilibrium with photo-heating from the UV background and/or local sources. While this could be determined on-the-fly, for simplicity, in this work we follow FB22 and assume that the clouds have a temperature of 10 4 K. Thus, the density of the clouds, 𝜌 cl , can be determined from equation 15. For clarity, we emphasise that 𝜌 cl is the density of the material inside the identical subgrid clouds hosted within the cloud particle, not the total mass of clouds divided by the cell volume or any similar property.  \nThe characteristic radius of any of the subgrid clouds hosted by the cloud particle is  \n𝑟 cl = GLYPH<18> 3 𝑚 part 4 𝜋𝜌 cl 𝑁 cl GLYPH<19> 1 3 , (16)  \nwhere we have used 𝑁 cl to relate the particle mass to the mass of an individual cloud. The turbulent velocity within the TRML is taken to be  \n𝑣 turb = 𝑓 turb 𝑣 rel . (17)  \nFollowing FB22 we adopt 𝑓 turb = 0 . 1. We also need to determine the cooling time within the TRML. Using 'mix' to denote gas properties within the TRML, this is  \n𝑡 cool ≡ 𝑢 mix 𝑛 2 H , mix Λ ' mix . (18)  \nWith hydrogen mass fraction 𝑋 mix the hydrogen number density is 𝑛 H , mix = 𝑋 mix 𝜌 mix / 𝑚 p . Λ ' mix = max [ Λ ( 𝑢 mix , 𝜌 mix , 𝑋 mix , 𝑍 mix ) , 0 ] where Λ is the net cooling rate calculated using A/r.pc/e.pc/p.pc/o.pc's cooling routines. Λ ' mix = 0 in the event of net heating (i.e. Λ < 0 . 0), leading to an infinite cooling time. 6 This does not ever occur in this work, but might for lower choices of 𝑢 cl and/or in the presence of a strong radiation field. Following Begelman & Fabian (1990) and Gronke & Oh(2018), and with the TRML in pressure equilibrium with the cloud andtheambientmediumwehave 𝑢 mix = √ 𝑢 cell 𝑢 cl , 𝜌 mix = √ 𝜌 cell 𝜌 cl , 𝑋 mix = √ 𝑋 cell 𝑋 cl and 𝑍 mix = √ 𝑍 cell 𝑍 cl . As noted in FB22, what exactly sets 𝑡 cool in the TRML is still an open area of research, so other choices exist (see e.g. Abruzzo et al. 2022).  \nWe can then calculate the quantity  \n𝜉 = 𝑟 cl 𝑣 turb 𝑡 cool , (19)  \nwhich compares the relative strengths of turbulent mixing and radiative cooling. Referring the reader to FB22 for the derivation, the mass growth and loss rates are, respectively,  \n/ 𝑚 grow = / 𝑚 0 𝜉 𝛼 (20)  \nand  \n/ 𝑚 loss = / 𝑚 0 , (21)  \nwhere  \n/ 𝑚 0 = 3 𝑓 mix 𝑚 part 𝑣 turb 𝜒 1 / 2 𝑟 cl , (22)  \nwith 𝑓 mix = 1 / 3 our fiducial choice, and  \n𝛼 = ( 1 / 4 𝜉 ≥ 1 1 / 2 𝜉 < 1 . (23)  \nThese mass exchanges give rise to a momentum exchange, as expressed in equation 6. Additionally, we define the drag term (which enters into equation 5) as  \n/ p drag = -1 2 𝐶 drag 𝑁 cl 𝜋𝑟 2 cl 𝜌 cell 𝑣 rel v rel , (24)  \nwhere, as in FB22, we take 𝐶 drag = 1 / 2.  \nThe mass exchanges also give rise to an exchange of enthalpy between the particle and the cell, such that the net flow of thermal energy from cell to particle is  \n/ 𝐸 th = 𝛾 GLYPH<0> 𝑢 cell / 𝑚 grow -𝑢 cl / 𝑚 loss GLYPH<1> . (25)  \nAs previously mentioned, the temperature of the cloud particle is fixed; the specific thermal energy brought in by accreted material is radiated away in the TRML (this condition has set / 𝑚 grow in the first place) so the cloud particle does not gain the energy described by equation 25. However, the gas cell loses this energy and also experiences an additional heating rate as the kinetic energy lost by the particle in the rest frame of the cell is thermalised. Thus, the thermal energy of the cell changes by the following rate:  \n/ 𝐸 cell , th = 1 2 / 𝑚 loss 𝑣 2 rel -/ 𝐸 th . (26)  \nNote that the total energy of the cell also changes as a result of the momentum transfers (including the drag force).", '2.2.3 Validation of the FB22 implementation': "In order to validate our implementation, we present a simple numerical experiment. A single cloud particle is placed in a hot wind tunnel (similar in principle to many cloud-crushing simulations) and its subsequent evolution is compared to the FB22 solution. We choose the same configuration as FB22, section 3.2.3 (the results of which are shown in their fig. 4). The background wind has a velocity of 1000 km s -1 and a temperature of 10 7 K. Both the wind and clouds have a pressure of 𝑃 / 𝑘 B = 10 4 Kcm -3 and we use 𝜒 = 1000. To avoid having to implement A/r.pc/e.pc/p.pc/o.pc's cooling functions into our separate routines for directly integrating the FB22 model (or vice-versa), for this specific test we choose to adopt a fixed value of 𝑡 cool = 1 Myr (whichis of the correct order of magnitude for the cooling functions in the regime probed) to enable an exact evaluation of the performance of A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc with respect to the FB22 solutions. The background gas is not permitted to radiatively cool (other than onto the cloud particle) as a constant background is necessary for comparison to the idealised solutions. Likewise, we do not include any self-gravity. We perform several simulations with different initial values of 𝜉 for the cloud particles, each of which specifies an initial cloud mass (or, equivalently, an initial cloud radius).  \nWe choose a target gas cell resolution of 8 × 10 4 M ⊙ for the background wind. The domain has dimensions of 100 kpc × 25 kpc × 25 kpc with periodic boundary conditions. Initial positions of the mesh generating points are drawn from a low discrepancy sequence (the 𝑅 3 sequence of Roberts 2018) in order to reduce Poisson noise while avoiding a highly structured configuration as much as possible. Cell masses are initialised by multiplying the target density with their initial volume so that any remaining noise in the initial mesh generating point configuration manifests as a scatter in cell  \nFig. 1 shows the results of these simulations. We show the evolution of the cloud mass and velocity both in absolute terms or normalised to the initial mass or wind velocity, respectively. The solid, coloured lines show the output of the simulations while the dotted lines show the direct integral of the FB22 equations using a standard RK45 integration scheme. 7 It can be seen that for a range of initial values of 𝜉 = 𝜉 0 , the simulations have essentially perfect agreement with the expected evolution, confirming that our implementation correctly integrates the cloud evolution mass, momentum and energy transfer rates. For a detailed physical interpretation of the behaviour of the  \n<!-- image -->  \nFigure 1. The evolution of cloud particles placed in a wind tunnel, for various initial values of 𝜉 = 𝜉 0 (corresponding to different initial cloud masses). The wind has a velocity of 1000 km s -1 , temperature of 10 7 K and pressure of 𝑃 / 𝑘 B = 10 4 Kcm -3 . The density contrast between the clouds and the wind is 1000. To enable a one-to-one comparison to the FB22 solutions we adopt a fixed value of 𝑡 cool = 1 Myr. Solid lines indicate the results from the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc simulations while the dotted lines indicate the directly integrated solutions from FB22. We run the simulations for 100 initial cloud crushing times ( 𝑡 cc ) or until the particles lose 99 per cent of their mass. The latter case is indicated by a cross marker. Top left : cloud mass normalised to the initial cloud mass as a function of time normalised to 𝑡 cc . Bottom left : cloud velocity normalised by the wind velocity as a function of time normalised to 𝑡 cc . Top right : cloud mass as a function of time. Bottom right : cloud velocity as a function of time. It can be seen that the implementation of the FB22 cloud evolution model into A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc results in essentially perfect agreement with the directly computed solutions across a wide range of parameters.  \n<!-- image -->  \nmass rather than density. This configuration is then evolved for several domain sound crossing times with A/r.pc/e.pc/p.pc/o.pc's standard refinement, de-refinement and mesh regularisation schemes enabled, yielding a completely uniform density medium represented with an unstructured mesh and well-rounded cells which all have masses within the standard A/r.pc/e.pc/p.pc/o.pc tolerance of a factor of 2 of the target resolution. The cloud particle is given an initial mass of 800 M ⊙ (not to be confused with the initial subgrid cloud masses), placed in the domain and given an initial velocity of 1000 km s -1 (thus, in the rest frame of the particle, it is experiencing a 1000 km s -1 wind). We run the simulations for 100 cloud crushing times, where 𝑡 cc = 𝜒 1 / 2 𝑟 cl / 𝑣 rel , or until the cloud particle has lost 99 per cent of its initial mass.  \ncloud evolution model, we refer the reader to FB22. However, the experiment highlights the most salient points. Larger values of 𝜉 0 correspond to initially larger and therefore (for fixed 𝜒 ) more massive clouds. Large clouds (with 𝜉 0 ≫ 1) grow efficiently and are accelerated to a significant fraction of the wind velocity within a few 𝑡 cc . Small clouds (with 𝜉 0 ≪ 1) lose mass quickly. Small clouds evolve much faster than the larger clouds in real terms, undergoing a more rapid acceleration (despite the acceleration time-scale being longer relative to their own 𝑡 cc compared to larger clouds). Intermediate clouds with 𝜉 0 ≲ 1 lose mass to begin with, but the instantaneous 𝜉 becomes greater than unity as the cloud is entrained in the wind ( 𝑣 rel decreases) leading to subsequent net growth.", '2.3 Launching and recoupling of cloud and wind particles': 'In this work, we only consider the generation of clouds from within the ISM (rather than via precipitation from the CGM for example, although the scheme can also be applied in that context). In Paper I we described the launching of wind particles from the ISM in order to drive the hot phase of a galactic wind - the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc-H/o.pc/t.pc scheme. Cloud particles are launched following the same approach. Wetherefore refer the interested reader to Paper I, but summarise the details here, focusing primarily on the addition of the cloud particles to the scheme.  \nThe star formation rate of a cell, / 𝑚 ★ , is calculated in this work with the model of SH03, though our scheme is insensitive to this choice. The input wind and cloud mass loading factors, 𝜂 w0 𝑀 and 𝜂 cl0 𝑀 , relate the rate at which wind and cloud material is injected relative to the  \nstar formation rate:  \n/ 𝑚 w = 𝜂 w0 𝑀 / 𝑚 ★ , (27)  \n/ 𝑚 cl = 𝜂 cl0 𝑀 / 𝑚 ★ , (28)  \nWe sample / 𝑚 ★ , / 𝑚 w and / 𝑚 cl to stochastically generate star, wind and cloud particles from star forming gas cells. Star particles have the same mass resolution as the gas mass resolution of the simulation (so gas cells are typically converted completely into star particles), whereas we spawn wind and cloud particles that have a mass that is a factor 𝑓 w 𝑚 and 𝑓 cl 𝑚 , respectively, smaller than the gas mass resolution.  \nWind particles inherit the velocity of their parent gas cell, but are given a kick  \nΔ 𝑣 w = v u t 2 𝜂 wkin0 𝐸 𝜂 w0 𝑀 𝑢 ★ , (29)  \nwhere 𝜂 wkin0 𝐸 is the input wind kinetic energy loading and 𝑢 ★ is the characteristic specific energy associated with stellar feedback. Since the majority of the energy driving the winds originates from supernovae (while radiation, stellar winds etc. have a more local effect around massive stars), we adopt 𝑢 ★ = 5 . 26 × 10 5 ( kms -1 ) 2 , as in Paper I, corresponding to one SN of 10 51 erg for every 95.5 M ⊙ of stellar mass formed (consistent with the value used in Kim et al. 2020a). 8 The kick can be applied in some preferred direction (e.g. vertically out of the disc plane) or isotropically. This choice is discussed in detail in Paper I. For this work we choose to apply the kicks vertically. The wind particle receives a specific internal energy  \n𝑢 w = 𝜂 wth0 𝐸 𝜂 w0 𝑀 𝑢 ★ , (30)  \nwhere 𝜂 wth0 𝐸 is the input wind thermal energy loading.  \nIn this work we choose to set the initial velocity and internal energy of the cloud particles directly (i.e. this means they are independent of the mass loading). The specific internal energy of the clouds are set by the cloud model. As described above, for the FB22 model adopted in this work, this means that we give the clouds an initial temperature of 10 4 K.Notethatwecannotresolve the acceleration within the ISM, so the choice of the initial launch velocity is intended to compensate for this. In the future, we will derive this value more rigorously from high resolution simulations in tandem with an appropriate mass loading for cool gas that is able to escape the ISM (possibly leading to velocities that depend on halo properties). However, in this work, we choose a value of Δ 𝑣 cl = 100 kms -1 which we have determined empirically to be high enough to allow all cloud particles to escape the ISM in our idealised setup while still providing plenty of headroom for the interaction with the hot wind to accelerate them to higher velocities. We will demonstrate later in this work that this initial velocity kick is negligible compared to the subsequent cloud acceleration by the hot wind.  \nBoth wind and cloud particles are initially hydro-decoupled, meaning that they only experience gravity and do not participate in any hydrodynamical interactions (either explicitly or via a subgrid model). Once a wind particle finds itself in gas with a density lower than a threshold, 𝜌 rec , it recouples, depositing its mass, energy, momentum and metals into the gas of the ISM/CGM transition region and drives  \n8 It should be noted that 𝑢 ★ plays the role of a normalising reference value and is degenerate with the energy loadings. Its exact value is not particularly important, except when comparing energy loadings between works.  \na hydro-resolved hot wind. We choose 𝜌 rec = 0 . 1 𝜌 SF . The details of the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc-H/o.pc/t.pc wind recoupling model are described in detail in Paper I. The most salient point, where we have made significant improvements over previous similar models, is that we are able to ensure that we have sufficient resolution to resolve a low density, high specific energy outflow with the use of a novel refinement technique (\'displacement recoupling\'). The wind is initially kept at the same mass resolution as the wind particles (i.e. a factor 𝑓 𝑚, w better than the base resolution of the simulation). Wind material subject to the refinement criteria is identified with a passive scalar "dye" injected as the wind particle recouples. In this work, we gradually relax the refinement of the wind such that the mass resolution increases linearly from 0 . 1 𝑟 200 from the galaxy centre until it reaches the base resolution of the simulation at 0 . 5 𝑟 200 . A more physically motivated version of this distance-refinement level relation implemented for use in cosmological simulations will be presented in Bennett et al. in prep.  \nIf the ambient density around a cloud particle is greater than 𝜌 rec and the density contrast between the clouds and the ambient medium is smaller than a threshold value 𝜒 dec , it also remains completely hydro-decoupled. The first time that a cloud particle finds itself in gas such that 𝜌 cell ≤ 𝜌 dec and 𝜒 ≥ 𝜒 dec , the full cloud evolution model is enabled and exchanges of mass, momentum, energy and metals between the cloud and the host cell begin, as described in the previous sections. Once the cloud model has been enabled, if the density contrast ever falls below a second threshold value, 𝜒 rec , or its mass drops below a factor 𝑓 𝑚, rec of its initial value, it is fully recoupled into the cell. This proceeds via the \'standard recoupling\' scheme described in Paper I; all conserved quantities of the cloud particle are injected into the host cell. We adopt 𝑓 𝑚, rec = 0 . 1. We choose 𝜒 rec = 10, causing the cloud particle to recouple if it has encountered gas with densities not too dissimilar to the subgrid internal cloud density. We use a slightly higher value to enable the cloud evolution model, choosing 𝜒 dec = 1 . 1 𝜒 rec = 11, such that cloud particles are sufficiently clear of the recoupling threshold to avoid \'false starts\'.  \nAs described below, in this work all cloud particles are launched with sufficient velocity to exit the ISM and meet the conditions to turn on the cloud evolution model. In other settings (such as cosmological simulations), it is in principle possible for particles to remain trapped in the ISM in which case we will fall back to a maximum lifetime before recoupling. However, because our model is designed to treat winds that have left the ISM, the launching of cloud particles with low velocities that will always remain fully hydro-decoupled should be avoided as much as possible. When measuring outflow properties from high resolution simulations for use as inputs to A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc, the mass loading of the cool phase should only consider contributions from material that has a high enough velocity to leave the ISM.', '3 DEMONSTRATION OF THE MODEL': 'We will now apply the model using idealised non-cosmological simulations. The aim of this work is to demonstrate the general behaviour of the model in various regimes, rather than to advocate for a particular set of loading factors. A more comprehensive study of the effect of different loading factors will be the focus of future works in a cosmological context.', '3.1 Initial conditions and model setup': 'We use the fiducial initial conditions from Paper I and refer the interested reader there for the finer details of their preparation. The  \nTable 1. The parameters defining the initial conditions used in this work. For details about the models adopted and the definition of the symbols, see the main text. Note that the quantities reported in this table are all input parameters with the exception of 𝑟 200 , which is derived from the halo mass, concentration and cosmology, and the CGM mass inside 𝑟 200 , which is a derived quantity of the cooling flow solution given the other constraints.  \nsetup is an isolated system comprised of dark matter, a disc and bulge of pre-existing stars, a gas disc and a realistic CGM. The latter component is often omitted from idealised galaxy simulations but is crucial for this work, allowing us to study accretion onto the galaxy, the interaction of winds with the CGM and the balance of ejective and preventative feedback. The parameters describing the initial conditions are given for reference in Table 1. The dark matter is modelled as a spherically symmetric, static background potential and includes both an inner and outer halo component. The inner component follows a Navarro-Frenk-White (NFW, Navarro et al. 1997) profile with 𝑀 200 = 10 11 M ⊙ and a concentration of 10. Taking a Planck Collaboration et al. (2020) cosmology, this gives 𝑟 200 = 97 . 9 kpc (where we define 𝑀 200 and 𝑟 200 relative to the critical density at 𝑧 = 0). The outer halo component is modelled following Diemer & Kravtsov (2014), with the parameters 𝑠 𝑒 = 1 . 5 and 𝑏 𝑒 = 1. The stellar disc and bulge, and the gas disc are generated using the code M/a.pc/k.pc/e.pcN/e.pc/w.pcD/i.pc/s.pc/k.pc (Springel et al. 2005). The stellar and gas discs have exponential surface density profiles with a scale length of 𝑅 s = 2 . 5 kpc. The stellar disc has a Gaussian vertical density profile with a scale height of 𝑧 s = 0 . 25 kpc. The gas disc has a vertical density profile that is set to produce hydrostatic equilibrium at its initial temperature of 𝑇 0 = 10 4 K. We truncate the gas disc beyond five scale lengths and five scale heights. The disc has an  \ninitial metallicity of 𝑍 0 = 1 𝑍 ⊙ (where we adopt 𝑍 ⊙ = 0 . 0127). The stellar bulge is spherically symmetric and follows a Hernquist (1990) density profile with a scale length of 𝑟 s = 0 . 25 kpc. We use 𝑀 disc ,★ = 8 × 10 9 M ⊙ , 𝑀 bulge ,★ = 10 8 M ⊙ and 𝑀 disc , gas = 2 × 10 9 M ⊙ .  \nThe CGM gas is initialised to a steady state rotating cooling flow configuration for our total potential, the full details of which can be found in Stern et al. (2024). We choose a cooling flow solution with a sonic radius of the flow of 𝑟 sonic = 2 kpc and a circularisation radius of 𝑟 circ = 2 . 5 kpc. This means that the flow remains in the subsonic limit. We choose an initial CGM metallicity of 0 . 1 𝑍 ⊙ . In combination with the other constraints, this yields an initial CGM mass inside 𝑟 200 of 2 . 96 × 10 9 M ⊙ . The predicted steady-state mass flux of the cooling flow (in the absence of feedback) for these initial conditions is 0 . 2 M ⊙ yr -1 , though in practice the addition of a galaxy to the centre of the halo and minor discrepancies between the cooling functions used in A/r.pc/e.pc/p.pc/o.pc and those assumed when calculating the solution leads to the emergent inflow rate being a factor of a few higher. In order to provide a reservoir of mass that can flow into the halo and maintain the cooling flow, we initialise the CGM gas out to a distance of 600 kpc ( ∼ 6 𝑟 200 ). However, we do not require our full resolution far outside the halo. We therefore degrade the mass resolution smoothly outside 200 kpc by a factor of 3 every √ 2 × 200kpc, making appropriate modifications to the (de)refinement scheme to maintain this resolution structure during the simulation.  \nWe use a target gas cell mass resolution of 𝑚 g , tar = 8 × 10 4 M ⊙ . Recall that, as described above, a given gas cell may be at a finer resolution if subject to the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc hot wind refinement scheme or a coarser resolution if beyond 200 kpc from the system centre. Star particles, either present in the initial conditions or created during the simulation, also have a mass of 𝑚 ★ = 8 × 10 4 M ⊙ . Gas cells have an adaptive gravitational softening of 2.5 times the cell radius with a minimumvalueof50pc.Collisionlessparticles have a fixed softening of 195 pc.  \nFor the background hot wind, we take the input loading parameters for the high specific energy wind presented in Paper I, 𝜂 w0 𝑀 = 0 . 32, 𝜂 wkin0 𝐸 = 0 . 321, 𝜂 wth0 𝐸 = 0 . 579. We use 𝑓 w 𝑚 = 0 . 01 i.e. the wind particles and refined wind gas are 100 times less massive than the base resolution of the simulation. For the cool cloud component, we perform simulations with three different input mass loadings, 𝜂 cl0 𝑀 = 0 . 1, 1 and 5, which illustrate different regimes of the model. As stated earlier, in this work we give cloud particles an initial velocity kick of 100 km s -1 and temperature of 10 4 K. For the three mass loadings above, this corresponds to total input cloud energy loadings of 𝜂 cl0 𝐸 = 9 . 7 × 10 -4 , 9 . 7 × 10 -3 and 0.048, respectively. Note that these are negligible compared to the input energy loadings of the hot wind. We use 𝑓 cl 𝑚 = 0 . 01 i.e. cloud particles are initially 100 times lower mass than the base resolution of the simulation (but the same mass resolution as the wind particles and refined wind gas). Newly created cloud particles are assigned a randomly drawn cloud number, 𝑁 cl , such that the initial distribution of (subgrid) cloud masses across the whole population of cloud particles follows a mass function d 𝑁 / d 𝑚 ∝ 𝑚 -2 , motivated by the mass functions seen in the high resolution simulations of Tan & Fielding (2024). We also perform simulations with no cool cloud component (i.e. 𝜂 cl0 𝑀 = 0) and no wind at all (i.e. 𝜂 w0 𝑀 = 𝜂 cl0 𝑀 = 0).  \nFor the sake of simplicity, in all simulations in this work, wind and cloud particles inherit the metallicity of the ISM gas cell from which they were spawned. We will include differential levels of enrichment in future work (derived from high resolution simulations such as Kim et al. 2020a), which we expect to have significant impact on the metallicity of different CGM phases.', '3.2.1 Morphologies': "In Fig. 2 and 3 we show visualisations of the simulations with cloud particles after 0.75 and 2.25 Gyr, respectively. In the top row, we show a gas density slice 9 orientated vertically through the disk plane. The underlying structure of the Voronoi mesh is apparent in these slices; recall that these simulations are carried out at a resolution that is feasible for cosmological volume simulations, not just for individual galaxies. Below this, we plot the locations of cloud particles coloured by their current subgrid cloud mass. The particles shown are all those that are located inside the cells which are shown in the gas density slice. The third and fourth rows show the radial velocity of the gas and cloud particles, respectively. The bottom row shows a slice of the intra-cell cloud particle 3D velocity dispersion, 𝜎 cl . For each cell that contains cloud particles, we calculate 𝜎 cl as the mass-weighted velocity dispersion of all of their cloud particles.  \nBroadly speaking, Fig. 2 and 3 neatly illustrate three distinct scenarios for the impact of the cloud population on the background wind as the input cloud mass loading is varied: 1) there are relatively few clouds in the outflow so the background wind is barely impacted by the cloud component, 2) there are large amounts of cloud material in the outflow but the background wind is still coherent and energetic, and 3) a high cloud mass loading significantly disrupts the background wind which in turn limits the spread of cloud material within the CGM. We will now explore these results in more detail.  \nStarting with the 𝜂 cl0 𝑀 = 0 . 1 simulation (left-hand columns), in Fig. 2 (at 0.75 Gyr), the outflowing wind can be clearly seen in the density and radial velocity slices. The outflow forms a coherent bicone that is significantly under-dense compared to the CGM. The outflow is fastest up the centre of the bicone. The outflow covers a significant portion of the CGM, restricting inflows (and generally regulating the galaxy evolution as described in Paper I). By 2.25 Gyr (see Fig. 3) the opening angle of the outflow has reduced but is still a coherent flow. This effect was seen and discussed in detail in Paper I in the absence of cloud material. As the SFR of the galaxy drops (as we will examine below), the power of the wind similarly drops. This allows the inflowing CGM to confine the outflow, leading to a narrowing of the wind bicone. In both Fig. 2 and 3, a handful of cloud particles can be seen with a wide variety of cloud masses. At 0.75 Gyr, these particles largely trace the wind region and are travelling outwards at many hundreds of km s -1 (we will examine the outflow velocities quantitatively later in the work), much faster than their initial launch velocity. Cloud material is able to make it to 𝑟 200 and beyond. It can be seen that there is a slight under-density in cloud particles in the centre of the outflow, where the velocity of the wind is highest. By 2.25 Gyr, there are fewer cloud particles and they are more preferentially located closer to the boundary of the wind. In fact, some are actually located outside of the wind. This arises because cloud particles are not perfectly 'pinned' to the wind material. The wind velocity field can fluctuate faster than the clouds can respond, due to their inertia, so clouds can fall out of the flow. Furthermore, the wind is in general faster than the clouds so as the wind bicone narrows, cloud particles can find themselves outside of the wind. If the clouds survive the sudden change in relative velocity as they pass into oncoming CGM material, they can then fall back towards the galaxy, as can be seen for some particles.  \nAt 0.75 Gyr (Fig. 2), it can be seen that the hot wind in the  \n𝜂 cl0 𝑀 = 1 simulation (centre columns) is similar to the 𝜂 cl0 𝑀 = 0 . 1 case, albeit slower. It still occupies a large fraction of the CGM. Many more cloud particles are present, spanning the wind bicone. Again, there is a wide variety of cloud masses (we will examine the mass distribution in a quantitative manner below) and cloud material is able to be accelerated out of the halo. Already at this time, small amounts of inflowing cloud material is apparent. These are cloud particles launched at earlier times that have ended up outside of the wind. By 2.25 Gyr (Fig. 3), the wind has begun to be disrupted by regions of dense material that have formed out in the CGM. These divert the flow, leading to a complex velocity field. Accordingly, the morphology and velocity field of the ensemble of cloud particles is more complex, with outflowing and inflowing regions (that often overlap).  \nAt 0.75 Gyr, the wind in the 𝜂 cl0 𝑀 = 5 simulation (right-hand columns) is already very different to the other simulations. The wind bicone is much narrower and a dense shell of slow moving gas has formed on the edge of the outflow. Much of the CGM is dominated by inflowing material. There are a large number of cloud particles in the immediate vicinity of the galaxy that have not been efficiently accelerated away. In fact, there are a handful of particles falling inwards close to the disc. The particles that have been accelerated away are limited to the region on the edge of the wind. By 2.25 Gyr, the quantity of dense material out in the CGM has increased, significantly diverting the wind. This also leads to a complex spatial distribution of cloud particles as well as inflowing components. Very close to the disc, there are a mixture of outflowing and inflowing cloud particles.  \nTurning our attention to the maps of 𝜎 cl in Fig. 2 and 3, it can be seen that the velocity dispersion of cloud particles contained within a single cell can reach over 100 km s -1 . Thus, on the spatial scale of our gas cell resolution, the motion of the cloud material is not well described by a single valued velocity field. There are two main causes of this non-zero velocity dispersion. Firstly, there are regions where populations of outflowing cloud particles intersect populations of inflowing cloud particles. This can occur out in the CGM (for example, this can be seen at 2.25 Gyr for the 𝜂 cl0 𝑀 = 1 simulation) or close to the galaxy in a low altitude fountain flow (this is seen especially in the 𝜂 cl0 𝑀 = 5 simulation). The second, often more dominant, cause of this velocity dispersion is the differential acceleration of clouds of different masses when exposed to the same background wind. As we demonstrated in Fig. 1, lower mass clouds are accelerated faster than higher mass clouds. The values of 𝜎 cl of ∼ 50 -120 km s -1 seen within the wind (particularly apparent in Fig. 2 for the 𝜂 cl0 𝑀 = 1 simulation, but present in all runs) occur because lower mass clouds overtake higher mass clouds. Both of these effects mean that the dominant component of the dispersion is in the radial direction for most cells (though a more realistic, less spherically symmetric setup might alter this).", '3.2.2 SFRs, inflow/outflow rates and loading factors': 'In Fig. 4 we plot the SFR as a function of time for the three simulations presented above, as well as for the simulation with no wind ( 𝜂 w0 𝑀 = 𝜂 cl0 𝑀 = 0) and the simulation with a wind but no cloud component ( 𝜂 cl0 𝑀 = 0). We derive the SFR from the mass of stellar material formed within 20 Myr windows. Equivalents of the no wind and 𝜂 cl0 𝑀 = 0 simulations have been presented in Paper I and we refer the reader to that work for a detailed analysis. In the no wind simulation, the consumption of ISM gas is balanced by inflows from the CGM, resulting in a steady SFR close to 1 M ⊙ yr -1 . With  \n0', 't = 0 . 75 Gyr': '1  \nFigure 2. Visualisations of the fiducial simulations with cloud particles after 0.75 Gyr. The left, centre and right columns show the 𝜂 cl0 𝑀 = 0 . 1, 1 and 5 simulations, respectively. The top row shows slices of gas density, orientated vertically through the disc plane. The second row shows the locations of all cloud particles contained within cells intersected by the slice shown in the top row, coloured by their current subgrid cloud mass. The third and fourth rows show the radial velocity for gas and cloud particles, respectively. The fifth row shows slices of the intra-cell cloud particle 3D velocity dispersion (see the text for details). On all panels, we overlay dashed circles corresponding to 0.1, 0.2, 0.5 and 1 times 𝑟 200 . These are the reference surfaces used later in the work to measure mass and energy fluxes. On the top row, we also indicate (again with dashed lines) the region used when measuring various wind property profiles later in the work. The three simulations shown demonstrate different scenarios. On the left, a low cloud mass loading means that the background wind experiences a minimum of disruption, but there is only a small population of cloud particles. In the centre, a higher mass loading results in only a slight impact on the wind, but the wind is now filled with cloud particles. On the right, a mass loading of 5 has lead to significant disruption of the background wind, which has in turn limited the spatial extent covered by cloud particles. In the bottom row, we can see that the different trajectories of cloud particles means that the cloud velocity field would be poorly described by a single value per cell (which would enforce 𝜎 cl = 0).  \n<!-- image -->  \n1  \n0  \nFigure 3. Visualisations of the fiducial simulations with cloud particles after 2.25 Gyr. See the caption for Fig. 2 for details on what is displayed. The bicone of the outflow in the 𝜂 cl0 𝑀 = 0 . 1 simulation (left) has narrowed but is largely undisrupted. Both the 𝜂 cl0 𝑀 = 1 and 5 (centre and right, respectively) show significant disruption to the wind bicone, with the outflow being diverted by denser gas structures. This is more severe in the 𝜂 cl0 𝑀 = 5 case. The complex velocity field of the background wind (which also varies significantly in time) leads to cloud particles falling out of the flow. There are therefore populations of inflowing cloud particles.  \n<!-- image -->  \nFigure 4. SFR as a function of time for simulations without a wind and for different input cloud mass loading factors. We derive the SFR as the mass of stellar material formed within a 20 Myr window. Without a wind, the SFR remains steady. With the inclusion of a wind, inflows from the CGM are suppressed, leading to a reduction in the SFR as the ISM is consumed by star formation. Increasing the cloud mass loading factor results in further reductions in the SFR as ISM material is ejected from the galaxy. However, for the 𝜂 cl0 𝑀 = 5 simulation, the preventative feedback of the wind becomes ineffective at around 0.8 Gyr, leading to a return to steady-state star formation.  \n<!-- image -->  \na wind (but no cloud component), the SFR starts off close to the no wind case, but drops steadily, reaching 0 . 1 M ⊙ yr -1 by 2 Gyr. As shown in Paper I, the reduction of the SFR occurs because the wind suppresses the resupply of gas from the CGM. The wind is therefore a form of preventative feedback. The 𝜂 cl0 𝑀 = 0 . 1 simulation produces essentially the same SFR evolution, but the 𝜂 cl0 𝑀 = 1 simulation produces a marginally lower SFR. This is because, in addition to the preventative feedback of the wind, the launching of clouds provides an element of ejective feedback, removing star forming gas from the disc. This is strongest in the 𝜂 cl0 𝑀 = 5 simulation. It can be seen that the SFR initially declines much faster than the other simulations. However, around 0.8 Gyr the SFR suddenly increases and then reaches a steady state at around 0 . 5 M ⊙ yr -1 . This indicates that the preventative feedback of the wind has suddenly become ineffective, as we shall explore when examining the properties of the outflows.  \nIn Fig. 5 we plot the mass inflow and outflow rate through four reference spheres around the galaxy with galactocentric radius 𝑟 of 0.1, 0.2, 0.5 and 1 𝑟 200 for the simulations shown in Fig. 4. We also plot the SFR of each simulation for reference. We calculate and plot the mass flow rates for the gas cells (the \'hydro\' component) and cloud particles separately, as well as plotting their sum. For the gas cells, following Paper I, we determine the mass flow rates by considering the fluxes through a sphere at the reference radius. This is achieved by discretizing the reference sphere into 𝑁 pix equal area pixels using the HEALP/i.pc/x.pc library (Górski & Hivon 2011). In order to guarantee that the inter-pixel spacing is finer by a factor of a few than the diameter of the smallest cells intercepted by the sphere, we adopt 𝑁 pix = 786432. For each pixel centre, we search for the nearest Voronoi mesh-generating point and hence the gas cell within which the pixel is located. Gas cell properties are then mapped on the pixel. The mass flux per unit area through the pixel is then  \nF 𝑀 = 𝜌𝑣 r , (31)  \nwhere 𝜌 is the cell density and 𝑣 r is the radial velocity (i.e. normal to the spherical reference plane). We can then select only pixels with  \npositive (negative) 𝑣 𝑟 to compute the mass outflow (inflow) as  \n/ 𝑀 hydro out ( in ) = 𝐴 ∑︁ F 𝑀 , (32)  \nwhere we sum over the selected pixels, each with equal area 𝐴 = 4 𝜋𝑟 2 / 𝑁 pix . For cloud particles we obtain the mass outflow (inflow) rates as  \n/ 𝑀 cl out ( in ) = " 𝑚𝑣 r Δ 𝑟 , (33)  \nwhere the sum runs over all cloud particles with positive (negative) radial velocity located with in a spherical shell of thickness Δ 𝑟 . We adopt Δ 𝑟 = 2 kpc, 4 kpc, 10 kpc and 20 kpc at 𝑟 = 0.1, 0.2, 0.5 and 1 𝑟 200 , respectively. In contrast to the cloud particles, all wind particles have already recoupled well inside 0 . 1 𝑟 vir .  \nAs shown in Paper I, in the simulation with no wind (top row), inflows are completely dominant at all four radii. The inflow rate through 0.1 𝑟 200 remains constant at about 1 M ⊙ yr -1 over 2.5 Gyr. The SFR is essentially identical to this as it is regulated by resupply of material from the CGM to the ISM. Inflows at larger radii are not quite as constant, due to instabilities arising from our idealised setup. A small outflow can be seen at 0.1 𝑟 200 , but this only arises from disc material crossing the reference sphere. These should be compared against the other simulations in order to assess the level of \'true\' wind driven mass outflows. A very small outflow can be seen at late times crossing 0.5 and 1 𝑟 200 . This is an outwardly propagating sound wave arising from the idealised setup and has negligible impact. With the addition of the high specific energy wind component, but no cloud material ( 𝜂 cl0 𝑀 = 0, second row), inflows are efficiently suppressed at all radii. At 0.1 𝑟 200 , the outflow and inflow rates are roughly equal. Further out, outflow rates exceed inflow rates for the first 1 - 1.5 Gyr. The suppression of inflows acts as a form of preventative feedback, as described above, causing the SFR to drop as the ISM gas supply is used up. As examined in detail in Paper I, the dropping SFR causes a gradual decrease in wind power in absolute terms, resulting in inflows beginning to reassert themselves at large radii in the latter half of the simulation. It is worth noting that the outflow rate exceeds the SFR at all times and at all measured radii despite the input wind mass loading, 𝜂 w0 𝑀 = 0 . 32, being less than unity. In other words, the outflow contains more mass than is being ejected from the ISM, indicating entrainment of CGM material; this will be discussed further below.  \nWhen cool clouds are introduced, but at a small mass loading ( 𝜂 cl0 𝑀 = 0 . 1, third row), the results are very similar to the 𝜂 cl0 𝑀 = 0 simulation; the cloud phase is subdominant to the resolved hydrodynamic phase, so has little impact. However, it can be seen that there is outflowing cloud material at all radii, albeit at a much lower level than the rest of the gas. Increasing the input mass loading of the cloud material to 𝜂 cl0 𝑀 = 1 (fourth row) results in a similar picture. The hot wind continues to efficiently suppress inflows. A larger fraction of the total outward flowing mass flux is now in the form of cloud particles, scaling approximately linearly with the increase in the input mass loading. However, it is still the subdominant component, despite having over three times the input mass loading than the hot wind. This is because the hot wind is able to entrain significant amounts of CGM material, boosting the total mass of outflowing material well above the injected hot wind mass.  \nFurther increasing the cloud input mass loading to 𝜂 cl0 𝑀 = 5(bottom row) results in a change in behaviour. Initially, inflows are suppressed in common with the other simulations that include a wind. The total mass outflow rate is marginally higher through 0.1 𝑟 200 than the other cases, primarily (as we shall examine later) due to an increase in the amount of stripped cloud material in the wind. The SFR drops  \nFigure 5. SFR, mass inflow and outflow rates through spherical reference surfaces as a function of time. The different columns correspond to measurements taken at 0.1, 0.2, 0.5 and 1 𝑟 200 (the SFR is repeated in each column for reference). The locations of these reference surfaces are indicated in Fig. 2 and 3. Each row corresponds to a different simulation, with various input loadings. We plot the total inflow/outflow rates with solid lines and use dotted and dashed lines to indicate the contribution from the resolved hydro and cloud particles, respectively. Without any wind, the SFR remains high, due to CGM inflows. When a hot wind is included, unless 𝜂 cl0 𝑀 is very large (bottom row), strong preventive feedback reduces accretion, with the SFR remaining lower than the hydrodynamic outflows through all the measurement surfaces. At all radii and in all simulations, the cloud mass fluxes remain well below the hydro mass fluxes, even when 𝜂 cl0 𝑀 ≫ 𝜂 w0 𝑀  \n<!-- image -->  \n.  \ndue to preventative feedback, with the addition of a slight additional reduction due to the higher input mass loading causing more ISM material to be ejected. However, from around 0.5 Gyr onwards the inflow rates surpass the outflow rates in the inner radii, leading to an uptick in the SFR around 0.8 Gyr as the impact of the restored supply of gas to the ISM is felt. For the rest of the simulation, the wind is not able keep inflow rates below outflow rates at all radii, resulting in the resumption of a steady state SFR, albeit at a reduced level than the no wind case.  \nThe relative behaviour of the different simulations can be further understood by examining the emergent mass and energy loadings of the outflows. For both the resolved hydrodynamic and cloud particle outflows, in an analogous manner to the definition of the input loadings given in Section 2.3, we define the emergent mass loading as the ratio of the instantaneous 10 mass outflow rate through the reference sphere to the SFR averaged over the previous 20 Myr:  \n𝜂 hydro M = / 𝑀 hydro out / 𝑀 ★ , (34)  \n𝜂 cl M = / 𝑀 cl out / 𝑀 ★ . (35)  \nWecansimilarly define energy loadings by examining the ratio of the energy outflow rate to the energy input associated with star formation:  \n𝜂 hydro E = / 𝐸 hydro out 𝑢 ★ / 𝑀 ★ , (36)  \n𝜂 cl E = / 𝐸 cl out 𝑢 ★ / 𝑀 ★ . (37)  \nThe energy outflow rate for the resolved hydro phase is computed in a similar manner to the mass outflow rate. The energy flux through each HEALP/i.pc/x.pc pixel is  \nF 𝐸 = 𝜌𝑣 𝑟 GLYPH<18> 1 2 𝑣 2 + 1 𝛾 -1 𝑐 2 s GLYPH<19> , (38)  \nfor magnitude of the total velocity, 𝑣 , and sound speed, 𝑐 s = √︁ 𝛾𝑃 / 𝜌 . The total energy outflow rate is then  \n/ 𝐸 hydro out = 𝐴 ∑︁ F 𝐸 , (39)  \nwhere the sum runs over all pixels with positive 𝑣 𝑟 . The equivalent quantity for the cloud particles is  \n/ 𝐸 cl out = " 𝑚𝑣 r GLYPH<16> 1 2 𝑣 2 + 𝑢 GLYPH<17> Δ 𝑟 , (40)  \nwhere 𝑢 is the specific internal energy of the cloud and the sum runs over all cloud particles with positive radial velocity located with in a spherical shell of thickness Δ 𝑟 .  \nWe plot the mass and energy loadings for the simulations (except the no wind case, which has no appreciable outflows) in Fig. 6. For the simulation with no cloud material, the mass loading is between  \n10 For simplicity, we make the comparison between the instantaneous outflow rates and the very recent SFR. In order to more carefully examine the link between star formation and emergent outflows, one could attempt to account for the travel time of the outflow from the galaxy to the reference sphere. However, the outflows lack a single characteristic velocity, which would complicate this procedure.  \n∼ 1 -8 at all times and at all radii. Given the input mass loading of the hot wind is only 0.32, this indicates that a significant amount of inflowing CGM material is being turned around. The energy loading is relatively constant at all radii close to the input value of 0.9, indicating minimal losses. The 𝜂 cl0 𝑀 = 0 . 1 and 1 simulations do not substantially differ from the simulation without cloud material with regards to the (hydro) mass loading. There is a very slight enhancement overall which originates from material stripped from cloud particles. For both the 𝜂 cl0 𝑀 = 0 . 1 and 1 cases, the mass loading is a factor ∼ 3 to 10 below the injection value, decreasing in time. The missing mass is a result of a combination of shredded clouds (which boosts the hydro mass loading) and clouds which fail to reach the reference sphere. The latter contribution increases with time as the opening angle of the hot wind decreases (as described above), leading to more cloud particles falling out of the flow at low altitude or never being entrained to begin with. The input energy loading in the cloud phase is very small, corresponding to the initial launch velocity of 100 km s -1 and temperature of 10 4 K; we plot these input loadings with a horizontal dashed line on Fig. 6. It can be seen that the emergent energy loading of the clouds, 𝜂 cl E , is initially ∼ 5 times larger than the input at all radii. This is due to the efficient acceleration of clouds by the hot wind. The signature of this transfer of energy from the high specific energy wind to the clouds can be seen in the emergent energy loadings of the resolved hydro phase, 𝜂 hydro E ; there is a corresponding deficit relative to the case with no cloud material. 𝜂 cl E decreases with time in a similar manner to 𝜂 cl M as the mass of cool clouds being entrained in the wind drops.  \nThe simulation with a higher input cloud mass loading of 5 exhibits different behaviour. For the first ∼ 1 Gyr of the simulation, the emergent mass loading in the hydro phase is up to a factor of a few larger than the other simulations at 0.1 𝑟 200 . Contrastingly, the mass loading in the cloud phase is much lower than the input. As we shall confirm later, this indicates that the enhancement in the resolved hydro phase comes primarily from shredded cloud material. However, the emergent energy loading in the hydro phase is almost an order of magnitude lower than the input value. The high specific energy wind loses kinetic energy to the clouds. None the less, it can be seen that the emergent energy loading in the clouds is actually lower than the input value and does not explain the deficit (as is the case with lower input cloud mass loadings). The hot wind also experiences enhanced radiative cooling losses in two forms. Firstly, the increased mass of clouds in the wind results in greater losses via the (subgrid modelled) TRML of clouds. Unlike the kinetic energy transfer, this energy is not gained by the clouds but is instead lost from the wind. Secondly, the large amounts of mass stripped from the clouds into the hot wind increases the radiative cooling within the hot phase itself by increasing its density and dropping its temperature towards higher valued regions of the cooling function. Thus, the increased input mass loading of clouds results in a \'poisoning\' of the wind, leading to a sudden drop in the total mass loading at 0.2 𝑟 200 and beyond from around 0.5 Gyr onwards. The wind\'s capacity to perform preventative feedback is blunted, resulting in the restoration of inflows and subsequent star formation described above.', '3.2.3 Wind profiles': "In Fig. 7, we show various radial profiles of the wind region. These are calculated in a truncated cone expanding upwards from the disc with a cylindrical radius of 𝑅 = 2.5 kpc (equal to the disc scale length) in the disc mid-plane and an opening angle of 10 degrees. This region is illustrated in the density slices shown in Fig. 2 and 3.  \nFigure 6. Mass loading and energy loadings through spheres at various radii for simulations with different input cloud mass loadings. We split the contributions to the loadings into the resolved hydro component (first and third rows) and cloud particle component (second and fourth rows). Horizontal dashed lines indicate the input values (for the energy loadings, this is the combined kinetic and thermal). The simulation without a cool cloud component shows a mass loading that is around an order of magnitude larger than the input value at all radii, indicating that the outflow is entraining a substantial amount of CGM gas. The energy loading remains close to the input value. Adding clouds leads to some additional mass loading in the hydro component at the expense of the cloud mass loading, as material is stripped. Likewise, the cloud energy loading is increased from the input value by taking energy from the hydro phase. When the input cloud mass loading is increased to 5, this results in too much energy being removed from the hot wind to maintain a strong outflow.  \n<!-- image -->  \nFor the gas cells, we calculate the profiles in a similar way to the calculation of the outflow fluxes in order to mitigate fluctuating spatial resolution; at each logarithmically spaced distance from the origin, we discretize the surface into HEALP/i.pc/x.pc pixels and determine the quantity at the distance as the mean across all the (equal area) pixels. The profiles are therefore in some sense volume weighted. For the cloud particles, we find all the particles within shells and determine the properties as a mass weighted average. We plot the simulations at three different times. In the top row we plot the density profile of the hydrodynamically resolved wind (dashed lines), as well as the mean density of cloud material (i.e. total mass of clouds in a volume, not to be confused with the subgrid internal density of the clouds themselves, solid lines). In the middle row, we plot pressure profiles of the resolved wind. In the bottom row, we plot radial velocities of the wind and clouds. This radial velocity is defined as the velocity parallel to the lines of constant polar angle in the truncated cone. 11  \nFor the clouds, we also indicate the 1 𝜎 scatter in the cloud velocities with a shaded band.  \nAs shown in Paper I, the wind without a cloud component has radial properties close to the expectations from theories of energy driven winds (e.g. Chevalier & Clegg 1985), particularly at early times. Away from the ISM, the density and pressure profiles decline approximately as 𝑟 -2 and 𝑟 -10 / 3 , respectively. The velocity increases rapidly within 10 kpc as thermal energy is converted to kinetic, before flattening out at a peak of ∼ 1300 kms -1 . As described above, later in the simulation, as the wind power decreases due to declining SFR, inflowing CGM gas begins to impact the wind. Thus, at late times, the density profile flattens at larger radii and this transition can also be seen in the pressure and velocity profiles.  \nAdding in cloud material at a mass loading of 0.1 has very little impact on the resolved wind. At 0.75 Gyr and 1.5 Gyr, the peak velocity of the wind is marginally reduced. The cloud density profile is approximately an order of magnitude lower than the resolved wind  \nradii as the gas flow is initially directed vertically out of the disc plane and therefore perpendicular to the vector to the origin for most of the disc. The two definitions of radial velocity converge after a few kpc.  \nFigure 7. Radial profiles of various quantities in the wind region, defined as a truncated cone expanding upwards from the disc with a cylindrical radius of 2.5 kpc (the disc scale length) in the disc mid-plane and opening angle of 10 degrees (this region is indicated in Fig. 2 and 3). The top row shows density profiles for hydro (solid) and cloud particle (dashed) components. Note that for the cloud particles, this is the average density of cloud mass in the region not the internal density of the clouds themselves. The middle row shows pressure profiles for the hydro component, though due to our assumption of pressure equilibrium this is also the internal pressure of the clouds. The bottom row shows velocity profiles of hydro (dashed) and cloud (solid) components. For the clouds, the shaded band indicates the 1 𝜎 scatter. Increasing the input cloud mass loading tends to reduce the normalisation of the velocity profiles for both components as the energy of the hot, fast wind is spread around more mass. This becomes less clear at late times due to the more messy velocity field (see Fig. 3).  \n<!-- image -->  \nat all radii. At late times, the uptick in density at large radii seen in the resolved wind is not reflected in the cloud population. This is because there are not significant inflows of cloud material from outside, unlike the resolved gas. At 0.75 Gyr, the cloud velocity increases gradually as a function of radius out to 𝑟 200 , indicating a consistent acceleration by the background wind. That said, the clouds are always much slower than the wind, reaching about half the velocity of the background by 𝑟 200 . At later times, the turnover in the radial profile of the velocity in the wind leads to a corresponding reduction in the acceleration of the cloud particles.  \nWhen the input mass loading of the cloud particles is increased to 1, there is a corresponding increase in the density profile of cloud material. At 0.75 Gyr, the gas and cloud particle density profiles  \nhave similar amplitudes and slopes throughout much of the halo (we examine this in more detail below). The density and pressure profiles of the gas remain similar to the equivalents in the 𝜂 cl0 𝑀 = 0 and 0.1 simulations. However, the wind does not reach such a high peak velocity, as momentum is transferred to the cloud particles. This also leads to the velocity profile of the cloud particles being shallower than the 𝜂 cl0 𝑀 = 0 . 1 case. Eventually, at later times, the slight reduction of the preventative feedback ability of the hot, fast wind leads to some re-establishment of inflows, which flattens the density profiles and turns over the velocity profiles.  \nFinally, we examine the 𝜂 cl0 𝑀 = 5 simulation. At 0.75 Gyr the density of both the gas and cloud material are higher than the other simulations. The acceleration of the wind is blunted even further,  \nFigure 8. Contribution to the total density profiles (i.e. mass fraction) in the wind, shown in Fig. 7, by pristine hydro (gas that has never been inside a cloud, solid line), processed hydro (gas that has been inside a cloud at some point in the past, dotted line) and clouds (dashed line). The columns show three different times, while the top row, middle and bottom rows show the 𝜂 cl0 M = 0 . 1, 1 and 5 simulations, respectively. As the input cloud mass loading is increased, the contribution of the cloud and processed components increases.  \n<!-- image -->  \nreaching a peak velocity less than half that of the simulation without cloud material. Likewise, the cloud particles are also accelerated more slowly. This reduction in wind velocities, along with the resulting re-establishment of inflows (as discussed above) leads to a significant enhancement of the density in both phases at late times. It can be seen in Fig. 3, however, that there are significant structures in both the density and velocity slices which makes a spatially averaged profiles an oversimplification.  \nIn Fig. 8, we plot the relative mass fraction of cloud and wind material in the outflows as a function of radius, using the same spatial cut as the radial profiles shown in Fig. 7. We consider mass in cloud particles, 'pristine' gas, that has never been inside a cloud particle, and 'processed' gas, that has been inside a cloud particle at some point in the past. To enable the quantification of the pristine and processed mass fractions, we make use of a passive scalar 'dye'. The mass fraction of the dye is initialised to zero in all gas at the beginning of the simulation. Mass that is transferred from a cloud particle to a gas cell is assigned a dye mass fraction of unity. The dye is advected by the hydro scheme. It can be seen that for the simulation with an input cloud mass loading of 0.1 (top row) the mass in the outflow region is dominated by pristine gas at all times and distances. For the 𝜂 cl0 𝑀 = 1 simulation (middle row) at 0.75 Gyr, once outside of the ISM, the cloud particles are the dominant mass component until ∼ 10 kpc, with a peak at 2 kpc, followed by a decline. This decline is matched by an increase as function of radius of the processed gas component. This occurs as the clouds are gradually shredded into the wind as they travel. At later times, the reduction in the absolute outflow rates and the increase in inflows leads to the pristine gas component becoming more dominant. The 𝜂 cl0 𝑀 = 5 case  \n(bottom row) has a similar, but more extreme behaviour. As discussed previously, the wind struggles to propel such a large mass of clouds to large distances. Thus, at 0.75 Gyr the clouds only dominate over the other mass components inside 0.1 𝑟 200 . However, the result is the wind is then overwhelmingly dominated by material stripped from clouds. This is true even at later times when inflows begin to become important.", '3.2.4 Cloud evolution': 'In the preceding figures, we have examined the properties of the entire population of clouds in the simulation. However, the Lagrangian nature of our approach to modelling the clouds enables us to follow the evolution of individual clouds over time, as we demonstrated in Fig. 1 for our simple wind tunnel tests. In Fig. 9, we plot some examples of evolutionary tracks for clouds extracted from the 𝜂 cl 𝑀 = 1 simulation. At 0.5 Gyr and 1.5 Gyr, we find all of the cloud particles created in the preceding 5 Myr. We plot their cloud mass, galactocentric distance and radial velocity as a function of time. For the 0.5 Gyr sample, we only plot half of the tracks (selected randomly) to avoid overcrowding the plot, while for the 1.5 Gyr sample we plot all tracks. 12 The colour used to plot the tracks reflect the initial cloud mass.  \nThere are general trends that are qualitatively consistent with the  \n]  \n1  \n[  \nFigure 9. Time evolution of cloud mass, galactocentric radius and radial velocity for a selection of cloud particles created in the 5 Myr preceding 0.5 and 1.5 Gyr in the 𝜂 cl0 𝑀 = 1 simulation. For the 0.5 Gyr selection (left) we plot half of the particles created in that time window (selected randomly) to avoid overcrowding the plot. For the 1.5 Gyr sample, we plot all tracks. Lines are coloured by their initial cloud mass. Clouds that are initially low mass are accelerated rapidly, but tend to be shredded. More initially massive clouds maintain or grow their mass, but are accelerated relatively slowly. Some clouds leave the hot wind flow and fall back towards the galaxy.  \n<!-- image -->  \nresults of our simple wind tunnel test. The clouds with the lowest initial masses are accelerated fastest but often experience a net mass loss at all times. More massive clouds experience a smaller acceleration but are more likely to grow or keep a stable mass. However, the spatially and temporally evolving background wind in these simulations produces additional complexities that are not present in a uniform wind tunnel. As demonstrated above, in the period spanned by our first set of tracks (0.5 - 1.5 Gyr), the wind structure is more well behaved than during the second period (1.5 - 2.5 Gyr), in the sense that there are fewer substructures and inflowing features. In the first set of tracks, it can be seen that clouds with initial masses below ∼ 100 M ⊙ are rapidly accelerated up to the wind velocity (see Fig. 7), at which point their velocity track suddenly flattens out. The scatter in the peak velocity reached is predominantly a result of the variation of the wind velocity with polar angle; the flow is faster closer to the  \ncentre of the wind bicone than at the edges. This can be seen in Fig. 2, while Paper I examines this in more detail. All of these clouds have been destroyed (in our implementation, they have lost more than 90% of their initial mass and are then fully recoupled to the gas) within 250 Myr. Many are destroyed on much shorter time-scales. Several of the cloud particles reach or exceed 𝑟 200 (97.9 kpc). The wind velocity drops suddenly in this region at these times, leading to a sudden deceleration of the clouds. A few of the low mass tracks show an earlier and more gentle deceleration; these are particles that are close to the edge of the wind and stray into regions of low wind velocity. All of these are destroyed before they fall out of the outflow region. For more massive initial clouds, we see similar trends, but the evolution is slower. Clouds experience lower accelerations, so they are not as closely coupled to the background gas. There is a large scatter in the tracks even with clouds of similar initial masses, seeded  \nFigure 10. Radial profiles at 0.75 Gyr in the 𝜂 cl0 𝑀 = 1 simulation within the wind region (defined as in Fig. 7) logarithmically binned by current cloud mass (left) and initial cloud mass (right). In the top row, we plot the relative contribution of clouds in a given bin to the total cloud mass density. The unmodified input cloud mass function ( 𝑑𝑁 / 𝑑𝑚 ∝ 𝑚 -2 ) would result in an equal contribution from the four logarithmic bins in the range 10 -10 5 M ⊙ ; this value is indicated for reference with the horizontal dashed line. The initial mass distribution of clouds is modified with increasing radii as low mass clouds shred and higher mass clouds grow. In the bottom row, we plot mass weighted radial velocity profiles for each cloud mass bin. Lower mass clouds are accelerated much faster than higher mass clouds, giving rise to a difference of around 200 km s -1 between the lowest and highest mass clouds at a given location.  \n<!-- image -->  \nby variations in the local properties of the background wind and then amplified with time. Some clouds grow at all times, others lose mass at later times.  \nFor the particles selected at 1.5 Gyr, the general trends remain, but more complex trajectories are evident. The properties of the background wind fluctuate significantly in space and time. The background velocity field experienced by a cloud can change suddenly as the opening angle of the wind changes as a function of radius due to collimation by the CGM. It can also change due to fluctuations originating upstream (caused, for example, by the changing SFR or substructures close to the disc) catching up with a cloud particle moving significantly slower than the wind. The lowest mass clouds are closely coupled to the background flow, but the lower accelerations experienced by more massive clouds mean that they are slow to react to these changes. This leads to an enhanced fraction of clouds falling out of the flow. Many clouds turn around and accelerate back towards the galaxy. These can be seen in Fig. 3, flowing down the edges of the outflow. Some fall back into the ISM, their density contrast drops below our threshold value of 10 and are recoupled. Others re-enter outflows as they near the galaxy and are accelerated away from the galaxy once again. The more massive a cloud, the more likely it is to survive by growing in response to the sudden increase of relative velocity as it enters the outflow. On the other hand, if it is too massive, it takes too long to arrest its fall before it reaches the ISM.  \nIn Fig. 10, we once again show radial profiles (in the outflow region, defined as in Fig. 7) of the cloud population from the 𝜂 cl0 𝑀 = 1 simulation at 0.75 Gyr, but binned by cloud mass. In the left panels  \nwe split the population by the current cloud mass while in the right panels we split by the initial cloud mass. We show the contribution of a given mass bin to the total cloud density as well as the mass weighted mean radial velocity for each mass bin. Our mass bins are logarithmically spaced. In combination with our adopted initial cloud mass function of d 𝑁 / d 𝑚 ∝ 𝑚 -2 , this means that the mass bins within our initial cloud mass range (10 -10 5 M ⊙ ) contribute equally to the total cloud mass at small distances from the galaxy. However, as lower mass clouds are shredded while larger mass clouds remain stable or grow (as shown previously), this means that the cloud mass at larger radii is dominated by the more massive clouds (or, indeed, the clouds that started off as the most massive). Turning to the radial velocity profiles, it can be seen that there is a wide spread in normalisations of the profiles. As previously shown, the more massive the cloud, the lower acceleration it experiences. This leads to a spread of around 200 km s -1 between the least and most massive bins. This is a major contributor to the cloud-to-cloud velocity dispersions within cells shown in Fig. 2 and 3, as more massive clouds are overtaken by less massive clouds.', '3.2.5 Ram pressure is not the dominant source of acceleration': 'Finally, in order to demonstrate that it is the accretion of high momentum material into the clouds that provides their high velocities, we perform two resimulations of the 𝜂 cl0 𝑀 = 1 case. We set the cloud mass growth rate, / 𝑚 grow , to zero. In the first simulation, we also set the drag coefficient, which appears in equation 24, to zero. This  \n0 . 75 Gyr  \n<!-- image -->  \n-100  \nFigure 11. Slices as in Fig. 2 for simulations with 𝜂 cl0 𝑀 = 1 at 0.75 Gyr. In the simulation shown in the left column, we disable both cloud growth and acceleration by ram pressure. In the middle column, we disable cloud growth but allow acceleration by ram pressure. The right column contains the fiducial physics simulation. Without a non-zero mass growth term (which also brings momentum into the cloud), clouds are not efficiently accelerated away from the galaxy and are rapidly destroyed.  \nmeans that the clouds will travel ballistically after being launched with the initial velocity of 100 km s -1 , only experiencing acceleration due to gravity. They still lose mass due to shredding, however. In the second re-simulation, we set the drag coefficient back to our fiducial choice of 1/2 to see if this source of acceleration is enough to propel clouds away from the galaxy. In Fig. 11, we show slices (equivalent to Fig. 2) of these simulations after 0.75 Gyr, alongside the fiducial full model run. It can be seen that only the fiducial simulation including cloud mass growth and associated momentum gain from the hot, fast wind leads to clouds being accelerated throughout the halo. The two re-simulations only show evidence of a small number of clouds in an x-shaped distribution (as seen in the slice). These are the only clouds that have been able to survive. As a reminder, in this work, we give the cloud particles an initial vertical kick of 100 km s -1 away from their birth site. However, this is added to the initial velocity of the star forming gas from which they were  \nborn, so the cloud trajectory is not exactly vertical. Without the mass growth and acceleration terms, cloud particles that end up in the hot outflow are rapidly destroyed before they can be accelerated by the drag term (in the simulation that includes it). Only those cloud particles launched exactly along the edge of the wind bicone (where the relative velocity is smallest) manage to survive, giving rise to the x-shaped distribution in the slice. Even then, these are only the most massive clouds. Lacking an efficient acceleration source, this very small amount of cloud mass quickly turns around and fall back to the galaxy. This confirms, for the full model, that the acceleration due to momentum-rich mass being accreted into the cloud (the first term of equation 5) dominates over the drag term and that our initial kick of 100 km s -1 is negligible compared to the velocities the clouds reach when accelerated by the hot wind.', '4.1 Model applications, interpretation and extensions': "In this work, we have demonstrated A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc in idealised noncosmological tests in order to highlight the main features of our scheme in a clean manner. In practice, the model is intended to be deployed in large-scale cosmological simulations. Bennett et al. in prep. will present the first application of the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc-H/o.pc/t.pc part of the scheme (introduced in Paper I) in cosmological simulations. The parts of the model presented in this article will be included in a cosmological setting in another future work. None the less, we stress that because the cloud-wind interaction terms depend solely on local properties (i.e. without reference to larger scale galaxy or halo properties), this aspect of the scheme is already fully deployable in cosmological simulations.  \nThe highly complex small-scale interactions of stellar feedback with the ISM that give rise to galactic winds are intentionally omitted in our model (via the use of initially hydro decoupled particles) because it is impossible to capture them properly at the resolution accessible in large-scale cosmological simulations. Attempting to launch galactic winds by direct injection of feedback energy into the ISM at these resolutions provides no predictive power. In fact, such an approach actively hinders the interpretation of results because the necessarily unphysical ISM behaviour and emergent wind properties are likely to be over-interpreted. Instead, in our scheme, the input mass and energy loading, along with the initial cloud mass function, encodes the properties of the wind as it leaves the ISM, distilling all of the unresolvable small-scale physics into an interpretable and numerically robust parametrisation.  \nIn this work, we used fixed values of mass and energy loadings for both wind phases for simplicity. The values were chosen to demonstrate various regimes of cloud-wind co-evolution. In practice, as discussed in Paper I, it is more appropriate for the loading factors to be adaptively selected during the simulation. There are essentially two approaches that can be adopted. Firstly, the loadings can be scaled with properties of the dark matter halo (e.g. local dark matter velocity dispersion, as in Oppenheimer & Davé 2008; Okamoto et al. 2010) or the galaxy (e.g. stellar mass, as in Davé et al. 2016). These scalings can be tuned to match some observational constraints on galaxy properties (e.g. Vogelsberger et al. 2013; Pillepich et al. 2018) or to capture a link between galaxy and wind properties seen in higher resolution simulations (e.g. Davé et al. 2016, 2019). Scaling with large-scale properties permits the placing of a wind with particular properties into a particular galaxy/halo, enabling controlled experiments to be performed or simply the recovery of certain galaxy properties (in the case of tuned scalings). An alternative approach is to scale the loadings with local ISM properties based on very high resolution simulations that are more able to resolve the relevant physics, providing more predictive power. One can go one step further and allow wind/cloud particles to have their initial velocity and temperature drawn from a distribution (see e.g. Kim et al. 2020b) rather than being single valued. This means that the winds (and their impact) are also more sensitive to other aspects of the entire galaxy formation model (the ISM model in particular), which can be an advantage or disadvantage depending on the aims of a given application. Therefore, the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc implementation is compatible with either method of scaling the loadings.  \nIn this work, we only tested one model of cloud-wind interaction. However, the underlying A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc methodology permits the FB22 model to be easily amended or entirely replaced. As described in Section 2.2, all that is needed is to provide a functional form of the mass, momentum, energy and metal exchange rates between  \nthe cloud and background wind. Physics not included in this work, for example the impact of thermal conduction and magnetic fields, can be included by changing how these various exchange terms arise from the current properties of the clouds and the local ambient wind. Likewise, a more detailed treatment of the internal phase structure of the cloud (the division between ionised, neutral and molecular gas) could be adopted with a suitable analytic prescription (e.g. Vijayan & Krumholz 2024). Processes leading to cloud splitting can be included by relaxing the assumption that the cloud number, 𝑁 cl , is an immutable property of the cloud particle. 13 Furthermore, any additional physics that impacts the properties of the wind as it leaves the ISM can be folded into the determination of the input mass and energy loadings. As an example, we can consider the potential impact of cosmic rays. Cosmic rays may increase the mass and energy loadings of cool gas leaving the ISM (see e.g. Rathjen et al. 2023; Armillotta et al. 2024). This effect can be included by altering the input mass and energy loadings, ideally by measuring them in a high resolution simulation. If one wished to include the effects of cosmic rays on the acceleration of clouds further away from the ISM, this could be achieved by modifying the cloud-wind interaction terms to account for this by providing an additional acceleration term (with an appropriate backreaction).  \nWhile we have primarily discussed the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc cloud model in the context of stellar feedback driven winds, the scheme can be used to model the behaviour of unresolvable cool clouds in other contexts. For example, the scheme could be used to model the evolution of cool material entrained in AGN outflows. Furthermore, the creation of clouds is not limited to ejection in winds. With the inclusion of suitable cloud creation criteria, A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc could be used to model the formation of clouds via thermal instability in the CGM/IGM and their subsequent precipitation (see e.g. Field 1965; McCourt et al. 2012; Sharma et al. 2012; Voit & Donahue 2015).  \nFinally, we emphasise that a key motivation behind the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc model is to enable us to make predictions of the observational signature of cool material in winds and the CGM and, importantly, its dependence on the assumed small-scale physics. Our model makes predictions for both the quantity, spatial distribution and small-scale properties of the clouds. We intend to translate these into mock observations of the clouds in absorption and emission. We leave a detailed discussion of this topic for a future work, barring a few brief remarks. As explained in Section 2, cloud particles represent the distribution of cloud material in a statistical sense; using one particle for one cloud would require a prohibitively expensive number of resolution elements. The particles are point masses with no extent, so one cannot compute an absorption spectra by simply ray-tracing through the gas distribution. Hummels et al. (2023) present a tool, C/l.pc/o.pc/u.pc/d.pcF/l.pc/e.pc/x.pc, that makes predictions for absorption-line signatures of cool material in the CGM distributed into complex structures of clouds. C/l.pc/o.pc/u.pc/d.pcF/l.pc/e.pc/x.pc places clouds in a CGM in a Monte Carlo fashion, taking as input various free parameters governing the distributions of cloud positions and properties. Once it has been populated with clouds, sightlines can be sent through the CGM and absorption spectra can be computed based on the clouds intersected. A similar approach could be used  \n13 The only constraint is that cloud particles can only contain a population of identical clouds because clouds of different masses follow different trajectories, precluding our approximation of the population as a point mass. For example, one could split all clouds in a particle into two equal clouds, resulting in a doubling of 𝑁 cl . However, any uneven splitting of clouds requires a splitting of the cloud particle itself into two (or more) particles that each contain homogeneous cloud populations. This functionality can be added in the future if needed.  \nto forward-model A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc outputs, but with the distribution of clouds constrained by the cloud particle locations. Thus we can make predictions for the observational signature of 'low ions' (e.g. C /i.pc/i.pc, C/i.pc/i.pc/i.pc, Si /i.pc/i.pc, Si /i.pc/i.pc/i.pc, Mg/i.pc/i.pc) and H /i.pc. In addition to absorption from the cool material within the cloud, we expect absorption from higher ions (e.g. Si /i.pc/v.pc, N /v.pc, O /v.pc/i.pc) within the intermediate temperature ( ∼ 10 5 K) gas of the TRML. This will require a model for the temperature and ionisation structure within the TRML (see e.g. Ji et al. 2019; Tan & Oh 2021; Chen et al. 2023) which is not a direct prediction of A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc. The TRML radiates energy away as material cools onto the cloud, primarily in Ly 𝛼 emission. As we know / 𝑚 grow for our cloud particles, we can relate this to a total luminosity (see e.g. Gronke et al. 2022) which is is a good proxy for Ly 𝛼 luminosity. As with the TRML absorption signature, predicting the emissivity of metal lines requires a more detailed model of the TRML (see e.g. Tan & Oh 2021; Chen et al. 2023).", '4.2 Comparison to other schemes': "Huang et al. (2020) present a subgrid model, named P/h.pcEW, in the meshless finite mass (MFM) code G/i.pc/z.pc/m.pc/o.pc (Hopkins 2015), where wind particles are treated as collections of clouds in an analogous manner to our scheme. Their model is inspired by the results of cloud crushing simulations with thermal conduction (Brüggen & Scannapieco 2016). In this framework, there is a competition between the tendency of thermal conduction to lead to the evaporation of clouds versus its ability to extend the lives of clouds by suppressing instabilities (e.g. Kelvin-Helmholtz) that would otherwise disrupt the cloud. The model does not include the effects of radiative cooling onto the cloud (as in the FB22) model, so there is no possibility for clouds to gain mass or momentum from the ambient wind and the only form of acceleration is ram pressure. Particles can only deposit mass into background gas, rather than needing to handle simultaneous deposition and accretion of material as in our case. Huang et al. (2022) present an application of the P/h.pcEW model in cosmological simulations. They launch all galactic wind material in the form of cloud particles. In other words, P/h.pcEW does not treat winds as multiphase, but as single phase outflows composed entirely of cool clouds. Unlike our model, there is no background hot, fast wind that can accelerate cool clouds once they leave the ISM. Cloud particles are therefore ejected from the ISM (with a brief hydro decoupled phase, as in our scheme) and are then decelerated by ram pressure as they impact the CGM, with the cloud evolution model determining how far they can travel before being completely destroyed.  \nThe impact of P/h.pcEW is compared to simulations where the wind particles deposit all of their mass into local gas once they meet the recoupling criteria. There are two primary impacts of the model. Firstly, the simulations with the P/h.pcEW model are in better agreement with low-redshift galactic stellar mass functions (GSMFs) for 𝑀 ★ < 10 11 M ⊙ due to increased wind recycling. This increase in recycling arises because the P/h.pcEW model permits a more gentle and spatially extended deposition of wind particle material, rather than a single point injection. However, we note that schemes similar to the nonP/h.pcEW model (e.g. Pillepich et al. 2018; Davé et al. 2019) can also achieve good agreement with the low mass end of the GSMF. The other major conclusion is that the P/h.pcEW model leads to a very different distribution of metals within CGMs. However, G/i.pc/z.pc/m.pc/o.pc is a completely Lagrangian scheme in the sense that there are no mass fluxes between resolution elements. This means that, in the absence of an additional metal mixing scheme (for discussions of G/i.pc/z.pc/m.pc/o.pcspecific implementations see e.g. Su et al. 2017; Rennehan 2021)  \nthere can be no metal fluxes between resolution elements. Huang et al. (2022) report that this drives the differences between the nonP/h.pcEW and P/h.pcEW simulations. In both cases, metals stay locked up in the gas particles into which they are deposited; the P/h.pcEW model spreads out this injection of metals over many more particles. It is therefore unclear whether this difference has physical significance or is compensating for a lack of metal mixing in the underlying scheme.  \nBoth A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc and P/h.pcEW model the cloud component of the fluid in a Lagrangian manner. However, an alternative approach is to use an Eulerian discretisation, known as a 'two fluid' or 'multifluid' approach. Weinberger & Hernquist (2023) introduce a generalised scheme along these lines in A/r.pc/e.pc/p.pc/o.pc, based on the stratified flow model of Chang & Liou (2007). Each finite volume resolution element represents a mixture of two different fluid phases, each with their own state vector. One can assume each fluid behaves in a similar manner to a single fluid, except with the addition of terms governing the exchange of mass, momentum and energy. In practice, this means solving three Riemann problems at each interface: the interactions between the first fluid on each side of the interface, between the second fluid and between the first and second fluid. Then, additional source/sink terms can be included to capture any additional interaction physics (such as a cloud-wind interaction) in an analogous manner to our implementation. However, in the case of cool clouds embedded in a wind, cloud - cloud collisions should not occur, due to their low volume filling fractions. For this reason, the two-fluid scheme of Butsky et al. (2024), implemented in the adaptive mesh refinement (AMR) code E/n.pc/z.pc/o.pc, models the cool component as being pressureless (although the clouds themselves have an internal pressure), which also has the advantage of simplifying the problem from the generalised case presented in Weinberger & Hernquist (2023). A benefit of multifluid schemes is that they do not require a second set of resolution elements (though the memory and computational requirements per element still increases with each additional fluid) and the fractional volume density of each component is well defined throughout the domain. Unlike our cloud particle based approach, multifluid schemes do not have to perform neighbour searches to associate the cool material with its neighbouring hot gas. That said, in A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc these searches only take up ∼ 2% of the total computational expense. Weinberger & Hernquist (2023) demonstrate a compelling application for the multifluid approach, replacing the SH03 eEoS model for the ISM, which assumes a two-phase medium in each cell with the same velocity, with two fluids that can have relative motion (but using the same source/sink terms). The general behaviour of the eEoS and multifluid simulations (idealised and cosmological) are the same on large scales, but the multifluid scheme allows for the production of a thin disc of cool material and a thicker disc of warm material, as might be intuitively expected.  \nHowever, in the context of cloud-wind interactions, multifluid schemes have several drawbacks which limit their applicability for modeling galactic outflows. Firstly, each phase in the resolution element can have its own independent velocity, but the velocity field for a given phase is single valued. As we demonstrated in this work, the velocity dispersion between cool clouds on our resolved scale is typically non-zero. Clouds with different properties experience differential acceleration even when exposed to the same background wind, leading to smaller clouds overtaking more massive clouds. Likewise, we demonstrate fountain flow behaviour with inflowing clouds passing outflowing clouds. In a multifluid scheme, regardless of whether the cool phase is modelled as a pressureless fluid or not, such behaviour is not possible. Intersecting flows of clouds 'collide' and result in a single, physically meaningless bulk velocity. This can only be avoided if the spatial resolution in the simulation is  \nhigh enough to resolve the inter-cloud separation, in which case the advantages of the multifluid scheme are lost. A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc occupies the opposite limit; the cloud particles are collisionless with respect to each other. The cloud velocity field can therefore be multivalued on any spatial scale. As currently implemented, our scheme is inappropriate for regimes where collisions between clouds are likely (e.g. the ISM), but this is not a concern for our intended applications.  \nSecondly, it is difficult to track the time evolution of properties of distinct cloud populations with a multifluid scheme. Only properties averaged across the volume element are known. For example, one can know the total mass of clouds in a cell but cloud scale properties (such as cloud masses or equivalently cloud radii) are not explicitly tracked. If a scheme could be constructed such that these properties were advected around with the fluid flow they would become poorly-defined when the cool fluid mixes with itself. For example, if there was a flux of clouds with a particular mass into a cell that contained clouds of another mass, the information of the two distinct cloud populations is necessarily lost. One could use an additional fluid for every possible configuration of cloud but this would rapidly become numerically intractable. Thus, cloud models have to make the assumption that either all clouds everywhere in the simulation domain are identical or that a cloud's properties can be instantly inferred from the local properties of a cell, discarding all prior evolution. One could imagine assuming that the cloud masses in a cell follow some statistical distribution (either assumed to be universal or determined from local properties), but this then falls foul of the single valued velocity problem described in the previous paragraph; clouds with different masses follow different trajectories, so cannot be represented by a single fluid. Additionally, we have shown that the cloud distribution at a given point in the wind (e.g. Fig. 10) emerges from the intersection of the trajectories of different clouds, launched at different times, that may have been stripped or grown on their way to that location, making an estimate of the distribution from purely local and instantaneous properties difficult.  \nFinally, we note a potential numerical issue with the multifluid scheme that could arise in our context. Weinberger & Hernquist (2023) report that with a fixed mesh their scheme is very diffusive, due to numerical mixing of the phases from advection errors. Once A/r.pc/e.pc/p.pc/o.pc is allowed to operate in its usual mode, with the mesh generating points moving with the local fluid velocity, these errors are reduced to a completely negligible level because the fluxes between the cells are significantly reduced. In the idealised test they use to demonstrate this behaviour, the two fluids have the same velocity. However, in the general case where there is a non-zero relative velocity between the two fluids, it is not possible to choose a mesh-generating point velocity such that both fluids are at rest with respect to the mesh. One could choose to move the mesh with the velocity of either fluid or the (weighted) bulk velocity of the combination, but it will lead to a relative mesh - fluid motion in either one or both of the fluids. When the relative velocity of both of the phases is small, the increase in diffusivity may be negligible. However, as we demonstrate in this work (e.g Fig. 7), the relative velocity between the wind and the clouds can be on the order of 600 -800 kms -1 (or equivalently, M ∼ 6 -8 in the cloud phase). Thus, any choice for the motion of the mesh would necessarily result in a large relative mesh - fluid velocity which may promote advection errors. In a static mesh code (such as AMR), the choice is not available.  \nIn addition to describing a general framework for a two fluid scheme, Butsky et al. (2024) also implement a cloud evolution model inspired by FB22, so it is instructive to make a direct comparison to our implementation. As described above, the choice to use a two fluid scheme removes the ability to take the history of a cloud in  \naccount when determining its future evolution. Butsky et al. (2024) set the cloud radius based on instantaneous and local properties of the cell, inspired by the mist model of McCourt et al. (2018). The implication of this is that as subgrid clouds gain mass they instantly fragment, increasing the number of clouds in the cell, in order to maintain a constant cloud radius. Likewise, as clouds lose mass, they must be instantaneously merging with other clouds, reducing the cloud number, in order to preserve the mandated cloud radius. While it may be appropriate for a fine mist, this property of the scheme is incompatible with the FB22 framework, where the current state of an individual cloud arises from the integration of its past evolution, an essentially Lagrangian behaviour.  \nIn another difference between our two implementations, Butsky et al. (2024) use a simplified form of the momentum transfer terms; compare our equation 6 to their equation 14. In their implementation, momentum transfer follows the net mass flow. In the FB22 model clouds gain mass via the TRML and lose mass through shredding simultaneously, meaning that the sign of the net momentum transfer is not constrained to have the same sign as the net mass transfer. As can be seen in our equation 7, the acceleration of the cloud due to mass transfer always reduces the relative velocity of the cloud and wind and its magnitude only depends on the mass growth rate, not the mass loss rate. This means that even if a cloud is losing mass overall (or its mass is not changing), it is still accelerated by accreted wind material as long as / 𝑚 grow ≠ 0. This can be seen in our Fig. 1 and 9. In the Butsky et al. (2024) scheme, only clouds that have a net mass growth ( / 𝑚 = / 𝑚 grow - / 𝑚 loss > 0) can be accelerated by the accretion of hot wind material and, even then, this acceleration is underestimated. A related difference between our schemes arises in the form given for the energy transfer from hot phase to cool phase (their equation 18). Their prescription allows energy transfer based on the magnitude and sign of the net mass transfer, rather than considering the mass growth and loss rates independently. The net energy transfer should not depend purely on the net mass transfer rate, but on the growth and loss rates independently (see our equation 25). Additionally, the energy transfer determined via their equation 18 enforces the specific energy of the mass being transferred to always be that of the hot phase, regardless of the sign of the net mass transfer. The consequence is that as cool gas is mixed into the hot phase, the specific energy of the resulting mixed gas is the same as the hot phase before mixing. Thus, the cool gas has been heated to the temperature of the hot phase, violating energy conservation.", '5 CONCLUSION': "A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is a novel scheme implemented within the A/r.pc/e.pc/p.pc/o.pc code that allows the modelling of multiphase stellar feedback-driven galactic winds within coarse resolution cosmological hydrodynamic simulations of galaxy formation. In Paper I, we demonstrated aspects of the scheme that allow high specific energy (i.e. hot and fast) winds to be accurately modelled. This is particularly important to properly capture the operation of such winds as a form of preventative feedback. In this work, we have presented the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc cloud particle treatment, which permits the inclusion of a population of cool clouds embedded in the hot, fast wind. Resolving the relevant interactions between cool clouds and the ambient medium in a cosmological galaxy formation simulation is intractable because of the vast dynamic range in spatial scales that must be captured. However, their inclusion is necessary in order to model a truly multiphase wind, while preceding analytic work has demonstrated the importance of their impact on high specific energy outflows.  \nWemodel cool clouds using collisionless N-body particles (which we term cloud particles) that can move relative to the gas resolved by A/r.pc/e.pc/p.pc/o.pc's hydrodynamic scheme. A cloud particle represents an unresolved population of cool clouds. In addition to feeling gravity, the cloud particle exchanges mass, momentum, energy, and metals with the ambient gas with which it is co-located. Models of cloud evolution and cloud-wind interactions, derived from high resolution simulations of 'cloud crushing' and turbulent radiative mixing layers, can be distilled into the various exchange terms, as well as into the subgrid evolution of the cloud population. A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc is agnostic as to the choice of cloud-wind interaction model, granting substantial flexibility. In addition to describing the underlying scheme, in this workwepresentedanimplementationoftheFB22cloud-windmodel into the A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc framework.  \nWe demonstrated the scheme using idealised non-cosmological simulations carried out at a coarse resolution appropriate for a large volume cosmological simulation. The setup duplicated that used in Paper I, featuring a 𝑀 200 = 10 11 M ⊙ system with a cooling flow CGM. We included a high specific energy wind, then performed simulations with varying input mass loadings for the cool cloud phase, enabling us to illustrate various characteristic regimes. The main findings are as follows:  \n- (i) As expected from preceding analytic work, the FB22 cloudwind interaction model enables the efficient acceleration of clouds embedded in the hot wind. All clouds are accelerated, but their fate depends on their initial mass. Initially low mass clouds are rapidly accelerated, but are relatively short lived. Higher mass clouds take longer to be accelerated up to the background wind velocity, but are more likely to survive.\n- (ii) The strong dependence of cloud trajectories and mass evolution as a function of initial cloud mass and current local wind conditions gives rise to a complex distribution of cloud masses throughout the wind. Importantly, this means that the distribution cannot be captured by instantaneous local properties, making it impossible to capture with a simpler statistical model.  \n(iii) As shown in Paper I, the properties of the background wind vary in both space and time. There is a variation in the wind properties with polar angle, the wind is collimated by the CGM (which changes the opening angle of the wind with height) and changes of the SFR of the galaxy over time result in fluctuations in wind properties. Clouds take time to respond to changes in the local properties of the background wind, with less massive clouds having a stronger coupling to the flow than more massive clouds. This leads to a complex range of cloud trajectories, with some clouds moving out of the wind and falling back towards the galaxy. Some of these return to the ISM, but others can be re-accelerated before they reach the galaxy. This adds further complexity to the spatial distribution of clouds in a way that cannot be captured in a simple 1D model.  \n(iv) The emergent mass outflow rate of cool clouds at some distance from the galaxy is a non-monotonic function of the input cool cloud mass loading factor because of the back-reaction on the high specific energy wind. When the input cool cloud mass loading is negligible, the impact on the wind is similarly negligible, allowing it to continue to efficiently regulate the galaxy's SFR by suppressing CGMinflows (i.e. preventative feedback). When we used a relatively high input cloud mass loading, this lead to an initially high flux of cool clouds and a noticeable drop in the SFR because of the increase in the amount of mass being removed from the ISM (i.e. ejective feedback). However, increasing the mass of clouds in the high specific energy wind necessarily results in an increase in the transfer of energy from the wind to the clouds. This resulted in an eventual  \n'poisoning' of the wind, the resumption of strong inflows from the CGM and a return to higher SFRs. We demonstrated a simulation with an input cloud mass loading between those two extremes. This was able to fill a substantial portion of the halo with cool clouds carried by the wind without overly impacting the ability of the hot wind to perform preventative feedback. The exact boundaries of these different regimes will depend strongly on properties of the galaxy, CGM, halo, and mass and energy loading of the high specific energy wind.  \n- (v) Several phenomena illustrated in our simulations indicate that two-fluid/multifluid schemes, an alternative approach to that adopted by A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc, are not as well suited for modelling cool clouds entrained in a wind. While allowing a relative velocity between the hot and cool phases, these schemes enforce a single valued velocity field for the cool material. We have demonstrated that on the scales resolvable in this class of simulation the cool cloud velocity field is multivalued because of the differential acceleration of clouds with different masses and the ubiquitous presence of intersecting cloud trajectories (e.g. in fountain flows). Additionally, as noted in point (ii), the distribution of cloud properties at a given location results from an integration of their past trajectory and mass evolution in a manner that cannot be derived from instantaneous or local properties. It is difficult to conceive of a method of tracking this evolution with an Eulerian approach as used in two-fluid/multifluid schemes.  \nIn this work, we demonstrated A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc with a particular model for cloud-wind interactions (FB22), but the scheme is agnostic to this choice. We can therefore use this approach to study the large-scale consequences of other theories of cloud evolution, which may include additional physics such as thermal conduction and magnetic fields. We adopted constant mass and energy loading factors to provide a clean numerical experiment, but in future work these will be varied according to halo, galaxy and/or ISM properties. We performed idealised non-cosmological simulations in this work, but the scheme has been specifically designed to work in cosmological volume simulations; Bennett et al. in prep. will demonstrate some first cosmological applications of A/r.pc/k.pc/e.pc/n.pc/s.pc/t.pc/o.pc/n.pc/e.pc.", 'ACKNOWLEDGEMENTS': "We are grateful to Rüdiger Pakmor for helpful comments. This work was supported by the Simons Collaboration on 'Learning the Universe.' GLB acknowledges support from the NSF (AST-2108470, XSEDE grant MCA06N030), NASA TCAN award 80NSSC21K1053, and the Simons Foundation (grant 822237). CGK and ECO acknowledge grant 10013948 from the Simons Foundation to Princeton University, to support the Learning the Universe Collaboration. Computations were performed on the HPC systems Raven and Freya at the Max Planck Computing and Data Facility (MPCDF). The following open source software packages were used in this work: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), Matplotlib (Hunter 2007), nanoflann (Blanco & Rai 2014), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020).", 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding author.', 'REFERENCES': 'This paper has been typeset from a T E X/L A T E X file prepared by the author.'}