David A Bonhommeau
email: david.bonhommeau@univ-reims.fr

A practical law to predict the appearance sizes of multiply charged rare-gas and molecular clusters

Keywords: Rare-gas clusters, molecular clusters, critical sizes, Lennard-Jones interactions, liquid drop model

A dimensionless law depending on cluster size N is derived from the Rayleigh limit expressed in reduced Lennard-Jones (LJ) units to predict the critical sizes n c (z) of clusters carrying z positive charges. This relationship provides suitable estimates of n c (z) (z = 2-4) for rare-gas clusters, including neon clusters whose predicted critical sizes deviate from experimental expectations by less than 12% for different choices of LJ parameters. An extension to 11 nonpolar and 15 polar molecular clusters, from dimers to aromatic hydrocarbons, is achieved that demonstrates the broad applicability of the formula despite inaccuracies for highly polar systems.

Introduction

The search for critical (or appearance) sizes n c (z) above which clusters with z charges are stable is a long-term endeavour. The existence of highly charged metal clusters (z ≤ 14) has been first observed on sodium clusters by Näher et al. [START_REF] Näher | Observation of highly charged sodium clusters[END_REF]. Although the charge is expected to be distributed over the surface of alkali metal clusters [START_REF] Näher | Fission of Metal Clusters[END_REF], molecular dynamics simulations of uniformly charged sodium clusters following a tight-binding approach [START_REF] Bulgac | Thermal properties of Na 8 microclusters[END_REF][START_REF] Calvo | Geometric size effects in the melting of sodium clusters[END_REF] yielded n c (2) = 27 and n c (3) = 63 [START_REF] Li | Molecular dynamics study of neutral and multiply charged sodium clusters[END_REF], in close agreement with experimental data [START_REF] Näher | Observation of highly charged sodium clusters[END_REF]. However, the critical size obtained for z = 4, namely 110, slightly underestimated the experimental values, derived from ion bombardment experiments, that ranged from 122 ± 2 (H + at 5 keV) to 115 ± 5 ( 40 Ar 8+ at 40 keV) [START_REF] Chandezon | Critical Sizes against Coulomb Dissociation of Highly Charged Sodium Clusters Obtained by Ion Impact[END_REF].

Uniform charge distribution has also been invoked to reproduce the critical sizes of van der Waals and hydrogen-bonded clusters [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF]. Clusters are then modeled as continuous media of dielectric constant ǫ and their total energy is written as the sum of a volume term, a surface term, and a Coulomb term by analogy to liquid drop models (LDMs). Based on the work of Briant and Burton [START_REF] Briant | Molecular dynamics study of the structure and thermodynamic properties of argon microclusters[END_REF] devoted to argon microclusters, Echt et al. expressed the surface energy E s as E s = (AN 2/3 + BN 1/3 )V LJ [START_REF] Näher | Observation of highly charged sodium clusters[END_REF] where N is the cluster size, V LJ is the corresponding Lennard-Jones (LJ) pair well depth, and A = 10.3 and B = -2.7 are dimensionless parameters. This model properly predicts the critical sizes of doubly charged krypton, hydrogen sulfide, or 1,1-difluoroethane clusters but more significant disagreements arose for neon and nitrogen clusters [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF].

The same expression of surface energy has been recently employed to ascertain the values of the maximum charge a multiply charged cluster can sustain before becoming unstable [START_REF] Bonhommeau | Rayleigh limit and fragmentation of multiply charged Lennard-Jones clusters: Can charged clusters provide clues to investigate the stability of electrospray droplets?[END_REF]. This critical charge, commonly called Rayleigh limit and denoted as Q R , was evaluated from Monte Carlo (MC) simulations performed on uniformly charged clusters for sizes ranging from N = 55 to N = 1000 at temperatures close to the melting point of the neutral clusters. Uniform charge distribution was achieved by attaching the same reduced charge to each particle, then producing a set of N identical particles. The Rayleigh limit was estimated by tracking the change in quadrupole moments as a function of the charge carried by the clusters, noting that elongation is related to prolate, that is quadrupole, deformations. The qualitative trend of MC data has been properly modeled by defining the reduced charge Q * LDM as

Q * LDM = Q LDM √ 4πǫ 0 σ LJ V LJ = 48 π 1/6 AN + BN 2/3 (2)
where ǫ 0 is the permittivity of free space, σ LJ is the LJ pair diameter, and Q LDM is the estimate of the Rayleigh limit Q R , the 'LDM' subscript reminding that E s is the surface energy used by Echt et al. in their LDM. However, Q * LDM overestimated MC data and the constant parameters A and B had to be divided by a factor of about 2 to reach more quantitative agreement with bulk properties expected for a set of eleven atomic and molecular droplets. This recent work was focused on large clusters and, consequently, no evaluation of the suitability of Eq.( 2) to account for the values of critical sizes of multiply charged clusters had been undertaken. The present letter aims to fill this gap while amending Eq.( 2) to free it from the division by an empirical constant factor. The values of bulk surface tensions are first compared to available data before discussing the critical sizes of rare-gas (Rg) clusters and 26 molecular clusters.

Method

Eq.( 2) is based on Rayleigh's formula that relates droplet instability to its ability to overcome the cohesion at surface [START_REF] Rayleigh | On the Equilibrium of Liquid Conducting Masses charged with Electricity[END_REF]

Q R = π 8ǫ 0 γD 3 ( 3
)
where γ is the surface tension and D is the droplet diameter. γ can be identified to the surface energy by unit area on the basis of the expression of E s reported in Eq.( 1) [START_REF] Bonhommeau | Rayleigh limit and fragmentation of multiply charged Lennard-Jones clusters: Can charged clusters provide clues to investigate the stability of electrospray droplets?[END_REF]. Although originally devised for argon microclusters by reckoning the difference between the energy of clusters containing N atoms and the energy of N atoms in the bulk phase [START_REF] Briant | Molecular dynamics study of the structure and thermodynamic properties of argon microclusters[END_REF], this expression of surface energy was proved to be well suited to model the stability of uniformly charged continuous media, as done by Echt et al. in their LDM aiming to confirm the critical sizes deduced from experimental mass spectra. However, the fissibility parameter X = E(Coulomb)/2E(surface) [START_REF] Last | Fragmentation channels of large multicharged clusters[END_REF] that would arise from the ratio between the Coulomb energy of MC calculations performed in free space [START_REF] Bonhommeau | Rayleigh limit and fragmentation of multiply charged Lennard-Jones clusters: Can charged clusters provide clues to investigate the stability of electrospray droplets?[END_REF] and 2E s would underestimate the fissibility parameter representative of the LDM proposed by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF] since the Coulomb energy is decreased by the dielectric constant ǫ in continuous media whereas the surface energy has the same expression in both cases. To counterbalance this underestimation of the fissibility parameter, E s is therefore divided by the dielectric constant ǫ. Moreover, droplets are assumed spherical within Rayleigh's approximation and their radii can be calculated from R = (3N/4π) 1/3 σ LJ in first approximation. Axially symmetric cluster deformations and temperature dependence of cluster radii can be taken into account by analyzing parallel tempering simulations performed on neutral LJ clusters [START_REF] Bonhommeau | Charge localization in multiply charged clusters and their electrical properties: Some insights into electrospray droplets[END_REF][START_REF] Bonhommeau | Rayleigh limit and fragmentation of multiply charged Lennard-Jones clusters: Can charged clusters provide clues to investigate the stability of electrospray droplets?[END_REF] composed of 55, 147, 309, and 1000 particles. The average cluster radii are evaluated from the expression of the ellipsoidal volume by determining the location of surface particles with a geometric model [START_REF] Miller | Structure and stability of charged clusters[END_REF], diagonalizing the inertia matrix of the cluster, and identifying the corresponding semi-principal axes to the maximum values taken by the projected surface particle coordinates along the inertia axes. As illustrated in Figure S1 of the supplementary material, the reduced radii R * = R/σ LJ of the four selected clusters roughly evolve according to the same linear equation, namely R * (T * ) = (3N/4π) 1/3 (α R T * + β R ), at temperatures above the melting point. The dimensionless parameters α R and β R hardly rely on N and averages over the four cluster sizes yield α R ≈ 0.526 and β R ≈ 0.811.

The previous approximations on surface energy and cluster radii then lead to the following expressions for the reduced surface tension γ * LDM,ǫ and the reduced critical charge

Q * LDM,ǫ γ * LDM,ǫ ≈ E * s 4πR * (T * ) 2 ǫ = 1 (36π) 1/3 ǫ A + BN -1/3 (α R T * + β R ) 2 (4a) Q * LDM,ǫ = Q LDM,ǫ √ 4πǫ 0 σ LJ V LJ = 48 π 1/6 (AN + BN 2/3 )(α R T * + β R ) ǫ (4b)
where Q LDM,ǫ = Q LDM / √ ǫ and ǫ is the dielectric constant calculated from the Clausius-Mossotti (CM) formula on the basis of the polarizability [START_REF]NIST Computational Chemistry Comparison and Benchmark DataBase[END_REF] and density [15] of the cluster constituents at temperatures close to the melting and boiling points (the complete set of temperatures, densities, and polarizabilities used in this work are supplied in Table S1 of the supplementary material). When data were not available in the NIST databases, they were taken from the work by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF] that mainly focused on physical properties near the boiling point. Note that the CM formula is only valid for gases and nonpolar liquids, and the treatment of polar liquids would require a description of dipole moments beyond the Langevin-Debye equation [START_REF] Debye | Some results of kinetic theory of isolators[END_REF] whose applicability is restricted to extremely dilute polar molecules in nonpolar solvents. As an example, more sophisticated approaches, like Kirkwood's model [START_REF] Kirkwood | The Dielectric Polarization of Polar Liquids[END_REF] that improves Onsager's theory [START_REF] Onsager | Electric Moments of Molecules in Liquids[END_REF] by taking into account the effect of dipole associations, should be considered for liquid clusters composed of polar molecules [START_REF] Wilson | The dielectric constants of polar liquids[END_REF]. However, the inclusion of such effect would yield unrealistic critical sizes since the expression of E s is not expected to properly take into account hydrogen bonding and dipole-dipole interactions. The omission of dipole effects in surface energies and dielectric constants will be shown to yield satisfactory predictions for the critical sizes of polar clusters with moderate permanent dipole moments (µ 0.5 D). Finally, reminding that γ * = γσ 2 LJ /V LJ and T * = kT /V LJ where k is the Boltzmann constant, the surface tension γ LDM,ǫ and the critical charge Q LDM,ǫ can be written as

γ LDM,ǫ = 1 (36π) 1/3 ǫσ 2 LJ A + BN -1/3 (α R kT /V LJ + β R ) 2 V LJ (5a) Q LDM,ǫ = 3072 1/6 π 1/3 ǫ 0 σ LJ ǫ (AN + BN 2/3 )(α R kT + β R V LJ ) (5b)
These quantities might depend on the particular LJ parameterization although an increased pair diameter may be compensated by a shallower well. For Rg atoms, we selected LJ parameters commonly used to model liquids or solids [START_REF] Herrero | Isotope effects in structural and thermodynamic properties of solid neon[END_REF][START_REF] Ramirez | Quantum path-integral study of the phase diagram and isotope effects of neon[END_REF][START_REF] Rowley | Monte Carlo Grand Canonical Ensemble Calculation in a Gas-Liquid Transition Region for 12-6 Argon[END_REF][START_REF] Ohtori | Explicit expression for the Stokes-Einstein relation for pure Lennard-Jones liquids[END_REF]. A particular emphasis is put on multiply charged neon clusters whose critical charges were poorly predicted by LDMs [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF][START_REF] Nakamura | Shell effects on stability and fragmentation of multiply charged rare gas clusters[END_REF] and more sophisticated approaches such as path integral molecular dynamics simulations [START_REF] Calvo | Does Vibrational Delocalization Stabilize Multiply-Charged Neon Clusters?[END_REF].

For instance, LJ parameters based on a neon pair potential employed to model the photodissociation dynamics of hydrogen halides in floppy neon clusters [START_REF] Slavicek | Pickup and Photodissociation of Hydrogen Halides in Floppy Neon Clusters[END_REF] are considered to evaluate the influence of zero-point energy (ZPE) on the prediction of critical sizes. A similar approach, based on an effective He-He potential built to perform the molecular dynamics of Ne + 4 He 100 [START_REF] Bonhommeau | Fragmentation of ionized doped helium clusters: Theoretical evidence for a dopant ejection mechanism[END_REF], is considered to estimate the critical sizes of multiply charged helium clusters although no experimental data is available for these systems. For molecular van der Waals clusters, we followed the procedure advocated by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF] and approximate LJ parameters are deduced from the relations

V LJ = 0.77 T c (6a) σ LJ = 2.44 T c p c 1/3 (6b)
where T c (in K) and p c (in atm) are the critical temperature and pressure extracted from the NIST Chemistry WebBook [15]. Values of T c and p c are reported in Table S1 of the supplementary material and the corresponding LJ parameters are presented in Tables 1,2, and 3 for Rgs, nonpolar molecules, and polar molecules, respectively. For H 2 , another set of LJ parameters including ZPE effects is also considered that makes use of the work by Sterling et al. [START_REF] Sterling | Simulations of quantum crystals by classical dynamics[END_REF] that built an effective potential for solid H 2 based on the Silvera-Goldman pair potential [START_REF] Silvera | The isotropic intermolecular potential for H 2 and D 2 in the solid and gas phases[END_REF]. The appearance sizes n c (z) are finally extracted from Eq.(5b) by determining the cluster size where Q LDM,ǫ exceeds z elementary charges. It is worth noting that this approach somewhat differs from mass spectrometry measurements where the Rayleigh limits of multiply charged clusters rather manifests themselves as intervals of instability extending from n c (z) to n st (z), the size at which all the clusters with z charges are effectively stable [START_REF] Chen | Coulomb Fission in Dielectric Dication Clusters: Experiment and Theory on Steps That May Underpin the Electrospray Mechanism[END_REF].

Bulk surface tensions

The reliability of Eq.(5a) can be ascertained by comparing its predictions as N tends to infinity with bulk surface tensions derived from experimental data or empirical models. In particular, Ferguson and Kennedy [START_REF] Ferguson | Free and total surface energies and related quantities[END_REF] proposed the following equation, sometimes called Ferguson's equation (FE), to evaluate bulk surface tensions of cryogenic liquids

γ = γ 0 1 - T T c p ( 7
)
where T c (in K) is the critical temperature, γ 0 is a parameter characteristic of the liquid, and p is a parameter deemed as universal by Ferguson and Kennedy. In the present work, FE is applied to 4 rare gases (He, Ne, Ar, and Xe), three homonuclear diatomic molecules (H 2 , N 2 , and O 2 ), carbon monoxide and methane [START_REF] Sprow | Surface Tensions of Simple Liquids[END_REF][START_REF] Smith | Surface Tension and Energy of Liquid Xenon[END_REF] (the corresponding parameters are collected in Table S2 of the supplementary material). Note that the NIST values of T c slightly differ from those advocated for the adjustment of FE, the larger discrepancy arising for O 2 where the critical temperature in FE is overestimated by about 3% [START_REF] Sprow | Surface Tensions of Simple Liquids[END_REF]. These data are supplemented by bulk surface tensions calculated from Fuks and Bellemans' linear equation for krypton [START_REF] Fuks | The surface tension of krypton, methane and their mixtures[END_REF], namely γ Kr (mN/m) = 16.344 -0.2089(T -116) [START_REF] Briant | Molecular dynamics study of the structure and thermodynamic properties of argon microclusters[END_REF] and experimental data from the literature for water [START_REF]CRC Handbook of Chemistry and Physics, Internet Version 2005[END_REF]. Bulk surface tensions of these 11 species are compared to predictions from Eq.(5a) in Figure 1. The qualitative trends, like the steady increase of surface tension as the mass of Rgs increases and the lower surface tension at boiling, are preserved by Eq.(5a). Apart from helium, neon and water, the deviation between predicted surface tensions and reference data does not exceed ∼ 10% at the boiling point, the agreement being particularly satisfactory for the heavier Rgs and diatomic molecules. At the melting point, the relative error rises to about 20-30% for molecular clusters but remains low (< 10%) for Rgs with classical behaviors. However, it is worth noting that these deviations to expectations will be damped in expression (5b) of the critical charge that relies on the square root of surface tension. Considering an alternative definition for the cluster radius as the mean distance between surface particles and the cluster center of mass is found to worsen the agreement with the reference data.

Critical sizes of rare-gas clusters

Critical sizes n c (z) (z = 2 -4) of Rg clusters predicted from Eq.(5b) are collected in Table 1 and compared to available experimental and theoretical data. Eq.(5b) yields n c (2) = 278 and n c (3) = 620 at the melting point when using LJ parameters employed in quantum path-integral Monte Carlo simulations to investigate the isotopic effects in structural and thermodynamic properties of neon at low temperature (T ≤ 24 K) [START_REF] Herrero | Isotope effects in structural and thermodynamic properties of solid neon[END_REF], in close agreement with the experimental values of Mähr et al. [START_REF] Mähr | Multiply Charged Neon Clusters: Failure of the Liquid Drop Model?[END_REF], namely n c (2) = 284 and n c (3) = 656. Considering neon properties at temperatures close to the boiling point [START_REF]The melting and boiling points of rare gases are very close to each other at p=0.1 MPa and significant temperature effects between these two limits are unlikely[END_REF] or other LJ parameterizations specific to solid or liquid neon [START_REF] Ramirez | Quantum path-integral study of the phase diagram and isotope effects of neon[END_REF] do not significantly alter the results, deviations to experiments never exceeding 12%. Interestingly, the microscopic attempts to determine the first critical sizes of neon clusters improved the results of Echt et al., that is n c (2) = 868 and n c (3) = 2950, but they failed to confirm the experimental data. Nakamura described the total energy of charged Rg clusters as a neutral and a Coulomb contribution [START_REF] Nakamura | Shell effects on stability and fragmentation of multiply charged rare gas clusters[END_REF], the neutral term being either represented by the clusters global minimum energy or by a power law in N stemming from LDMs. Despite the implicit inclusion of shell effects in global minimum energy, the best estimate of n c (2) was 635 at T = 0 K, more than twice the expected value. Calvo undertook to evaluate the influence of vibrational delocalization of neon atoms by performing path integral molecular dynamics simulations [START_REF] Calvo | Does Vibrational Delocalization Stabilize Multiply-Charged Neon Clusters?[END_REF] but the quantum result at the melting point (n c (2) = 626) only slightly improved the classical result (n c (2) = 703) then making the author conclude that a full quantum modeling of both electrons and nuclei should presumably be required to get accurate results. This weak influence of vibrational delocalization of neon atoms is also ascertained from Eq.(5b) when considering LJ parameters deduced from an effective Ne-Ne interatomic potential that includes ZPE effects [START_REF] Slavicek | Pickup and Photodissociation of Hydrogen Halides in Floppy Neon Clusters[END_REF]. Critical sizes remain in the vicinity of experimental values (e.g., n c (2) = 274 and n c (3) = 610 close to the melting point), underestimating them by less than ∼ 10%.

In contrast, vibrational delocalization might have a significant effect on the critical sizes obtained for helium clusters whose effective He-He potential only retains 20% of the He 2 dissociation energy [START_REF] Janzen | An accurate potential energy curve for helium based on ab initio calculations[END_REF][START_REF] Bonhommeau | Fragmentation of ionized doped helium clusters: Theoretical evidence for a dopant ejection mechanism[END_REF], to be compared with the 82% specific to the aforementioned effective Ne-Ne potential [START_REF] Aziz | The Ne-Ne interatomic potential revisited[END_REF][START_REF] Slavicek | Pickup and Photodissociation of Hydrogen Halides in Floppy Neon Clusters[END_REF]. For instance, the neglect of ZPE in the He-He interatomic potential to model the dynamics of doped helium clusters is well known to significantly overestimate the cage effect exerted on the dopant [START_REF] Bonhommeau | Fragmentation of ionized doped helium clusters: Theoretical evidence for a dopant ejection mechanism[END_REF]. The appearance size of helium dications making use of such an effective potential only reaches 2563 at temperatures close to the melting point, to be compared with the value of 17092 predicted by the LDM of Echt et al. (LDME). Deviations between the two models still increase for the higher order critical sizes n c (3) and n c (4). The result of Eq.(5b) should however be considered with caution since the effective He-He potential was significantly harder to converge than the Ne-Ne potential due to the small binding energy of He 2 (∼ 1.3 × 10 -3 K [START_REF] Toennies | Spectroscopy of atoms and molecules in liquid helium[END_REF]). Thus, the wavefunctions centered on He atoms had to be truncated to 1.6 Å during the iterative process leading to convergence [START_REF] Bonhommeau | Modeling the fragmentation dynamics of ionic clusters inside helium nanodroplets: the case of He 100 Ne + 4[END_REF]. This artifice was not required to generate the effective Ne-Ne potential that converged after as few as 4 iterations [START_REF] Slavicek | Pickup and Photodissociation of Hydrogen Halides in Floppy Neon Clusters[END_REF]. Consequently, the "true" effective He-He potential of helium clusters might be more weakly bound and the results provided for helium in Table 1 should be regarded as lower limits of the real critical sizes. Despite the lack of information on multiply charged helium clusters, critical sizes can also be estimated from the knowledge of helium bulk surface tension derived from FE (see Eq.( 7) and Table S2 of the supplementary material) by expressing the Rayleigh limit reported in Eq.(3) as

Q bulk R = 4 3πǫ 0 γM N A ρ N (9)
where N A is Avogadro's number, N is the cluster size, ρ is the density of helium, and M is the molar mass of helium. The corresponding critical sizes are significantly larger than that predicted by Eq.(5b) since n c (2) = 5744 at the melting temperature and n c (2) = 15670 at the boiling temperature (see Table S3 of the supplementary material), the latter value being in good agreement with the prediction of the LDME. Knowing that the replacement of surface tension by its bulk value is generally valid for clusters composed of thousands of atoms or molecules we might expect Eq.( 9) to provide a suitable estimate of helium critical sizes although we cannot exclude more exotic behaviors related to the highly quantum nature of helium clusters.

The relative failure of microscopic simulations to properly estimate neon critical sizes has to be tempered since excellent agreement has been obtained by Calvo for argon clusters [START_REF] Calvo | Does Vibrational Delocalization Stabilize Multiply-Charged Neon Clusters?[END_REF], especially at T = 0 K where classical and quantum simulations yielded n c (2) = 93 and n c (2) = 94, respectively. The application of Eq.(5b) leads to a critical size for argon dications of similar accuracy (n c (2) = 89), in much better agreement with experimental results than the LDME (n c (2) = 122). Although the quality of the prediction drops for n c (3) (deviation of 15% at the boiling temperature), it remains similar to that of Nakamura's calculations [START_REF] Nakamura | Shell effects on stability and fragmentation of multiply charged rare gas clusters[END_REF]. For krypton and xenon clusters, the critical sizes predicted from Eq.(5b) for LJ parameters characteristic to liquids [START_REF] Ohtori | Explicit expression for the Stokes-Einstein relation for pure Lennard-Jones liquids[END_REF] slightly underestimate the experimental values. In particular, n c (2) = 66 and n c (3) = 145 for liquid krypton clusters, these two values underestimating by ∼ 7% the experimental appearance sizes. On the contrary the LDME leads to the correct value for n c (2) but it overestimates n c (3) by 13%. This approach also provides results of same accuracy as Eq.(5b) for xenon critical sizes when z = 2 and z = 3.

Therefore, Eq.(5b) yields critical sizes of similar or higher accuracy than previous LDMs and microscopic simulations for Rg clusters. Apart from helium whose treatment is made involved by the lack of experimental data, results at the melting point are in better agreement with experimental expectations than those at the boiling point. The predictions of Eq.(5b) never deviate from the experimental appearance sizes by more than 15% and the minor influence of ZPE on the predicted neon critical sizes indicate that the lower polarizability of neon atoms should play a major role in accounting for the larger critical sizes of multiply charged neon clusters. Rg clusters should be ideal systems to evaluate the applicability of Eq.(5b) due to their zero permanent dipole moment and the relative adequacy of LJ pair potentials to describe their interatomic interactions. The differences with experimental data may be essentially attributed to inaccuracies in the definition of surface tension and to the assumption of uniform charge distribution.

Extension to molecular clusters

Nonpolar molecules

The critical sizes n c (z) (z = 2 -4) of nonpolar molecular clusters are presented in Table 2. The first striking point is the noticeable effectiveness of Eq.(5b) to predict the experimental appearance sizes of N 2 (n c (2) = 99 and n c (3) = 215), methane (n c (2) = 70), and benzene (n c (2) = 22 -26, n c (3) = 52, and n c (4) = 92) when the results are considered in the vicinity of the melting point. The excellent agreement obtained for N 2 clusters contrasts with the apparent failure of the LDME that yields n c (2) = 145 for these clusters and similar disagreements when z > 2. The experimental critical sizes of other small nonpolar molecules like O 2 (n c (2) = 92) and ethylene (n c (2) = 51) lie in between the values obtained at the melting and boiling temperatures with Eq.(5b), namely n c (2) = 82 -99 for O 2 , and n c (2) = 45 -56 for ethylene with alike accuracies obtained for z ≥ 3. Predictions from Eq.(5b) and LDME provide the same result for SF 6 (n c (2) = 44) in reasonable agreement with the experimental value (n c (2) = 39). Two LJ well depths are considered for CO 2 , V LJ = 234.2 K is based on Eq.(6a) whereas V LJ = 290 K was adjusted by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF] to improve the agreement between their LDM and experimental results, noting that this LJ well depth allowed to better fit the bulk sublimation energy of CO 2 . The latter parameterization is found to improve the quality of both theoretical models. As an example, the critical size of CO 2 dications drops from n c (2) = 61 to n c (2) = 48 when Eq.(5b) is employed, to be compared with the experimental value of n c (2) = 43 -45. Apart from n-butane (n c (2) = 24 -32 and n c (3) = 51 -70), the agreement with experimental data seems to decrease for large hydrocarbons like 2,2-Dimethylpropane and (E)-2-Butene. Beside the possible influence of the uniform charge distribution and surface tension definition discussed for Rgs, this discrepancy might also stem from the spherical symmetry of particles used in benchmark MC simulations [START_REF] Bonhommeau | MCMC 2 : A Monte Carlo code for multiply-charged clusters[END_REF], spherical symmetry becoming less and less representative of linear systems as their carbon chain lengthens. Note that the neglect of the temperature dependence in the expression of R * (T * ) (see Eq.(4a)) is found to improve the quality of the predictions for the larger hydrocarbons, but the relative error between experimental and theoretical (see Eq.(5a)) bulk surface tensions is increased and the predicted appearance sizes of small molecules like N 2 , O 2 , and CH 4 are less accurate.

For H 2 , the spectacular disagreement between the LDME (n c (2) = 863) and Eq.(5b) (n c (2) = 317 -360) reminds the results obtained for neon and helium clusters. Replacing H 2 LJ parameters derived from Eqs.(6a) and (6b) by parameters including ZPE effects does not drastically modify the critical sizes predicted by Eq.(5b) (n c (2) = 314 -364) although more than one third of the dissociation energy is lost compared to the well depth of the Silvera-Goldman pair potential. In close analogy with the results obtained for helium clusters, critical sizes of dications predicted from Eq.( 9) at the boiling temperature (n c (2) = 846) are in good agreement with the LDME when z = 2 (see Table S3 of the supplementary material). However, the lack of experimental data for H 2 forbids the drawing of any unambiguous conclusion about the relative effectiveness of the different approaches.

Polar molecules

Strictly speaking, LDMs and Eq.(5b) should not be appropriate to model the critical sizes of polar liquids since the dielectric constant is calculated from the Clausius-Mossotti relationship. However, we may wonder if the experimental order of magnitude is reachable under some conditions. Critical sizes of polar liquids with permanent electric dipole moments ranging from µ(CO) = 0.112 D to µ(CH 3 CN) = 3.92 D are reported in Table 3. For small molecules with µ < 0. 5b) is more subject to debate for NO since no reference data is available for this molecule. Moreover, the critical sizes predicted for larger dications like those containing (Z)-2-Butene (n c (2) = 37 and n c (3) = 82) seem less accurate, reminding the results obtained with (E)-2-Butene, and their reliability decreases for n c (3) whereas the LDME provides more reliable results for n c (3) than n c (2). The critical sizes of toluene (µ = 0.332 D) resemble those of benzene as expected from the similarities between these two molecules. More interestingly, the results for the highly polar pyridine (µ = 2.19 D and n c (2) = 23) compares well with experimental sizes (n c (2) = 12-27) although this agreement may be fortuitous because of the broad range of stability of pyridine. Satisfactory predictions are also reached for some highly polar molecules composed of fluorine, namely CHF 3 (n c (2) = 34 -41, n c (3) = 74 -91, and n c (4) = 130 -160) and CH 3 CHF 2 (n c (2) = 26 -32, n c (3) = 57 -71). This unexpected behavior cannot be generalized to other highly polar species since Eq.(5b) and the LDME both fail to predict the appearance sizes of water and ammonia. Furthermore, the LDME provides better predictions than Eq.(5b) for H 2 S (n c (2) = 37 and n c (3) = 85) and it is the only model that estimates properly the critical sizes of CH 3 OCH 3 (n c (2) = 31 and n c (3) = 77) although it fails to supply reliable results for CHF 3 . The failure of Eq.(5b) is also prominent for acetonitrile since this equation predicts n c (2) = 20 whereas the experimental interval of stability for acetonitrile dications extends from 28 to 32 molecules [START_REF] Chen | Coulomb Fission in Dielectric Dication Clusters: Experiment and Theory on Steps That May Underpin the Electrospray Mechanism[END_REF]. Surprisingly, the critical sizes derived from Eq.( 9) by replacing the surface tension by its bulk value leads to the correct critical sizes for water (n c (2) = 34 -42) and acetonitrile (n c (2) = 32 at T = 20 • C) whereas surface tension should not match its bulk value for clusters as small as 30-40 molecules (see Table S3 of the supplementary material). However, Eq.( 9) cannot be blindly recommended since it provides poor predictions for neon (n c (2) = 505 -573) and argon (n c (2) = 123 -130) clusters. Consequently, neither the LDME nor Eq.(5b) supply reliable results for all the polar species, as expected from the neglect of polarization in both models, and compensation of errors probably account for the good agreement obtained for some highly polar molecules.

Conclusion

In this letter, a simple law depending on cluster size N has been proposed to evaluate the critical charge, and subsequently the appearance sizes, above which 31 atomic and molecular clusters become unstable. Three parameters are only required to get an estimate of this charge: the dielectric constant of the medium, the LJ pair diameter of the species, and its LJ pair well depth. The temperature-dependence of cluster radii has been adjusted on MC data above the melting point of neutral LJ clusters to improve the connection of the model to bulk properties, chief among which is surface tension. The critical sizes n c (z) (z = 2 -4) of multiply charged clusters are satisfactorily predicted as confirmed by their comparison with experimental data whose reliability is supported by their reproducibility for a broad range of van der Waals and hydrogen-bonded clusters [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF][START_REF] Gotts | Observations of Coulomb explosion in doubly charged atomic and molecular clusters[END_REF][START_REF] Chen | Coulomb Fission in Dielectric Dication Clusters: Experiment and Theory on Steps That May Underpin the Electrospray Mechanism[END_REF]. In particular, a remarkable agreement is obtained for multiply charged neon clusters with several LJ parameterizations and the results are generally of similar or higher accuracy than those provided by LDMs for nonpolar molecular systems. Although Eq.(5b) yield proper estimates of appearance sizes for most polar molecules with moderate permanent dipole moments (µ = 0.5 D) investigated in the present work, the excellent agreement obtained for some highly polar molecules can be presumably attributed to compensation of errors. Despite its effectiveness, the current model assumes uniform charge distribution, neglects polarization effects, and the dependence on temperature is only taken into account through the CM dielectric constants and approximate radii of deformed LJ clusters. The present work is also focused on simple multiply charged clusters composed of a unique type of atoms or molecules interacting with each other by LJ pair potentials and Coulomb interactions. Extensions to multicomponent clusters, widely used in electrospray ionisation experiments, is also posed knowing that LJ parameters can be built by means of combining rules.

Supplementary material

Estimated radii for LJ clusters composed of 55, 147, 309 and 1000 particles are depicted in Figure S1 (Section S1). General properties of Rg and molecular clusters, including the molar masses, the densities at the melting and boiling points, the critical pressures and temperatures, and the polarizabilities of the cluster constituents are listed in Table S1 (Section S2). Section S3 collects bulk surface tensions of some cryogenic liquids (Table S2) and critical sizes for these liquids and some electrospray compatible solvents (water and acetonitrile). TABLE 1: Critical sizes n c (z) (z = 2 -4) of multiply charged Rg clusters at temperatures close to the melting and boiling points. Values at the boiling point are written in parentheses. Unless otherwise specified, benchmark theoretical results (columns "Calc.") and experimental data (columns "Exp.") refer to the LDM and the experimental results by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF]. TABLE 1: Critical sizes n c (z) (z = 2-4) of multiply charged nonpolar molecular clusters at temperatures close to the melting and boiling points. Values at the boiling point are written in parentheses. Unless otherwise specified, benchmark theoretical results (columns "Calc.") and experimental data (columns "Exp.") refer to the LDM and the experimental results by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF]. TABLE 3: Critical sizes n c (z) (z = 2 -4) of multiply charged polar molecular clusters at temperatures close to the melting and boiling points. Values at the boiling point are written in parentheses. Unless otherwise specified, benchmark theoretical results (columns "Calc.") and experimental data (columns "Exp.") refer to the LDM and the experimental results by Echt et al. [START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF].

FIG. 1: Comparison between bulk surface tensions from reference data [START_REF] Sprow | Surface Tensions of Simple Liquids[END_REF][START_REF] Smith | Surface Tension and Energy of Liquid Xenon[END_REF][START_REF] Fuks | The surface tension of krypton, methane and their mixtures[END_REF][START_REF]CRC Handbook of Chemistry and Physics, Internet Version 2005[END_REF] and Eq.(5a). (a) Comparison at the melting point. (b) Comparison at the boiling point. (c) Relative errors at the melting and boiling points. For He, Ne, and H 2 , LJ parameters that include ZPE effects [START_REF] Bonhommeau | Fragmentation of ionized doped helium clusters: Theoretical evidence for a dopant ejection mechanism[END_REF][START_REF] Slavicek | Pickup and Photodissociation of Hydrogen Halides in Floppy Neon Clusters[END_REF][START_REF] Sterling | Simulations of quantum crystals by classical dynamics[END_REF] have been considered.

a In each cited reference the value with best agreement to experimental data has been selected disregarding the temperature of the calculation. a In each cited reference the value with best agreement to experimental data has been selected disregarding the temperature of the calculation. b For acetonitrile, the dielectric constant is given at 20 • C instead of the melting temperature. c For pyridine, the dielectric constant is given at 20 • C, in between the melting and boiling points.

  5 D like CO (n c (2) = 85 -93) and N 2 O (n c (2) = 45, n c (3) = 98 -99, and n c (4) = 172 -173), a reasonable agreement with experimental results is obtained, the LDME leading to better predictions for N 2 O but more dubious for CO. Critical sizes of CF 3 Cl (n c (2) = 33 -45 and n c (3) = 73 -98) and propene (n c (2) = 31 -44 and n c (3) = 69 -97) are also in good agreement with experimental expectations. The relative effectiveness of the LDME and Eq.(

Figure 1 :

 1 Figure 1: Comparison between bulk surface tensions from reference data [32, 33, 34, 35] and Eq.(5a). (a) Comparison at the melting point. (b) Comparison at the boiling point. (c) Relative errors at the melting and boiling points. For He, Ne, and H 2 , LJ parameters that include ZPE effects [27, 26, 28] have been considered.

Table 1 :

 1 Critical sizes nc(z) (z = 2 -4) of multiply charged Rg clusters at temperatures close to the melting and boiling points. Values at the boiling point are written in parentheses. Unless otherwise specified, benchmark theoretical results (columns "Calc.") and experimental data (columns "Exp.") refer to the LDM and the experimental results by Echt et al.[START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF].

	Species	LJ parameters	ǫ		nc(2)			nc(3)			nc(4)
		σ LJ ( Å) V LJ (K)		This work	Calc. a	Exp.	This work	Calc. a	Exp.	This work	Calc. a	Exp.
	He[27]	3.754	2.209	1.059	2563	17092	-	5740	76895	-	10178	198966	-
				(1.050)	(1865)			(4174)			(7398)
	Ne[20]	2.782	35.79	1.190	278	868, 626[25], 635[24]	284[36]	620	2950, 2259[24] 656[36]	1095	6424	-
				(1.183)	(268)			(598)			(1056)
	Ne[21]	2.782	37.29	1.190	271			602			1065
				(1.183)	(261)			(581)			(1027)
	Ne[26]	2.895	34.53	1.190	274			610			1079
				(1.183)	(264)			(588)			(1039)
	Ar[22]	3.405	119.8	1.525	89	122, 92[24], 93[25]	91, 89[43]	196	333, 251[24]	226	345	648	-
				(1.516)	(87)			(192)			(339)
	Kr[23]	3.670	167.0	1.676	66	71, 54[24]	71	145	177, 135[24]	156	254	331	-
				(1.667)	(65)			(142)			(250)
	Xe[23]	3.924	257.4	1.918	47	46, 40[24]	51	104	107, 87[24]	114	183	196, 151[24]	208
				(1.908)	(47)			(103)

Table 2 :

 2 Critical sizes nc(z) (z = 2 -4) of multiply charged nonpolar molecular clusters at temperatures close to the melting and boiling points. Values at the boiling point are written in parentheses. Unless otherwise specified, benchmark theoretical results (columns "Calc.") and experimental data (columns "Exp.") refer to the LDM and the experimental results by Echt et al.[START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF].

	Species	LJ parameters	ǫ		nc(2)		nc(3)			nc(4)
		σ LJ ( Å) V LJ (K)		This work Calc. a	Exp.	This work Calc. a	Exp. This work Calc. a	Exp.
	H 2	3.349	25.55	1.246	360	863	-	803	2824	-	1420	6014	-
				(1.225)	(317)			(707)			(1249)
	H 2 [28]	3.891	20.46	1.246	364			810			1433
				(1.225)	(314)			(700)			(1238)
	N 2	3.795	97.17	1.462	96	145	99	212	411	215	374	814	-
				(1.425)	(88)			(194)			(342)
	O 2	3.560	119.0	1.575	99	120	92	218	330	-	384	646	-
				(1.491)	(82)			(180)			(317)
	CO 2	3.929	234.2	(1.611)	(61)	54	45, 43[43]	(134)	143	109	(236)	274	216
		4.100	290.0	(1.611)	(48)	43	45, 43[43]	(105)	112	109	(184)	216	216
	SF 6	4.998	245.4	(1.522)	(44)	44	39	(96)	121	-	(169)	234	-
	CH 4	3.933	146.8	1.631	70	78	70[44]	154	200	-	272	379	-
				(1.583)	(64)			(141)			(248)
	C 2 H 4	4.348	217.5	1.982	56	47	51	123	115	108	216	214	192
				(1.816)	(45)			(99)			(173)
	CH 3 CH 2 CH 2 CH 3	5.481	327.3	2.031	32	26	31	70	63	78	122	116	-
				(1.795)	(24)			(51)			(90)
	E-CH 3 CHCHCH 3	5.359	330.0	(1.855)	(37)	26	31	(82)	62	66	(143)	114	-
	C(CH 3 ) 4	5.844	334.0	1.870	23	24	30	50	58	-	88	108	-
				(1.823)	(22)			(48)			(83)
	C 6 H 6	5.531	432.7	2.210	24	18	23[7, 43], 22-26[30]	51	38	52	90	69	92
				(2.064)	(21)			(45)			(78)

Table 3 :

 3 Critical sizes nc(z) (z = 2 -4) of multiply charged polar molecular clusters at temperatures close to the melting and boiling points. Values at the boiling point are written in parentheses. Unless otherwise specified, benchmark theoretical results (columns "Calc.") and experimental data (columns "Exp.") refer to the LDM and the experimental results by Echt et al.[START_REF] Echt | Dissociation channels of multiply charged van der Waals clusters[END_REF].

	Species	µ a (D)	LJ parameters	ǫ		nc(2)		nc(3)			nc(4)
			σ LJ ( Å) V LJ (K)		This work	Calc.	Exp.	This work	Calc.	Exp. This work Calc.	Exp.
	CO	0.112	3.839	103.5	1.527	93	131	98	205	364	-	361	715	-
					(1.487)	(85)			(189)			(332)
	NO	0.153	3.440	138.8	(1.684)	(121)	83	-	(267)	195	-	(471)	358	-
	N 2 O	0.167	3.978	238.4	1.810	45	44	51	99	104	105	173	191	184
					(1.805)	(45)			(98)			(172)
	H 2 O	1.855	3.508	498.2	1.798	27	28	35, 36[43], 35-40[30]	60	69	-	104	130	-
					(1.757)	(25)			(53)			(93)
	H 2 S	0.970	3.942	287.5	2.091	46	37	39	100	85	84	176	156	-
					(2.028)	(43)			(94)			(164)
	NH 3	1.470	3.752	312.2	1.888	41	38	51, 50-55[30]	89	91	121	156	169	-
					(1.810)	(37)			(80)			(141)
	CH 3 CN b	3.920	5.482	419.7	1.777	20	-	28-32[30]	42	-	-	74	-	-
					(1.690)	(18)			(38)			(66)
	CHF 3	1.650	4.505	230.3	1.622	41	52	41	91	142	90	160	277	168
					(1.513)	(34)			(74)			(130)
	CF 3 Cl	0.500	4.854	232.4	1.794	45	44	43	98	114	87	172	218	-
					(1.619)	(33)			(73)			(128)
	CH 3 CHF 2	2.270	5.012	297.6	1.809	32	34	33	71	88	69	124	167	-
					(1.659)	(26)			(57)			(99)
	CH 3 OCH 3	1.300	4.781	308.8	(1.786)	(43)	31	35	(94)	77	72	(164)	145	-
	CH 3 CHCH 2	0.363	4.896	281.2	2.116	44	33	39	97	79	78	170	146	-
					(1.840)	(31)			(69)			(120)
	Z-CH 3 CHCHCH 3	0.257	5.340	335.3	(1.881)	(37)	25	31	(82)	60	63	(143)	109	-
	C 6 H 5 CH 3	0.332	5.971	456.6	2.390	25	-	-	55	-	-	96	-	-
					(2.017)	(18)			(38)			(67)
	C 5 H 5 N c	2.190	5.440	476.6	2.269	23	-	12-27[30]	49	-	-	86	-	-

a Permanent electric dipole moment of cluster constituents.

Acknowledgments

DAB is indebted to CNRS (Centre National de la Recherche Scientifique) for the award of a Chaire d'excellence and warmly acknowledges the supercomputing centre of Champagne-Ardenne (ROMEO) for computational resources and technical support.