Interaction between particles with inhomogeneous surface charge distributions: Revisiting the Coulomb ﬁssion of dication molecular clusters

An analytical solution describing the electrostatic interaction between particles with inhomogeneous surface charge distributions has been developed. For particles, each carrying a single charge, the solution equates to the presence of a point charge residing on the surface, which makes it particularly suitable for investigating the Coulomb ﬁssion of doubly charged clusters close to the Rayleigh instability limit. For a series of six separate molecular dication clusters, centre-of-mass kinetic energy releases have been extracted from experimental measurements of their kinetic energy spectra following Coulomb ﬁssion. These data have been compared with Coulomb energy barriers calculated from the electrostatic interaction energies given by this new solution. For systems with high dielectric permittivity, results from the point charge model provide a viable alternative to kinetic energy releases calculated on the assumption of a uniform distribution of surface charge. The equivalent physical picture for the clusters would be that of a trapped proton. For interacting particles with low dielectric permittivity, a uniform distribution of charge provides better agreement with the experimental results.


I. INTRODUCTION
The behaviour of multiply charged clusters and nanodroplets close to the Rayleigh instability limit is a phenomenon of broad relevance and practical importance for a wide range of processes, including electrospray ionization 1,2 and the separation of carbon nanotubes 3 . Identifying the point at which highly charged, finite collections of atoms and molecules become unstable has been the subject of numerous experimental and theoretical studies, [4][5][6][7][8][9][10][11][12][13][14] but it is only recently that detailed patterns of charge separation close to the Rayleigh instability limit have been established. 15,16 Molecular systems for which specific fragmentation patterns have been mapped include: (H 2 O) 2+ n , (NH 3 ) 2+ n , (CH 3 CN) 2+ n , (C 5 H 5 N) 2+ n , (C 6 H 6 ) 2+ n ; 16 (C 6 H 6 ) z+ n , (CH 3 CN) z+ n and (C 4 H 8 O) z+ n with z = 3 and z = 4; 17 and (NH 3 ) z+ n clusters with z ≤ 8, 18 . Precise quantitative data are available for a range of doubly charged clusters, where Coulomb fission has been shown to result in asymmetric fragmentation and the formation of two singly charged clusters. 15,16 For these latter examples, an accurate determination of the centre-of-mass kinetic energy release can provide a measure of the Coulomb repulsion experienced by the two fragments as they separate. The exact mechanisms leading to charge separation is still the subject of intense interest, but experiments show that the outcome of the fragmentation step depends on the size and, to some extent, the composition of a cluster. 16 This combined experimental and theoretical study of cluster dication fragmentation, 16 showed that the barrier to charge separation depends strongly on the dielectric constant (or polarisability) of a cluster. In these examples, Coulomb fission has been a) Electronic mail: fav@triniti.ru. b) Electronic mail: Anthony.Stace@nottingham.ac.uk c) Electronic mail: Elena.Besley@nottingham.ac.uk modelled using a theory of electrostatic interactions that has been developed for charged particles of dielectric materials 19 . The theory assumes that, for each fission product, the single free charge is uniformly distributed over the surface of the latter; an approach that has been more widely adopted in the literature (see papers [20][21][22][23] and references therein). A more general solution to this problem has been proposed by Munirov and Filippov 20 .
Whether the assumption of a uniform positive or negative charge distribution over the surface of a particle is an appropriate description is open to speculation. The alternative is to assume that charge is localized, possible in the form of, for example, a proton in the case of a positive charge or as a surface quantum state in the case of an excess electron. There are a number of examples where experimental data have been interpreted in terms a localized charge. For example, the anomalous behaviour of (H 2 O) 2+ 21 has been interpreted as H 3 O + localized at the centre of a water cluster [24][25][26][27][28][29] and a similar picture of a localized central charge has been proposed to explain optical spectra recorded from xenon cluster ions, Xe + n . 30 For these examples, which are essentially spherical clusters with point charges at their centres, Gauss's law states that such a configuration generates a potential equivalent to that of the charge being uniformly distributed across the surface. In the case of singly charged anionic clusters, a combination of separate solvated and surface-bound states has been used to account for experimental electron photodetachment spectra. [31][32][33] For multiply charged clusters, theoretical studies of the stability of rare gas clusters have been model on collections of surface-bound charges each occupying a position determined predominately by Coulomb repulsion. 10,14 In contrast, molecular dynamics simulations of mixed water/ammonia/methanol clusters show that charge carriers in the form of NH + 4 are distributed throughout. 34 As an alternative to the assumption of a homogeneous distribution of free charge over the surface of the interacting par-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
ticles, this paper considers the case where for two charged, dielectric, spherical particles originating from Coulomb fission there are inhomogeneous distributions of free charge on their surfaces. The charge distributions are described by δfunctions of angular variables and a general solution is derived for a point charge of one elementary unit residing on the surface of each particle of arbitrary radius. The proposed solution is ideally suited to the study of processes involving the coalescence and Coulomb fission of charged clusters, liquid droplets and other nanoscale particles, where the particles concerned carrying one or few elementary charges (in absolute value). It may also be relevant for the study of electrostatic interactions between patchy colloidal particles 35 . The model is tested against fragmentation data recorded of a series of doubly charged molecular clusters composed of benzene, tetrahydrofuran (THF), pyridine, ammonia, acetonitrile, and water molecules. In particular, values predicted for the Coulomb repulsion energy between separating fragments are compare with accurate kinetic energy release measurements extracted from experimental data.

II. ELECTROSTATIC INTERACTION ENERGY OF INHOMOGENEOUSLY DISTRIBUTED FREE CHARGES
To account for an inhomogeneous distribution of free charge residing on the surface of the interacting particles, we assume a δ -function distribution of charge which is dependent on angular variables and the problem is solved following an approach proposed by Munirov and Filippov 20 .
A. Geometry of the problem and expansion of the electrostatic potential We consider two spherical particles with radii a 1 and a 2 , and charges q 1 and q 2 , which are generally non-uniformly distributed over their surfaces; the particles have dielectric constants ε 1 and ε 2 and they are placed in a uniform dielectric medium with permittivity ε. We introduce a Cartesian coordinate system such that the z-axis is directed along a line connecting the centres of the particles (see Fig. 1). The choice of the plane xz remains arbitrary.
In bispherical coordinates, [36][37][38] denoted as (ξ , η, ϕ) (see Fig. 1): where a = a 1 sinh ξ 1 = a 2 sinh ξ 2 , ξ 1 and ξ 2 , are the coordinate surfaces, which coincide with the surfaces of the particles, such that where R is the distance between centers of the particles: R = z 1 − z 2 = a (coth ξ 1 + coth ξ 2 ) = a 1 cosh ξ 1 + a 2 cosh ξ 2 , z 1 and FIG. 1. Geometry describing the interaction between two particles with radii a 1 and a 2 in a bispherical coordinate system (ξ , η, ϕ). z 2 are z-coordinates of centers of the particles: The Lame coefficients h ξ , h η and h ϕ in bispherical coordinates are determined from the expressions 38 : The electrostatic interaction between particles in a uniform dielectric is determined by the Laplace equation ∆φ = 0, which in bispherical coordinates can be solved by separating variables with the introduction of a new quantity φ (ξ , η, ϕ) = ψ(ξ , η, ϕ) cosh ξ − cos η, here φ is the electrostatic potential. Bound solutions of the Laplace equation in bispherical coordinates for regions inside (ψ I and ψ II ) and outside the particles (ψ III ) can be represented as follows This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
where P m (cos η) is the Legendre polynomial function, and the potential expansion coefficients A ±m , B ±m , C ±m and D ±m are determined in Appendix A. The electrostatic potential satisfies the following boundary conditions 39 where σ 1 , σ 2 are the surface densities of free charge residing on the particles, which are generally functions of η and ϕ.
The distribution of surface charge is usually given in a spherical coordinate system with an origin at the centre of the i-th particle as Here, θ i is the polar angle of a point on the surface of the i-th particle with a pole at its center, i.e. the angle between the radius vector to the point in question and the positive direction of the z-axis as denoted by θ 10 and θ 20 in Fig. 1, and ϕ is the azimuthal angle of this point. Decomposition of the surface charge represented by Eq. (7) in the bispherical coordinate system has been proposed by Munirov and Filippov 20 : where and the decomposition coefficients b nm i, are defined by the following expression In the case of a uniform surface charge distribution in Eq. (7), n = 0 and using Eq. (10) we obtain For the case of an axially symmetric charge distribution, m = 0, such that and from Eq. (10) we obtain B. Point charges represented by δ -functions on the surfaces of particles In a bispherical coordinate system, charge located on the surface of a sphere at the position defined by the coordinates η i0 and ϕ i0 (i = 1, 2) can be expressed as The value of σ i0 is determined by integration of the charge distribution (14) over the surface of the ith particle: (15) leading to the following expressions and Equation (17) can be expressed in the form of Eq. (8). Taking and using (8), S i (η, ϕ) can be represented as which, using (17) gives Multiplying Eqs. (19) and (20) by cos (mϕ) P m n (cos η) or sin (mϕ) P m n (cos η) and integrating over sin ηdηdϕ using This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
Eqs. (C5) and (C6) we obtain from (19): cos mϕ sin mϕ P m n (cos η) sin ηdηdϕ = = 2π 2n + 1 1 2 )|ξi| (21) and from (20) 2π 0 π 0 S i (η, ϕ) cos mϕ sin mϕ P m n (cos η) sin ηdηdϕ = = q i a 2 (cosh ξ i − cos η i0 ) 3/2 cos mϕ i0 sin mϕ i0 P m n (cos η i0 ). (22) Eqs. (21) and (22) can be combined to give For η = 0, π and ϕ = 0, 2π we should consider the limits of η → 0 +0 , π −0 and ϕ → 0 +0 , 2π −0 . Since choice of the xz plane is arbitrary, it is chosen such that ϕ 10 = 0. In this case, for decomposition of the surface charge of the first particle, terms with negative index "−m" are equal to zero. Note that for fixed η i0 , the charge moves across the surface as the interparticle distance changes. To avoid this, it is necessary to specify an angle in a spherical coordinate system with a pole at the centre of each particle. Therefore, we have cos η 10 = cosh ξ 1 − sinh 2 ξ 1 cosh ξ 1 + cos θ 10 , cos η 20 = cosh ξ 2 − sinh 2 ξ 2 cosh ξ 2 − cos θ 20 (25) where θ 10 , θ 20 are polar angles in a spherical coordinate system with poles at the centres of the first and second particles, respectively. For the case of several point charges present on the surface, it is necessary to sum expression (23) for each charge. From Eqs. (8) and (23) we obtain the following expression for the distribution of density of free charge cos mϕ i0 cos mϕ + sin mϕ i0 sin mϕ P m (cos η)P m (cos η i0 ) , i = 1, 2. (26) In the case of an axially symmetric system, we proceed from the relationship and find that Therefore, recalling that the δ -function is a functional, we give the distribution of surface charge a convenient form From Eqs. (8) and (29), proceeding with similar steps to those used to obtain Eq. (26), we derive the distribution of axially symmetric free charge Taking into account Eq. (C1), this leads to Eq. (29).
An alternative method of representing the surface charge density in the presence of point charges defined by a δfunction is shown in Appendix B. This method has been derived in spherical coordinates. Our numerical tests show that this alternative approach is less effective and requires significantly more computational resources. Therefore, in this paper all calculations have been performed using the method described in this section. In numerical solutions to the problem, a finite number of terms is considered up to a given value: = max , where the accuracy of representation (29) using finite terms can be examined using expression (C3).

C. Density of the total surface charge
The distribution of the total charge density on the surface of particles is determined from a discontinuity of the normal component of the electric field From the boundary condition (5) it follows that Using Eqs. (2), (3), (4) and (32), we obtain from Eqs. (31) (33) This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
For the axially symmetric problem, Eqs. (33) and (34) are reduced to Note from Eqs. (33)(34)(35) that the distribution of total charge includes the quantities D e ( + 1 2 )ξ1 and C e ( + 1 2 )ξ2 , and this is unlike the expression for the electrostatic force derived in the next section, where the force is defined simply from products of the expansion coefficients of the electrostatic potential. In the case of a uniform distribution of free surface charge, this does not present any difficulties, as F ±m used in expression (A1) and defined by Eqs. (A5), decrease exponentially with , and accordingly D e ( + 1 2 )ξ1 and C e ( + 1 2 )ξ2 also decrease exponentially. In the case of a δ -like distribution of free charge, F ±m in (A5) is expressed as (the upper line in curly brackets refers to the index "+m", the lower line refers to "−m") the values are of order O(1) for all . We next consider the simple case of an axially symmetric system. In this case, from (A1) -(A5) we obtain (denoting where For sufficiently high values of , the diagonal elements of the matrices A , B and C become negligible because, even with δ -localized charges, they decrease with as e − ξ 1 or e − ξ 2 . Therefore, at sufficiently large ≥ L max the expansion coefficients of the potential for each particle cease to depend on their own charge. In this case, if we introduce new variables for ≥ L max systems, equation (37) takes the form For a uniform distribution of free surface charge, Eq. (41) reduces to and for δ -type charges localised along the z-axis gives In the case of δ -localized charges, one can see from Eqs. (43) and (45) that with increasing the new coefficients c and d remain of the order of O(1), and calculation of the total charge density distribution is not trivial. For ≥ L max , matrices (38), (39) and (40) take the following form (for new variables) (C ) 11 = 0, This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
Using the matrix sweep method described in Appendix A requires the introduction of a direct sweep matrix α +1 (A7), which for ≥ L max also takes on a diagonal form In the limit of → ∞ Eq. (49) leads to (50) and from the condition α +1 = α we obtain For the vector β +1 from Eq. (A7) one can obtain which in the limit of → ∞ gives These expressions show that in the case of a uniform distribution of surface charge, β +1,1 = β ,1 e −ξ 2 and β +1,2 = β ,2 e −ξ 1 , i.e. as increases, they tend to zero. For the case of a point (or circumferential distribution) of δ -localized charge lying on the axis, we have that For large , when δ -localized charges are located on the zaxis (i.e., at η i0 = 0 or π) coefficients β +1,i tend to a constant value determined by the following relations We can now write the asymptotic solution of the system (A7) (which we denote as c and a ) in the following form For a sufficiently large number = N (N can be equal to infinity), solution (57) takes the form Analysis of the asymptotic expression (C4) shows that constancy (in absolute value) of the ratio of two neighboring Legendre polynomials P /P −1 is only possible if θ i0 = 0 and θ i0 = π (repeated cycles are at θ i0 = π p/q, where p and q are integers).
Finally, we consider the case of point charges lying on the z axis. In this case, cos η i0 = 1, P (cos η i0 ) = 1, or cos η i0 = −1, P (cos η i0 ) = (−1) , and the coefficients b i, cease to depend on and become constant. From Eq. (58) we obtain and from Eq. (56) we obtain These expressions show that for sufficiently large > L max , the coefficients b 1 and b 2 do not depend on . For high multipole expansion terms of the total surface charge shown in (35) we find This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

D. Electrostatic interaction force
The Maxwell tensor of tensions can be used 39 to calculate the electrostatic force applied to a dielectric body. In reference to Fig. 1, we perform calculations of the electrostatic force for a particle located on the positive side of the z axis, such that the repulsion force due to the presence of a second particle is positive and the attraction is negative where e are orthonormal basis vectors. For Cartesian components of the interaction force, Munirov and Filippov 20 found Here it is assumed that C −m = C m and D −m = D m for m = 0. For the case of uniformly charged particles, only the zcomponent of the electrostatic force is non-zero, as follows from Eqs. (67)-(69), and has the following form

III. EXPERIMENTAL SECTION
As discussed in earlier publications [16][17][18] , observations on the fragmentation patterns of multiply charged molecular clusters have been undertaken on an apparatus that combines a high resolution reversed geometry mass spectrometer (VG Analytical ZAB-E) with a pulsed supersonic cluster source. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
Neutral molecular clusters from each of the liquids were ionized by 70 eV electrons and the ion beam extracted from the ion source at a potential of between +5 and +7 kV into the flight tube of the mass spectrometer. Doubly charged cluster ions, M n 2+ , with values of n known to be close to the Coulomb instability limit have been mass-selected using a magnet, and the ionic products of Coulomb fission in the second field-free region (2ffr) between the magnet and an electrostatic analyser (ESA) were identified by scanning the voltage on the latter. This linked-scan procedure provides a massanalysed ion kinetic energy (MIKE) spectrum 40 , which can be used to identify ionic fragments according to their laboratoryframe kinetic energy. In addition, the energy spread in a peak can be related to the centre-of-mass kinetic energy released during Coulomb fission 40 ; in effect, this is measure of the repulsion experienced by the two positively charged ions as they separate. To detect the principal charged fragment from the fission of a doubly charged cluster, the ESA was scanned to record ionic fragments with laboratory-frame kinetic energies of between 5-7 and 10 keV. Within this energy range, there are no background ion signals from other processes, such as the loss of neutral molecules, which means the very weak signals that arise from Coulomb fission can be recorded without interference. However, this approach does mean that only the largest of the charged fragments is detected and therefore, the size of the smaller fragment has to be inferred from mass and charge balance. Figure 2 gives an example of a kinetic energy spectrum recorded following the Coulomb fission of [(H 2 O) 37 ] 2+ , and as noted above, the key measurement in these experiments is the kinetic energy released as a consequence of Coulomb repulsion between the two like-charged fragments as they separate. Due to a combination of high kinetic energy release and comparatively light mass, the peak profiles recorded from the fragmentation of molecular dications show the presence of a sequence of pronounced dish-shaped peak profiles 16,17 . These arise from the preferential transmission of fragment ions that are either strongly forward or backward scattered with respect to the laboratory frame of reference. Whilst it is possible to extract a value for the average kinetic energy release (KER) from the full width at half maximum (FWHM) for each peak assigned to a particular fission process 40 , this approach does have limitations. It relies heavily on the quality of the experimental data and is very sensitive to how accurately the width (∆E) of a profile can be measured (KER is ∝ ∆E 2 ) 40 ; hence, poor signal-to-noise levels on the edges of peaks can lead to errors in KER values.
To improve the accuracy of energy release measurements, a method for calculating peak profiles proposed by Beynon and co-workers 41,42 has been adopted in the form of a Monte Carlo simulation 17 . From random values for the kinetic energy release, centre-of-mass velocities for the fragments are calculated on the assumption that, in the centre of mass frame, the scattering of ions is equally probable in all direction. These velocities are then transformed to the laboratory-frame as two components, v z , which determines whether or not a fragment ion will pass through the final slit on the mass spectrometer, and v xy , which determines how rapidly a fragment ion will pass through the electrostatic analyzer and reach the detector 18,41 . Since the position in the flight tube where fission occurs also influences the probability of an ion passing through the final slit, it is assumed that ions have equal probability of fragmenting at all positions along the 2ffr. A total of 10 6 simulations were run for each set of conditions and for those ions calculated to have reached the detector, their centre-of-mass kinetic energies were transformed into a laboratory-frame peak profile. The simulation of peak profiles was found to be sensitive to two centre-of-mass kinetic energy variables: the minimum kinetic energy, E min , which primarily contributes to the centre of a profile and the maximum allowed kinetic energy, E max , which determines the width of the profile and the separation between the wings. The results shown in Fig. 2 were all calculated using E min = 0.7 eV together with a uniform spread in kinetic energies up to E max = 1.1 eV; from visual inspection these values are probably accurate to ±0.05 eV. The calculated average kinetic energy (KER) is 0.89 eV. The most significant results from these experiments are the values determined for E max , as these represent the maximum repulsion energy experienced by singly charged fragments as they separate. In terms of the calculated interaction energies, E max should be less than or equal to the height of a Coulomb barrier. A value for E max less than the height of a particular barrier is possible if some of the kinetic energy is assumed to have been dissipated into internal energy as the fragments separate. Such behaviour would result in the evaporation of single molecules and further scatter the fragment ions, leading to an in-filling of the dish shaped peak profiles.
A series of measurements have been taken for six separate molecular systems, and these fragmentation processes are summarised in Table I. In each case a single fragmentation step is taken to be representative of all the profiles measured for a particular, mass-selected molecular dication. Although there are subtle differences between profiles, the accuracy of both the measurements and the simulations do not warrant further detail. For several examples, kinetic energy data for the same molecular system have been recorded using measure-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
ments undertaken at two ion source potentials, 5 kV and either 6 or 7 kV, and in each case the peak profiles could be accurately simulated using the E min and E max information given in Table II. IV. NUMERICAL RESULTS, COMPARISON WITH EXPERIMENT, AND DISCUSSION Table II presents experimental data from measurements on the kinetic energy associated with Coulomb fission of doubly ionised clusters, together with computational predictions for the corresponding electrostatic interaction energy barriers. Also shown are the dielectric permittivities of the fluids forming the molecular clusters and the sizes of the fragments into which the clusters decay. The latter have been determined from the masses and densities of the fluids concerned. Pairs of particles are subsequently identified by the numbers given in the first column of Table I.
The reduced electrostatic interaction force between two particles has been calculated as a function of surface-tosurface separation and has been scaled by the Coulomb interaction force between two point charges located at the centreto-centre inter-particle distance R. Five variants are considered which differ in the type and location of the surface free charge: variant I with the homogeneous distribution of surface charge and variants II-V with the δ -charges located along the z-axis as shown in Fig. 3. Figures 4(a)-(d) show the reduced electrostatic force calculated, with reference to Table I, for pairs of charged particles involved in the following fragmentation processes 1 (benzene), 9 (acetonitrile), 11 (water), and 12 (water). These cases differ in the dielectric permittivities of the particles (the corresponding liquids) and the sizes of fragmented pairs. Note that in the case of a uniform distribution of surface charge (variant I) no attraction between fragments is observed at all surface-to-surface separation distances if the fragments have low dielectric permittivities, such as benzene (see Figure 4(a)), and/or the fragments are not too dissimilar in size, as in process 12 (see Figure 4(d)). For variant I, at very short separation distances of L ≈ 10 −2 nm (L = R − a 1 − a 2 ) the attraction between fragments increases with the increase in the ratio of fragment radii a 1 /a 2 . Comparison of Figures 4(b) and 4(c) shows the stronger attraction in process 9 (a 1 /a 2 =1.60) than in process 11 (a 1 /a 2 =1.39).
It was found in 20 that generally the dependence of the interaction force on the charges and radii of particles can be described by the following formula: where f is an unknown function of these arguments, which varies from 1 for large separation distances, L, to 0 at the maximum point of the interaction potential. According to this formula, for a pair of identical molecules with constant charges the interaction force increases, at large separations, if the radius of one particle (a 2 ) is decreased whilst keeping the radius of another particle (a 1 ) constant. However, at short separations, the increase in the ratio of particle radii leads to stronger attractive charge-induced polarisation forces thus making the effect of particles sizes on the interaction energy barrier difficult to predict. These considerations are supported by data presented in Table II showing that with an increase in the ratio of particle radii, a 1 /a 2 , the height of the barrier (the value of U max ) decreases in pairs 4 and 5, 6 and 8, 9 and 10, but goes up in pairs 1 and 2, 11 and 12 (note that in the latter pair both radii are changed).
If the point charges located on the surfaces face one another (variant IV), all pairs experience repulsion. If the point charges are located on opposite sides (variant II), the interaction at very short distances is at its weakest and, for sufficiently large values of the dielectric permittivity, the interaction force is very close in value to the case of a uniform distribution of surface charge (Figs. 4(b)-(d)). For an arrangement of δ -point charges localised on the same side of the particles (either on the right as in variant III or on the left as in variant V) a repulsive barrier becomes attractive for all pairs at short separation.
For pairs of particles involved in fragmentation processes 1 (benzene) and 9 (acetonitrile), the dependence of the elec-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  trostatic interaction energy on the surface-to-surface separation is shown in Fig. 5, and values for the Coulomb barriers (U max ) are summarised in Table II. If the electrostatic interaction energy passes through a maximum as the separation distance decreases (as, for example, in variants III and V), its value is shown in Table II in normal font, however, if the maximum value for the interaction energy is reached at the shortest separation distance it is highlighted in bold.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347
To establish why there is a switch in electrostatic behaviour from repulsion to attraction in particles carrying the same sign of charge, we consider the distribution of total surface charge at short separation distances L. Figure 6 presents the case for a uniform distribution of free surface charge at L = 10 −1 nm for pairs of particles involved in fragmentation processes 1 and 9. It can clearly be seen that for the charged clusters involved in process 1 (benzene, low dielectric constant of ε = 2.28) the density of the total surface charge is positive at all angles (i.e. everywhere on the surface of the fragments), such that for this pair, attraction is not observed at any of the separation distances under consideration. However, for the pair involved in process 9 (acetonitrile, ε = 37.5), a significant fraction of negative total surface charge is accumulated on the larger par-ticle in proximity to the second particle. It is the interaction between this negative charge and the neighbouring particle, which has positive total charge density everywhere on its surface, that leads to attraction between the cationic fragments. For cases where each particle has a δ -like positive point free charge, the distributions of total surface charge at L = 10 −1 nm for fragmentation processes 1 and 9 for variants II and IV are shown in Fig. 7, and for variants III and V are shown in Fig. 8. From these plots the following interesting observations can be made.
• In variant II, the surface point charges are furthest apart, and for low values of ε 1 and ε 2 the sign of the total surface charge is positive at all angles, therefore, in Fig. 4(a) we do not observe any attraction. If the interacting particles have high values of ε 1 and ε 2 , negative charge is acquired on that part of the surface of the larger particle that is closest to the second particle, therefore, in Fig. 4(b) we observe attraction at short separation distances.
• In variant IV, the positive point charges face one another and no attraction is observed for any of the pairs, This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  Table I. The calculated force is divided by the corresponding Coulomb interaction force between two point charges. The curves are labeled by the variant number of the free charge distribution.

PLEASE CITE THIS ARTICLE AS
despite the fact that near the point charges there exists an area of induced negative charge. This more delocalised region of negative charge resides at a greater distance from the point charge of the neighbouring particle than the corresponding point charge, and the overall effect is that the negative charge only slightly weakens any inter-particle repulsion. This effect leads to the appearance of a clearly defined minimum in the dependence of the electrostatic force on separation distance L for high values of ε 1 and ε 2 (Figs. 4(b)-(d)). Note that in fragmentation process 9, the total charge on the smaller particle always remains positive.
• For the case of point charges facing the same direction, to the right (variant III) or to the left (variant V), there is for all pairs a strong attraction at short distances. Here, negative charge can be induced on both particles, if the values of ε 1 and ε 2 are sufficiently high. This is the case shown in Fig. 8(b) for variant V, where the smaller particle lies between the point charges,but it also applies when either particle lies between the point charges.
• At the location of the point charge, an increase in the density of total charge of the same sign as that of the point charge is observed in almost all cases. The only exception is variant IV where the point charges are at the shortest separation from each other. In variant IV, if the dielectric permittivities of the interacting particles are small, the sign of the total charge (at the location of the point charge) is opposite to the sign of the point charge for both particles; if the dielectric permittivities are large, this only occurs for the particle with larger radius (see Fig.7(b)).
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The latter statement refers to a rather unusual, and not selfevident circumstance, which we will now consider in further detail. We calculate the angular distribution of the total surface charge for different values of inter-particle separation, L, as shown in Figs. 9-11. For a uniform distribution of surface free charge, the angular distribution of the total charge as a function of inter-particle distance (Fig. 9) is well-known. We only note that on the smaller particle, in the region closest to the neighboring particle, an increase in the density of total charge of the same sign as the point charge is observed ( Fig. 9(b)). Figures 10 and 11 show that at the location of the point charge an increase in the density of total surface charge of the same sign persists even at large inter-particle separation where their mutual influence is negligible. Consider the case of a single particle containing a point charge on its surface. From Eq. (C9) in the limit of R → a 1 , the total charge distribution on a sphere in the case of δ -localised charge is given by

PLEASE CITE THIS ARTICLE AS
where θ is the polar angle in spherical coordinates with a pole at the centre of the dielectric sphere and an axis directed to the point charge, i.e. the point charge is located at θ 0 = 0. Using Eq. (C1) we transform Eq. (72) to We see from Eq. (73) that in the limit of θ → θ 0 (note that for θ 0 = θ , δ (cos θ 0 − cos θ ) = 0 and P n (cos θ ) > 0 for all n at θ < π/2) the sign of σ t,δ is the same as the sign of q. This explains the appearance of an additional contribution to the total charge that has the same sign as the point charge close by. From the last term on the r.h.s. of this expression it is also evident that the environment reacts to the presence of the surface point charge such that the "height" of the δ -function decreases.

V. CONCLUDING REMARKS
A general solution to the problem of calculating the electrostatic interaction between particles with inhomogeneous dis-This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. tributions of free surface charge has been presented and compared with the more established case of a uniform distribution of free surface charge. Point charges on each of the particles have been described by a δ -function of angular variables, and electrostatic energy barriers have been calculated over the range of particle-particle distances, 10 −2 ≤ L ≤ 10 3 nm, regions where the effects of induced interactions are found to be at their greatest. When compared with results for a uniform distribution of charge, the calculations reported in Table II show that the exact location of the point charge has a marked influence on the value of any Coulomb barrier in the electrostatic energy. These barriers can affect the behaviour of charged particles in two ways: first, there is the coalescence of charged particles and which of the particle orientations offers the most facile pathway to facilitate that process. A second aspect to the work concerns the fragmentation of multiply charged particles and here the barrier will contribute to the process of charge separation. In this work we have addressed the latter topic with examples taken from the fragmentation of dication molecular clusters. As has been discussed in earlier publications on the fragmentation of dication clusters, [16][17][18] there are other contributions to the energy barrier experienced by a dication at the onset to fragmentation, and these will come from the breaking of physical and/or This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The total surface charge density on the larger particle (a) and the smaller particle (b) as a function of the polar angle θ 1 or θ 2 and the surface-to-surface separation distance L for the fragmentation process 9 (acetonitrile, ε = 37.5). Free charge is uniformly distributed on the surface of each particle.

PLEASE CITE THIS ARTICLE AS
(a) (b) FIG. 10. The total surface charge density on the larger particle (a) and the smaller particle (b) as a function of the polar angle θ 1 or θ 2 and the surface-to-surface separation distance L for the fragmentation process 1 (benzene, ε = 2.28). Each particle has a δ -like positive point free charge located at θ 10 = θ 20 = 0 (variant III).
chemical bonds, which in the examples given will most probably take the form of either van der Waals or hydrogen bond interactions. The calculations discussed here are only concerned with events taking place after these bonds have been broken and the individual particles start to separate; at that point they are subject to the influence of a purely repulsive Coulomb interaction, which depending on how polarizable the particle are, can be moderated by the effects of short-range, attractive charge-induced multipole interactions.
Accurate experimental measurements of the Coulomb barrier in the form of values for the kinetic energy release following charge separation have been presented for six molecular dications. Table II summarises these data together with calculated values of the barrier for different surface locations of the two point charges as the particles separate (columns II-V). Also given are barriers calculated for a uniform distribution of free surface charge. As noted earlier, these calculated barrier heights are to be compared with experimental results determined for E max , which represent maxima in the Coulomb energy as recorded from peak profiles of the type shown in Fig. 2; the value of E max should be less than or equal to U max . As can be seen from Table II, for weakly polarisable clusters composed of either benzene or THF molecules, this latter requirement is most closely met with a uniform distribution of surface charge (variant I). As ε increases in value, the point charge orientation θ 10 = θ 20 = π (variant V) becomes a realistic alternative, particularly when the ±0.05 eV experimental error limit is taken into consideration. However, for the most This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. FIG. 11. The total surface charge density on the larger particle (a) and the smaller particle (b) as a function of the polar angle θ 1 or θ 2 and the surface-to-surface separation distance L for the fragmentation process 9 (acetonitrile, ε = 37.5). Each particle has a δ -like positive point free charge located at θ 10 = θ 20 = 0 (variant III). polar of the clusters (water), almost any of the variants could be appropriate. One reason for the change in behaviour is that the charge could become localised in the form of a proton or similar charge carrier, 43 and for at least two of the systems, (H 2 O) 2+ n and (NH 3 ) 2+ n , this is a distinct possibility.

PLEASE CITE THIS ARTICLE AS
If we examine what a inhomogeneous charge distribution might mean for the outcome of particle -particle collisions, for example when patchy colloidal particle coalesce, 35 then under those circumstances pathways with the lowest energy barriers might be expected to be more favourable. Using the data in Table II as an illustration, it can again be seen that there are differences depending on the dielectric constants of the materials involved. For low values of ε, variant II offers the lowest Coulomb barrier to coalescence, which bearing in mind this figure is for the gas phase, will be of the order of a few meV when the dielectric constant of any solvent is taken into consideration. Once the dielectric constant of the particles has a value of 25 or more, there is, as before, little to distinguish between barriers presented by a uniform charge distribution and variants II, III, and V. What might then influence events is the ability of the short-range attractive barrier to hold the particles in place.
Central to any discussion of charge orientation is the timescale over which competing events take place. For particles moving apart following fragmentation, those events are going to be the rotational period and the time taken for the particles to separate to a distance where the Coulomb potential no longer has a significant influence. Taking values for size, density and kinetic energy that are representative of the clusters given in Table II, the rotational period is estimated to be about 10 −10 s, and the time taken to separate to a distance of 10 nm is approximately 10 −11 s; thus, starting with a particular charge orientation (any of variants II-V), this would not be expected to change during the time taken for the particles to separate. From potential continuity conditions (5) we express the coefficients A ±m and B ±m in C ±m and D ±m . Then, using conditions (6) and the properties of the associated Legendre functions 44 , after simple algebra we obtain the following equations for determining the potential expansion coefficients 20 : where y = C ±m , D ±m T , This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5119347 and This series of equations is a system with a block threediagonal matrix.
Taking into account the fact that for = m from (A2) we have A ±m = 0 and putting y L max = 0 for a sufficiently large = L max , we will reduce the system of equations (A1) to the following form (N = L max − 1) C ±m m y m − B ±m m y m+1 = F ±m m , −A ±m y −1 + C ±m y − B ±m y +1 = F ±m , −A ±m N y N−1 + C ±m N y N = F ±m N , = m + 1, m + 2, . . . , N − 1.
(A6) Equations (A6) are solved by using the matrix sweep method 48 for each given value m = 0, 1, 2, . . . , N. The algorithm of the matrix sweep method has the form 48 (we omit the "±m" superscripts; the "−1" superscript denotes the inverse matrix): and the solution in this case will be: y N = β N+1 , y N−1 = α N y N + β N .

(A9)
Appendix B: Alternative method of representing the point charges on the surface of particles Let us consider the case of point charges located on the surfaces of spheres at points with coordinates θ i0 and ϕ i0 : The value of σ i0 is determined by an integral over the surface of the i-th particle (µ i = cos θ i ): From Eqs. (B1) and (B2) we have Using the relation (C2) for the expansion coefficients of Eq. (7) in spherical coordinates with a pole at the centre of the i-th particle we find Then, from (9) and (10) we find expansion coefficients, which are already represented in the bispherical coordinate system. In the case of an axially symmetric arrangement of point charges, we can proceed from the expression: then σ i0 are defined by expressions and we obtain Finally, we use Eq. (C1) and for the expansion coefficients of Eq. (7) in spherical coordinates with a pole at the centre of the i-th particle we obtain σ i,n = q i 4πa 2 i (2n + 1) P n (cos θ i0 ) . (B8) Note, that this relation can be obtained from Eq. (B4) by assuming m = 0. Furthermore, we find the required expansion coefficients in the bispherical coordinate system from Eqs. (9) and (13). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. A solution to the electrostatic problem for a system consisting of an uncharged dielectric ball and a point charge is defined by the equations (see 38  Here a 1 is the ball's radius, ε 1 is the ball's dielectric permittivity, q is the point charge, ε is the dielectric permittivity of the medium, R is the distance between the center of the ball and the point charge, r and θ are the radius and polar angle in spherical coordinates with a pole at the centre pf the ball and an axes directed to the point charge. For the total charge distribution on the ball from Eq. (31) we find σ 1,t = q εR 2 ∞ ∑ n=0 a 1 R n−1 n (2n + 1) (ε − ε 1 ) nε 1 + (n + 1) ε P n (cos θ ) . (C9)