Romain Guérout

Pavel Soldan

Mireille Aymar

Johannes Deiglmayr

Olivier Dulieu
email: olivier.dulieu@lac.u-psud.fr

Olivier Dulieu Core

P Soldán

Core repulsion effects in alkali trimers

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

Introduction: a brief overview on cold molecules

The research field on translationally cold (≈ 1K and lower) and ultracold (≈ 1mK and lower) molecules is continuously expanding in many directions, involving an increasing number of groups throughout the world. The availability of gaseous samples of cold neutral molecules opens entirely new avenues for fascinating researches, offering the possibility to control all degrees of freedom of a quantum system. Molecular ions can also be trapped and cooled down inside atomic ion traps behaving like ionic crystals, where each ion is kept at a specific position with a residual motion equivalent to a temperature smaller than 1 K [START_REF] Gerlich | The study of cold collisions using ion guides and traps[END_REF]. Cold neutral molecules brought new perspectives in high-resolution molecular spectroscopy [START_REF] Stwalley | [END_REF]3,4,5]. The expected accuracy of the envisioned measurements with ultracold molecules makes them appear as a promising class of quantum systems for precision measurements related to fundamental issues: the existence of the permanent electric dipole moment of the electron [6,7,8,9] related to CP-parity violation [10,11,12], and the time-independence of the electronto-nuclear and nuclear-to-nuclear mass ratios [13,14,15,16,17,18], or of the finestructure constant [19]. Various proposals have been suggested for achieving quantum information devices [20,21,22,23] based on cold polar molecules (i.e. molecules exhibiting a permanent electric dipole moment). Elementary chemical reactions at very low temperatures could be manipulated by external electric or magnetic fields, therefore offering an extra flexibility for their control [24]. The large anisotropic interaction between cold polar molecules is also expected to give rise to quantum magnetism [25], and to novel quantum phases [26,27]. The achievement of quantum degeneracy with cold molecular gases [28,29,30,31,32,33,34,35] together with the mastering of optical lattices [36] built a fantastic bridge between condensed matter and dilute matter physics. Indeed, the smooth crossover between the Bose-Einstein condensation (BEC) of fermionic atomic pairs and the Bardeen-Cooper-Schrieffer (BCS) delocalized pairing of fermions related to superconductivity and superfluidity, has been observed experimentally [35,37,38,39].

Cold and ultracold neutral molecules can be formed along two paths:

• Existing molecules in their vibrational ground state can be slowed down and cooled by various means, such as interactions with external electric and magnetic fields, or by collision with surrounding particles (see for instance the review articles in refs. [40,24,41]. This approach concerns a broad variety of small molecules (OH, NH, NO, SO 2 , ND 3 , CO, YbF, C 7 H 5 N, H 2 CO, LiH, CaH,...), but is currently limited to the production of molecules with a translational temperature around 10 mK.

• Pairs of trapped ultracold atoms can be associated using laser fields (photoassociation, or PA) [3], or time-varying magnetic fields (magneto-or Feshbach-association) [42]. Translational temperatures as low as a few tens of microKelvins can be reached, and even lower when quantum degeneracy is achieved. In contrast to the previous case, these results are up to now limited to alkali diatomic molecules which are most often formed in highly-excited vibrational levels, i.e. with high internal energy.

In the latter case, a breakthrough occurred in 2008, when the possibility to create ultracold bialkali molecules in the lowest vibrational level of their ground state or their lowest metastable triplet state, has been demonstrated. Caesium dimers have been created in their ground state v = 0 level after a PA step followed by spontaneous emission and further vibrational pumping, using a sequence of shaped laser pulses [43]. Dipolar molecules, namely LiCs, have been observed in their absolute ground state rovibrational level v = 0, J = 0 after PA and spontaneous emission [44]. Samples close to the degeneracy regime of ultracold molecules in the v = 0 level of their ground state have been observed for Cs 2 [45,46], and KRb [47], and in the v = 0 level of their lowest triplet state for Rb 2 [48], using the STIRAP technique (Stimulated Rapid Adiabatic Passage) to transfer the population from initial high-lying vibrational levels.

Motivation of the present work

The growing availability of ultracold samples of alkali diatomics brings the possibility to observe and study their interactions with neighboring ultracold atoms and molecules. In two independent -and almost simultaneous-experiments, inelastic rate constants for atom-molecule and molecule-molecule collisions have been extracted, using optically trapped Cs 2 molecules created by PA and spontaneous emission [49,50]. The rate constants have been found independent from the initially populated rovibrational level. A resonant feature observed in the loss rate after collisions between trapped ultracold Cs 2 molecules created by Feshbach association has been interpreted as Cs 2 -Cs 2 bound states [51]. Inelastic atom-molecule collisions in a trapped sample of RbCs molecules have been studied for both RbCs-Cs and RbCs-Rb cases [52]. No systematic dependence on the internal state of the molecule has been probed within the experimental precision.

Theoretical knowledge of the electronic structure of alkali triatomic systems is strongly needed for theoretical dynamical studies, which would support these experiments. There has been a significant research activity concerning the theoretical study of the structure and dynamical properties of homonuclear alkali-metal trimers. Non-additive effects in spin-polarized alkali-metal trimers were studied by Soldán et.al. [53] and the three-dimensional potential energy surfaces for the lowest quartet states were constructed for Li 3 [54,55], Na 3 [56,57], K 3 [58] and Rb 3 [59]. These were then used to study the corresponding spin-polarized reactive atom-dimer collisions at very low temperatures [55,58,60,61,62,63,64,65,66,67]. General trends for the lowest quartet state of homonuclear alkali trimers -correlated to three spin-polarized alkali atoms -have been described by Hutson and Soldán [68]. One of the main characteristic of these systems is the importance of non-additive three-body forces, which are predicted to be large especially close to the equilibrium geometries [53]. Much less has been done in the case of heteronuclear alkali-metal trimers. One of us has actually started with a systematic study the lowest quartet states of Li 2 A mixed systems (with A=Na, K, Rb, Cs) [69]. He concluded that the single-reference coupled-cluster approach, which has been successfully used for homonuclear alkali-metal trimers [58,69], would be for various reasons very difficult to employ for calculations of the potential energy surfaces of heteronuclear alkali-metal trimers, and that an alternative approach to the problem should be sought.

In the present paper, we present such an alternative approach. We extend our previous works on effective two-electron diatomic molecules like alkali dimers [70,71,72] and alkali hydrides [73] to alkali trimers. We propose a preliminary study of the potential energy surfaces of the heaviest alkali trimer, Cs 3 , modeled as an effective three-electron molecule. Our approach is based on large-core effective core potentials (ECP) with core polarization potentials (CPP). We focused our work on two aspects: (i) the estimation of repulsion effects between the three Cs + ionic cores at short distances within the present ECP+CPP approach; (ii) the comparison of the results with the results obtained by the MOLPRO package for ab-initio calculations.

Method of calculation

Ab initio methods and basis sets

The potential energies were calculated making use of the CIPSI package (Configuration Interaction by Perturbation of a multiconfigurational wave function Selected Iteratively) [74]. As in our previous studies on alkali dimers, the alkali atom A is described by an ℓ-dependent Effective Core Potential (ECP) for the ionic core A + [75,76] including effective core polarization potential (CPP) terms [77,78], and by a large set of uncontracted Gaussian functions for the valence electron. The atom is modeled as an effective one-electron system, permitting us to perform Full Configuration Interaction (FCI) calculations for any alkali dimer or trimer, considered as an effective two-electron or three-electron system, respectively.

As the distance between different cores becomes smaller than the equilibrium distance of the system, it is well known that their electrostatic repulsion is not accounted for by this effective large-core potential. At the same time, the addition of the CPP term leads to an non-physical attractive behavior at short distances. For diatomic molecules, an empirical repulsion term is usually added to the calculation (see e.g. Refs. [79,80,81,82]). If this short-range repulsive interaction in the triatomic molecule can be treated as the sum of pairwise interaction terms, computation of the potential surfaces would be greatly simplified. In order to check if such an approximation could be justified, we calculated the core-core repulsion at the Hartree-Fock level making use of the MOLPRO 2006.1 Quantum Chemistry Package [START_REF] Werner | Molpro, version 2006.1, a package of ab initio programs[END_REF]. For lithium and sodium, the all-electron correlation-consistent polarized core-valence quadruple-ζ cc-pCVQZ basis sets [START_REF] Oren | [END_REF] were used. For potassium, rubidium, and caesium, the small-core scalar relativistic effective core potentials ECP10MDF, ECP28MDF, and ECP46MDF, respectively, together with the corresponding uncontracted valence basis sets [85] were employed.

For the caesium dimer and trimer, the counterpoise-corrected dimer interaction energies were optimized using an algorithm implemented in MOLPRO. In this algorithm the interaction energies were calculated at the coupled-clusters level making use of a single-reference restricted open-shell variant [86] of the coupled cluster method [87] with single, double and non-iterative triple excitations [RCCSD(T)]. The small-core scalar relativistic effective core potentials ECP46MDF was employed together with its uncontracted valence basis set, which was augmented by one set of diffuse function in an even-tempered manner (aug-ECP46MDF). All electrons from the "outer-core" 5s5p orbitals were included in the RCSSD(T) calculations. The full counterpoise correction of Boys and Bernardi [88] was applied to all the interaction energies in order to compensate for the basis set superposition errors.

Core Polarization Potential for a triatomic molecule

The core polarization potential (CPP) is given by [77]:

V CP P = - 1 2 c α c f c • f c (1)
The electric field f c is the one produced at point r c by all the other electrons and the ionic cores with static dipole polarizability α c :

f c = i r ci r 3 ci h(r ci , ρ c ) - c ′ =c R cc ′ R 3 cc ′ Z c ′ (2)
where r ci is the position vector from core c to electron i, R cc ′ the one between cores c and c ′ . The cutoff function h(r, ρ) ensures the convergence of the integral over electronic coordinates. The V CP P term contains a purely geometrical part v nn depending only on the relative positions of the nuclei:

v nn = - 1 2 c α c c ′ ,c ′′ =c R cc ′ • R cc ′′ R 3 cc ′ R 3 cc ′′ Z c ′ Z c ′′ (3)
For a diatomic molecule, v nn is given by the familiar formula:

v nn = - Z 2 2 α 1 +Z 2 1 α 2 2R 4
. For a triatomic system, we have:

v nn = - Z 2 2 α 1 + Z 2 1 α 2 2R 4 12 - Z 2 3 α 1 + Z 2 1 α 3 2R 4 13 - Z 2 3 α 2 + Z 2 2 α 3 2R 4 23 - Z 2 Z 3 α 1 R 12 • R 13 R 3 12 R 3 13 - Z 1 Z 3 α 2 R 12 • R 23 R 3 12 R 3 23 - Z 1 Z 2 α 3 R 13 • R 23 R 3 13 R 3 23 (4)
In those formulas, Z i and α i are, respectively, the net charge and static polarizability of ion core i. This term is independent of the ab initio calculation itself and is added a posteriori.

Volume effect for triatomic systems

The part of the interaction between two nuclei other than the point charge repulsion has been dubbed "volume effect" by Jeung in his 1997 paper on alkali diatomics [89]. Here, we are interested in this volume effect for alkali triatomic systems. More precisely, we investigate whether it can be reasonably approximated as a sum of two-body terms.

The two-body repulsion energies V cc are calculated for all the alkali homonuclear diatomic molecules.

The three-body repulsion energies V ccc are calculated for particular geometries, namely for D ∞h and D 3h symmetries, to provide a representative investigation. Calculations ar done at the Hartree-Fock level of theory starting with the orbitals of the free atoms. Those orbitals are kept frozen as the internuclear distance is shortened. In this respect, our calculations for the two-body term are similar to those of the "A" column from the tables of ref. [89]. This comparison is shown in Figure 1 in logarithmic scale. Apart from a slight discrepancy in the case of Na 2 , there is a good agreement between both calculations, and it is seen that the two-body core repulsion term could be well fitted by an exponential form, as already stated in ref. [80]. Three-body core repulsion energies have been calculated for all homonuclear alkali trimers at linear symmetric and equilateral geometries according to the schemes of Figure 2. In Figure 3, we compare the three-body core repulsion term for linear geometries V lin ccc (R) (the internuclear distance R being defined in Figure 2), with the sum of two-body term 2V cc (R) + V cc (2R). In the same manner, we compare in Figure 4 the three-body repulsion term for equilateral geometries V eq ccc (R) to 3V (2) (R). In both cases, it is seen that the three-body core repulsion term is well described by a pairwise additive lemma, up to the region where it becomes negligible. Therefore, just like for alkali dimers, the V ccc term can be added a posteriori to the potential surface calculations, whatever the chosen grid is for them.

R R R θ R' r' R R Figure 2
. Linear (D ∞h ) and equilateral (D 3h ) geometries at which the three-body core repulsion term have been calculated. The Jacobi coordinates {r ′ , R ′ , θ} used in the calculation of the potential energy surfaces of Cs 3 are also presented. Blue : Three-body term V lin ccc (R). Red dashed : sum of the relevant two-body terms 2V cc (R) + V cc (2R) (see text).

Preliminary application to caesium trimer

The inclusion of the three-body core repulsion term is considered in preliminary calculations concerning the 4 A ′ 2 state of Cs 3 which correlates to three spin-polarized caesium atoms. The used basis set is kept relatively small to maintain reasonable size for the FCI, reaching 152 532 Slater determinants in C s symmetry in which the calculation is performed. The basis set (5s3p4d/4s3p3d) is detailed in Table 1. We use ℓ-dependent CPP's [78] which allow us to adjust the calculated atomic energies to the experimental ones for the 6s 2 S, 6p 2 P and 5d 2 D atomic states of caesium. This is achieved by tuning the cutoff radii ρ ℓ of Eq. 2 to the values displayed in Table 2.

The triplet state of Cs 2

With these basis set and cut-off radii, we first calculate the lowest 3 Σ + u potential energy curve for Cs 2 dissociating into Cs(6s)+Cs(6s) in the framework of our ECP+CC-FCI Blue : Three-body term V eq ccc (R). Red dashed : sum of the relevant two-body terms 3V (2) (R) (see text). method. We found an equilibrium internuclear distance R e = 11.82 a 0 and a well depth of D e = 380.6 cm -1 . These numbers can be compared to our optimized results obtained at the RCCSD(T)/aug-ECP46MDF level R e = 12.19 a 0 and D e = 256.4 cm -1 and the results from Ref. [53,68], who obtained R e = 12.44 a 0 and D e = 246.8 cm -1 at the RCCSD(T)/ECP46MWB level of theory. The MOLPRO calculation yields a potential well which is deeper by is a more than 140 cm -1 than the one of the ECP+CC-FCI calculation. Apparently, the small-core ECP in combination with a rich valence basis set in the RCCSD(T) calculations provides more realistic description of the interaction energy than the large-core ECP in combination with a smaller valence basis set used in the FCI calculations (see also the discussion in Section 6). The potential energy well of the triplet state results from the competition between the exchange energy at short distances and the attractive dispersion forces at long distances and is very sensitive to the quality of the basis. For this 2. ℓ-dependent cutoff radii ρ ℓ (in units of a 0 ) used in the core polarization potential. The polarization for Cs + ionic core of 16.33 a 3 0 is taken from Ref [90].

term does not change the characteristics of the well but ensures a correct exponential behavior at short internuclear distances (Figure 5).

Quartet state of Cs 3

We then calculate the potential energy surface of the lowest quartet state of Cs 3 . Calculations are carried out in Jacobi coordinates {r ′ , R ′ , θ} (depicted in Figure 2) with the angle θ fixed at π/2 in order to explore C 2v nuclear geometries, as a representative geometry. We have computed 984 points on a grid where 4.4 a 0 < r ′ < 40 a 0 and 4 a 0 < R ′ < 40 a 0 . The minimum of the surface is found at a D 3h geometry with a bond length of R ′ e = 10.72 a 0 and a well depth of D e = 1518 cm -1 with the inclusion of the pairwise three-body core repulsion term, which ensures a realistic short-range repulsive wall. These numbers can be compared to our optimization results obtained at the RCCSD(T)/aug-ECP46MDF level R ′ e = 11.18 a 0 and D e = 1208 cm -1 and to the results from Ref. [53,68], who obtained R ′ e = 11.33 a 0 and D e = 1139 cm -1 at the RCCSD(T)/ECP46MWB level of theory. It is worth noting that as the minimum of the quartet potential well is located at shorter distances than the minimum of the Cs 2 u state of Cs 2 . Black plain curve : ECP+CPP-FCI calculation with the basis detailed in Table 1. Red dashed curve : obtained at the RCCSD(T)/aug-ECP46MDF level. Green chain dotted curve : ECP+CPP full CI result obtained with basis set "A" from Ref. [70] (see the discussion in Section 6). Blue dotted curve : Multiparameter Morse Long Range fitting to spectroscopic data measured by Xie et al [91].

triplet state, the V ccc term indeed contributes to the depth of the well, and not only in the region of the repulsive wall. However, the magnitude of V ccc is still small, as the equilibrium internuclear distance is R ′ e = 10.7 a 0 and the well depth is D e = 1522 cm -1 if we neglect it. As expected from the diatomic calculations above, we obtain a difference of about 300 cm -1 on the well depth of the quartet state of Cs 3 compared to the MOLPRO calculations. This confirms the limited quality of our basis set at the current level of our computations.

Doublet states of Cs 3

The three lowest potential surfaces for Cs 3 are presented in Figure 6. The 2 B 2 and 2 A 1 states are the two components of a 2 E ′ state at D 3h symmetry subjected to Jahn-Teller effect. The left (resp. right) column shows the surfaces without (resp. with) the addition of the core repulsion term. Note the strongly unphysical attractive behavior in the region r ′ < 5 a 0 , which is removed when V ccc is added. The global ground state 2 B 2 is characterized by a well depth of 5437.1 cm -1 at the geometry {r ′ = 10.58 a 0 , R ′ = 7.36 a 0 } with the inclusion of the core repulsion term. This changes to a well depth of 5489 cm -1 at {r ′ = 10.57 a 0 , R ′ = 7.3 a 0 } without the core repulsion term. Once again, the effect of V ccc is larger than in the dimer, as the equlibrium distance is shorter. The other component 2 A 1 is characterized by a well depth of 5258.2 cm -1 at {r ′ = 8.75 a 0 , R ′ = 9.28 a 0 } with the core repulsion term (5312 cm -1 at {r ′ = 8.66 a 0 , R ′ = 9.26 a 0 } without it). We present in Figures 7 and8 the cut through the D 3h geometry of the previous surfaces. We recall that in the D 3h symmetry point group, the two doublet states 2 B 2 and 2 A 1 discussed before correlates to the two components of a doubly degenerate 2 E ′ state. We clearly see in these figures the abrupt attractive behavior which occurs at short distances and that the addition of the core repulsion term gives us back a realistic repulsive wall of potential.

Internuclear distance (a₀) Potential energy (cm¹) Plain black (resp. red dashed) curve : with (resp. without) the core repulsion term V eq ccc (R). Both components of the 2 E ′ state are represented.

Future prospects

Preliminary potential energy surfaces for doublet and quartet states of the Cs 3 molecule have been computed in the framework of a ECP+CPP-FCI quantum chemistry approach. The core repulsion interaction has been taken into account as an additional empirical correction as it is routinely done for diatomic molecules calculations. We showed that this term is well described by a sum of two-body terms. As mentioned before in ref. [80], it can be fitted by an exponential formula. Therefore, this term needs not to be calculated at each point of the three dimensional grid spanned by the degrees of freedom of a triatomic molecule on which the ab initio calculations are performed. It can easily be added a posteriori to cancel any non-physical behavior which might occur at short internuclear distances.

To estimate the quality of this basis set, we have compared our calculations to high level ones which uses a different method than ours, with the MOLPRO package. Such a comparison between different methods is crucial to assess their accuracy and their consistency, as no spectroscopic data are available for the Cs 3 system, as well as for most of the triatomic alkali systems. The calculations presented here are exploratory in nature: the small basis set used (Table 1) makes the calculation time short enough to generate over a thousand ab initio points. Figure 5 shows the limitation of this basis set where the well depth of the triplet state of Cs 2 is overestimated by more than 100 cm -1 . For the quartet state of Cs 3 , our calculated well depth is consequently 300 cm -1 lower than the RCCSD(T) calculations. If we employ the basis set previously used for the caesium dimer [70] in the ECP+CPP-FCI method, this yields a very good agreement with the RCCSD(T) calculations, as the difference is now found at about 20 cm -1 between them (Figure 5). The agreement is even more spectacular with the potential curve extracted from the spectroscopy of the lowest triplet state of Cs 2 [91].

We are currently calculating a comprehensive set of ab initio points with the ECP+CPP-FCI with this extended basis for the Cs 3 molecule which will be the subject of a future paper. Such FCI calculations now involve more than 500,000 Slater determinants. In Figure 9, we display the first calculated points for the D ∞h symmetry of the lowest Cs 3 quartet state, compared to the MOLPRO calculation as described in Section 5. The figure suggests that the results are in good agreement with each other as the position of the minimum looks similar. The point calculated around 12.3 a 0 with the ECP+CCP full CI method is deeper than the minimum of the MOLPRO curve by about 40 cm -1 , while the difference between corresponding triplet potential curves in Figure 5 is seen to be about 20 cm -1 This suggests that pairwise additivity of the forces is a reasonably good approximation in this case. Green chain dotted curve : calculation at the RCCSD(T)/aug-ECP46MDF level of theory. Plus signs : preliminary calculations with the ECP+CPP-FCI approach using basis set "A" from Ref. [70].

Science Foundation. J.D. acknowledges partial support of the French-German University (http://www.dfh-ufa.org), and R.G. support from Institut Francilien de Recherches sur les atomes froids -IFRAF (http://www.ifraf.org).

Figure 1 .

 1 Figure 1. Two-body core repulsion energies for alkali dimers. Blue : our work. Red dashed : Ref. [89]. Note the logarithmic scale used for the potential energy axis.

Figure 3 .

 3 Figure 3. Additivity of the three-body core repulsion term for linear geometries. Blue : Three-body term V lin ccc (R). Red dashed : sum of the relevant two-body terms 2V cc (R) + V cc (2R) (see text).

Figure 4 .

 4 Figure 4. Additivity of the three-body core repulsion term for equilateral geometries. Blue : Three-body term V eq ccc (R). Red dashed : sum of the relevant two-body terms 3V(2) (R) (see text).

Figure 5 .

 5 Figure 5. Calculated potential energy curves for the lowest3 Σ + u state of Cs 2 . Black plain curve : ECP+CPP-FCI calculation with the basis detailed in Table1. Red dashed curve : obtained at the RCCSD(T)/aug-ECP46MDF level. Green chain dotted curve : ECP+CPP full CI result obtained with basis set "A" from Ref.[70] (see the discussion in Section 6). Blue dotted curve : Multiparameter Morse Long Range fitting to spectroscopic data measured by Xie et al[91].

Figure 6 . 3 2

 63 Figure 6. Lowest three calculated surfaces for Cs 3 in Jacobi coordinates with fixed θ = π/2. Left (resp. right) column : without (resp. with) the inclusion of the cores repulsion term V ccc . The red dashed curve locates the equilateral D 3h geometry R ′ = √ 3 2 r ′ . Contours are distant of 500 cm -1 and the dashed contour is the energy of the dissociation Cs( 2 S) + Cs( 2 S) + Cs( 2 S).

Figure 7 .Figure 8 .

 78 Figure 7. Cut through the D 3h geometry for the lowest quartet states of Cs 3 . Plain black (resp. red dashed) curve : with (resp. without) the core repulsion term V eq ccc (R). Green chain dotted curve : quartet state on the RCCSD(T)/aug-ECP46MDF level of theory (see the discussion in section 6.

Figure 9 .

 9 Figure 9. Quartet 4 Σ + u state of Cs 3 in D ∞h symmetry.Green chain dotted curve : calculation at the RCCSD(T)/aug-ECP46MDF level of theory. Plus signs : preliminary calculations with the ECP+CPP-FCI approach using basis set "A" from Ref.[70].

Table 1 .

 1 Gaussian basis set used on each caesium atom in the ECP+CPP-FCI approach.

	Angular momentum Exponent Contraction coefficients
	s	0.347926	0.411589
	s	0.239900	-0.682422
	s	0.050502	1.
	s	0.036900	1.
	s	0.00515	1.
	p	0.1837	1.
	p	0.0655	1.
	p	0.0162	1.
	d	0.2106894	0.18965
	d	0.065471	0.22724
	d	0.021948	1.
	d	0.011200	1.
		ℓ	ρ ℓ
		s 2.6248
		p 1.87
		d 2.8111

3

Σ + u curve, the inclusion of the two-body core repulsion

Table

Acknowledgments

This work is performed in the framework of the network "Quantum Dipolar Molecular Gases" (QuDipMol) of the EUROCORES/EUROQUAM program of the European


Ilias Rmouque

Maxime Devanne

Jonathan Weber

Germain Forestier

Cédric Wemmert
email: wemmert@unistra.fr

Ensemble clustering for histopathological images segmentation using convolutional autoencoders

Keywords: Digital Pathology, Deep Learning, Autoencoders, Histopathology, Ensemble Learning, Segmentation

Unsupervised deep learning using autoencoders has shown excellent results in image analysis and computer vision. However, only few studies have been presented in the field of digital pathology, where proper labelling of the objects of interest is a particularly costly and difficult task. Thus, having a first fully unsupervised segmentation could greatly help in the analysis process of such images. In this paper, many architectures of convolutional autoencoders have been compared to study the influence of three main hyperparameters:

(1) number of convolutional layers, (2) number of convolutions in each layer and (3) size of the latent space. Different clustering algorithms are also compared and we propose a new way to obtain more precise results by applying ensemble clustering techniques which consists in combining multiple clustering results.

INTRODUCTION

Pathology is essential for the diagnosis evaluation and understanding of many underlying biological and physiological mechanisms. It is usually a visual evaluation by pathologists of a tissue sample using a microscope to identify its structural properties. Currently, the visual evaluation of microscopic specimens is largely an unassisted process, and the pathologist's accuracy is established through extensive training, comparative analysis, peer quality control and personal experience. However, this field has undergone several technological revolutions in recent years with the advent of virtual microscopy (conversion of glass slides into high-resolution images called Whole Slide Images -WSI), often referred to as "digital pathology". Thus, major efforts have been made to design image analysis tools, for example to identify basic biological structures (stroma, immune cells, tumour, etc.), in order to make it easier for doctors to (semi-)automate the interpretation of slides. Meanwhile, automatic image analysis algorithms have recently made extraordinary progress, particularly with the advent of the deep learning methods introduced by Lecun et al. [START_REF] Lecun | Deep learning[END_REF]. Indeed, the performances of these methods have exploded in recent years, allowing the detection, classification and segmentation of objects of interest in images and particularly in medical images with high precision [START_REF] Bukala | Classification of histopathological images using scaleinvariant feature transform[END_REF][START_REF] Ando | Cell image segmentation by feature random enhancement module[END_REF]. But most of these approaches operate in supervised mode, i.e. they require many examples in order to provide an effective model. However, obtaining quality annotations on histopathological images remains very costly. For example, in the field of colorectal cancer WSI segmentation, Qaiser et al. proposed a method based on persistent homology to classify tumour and non-tumour patches from Hematoxylin & Eosin stained histology images [START_REF] Qaiser | Persistent homology for fast tumor segmentation in whole slide histology images[END_REF]. To train their system, more than 18,000 annotated patches were needed. At the same time, unsupervised approaches have shown their interest in many applications for image analysis, such as remote sensing [START_REF] Liang | Remote sensing image classification based on stacked denoising autoencoder[END_REF][START_REF] Mei | Unsupervised spatial-spectral feature learning by 3d convolutional autoencoder for hyperspectral classification[END_REF]. Recently, they have also been applied to histopathological WSI analysis for cells segmentation [START_REF] Junior | Unsupervised segmentation of leukocytes images using particle swarm[END_REF] or regions of interest classification [START_REF] Figueira | Adversarialbased domain adaptation networks for unsupervised tumour detection in histopathology[END_REF]. In particular for cancer, the authors in [START_REF] Yamamoto | Automated acquisition of explainable knowledge from unannotated histopathology images[END_REF] describe an unsupervised approach for extracting interesting information from WSI that obtains better accuracy than human for prognostic prediction of prostate cancer recurrence.

In this paper, we are interested in automatic segmentation in order to quickly extract regions of interest (tumours for example) to make a more precise analysis of these areas only. However, only few approaches on fully unsupervised segmentation of WSI have been proposed. The first attempt to segment regions of interest from WSI without any prior information or examples has been performed in [START_REF] Khan | Hymap: A hybrid magnitude-phase approach to unsupervised segmentation of tumor areas in breast cancer histology images[END_REF]. The authors highlight tissue morphology in breast cancer histology images by calculating a set of Gabor filters to discriminate different regions. In [START_REF] Fouad | Unsupervised morphological segmentation of tissue compartments in histopathological images[END_REF], the authors use mathematical morphology to extract 'virtual-cells' (e.g. superpixels), for which morphological and colour features are calculated to then apply a consensus clustering algorithm to identify the different tissues in the image. More recently, a similar approach has been presented in [START_REF] Landini | Epithelium and stroma segmentation using multiscale superpixel clustering[END_REF], adding a semi-supervised self-training classifier to the previous techniques that enhances the results at the cost of partial supervision. All these approaches propose to cluster the image based on predefined features. However, deep learning approaches, Elsevier particularly via autoencoding architectures, make it possible to avoid manual definition of features by calculating a condensed representation of the image in a latent space by applying convolutional filters. Unfortunately, as stated in [START_REF] Raza | A tour of unsupervised deep learning for medical image analysis[END_REF], most applications of autoencoders in digital pathology were developed to perform cell segmentation or nuclei detection [START_REF] Xu | Stacked sparse autoencoder (ssae) for nuclei detection on breast cancer histopathology images[END_REF][START_REF] Hou | Sparse autoencoder for unsupervised nucleus detection and representation in histopathology images[END_REF], or stain normalisation [START_REF] Janowczyk | Stain normalization using sparse autoencoders (stanosa): application to digital pathology[END_REF]. Therefore, we propose here to study the potential of these approaches for WSI tissue segmentation. The aim is to try to automatically identify clusters corresponding to each type of tissue in the WSI that could then be labelled by pathologists.

In this paper, we present a study on how convolutional autoencoders perform on WSI segmentation by comparing different approaches. First, different autoencoders architectures are compared to quantify the importance of hyperparameters of interest (number of convolutional layers, number of convolutions by layer and size of the latent space). Then, a multi-resolution approach using an ensemble clustering framework is evaluated, to see if such ensemble techniques could provide more accurate results.

METHODS

Convolutional autoencoders

In this section, we explore of the use of convolutional autoencoders to cluster WSI histopathological images. For this, we present several experiments to evaluate the importance of each hyperparameter.

As shown in Figure 1, a Convolutional AutoEncoder (CAE) is a deep convolutional neural network composed of two parts: an encoder and a decoder. The main purpose of the CAE is to minimise a loss function L, evaluating the difference between the input and the output of the CAE (usually Mean Squared Error). Once this function is minimised, we can assume that the encoder part builds up a suitable summary of the input data, in the latent space, as the decoder part is capable of reconstructing an accurate copy of it from this encoded representation.

The encoder is first constituted of the input layer (having the size of the input image) which is connected to N convolutional layers of diminishing size, up to an information bottleneck of size Z, called the latent space. The bottleneck is connected to a series of N convolutional layers of increasing size, until reaching the size of the input. This second part is called the decoder. Each convolution layer is composed of C convolutions and is followed by three other layers: a batch normalisation, an activation function (ReLU) and a max pooling of size (2,2).

To perform the clustering, a trained CAE is used to encode each patch of the whole image. Then, this encoded representation of the patch (in the latent space) is given as the input of a clustering algorithm and a cluster is assigned to the patch.

We decided to evaluate the influence of the three hyperparameters N, Z and C. For each one, different values were tested while fixing the two others (N = 2, Z = 250, C = 10). To evaluate the quality of the results, the Adjusted Rand Index (ARI) is calculated to compare the obtained clustering to the annotations of the expert. The Rand Index computes a similarity measure between two clusterings by considering all pairs of samples and counting pairs that are assigned in the same or different clusters in the predicted and true clusterings. The score is then normalised into the ARI score by:

ARI = (RI -Expected RI) (max(RI) -Expected RI) (1)
Values of the ARI are close to 0 for random labelling independently of the number of clusters and samples, and exactly 1 when the clusterings are identical (up to a permutation).

Each CAE was trained over a set of 10,000 different patches randomly selected. As the result of both the clustering and the training of the CAE are non-deterministic, due to a high sensitivity to the initial conditions, 10 autoencoders were trained and the results averaged for each hyperparameter value.

We also investigated the performance of several clustering algorithms, i.e Kmeans, Agglomerative clustering (AggCl), Gaussian mixture (GM) and also the not too deep clustering method (N2D) exposed in [START_REF] Mcconville | N2d: (not too) deep clustering via clustering the local manifold of an autoencoded embedding[END_REF]. A clustering performed directly with the Kmeans algorithm on the raw data (without any data reduction by the CAE) has been calculated as a baseline to evaluate the benefit of encoding the data with the CAE.

Ensemble clustering

As exposed in [START_REF] Yamamoto | Automated acquisition of explainable knowledge from unannotated histopathology images[END_REF], both micro-structures and macro-structures give different information. Pathologists also agree that identifying a single cell is way more difficult without its surrounding context and they always look at the WSI at lower magnification (to better capture the context) before zooming in at high magnification. Furthermore, in [START_REF] Alsubaie | A multi-resolution deep learning framework for lung adenocarcinoma growth pattern classification[END_REF] an example of multi-resolution lung cancer adenocarcinoma classification using deep learning shows improvements in the overall accuracy.

Thus, we explored a way to improve the results by using an ensemble of clustering methods, each focusing on a different resolution. The objective is to merge low level information (context) with high level information (shape of the cells, etc.). For this, the consensus method proposed in [START_REF] Wemmert | A multi-view voting method to combine unsupervised classifications[END_REF] was used. This method is based on a the evaluation of the similarity between different clusterings and the definition of corresponding clusters. Then, a multi-view voting approach is computed to produce a single result representing all clusterings. An example of the architecture of the approach is depicted in Figure 2.

We explored different configurations, but we only present the two most representative which highlight how the quality of the results can be improved by using ensemble clustering. The first configuration, E multires is composed of three clustering algorithms (Kmeans) working on the latent space representation of the image obtained by different CAE trained at different resolutions: 10× with 8 clusters, 5× with 6 clusters and 5× with 8 clusters. As the reconstructed image from the autoencoder seems to focus more on colour intensity than real structures, a second ensemble configuration has been tested. To add diversity and to force the final result to focus its attention more on the structure of the objects, a clustering working on a binary image (by thresholding the intensity of the initial image) has been computed. Thus, the second configuration (E struct ) is composed

EXPERIMENTS AND RESULTS

Our study was performed on 8 WSI of Haematoxylin Eosin Saffron (HES) stained tissue extracted from a cohort of patients built within the scope of the AiCOLO project (INSERM/Plan Cancer) studying colon cancer. The images have been provided by Georges Franc ¸ois Leclerc Centre (Dijon, France) and acquired from two different centres. An example is given in Figure 3a. HES stain distinguishes cell nuclei in purple, from extracellular matrix and cytoplasm in pink.

All images have been acquired at 20× magnification (corresponding to 0.5 µm/pixel) but stored at several resolutions in a pyramidal format. The size of each image is around 90, 000 × 50, 000 pixels. To train autoencoders, 10, 000 patches of size 128 × 128 pixels were randomly extracted at 10× resolution from all images (and 5× for the ensemble approach), as this seems to be the minimal amount of information required by human expert to classify the tissue. Meanwhile, sparse manual annotations of the five classes of tissue, tumour, stroma, outer layer mucosa (crypts of Lieberkuhn and connective tissue), immune cells, and necrosis, and two classes for background and artefacts (ink marks, etc.) have been performed by pathologists on the images (using Cytomine [START_REF] Marée | Collaborative analysis of multi-gigapixel imaging data using Cytomine[END_REF]), to be able to evaluate the relevance of the clustering.

Evaluation of all hyperparameters of the CAE

First, results obtained without using the latent space representation (see Table 1) are worse than all those obtained when clustering the encoded data. This confirms the interest of using a CAE for WSI clustering. As shown in Figure 4a, it appears that the number of convolutions in each layer of convolutions (hyperparameter C) does not greatly affects the quality of the autoencoder as only a apart from a slight narrowing of the variability of the results. It's quite easy to figure out why: passed a certain number, additional convolution brings to few complementary information. Figure 4b shows the evaluation of the ARI with different number of convolution layers in the architecture. We can notice an increase of the quality index up to 4 layers and then a brutal drop at 5. This indicates clearly that too many convolutions (and poolings that downsample the information) reduce the information that can further not be properly processed.

Nonetheless, as seen in Figure 5, the latent space size Z, seems to greatly influence the pertinence of the CAE. Indeed, the ARI clearly grows as there is more space to encode the latent representation, as a more precise information can be stored. Also, the more information is present in the latent representation, the more classes can easily be differentiated. However, it is also clear that a too large latent space will not be able to summarise efficiently the information, and thus, will not help the clustering algorithm to discriminate the different tissues. Moreover, the larger the latent space, the more memory and time are needed to train the network.

Comparison of the CAE with the ensemble approach

As seen in the previous experiment, the ARI tends to give low scores because we only have very few annotations on each class of interest. So we decided to compute a second evaluation criterion based on the ability of the clustering to detect tumours areas in the image, as it is the main class of interest in our project. To associate the tumour class to a cluster, we calculated its tumour density (number of labelled tumour pixels / number of total labelled pixels in the cluster). All clusters having a density over 50% are kept as 'tumour', the others are labelled as 'not tumour'. Thus, two evaluation criteria have been calculated on the results and are presented in Table 1: the ARI as in the previous experiment (see Eq.1) and the FScore on the two-classes problem (tumour vs. not tumour) [START_REF] Van Rijsbergen | Information retrieval: theory and practice[END_REF].

Discussion

Classical methods applied on the latent space representation of the CAE tend to show acceptable results. However, both ensemble clustering configurations seem to be more efficient in finding coherent clusters corresponding to the classes of interest defined by the pathologists.

Among all the exposed methods, E struct seems to give the best results. It tends to confirm the importance of the shape of the objects on histopathological images. Furthermore, it shows that even if convolutional autoencoders aim at automatically finding the best features to encode images, they can also take advantage of pre-computed features for some specific tasks.

CONCLUSION

In this paper, we compared different configurations of convolutional autoencoders in the field of unsupervised learning for WSI histopathological image segmentation. For this, different CAE architectures have been compared to try to find the best configuration and to study the influence of each hyperparameter. Then, we proposed a new approach that uses ensemble clustering technique to take advantage of multiresolution information and structural features in the image. This confirms the importance of having diversity in an ensemble learning framework and that working at different resolutions at the same time can really improve the quality of the results.

Figure 1 :

 1 Figure 1: Architecture of a CAE with N = 2, C = 10 and Z = 50

Figure 2 :

 2 Figure 2: Architecture of the first ensemble configuration E multires : two CAE trained at different resolutions (10× and 5×) produce different latent representations that are clustered. The three resulting clusterings are then merged through the multi-view voting algorithm proposed in [19].

Figure 3 :

 3 Figure 3: Example of (a) a raw WSI and (b) a clustering result of this image.

Figure 4 :

 4 Figure 4: Evaluation of the ARI for the two main hyperparameters of the convolutions of the CAE comparing Kmeans clustering on 7, 8, 9 and 10 clusters.

5 :

 5 Evaluation of the ARI with different latent space sizes, comparing Kmeans clustering on 7, 8, 9 and 10 clusters.

Table 1 :

 1 Evaluation of the ARI and FScore of all clustering results obtained with the different methods.

	Raw data			Encoded data
	Kmeans	Kmeans	AggCl	GM	N2D	E multires	E struct
	Image 1 0.39 0.89 0.48 0.96 0.38 0.94 0.27 0.73 0.43 0.97 0.47 0.96 0.42 0.88
	Image 2 0.27 0.68 0.33 0.63 0.29 0.62 0.19 0.43 0.29 0.68 0.31 0.66 0.46 0.62
	Image 3 0.25 0.76 0.39 0.87 0.35 0.85 0.22 0.76 0.31 0.88 0.37 0.87 0.45 0.91
	Image 4 0.08 0.48 0.08 0.50 0.13 0.61 0.05 0.51 0.12 0.55 0.08 0.54 0.08 0.75
	Image 5 0.11 0.65 0.11 0.64 0.10 0.60 0.10 0.62 0.11 0.65 0.12 0.65 0.17 0.72
	Image 6 0.37 0.68 0.52 0.75 0.51 0.72 0.49 0.67 0.43 0.76 0.51 0.77 0.57 0.75
	Image 7 0.28 0.68 0.35 0.73 0.33 0.75 0.14 0.61 0.37 0.74 0.41 0.76 0.36 0.84
	Image 8 0.33 0.63 0.44 0.71 0.42 0.70 0.07 0.44 0.37 0.69 0.44 0.75 0.45 0.75
	Mean 0.26 0.68 0.34 0.72 0.31 0.72 0.19 0.60 0.30 0.74 0.34 0.75 0.37 0.78
	Stdv 0.11 0.12 0.16 0.14 0.14 0.12 0.14 0.13 0.13 0.13 0.16 0.13 0.16 0.09

Acknowledgement

This work was supported by the AiCOLO project funded by INSERM/Plan Cancer (MIC 2019).