Binbin Xu

Sabir Jacquir
email: sjacquir@u-bourgogne.fr

Gabriel Laurent

Jean-Marie Bilbault

Stéphane Binczak

A Hybrid Stimulation Strategy for Suppression of Spiral Waves in Cardiac Tissue

Atrial fibrillation (AF) is the most common cardiac arrhythmia whose mechanisms are thought to be mainly due to the self perpetuation of spiral waves (SW). To date, available treatment strategies (antiarrhythmic drugs, radiofrequency ablation of the substrate, electrical cardioversion) to restore and to maintain a normal sinus rhythm have limitations and are associated with AF recurrences. The aim of this study was to assess a way of suppressing SW by applying multifocal electrical stimulations in a simulated cardiac tissue using a 2D FitzHugh-Nagumo model specially convenient for AF investigations. We identified stimulation parameters for successful termination of SW. However, SW reinduction, following the electrical stimuli, leads us to develop a hybrid strategy based on sodium channel modification for the simulated tissue.

Introduction

Although human data support the idea that rotors are a crucial mechanism for fibrillation maintenance in both atria and ventricles, there are clear differences between the 2 chamber types. Ventricular fibrillation (VF), responsible for sudden cardiac death, can be considered as a lethal "three-dimensional" cardiac arrhythmia with multiple micro re-entries located deep inside the endocardial muscle [START_REF] Nanthakumar | Optical mapping of langendorff-perfused human hearts : establishing a model for the study of ventricular fibrillation in humans[END_REF]. VF has been widely studied [START_REF] Wellens | Ventricular arrhythmias : Mechanisms and actions of antiarrhythmic drugs[END_REF][START_REF] Moreno | Ionic currents and ventricular fibrillation dynamics[END_REF][START_REF] White | Ventricular fibrillation and defibrillation : contemporary understanding of mechanisms, Anaesthesia, Pain, Intensive Care and Emergency A[END_REF], but because of electrophysiological characteristics, the only treatment available is indeed an electrical defibrillation via automated external defibrillators (AED), or via implantable cardiac defibrillators (ICD). Conversely, most of "two-dimensional" macro reentry ventricular tachycardia (VT) can be stopped by rapid burst of stimulations delivered by an ICD via a single right ventricular lead [START_REF] Mirowski | The automatic implantable cardioverter-defibrillator : An overview[END_REF].

As opposed to VF, atrial fibrillation (AF) is associated with a better outcome but it is by far the most frequently diagnosed cardiac arrhythmia affecting, for example, approximately 1% of the adult population of the United States.

Recent clinical trial evidence suggests that the presence of AF is an independent predictor of morbidity and mortality. One of the primary goals of paroxysmal and persistent AF treatment is to relieve symptoms by a rhythm control strategy [START_REF] Patel | Atrial fibrillation : Pharmacological therapy[END_REF].

Endocardial walls of the atria are much thinner than their ventricular counterparts, therefore AF may be considered as a "two-dimensional" cardiac arrhythmia. Specifically, therapy for restoring and maintaining a normal sinus rhythm is multidimensional with treatment options spanning from pharmacologic therapy to electrical cardioversion and catheter-based radiofrequency ablation [START_REF] Camm | Guidelines for the management of atrial fibrillation[END_REF]. But even combined, these methods have limited efficacy as the AF recurrence rate is high, and are limited by numerous side effects [START_REF] Weerasooriya | Catheter ablation for atrial fibrillation : Are results maintained at 5 years of follow-up ?[END_REF]. To date, there are two known concepts to restore a normal sinus rhythm by using electrical stimulation. The first one is to "reset" the whole heart by applying high voltage via a dedicated defibrillation device (internal or external) in order to extinguish most of the fibrillatory wavefronts. External cardiac defibrillation cannot be rea-lized very often as it has to be done under general anesthesia. Maintenance of normal sinus rhythm using repeated internal cardioversion shocks via implanted devices has also been shown to be effective and safe, however increase in defibrillation thresholds and patient intolerance due to the pain were the main reasons for discontinuation of this therapy [START_REF] Geller | Treatment of atrial fibrillation with an implantable atrial defibrillator -long term results[END_REF]. The second way of suppressing AF is based on Allessie's concept asserting that there is a critical number of micro-reentries below which the cardiac arrhythmia can self terminate [START_REF] Allessie | Experimental evaluation of moe's multiple wavelet hypothesis of atrial fibrillation, Cardiac Electrophysiology and Arrhythmias[END_REF].

Micro-reentries can be observed experimentally or numerically as spiral waves (SW). The rotors of SW, defined as rotating motors that give rise to the spiraling wavefronts, gyrate at relatively high frequencies [START_REF] Jalife | Rotors and spiral waves in atrial fibrillation[END_REF]. From such rotors, SW propagate throughout the myocardium in very complex ways [START_REF] Moreno | Ionic currents and ventricular fibrillation dynamics[END_REF], for example SW meandering and drifting, whose dynamics have been well studied [START_REF] Fenton | Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity[END_REF][START_REF] Hakim | Spiral wave meander in excitable media : The large core limit[END_REF][START_REF] Henry | Spiral wave drift in an electric field and scroll wave instabilities[END_REF].

It is known that, under certain conditions, electrical stimuli can modify SW dynamics [START_REF] Krinsky | Autowave approaches to cessation of reentrant arrhythmias[END_REF][START_REF] Davidenko | Effects of pacing on stationary reentrant activity : theoretical and experimental study[END_REF][START_REF] Kamjoo | Importance of location and timing of electrical stimuli in terminating sustained functional reentry in isolated swine ventricular tissues : evidence in support of a small reentrant circuit[END_REF]. Moreover, they can terminate SW either by preventing rotor formation and maintenance, by re-driving / pacing SW [START_REF] Stamp | Suppressing arrhythmias in cardiac models using overdrive pacing and calcium channel blockers[END_REF] or by feedback control [START_REF] Panfilov | Elimination of spiral waves in cardiac tissue by multiple electrical shocks[END_REF]. This hypothesis could explain the global mechanism of defibrillation by electrical stimulation. Namely, at cellular level, the propagation of action potential in excitable medium is controlled by voltage-gated ion channels in cells.

These channels are shut when the membrane potential is near the resting potential. They rapidly begin to open if the membrane potential increases to the threshold value. However, as opposed to VT, AF can't be stopped by single site burst stimulations via a pacing lead [START_REF] Camm | Guidelines for the management of atrial fibrillation[END_REF], but a multifocal and simultaneous stimulation may help to extinguish most or all fibrillatory wavefronts in a selected area by resetting cellular membrane potential under its threshold value. This approach has already been investigated in previous studies. Allessie has shown that, during electrically induced AF in chronically instrumented conscious dogs, it was possible to control the local fibrillatory process by rapid focal pacing through a grid of electrodes at the left atrial appendage [START_REF] Allessie | Regional control of atrial fibrillation by rapid pacing in conscious dogs[END_REF].

However, this method failed to terminate AF as the area under control was very limited. Other authors have used an array of control points in order to terminate spatiotemporal chaos in excitable media and focused mainly on the stimulation amplitude [START_REF] Sridhar | Controlling spatiotemporal chaos in excitable media using an array of control points[END_REF]. However, we strongly believe that there are many other parameters to be taken into account.

Other methods to suppress SW based on boundaries conditions or spatial inhomogeneity [START_REF] Osipov | Controlled movement and suppression of spiral waves in excitable media[END_REF][START_REF] Zhang | Suppression of spiral waves and spatiotemporal chaos by generating target waves in excitable media[END_REF] are either difficult in controlling the inhomogeneity of a real cardiac tissue, or are associated with the generation of new waves. In [START_REF] Shajahan | Spiral-wave turbulence and its control in the presence of inhomogeneities in four mathematical models of cardiac tissue[END_REF], a spatially extended mesh which could be more successful in the elimination of spiral turbulence has been mentioned. But the optimal parameters of the controlling signals have not been defined. Moreover, these techniques have strong limitations as they not only eliminate SW, but also reinduce new ones.

In fact, the ideal atrial stimulation/defibrillation method should suppress enough SW in a dedicated region for the whole heart to stop fibrillating, but without generating new ones.

In order to address this major concern, we studied two different options by using 2D numerical simulations of the model, as described in section 2. Then, we optimized the stimulation parameters (section 3). For example, we defined the number of required stimulated sites via an electrode array ; we optimized the frequency and duty cycle of the stimulation, and the output amplitude which has to be set as low as possible to prevent from new SW and to protect the whole heart from arrhythmia re-introduction. But this minimal amplitude being still too large, we investigated then a second option, in developing a hybrid strategy based on a sodium channel modification of the simulated tissue. By increasing the threshold of new SW generation, we showed that it is possible to protect the tissue from arrhythmia re-induction. We sought to identify the best balance between cells number needed to be stimulated compared to the number of total cells simulated and the SW generation threshold. In the last section 4, we concluded this paper and gave some future prospects of this work.

Model and methods

In literature, there are many cardiac excitation models, for example, Luo-Rudy 1 and 2 models [START_REF] Luo | A model of the ventricular cardiac action potential. depolarization, repolarization, and their interaction[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes[END_REF], Ten Tusscher-Noble-Panfilov model [START_REF] Tusscher | A model for human ventricular tissue[END_REF]. These models represent an interesting way to investigate the heart activity at cellular level, but the control of their parameters is quite complex. The FitzHugh-Nagumo (FHN) model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Nagumo | An active pulse transmission line simulating nerve axon[END_REF] remains simple and keeps the essential mathematical properties of excitation and propagation from the electrochemical properties of sodium and potassium ion flows. FHN model is then used to study the ionic current of nerve or cardiac membrane [START_REF] Scott | Neuroscience : a mathematical primer[END_REF][START_REF] Kogan | Introduction to Computational Cardiology[END_REF], completed by adding terms to a non-linear "relaxation oscillators" dimensionless equation. We start from the generic ionic equation of excitable medium with realistic units,

C m ∂V ∂t = C m d∆V (x, y, t) -i N a -i other + i stim , (1)
where C m is the membrane capacitance ; V (x, y, t) is the transmembrane potential at location (x, y) at time t ; d is the diffusion parameter ; ∆ is the continuous Laplace operator ; i N a is the sodium current ; i other is the sum of other ionic currents and i stim is the stimulus current.

The dynamics of the currents in FHN model are characterized as Fast-Slow ones, where the sodium one is the fast part and is expressed as

i N a = k • g N a • V (V -V a )(V -V b ). (2)
Biologically, k intervenes in the control of the sodium current dynamics and kinetics of sodium channels. In this model, changing k results in modifying the dynamics of sodium current, then the value and the position of nullclines in the phase space. g N a is the maximum sodium conductance ; V a and V b are respectively the activation and inactivation potential ;

Substituting Eq.( 2) into Eq.( 1), we obtain :

C m ∂V ∂t = C m d∆V -k • g N a • V (V -V a )(V -V b ) -i other + i stim . ( 3
) Setting that v = V V b , α = V a V b
and converting the continuous Laplace operator to discrete Laplace operator with spatial step δh, yields :

∂v ∂τ = D∆v -k • g N a • v(v -α)(v -1) -I other + I stim , (4)
τ = g N a • V 2 b C m • t (5)
D = C m 1 g N a • V 2 b • δh 2 • d (6)
I stim = 1 g N a • V 3 b • i stim ( 7
)
where I stim , the normalized stimulus current, is defined by :

               I stim (x, y, τ ) = A(x, y)rect F θ τ * N n=0 δ(τ - n F ) A(x, y) =      A, if cell (x, y) is stimulated ; 0, otherwise. (8)
A(x, y) • rect( F θ τ ) is a rectangular pulse with amplitude A(x, y), frequency F , duty cycle θ, duration θ F , and centered in time τ = 0. The function

N n=0 δ(t - n F )
is a Dirac comb corresponding to a series of Dirac delta functions spaced at intervals of n F .

The current I other in FHN model is introduced by a recovery variable w which indicates the capacity of the medium to revert to its resting state after the propagation of impulsions :

∂w ∂τ = ε(v -γw) (9)
where ε and γ are dimensionless constant parameters.

The above equations allow us to obtain the realistic values from numerical ones and vice versa. Some typical values for the realistic parameters are used here : [START_REF] Scott | Neuroscience : a mathematical primer[END_REF][START_REF] Clayton | A guide to modelling cardiac electrical activity in anatomically detailed ventricles[END_REF]), d = 0.026 cm 2 s -1 . However, considering the difference of action potentials (AP) lengths 1 and the different properties of cell models, this conversion can only provide some qualitative features.

C m = 10 -6 F cm -2 , g N a = 23 mS V -2 cm -2 , V b = 122 mV ([
The method used to suppress SW is obtained by sending a train of monophasic stimuli from external electrodes distributed periodically in the system. The cardiac tissue is covered by an electrodes network. The size of each electrode is set to be the same as a cardiac cell's. The electrodes network has a square grid form but other patterns and sizes could be considered.

Optimization of parameters and Results

Simulations

Global stimulation

In the global stimulation case, the effects of the parameters of the stimulation signal on the suppression of SW are investigated. Under the hypothesis that there are SW evolving on the tissue, the whole tissue is stimulated. The results are presented by the relationship between the SR and I stim parameters (amplitude, frequency and duty cycle), as shown in Figs. 2, 3 and4.

There is a decreasing relationship between SR and each of these stimulation parameters. Indeed, the larger the applied energy is, the smaller the number of stimulated cells can be. However, there exists asymptotically a minimal stimulated cells number needed in the SW suppression process, whatever are the amplitude, frequency or duty cycle of the stimulus, which is in agreement with [START_REF] Allessie | Regional control of atrial fibrillation by rapid pacing in conscious dogs[END_REF]. The slope is very important for smaller amplitudes, lower frequencies and smaller duty cycles. This means that in these regions, small variations of amplitude or frequency could greatly change the result of suppression process of SW.

Furthermore, decreasing too much both parameters would increase SR which is contradictory to our main goal : stimulating less cells. As a consequence, a reasonable compromise should be taken. In following simulations, we take A = 83.53 µA cm -2 , F = 21.40 Hz, θ = 50% as our benchmark parameters.

Observing these results, it seems possible to suppress SW by choosing the appropriate parameters, but this is valuable only if the stimulation grid covers the whole tissue. This global stimulation is not possible in the clinical case, because it is difficult to stimulate the whole cardiac tissue. If this grid covers just a localized area of the tissue, SW can actually be suppressed in this area, but outside this area the stimulation provokes other waves (see wave fronts in Fig. 5c and5d, propagating with a growing rough square shape outwards the stimulated area). In certain conditions, for instance : with local inhomogeneity, presence of obstacles or collision of waves, the new generated wave could be transformed into spiral waves. This implies that the amplitude of I stim is above the threshold of new wave generation. Consequently this kind of stimulation cannot be considered as a real suppression of SW. As a result, it is necessary to determine this threshold as described in the following subsection.

(a) (b) (c) (d) (e) (f)

No Suppression of Spiral Waves

Suppression of

Wave generation threshold

Let us consider a tissue initially in a resting state and stimulated by an electrodes grid with fixed frequency and duty cycle values and let us investigate the influence of the variation of k on the generation of waves. In order to investigate the influence of k on wave generation, stimulations with a single electrode injected to the tissue initially in resting state have been performed. As shown in Fig. 7, the minimum amplitude of I stim leading to wave generation can reach A = 2.09 mA cm -2 when k = 0.75, which is 238 times larger than what is needed for k = 1 (A = 8.77 µA cm -2 ). So, in the next subsection, the effect of the modification of the parameter k combined to the SW suppression is studied.

No wave generation

Wave generation 0.75 0.8 0.85 0.9 0.95 1

Hybrid stimulation

From precedent results, we know that suppressing locally SW by electrical stimuli without generating new waves may be obtained by modifying the value of the parameter k. Biologically, it could correspond to use pharmacological agents in order to reduce the conductance of the sodium channels. This yields a hybrid method. First, our simulation results show that when k < 0.73, it is even not possible to excite the tissue. The cells have only sub-threshold oscillations even if the stimulation strength becomes substantial. Namely, considering the fact that FHN model is a fast-slow dynamical system, the velocity of traveling wave depends mainly on the fast dynamics i.e on the depolarization phase of the action potential. In this case, the front velocity value is mainly a function of the nonlinear threshold α of the sodium conductance and the diffusion parameter D [START_REF] Binczak | Reaction-diffusion network for geometric multiscale high speed image processing[END_REF][START_REF] Binczak | Experimental propagation failure in a nonlinear electrical lattice[END_REF] and depends also on the value of k. There exists a velocity threshold, below which the system is unstable and waves fail to propagate through the tissue [START_REF] Scott | Neuroscience : a mathematical primer[END_REF]. Since k controls the dynamics of sodium current, it is clear that a relatively small value of k will result in a too low excitability of stimulated tissue and the collapse of the waves. waves and the tissue is still excitable after the suppression of SW. This region is interesting in the optimization of the stimulation process.

Conclusion and Discussion

It is well known that atrial fibrillation (AF) begets AF [START_REF] Wijffels | Atrial fibrillation begets atrial fibrillation : A study in awake chronically instrumented goats[END_REF]. However, according to Allessie's concept, AF can also self terminate if the spiral waves (SW) number is not enough to self maintain this cardiac arrhythmia [START_REF] Allessie | Experimental evaluation of moe's multiple wavelet hypothesis of atrial fibrillation, Cardiac Electrophysiology and Arrhythmias[END_REF]. In fact, the ideal atrial stimulation/defibrillation method should be able to suppress enough SW in a dedicated region to stop fibrillating for the whole heart, but without generating new ones. In this intend, we studied here a simulated tissue of cardiac cells using FHN model. Compared to other some realistic ionic models, the FHN model is usually considered as a simple one and has its shortcomings.

However, we adopted it because it seems to be well adapted for modeling the main dynamics of ion channels we needed in our study.

In this numerical model study, we showed that it may be possible to suppress spiral waves (SW) inside a local area by using a hybrid strategy. This strategy combined the optimization of multifocal stimulations delivered via a grid of electrodes, and simulated sodium channel modification of the cardiac tissue. Although this strategy may only be a local and temporary solution (the confined area may be colonized by outside wandering SW), we believe that it could be of interest to restore the normal sinus rhythm under certain conditions considering for example Allessie's concept. In our study, we took into account physical parameters like the size and geometry of the electrode network as well as the number of the stimulating electrodes needed. We also identified some ideal pacing parameters in order to prevent from SW reinduction. However, in order to address this last issue, one of the options available is to modify the action potential dynamics (we choose here the sodium channel modification) in a simulated cardiac tissue.

In clinical practice, it is well known that antiarrhythmic drugs can optimize atrial defibrillation success rate. N a + channel blockers are commonly used to control the cardiac rhythm, and especially during AF [START_REF] Camm | Guidelines for the management of atrial fibrillation[END_REF][START_REF] Goldin | Resurgence of sodium channel research[END_REF][START_REF] Kneller | Mechanisms of atrial fibrillation termination by pure sodium channel blockade in an ionically-realistic mathematical model[END_REF][START_REF] Burashnikov | Atrium-selective sodium channel block as a strategy for suppression of atrial fibrillation : Differences in sodium channel inactivation between atria and ventricles and the role of ranolazine[END_REF]. Whereas, their use can be limited by potentially serious side effects like other antiarrhythmic drugs, they remain at the first line therapy in AF treatment [START_REF] Trial | preliminary report : effect of encainide and flecainide on mortality in a randomized trial of arrhythmia suppression after myocardial infarction[END_REF].

In summary, we showed that a hybrid strategy applied to the FitzHugh-Nagumo model helped to identify specific regions where it was possible to suppress SW without provoking new reentries. In order to confirm these encouraging results, it would be of interest to apply the same strategy to more realistic ionic models. The next step could also be to use in vitro culture of cardiomyocytes [START_REF] Jacquir | Cardiac arrhythmias induced by an electrical stimulation at a cellular level[END_REF][START_REF] Jacquir | In vitro simulation of spiral waves in cardiomyocyte networks using multielectrode array technology[END_REF][START_REF] Jacquir | Investigation of micro spiral waves at cellular level using a microelectrode array technology[END_REF] . The MEA (Microelectrode Array) technology could be a useful tool in testing the role of multifocal stimulation strategies on induced SW. As a

Fig. 1 with the parameter k = 1 .

 1 Fig. 1 with the parameter k = 1. Three categories of simulations have been carried out to optimize the suppression of SW : global stimulation, wave generation threshold by stimulation and hybrid stimulation.

Figure 1 :8Figure 2 :

 12 Figure 1: Suppression of SW at different times : (a) t = 0 ms ; (b) t = 26.29 ms ; (c) t = 181.11 ms ; (d) t = 219.09 ms ; (e) t = 265.82 ms ; (f) t = 391.43 ms. (Tissue size 200 × 200 cells grid, k = 1, A = 375.88 µA cm -2 , F = 21.40 Hz, θ = 50%, SR = 2.1%).

Figure 3 :8Figure 4 :

 34 Figure 3: SR vs. frequency of stimulus (A = 83.53 µA cm -2 , θ = 50%, k = 1).

Figure 5 :

 5 Figure 5: New waves generated outside the stimulated region in the SW suppression process. Tissue size : 350 × 350, stimulated region 150 × 150, spiral waves in the center 100 × 100. (a) t = 0 ms ; (b) t = 27.00 ms ; (c) t = 261.73 ms ; (d) t = 1123.80 ms. (k = 1, A = 83.53 µA cm -2 , F = 21.40 Hz, θ = 50%, SR = 1.96% in spiral waves region).

Fig. 6 10Figure 6 :

 66 Figure 6: Comparison of thresholds of wave generation for k = 1 and k = 0.75. (F = 21.40 Hz, θ = 50%).

3 k 2 )Figure 7 :

 327 Figure 7: Amplitude of stimulus vs. Parameter k characterizing the threshold of waves generation using one stimulating electrode. (F = 21.40 Hz, θ = 50%).

10 1 Figure 8 :Amplitude of Istim (µA cm - 2 )Figure 9 :

 1829 Figure 8: Comparison of the waves generation threshold and the SW suppression threshold with F = 21.40 Hz, θ = 50%. k = 1.

12Figure 10 :

 10 Figure 10: Comparison of the waves generation threshold and the SW suppression threshold with F = 21.40 Hz, θ = 50%. k = 0.75.

Sodium-based AP lasts for less than 1 ms, calcium-based AP lasts for 100 ms or longer and normal human atria1 AP lasts usually for 383 ± 103ms[START_REF] Dawodu | The shape of human atrial action potential accounts for different frequency-related changes in vitro[END_REF].

clinical perspective, we may infer that resetting one or more local regions of the left atrium via implanted array of electrodes connected to a dedicated pacing device may help to eradicate enough SW to get back to normal sinus rhythm.


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

Rayleigh limit and fragmentation of multiply charged Lennard-Jones clusters: Can charged clusters provide clues to investigate the stability of electrospray droplets?

Keywords: Monte Carlo, charged clusters, Rayleigh limit, fission, electrospray ionization

The fragmentation of multiply charged clusters composed of N ≤ 1000 Lennard-Jones particles augmented with electrostatic interactions is explored by classical Monte Carlo and molecular dynamics simulations with the stated goal of establishing possible analogies with electrospray droplets. Clusters with few charge carriers are shown to be only subject to particle ejection and their Rayleigh limit can be estimated by quantifying the loss of charged particles. On the contrary, uniformly charged clusters can both evaporate particles and undergo fission, making them better candidates to model electrospray droplets. Critical charges delimiting regions of instability of these clusters are defined from the calculation of lower order multipole moments and asymmetry parameters based on the knowledge of moments of inertia. The first discontinuity of quadrupole moments and asymmetry parameters is related to cluster elongation before twofold fission and the corresponding charge is deemed to be a good estimate of the Rayleigh limit. Octopole moments are negligible about this charge, their discontinuities arising at higher charges when threefold fissions are allowed. The size dependence of these critical charges is qualitatively predicted from Rayleigh's formula and the expression of surface energy advocated in liquid drop models. Deviations below 15% are commonly achieved when comparing Rayleigh limits extracted from experimental data with theoretical predictions based on Monte Carlo simulations or liquid drop models for a set of eleven atomic and molecular liquid clusters. Although manifold fission of uniformly charged clusters is unlikely close to the Rayleigh limit, successive asymmetric fissions are found to occur in conjunction with other fragmentation mechanisms, including the expansion of ring-shaped structures, at charges more than twice as large as the Rayleigh limit.

I. INTRODUCTION

The stability of multiply charged droplets dates back to the nineteenth century when Lord Rayleigh defined a critical charge Q R = π 8ǫ 0 γD 3 , the so-called Rayleigh limit, above which a charged, incompressible spherical droplet of diameter D and surface tension γ becomes unstable. 1 From then, this idea was extended to a myriad of species from atomic nuclei, 2 to metal clusters, 3 charged van der Waals clusters, 4 and electrospray droplets. 5 In fact, the stability of multiply charged van der Waals clusters is not generally defined by evaluating their Rayleigh limit but rather by determining a critical size n c (z) above which clusters with z charges start to stabilize and can therefore be detected by mass spectrometry measurements. For carbon dioxide cations, the critical sizes are n c (2) = 45, n c (3) = 109, and n c (4) = 216. 4 Dication clusters (CO 2 ) 2+ n (n > ∼ 44) do not feature any fission but have a significant propensity to evaporate one or two neutral monomers. On the contrary, both evaporation of neutral monomers and asymmetric fission are relevant to account for the formation of doubly charged fragments from (CO 2 ) 3+ n (109 < ∼ n < ∼ 123), the products of fission being a small (∼ 8% of the parent ion mass) singly charged cluster and a much larger (∼ 92% of the parent ion mass) doubly charged cluster. 4 This competition between evaporation and fission processes is not rare among van der Waals clusters and it has recently been invoked to interpret the mass spectra obtained upon Coulomb fission of (H 2 O) 2+ n and (NH 3 ) 2+ n . 6 The experimental critical sizes of nonpolar and heavy rare gas (Kr and Xe) clusters are also properly reproduced by theoretical calculations based on the liquid drop model (LDM) assuming a uniform distribution of the z charges. 4,7 The significant discrepancy obtained for the lighter neon clusters, namely an overestimation of n c (2) by a factor of at least 2.2, is still subject to debate since neither the inclusion of shell effects 8 nor that of vibrational delocalization 9 were able to significantly improve the agreement of the LDM to experiment. This result was all the more surprising that the prediction from the Rayleigh model, typically valid for small but macroscopic droplets, only underestimated the observed critical size of n c (2) = 284 by 10%. 10 The Rayleigh model is also of practical interest in proteomics 11 and atmospheric physics 12 to reckon the maximum charge an electrospray droplet can sustain, and two main mechanisms have been postulated to account for the distribution of fragments in electrospray ionization experiments coupled to mass spectrometric measurements (ESI-MS): the ion evapo-ration model (IEM) 13 and the charged residue model (CRM). 14 The first stage of both models is identical. The parent droplet, whose charge is about 70% of the Rayleigh limit, 15 shrinks by evaporating solvent molecules which mechanically increases the electrostatic repulsion between charge carriers and eventually leads to the asymmetric fission of the droplet into a cloud of "offspring" droplets that represent, as a whole, about 15% of the parent droplet charge and 2% of its mass. 16,17 The charge of the second-generation parent droplet (about 85% of the initial charge) then lies below the Rayleigh limit and additional solvent evaporations may induce new disintegration events, depending on the nature of the fragmentation mechanism. IEM and CRM differ in the way the ionized analytes are released into the gas phase after the production of the first offspring droplets. IEM assumes that an ion, with its solvation shell, can be desorbed from charged droplets with radii of about 10 nm, although the Rayleigh limit is not reached for such sizes, whereas CRM postulates that disintegration cycles should proceed until droplets can contain only one analyte. Despite debates, it is commonly accepted that IEM should hold for small ions and CRM for large systems such as globular proteins. 18 In the same way as electrospray droplets, multiply charged clusters can be prone to evaporation and fission events where the propensity for these decay processes is due to their activation barriers. A dimensionless fissility parameter X = E Coulomb /2E surf ace that relates the repulsive Coulomb energy to the cohesive surface energy is sometimes introduced to quantify the barrier height to fission of multiply charged clusters such as Morse 19 or metal clusters. 3,20 When X < 1 the shape of a uniformly charged droplet at the saddle point is prolate, Rayleigh's predictions should not be reliable, and fission can be observed if the barrier height remains smaller than the activation energy of other decay processes. At X = 1 the saddle point structure is spherical, the barrier to fission vanishes, and the Rayleigh model holds. When the Coulomb energy exceeds twice the surface energy (X > 1) homogeneously charged droplets as well as surface-charged metal clusters become unstable and undergo fission. In particular, molecular dynamics simulations performed on metal clusters (X > ∼ 1) with up to 923 atoms revealed a significant propensity for symmetric fission as the charge was uniformly distributed within the cluster, the fission process gaining in asymmetry as the cluster charge increased. 20 In the present work, we adopt a coarse-grained point of view where charged droplets are modeled as clusters composed of N Lennard-Jones (LJ) particles (N ≤ 1000) augmented with electrostatic interactions. Throughout the paper a particular emphasis will be devoted to the largest clusters (N = 1000) with uniform distribution of charges since these systems are subject to evaporations and fissions as it will be discussed in Section III A 1. This work first aims to evaluate the relevance of models based on multiply charged LJ clusters to get some insight into the stability of charged droplets, and especially electrospray droplets.

In particular, we were concerned with the devising of non-parameterized methods able to unambiguously characterize electrostatic instabilities responsible for fission events in large charged droplet-like systems. As detailed in Sections II A and III A 2, Monte Carlo (MC) simulations coupled to the calculation of multipole moments or asymmetry parameters are found to provide a simple way to identify some critical charges, chief among which is the Rayleigh limit of the charged clusters. A simple model based on Rayleigh's formula and the liquid drop model is then described in Section III B to provide a qualitative understanding of the size dependence of these critical charges. The suitability of this model together with the accuracy of MC calculations are evaluated on a set of eleven atomic and molecular clusters in the same section. These results are supplemented by molecular dynamics simulations performed at charges above the Rayleigh limit in Section III C to quantify the propensity for successive asymmetric fissions, fragmentation processes not detected in former studies dedicated to the fragmentation dynamics of 309-particle charged clusters. 21

II. METHOD

A. Monte Carlo simulations

The melting temperature of charged clusters, especially at the Rayleigh limit, is expected to be significantly lower than the melting temperature of the corresponding neutral clusters due to the electrostatic repulsion between charge carriers. The melting temperature of neutral clusters can thus serve as a reference temperature for Monte Carlo simulations on liquid charged clusters. This choice is motivated by the unambiguous definition of the melting temperature as maximum of the heat capacity curve for all the cluster sizes and by our will to minimize the influence of temperature on dissociation events while enabling the charged clusters, particularly at the Rayleigh limit, to be liquid. MC simulations are therefore divided into two stages: (i) parallel tempering (PT) simulations on neutral clusters to evaluate the position of their melting peak and (ii) parallel charging (PC) simulations 22 of charged clusters at the melting temperature of the neutral cluster to identify the critical charges delimiting regions of cluster instability.

Parallel tempering: melting temperatures of neutral clusters

Parallel tempering (PT) simulations 23,24 were performed in the NV T ensemble on a selection of neutral clusters with sizes ranging from N = 55 to N = 1000. The potentialenergy surface (PES) describing these clusters is a sum of LJ interactions

E neut = 4u N i<j σ r ij 12 - σ r ij 6 , (1)
where u and σ are the well depth and pair diameter of the LJ potential, and r ij is the distance between particles i and j. In typical LJ reduced units, namely E * neut = E neut /u and r * = r/σ, this PES writes

E * neut = 4 N i<j 1 r * ij 12 - 1 r * ij 6 . (2)
We used as input configurations the global minimum structure of LJ clusters reported in the Cambridge Energy Landscape Database. [25][26][27] Additional computational details on the method are reported in Section S1A of the Supplementary Material.

The melting temperatures of clusters composed of 55, 100, 147, 309, 561 and 1000 LJ particles are reported in Figure 1 as a function of N -1/3 . For N = 1000, the uncertainty on the melting peak position is probably larger than for smaller clusters, in part due to the reduced number of MC sweeps at the accumulation stage. However, the convergence of the melting curves of large neutral LJ clusters (N > ∼ 309) especially requires a significant number of MC sweeps (ie, more than 10 8 ) when the position and the magnitude of a premelting peak are sought, 28 which is not the object of the present work. Furthermore, the linearity of the curve in Figure 1 has already been pointed out by Pahl et al. 29 in their analysis of clusters binded by extended Lennard-Jones potentials. 30 This linear trend is expected to enable the extrapolation of the melting temperatures of large clusters up to the bulk, then circumventing the drawbacks of time-consuming MC simulations. However, the primary goal of the fitting procedure is here to assess the reliability of the full range of melting temperatures used as inputs in the parallel charging simulations of clusters with N ≥ 147 although we are aware that power laws in N -1/3 are generally more suitable for representing the properties of clusters further from magic numbers (N > ∼ 1000). As a complement, bulk melting points of rare gases extrapolated from our fitted curve are reported in Section S1B of the Supplementary Material and the good agreement to available experimental data and to the theoretical results from Pahl et al 29 is briefly discussed.

Parallel charging: Rayleigh limit of liquid charged clusters

Parallel charging (PC) 31,32 simulations performed to investigate the properties of multiply charged clusters are based on the same philosophy as PT simulations except that the replicas differ in the value of the charge carried by the charged particles rather than the cluster temperature, and configurations of the neutral cluster at its melting point are used as input configurations. This subsequently modifies the acceptance probability for swapping two configurations A and B according to

acc(AB → BA) = min {1, exp ([E * A (Q * A ) + E * B (Q * B ) -E * A (Q * B ) -E * B (Q * A )]/T * )} (3)
where E * i (Q * j ) (i, j ∈ {A, B}) is the reduced potential energy of configuration i with total charge Q * j = Q j / √ 4πǫ 0 σu, T * = k B T /u is the reduced temperature with Boltzmann's contant k B , and ǫ 0 is the vacuum permittivity. The potential energy of multiply charged clusters is then obtained by adding electrostatic charge-charge interactions to Eq. ( 2),

E * = 4 N i<j 1 r * ij 12 - 1 r * ij 6 + N i<j q * i q * j r * ij (4)
where q * i is the reduced charge on particle i. For the sake of simplicity, the charged particles will always be assigned the same charge, q * i = q * for all i. PC simulations have been performed on clusters composed of N = 55, 100, 147, 309 and 1000 particles. However, our discussion will mainly focus on LJ clusters with N ≥ 147 that will be shown to represent good candidates for discriminating between stability patterns of small and large clusters.

We will investigate the influence of low (∼ 10%) and high (100%) proportions of charge carriers for a broad range of total cluster charges extending from zero to charges above the Rayleigh limit. Additional information on the PC simulations is provided in Section S1C of the Supplementary Material.

Several approaches afford to provide some estimate of the Rayleigh limit of multiply charged clusters from PC simulations without requiring the knowledge of the cluster surface tension. The Rayleigh limit can be defined as the critical charge that entails partial evaporation in more than 50% of equilibrium structures, 33,34 where a particle can be assumed to evaporate as soon as it reaches a spherical layer of 0.5σ from the container edge. This criterion is probably suitable for small clusters where thermal evaporation of neutral particles is negligible but it is not reliable for large clusters due to their much higher melting temperatures and their propensity for fissions. In such a case, an analysis of fragment distribution can supply additional information on charge-driven instabilities as it will be discussed in Section III A 1. However, a more astute way to evaluate the Rayleigh limit is probably to calculate the multipole moments of the charged clusters, these quantities depending only on the cluster positions and charges. Unlike the aforementioned methods, no parameterization is needed and multipole moments intrinsically contain some information on cluster deformations. As an example, the elongation of a sphere along one direction is represented by an increase of its quadrupole moment. This concept is exploited in nuclear physics to investigate the deformations of spheroidal atomic nuclei that arise from the decomposition of the length of the radius vector pointing from the nucleus origin to its surface in terms of Y lm (θ, φ) spherical harmonics. When nuclei are axially symmetric φ can be discarded from variables since m = 0 and spherical harmonics become real functions of cos θ. The length of the radius vector then writes 2,35

R(θ) = c(α)R 0 1 + lmax l=2 β l 2l + 1 4π P 0 l (cos θ) (5)
where c(α) is determined to ensure volume conservation condition, R 0 = r 0 A 1/3 with r 0 the estimated radius of nucleons excluded spherical volume (eg, r 0 = 1.4 fm 35 or r 0 = 1.2 fm 36 ) and A the atomic mass number, β l are the deformation parameters, and P 0 l (cos θ) are Legendre polynomials of order l.

The surface particles, and therefore the length of the radius vector, of charged clusters could be easily determined, 34 and it would be possible to calculate deformation parameters from least-square fitting of Eq.( 5). However, their values would be meaningless for dissociated equilibrium structures due to the decomposition in spherical harmonics. In the rest of this paper we therefore limited our inspection of the properties of charged clusters to the relationship between the Rayleigh charge and the values of the lower order multipole moments, especially quadrupole and octopole moments, without attempting to derive any deformation parameter from these multipole moments although the concepts behind deformation parameters will be used to rationalize the early stages of the fragmentation processes. The components of the multipole moment of order l are derived from the general expression 37

ξ(l) αβ•••ν = (-1) l l! n i=1 q i r 2l+1 i ∂ ∂r iν • • • ∂ ∂r iβ ∂ ∂r iα 1 r i (6)
where α, β, • • • , ν ∈ {x, y, z}, q i is the charge on particle i, r i the norm of the position vector of particle i, and n the number of charge carriers. The quadrupole moments (l = 2) and octopole moments (l = 3) then write ξ( 2)

αβ = Θ αβ = 1 2 n i=1 q i (3r iα r iβ -r 2 i δ αβ ), (7a)
ξ(3) αβγ = Ω αβγ = 1 2 n i=1 q i [5r iα r iβ r iγ -r 2 i (r iγ δ αβ + r iβ δ αγ + r iα δ βγ )], (7b)
where their reduced counterparts are easily deduced from ξ(l) *

αβ•••ν = ξ(l) αβ•••ν / σ l √ 4πǫ 0 σu .
Dipole moments are not considered here since they are null whatever the cluster size and charge, in close analogy with multipole deformations of atomic nuclei as clear from the absence of l = 1 in Eq.( 5). In order not to overload the figures reported in Section III, we will mainly focus our discussion on multipole moments ξ(l) αα•••α (l = 2 -3, α ∈ {x, y, z}) but the conclusions we will draw would remain valid for other components of the multipole moments.

Moments of inertia can also be considered as valuable quantities to evaluate cluster deformations. In particular, overall parameters such as the asymmetry parameter commonly used to investigate the shape of asymmetric top molecules in rotational spectroscopy 38 may be practical alternatives to multipole moments. This parameter is defined as

κ = 2B -A -C A -C (8)
where

A = I -1 a , B = I -1 b , and C = I -1 c , I x (x ∈ {a, b, c}) being the three moments of inertia. It evolves from -1 (prolate symmetric top, B = C) to +1 (oblate symmetric top, A = B),
and it is zero for asymmetric clusters (eg, A = 3C and B = 2C). However, it is not defined for spherical or icosahedral systems (A = B = C) which may harm data interpretation when such highly symmetric structures are explored by MC simulations. For this reason we will preferentially discuss results related to multipole moments, then referring occasionally to moments of inertia and asymmetry parameters to complete the analysis of our results.

B. Dynamics of uniformly charged clusters

The Monte Carlo simulations of uniformly charged 1000-particle clusters are supplemented by molecular dynamics simulations at constant energy to describe more precisely their relevant fragmentation mechanisms. Six independent batches of 160 trajectories are run at a total charge of Q * = 180 to evaluate the propensity for successive asymmetric fissions. The input configurations come from preliminary PT simulations performed in relatively small containers (R c = 8σ) at T * = 0.53 to prevent any fission of the clusters before running the dynamics. Newton's equations of motion are integrated over the reduced time

t * = t u mσ 2
, where m is the particle mass, via a standard velocity Verlet algorithm by making use of the potential-energy surface presented in Eq.( 4). [START_REF] Bonhommeau | [END_REF] Each trajectory lasts 10 6 reduced time steps of value 10 -4 , and information on the trajectories is collected every 10 3 steps to enable subsequent analyses of fragment distributions. The parameterized procedure employed to follow the production of fragments along the dynamics and to characterize the Monte Carlo equilibrium structures is detailed in Section S2 of the Supplementary Material.

III. RESULTS

A. Rayleigh limit of charged clusters

Estimation of the Rayleigh limit from the analysis of fragment sizes

Two mechanisms of decay have been identified on 309-particle LJ clusters at charges close to the Rayleigh limit, 21 namely particle ejection (sometimes called evaporation here) and twofold fission. The latter process only occurs when the number of charge carriers is large (n ≈ N) whereas the former persists at all n. The same behavior is encountered on smaller (N = 147) and larger (N = 1000) clusters as evidenced by Figure 2 that reports the average number of monatomic fragments, the average fragment size when excluding monomers and the largest fragment, and the average size of the largest fragment relative to the parent cluster size, as a function of the charge ratio Q * max (1000, 100) = 90.0, and Q * max (1000, 1000) = 170. The size of clusters poor in charge carriers (n = N/10) hardly decreases by a few percents due to an increasing loss of monomers when the cluster charge increases, as illustrated in Figures 2a and2c. The average size of other fragments never exceeds 2 (see Figure 2b) which indicates that equilibrium structures contain one dimer at most. In fact, trimers may also be formed with a negligible abundance that roughly reaches 0.2 trimers per equilibrium structure when focusing on the dissociation of the largest parent cluster in the highest charge state (N = 1000 , Q * = 90). Progressive ejection of monomers is also the main source of particle loss for uniformly charged clusters at 2d and2e. At higher charges, the dissociation process yields a decrease of the production of smaller fragments, such as monomers, parallel to the formation of larger fragments as evidenced by the sharp increase of the average fragment size in Figure 2e. Two large fragments are generally produced (see Figure S5a given in the Supplementary Material), the largest one being roughly half as large as the parent cluster (see Figure 2f). Above this threshold for twofold fissions, fragmentations into three or four fragments can also occur, the largest fragment then containing about one third or one fourth of the parent cluster particles. The Rayleigh limit of uniformly charged clusters can be presumably identified as the charge above which fission proceeds (ie, Q * ≈ (0.7-0.8) Q * max in our simulations) but such a criterion cannot apply to clusters poor in charge carriers due to their steady loss of small fragments by evaporation. For such clusters, the Rayleigh limit can be estimated by evaluating the charge Q * evap for which 50% of the parent clusters have lost one charge carrier, this loss being measured by the propensity the particle has to reach a layer of width 0.5σ from the container edge, or by estimating the charge that maximizes the fraction of charged particles at surface. 33,34 As an example, the percentage of fragmented parent clusters, presented in Figure 3, can help with the determination of Q * evap . When n = N/10, the fragmentation dynamics of small charged clusters (N < ∼ 309) is mainly governed by the evaporation of monatomic fragments, and especially charge carriers. 21 Although the evaporation of neutral particles is not negligible for the larger clusters, its amount remains roughly constant as the total cluster charge increases. The ejection of neutral particles is then probably more influenced by the value of temperature, that is a fixed parameter in PC simulations, than possible cluster rearrangements stimulated by the ejection of the charged particles, whose rate differs between replicas. Therefore, the loss of charge carriers can be used as a signature of electrostatic instability to define the Rayleigh limit of these clusters. The charge that leads to 50% of equilibrium structures subject to evaporation is

Q * /Q * max (N, n). Q * max (N, n), referred as Q * max in
Q * < ∼ (0.7 -0.8) Q * max as presented in Figures
Q * evap (147) ≈ 24, Q * evap (309) ≈ 34
, and Q * evap (1000) ≈ 56, in good agreement with recent Rayleigh limit estimates performed on 309-particle clusters at three different pressures that lead to Q * evap (309) = 37 ± 6. 21 However, we must keep in mind that evaporation is favored in larger containers where fragments have lower electrostatic interactions and Q * evap could be an upper bound of the Rayleigh limit. This comment should also hold true for uniformly charged clusters.

Alike results could have been obtained by plotting the charge density of these clusters as a function of their charge as depicted in Figure S1 of the Supplementary Material. The charge density increases linearly while the cluster charge remains below the Rayleigh limit due to the infrequent particle ejections that cannot counterbalance the gain in charge brought by the charge carriers. At the Rayleigh limit, evaporation accelerates significantly and the charge density starts to level off. This change of regime of the charge density could be exploited to locate the Rayleigh limit but the percentage of fragmented parent clusters seems easier to interpret. It is worth noting that the stability of small uniformly charged clusters could also be defined in terms of percentages of fragmented parent clusters, as shown in Figure 3b.

However, the significant loss of particles at low charge combined with the mixing of thermal evaporation and fission makes this criterion ineffective for large clusters.

Non-parameterized method to estimate the Rayleigh limit

The aforementioned criteria are practical methods to estimate the Rayleigh limit of multiply charged clusters although they require the setting of some parameters: the distance above which two particles are considered to belong to different fragments for the method based on fragment detection, the spherical layer width and the percentage of trajectories above which the limit of instability is assumed to be reached for clusters only prone to evaporation. In contrast, multipole moments are a non-parameterized alternative to track the charge-driven instability of charged clusters. Indeed, these quantities contain the information on cluster deformations related to changes in charge distribution, and discontinuities specific to fission events should arise at the vicinity of the Rayleigh limit. Furthermore, evaporation is isotropic on average whereas fission occurs along preferential directions, and multipole moments should therefore enable a discrimination between these two decay pro-cesses. Quadrupole and octopole moments of clusters poor in charge carriers oscillate about zero whatever the cluster size, as illustrated in Figure S2 of the Supplementary Material.

The magnitude of oscillations can be reduced by simply increasing the number of configurations included in the averages of multipole moments (see Figure S3 of the Supplementary Material). Weak multipole moments are compatible with a spherical distribution of charges, expected for isotropic evaporations. No peculiar charge-driven marker of instability is noticeable on the plots and Q * evap cannot be deduced from the fluctuations of multipole moments. The same trend can be observed when plotting the asymmetry parameter κ as a function of the cluster charge, κ remaining almost constant at a small value close to that obtained for the neutral cluster (see Figure S4 of the Supplementary Material). Moreover, the slope of moments of inertia slightly changes at Q * evap but the location of the inflection point is not as precise as in Figure 3. On the contrary, clear discontinuities are visible on quadrupole and octopole moments of clusters with uniform charge distribution, as depicted in Figure 4 for N = 147, 309, and 1000. Three regions emerge from the analysis of these plots that can be easily related to the analysis of fragments discussed in the previous section and to Figures S5b-d supplied in the Supplementary Material that report the average number of fragments containing more than 5 and more than N/10 particles with respect to the parent cluster charge. In region I evaporation of charged particles is the unique fragmentation process whereas twofold and manifold fissions (more than two fragments formed by fission) characterize regions II and III, respectively. Fissions in regions II and III are also accompanied by evaporation, as revealed in Figures 2e and2f by the decrease of the fragment sizes as a function of the cluster size in these regions, but evaporation can be disregarded to characterize the charge-driven instability of these clusters.

As the charge increases, enhanced Coulomb repulsions tend to destabilize the cluster by altering the balance between cohesive LJ and repulsive forces. Uniformly charged clusters can respond to this instability by elongating to maximize the number of particles at surface which decreases the value of the fissility parameter compared to the spherical shape. Fission occurs when the cluster cannot accommodate the charge any longer. Therefore, the charge separating regions I and II, denoted Q * I , can be identified as the Rayleigh limit of the cluster, or at least an upper bound of this critical charge due to the influence of the fixed-volume container. At Q * I , all the clusters present discontinuities of their quadrupole moments, and fission into two fragments occurs. Octopole moments are negligible in region II but they significantly increase after passing through Q * II , the limit charge between regions II and III. This difference between quadrupole and octopole moments can be simply comprehended by reminding that spheroidal deformations can be evaluated by decomposing the length of the radius vector of any compact cluster in spherical harmonics. Axially symmetric multipole deformations (l = 1 to 4) are plotted in Figure S6 of the Supplementary Material for β l = 1 as examples of typical spheroidal deformations. Quadrupole deformations (l = 2) correspond to elongations, prolate deformations along the axis of minimum moment of inertia, prior to twofold fissions whereas octopole deformations (l = 3) yield a three-lobe shape that heralds future threefold fissions. However, we must not conclude that higher order multipole moments would provide accurate estimates of critical charges specific to manifold fissions implying more than three fragments. In such a case, we would expect hexadecapole moments to exhibit a discontinuity at Q * > Q * II whereas their magnitude is non-zero at any charge above Q * I according to Figure S7 given in the Supplementary Material. As a consequence, Q * I and Q * II might be the only critical charges whose values can be predicted from the analysis of multipole moments. This is not a significant restriction of the method to attempt comparisons with electrospray droplets since the dynamics of these droplets is governed by the deviation between their charge and the Rayleigh limit, whose estimate is supplied by

Q * I .
It is also worth noting that discontinuities are also visible on moments of inertia and asymmetry parameters presented in Figure S8 of the Supplementary Material. At all sizes, the first discontinuity can be easily interpreted as an estimate of the Rayleigh limit since the asymmetry parameter reaches a value of -1 specific to prolate shapes. However, the origin of the second discontinuity can only be ascertained by plotting the octopole moments presented in Figure 4.

B. Rayleigh limit and liquid drop model

As pointed out in the introduction, the Rayleigh limit depends on surface tension γ and droplet diameter D according to Q R = π 8ǫ 0 γD 3 . When posing the reduced surface tension γ * = γσ 2 /u and the reduced diameter D * = D/σ the Rayleigh limit then writes

Q * R = 2πγ * D * 3 . (9)
A precise evaluation of the Rayleigh limit from this formula would require the knowledge of surface tension, that is the surface Gibbs energy per unit area. In order to rationalize the trend observed on Rayleigh charges reported in Figure 5a without entering the complexity of an accurate evaluation of the surface Gibbs energy, a crude estimate of surface tension is obtained by dividing the surface energy of the droplet by its surface area. The reduced surface energy is derived from a liquid drop model introduced by Kreisle and coworkers to investigate the dissociation channels of multiply charged CO 2 clusters 40

E * s = AN 2/3 + BN 1/3 ( 10
)
where the dimensionless parameters A = 10.3 and B = -2.7 are extracted from a work by Briant and Burton devoted to the properties of argon microclusters. 41 Including this expression together with the cluster radius R * = (3N/4π) 1/3 yields γ * LDM and Q * LDM , our estimates of γ * and Q * R based on the liquid drop model, according to

γ * LDM = 1 (36π) 1/3 A + BN -1/3 (11a) Q * LDM = 144 3π 1/6 AN + BN 2/3 1/2 . (11b)
At large N values, when the cluster approaches bulk properties, γ * LDM = A/(36π) 1/3 ≈ 2.13 and the leading power of Q * LDM is N 1/2 as expected from Eq.( 9) at constant surface tension. Moreover, γ * LDM converges to the bulk value at N ≈ 500, as illustrated in Figure 5b, in good agreement with results obtained on water droplets where surface tension is expected to be roughly constant for droplet radii exceeding 1-2 nm. 42 The overall trend of Q * LDM is in agreement with Q * I calculated from Monte Carlo simulations despite an overestimation by 75% for the smallest cluster (N = 55) and 35% for the largest one (N = 1000). This deviation is attributed to the crude description of surface tension considered in Eq.(11a) that also accounts for the better agreement of Q * LDM with Q * II , a critical charge above the Rayleigh limit.

The overestimation of surface tension by γ *

LDM is evaluated on a set of 11 multiply charged clusters composed of rare gases, homonuclear dimers and polyatomic molecules whose properties are reported in Table I. [START_REF] Lide | CRC Handbook of Chemistry and Physics, Internet Version[END_REF][START_REF] Jacobsen | Thermodynamic Properties of Cryogenic Fluids[END_REF][START_REF] Won | [END_REF] Reference surface tensions are mainly derived from Ferguson's equation 46

γ = γ 0 (1 -T /T c ) p (12)
where T c is the critical temperature, γ 0 is a parameter characteristic of the liquid and p is deemed to be a universal parameter. 47,48 Such a parameterization is not available for water and ethanol whose surface tensions have been taken from the literature, [START_REF] Lide | CRC Handbook of Chemistry and Physics, Internet Version[END_REF] and for krypton whose surface tension has been estimated from the relationship introduced by Fuks and Bellemans 49

γ Kr = 16.344 -0.2089(T -116) (13) expressed in dyne/cm. It is also worth noting that properties of ethanol at 25 • C have been considered in place of data at temperatures close to its melting point due to the absence of surface tension values at low temperatures for this molecule. γ LDM are deduced by converting γ * LDM (∞) = A/(36π) 1/3 to atomic units from the knowledge of LJ pair diameters and well depths available for the 11 species. 4,[50][51][52] The analysis of the ratios γ LDM /γ reveals that surface tension is typically overestimated by an average factor of 2.04. Thus, a simple way to improve the quality of the correspondence between Q * LDM and Q * I is to divide γ * LDM by 2.04. The scaled surface tension yields much better agreement between Q * LDM and Q * I , the deviation dropping to 3 -5% for clusters with sizes between N = 100 and N = 1000, as illustrated by the blue curve plotted in Figure 5a.

These two theoretical estimations of the Rayleigh limit can be compared to the critical charge Q R that would arise from Rayleigh's formula on the basis of experimental densities and surface tensions of realistic droplets. Table I Q LDM (dashed curve) and Q R (solid curves) depicted in Figure 6 is always very similar. The agreement is especially remarkable for rare gases, water and ethanol, then demonstrating that a suitable approximation of Rayleigh limits of homogeneous clusters and droplets can be obtained by performing MC calculations, and in particular PC simulations, on uniformly charged LJ clusters. The prediction of Rayleigh limits associated to more realistic multicomponent electrospray droplets is however not possible with single LJ interactions. Binary or even ternary LJ interactions are required to evaluate the suitability of such MC approaches for determining the Rayleigh limits of mixtures.

C. Fragmentation of uniformly charged clusters

Successive asymmetric fissions are commonly invoked to account for the fragmentation of electrospray droplets into streams of tiny daughter droplets. Although uniformly charged clusters can undergo fissions, series of fissions are unlikely close to the Rayleigh limit. For instance, a twofold symmetric fission of 1000-particle clusters at Q * I (1000) ≈ 117 should roughly lead to two 500-particle fragments with a total charge of 58.5 well below their Rayleigh limit of Q * I (500) ≈ 78 derived from Figure 5a. The fragments are therefore stable with respect to fission and particle ejection should be the only fragmentation process in play at later times. In contrast, the fragments may remain unstable after the first fission if the parent cluster charge exceeds 156, a charge that enables threefold fissions (region III in Figure 4). Constant energy MD simulations have therefore been performed on 1000particle uniformly charged clusters at Q * = 180 > 2 Q * I (500) to evaluate the propensity for successive fissions and their characteristic time scales. Meaningful mass spectra obtained along the dissociation dynamics of these clusters are depicted in Figure 7. Evaporation is the main fragmentation process in play at the earlier times of the dynamics and most of the parent clusters have lost a few tens of particles at t * < ∼ 30. At t * = 40, a large amount of parent clusters have already undergone symmetric or asymmetric fission which produces a broad distribution of fragments extending up to M * f rag ≈ 750. Interestingly, the mass spectra are steadily moved to lower masses as time proceeds until yielding a distribution centered at M * f rag ≈ 250, one fourth of the parent cluster mass, that should be the result of additional fissions. This indicates that the fragmentation dynamics of highly charged 1000-particle clusters can be sequential. These behaviors observed at four selected times of the dynamics are confirmed by the analysis of the average production of fragments presented in Figure 8a.

The disappearance of the parent clusters is complete at about t * = 10 due to a progressive ejection of charge carriers that results in the exponential decay of the number of parent clusters. Significant losses of particles by symmetric or asymmetric fissions only start to occur at t * ≈ 30. These primary fission events can be followed by additional fissions when the fragments are composed by more than 400 charge carriers (Q * > 72) since their charge then exceeds their Rayleigh limit, Q * I ≈ 70. However, we must emphasize that a few fragments with a bit less than 300 particles (Q * = 54 and Q * I ≈ 60) have also been found to dissociate. This deviation to expectations could be related to the high-temperature structure of these fragments since we must remind that the temperature of MD simulations is initially fixed to the melting temperature of the neutral parent cluster (T * = 0.53) and the temperature of the system can still increase during the dissociation process. Therefore, the values of Q * I are unambiguously upper bounds of the Rayleigh limits specific to charged fragments. After a simulation time of t * = 100 the production of fragments is much slower but it is not completely ended and mass spectra could still shift to lower masses (M * f rag < ∼ 250). Note that the significant abundance of small fragments (less than 50 particles in Figure 8a) is especially governed by the ejection of monatomic ions and it is not a signature of manifold fissions. Evaporation also accounts for the existence of about 12.4% of partially evaporated parent clusters at the beginning of the dynamics. This initial source of particle loss could have been avoided by reducing the radius of the spherical containers used in the preliminary PT simulations but it is not expected to alter the conclusions of our simulations.

A close inspection of the 960 trajectories proves that uniformly charged clusters dissociate according to two types of fragmentation patterns where asymmetric fission prevails: (i) an elongation leading to primary fissions that can be followed by secondary, and even tertiary, fissions in 58% of the trajectories and (ii) the expansion of a ring-shaped structure leading to fragments moving apart in opposite directions, the number of fragments produced depending on the ring geometry, in 36% of the trajectories. 6% of the trajectories correspond to nascent rings that rapidly break up because of cluster elongation. These trajectories have not been included in one of the two families aforementioned. Percentages of trajectories provided in the rest of the text are always calculated with respect to the total number of trajectories to avoid any ambiguity. Four groups of fragmentation patterns have been identified to characterize the 58% of manifold fissions initiated by cluster elongation. Two thirds (∼ 43%) correspond to series of asymmetric twofold, and to a lesser extent threefold, fissions. In ∼ 4% of the trajectories, successive asymmetric fissions may remind the jet of tiny droplets ejected by electrospray droplets when their charge reaches the Rayleigh limit, as illustrated in Figure 9a, and the proximity in time of several asymmetric fission events can sometimes be interpreted as fourfold or fivefold fissions (< 3%). The second group (∼ 12%) is very similar to the first one except that symmetric fissions, defined in the present work as fissions where the fragment sizes differ from each other by less than 10%, can occur in conjunction with asymmetric fissions (see Figure 9b). Trajectories only subject to symmetric fissions are extremely rare (< 1%). A more exotic mechanism is the trigonal elongation of the cluster (∼ 3%) that leads to threefold fissions when the three lobes are not deformed (see Figure 9c). In the contrary case, four fragments can be produced.

Beside fragmentations related to cluster elongations, manifold fissions of multiply charged clusters at Q * = 180 are also found to dissociate subsequent to the expansion of three kinds of ring-shaped structures. The two most abundant are quadrilateral (∼ 15%) and triangular (∼ 14%) shapes that are depicted in Figures 9e and9f, respectively. The expansion of circular structures (see Figure 9d) is found half as frequent as the previous mechanisms (∼ 7%). The pictures in Figure 9 represent ideal cases and it is not rare to identify a nearly triangular ring that finally deforms into a more quadrilateral or circular shape then producing more than three fragments, all the more since large fragments can undergo additional fissions if their charge is above their Rayleigh limit. However, fissions coming from the expansion of triangular rings will yield less fragments than those originating from the breaking up of quadrilateral or circular rings, on average. Fragments produced by the fission of triangular rings are also less symmetric than those coming from trigonal elongations. Finally, it is worth noting that no star-shaped clusters, theoretically predicted by force field MD simulations on solvated macroions, 53 have been detected in our simulations.

On average, typical trajectories form 3 to 5 fragments composed of more than 100 particles, as illustrated in Figure 8b, which is in agreement with PC simulations that predict the formation of at least three fragments at Q * = 180 > Q * II ≈ 154. The large fragments can also be surrounded by some smaller fragments (less than 100 particles) and a cloud of about 60 monatomic ions, as confirmed by Figure S9 of the Supplementary Material.

IV. CONCLUSION

The search for similarities between charged clusters and realistic systems is a long-term endeavour. Last et al. 19,54 Electrospray droplets are commonly assumed to fragment by asymmetric fission cycles releasing streams of tiny "offspring" droplets carrying 1-2.3% of the mass and 10-18% of the charge of the parent droplet. 17 Unfortunately, no multiply charged LJ cluster can dissociate in this way. Clusters poor in charge carriers are only prone to evaporations and uniformly charged clusters, though both subject to evaporations and fissions, cannot increase their charge density by losing particles after the first fission event, which makes successive asymmetric fissions impossible at charges close to the Rayleigh limit. As demonstrated on 1000-particle clusters, series of asymmetric fissions can only occur at charges well beyond the Rayleigh limit where twofold fissions have a significant probability to produce unstable fragments, namely fragments whose charges exceed their Rayleigh limit. However, at these high charges other fragmentation channels are also open, like the formation of ring-shaped structures that fragment by expansion of the ring. Despite the limited analogies between electrospray droplets and charged LJ clusters, these systems can be used to devise approaches aiming to investigate the stability of charged droplets. In particular, the Rayleigh limit of clusters poor in charge carriers can be estimated by quantifying the ejection rate of charged particles, by reckoning the charge density of the cluster or by identifying the charge that maximizes the fraction of charged particles at surface. 34 The stability of systems prone to fissions, such as uniformly charged LJ clusters, is better characterized by the calculation of multipole moments or asymmetry parameters that enable a clear discrimination between evaporation, isotropic on average, and fission processes. Regions of instability can be defined on the basis of multipole moments discontinuities that arise at some critical charges.

When identifying the lowest of these charges, called Q * I , as the Rayleigh limit of the cluster, good agreement is obtained with previous MC simulations in the NP T ensemble. 21 The size dependence of Q * I is also properly reproduced by Rayleigh's formula when expressing surface tension as the LDM surface energy by unit area. Comparisons with Rayleigh limits predicted for charged droplets composed of rare gas atoms, homonuclear dimers and polyatomic molecules yield deviations that rarely exceed 15%. Although any extensive comparison between charged clusters and multicomponent electrospray droplets would be risky, the stability of charged droplets described by all-atom force fields could presumably take advantage of parallel charging (PC) simulations coupled to the calculation of multipole mo-ments. In our PC simulations replicas differ in their reduced charge q * = q/ √ 4πǫ 0 σu and the force field-based analogs could consider replicas that differ in their charge or LJ parameters. This approach might provide an effective way to get some insight into the stability of charged droplets to guide MD simulations at lesser computational cost. For instance, MD simulations performed in the present work on 1000-particle clusters at Q * = 180 were 2.2 times as long as PC simulations performed at 24 different charges on clusters of same size.

Practically speaking, charge-driven instabilities of electrospray droplets are however mainly tackled by molecular dynamics simulations making use of classical force fields to model solvated small ions, 55 the formation of salt clusters 56 or the release of proteins into the gas phase. 57 Beside the questions related to the diversity of analytes produced by ESI, Consta et al. 58 introduced a transfer reaction coordinate, based on the evaluation of the probability a fictitious particle located on an ion has to reach another ion after jumping on intermediate solvent molecules, to track the fragmentation of charged droplets. This kind of reaction coordinate is not applicable to coarse-grained uniformly charged clusters composed of identical particles. However, it might help to follow more precisely the dynamics of heterogeneous clusters coming from the replacement of simple LJ interactions by binary or ternary LJ potentials, then enabling the modelization of the properties of droplet-like systems without entering the complexity of all-atom simulations.

V. SUPPLEMENTARY MATERIAL

Complements on Monte Carlo simulations are collected in Section S1. Computational details on parallel tempering and parallel charging simulations are given in Sections S1A and S1C. The extrapolation of bulk melting points of rare gases from the melting temperatures of LJ clusters are reported in Table S1 and discussed in Section S1B. Section S2 describes the algorithm used in MC and MD simulations for the detection of fragments.

Section S3 provides supplementary data on cluster deformations for clusters poor in charge carriers (Section S3A) and uniformly charged clusters (Section S3B). Section S3a collects figures of charge densities (Figure S1), multipole moments (Figure S2), test of convergence of quadrupole moments (Figure S3), moments of inertia and asymmetry parameters ( Species General properties [START_REF] Lide | CRC Handbook of Chemistry and Physics, Internet Version[END_REF][START_REF] Jacobsen | Thermodynamic Properties of Cryogenic Fluids[END_REF][START_REF] Won | [END_REF] FE parameters [46][47][48] Surface tensions [START_REF] Lide | CRC Handbook of Chemistry and Physics, Internet Version[END_REF]49 LJ parameters Symbols and curves with same color correspond to the same species.

  justified their work on the dynamics of fission and Coulomb explosion of multiply charged Morse clusters by analogies with Coulomb instability in atomic nuclei or charged droplets, and pseudo-Coulomb instability in optical molasses. Calvo 20 compared the effect of uniform and surface distribution of charges on the fragmentation dynamics of sodium clusters as first candidates to better comprehend the fragmentation patterns of metal clusters. In the present work, we were concerned by the stability of multiply charged LJ clusters (N ≤ 1000) and possible analogies with the fragmentation of electrospray droplets.

TABLE 1 :FIG. 1 :FIG. 2 :FIG. 4 :FIG. 6 :FIG. 7 :

 112467 FIG. 1: Melting temperatures T *m of neutral LJ clusters (circles) as a function of N -1/3 . The linear fit of the data ranging from N = 100 to N = 1000 is plotted in red.

  TABLE I. Parameters of LJ potentials (pair diameter σ and pair well depth u), parameters of Ferguson's equation (FE) and correspondingsurface tensions γ, LDM bulk surface tension γ LDM from Eq.(11a), LDM Rayleigh limit Q LDM from Eq.(11b), Rayleigh limit Q I extracted from PC simulations, and Rayleigh limit Q R predicted by Rayleigh's formula, for rare gas clusters and a selection of clusters composed of dimeric and polyatomic molecules at temperatures close to their melting point. The uncertainty on Q I comes from the uncertainty on Q * I due to the charge difference between replicas surrounding the discontinuity of quadrupole moments. δQ I (resp. δQ LDM ) is the relative deviation between Q I (resp. Q LDM ) and Q R .

FIG. 1 .FIG. 3 .

 13 FIG. 1. Melting temperatures T * m of neutral LJ clusters (circles) as a function of N -1/3 . The linear fit of the data ranging from N = 100 to N = 1000 is plotted in red.

  the rest of the text, is the charge of the highly charged replica used in PC simulations of LJ clusters composed of N particles and n charge carriers, namely Q

* max (147, 15) = 37.5, Q * max (147, 147) ≈ 52.9, Q * max (309, 31) = 62.0, Q * max (309, 309) ≈ 81.6,

  collects the Rayleigh limits obtained for charged clusters composed of 1000 particles. The deviation between Q LDM or Q I and Q R lies below ∼ 5% for rare gases, remains below 15% for N 2 , CO, CH 4 , H 2 O and EtOH, and hardly reaches 20% for H 2 and O 2 . More interestingly, the size dependence of Q I (open symbols),