Laurent Jutier

Céline Léonard

Dr Céline Leonard
email: celine.leonard@univ-mlv.fr

Ab initio study of the C 2 O + cation

Keywords: CCO, ab initio spectroscopy, electronic potentials, rovibrational levels, spin-orbit coupling

The molecular cation C 2 O + plays an important role in the chain reactions which take place in the interstellar clouds involving the monoatomic fragments C, C + , O, O + and diatomic C 2 , C + 2 , CO, CO + , O 2 , O + 2 from which C 2 O + can be formed. In particular, C 2 O + can be produced by the ion-neutral reaction of CO [1]. Furthermore, this cation is a link towards molecular cations of C n O + type, with n > 2. Some of the occuring mechanisms have already been identified twenty five years ago, using mass spectroscopy, by Schildcrout and Franklin [2].

In 1982, Haese and Woods [3] performed ab initio calculations on the electronic states of C 2 O + , and established that the most stable isomer is the linear [C-C-O] + . They presented a 2 Π for the ground state, which is in agreement with the more recent works performed by Maclagan and Sudkeaw in 1993 [4], who also calculated the harmonic vibrational frequencies of the 2 Π state.

In 1995, Cao and Tian [5] studied theoretically the reactions {C + 2 + CO} and {C + 2 +O}. They found that the reaction product is mostly C 2 O + in the 4 Σ -state in both cases and they computed harmonic frequencies for this state. But there is still an uncertainty about the nature of the ground electronic state.

The objective of the present ab initio study is to obtain an accurate relative position of the 2 Π and the 4 Σ -energies. In a first part, all the electronic states associated to the first dissociation asymptotes are studied, showing that the minima of the first 2 Π and 4 Σ -states are close to each other and are much lower than the other states. Then, the second part of this work is devoted to accurate electronic calculations of these both states in order to confirm the nature of the ground state. Moreover, using new ab initio potential energy functions, the rovibronic energies corresponding to these lowest electronic states are computed variationaly. The study is completed by an investigation of the spin-orbit coupling between the 2 Π and 4 Σ -states.

2 Study of the first electronic states

Asymptotes of dissociation

All the molecular electronic states belonging to the first three asymptotes of dissociation are studied. Relative energies of these asymptotes are already known from Moore [START_REF] Moore | Atomic Energy Levels[END_REF] for atomic fragments, Huber and Herzberg [START_REF] Huber | Constants of Diatomic Molecules[END_REF] for CO and CO + and Wang, Zhu and Yang [START_REF] Wang | [END_REF] for C + 2 . The corresponding states for linear geometries are reported in Table 1. The dissociation along CO shows a high density of states: 6 doublets, 6 quartets and 2 sextets states are distributed in the two first asymptotes separated by less than 8000 cm -1 . These asymptotes are associated to the first electronic states of C + 2 combined with O in its electronic ground state 3 P g . For this molecule, the ambiguity concerning the nature of the electronic ground state was raised very recently by Wang et al. [START_REF] Wang | [END_REF]. They found that the electronic ground state is X 4 Σ - g with an equilibrium bond length of 2.672 bohr and the first excited electronic state is a 2 Π u with R eq CC at 2.497 bohr. The dissociation along CC is by far less complicated. The first asymptote is due to the combination of C + and CO in their respective electronic ground states 2 P u and X 1 Σ + with R eq CO = 2.132 bohr [START_REF] Huber | Constants of Diatomic Molecules[END_REF]. The second and the third ones combine CO + in X 2 Σ + with R eq CO = 2.107 bohr [START_REF] Huber | Constants of Diatomic Molecules[END_REF] with C in 3 P g and

Linear cuts

Using the MOLPRO software [9], the energies of the molecular electronic states given in Table 1 are obtained at the Multi-Configuration Self Consistent Field (MCSCF) level [10,11] with the correlation consistent polarized Valence Triple Zeta (cc-pVTZ) basis set [12] leading to a total of 90 contracted s, p, d and f gaussian functions. The active space comprised all the valence molecular orbitals (OMs) and all the states associated to a given spin multiplicity are averaged alltogether with one set of optimised OMs. Figures 1 and2 present the linear cuts of the electronic states of Table 1 dissociating along CC and CO respectively. The computations were done in the C 2v symmetry group and the non-dissociating part of the molecule was not relaxed. In the first case, R CO is fixed arbitrarily at 2.15 bohr and in the second one R CC at 2.60 bohr, these values being close to the isolated diatomic bond lengths. The density of states in Figure 2 is as large as expected from Table 1 and the avoided-crossings occur for large values of R CO . In Figure 1 both states in the CC dissociation for large R CC is certainly linked to the interactions between the different fragments. In the case of the 4 Σ -state, the interaction between CO + (X 2 Σ + ) and C corresponds to an ion-quadrupole interaction which behaves like ∼ -1/R 3 CC . Whereas for the 2 Π state, CO + C + corresponds to a ion-polarisability interaction which behaves like

∼ -1/R 4
CC since CO is almost apolar [13].

Linear equilibrium geometries

The linear equilibrium geometries of the first doublet and quartet states have been optimised

and are reported in Table 2 with the corresponding energies and the predominant electronic configuration. The values confirm that the lowest 4 Σ -state seems to be the ground electronic state separated by only 0.011645 a.u., i.e. 0.32 eV from the lowest 2 Π state. This value is relatively close to the separation of 0.93 eV between the X 4 Σ - g and the 2 Π u states in C + 2 . Table 2 and Figures 1 and2 show also that the 2 Π state has a floppier minimum along R CC than the 4 Σ -state associated with a different value of R eq CC (in agreement with the equilibrium bond lengths of the correlating electronic states of the C + 2 fragment [START_REF] Wang | [END_REF]), whereas both states present similar values of R eq CO and similar steepness around the minimum along R CO . These equilibrium geometries are in agreement with previous ab initio values [4] and [5]. The cut along the CC dissociation presents a crossing of these two states for a value of R CC intermediate between the equilibrium positions of 4 Σ -and 2 Π. This behaviour seems typical of the CC bond length since it has been already found in the diatomic C + 2 by Wang et al. [START_REF] Wang | [END_REF] and in the larger system [HCCO] + by Papakondylis and Mavridis [14].

For each spin multiplicity, there is a difference between the geometry of the lowest state and the geometries of the excited states which have similar equilibrium geometries.

Bent cuts

The variation of the potential of the lowest doublet and quartet electronic states with the bending angle was investigated. The doublet and the quartet states have been computed independently using the C s symmetry for all the points. The variation of the energy of the doublet states with the bending angle is shown in Figure 3 and in Figure 4 for the quartet states.

The bond lengths were fixed at the equilibrium geometry of the lowest 2 Π state and of the lowest 4 Σ -state for the doublet and the quartet states respectively. The two Renner-Teller From the results of the previous section, the inclusion of more electronic correlation in the electronic treatment is necessary in order to confirm if the 4 Σ -state is the electronic ground state of C 2 O + . After a MCSCF step producing a reference wavefunction, the Internally Contracted Multi-Reference Configuration Interaction (MRCI) method [15,16] including the multireference Davidson correction [17] (MRCI+Q) is applied to the lowest 4 Σ -state and to each electronic component of the lowest 2 Π state separately. For the 2 Π state, the two Renner-Teller components are averaged together at the MCSCF level. The equilibrium geometry of these states was obtained using different basis sets. Table 3 shows that the equilibrium bond lengths and energies of both states decrease with the size of the basis set. But Table 3 shows also the variation of the energy difference ∆E = E eq ( 4 Σ -) -E eq ( 2 Π) with the basis. This difference is positive for all the basis set and seems also to converge at 2388 cm -1 for the cc-pV5Z basis set [12] (giving a total of 273 contracted s, p, d, f , g and h gaussian functions). This basis set will be used in order to determine the potential energy functions (PEFs) of both states in the following part of the study. the static correlation must be included in order to place correctly the 2 Π state with respect to the 4 Σ -. We checked that this relative order is the same at the MRCI level.

Determination of the PEFs

The MRCI+Q / cc-pV5Z electronic energy was computed for 33 geometries around the respective equilibrium geometry for both the a 4 Σ -and the X 2 Π electronic states. All the points (24 linear and 9 non-linear) were performed in the C s symmetry group.

For the doublet state, the geometries spread over a range such as: 2.56 ≤ R CC ≤ 3.50 bohr, 1.91 ≤ R CO ≤ 2.41 bohr and 150 • ≤ θ ≤ 180 • . For the quartet state: 2.17 ≤ R CC ≤ 3.10 bohr,

1.93 ≤ R CO ≤ 2.43 bohr and 150 • ≤ θ ≤ 180 • .
These points were fitted by a polynomial expansion up to the sixth order of R CC and to the fourth order of R CO and θ:

V (R CC , R CO , θ) = ijk C ijk (R CC -R eq CC ) i (R CO -R eq CO ) j (θ -π) k (1)
with 0 ≤ i ≤ 6, 0 ≤ j, k ≤ 4, 0 ≤ i+j +k ≤ 4 and k must always be even in order to respect the symmetry of the bending mode. The reference is chosen as the equilibrium geometry (defined as the geometry for which the first derivatives of the PEF are zero) of the given state. The 24 expansion coefficients, optimised in a least square fitting with a Root Mean Square (RMS) of 0.4 cm -1 , are given in Table 4 with the reference geometries. The coefficient C 002 , associated to the square of the bending displacement, of the two linear-linear Renner-Teller components of the X 2 Π state confirms that the lower component is 2 A ′′ as shown in Figure 3 for large values of θ.

The norm of the equilibrium dipole moment, the reference being the center of mass of the molecule, is given for both states in Table 4 as Figure 5 shows the 2D crossing between the X 2 Π and a 4 Σ -states for the linear geometries.

The crossing surface is almost independent of R CC and the minimum is at R CC = 2.574 bohr and R CO = 2.140 bohr, 2434 cm -1 above the X 2 Π minimum and very close to the equilibrium of the a 4 Σ -state (cf Table 4).

The harmonic and anharmonic constants are obtained from the PEFs derivatives at the minima [18,[START_REF] Mills | Molecular Spectroscopy[END_REF], and are given in Table 5 with the equilibrium rotational constants for both electronic states. The harmonic wavenumbers predict a Fermi coupling in the X 2 Π state between the modes 1 and 2. This Renner-Teller state presents also strong anharmonic constants when the bending is involved in the A ′ component. This could be linked to the avoided crossing occuring in this component for θ ≈ 95 • (cf. Figure 3). Moreover, the values of ω 1 and ω 2 are approximatively twice in the quartet state than in the doublet state, whereas ω 3 are identical in both states.

The definitions of the stretching normal modes are given in Table 6 showing that, for both electronic states, Q 3 corresponds almost completely to the stretching of the CO bond whereas Q 1 is mostly associated to the CC stretch with a small mixing with the CO stretch. The mixing between the two bond stretches is slightly stronger for the a 4 Σ -state than for the X 2 Π state in relation with the fact that the harmonic wavenumbers ω

Spin-orbit coupling

The spin-orbit interaction was calculated between the X 2 Π and a 4 Σ -states and between the 4 spin-electronic components of the X 2 Π state. The spin-orbit coupling can be treated perturbatively and the coupling terms are the matrix elements of the non-relativistic part of the Breit-Pauli operator ĤBP SO [START_REF] Berning | [END_REF]21] on the basis of the spin-electronic wavefunctions. The wavefunctions were obtained with the MCSCF method and the uncontracted s, p and d functions of the cc-pV5Z basis set. The z axis corresponds to the linear molecular axis and x-axis bisects the bond angle θ. For the X 2 Π state, the couplings between the different components can be determined from:

LS z = A ′ φ σ=1/2 S=1/2 | ĤBP SO |A ′′ φ σ=1/2 S=1/2 = iA (2)
A equals -22.52 cm -1 at the equilibrium geometry leading to a spin-orbit constant A SO = -2A = 45.05 cm -1 (cf. ref [21]). This positive value is in agreement with a normal multiplet splitting and is comparable with the spin-orbit splitting of 43.5 cm -1 between the 3 P g0 and 3 P g2 states of the carbon atom [START_REF] Moore | Atomic Energy Levels[END_REF]. The LS z integral has been calculated for all the geometries of the grid used for the X 2 Π state and the corresponding A values have been fitted using the polynomial expansion of equation 1 with the same restrictions on i, j and k. The corresponding coefficients are given in Table 7. The fitted function has been obtained with a RMS of 10 -8 a.u. The variations of A with the geometry are shown in Figure 6. The coupling terms between the X 2 Π and a 4 Σ -states are given by:

LS 1 x = A ′ φ σ=1/2 S=1/2 | ĤBP SO |A ′′ φ σ=1/2 S=3/2
(3)

LS 1 y = A ′′ φ σ=1/2 S=1/2 | ĤBP SO |A ′′ φ σ=3/2 S=3/2 (4)
with

LS 1 x = iLS 1 y . LS 2 x = A ′ φ σ=1/2 S=1/2 | ĤBP SO |A ′′ φ σ=-1/2 S=3/2 = -LS 1 x √ 3 ( 5
)
LS 2 y = A ′′ φ σ=1/2 S=1/2 | ĤBP SO |A ′′ φ σ=-1/2 S=3/2 = LS 1 y √ 3 ( 6
)
with also LS 2 x = iLS 2 y . As an example, a zoom of the region of the crossing between the X 2 Π and a 4 Σ -states including the spin-orbit coupling is shown in Figure 7 for R CO fixed at 2.123 bohr ≈ R eq CO (X 2 Π). The first avoided crossing occurs at R cc = 2.6402 bohr with a splitting of 5.77 cm -1 between the Ω = 1/2 components and the second avoided crossing at R cc = 2.6424 bohr, the Ω = 3/2 components being separated by 9.82 cm -1 . |LS 1 y | = 3.47 cm -1 for both crossings. In conclusion, the spin-orbit coupling between the two electronic states is small and we checked that this coupling is almost not geometry dependent. where 2S + 1 is the spin-multiplicity and {K} is the notation for the different values of K: K = 0, 1, 2, 3 gives Σ, Π, ∆, Φ vibronic levels.

Using the code RVIB3 [22], the rovibronic levels of the X 2 Π state on one hand and the rovibrational levels of the a 4 Σ -state on the other hand have been determined up to 2200 cm -1

and 5000 cm -1 respectively using the PEFs presented in Table 4. These two sets were treated independently.

For the doublet state, the internal spin-orbit interaction is included via the geometry dependant analytical representation of A (see Table 7) [21] and the rovibronic levels have been obtained

for J = P = 1/2, 3/2, 5/2 and 7/2, where J is the quantum number associated with the projection of the total angular momentum on the molecular axis. For the quartet state, the internal spin-vibrational coupling is not taken into account at the present time in the code RVIB3 and the rovibrational levels have been obtained for J = K = 0, 1, 2, 3 and 4.

For both states, the contraction scheme involves 25 harmonic oscillator eigenfunctions for each stretching mode and 41 Legendre polynomials for the bending mode as basis functions. The primitive set of integration points comprised 46 and 52 points for the stretching and the bending modes respectively. The resulting energies are given in Tables 8 and9. The energy reference is the zero rovibrational energy of the corresponding electronic state. The density of levels is high for the X 2 Π state due to the spin-orbit and the Renner-Teller couplings, and the assignment of the levels above 2200 cm -1 becomes hazardous. The level noted * are tentatively assigned.

Such levels are more numerous among the κ levels (associated with the A ′ component) than among the µ levels (associated with the A ′′ component) ; this results from the Fermi coupling and the strong anharmonicity for the levels associated to the A ′ component as predicted above.

The vibrational levels of the quartet state are more separated in energy and they have been assigned up to 5000 cm -1 . Several levels are only tentatively labeled due to the Fermi coupling eventhough the harmonic wavenumbers do not emphasize this behaviour in the quartet state .

Since the minimum of the crossing surface (R CO , R CC , θ = 180 • ) of both states occurs at 2434 cm -1 above the X 2 Π minimum, the presented levels of the doublet state should not be perturbated by the presence of the quartet state. Moreover, the rovibrational levels of the a 4 Σ - state have been computed assuming that the X 2 Π and a 4 Σ -states do not interfer by spin-orbit coupling. This assumption is confirmed by the following development. The spin-orbit coupling between vibrational levels v and v ′ corresponding to different spin-electronic states se and s ′ e ′ can be seen, in the Born-Oppenheimer approximation, as the quantity:

SO vv ′ = Ψ se Ψ e v | ĤSO |Ψ s ′ e ′ Ψ e ′ v ′ = Ψ e v |Ψ e ′ v ′ Ψ se | ĤSO |Ψ s ′ e ′ (8)
The first integral is the Franck-Condon factor. This factor is assumed to be rather large (but never larger than 1) for the energy levels of the a 4 Σ -state situated below 5000 cm -1 from the zero rovibrational energy due to the fact that the X 2 Π and a 4 Σ -states intersect close to the equilibrium of the a 4 Σ -state (cf. Figure 5). However, the second integral corresponds to quantities computed in equations 3 to 5 and is associated with small values (≈ 3.5 cm -1 ).

Moreover, the zero vibrational energy of the quartet state lies 2864 cm -1 above the zero vibrational energy of the doublet state, reaching a high density of vibronic levels region in the doublet state. It is very difficult to predict the accurate position of the levels in such a situation.

Consequently the resonances between the X 2 Π and the a 4 Σ -states rovibrational levels have

Electronic transitions

In order to predict the possible transitions in UV absorption spectroscopy, the vertical transition moments and energies were computed at the equilibrium geometries of the X 2 Π and 4 Σ - states. The results are presented in Table 10. The first transitions start from the X 2 Π state and reach the first 2 ∆, 2 Σ + and 2 Σ -electronic states. These transitions correspond to relatively large transition dipole moments. They should be observed between 450 and 300 nm. They should be caracterized by long progressions involving the stretching mode which is associated to the variation of CC bond length, i.e. ν 1 due to the strong difference of the CC equilibrium bond lengths between the X 2 Π and the excited doublet states (cf. Table 2). The first transition associated with quartet states involves the a 4 Σ -and the first 4 Π electronic states and should be mesured at 232 nm, but the corresponding transition dipole moment is small (≈ 0. In conclusion, at room temperature, 3 2 kT = 308 cm -1 if T = 295 K, the absorption spectra should only show transitions from the X 2 Π state since the a 4 Σ -state is 2388 cm -1 above and the spin-orbit coupling is too small to couple these states efficiently. This study has given a general view about the spectroscopy of the C 2 O + cation where an ambiguity remained about the nature of the electronic ground state. The excited electronic states have been investigated showing that the lowest states are X 2 Π and a 4 Σ -. High level of configuration interaction calculations have allowed to guaranty that the electronic ground state is the Renner-Teller 2 Π state and that the first excited state is the 4 Σ -state with an equilibrium transition energy, T e , of 2388 cm -1 . In order to complete the study, the rovibrational levels of both states have been computed. It has been shown also that the spin-orbit coupling between these electronic states can be neglected. The knowledge of the electronic ground state has allowed to determine the vertical electronic transitions. These predictions will be useful for further experimental investigations of the spectroscopy of this cation.

Note : the ab initio points are available upon request. The origin of energy is the minimum of the X 2 Π state, and the step among curves is 1000 Table 1: First asymptotes of dissociation in both cases : {CO + C} + or {C 2 + O} + , and resulting molecular states. The energy reference is taken at the lowest energetic asymptote : { CO[ 1 Σ + ]+C + [ 2 P u ]}. The fragment energies are taken from ref. [START_REF] Moore | Atomic Energy Levels[END_REF], [START_REF] Huber | Constants of Diatomic Molecules[END_REF] and [START_REF] Wang | [END_REF]. R CC → ∞ : a : Maclagan and Sudkeaw [4] using MP2 / 6-31G** b : Cao and Tian [5] using UHF / 6-31G* Table 3: Equilibrium geometries (in bohr) and corresponding energies of the 2 Π and 4 Σ - electronic states calculated with MRCI+Q and cc-pVnZ basis sets (n = 2, 3, 4, 5 (D, T, Q, 5)) [12]. ) of the first rovibronic levels of the X 2 Π electronic state given with respect to the zero rovibrational energy at 1679.9 cm -1 . 9: Energies of the first rovibrational levels of the a 4 Σ -electronic state given with respect to the zero rovibrational energy at 2155.7 cm -1 .

cm -1 .
Energy (cm -1 ) Fragments States for the linear molecule [CCO] + 32345 CO + [ 2 Σ + ] C[ 1 D g ] 2 Σ + , 2 Π , 2 ∆ 22151 CO + [ 2 Σ + ] C[ 3 P g ] 2 Σ -, 2 Π , 4 Σ -, 4 Π 0 CO[ 1 Σ + ] C + [ 2 P u ] 2 Σ + , 2 Π R CO → ∞ : Energy (cm -1 ) Fragments States for the linear molecule [CCO] + 47821 CC + [ 2 Π u ] O[ 3 P g ] 2 Σ + , 2 Σ -, 2 Π , 2 ∆ , 4 Σ + , 4 Σ -, 4 Π , 4 ∆ 40312 CC + [ 4 Σ - g ] O[ 3 P g ] 2 Σ + , 2 Π , 4 Σ + , 4 Π , 6 Σ + , 6 Π
Basis 2 Π 4 Σ - R eq CC R eq CO E (u.a.) R eq CC R eq CO E (u.a.) ∆E (cm -1 ) cc-pVDZ 2.
2 Σ 1/2 µ 2 Σ 1/2 κ 2 Π 1/2 2 Π 3/2 v 1 v 2 v 3 Energy v 1 v 2 v 3 Energy v 1 v 2 v 3 Energy v 1 v 2 v 3
2 ∆ 3/2 2 ∆ 5/2 2 Φ 5/2 2 Φ 7/2 v 1 v 2 v 3 Energy v 1 v 2 v 3 Energy v 1 v 2 v 3 Energy v 1 v 2 v 3
4 Σ - 4 Π 4 ∆ 4 Φ v 1 v 2 v 3 Energy v 1 v 2 v 3 Energy v 1 v 2 v 3 Energy v 1 v 2 v
Σ - 2 Π 2 ∆ ↑ 2 Σ + ↑ 2 Σ -↓ 4 Π 2 Π↓ 2 Π 2 Σ + CO + [ 2 Σ + ] + C[ 1 D g ] CO + [ 2 Σ + ] + C[ 3 P g ] CO[ 1 Σ + ] + C + [ 2 P u ]
Σ - 2 Π 2 ∆ 2 Σ + ↑ 2 Σ -↓ 4 Π 2 Π 4 ∆ 4 Σ + 2 Σ + ↑ 4 Σ + 4 Π 6 Π 6 Σ + → CC + [ 2 Π u ] + O[ 3 P g ] CC + [ 4 Σ - g ] + O[ 3 P g ]

  Figures 1 and 2 near here **********************

*

  ********************* Table 2 near here **********************

  Figures 3 and 4 near here **********************

*

  ********************* Table 3 near here ********************** Contrary to the global study of the large set of electronic states with averaged molecular orbitals, at the MRCI+Q level, the electronic ground state is always the 2 Π state, showing that

  additional data. ********************** Table 4 and Figure 5 near here **********************

  1 and ω 3 are closer in the quartet state. This mixing is also caracterised by a large value of x 13 in the a 4 Σ -state. ********************** Tables 5 and 6 near here **********************

Figure 6

 6 Figure 6 near here **********************

Figure 7

 7 Figure 7 near here **********************

  not been investigated. ********************** Tables 8 and 9 near here **********************

  04 D). ********************** Table 10 near here **********************

  Figure captions:

Figure 1 :

 1 Figure 1: linear dissociation of CCO + along R CC (bohr) towards the first three asymptotes with MCSCF/cc-pVTZ computations.

Figure 2 :

 2 Figure 2: linear dissociation of CCO + along R CO (bohr) towards the first three asymptotes with MCSCF/cc-pVTZ computations.

Figure 3 :

 3 Figure 3: variation of the energy of the first doublet electronic states with the bending coordinate from MCSCF/cc-pVTZ computations.

Figure 4 :

 4 Figure 4: variation of the energy of the first quartet electronic states with the bending coordinate from MCSCF/cc-pVTZ computations.

Figure 5 :

 5 Figure 5: 2D contour plots of the MRCI+Q / cc-pV5Z PEFs for the X 2 Π and a 4 Σ -states for linear geometries. The full lines describe the a 4 Σ -state and the dashed lines the X 2 Π state.

Figure 6 :

 6 Figure 6: 1D cuts (a, b, c) of the variation of the spin-orbit coupling term A (see eq. 2) of the X 2 Π electronic state with R CC , R CO and θ. The other coordinates are respectively fixed to their equilibrium values: R CC = 2.8944 bohr, R CO = 2.1235 bohr and θ = 180 • .

Figure 7 : 14 Page

 714 Figure 7: Crossings between the X 2 Π and a 4 Σ -electronic states including the spin-orbit coupling. The Ω = 1/2 components are in full-lines and the Ω = 3/2 components are in dashed lines where Ω = |Λ + σ| (Λ and σ are the quantum numbers associated with the projection along z of the total electronic orbital and spin angular momenta respectively).

4

 4

R

  Figure 1:

4

 4

Figure 2

 2 Figure 2:

Figure 3

 3 Figure 3:

Figure 4 8 R

 48 Figure 4:

Figure 5

 5 Figure 5:

  1 D g . The density of states is smaller in that case with only 7 doublets and 2 quartets.

	**********************
	Table 1 near here
	**********************

  + is the 4 Σ -state. The different behaviour of

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

4 Σ -state is correlated to the { CO + [ 2 Σ + ]+C[ 3 P g ]} and the {C + 2 [ 2 Π u ]+O[ 3 P g ]} asymptotes and the 2 Π state to the {CO[ 1 Σ + ]+C + [ 2 P u ]} and the {C + 2 [ 4 Σ - g ]+O[ 3 P g ]} asymptotes. At this level of computation, the electronic ground state of C 2 O

Table 2 :

 2 Dominant configurations, energies and optimized linear geometries for the first electronic states of C 2 O + calculated by CASSCF method and cc-pVTZ basis set (cf. text).

		2.96	2.11	-150.320846
		2.984 a	2.132 a
	2 ∆ [...(1π 4 )(7σ 1 )(2π 2 )]	2.49	2.20	-150.240009
	2 Σ + [...(1π 4 )(7σ 1 )(2π 2 )]	2.47	2.23	-150.216771
	2 Σ -[...(1π 4 )(7σ 1 )(2π 2 )]	2.45	2.24	-150.211309
	4 Σ -[...(1π 4 )(7σ 1 )(2π 2 )]	2.57	2.13	-150.332491
		2.551 b	2.098 b
	4 Π [...(1π 3 )(7σ 2 )(2π 2 )]	2.64	2.47	-150.152082
	4 Σ -[...(1π 3 )(7σ 1 )(2π 3 )]	2.56	2.42	-150.113733

State / Configuration R CC (bohr) R CO (bohr) Energy (hartree) 2 Π [...(1π 4 )(7σ 2 )(2π 1 )]

Table 4 :

 4 Polynomial expansion coefficients (in u.a.) of the PEFs, equilibrium geometries and equilibrium dipole moment of the X 2 Π state and of the a 4 Σ -state. the reference is taken at the center of mass.

	Page 15 of 26
	954 2.148 -150.505063 2.580 2.166 -150.493383	2569
	cc-pVTZ 2.908 2.134 -150.621250 2.553 2.152 -150.610359	2393
	cc-pVQZ 2.897 2.127 -150.657942 2.546 2.145 -150.647108	2378
	cc-pV5Z 2.895 2.125 -150.669238 2.544 2.144 -150.658356	2388
	URL: http://mc.manuscriptcentral.com/tandf/tmph

a

Table 5 :

 5 Spectroscopic constants (in cm -1 ) of X 2 Π (A";A') and a 4 Σ -electronic states obtained from PEFs given by equation 1.

	X 2 Π

Table 6 :

 6 Decomposition of stretching normal modes Q 1 and Q 3 on displacements of internal coordinates ∆R CC = R CC -R eq CC and ∆R CO = R CO -R eq CO from the analysis of PEFs given in

Table 4 .

 4 State Normal mode ∆R CC ∆R CO

	X 2 Π	Q 1	0.8599 0.5105
		Q 3	-0.0277 0.9996
	a 4 Σ -	Q 1	0.8031 0.5959
		Q 3	-0.2310 0.9730

Table 7 :

 7 Polynomial expansion coefficients of the spin-orbit function A of the 2 Π state. The reference is the equilibrium geometry of the X 2 Π electronic state.

	Coefficient (10 -4 u.a.)
	C 000	-1.026284 C 100 0.568288 C 010 -0.478982
	C 200	-0.500797 C 110 0.818086 C 020 -0.264850
	C 002	0.213810 C 300 -0.021868 C 210 -0.141058
	C 120	0.147772 C 030 0.002741 C 102 -0.207771
	C 012	0.092577 C 400 0.240398 C 310 -0.653277
	C 220	0.594445 C 130 -0.364503 C 040 0.060324
	C 202	0.035374 C 112 0.075825 C 022 -0.117369
	C 004	-0.015375 C 500 -0.029922 C 600 -0.079419

Table 8 :

 8 Energies (in cm-1

Table 10 :

 10 First MRCI+Q / cc-pV5Z vertical transitions at the equilibrium geometry of the X 2 Π state for the doublet states and at the equilibrium geometry of the a 4 Σ -for the quartet states.

	3 Energy

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics F o r P e e r R e v i e w O n l y

Acknowledgments

The authors would like to thank P. Rosmus, G. Chambaud and F. Le Quéré for stimulating discussions.


Nicole Feautrier
email: nicole.feautrier@obspm.fr

Annie Spielfiedel

Christian Balança

Fabrice Dayou

Francois Lique

Maria Luisa Senent

A Spielfiedel †

C Balanc ¸a

M.-L Senent

Molecules as Diagnostic Tools in the Interstellar Medium

Keywords: Potential energy surface, inelastic collision, rotational excitation, interstellar medium Potential energy surface, inelastic collision, rotational excitation, interstellar medium

Introduction

Over the next few years, ground-based and space-based missions will open up the universe to high spatial and spectral resolution studies at infrared and submillimeter wavelengths. This will allow to study, in much greater details, the molecular composition of a variety of astrophysical regions, and, subsequently, to trace the origin and evolution of molecules in space. These new missions are expected to lead to the detection of many thousands of new spectral features. Identification, analysis and interpretation of these features, in terms of the physical and chemical characteristics of the astronomical sources, will require detailed astronomical tools supported by both laboratory measurements and theoretical calculations of spectroscopical data, reaction rate coefficients and collisional excitation rates for species of astrophysical relevance.

Whereas a proper modeling of interstellar molecular emission requires excitation calculations using collisional rate coefficients, to date, the determination of such data by means of the most accurate methods concern only few molecular systems of astrophysical interest. The collisional rate coefficients with the most abundant interstellar species He and H 2 as colliders are needed given that the steady-state populations in each ro-vibrational level of molecules are driven by the competition between collisional and radiative processes. Taking into account the sensitivity of the radiative transfer models to the collisional rate coefficients, and the influence of careful calculation of the excitation rate coefficients on the interpretation of observed spectra(1; 2; 3), it is thus of a great importance to provide accurate values of these rate coefficients for the molecules of astrophysical interest.

The present paper focuses in a first part on recent theoretical studies of collisional excitation processes at low temperatures for several molecular systems with He as a collider. In cold environments such as the interstellar medium (ISM), collisions with H 2 are the most important due to its high abundance.

Nonetheless, because of the near-sphericity of H 2 and its large rotational constant, collisions with He can be considered as a model for collisions with para-H 2 . In a second part of the paper, we make use of the newly obtained rate coefficients to analyse their impact on astrophysical modeling through radiative transfer models.

2 Molecular Collisions in the Interstellar Medium

Critical densities

The steady-state population in each quantum state of a molecular system depends on the competition between radiative and collisional processes. Their determination in the optically thin approximation relies on the statistical equilibrium equations(4), written for each level j as:

0 = l>j [n l A lj + (n l B lj -n j B jl )J] (1)
-

l<j [n j A jl + (n j B jl -n l B lj )J] + l n P [n l k lj -n j k jl ]
where the A and B terms are the Einstein coefficients for spontaneous and stimulated emission and absorption associated with the quantum states l and j, for the incident radiative field J. The k jl (k lj ) terms are the collision rate coefficients from j to l (l to j) and n P is the density of colliding particles (H 2 or He). The collision rate coefficients at a temperature T are related to the collisional cross sections through:

k ij (T ) = 8k B T πµ 1 2 1 k B T 2 × ∞ 0 E k σ ij (E k )exp(-E k /k B T )dE k (2)
The knowledge of the collision rate coefficients allows one to define for each level j a critical density n crit j such that(5):

n crit j = A j→j-1 i k ji (3)
where the sum runs over all reachable states. We may consider that the so-called Local Thermodynamical Equilibrium is fulfilled if the density n P is much larger than n crit j , since collisions tend to equilibrate the populations of the considered species in such a case. Conversely, the radiative processes are dominant for densities much lower than n crit j . Between these two extreme cases, both processes compete and accurate collision rates with the most abundant species are thus needed for the precise determination of the steadystate populations, and subsequently the physical conditions of the astronomical source. It is clear that studying different molecules and different ro-vibrational levels of each of them allows one to explore different regions of the studied astrophysical object. This point will be discussed in more details in the last section.

Potential energy surfaces

Within the Born Oppenheimer approximation, the determination of the collision rate coefficients requires two steps: the calculation of the interaction between nuclei and electrons for fixed nuclear geometries, leading to the so-called potential energy surface (PES), and the determination of the scattering cross sections associated with the molecular fragments for each quantum state of interest. In the context of observations in the mm and submm ranges, the collisional processes to be studied concern astrophysical regions characterized by relatively low temperatures (5K < T < 1500K). Moreover, for H 2 or He as the colliding particles, the molecular systems considered are generally of a van der Waals type. Under such conditions, electronic excitations and reactive channels are quite unlikely, and the collisional process can be treated by means of a single PES.

Recent progress in ab initio quantum chemistry facilitate very accurate determinations of the PES of small non-reactive systems. For such molecular systems, the coupled cluster theory including singles, doubles and perturbative triples excitations (CCSD(T)) is generally assumed to be very accurate, as long as only one configuration contributes significantly to the electronic state. It is generally the case for collisional excitation of rotational states in the interstellar medium, given that the temperature is very low and the relevant species are in their ground vibrational state. In that case, the geometries of the two partners are frozen at their equilibrium values, or their average values within the ground vibrational state. For collisions at higher energies, leading to possible vibrational excitation, this point must be checked out carefully as more than one configuration may contribute to the ground electronic PES when the interatomic distances of the studied molecule vary. In the case that a single electronic configuration is no more fulfilled, multiconfigurational approaches, such as the CASSCF-MRCI, may thus be needed to properly describe the PES.

In the following we discuss the results obtained for the CS, SiO and SO molecules with He as the collider. For each of these van der Waals systems, the determination of the PES relied on the use of the coupled-cluster method at the CCSD(T) level of theory. We employed the augmented correlation-consistent atomic basis sets of quadruple zeta quality (aug-cc-pVQZ) (7). To improve further the description of the wavefunction between the interacting fragments, we have added to this set the (3s3p2d2f1g) bond functions optimized by Cybulski and Toczylowski( 8), placed at the midpoint betwen the He atom and the studied diatom. The calculated interaction energies were corrected from the basis set superposition error (BSSE) at all geometries using the Boys and Bernardi(6) counterpoise procedure. All the calculations were performed with the MOLPRO suite of programs(9).

Dynamics

The most reliable approach to compute inelastic scattering cross sections is to perform quantum calculations by mean of the close-coupling method, first introduced by Arthurs & Dalgarno (10). The main limitation of such a theoretical approach comes from the large number of coupled channels that increases rapidly with energy, in relation with the number of available levels to be accounted for. While the exact close coupling approach, including all the angular momentum couplings, is feasible for rotational excitation of small molecules at low energies, it becomes unrealistic if the number of channels becomes too large, particularly if vibrational excitation has to be included. There exist various approximations (quantal or semi-classical) to the full close-coupling approach which were developed to circumvent such a limitation. Among them, the fast and popular Infinite Order Sudden (IOS) quantum approximation usually converges well whenever the relevant ∆j rotational energy spacings are small compared to the available collisional energy. When the energy is increased this condition is fulfilled for the lower j levels, but not for the highest j levels. In that respect, the extension of the present scattering calculations to the obtention of inelastic rates for warm gas or shocked regions may be problematic if accurate results are sought. All the results presented here were obtained within the full close-coupling approach using the MOLSCAT code [START_REF] Hutson | MOLSCAT computer code, version 14[END_REF]. the angular grid varied from 0 to 180 degrees by steps of 15 degrees. For the CS-He surface, the radial coordinate R was assigned values from 3.5 bohr to 16.0 bohr by steps of 0.25 bohr. The global minimum of the CS-He interaction energy was found to be -21.8051 cm -1 (R=7.26 bohr, θ = 75 degree). For the SiO-He surface, R ranged from 4.0 to 16.0 bohr by steps of 0.25 bohr and, at short atom-diatom separations, additional geometries were computed to insure that the interaction energy becomes at least 3500 cm -1 at any angle. The global minimum of the SiO-He interaction energy was found to be -26.596 cm -1 at the linear He-OSi geometry (R=7.993 bohr). The calculated interaction energies were fitted by means of the procedure described by Werner et al. (1988) in order to obtain the V (r = r e , R, θ) numerical expansion required for scattering calculations.

The new CS-He and SiO-He ab initio PESs were employed to perform full close-coupling calculations of the rotational (de)excitation cross sections for energies up to 3000 cm -1 (1; 12). The leading rotational de-excitation cross sections are decreasing functions of increasing kinetic energies, and display oscillations at low collision energies corresponding with resonances (of shape or Feshbach type) due to the presence of an attractive van der Waals potential well in the PES. Thermal average of the scattering cross sections yields the collision rate coefficients. We report in Fig. 1 the de-excitation rate coefficients of both CS and SiO molecules for a given initial rotational level j = 7 as a function of ∆j = jj ′ . It is seen that the de-excitation rate coefficients decrease with increasing ∆j, and propensity rules arise, namely: even ∆j rate coefficients larger than odd ∆j ones for CS and the reverse for SiO. Such behaviour is observed for all initial j level and all temperatures. Whereas the magnitude of the inelastic collision rate coefficients is driven by the anisotropy of the PES, almost similar for both molecular systems, the opposite propensity rules observed relate with the asymmetric character of each PES (14). In the case of CS-He collisions, the PES displays a shallow potential well of about 22 cm -1 close to a T-shape orientation of the fragments. Under such conditions, propensity rules almost similar to that of symmetric molecules are obtained (the excitation corresponding with odd ∆j transitions is forbidden in such a case). Conversely, for the SiO-He PES a shallow potential well of about 27 cm -1 is also obtained but at the linear He-OSi geometry, mainly due to the large differences between the atomic polarizabilities of the Si and O atoms (α Si ≫ α O ), that strongly favour odd ∆j transitions over even ∆j ones.

The rate coefficients for collisions with para-H 2 (j = 0) are expected to be larger than those associated with He as a collider, owing to the smaller reduced mass of the two collision partners in the case of para-H 2 . It is thus often applied a scaling factor to the collision rate coefficients with He to get a first estimate of the rates with para-H 2 , the scaling factor being close to 1.4 for molecules much heavier than H 2 . In such a case the underlying approximation is to consider identical cross section values for the two colliding systems. Another widely used approximation is to consider the interaction potential between a given molecule and He identical to the one associated with para-H 2 , and then perform scattering calculations using the reduced mass corresponding with para-H 2 to get the associated rate coefficients. The previous results for the CS-para H 2 and SiO-para H 2 colliding systems were estimated in such a way by Turner et al. (15), using CS-He and SiO-He PESs derived from the electron gas model. To assess the effect of the new ab initio PES on collision rate coefficients, we employed an identical approximation. For the two colliding systems, the comparison with the results of Turner (15) shows quite important disagreements, from a factor 2 up to 5 depending on the temperature and ∆j, that clearly point out the major role played by the PES.

2.4.2

SO. The previously available rate coefficients were provided by Green( 16) for the SO-H 2 system at temperatures ranging from 50K to 350K. The author used a CS-H 2 PES adapted from an electron gas model of CS-He. The rate coefficients among fine levels were obtained with a IOS scaling relationship (17). A new ab-initio PES was calculated using the CCSD(T) method and the basis sets described in paragraph 2.2. The SO r-distance was frozen at its experimental minimum energy distance (r e = 2.80 bohr). The values of the radial coordinate R ranged from 4 bohr to 16 bohr in steps of 0.25 bohr. The angular grid was uniform with a 15 degree spacing from 0 to 180 degree. At some angles, it was necessary to include smaller values of R to ensure that the calculations extended into the short-range repulsive region of the potential. The global minimum was found to be -34.921 cm -1 (R= 7.29 bohr, θ=180 degrees). The fitting procedure of Werner et al. (1988) was used.

From the new ab-initio PES, we have carried out full close coupling calculations with the exact energy levels including the fine structure interaction. Rate coefficients were obtained for temperatures ranging from 5K to 50K (18). The results displayed the expected propensity rule ∆j = ∆N .

The use of recoupling techniques from spin-free cross sections was investigated at low energy. This approximation is not valid for excitation between the first levels (N ≤ 5) where the fine structure splitting is large compared to the rotational splitting but it provides a reasonable estimate of rates for high-J levels at large temperatures. A comparison with the previous rate coefficients obtained by Green(16) for SO-H 2 , and transformed into SO-He coefficients by the ratio of masses, are on average a factor of 2-3 below the recent ones for the ∆j = ∆N transitions. As previously discussed for CS-He and SiO-He collisions, the precision of the PES appears to be also critical for the SO-He system. However, the differences may be as large as a factor of 10 for transitions involving the low-energy levels, in particular for the 0 1 -1 0 transition. This is due to the fine structure interaction (not explicitly included by Green) which leads to the 1 0 level energetically lower than the 0 1 level. We thus consider that the new collision rate coefficients are more reliable than the previous ones.

Astrophysical consequences

The population of the ro-vibrational levels of interstellar molecules is determined by both the density and by the radiation field (statistical equilibrium equations) whereas the radiation field is also determined by the population of the ro-vibrational levels of interstellar molecules. Indeed, we use radiative transfer codes to solve simultaneously the statistical equilibrium equations and the radiative transfer equation. These codes allow the determination of the physical conditions of molecular clouds (temperature, density, abundance of molecule ...) since, using them, we can simulate, for a set of physical conditions of a molecular cloud, the spectral line shapes of ro-vibrational transitions of observed molecules and compare them to the observations in order to validate the adopted model of the cloud.

The consequence on the determination of ISM physical conditions of using one set of collision rate coefficients or another is direct as the rate coefficients contribute through the product k ij n P to the statistical equilibrium equations (equation 1): increasing the rate coefficients k ij is equivalent to decreasing the density n P .

Taking as example the SO molecule, we have studied the astrophysical effect of the new rates by modeling the excitation of SO through a Large Velocity Gradient (LVG) radiative transfer code [START_REF] Sobolev | Theoretical Astrophysics[END_REF]. We have compared(2) the intensities I ν (expressed in terms of brightness temperature T B = (c 2 /2kν 2 )I ν ) obtained for observed 0 1 -1 0 , 2 2 -1 1 , 2 3 -1 2 , 6 5 -5 4 lines by using the rates of Lique et al. (18) and those of Green (16). The corresponding ratio is plotted in Fig. 2 as a function of the density for SO column densities ranging from 3.10 12 to 3.10 15 cm -2 by step factors of 10. )) for 0 1 -1 0 , 2 2 -1 1 , 2 3 -1 2 and 6 5 -5 4 lines of SO (from top to bottom). The H 2 volume density varies between 10 2 and 10 8 cm -3 and the SO column density from 3 10 12 and 3 10 15 cm -2 by a step factor of 10 (blue, green, red and black lines, respectively). From left to right the kinetic temperature corresponds to 50, 100, and 150 K, respectively.

We can see that the brightness temperature ratio is really different from one and that differences are important especially for the frequently observed 0 1 -1 0 line. In addition, the ratio is larger than one for the 0 1 -1 0 line whereas the situation is opposite for the other lines. The derivation of density or abundance of SO from observations will hence be dependent of the transitions used if the rate coefficients of Green (16) are used : if the 0 1 -1 0 line is used, the derived abundance of SO (or density) will be underestimated whereas an overestimation of the abundance of SO (or density) will be found if the other transitions are used.

We have applied these results to modeling of observations in TMC-1. We have used the observed line profiles of the transitions 0 1 -1 0 , 2 3 -1 2 , 3 4 -2 3 and 2 2 -1 1 in order to derive physical conditions and the SO abundance in this molecular cloud. We have used a non-local radiative transfer code and by comparing observations and synthetic spectra, we have built a model which allows a good fit of the observed spectra as shown in Fig. 3 (see (2) for details).

The structure of TMC-1 is well described by a cylindric filament containing several cores (T K = 8 K, n(H 2 ) ∼ 3 10 4 cm -3 ) along the ridge and an envelope (T K = 10 K, n(H 2 ) ∼ 6 -8 10 3 cm -3 ) surrounding these condensations.

Compared to previously published determinations of the SO abundance in this molecular cloud using Green's rate coefficients (16), the new results on SO appear to be more robust, providing better diagnostics to derive the physical conditions of cold molecular clouds. In particular, the use of previous rate coefficients could result in overestimating the abundance of SO by a factor 2 to 10 in cold molecular clouds.

Concluding remarks

From these few examples, it appears clearly that reliable astrophysical modeling needs accurate molecular data and in particular accurate collisional rate coefficients if the Local Thermodynamical Equilibrium is not fulfilled. This will be of a particularly great importance for the interpretation of the data that extremely powerful instruments such as Herschel and ALMA will provide in the next years.

Accurate predictions of rotational excitation cross sections can be obtained for a large set of small molecules in interaction with He and H 2 as long as only one electronic configuration contributes to the potential energy surfaces. The accuracy of the PES can reach a few cm -1 and close-coupling calculations can be performed. This fact is actually true for low temperatures but may be more questionnable for larger temperatures. In fact, the best level of accuracy is out of reach of current approximate approaches for scattering. Accordingly, the main source of uncertainties may come in the future from the scattering calculations instead of the PES, especially for highly excited rotational states or vibrationally excited molecules, for which fully converged close-coupling calculations become very difficult.

Acknowledgments

We are very pleased to thank Pavel Rosmus for his continuous and friendly help and interest in our work.

2

 2 Collisional excitation of molecules in the Interstellar Medium

  molecules in the Interstellar Medium 3

2. 4 Figure 1 .

 41 Figure 1. Rotational de-excitation rate coefficients for the CS(j = 7) (left panel) and SiO(j = 7) (right panel) molecules in collisionwith He, as a function of ∆j = jj ′ and at a given kinetic temperature.

  molecules in the Interstellar Medium 5

Figure 2 .

 2 Figure2. Brightness temperature ratios (T B (Green[13] rate coefficients) / T B (Lique et al.[15])) for 0 1 -1 0 , 2 2 -1 1 , 2 3 -1 2 and 6 5 -5 4 lines of SO (from top to bottom). The H 2 volume density varies between 10 2 and 10 8 cm -3 and the SO column density from 3 10 12 and 3 10 15 cm -2 by a step factor of 10 (blue, green, red and black lines, respectively). From left to right the kinetic temperature corresponds to 50, 100, and 150 K, respectively.

Figure 3 .

 3 Figure 3. Theoretical line profile (dashed line) compared to SO observations in TMC-1 (full line).

  ://mc.manuscriptcentral.com/tandf/


Martial Mazars
email: martial.mazars@th.u-psud.fr

L P T -Orsay

Yukawa potentials in systems with partial periodic boundary conditions I : Ewald sums for quasi-two dimensional systems

Keywords: Yukawa potentials, Ewald sums, Monte Carlo, electroneutrality, colloids

Yukawa potentials in systems with partial periodic boundary conditions I : Ewald sums for quasi-two dimensional systems.

Introduction

Yukawa interaction potentials, named likewise after the achievement of the meson theory by Hideki Yukawa in 1935 [START_REF] Yukawa | Meson theory in its developments[END_REF], were introduced by Debye and Hückel in 1923 as a mean field approximation in the Poisson-Boltzmann equation for ionic solution [START_REF] Debye | [END_REF] ; these potentials, in their static form, are given by

V (r) = Q exp(-κ | r |) | r | ( 1
)
where κ is the inverse screening length, and Q an effective Yukawa charge. These quantities are related to some physical parameters of systems ; some examples are given in Table 1.

The Yukawa potential is solution of the Helmholtz equation

(∆ -κ 2 )V = -4πQδ(r) (2)
with the boundary condition V (∞) = 0.

In theories of ionic systems as liquids [3,[START_REF] Toda | Statistical Physics I : Equilibrium Statistical Mechanics[END_REF], colloids [START_REF] Verwey | Theory of the Stability of Lyophobic Colloids[END_REF][START_REF] Hunter | Foundations of Colloid Science[END_REF][START_REF] Russel | Colloidal Dispersions[END_REF][START_REF] Belloni | [END_REF], plasmas [9,10], etc., κ -1 is known as the Debye length, this length gives an estimation of the radius of the screening sphere and Q is an effective charge that may be related to the surface potential ϕ S of particles with hard-core interaction. More precisely, for hard spheres systems we have

Q = ϕ S a exp(κa)
with a the radius of the hard spheres (cf. Table 1).

In Yukawa's meson theory and according to the general principle of the quantum theory, κ is related to the mass m π of the meson by m π = κh/c and Q is an effective coupling constant for nuclear interactions [START_REF] Yukawa | Meson theory in its developments[END_REF].

On general grounds, Yukawa interaction potentials can be considered as a reasonable approximation of effective interaction potentials between particles when some microscopic degrees of freedom may be approximated to a continuous background that screens the direct interaction between particles, while the spherical symmetry of the interaction is preserved [11]. Despite the apparent simplicity of the analytical form of the Yukawa potentials in Eq.( 1), the physical mechanisms by which such effective potentials may be derived from the microscopic degree of freedom can be quite complicated and these mechanisms depend strongly on the properties of studied system. For instance, one may compare the difference between the charging mechanisms of dust grains in dusty plasma systems [9,10,15,16,17] and the equilibrium Poisson-Boltzmann approximation or the other approximations schemes in colloidal systems [START_REF] Verwey | Theory of the Stability of Lyophobic Colloids[END_REF][START_REF] Belloni | [END_REF]11,12,13,14].

In computer simulations, one uses periodic boundary conditions to avoid irrelevant surface bias.

As the lattice sums of Yukawa potential are absolutely convergent for any value of the inverse

Debye length κ, one may safely use a truncation of the potential and minimum image convention, as long as κ is large enough. But, in the Coulomb limit (κ → 0), truncation of the potential introduces important bias and errors [18,19,20] ; therefore, there is a domain with κ = 0 where one has to handle lattice sums of Yukawa potential with caution by using convenient algorithm as Ewald methods [21,22]. As outlined before by authors of refs. [21,22], truncation of Yukawa potentials for bulk systems may safely be used as long as exp(-κL)/L is small. Typically, for systems with number density ρ * 1 with N 1000 Yukawa particles, truncation of the Yukawa potential can be used for κ * = κσ ≥ 2 (cf. Table 1) [23,24,25,26]. It is also worthwhile to note that a code with Ewald sums for Yukawa potential, for any value of κ, can easily be obtained, with minimal modifications, from a code with Ewald sums for Coulomb potential [21,22].

Many interesting experimental systems in plasmas physics [27,28] and in colloids science are quasi-two dimensional systems [30,31], or even quasi-one dimensional systems [32]. To study such systems with numerical simulations, one have to use partial periodic boundary conditions in one or two directions. The purpose of the present work is to provide a detailed derivation of the Ewald method for Yukawa potentials in quasi-two dimensional systems. It seems that Ewald method for Yukawa potentials in quasi two dimensional systems have already been used in some works of refs. [28], as outlined in ref. [29] ; however, the analytical details of the Ewald sums for Yukawa potentials are not gven in these works. In a paper numbered II, the Lekner sums [33,35] for Yukawa potentials will also be given and in the paper III, we give numerical results for a bilayer system of particles interacting via a Yukawa potential ; some preliminary results for this Yukawa bilayer system are given in the section 4 of the present paper. The analytical form of Ewald sums for Yukawa potential obtained in the present work shows that, as for bulk systems, a code can be easily written from one code for Coulomb potential with minimal modifications. In this paper, a special attention is paid to the Coulomb limit and its relation with electroneutrality. From this analysis we obtain the constants useful to compute the total energy. These constants depend on the geometry of systems and on the way the electroneutrality is achieved.

The present work is organised as follow. In section 2, we derive the Ewald sums of Yukawa interaction potential for quasi-two dimensional systems by using a method similar to the one used by Grzybowski and co-workers for the Coulomb interaction potential [36]. At the end of section 2, the particle-particle interaction energy for a system of N particles in a quasi-two dimensional system, this result does well agree with a recent result [37] obtained from both the Ewald sums for Yukawa potential in three dimensional systems [22] and the Parry method [38]. In section 3, the small screening limit and its relations with electroneutrality of systems are outlined. In particular in this section, we show explicitly how singular contributions are cancelled by electroneutrality for three different systems ; these results allow to link the total energy of systems interacting with Yukawa potential to the energy of systems interacting with Coulomb potential.

Section 4 is devoted to a discussion on the practical use of the Ewald method for Yukawa potential and to numerical tests. For completeness two appendices have been added at the end of the paper. In appendix A, we provide a full derivation of the Ewald sums for Yukawa potentials in three dimensional systems with the method by Grzybowski and co-workers [36] ; the results of this appendix agree with the derivation done by Salin and Caillol who have computed the Ewald sums from the Green's function of the Helmholtz equation [22]. In appendix B, we compute the forces obtained with the Ewald method for molecular dynamics implementations.

2 Ewald sums for Yukawa interaction in quasi-two dimensional systems.

We consider a system made of N particles interacting via Yukawa potentials. The simulation box has partial periodic boundary conditions, with spatial periodicity L x and L y , applied respectively to directions x and y, whereas no periodic boundary conditions are taken in the third direction z, parallel to the unitary vector êz . In the simulation box, the position of the particle i is defined

r i = s i + z i êz (3)
The particle-particle interaction energy is given by

E cc (Yukawa; κ) = 1 2 N i=1 N i =j Q i Q j Φ(r ij ) + 1 2 N i=1 Q 2 i Φ 0 (4)
with Φ(r) defined by

Φ(r) = n exp(-κ | r + n |) | r + n | and Φ 0 = n =0 exp(-κ | n |) | n | (5)
In previous equations, we have used the condensed notations

| r + n |= (x + n x L x ) 2 + (y + n y L y ) 2 + z 2 and | n |= n 2 x L 2 x + n 2 y L 2 y
with n x and n y integer numbers associated with periodic images of the box. In appendix A, one would find the computation for periodic images in all the three directions of the simulation box.

The computation starts with the identity

exp(-κ | r + n |) | r + n | = 1 √ π ∞ 0 dt √ t exp(- κ 2 4t -| r + n | 2 t) (6)
that will allow to use later the Poisson-Jacobi identity.

With Eq.( 6), Φ(r) may be written as

Φ(r) = 1 √ π ∞ 0 dt √ t exp(- κ 2 4t ) n exp(-| r + n | 2 t) = 1 √ π n ∞ α 2 dt √ t exp(- κ 2 4t ) exp(-| r + n | 2 t) + 1 √ π α 2 0 dt √ t exp(- κ 2 4t -z 2 t) n exp(-| s + n | 2 t) = I 1 (r; α; κ) + I 2 (r; α; κ) (7)
where the integral has been split into a short ranged contribution I 1 (r; α; κ) and a long ranged contribution I 2 (r; α; κ). In Eq.( 7), α is the Ewald damping parameter and, for practical applications, its value is chosen to balance efficiency and accuracy ; a comparison between the results of appendix A and the results in ref. [22] shows that one may interpret the parameter α in Eq. [START_REF] Russel | Colloidal Dispersions[END_REF] as the gaussian width of a charge distribution that screens the central charge in the short ranged contribution [Note : In ref. [22], the parameter α is denoted by β, and κ by α]. As remarkably outlined in ref. [22], for three dimensional systems, the optimum choice for α is independent of κ and is the same as for Ewald sums for Coulomb interaction, i.e. as κ → 0.

With the integral relation

∞ α 2 dt √ t exp(- κ 2 4t -a 2 t) = ∞ 2α dt exp(- a 2 t 2 4 -κ 2 t 2 ) = √ π 2a [ exp(κa)erfc(αa + κ 2α ) + exp(-κa)erfc(αa - κ 2α )] (8)
the short ranged contribution to the potential is given by

I 1 (r; α; κ) = 1 2 n D(| r + n |; α; κ) | r + n | (9)
where

D(r; α; κ) = exp(κr)erfc(αr + κ 2α ) + exp(-κr)erfc(αr - κ 2α ) (10)
On Figure 1, we show the function D(r; α; κ) versus κr for several values of the ratio λ = κ/α ; on this figure, the function C(r; α; κ), used to compute the short ranged contribution of forces and defined by Eq.(B.4) in appendix B, is also represented.

The long ranged contribution is computed by using the Poisson-Jacobi identity

n exp(-| s + n | 2 t) = 1 A π t G e iG.s exp(- G 2 4t ) (11)
where A = L x L y is the surface of the periodic simulation cell ; thus I 2 (r; α; κ) becomes

I 2 (r; α; κ) = √ π A α 2 0 dx x 3/2 exp(- κ 2 4x -z 2 x) + √ π A G =0 e iG.s α 2 0 dx x 3/2 exp(- G 2 + κ 2 4x -z 2 x) (12)
With the help of Eq.( 8) and after some algebra, we obtain

α 2 0 dx x 3/2 exp(- k 2 4x -z 2 x) = √ π k [ exp(kz)erfc( k 2α + αz) + exp(-kz)erfc( k 2α -αz)] = √ π F(k; z; α) (13
I 2 (r; α; κ) = π A [F(κ; z; α) + G =0 e iG.s F( G 2 + κ 2 ; z; α)] (14)
Therefore, Φ(r) is given by

Φ(r) = 1 2 n D(| r + n |; α; κ) | r + n | + π A [F(κ; z; α) + G =0 e iG.s F( G 2 + κ 2 ; z; α)] (15)
with D(| r + n |; α; κ) and F(k; z; α) defined respectively by Eqs. (10) and (13).

The lattice sums in Eq.( 5) are absolutely convergent, as a consequence, as long as κ = 0, there is no diverging contribution in the analytical decomposition of Φ(r) into short and long ranged contributions. As κ → 0, a 1/κ-asymptotic behaviour is found ; this stems from the conditionally convergent property of the Coulomb lattice sums. This diverging behaviour of the Coulomb lattice sums is handled by taking into account the boundary conditions as | r |→ ∞ and the electroneutrality of the simulation cell [36]. The analysis of the asymptotic behaviour of Yukawa-Ewald sums as κ → 0 and their links with the Coulomb-Ewald sums will be discussed extensively in the next section.

To compute the energy E cc of the interaction between the particles, the self contribution Φ 0 is needed. This contribution is computed as for Φ(r). We have

Φ 0 = n =0 exp(-κ | n |) | n | = 1 √ π n =0 ∞ α 2 dt √ t exp(- κ 2 4t ) exp(-| n | 2 t) + 1 √ π α 2 0 dt √ t exp(- κ 2 4t -z 2 t) n exp(-| n | 2 t) - 1 √ π α 2 0 dt √ t exp(- κ 2 4t -z 2 t) (16)
where the last integral has been added and substracted to the long ranged contribution to permit the use of the Poisson-Jacobi identity. Then

Φ 0 = I (0) 1 + I (0) 2 (17)
and after simple algebra we find 4), ( 15) and ( 18), the particle-particle interaction energy is given by

                   I (0) 1 = 1 2 n =0 D(| n |; α; κ) | n | I (0) 2 = 2π A G erfc( √ G 2 + κ 2 /2α) √ G 2 + κ 2 -2[ α √ π exp(- κ 2 4α 2 ) - κ 2 erfc(κ/2α)] (18
E cc (Yukawa; κ) = 1 4 N i,j Q i Q j n D(| r ij + n |; α; κ) | r ij + n | + π 2A N i,j Q i Q j G =0 e iG.s ij F( G 2 + κ 2 ; z ij ; α) + π 2A N i,j Q i Q j F(κ; z ij ; α) -( N i=1 Q 2 i )( α √ π exp(- κ 2 4α 2 ) - κ 2 erfc(κ/2α)) (19)
where the prime in the summation over the image indicates that the term i = j has to be omitted for n = 0.

The total particle-particle interaction energy is given by E cc (Yukawa; κ) and may be used in Monte-Carlo simulations ; Eq.( 19) agrees with the results obtained by using another method in ref. [37]. For Molecular Dynamics simulations, the computation of forces is needed, this derivation is done in appendix B from the results of section 2 and 3.

For repulsive Yukawa interaction between particles, one needs to compute the particle-background and the background-background energies to obtain the total energy of the system and thus to be able to implement correctly a Monte-Carlo sampling of the phase space. These computations are done in the next section for three different neutralizing backgrounds, and their limits as κ → 0 to a non-screened Coulomb system are outlined.

3 The Coulomb limit and electroneutrality.

For Coulomb interactions in systems with periodic boundary conditions, the electroneutrality property of the system is important to allow to define the electrostatic potential. Electrostatic Yukawa interactions are equivalent to non-screened Coulomb interactions as the screening length tends to infinity (or the inverse screening length κ → 0) ; therefore, in the Coulomb limit, the Yukawa interaction energy of the system, computed with Eqs.(4-5), has to be equal to the Coulomb interaction energy. The other limit, κ → ∞, is either a hard sphere system (a, σ, σ ij = 0) or an ideal gas limit (a, σ, σ ij = 0).

The series in Eq.( 5) are absolutely convergent when κ > 0 (Yukawa) and conditionally convergent when κ = 0 (Coulomb). For Coulomb interactions, the conditional convergence property of the series leads to a diverging behaviour and this diverging behaviour is cancelled by using the electroneutrality of the system ; thus some constants, that depend on the geometry of the simulation box and on the manner in which electroneutrality is achieved, arise. As long as κ = 0, the series are absolutely convergent and no diverging behaviour is found in the total interaction energy. To achieve consistency between Yukawa and Coulomb potentials in the Coulomb limit, one has to recover from Eqs.( 19) and ( 5) a diverging behaviour corresponding

to Coulomb potential and one has to be able to cancel it by using the electroneutrality of the system. This section is devoted to the examination of the Ewald sums consistency between Yukawa and Coulomb potentials for quasi-two dimensional systems.

In subsection 3.1, the diverging terms appearing in the particle-particle interactions (Eq.19) as κ → 0 are displayed and in subsections 3. 3.1 The Particle-Particle interaction energy in the Coulomb limit.

According to Table 1, if temperatures are high and/or if counterions are low, the inverse Debye length may become quite small. In the special case κ = 0, one should recover the properties of One Component Plasma systems, thus, as κ → 0, Coulomb interaction energy has to be recovered from Eq.( 19).

As κ → 0, we have we have also

               D(| r |; α; κ) = 2 erfc(α | r |) + [erfc(α | r |) -exp(-α 2 | r | 2 ) √ πα | r | ] | r | 2 κ 2 + o(κ 4 ) κ F(κ; z; α) = 2 -2[z erf(αz) + 1 α √ π e -α 2 z 2 ]κ + κ 2 + o(κ 3 ) (20
G F( √ G 2 + κ 2 ; z; α) = [ exp(Gz)erfc( G 2α + αz) + exp(-Gz)erfc( G 2α -αz)] + κ 2 2G [(1 -Gz)e Gz erfc( G 2α + αz) -(1 + Gz)e -Gz erfc( G 2α -αz) - 2G α √ π exp(- G 2 4α 2 -α 2 z 2 )] +o(κ 4 ) (21)
and

α √ π exp(- κ 2 4α 2 ) - κ 2 erfc( κ 2α ) = α √ π - κ 2 + 3κ 2 4α √ π + o(κ 3 ) (22)
Therefore, we have

E cc (Yukawa; κ → 0) = E cc (Coulomb) + π A ( i Q i ) 2 1 κ + ( i Q 2 i ) κ 2 + o(κ 2 ) ( 23
)
with

E cc (Coulomb) = 1 2 N i,j Q i Q j n erfc(α | r ij + n |) | r ij + n | + π 2A N i,j Q i Q j G =0 e iG.s ij F(G; z ij ; α) - π A N i,j Q i Q j [z ij erf(αz ij ) + 1 α √ π e -α 2 z 2 ij ] -α √ π ( i Q 2 i ) (24)
the Coulomb particle-particle interaction energy for a system with quasi-two dimensional geometry [36]. In Eq.( 23), the singular contribution as 1/κ appears explicitly. This singular term is cancelled if the system of particles has the electroneutrality property (i.e. i Q i = 0), then the small screening limit (κ → 0) of the energy of the system, E cc (Yukawa; κ) in Eq.( 19), is exactly 24), the energy of the system of particles in interaction via a Coulomb potential. The electroneutrality implies that some pairs of particles in the system have attractive Yukawa interaction as, for instance, in the restricted primitive model of electrolytes (RPM)

E cc (Coulomb) in Eq.(
with Coulomb poyentials.

However, for applications to plasmas or colloids, one rather uses the Yukawa One Component directions of the space, the electroneutrality is obtained by embedding the particles in a uniform neutralizing background [22] (cf. Appendix A). For quasi-two dimensional systems, such uniform volume neutralizing background can not be used since boundary conditions in the zdirection depend on the model studied. In the next three subsections, the small screening limits of three different neutralizing backgrounds, adapted to the quasi-two dimensional geometry, are examined. These backgrounds are monolayer, bilayer and slab neutralizing backgrounds, they are schematically represented on Figure 2.

The monolayer neutralizing background

The monolayer neutralizing background represented on Figure 2 (a) is a plan with an uniform surface charge density localised at z = 0, particles may be localised on both sides of the plan or their location may be restricted also to only one half-space. Hereafter, system (a) designs the system made by the monolayer and particles ; for this system, the charge density in the right hand side of the Helmholtz equation, Eq.( 2), is given by

ρ a (r) = i Q i δ(r -r i ) + σδ(z) (25)
where δ(z) is the Dirac distribution. Then, the electroneutrality reads as

N Q + Aσ = 0 ( 26
)
The energy of the system (a) is given by

E (a) (Yukawa; κ) = E cc (Yukawa; κ) + E cB (κ) + E BB (κ) (27)
with E cB (κ) is the interaction energy of the particles with the background, while E BB (κ) is the self energy of the background.

                 E cB (κ) = - N Q 2 A N i=1 S ds +∞ -∞ dzδ(z) n exp(-κ | r i -r + n |) | r i -r + n | E BB (κ) = N 2 Q 2 2A 2 S ds +∞ -∞ dz δ(z ) S ds +∞ -∞ dzδ(z) n exp(-κ | r -r + n |) | r -r + n | (28
After elementary algebra, we find

                 E cB (κ) = - N Q 2 A √ π N i=1 ∞ 0 dt t 3/2 exp(-z 2 i t - κ 2 4t ) E BB (κ) = N 2 Q 2 2A √ π ∞ 0 dt t 3/2 exp(- κ 2 4t ) (29)
and thus we have

E (a) (Yukawa; κ) = E cc (Yukawa; κ) -2π N Q 2 A 1 κ N i=1 exp(-κ | z i |) + π N 2 Q 2 A 1 κ (30)
In the Coulomb limit, all singular contributions as 1/κ are cancelled by taking into account the electroneutrality (cf.Eq.( 26)) ; in this limit we have

E (a) (Yukawa; κ → 0) = E cc (Coulomb) + 2π N Q 2 A N i=1 | z i | +o(κ) (31)
The second contribution in Eq.( 31) stems from E cB (κ → 0). The energy of the Coulomb system is obtained by taking κ = 0. The energy E (a) (Yukawa; κ) computed with the Ewald sums is therefore defined for all value of κ and the limit as κ → 0 is well defined : it corresponds to a quasi-two dimensional coulombic system made of N particles carrying charge Q (One Component Plasma), and a monolayer with a surface charge density σ to achieve electroneutrality, Eqs. (25) and (26).

The bilayer neutralizing background

The bilayer neutralizing background is represented on Figure 2 (b), this background is made of two plans, h is the separation between both plans and σ the surface charge density in each plan.

Particles may be localised between the plans or outside the slab. We design by system (b) the system made by the bilayer and the particles. For this system, the charge density appearing in the Helmholtz equation is

ρ b (r) = i Q i δ(r -r i ) + σδ(z - h 2 ) + σδ(z + h 2 ) ( 32
)
12

N Q + 2Aσ = 0 (33)
For system (b), the energy is given by

E (b) (Yukawa; κ) = E cc (Yukawa; κ) + E cB (κ) + E BB (κ) (34)
with the same definitions as for system (a). With simple computations, one finds

                                   E cB (κ) = - N Q 2 A √ π( N i=1 ∞ 0 dt t 3/2 exp(-(z i - h 2 ) 2 t - κ 2 4t ) + N i=1 ∞ 0 dt t 3/2 exp(-(z i + h 2 ) 2 t - κ 2 4t
))

E BB (κ) = N 2 Q 2 8A √ π(2 ∞ 0 dt t 3/2 exp(- κ 2 4t ) + ∞ 0 dt t 3/2 exp(-h 2 t - κ 2 4t )) (35)
then

E (b) (Yukawa; κ) = E cc (Yukawa; κ) -π N Q 2 A 1 κ N i=1 ( exp(-κ | z i - h 2 |) + exp(-κ | z i + h 2 |)) +π N 2 Q 2 2A 1 κ (1 + exp(-κh)) (36)
Similarly to system (a), the singular contributions in system (b) are cancelled by the electroneutrality. As h → 0, we recover

E (b) (Yukawa; κ) → E (a) (Yukawa; κ) (37)
and as κ → 0 and h = 0, one has

E (b) (Yukawa; κ → 0) = E cc (Coulomb)+π N Q 2 A N i=1 (| z i + h 2 | + | z i - h 2 |)-π N 2 Q 2 A h+o(κ) (38)
The same analysis may be done easily for multilayer neutralizing backgrounds.

If N 0 = N/2 particles are confined in each plan, then we have z i = -h/2 or z i = h/2 and the energy of the system is then E (b) (Yukawa; κ → 0) is exactly the energy of the bilayer Wigner system studied in refs. [35,39]. This will provide an useful reference point.

E (b) (Yukawa; κ → 0) = E cc (Coulomb) + 2πL 2 σ 2 h + o(κ) (39

The slab neutralizing background

The slab neutralizing background is represented on Figure 2 (c), the slab between the plans located at z = -h/2 and z = h/2 is filled with an uniform volume charge density ρ 0 . Outside the slab, the volume charge density is null. Particles may be localised inside or outside the slab.

We design by system (c) the system made by the slab and the particles. For this system, the charge density appearing in the Helmholtz equation is

ρ c (r) = i Q i δ(r -r i ) + ρ 0 (θ(z + h 2 ) -θ(z - h 2 )) (40)
with θ(z) the Heaviside distribution. For systems (c), electroneutrality reads as

N Q + Ahρ 0 = 0 (41)
Similarly to systems (a) and (b), the total energy of system (c) is given by

E (c) (Yukawa; κ) = E cc (Yukawa; κ) + E cB (κ) + E BB (κ) (42)
with

                 E cB (κ) = - N Q 2 Ah √ π N i=1 h 2 -h 2 dz ∞ 0 dt t 3/2 exp(-(z i -z) 2 t - κ 2 4t ) E BB (κ) = N 2 Q 2 A 2 h 2 4π κ 2 S ds h 2 -h 2 dz (1 -e -κh/2 cosh(κz )) (43)
To compute E cB (κ), one has to consider the set of particles that are confined in the slab and the set of particles being outside the slab, therefore we define

I = { particles i :| z i |> h 2 } and J = { particles j :| z j |< h 2 }
and we note

               E cB (κ) = - N Q 2 Ah √ 4π κ 2 [ sinh( κh 2 ) I e -κ|z i | + N J -e -κh/2 J cosh(κz j )] E BB (κ) = N 2 Q 2 Ah 2π κ 2 [1 -2 κh e -κh/2 sinh( κh 2 )] (44)
Again, as κ → 0, the singular contribution is cancelled by the electroneutrality and we have

E (c) (Yukawa; κ → 0) = E cc (Coulomb) + π N 2 Q 2 h 2L 2 [ N J N - 1 3 + 4 N h I | z i | + 4 N h 2 J z 2 j ] + o(κ) (45)
As h → 0 and with the constraint ρ 0 h = σ, we find easily

E (c) (Yukawa; κ) → E (a) (Yukawa; κ) (46)

Discussion

A comparison between Eq.( 19), that gives the total particle-particle interaction energy for Yukawa potential, and Eq.( 24 When all constants arising from the backgrounds contribution and self energies are included, then we may obtain a code for Yukawa potential for any value of the inverse screening length κ, and this code will be consistent with the Coulomb limit. One should note also that the complexity of the long ranged contribution is the same as for Coulomb potential : the sum over the pairs and reciprocal lattice vectors G may not be simplified into summations over particles (with exceptions for some particular systems as the bilayers studied in refs. [39,35]).

As already noted by authors of refs. [21,22,24], the use of Ewald sums for Yukawa potentials is only necessary when κ is small, otherwise a brute truncation of the potential and the minimum image convention are sufficient to perform accurate computations. To be more specific, we may estimate a lower bound κ cut for κ that allows to obtain a relative accuracy of the order δ in the computation of the energy. The criterion we use to estimate κ cut is as follow.

Consider a system of Yukawa particles with hard sphere diameter a, then the interaction energy at contact is E c = Q 2 exp(-κa)/a, while, if the separation between particles is L/2, with L the smallest spatial periodicity of the simulation box, the energy is

E L = 2Q 2 exp(-κL/2)/L.
Then, a direct truncation of the potential will allow to obtain a relative accuracy of the order δ if

E L /E c < δ.
Then, we obtain an estimation of the lower bound κ cut that allows safe truncations of Yukawa potentials for distance greater than L/2, as

κ * = κa > κ * cut = 2 L * -2 ln( 2
δL * ) = σ * N √ N -σ * N ln( 4σ * N δ 2 N ) ( 47
)
where we use the reduced units L * = L/a, κ * = κa, σ * N = a 2 σ N and σ N is the surface coverage of the simulation box defined by σ N = N/L 2 , N being the number of Yukawa particles.

On table 2, we give some value of κ * cut for different values of the surface coverage and the number of particles, for relative accuracy δ = 10 -4 and 10 -6 .

In Table 3, we report some numerical tests of the Ewald sums for Yukawa potentials. For these tests, we consider a pair of Yukawa particles with effective charge Q 1 = +1 and Q 2 = -1 in a box with L x = L y = 1 (in these reduced units, we have κL ≡ κ ≡ κ * L * ) ; for this system, one has Q 1 + Q 2 = 0, thus we can use straightforwardly Eq.( 19) to compute the energy of the system and as κ → 0, Eq.( 23) is well behaved. The five configurations reported in Table 3 have already been considered before to test Lekner sums for Coulomb potentials [34,33,[START_REF] Tildesley | The Molecular Dynamics Method[END_REF][START_REF] Liem | [END_REF] ; the values reported in the Column 'Coulomb Limit' have been computed with analytical, Lekner [33] and/or the Hautman-Klein methods [42,[START_REF] Tildesley | The Molecular Dynamics Method[END_REF]. For κ = 10 -6 , the results obtained with Eq.( 19) perfectly agree with the Coulomb limit for Ewald damping parameters α = 6 and α = 12 and with k × k = 16 × 16 for the number of vectors G in the reciprocal space. The results obtained with a direct truncation of the Yukawa potential do not allow to reproduce the correct Coulomb limit and energies computed with a direct truncation for κ = 1 are very different than the value obtained with Ewald sums. On Figure 3, we represent, for three configurations of the pair, the particle-particle energy E cc (Yukawa; κ) as a function of κ (symbols) and the energy computed with a direct truncation of the Yukawa potential (lines) ; the Coulomb limit for each configu- 2 (κ ≡ κ * L * ). One may also note on Fig. 3 that, for the configuration (0.1,0.1,0.1) of the pair, by chance the energy computed with Ewald sums agrees quite well with the one computation using a direct truncation of the potential ; however an inspection of the numerical values shows that the difference is greater than 10 -4 .

On Figure 4, we represent ∆E cc (κ, α) = E cc (κ, α) -E cc (κ, α = 6.) as a function of the Ewald damping parameter α for the three configurations of Figure 3. The contributions of the reciprocal space are computed with k × k = 16 × 16 for the number of vectors G. These data show the stability of the Ewald sums for Yukawa potentials versus the arbitrary choice of the damping parameter α ; more precisely, an inspection of the numerical data shows that, if 4. ≤ α < 15., then | ∆E cc (κ, α) |≤ 10 -6 . Therefore, a choice αL 5. -6. is convenient to use Ewald sums for Yukawa potentials, as already noted for Coulomb potentials [22]. On the one hand, if α is too small, then the minimum image convention is insufficient and additional images have to be included in the real space contribution ; on the other hand, if α is too large additional G-vectors have to be included in the reciprocal space contributions (cf. Fig. 4).

On Table 4 and Figure 5, we report some preliminary results from Monte Carlo computations on Yukawa Bilayer systems. This system is the same as the one studied in refs. [35,39] but the particles interact with a repulsive Yukawa potential and the electroneutrality is achieved with the background of subsection 3.3. The geometrical and physical parameters of the bilayer correspond to run a of ref. [35] ; more precisely, we have h = 1, N = 512, L = 28.36 and the coupling constant Γ = Q 2 /kT a = 196. Intralayer and interlayer energies and correlation functions are defined as in ref. [35]. On Table 4, data for Coulomb interaction are extracted from ref. [35], the results for Yukawa interaction have been obtained with the Ewald sums given by Eqs.( 19) and (26). For κ = 10 -4 , energies and correlation functions, on Fig. 5 (a,b), agree perfectly ; this shows that the Coulomb limit is correctly obtained with the results of section 2 and 3. For larger values of κ, the results obtained with the Ewald sums are compared to results obtained with a direct truncation of the Yukawa potential. For κ = 0.1, energies and correlation functions obtained with a direct truncation differ significantly from the results obtained with Ewald sums. For κ = 0.25, energies obtained with a direct truncation of the Yukawa potential agree with energies computed with Ewald sums, but, as it appears on Fig. 5 (c,d), correlation functions disagree. For larger values of κ, both energies and correlation functions computed with Ewald sums or with a direct truncation of the interaction potential are in agreement with a good accuracy. More details and properties on Yukawa bilayer systems will be given in the forthcomming paper numbered III.

There are some deviations from Yukawa potentials [3,9,43] ; this is particularly true for quasitwo dimensional systems since counterions and coions density profiles in the z-directions are non-uniform. However, if, as in the theoretical analysis done in ref. [14], one defines an effective screening length that depends only on the z-coordinates of the pair of particles, then one may use the Yukawa-Ewald sums, Eq.( 19), in replacing the uniform screening parameter κ by a non-uniform κ(z, z ) (cf. Eqs. (25,26) of ref. [14]). Otherwise, if the screening of the Coulomb interaction is more complicated, one should come back to the reliability of the Yukawa interaction potential for the system considered.

As a final comment to this paper, we would like to outline the fact that Coulomb limits may be obtained accurately for Yukawa potentials with the use of Ewald sums ; but, as it is shown in Section 3, the effective Coulomb limit of a system is depending on the manner the electroneutrality is restored from the background. To compute the particle-particle interaction energy of the system with periodic boundary conditions in three dimensions, we follow the same method as in section 2 for quasi-twodimensional systems. Potentials are still defined by

Φ(r) = n exp(-κ | r + n |) | r + n | and Φ 0 = n =0 exp(-κ | n |) | n | (A.1)
but summations over the periodic images are taken along the three dimensions of the space. With periodic boundary conditions in all three directions of the space, the Poisson-Jacobi identity reads as

n e -|r+n| 2 t = 1 V ( π t ) 3/2 G e iG.r exp ( - G 2 4 1 t ) (A.2)
From Eqs.(A.1) and ( 6), we have

Φ(r) = 1 √ π n ∞ α 2 dt √ t exp(- κ 2 4t ) exp(-| r + n | 2 t) + 1 √ π α 2 0 dt √ t exp(- κ 2 4t -z 2 t) n exp(-| r + n | 2 t) = I (3D) 1 (r; α; κ) + I (3D) 2
(r; α; κ)

(A.3)
and with Eq.( 8), we find easily

I (3D) 1 (r; α; κ) = 1 2 n D(| r + n |; α; κ) | r + n | (A.4)
with D(| r + n |; α; κ) given by Eq. (10).

I (3D) 2
(r; α; κ) is computed by using the Poisson-Jacobi Identity (A.2), we have then Without any difficulties, Φ 0 is also computed as follows

I (3D) 2 (r; α; κ) = 4π V e -κ 2 /4α 2 κ 2 + 4π V G =0 e iG.r exp ( -(G 2 + κ 2 )/4α 2 ) (G 2 + κ 2 ) (A.5)
Φ 0 = 1 2 n =0 D(| n |; α; κ) | n | + 4π V G =0 exp ( -(G 2 + κ 2 )/4α 2 ) (G 2 + κ 2 ) + 4π V e -κ 2 /4α 2 κ 2 -2[ α √ π exp(- κ 2 4α 2 ) - κ 2 erfc( κ 2α )] (A.6)
Therefore the particle-particle interaction energy is given by

E (3D) cc (Yukawa; κ) = 1 4 N i,j Q i Q j n D(| r ij + n |; α; κ) | r ij + n | + 2π V N i,j Q i Q j G =0 e iG.r ij e -(G 2 +κ 2 )/4α 2 (G 2 + κ 2 ) + 2π V ( N i=1 Q i ) 2 e -κ 2 /4α 2 κ 2 -( N i=1 Q 2 i )[ α √ π exp(- κ 2 4α 2 ) - κ 2 erfc( κ 2α )] (A.7)
The singular term as κ → 0 in Eq.(A.7) is as 1/κ 2 , in the Coulomb limit we have

E (3D) cc (Yukawa; κ → 0) = E (3D) cc (Coulomb) + 2π V ( i Q i ) 2 1 κ 2 + ( i Q 2 i ) κ 2 + o(κ 2 ) (A.8) with E (3D) cc (Coulomb)
given by

E (3D) cc (Coulomb) = 1 2 N i,j Q i Q j n erfc(α | r ij + n) | r ij + n | + 2π V N i,j Q i Q j G =0 e iG.r ij e -G 2 /4α 2 G 2 - α √ π ( i Q 2 i ) (A.9)
The singular term in Eq.(A.8) has to be cancelled by the electroneutrality of the system. Assuming that for all particles we have Q i = Q > 0 and that the electroneutrality is achieved by an uniform background ρ 0 , then the charge density in the simulation box is

ρ (3D) (r) = i Q i δ(r -r i ) + ρ 0 (A.10)
and the electroneutrality reads as Thus, for the three dimensional system, the total energy is

N Q + V ρ 0 = 0 (A.11) 20
E (3D) (Yukawa; κ) = E (3D) cc (Yukawa; κ) + E (3D) cB (κ) + E (3D) BB (κ) (A.12) with                  E cB (κ) = - N Q 2 V N i=1 V dr n exp(-κ | r i -r + n |) | r i -r + n | E BB (κ) = N Q 2 2V V dr V dr n exp(-κ | r -r + n |) | r -r + n | (A.13)
After a simple computation of gaussian integrals, we obtain

             E cB (κ) = -4π N Q 2 V 1 κ 2 E BB (κ) = 2π N Q 2 V 1 κ 2 (A.14)
and therefore, for a system with an uniform neutralizing background, we have For molecular dynamics implementations, the computation of forces is needed ; in the following we give the forces obtained for Yukawa potentials in quasi-two dimensional systems.

E (3D) (Yukawa; κ → 0) = E cc (Coulomb) + o(κ) (A.
The force acting on particle i is given by

F i = F (cc) i + F (cS) i = -∇ r i E (cc) i -∇ r i E (cS) i = j F (cc; short) j/i + j F (cc; long) j/i + F (cS) S/i (B.1)
where

E (cc) i and E (cS) i
are respectively interaction energies of particle i with other particles (cf.

Eq.( 19)) and with the neutralizing background (cf. Eqs. (30,36) or ( 43)). In Eq.(B.1), the force acting on particle i due to other particles is split into a short and a long ranged contribution.

With the chain rule

∂f (| r ij + n |) ∂x i = ∂ | r ij + n | ∂x i × df (r) dr | r=|r ij +n| (B.2)
and Eq.( 19), the short range contribution of the force acting on particle i due to particle j is given by

F (cc; short) j/i = -∇ r i E (cc; short) j/i = 1 4 Q j Q j n C(| r ij + n |; α; κ) | r ij + n | 3 [(x ij + n x L x )ê x + (y ij + n y L y )ê y + z ij êz ] (B.3) with x ij = x i -x j and F (r; α; κ) defined by C(r; α; κ) = (1 -κr) exp(κr)erfc(αr + κ 2α ) + (1 + κr) exp(-κr)erfc(αr - κ 2α ) (B.4)
Similarly, the long range contribution to particle-particle forces is given by

F (cc; long) j/i = -∇ r i E (cc; long) j/i = -∇ s i E (cc; long) j/i - ∂ ∂z i E (cc; long) j/i êz (B.5)
From Eq.( 19), we find

-∇ s i E (cc; long) j/i = - π 2A Q i Q j G =0 sin(G.s ij ) F( G 2 + κ 2 ; z ij ; α) G (B.
j/i = - π 2A Q i Q j G e iG.s ij B( G 2 + κ 2 ; z ij ; α) (B.7) with B(k; z; α) = exp(kz)erfc( k 2α + αz) -exp(-kz)erfc( k 2α -αz) (B.8)
The force induced by the background on particle i depends on the geometrical parameter of the background. For systems (a), (b) and (c) described in section 3, there are some discontinuities in spatial distributions of the neutralizing background, hence forces are also discontinuous at those locations. Such behaviours may induce some complicated bias in molecular dynamics computations, especially when particles cross the surfaces of discontinuity. To implement molecular dynamics codes for systems (a), (b) and (c), one would have to estimate numerically the influence of these discontinuities on the trajectories of particles.

Forces induced by backgrounds on particle i may be easily computed with

F (cS) (a,b,c)/i = - ∂ ∂z i E (cS)
(a,b,c)/i êz (B.9)

For system (a), in the half-space z i > 0, we have

F (cS) (a)/i = - 2πN Q 2 A exp(-κz i ) êz (B.10)
for system (b), in the slab -h/2 < z i < h/2,

F (cS) (b)/i = 2πN Q 2 A e -κh/2 sinh(κz i ) êz (B.11)
and for system (c), in the slab -h/2 < z i < h/2, List of Tables Table 1 : Definitions of the Yukawa parameters Q and κ in Eqs.(1,2) as functions of physical parameters and thermodynamical states of systems with screened-Coulomb interactions. For all systems, k B is the Boltzmann constant, T the temperature, e the elementary charge, 0 the vacuum dielectric constant and the dielectric permittivity of the medium. In the Debye-Hückel approximation of electrolytes, q i and ρ i are respectively the charge and the number density of ions of species i, and for the Primitive Model of electrolytes, σ ij = (σ i + σ j )/2, is the mean ionic diameter of ions i and j [3,[START_REF] Toda | Statistical Physics I : Equilibrium Statistical Mechanics[END_REF]. For colloidal systems, Z is the charge number and σ the hard-core diameter of macroions, ρ is the number density of the monovalent counterions and λ B = e 2 / k B T is the Bjerrum length [START_REF] Verwey | Theory of the Stability of Lyophobic Colloids[END_REF][START_REF] Hunter | Foundations of Colloid Science[END_REF][START_REF] Russel | Colloidal Dispersions[END_REF][START_REF] Belloni | [END_REF]. In Plasmas and Dusty Plasmas, the parameters for screening-Coulomb potentials are Z and a the charge number and the radius of the impurity (dust particle), n 0 is the plasma density, T i and T e are the ion and electron temperatures [9,10].

F (cS) (c)/i = √ 4π N Q 2 A e -κh
Table 2 : Values of κ * cut for several values of the surface coverage and number of particles in the simulation box. If κ * = κa ≥ κ * cut , then one may use a direct truncation of Yukawa potential for distances greater than L/2 ; otherwise the long range of the potential has to be handled with a convenient algorithm.

Table 3 : Numerical tests of the Coulomb limit of the Ewald sums for the Yukawa potentials. The configuration of the pair of particles is defined by (x 12 , y 12 , z 12 ) particle 1 with a charge Q 1 = +1 located at (0,0,0) and particle 2 with a charge Q 2 = -1 located at (x 12 , y 12 , z 12 ) in a box with L x = L y = 1. The values of particle-particle Coulomb energies are extracted from original works by others, in the original works the computations of energies have been done analytically or with the Lekner method or with the Hautmann-Klein method (HK) [42]. Evaluations of particle-particle interaction energies in columns Table 4 : Average energies computed with the Monte Carlo algorithm for a Yukawa bilayer system. The geometrical and physical parameters of the bilayer correspond to run a of ref. [35] ; the coupling constant is Γ = 196, N = 512 is the number of particles and the spatial periodicity L = 28.36. βU/N is the average total energy per particle, βE intra /N and βE inter /N are respectively the average intralayer and interlayer energies per particle, βσ U /N is an estimation of the statistical fluctuation σ U of the total energy. Data for the Coulomb interaction (κ = 0) are extracted from ref. [35] The particle 1 carries a charge Q 1 = +1 and is located at (0,0,0) ; the particle 2 carries a charge Q 2 = -1 and three positions of the particle 2 are considered in a box with periodic boundary conditions with L x = L y = 1.. The three positions of particle 2 are : a : (0.1, 0.1, 0.1) ; b : (0., 0., 0.25) and c : (0.5, 0.5, 0.). Symbols refer to evaluations done with Ewald sums while lines represent evaluations done with a direct truncation at L/2. of the Yukawa potentials. For each configuration, the value of the Coulomb limit is indicated by a thick horizontal line and limiting values are explicitly given. These configurations have already been considered in some previous works [33,34,[START_REF] Tildesley | The Molecular Dynamics Method[END_REF][START_REF] Liem | [END_REF] (see also Table 3). 4. The correlation functions are defined in ref. [35]. For the Coulomb interaction, the correlation functions have been obtained from the ref. [35] where the computations have been done with Ewald, Hautmann-Klein and Lekner methods.

  2 to 3.4, we provide the results obtained for several backgrounds necessary to fulfill electroneutrality : a monolayer in subsection 3.2 (Fig 2.a), a bilayer (subsection 3.3, Fig 2.b) and a slab (subsection 3.4, Fig 2.c). For any other neutralizing backgrounds, one may also compute the particle-background and the background-background energies according to similar methods.

) 9 Page

 9

10 Page

 10 Plasma model (YOCP) in which all particles have effective charge of the same sign and interactions between particles are purely repulsive ; therefore the electroneutrality can not be fulfilled by summation over the particles. For YOCP with periodic boundary conditions in the three

) 11 Page

 11

) 13 Page

 13

CardI = N I and CardJ = N J 14 Page

 14 we have N I + N J = N . With these definitions of sets I and J, we find

  ) giving the energy for Coulomb potential, shows that a code for Yukawa potential in quasi-two dimensional systems may easily be built from existing computer codes for Coulomb potential. Apart from constants, the main modifications to perform are : to change the function erfc(α | r ij + n |) in the short range contribution to the Coulomb energy by the function D(| r ij + n |; α; κ)/2 ; in the long ranged contribution, for G = 0, to change G in the function F(G; z ij ; α) by √ G 2 + κ 2 and to modify the contribution for G = 0.

16 Page

 16 by a thick red line. It clearly appears on figure3, that if κ < 10 the use of Ewald sums is necessary ; this agrees with the results of Table

17 Page

 17

  Appendix A: Ewald sums for Yukawa potential in three dimensional systems.

15 ) 21 Page

 1521

28

 28

  Ewald , refer to evaluations performed by using straightforwardly the Ewald sums for Yukawa potentials given by Eqs.(19), and with k × k = 16 × 16 for the number of vectors G in the reciprocal space and α as indicated in the Table.The term Truncation refers to evaluations of the energies of pairs with a direct truncation at L/2. of the Yukawa potentials.

Figure 1 :

 1 Figure 1: Representation of D(r; α; κ) and C(r; α; κ) as functions of κr for several values of the ratio λ = κ/α. Function D and C are respectively defined by equations (10) and (B.4), the function D is used to compute the short ranged part of the energy while the function C is used to compute the short ranged part of the force between particles.As κr → 0, we have D(r; α; κ) 2 and C(r; α; κ) 2, this corresponds to a non screened Coulomb interaction between particles.

Figure 2 :

 2 Figure 2: Schematic representation of the neutralizing backgrounds of three systems : (a) monolayer with surface charge density σ = -N Q/A, (b) bilayer with surface charge density σ = -N Q/2A in each layer, and (c) slab with a volume charge density ρ 0 = -N Q/Ah.

Figure 3 :

 3 Figure 3: Representation of E cc (Yukawa; κ) for three configurations of a pair of Yukawa particles as functions of κ. E cc (Yukawa; κ) are computed by using equation (19), with α = 6.0 and k × k = 16 × 16 for the number of vectors G in the reciprocal space.

34 PageFigure 4 :

 344 Figure 4: Dependence of the particle-particle interaction energy on the Ewald parameter α. The configuration of the pair of particles considered are the same as those on Figure 3, with Q 1 = +1 and Q 2 = -1. The quantity ∆E cc (κ, α) = E cc (κ, α) -E cc (κ, α = 6.) is computed with k × k = 16 × 16 for the number of vectors G in the reciprocal space and (a) κ = 10 -6 ; (b) κ = 1. Inspection of the numerical values shows that, if 4. ≤ α < 15., then | ∆E cc (κ, α) |≤ 10 -6 .

Figure 5 :

 5 Figure 5: Representation of the intralayer g 11 (s) and interlayer g 12 (s) correlation functions obtained for the Yukawa bilayer system with the same physical parameters as in Table4. The correlation functions are defined in ref.[35]. For the Coulomb interaction,

35 PageFigure 2 :

 352 Figure 1: M. Mazars, Yukawa I : Ewald sums.

37 Page

 37 Figure 3: M. Mazars, Yukawa I : Ewald sums.

Figure 4 : 39 Page

 439 Figure 4: M. Mazars, Yukawa I : Ewald sums.

Figure 5 : 40 Page

 540 Figure 5: M. Mazars, Yukawa I : Ewald sums.

) 6 Page 6 of 45 URL: http://mc.manuscriptcentral.com/tandf/tmph

	Molecular Physics
	and thus
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

) 7 Page 7 of 45 URL: http://mc.manuscriptcentral.com/tandf/tmph

	Molecular Physics
	With Eqs.(
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

19 Page 19 of 45 URL

	Molecular Physics
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

: http://mc.manuscriptcentral.com/tandf/tmph

Page 29 of 45 URL: http://mc.manuscriptcentral.com/tandf/tmph

  and data for the Direct Truncation have been computed with a direct truncation of the Yukawa potential for distance greater than L/2.

		Molecular Physics
	Systems and/or		κ		Q	physical
	approximations	(inverse screening length) (effective charge) parameters
	Debye-Hückel F (Primitive Model o r of electrolytes)	κ 2 =	4π k B T i	ρ i q 2 i	q i q i exp(κσ ij ) (1 + κσ ij )	q i , ρ i and σ ij
	P
	e e κ 2 = 4πλ B ρ r (monovalent counterions) Colloids	Ze exp(κσ/2) (1 + κσ/2)	Z, λ B , σ and ρ
	Plasmas and Dusty Plasmas (a = 0) (a = 0)	R e v i e κ 2 = 4πn 0 e 2 0 k B ( 1 T e + 1 ) T i w O Ze 0 n Ze exp(κa) 0 (1 + κa)	Z, n 0 , T i , T e and a
					l y
			29

Table 1 M

 1 . Mazars, Yukawa I : Ewald sums.

	30

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 2

 2 M. Mazars, Yukawa I : Ewald sums.

	31

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 3

 3 M. Mazars, Yukawa I : Ewald sums.

	32

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 4

 4 M. Mazars, Yukawa I : Ewald sums.

	33

URL: http://mc.manuscriptcentral.com/tandf/tmph

URL: http://mc.manuscriptcentral.com/tandf/tmph

ACKNOWLEDGEMENTS

I acknowledge the computation facilities of the Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) under the project 72104. The computations have been done on IBM Regatta Power 4. I am very grateful to D. Levesque, J.-J. Weis and J.-M. Caillol for useful and interesting discussions ; and to V. Huet for her help in the preparation of the manuscript.