Stefan Werner Hieke
email: s.hieke@mpie.de

Marc-Georg Willinger
email: willinger@fhi-berlin.mpg.de

Zhu-Jun Wang

Gunther Richter
email: richter@is.mpg.de

Dominique Chatain
email: chatain@cinam.univ-mrs.fr

Gerhard Dehm
email: g.dehm@mpie.de

Christina Scheu
email: c.scheu@mpie.de

On pinning-depinning and microkink-flow in solid state dewetting: Insights by in-situ ESEM on Al thin films

Keywords: In situ, Solid state dewetting, Thin films, Capillarity phenomena, Interface diffusion

The dynamics of solid state dewetting phenomena of a 50 nm thick, mazed bicrystalline Al film on single crystalline α-Al 2 O 3 (sapphire) substrates was studied in-situ using an environmental scanning electron microscope (ESEM). The bicrystalline Al thin films served as a model system where the influence of grain boundaries and texture effects are well determined compared to polycrystalline films. The experiments were performed in controlled oxidizing and reducing atmospheres at 773 K and 823 K, respectively, to shed light on the differences in dewetting mechanisms and dynamics. While the reducing atmosphere led to spontaneous dewetting at 823 K after an incubation time of a few minutes, a hierarchical dewetting process was observed for the sluggish dewetting under oxidizing conditions. Voids initiated at (substrate or surface) defects and expanded trying to maintain a hexagonal shape. Pinning and depinning processes led to a discontinuous void growth and irregular void shapes including finger instabilities. As a consequence, the void growth followed a variety of power law exponents between 0.10 and 0.55. A new microkink-flow mechanism was discovered at the terminating Al planes at the void.

Introduction

Interpretation of dynamics on post-mortem samples is challenging due to the loss of time resolution which may obscure the actual interplay of the dominating mechanisms [START_REF] Jacquet | Grain growth: The key to understand solid-state dewetting of silver thin films[END_REF][START_REF] Niekiel | Texture evolution and microstructural changes during solid-state dewetting: A correlative study by complementary in situ TEM techniques[END_REF].

In-situ techniques allow the direct observation of time-dependent phenomena, e.g. layer growth, grain growth, oxidation or phase transformation and reveal the details of the involved mechanisms [START_REF] Jacquet | Grain growth: The key to understand solid-state dewetting of silver thin films[END_REF][START_REF] Dehm | -situ Electron Microscopy: Applications in Physics[END_REF][START_REF] Oh | Oscillatory Mass Transport in Vapor-Liquid-Solid Growth of Sapphire Nanowires[END_REF][START_REF] Wang | Direct Observation of Graphene Growth and Associated Copper Substrate Dynamics by in Situ Scanning Electron Microscopy[END_REF]. In the present work in-situ environmental scanning electron microscopy (ESEM) was used to investigate the dynamics of the solid state dewetting behavior of mazed bicrystalline Al thin films on sapphire under well-defined atmosphere and temperature conditions.

Solid state dewetting is usually controlled by surface diffusion and capping layers can be added on thin films in order to suppress dewetting [START_REF] Thompson | Solid-State Dewetting of Thin Films[END_REF][START_REF] Cao | Preventing Thin Film Dewetting via Graphene Capping[END_REF]. However, material transport via short-circuit pathways below the passivating layer can occur [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF] and lead to nanowire growth [START_REF] Kosinova | Whiskers growth in thin passivated Au films[END_REF]. Void formation despite the presence of capping layers is known as a degradation phenomenon for Al conductor lines on SiO 2 passivated Si [START_REF] Yost | Voiding due to thermal stress in narrow conductor lines[END_REF]. Thermally induced stresses were identified as driving forces for dewetting with grain boundary diffusion as the major diffusion path [START_REF] Yost | Stress-driven diffusive voiding of aluminum conductor lines[END_REF][START_REF] Kristensen | Hole formation in thin aluminium films under controlled variation of strain and temperature[END_REF]. In a recent work, we have studied post-mortem the formation of annealing-induced voids within <111> textured (several grain orientations) 250 nm thick Al thin films on single crystalline sapphire [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. Voids formed via retraction of the Al thin film underneath the continuous oxide layer. Grain boundary and interface diffusion were the dominating diffusion pathways while surface diffusion was suppressed [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. Irregular void shapes indicated a discontinuous void growth mechanism [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. Temperature induced void formation of Al thin films covered with a capping layer was observed in literature for Al on passivated Si [START_REF] Kristensen | Hole formation in thin aluminium films under controlled variation of strain and temperature[END_REF][START_REF] Korhonen | Mechanisms of thermal stress relaxation and stress-induced voiding in narrow aluminum-based metallizations[END_REF] and Al on single crystalline sapphire [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Dutta | Sub-Surface Oxidation at the Aluminum-Sapphire Interface During Low-Temperature Annealing[END_REF][START_REF] Dutta | Solid State Annealing Behavior of Aluminum Thin Films on Sapphire[END_REF]. All mentioned studies analyzed the solid state annealing behavior indirectly post-mortem without direct observation of the involved processes, except for a study of Sugano et al. [START_REF] Sugano | In-situ observation and formation mechanism of aluminum voiding[END_REF]. They observed the annealing induced void formation in Al thin films by optical microscopy in-situ on the µm level [START_REF] Sugano | In-situ observation and formation mechanism of aluminum voiding[END_REF].

Solid state dewetting of thin films without capping layers was observed by various insitu techniques including electrical resistance measurements [START_REF] Saxena | Kinetics of voiding and agglomeration of copper nanolayers on silica[END_REF], Auger electron spectroscopy [START_REF] Beszeda | Morphological evolution of thin gold films studied by Auger electron spectroscopy in beading conditions[END_REF], low-energy electron microscopy [START_REF] Cheynis | Dewetting dynamics of silicon-on-insulator thin films[END_REF], laser light transmission [START_REF] Jiran | Capillary instabilities in thin films[END_REF][START_REF] Jiran | Capillary instabilities in thin, continuous films[END_REF], confocal laser microscopy [START_REF] Schaefer | Stability and dewetting kinetics of thin gold films on Ti, TiO x and ZnO adhesion layers[END_REF][START_REF] Jahangir | In-situ investigation of thermal instabilities and solid state dewetting in polycrystalline platinum thin films via confocal laser microscopy[END_REF], atomic force microscopy [START_REF] Mizsei | In situ AFM, XRD and Resistivity Studies of the Agglomeration of Sputtered Silver Nanolayers[END_REF], X-ray diffraction [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Müller | Microstructure evolution during dewetting in thin Au films[END_REF][START_REF] Amram | Phase Transformations in Au-Fe Particles and Thin Films: Size Effects at the Micro-and Nano-scales[END_REF], time resolved differential reflectometry [START_REF] Luber | Solid-state dewetting mechanisms of ultrathin Ni films revealed by combining in situ time resolved differential reflectometry monitoring and atomic force microscopy[END_REF], scanning electron microscopy (SEM) [START_REF] Jacquet | Grain growth: The key to understand solid-state dewetting of silver thin films[END_REF][START_REF] Jacquet | On the solid-state dewetting of polycrystalline thin films: Capillary versus grain growth approach[END_REF] and transmission electron microscopy (TEM) [START_REF] Niekiel | Texture evolution and microstructural changes during solid-state dewetting: A correlative study by complementary in situ TEM techniques[END_REF][START_REF] Jang | Crystallographic analysis of the solid-state dewetting of polycrystalline gold film using automated indexing in a transmission electron microscope[END_REF][START_REF] Niekiel | The process of solid-state dewetting of Au thin films studied by in situ scanning transmission electron microscopy[END_REF]. High spatial resolution is only achieved by SEM or TEM, but for these methods the observed processes are usually limited to (ultra)high vacuum conditions. This work reports direct observations of the dynamics of the solid state dewetting behavior of bicrystalline, nominal 50 nm thick Al films on sapphire in oxidizing atmosphere using an ESEM. The Al thin film possesses a heteroepitaxial orientation relationship with the single crystalline substrate with the close packed planes and directions parallel to each other, which will be termed OR a/b in this article: OR a/b {111} || (0001) (± indicates the symmetrically equivalent twin variants a and b obtained by a 180° rotation around the <111> surface normal). The mazed bicrystalline Al thin films serves as a defined model system to reduce the interference of texture evolution and different grain boundary types as observed for oligocrystalline or polycrystalline thin films [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Müller | Microstructure evolution during dewetting in thin Au films[END_REF]. The mechanisms of annealing induced void formation, growth and faceting as well as finger instabilities and their evolution with time were tracked in-situ and further analyzed. A microkink-flow mechanism is observed for the first time in-situ. Dewetting under reducing atmosphere was also investigated to emphasize the impact of the surface oxide on the dynamics.

Experimental

Thin film growth

Thin films were grown in a molecular beam epitaxy system (DCA instruments, Finland) with a base pressure of 4•10 -7 Pa. (0001) oriented α-Al 2 O 3 two inch wafers (miscut <0.25°, one side polished, thickness 330±25 µm, CrysTec Kristalltechnologie, Berlin, Germany) were used as substrates. Ar + -sputtering (200 eV, 5 min) was applied within the molecular beam epitaxy system to remove contaminations and structural defects were healed by subsequent annealing at 1273 K for 60 min. The nominal film thickness and deposition rate were 50 nm and 0.03 nm/s, respectively. During deposition, the sapphire substrate was rotated and kept at room temperature.

In-situ annealing and microstructural characterization

In-situ isothermal annealing experiments were performed in a field emission gun (FEG) ESEM FEI Quanta 200 equipped with a heating stage and a gas supply unit (mass flow controllers from Bronkhorst) [START_REF] Wang | Stacking sequence and interlayer coupling in few-layer graphene revealed by in situ imaging[END_REF]. Four short time (<1 h) and a single long-term (~15. Prior to the in-situ experiment and post-mortem, the samples were characterized by SEM and electron backscatter diffraction (EBSD) using FEG SEMs (Jeol JSM 6500F and Zeiss Merlin) or focused-ion-beam (FIB) systems (FEI Helios Nanolab 600 or 600i dual beam SEM/FIB). The microstructure of the as-deposited film was also characterized using X-ray diffraction (XRD) in a GE Seifert diffractometer equipped with a Cu X-ray source (K α = 1.54•10 -10 m). The grain size was estimated by the line intercept method with horizontal straight lines using the TSL OIM software applied on post-mortem EBSD data. This is an apparent grain size (mean value ± standard deviation) due to the mazed microstructure of interpenetrating twinned grains. For the as-deposited film the grain size was determined from cross-sectional TEM samples.

Statistical void size data was obtained from post-mortem plan-view SEM micrographs and is presented by a mean value ± standard deviation of the equivalent void radius (only voids with an area >10 µm² considered). Individual voids were tracked during the in-situ experiment. The perimeter and the equivalent void radius as a function of time were extracted from the in-situ SEM datasets by image processing using ImageJ [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF] including the plugin suite FeatureJ [START_REF] Meijering | FeatureJ: An ImageJ Plugin Suite for Image Feature Extraction[END_REF] and the plugin Template Matching and Slice Alignment [START_REF] Tseng | A new micropatterning method of soft substrates reveals that different tumorigenic signals can promote or reduce cell contraction levels[END_REF].

The velocities of facets and microkinks were determined by measuring the position with respect to a reference position (temporal resolution 17.7 s for the short-term dataset shown in Figure 2, 88.5 s for the long-term dataset shown in Figure 6).

The error of the displacement x is estimated with two times the pixel size. For facets, the normal was constructed for the analyzed facets for a reference frame and the position along the normal was tracked. Diffusion constants D for Al were derived by analyzing the mean square displacement x 2 of facets/microkinks as a function of time t ( . The measured values are apparent diffusion constants as discussed in section 0.

Site-specific cross-sectional specimens for TEM analysis were prepared by FIB using a modified procedure of the lift out technique reported in literature [START_REF] Giannuzzi | A review of focused ion beam milling techniques for TEM specimen preparation[END_REF][START_REF] Schaffer | Sample preparation for atomic-resolution STEM at low voltages by FIB[END_REF][START_REF] Langford | In situ lift-out: Steps to improve yield and a comparison with other FIB TEM sample preparation techniques[END_REF]. More details and the parameters of the applied method, briefly described here, can be found in [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. A Pt layer was deposited to protect the specimen surface from Ga ion milling. After the transfer of 1 -2 µm thick TEM lamellae to a TEM grid, the sample was thinned to electron transparency with accelerating voltages from 30 kV down to 5 kV to minimize beam damage. The final ion beam current at 5 kV was set to 15 pA.

TEM characterization of all films was carried out using a FEI Titan Themis 60-300 X-FEG S/TEM instrument operated at 300 kV, equipped with a probe corrector for spherical aberration. For scanning TEM (STEM) measurements a spot size of 0.15 nm and a convergence semi-angle of the electron probe of 23.8 mrad were used. Probe currents of ~70 pA were applied for imaging and up to 400 pA for energy-dispersive X-ray spectroscopy (EDS) measurements. EDS experiments, were carried out in STEM mode using a FEI Super-X windowless EDS system with 4 synchronized silicon drift detectors attached to the Titan Themis. A high angle annular dark field (HAADF) detector covering the angular range of collection semi-angles from 73 to 352 mrad was used for imaging. High signal to noise ratios for the atomically resolved STEM HAADF micrographs, were obtained by stacking of up to 60 frames using the Smart Align software [START_REF] Jones | Smart Align-a new tool for robust non-rigid registration of scanning microscope data[END_REF] 3. Results

In-situ microstructural evolution

Discontinuous void growth mechanism

The initial, dense, homogenous Al film possessed defects like scratches and surface contamination in the form of aggregates. They were visible at a few locations and were used as markers during the in-situ experiments. These side facets of the void which do not interfere with the scratch were tracked (see facets 2, 3 and 4 in Figure 2a) and their position changes are shown as colored lines in Figure 2a. The positions along the facet normal of a reference frame (t 0 + 283 s) are plotted in Figure 2b. The irregular void growth and shift of the centroid prevented the tracking of all facets from the first frames and the identification of a unique starting point. Thus, the relative distance of each track is shown, whereas the first detected point serves as the relative starting point. Pinning events and the subsequent depinning at later stages can clearly be observed for facets 3 and 4 ( Figure 2b). At the beginning of the void growth, a faster diffusion was observed. For facet 3, the shortest facet, a maximum value of ~8 •10 -15 was measured for the first frames. For the longer facets 2 and 4, the initial points were not tracked, but estimated from a straight line from the point where the void initiated to the facet position at t 0 + 283 s. The diffusion constant for them was in the range ~4•10 -15 . An average diffusion constant of ~1 -2•10 -

15

was found for all moving facets.

Void growth kinetics

The velocity of the retracting facets decreases as the hole grows. For example, in the case of void 1 ( Figure 1), the average velocity of the moving facets 2 -4 is in the range of ~0.5 -1 for the main growth regime. In the initial stages, facet velocities up to ~15 were observed depending on the pinning sites (e.g. Void 1 in Supplementary figure 1a). The equivalent void radius r void can be described by a power law function of the annealing time t: r void ~ t x . The power law fit applies for the main growth regime. Different power law exponents were extracted for different voids tracked in the same dataset: r void 1 ~ t 0.36 , r void 2 ~ t 0.26 and r void 3 ~ t 0.11 and they stay constant for short-term experiments (Supplementary figure 1). Some voids are already present when reaching the target temperature of 773 K indicating that they form during heating up (e.g. voids 5 and 6 in Figure 3). In addition, new voids form during the isothermal treatment after various incubation times (e.g. void 4). In the long-term experiment, the growth of individual voids slows down (void 4 and 5) or even stops. This can be seen in the change of the slope in the log-log-plot in Figure 3c. In summary, the power law exponent x of the individual voids varies from 0.10 and 0.55 depending on the pinning of the retracting film at obstacles.

Formation of microkinks and fingers

Varying boundary conditions during void growth causes the formation of complex void morphologies. Faceting instabilities or microfaceting, finger instabilities, and pinch-off processes occur as a result of pinning as shown in Figures 34567and Supplementary movie 2 for the long-term experiment (~15.3 h). Pinning and depinning of a retracting Al front during the dewetting process is shown for void 4 in the SEM image sequence (only selected frames)

and the perimeter evolution of the void is displayed in Figure 4 (see also Supplementary Figure 5 shows an unexpected microkink-flow diffusion mechanism along void facets resulting in film retraction and (micro)facets. Supplementary figure 2 shows the size of the microkinks into the void which is in the range from ~20 nm to 300 nm (max. ~800 nm) in post-mortem SEM plan-view micrographs. The position of the microkink along the facet was tracked as a function of time and is shown as overlay on the SEM micrograph in Figure 6a where the first and last frame, where the microkink was observed, were merged. Within the experiment the microkink moved back and forth, once each in this sequence, removing two "steps/rows" of the facet. The relative distance to the intersection of the facet with the adjacent one at t 0 + 351 min (reference position) as a function of time t is shown in Figure 6b. (Supplementary movie 2).

In-situ annealing in reducing atmosphere

In-situ annealing at ~823 K in reducing atmosphere (67% H 2 /33% N 2 , ~50 Pa) led to a different solid state dewetting behaviour compared to the one observed when annealing in an oxidizing atmosphere. After dwell times in the range of minutes once the desired temperature has been reached, the film breaks-up into isolated, µm-sized features within the Al film and regions with exposed substrate. The process is relatively fast: at 823 K it takes less than 35 s. Supplementary figure 3 shows the resulting, complex microstructure which is beyond the scope of this article. However, these experiments in reducing conditions confirm that the sluggish void formation, diffusion-limited dewetting and pinning induced discontinuous void growth, as observed in an oxidizing atmosphere, is due to the presence of a stable and continuous surface oxide discussed for the ex-situ annealing experiment [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF].

Comparison of the microstructure before and after annealing

Initial microstructure after deposition

Microstructure after in-situ annealing

Several in-situ annealed samples (O 2 , ~773 K) were investigated post-mortem. Postmortem EBSD analysis of short-term as well as of long-term in-situ annealed samples reveal no change in texture (Supplementary figure 4). After annealing, the Al twin growth variants OR a/b are still arranged in a mazed structure, as reported for some fcc metallic thin films on single crystalline (0001) α-Al 2 O 3 for certain growth conditions [START_REF] Dehm | Growth and microstructural stability of epitaxial Al films on (0001) α-Al 2 O 3 substrates[END_REF][START_REF] Amram | Grain boundary grooving in thin films revisited: The role of interface diffusion[END_REF][START_REF] Amram | On the role of Fe in the growth of single crystalline heteroepitaxial Au thin films on sapphire[END_REF][START_REF] Ramanathan | Microstructural study of epitaxial platinum and Permalloy/platinum films grown on (0001) sapphire[END_REF][START_REF] Radetic | Mechanism and dynamics of shrinking island grains in mazed bicrystal thin films of Au[END_REF]. Upon annealing the grain size increases to 1.7±1.4 µm (short-term, ~40 min) and 1.9±1.5 µm (long-term, ~15.3 h) which is a factor of ~10 compared to the as-deposited thin films. The stagnation of grain growth is in accordance to literature observations [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. Plan-view SEM micrographs after ~15.3 h in-situ annealing at 773 K reveal that the voids occupy ~3.9% of a ~381 x 286 µm large area (Supplementary figure 5).

For TEM analysis, selected regions of the short-term (~40 min) in-situ annealed sample containing voids as well as finger instabilities were prepared in cross-section by FIB (Figure 9a, Supplementary figure 6). Both, the Al film and the void, are capped by a surface oxide layer 5.5±0.4nm thick. The voids underneath a continuous surface oxide membrane appear dark in the SEM micrographs. After a collapse of the surface membrane, the oxide layer is in contact with the substrate as marked in Figure 9a and voids appear bright in plan-view SEM micrograph due to charging effects of the sapphire substrate. This was only rarely observed in-situ, but for the majority of voids post-mortem after the in-situ experiment. The collapse probably occurred after the in-situ experiment due to the small film thickness and the pressure difference between the microscope environment and atmosphere. Figure 9b shows the Al|void interface region without rim at the surface triple junction. The total thickness (metal plus oxide films) remains almost constant indicating that Al redistributes over large areas [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. The STEM micrograph in Figure 9c Al|substrate interface is observed as for the as-deposited thin film TEM sample (Figure 9e).

The high contrast at the interface is a FIB artefact as described above.

Discussion

Comparison of the solid state dewetting behavior in oxidizing and reducing conditions

This in-situ study of 50 nm passivated mazed bicrystalline Al films allowed the direct observation of the dewetting behavior, as the initiation and growth of voids in oxidizing conditions. During the in-situ annealing experiments in reducing atmosphere, we speculate that cracks form in the covering oxide membrane the Al film expands upon heating relative to the more rigid sapphire substrate. The lack of oxygen prevents healing of the cracks [START_REF] Laurent | Wettability of monocrystalline alumina by aluminium between its melting point and 1273 K[END_REF][START_REF] Levi | Oxygen induced interfacial phenomena during wetting of alumina by liquid aluminium[END_REF] and outward Al diffusion becomes possible. Surface diffusion is then the dominant diffusion mechanism for Al redistribution resulting in breaking-up of the film in µmsized features on the sapphire substrate (Supplementary figure 3). In contrast, for oxidizing conditions the surface oxide is stabilized, e.g. in O 2 atmosphere in the ESEM in this in-situ study or by O 2 or H 2 O as residual gases during ex-situ annealing studies [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Dutta | Solid State Annealing Behavior of Aluminum Thin Films on Sapphire[END_REF]. This is also indicated by an increase in thickness of the surface oxide while still continuously covering the film and the voids in accordance to our previous ex-situ annealing study [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF].

The diffusion processes involved in void growth have been discussed in our earlier study in depth [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. In brief, Al diffusion from the Al|void|sapphire triple junction towards the Al|void|surface oxide and subsequent outward Al diffusion at the Al|surface oxide interface was assumed as major Al diffusion pathways. Still, interfacial diffusion at the Al|sapphire interface cannot be excluded.

New voids formed after different incubation times constantly as also reported in literature [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Shaffir | Void formation in gold films on yttrium-doped zirconia in the initial stage of de-wetting[END_REF][START_REF] Shaffir | The mechanism of initial de-wetting and detachment of thin Au films on YSZ[END_REF][START_REF] Müller | Dewetting of Au and AuPt alloy films: A dewetting zone model[END_REF]. In the ESEM, void initiation at surface scratches was observed in-situ. However, voids were also observed to a minor degree within regions without visible surface defects (compare Supplementary movie 2). Void formation can also start at defects at substrate imperfections at the interface to the films, which could be e.g. contaminants or surface steps being present prior Al film deposition. 2D defects like Σ3 twin boundaries might also act as a starting point for void initiation. Müller and Spolenak reported that pores initially form at the film|substrate interface, grow preferentially along grain boundaries in Au thin films and evolve to voids [START_REF] Müller | Dewetting of Au and AuPt alloy films: A dewetting zone model[END_REF]. In addition to the void initiation, all mentioned defects are involved in pinning events determining the evolution of the void with time as described in the next section.

In-situ observed dynamics, pinning and microkink-flow mechanism

Interface and grain boundary diffusion are important diffusion pathways in the present sapphire substrate|Al film|surface oxide system as discussed in [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. Al self-diffusion constants at 773 K are reported in the range D ~4•10 -14 [START_REF] Lundy | Diffusion of Al 26 and Mn 54 in Aluminum[END_REF][START_REF] Mehrer | .13 Aluminum group metals: Datasheet from Landolt-Börnstein -Group III Condensed Matter[END_REF][START_REF] Brandes | Smithells Metals Reference Book[END_REF] which is higher than the measured ones (~1 -2•10 -15 for void 1 in Figure 1). The determined ones are only apparent diffusion constants as no (local) increase in film thickness was measurable. To calculate the average characteristic length L a time t of 40 min and the D volume = 4•10 -14 from literature [START_REF] Lundy | Diffusion of Al 26 and Mn 54 in Aluminum[END_REF][START_REF] Mehrer | .13 Aluminum group metals: Datasheet from Landolt-Börnstein -Group III Condensed Matter[END_REF][START_REF] Brandes | Smithells Metals Reference Book[END_REF] were inserted in and a value of ~13 µm was found for the short-time experiment ( Figure 1, Supplementary movie 1). This value is higher than the maximum change of position for the tracked Al|void interface of ~3 µm. Interface diffusion is supposed to be the dominating diffusion process [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF] resulting in even higher characteristic lengths L compared to volume diffusion. This indicates that Al could be redistributed to the surrounding remaining Al film. A simple calculation shows that an increase of 1 nm in film thickness for the Al film surrounding the void is sufficient to distribute all Al contained in the volume of void 1 (Supplementary figure 1b, additional radius: ~11µm). However, the obtained, apparent diffusion constant can be considered as lower-bound estimate and differ for the void growth just by one order of magnitude. For a more detailed model, the influence of the newly formed surfaces and interfaces, capillary forces and especially pinning events should be considered as addressed by Amram et al. [START_REF] Amram | Anisotropic hole growth during solid-state dewetting of single-crystal Au-Fe thin films[END_REF]. Driving forces and diffusion paths may locally vary as visible in the pinning events.

The in-situ observed void size evolution is represented by the change of the equivalent radius with annealing time t (see Supplementary figure 1c). It can be described by one power law r void ~ t x for the short-term experiments (~40 min) where r void is the equivalent void radius and the individual exponent x for each void varies in the range from 0.10 to 0.55.

The results of the long-time (~15.3 h) annealing experiment can only be explained by using two power law functions with two different exponents for one void which might be a result of a decreasing driving force and/or varying diffusion path. The void size evolution for ex-situ annealed 250 nm Al thin films containing two twinned ORs (rotated by 30° relative to each other) was described by a power law r void ~ t 0.23 [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF]. The different exponent of the ex-situ study can be caused by statistics (many voids, but poor time resolution for the ex-situ annealing study and vice versa for an in-situ study) as well as by other parameters such as film thickness, texture, grain size (evolution), and distribution of obstacles. Especially the film thickness is known to strongly influence the solid state dewetting kinetics [START_REF] Thompson | Solid-State Dewetting of Thin Films[END_REF].

This in-situ study reveals discontinuous void growth mechanisms due to pinning: void facet pinning, finger formation and microkink motion (see sketch in Figure 10). Pinning of the triple line can explain the irregular shaped voids observed in exsitu studies which deviate from the initial hexagonal void shape ( Figure 10a) [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF][START_REF] Dutta | Solid State Annealing Behavior of Aluminum Thin Films on Sapphire[END_REF]. Depinning can additionally occur as shown in Figure 4. The degree of complexity to describe the dewetting phenomena is increasing with annealing time due to varying boundary conditions affecting the void growth and also leading to additional microstructural features, e.g. fingers and microfacets.

Pinning induces the formation of Al fingers as schematically shown in

Figure 10b. While the front around the finger recedes, the finger gets longer and may even pinch-off as observed in-situ (Figure 7). Fingers are observed in most cases parallel to the <110> direction and confirm anisotropic film retraction [START_REF] Ye | Anisotropic edge retraction and hole growth during solid-state dewetting of single crystal nickel thin films[END_REF]. The formation of fingers and pinch-off during solid state dewetting is found for various systems in literature [START_REF] Thompson | Solid-State Dewetting of Thin Films[END_REF][START_REF] Jiran | Capillary instabilities in thin films[END_REF][START_REF] Bussmann | Dynamics of solid thinfilm dewetting in the silicon-on-insulator system[END_REF][START_REF] Zucker | A model for solid-state dewetting of a fully-faceted thin film[END_REF] and is commonly described to be caused by an unstable, retracting edge. A finger is a feature of very high curvature which instability increases as its length increases. The increasing tensile stress might cause the transverse contraction leading to a concave shape in contrast to the convex Al sidewall at the Al|void interface (Figure 9b,c). Freestanding fingers were reported as the origin of adjacent faceted particles after dewetting without capping layer [START_REF] Müller | Microstructure evolution during dewetting in thin Au films[END_REF]. Instead of pinch-off, depinning may occur as shown in Figure 4. Thus film retraction continues, the finger regresses and a facet evolves. Fingers can also completely disappear by Ostwald ripening as observed in-situ by Niekiel et al. [START_REF] Niekiel | The process of solid-state dewetting of Au thin films studied by in situ scanning transmission electron microscopy[END_REF]. Supplementary figure 2 shows faceting as well as finger coincidentally within the same void. In the previous ex-situ study [START_REF] Hieke | Microstructural evolution and solid state dewetting of epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF] annealing of a 250 nm thick Al film, faceting and film retraction adjacent to a pinned position were also observed, while almost no finger formation occurred.

With increasing annealing time faceting instabilities were observed. Faceting instabilities are described in literature as a result of anisotropic edge retraction [START_REF] Ye | Anisotropic edge retraction and hole growth during solid-state dewetting of single crystal nickel thin films[END_REF] After reaching an intersecting obstacle, a straight facet has been formed again, as shown in the SEM sequences in Figure 4 and Figure 5, until the microkink flow mechanism turns back. The initial stages of the microkinkflow mechanism seem to occur faster than the SEM frame time in the in-situ experiment. The obtained, apparent diffusion constants are one or two magnitudes lower compared to the facet movement. The reason for the different diffusion constants is not clear, but even at this time scale pinning events were observed. Different microkink heights (~20 -800 nm) were found, whereas for higher microkinks faceting could be resolved by SEM (Supplementary figure 2). The combination of a pinned microkink with one or several following microkinks pinned at the same obstacles results in an increase in microkink height and finally faceting instabilities (observed in-situ for the long-term dataset in Supplementary movie 2). Postmortem EBSD revealed that one void is usually surrounded by several twin boundaries which could act as pinning obstacle (Supplementary figure 4).

Conclusion

The solid state dewetting behavior of a self-passivated mazed bicrystalline Al thin film model system was investigated in-situ in an ESEM. In-situ annealing experiments in reducing atmosphere at ~823 K led to rapid solid state dewetting and the formation of µm-sized features. In contrast, in oxidizing atmosphere the surface oxide is stabilized, surface diffusion suppressed and dewetting slowed down, which takes place by formation and growth of voids In summary, the in-situ annealing ESEM experiments confirmed the proposed discontinuous void formation mechanism for surface-passivated Al thin films on sapphire [START_REF] Hieke | Annealing induced void formation in epitaxial Al thin films on sapphire (α-Al 2 O 3 )[END_REF] and resolved unexpected mechanisms and kinetics. This study shows the potential of in-situ annealing experiments in the ESEM under well-defined atmospheres and controlled temperature conditions for resolving time-dependent details of solid state dewetting.

  3 h) annealing experiments have been carried out (full datasets can be found in Supplementary movie 1 and Supplementary movie 2). The solid state dewetting experiments were performed in oxidizing atmosphere (O 2 , total pressure ~30 Pa) at 773 K and at 823 K in reducing atmosphere (67% H 2 /33% N 2 , total pressure ~50 Pa) atmosphere. All gases possessed a purity of 99.99%. The base pressure within the ESEM chamber prior annealing was below 5•10 -3 Pa. Temperatures were measured with a K-type thermocouple, that was mechanically clamped onto the sample. Secondary electron micrographs were acquired at an acceleration voltage of 5 kV using an Everhart-Thornley detector. Heating of the sample to the target temperatures of 773 K and 823 K as well as the selection of the region of interest for isothermal annealing experiments took up to 15 min. The areas displayed in this article are representative as controlled by additional micrographs recorded after annealing at different locations. Comparing irradiated to non-irradiated regions, no impact of the electron beam was found. The in-situ setup allowed fast cooling (estimated cooling time down to 473 K ~30 seconds) and hence quenching of the microstructure.

Figure 1 Figure 1 Figure 1 (

 111 Figure 1 Left: Void evolution with annealing time observed in-situ in the ESEM (~40 min, O 2 , ~773 K. A part of the micrographs is missing due to sample drift. The full region of interest can be seen in Supplementary figure 1 and the complete dataset is shown in the

Figure 2 a

 2 Figure 2 a) Overlay of the positions along the normal of the facet as a function of time on the reference frame at t 0 + 283 s. b) The relative distance as a function of time for the facets 2,3 and 4 reveals a discontinuous Al film retraction. Inset shows a zoom-in between 0 and 900 s.

Figure 3 a

 3 Figure 3 a) Microstructure of the region of interest for the long-term in-situ annealing experiments after reaching the final temperature ~773 K. b) Temperature induced microstructural changes after ~15.3 h annealing (O 2 , ~773 K, see text for details, the full dataset can be seen in the Supplementary movie 2). The evolution of different features with time is indicated by the white rectangles. c) Equivalent void radius evolution with annealing time for the voids 4 to 6 (Note: for void 4 the time of the initial void formation was set as t 0 , for void 5 and 6, the absolute time was used because the initial void formation was not observed).

movie 2 )

 2 Figure 4.

Figure 4 Figure 3 )

 43 Figure 4 Left: Sequence of SEM micrographs of void 4 (indicated in

For

  both motions, forward and backward, the first steps could not be resolved because of limits in spatial and/or time resolution. During the turning back, a jump of the microkink occurred (~620 nm in 885 s). The motion of the microkink is discontinuous on a local scale due to pinning. The average backward motion occurs faster ( ~0.16 ) compared to the forward motion ( ~0.06 ). During the jump, the extrapolated velocity ~0.70 is a magnitude higher. The analysis of the backward motion (jump excluded) revealed an average Al|void interface diffusion constant of ~2•10 -16 , whereas the apparent diffusion constant varied from ~4•10 -17 to 3•10 -16 on a local time scale. For the forward motion an average apparent diffusion constant of ~8•10 -17 with local variations from ~2•10 -17 to 3•10 - 16 was calculated.

Figure 5

 5 Figure 5 The solid state dewetting mechanism of a large void via a microkink-flow process along the Al sidewall facet is shown in the sequence of SEM micrographs of the marked

Figure 6 a

 6 Figure 6 a) Overlay of the positions of a microkink moving back and forth along the facet as a function of time on the merged, color-coded SEM frames at t 0 +351 min (Al|void interface in white) and t 0 +916 min (Al|void interface in intense blue). b) The relative distance as a function of time for both microkink motion reveal a discontinuous motion on a local scale affected by pinning and different local velocities. Note the different magnitude of time compared to Figure 2.

Figure 7 Figure 3 .

 73 Figure 7 The formation of an isolated finger is revealed in the sequence of SEM micrographs of the marked region of void 5 in

Figure

  Figure 8a and d show the twin structure with vertical grain boundaries and an average width of the grains in the range from 50 -200 nm for the as-deposited thin film. The {111} Al pole figure in

Figure 8b ,

 8b Figure 8b, obtained by XRD, revealed a ~1:1 ratio between the two twinned in-plane orientations OR a/b in accordance to literature [40]. The Al thin film is covered by an amorphous and continuous native, 3.6±0.4 nm thick Al 2 O 3 oxide layer (

Figure 8

 8 Figure 8 The initial microstructure of the as-deposited Al thin film. a) Cross-sectional STEM BF micrograph revealing the columnar grain structure of the Al thin film. b) {111} XRD pole figure with the two twin-related growth variants OR a/b marked. c) Combined Al-and O-EDS map showing a uniformly thick Al film covered by a continuous surface oxide layer. g) Atomic resolved STEM HAADF micrograph of a triple phase boundary at the Al|substrate interface with two twin-related Al grains with OR a respectively OR b taken in [ ] Al zone axis, [ ] Al zone axis and [ ] Al 2 O 3 zone axis.

  Figure 9d, confirm that the finger consists of pure Al. Neither interfacial pores nor voids are observed in the HRSTEM micrographs of regions with intact Al film. An atomically sharp

Figure 9

 9 Figure 9 Annealing induced microstructural changes revealed by post-mortem analysis (~40 min, O 2 , ~773 K, corresponding in-situ movie: Supplementary movie 1). a) Crosssectional SEM micrograph of the region marked in Supplementary figure 6. Voids, a finger and a region with a collapsed surface membrane are observed. Cross-sectional STEM HAADF micrographs showing b) the Al|void interface region with a surface membrane, continuously covering both, the Al thin film as well as the void and c) the marked finger. d) Combined EDS map of the finger region. The finger consists of pure Al covered with an intact Al oxide membrane. d) Atomic resolved STEM micrograph of a region adjacent to the Al|void interface confirming OR a/b. The twin variants cannot be distinguished in the present Al [ ] zone axis which is parallel to the [ ] Al 2 O 3 zone axis (Al and O are marked red and white, respectively).

Figure 10

 10 Figure 10 Sketch of pinning affected void growth mechanisms in plan-view: a) facet pinning, b) finger formation and c) microkink-flow (Al thin film in grey, void in white, pinned triple line in red, growth directions in white, further growth steps in dashed white).

  Figure 10c. Subsequently, Al transport occurs along the facets in <110> directions by a

  in the Al film underneath the oxide capping layer via interface diffusion. Faceted void formation and growth were observed at ~773 K in 30 Pa O 2 atmosphere and the bicrystalline {111} || (0001) 50 nm thick Al film kept its orientation. Void growth dynamics is influenced by the pinning events of the void|film|substrate triple line. It can be described by a power law r void ~ t x with an exponent x varying in the range of 0.10 to 0.55. The void growth dynamics is discontinuous due to pinning which leads to two mechanisms: discontinuous film retraction normal to facets and microkink-flow. The film retraction can be described by an apparent diffusion constant in the range of ≤ 2•10 -15 .Fingers may form at pinning sites. With time, depinning occurs and dewetting continues or the fingers pinch-off and form elongated islands as observed in some of the dewetted areas.Such events affect the growth of the void shapes since depinning requires sufficient incubation time. As a consequence of local pinning, the void shapes deviate from an ideal hexagon.Microkink-flow was observed as a new mechanism and occurs along the facets of the void causing a stepwise film retraction normal to the facet. Microkinks are related to faceting instabilities which are intermediate steps in cases where straight facets again evolve. They also result from the combination of multiple microkinks and pinning.

Acknowledgement

The authors thank Ilsemarie Lakemeyer for experimental support with the thin film growth.

Helpful discussions with Wayne Kaplan are highly appreciated.

Appendices

A1 Supplementary figures

Supplementary figure 1: a,b) The full region of interest observed in-situ containing three voids after selected annealing times (~40 min, O 2 , ~773 K. The sequence for void 1 is shown in Figure 1 and the full dataset can be seen in the Supplementary movie 1). c) Increase of the equivalent void radius with annealing time for the three analyzed voids 1-3 of the dataset and the corresponding power law fits.


E Giner

G David

A Scemama

J P Malrieu

A simple approach to the state-specific MR-CC using the intermediate Hamiltonian formalism

This paper presents a rigorous state-specific multi-reference coupled cluster formulation of the method first proposed by Meller et al. [J. Chem. Phys. 104, 4068 (1996)]. Guess values of the amplitudes of the single and double excitations (the T operator) on the top of the references are extracted from the knowledge of the coefficients of the Multi-Reference Singles and Doubles Configuration Interaction (MR-CISD) matrix. The multiple parentage problem is solved by scaling these amplitudes from the interaction between the references and the singles and doubles. Then one proceeds to a dressing of the MR-CISD matrix under the effect of the triples and quadruples, the coefficients of which are estimated from the action of exp( T ). This dressing follows the logic of the intermediate effective Hamiltonian formalism. The dressed MR-CISD matrix is diagonalized and the process is iterated to convergence. As a simplification, the coefficients of the triples and quadruples may in practice be calculated from the action of T 2 only, introducing 5th-order differences in the energies. The so-simplified method is tested on a series of benchmark systems from Complete Active Spaces (CASs) involving 2-6 active electrons up to bond breakings. The comparison with full configuration interaction results shows that the errors are of the order of a few millihartree, five times smaller than those of the CAS-CISD, and the deviation to strict separability is lower than 10 µ hartree. The method is totally uncontracted, parallelizable, and extremely flexible since it may be applied to selected MR and/or selected CISD. Some potential generalizations are briefly discussed.

I. INTRODUCTION

In the domain of molecular physics and quantum chemistry, the many-body problem is perfectly clear as long as it is formulated from a single reference. The perturbative expansion of the wave operator and its diagrammatic transcription offers a guide to understanding the relation between the multiplicative structure of the wave function and the additive structure of the energy. The linked cluster theorem 1 clarifies the questions of the size consistency and of the strict separability into fragments. The defects of truncated Configuration Interaction (CI) are well understood and algorithms have been proposed to respect approximately (CEPA-0, 2,3 CEPA-n [4][5][6] ) or strictly ((SC) 2 CI) 7 the cancellation of unlinked diagrams. The Coupled Cluster (CC) method [8][9][10] is definitely the most elegant formalism and can be considered as the standard treatment in its Coupled Cluster Singles and Doubles (CCSD) version, or in the CCSD(T) version, 11 which incorporates the fourth-order effect of the triply excited determinants. But all these approaches fail when one cannot expect that a single determinant will represent a reliable starting point to conveniently generate the wave function. This is precisely the situation in many domains. The excited states present an intrinsic multi-determinantal character, and frequently a multi-configurational character. a) Electronic mail: scemama@irsamc.ups-tlse.fr Such are the magnetic systems in their low energy states, and the treatment of chemical reactions, in which chemical bonds are broken, also requires to consider geometries for which a single determinant picture is not relevant. A generalized linked cluster theorem has been established by Brandow, 12 which gives a conceptual guide, but the conditions that must be fulfilled for its demonstration (Complete Active Space (CAS) as reference space, mono-electronic zero-order Hamiltonian) would lead to strongly divergent behaviors of the corresponding perturbative expansion in any realistic molecular problem. Practical computational tools have been proposed, most of them being state-specific. One may quote second order perturbation expansions based on determinants from selected references (Configuration Interaction Perturbing a multi-configurational zeroth-order wave function Selected Iteratively, CIPSI), 13,14 intermediate Hamiltonian dressing, 15 in the so-called shifted B k technique. 16 These methods are not strictly size-consistent, and the conditions to satisfy the strict separability of determinant-based expansions require to define sophisticated zero-order Hamiltonians. 17 Contracted perturbative expansions, which perturb the multideterminant zero-order wave function under the effect of linear combinations of outer-space determinants, have also been proposed. One may quote the CASPT2 method 18,19 which uses a mono-electronic zero-order Hamiltonian, faces intruder state problems, and is not size consistent; the NEVPT2 method [20][21][22] which uses a bi-electronic zero-order Hamiltonian (the Dyall's one 23 ) and is size consistent and intruder-state free; the method from Werner; 24 as well as the perturbation derived by Mukherjee et al. [25][26][27] from their multi-reference coupled cluster (MR-CC) formalism.

If the CASSCF wave function is considered as the counterpart of the single-determinant reference, the CAS-CISD is the counterpart of the CISD, with the same sizeinconsistence defect, and the research of MR-CEPA and MR-CC has been the subject of intense methodological researches for about 20 years, without evident success. The cancellation of all unlinked terms in the MR expansion (i.e., a MR-CEPA or MR-(SC) 2 -CI) formalism is not an easy task. 28,29 If one lets aside the MR-CC methods that attribute a specific role to a single reference, 30 a few state-specific strictly multi-reference CC methods have been proposed: one suggested by one of the authors and collaborators, 31 another one by Mukherjee and co-workers 32 which has been intensively tested, and a third one in a Brillouin Wigner context. [START_REF] Mášik | Advances in Quantum Chemistry[END_REF][START_REF] Hubač | [END_REF] We return here on the first proposal which was only tested on a single problem. We shall present an improved formulation of the method in Section II, then the principle of its implementation (Section III), followed by numerical illustrations of its accuracy (Section IV). Section V comments on the properties of the method and Section VI discusses the advantages of this formalism, its flexibility, and possible extensions.

II. METHOD

A. Principle

Let us call |I⟩ the reference determinants, the number of which will be called N. The reference space may be a CAS, but this is only compulsory if one wants to satisfy the strict-separability property. If not the method is applicable to incomplete model spaces as well. The projector on the model space is

P0 =  I |I⟩⟨I|. (1)
Let us consider a zero-order wave function restricted to the model space,

|Ψ m 0 ⟩ =  I c m I |I⟩. (2)
This function may be either the eigenfunction of P0 Ĥ P0 ,

P0 Ĥ P0 |Ψ m 0 ⟩ = E m 0 |Ψ m 0 ⟩, (3)
or the projection of the eigenvector of the CAS-CISD on the model space,

|Ψ m 0 ⟩ = P0 |Ψ m CAS-CISD ⟩. (4)
The CAS-CISD wave function is written as

|Ψ m CAS-CISD ⟩ =  I c m I |I⟩ +  i c m i |i⟩, (5)
where |i⟩ are the singles and doubles (the determinants of the CAS-CISD space which do not belong to the reference space). We want to follow the Jeziorski-Monkhorst 35 expression of the wave operator Ωm which is supposed to transform the zero-order wave function into the exact one,

Ωm |Ψ m 0 ⟩ = |Ψ m ⟩, (6)
as a sum of reference-dependent operators,

Ωm P0 =  I Ωm I |I⟩⟨I|. (7)
Each of the Ωm I 's will take an exponential form

Ωm I = exp( T m I ), (8)
and each operator T m I will be truncated to the single and double excitations, as one does in the CCSD formalism.

B. The multi-parentage problem and the extraction of guess values of the excitation amplitudes from the CAS-CISD eigenvector

One may easily recognize that there exist some degrees of freedom in the determination of the wave operators. In the single reference CCSD expansion, one searches for the amplitudes of the excitations sending from the reference Φ 0 to the singly and doubly excited determinants. One evaluates the amplitudes of the triples and quadruples as given by the action of T 2 on Φ 0 , and the eigenequation is projected on each of the singles and doubles. If the number of singles and doubles is n, one may write a set of n coupled quadratic equations on the amplitudes. But it may be convenient to guess a first evaluation of these amplitudes from the coefficients of the singles and doubles in the CISD matrix, which may be done in a unique manner. From these amplitudes one may obtain a guess of the coefficients of the triples and quadruples and it is convenient (ensuring for instance a better convergence than solving coupled biquadratic equation) to write the process as an iterative dressing of the CISD matrix, in the spirit of Intermediate Effective Hamiltonian formalism. 15 This is a procedure to solve the set of non-linear equations and the results are identical to those obtained with the traditional formulation of the CCSD method.

In the multi-reference context one faces a genealogical problem, sometimes called the multiple-parentage problem. Actually for a state-specific formalism, one has only one coefficient for each of the singly and doubly excited determinants |i⟩. In principle one may decide that this determinant is obtained from each of the references and one would then write

c m i =  I dm I i c m I , (9)
but one must find a criterion to define the N dIi amplitudes from the knowledge of a single coefficient. Returning to a perturbative estimate of the coefficients of the singles and doubles starting from Ψ m 0 , the first-order expression of these coefficients,

c m(1) i = ⟨Ψ m 0 | Ĥ |i⟩ E m 0 -⟨i| Ĥ |i⟩ =  I c m I ⟨I| Ĥ |i⟩ E m 0 -⟨i| Ĥ |i⟩ , ( 10
)
suggests that the amplitudes of the excitation operators from the references to the singles and doubles might satisfy dm

I i dm J i = ⟨I| Ĥ |i⟩ ⟨J| Ĥ |i⟩ . ( 11
)
This scaling had been proposed in Ref. 31. This condition may be expressed as

dm I i = λ m i ⟨I| Ĥ |i⟩, (12)
where the quantity λ m i is the inverse of an energy. Re-injecting this expression in Eq. ( 9) leads to

c m i = λ m i  I c m I ⟨I| Ĥ |i⟩, (13)
which defines λ m i as

λ m i = c m i ⟨Ψ m 0 | Ĥ |i⟩ . ( 14
)
These are the key equations which define guess values of the amplitudes of the excitations leading from the references to the singles and doubles. Notice that we only consider amplitudes for the excitations which correspond to physical interactions, and since Ĥ is at most bi-electronic, one only introduces single-and double-excitation operators. Finally, we can re-express the CAS-CISD wave function as

|Ψ m CAS-CISD ⟩ =  I c m I 1 +  i dm I i TIi |I⟩, (15)
where

TIi |I⟩ = |i⟩. (16)
In order to reach the coupled cluster formalism, it is necessary to subtract the contribution of the product of single amplitudes from the amplitudes dm I i , The single excitation amplitudes s m I i are simply given by s m I i = dm I i , and the double excitation amplitudes are obtained from Eq. (17),

c m i =  I dm I i c m I =  I ∈Singles(i) s m I i c m I +  I ∈Doubles(i) d m I i +  (k,l)∈I → i ±s m I k s m I l c m I , (17)
d m I i = dm I i -  (k,l)∈I → i ±s m I k s m I l . (18)

C. Coefficients of the triples and quadruples

One may then generate the triples and quadruples |α⟩. Among them only those which interact with the singles and doubles (i.e., which are generated by the action of Ĥ on the singles and doubles and which do not belong to the CAS-CISD space) have to be considered. One may find the references with which they present either 3 or 4 differences in the occupation numbers of the molecular orbitals (MOs). These reference determinants may be called the grand-parents of |α⟩. The comparison between |α⟩ and each of its grand-parents |I⟩ defines the excitation operator from |I⟩ to |α⟩ as a triple or quadruple excitation,

TIα |I⟩ = |α⟩, (19)
which may be expressed in second quantization as the product of 4 (or 3) creation operators and 4 (or 3) annihilation operators,

TIα = a † q a † p a † n a † l a e a f a g a h . (20)
The creations run on active and virtual MOs, the annihilations run on active and inactive occupied MOs but the number of inactive indices among the creation and/or among the particles must be equal to 3 or 4, otherwise the determinant would belong to the CAS-CISD space. Knowing the operator, it may be factorized as the product of all complementary double or single excitation operators in all possible manners (each excitation keeping untouched the M s value). Then we may write the contribution to the coefficient of |α⟩ differing by 3 orbitals from |I⟩ as

d m I α =  (k,l)∈(I → α) ±d m I k s m I l +  ( j,k,l)∈(I → α)
±s m I j s m I k s m I l (21) and similarly if they differ by 4 orbitals,

d m I α =  (k,l)∈(I → α) ±d m I k d m I l +  ( j,k,l)∈(I → α) ±d m I j s m I k s m I l +  (i, j,k,l)∈(I → α) ±s m I i s m I j s m I k s m I l . (22)
Finally, one may write the coefficient c m α as

c m α =  I d m I α c m I . (23)

D. Dressing of the CAS-CISD matrix

If one considers the eigenequation relative to ⟨i|,

⟨i| Ĥ |i⟩ -E m c m i +  I ⟨i| Ĥ |I⟩c m I +  j i ⟨i| Ĥ | j⟩c m j +  α ⟨i| Ĥ |α⟩c m α = 0, (24)
one may decompose the last term,

 α ⟨i| Ĥ |α⟩c m α =  α ⟨i| Ĥ |α⟩  I d m I α c m I =  I  α d m I α ⟨i| Ĥ |α⟩ c m I . ( 25
)
Introducing the quantities,

⟨i| ∆m |I⟩ =  α d m I α ⟨i| Ĥ |α⟩, (26)
one may write the eigenequation (24) as

⟨i| Ĥ |i⟩ -E m c m i +  I ( ⟨i| Ĥ |I⟩ + ⟨i| ∆m |I⟩ ) c m I +  j i ⟨i| Ĥ | j⟩c m j = 0 ( 27
)
which suggests to treat the effect of the triples and quadruples as a column dressing of the CAS-CISD matrix. A similar idea has been exploited in the single-reference CCSD context, which may be presented and managed as an iterative dressing of the column between the reference and the singles and doubles. 36 The coupled cluster dressed CAS-CISD Hamiltonian may be written as PCAS-CISD ( Ĥ + ∆m ) PCAS-CISD , which is non-Hermitian. Defining the projector on the singles and doubles as

PSD = PCAS-CISD -P0 , ( 28
) PCAS-CISD ∆m PCAS-CISD = PSD ∆m P0 , (29)
one may define an equivalent Hermitian dressing ∆m′ in the case where one considers the Hermitization of the dressed CAS-CISD matrix to be desirable,

⟨i| ∆m′ |I⟩ = ⟨I| ∆m′ |i⟩ = ⟨i| ∆m |I⟩ (30)
provided that one introduces a diagonal dressing of the CAS-CISD matrix,

⟨I| ∆m′ |I⟩ = - 1 c m I  i ⟨I| ∆m′ |i⟩c m i . ( 31
)
The diagonalization of the matrices

PCAS-CISD ( Ĥ + ∆m ) PCAS-CISD and PCAS-CISD ( Ĥ + ∆m′ ) PCAS-CISD
will give the same desired eigenenergy and eigenvector

PCAS-CISD ( Ĥ + ∆m ) PCAS-CISD |Ψ m CC ⟩ = E m CC PCAS-CISD |Ψ m CC ⟩, ( 32
) PCAS-CISD ( Ĥ + ∆m′ ) PCAS-CISD |Ψ m CC ⟩ = E m CC PCAS-CISD |Ψ m CC ⟩. ( 33
)
Of course the process has to be iterated, and the resulting eigenvector defines new coefficients on both the references and the singles and doubles, which lead to new amplitudes, new evaluations of the coefficients of the triples and quadruples, and new dressings. Since the eigenvectors of the dressed matrices are identical, the two formulations, Hermitian or non-Hermitian, converge to the same solution. The converged solutions are the MR-CCSD energy and the MR-CCSD amplitudes, which define the exponential wave operator. One should keep in mind that the building of the dressed CAS-CISD matrix is not the only way to solve the non-linear equations of the MR-CC formalism, and more standard techniques like the Newton method or DIIS might be applied here.

E. A convenient approximation

The coupled cluster expansion can be related to a perturbation expansion, for instance, in terms of diagrammatic development of the wave function. Considering only the second power of T in the evaluation of the coefficients of the triples and quadruples ensures that these coefficients are correct to the second-order and that the energy is correct to the fourth-order. Hence, the second sum of ( 21) and the two last sums of ( 22) introduce only contributions of the third and fourth orders to the wave function, and of the fifth and sixth orders to the energy. Consistently, if these terms are neglected, one obtains d m I i = dm I i and a simplified expression of d m I α ,

d m I α =  (k,l)∈(I → α) ± dm I k dm I l . (34)
Splitting the excitation operators T into single and double excitation operators, T1 and T2 respectively, the terms . The price to pay is a certain violation of the strict separability, i.e., of the additivity of the energies of two subsystems A and B when treating the AB system at infinite distances. This will be illustrated in Section V B.

( T1 ) 2 , T1 T2

III. IMPLEMENTATION

The proposed algorithm was implemented in the Quantum Package, 37 an open-source series of programs developed in our laboratory. The bottleneck of this algorithm is the determinant comparisons needed to determine the excitation operators and phases during the reconstruction of the genealogy of the |α⟩'s. This was made possible, thanks to a very efficient implementation of Slater-Condon's rules. 38

A. General structure

At each iteration step, one first assigns the values of the λ m i parameters obtained from the eigenvector of the (dressed) CAS-CISD matrix according to Equation (14). From these parameters λ m i , the amplitudes of the single and double excitations are uniquely defined. Then one loops on the singles and doubles |i⟩. On each of them, one reapplies the excitation operators to generate the |α⟩'s. Those which belong to the CAS-CISD space are eliminated. The parents of |α⟩ (that are all the singles and doubles |k⟩'s such that ⟨α| Ĥ |k⟩ 0) are generated. If one of the |k⟩'s has already been considered in the loop on the |i⟩'s (k < i), this |α⟩ has been already generated and taken into account and must not be double counted. While generating the parents of |α⟩, its interactions with them, ⟨k| Ĥ |α⟩, are stored. At this step, the reference grand-parents |I⟩ are identified as having 3 or 4 differences with |α⟩. Then, the excitation operator leading from |I⟩ to |α⟩ is expressed in all possible manners as products of two complementary single or double excitations. For each couple of complementary excitations, the product of the amplitudes is accumulated to compute d m I α according to Equation (34). Finally, the product of ⟨k| Ĥ |α⟩ with d m I α is accumulated in ⟨k| ∆m |I⟩ for each parent |k⟩ of |α⟩ according to Equation (26). Once the loop on the |i⟩'s is done, all the |α⟩'s have been generated, and the column dressing is completed. Then, in order to fit with a symmetric diagonalization technique, one symmetrizes the dressing as mentioned in Sec. II (Equations ( 30) and ( 31)). The dressed CAS-CISD matrix is diagonalized and the process is repeated up to convergence of the calculated dressed energy. From the computational point of view, this process requires the storing of the dressing columns which scales as N × (n occ n virt ) 2 where N is the number of determinants in the reference, and n occ and n virt are, respectively, the number of occupied and virtual MOs. This amount of memory is reasonable and does not represent a bottleneck for the present applications. Regarding the CPU time, the costly part concerns the handling of the |α⟩'s which scales as (n occ n virt ) 4 . Nevertheless, the process is perfectly parallelizable as all the work done with the |α⟩'s generated from |i⟩ does not depend on any other |α⟩.

B. Practical issues

The definition of λ m i can lead to numerical instabilities when ⟨Ψ m 0 | Ĥ |i⟩ is small. Nevertheless, in such cases the contribution of |i⟩ to the post-CAS correlation energy is also small, suggesting that one might use a perturbative estimate of λ m i . In practice, we use the perturbative λ i according to two different criteria. The first one concerns the ratio of the variational coefficient c m i (obtained at a given iteration) over its perturbative estimate (see (10)). If = 1

E m 0 -⟨i| Ĥ |i⟩ . ( 35
)
In such situations, the coefficient c m i is not determined by its interaction with the reference determinants, but comes from higher-order effects. The second criterion concerns the absolute value of each of the d m I i defined according to (12). If any of these terms calculated with the λ m i obtained from the variational calculation (see (14)) is larger than 0.5, the perturbative λ . This condition avoids numerical instabilities occurring when both c m i and ⟨ψ m 0 | Ĥ |i⟩ are small and allows us the control of the maximum value of the amplitudes. As soon as along the iterations one of the |i⟩'s fulfills one of these criteria, it will be treated perturbatively in the following iterations. This precaution avoids significant oscillations due to back and forth movements from perturbative to variational treatment of the λ m i . The numerically observed residual oscillations are of the order of magnitude of 10 -6 E h , which may certainly be attributed to the non linear character of the numerical algorithm. Nevertheless, the order of magnitude of the residual oscillations is much smaller than the chemical and even spectroscopic accuracy. Similar instability problems have been noticed in the Mk-MRCC treatment. 39,40

IV. NUMERICAL TEST STUDIES

We decided to test the accuracy and robustness of the method on a series of benchmarks, some of which have been used in the evaluations of other MR-CC proposals and of alternative MR approaches. They essentially concern model problems, especially bond breaking problems or the treatment of degenerate situations. They require to use a CAS with up to six electrons in six MOs. In all cases, the method converged in a few iterations. A systematic comparison is made with full configuration interaction (FCI) estimates, either taken from the literature or obtained from a CIPSI type variation+perturbation calculation 13,14 where the perturbative residue is about -6 mE h . Of course, the CAS-CISD is already a rather sophisticated treatment, which takes into account, although in a size-inconsistent manner, the leading correlation effects, both the non-dynamical part in the CAS and the dynamical part in the CISD step. One may expect that the improvement brought by the MR-CC treatment will be significant when the number of important inactive double excitations is large.

In order to have a global view of the performance of the here-proposed algorithm, we report potential energy curves and the error to FCI estimate of our MR-CCSD algorithm together with the CAS-CISD. Tables showing the error with respect to the FCI estimate of the MR-CCSD and CAS-CISD are also presented, complemented by the total energies of the FCI estimate. The non-parallelism error (NPE) is here calculated as the difference between the minimum and maximum error to the FCI estimate. The spectroscopic constants are obtained from an accurate fit of the obtained potential energy curves with a generalized Morse potential representation. The spectroscopic constants reported here are the equilibrium distance R eq in Å, the frequency k eq in E h /Å 2 , and the atomization energy D e in kcal/mol. CAS-SCF orbitals were obtained with the GAMESS(US) program 41 and all Multi-Reference Singles and Doubles Configuration Interaction (MR-CISD), CIPSI and MR-CCSD calculations were performed using the Quantum Package. 37

A. Single-bond breakings

The treatment of the breaking of a single bond in principle requires only a CAS-SCF zero-order treatment including two electrons in two MOs. We have considered three problems of that type.

Bond breaking of the F 2 molecule

The F 2 molecule is a paradigmatic molecule since it is a case where the dynamical correlation brings a crucial contribution to the bonding. Despite the closed shell character of the wave function in the equilibrium region, the single reference Hartree-Fock (HF) solution is unbound (by 18 kcal/mol) with respect to the restricted open shell HF solution of the fluorine atoms. The 2-electron in 2-MO CASSCF treatment binds the molecule by 18 kcal/mol, but the experimental binding energy is much larger (39 kcal/mol). Going to a full valence CASSCF (14 electrons in 10 MOs) does not bring any improvement. The role of the dynamical correlation has been extensively studied and may be seen as a dynamic response of the lone pair electrons to the fluctuation of the electric field created by the two electrons of the σ bond. 42,43 The concept of orbital breathing has been proposed to express the fact that the orbitals of the lone pairs tend to become more diffuse on the negative center and more contracted on the positive center in the ionic valence-bond (VB) components of the CAS. These dynamic relaxation processes can only take place if one uses non-minimal basis sets.

The calculations of F 2 were made in the cc-pVDZ basis set 44 keeping the 1s electrons frozen, and accurate FCI estimates are taken from the work of Bytautas et al. 45 Figure 1 shows an exponential convergence of the energy along the

The C-C bond breaking in ethane

This calculation is performed in the 6-31G basis set, keeping the 1s electrons frozen. The geometrical parameters are given in Table II. The potential energy curves are reported in Figure 4 and the errors with respect to the FCI estimate appear in Figure 5. These data show that the error with respect to the FCI energy (<3.5 mE h ) is greatly reduced, by a factor 6 in average, with respect to the CASSDCI values. According to Table III, the NPE goes from 2.01 mE h to 1.32 mE h for, respectively, the CAS-CISD and MR-CCSD approaches. Concerning the spectroscopic constants, the impact of the CC treatment is modest but goes in the right direction.

The rotation of the ethylene molecule around its C-C bond

This twisting breaks the π bond. The calculation is performed in the 6-31G basis set at the geometry given in Table II, keeping the 1s electrons frozen. The occupied MOs in the inactive space involve 10 electrons, and despite the modest size of the basis set, one may expect a significant size-consistence defect of the CAS-CISD results, since they miss the repeatability of inactive double excitations on the SD determinants. The potential energy curve along the angle of rotation is reported in Figure 6 and the error to the FCI estimate is reported in Figure 7. From these data, it appears that the global shape of the potential energy curve obtained using the CC treatment is more parallel to the FCI curve than using the CAS-CISD approach. From Table IV, one observes that the error with respect to the FCI estimate (about 3 mE h ) is reduced by a factor of 6 when going from CAS-CISD to MR-CCSD. Also, the NPE is also reduced from 1.6 mE h to 0.3 mE h .

B. Insertion of Be in H 2

The Be + H 2 system is a popular benchmark for MR-CC methods. [46][47][48][49][50][51][52][53][54] This chemical reaction is somewhat more complex than the single bond breaking and the use of CAS(2,2) is questionable but this reference space has been used here for comparison with the literature. The geometries

where the berylium atom is at the origin and the hydrogen atoms at the coordinate (x, 0, ±z). As in previous calculations, the 1s electrons are kept frozen. The two active orbitals result from the CASSCF(2,2) optimization starting with the HOMO and LUMO orbitals of the RHF calculation. The results are given in Table V. From this table, it appears that the MR-CCSD reduces both the absolute error to the FCI energy and the NPE by a factor between 2 and 3, with respect to the CAS-CISD calculations. The NPE with our proposal is of 1.297 mE h in the cc-pVDZ basis set. This value compares ), which falls to 1.355 mE h after a perturbative inclusion of the triple excitations, 53

C. Two-bond breakings

Three systems have been treated using a CAS with four electrons in 4 active MOs. Two of them concern the simultaneous breaking of two bonds.

Breaking of the C= =C double bond of ethylene

The dissociation of the ethylene molecule by breaking the double bond was studied in the 6-31G basis set, with the geometry given in Table II. We report the error (<3 mE h ) with respect to the FCI estimate in Figure 8. The corresponding values appear in Table VI. Again, the error to estimated FCI energy is reduced by a factor of 4, but the NPE is reduced only by 20% with the CC treatment.

Two-bond breaking in H 2 O

This is a rather well known test problem for MR-CC methods. The calculation is done with the cc-pVDZ basis set at five different geometries obtained from the equilibrium geometry (R e = 1.843 45 Å and ∠ HOH = 110.6 • ), in order to compare with the values of the literature. 40 The results appear in Table VII. The benefit of the MR-CCSD with respect to the CAS-CISD treatment is significant: the maximum error is 1.4 mE h , better than the 6.4 mE h given by the Mk-MR-CC treatment. This improvement may be due to the hereproposed treatment of the amplitudes responsible for potential divergences. The NPE goes from 2 mE h to 0.7 mE h when the CC treatment is applied.

D. Triple-bond breaking in N 2

The dissociation of N 2 is a popular benchmark for multireference methods. It requires a 6-electron in 6-MO CAS. Calculations were done with the cc-pVDZ basis set at the six inter-atomic distances taken from Ref. 40. The results are reported in Table VIII. As the number of inactive electrons is rather small, the CAS-CISD results are already quite good: the NPE is 1.088 mE h and the error is in the order of 8-9 mE h . The CC treatment reduces the mean error to ∼2.5 mE h and the NPE decreases to 0.45 mE h . These results are much better than those of the Mk-MR-CCSD method and of a quality similar to that obtained by the Single Reference-Multi-Reference Coupled Cluster (SR-MR-CC) method of Oliphant and Adamowicz. 30

V. PROPERTIES

A. Internal decontraction

The method is internally decontracted. The coefficients of the references as well as those of the singles and doubles change along the iterations. If the reference space is a valence CAS, treating the non-dynamical correlation effects, the method takes care of the impact of the dynamical correlation on the non-dynamical part. The phenomenon is especially important in magnetic systems where the dynamical charge polarization effects increase dramatically the weight of the ionic valence bond components, diminishing severely the effective energy of these components. 43 This effect is already present in the CAS-CISD calculation but the MR-CCSD treatment eliminates the size consistency defect and slightly improves the quality of the projection of the wave function on the CAS.

B. Size consistence and separability

The size consistence and strict separability properties are crucial issues of multi-reference methods. 57 The rigorous formulation of the method does not introduce any unlinked diagram and is therefore size-consistent. However, the fifthorder approximation introduced by simply considering the products of dm I k (see Section II E) breaks the strict separability. The strict separability requires that in the splitting into two subsystems A and B, the active and inactive MOs are localized on one of the two subsystems A or B. Actually, as occurs for the Mk-MR-CC formalism, the method is not invariant with respect to the unitary transformation of the MOs in their class (inactive occupied, active, inactive virtual). This dependence will be studied in a future work, but the error to FCI being small, we do not expect a strong dependence on the MO definition. As was shown in the study of bond breakings, the asymptotic size-consistency error (which is demonstrated to be zero when localized MOs are used) is negligible in a basis of symmetry-adapted MOs.

Numerical tests were performed in the cc-pVDZ basis set on a Li + F problem at R = 16 Å using a 2-electron in 2-orbital CAS-CISD and MR-CCSD. At the CAS-CISD level, the nonadditivity error is 6 × 10 -5 E h , and it is reduced by a factor of 10 at the MR-CCSD level. The Be + H 2 system was also studied for R = 20 Å using both a CAS-CISD with the two electrons of the Be atom in the 2s and 2p orbitals of Be, and also adding to the CAS the 2 electrons of H 2 in the bonding and anti-bonding MOs of the molecule. The CAS-CISD energies deviate from additivity by, respectively, 1.8 mE h and 0.9 mE h . The additivity is recovered by the MR-CCSD at a precision 5 × 10 -6 E h . The total energies are given in Table IX.

This small deviation is due to the approximation introduced in Section II E (d m I i ≈ dm I i ): in the case where the interaction H I i is zero, our approximation of the amplitudes, based on this matrix element, cancels the corresponding d m I i , while the amplitudes of the single excitations contained in the dm I i may be non-zero. This is the case at the separation limit into two free radicals A• and •B, and the single excitations s I, Ti A → r A I on A and s I, T j B → s B I on B may be non-zero (as will be the case for instance for the spin polarization of the cores), while ⟨I|H |a † r A a † sB a jB a i A I⟩ = 0, since it is a dispersion type double excitation. Our approximation omits the action of the products s I, Ti A → r A I × s I, T j B → s B I on the singles and doubles, the dominant part concerning intra-system excitations. In our opinion, the computational benefit of the neglect of some high powers of the single excitations legitimates this approximation in the MR-CC treatment.

C. Eigenfunction of Ŝ2

The here-proposed method does not provide an eigenfunction of Ŝ2 as we consider only the determinants that are connected by an application of Ĥ to the determinants belonging to the CAS. This treatment does not include higher excitations which will generate the full space associated with a given space part. Along all the performed calculations on singlet states, the order of magnitude of the expectation value of Ŝ2 calculated on both the CAS-CISD and the projected MR-CCSD wave functions never exceeded 10 -3 . A solution working with the same restricted space but providing a strict eigenfunction of Ŝ2 is presented in Sec. VI.

VI. PROSPECTS

A. Computational cost

The method is extremely flexible, either on the choice of the reference space and/or of the excitations from it. As the theory is determinant based, one can take advantage of this flexibility to realize a CIPSI-like selection of the dominant contributions of both the references and the single and double excitations. Further works will investigate the various possibilities such as the combination of MR-CC with perturbation theory in order to target more realistic systems.

B. Excited states

The method is applicable to excited states using several approaches. The formalism being state specific, the dressing technique of the CI matrix can be applied to any state dominated by the reference determinants, as long as a state following procedure is applied. For states belonging to the same symmetry, the resulting eigenvectors will not be strictly orthogonal but might be orthogonalized a posteriori. Another possibility consists of a state average procedure where the amplitudes are obtained from the values of the quantities λ i averaged over all desired eigenstates,

λ i =  m λ m i (c m i ) 2  m c m i 2 . ( 38
)
If one refers to the perturbative expression of the first order coefficients,

λ m i ≈ 1 E m 0 -⟨i| Ĥ |i⟩ , (39)
this approximation should be relevant when the states are close enough in energy.

A recent paper 58 has proposed a generalization of this approach to the simultaneous treatment of several states of the same symmetry. The basic ideas are the same, except for the fact that the extraction of the amplitudes is more complex.

The method requires to partition the reference space into a main and an intermediate model spaces, in the spirit of the intermediate Hamiltonian formalism. 15 This proposal will be tested in a further work.

C. Spin multiplicity conservation

As already mentioned the space of the singles and doubles is limited here to the determinants which interact with at least one of the references. If the references already have several open shells, this condition drastically reduces the size of the CASSD with respect to the space involving all the determinants of the configurations which interact with the references through at least one of their determinants. The price to pay is some minor deviation of the eigenvector from an Ŝ2 eigenvector. If one wants to ensure this property, one may follow one of two following routes. If the one calls i 1 ,i 2 , . . . ,i k , . . . the determinants which belong to the same space configuration and interact with at least one reference, they may receive various values of the effective inverse energy λ m ik . If one defines a unique value λ m i for all the determinants of the configuration, for instance, the mean value of the λ m ik 's,

λ m i =  k λ m ik (c m ik ) 2  k (c m ik ) 2 , (40)
the vector

|i⟩ = λ m i  k ⟨Ψ m 0 | Ĥ |i k ⟩|i k ⟩ (41)
is an eigenfunction of Ŝ2 . This is due to the fact that Ŝ2 does not change the space part of the determinants and that it commutes with Ĥ. Consequently, the excitation operator TIi

d m I i = λ m i  k ⟨I| Ĥ |i k ⟩ TIk (42)
commutes with Ŝ2 . The products of such operators, which have been used to generate the coefficients of the triples and quadruples, also commute with Ŝ2 , and the CC process should keep the exact spin multiplicity. The other way to ensure that property, which is by far more expensive, would consist of starting from the complete CASSD, including the determinants i ′ h which do not interact with the references. The eigenvector is then an eigenfunction of Ŝ2 . A parentage must be defined between the references and the determinants i ′ , possibly by considering the interactions between the determinants of type i ′ and the determinants of type i. But this strategy would introduce many-body (more than 2-body) operators, which would make the algorithm much more complex and costly.

D. Possible extensions

The method is extremely flexible. One may, for instance, define a subset of inactive MOs (occupied and virtual) which have the major contribution to the dynamical correlation, and divide the CISD subspace into two parts, the "internal" excitations involving only the MOs of the selected subset would define a selected fraction of the CAS-CISD, and the effect of the triples and quadruples issued from the action of the "internal" double excitations on the selected CAS-CISD would lead to the dressing of the interactions between the references and the selected singles and doubles, while the rest of the CAS-CISD matrix would be untouched. This would correct the major size-inconsistency defect of the CAS-CISD treatment and introduce the major fourth-order contributions. Combinations of MR-CC and perturbations will be considered as well.

In order to reduce the computation time, the purely inactive double excitations, which are the most numerous, may be treated in a specific manner. Let us consider a composite system A • • • B where the active space is localized on A, and where the interactions between A and B vanish, and suppose that one is interested in the spectrum of A, B being in its ground state |0 B ⟩. The model space determinants have the form

|I ′ A ⟩ = |I A 0 B ⟩.
The effect of the double excitations on B should be a simple shift of the energies. This is the weak separability condition. It is not satisfied at the CAS-CISD level, due to spurious normalization effects. The inactive double excitations Tl B on B are possible on all references |I A ⟩ and the coefficient of the resulting determinant should be

c Tl B I A = c l B C I A , (43)
where c l B is the coefficient relative to the double excitation Tl B on B. The correlation energy brought by the double excitation on B is

ϵ l B = ⟨ Tl B 0 B | Ĥ |0 B ⟩c l B , (44)
and this effect may be introduced as an energy shift on the diagonal energies of all the determinants of the CAS-CISD matrix, as occurs for the (SC) 2 dressing in the single reference context. At convergence this dressing ensures the weak separability property. We suggest to extend this remark to all the inactive double excitations whatever the interaction between the active and inactive orbitals. Considering the CAS-CISD, one may define a first evaluation of the contribution of any inactive double excitation as the averaged contributions of the inactive double excitation Tl on the various references,

ϵ l =  I ⟨ Tl I| Ĥ |I⟩c Tl I C I  I C I 2 . ( 45
)
Then the effect of all the inactive double excitations should be summed and introduced as a shift of the diagonal energies of all the SD determinants in the dressed CAS-CISD matrix. The inactive double excitations should no longer be considered in the generation of the triples and quadruples determinants |α⟩ in the MR-CC treatment. Another extension will concern the Difference Dedicated Configuration Interaction (DDCI), which is a standard method in the treatment of magnetic systems. 59 The method is based on a CAS and eliminates the purely inactive double excitations, which are the most numerous and do not contribute (at least to the second-order) to the vertical energy differences. Due to this restriction, the error of size consistence is reduced but it exists and may be non-negligible when the number of active electrons increases. It is perfectly possible to couple-clusterize the DDCI method by repeating the effect of the semi-active excitations only on the determinants appearing in the DDCI space.

VII. CONCLUSION

This work shows the relevance of a solution previously proposed to the problem of the multiple parentage faced by all Multi-reference treatments, as soon as the number of targeted vectors is lower than the number of references. The proposed MR-CC algorithm is simple. It only introduces two-body excitation operators and the number of amplitudes to be determined is reduced to the very minimum. It proceeds through an iterative dressing of the MR-CI matrix, formulated in terms of a standard eigenvalue equation. It is parallelizable in the most expensive step (the generation of the coefficients of the triples and quadruples). It is entirely decontracted and may be applied to excited states. A convenient simplification introducing a fifth-order deviation on the energy was introduced. For the list of benchmark studies we have performed, the results are extremely encouraging. The present version is state-specific but the principles of extension to a multi-root version have been formulated. 58 This work actually opens into several directions, which have to be explored in the future. The reduction of the computational cost might be done using several approximations involving the selection of the references and/or singles and doubles according to various criteria. Furthermore, the excited states may be treated either in the present state-specific formalism or using its multi-state generalization.

On a different perspective, the multiple parentage problem, which was faced here with the purpose of building a logically consistent computational tool to go in the direction of the exact solution, also concerns the building of rational valence-only effective Hamiltonians. In such an approach, the idea is to map the information coming from a sophisticated treatment into a minimal effective Hamiltonian, the parameters of which should be as physically meaningful as possible. 60 We believe that the solution we proposed to the multiple parentage problem offers a rational solution to this reduction of information. This remark illustrates the intrinsic link between the two main tasks of quantum chemistry, namely, the production of physically grounded interpretative models on one hand and the conception of rigorous computational tools.

where • s m I i are

  the amplitudes of the single excitation operators between the reference |I⟩ and |i⟩, • d m I i are the amplitudes of the double excitation operators between the reference |I⟩ and |i⟩, • (k, l) ∈ (I → i) denotes the couples (k, l) for which TIk TIl = ± TIi , the sign being governed by the permutation logics.

  0.5] then the amplitudes involving |i⟩ are determined using the perturbative λ m(pert) i defined as λ m(pert) i

  the working amplitudes d m I i are set to d m(pert) I i

FIG. 1 .FIG. 2 .

 12 FIG. 1. F 2 molecule at R = 1.45 Å, cc-pVDZ basis set. Convergence of the MR-CCSD energy along the iterations.

4 .FIG. 5 .

 45 FIG. 5. Ethane molecule, 6-31G basis set. Errors with respect to the FCI estimate as a function of the C-C distance.

7 .

 7 Ethylene molecule, 6-31G basis set. Errors with respect to the FCI estimate as a function of the rotation around the C-C bond.

TABLE I .

 I F 2 molecule, cc-pVDZ basis set. Total energies are given in E h , and the energy differences are given in mE h .

	R (Å)	E CAS-CISD -E FCI	E MR-CCSD -E FCI	FCI estimate
	1.14	19.223	3.726	-199.007 18
	1.20	19.495	3.823	-199.048 11
	1.30	19.825	4.102	-199.085 10
	1.36	19.829	4.045	-199.095 17
	1.411 93	19.721	3.920	-199.099 20
	1.50	19.518	3.830	-199.099 81
	1.60	19.094	3.466	-199.095 10
	1.80	18.038	2.843	-199.080 90
	2.0	16.850	2.034	-199.068 82
	2.2	16.280	1.783	-199.061 65
	2.40	16.225	1.936	-199.058 23
	2.80	16.055	1.794	-199.055 77
	8.00	16.241	1.893	-199.055 45
	R eq	1.466	1.465	1.460
	k eq	0.730	0.739	0.795
	D e	26.01	26.91	28.31
	NPE	3.774	2.319

MR-CC iterations. The here-reported calculation, performed in a medium size basis set, does not afford a sufficient flexibility to reach the experimental binding energy (the estimated FCI binding energy in this basis is D e = 28.3 kcal/mol). The potential energy curves and the error to FCI estimate are reported, respectively, in Figures

2 and 3

, and the estimated FCI values together with the error of the MR-CCSD and CAS-CISD calculations appear in Table

I

. The average error, lower than 5 mE h , is reduced by a factor close to 6, and the NPE is only reduced by 40% by the MR-CCSD calculations (3.7 mE h for the CAS-CISD and 2.3 mE h for the MR-CCSD).

TABLE II .

 II Geometries used for ethane and ethylene.

	Geometrical parameters	C 2 H 6	C 2 H 4
	C-H	1.103 Å	1.089 Å
	C-C	1.550 Å	1.335 Å
	H-C-C	111.2 •	120.0 •
	H-C-H	107.6 •	120.0 •
	H-C-C-H	180.0 •	180.0 •

TABLE III .

 III Ethane molecule, 6-31G basis set. The FCI estimate is the CIPSI calculation. Total energies are given in E h , and the energy differences are given in mE h .

				FIG. 6. Ethylene molecule, 6-31G basis set. Energy as a function of the
				rotation around the C-C bond.
	R (Å)	E CAS-CISD -E FCI	E MR-CCSD -E FCI	FCI estimate
	4.00	15.508	1.834	-79.253 166
	3.50	15.770	1.915	-79.256 574
	3.00	16.379	2.074	-79.268 617
	2.60	17.037	2.622	-79.293 972
	2.30	17.402	3.009	-79.326 999
	2.00	17.519	3.134	-79.370 376
	1.90	17.510	3.150	-79.385 598
	1.80	17.482	3.152	-79.399 969
	1.70	17.442	3.106	-79.412 107
	1.65	17.419	3.035	-79.416 695
	1.60	17.395	3.046	-79.419 813
	1.55	17.371	3.055	-79.420 987
	1.50	17.347	3.062	-79.419 613
	1.45	17.326	3.083	-79.414 941
	1.40	17.306	3.099	-79.406 030
	1.35	17.291	3.135	-79.391 701
	1.30	17.284	3.153	-79.370 480
	R eq	1.549	1.550	1.550
	k eq	1.018	1.017	1.015
	D e	104.52	104.99	105.75
	NPE	2.011	1.319

TABLE IV .

 IV Rotation of the ethylene molecule, 6-31G basis set. The FCI estimate is the CIPSI calculation. Total energies are given in E h , and the energy differences are given in mE h .

	Angles (deg)	E CAS-CISD -E FCI	E MR-CCSD -E FCI	FCI estimate
	0	13.255	2.935	-78.216 340
	10	13.132	2.935	-78.214 241
	20	13.196	2.938	-78.208 391
	30	13.331	2.955	-78.198 732
	40	13.513	2.991	-78.185 373
	50	13.750	3.035	-78.168 619
	60	14.043	3.120	-78.149 094
	70	14.368	3.212	-78.128 205
	80	14.631	3.258	-78.109 498
	85	14.694	3.237	-78.103 326
	90	14.605	3.227	-78.100 966
	NPE	1.562	0.323
	are given by the relation

z = 2.54 -0.46x (a.u.),

TABLE V .

 V BeH 2 molecule, cc-pVDZ basis set. Total energies are given in E h , and the energy differences are given in mE h .

				FIG. 8. Potential energy curves of the ethylene molecule, 6-31G basis set.
				Errors with respect to the FCI estimate as a function of the C-C distance.
				favourably with the values reported in the literature using a
				smaller basis set, namely, the Mk-MRCCSD (2.238 mE h	49,53
	X (a.u.)	E CAS-CISD -E FCI	E MR-CCSD -E FCI	FCI estimate
	0	1.212	0.470	-15.835 475
	0.2	1.202	0.471	-15.833 837
	0.4	1.204	0.468	-15.828 821
	0.6	1.217	0.433	-15.821 276
	0.8	1.245	0.458	-15.812 169
	1	1.291	0.481	-15.802 262
	1.2	1.358	0.518	-15.791 845
	1.4	1.450	0.438	-15.780 681
	1.6	1.567	0.615	-15.768 174
	1.8	1.679	0.655	-15.753 649
	2	1.648	0.761	-15.736 620
	2.2	1.807	0.837	-15.716 964
	2.4	2.107	1.032	-15.695 068
	2.6	2.538	1.464	-15.672 378
	2.8	3.333	1.767	-15.655 899
	3	4.093	1.276	-15.666 560
	3.2	3.880	1.141	-15.690 150
	3.4	3.546	0.958	-15.714 088
	3.6	2.623	0.983	-15.735 526
	3.8	2.492	0.786	-15.752 272
	4	2.681	0.569	-15.760 567
	NPE	2.891	1.297

  1.663 mE h of the exponential wave function Ansatz of Hanrath. 51

TABLE VI .

 VI Dissociation of the ethylene molecule, 6-31G basis set. The FCI estimate is the CIPSI calculation. Total energies are given in E h , and the energy differences are given in mE h .

	R (Å)	E CAS-CISD -E FCI	E MR-CCSD -E FCI	FCI estimate
	1.20	9.962	2.635	-78.179 508
	1.25	9.960	2.614	-78.200 805
	1.30	9.961	2.578	-78.212 507
	1.35	9.971	2.578	-78.216 869
	1.40	9.985	2.549	-78.215 666
	1.45	10.006	2.558	-78.210 313
	1.50	10.028	2.539	-78.201 918
	1.55	10.055	2.517	-78.191 365
	1.60	10.085	2.507	-78.179 345
	1.65	10.121	2.500	-78.166 404
	1.70	10.158	2.498	-78.152 957
	1.80	10.238	2.528	-78.125 775
	1.90	10.313	2.545	-78.099 562
	2.00	10.368	2.635	-78.075 294
	2.30	10.310	2.411	-78.017 944
	2.60	9.890	2.228	-77.983 952
	3.00	9.358	1.863	-77.964 220
	3.50	8.321	1.011	-77.956 421
	4.00	8.616	1.237	-77.955 127
	R eq	1.362	1.362	1.362
	k eq	2.043	2.039	2.042
	D e	163.48	163.54	164.47
	NPE	2.047	1.624

TABLE VII .

 VII Symmetric dissociation of the water molecule, cc-pVDZ basis set. The FCI total energy55 is given in E h , and the deviations to this reference are given in mE h . Comparison with Mukherjee's state specific MR-CC values (E Mk-MR-CCSD -E FCI ) obtained from Ref. 40.

	R (R e )	E CAS-CISD -E FCI	E Mk-MR-CCSD -E FCI	E MR-CCSD -E FCI	FCI
	1 R e	4.923	2.909	1.407	-76.241 860
	1.5 R e	4.674	4.817	1.248	-76.072 348
	2.0 R e	3.665	6.485	0.855	-75.951 665
	2.5 R e	3.097	5.672	0.763	-75.917 991
	3.0 R e	2.959	3.987	0.845	-75.911 946
	NPE	1.964	3.576	0.644

TABLE VIII .

 VIII Dissociation of the N 2 molecule, cc-pVDZ basis set. The FCI total energy 56 is given in E h , and the deviations to this reference are given in mE h . Comparison with Mukherjee's state specific MR-CC values (E Mk-MR-CCSD -E FCI ) and SR-MR-CCSD values (E SR-MR-CCSD -E FCI ) obtained from Ref. 40.

	R (a 0 )	E CAS-CISD -E FCI	E Mk-MR-CCSD -E FCI	E SR-MR-CCSD -E FCI	E MR-CCSD -E FCI	FCI
	2.118	8.320	6.916	1.515	2.403	-109.278 339
	2.4	8.774	8.785	1.528	2.472	-109.238 397
	2.7	9.163	11.071	1.618	2.589	-109.160 305
	3.0	9.408	13.506	1.789	2.671	-109.086 209
	3.6	9.205	18.413	2.187	2.485	-108.994 906
	4.2	8.485	18.950	2.246	2.220	-108.966 950
	NPE	1.088	12.033	0.731	0.451

TABLE IX .

 IX Total energies (cc-pVDZ) for the numerical check of the separability of Li • • •F and Be • • •H 2 . (The H-H distance was set to 1 Å to magnify the correlation on H 2 .)

		CAS-CISD	MR-CCSD
	Li	-7.433 463	-7.433 463
	F	-99.528 807	-99.531 870
	Li + F	-106.962 270	-106.965 334
	Li • • •F	-106.962 209	-106.965 342
	Be	-14.618 519	-14.618 536
	H 2	-1.140 073	-1.140 073
	Be + H 2	-15.758 592	-15.758 609
	Be • • •H 2 CAS(2,4)	-15.756 791	-15.758 604
	Be • • •H 2 CAS(4,6)	-15.757 693	-15.758 605