Effective viscosity from cloud-cloud collisions in three-dimensional global SPH simulations

Analytic estimates of the viscous time-scale due to cloud-cloud collisions have been as high as thousands of Gyr. Consequently, cloud collisions are widely ignored as a source of viscosity in galactic disks. However, capturing the hydrodynamics of discs in simple analytic models is a challenge, both because of the wide dynamic range and importance of 2D and 3D effects. To test the validity of analytic models we present estimates for the viscous time-scale that are measured from three dimensional SPH simulations of disc formation and evolution. We have deliberately removed uncertainties associated with star-formation and feedback thereby enabling us to place lower bounds on the time-scale for this process. We also contrast collapse simulations with results from simulations of initially stable discs and examine the impact of numerical parameters and assumptions on our work, to constrain possible systematics in our estimates. We find that cloud-collision viscous time-scales are in the range of 0.6-16 Gyr, considerably shorter than previously estimated. This large discrepency can be understood in terms of how the efficiency of collisions is included in the analytical estimates. We find that the viscous time-scale only depends weakly on the number of clouds formed, and so while the viscous time-scale will increase with increasing resolution, this effect is too weak to alter our conclusions.


INTRODUCTION
Arguably, the most successful model for the formation of disc galaxies is the ΛCDM model, in which galaxies are formed from the dissipational collapse of baryonic gas within a dark matter halo (White & Rees 1978;Blumenthal et al. 1984;Davis et al. 1985;White & Frenk 1991;Kauffmann et al. 1993;Cole et al. 1994;Baugh 2006;Benson 2010). While the physical viscosity of the baryonic gas is not anticipated to have a strong influence on gas evolution except in magnetized or hot environments such as a galaxy cluster (Sijacki & Springel 2006), effective kinematic viscosities could in principle impact disc evolution. Simulations by Lin & Pringle (1987) with a viscous time-scale close to the star-formation time-scale showed that viscous evolution with infall can reproduce the ubiquitous exponential density profile from a range of initial conditions. In this work the viscosity was assumed to be caused by large-scale turbulent motions dissipating kinetic energy and transporting angular momentum.
Feedback from supernovae can be a source of viscosity E-mail: williams@ap.smu.ca by feeding this turbulence (Vollmer & Beckert 2003). Additionally, the self-gravity of the gaseous disc can provide an effective viscosity (Vollmer & Beckert 2002). This can take the form of large-scale instabilities (Rafikov 2009;Gammie 2001), or of interactions between giant molecular clouds (Vollmer & Beckert 2002). These cloud interactions potentially generate viscosity through two different mechanisms. Firstly, gravitational scattering can increase the velocity dispersion of the cloud population, converting orbital energy into large-scale turbulence (Gammie et al. 1991;Fukunaga & Tosa 1989;Agertz et al. 2009). Secondly, during inelastic collisions between clouds, shocks convert orbital energy into turbulence and heat within the colliding clouds (e.g. Gittins et al. 2003;Kitsionas & Whitworth 2007;Anathpindika 2009). Radiative processes contribute to the dissipation of kinetic energy during these collisions, and are also important for dissipating turbulent energy that has cascaded into thermal energy. These processes are significant even in the absence of star-formation: the observations compiled by Dib et al. (2006) show that the velocity dispersion of HI gas does not strongly depend on the star-formation rate below a certain threshold, and the AMR simulations of Agertz et al. (2009) suggest that a 'baseline' turbulence is caused by interactions between clouds, and that this is only supplemented by supernova feedback at high star-formation rates. It has been argued (Vollmer & Beckert 2002;Bell 2002, hereafter B02) that cloud-collisions are not an efficient source of viscosity. In particular, in B02 the time-scales for viscosity due to cloud collisions are estimated to be on the order of tν ∼ 1000 Gyr in most local spiral galaxies, although the time-scales might be considerably lower in earlier gasrich galaxies or in galaxies where the velocity distribution of GMCs has been stirred up by some mechanism (such as galaxy interaction e.g. Hernquist & Mihos 1995). Vollmer & Beckert (2002) argues that because molecular clouds evaporate at an age of ∼ 10 7 yr, and this is less than the time between collisions (∼ 10 8 yr), cloud collisions are very rare. However, cloud formation times, assuming that collapse and formation of H2 are the dominant factors in forming a cloud, appear to equally short (Glover & Mac Low 2007). This leads to a scenario in which the number density of clouds is roughly constant, although the short life-time may affect the velocity dispersions of molecular clouds as they have less time to build up a large deviation from circular velocity through scattering events with other clouds. In this steady state, the effective collision time-scale should remain similar.
It has also been argued that physical collisions between clouds have a smaller effect than gravitational scattering (Jog & Ostriker 1988), particularly as if magnetic fields are taken into account, cloud interaction cross-sections may be underestimated (Ozernoy et al. 1998). On the other hand, Das & Jog (1996) modelled a system of cloud particles, finding that cloud collisions rather than local gravitational interactions (scattering events) dominate the mass distribution and velocity dispersion of molecular clouds, suggesting that cloud collisions may indeed be important. However, as far as we are aware, the effective viscosity of direct cloud-cloud collisions has not yet been examined in global three dimensional numerical hydrodynamic models.
Most simulations of cloud formation and the associated disc dynamics have been performed in two dimensions and/or on a small scale using shearing-box studies (e.g. Kim & Ostriker 2007). However, increased computing power and the availability of locally adaptive algorithms have recently enabled galaxy-scale simulations with sufficiently high resolution to resolve cloud-formation in discs. Numerical experiments have been performed using both AMR (Tasker & Tan 2009;Tasker 2011;Agertz et al. 2009) and SPH (Robertson & Kravtsov 2008) with resolutions as fine as 6 pc. The nontrivial cooling processes and chemistry make these simulations a significant technical challenge. Agertz et al. (2009) and Tasker & Tan (2009) ran suites of high resolution AMR simulations of Milky-Way and M33-like disc galaxies, and reported on the properties of the clouds generated by their models, including cloud-cloud velocity dispersion. However, neither study has provided an estimate of the viscous timescale due to cloud-cloud collisions. Furthermore, the discs of Tasker & Tan (2009) are much more stable than the Milky Way, with a density distribution chosen to give a constant value of the Toomre Q parameter (Toomre 1964), and a static dark matter and stellar component, which may inhibit some of the instabilities important to cloud formation.
In this paper we revisit the calculations of B02 with full three dimensional SPH models. This is not entirely triv-ial since there is no universally agreed upon cloud finding process. However, the use of a particle method enables the Friends-of-Friends (Davis et al. 1985) group finding methodology and we adapt that to our simulations. Hence given our cloud population our primary goal is to see whether the analytical calculations are supported, and if not what the implications are. It is important to note that the results of such simulations could highlight non-physical evolution in numerical schemes with artificial viscosities, of which SPH is a notable example (Valdarnini 2011). We also investigate the issue of numerical artefacts in our calculated results. This is a key issue since structure formed within simulations starting from smooth initial conditions is inevitably the result of amplification of noise in the initial conditions. While a full calculation in the cosmological context (e.g. Katz 1992;Thacker & Couchman 2001;Governato et al. 2007;Stinson et al. 2010;Scannapieco et al. 2009;Brook et al. 2008) is beyond the scope of this paper, primarily due to resolution limitations, we instead consider two classes of isolated models. We examine an equilibrium system with similar parameters to the Milky Way consisting of a gas disc, a stellar disc and bulge, and a dark matter halo. Here the gas disc is stabilized by the other components which dominate the system's mass. We also consider the dissipational collapse problem that has been used extensively elsewhere (e.g. Gott & Thuan 1976;Carlberg 1984;Katz & Gunn 1991;Brook et al. 2004;Kaufmann et al. 2006). By contrast with the Milky Way model, this collapse produces a very unstable disc, and so we investigate both high-stability gas-poor systems and low-stability gas-rich systems. These models include hydrodynamics, gravitational interactions, and cooling with a dynamic temperature floor. By removing the numerous unknowns associated with star-formation and feedback (as discussed in numerous places e.g. Thacker & Couchman 2000;Ceverino & Klypin 2009;Christensen et al. 2010) we hope to isolate the impact of cloud-cloud interactions and place lower bounds on the viscous time-scale.
The layout of the paper is as follows: in section 2 we outline the details of our simulation code. We also discuss the initial conditions, our cloud identifying approach and also the underlying theory of the effective viscosity. Results are presented in section 3 followed by a brief conclusion.

Simulation Code
We model the dark matter, stars and gas using a specially adapted version of the OpenMP n-body AP 3 M (Couchman 1991) SPH (Monaghan 1992) code HYDRA (Thacker & Couchman 2006). Our modifications are: (i) The cooling curve has been extended to down to 10 K using the cooling function (Λ) of Wada & Norman (2001), although we set our fiducial temperature floor to 300K to make our results more comparable with Tasker & Tan (2009), except in cases where we investigate the effect of a lower floor. The earlier cooling curve of Sutherland & Dopita (1993) is retained for T > 10 4 K. The combined cooling curve is plotted in Fig. 1.
(iii) The parallelisation algorithm has been altered so that during the particle-particle gravity and SPH calculation regions containing a large number of particles are split over all processors, instead of each processor receiving a single region. This greatly improves load-balance in simulations containing many dense clumps of particles.

Dynamic Temperature Floor
We use a method similar to Robertson & Kravtsov (2008) to ensure that the Jeans mass is resolved in our simulations. This is to satisfy the Truelove et al. (1997) criterion and avoid artificial fragmentation -crucial in simulations of cloud formation. The method is in the form of a dynamic pressure floor. The Jeans mass (Jeans 1902) is defined as Bate & Burkert (1997) noted each particle should satisfy 2N neigh mgas < mJeans (where N neigh is the number of SPH neighbours for the particle and mgas is the gas particle mass) to avoid artifical fragmentation. Defining the local ratio of the Jeans mass to the SPH kernel mass as hJeans, then this requirement is equivalent to the statement that with NJeans being the required factor by which the Jeans mass must be resolved, which in the Bate and Burkert case is set to 2. In an ideal gas cs ∝ √ u, so we can fulfill this criterion by applying whenever hJeans < NJeans. We found spurious (sub-resolution) string-like structures forming within clouds for low values of NJeans, and found that NJeans = 50 removed these structures and resulted in a more homogeneous interior for clouds.

Milky Way Model
We produce our Milky Way model using the GALACTICS package (Kuijken & Dubinski 1995;Widrow & Dubinski 2005;Widrow et al. 2008) with the parameters in Table 2 of Widrow et al. (2008). Through an iterative process, this package produces an equilibrium system consisting of an exponential stellar disc, a stellar bulge and a dark matter halo. The disc is exponential radially, follows sech 2 vertically, and has a radial dispersion profile of where we set Rσ = R d for simplicity. We generate the gas disc by copying the disc star particle positions and flipping the coordinates across the x = y plane to prevent particles having coincident positions. Bulge particles are not copied. The masses of the gas and star particles are scaled to give the appropriate mass ratio. The gas disc is given a dispersionless velocity profile output by GalactICs and is initially isothermal at 10 4 K. The disc scale length is 2.81 kpc, truncated at 30 kpc with a truncation scale-length of 0.1 kpc. The scale height is initially 0.36 kpc, and the total disc mass is 5.2 × 10 10 M . The halo density profile is where a h is the halo scale parameter, r h is the cutoff radius, δr h is the scale length for the truncation, γ is the 'cuspiness' parameter (equal to unity for an NFW profile), and σ h is a velocity parameter that sets the mass of the halo. Halo parameters are given in Table 1. We also ran a test with a static analytic halo potential, to explore if the discretization of the halo has any effect on cloud formation.
The bulge density profile is where p = 1 − 0.6097/n + 0.05/n 2 gives a Sérsic profile with n the Sérsic index. Re is the radial scale parameter, and in GALACTICS ρ b is parametrized by the velocity parameter We set these parameters to n = 1.31, σ b = 272 km s −1 , and Re = 0.64 kpc.
We have named our fiducial run LowSoftMW. To test the effects of change the resolution, softening length, temperature floor, gas mass fraction, and artificial viscosity we investigate a total of ten different runs, summarized in Table 2. Both MidSoftMW and HighSoftMW have higher gravitational softening lengths; MedGasMW and HighGasMW have higher gas mass fractions; LowFloorMW has a lower temperature floor; LowViscMW has lower artificial viscosity parameters (α, β); and LowResMW has a lower resolution. In addition, as a convergence check we ran a higher resolution simulation (HighResFlatMW) with a total of 3.5 × 10 6 particles and a softening length of 45 pc, although we do not consider this our fiducial run as the simulation time did not reach a full Gyr. We found when running a simulation of this high resolution with identical initial conditions to LowSoftMW that the disc was dominated by a strong ring-shaped shock propagating outwards. This is caused by a combination of the rapid vertical collapse of the disc as it initially cools, and the rotation curve not being quite precise enough as GALACTICS is intended for collisionless mechanics and does not take into account gravitational softening or the pressure gradient of the gaseous disc. At the lower resolutions this shock is not well captured, and the disc quickly returns to equilibrium, so this is only a problem at our highest resolution. To prevent the shock becoming a problem it is necessary to start the simulation from an initially flattened state akin to the later evolution of the cooled disks. We therefore flattened the gas disc to a similar scale height as the cooled disks, which is a factor of 10 smaller. Circular velocities (vcirc) were then set up using radial accelerations (a rad ) generated from a single iteration of the HYDRA code, and explicitly setting a rad = v 2 circ /R for each gas particle, where R is the radial coordinate of the particle. We also performed a simulation (FlatMW) with these initial conditions but at our fiducial (moderate) resolution, for a fair comparison of the effects of resolution.

Monolithic Collapse Model
This model consists of a spherically symmetric distribution of gas within an equilibrium NFW dark matter halo. We generate the halo using GALACTICS according to the parameters in Table 1, giving a halo with M = 1.1 × 10 12 M .
For the gas we use the 'high-entropy' profile of Kaufmann et al. (2009), which was produced from equation 1 of Kazantzidis et al. (2004), setting c = 1, α = 1, β = 3, and γ = 0. Kaufmann noted that a gas density profile that is shallower than the NFW profile (as expected in models with pre-heating feedback e.g. Mo & Mao 2002) produces an angular momentum distribution in the final object that better fits observations. In this model, the gas collapses into clumps which combine to form an unstable disc.
As in Kaufmann et al. (2007), the initial temperature profile is calculated to provide hydrostatic equilibrium according to where µ is the mean molecular weight of the gas (taken as its primordial value, µ ≈ 0.59mH ), kB is the Boltzmann constant, ρG is the initial gas density, and Mtot(r) is the total mass (gas and dark matter) within a sphere of radius r. We give the gas a flat velocity profile. The positions of the gas particles in our initial conditions are simply the generated positions of the dark matter particles flipped as in section 2.2.1.
To set up a rotating halo, GALACTICS swaps a fraction of the dark matter particles' velocities over the radial axis to increase the number of particles rotating in the same direction. We assume the gas and the dark matter have the same specific rotational momentum, i.e.
so that the spin parameter (Binney & Tremaine 2008) of the gas is equal to the spin parameter of the halo, We used a spin parameter of λG = 0.038, close to the median value observed in simulations (Bullock et al. 2001;Barnes & Efstathiou 1987). After each gas/DM halo is produced, it is evolved for 0.5 Gyr with cooling switched off to ensure the ICs are stable. Our first model (HighSoftC)  We also investigated models with lower softening lengths and temperature floors to see if smaller clouds were resolved. A low resolution run was performed as a convergence check, and finally we performed a model with a low gas fraction to see the effect of increasing the disc's stability. These models are summarised in Table 3.

Identification algorithm
We identified clouds using a two-step process. A density threshold was applied and then particles above this threshold were linked into groups using the friend-of-friends algorithm. Using a density threshold partially avoids the notorious 'string of pearls' effect that that may lead to spurious filamentary structures or the merging of many smaller structures into one larger one. We set the density threshold at a density of n = 7 M pc −3 , and used a linking length of 50 pc. The high density threshold ensures that only dense cloud-like objects are selected, while the linking length is close to the size of the softening ensuring small fluctuations below this threshold can be skipped over. We set the minimum cloud size to 30 particles, giving minimum cloud mass of 1.6×10 5 -8.5×10 6 M , depending on resolution. We found these parameters largely selected dense, cool clumps while excluding other objects such as filaments from cloud encounters. Because our lower limit for cloud size is a number of particles and not a mass, we include clouds of increasingly small mass as we increase the resolution of our simulation. This may create a resolution dependence until individual molecular clouds are resolved -a level that we do not achieve in this work, although as noted in section 3.1.2 the major axisymmetric instabilities appear to be resolved.
Clouds are tracked from output to output by examining the particles resident in each cloud. If the cloud A at time ti contains at least half of the particles contained by cloud B at the time of the following output ti+1, then A is a parent  Table 2. Summary of Milky Way runs. l sof t is the minimum softening length, T f loor is the temperature floor, n * , ng, and n DM are the numbers of star, gas and dark matter particles, mg/m * is the gas/star mass ratio for the disc, t end is the total simulation time, h disc is the scale height of the disc, and α and β are the artificial viscosity parameters. Name HighSoftC 514 3 × 10 4 5 × 10 5 1 × 10 5 0.148 4.5 MidSoftC 200 3 × 10 4 5 × 10 5 1 × 10 5 0.148 3.9 LowSoftC 60 3 × 10 4 5 × 10 5 1 × 10 5 0.148 3.3 LowSoftFloorC 60 300 5 × 10 5 1 × 10 5 0.148 3.7 LowResC 60 300 1 × 10 5 1 × 10 5 0.148 4.6 LowMassC 512 3 × 10 4 5 × 10 5 1 × 10 5 0.030 7.8 Table 3. Summary of collapse runs. l sof t is the minimum softening length, T f loor is the temperature floor, ng and n DM are the numbers of gas and dark matter particles, mg/m DM is the gas/dark matter mass ratio, and t end is the total simulation time.
of B. If B contains at least half of the particles contained by cloud A, then B is a child of A. If B has several parents, then a merger has occurred. If A has several children, then a separation has occurred. If A is the only parent of B, and B is the only child of A, then B is identified as the same cloud as A. This categorization allows for multiple parents to join in a merger and it is also possible for a parent to split into into multiple children. During simulations we observed that mergers can be complex with clouds merging and separating several times before settling into a single cloud, or in some cases while no longer interacting. This means our statistics are perhaps better thought of as recording 'interaction' rates (including 'self-interaction') rather than cloud collision rates.

Treatment of cloud energy and interactions
To quantify the energy loss due to interactions, we compare the kinetic and potential energies of clouds in sequential outputs across a separation or merger event. For the 'combined' stage of the interaction (the earlier output for a separation, the later output for a merger) we calculate the centre of mass kinetic energy, where P is the set of all particles 'involved' in the interaction (defined below), and v i and mi are the velocity and mass of particle i. This is compared with the sum of the centre of mass kinetic energies of the clouds during the 'separated' stage of the interaction, where pj is the set of all particles in cloud j and C is the set of all clouds involved in the interaction during the separated stage.
The merged or unseparated cloud does not typically contain all of the particles from the clouds that formed it or that separated from it; there are always a number of particles that are expelled during the interaction, while other particles may accrete on to the clouds during the interaction. These particles carry kinetic energy, and so to ensure that we are measuring a real loss of energy from the system and not just an apparent energy loss from particles leaving the cloud phase, we define P by The final step in the energy budget calculation is to ensure that the energy change due to clouds moving in the potential well of the dark matter and baryons is accounted for. We take out the effects of these gravitational interactions by calculating the potential between P and all other particles during both dumps, and subtracting the difference from the kinetic energy that was calculated.

Energy analysis
From the cloud energy budget we can obtain an estimate for the total time-scale for dissipation of kinetic energy from cloud-cloud interactions. By analogy with the star-formation time-scale, typically defined as tSF R = Σgas/(dΣ * /dt), we define the viscous time-scale due to cloud-cloud collisions to be the dissipative time-scale, where K is the total rotational kinetic energy of the gas, and −dKC /dt is the rate at which this kinetic energy is dissipated due to collisions. This dissipation rate can be written as where C is the interaction rate (determined by counting the number of interactions that occur within a time period) and ∆K col is the energy lost per interaction. In practice, we average over n col interactions so that where ∆t is the time period that the n col interactions occurred over (and hence C = n col ∆t ), ∆Ki is the kinetic energy lost in a particular interaction i and K(ti) is the total kinetic energy in gas at the time of that interaction.
It is important that we connect this method of measuring the dissipative time-scale in our models with definitions used elsewhere. It is commonly argued (e.g. Bell 2002) that the form of the viscous time-scale is where R is the radial coordinate and ν the (effective) viscosity. To see how this form arises in our measurements, consider the following argument: If we neglect radial velocity, then the kinetic energy per unit volume of a component of fluid in a rotating disc is k = ρ(RΩ) 2 /2. We can convert the rate of viscous dissipation for a generic fluid (Φ, the energy lost per unit volume per unit time) from Mihalas & Weibel Mihalas (1984) into cylindrical coordinates and again assume angular velocity dominates, simplifying it to: where the prime indicates a radial derivative. We can substitute these values into our definition for tν col because Φ = −dkC /dt, so If we then take a power law for rotation Ω ∝ R −α then which agrees with R 2 /ν within a factor of 1/2α 2 . For a flat rotation curve, α = 1 and this factor is merely 1/2 -hence the dissipative time-scale is of the order of the traditional viscous time-scale. Note, Lin & Pringle (1987) give a different prefactor -(2 − α)/(α). However, these values all agree within an order of magnitude, provided α is not extremely large or small. Although our viscous time-scales are calculated over the whole disc to ensure sufficient numbers of interactions are measured, and the analytical R 2 /ν is a local value at a specific radius, we should not expect this to have an effect beyond an order of magnitude, assuming analytical viscous time-scales have been calculated at a representative radius.

General Evolution
The evolution of all models excluding HighSoftMW are similar 1 . In these models the gas disc is initially close to equilibrium. However, the gas rapidly cools and becomes unstable, collapsing vertically (except in FlatMW and HighRes-FlatMW, which are produced from already-collapsed initial conditions), and forming spiral instabilities which fragment into large number of small (m ∼ 10 6 -10 7 M , R ∼ 100 pc) clouds.
After this epoch of rapid cloud formation, the clouds merge and continue to accrete material. The number of clouds drops, as illustrated in Fig. 2, while the total mass within clouds continues to increase until both reach a less dramatic stage from around 0.8-1.0 Gyr, where the number of clouds decays only gradually as the mass within clouds gradually increases. A face-on view of the evolution of Low-SoftMW is shown in Fig. 3, and a snapshot of HighRes-FlatMW is shown in Fig. 4. In HighSoftMW cloud collapse was quenched by the high softening length, and instead the disc was dominated by large scale instabilities (Fig. 5). The higher gas mass in HighGasMW and MedGasMW reduced the hydrodynamic time-step and so these simulations could only be run for ∼ 0.45 Gyr, while the increased computational load of the high resolution run HighResFlatMW also made a full simulation of 1.0 Gyr unfeasible, and so this simulation was evolved for ∼ 0.3 Gyr.
The gas disc separates into two phases: diffuse gas which retains a moderate temperature (∼ 10 3 to ∼ 10 4 K) though shock heating and a low cooling time, and dense gas whose temperature is tightly controlled by the Robertson-Kravstov dynamic temperature floor. It should be noted that while our models lack direct stellar feedback, the dynamic floor can heat the dense gas to temperatures as high as 3×10 4 K. This temperature is equivalent to a sound speed of ∼ 26 km s −1 , which is on the order of the velocity dispersion generated by various feedback mechanisms (Thacker & Couchman 2000;Governato et al. 2007;Ostriker & Shetty 2011). Hence while we expect implementing feedback would change our results, the difference may not be large. This is further supported by the findings of Shetty & Ostriker (2008) , who found that the properties of large clouds are not strongly sensitive to feedback. Tests were also performed with a higher cooling floor of 3 × 10 4 K, and no clouds were formed. This perhaps demonstrates that a static cooling floor is a worse approximation to feedback as it inputs energy into any cool region of gas regardless of density, impeding any collapse that would have actually formed stars, in contrast with a dynamic temperature floor which inputs energy only into dense star-forming gas.

Cloud formation & numerical issues
We now draw attention to the differences between the simulations illustrated in Fig. 2. While LowResMW produces clouds at the same time as LowSoftMW (top left), it produces fewer of them as the mass spectrum is truncated. Simlarly, FlatMW produces clouds at the same time as High-ResFlatMW, but in smaller numbers (bottom right). Hence there is a trend of producing more clouds with increasing resolution. Overall, the flat initial conditions of FlatMW and HighResFlatMW produced clouds earlier and in greater numbers than in LowSoftMW. LowViscMW appears identical to LowSoftMW, suggesting that numerical artefacts due to artificial viscosity are not a significant effect (top right). LowFloorMW produced more clouds than LowSoftMW as the lower cooling floor allows the disc to become more unstable to cloud formation from Toomre instabilities. We also found that clouds formed earlier and were more numer-ous with increasing gas fraction, as demonstrated by High-GasMW and MedGasMW (bottom left).
We found that replacing the halo with a static potential did nmot have a significant effect -the mass spectra and number of clouds formed over 430 Myr of evolution were almost identical (Fig. 6).
As expected, the gravitational softening parameter has a significant effect on cloud formation. With a softening of 60 pc (LowSoftMW), a maximum of ∼ 300 clouds were formed at a time of 0.02 Myr, while with a softening of 200 pc (Mid-SoftMW), half as many were formed (∼ 150), and the peak number was achieved later (0.04Myr). It should be noted though, that both models have a similar fraction of mass in clouds (∼ 80%). Increasing the softening yet further to 500 pc (HighSoftMW), leads to almost no clouds forming other than a few clouds in the centre of the galaxy after about a Gyr of evolution (not shown in Fig. 2). These results match what would be expected on theoretical grounds. Increasing the softening length delays cloud formation and produces fewer, more massive clouds, unless the softening length is increased above a certain threshold, beyond which cloud formation is prevented.
It seems most likely that this threshold softening length is related to the wavelength of the unstable mode that causes cloud formation. We can calculate this using the two-fluid (gas/star) Qgs stability parameter from Jog & Solomon (1984), Rafikov (2001), and Li et al. (2005). The individual Q parameters for stars and gas are defined as where Σs and Σg are the gas and stellar surface densities, σs the stellar radial velocity dispersion, cg the gas sound speed, and κ the epicyclic parameter. Note that Qs differs from Toomre (1964)'s definition of Q for a collisionless system by a factor of 3.36/π. If we define q = 2πσs/(κλi), f = cg/σs, where λi is the wavelength of a particular mode of instability, and treat the stars as a fluid with sound speed equal to σs as in Rafikov (2001) with a stability condition of Qgs < 1. We calculate Qgs by using azimuthal means of Ω, Σ, κ, cg and σs, and setting λi to λmin, the wavelength that minimizes Qgs. It is worth cautioning that these parameters are derived from linear perturbation theory and may not adequately describe the system once clouds have formed. Nevertheless, λmin does not rapidly vary from t = 1 Myr to t = 200 Myr for LowSoftC as shown in Fig. 7. λmin is fairly small (< 1kpc) until a radius of 10 kpc at which point it triples in size. This jump is due to the small wavelength gas instabilities starting to dominate over the large wavelength stellar instabilities. A comparison with the face-on density plots (e.g. Fig. 3) shows that clouds form within 10 kpc. In this region λmin is of the order of 100s of pc. The 'threshold' resolution for cloud formation (assuming 4 to 5 softening  lengths are required) in our models lies somewhere between 200 pc and 500 pc, and is consistent with this range. This quantifies an often quoted caveat for galaxy models -if the gravitational softening length is larger than the wavelength of the most unstable modes, then fragmentation is artificially frustrated.
The size of the unstable perturbations can be used to crudely estimate the masses of clouds. Assuming that the disc fragments into clumps of mass ∼ πΣλ 2 min , then for the LowSoftMW simulation (for example) the typical cloud masses should be the order of several 10 6 M , which is admittedly significantly larger than average molecular cloud masses and actually much closer to giant molecular cloud complex masses. Nonetheless, this value is broadly consistent with our spectrum of cloud masses (e.g. Fig. 8). However, we caution against over interpretation as the mass spectrum convolves together an initial spectrum and its subsequent evolution. If this simple approach to calculating ini- tial cloud masses were accurate we would not expect a higher resolution model to produce smaller clouds from this mode of instability, although non-azimuthally symmetric modes which may produce smaller scale instabilities have been excluded from this analysis. Smaller clouds could also be produced in a higher resolution Milky Way model by changing the initial conditions, or if these giant clouds undergo further fragmentation.

Cloud Mass Functions
The mass functions of our clouds (Fig. 8)    feedback. Resolution could potentially also be an issue: although our mass function does not greatly vary between our low and moderate resolution models in our fiducial simulations, our high resolution flat model produced lower mass clouds than the moderate resolution flat model (Fig. 10).

Viscous time-scales
The viscous time-scale is calculated using the method described in section 2.3 and is plotted in Fig. 11. Each point is calculated from 600 collisions. There is a general trend toward lower time-scales as the simulation evolves, and the final time-scales are generally below 10 Gyr, with many approaching 1 Gyr. This decreasing trend coincides with a trend of the number of clouds lowering and the mass of individual clouds increasing. The time-scales are less than a Hubble Time, and so should have some significant effect on the evolution of a galaxy, contrary to the predictions of B02.
This energy loss is seen in the mean specific kinetic energy of the gas in LowSoftMW, dropping from 1.9 × 10 14 erg/g at t = 170 Myr to 1.1 × 10 14 erg/g at t = 1010 Myr. As expected, this loss is primarily in the clouds -the diffuse gas only drops from 1.9×10 14 erg/g to 1.7×10 14 erg/g in the same time period. We can make an additional crude estimate of the total viscous time-scale, tν = (t2 − t1)k1/(k1 − k2), where ti and ki are the time and specific kinetic energy at these two outputs. This gives a time-scale of tν ∼ 2 Gyr, around half of the value of tν = 4.5 Gyr from our method in section 2.3. This is either because additional energy is being dissipated through internal processes in clouds, or because interactions are being missed by our interaction finding procedure. The energy loss of the gas over this time period in LowViscMW is the same to a precision of less than 0.1×10 14 erg/g, suggesting that this additional energy loss is 'physical' and not dominated by artificial viscosity. Nevertheless, we can not evaluate how much of this additional energy loss is directly due to cloud-cloud collisions, and so it is more informative to use the procedure in 2.3 to calculate the viscous time-scale. The mean viscous time-scales from all interactions over each entire simulation for both the Milky Way and collapse models are tabulated in Table 4. Despite the variation of parameters, many of the time-scales are within a narrow range, from 3-5 Gyr. Modifying the artificial viscosity (LowVis-cMW) did not appear to significantly change the viscous time-scale. The softening length in HighSoftMW (600 pc) was large enough to completely quench cloud formation, except for a few clumps that formed within the central bar instability. We do not include a viscous time-scale here as the mechanisms for formation and interaction are different to those of molecular clouds in nearly circular orbits. Feedback processes from star formation and AGN would also be more important here than in the other models. However, lowering the softening length from 100 pc to 60 pc (MedSoftMW to LowSoftMW), while increasing the number of clouds produced, did not significantly alter the viscous time-scale.
HighGasMW has a significantly shorter viscous timescale at 0.6 Gyr, and indeed there appears to be a trend of decreasing viscous time-scale with increasing gas fraction. This is clearer if we compare the models over the same time period. The viscous time-scale over the first 430 Myr is 7.1 Gyr for LowSoftMW, 1.5 Gyr for MedGasMW, and 0.6 Gyr for HighGasMW. Increasing the gas fraction increases the mass of the cloud population (Fig. 9), which increases the frequency and dissipative efficiency of collisions. HighResFlatMW is our highest resolution simulation, but has different initial conditions to LowSoftMW due to the more stringent stability requirements at high resolution (detailed in section 2.2.1). The flat discs of FlatMW and High-ResFlatMW caused cloud formation to occur earlier than in LowSoftMW. A resolution dependence is also evident: The 2.5x increase in mass resolution from FlatMW to HighRes-FlatMW caused a 1.4x increase in the viscous time-scale, and the 5x increase in mass resolution from LowResMW to LowSoftMW caused a 1.8x increase in the viscous time-scale.  Figure 10. Cumulative cloud mass spectra from flat initial conditions, including our highest resolution model.

Monolithic Collapse Model
In all models the gas collapse proceeds as soon as cooling is turned on, thus breaking the hydrostatic equilibrium. The High-S profile slowed the collapse sufficiently for the infalling gas to fragment into clouds at a large radius, although these clouds are are too diffuse to be found by the cloud identification algorithm. As the simulation progresses, these clouds start to merge (from t ∼ 3 Gyr in all runs except for Low-MassC), and reach the effective threshold density of our cloud-finder. The number of clouds quickly reaches a maximum (see Fig. 15). These clouds combine to form a disc. The number and size of clouds these discs fragment into varies greatly between our models.
In HighSoftC, MidSoftC, LowSoftC and LowResC, the disc is extremely unstable, collapsing into ∼ 7 massive (several times 10 9 M in mass) clumps (Fig. 12). These are not small-scale GMC-style clumps as found in the Milky Way simulations, and perhaps this level of collapse is more anal-  Figure 11. Viscous time-scales for disc models that ran for > 800 Myr (left) and 800 Myr (right). At early times, some models give negative time-scales, but as these values are large, they are not as dynamically important and are not plotted.
ogous to the gas-rich clump-cluster galaxies found at high redshift (e.g. Elmegreen & Elmegreen 2005). In the simulations of Bournaud et al. (2007) and Dekel et al. (2009), the large clumps in clump-cluster galaxies coalesce into a central bulge, forming a more stable disc. These simulations differ to ours particularly in that they include star-formation and feedback. With infalling material, Dekel et al. (2009) finds the clumpy phase can last for several Gyr.
The heavy clustering in this discs dictated that they could only be evolved for < 1 Gyr after formation (which takes ∼ 3 Gyr) due to problems with the SPH solver. The high densities cause a large increase in the number of particles with smoothing lengths at the minimum allowed which contributes to an O(n 2 ) slowdown.
The simulations of Kaufmann et al. (2009), while including star-formation (but not explicit feedback), also produce a disc with large-scale gravitational instabilites. Both our and Kaufmann's models have a temperature floor of 3 × 10 4 K, as a very crude form of feedback. Including starformation and more self-consistent feedback method could produce a stable disc (Stinson et al. 2006;Christensen et al. 2010).
In LowSoftFloorC, the low temperature floor allows the halo clouds to condense into dense (n ∼ 10 4 − 10 5 cm −3 ) clumps (Fig. 13). Their low cross-section means that their coalescence has properties of a collisionless collapse. So in addition to an unstable disc, there exists a swarm of clumps with a half-mass height of 7.8 kpc. Their ellipsoidal distribution and high densities are reminiscent of globular clusters, but the inclusion of feedback would definitely increase the cloud cross-sections and produce a more dissipated and flattened disc.
LowMassC is the only run that produces a disc that does not collapse into large clumps (Fig. 14), although it took considerably longer to form (∼ 4.5 Gyr) and the disc is still dominated by spiral instabilities. Discs are unstable to bar formation when the disc mass fraction is greater than the spin parameter (m d > λG) (Efstathiou et al. 1982;Foyle et al. 2008), so a lower mass disc is more stable. If the bar is too strong, it may fragment into large clumps. This in- Figure 14. LowMassC at t = 6.0 Gyr. The disc undergoes spiral instabilities but does not fragment into clumps as the other collapse models do. stability may well drive the infalling clouds into a few large clumps in the higher mass models.
As seen in Table 4, the viscous time-scales for the collapse runs trend toward lower values than the Milky Way simulations -around 1 − 2 Gyr. Though the number of interactions is not as large as in the Milky Way models, they occur over a short period (e.g. all 566 interactions in Low-SoftC are within ∼ 500 Myr). The number of clouds is small, so each cloud undergoes many collisions, producing a short viscous time-scale.

Comparison with Analytical Model
B02 argued that while cloud collisions are not uncommon (occurring > ∼ 1 time per orbit), the low efficiency of cloud collisions produces a long viscous time-scale. This efficiency is measured with a parameter η, equal to the fraction of a cloud's energy that is lost in a collision (not entirely dis- Figure 12. Impact of varying the softening length and resolution in collapse runs at t=3.5 Gyr. Top left is HighSoftC (514 pc, 3 × 10 4 K), top right is MidSoftC (200 pc, 3 × 10 4 K), bottom left is LowSoftC (60 pc, 3 × 10 4 K) and bottom right is LowResC (60 pc, 300 K). Although HighSoftC, MidSoftC and LowSoftC produce different numbers of clouds initially (more clouds for a shorter softening length), after ∼ 500 Myr of collisions all three models have ∼ 7 large clumps. Despite the low temperature floor, the limited resolution of LowResC produces an unstable disc, instead of a swam of dense clumps as in LowSoftFloorC. similar from a coefficient of restitution). When two clouds merge completely, the fraction of kinetic energy lost is well approximated by η = (v rel /vrot) 2 , where v rel is the relative velocity of the clouds, and vrot is their rotational velocity which is roughly constant for a galaxy. This is consistent with our numerical results. The analytical model of B02, finds that η 10 −2 for a Milky-Way-like model, concluding that cloud-cloud collisions are not an efficient sink of energy, with tν ∼ 1000-2000 Gyr.
The complex interactions that occur between clouds in our simulation mean that it is not straightforward to determine values for η. Several of our merger and separation events can take place within what is really a single extended interaction, which lowers the average time between interactions significantly. Indeed, we find the interaction rates are on the order of one separation or merger event per cloud every 50 − 60 Myr for LowSoftMW, Med-SoftMW, LowFloorMW, LowViscMW, HighSoftC and Low-MassC. The greatest interaction time-scale was in , and the smallest was in LowSoftC (14 Myr).
It is difficult to track the number of interactions over a full merger process, as additional clouds often interact with the merging clouds. We carefully examined examined a span of time around each of a sample of 10 recorded interactions in LowSoftMW on an iteration-by-iteration basis to determine the number of recorded interactions per 'real' interaction. These interactions were selected so that they were evenly distributed across the simulation (∼ 2 every 5000 iterations). We initially examined a period of ±800 iterations around the interaction, and if no 'real' interaction was observed during To smooth the data, each plotted point is an average of the 29 data points centred on it. Being very unstable, these systems formed a few large clumps rather than many small clumps.
this time, this was extended to ±2000 iterations. Several different behaviours were observed: • In two cases, no real interactions were observed; outer parts of a cloud were attaching and detaching to the main cloud, and dissolving and condensing across the cloud density threshold, causing a number of recorded interactions which did not correspond to any clear long-term merger, scattering or separation event.
• Three events were 'messy' interactions with 6, 7 and 16 recorded interactions per real event; the event consisting of 16 recorded split and merge events was a scattering event where the clouds passing by each other several times before separating for a final time.
• The last event was a series of mergers in rapid succession -3 recorded and 3 'real' mergers.
Overall, there was a mean of 4.9 recorded interactions per examined period, with a standard deviation of 4.3. A total of 11 'real' interactions were observed, giving 4.5 recorded interactions per real interaction. This increases our interaction time-scale to one event per ∼ 250 Myr for the LowSoftMW-like models. This is approximately once per orbit at a solar radius. The analytic estimate in B02 of the cloud-cloud collision rate is ∼ 100 Myr, which is of similar order.
We can estimate an η for the interactions in our models by where ∆K and ∆φ are the change in kinetic and potential energy of a cloud, and Kc is the total kinetic energy of both clouds before collision. η can be negative, as energy is converted from internal motions into orbital kinetic energy during separations. The clouds all have similar velocity because of the flat rotation curve, so the total energy lost is primarily dependent on η and the cloud masses. We find for most interactions |η| is on the order of ∼ 0.002 (Fig. 16). If we separate our η values into two sets, η− for η < 0 and η+ for η > 0, we find that the median value of |η−| is greater than the median value of |η+|, even though the viscous timescale is positive. This is because although η, the relative energy change is larger for interactions which increase orbital energy (η− < 0) than those which decrease orbital energy (η+ > 0), the absolute change in energy is larger for interactions which decrease orbital energy than increase it -i.e. these interactions tend to occur between clouds with greater mass. Although it is not apparent on these plots, there are several collisions for which η is very large, with η > 0.1. These interactions occurred within the 1 kpc of galaxy centre, and only after ∼ 400 Myr. These are clouds that have been strongly scattered by interactions and fallen down the potential well, colliding with speeds of > 100 km s −1 . Our interactions are no more efficient at removing energy than in B02, and are no more common, yet the B02 model predicts tν ∼ 1000-2000 Gyr, while our simulations have tν < 10 Gyr. Our simulated discs are more energetic than standard Milky Way models: the velocity dispersion in LowSoftMW is ∼ 20 km s −1 at 7.5 kpc, more than triple the standard Milky Way value used in B02 (6 km s −1 ). However, this is not the cause of the large difference between the model of B02 and our own. Here we derive our own model for η, and contrast this with the model in B02 to find the source of this disparity.
We can split the velocity components of v rel into tangential and radial components to give If we make the epicyclic approximation (Binney & Tremaine 2008), that the deviation from a circular orbit is small compared to the radius of the orbit (R = Rg +x, where Rg is the 'guiding centre' of the orbit, and x R is the radial excursion), thenṘ =ẋ = Xκ cos(κt + α), (where X is the maximum radial excursion of a cloud, kappa the epicyclic frequency, and α is a phase parameter) andφ = Rgvrot/R 2 from conservation of momentum in a flat rotation curve. Hence R 2 (φ1 −φ2) 2 = (v 2 rot /R 2 )(Rg,1 − Rg,2) 2 -the tangential component of the difference in velocity depends only on the radial distance between the clouds' guiding radii.
The radial component is more difficult to calculate, as it depends on the phase of the interaction. We can estimate the maximum η by assuming the clouds are perfectly out of phase -that is, as R ∼ Rg. For clouds to collide precisely out of phase, they must have the same guiding radius, and so Rg,1 − Rg,2 = 0. Hence, if X1 ∼ X2 ∼ X, then ηmax,r = 8X 2 /R 2 . If the clouds are at their maximum deviation when they collide, then their radial velocities are zero, but their relative φ velocities are maximised, that is,φ1 −φ2 = 2X, and so η max,φ = 4X 2 /R 2 . These coefficients give the maximum η, but we should nevertheless expect η ∼ X 2 /R 2 , i.e. η depends on the radial excursion of clouds.
This can also be expressed in terms of a velocity dispersion. We can calculate the velocity dispersion by Assuming a flat rotation curve and that X and κ are more or less constant within the region of interest, the radial component isẋ = Xκ cos(κt + α) , hence and the tangential component is R(φ − Ωg) = −2XΩg sin(κt + α), hence This gives v 2 s = 3v 2 rot (X 2 /R 2 ), and so From these expressions for η we can determine the dissipative time-scale from tν = tc/η.
We next summarise the model of B02. In the limit of rapid collisions, the kinematic viscosity due to cloud-cloud collisions can be modelled as a Reynolds stress and expressed as ν ∼ λ d vs (Faber 1995), where vs is the velocity dispersion, and λ d is the mean free path. The mean free path is λ d = vstc, where tc is the typical time between collisions. Similarly to our result, B02 states η ∼ ∆R 2 /R 2 , where ∆R is the radial distance between collisions. For the case of very rapid collisions, ∆R ∼ λ d , so η ∼ λ 2 d /R 2 . This gives a viscous time-scale that should be equal to the dissipative time-scale, Hence if we follow the description given in B02, the results should be equivalent to ours. Continuing to follow B02, we can set tc = M cloud h/(Σgvsπr 2 cloud ), so We can evaluate this using the Milky Way parameters of B02, that r = 7.5 kpc, vrot = 220 km s −1 , vs = 6 km s −1 , Σg = 50M pc −2 , M cloud = 10 5 M , h = 100 pc, and r cloud = 10 pc to result in tν = 14 Gyr. However, B02 states tν ∼ 2000 Gyr. This disagrees by a factor of 1/η. It appears that B02 includes an additional factor of η ∼ 0.008 in the denominator -i.e. t ν,Bell ∼ R 2 /(ην). This η is not necessary, as it is already included in the radial excursion or velocity dispersion, and as is clear from equation 31, the expression tν ∼ R 2 /νη is not equivalent to the dissipative time-scale.
In B02 's rare collision case, ν ∼ vs∆R(tκ/tc), where tκ = 2π/κ is the epicyclic time-scale. For a flat rotation curve κ = √ 2Ω ∼ v0/R. The excursion ∆R is on the order of the radial excursion of the epicyclic motion of the clouds. B02 state ∆R ∼ vs/κ ∼ vsR/vrot, which is consistent with our result in equation 29. Putting this together gives i.e. η ∼ 2πv 2 s /(v 2 rot ) ∼ 0.023. B02 uses a low surface brightness galaxy in this case, with Σg = 10M pc −2 and vrot = 100 km s −1 , which results in tν ∼ 23 Gyr. Again, the value in B02 is much larger, tν ∼ 1000 Gyr, which again is higher than our calculated value by a factor of approximately 1/η. These model are intended to apply in the limits of very frequent or very infrequent collisions where tcΩ 1 or tcΩ 1. In our simulations, we found that clouds collide about once per orbit, i.e. Ωtc ∼ 1. However, we can contrast these results with those of Goldreich & Tremaine (1978), who solve the Boltzmann equation for a system of inelastically colliding particles in a disc, and find for arbitrary Ωtc that the viscosity is of order after we make the substitution that λ d ∼ vstc. For Ωtc = 1, ν = 1/2(λ d vs). The frequent collision case of B02, ν ∼ vsλ d , is accurate to this within an order of magnitude if we exclude the erroneous factor of 1/η. Substituting our typical cloud and disc parameters at 7.5 kpc (h ∼ 25 pc, Σg ∼ 100M /pc −2 , vs ∼ 20 km s −1 , r cloud ∼ 35 pc, and M cloud ∼ 10 7 M ) for LowSoftMW at t = 1 Gyr into this model gives a viscous time-scale of 1.1 Gyr. This values somewhat underestimates our numerical results for the Milky Way models in Table 4, for most of which tν 4.0 Gyr. The unstable disc of LowSoftC, forming from a collapse without stars, has very different properties at R = 7.5 kpc, with h ∼ 250 pc, Σg ∼ 5000M /pc −2 , vs ∼ 100 km s −1 , r cloud ∼ 100 pc, and M cloud ∼ 10 9 M . This gives tν = 0.35 Gyr, which agrees with our simulation result (0.8 Gyr) within a factor of ∼ 2. The analytical expression for tν was evaluated from order-of-magnitude arguments and assumptions that may not be entirely valid in our simulations -particularly in models with very few clouds, such as LowSoftC. Numerical factors also vary our simulation results by a factor of ∼ 4. Given these issues, it is not surprising that the agreement is not exact.
Interestingly, despite the different disc properties, Low-SoftC and LowSoftMW have similar viscous time-scales in Correlation between peak number of clouds (N cloud,max ) and viscous time-scales (tν ) for all models whose time-scale is given in Table 4. The plotted fit is tν = (0.67 Gyr)(N cloud,max ) 0.39 . both our numerical simulations and in this analysis. This is to be expected from equation 32. We should expect the typical cloud mass to increase with the gas density and typical cloud radius, and so M cloud /(πr 2 cloud Σg) should vary only weakly. Hence the viscous time-scale will primarily depend primarily on h and vs. This suggests that timescales will not vary greatly for models beyond those simulated here -perhaps even of higher resolution. To quantify this, we note that there appears to be a correlation between the maximum number of clouds formed (N cloud,max ) and the viscous time-scale (Fig. 17). Performing a fit to a power-law tν ∝ (N cloud,max ) m , we find a power index of m = 0.39 ± 0.19. This predicts a viscous time-scale of tν ∼ 23 Gyr for N cloud,max = 10 4 , and tν ∼ 60 Gyr for N cloud,max = 10 5 , although we caution that this is a purely empirical fitting and is not likely to be very accurate.

CONCLUSIONS
Previous estimates of the viscous time-scale suggest that the viscous time-scale for cloud-cloud collisions in a Milky-Way-like galaxy is large, with tν > 1000 Gyr. To test the hypothesis that the viscous time-scale is long, we performed simulations using the AP 3 M SPH code HYDRA with cooling down to 10 K and a dynamic temperature floor. The simulations fell into two sets of models: initially stable gaseous discs within a dark-matter haloes and stellar discs, and a gaseous spheres collapsing inside dark-matter haloes. These two models were chosen to bracket a wide range of stability. The viscous time-scale was measured by tracking clouds with a friends-of-friends algorithm, and determining the energy loss when clouds collided.
Although our cloud masses are larger than those found in other simulations, potentially due to insufficient resolution, a simple analysis suggests that we are resolving the wavelength of the most unstable mode. However, further instabilities (in particular, non-axisymmetric instabilities that we exclude from our analysis) may appear at higher resolutions, and while the inclusion of energy input from stellar feedback may not greatly alter the properties of clouds, it may contribute to cloud evaporation and affect their collisional behaviour by increasing their cross-section through heating.
Identifying clouds and interactions between clouds is still a difficult task, as clouds have complex structures and dynamics. The friends-of-friends algorithm often identifies clouds as merging and separating several times over a period that upon visual inspection appears to be a single interaction. Through a detailed examination of 10 interaction events, we determined that each 'real' interaction corresponds to ∼ 4.5 interactions found by our algorithms. This also complicated our estimates for η = ∆K cloud /K cloud , the efficiency of energy loss per cloud interaction. We found that despite our low viscous time-scales, η was not large, with η ∼ 0.002 per recorded interaction.
Most models from both sets of initial conditions collapsed into discs dominated by clumps of gas. The Milky Way models produced a more stable disc with a large number of small clouds, while the collapse models produced a highly unstable disc consisting of a small number of massive clumps. Despite this large disparity, the viscous timescales were similar, with tν = 4.5 Gyr for LowSoftMW, and tν = 0.8 Gyr for LowSoftC. These values are much smaller than estimates using the formulation of B02, which overestimate the viscous time-scale by appearing to erroneously include inefficiency of cloud collisions twice. Removing this factor gives analytic estimates of tν = 1.1 Gyr for Low-SoftMW and tν = 17 Gyr for LowSoftC. These values do not exactly coincide with our measured values as they are based on simple arguments that are particularly inaccurate for LowSoftC. However, they all agree with the general statement that viscosity due to cloud-cloud collisions is not negligible.
The scatter of tν across our models was quite small (0.6-16.0 Gyr), despite the range of cloud properties. Hence our viscous time-scales are applicable for a wider range of galaxies than those modelled here, although viscous timescales will likely increase somewhat as resolution improves. For a simulation capable of resolving 10 5 clouds, we predict a viscous time-scale of around 60 Gyr, admittedly making the effect comparatively weak within a Hubble time, but nonetheless over an order of magnitude faster than previous estimates.
These results suggest that viscosity due to cloud-cloud collisions, while not dominant, does not have a completely negligible effect on the evolution of a galaxy. Although our models may underestimate the viscous time-scales due to resolution effects, it still appears that cloud-cloud viscosity is more significant than previously estimated. While numerical models of galaxies may be able to model this directly (as we do in this work), it may be necessary to include a cloudcloud viscous term in analytical and semi-analytical models of disc evolution.