On the CCN [de]activation nonlinearities

We take into consideration the evolution of particle size in a monodisperse aerosol population during activation and deactivation of cloud condensation nuclei (CCN). The phase portrait of the system derived through a weakly-nonlinear analysis reveals a saddle-node bifurcation and a cusp catastrophe. An analytical estimate of the activation timescale is derived through estimation of the time spent in the saddle-node bifurcation bottleneck. Numerical integration of the system portrays two types of activation/deactivation hystereses: one associated with the kinetic limitations on droplet growth when the system is far from equilibrium, and one occurring close to equilibrium and associated with the cusp catastrophe. The hysteretic behaviour close to equilibrium imposes stringent time-resolution constraints on numerical integration, particularly during deactivation.

Abstract. We take into consideration the evolution of particle size in a monodisperse aerosol population during activation and deactivation of cloud condensation nuclei (CCN). The phase portrait of the system derived through a weaklynonlinear analysis reveals a saddle-node bifurcation and a cusp catastrophe. An analytical estimate of the activation timescale is derived through estimation of the time spent in the saddle-node bifurcation bottleneck. Numerical integration of the system portrays two types of activation/deactivation hystereses: one associated with the kinetic limitations on droplet growth when the system is far from equilibrium, and one occurring close to equilibrium and associated with the cusp catastrophe. The hysteretic behaviour close to equilibrium imposes stringent time-resolution constraints on numerical integration, particularly during deactivation.

Background
Atmospheric clouds are visible to human eye for they are composed of water and ice particles that effectively scatter solar radiation. The multi-micrometre light-scattering cloud droplets form on sub-micrometre aerosol particles (cloud condensation nuclei, CCN) in a process referred to as CCN activation. The concentration (100-1000 cm -3 ) and mean size (1-10 µm) of activated particles can vary by an order of magnitude depending on the size spectrum and composition of CCN. On one hand, CCN physicochemical properties are influenced by anthropogenic emissions of particles into the atmosphere. On the other hand, the resultant size spectrum of cloud droplets determines how effectively the clouds interact with solar radiation and how effectively they produce precipitation (see e.g. a recent NPG paper by Feingold and Koren, 2013, for a discussion of aerosol-cloud-precipitation interaction chain, unconventionally modelled as a predatorprey problem). CCN activation is thus the linking process between the microscopic human-alterable atmospheric composition and the macroscopic climate-relevant cloud properties. As once aptly stated, "there is something captivating about the idea that fine particulate matter, suspended almost invisibly in the atmosphere, holds the key to some of the greatest mysteries of climate science" (Stevens and Boucher, 2012). This has certainly contributed to the wealth of literature on the subject published since the first studies of the 1940-ties (Howell, 1949;Tsuji, 1950), for a thorough list of references see e.g., Khvorostyanov and Curry (2014, chpt. 7).
Deactivation is the reverse process in which cloud droplets evaporate back to aerosol-sized particles. The process is also referred to as aerosol regeneration, aerosol recycling, dropto-particle conversion or simply droplet evaporation (see section 1 in Lebo and Seinfeld, 2011, for a review of modelling studies). Both activation and deactivation are particular cases of particle condensational growth which, in context of cloud modelling, is generally regarded as reversible to contrast the irreversible collisional growth (see e.g., Grabowski and Wang, 2013). The reversibility of condensational growth is a sound (and often a constituting) assumption for cloud models for which activation and deactivation are subgrid processes, both in terms of time-and length-scales. Yet, when investigated in short-enough timescales, condensation and evaporation exhibit a hysteretic behaviour in an activationdeactivation cycle. The hysteresis can be associated with the kinetic limitations in the vapour and heat transfers to/from the droplets (Chuang et al., 1997) and has been previously depicted in the studies of Korolev and Mazin (2003, discussion of Fig. 1) and Korolev et al. (2013). As we point out in this note, the system can exhibit a hysteretic behaviour also in a close-to-equilibrium régime where the kinetic limitations do not play a significant role.

Droplet growth laws in a nutshell
The key element in the mathematical description of CCN activation/deactivation is the equation for the rate of change of particle radius r w (so-called wet radius) due to water vapour transfer to/away from the particles. It is modelled by a diffusion equation in a spherical geometry: where ρ w is the liquid water density, ρ v is the ambient vapour density (away from the droplet), ρ • is the equilibrium vapour density at the drop surface and the D eff = D eff (T, r w ) is an effective diffusion coefficient in which the temperature dependence stems from an approximate combination of the Fick's first law and Fourier's law (latent heat release) into a single particle-growth equation (so-called Maxwell-Mason formula), while the radius dependence stems from corrections limiting the diffusion efficiency for smallest particles. For derivation and discussion see Khvorostyanov and Curry (2014, section 5.1.4 1 + A/r w -ϰ r d 3 / r w 3 c 0 + c 2 ξ c 2 Figure 1. Köhler curve for CCN with rd = 0.05µm, κ = 1.28 (NaCl) and its Taylor expansions at rc and at infinity.
respect to plane surface of pure water) and the equilibrium relative humidity RH eq = ρ • /ρ vs , the drop growth equation is given by: The crux of the matter is the dependence of RH eq on r w due to the droplet curvature and due to the presence of dissolved substances. It is given by the so-called Köhler curve, the leading terms of its common κ-Köhler form can be approximated with (for r d r w what is an reasonable assumption in context of activation/deactivation): where A ∼ 10 −3 µm is a temperature-dependant coefficient related with surface tension of water, while the dry radius r d and the solubility parameter κ (in general, 0 < κ < 1.4, see Petters and Kreidenweis, 2007) are proxy variables depicting the mass and chemical composition of the substance the CCN are composed of. The ∂ ∂r w RH eq derivative has an analyticallyderivable root corresponding to the maximum of the Köhler curve at (r c , RH c ) where r c = 3κr 3 d /A is the so-called critical radius and RH c = 1+ 2A 3rc is the critical relative humidity. Neglecting the details of the activation kinetics, the CCN can be considered activated (deactivated) when they grow beyond (shrink below) their critical radius.

Saddle-node bifurcation at Köhler curve maximum
Rewriting equation 2 in terms of ξ = r 2 w + C (where C is an arbitrary constant) gives: Rewriting RH eq in terms of ξ c = r 2 w − r 2 c and Taylorexpanding it around ξ c = 0 gives: where c 0 = RH c , c 1 is zero as we are expanding around the root of ∂ ξc RH eq and c 2 = − A 4 r −5 c is negative. Combining equations 5 and 6 gives: which is the normal form of the saddle-node bifurcation (Strogatz, 2014, section 3.1). The phase portrait of the system can be recognised in a standard cloud-physics Köhler curve plot given in Figure 1. In the nearest vicinity of the critical radius, when RH approaches RH c from below, there are two fixed points in the system: one stable (for r w < r c , unactivated CCN) and one unstable (for r w > r c , activated CCN). The bifurcation occurs when RH = RH c when the fixed points coalesce into a half-stable fixed point. When RH > RH c , there are no fixed points at all.

Activation timescale estimation
Interestingly, the analysis of the CCN activation/deactivation in terms of saddle-node bifurcation provides a way to estimate the timescale of activation. Following Strogatz (2014, section 4.3), the coalescence of the fixed points is associated with a passage through a bottleneck. The key observation is that for the parabolic normal form of the saddle-node bifurcation, the time of the passage through the bottleneck dominates all other timescales. Thus, the timescale of the process can be estimated by integrating ξ c from −∞ to ∞: The activation timescale τ act given by eq. 8, plotted as a function of RH and r d (and substituting r c and RH c by their analytic formulae given in the preceding section) is presented in Figure 2. It matches remarkably the data obtained through numerical calculations presented in Hoffmann (2016). The white region in the plot corresponds to situation where activation does not happen. The range of RH depicted in the plot is chosen to match the one of Figure 2 in Hoffmann (2016), while in principle the presented weakly non-linear analysis of the system is applicable only close to the equilibrium (i.e., close to the edge of the white region in the plot).

Cusp catastrophe of the RH-coupled system
The key limitation of the preceding analysis is that the evolution of particle size is not coupled with the evolution of ambient heat and moisture content, and hence the relative humidity. Limiting the analysis to a monodisperse population, the coupling efficiency is determined by the total number of particles in the system. The so-far assumed constant RH approximates thus the case of small number of droplets.
To at least partially lift the constant-RH assumption, while still allowing for a concise analytic description of the system dynamics, let us consider a simple representation of the moisture budget in the system under a temporary assumption of constant temperature and pressure (and hence constant volume, constant ρ vs , A and D eff ). The rate of change of the ambient relative humidityṘH can be expressed then as a function of the droplet volume concentration N : where the form of α stems from defining the density of liquid water in the system as N ρ w 4 3 πr 3 w . Integrating in time gives: which combined with eq. 2 and expressed in terms of ξ with C = 0 leads to the following phase portrait of the RHcoupled system (assuming r w r d ): where the group of terms labelled as f can be intuitively thought of as corresponding to the Köhler curve with an additional term representing the simplified RH dynamics. Figure 3 depicts the dependence of f on the droplet radius r w = √ ξ and the droplet concentration N . To facilitate analysis, the zero-crossings of the first derivative of f with respect to r w are plotted as well using analytically derived formula: For N = 0, f has the Köhler-curve shape depicted in Fig. 1 which, as discussed in the preceding sections, implies a saddle-node bifurcation. With N greater than zero but less than ca. 50/cm 3 , a second saddle-node point appears as the αN ξ 3 2 term causes f to have a local minimum above the critical radius. At ca. N=50/cm 3 , both the first and second derivatives of f vanish implying a cusp point in the f surface. For larger N , f is monotonic, hence the saddle-node bifurcation ceases to exist in the system. This phase portrait reveals a cusp bifurcation with the plotted zero-crossing line depicting its stability diagram. The cusp bifurcation (Kuznetsov, 2004, chpt. 8.2), an imperfect supercritical pitchfork bifurcation (Strogatz, 2014, chpt. 3.6), features a cusp catastrophe what allows to envision a "catastrophic" jump from one equilibrium to another and a hysteretic behaviour of the system when approaching (in therm of r w ) the local minimum of f from below (activation) and from above (deactivation) for small enough N .

Adiabatic vertically-displaced air parcel system
In order to lift the assumptions of constant temperature and pressure, the system evolution can be formulated by supplementing the drop growth equation with two equations representing the hydrostatic balance and the adiabatic heat budget. This leads to a commonly used so-called air-parcel framework depicting behaviour of a vertically displaced adiabatically isolated mass of air: where ρ d and p d are the dry-air (background state) density and pressure, w is the vertical velocity of the parcel, q = ρ v /ρ d is the water vapour mixing ratio, c pd is the specific heat of dry air, l v is the latent heat of vapourisation and g is the acceleration due to gravity. As discussed in section 5, for a monodisperse population of N particles, the changes in the mass of liquid water in the system are proportional to the particle concentration, henceq ∼ N . Consequently, the analysis of the activation/deactivation dynamics presented in sections 3-4 under the assumption of constant RH corresponds to the behaviour of the air-parcel system defined by eq. 13 in the limit of: w → 0 (and henceṗ d ≈ 0) i.e., slow, close-toequilibrium evolution of the system relevant to fixedpoint analysis (by some means pertinent to the formation of non-convective clouds such as fog) and -N → 0 (and henceṙ ≈ 0) i.e., weak coupling between particle size evolution and the ambient thermodynamics (pertinent to the case of low particle concentration).

Numerical simulations
Since the system defined by (13) is less susceptible to a simple analytic analysis, we proceed with numerical integration. Furthermore, employing numerical integration allows to evaluate the Köhler curve in unapproximated form (4) to corroborate the findings obtained with the assumption of r d r w . To this end, a numerical solver was implemented using the libcloudph++ library (Arabas et al., 2015) and the CVODE adaptive-timestep integrator (Hindmarsh et al., 2005). The solver code is free and open-source and is available as an electronic supplement to this note. In order to depict an activation-deactivation cycle, the vertical velocity w was set to a sinusoidal function of time t such that the maximal displacement is reached at t = t hlf and the average velocity is <w>: w = <w> π 2 sin π t t hlf (14) Figure 4 summarises results of nine simulations in three types of coordinates: displacement vs. supersaturation (the top row), supersaturation vs. wet radius (the middle row, same coordinates as in Fig. 1) and displacement vs. wet radius (bottom row). The nine model runs correspond to three sets of aerosol parameters (left, middle and right columns) and three values of mean vertical velocity (depicted by line thickness). The varied aerosol input parameters are the concentration (N STP of 50 and 500 cm −3 , STP subscript corresponding to the values at standard temperature and pressure) and the dry radius (r d of 0.1 and 0.05 µm). In all panels, black lines correspond to air-parcel ascent (activation) and orange lines correspond to the descent (deactivation). Besides integration results, the panels in the middle row feature the Köhler curve plotted with thick grey line in the background.
The plots depict that for mean velocities of 100 cm/s and 50 cm/s activation and deactivation are not symmetric and happen far from equilibrium (the Köhler curve). This type of hysteresis corresponds to the kinetic limitations on the transfer of water molecules to/from the droplet surface what prevents the droplets from attaining equilibrium under rapidly changing ambient conditions.
At much lower velocity of 0.2 cm/s, the processes are symmetric and match the equilibrium curve, but only for the N = 500 cm −3 and r d = 0.1 µm (middle column). A twofold decrease of the dry radius (right column) as well as a tenfold decrease of particle concentration (left column) both cause the system to exhibit a hysteretic behaviour also at the lowest considered velocity. This hysteresis is characterised by a "jump" in the wet radius that qualitatively matches the envisioned catastrophic behaviour associated with the cusp bifurcation. This behaviour is robust to further reduction in the vertical velocity (not shown) confirming a close-toequilibrium régime was attained.
The adaptive-timestep solver statistics (not shown) reveal that regardless of the chosen relative accuracy, for all considered input parameters, there are two instants for which the solver needs to significantly reduce the timestep: when resolving the supersaturation maximum during activation and when resolving the "jump" back to equilibrium during deactivation. It is a robust feature that deactivation requires roughly an order of magnitude shorter timestep as compared to activation (ca. 0.01 s vs. 0.1 s for a relative accuracy of 10 -6 ). The only exception from this rule is the symmetric case which does not feature the (catastrophic) "jump" back onto the equilibrium curve. This confirms that it is the hysteretic behaviour that imposes the tightest constraints on the timestep.

Monodisperse system: limitations and applicability
The key advantage of the embraced monodisperse simulation is simplicity -in terms of model formulation, result analysis but also integration. Simulations of the particle size spectrum evolution during activation are prone to numerical difficulties both due to the stiffness of the system as well as due to the sensitivity to the size spectrum discretisation (Arabas and Pawlowska, 2011).   The key inherent limitations for applicability of monodisperse simulations are the lack of description of the cloud droplet size spectrum shape and the lack of differentiation between activated and unactivated CCN. The fact that all considered particles activate also implies that the model is unable to resolve the noise-induced excitations to which the system is susceptible. The system exhibits an excitable behaviour if subject to fluctuations in the forcing terms (e.g. in the cooling rateṪ , see Hammer et al., 2015, discussion of Fig. 10-11 and other studies referenced therein). The excitations influence the partitioning between activated and unactivated CCN, and decay when the characteristic timescale (period) of fluctuations is largely longer or shorter than the activation timescale discussed in section 4.
These limitations certainly restrain the relevance of the presented calculations to real-world problems. Yet, let us underline that both the monodisperse spectrum and even the no-RH-coupling assumption are in fact contemporarily used in atmospheric modelling in the recently popularised particlebased (Lagrangian, super-droplet) techniques for representing aerosol, cloud and precipitation particles in models of atmospheric flows (see Shima et al., 2009;Arabas and Shima, 2013, and works referred therein). In these models, in the spirit of the particle-in-cell approach, the liquid water is represented with computational particles, each representing a multiplicity of real-world particles with monodisperse size. In such models, the particles undergo repeated activationdeactivation cycles. Consequently, the close-to-equilibrium catastrophic hysteresis observed in the presented simulations, even if of no foreseeable relevance to the macroscopic behaviour of the large-scale cloud systems modelled with the particle-based techniques, has to be taken into account when developing numerical integration schemes.

Concluding remarks
With this note we intend to bring attention to the presence of nonlinear peculiarities in the equations governing CCN activation and deactivation, namely a saddle-node bifurcation and a cusp catastrophe. We have shown that conceptualisation of the process in terms of bifurcation analysis yields a simple yet practically-applicable description of the system allowing analytic estimation of the timescale of activation. Both through weakly-nonlinear analysis and through numerical integration, we have depicted the presence of a cusp catastrophe in the system and the corresponding hysteretic behaviour near equilibrium (i.e., for very small air-parcel velocities). The near-equilibrium hysteresis was observed to determine the timestepping constraints for numerical integration when simulating an activation-deactivation cycle of a monodisperse droplet population. It is a finding of interest for cloud modelling community since monodisperse activation/deactivation models of the studied type play a constituting role in the more-and-more widespread particle-based models of aerosol-cloud interactions.