Djamal Rabli
email: djamal.rabli@gmail.com

Ronald Mccarroll
email: ronald_mccarroll@upmc.fr

The low-lying adiabatic states of the K + 2 alkali dimer

Keywords: model potential, K + 2, ab-initio, adiabatic energies, molecular spectroscopic constants

 with those of Magnier and Frécon [2] based on a model potential approach reveals a number of serious disagreement concerning several excited states. In particular, the 5

Π g and 2 2 ∆ u states which are found to be repulsive by Magnier and Frécon, but attractive when ab-initio techniques are employed.

To clarify the origin of this disagreement, the adiabatic energies and spectroscopic constants are re-computed for the low-lying states of the K + 2 alkali dimer within a model potential framework. Contrary to the claims of Magnier and Frécon, the new results based on a model potential approach agree well with the ab-initio ones.

Introduction

During the past few decades, several studies in quantum chemistry have been developed to obtain an accurate description of potential energy surfaces of excited electronic states [3,4,5,6]. Such surfaces provide the starting point for the investigation either of non-adiabatic effects via avoided crossings or of the dissociation into different fragments along a minimum energy path in presence of an external electromagnetic field.

In the particular case of diatomic molecules, potential energy surfaces simply become potential energy curves, which can be computed by solving the electronic timeindependent Schrodinger equation for different internuclear separations in the Born-Oppenheimer approximation [7]. Morever, for diatomic molecular ions formed by a neutral atom and an ion, it is well known that the corresponding potential energy curves generally exhibit a certain number of avoided crossings which can be exploited to obtain cross sections of charge transfer and excitation processes which may occur between different molecular adiabatic states.

Over the past few years, potential energy curves of several diatomic molecules have been investigated with a special focus both on neutral and ionized alkali-dimer systems [8,9,10,11]. This particular interest is generally due to the recent development in the cold collision dynamics involving alkali diatomic systems.

In the present paper, we are concerned with the K + 2 alkali-dimer molecular cation, which has been the subject of several theoretical and experimental studies. From the experimental perspective, according to the measurements of Stwalley and Wang [12] using a spectroscopic technique of the photoassociation of ultra cold atoms, the estimated ground state dissociation energy D e of the K + 2 ion was found to be 6670 cm -1 , while in an experimental work using a supersonic alkali-metal beam with a sequential two photon ionization technique, developed by Leutwyler et al. [13], the dissociation energy was found to be equal to 6404 cm -1 . On the other hand, according to Broyer et al. [14] who employed a technique of resonant two and three-photon ionization, the K + 2 ground state 1 2 Σ + g , has harmonic and anharmonic vibration constants, ω e =73.40 cm -1 and ω e χ e = 0.2 cm -1 , respectively.

From the theoretical point of view, multi-electron ab-initio calculations have been performed by Konowalow and Rosenrantz [15] to compute the electronic structure of a few light alkali diatomic molecules and their molecular cations, such as Li + 2 and Na + 2 . Their calculated data was then scaled to predict qualitative results for the K + 2 molecular ion. On the other hand, Ilyabaev and Kaldor [16] used a Fock-space open-shell coupled cluster approach to compute the potential energy of the three lowest states, while Valance [17] employed a Hellmann-type pseudopotential to generate the potential energy curves of the six lowest Σ + states. In a recent study, electron affinity calculations based on the equation of motion coupled cluster method has been employed by Skupin et al. [18] to accurately describe the six (four Σ + and two Π) lowest lying electronic states.

A theoretical investigation of the first 58 adiabatic states including spin-orbit effects has been performed by Jraij et al. [19], whereas in an extensive study, Berriche et al. [1] employed an ab-initio approach to compute the adiabatic energies of states dissociating up to K + +K(6s) limit, for a wide range of internuclear separations (2 ≤ R ≤ 92) a.u.

In addition to the ab-initio and pseudo potential methods which provide most available data describing the K + 2 ion, model potential techniques were employed by Henriet and Masnou-seeuws [20] to generate the spectroscopic constants of the ground state. In an alternative study also based on model potential approach, Henriet [21] has been successful in providing the adiabatic energies of the first excited states for short and intermediate internuclear separations up to 20 a.u. But to the best of our knowledge, the most extensive work exploiting model potential techniques remains the one by Magnier and Frécon [2] who used a potential of Klapisch [22] to compute adiabatic energies for states dissociating up to K + +K(6s) limit for short and intermediate internuclear distances, 4 ≤ R ≤ 35 a.u.

While for the ground state studies based on different approaches indicate a global agreement, for several excited states an important disagreement is observed between theoretical results. For example, an important disagreement is observed when comparing the ab-initio results of Berriche et al. [1] and those of Magnier and Frécon [2] based on model potential techniques. In particular, the 5

2 Σ + u , 6 2 Σ + g , 6 2 Σ + u , 7 2 Σ + g , 7 2 Σ + u , 3 2 Π g , 4 2
Π g and 2 2 ∆ u states which are found to be repulsive by Magnier and Frécon, whereas these states exhibit a well pronounced energy minimum when ab-initio [1] techniques are employed.

As shown by Rabli and McCarroll [23,24] in previous studies treating the Li + 2 and NaLi + molecular ions, in principle model potential methods should be suitable for the 2 system within a model potential framework to generate the adiabatic energy curves and the molecular spectroscopic constants of molecular states dissociating up to the K + + K(6s) limit. Our aim is to clarify whether this disagreement between the abinitio and model potential methods originates from some conceptual defect of the model potential technique or whether there is some source of error in the calculations of Magnier and Frécon [2].

Except where stated, atomic units will be used throughout.

Theory

Model potential methods have been widely employed in describing bound states of alkaline atoms. In the particular case of the K + 2 molecular ion, the interaction of the valence electron with the closed shell K + (3p 6 ) is described by a model potential proposed by Klapisch [22] which has been used in a previous work by Magnier and Frécon [2]:

V model (r) = - q r - Z -1 r (e -ar + cre -br ), (1)
where r denotes the radial distance of the electron form the nucleus, q = 1 the ionic core charge and Z = 19 the atomic number. The parameters (a,b,c) shown in Table 2 are in units of inverse length au -1 , and have been optimised such that the eigenvalues of the model atomic Hamiltonian

H model = T + V model , (2)
where T is the kinetic energy operator, reproduce the experimental energies of the Potassium atom.

The description of the energies in the atomic limit, is indeed crucial, and has already been adequately discussed in much of the earlier work on model potentials notably by McCarroll and Valiron [31,32]. Other aspects of the model potential approach were treated by several authors [33,34,35] on the introduction of other non-local constraints, such as the orthogonality to the inner orbitals. However, in the present study, we mainly recall that the aim of the model potential is to represent in a physically meaningful way the ground and excited states with a reasonable precision. Since spin-orbit interactions q Z a b c 1 19 3.474 1.725 0.588333

Table 2: Klapisch model potential parameters [22]. a,b and c are given in units of inverse length (1/au).

e - t i E B( K + ) A( K + ) x x r A r B R Figure 1:
Relative coordinates for the active electron.

have been neglected (Hamiltonian in Eq (3) does not contain any spin dependence), the precision of our computed energies in the asymptotic limit and which are shown in Table 1, is adequate for the particular applications of this work. Fixing the two Potassium cores at A and B positions respectively as shown in Fig. 1, the K + 2 molecular cation is then treated as a set of two closed shell K + (3p 6 ) cores and a single active electron moving in the field created by the two ionic cores. The electronic Hamiltonian is given by

H e = T + V model (r A ) + V model (r B ) + U core , (3)
where T denotes the electronic kinetic energy operator, r A and r B the radial distances of the active electron from the cores K + respectively. U core is mainly The nuclei repulsion U core term is modified by a polarization term, namely:

U core = 1 R - α d R 4 - α q R 6 , (4)
where α d and α q are the dipole and quadrupolar polarizabilities, respectively, which have been used in the work of Magnier and Frécon [2], α d = 5.354 a.u [25] and α q = 4.69 a.u [START_REF] Baskkin | Atomic energy levels and Grotrian diagrams[END_REF].

The spectrum of the effective Hamiltonian H e defined in Eq (3) is determined by conventional variational techniques, using a basis set of Slater-type orbitals ψ sto , expressed in prolate spheroidal coordinates

λ ∈ [1, ∞], µ ∈ [-1, 1] and , φ ∈ [0, 2π] where λ = 1 2R (r A + r B ), µ = 1 2R (r A -r B ), φ, (5)
is the azimuthal angle.

In this prolate spheroidal coordinates system, the basis functions of Slater-type orbitals take the form:

ψ sto (λ, µ, φ) = N [(λ 2 -1)(1 -µ 2 )] Λ 2 λ l µ k e -R 2 ζ(λ+ µ) e imφ , (6)
where N denotes a normalization factor, a parameter that takes the values ±1 depending on whether the orbital is centered on nucleus A or nucleus B, (l, k) integers, ζ a variational parameter describing the decay of the orbital, m the projection of the electronic angular momentum on the z-axis (coinciding with R) and Λ =| m |=0,1,2,. . . for Σ + , Π and ∆ molecular states respectively.

Results

3.1. Σ + states
We generated the potential energy curves of the 2 Σ + states dissociating up to the K + +K(6s) limit. Corresponding graphics are plotted in Fig. 2, while the spectroscopic constants are listed and compared to available values issued from different approaches in Table 3.

For the 1 2 Σ + g ground state correlated to K + +K(4s) in the asymptotic limit, present work yields a dissociation energy D e =6630 cm -1 , an equilibrium distance R e =8.33 a.u, vibration constants ω e =73.04 cm -1 and ω e χ e = 0.201 cm -1 . Those values may be compared to the ab-intio results [1], D e =6454 cm -1 , R e =8.58 a.u, ω e =73.05 cm -1 and ω e χ e = 0.14 cm -1 . On the other hand, as indicated in Table 3, it is interesting to note that our dissociation energy compares well to the experimental value D e =6670 cm -1 , obtained by Stwalley et al. [12].

For the first excited state 1 2 Σ + u also correlated to K + +K(4s) in the asymptotic limit, the analysis of our potential energy curve indicates the existence of an energy minimum occurring at an equilibrium distance R e = 22.18 a.u, with a corresponding dissociation energy D e =80 cm -1 . As shown in Table 3, these values compare well to the ones issued from different approaches [1,2,18]. Our vibration constants ω e =5.72 cm -1 and ω e χ e =0.102 cm -1 are in perfect agreement with those based on ab-initio method [1], ω e =5.72 cm -1 and ω e χ e =0.13 cm -1 .

A good agreement is also observed between our 2 2 Σ + g state spectroscopic constants and those resulting from previous studies [1,2]. Present study yields a dissociation energy D e =3070 cm -1 and an internuclear equilibrium distance R e =15.96 a.u. These values may be compared to the ab-initio ones [1], D e =3038 cm -1 , R e =16.02 a.u and to the ones resulting from a model potential approach [2], D e =3048 cm -1 and R e =16.00 a.u.

For the 2 2 Σ + u state correlated to K + +K(4p) in the asymptotic limit, an energy minimum is observed at R e = 29.33 a.u, with a corresponding dissociation energy D e =199 cm -1 . Our vibration constants ω e =5.83 cm -1 and ω e χ e =0.043 cm -1 compare favorably with the ab-initio results [1], ω e =5.83 cm -1 and ω e χ e =0.05 cm -1 .

The 3 2 Σ + u and 3 2 Σ + g states are correlated to K + +K(5s) in the asymptotic limit. As predicted by previous studies, present work shows that both states are attractive with spectroscopic constants listed in Table 3, that compare favorably with values resulting from other different theoretical approaches [1,2].

For both states 4 2 Σ + g and 4 2 Σ + u which are correlated to K + +K(3d) in the asymptotic limit, we have found an equilibrium distances R e = 19.04 a.u and R e =43.9 a.u respectively. A good agreement is also observed between our 4 2 Σ + u and 4 2 Σ + g spectroscopic constants and those based on an ab-initio approach [1].

The potential energy curve of our 5 2 Σ + g state correlated to K + +K(5s) in the asymptotic limit, exhibits an energy minimum at R e =39.74 a.u with a dissociation energy D e =1211 cm -1 . As shown in Table 3, spectroscopic constants obtained in the present work are in good agreement with values resulting from other different theoretical approaches.

According to the ab-initio studies of Berriche [1], the 5 2 Σ + u state is attractive, whereas Magnier and Frécon [2] employing model potential techniques described a such state as completely repulsive. In the present study, our potential energy curve has an energy minimum occurring at an equilibrium internuclear distance R e =56.04 a.u, with a corresponding dissociation energy D e =169 cm -1 . Those may be compared to the ab-initio values [1], R e =56.35 a.u and D e =232 cm -1 . Equally, our vibration constants ω e =2.91 cm -1 and ω e χ e =0.012 cm -1 compare well to the ab-initio results [1], ω e =3.04 cm -1 and ω e χ e =0.01 cm -1 .

Another state which has been the subject of disagreement concerns the 6 2 Σ + g state correlated to K + +K(4d) in the dissociation limit. The analysis of our potential energy curve shows the existence of two minima. Therefore, contrary to Magnier and Frécon claims, as predicted by other ab-initio methods [1], our model potential approach confirms the attractive character of the 6 2 Σ + g state. As indicated in Table 3, our 6 2 Σ + g state has an energy minimum in the inner region, occurring at R e =31.39 a.u, with a corresponding dissociation energy D e =1200 cm -1 , vibration constants ω e = 9.56 cm -1 and ω e χ e = 0.019 cm -1 . Morever, for the minimum located in the outer region, we obtained an equilibrium distance R e =52.48 a.u with a well depth D e =1272 cm -1 . Those to be compared to the ab-initio values R e =53.17 a.u and D e =1130 cm -1 .

For the states 7 2 Σ + g and 7 2 Σ + u contrary to Magnier and Frécon [2] who concluded that such states are repulsive, our model potential results clearly indicate that 7 2 Σ + g and 7 2 Σ + u are attractive as predicted by ab into approach [1]. We found an equilibrium distance R e = 42.74 a.u with an dissociation energy D e = 562 a.u for the 7 2 Σ + g state, while for the 7 2 Σ + u state our equilibrium distance and dissociation energy are respectively R e = 67.74 a.u and D e = 26 a.u.

Π states

Calculations have also been carried out for the first eight 2 Π states for internuclear separations up to 80 a.u. Corresponding graphics are plotted in Fig. 3, while the spectroscopic constants are listed and compared to available values issued from different approaches in Table 4.

Our potential energy curve for the 1 2 Π g state does not exhibit any minimum. Therefore, as predicted by previous studies based on different methods, our 1 2 Π g state correlated to K + +K(4p) in the dissociation limit, is purely repulsive.

On the other hand, for the 1 2 Π u state also correlated to K + +K(4p) in the dissociation limit, present study yields a dissociation energy D e =2652 cm -1 and an equilibrium distance R e = 9.58 a.u. Those may be compared to the ab-initio results [1], D e =2367 cm -1 and R e = 9.96 a.u and to those employing model potential approach [2], D e =2682 cm -1 and R e = 9.5 a.u. 8 The 2 2 Π g and 2 2 Π u states are correlated to K + +K(3d) in the asymptotic limit. In the present study, as indicated in Table 4, both states have been found to be attractive with spectroscopic constants which compare well with other values issued from previous works .

Another disagreement between the ab-initio results [1] and those based on model potential approach [2] concerns the 3 2 Π g state, and which has been sorted out in the present study. Contrary to the conclusion of Magnier and Frécon [2], present results agree with the ab-initio ones. Our potential energy curve has an energy minimum occurring at R e = 58.03 a.u with a very weak dissociation energy D e =2 cm -1 , while the ab-initio values [1] are R e = 57.92 a.u. and D e =5 cm -1 .

For the 3 2 Π u state, present work yields a dissociation energy D e =691 cm -1 and an equilibrium distance R e =32.00 a.u. These values may be compared to the ab-initio ones [1], D e =615 cm -1 and R e =32.22 a.u, and to the previous ones based on model potential approach [2], D e =681 cm -1 and R e =32 a.u. As indicated in Table 4, our vibration constants ω e =9.42 cm -1 and ω e χ e = 0.032 cm -1 , are slightly different from the ab-initio values [1], ω e =8.32 cm -1 and ω e χ e = 0.01 cm -1 The 4 2 Π g correlated to K + +K(3d) in the asymptotic region, has also been a subject of disagreement between ab-initio results [1] and model potential ones [2]. Present investigation clearly indicates that 4 2 Π g potential energy curve has an energy minimum occurring at R e = 59.34 a.u, with a corresponding dissociation energy D e = 210 cm -1 . These values are in agreement with ab-initio results [1], R e =60.69 a.u and D e =198 cm -1 .

As indicated in Table 4, our equilibrium distance R e = 41.75 a.u and dissociation energy D e =1202 cm -1 for the 4 2 Π u state, compare well to previous results particularly to the ones by Magnier and Frécon [2] giving, D e =1216 cm -1 and an equilibrium distance R e =41.75 a.u. Present study also gives vibration constants ω e =6.96 cm -1 and ω e χ e = 0.010 cm -1 , which compare favorably with the ab-initio values [1], ω e =7.4 cm -1 and ω e χ e =0.01 cm -1 .

∆ states

We generated the adiabatic energies of the 2 ∆ states for internuclear distances up to 80 a.u. Our energy curves are plotted in Fig. 4, while the corresponding spectroscopic constants are listed and compared to available data resulting from different approaches in Table 5.

Our potential energy curve relative to the 1 2 ∆ g state exhibits an energy minimum occurring at an internuclear distance R e = 16.15 a.u, with a corresponding dissociation energy D e = 164 cm -1 . While these values compare favorably to those using a model potential approach of Magnier and Frécon [2] and Henriet [21], a difference of 140 cm -1 is observed when comparison is made with ab-initio results [1] which seem to underestimate the dissociation energy value.

On the other hand, while our 1 2 ∆ u state is repulsive as predicted by previous studies, our potential energy curve for the 2 2 ∆ g state, exhibits an energy minimum occurring at an equilibrium internuclear distance R e =28.79 a.u, with a corresponding dissociation energy D e =1189 cm -1 . Those to be compared to the ab-initio results [1], R e = 30.65 a.u and D e =698 cm -1 .

As previously mentioned, the 2 2 ∆ u state correlated in the asymptotic limit to K + +K(4d) has been the subject of a disagreement, when comparing ab-initio results [1] and previous model potential [2] ones. Contrary to the previous model potential results [2], our study clearly indicates that the 2 2 ∆ u state is attractive with a dissociation energy D e =67 cm -1 and an internuclear equilibrium distance R e =53.50 a.u. These values compare favorably to the ab-initio ones [1], R e =55.82 a.u and D e =40 cm -1 .

Avoided crossings

The obtained potential energy curves exhibit several avoided crossings between adjacent adiabatic states which belong to the same molecular symmetry. In Table 6, we give the internuclear distance R c of the pseudo-crossing and the corresponding energy gap ∆E c . Such quantities are required to evaluate cross sections of charge transfer process when Landau-Zener [30] approximation is used. Except the avoided crossing 2 2 Π g -4 2 Π g , our internuclear distances R c of the pseudo-crossing compare favourably to the ones issued from an ab-initio approach [1].

Conclusion

We have investigated the K + 2 alkali-dimer using a model potential approach to compute the adiabatic potential energies and the molecular spectroscopic constants for the first low-lying states (fourteen 2 Σ + , eight 2 Π and four 2 ∆), dissociating up to the K + +K(6s) limit.

Present results compare well to those based on ab initio [1] methods, which is a clear indication that model potential approach remains an efficient way of providing an accurate description of molecular systems involving a single active electron.

On the other hand, it is interesting to underline that previous existing disagreement between ab-initio results [1] and those based on a model potential approach [2], concerning the 5 2 Σ + u , 6 2 Σ + g , 6 2 Σ + u , 7 2 Σ + g , 7 2 Σ + u , 3 2 Π g , 4 2 Π g and 2 2 ∆ u states, has been sorted out in the present work. Contrary to the conclusion of Magnier and Frécon [2], our model potential results show that the potential energy curves relative to these states are attractive with a well pronounced energy minimum. This is in a perfect agreement with results predicted by studies exploiting ab-initio techniques.

Figure 2 :

 2 Figure 2: K +2 potential energy curves given with respect to the ground state dissociation limit. Σ + g in solid lines, Σ + u in dashed lines.

Figure 3 :

 3 Figure 3: K + 2 potential energy curves given with respect to the ground state dissociation limit. Πu in solid lines, Πg in dashed lines.

Figure 4 :

 4 Figure 4: K +2 potential energy curves given with respect to the ground state dissociation limit. ∆g in solid line, ∆u in dashed lines.

Table 1 :

 1 Comparison of calculated and experimental energy levels in cm -1 of different states in the dissociation limit.

	Level	Present work	[2]	[1]	Exp [28]
	K + + K(4s)	0.00	0.00	0.00	0.00
	K + + K(4p)	13092.54	13077.61 13022.99 13042.90
	K + + K(5s)	21010.53	20980.67 21018.02 21026.55
	K + + K(3d)	21665.98	21575.01 21535.54 21536.99
	K + + K(5p)	24712.41	24670.70 24750.85 24720.14
	K + + K(4d)	27438.61	27366.94 27417.26 27398.15
	K + + K(6s)	27461.49	27424.01 27451.93 27450.71
	description of the K + 2 system without difficulties. In the present work we re-investigate the K +

Table 3 :

 3 Spectroscopic constants for 2 Σ + states, R e in a.u, D e , ω e , ω e χ e , B e and T e in cm -1 .

	State	R e	D e	T e	ω e	ω e χ e	B e	Reference
	1 2 Σ + g	8.33	6630	0	73.04	0.201	0.0444	Present work
		8.58	6454	0	73.05	0.14	0.041818	[1]
		8.3	6685	0	73.2			[2]
					5

Table 4 :

 4 Spectroscopic constants for 2 Π states, R e in a.u, D e , ω e , ω e χ e , B e and T e in cm -1 .

	State	R e	D e	T e	ω e	ω e χ e	B e	Reference
	1 2 Π g Repulsive						Present work
		Repulsive						[1]
		Repulsive						[2]
		Repulsive						[18]
	1 2 Π u	9.58	2652 17072 36.42 0.125	0.0336	Present work
		9.96	2367 17113 35.16 0.29 0.031017	[1]
		9.5	2682 17066 37.86			[2]
		10.41	2218 16838 31.76 0.00		[16]
		9.552	2463 17089 40.1	0.09		[18]
	2 2 Π g	37.92	89	28154 2.56 0.019	0.0021	Present work
		38.47	81	27910 2.94	0.03 0.002081	[1]
		38.00	88					[2]
	2 2 Π u	20.68	1511 26730 14.28 0.034	0.0072	Present work
		20.9	1376 26619 16.11 0.03 0.007047	[1]
		20.75	1504					[2]
	3 2 Π g	58.03	2	31346 0.81 0.084	0.0009	Present work
		57.92	5	31206 1.11	0.04 0.000918	[1]
		Repulsive						[2]
	3 2 Π u	32.00	691 30654 9.42 0.032	0.0030	Present work
		32.22	615 30595 8.32	0.01 0.002967	[1]
		32	681					[2]
	4 2 Π g	59.34	210 33809 2.44 0.007	0.0009	Present work
		60.69	198 33674 2.76	0.01 0.000836	[1]
		Repulsive						[2]
	4 2 Π u	41.75	1202 32817 6.96 0.010	0.0018	Present work
		42.37	983 32889	7.4	0.01 0.001716	[1]
		41.75	1216					[2]

Table 5 :

 5 Spectroscopic constants for 2 ∆ states, R e in a.u, D e , ω e , ω e χ e , B e and T e in cm -1 .

	State	R e	D e	T e	ω e	ω e χ e	B e	Reference
	1 2 ∆ g	16.15	164 28051 10.90 0.181	0.0118	Present work
		17.92	25	27971 9.43	0.1	0.009593	[1]
		16.5	162					[2]
		16	274					[21]
	1 2 ∆ u Repulsive						Present work
		Repulsive						[1]
		Repulsive						[2]
	2 2 ∆ g	28.79	1189 32859 9.38 0.019	0.0037	Present work
		30.65	698 33174 7.83	0.01 0.003278	[1]
		28.75	1174					[2]
	2 2 ∆ u	53.50	67	33968 3.58 0.047	0.0011	Present work
		55.82	40	33832 1.77	0.01 0.000988	[1]
		Repulsive						[2]
					11

Table 6 :

 6 Avoided crossing parameters R c and ∆E c

	Crossing	R c (a.u) ∆E c (a.u)	Reference
	2 2 Σ + g -3 2 Σ + g	5.71	5.57 10 -3 Present work
		5.65		[1]
	3 2 Σ + g -4 2 Σ + g	18.28	1.74 10 -3 Present work
		18.19		[1]
	5 2 Σ + g -6 2 Σ + g	10.28	1.16 10 -2 Present work
		10.38		[1]
	6 2 Σ + g -7 2 Σ + g	43.02	5.05 10 -3 Present work
		40.91		[1]
	3 2 Σ + u -4 2 Σ + u	29.58	1.79 10 -3 Present work
		30.07		[1]
	2 2 Π g -4 2 Π g	9.90	1.99 10 -2 Present work
		6.74		[1]


Raphaël Pile

Guillaume Parent

Emile Devillers

Thomas Henneron

Yvonnick Le Menach

Jean Le Besnerais

Jean-Philippe Lecointe

Application Limits of the Airgap Maxwell Tensor

Keywords: Maxwell tensor, Magnetic stress, Electrical machines, Magneto-mechanical, Vibroacoustic

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.

I. INTRODUCTION

M AXWELL tensor (MT) is used to accurately compute overall forces on a given body. In particular, electrical machines designers are used to apply MT technique in the middle of the airgap to calculate the electromagnetic torque experienced by the rotor. This technique has been naturally extended to the study of magnetic force waves experienced by the outer structure -generally the stator -for vibroacoustic analysis. The method is generally based on a circular path in the airgap to ease physical interpretations with the Fourier transform [START_REF] Gieras | Noise of polyphase electric motors[END_REF][START_REF] Xu | Review of electromagnetic vibration in electrical machines[END_REF]. In that case, one may look for a circular path as close as possible from the stator bore radius to accurately estimate forces experienced by the stator. Indeed, for magnetic linear media [START_REF] Henrotte | Computation of electromagnetic force densities: Maxwell stress tensor vs. virtual work principle[END_REF], the theoretical magnetic surface force density at the air-ferromagnetic interface is exactly given by MT expression.

However, the radius of the circular path has a strong influence on calculated stress harmonics [START_REF] Mizia | Finite element force calculation: comparison of methods for electric machines[END_REF]- [START_REF] Hallal | Study of electromagnetic forces computation methods for machine vibration estimation[END_REF]. In Finite Element Analysis (FEA) simulations, the mesh is generally the finest in the airgap so that it is recommended to evaluate forces in the middle of the airgap. Besides, 2D FEA based on scalar potential cannot ensure high accuracy of both radial and tangential flux density at the air-ferromagnetic interface, which are both involved in MT expression. Then the numerical application of the MT at the interface is source of numerical errors [START_REF] Ren | Comparison of Different Force Calculation Methods in 3D Finite Element Modelling[END_REF].

Therefore, it is interesting to know how to transfer surface forces based on MT calculated in the airgap to another radius, in particular at tooth tip radius. Such a transformation would also be useful when using analytical electromagnetic models such as permeance magneto-motive force models, which focus on the accurate calculation of the flux density distribution at the airgap radius.

Therefore, this article compares the deviation between the application of airgap surface force based on MT and the surface forces at stator bore radius. Due to numerous numerical artifacts in practical electrical machines, only a slotless stator and rotor with constant permeability is studied in this article. It allows to obtain the analytic expression of the airgap potential and resulting surface forces at any radius. Thus the surface force transfer deviation is quantified by analytically-derived transfer coefficients (20) and (24) which depend on airgap geometry and magnetic surface force wavenumber.

II. PROBLEM DEFINITION

A. Maxwell Tensor

According to [START_REF] Henrotte | Computation of electromagnetic force densities: Maxwell stress tensor vs. virtual work principle[END_REF][START_REF] Rosensweig | Ferrohydrodynamics[END_REF], for an incompressible linearly magnetizable media the magnetic flux density B is related to the magnetic field H ∀x ∈ Ω by B(x) = µ(x, B)H(x) such that the magnetic stress tensor reduces to the following form:

T m = -µ H • H 2 I + µHH (1)
with I the identity tensor. In particular, the Vacuum Maxwell stress tensor is obtained for µ = µ 0 . The theoretic application of this magnetic stress tensor for the computation of magnetic surface force density P on a contour which exactly corresponds to the air-ferromagnetic interface (defined by its local normal vector n) leads to [START_REF] Henrotte | Computation of electromagnetic force densities: Maxwell stress tensor vs. virtual work principle[END_REF][START_REF] Ren | Comparison of Different Force Calculation Methods in 3D Finite Element Modelling[END_REF]:

P = 1 2 1 µ 0 - 1 µ B 2 n - µ 0 -µ 2 H 2 t n (2)
where B n and H t are the magnetic flux density locally normal (resp. magnetic field locally tangential) to the interface.

B. Airgap Surface Force

One of the most used method to compute magnetic surface forces based on MT is to use Stoke's theorem [START_REF] Landau | The classical theory of fields[END_REF] along a Rotor Stator Γ R s R ag Fig. 1. Slotless electrical machine used to compare magnetic surface force on the stator and in the airgap circular closed boundary Γ in the airgap at a given radius R ag as illustrated in Fig 1 . The result gives the total magnetic force F m acting on either the stator or rotor [START_REF] Gieras | Noise of polyphase electric motors[END_REF][START_REF] Ren | Comparison of Different Force Calculation Methods in 3D Finite Element Modelling[END_REF]:

F m = Γ - µ 0 2 H 2 n + µ 0 H n H dΓ (3)
Then the airgap MT approximation is to assume the term under the integral to be the magnetic surface force density in the polar referential:

P r (R ag ) ≈ 1 2µ 0 B r (R ag ) 2 - µ 0 2 H θ (R ag ) 2 P θ (R ag ) ≈ B r (R ag )H θ (R ag ) (4)
with P r the radial magnetic surface force density, P θ the tangential magnetic surface force density, B r the radial magnetic flux density, and H θ the tangential magnetic field.

Neglecting the ferromagnetic permeability µ contribution into (2) gives an expression very similar to (4) justifying the use of airgap surface force approximation. The application of (2) at the air-ferromagnetic interface remains source of numerical errors [START_REF] Ren | Comparison of Different Force Calculation Methods in 3D Finite Element Modelling[END_REF]: in the case of numerical model based on FEA, the continuity of both normal magnetic flux B n density and tangential magnetic field H t cannot be ensured on the interface between materials. As a consequence, publications often use airgap approximation (4) for the vibroacoustic design of electrical machines because the middle of the airgap is better meshed, it is less sensitive to variational formulation, and it allows the physical understanding of each magnetic force wave (for instance with analytical models) [START_REF] Gieras | Noise of polyphase electric motors[END_REF][START_REF] Xu | Review of electromagnetic vibration in electrical machines[END_REF].

C. Studied Slotless Case

Understanding the sources of airgap surface force (4) variations is a difficult task because of numerous artifacts such as slotting effect, sharp geometries, interference between the wavenumbers, etc. To avoid these main artifacts, a specific case is defined:

• Slotless stator and rotor to avoid sharp geometries and slotting effect. • Single-wave excitation to avoid interference between different wavenumbers. • Analytical solving to avoid meshing and numerical errors. The case presented in Fig. 1 proposes to fulfill all these constraints. The stator bore radius is R s and the middle airgap radius is R ag . In order to have only one magnetic wavenumber, the magnetic potential z-component A z is imposed at radius R ag such that ∀θ ∈ [0, 2π]:

A z (R ag , θ) = β n sin(nθ + ψ n ), n ∈ N * (5)
The continuity of tangential magnetic field can be expressed with the magnetic potential at the interface between the air (permeability µ 0 ) and ferromagnetic media (permeability µ) ∀θ ∈ [0, 2π]:

1 µ 0 ∂A z ∂r (R - s , θ) = 1 µ ∂A z ∂r (R + s , θ) (6)
Thus the boundary condition at r = R s with µ r → ∞ becomes:

∂A z ∂r (R s , θ) = 0 (7)
For these boundary conditions, Poisson's equation is solved for the 2D magnetic vector potential in polar coordinates ∀θ ∈ [0, 2π] and ∀r ∈ [R ag , R s ]:

1 r ∂ ∂r 1 r ∂A z ∂r + ∂ 2 A z ∂θ 2 = 0 (8)
This partial derivative problem can be solved analytically [START_REF] Lubin | Exact analytical method for magnetic field computation in the air gap of cylindrical electrical machines considering slotting effects[END_REF]. Thus the magnetic potential, flux density and field can be analytically computed in the studied domain. The next section provides the calculation steps which lead to the analytical expression of flux density (15) in the slotless case.

D. Analytic Magnetic Solving in Slotless Case

A solution exists and is unique for the previous system composed of [START_REF] Landau | The classical theory of fields[END_REF], [START_REF] Rosensweig | Ferrohydrodynamics[END_REF] and (5) according to [START_REF] Landau | The classical theory of fields[END_REF]. Then a method consists in stating a function A z and to check if it fulfills the boundary conditions: a solution similar to [START_REF] Lubin | Exact analytical method for magnetic field computation in the air gap of cylindrical electrical machines considering slotting effects[END_REF] is searched ∀θ ∈ [0, 2π] and ∀r ∈ [R ag , R s ]:

A z (r, θ) = γ n E n (r, R s ) + α n E n (R ag , r) E n (R ag , R s ) sin(nθ + ψ n ) (9)
with E n a polynomial function defined by:

E n : (x, y) ∈ R 2 → (x/y) n -(y/x) n (10)
Under this form, the vector A z is satisfying Poisson's equation [START_REF] Landau | The classical theory of fields[END_REF]. Then the coefficient γ n , α n , have to be determined in order to satisfy the boundary conditions. If a correct set of these coefficients is found, then A z would be the unique solution of the problem. First, satisfying the boundary condition [START_REF] Ren | Comparison of Different Force Calculation Methods in 3D Finite Element Modelling[END_REF] in the airgap leads to:

γ n = β n (11)
The second boundary condition [START_REF] Rosensweig | Ferrohydrodynamics[END_REF] leads to:

α n = 2β n /F n (R ag , R s ) (12)
with F n a polynomial function defined by:

F n : (x, y) ∈ R 2 → (x/y) n + (y/x) n (13)
Thus, the unique solution of the system is entirely defined with the geometrical parameters and excitation's wavenumber. The flux density is derived from the magnetic potential:

B = curl(A) = 1 r ∂A z ∂θ e r - ∂A z ∂r e θ = B r e r + B θ e θ (14)
With ( 14) and ( 9) both radial and tangential magnetic flux can be analytically expressed for this slotless machine:

B r (r, θ) = nβ n r E n (r, R s ) + 2En(Rag,r) Fn(Rag,Rs) E n (R ag , R s ) cos(nθ + ψ n ) B θ (r, θ) = - nβ n r F n (r, R s ) + 2Fn(Rag,r) Fn(Rag,Rs) E n (R ag , R s ) sin(nθ + ψ n ) (15)

III. TRANSFER OF SURFACE FORCES IN SINGLE MAGNETIC WAVE CASE

This section aims to compare the airgap surface forces (4) calculated at radius R ag with the exact value calculated at R s . Then an airgap surface radial force function is defined:

P r : [R ag , R s ] × [0, 2π] → R (r, θ) → 1 2µ0 B 2 r (r, θ) -B 2 θ (r, θ) (16
) In order to better understand the airgap transfer of the forces, a Fourier decomposition is performed:

P r (r, θ) = P r,0 (r) + P r,2n (r) cos(2nθ + 2ψ n ) (17)
The previous result (15) allows to express the magnetic flux under the form:

B r (r, θ) = B r,n (r) cos(nθ + ψ n ) B θ (r, θ) = B θ,n (r) sin(nθ + ψ n ) (18)
Then the airgap radial surface force can be decomposed under the form:

P r,0 (r) = B 2 r,n (r) -B 2 θ,n (r) /2 P r,2n (r) = B 2 r,n (r) + B 2 θ,n (r) /2 (19)
This is a classical result: a magnetic wavenumber n in the airgap is recomposed into magnetic force wavenumbers 0 and 2n. For the studied slotless case ∀θ ∈ [0, 2π], B θ (R s , θ) = 0.

In order to avoid θ dependency, the amplitude of each harmonics are compared independently. To this purpose, a function is defined ∀r ∈ [R ag , R s ], for k ∈ N:

R k (r) = P r,k (r)/P r,k (R s ) (20)
The next step is to introduce the analytic solution of the magnetic flux (15) into (20) in order to get the analytic expression of each function R 0 and R 2n . After some calculations, the expression of these two functions simplify remarkably well as follows for k = 0 or k = 2n:

R k (r) = (R s /r) 2 F k (r, R s )/2 (21)
Note that ∀r, F 0 (r, R s ) = 1. Then combining (17) with ( 21), the relation between the airgap radial surface force and the exact radial surface force at position r = R s can be established:

P r (R s , θ) = r 2 R 2 s P r,0 (r) + 2P r,2n (r) F 2n (r, R s ) cos(2nθ + 2ψ n ) (22)
The differences between exact radial surface force density and airgap radial surface force depends on the airgap width, the radius of application and the wavenumber, but not on the magnetic excitation β n . If the airgap surface forces are equivalent to the exact surface forces at R s , then both R 0 0.9 0.92 0.94 0.96 0.98 1 0.8 0.9

1

Adimensional radius r/Rs R 2n (r) n = 1 n = 2 n = 4 n = 8
Fig. 2. Comparing stator and airgap Maxwell Tensor magnetic surface force for the harmonic 2n

and R 2n should be very close to 1 for every wavenumbers and for any radius. Fig. 2 shows that these coefficients are not always negligible and they can strongly affect high magnetic surface force wavenumbers. Nevertheless, the vibro-acoustic of electrical machine generally considers only low wavenumbers such that the airgap surface force stays relatively accurate in the radial direction with narrow airgaps. A similar method is applied to estimate the impact of the airgap transfer on the tangential surface force, with the following function:

P θ : [R ag , R s ] × [0, 2π] → R (r, θ) → 1 µ0 B r (r, θ)B θ (r, θ) (23)
Then the next result can be easily demonstrated following the same methodology:

P θ (r, θ) = - R 2 s r 2 E 2n (r, R s ) 2 P r,2n (R s ) sin(2nθ + 2ψ n ) (24)
In order to estimate the relative error ∆, the airgap width is defined as g with g = 2 (R ag -R s ) R ag . Then a Taylor approximation on (24) leads to:

∆ ∝ 2ng/R ag (25)
However on the studied case, the tangential surface force density should be null. A new criterion can be defined: ∆ 1 means that airgap surface forces are accurate. Note that the resultant torque is constant for any radius -2π 0 P θ (r, θ)rdθ = 0 -but it is not true for the total radial force because it is varying linearly with the radius of application: P r (r, θ)dθ = Rs r P r (R s , θ)dθ. Nevertheless, the transfer coefficient does not change the unbalanced magnetic forces F x and F y in the Cartesian referential.

IV. GENERALIZATION TO MULTI-HARMONIC CASE

The previous section demonstrates a new transfer coefficient for single magnetic excitation wave. In this section, the coefficient (21) is validated to any combination of wavenumbers.

Pr is the Fourier transform of P r and if P (r, k) is the k th harmonic of the airgap surface force (4) :

P r (r, θ) = k P (r, k) cos(kθ + φ k ) (26)
From this airgap surface force, a transfer model is proposed by assuming that -despite the recombination of several flux

P r (R s , θ) = k 2 P (r, k) F k (r, R s ) cos(kθ + φ k ) (27)
The validity of the assumption is verified for the same slotless case as in Fig. 1 with infinite ferromagnetic permeability. Thus the flux density is entirely solved analytically using a linear combination of (15). The numerical inputs of the validation case are presented in Table I. The results are presented in Fig. 3: P A r denotes the application of ( 16) at r = R ag , P P r denotes the application of (27), and P S r denotes the exact value of magnetic surface forces obtained with (16) at r = R s . Comparing the L 2 -norm of the two signals with respect to the stator surface forces leads to:

P P r -P S r L 2 ≈ 10 -5 P A r -P S r L 2 (28)
A perfect match is observed between the P P r and the P S r , which means that the intuited formula ( 27) is correct at least for the slotless geometry. A similar result can be found with the tangential component.

This last result allows the partial use of these coefficients for slotted electrical machine: all the previous results hold when the angular domain is restricted. Thus supposing the previous reasoning is applied on the angular opening of an electrical machine stator's tooth tip -with infinite ferromagnetic permeability -then the transfer method (27) can be applied but only for this restricted angular opening. Nevertheless, the transfer of the airgap surface force based on MT in front of the slotted area might not be properly predicted with these coefficients.

V. CONCLUSION

This paper objective was to understand to what extent the surface force density computed at the airgap radius with Maxwell Tensor differs from the exact magnetic force density. For this purpose, a simplified slotless academic machine was defined in order to find analytical surface force transfer coefficients. The study of a single flux density excitation demonstrates the existence of a a new transfer coefficient. The final part of the paper shows that these transfer coefficients remain unchanged with a multi-harmonic magnetic excitation.

The importance of these transfer coefficients depends on the airgap size and force wavenumbers to be taken into account for vibroacoustic analysis. Only these two factors were identified. The transfer coefficients have little influence on low magnetic force wavenumbers but significant influence on high force wavenumbers. Consequently, these coefficients can be used to estimate the deviation between the surface force calculated in the air and the surface force calculated on the stator bore radius. Considering infinite permeability, the transfer coefficients could be applied to correct the radial magnetic surface forces from the airgap Maxwell Tensor with a real electrical machine, but only in front of a tooth tip.

Future research work will address the generalization of these coefficients to slotted areas of electrical machines and to tangential surface forces. The authors expect to use this slotless machine as a reference case in future development of magneto-mechanical coupling methods: the simplified geometry allows to solve analytically the magnetic and mechanical equations such that it could be compared to numerical and experimental results in future research work.

TABLE I INJECTED

 I MAGNETIC POTENTIAL SPECTRUM FOR THE VALIDATION CASE WITH RAG = 42.5 mm AND RS = 45.0 mm.

	Wavenumber n	5	7	15		17	19
	Amplitude βn [T.m] 5.2615 0.6861 0.2119 0.2725 0.1926
	Phase ψn [rad]	3.8893 0.7100 0.4700 0.6849 0.5410
		•10 5
		0
	]
	Surface force [N/m 2	-1 -2	P A r P P r P S r
		0	0.5	1	1.5	2	2.5	3
				Angular position [rad]
		•10 4
	] Surface force [N/m 2	3.7 8.6					P A r (Airgap) P P r (New method) P S r (Reference)
		0 2 4	8 10 12 14	20 22 24 26	30 32 34 36 38
				Angular wavenumber
	Fig. 3. Comparison of several methods for surface force density using
	Maxwell Tensor
	density harmonics into one force harmonic -the coefficient
	(21) stays valid:


domain: math info

text:
Qinmeng Zou

Frédéric Magoulès
email: frederic.magoules@hotmail.com

Convergence Detection of Asynchronous Iterations based on Modified Recursive Doubling

Keywords: recursive doubling, asynchronous iterations, convergence detection, message passing interface

This paper addresses the distributed convergence detection problem in asynchronous iterations. A modified recursive doubling algorithm is investigated in order to adapt to the non-power-of-two case. Some convergence detection algorithms are illustrated based on the reduction operation. Finally, a concluding discussion about the implementation and the applicability is presented.

Introduction

Consider the following linear system Ax = b, where A ∈ R n×n and b ∈ R n . A splitting A = M -N, yields an iterative scheme

x k+1 = T x k + c,
where k ∈ N, T = M -1 N and c = M -1 b. Generally, this scheme is well suited for parallel computing of the form

x k+1 i = T i x k + c i , i ∈ {1, . . . , p},
where p is the number of processors. Here, x i and c i might be two values or two smaller vectors. Similarly, T i might be a row vector or a smaller matrix. They are distributed in different processors. However, a specific point is required at the end of each iteration to synchronize between the processors. The waste of time may be significant in the case of unbalanced working load and node failure, which gives rise to the asynchronous iterative methods. The asynchronous iterative scheme has been proposed by Chazan and Miranker [START_REF] Chazan | Chaotic relaxation[END_REF] for the solution of linear equations and generalized by several researchers (see, e.g., [START_REF] Miellou | Algorithmes de relaxation chaotique à retards[END_REF][START_REF] Baudet | Asynchronous iterative methods for multiprocessors[END_REF][START_REF] Tarazi | Some convergence results for asynchronous algorithms[END_REF][START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]) for the general problem

x k+1 = f x k ,
where f is a fixed point mapping. The asynchronous iterative scheme is shown as follows

x k+1 i = f i x τ i,1,k 1 , . . . , x τ i,p,k p , i ∈ P k , x k i , i / ∈ P k ,
where τ i,j,k ≤ k is a sequence of iterations with retards for each element j in each processor i, and P k ⊂ {1, . . . , p} is a sequence of subsets of processor numbers. In this case, processors are not required to wait for receiving all messages and allowed to keep on their own pace. We often add the following conditions to better investigate the chaotic process

k ∈ N | i ∈ P k = +∞, ∀i ∈ {1, . . . , p}, lim k→+∞ τ i,j,k = +∞, ∀i, j ∈ {1, . . . , p},
where |.| is the cardinality of a set that measures the number of elements. It means that no processors should be abandoned forever and more and more recent values should be used. Asynchronous iterative algorithms must terminate after a finite number of iterations, as suggested in [START_REF] Savari | Finite termination of asynchronous iterative algorithms[END_REF]. Thus, a practical implementation involves a set of admissible solutions S, such that

x * ∈ S, where x * is a solution vector. We would like to find a vector x established by the components from each processor, and we have to evaluate x ∈ S. If true, then x * = x; otherwise, continue the computation, as well as the evaluation. Thus, the termination condition can be expressed by a residual evaluation

f (x) -x < , > 0,
where . is a norm, is a well-chosen threshold. x is given as an arbitrary combination of local components

x = x k1 1 , . . . , x kp p , k 1 , . . . , k p ∈ N.
The major problem of termination detection is how to collect x ki i and execute the evaluations. Recently, several developments for the asynchronous iterations have been proposed in different domains, such as domain decomposition methods [START_REF] Magoulès | Asynchronous optimized Schwarz methods with and without overlap[END_REF][START_REF] Magoulès | Asynchronous iterative sub-structuring methods[END_REF][START_REF] Magoulès | Asynchronous iterations of Parareal algorithm for option pricing models[END_REF], convergence detection methods [START_REF] Miellou | A new stopping criterion for linear perturbed asynchronous iterations[END_REF][START_REF] Magoulès | Distributed convergence detection based on global residual error under asynchronous iterations[END_REF], and programming libraries [START_REF] Magoulès | JACK: an asynchronous communication kernel library for iterative algorithms[END_REF][START_REF] Magoulès | JACK2: An MPI-based communication library with non-blocking synchronization for asynchronous iterations[END_REF]. In this paper, we propose a modified recursive doubling algorithm applied to the non-blocking collective communication, which is addressed in the next section. In Section 3, we present some convergence detection strategies based on our new method. Finally, further discussion is given in Section 4 about the implementation and performance.

Modified Recursive Doubling

Recall that p is the number of processors. We define a µ 0 ∈ N such that

p 0 = 2 µ0 ≤ p < 2 µ0+1 ,
where p 0 denotes a pivot. Note that for the parallel iterative methods, the Allreduce function is very useful because we need to collect residual values from different processors. For the asynchronous iterations, however, the collective operations should be not only efficient, but also performed in a non-blocking way.

Traditional recursive doubling (see, e.g., [START_REF] Thakur | Optimization of collective communication operations in MPICH[END_REF]) is one of the possible algorithms for the Allreduce function. It involves a power-of-two number of processors, which adapts only to some special situations. We consider a modified version of recursive doubling, in which a backward shift and a forward shift are required in the general case. The first step is sending data from the extra processors to the first several processors, called backward shift, illustrated in Figure 1. During this process, we proceed the corresponding arithmetical operations, such as summation, maximization, and minimization. Then, the recursive doubling algorithm is proceeded only within the power-of-two processors to exchange data and execute reduction operation as shown in Figure 2. Finally, a forward shift is proceeded to send back the final data to the extra processors as illustrated in Figure 3, which is indeed the inverse operation of the first shift.

Asynchronous iterations require non-blocking communication, which can be implemented is several ways. For example, we might prefer to create a new thread for a desired collective

P0 P1 P2 P3 P9 P8 P7 P6 P5 P4 Figure 1: Backward shift P0 P1 P2 P3 P9 P8 P7 P6 P5 P4 Figure 2: Recursive doubling P0 P1 P2 P3 P9 P8 P7 P6 P5 P4
Figure 3: Forward shift function, and then design its behavior by some external interface functions; we could also create a state-based interface that should be invoked in user applications, in which some lightweight functions act as different states in the life cycle of a collective operation. Here, we adopt the latter and give an example of state diagram depicted in Figure 4, which implements a non-blocking Allreduce function. From the picture, we can see that each cycle begins with the backward shift operation. If the rank of processor belongs to the extra range, it sends data and enters into the forward shift state. In other cases, the relative processors must enter into the updating loop, which involves all the processors having a rank smaller than the pivot p 0 . Finally, a forward shift process is executed to gather the final results to the non-power-of-two processors. Notice that the first several processors within the exponential area engage as well in the shift subroutines.

The amount of data exchanged by each processor depends on the way of collecting residual values. We need exactly log p 0 + 2 steps to finish a cycle in the synchronous case. If there is only a floating point residual value being exchanged in each processor, then totally p 0 log p 0 + 2(p -p 0 ) data are exchanged in each cycle. In asynchronous mode, this number keeps the same. However, processors wait no longer the others and conduct iterations on their own pace.

Convergence Detection Algorithms

We could also develop other collective operations based on the backward-forward recursive doubling algorithm. In practice, however, these functions are rarely used in the context of asynchronous iterations because we expect to exploit the most recent values as much as [ else ]

[ i >= p 0 ] / Send shift msg. when ( shift msg. received ) / Send recursive msg.

[ else ]

[ j < μ 0 ]

[ i < p -p 0 ] / Send shift msg.

[ else ] We consider first an inexact residual collection strategy that involves only the Allreduce function, depicted as follows.

res_norm = res_thresh res_loc = res_thresh res_glb = res_thresh flag = 1 while res_norm ≥ res_thresh do

z i = x i Compute(x i , A i , b i , x) Send(x i ) Allreduce(res_glb, res_loc, flag) if flag then res_norm = res_glb res_loc = x i -z i ∞ flag = 0 end if Recv(x) end while
We mention here that although such algorithm is not exact, it might be efficient due to the simplicity and still has an acceptable precision. In the algorithm, we take the maximum norm as an example to compute residual and omit some function parameters, e.g., the arithmetic operation of reduction. Unlike the message-passing standard [START_REF]Message Passing Interface Forum. MPI: A message passing interface standard[END_REF], our implementation is based on the state that requires the function invocation repeatedly, not just a request handler. The Compute function could be any appropriate iterative algorithm, such as Jacobi method or gradient method (see, e.g., [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]). This is inexact because res_loc might not be monotone all over the iterations. Sometimes global residual indicates a convergence signal but local residual rises instead due to the retard term.

Now we give a second version that leads to an exact solution in view of the residual collection, shown as follows.

res_norm = res_thresh sflag = 1 eflag = 0 while res_norm ≥ res_thresh do Compute(x i , A i , b i , x) Send(x i ) if sflag then Snapshot(x, x i , eflag) if eflag then z i = xi Compute(x i , A i , b i , x) res_loc = xi -z i ∞ sflag = 0 eflag = 0 end if else Allreduce(res_glb, res_loc, eflag) if eflag then res_norm = res_glb sflag = 1 eflag = 0 end if end if Recv(x)
end while This algorithm involves a distributed snapshot process that generates a consistent solution buffer [START_REF] Chandy | Distributed snapshots: Determining global states of distributed systems[END_REF] (see also, e.g., [START_REF] Savari | Finite termination of asynchronous iterative algorithms[END_REF][START_REF] Magoulès | Distributed convergence detection based on global residual error under asynchronous iterations[END_REF]). The snapshot algorithm first sends x i to the processors that depend on x i . In this situation, we call them dependent neighbors; then, processor i begins to wait for the necessary data from some other processors, which are called essential neighbors [START_REF] Bertsekas | Parallel and Distributed Computation: Numerical Methods[END_REF]; finally, it returns a collection of essential data that are used for the residual computation. Here we simplify the process by assuming that the communication follows an "all-to-all" pattern, which implies that both dependent neighbors and essential neighbors are all the other processors so that they are the same. For the general case, the algorithm would be similar. We first set sflag = 1 to enable snapshot process. Then, we compute res_loc when snapshot finishes and set all sflag = 0, which enables the reduction process. Finally, Allreduce is called repeatedly that is exactly the first algorithm, except that this time we keep a set of consistent data that provides an exact result.

Further Discussion

In this section we first discuss the implementation issue of the convergence detection algorithms. Here we take the message passing interface (MPI) standard as an example. Notice that in order to implement a non-blocking function, we should execute the relative instructions in an independent thread, which involves an explicit construction or an implicit In the next version, these functions will be redefined in order to provide a more flexible interface. Notice that if the number of processors falls on the power-of-two case, the iterations in Figure 4 will jump over all the shift steps appropriately. Such case has been proven very efficient in several situations [START_REF] Thakur | Optimization of collective communication operations in MPICH[END_REF], whereas our algorithm can benefit from it as well. On the other hand, our Send operation is implemented in a blocking mode because it causes rarely a negative impact in practice on the efficiency. We could avoid wasting time by switching it to non-blocking mode without changing so much codes.

Finally, we mention here that our algorithm is suitable for a relatively "close" distributed environment; otherwise, there might be a great deal of communication operations exchanging data between long-distance nodes, which increases the transfer time. In such case, a tree-based algorithm is preferred. However, asynchronous iterations may not exhibit advantages in a completely local cluster, even perform sequentially like synchronous scheme with much more ongoing messages. Consider a two-point boundary value problem with an asynchronous relaxation solver [START_REF] Chazan | Chaotic relaxation[END_REF]. We implement the mathematical operations by Alinea [START_REF] Magoulès | Alinea: An advanced linear algebra library for massively parallel computations on graphics processing units[END_REF] and the asynchronous iterations by JACK [START_REF] Magoulès | JACK2: An MPI-based communication library with non-blocking synchronization for asynchronous iterations[END_REF], which have been proven very efficient for the large-scale scientific computing [START_REF] Magoulès | Auto-tuned Krylov methods on cluster of graphics processing unit[END_REF][START_REF] Magoulès | Optimized Schwarz method without overlap for the gravitational potential equation on cluster of graphics processing unit[END_REF][START_REF] Magoulès | Green computing on graphics processing units[END_REF]. The finite difference scheme is adopted for the discretization. The matrix dimension n = 10000 and b is chosen arbitrarily from -10 to 10. Results are shown in Figure 5. We observe that the iteration curve shows synchronous behavior that exists a bottleneck within a specific range of processors. The experiment was performed on a cluster of Intel Xeon CPU E5-2670 v3, connected by FDR Infiniband network with 56 Gbit/s, which is concentrated and favors synchronous iterations. Furthermore, an "all-to-all" algorithm generates huge amounts of messages in the asynchronous mode, which makes the network too messy to be efficient. In this case, we prefer the traditional synchronous iterative scheme, even for large-scale parallel computing.

Figure 4 :

 4 Figure 4: Statechart of a non-blocking Allreduce function

Figure 5 :

 5 Figure 5: Asynchronous iterations in a concentrated environment

Acknowledgment

This work was supported by the French national programme LEFE/INSU and the project ADOM (Méthodes de décomposition de domaine asynchrones) of the French National Research Agency (ANR).


text:
Qinmeng Zou

Frédéric Magoulès
email: frederic.magoules@hotmail.com

A New Cyclic Gradient Method Adapted to Large-Scale Linear Systems

Keywords: gradient methods, linear systems, large-scale problems, Barzilai-Borwein methods

This paper proposes a new gradient method to solve the large-scale problems. Theoretical analysis shows that the new method has finite termination property for two dimensions and converges R-linearly for any dimensions. Experimental results illustrate first the issue of parallel implementation. Then, the solution of a large-scale problem shows that the new method is better than the others, even competitive with the conjugate gradient method.

Introduction

We are interested in investigating new gradient methods for the solution of linear system

Ax = b, (1)
where A ∈ R n×n is symmetric positive definite (SPD) and b ∈ R n . This problem is equivalent to the minimization of a convex quadratic function

f (x) = 1 2 x Ax -b x. (2)
Gradient methods generate a sequence of the form

x k+1 = x k -α k g k , k = 0, 1, . . . , (3)
where g k = Ax k -b. It is well known that the steepest descent (SD) method [START_REF] Cauchy | Méthode générale pour la résolution des systèmes d'équations simultanées[END_REF] performs poor in most cases, where the steplength can be written as follows

α SD k = g k g k g k Ag k . (4)
The iterates generated tend to asymptotically alternate between two directions [START_REF] Akaike | On a successive transformation of probability distribution and its application to the analysis of the optimum gradient method[END_REF]. In contrast, the conjugate gradient (CG) method [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] is often the method of choice that will terminate in at most n iterations. It is very attractive because of its high efficiency and low storage requirement. Nonetheless, CG iteration depends strongly on the search of direction calculation, i.e., any derivation such as round-off errors can seriously degrade performance [START_REF] Fletcher | On the Barzilai-Borwein method[END_REF].

In the past several decades, a renewed interest for gradient methods has appeared since Barzilai and Borwein [START_REF] Barzilai | Two-point step size gradient methods[END_REF] proposed two efficient nonmonotone steplengths

α BB1 k = g k-1 g k-1 g k-1 Ag k-1 . ( 5
) α BB2 k = g k-1 Ag k-1 g k-1 A 2 g k-1 . (6)
The motivation consists in approximating the Hessian and imposing some quasi-Newton properties. Some theories and experiments have shown that BB methods have good performance and are competitive with CG methods when low accuracy is required or small perturbation exists [START_REF] Fletcher | On the Barzilai-Borwein method[END_REF]. The convergence has been proven by Raydan [START_REF] Raydan | On the Barzilai and Borwein choice of steplength for the gradient method[END_REF]. Furthermore, Friedlander et al. [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF] provided a general framework under the name of "gradient method with retards" that SD and BB both belong to it, as well as several alternate methods proposed later [START_REF] Dai | Alternate step gradient method[END_REF][START_REF] Dai | Alternate minimization gradient method[END_REF][START_REF] Zhou | Gradient methods with adaptive step-sizes[END_REF]. Motivated by the two-dimensional finite termination property, Yuan [START_REF] Yuan | A new stepsize for the steepest descent method[END_REF] provided a somewhat complicated steplength

α Y k = 2 1 α SD k-1 -1 α SD k 2 + 4g k g k s k-1 s k-1 + 1 α SD k-1 + 1 α SD k , (7)
where

s k-1 = x k -x k-1 .
He gave two algorithms and some variants were investigated further by Dai and Yuan [START_REF] Dai | Analysis of monotone gradient methods[END_REF]. Among these methods, the second variant (DY) is the most efficient one according to the experiments in [START_REF] Dai | Analysis of monotone gradient methods[END_REF], where iterates are generated of the form

α DY k = α SD k , k mod 4 < 2, α Y k , otherwise. (8)
In this paper, we address the properties of cyclic gradient methods, especially their parallel behavior. We propose a new algorithm based on the Yuan steplength, which has also the two-dimensional finite termination property. In the next section, we introduce the cyclic gradient methods and propose our new steplength. In Section 3, we give the convergence results of the new method. Some numerical results are presented in Section 4. Finally, a concluding remark is shown in Section 5.

Cyclic Gradient Methods

Friedlander et al. [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF] proposed an ingenious framework that gives rise to a great number of potentially efficient algorithms. Firstly, assume that m ∈ N represents retard that allows to employ the information from previous iterations. Let

k = max{0, k -m}, (9)
then a collection of possible choices of steplength can be set as follows

α GMR k = g τ (k) A ρ(k) g τ (k) g τ (k) A ρ(k)+1 g τ (k) , (10)
where

τ (k) ∈ k, k + 1, . . . , k -1, k , (11)
and ρ(k) ∈ {q 1 , . . . , q m } , q j ≥ 0,

where k ∈ N. The next theorem summarizes the convergence result in [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF].

Theorem 1 [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF]. Consider the linear system (1) with A ∈ R n×n is SPD and b ∈ R n , where x * = A -1 b is the exact solution. Consider the gradient method (3) being used to solve (1) and the steplength α k given by [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF]. Then the sequence {x k } converges to x * starting from any point x 0 .

For a proof of the above theorem, see [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF]. Incidentally, several potential algorithms were provided therein, including the first cyclic gradient method under the name of cyclic steepest descent (CSD) as suggested in [START_REF] Dai | Alternate step gradient method[END_REF], which can be summarized as follows

α CSD k = α SD k , k mod m = 0, α k-1 , otherwise. ( 13
)
Notice that if we choose ρ(k) = 0 and τ (k) = k + 1, . . . , k -1, k, then (10) becomes CSD method, which satisfies the Theorem 1. On the other hand, Dai [START_REF] Dai | Alternate step gradient method[END_REF] proposed a variant called cyclic Barzilai-Borwein (CBB) method. They suggested that

α CBB k = α BB1 k , k mod m = 0, α k-1 , otherwise. (14)
Similarly, if we choose ρ(k) = 0 and τ (k) = k, k + 1, . . . , k -1, then (10) becomes CBB method.

Although these methods greatly speed up the convergence, their motivation is too straightforward to further accelerate the iterations, which relies on the nonmonotone property to search the whole space without sink into any lower subspace spanned by eigenvectors [START_REF] Fletcher | On the Barzilai-Borwein method[END_REF]. This allows to reduce the gradient components more or less in the same asymptotic rate [START_REF] Dai | Analysis of monotone gradient methods[END_REF].

The recent literature showed that Yuan steplength may lead to efficient algorithms [START_REF] Yuan | A new stepsize for the steepest descent method[END_REF][START_REF] Dai | Analysis of monotone gradient methods[END_REF]. All methods therein have two-dimensional finite termination property, i.e., if ( 8) is applied to a linear system in two-dimensional space, then the algorithm will terminate in at most 3 iterations. In general, such property seems not attractive in practice. However, experiments showed that they perform well in higher dimensions and are competitive with BB methods for large-scale problems [START_REF] Dai | Analysis of monotone gradient methods[END_REF].

Inspired by the Yuan steplength, we suggest a simple way of modifying steepest descent model to a cyclic gradient method. Consider a steplength of the form

α YB k = α SD k , k mod 3 = 0 or 2, α Y k , k mod 3 = 1. (15)
Here we modify the order of SD and Y compared to the original YB formula, which is useful for the development of the new algorithm. Apart from this change, ( 15) is indeed the second algorithm propose by the pioneering work of Yuan [START_REF] Yuan | A new stepsize for the steepest descent method[END_REF]. It keeps the two-dimensional finite termination property that performs as well as BB for large-scale problems and better for small-scale problems. We could introduce simply the cyclic behavior based on (15

) of the form ∀m ∈ N, if k mod (3 + m) > 2, then α k = α k-1 . (16)
Besides, we find that De Asmundis et al. [START_REF] De Asmundis | On spectral properties of steepest descent methods[END_REF] gives an interesting view about the iterations of SD method, where the technique of alignment was proposed therein to force the gradients into one-dimensional subspace and avoid the zigzag pattern. Notice that the inverse of constant Rayleigh quotient such as SD and BB steplengths has the similar property. Thus, constant SD with retards can also give rise to the alignment behavior and keep the nonmonotone benefit. To achieve this goal, we need to impose a repeat time to the zigzag process. Meanwhile, we want to keep the process based on Yuan steplength in the first several iterations. These motivations lead to a new method of the form

α CY k =      α Y k , k mod (l + m + 2) = 1, α SD k , k mod (l + m + 2) < l + 2, α k-1 , otherwise, (17)
where l ≥ 1 and m ≥ 1. Such formula seems complicated, but indeed easy to understand. There are three components consisting in [START_REF] Magoulès | Asynchronous optimized Schwarz methods with and without overlap[END_REF]: the first SD and Y are used to insure the finite termination property; the parameter l acting on the second part of SD is used to keep several zigzag iterations; finally, the retard term m induces alignment and provides nonmonotone behavior to leap from the lower subspace.

Convergence Analysis

By the invariance property under any orthogonal transformation, we can assume without loss of generality that

A = diag(λ 1 , . . . , λ n ), (18)
where

1 = λ 1 ≤ • • • ≤ λ n . (19)
We follow the convergence framework established by Dai [START_REF] Dai | Alternate step gradient method[END_REF] and adapt it to our method. Let

G(k, µ) = µ i=1 g 2 i,k , (20)
where g i,k is the ith component of g k . A preliminary property is defined as follows.

Definition 2 (Property A). Suppose that matrix A has the form [START_REF] Magoulès | Asynchronous iterative sub-structuring methods[END_REF] with condition [START_REF] Raydan | On the Barzilai and Borwein choice of steplength for the gradient method[END_REF] holds. If ∃ξ ∈ N, ∃M 1 , M 2 > 0, such that ∀µ ∈ {1, . . . , n -1}, ∀ > 0, ∀j ∈ {0, . . . , min{k, ξ} -1},

• λ 1 ≤ α -1 k ≤ M 1 ; • if G(k -j, µ) ≤ and g 2 µ+1,k-j ≥ M 2 , then α -1 k ≥ 2 3 λ µ+1 , then the steplength α k has Property A.
The convergence framework of Dai can be deduced from Property A, stated as follows.

Theorem 3 [START_REF] Dai | Alternate step gradient method[END_REF]. Consider the linear system (1) with A ∈ R n×n of the form [START_REF] Magoulès | Asynchronous iterative sub-structuring methods[END_REF] and b ∈ R n . Consider the gradient method (3) being used to solve [START_REF] Akaike | On a successive transformation of probability distribution and its application to the analysis of the optimum gradient method[END_REF]. If the steplength α k has Property A, then the sequence { g k } converges to 0 R-linearly for any starting point x 0 .

For a proof of the above theorem, see [START_REF] Dai | Alternate step gradient method[END_REF]. Many gradient methods have Property A as mentioned in [START_REF] Dai | Alternate step gradient method[END_REF], e.g., the gradient method with retards [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF]. Inspired by the demonstration therein, we now develop a convergence result for the CY method. Proof. Note that (17) has three alternate steplengths, whereas the SD updating process and the constant process using the last SD steplength both follow the framework [START_REF] Friedlander | Gradient method with retards and generalizations[END_REF], which has been proven to have the Property A [START_REF] Dai | Alternate step gradient method[END_REF]. Therefore, we only investigate the Yuan steplength.

Recall that Yuan steplength has the following property

1 α SD k-1 + 1 α SD k -1 < α Y k < min α SD k-1 , α SD k , (21)
which is given in [START_REF] Yuan | A new stepsize for the steepest descent method[END_REF]. Hence,

λ 1 ≤ 1 α SD k < 1 α Y k < 1 α SD k-1 + 1 α SD k ≤ 2λ n . (22)
Then the first condition of Property A holds by setting M 1 = 2λ n . For the second one, let M 2 = 2 and ξ = 1, which yields j = 0. Suppose that

G(k, µ) ≤ , g 2 µ+1,k ≥ M 2 , (23)
for all µ ∈ {1, . . . , n -1}, and > 0. Hence, the inverse of Yuan steplength becomes

1 α Y k > 1 α SD k = g k Ag k g k g k = n i=1 λ i g 2 i,k n i=1 g 2 i,k ≥ λ µ+1 n i=µ+1 g 2 i,k µ i=1 g 2 i,k + n i=µ+1 g 2 i,k ≥ λ µ+1 µ i=1 g 2 i,k g 2 µ+1,k + 1 ≥ λ µ+1 2 + 1 = 2 3 λ µ+1 (24)
Hence, the second condition of Property A is satisfied, which completes the proof.

Numerical Results

We first address the issue of parallel implementation. The dot product is engaged in the computation of steplength, which is the major obstacle of parallelization. Here we have two strategies to realize this goal. Let A i be the band matrix stored in the ith processor. The first one (Gather Algorithm, GA) is to gather the vector q i = A i * g and execute dot product with global vectors, shown as follows Allgatherv(q, q i ) α = Dot(g, g) / Dot(g, q) Then, the second one (Reduce Algorithm, RA) consists in computing the dot product locally, shown as follows c i = Dot(g i , q i ) Allreduce(c, c i , SUM) α = Dot(g, g) / c Besides, we can see that global gradient vector is used in each iteration that we must proceed another Allgatherv function to communicate with other processor. Let p be the number of processors. The two experiments are proceeded by Alinea [START_REF] Magoulès | Alinea: An advanced linear algebra library for massively parallel computations on graphics processing units[END_REF] (see also, e.g., [START_REF] Magoulès | Asynchronous optimized Schwarz methods with and without overlap[END_REF][START_REF] Magoulès | Asynchronous iterative sub-structuring methods[END_REF][START_REF] Magoulès | Asynchronous iterations of Parareal algorithm for option pricing models[END_REF]) and JACK [START_REF] Magoulès | JACK: an asynchronous communication kernel library for iterative algorithms[END_REF][START_REF] Magoulès | JACK2: An MPI-based communication library with non-blocking synchronization for asynchronous iterations[END_REF] (see also, e.g., [START_REF] Magoulès | Distributed convergence detection based on global residual error under asynchronous iterations[END_REF]) and results are illustrated in Figures 1 and2. We can see that generally the results are not good because the first one imposes so much computation and communication load, while the second one causes indeed the problem of loss of precision. These problems exist in all projection methods and by now we have not yet managed to find a solution.

The second experiments are proceeded by Matlab R2017b with a large-scale problem provided by The SuiteSparse Matrix Collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF], with n = 50000 and 349968 non-zero values. All parameters are chosen under a training problem given in [START_REF] Barzilai | Two-point step size gradient methods[END_REF], such that l = 4, m = 3 for CY, m = 3 for CSD, and m = 4 for CBB. Average results are shown in Table 1. We use bold numbers indicating the most efficient algorithms under each residual threshold. Backslash represents a number of iterations more than 10000. From Table 1, we see that the CY method performs better than other methods except CG. Although we do not yet beat CG in high precision, our method is still competitive in most cases. Specifically, we can see that CY is stable throughout the iterations and much better than CG when low accuracy is required.

Conclusions

In this paper, we have proposed a new gradient method and shown that it is very competitive with CG method and better than the others for large-scale problems. However, it is still lack of theoretical evidence supporting such results. It is also important to find a better parallelization strategy. Therefore it still remains to study the properties and high performance implementation of gradient methods.

Theorem 4 .

 4 Consider the linear system (1) with A ∈ R n×n of the form (18) and b ∈ R n . Consider the gradient method (3) with steplength (17) being used to solve (1). Then the steplength α CY k has Property A.

Figure 1 :

 1 Figure 1: Parallel CY method with GA implementation

Figure 2 :

 2 Figure 2: Parallel CY method with RA implementation

Table 1 :

 1 Gradient methods with different residual threshold 10 -1 10 -2 10 -3 10 -4 10 -5 10 -6

	CG	58	735	2617 4251 5786 7535
	CY	13	208 1153 4275 6181	\
	CSD	27	212 1470 4357 6713	\
	CBB	31	241	1965 7534	\	\
	DY	15	200	1595 6415	\	\
	BB1	356	984	2494 4612 8633	\
	SD	61	5773	\	\	\	\

Acknowledgment

This work was supported by the French national programme LEFE/INSU and the project ADOM (Méthodes de décomposition de domaine asynchrones) of the French National Research Agency (ANR).