Keywords: Conditions aux limites absorbantes, coordonées sphéroïdales prolates, opérateur Dirichlet-to-Neumann, problèmes de scattering

Construction et analyse de performance de conditions aux limites absorbantes de type DtN localisé pour des problèmes extérieurs d'Helmholtz. Partie II : Frontières sphéroïdales prolates Résumé : Nous proposons une nouvelle classe de conditions aux limites absorbantes de type DtN localisé applicables à des frontières extérieures de forme ellipsoïdale lors de la résolution de problèmes de scattering d'ondes acoustiques par des obstacles allongés. Ces conditions sont : (a) exactes pour les premiers modes, (b) faciles à implémenter et à paralléliser, (c) compatibles avec la structure locale d'un schéma numérique de type éléments finis, et (d) applicables à des frontières extérieures de forme ellipsoïdale ce qui est plus approprié en terme de coût de calcul lorsque l'on entoure un obstacle allongé. De plus, ces conditions coïncident avec la condition DtN localisée classique construite pour des frontières de forme sphérique. Nous étudions analytiquement et au travers de tests numériques l'effet du régime de fréquence et de l'allongement de la frontière sur la performance de ces conditions lors de la résolution des problèmes de radiateur et de scattering. Nous comparons aussi leurs performances à celles de la condition aux limites absorbante d'ordre deux (BGT2) construite par Bayliss, Gunzburger et Turkel lorsqu'elle est exprimée en coordonnées sphéroïdales prolates. Cette analyse révèle que, en régime basse fréquence, la condition d'ordre deux DtN (DtN2) admet un très bon niveau de performance quel que soit l'allongement de la frontière. De plus, la condition aux limites DtN2 est plus performante que la condition BGT2. Cette supériorité apparaît clairement pour de grandes excentricités.

Introduction

Exterior Helmholtz problems are classical mathematical models for studying scattering problems arising in many applications such as sonar, radar, geophysical exploration, nondestructive testing, etc... Despite their simplicity, this class of problems is not completely solved particularly from a numerical point of view. For example, the computation of the solutions of these problems requires first to limit it to a finite domain. This is often achieved by surrounding the given scatterer(s) (or radiator) by an artificial boundary that is located at some distance (measured in multiples of wavelength of interest) from its surface. A so-called "nonreflecting" boundary condition is then prescribed on the artificial boundary to represent the "far-field" behavior of the scattered field. The challenge here is the development of a simple but reliable as well as cost-effective computational procedure for representing the far-field behavior of the scattered. The quest for such conditions is ongoing (see, e.g., the recent review by Turkel in the book [START_REF] Turkel | Iterative methods for the exterior Helmholtz equation including absorbing boundary, conditions[END_REF]).

We propose in this work new three-dimensional approximate local DtN boundary conditions to be employed on prolate spheroid boundaries that are primary candidates for surrounding elongated scatterers. The idea for constructing such conditions is driven by several considerations chief among them the following two reasons. First, the widely-used second order absorbing boundary condition (BGT2) designed by Bayliss, Gunzburger and Turkel for spherical-shaped boundaries [START_REF] Bayliss | Boundary conditions for the numerical solution of elliptic equations in exterior regions[END_REF] performs poorly when it is expressed in elliptical coordinates and applied to prolate spheroid boundaries in the low frequency regime [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]. The accuracy deteriorates significantly for large eccentricity values of the boundaries as observed in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]. The damping effect introduced to this condition [START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF] improves the performance for small eccentricity values. However, the modified BGT2 still performs poorly for eccentricity values larger than 0.6 in the (relatively) low frequency regime (see Figure 15 in [START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF]). Hence, there is a need for constructing local absorbing boundary conditions (ABC) that extend the range of satisfactory performance. Second, the three-dimensional approximate local DtN conditions designed for spherical-shaped boundaries [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF] perform very well for low wavenumber values as reported in [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. We recall that, in R 3 , the second-order DtN boundary condition and the BGT2 condition are identical [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. Nevertheless, using these conditions on spherical-shaped exterior boundaries when solving scattering problems by elongated scatterer often leads to larger than needed computational domains, which hampers computational efficiency. This suggests that approximate local DtN boundary conditions designed for prolate spheroidal-shaped boundaries is an attractive alternative for improving the computational performance.

Given that, this work is devoted to the construction of theses conditions and to the assessment of their performance when employed on prolate spheroidal-shaped boundaries for solving three-dimensional radiator and scattering problems. The idea of constructing three-dimensional approximate local DtN boundary conditions is not new. Indeed, as stated earlier, such conditions have been already derived for spherical-shaped boundaries [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF]. The construction procedure adopted in [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF] is based on the localization of the truncated global DtN boundary condition [START_REF] Keller | Exact nonreflecting boundary conditions[END_REF]. The key ingredient of this procedure is the trigonometric identities that express high order derivatives of sine and cosine functions. However, this property is not satisfied by the angular spheroidal wave functions [4]. Consequently, the procedure used in [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF] is no longer applicable to the truncated global DtN boundary operator when expressed in elliptical coordinates [START_REF] Grote | On nonreflecting boundary conditions[END_REF][START_REF] Givoli | Exact representations on artificial interfaces and applications in mechanics[END_REF]. Hence, the construction methodology we propose for deriving the class of approximate local DtN boundary conditions in prolate spheroidal coordinates can be viewed as an inverse-type approach. More specifically, we start from a Robin-type boundary condition with unknown coefficients. Unlike the case of polar coordinates, these coefficients depend on the angles (ϕ, θ) of the prolate spheroidal coordinates. Such dependence is necessary to preserve the symmetry and local nature of the resulting boundary conditions. Then, we require that the considered condition to be an exact representation of the first modes. Consequently, the coefficients are the unique solution of a linear algebraic system.

We assess mathematically and numerically the performance of the constructed approximate local DtN boundary conditions. More specifically, we analyze the effect of low wavenumber and the eccentricity on the performance of these conditions in the case of three-dimensional radiator and scattering problems. We adopt the on-surface radiation condition formulation (OSRC) [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering using an on-surface radiation boundarycondition approach[END_REF] in order to perform analytically this investigation. We note that such formulation is not appropriate for high frequency regime as observed previously in [START_REF] Antoine | Fast approximate computation of a time-harmonic scattered field using the onsurface radiation condition method[END_REF]. The main interest in the following analyses is to evaluate the performance of the proposed approximate local DtN conditions at low wavenumber to see if relatively small computational domains can be employed in order to avoid excessive computational cost. The OSRC formulation must be viewed as an extreme case while an exterior elliptical-shaped boundary surrounding an elongated scatterer would be less "demanding" on the boundary condition. The analysis herein shows that, in the case of the radiator, we can observe the DtN2 condition performs similarly to the BGT2 condition for values of e smaller than 0.6 and it outperforms BGT2 when e is bigger than 0.6. In the case of the scattering problem, the DtN2 condition outperforms the BGT2 condition even for small values of e, the DtN2 solution being clearly better for e ≥ 0.6.

The reminder of this paper is as follows. We specify in Section 2 the nomenclature and assumptions adopted in this paper. Then, we present the class of new local DtN boundary conditions in prolate spheroidal coordinates in Section 3, along with their properties. We show that these conditions are exact for the first modes. In section 4, we perform the mathematical and numerical analysis to assess computational performance of these conditions when applied for solving radiator and scattering problems. We conclude this work in Section 5.

Preliminaries

Throughout this paper, we use the prolate spheroidal coordinates (ξ, ϕ, θ) which are lrelated to the cartesian coordinates (x, y, z) by the following transformation:

x = b sin ϕ cos θ, y = b sin ϕ sin θ, z = a cos ϕ (1)
INRIA where ϕ ∈ [0, π), θ ∈ [0, 2π). The parameters a and b respectively represent the major and the minor axes and are given by:

a = f cosh ξ, b = f sinh ξ (2)
where ξ is a strictly positive real number and f is the interfocal distance and is defined by:

f = a 2 -b 2 (3)
We also define the eccentricity e on a prolate spheroid at ξ = ξ 0 , by:

e = 1 cosh ξ 0 = 1 - b 2 a 2 (4)
The eccentricity e characterizes the slenderness of the surface. It satisfies 0 < e < 1. Note that when e → 0, the prolate spheroid degenerates into a sphere and when e → 1, the spheroid degenerates into a line with length 2f on the z-axis.

We recall that the mn-th prolate spheroidal mode u mn (m ≤ n) is given by [START_REF] Senior | Scalar diffraction by a prolate spheroid at low frequencies[END_REF]:

u mn = R (3) mn (kf, cosh ξ)S mn (kf, cos ϕ) cos mθ (5)
The functions R

mn (kf, cosh ξ) are the radial spheroidal wave functions of the third kind and S mn (kf, cos ϕ) denote the angular spheroidal wave functions which are solution to the differential equation [4]:

∂ ∂ϕ sin ϕ ∂S mn ∂ϕ + sin ϕ λ mn -(kf ) 2 cos 2 ϕ - m 2 sin 2 ϕ S mn (kf, cos ϕ) = 0 (6)
where λ mn := λ mn (kf ) is called the prolate spheroidal eigenvalue (see p.11 in Ref. [4]).

3 The new approximate local boundary conditions and their derivation

Next, we construct the new local DtN absorbing boundary conditions in prolate spheroidal coordinates, and describe their properties.

The approximate local DtN boundary conditions in Prolate spheroidal coordinates

The three-dimensional first-(DtN1) and second-order (DtN2) local Dirichlet-to-Neumann boundary conditions, defined on the prolate spheroid-shaped surface ξ = ξ 0 , are given by:

DtN1 : ∂ u ∂ξ = √ 1 -e 2 e R 00 u (7)
DtN2 :

∂ u ∂ξ = √ 1 -e 2 (λ 01 -λ 00 ) e λ 01 R 01 -λ 00 R 00 -(R 00 -R 01 ) (eka) 2 cos 2 ϕ u + (R 00 -R 01 ) ∆ Γ u (8)
where the coefficients R mn depend on the radial spheroidal wave functions of the third kind R

(3) mn , the wavenumber ka, and the eccentricity e as follows:

R mn = ∂R (3) mn ∂ξ (eka, e -1 ) R (3) mn (eka, e -1 ) (9)
and ∆ Γ denotes the Laplace Beltrami operator, which reads in prolate spheroidal coordinates (ξ, ϕ, θ) as:

∆ Γ = 1 sin ϕ ∂ ∂ϕ sin ϕ ∂ ∂ϕ + 1 sin 2 ϕ ∂ 2 ∂θ 2 (10)
The following four remarks are noteworthy:

First, the Robin-type boundary condition (7) (resp. ( 8)) is, by construction, an exact representation of the first (resp. the first and second) prolate spheroidal mode(s) (see Eq.5).

Second, the boundary conditions ( 7), ( 8) are called local DtN conditions because they result from a localization process of the truncated global DtN boundary operator defined in [START_REF] Grote | On nonreflecting boundary conditions[END_REF][START_REF] Givoli | Exact representations on artificial interfaces and applications in mechanics[END_REF]. The local feature of these conditions is of a great interest from a numerical view point. Indeed, the incorporation of these conditions in any finite element code introduces only mass-and stiffness-type matrices defined on the exterior boundary. The coefficients λ mn and R mn can be computed once for all at the preprocessing level.

Third, it must be point out that when e = 0 (the prolate spheroid becomes a sphere), conditions ( 7) and ( 8) are identical to the three-dimensional DtN conditions designed for spherical shaped boundaries [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF][START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. This property can be easily established using the asymptotic behavior of the radial spheroidal wave functions of the third kind R

(3) mn and the prolate spheroidal eigenvalues λ mn .

Last, approximate local DtN conditions of order higher than two are inappropriate for conventional finite element implementations since they require regularity higher than C 0 . Consequently, we consider only the first (DtN1) and second order (DtN2) conditions in the following.

INRIA

The procedure for constructing the approximate local DtN conditions

As stated earlier in the introduction, the idea of constructing approximate local DtN boundary conditions is not new. Indeed, such conditions have been already derived for sphericalshaped boundaries in [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF]. The construction procedure adopted in [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF] is based on the localization of the truncated global DtN boundary condition [START_REF] Givoli | Exact representations on artificial interfaces and applications in mechanics[END_REF]. The key ingredient of this procedure is the trigonometric identities that express high order derivatives of sine and cosine functions (see for example Eqs. (C15-C21) p. 215 in [START_REF] Harari | Computational methods for Problems of Acoustics with Particular Reference to Exterior Domains[END_REF] for the 3D case). However, this property is not satisfied by the angular spheroidal wave functions (see Eq. ( 6)), and therefore, the procedure used in [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF] cannot be applied to the truncated global DtN boundary operator when expressed in elliptical coordinates [START_REF] Grote | On nonreflecting boundary conditions[END_REF][START_REF] Givoli | Exact representations on artificial interfaces and applications in mechanics[END_REF]. The approach we propose for constructing these approximate local DtN boundary condition in the case of prolate spheroid-shaped boundaries can be viewed as an inverse-type methodology. More specifically, we start from a Robin-type boundary condition with unknown coefficients. These coefficients depend on the prolate spheroidal eigenvalue.In addition, we require that the considered condition is an exact representation of the first modes. Consequently, the coefficients are the unique solution of a linear algebraic system. Next, we describe this procedure for construction the approximate local DtN1 and DtN2 conditions given by ( 7) and ( 8).

i. Construction of the first order DtN1 boundary condition

We consider a Robin-type boundary condition that is an exact representation of the first spheroidal mode u 00 given by Eq. ( 5). Hence, we set:

∂u ∂ξ = Au (11)
where the constant A satisfies (sinh ξ) ∂R

(3) mn ∂ξ (kf, cosh ξ)S 00 (kf, cos ϕ) = A R

(3) 00 (kf, cosh ξ)S 00 (kf, cos ϕ) (12)

Then, we substitute, at ξ = ξ 0 , Eqs, (2)-(4) into Eq. ( 12), and deduce that

A = √ 1 -e 2 e R 00 (13)
where R 00 is given by Eq. ( 9).

ii. Construction of the second order DtN2 boundary condition

The goal here is to construct a Robin-type boundary condition that is an exact representation of the first two modes u 00 and u 01 given by Eq. ( 5). Hence, we set::

∂u ∂ξ = C u + D (∆ Γ -(eka) 2 cos 2 ϕ)u (14)
where C and D are constant (independent of ϕ) to be determined. Note that, unlike DtN2 boundary condition for the spherical-shaped boundaries, the coefficients of this condition depend on the angular variable ϕ. Such dependence is necessary for constructing a symmetric boundary condition since the angular spheroidal wave functions satisfied differential equation given by Eq. ( 6). Indeed, it is easy to verify that, for all modes u mn , we have:

∆ Γ u mn = -λ mn + (eka) 2 cos 2 ϕ u mn (15)
In order to determine the constants C and D, we assume that, at ξ = ξ 0 , we have:

∂u mn ∂ξ = C u mn + D ∆ Γ -(eka) 2 cos 2 ϕ u mn ; m = 0 and n = 0, 1 (16)
Then, using Eq. ( 15), it follows that (C, D) is the unique solution of the folllowing 2 × 2 linear system:

         C -D λ 00 = √ 1 -e 2 e R 00 C -D λ 01 = √ 1 -e 2 e R 01 (17)
where the coefficients R mn are given by Eq. ( 9). The DtN2 boundary condition given by Eq. ( 8) is a direct consequence of solving the system [START_REF] Senior | Scalar diffraction by a prolate spheroid at low frequencies[END_REF] and substituting the expressions of (C, D) into Eq. ( 14).

The performance of the new approximate local DtN conditions

In the following, we assess analytically and numerically the performance of the approximate local DtN1 and 2 boundary conditions given by ( 7) and [START_REF] Grote | On nonreflecting boundary conditions[END_REF]. More specifically, we analyze the effect of low wavenumber ka and the eccentricity e on the performance of DtN1 and DtN2 in the case of both radiator and scattering problems. We adopt the on-surface radiation condition formulation (OSRC) [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering using an on-surface radiation boundarycondition approach[END_REF] in order to perform this investigation analytically.

We note that such formulation is not appropriate for high frequency regime as observed previously in [START_REF] Antoine | Fast approximate computation of a time-harmonic scattered field using the onsurface radiation condition method[END_REF]. The main interest in the following analyses is to evaluate the performance of the proposed approximate local DtN conditions at low ka to see if relatively small computational domains can be employed in order to avoid excessive computational cost. The INRIA OSRC formulation is an extreme case while an exterior elliptical-shaped artificial boundary surrounding an elongated scatterer would be less "demanding" on the boundary condition.

As in [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF][START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF][START_REF] Reiner | Improvement of the performance of the BGT2 condition for low frequency acoustic scattering problems[END_REF], we assess the performance of the ABCs DtN1 and DtN2 using the specific impedance introduced in [START_REF] Geers | Doubly asymptotic approximations for transient motions of submerged structures[END_REF][START_REF] Geers | Third-order doubly asymptotic approximations for computational acoustics[END_REF] as a practical tool for measuring the efficiency of ABCs in the context of the OSRC formulation. This non-dimensional quantity measures the effect of the truncated medium in physical terms. It provides a convenient indicator of the performance of a given approximate representation. In the elliptical coordinates system, the specific impedance can be expressed as follows:

Z = i √ 1 -e 2 ka u ∂u ∂ξ | ξ=ξ0 (18)
4.1 The performance of the local DtN conditions for single prolate spheroidal mode

Next, we assess analytically and numerically the performance of the approximate local DtN1 and DtN2 boundary conditions given by ( 7) and ( 8) when applied to three-dimensional radiator problems. Recall (see Eq. (64), p. 3645 in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]), that the exact specific impedance Z ex3 mn , for the mn-th prolate spheroidal mode, on the surface of a prolate spheroid at ξ = ξ 0 is given by:

Z ex3 mn = i √ 1 -e 2 kaR (3)
mn (eka, e -1 ) ∂R

(3) mn ∂ξ (eka, e -1 ) = ieka R mn (19)
where the coefficient R mn is given by Eq. ( 9). In addition, Its asymptotic behavior as ka → 0 (see Eq. (67), p. 3646 n [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]) is given by:

Z ex3 mn ∼ (ka) n+1 (n + 1)(2n -1)!! 2 -i ka n + 1 ; as ka → 0 ( 20
)
where the double factorial is defined in [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF][START_REF] Abramovitz | Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables[END_REF].

Analytical study

The following lemma states the expression of the approximate specific impedances for the mn-th prolate spheroidal mode, on the surface of a prolate spheroid at ξ = ξ 0 . This lemma results from substituting u = u mn in Eq. ( 18), and using the boundary conditions ( 7) and ( 8) to evaluate ∂u mn ∂ξ (see [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF] for the detailed proof). The approximate specific impedance (Z DtN2 mn ) corresponding to the second-order Dirichlet-to-Neumann absorbing boundary condition (DtN2) is given by:

Z DtN2 mn = λ 01 -λ 00 (λ 01 -λ mn ) 1 Z ex3 00 + (λ mn -λ 00 ) 1 Z ex3 01 (22)
where λ mn := λ mn (kf ) is called the prolate spheroidal eigenvalue (see p.11 in Ref. [4]).

Remark 4.2 Note that when e = 0, that is the prolate spheroid degenerates to sphere, the approximate DtN specific impedances given by Eqs. ( 21) and ( 22) are identical to the ones obtained in the case of spherical-shaped radiators [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]. In addition, when m = n = 0 (resp. (m, n) = (0, 1)), Z DtN2 mn coincides with Z ex3 00 (resp. with Z ex3 01 ), as expected.

The next proposition states the asymptotic behavior of the approximate specific impedance Z DtN2 mn as ka → 0.

Proposition 4.3 The asymptotic behavior of the DtN2 specific impedance of the mn th (m ≤ n) prolate spheroidal mode as ka → 0 is given by:

Z DtN2 mn ∼ -2 n(n + 1) -2 (n(n + 1) + 2) 2 (ka) 2 -i 2ka n(n + 1) + 2 ; n ≥ 2 (23)
Proof of Proposition 4.3. The asymptotic behavior of the DtN2 specific impedance given by Eq. ( 23) results from substituting the asymptotic behavior of the exact specific impedance given by Eq. (20) into Eq. ( 22) and from using the following property of the prolate spheroidal eigenvalues (see p.11 in Ref. [4]):

λ mn (0) = n (n + 1) (24)
Remark 4.4 First, we note that the asymptotic behavior of Z DtN1 mn (resp. Z DtN2 mn ; n ≤ 1) is identical to the behavior of the exact specific impedance Z ex3 00 (resp. Z ex3 mn ; n ≤ 1), which follows from the design of the DtN boundary conditions. Second, for n ≥ 2, the asymptotic behavior for DtN2 is independent of the eccentricity e, and is identical to the case of a sphere (see Eq. (92), p. 32 in [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]) as well as to the approximate specific impedance corresponding to the second-order BGT condition when applied on prolate spheroid radiator (see Eq. (69), p. 3646 in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]). Last, the asymptotic behavior given by Eq. (20) and Eq. (23) indicates that for higher modes (n ≥ 2), the real part of the exact impedance tends to zero, as ka → 0, faster than the DtN2 specific impedance. This result suggests that DtN2 boundary condition may not be appropriate for single higher modes.

Numerical investigation

We have performed several experiments to investigate numerically the effect of the wavenumber and the slenderness of the boundary on the performance of the second-order DtN boundary condition DtN2 given by Eq. ( 8) when computing single radiating modes. We have INRIA compared the results to the ones obtained with BGT2 conditions when applied on prolate spheroid boundaries (see Figs. [START_REF] Turkel | Iterative methods for the exterior Helmholtz equation including absorbing boundary, conditions[END_REF] to (21) in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]). We report in this section the results obtained when computing modes (m, n) = (0, 2) and (m, n) = (1, 1) for illustration. More numerical results can be found in [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF]. Note that for lower modes, DtN2 is exact, and therefore outperforms the BGT2 absorbing boundary condition. The numerical results presented here have been obtained for two low wavenumber values ka = 0.1 (see Figs. ( 1), ( 2), ( 5), and ( 6)), and for ka = 1 (see Figs. ( 3), ( 4), [START_REF] Givoli | Exact representations on artificial interfaces and applications in mechanics[END_REF], and ( 8)). Because of the dependence of the BGT2 specific impedance on ϕ ∈ [0, π) (see Eq. ( 69) in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]), we have reported the obtained results on the magnitude of the specific impedance in polar coordinates, and the relative errors as a function of ϕ. These results are obtained for three eccentricity values e = 0.1 (a prolate spheroid boundary very close to a sphere), e = 0.6 (a "regular" prolate spheroid boundary), and e = 0.9 (the prolate spheroid boundary is very elongated). As expected (see Eq. ( 20) and Proposition 4.3), DtN2 boundary condition performs overall poorly when computing higher modes regardless of the wavenumber and eccentricity values since the relative error is above 25% for most situations (see Figs ( 2), ( 6), (4), and ( 8)). Note that the level of accuracy delivered by DtN2 boundary condition seems to be better than by BGT2 boundary condition (except for e = 0.6 when computing mode (m, n) = (0, 2)), yet this level of accuracy remains very high for practitioners. Last, when the prolate spheroid tends to be close to a sphere e = 0.1, both conditions retain an acceptable level of accuracy ( the relative errors are less than 10%).

The performance of the local DtN conditions for prolate spheroidshaped scatterer

Next, we assess analytically and numerically the performance of the approximate local DtN1 and DtN2 boundary conditions given by ( 7) and ( 8) when applied to three-dimensional scattering problems by sound-soft prolate spheroid-shaped obstacle . Recall (see Eq. (85), p. 3649 in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]), that the exact specific impedance Z ex3 on the surface of a prolate spheroid at ξ = ξ 0 is given by:

Z ex3 = - i √ 1 -e 2 ka u inc ∂ ∂ξ (u scat )| ξ=ξ0 (25)
where u inc is the incident plane wave with incident angle ϕ 0 defined by u inc = e ikf cosh ξ(cos ϕ cos ϕ0+tanh ξ sin ϕ sin ϕ0 cos θ) (26)

u scat is the acoustic scattered field by a sound-soft prolate spheroid-shaped obstacle. u scat can be represented by the following series [START_REF] Senior | Scalar diffraction by a prolate spheroid at low frequencies[END_REF]

u scat = -2 ∞ m=0 ∞ n=m (2 -δ 0m )A mn i n R (3) mn (kf, cosh ξ)S mn (kf, cos φ) cos mθ (27)
where

A mn = R (1)
mn (eka, e -1 )S mn (eka, cos φ 0 )

N mn R (3) mn (eka, e -1 ) (28)
and N mn represents the normalization factor of the angular spheroidal wave functions [4] and δ 0m is the Kronecker delta symbol. Moreover, we also recall (see Eq. (91), p. 3655 in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]), that the asymptotic behavior of the exact specific impedance Z ex3 of the scattered field on the surface of a prolate spheroid as ka → 0 is given by:

Z ex3 ∼ Z ex3 00 ∼ (ka) 2 -ika (29)
where Z ex3 00 is given by Eq. ( 20) for m = n = 0.

Analytical study

The next lemma states the expression of the approximate specific impedances on the surface ξ = ξ 0 of a prolate spheroid sound-soft scatterer. This lemma results from substituting u = -u inc in Eq. ( 18), and using the boundary conditions ( 7) and ( 8), and the expression of the incident plane wave given by Eq. ( 26) to evaluate ∂u inc ∂ξ at ξ = ξ 0 (see [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF] for the detailed proof).

INRIA

Lemma 4.5 The approximate specific impedance (Z DtN1 ) corresponding to the first-order Dirichlet-to-Neumann boundary condition (DtN1) is given by:

Z DtN1 = Z ex3 00
(30) The approximate specific impedance (Z DtN2 ) corresponding to the second-order Dirichlet-to-Neumann boundary condition (DtN2) is given by:

Z DtN2 = λ 01 -λ 00 λ 01 Z ex3 00 - λ 00 Z ex3 01 -2iαka -(ka) 2 δ -(eka) 2 cos ϕ (31)
where

α = cos ϕ cos ϕ 0 + √ 1 -e 2 sin ϕ sin ϕ 0 cos θ, δ = ∂α ∂ϕ 2 + 1 sin 2 ϕ ∂α ∂θ 2 (32)
Remark 4.6 Note that when e = 0, that is the prolate spheroid degenerates to sphere, the approximate DtN specific impedances given by Eqs. ( 21) and ( 22) are identical to the ones obtained in the case of spherical-shaped scatterers [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF].

The next proposition states the asymptotic behavior of the second-order DtN approximate specific impedance Z DtN2 as ka → 0. Proposition 4.7 The asymptotic behavior of the approximate specific impedance Z DtN2 as ka → 0 is given by:

Z DtN2 ∼ (1 -α) (ka) 2 -ika (33)
Proof of Proposition 4.7. The asymptotic behavior of the DtN2 specific impedance given by Eq. (33) results from substituting the asymptotic behavior of the exact specific impedance given by Eq. (20) into Eq. (31) and from using the the property of the prolate spheroidal eigenvalues given by Eq. (24).

Remark 4.8 First, we note that the asymptotic behavior of Z DtN1 is identical to the behavior of the exact specific impedance Z ex3 (see Eq. (29) and Eq. ( 30)). Second, Eq. (33) indicates that the asymptotic behavior of Z DtN2 depends on the eccentricity as well as on the observation angle ϕ. This dependence is comparable to the asymptotic behavior of the BGT2 approximate specific impedance (see Eq. (93), p. 3657 in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]). Last, when e = 0, the asymptotic of both Z DtN1 and Z DtN2 are identical to the case of spheres (see Eq. (140) p. 47 in [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF]).

Numerical investigation

We have performed several experiments to investigate numerically the effect of the wavenumber and the slenderness of the boundary on the performance of the second-order DtN boundary condition DtN2 given by Eq. ( 8) when solving sound-soft scattering problems in the OSRC context. We have compared the results to the ones obtained with BGT2 conditions when applied on prolate spheroid boundaries (see Figs. ( 22) to (33) in [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF]). We report in here the results obtained for three different values of the incidence angle ϕ 0 = 0, π 4 and π 2 for illustration. More numerical results can be found in [START_REF] Saint-Guirons | Construction et analyse de conditions aux limites absorbantes pour des problèmes de propagation d'ondes[END_REF]. Moreover, since the impedences in this case depend on ϕ ∈ [0, π) and θ ∈ [0, 2π), we have reported 3D plots of obtained results in Figs. ( 9)-( 14). These results have been obtained for three eccentricity values e = 0.1 corresponding to a prolate spheroid boundary very close to a sphere (see Fig. [START_REF] Harari | Computational methods for Problems of Acoustics with Particular Reference to Exterior Domains[END_REF] and Fig. [START_REF] Keller | Exact nonreflecting boundary conditions[END_REF]), e = 0.6 corresponding to a "regular" prolate spheroid boundary (see Fig. [START_REF] Harari | Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior domains[END_REF] and Fig. ( 13)), and e = 0.9 (corresponding to a very elongated prolate spheroid boundary (see Fig. [START_REF] Harari | Analytical study of the effect of wave number on the performance of local absorbing boundary conditions for acoustic scattering[END_REF] and Fig. ( 14)). One can observe that unlike the radiator problem, the DtN2 boundary condition delivers an excellent level of accuracy in the low frequency regime for all eccentricity values. Indeed, the results depicted in Figs. [START_REF] Keller | Exact nonreflecting boundary conditions[END_REF], [START_REF] Kriegsmann | A new formulation of electromagnetic wave scattering using an on-surface radiation boundarycondition approach[END_REF], and [START_REF] Reiner | The performance of local absorbing boundary conditions for acoustic scattering from elliptical shapes[END_REF] indicate that the relative error on the approximate DtN2 specific impedance is always below 2%, while The BGT2 accuracy deteriorates significantly for e ≥ 0.6 (the relative error is larger than 40%).

INRIA 5 Conclusion

We have designed a new class of approximate local ABCs to be applied on elliptical-shaped exterior boundaries when solving acoustic scattering problems by elongated obstacles. These conditions are exact for the first radiation modes, they are easy to implement and to parallelize, and they preserve the local structure of the computational finite element scheme. The analysis reveals that in the case of the radiator, DtN2 boundary conditions, by construction, outperforms the BGT2 conditions while performing similarly for small values of e and the DtN2 boundary condition outperforms BGT2 for lower single modes. However, for higher modes both conditions performs poorly. In the case of the scattering problem, the situation is different. Indeed, the BGT2 boundary conditions performs poorly for eccentricity values e ≥ 0.6, while the DtN2 boundary condition delivers an excellent level of accuracy in the low frequency regime for all eccentricity values. We plan to investigate the effect of large values of the wavenumber on the accuracy of the DtN2 condition.

Lemma 4 . 1

 41 The approximate specific impedance (Z DtN1 mn ) corresponding to the first-order Dirichlet-to-Neumann absorbing boundary condition (DtN1) is given by:

m = 0 n = 2 kaFigure 1 :Figure 2 :Figure 3 :

 2123 Figure1: Absolute value of the specific impedance for the exact (solid), the DtN2 (crossed), the BGT2 (dashed).

Figure 4 :Figure 5 :

 45 Figure 4: Relative error of the specific impedance for the DtN2 (crossed), the BGT2 (dashed).

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Relative error of the specific impedance for the DtN2 (crossed), the BGT2 (dashed).

Figure 9 : 4 (Figure 10 : 4 (Figure 11 : 4 (

 94104114 Figure 9: Absolute value of the specific impedance for the exact (clear grey), the DtN2 (black), the BGT2 (dark grey) for ka = 0.1 (left) and ka = 1 (right) for the incident angles θ 0 = 0 (up), θ 0 = π 4 (middle), θ 0 = π 2 (down).

Figure 12 : 2 (

 122 Figure 12: Relative error of the specific impedance for the DtN2 (black), the BGT2 (dark grey) for ka = 0.1 (left) and ka = 1 (right) for the incident angles θ 0 = 0 (up), θ 0 = π 4 (middle), θ 0 = π 2 (down).

Figure 13 : 2 (

 132 Figure 13: Relative error of the specific impedance for the DtN2 (black), the BGT2 (dark grey) for ka = 0.1 (left) and ka = 1 (right) for the incident angles θ 0 = 0 (up), θ 0 = π 4 (middle), θ 0 = π 2 (down).

Figure 14 :

 14 Figure 14: Relative error of the specific impedance for the DtN2 (black), the BGT2 (dark grey) for ka = 0.1 (left) and ka = 1 (right) for the incident angles θ 0 = 0 (up), θ 0 = π 4 (middle), θ 0 = π 2 (down).

RR n 6395

Acknowledgements

The authors acknowledge the support by INRIA/CSUN Associate Team Program. Any opinions, findings, conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of INRIA or CSUN.

RR n 6395


Christophe Gourlaouen

Sergi Vela

Sylvie Choua

Mathilde Berville

Jennifer A Wytko

Jean Weiss

Vincent Robert

Pairing-up viologen cations and dications: a microscopic investigation of van der Waals interactions

The microscopic origin of van der Waals-and magnetic-interactions in 4,4 0 methyl viologen cationbased units (MV + and MV 2+ ) was inspected using wave function (variational DDCI and perturbative MP2, CASPT2) and density functional theory (DFT) calculations. The analysis deepens the comprehension of the magnetic behavior of experimental bis-viologen cyclophanes ([CYC] 2(+) ), in which the MV + units are connected through alkyl linkers of different lengths. The formation of the so-called long-multicenter bonds in such radical dimers, responsible for the quenching of the magnetic response, was analyzed in +) . Dynamical correlation effects, accessible from second-order perturbation corrections, were decisive in observing a bonding regime characterized by an equilibrium distance of 3.3 Å and a 45 kJ mol À1 dissociation energy. At larger intermolecular distances, our calculations on [MV 2 ] 2(+) indicate that the singlet and triplet states are energetically competing (i.e. weak exchange interactions, J AB ). Despite the absence of any clear bonding regime at the MP2 level, the puzzling association of two di-cations into [MV 2 ] 4+ is anticipated at 3.3 Å using weakly screened point charges (e = 1.5) to account for the Coulomb interactions between the solvated subunits. The main conclusion is that both dispersion interactions and environment effects are required to overcome the Coulomb repulsion associated with doubly-charged species. All these data provide some complementary insights into the nature and amplitude of interactions between cation and dication units, and their relevance in various experimental manifestations.

Funding information

Providing accurate funding information will enable us to help you comply with your funders' reporting mandates. Clear acknowledgement of funder support is an important consideration in funding evaluation and can increase your chances of securing funding in the future.

We work closely with Crossref to make your research discoverable through the Funding Data search tool (http://search.crossref.org/funding). Funding Data provides a reliable way to track the impact of the work that funders support. Accurate funder information will also help us (i) identify articles that are mandated to be deposited in PubMed Central (PMC) and deposit these on your behalf, and (ii) identify articles funded as part of the CHORUS initiative and display the Accepted Manuscript on our web site after an embargo period of 12 months.

Further information can be found on our webpage (http://rsc.li/funding-info).

What we do with funding information

We have combined the information you gave us on submission with the information in your acknowledgements. This will help ensure the funding information is as complete as possible and matches funders listed in the Crossref Funder Registry.

If a funding organisation you included in your acknowledgements or on submission of your article is not currently listed in the registry it will not appear in the table on this page. We can only deposit data if funders are already listed in the Crossref Funder Registry, but we will pass all funding information on to Crossref so that additional funders can be included in future.

Please check your funding information

The table below contains the information we will share with Crossref so that your article can be found via the Funding Data search tool. Please check that the funder names and grant numbers in the table are correct and indicate if any changes are necessary to the Acknowledgements text.

Funder name

Funder's main country of origin Funder ID (for RSC use only)

Award/grant number

Researcher information

Please check that the researcher information in the table below is correct, including the spelling and formatting of all author names, and that the authors' first, middle and last names have been correctly identified. Names will be indexed and cited as shown on the proof, so these must be correct.

If any authors have ORCID or ResearcherID details that are not listed below, please provide these with your proof corrections.

Please ensure that the ORCID and ResearcherID details listed below have been assigned to the correct author. Authors should have their own unique ORCID iD and should not use another researcher's, as errors will delay publication.

Please also update your account on our online manuscript submission system to add your ORCID details, which will then be automatically included in all future submissions. See here for step-by-step instructions and more information on author identifiers. For your information: You can cite this article before you receive notification of the page numbers by using the following format: (authors), Phys. Chem. Chem. Phys., (year), DOI: 10.1039/c8cp04543g.

Editor's queries are marked on your proof like this Q1, Q2, etc. and for your convenience line numbers are indicated like this 5, 10, 15, ...

Please ensure that all queries are answered when returning your proof corrections so that publication of your article is not delayed.

Query reference

Query Remarks

Q1

Please confirm that the spelling and format of all author names is correct. Names will be indexed and cited as shown on the proof, so these must be correct. No late corrections can be made.

Introduction

Organic-based radicals have attracted much attention due to their ubiquitous character, being capable to act as spin holders, 1 switchable materials 2,3 and binding ligands, 4 not to mention their ability to form three dimensional networks. 5 From their flexible electronic structure, unusual properties including conductivity, magnetism, spin crossover or charge ordering have been observed and stimulated intensive work. [6][7][8][9][10][11] Furthermore, changes in the spin states can be controlled by manipulating the unpaired electrons, taking advantage of their intrinsic degree of freedom. 12 Bistability in organic radical species is largely dominated by their tendency to dimerize, quenching the desired property of interest while forming a diamagnetic ground state. 13,14 In some cases, the (apparently) non-dimerized paramagnetic phase has also been found to be the result of a fast dimerization thus giving the impression that they are, indeed, not dimerized. 15 Such dimerization is also active in some charged (monooxidized or reduced) species, where the repulsive Coulomb contribution is overcome by interactions that have been attributed to overlaps between the p-SOMOs (singly-occupied molecular orbitals) centred on each radical cation. 16 This type of interaction presents spectroscopic similarities to covalent bonds but differs in their nature, strength, and equilibrium distance. [17][18][19][20] In this context, paraquat cations (e.g. 1,1 0dimethyl-4,4 0 -bipyridinium) within the viologen family, are particularly attractive platforms as the formation of so-called p-bonds is likely to reversibly generate dimeric units under redox stimulus. 21 The bipyridinium radical-cations are held together in sandwich-like arrangements displaying much shorter distances (ca. 3.0 Å) than the standard p-p stacking distance (3.6 Å). More recently, synthetic strategies have been developed to vary the number of alkyl linkers between viologen units (cyclophane structures [CYC] 2(+) leading to compounds 1, 2 and 3, see Fig. 1), and to ultimately modulate the host-guest properties. 21 In the cyclophane units, the cation radicals are linked by two alkyl chains with either five, six or seven carbon atoms (see Fig. 1).

The p-dimerization of these species was analyzed by means of electrochemistry, vis/NIR spectroscopy and spectroelectrochemistry. In the meantime, the electronic structures were tentatively inferred by EPR spectra. A well-resolved EPR signal could only be observed for 2 in fluid solution revealing the presence of a biradical with a strong antiferromagnetic interaction (ca. À734 cm À1 ). 22 Unfortunately, it was not possible to measure a clear EPR signal for compounds 1 and 3. Quite interestingly, the crystal structure of 2 was obtained for the tetra-cation [CYC] 4+ , 23 a regime where doubly-charged organic subunits are held together. Density functional theory (DFT) calculations were then performed on isolated cyclophanes in a model of water, stressing a marked propensity towards antiferromagnetism. 22 More recently it has been suggested that large dispersion forces dominate the bonding regime in viologen cation-radical p-dimer (so-called pimer). 24 The analysis was carried out from DFT calculations including dispersion corrections and concluded on the competition between covalent and dispersion attractions with Coulomb repulsions.

The heart of the current investigation is whether the strong antiferromagnetism and packing are a consequence of the screened charges in [CYC] 2(+) and [CYC] 4+ units or must be attributed to dispersion interactions. Therefore, we felt that a microscopic analysis would be instructive to rationalize both the p-dimerization and the magnetic properties of this family of cations and radical cations. Starting from experimental X-ray data on [CYC] 2(+) (see Fig. 1), viologen radical-based compounds were examined using DFT and wave function calculations. We first concentrated on the [MV 2 ] 2(+) radical model, with the aim of understanding the origin of the interaction at van der Waals distances. Based on the structure of the wave function, the relative importance of the so-called ionic and neutral forms is discussed in comparison to the C 2 H 6 molecule, taken as a prototype of a ''two electrons in two MOs'' covalentlybonded molecule. Then, a similar inspection was carried out on the tetracationic form [MV 2 ] 4+ . The puzzling binding of two doubly-charged MV 2+ units, as observed in the crystal structure of 2, deserves particular attention in order to understand its nature. Finally, the magnetic properties of [MV 2 ] 2(+) were examined based on calculations of magnetic exchange couplings on the optimized dimer geometries. Due to the flexible nature of the alkyl linkers in the cyclophanes, different conformations were expected to be present in solution. These inspections allow us to understand how the relative positions of the MV + units modulate the nature and amplitude of the magnetic interactions.

In this framework, our intention is to examine the dimerization process holding both viologen units at rather short distances in comparison with standard p-p stacking distances.

Technical background

For [MV 2 ] 2(+) dication species, results were first obtained from complete active space self-consistent field (CASSCF) calculations including two electrons in two molecular orbitals (MOs), namely CAS [2,2]. All these calculations were performed using the MolCAS 7.2 package. 25 The atoms were described with extended basis sets of ANO-RCC types, [4s3p2d1f] for C and N atoms, and [2s1p] for the H ones. 26 Based on the CAS[2,2]SCF singlet ground state wave function, complementary second-order perturbation calculations (CASPT2) were performed. Total energies were evaluated at the CAS[2,2]PT2 level with respect to r z , the distance between two viologen units, and the structure of the wave function was analyzed. Such a framework allows one to build the wave function and to extract the information relative to the important dynamical correlation effects. The active space was extended to CAS [6,6] with no quantitative change in the calculated CASPT2 binding energies (42 kJ mol À1 for CAS [6,6]PT2 and 45 kJ mol À1 for CAS[2,2]PT2). Calculations were restricted to MP2 type (Mo ¨ller-Plesset) for the closed-shell tetra-cation [MV 2 ] 4+ analogues.

Second, the nature of the radical interaction was examined through a configuration interaction (CI) restricted to single excitations (CAS[2,2]SCF + S). These CI calculations were carried out with the CASDI code 27 to retrieve the leading charge redistribution beyond the mean field CASSCF picture. For comparison with a prototype of covalent bond, we performed CAS[2,2]SCF + S calculations for the C 2 H 6 ethane molecule at equilibrium and stretched carbon-carbon distances (1.54 and 3.0 Å) using a [3s2p1d] and [2s1p] basis sets for carbon and hydrogen atoms, respectively. At 3.0 Å, the active space consists in the bonding and antibonding s-type MOs based on the CH 3 SOMOs. This set of MOs was then used as a starting set to converge the CAS [2,2] solution at a 1.54 Å carbon-carbon distance. From the presence of a center of inversion, the active MOs are the in-phase (u) and out-of-phase (g) linear combinations of the SOMOs localized on each individual unit. It is wellknown that the CAS[2,2]SCF + S singlet wave function mostly reads C = l|g% g| À m|u % u| where g and u are the canonical bonding and antibonding MOs (see Fig. 2a andb). A unitary transformation allows one to rewrite C in the local atomic basis set to evaluate the ratio between the so-called ionic and neutral form as r = (l À m)/(l + m). 28 The diradical character is featured by a reduction of r, and can be considered as a measure of covalency to be compared with the so-called p-dimerization mechanism. Third, in the presence of two unpaired electrons in a dimer of viologen cation radicals, low-lying singlet (total spin S = 0) and triplet (S = 1) states are expected to compete. The energy difference defines the so-called magnetic exchange coupling constant 2J AB = E S À E T ruling the Heisenberg spin-Hamiltonian H = À2J AB S 1 S 2 , where S 1 and S 2 are the 3D-spin operators acting on each subunit, and are associated with spin values 1 2 . Two numerical approaches were used to evaluate the exchange coupling constants. Difference dedicated configuration interaction (DDCI) calculations were carried out using the CASDI code. 27 Such methodology has proven to provide not only very accurate estimates of magnetic exchange coupling constants but also microscopic interpretations. 29,30 The atomic basis sets were safely reduced down to [3s2p1d] for carbon and nitrogen atoms, and [2s] for hydrogen ones based on similar inspections carried out on verdazyl organic radicals. 31 Throughout the DDCI calculations, the triplet CAS[2,2] MOs were used. We checked that any modification of the set of MOs leads to negligible changes in the singlet-triplet energy difference. To reduce the computational cost, the two alkyl chains of the cyclophane units were substituted by four methyl groups (yielding pairs of cation viologen units, referred to as [MV 2 ] 2(+) and [MV 2 ] 4+ , respectively). We confirmed that this simplification has a limited impact on the calculated energy differences. Quantitatively, the triplet-singlet energy difference at the DDCI level are 4697 and 4794 cm À1 for the full and simplified diradicals at 2.9 Å distance, respectively (i.e. less than 3% relative difference). This negligible modification is coherent with the fact that such alkyl chains are poor mediators of magnetic interaction in diradicals [CYC] 2(+) .

Taking advantage of the relative low-cost of the DFT calculations, different degrees of freedom (interplanar distance, slippages and rotations) were considered to derive magnetostructural analysis. The broken-symmetry (BS) approach 32,33 provides the energies of a fictitious open-shell state and a high-spin state. Different strategies were proposed to evaluate the exchange-coupling constant from these values. 34 The relative singlet-triplet energy differences between DDCI and DFT calculations are less than 5%. All calculations were performed using the B3LYP functional (restricted or unrestricted), and the atoms were described by TZVP basis sets as implemented in Gaussian 09 [START_REF] Frisch | Gaussian 09, Revision D.01[END_REF] in gas phase.

Results and discussion

Wave function theory based calculations were performed (i) to analyze the bonding mechanism, (ii) to trace the stability of dication [CYC] 2(+) and tetra-cation [CYC] 4+ cyclophane units, and (iii) to provide exchange coupling constant values in the diradical [CYC] 2(+) .

Nature of the interaction between two MV + units

Despite the well-accepted terminology ''pimer'', the interaction has been viewed as resulting, at least in part, from s-overlap between p systems. 16 To clarify the nature of the interaction between the viologen cation radicals MV + , the variations of CAS [2,2]SCF and CAS[2,2]PT2 total energies were followed with respect to the intra-planar distance (r z ), maintaining a perfect stacking between the viologen units. As expected, the active MOs are the in-phase (u) and out-of-phase (g) linear combinations of the SOMOs localized on each individual unit (see Fig. 2a andb). Fig. 2c displays the calculated energy of the ground state singlet performed at various levels (CAS[2,2]SCF and CAS[2,2]PT2). To compare the energy variations, the r z = 10 Å values was set to zero as a common reference. Despite the presence of two unpaired electrons, the interaction cannot be understood from a ''two electrons in two MOs'' picture, that might be qualitatively acceptable for neutral radicals. [START_REF] Koivisto | [END_REF] The absence of any binding regime at the CASSCF level suggests that the interaction does not fall in the traditional covalent bond scenario. The one-electron picture constructed on the SOMOs overlap is not appropriate to qualitatively account for the weak bond formation. As soon as second-order perturbation corrections are included, the interactions between units considerably change. A marked minimum is observed for r z = 3.3 Å (see Fig. 2), associated with a binding interaction of 45 kJ mol À1 . In the present case, the bonding character is reached with dynamical correlation effects, similarly to what is wellknown between closed-shell systems such as benzene rings (i.e. dispersion-like interactions). At this stage, one should mention that shorter equilibrium distances (ca. 3.0 Å) have been reported experimentally. 16 Solvent molecules in solution as well as counter-anions in the solid state are likely to screen the Coulomb repulsions, allowing closer inter-molecular contacts.

To gain further insight into the nature of this interaction, the ratio between the ionic and neutral amplitudes (see com- From the wave function analysis, this quantity amounts to 0.27 in [MV 2 ] 2(+) . This is significantly smaller than the values of 0.84 and 0.32 found in C 2 H 6 at equilibrium and stretched carboncarbon distances, respectively. From this inspection, the contribution of the ionic forms is reduced by a factor of three in the [MV 2 ] 2(+) cation as compared to the value in C 2 H 6 at equilibrium distance, and still smaller than the one in the dissociated regime. This observation is consistent with the prime role played by the electronic correlation effects. This highlights the radical character with vanishingly-small electron delocalization in the interaction between the two MV + units.

Binding of two MV 2+ units?

We then moved to the tetracationic derivatives [MV 2 ] 4+ which was reported in the crystal structure of 2. 21 Multiply-charged species are not so scarce in solid state chemistry. In such species, the chemical bonding in an ionic solid is dominated by electrostatic interactions. Even at the electrostatic level, the Coulomb interactions between units are screened thanks to the environment. 37,38 In addition, the intrinsic polarizabilities of the building blocks are likely to give rise to stabilizing interactions. From previous investigations, and in the light of our current results on the [MV 2 ] 2(+) cation, the dispersion forces significantly overcome the Coulomb repulsion since a 45 kJ mol À1 binding interaction is calculated in vacuum (see previous section). How strong are dispersion interactions in this class of compounds, where highly charged species have been observed? It was rather tempting to assess the ability of such contributions to generate a tetra-cation from two closed-shell MV 2+ units in the absence of any counter ion. Thus, secondorder Mo ¨ller-Plesset calculations (MP2) were carried out by varying the interplanar distance r z . In sharp contrast with the parent [MV 2 ] 2(+) , no binding regime is observed for [MV 2 ] 4+ : the Coulomb repulsion exceeds the dispersion interactions in the examined r z distance window in vacuum. Nevertheless, no matter the phase (condensed: liquid or solid), solvent molecules as well as counter ions act as mediators of the pdimerization. We can evaluate their participation through the introduction of a dielectric constant e. The electrostatic interactions are effectively reduced as a consequence of the charge screening. Therefore, the Coulomb repulsion was first corrected to account for these effects through the following term: (À4/r z ) (1 À 1/e), and e was varied starting from an e = 1 value (i.e., vacuum). Formally, this quantity corresponds to the change in the electrostatic repulsion between two 2+ charges upon going from vacuum to a dielectric medium characterized by e. As soon as a rather weak e = 1.5 value is reached, the MP2 potential energy curve suddenly exhibits a clear minimum for r z = 3.5 Å (Fig. 3, curve MP2 + e). The resulting binding energy is of the order of 40 kJ mol À1 , which is a rather large value considering the nature of the system. At this stage, one may first question the amplitude of the binding energy compared to the value found for the dication [MV 2 ] 2(+) . Evidently, the interplanar distance being very similar, the screened Coulomb correction is four times larger than it is in the dication. Thus, one may then argue that this correction may mask the contributions of the dispersion-like interactions. However, we checked that the potential energy curve remains dissociative at the Hartree-Fock level even if these contributions are included (Fig. 3, HF + e curve). Despite the questionable validity of such a description at rather short distances, the impact of the electrostatic continuum was also estimated by screening the correlation contributions, a second order effect. The energy differences between the MP2 and HF values were renormalized by 1/e 2 (see Fig. 3, HF + e + D curve). Such procedure underestimates the correlation effects since both inter-and intrasubunit contributions are screened. Still, a bonded geometry is observed for r z = 3.9 Å (and a non-negligible binding energy of ca. 6 kJ mol À1 ). Finally, our views were complemented using HF calculations in the presence of a PCM, followed by a perturbative treatment. The PCM polarizes the MOs and account for the MOs distortion induced by the solvent. Such procedure allows one to grasp part of the dispersion-induced charge rearrangements, the importance of which has already been stressed. 39 In water, a binding regime is observed, which confirms that both screening and dispersion are required to hold the 2+ units together. In other words, the viologen polarization remains decisive to establish a bonding regime, which cannot be reached using a mean-field description.

Magneto-structural analysis

Starting from the DFT optimized geometry of [MV 2 ] 2(+) , characterized by an inter-molecular distance r z of 3.6 Å and a co- facial arrangement (r x = 0), the J AB values were calculated for two translations, the interplanar distance (r z ) and the longitudinal slippage (r x ), and two rotations, a around the stacking (z axis) and b around the x axis (see Fig. 4). Other variables were found to be less relevant. We did not intend to provide a full potential energy surface of the low-lying singlet and triplet states, but rather to evaluate the relative changes in the exchange-coupling constant in the vicinity of the equilibrium geometry. As expected, the J AB intensity increases exponentially from ca. À90 to 0 cm À1 as the r z value is varied from 3.5 to 6 Å (see Fig. 4a). Such behaviour was analysed in similar compounds and can be attributed to the exponential variation of the SOMOs overlap and hopping integrals. 31,40 At this stage, let us mention that the available crystal structures for compounds 1 and 3 exhibit interplanar distances which are larger than 4.8 Å. Thus, one would expect quasi-degenerate singlet and triplet, as reflected by a negligibly small J AB value. The non-linear response of J AB with respect to r z indicates that thermal oscillations would result in a significantly different mean value for J AB due to dynamic effects, making it more antiferromagnetic than the value obtained using a static picture. 41 In compound 2, such phenomena would support the EPR measurements of a J AB value of À734 cm À1 . 22 In contrast, the J AB variations along the longitudinal slippage exhibit a double minimum curve (see Fig. 4b). The first one is observed for r x = 0 Å, with a J AB value reaching ca. À700 cm À1 .

For this ''on-top'' geometry, the SOMOs overlap is optimal through facing pairs of pyridinium rings. Since the superexchange mechanism grows with the square of the overlap the resulting AFM contribution is then maximized. When r x deviates from zero, the super-exchange mechanism is disrupted and J AB increases until r x reaches 4 Å. For this value, the overlap between units is restricted to a single pair of pyridinium rings (see Fig. 4b, right inset) and, thus, the resulting J AB value should be roughly a quarter of the initial one. In this case, it is À250 cm À1 .

Let us now examine the effects of rotations on the exchange coupling constant. Along the z axis (a parameter in Fig. 4e), the curve exhibits two wells for a = 01 and a = 901 associated with J AB values of À1400 to À800 cm À1 , respectively (see Fig. 4c), and a maximum for a = 401. As previously mentioned, these variations reflect the modulation of the overlap between the SOMOs. Due to the intrinsically delocalized nature of the latter (see Fig. 2a), a = 901 corresponds to a local minimum for the exchange coupling constant. Upon variation of b (see Fig. 4d), J AB follows an exponential law, being directly controlled by the superexchange decay along this distortion. 42,43 To conclude, the nature of the coupling between MV + units mostly falls in the antiferromagnetic regime, with rather large amplitude variations over the explored coordinates. It is only for r x = 2 Å and r x = 6 Å slippage values that the singlet and triplet states become quasi-degenerate.

Conclusions

In the present work, we shed some light on the electronic structure of [CYC] 2(+) architectures using quantum chemistry calculations. Leaving out the alkyl linkers, we first showed the persistence of the radical nature of the MV + units even when van der Waals interactions occur. This was done by computing and analyzing the CASSCF and CASPT2 dissociation curves. In particular, we stressed the negligible weight of the charge transfer electronic configurations (i.e. ionic forms) even at short distances, suggesting that the radical character of each [MV 2 ] 2(+) unit was preserved. Thus, the interaction in [CYC] 2(+) architectures cannot be restricted to a one-electron picture of overlapping SOMOs. In solvent conditions where additional screening effects should included, aggregation is observed at shorter distances as reflected by strong antiferromagnetic interactions. This led us to examine the tetracationic derivative [MV 2 ] 4+ . The appearance of a binding regime results from the competition between dispersion interactions and moderately screened Coulomb repulsions. The introduction of an e = 1.5 constant value to account for the presence of solvent molecules and/or anions significantly stabilizes a pair of MV 2+ cations. This observation supports that the combination of dispersion interactions and of screening are mandatory to allow the binding between the doubly-charged species.

From a magnetic point of view, the alkyl linkers are not relevant when computing the J AB values (DDCI calculations), justifying the transferability of the results obtained for [MV 2 ] 2(+) to the target [CYC] 2(+) architectures. Other deformations from the perfect stacking reduce the strength of J AB , with a dominating antiferromagnetic behavior. From available X-ray data, one can conclude that isolated dimers cannot create strong magnetic interactions. Even if the larger alkyl linkers provide the cyclophane more flexibility, the geometries that can be potentially adopted by these compounds are associated with weak magnetic interactions. Therefore, the only consistent explanation for the presence of strong antiferromagnetic couplings is that distances smaller than the calculated equilibrium distance of 3.3 Å can be reached between two different [CYC] 2(+) units. Despite the positive charges in both units, solvent molecules can be mediators of the p-dimerization, as extensively reported experimentally 13,14 and rationalized computationally. 44,45 In this class of compounds, screening and dispersion effects act hand-in-hand to generate unusual regimes.

First

  (given) and middle name(s) Last (family) name(s) ResearcherID ORCID iD Queries for the attention of the authors Journal: PCCP Paper: c8cp04543g Title: Pairing-up viologen cations and dications: a microscopic investigation of van der Waals interactions

Fig. 1

 1 Fig. 1 Cyclophane structure ([CYC] 2(+ ) ) consisting of two viologen cation radical units linked together through alkyl chains of variable size. Compounds 1, 2 and 3 correspond to n = 1, 2 and 3, respectively.

Fig. 2

 2 Fig. 2 Triplet canonical CAS[2,2]SCF MOs (a) bonding, g symmetry, and (b) antibonding, u symmetry of [MV 2 ] 2(+ ) . (c) CAS[2,2]SCF and CAS[2,2]PT2 [MV 2 ] 2(+ ) energies (kJ mol À1 ) with respect to the interplanar distance r z (Å). The energy reference is taken at r z = 10 Å, arbitrarily.

  putational details) was calculated from the CAS[2,2]SCF + S wave function at the equilibrium distance r z = 3.3 Å. The r value was then compared to the values calculated for the reference C 2 H 6 molecule at equilibrium (1.54 Å) and stretched (3.0 Å) carbon-carbon distances. The wave function provides the s/s* electronic distribution resulting from the CH 3 SOMOs overlap.

Fig. 3 4 |

 34 Fig. 3 Dissociation curves of [MV 2 ] 4+ at HF and MP2 level of theory and influence of a dielectric constant (e) on the electrostatic interactions considering two 2+ point charges and on the correlation contributions (D). For each curve, the reference energy is set to zero for r z = 5 Å.

Fig. 4

 4 Fig. 4 Evolution of J AB values upon a change in the interplanar distance (r z , a), the longitudinal slippage (r x , b), and a rotation along the stacking (a, c) and longitudinal (b, d) axes, defined in (e) computed in DFT level of theory. The interplanar distance is set at 3.6 Å for b, c and d, which is the equilibrium distance in these conditions.

6 |

 6 Phys. Chem. Chem. Phys., 2018, 00, 1À7 T h i s j o u r n a l i s c the Owner Societies 2018 c the Owner Societies 2018 Phys. Chem. Chem. Phys., 2018, 00, 1À7 | 7 PCCP Paper

Acknowledgements

We thank the HPC centre of Strasbourg for computing time. We are grateful to Dr Manuel Ruiz-Lopez for fruitful discussions.

Conflicts of interest

There are no conflicts to declare.


Mathilde Berville

Jimmy Richard

Monika Stolar

Sylvie Choua

Nolwenn Le Breton

Christophe Gourlaouen

Corinne Boudon

Laurent Ruhlmann

Thomas Baumgartner

Jennifer A Wytko

Jean Weiss

A Highly Stable Organic Radical Cation

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Herein, we describe two remarkable features of a viologen (1 2+ ) bearing methyl ester groups at the 2,2′,6,6′ positions. Not only do these four electron-withdrawing groups significantly facilitate both reduction potentials but they also stabilize the radical cation in organic solution, even in the presence of oxygen. Although viologen radical cations can be stabilized as dimers under argon at low temperature 4 or in host-guest systems, 14 and air-stable viologen radical cations encapsulated in catenanes have been reported by Stoddart,15 there are no other examples of simple viologen derivatives in a stable radical-cation state.

The tetraester viologen 1 2+ was synthesized in three steps s1 starting from chelidamic acid 16,[START_REF] Howath | A new Efficient Method for the Preparation of 2,6-Pyridinedimethyl Ditosylates from Dimethyl 2,6-Pyridinedicarboxylates[END_REF] 2 (Scheme 1), which was transformed to the methyl ester of the chloropyridine Single crystals of the dicationic species were obtained by slow diffusion of Et 2 O in MeCN. In the solid-state structure 20 f1 of 1 2+ (Figure 1), the bipyridinium rings are twisted with a 32°5 dihedral angle that is closer to the 36°twist of 1,1′diphenylviologen [START_REF] Porter | Isolation and Characterization of Phenyl Viologen as a Radical Cation and Neutral Molecule[END_REF] than to the nearly coplanar structures observed for most 1,1-dialkylviologen dications. The structure of 1 2+ also shows the possible formation of intramolecular hydrogen bonds between the ester groups and the N-methyl substituents (2.1-2.5 Å) or the pyridine protons (∼2.4 Å) (see Supporting Information).

The redox properties of 1 2+ were investigated by electrochemical methods. Cyclic voltammetry (CV) on a glassy carbon electrode in CH 3 CN solutions with tetrabutylammonium hexafluorophosphate (TBAPF 6 ; 0.1 M) as supporting f2 electrolyte showed three reversible reductions (Figure 2a). The first two one-electron reductions at -0.27 and -0.50 V versus Fc + /Fc are, respectively, attributed to the formation of the radical cation 1• + and the neutral 1 0 . The third reduction detected from the CV at -1.90 V versus Fc + /Fc is a twoelectron event that generates 1 2-. This reduction was almost absent in rotating disk electrode experiments (Figure 2b), suggesting that the neutral form (1 0 ) undergoes a relatively slow chemical reaction once formed. It is also possible that the generated 1 2-species is sparingly soluble in the 0.1 M TBAPF 6 CH 3 CN medium. Controlled-potential coulometry at -0.37 and at -0.60 V versus Fc + /Fc resulted in the respective exchange of 0.97 and 1.85 electrons per molecule of compound 1.

The first two reductions of 1 2+ are shifted by +560 and +740 mV, respectively, compared to those of MV 2+ (E red1,1/2 = -0.83 V and E red2,1/2 = -1.24 V vs Fc + /Fc) under the same experimental conditions. The significantly less cathodic values Electron paramagnetic resonance (EPR) studies confirmed 119 the paramagnetic, radical character of 1• + . The reaction of 1 2+ 120 with 1 equiv of ascorbic acid yielded a well-resolved EPR f5 121 spectrum (Figure 5) centered at g = 2.0032 corresponding to 122 two sets of four and six equivalent protons (A 1 ( 1 H) = 0.11 mT, 123 A 2 ( 1 H) = 0.43 mT), respectively, and one set with two 124 equivalent nitrogen atoms (A 1 ( 14 N) = 0.45 mT). The same 125 spectrum was obtained by electrolysis. The colorless solution 126 turned green upon one-electron reduction to generate the 127 radical cation 1• + . After the second reduction, the solution 128 turned red, and the intensity of the EPR spectrum decreased 129 gradually to give the silent neutral species 1 0 . 130 The air stability of the radical cation 1• + over time was also f6 131 investigated by EPR (Figure 6). After the radical cation was 132 generated by coulometry in MeCN, samples of the solution of 133 1• + were transferred to EPR tubes purged with argon or open 134 to air. The argon sample was stable for more than 40 h. Even 135 after three months, the intensity of the EPR spectrum of the 136 same solution remained constant. For samples under air, the intensity of the EPR signal of 1• + remained relatively stable over 40 h in two cases. In the third case (Figure 6, 1• + open air 3), the intensity of the EPR signal increased initially over the first 5 h and then decreased by ∼25% after 40 h. We assign the initial increase to the presence of some diamagnetic dimeric species (1• + ) 2 that dissociate into paramagnetic 1• + monomers within the first few hours. The presence of oxygen in the samples exposed to air was confirmed by a broad EPR signal (see Supporting Information). Under the same experimental conditions, the MV• + reference compound was considerably less stable under argon, as demonstrated by the rapid decay [START_REF]Although all kinetic experimental decays for MV• + showed the same trend, small variations were observed depending on the amount of residual oxygen in solution[END_REF] of the EPR signal within a few hours (Figure 6). In the presence of air, MV• + was very unstable, and no EPR signal was observed. Together, these data confirmed the high stability of

1• + .
Despite the stability of the radical cation, we were unable to grow single crystals of sufficient quality for X-ray diffraction.

Thus, optimizations of the geometries of 1 2+ and 1• + were performed without any symmetry restriction. For the dication 1 2+ , the twist between the pyridinium rings was reproduced (33°vs 32°), and the accord between experimental and calculated geometry remained quite satisfactory. A large highest occupied molecular orbital-lowest unoccupied molecular orbital (HOMO-LUMO) energy gap was calculated (4.60 eV) for 1 2+ , but the LUMO was stabilized by 0.7 eV compared to MV 2+ . These calculations corroborate with the electrochemical data.

When the dication 1 2+ is reduced to 1• + , the optimized geometrical parameters provide a more planar structure with a value of dihedral angle of 15°. This twist is higher than the 6°1 and 11°observed in the two independent molecules of the Xray structure of MV• + . 24 Comparison of these data suggests less electronic delocalization in 1• + than in MV• + . The calculated 1 H and 14 N isotropic constants (see Supporting Information) are in close agreement with the experimental values. As shown by the SOMO illustrated in Figure 3, the unpaired electron is essentially localized on the pyridine rings, on both N-methyl groups, and only slightly on the tetraester substituents. This observation supports the stabilizing effect of the four electronwithdrawing ester groups on the radical cation.

In conclusion, the tetraester derivative of methylviologen is a strong electron acceptor that converts easily and reversibly into its radical cation species. Because of the electron-attracting effect of the ester groups, the radical species is stable for months under argon and for hours in the presence of oxygen.

Figure 1 .

 1 Figure 1. Side views of the viologen 1 2+ in the solid state. Triflate anions and solvent molecules have been omitted for clarity.

Figure 2 .

 2 Figure 2. CV (a) and rotating disk electrode voltammetry (b) of 1 2+ and MV 2+ in CH 3 CN + 0.1 M TBAPF 6 : working electrode = glassy carbon, pseudoreference = Pt, counter electrode = Pt, internal reference = Fc. (a) v = 0.1 V/s; (b) v = 0.02 V/s, rotation: 1000 rpm.

Figure 3 .

 3 Figure 3. ADF view drawing of the Kohn-Sham SOMO of the radical cation 1• + .

  118

Figure 4 .

 4 Figure 4. UV-Vis-NIR monitoring of the reduction of 1 2+ . (inset) Corresponding CV at scan rate of 0.01 V s -1 . Solution (c = 1 × 10 -3 M) of 1 2+ in CH 3 CN + 0.1 M TBAPF 6 recorded with an optically transparent thin-layer (0.2 mm) electrochemical cell equipped with an Au mini-grid WE and CaF 2 optical windows, reference AgCl/Ag.

Figure 5 .

 5 Figure 5. X-band EPR spectrum of 1• + (a) in CH 3 CN at 10 -3 M at room temperature; (b) simulated spectrum.

Figure 6 .

 6 Figure 6. Double integrated EPR spectrum of 1• + and MV• + in CH 3 CN at 1 × 10 -3 M at room temperature under Ar and open to air.