W Ludwig
email: ludwig@insa-lyon.fr

J-Y Buffie

S Savelli

P Cloetens

Study of the interaction of a short fatigue crack with grain boundaries in a cast Al alloy using X-ray microtomography

Keywords: zusätzliche Analyse der Kornorientierungen mit Hilfe der Ru¨ckstreuelektronenbeugung (EBSD), ermöglicht die Diskussion und Interpretation der Ergebnisse hinsichtlich verschiedener Rissausbreitungsmechanismen Short fatigue crack, Microtomography, Grain boundary wetting, Gallium, PACS:

Zusammenfassung

Mit Hilfe von Synchrotronstrahlung-Röntgenmikrotomographie gelingt es die Kornstruktur einer Aluminium Guss-Legierung in der unmittelbaren Nachbarschaft eines Ermüdungsrisses dreidimensional abzubilden. Die Abbildung der Kornstruktur beruht auf der selektiven Benetzung der Aluminium Korngrenzen mit fl¨ussigem Gallium, welches als Kontrastmittel verwendet wird. Die komplizierte Form des Ermüdungsrisses im Inneren der Probe, sowie die an der Probenoberfläche beobachteten Verzögerungen in der Rissausbreitung ko¨nnen mit der Kornstruktur

Synchrotron Radiation X-ray microtomography is used to visualize and analyze simultaneously the three-dimensional shape of crystallographic grains containing a short fatigue crack in a cast Al alloy. The visualization of the grains is based on the decoration of Al grain boundaries by liquid Ga which serves as a selective contrast agent. The intricate three-dimensional shape of the fatigue crack, as well as the crack stops observed on the sample surface, are correlated to the grain structure of the material. Complementary measurements of the grain orientation on the sample surface by electron backscattering diffraction (EBSD) allow us to discuss and interpret the observations in terms of possible crack propagation mechanisms.

Re ´sume Ĺa

micro-tomographie utilisant le rayonnement X synchrotron est utilisée pour visualiser et analyser simultanément la forme tridimensionnelle de grains cristallographiques contenant une fissure courte de fatigue dans un alliage d'aluminium de moulage. La visualisation des grains repose sur la décoration des joints de grains de l'alliage d'aluminium par le gallium liquide qui sert de marqueur sélectif. La forme compliquée de la fissure de fatigue au sein du matériau ainsi que les arrêts de la fissure observées à la surface du matériau apparaissent corrélés à la structure des grains sous la surface de l'échantillon. La determination de l'orientation cristallographique des grains par diffraction des électrons rétro-diffusés (EBSD) permet de discuter et interpreter ces re´sultats en terme de mécanismes de propagation des fissures.

Introduction

The Paris law which relates the growth rate da/dN of a fatigue crack to the stress intensity factor ⌬K is extensively used to predict the fatigue life of engineering components. However, in many engineering metallic alloys, it has been observed that small cracks (with a size typically smaller than a few grains [START_REF] Miller | The short crack problem[END_REF]), do not show a "regular" growth with ⌬K (as described by the Paris law) but some accelerations and retardations due to their interactions with microstructural features [START_REF] Miller | The behaviour of short fatigue cracks[END_REF]. If neglected, this so called short crack phenomenon can lead to a serious over estimation of the fatigue life [START_REF] Suresh | Fatigue of materials 1st ed[END_REF]. It is therefore a crucial problem to tackle. Among all the microstructural features which can alter the propagation of short cracks, grain boundaries have been observed to play a key role [START_REF] Miller | The behaviour of short fatigue cracks[END_REF][START_REF] Zurek | The effect of grain size on fatigue growth of short cracks[END_REF][START_REF] Lankford | The growth of small fatigue cracks in 7075-T6 aluminium[END_REF][START_REF] Lankford | The influence of microstructure on the growth of small fatigue cracks[END_REF][START_REF] Zhang | Measurement of plastic zones associated with small fatigue cracks by selected area channeling patterns[END_REF][START_REF] Turnbull | The effect of grain size on the fatigue of commercially pure aluminium[END_REF][START_REF] Turnbull | The effect of grain size on fatigue crack growth in aluminium magnesium alloy[END_REF]. However, the experimental study of the interaction of short fatigue cracks with grain boundaries is greatly complicated by the fact that it is a three dimensional problem. A few studies have shown that the abnormal growth of fatigue cracks as observed from the surface of a sample can be accounted for by the presence of some sub surface microstructural features [START_REF] Clement | Short crack behaviour in nodular cast iron[END_REF][START_REF] Tokaji | Limitations of linear elastic fracture mechanics in respect of small fatigue cracks and microstructure[END_REF]. Those studies are based on serial polishing and optical observations of fatigue samples, a quite tedious and destructive experimental technique which strongly limits the number of observations and which is not free of artifacts (e.g. crack blurring in ductile materials). In the last five years, however, the availability of new powerful synchrotron X-ray sources has opened the way for three dimensional (3D) imaging of the interior of materials through the use of high resolution microtomography. This technique which gives direct images of internal defects in the micrometer range has been used, for example, to study fatigue crack closure in Al-Li alloys [START_REF] Guvelenir | Direct observation of crack opening as a function of applied load in the interior of a notched tensile sample of Al-Li 2090[END_REF] but also damage development in various heterogeneous materials by repeated inspection of the specimens at different deformation stages [START_REF] Buffie | Characterisation of internal damage in a MMC p using X-ray synchrotron phase contrast microtomography[END_REF][START_REF] Savelli | Analyse des me ´canismes et mode ´lisation de la dure ´e de vie en fatigue d'alliages de moulage AS7G03[END_REF][START_REF] Martin | Characterisation by X-ray micro-tomography of cavity coalescence during superplastic deformation[END_REF][START_REF] Babout | Characterization by X-ray computed tomography of decohesion, porosity growth and coalescence in model metal matrix composites[END_REF][START_REF] Owen | A synchrotron X-ray study of a TiSiC f composite during in situ straining[END_REF][START_REF] Everett | Spatial distri-bution of voids in HY-100 steel by X-ray tomography[END_REF]. The use of coherent X-rays has also been shown to improve the detection of cracks and microstructural features in strained metal matrix composites [START_REF] Cloetens | Observation of microstructure and damage in materials by phase sensitive radiography[END_REF].

In this work, high resolution microtomography has been used to visualize simultaneously the 3D shape of a short fatigue crack in a cast aluminium alloy and of the surrounding grains. The grain visualisation is achieved by applying liquid Ga to the sample surface. Ga penetrates along the grain boundaries into the bulk of the material and leads to the formation of microscopic grain boundary wetting layers which can be detected by X-ray absorption imaging [START_REF] Ludwig | Penetration of liquid gallium into the grain boundaries of aluminium: a synchrotron radiation microtomographic investigation[END_REF]. The 3D images obtained with this technique clearly show that the surface growth of short cracks is strongly correlated to sub surface microstructural features. Besides, the analysis of the crack shape with respect to the local grain boundaries also reveals the importance of the crack-plane tilt and twist for the propagation rate of stage I cracks as suggested by other authors [START_REF] Zhai | A crystallographic mechanism for fatigue crack propagation through grain boundaries[END_REF].

Experimental methods

A model AS7G03 cast alloy containing artificial pores has been used in this study. A complete description of the microstructure and fatigue deformation mechanisms of this material can be found elsewhere [START_REF] Buffie | Experimental study of porosity and its relation to fatigue mechanisms of model Al-Si7-MgO.3 cast Al alloys[END_REF][START_REF] Savelli | Analyse des me ´canismes et mode ´lisation de la dure ´e de vie en fatigue d'alliages de moulage AS7G03[END_REF], only the main details are given in what follows. The alloy, which chemical composition is given in Table 1, was chill mould cast in a silicon carbide mould. The artificial pores were obtained by introducing at 1093 K a mixture of hydrogen (H 2 ) and argon (Ar) gases in the liquid metal through a propeller. After cooling, the alloy was solution treated during 10 h at 813 K, quenched in cold water, and aged during 6 h at 433 K (T6 heat treatment). The resulting microstructure is shown in Fig. 2a andb. It consists in large grains with an average size of 300 µm (2D measurements), containing a mixture of primary α aluminium dendrites separated by the eutectic phase containing small silicon particles with an average size of 3 µm. The fatigue behaviour of the material has been studied during constant stress amplitude tests (120Յ⌬ σ Յ 240 MPa) at room temperature in air with a load ratio R ϭ 0.1 and a frequency of 10 Hz.

One fatigue sample has been prepared for tomographic investigation according to the sequence described below:

1. An interrupted fatigue test with a stress amplitude ⌬σ ϭ 160 MPa was carried out and the surface growth of a crack initiated on one sample's corner was monitored each 25.000 cycles during 375.000 cycles by optical microscopy. 2. After the test, the grain orientation at the surface around the crack was determined by the use of electron backscattering diffraction (EBSD, surface XZ defined in Fig. 1). 3. One needle shape tomography sample with a square cross section of 0.7 * 0.7 mm 2 and a length of 5 mm was carefully cut with a wire saw around the fatigue crack (see Fig. 1). 4. A first tomographic scan was performed in order to study the 3D shape of the crack in the bulk of the sample. 5. The sample was then exposed to liquid gallium and annealed for 1 h at 100 °C. 6. A second tomographic scan, revealing the grain boundary structure was recorded. 7. The (now fragile) sample has been embedded into epoxy resin, cut along the XY plane and polished in order to determine the crystallographic orientation of sub-surface grains via EBSD measurements.

Microtomographic imaging

The tomographic imaging experiments were carried out at the ID19 beamline of the European Synchrotron Radiation Facility (ESRF), Grenoble. The dedicated microtomographic set-up consists of a precision mechanic sample stage combined with a fast, high resolution detector system. The detector system itself is based on a 26 µm thick transparent fluorescent screen (LAG: Eu) [START_REF] Koch | X-ray camera for computed microtomography of biological samples with micrometer resolution using Lu 3 Al 5 O 12 and Y 3 Al 5 O 12 scintillators[END_REF] which transforms the X-rays into visible light and microscope optics to project the image on the Peltier cooled 1024 2 CCD camera. The latter has a dynamic range of 13 bits, fast readout (0.06 s/frame) and a low dark current (3 e-/s) [START_REF] Labiche | FRELON camera: Fast REadout LOw Noise[END_REF]. The spatial resolution was determined with the knifeedge method to be 1.7 µm (full width half maximum of the line spread function) for the optical magnification used in our experiments (effective pixel size: 0.95 µm). The need for ultimate spatial resolution in order to resolve cracks with micron and sub-micron widths restricts the field of view and consequently the sample cross-section 2 to approximately 1 mm. In order to obtain high quality 3D images the incoming 'white' synchrotron radiation was monochromatized to 15 keV using a synthetic Rb-B 4 C multilayer with large energy bandwidth (⌬λ / λ ϭ 2 × 10 Ϫ 2 ) [START_REF] Morawe | Design and performance of graded multilayers as focusing elements for X-ray optics[END_REF].

Typical scan times (1000 projections over 180 degrees, 13 bit dynamic range, 1 µm pixelsize) are of the order of 20 min. After some basic image processing (subtraction of dark current, flat field correction), the 3D distribution of the linear attenuation coefficient µ is calculated using the filtered backprojection algorithm [START_REF] Kak | Principles of computerized tomographic imaging[END_REF]. The resulting threedimensional dataset consists in a set of isotropic voxels 3 and can be visualised and analysed using dedicated 3D image processing software packages.

It should be noted that the coherence properties of third generation synchrotron beams can be used to perform phase sensitive X-ray imaging. By increasing the sample to detector distance to a few centimeters, sample defects and heterogeneities give rise to pronounced interference effects (Fresnel diffraction, in-line holography). The use of this phenomenon leads to a contrast enhancement which improves the sensitivity of the tomography technique for detecting cracks in materials by at least one order of magnitude [START_REF] Cloetens | Observation of microstructure and damage in materials by phase sensitive radiography[END_REF]. Part of the work presented in this article is based on this phase sensitive imaging technique.

Grain boundary imaging

Currently only a few experimental techniques allow us to characterize the three-dimensional arrangement of grains in the bulk of a polycrystalline material: serial sectioning, either mechanically or by the focussed ion beam technique (limited to grain sizes below 20 µm [START_REF] Inkson | 3D determination of grain shape in a FeAl-based nanocomposite by 3D FIB tomography[END_REF]), and X-ray tracking [START_REF] Poulsen | Three-dimensional maps of grain boundaries and the stress-state of individual grains[END_REF][START_REF] Lauridsen | Tracking: a method for structural characterization of grains in powders or polycrystals[END_REF]. The sectioning techniques have the disadvantage of being extremely laborious and destructive. X-ray tracking on the other hand is a nondestructive technique based on Laue diffraction and provides the shape and crystallographic orientation of the individual grains in the bulk of a polycrystalline material. However, the current spatial resolution of the latter technique does not exceed 20-40 µm and it does not allow to detect the existence of cracks.

Due to their small lateral width, natural grain boundaries can not be detected by the change in X-ray attenuation. However, for the case of Al alloys we recently pointed out the possibility to visualize the 3D grain structure by using liquid Ga as a contrast agent for X-ray absorption tomography [START_REF] Ludwig | Penetration of liquid gallium into the grain boundaries of aluminium: a synchrotron radiation microtomographic investigation[END_REF]. When applied to the surface of a polycrystalline Al sample, liquid Ga (T m ϭ 29.7 °C) penetrates rapidly (Х 10 µm/s) [START_REF] Ludwig | Development and Applications of synchrotron radiation microtomography[END_REF] along the grain boundaries where it tends to form liquid layers of up to one micrometer thickness. Due to the fact that Ga has a considerably higher X-ray attenuation coefficient than Al, these layers can be easily detected by means of high resolution absorption microtomography (see Fig. 5). In order to improve the visibility of the wetting layer, the sample is annealed for about one hour at 100 °C before the acquisition of the tomographic dataset. This annealing leads to volume diffusion of Ga into the neighboring grains and increases the width and the continuity of the wetting layers in the tomographic reconstruction. Note that only grain boundaries which fulfill the wetting condition ν gb Ͼ 2ν sl , where ν gb and ν sl denote the grain boundary and solid-liquid interface tensions respectively, are prone to penetration. Special grain boundaries and low angle grain boundaries are known to have lower surface tensions and might therefore not be detected by this technique. It should also be noted, that the penetration process leads to the complete loss of the mechanical strength of the material (grain decohesion by the liquid metal) and this technique has, therefore, to be considered as a destructive one which can only be applied at the end of the fatigue experiment.

In order to be able to characterize and visualize the individual grains in a convenient way, the tomographic dataset has to be segmented and labeled in order to extract the grain structure. The transformation of the tomographic raw data, containing the three-dimensional arrangement of the grain boundaries, into its segmented representation is shown in Fig. 5. All voxels belonging to the same crystallographic grain (a grain is defined by the closed grain boundary surface which delimits it) have been attributed the same colour label. Despite its apparent simplicity for the human eye, the computer implementation of this procedure is not straightforward and requires a considerable amount of image processing. In a first step the data set is binarised by applying an intensity threshold. In order to keep all voxels belonging to the grain boundaries, this threshold has to be chosen rather low and results in a relatively noisy representation of the grain boundary structure. This binarised dataset in turn has to be cleaned by a series of morphological filter operations before it can be passed to the actual segmentation routine (3D watershed algorithm) [START_REF] Gratin | Mathematical morphology in three dimensions[END_REF][START_REF] Meyer | Mathematical morphology: from two dimensions to three dimensions[END_REF]. The segmented volume is finally labeled and can now be visualized with adapted volume rendering software packages.

Results

Observations of the crack at the surface of the sample: the validation of the method

The results of the crack propagation experiment are summarized in Fig. 2. The optical micrograph on the left side (Fig. 2a) shows the crack path at the sample surface (plane XZ on Fig. 1). As already observed for other samples of this material, crack initiation occurred on a Si particle/matrix interface at the immediate vicinity of a pore. Fig. 2b) shows an EBSD map of the same zone and reveals the position and orientation of the grains along the crack path. The crack, as well as the silicon rich particles cannot be indexed in the map and appear as black points. Note, that the Si particles are not only located at the grain boundaries, but also inside grains. The crack growth chronology, as determined from surface observations, is reported in Fig. 2e. From this figure, it can be seen that the crack stopped at the surface for Х 50000 and Х 25000 cycles on two grain boundaries visible on the EBSD surface map (points 'X' and 'Z'). A shorter arrest of the crack (Х 20000 cycles) was also observed (point 'Y') but does not seem to be correlated with the presence of a grain boundary at the surface. The tomographic reconstructions of the sample surface before and after application of gallium are shown on Fig. 2c andd respectively. Although the presence of the sample surface deteriorates the quality of those images slightly 4it appears that the crack shape observed on reconstructed images is very close to that observed by optical microscopy. The white contrast in Fig. 2d is due to the enhanced absorption caused by the presence of liquid Ga which has penetrated the grain boundaries [START_REF] Ludwig | Penetration of liquid gallium into the grain boundaries of aluminium: a synchrotron radiation microtomographic investigation[END_REF]. Note that the gallium also wets the crack surface. 5 The rather large width of the Ga rich zone (about 10 µm) is a result of the volume diffusion from the wetted grain boundary into the adjacent grains during the mentioned heat treatment at 100°C. Again, a good match is obtained between the reconstructed gallium layers and the grain boundaries as observed by EBSD (Fig. 2b). From those comparisons it is concluded that both the crack and the grain shapes in the interior of the sample can be studied with a good accuracy from the tomographic images.

Observations of the crack in the bulk of the sample

The crack path in the interior of the sample is strikingly more complex than the relatively linear propagation observed at the surface. This is illustrated on Fig. 3a andc which show two tomographic slices parallel to the XY plane defined in Fig. 1. Three marked changes in the direction of the crack can be seen on Fig. 3a. In addition, when going further in the material along the same direction, small crack segments appear, apparently disconnected from the main crack front (Fig. 3d).

For the post-mortem EBSD analysis of the sample (# 7 in the experimental procedure) the sample was cut along the plane depicted in Fig. 3a in order to analyse the grain orientations in this particular section of the sample. This study revealed that the planar portions of the crack are close (within 10°) to crystallographic {111} planes.

A better understanding of the crack morphology can be obtained when looking directly at a three dimensional representation of the crack surface. However because of the dual phase microstructure of the material (aluminium matrix plus silicon particles) this 3D representation is not straightforward. Indeed, in the reconstructed images, the gray levels of voxels corresponding to the crack and to the Si particles are both lower than the gray levels of voxels corresponding to the Al matrix; a mere segmentation based on gray levels would therefore extract the crack and the particles simultaneously. To avoid this problem, a sequence of 3D morphological filter operations [START_REF] Gratin | Mathematical morphology in three dimensions[END_REF] has been used and allowed the selective extraction of the crack surface from the 3D tomographic data. The result of this operation is illustrated on Fig. 3b andd for the two slices described previously. It can be seen that the binary image, although slightly less accurate that the initial one, still provides a good description of the crack. It is thus possible to obtain a 3D rendering of the binary data set representing the crack. Such a rendering is shown on Fig. 4 where the crack surface intersection, examined in EBSD and optical microscopy is indicated by an arrow. A close examination of this figure reveals that, in fact, the apparently unconnected crack segments described in Fig. 3c do belong to the crack front. 6This illustrates clearly the difficulty of studying crack propagation from simple 2D observations. In the central part of the crack front (label #2), one can also distinguish the bifurcation number 3 shown previously in Fig. 3a. One striking feature of this 3D view is that the crack front does not appear as a regular circular line. It is rather formed by three independent main segments (labelled 1,2 and 3) around 200 µm in length separated by uncracked Al matrix invisible in this rendering. This suggests that the crack front propagates first in "easy" zones of the crystal which might be grains showing favourable crystallographic orientations. In order to check this assumption we can now use the information about the 3D shape of the grains along the crack path.

Correlation between the crack shape and the internal grain distribution along the crack front

A first approach for studying the relation between the grain distribution and the crack morphology is to look at 2D reconstructed slices obtained after the wetting of the metal by gallium. Fig. 5a shows the identical section like Fig. 3a after application of the liquid metal. It appears that two of the abrupt changes in the crack propagation (1& 3) are related to the passage of the crack into a new grain. Close to label 1 the crack deviates by Ϸ 45°when crossing the grain boundary between grains A and B. Then, for unknown reasons, the crack deviates by Ϸ 90°in the center of grain B. Finally, when crossing the grain boundary between grains B and C the crack deviates again by Ϸ75°.

Thanks to the 3D nature of our observations, it is possible to study how the crack goes from grain B into grain C all along the crack front by simply looking at the different slices parallel to the XY planes. However, as discussed previously, this is better achieved by considering directly 3D images of the grains. 3D grain images have been obtained from the initial 3D data by morphological segmentation of the data set (see section 2.2). Fig. 5 illustrates the transformation of the tomographic reconstruction into the segmented and labelled representation. It can be seen that the grain boundary structure after segmentation (Fig. 5b) is in good agreement with the tomographic raw data: all grains have been detected and only a few inaccuracies can be observed (e.g. the additional small grain in the lower part, arising from over-segmentation). As already mentioned, gallium also wets the crack surface. Thus, when the crack is transgranular, the corresponding grain is split in two distinct cells by the segmentation algorithm. In that case, parts belonging to the same grain were attributed identical colours (in the 3D images), but still exist as distinct cells (see for example cells A 1 and A 2 in Fig. 5b) in the numerical representation of the data set. Both data sets, the one containing only the crack and the labelled one containing the grains, were carefully aligned with respect to each other and finally merged for a simultaneous 3D visualization of both the crack and the grains. An example of such a visualization is shown in Fig. 6 which represents the 3D shape of two of the grains observed by EBSD in Fig. 2b. For a direct comparison the grains have been assigned the same colors in the EBSD map and in the 3D image. Part of grain G2 has been set transparent in the visualization: the smooth surface on the right side corresponds to a tomographic section close to the sample surface observed by optical microscopy and the rather rough surface (nearly parallel to the YZ plane) corresponds to one crack face in the grain.

It is interesting to note that the point Y7 corresponds well to the position, where the sub-surface crack front (illustrated schematically by the black lines in Fig. 6a had to cross a section of the grain boundary lying parallel to the crack front. A possible interpretation of this observation is given at the end of section 4.1. The bifurcation number 3 on Fig. 5a can also be studied with the help of the 3D rendering of the grains. The results of the analysis are shown in Figs. 7 and8. The first figure shows the 3D shape of the grain which is responsible for the crack bifurcation. It appears (Fig. 8a andb) that the crack enters a new grain (the one depicted in Fig. 7) along the curved lined A A". This line can be further divided in two sections: a first section AA' corresponding to a large advance (Х 100 µm) of the crack in the grain by tilting around an axis which is roughly the intersection of the crack front with the boundary of the grain; and a second section A'A" corresponding to a much smaller advance of the crack (Х 10-20 µm) in the new grain by twisting of the crack front. A possible explanation of the preferential advance of the crack along A-A' might be given in terms of a continuity condition for the crack surface, as discussed in section 4.2.

Discussion

To the best of our knowledge, this is the first time that a simultaneous 3D image of a short fatigue crack and of the surrounding crystallographic grains has been obtained. The combination of the 3D tomographic measurements with conventional 2D surface observations (crack growth history and EBSD mapping) provides interesting information about various aspects of the interaction of short fatigue cracks with the grain structure in polycrystalline material.

Crack arrest at grain boundaries

One may distinguish two types of models describing the influence of the grain structure on fatigue crack propagation in polycrystals:

1. Models based on the varying elastic properties between neighboring grains. 2. Models based on the dislocation mechanisms and crystal plasticity ahead of the crack tip.

Representative for the first class of models we consider here recent work by Ravichandran and Li It is likely that the preferential attack along A-A' is governed by compatibility with favorable crystallographic orientations inside the grain.

[33] which accounts for the influence of the varying elastic modulus in the grains due to the change in grain orientation. Assuming isostrain conditions, these authors have established a method to calculate the local values of the stress intensity factor K along a crack front in heterogeneously stressed material and their model predicts strong variations of the local K values in the case of microstructural short cracks.

In order to test the applicability of this model, we calculated the elastic modulus of the surface grains along the loading direction, using the EBSD mapping and the relationship [START_REF] Nye | Physical properties of crystals[END_REF]

E ϭ ͩ s 11 Ϫ2 ͩ s 11 Ϫs 12 Ϫ ͩ s 44 2 ͪͪ (l 2 1 l 2 2 ϩ l 2 2 l 2 3 ( 1
)
ϩ l 2 3 l 2 1 ) ͪ Ϫ1
where the terms s ij correspond to the elastic constants of Al and l i to the coordinates of the loading direction, expressed in terms of the vectors of the unit cell. The results are summarized in Table 2.

Due to the weak coefficient of anisotropy 8 of Al (1.2 compared to 2.5 in the case of Fe [START_REF] Cottrell | The mechanical properties of matter[END_REF]) the elastic moduli show only moderate variations (7% with respect to the mean value) and the above mentioned mechanism is not expected to play a dominant role. Indeed, two of the crack stops observed at the sample surface (points X and Z) occur at the transition from a grain of low modulus into a grain of high elastic modulus, whereas one would expect the opposite behavior (acceleration rather than retardation) of the crack propagation rate in this situation.

This, together with the experimental evidence of localized plasticity ahead of the crack tip [START_REF] Savelli | Analyse des me ´canismes et mode ´lisation de la dure ´e de vie en fatigue d'alliages de moulage AS7G03[END_REF] seems to indicate that the second class of mechanisms plays a preponderant role in the present case. The common point of these mechanisms is the postulate that some plasticity has to be initiated in the new grain before the crack can enter it [START_REF] Zhang | Measurement of plastic zones associated with small fatigue cracks by selected area channeling patterns[END_REF][START_REF] Turnbull | The effect of grain size on the fatigue of commercially pure aluminium[END_REF][START_REF] Lankford | The growth of small fatigue cracks in 7075-T6 aluminium[END_REF]. As a result the "waiting time" in the grain boundary should be a function of the misorientation between the two grains. Such stops are likely to happen frequently along the crack front, when the crack is small.

The quantitative analysis of the crack propagation data (i.e. the observed crack arrests at grain boundaries) in terms of such a dislocation based model, taking into account the local differences in the activation of the crystallographic slip systems, is under progress and will be reported elsewhere.

Irrespective of the underlying mechanism, the restriction of the observation of crack propagation on the sample surface may lead to erroneous Table 2 Elastic moduli along the loading direction, determined from the EBSD mapping (Fig. 2b interpretation of the growth kinetics as the threedimensional nature of the crack is neglected. There is experimental evidence that significant changes in the aspect ratio 2c/a (c surface length of the crack, a bulk length) may occur during the growth of short cracks [START_REF] Ravichandran | Three-dimensional crack-shape effects during the growth of small surface cracks under fatigue in a titanium-base alloy[END_REF]. The ideal characterization procedure should therefore include the continuous monitoring of the 3D crack shape during the fatigue experiment and the determination of the grain orientations in the bulk of the polycrystalline sample. The former task can be achieved by repeated microtomographic inspection during an interrupted fatigue test and the latter task can be performed with the help of a new non-destructive, synchrotron based diffraction technique [START_REF] Poulsen | Three-dimensional maps of grain boundaries and the stress-state of individual grains[END_REF][START_REF] Lauridsen | Tracking: a method for structural characterization of grains in powders or polycrystals[END_REF].

For the current study we do not dispose over such detailed information on the crack growth history in the bulk of the sample. Nevertheless it is interesting to note that the crack retardation observed at point Y on the sample surface corresponds well to the position where the sub-surface segment of the crack has to pass a section of the grain boundary which lies parallel to the crack front. A possible interpretation of the retardation at point Y can therefore be given in terms of a crack arrest below the sample surface, which in turn influences the crack growth at the surface. The dashed lines in Fig. 6 have been added to illustrate the possible evolution of the crack front according to the above mentioned scenario.

Tilt and twist

The 3D tomographic reconstruction gives direct access to geometrical and morphological parameters of the crack surface in the bulk of the sample. This allows us, for example, to determine in a quantitative way the mixed-mode character of short cracks which often consist of several segments growing in different directions with respect to the loading axis. Coupled with in situ monitoring of the crack growth during the fatigue experiment, this information could be used in the future to assess the growth kinetics in a full fracture mechanical context, for example, by including the effect of crack deflections on the stress intensity factor [START_REF] Suresh | Crack deflection: Implications for the growth of long and short fatigue cracks[END_REF].

However, interesting information can also be obtained from the analysis of the 3D crack morphology at a given state during fatigue life, i.e. without the knowledge of the crack history. The investigated crack consists of three segments with preferred growth directions. The analysis of the strongly inclined central segment with respect to the surrounding grain structure clearly revealed that the deflection of the crack is related to the passage of the crack into a new grain (Figs. 7 and8).

Two factors are to be taken into account for the propagation of the crack. The first one is the activation of some plasticity on a crystallographic well orientated plane in the next grain. Our observations show that in the studied case the crack plane is close to a {111} plane. The second factor that has to be considered is the geometry of the crack surface itself. It has been shown recently [START_REF] Zhai | A crystallographic mechanism for fatigue crack propagation through grain boundaries[END_REF] by EBSD surface observations of an Al-Li alloy that short fatigue cracks growing on crystallographic planes tend to propagate by tilting rather than twisting of the crack front.

The reason for this is illustrated schematically in Fig. 9. The propagation of the crack by simple tilt of the crack front is more favorable, in terms of energy, than a propagation by twist. In the second case, indeed, the fracture surface is larger.

Our observations of the crack penetrating the grain depicted in Fig. 7 are in favor of this mechanism: the propagation of the crack in this grain is larger in the region where it can be accommodated by simple tilt (line AA' on Fig. 8). On the contrary, the grain resists propagation when it is accompanied by a twist of the crack front (between A' and A").

The observation of significant changes in the growth directions of cracks in the same material by conventional post mortem fractographic analysis (Fig. 10) gives additional support for the above mentioned mechanism: rather than homogeneous propagation into all grains close to the crack tip, the crack seems to follow the perimeter of some of the grains until propagation into these grain starts in a lateral direction.

It should be emphasized that the analysis of a larger number of fatigue cracks, including the crystallographic characterization of the adjacent grains, is necessary to corroborate the validity of this tilting/twisting mechanism.

Conclusions

This paper presents a new technique for the simultaneous three-dimensional visualization of fatigue cracks together with the surrounding grain structure in the bulk of polycrystalline Al samples. Penetration of liquid Ga along the Al grain boundaries leads to the formation of microscopic wetting layers which in turn can be detected by high resolution X-ray absorption tomography.

The studied short fatigue crack appears to consist of three crack segments corresponding to the crack growth in different grains. The stop and go motion of the crack, observed by optical microscopy at the sample surface, is clearly related to the grain structure directly at and close to the sample surface. The tentative interpretation of these crack stops in terms of an elastic model supposing isostrain conditions fails, presumably due to the weak coefficient of anisotropy of Al. The three-dimensional analysis of the crack with respect to the grain structure indicates that the passage of the crack into a new grain occurs preferentially from regions on the grain perimeter, where the growth can be accommodated by tilting of the crack plane. On the contrary, the crack has difficulties in entering the new grain from regions where twisting of the crack plane is required and tends to grow along the grain perimeter.

Future experiments should include repeated microtomographic characterization of the 3D crack shape during (interrupted) fatigue tests. This, together with the characterization of the grain orientation in the bulk of the material via modern synchrotron diffraction techniques will provide the necessary information to assess the problem of fatigue crack propagation in terms of a full threedimensional fracture mechanical context and to analyse the physical mechanisms governing crack propagation.

Fig. 1 .

 1 Fig. 1. Schematic representation of the geometry of the fatigue sample. The indicated small parallelepiped (0.7×0.7×5 mm) containing the short fatigue crack was wire-cut from the original specimen and characterized by synchrotron radiation microtomography.

Fig. 2 .

 2 Fig. 2. a) Optical micrograph of the fatigue crack at the sample surface b) EBSD mapping of the same zone c) tomographic reconstruction of the sample surface before application of Gallium (energy: 15 keV, distance: 5 mm) d) tomographic reconstruction after application of Gallium e) plot of crack length versus number of cycles as observed by optical microscopy during the interrupted fatigue test.

Fig. 3

 3 Fig. 3. a,c): Tomographic slices perpendicular to the surface of observation revealing the complicated morphology of the crack in the bulk of the sample (energy: 15 keV, distance 42 mm). b,d): corresponding slices after segmentation by 3D morphologic reconstruction.

Fig. 4 .

 4 Fig. 4. Three dimensional rendering of the crack surface. Note the abrupt change in inclination of the central arc-shaped crack segment (#2). The finger-like crack extensions (see arrows) correspond to the crack fragments depicted in Fig. 3d.

Fig. 5 .

 5 Fig. 5. a) Tomographic reconstruction after application of Ga. The depicted slice corresponds to the same section shown in Fig. 3a (energy: 15 keV, distance: 5 mm). b) Result of 3D morphologic segmentation: each cell detected by the segmentation algorithm receives a distinct label, represented here by different colours. Cells belonging to the same grain (A 1 &A 2 and B 1 &B 2 ) have been attributed identical colours. c) The superposition of the tomographic and the segmented data set shows a good agreement between both representations.

Fig. 6 .

 6 Fig. 6. a) 3D rendering of the local grain structure close to the crack stops observed at the sample surface at points X and Y. Both stops might be attributed to the pinning of the crack front at the grain boundary between grains G1 and G2. b) EBSD mapping of the corresponding sample area.

Fig. 7 .

 7 Fig. 7. The abrupt deviation in the crack propagation direction coincides with the transition into a new grain. a) Volume rendering of the crack surface as in Fig. 3 but from a different viewpoint. b) Rendering of the crack surface and the grain that caused the sudden deviation.

Fig. 8 .

 8 Fig.8. a) 3D rendering of the crack inside the new grain (set to transparent). b) same scene with the surrounding grains. The crack entered the new grain between points A and A' by tilting. Between points A' and A" the crack followed the contour of the grain. It is likely that the preferential attack along A-A' is governed by compatibility with favorable crystallographic orientations inside the grain.

8

  Defined as the ratio of the elastic constants 2c 44 /c 11 -c 12 .

Fig. 9 .

 9 Fig. 9. The adaptation of a crack to a new plane may involve tilting and twisting. Note that twisting requires the creation of additional surfaces (black) in order to preserve the connectivity of the crack surface.

Fig. 10 .

 10 Fig. 10. SEM micrograph of the fracture surface of a sample tested at 120 MPa (N R =2075200). The characteristic crack growth patterns indicate significant changes in the propagation direction of the crack.

Table 1

 1

	Chemical composition (mass %) of the studied model cast Al
	alloy.
	Al Si	Fe	Mg	Ti	Sb	Cr rest
	92.5 6.6-7 0.09-	0.29-	0.11-	0.12-	0.2 Յ0.01
		0.14	0.34	0.14	0.14

Current address: Lab GEMPPM, INSA-Lyon, 69621 Villeurbanne Cedex, France.

For tomographic reconstruction the sample has to fit laterally entirely into the field of view.

3D equivalent to the 2D pixel; represents the smallest digital volume of the material from which the X-ray attenuation is determined.

On reconstructed slices in the bluk of the material, individual Si particles can be more clearly visualized, see for example Fig. 3a.

[START_REF] Lankford | The growth of small fatigue cracks in 7075-T6 aluminium[END_REF] It is always possible, however, to distinguish the crack from a grain boundary by comparing with the tomographic scan obtained before the application of gallium.

Details labeled "finger" in Fig.4.

Retardation of surface crack propagation for about 20,000 cycles.

Acknowledgements

We would like to thank P.H Jouneau (INSA Lyon) and P. Buffat (CIME, Lausanne) for the EBSD characterization. We are also much indebted to D. Jeulin and S. Bouchet (Ecole des Mines, Paris) who provided and implemented the 3D image segmentation routines. D. Bellet and Y. Brechet (INPG, Grenoble) are acknowledged for useful discussions on various aspects of the work.


Ivan Duchemin

Antoine Levitt

Computing photoionization spectra in Gaussian basis sets

We present a method to compute the photoionization spectra of atoms and molecules in linear response time-dependent density functional theory. The electronic orbital variations corresponding to ionized electrons are expanded on a basis set of delocalized functions obtained as the solution of the inhomogeneous Helmholtz equation with gaussian basis set functions as right-hand side. The resulting scheme is able to reproduce photoionization spectra without any need for artificial regularization or localization. We demonstrate that it is able to produce accurate spectra for semilocal exchange-correlation functionals even using relatively small standard gaussian basis sets.

Introduction

Time-dependent functional theory in the linear response regime (LR-TDDFT) is widely used to compute excitation energies of molecules. The traditional way to compute these properties involves solving the Casida equations [START_REF] Mark | Time-dependent density functional response theory for molecules[END_REF], which are the discretized LR-TDDFT equations in a basis set of localized functions (typically, gaussian basis sets). This is very convenient for the excitations lying below the ionization potential, for which the orbital variations are localized functions. This however fails to reproduce the photoionization spectrum, turning a continuous function into a series of infinitely sharp peaks (Dirac deltas). These peaks approximate the function in a weak sense (when integrated against smooth functions) in the limit of complete basis sets [START_REF] Dupuy | Finite-size effects in response functions of molecular systems[END_REF], but pointwise values are not easily obtained.

Various techniques have been used to remedy this. The simplest and most general technique is to use an artificial dissipation parameter η. This effectively adds a constant imaginary part to the Hamiltonian, broadening the infinitely sharp peaks and producing a continuous function. Another related technique is the use of complex absorbing potentials, where the imaginary part is only made to act away from the molecule, effectively attempting to broaden only the states that correspond to ionized electrons [START_REF] Jg Muga | Complex absorbing potentials[END_REF]. Another variation on the same idea is to add the imaginary part to the molecular orbitals energies directly [START_REF] Coccia | Ab initio lifetime correction to scattering states for time-dependent electronic-structure calculations with incomplete basis sets[END_REF]. When implemented in the time domain, these techniques can be understood as adding an artificial dissipation, resulting in states with finite lifetimes (complex energies).

These schemes shift the poles of the response function away from the real axis into the lower complex plane, resulting in response functions that can mathematically be expressed as sums of Lorentzians. Howerver, this does not respect the mathematical structure of photoionization spectra (infinitely sharp peaks before the ionization potential, continuous function with sharp changes around the ionization potentials as well as possible resonances). As a result, when the basis sets are small (as is the case for typical gaussian basis sets), spectra are often unsatisfactory.

Another class of methods attempts to recover the continuous photoionization spectrum from the discrete peaks by fitting rational functions to quantities that are computable on localized basis sets, such as moments of the oscillator strengths or values of the polarizability in the complex plane: examples include Stieljes imaging or Padé extrapolation (see [START_REF] Nunes | Molecular photoionization and photodetachment cross sections based on l 2 basis sets: Theory and selected examples[END_REF] for a review). These methods are intrinsically numerically unstable and need a non-trivial manual parameter selection.

More sophisticated schemes like the complex scaling and exterior complex scaling methods use the analytic continuation of the solutions to transform oscillatory tails into decaying ones, resulting in more accurate 1 spectra [START_REF] Cerioni | Accurate complex scaling of three dimensional numerical potentials[END_REF]. These schemes are however non-trivial to implement, and require significant fine-tuning of parameters (as do the dissipation-based schemes).

A more principled class of methods involve solving an approximate form of the equation exactly outside of the computational domain, resulting in a Dirichlet-to-Neumann map that acts as an effective boundary condition [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF]. These schemes are however hard to implement in systems without spherical symmetry.

In this paper, inspired by a method we recently developed to compute resonances in locally perturbed periodic media [START_REF] Duchemin | Efficient extraction of resonant states in systems with defects[END_REF], we propose a new scheme that does not rely on artificial dissipation or localization methods, and works for arbitrary molecules without symmetries. The method can be summarized as follows. The LR-TDDFT equations for the orbital variations δψ ± i are of the Helmholtz form

(±ω + i0 + + ε i + 1 2 ∆)δψ ± i = f
where f is a localized function (that depends self-consistently on all the δψ ± i ) and ε i < 0 are the Kohn-Sham one-particle energies. The i0 + in this equation means that the equation is to be solved for finite positive values iη, and η is to be taken to tend to zero after the calculation.

From the form of the Green function of the operator ±ω + i0 + + ε i + 1 2 ∆, it can be seen that, in the limit η → 0 + , for ω > -ε i , δψ + i will be delocalized. Physically, δψ + i represents an electron ionization; the positive values of η correspond to imposing an outgoing wave boundary condition. It is not appropriate to discretize the delocalized δψ + i on a localized basis set; in fact, doing so results in a singular photoionization limit that, in the limit η → 0 + , tends to a finite sum of Dirac masses located at the Casida excitation energies. These then have to be regularized to obtain a continuous spectrum. Instead, we perform the change of variable

δψ + i = (ω + i0 + + ε i + 1 2 ∆) -1 δφ + i ,
and discretize φ + i (which is localized) on a standard gaussian basis set. Equivalently, we discretize δψ + i on a basis set consisting of solutions of the Helmholtz equation with gaussian basis set functions as right-hand side.

By respecting the structure of the equations (localized excitation or delocalized ionization, depending on the value of ±ω + ε i ), this method allows us to compute photoionization spectra directly. Compared to the standard method of solving the Casida equations, it is much more accurate, resulting in accurate spectra using very moderate basis sets. It can be easily adapted to compute resonances, and is fully general, being suited to molecules as well as atoms. The flip side is an added computational cost, which, although formally having the same cubic scaling as usual TDDFT methods, has a higher prefactor, mostly due to the need for frequency-dependent integrals on a real space grid. These are however highly parallelizable.

Methods

Model

We consider a molecule with N spin-paired electrons modeled using time-dependent adiabatic Kohn-Sham density functional theory with a semilocal functional, and an electrostatic potential V nucl originating from the nuclei. The ground-state orbitals ψ i satisfy, for all i ∈ {1, . . . , N/2},

H[ρ]ψ i = ε i ψ i
where, in atomic units,

H[ρ] = -1 2 ∆ + V tot [ρ] V tot = V nucl + V Hxc [ρ] V Hxc [ρ](r) = ρ(r ) |r -r | dr + v xc (ρ(r))
The total density is ρ(r) = 2

N/2 i=1 |ψ i (r)| 2 .
We choose the ground-state orbitals to be real for simplicity, but the scheme extends trivially to complex orbitals.

Consider now a perturbing potential δV P . In the time-harmonic regime, the first-order response can be described by the Sternheimer equations [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF] (±ω + iη

+ ε i -H[ρ])δψ ± i -(f HXC δρ)ψ i = δV P ψ i
The variation in the Hartree-exchange-correlation kernel is given by

(f HXC δρ)(r) = δρ(r ) |r -r | dr + v xc (ρ(r))δρ(r).
and the variation in density by

δρ(r) = 2 N/2 i=1 ψ i (r)(δψ + i (r) + δψ - i (r)).
This defines the total variation in polarization δP = rδρ(r)dr.

When δV P (r) = -e • r, the linear relationship δP = α(ω + iη)e defines the polarizability tensor α, a 3 × 3 matrix. Finally, the photoionization cross-section is given by

σ(ω) = lim η→0 + 4πω c 1 3 Tr(α(ω + iη)),
with c the speed of light, and will be our main observable of interest.

Integral form and delocalization

The standard approach to discretizing these equations is to expand δψ i in a basis, which leads to the usual Casida equations. This is however inefficient when ω is greater than -ε i , at which point the electron becomes ionized and δψ + i is delocalized. To see this, we write the Sternheimer equations as

(±ω + iη + ε i + 1 2 ∆)δψ ± i = V tot δψ ± i + (f HXC δρ)ψ i + δV P ψ i (1)
Since V tot and ψ i are localized, the right-hand side is localized. Introduce the free Green's function G 0 (r, r ; z), the kernel of the inverse of the operator z + 1 2 ∆, well-defined when η > 0. We can then reformulate the Sternheimer equations in integral form

δψ ± i = G 0 (±ω + iη + ε i ) V tot δψ ± i + (f HXC δρ)ψ i + δV P ψ i
An explicit computation shows that the kernel of the Green function is given by

G 0 (r, r , z) = - 1 2π e ik(z)|r-r | |r -r |
where k(z) is the square root of 2z with positive imaginary part (so that G 0 is localized when η > 0). When z = ±ω + iη + ε i approaches the real axis, this Green function is oscillatory when Re(z) > 0, and decaying when Re(z) < 0. In the region ω > 0, -ω + ε i is always negative, and therefore δψ - i will be localized. However, the behavior of ψ + i depends on whether ω < -ε i or ω > ε i (below or above ionization threshold). In the case ω < -ε i , δψ + i will be localized; in the case ω < -ε i , it will be oscillatory. This explains why usual methods, based on the discretization of δψ i , have trouble reproducing the ionization region ω > -ε i .

Our method

Our approach is to use the change of variables

δψ + i = G 0 (±ω + iη + ε i )δφ + i when ω + ε i > 0
and to discretize δφ + i in a Gaussian basis set instead of δψ + i directly. This ensures automatically the correct asymptotic behavior for δψ + i . The δφ + i satisfy the integral equation

δφ + i -V tot G 0 (+ω + iη + ε i )δφ + i -(f HXC δρ)ψ i = δV P ψ i (2)
When ±ω + ε i ≤ 0, we discretize δψ ± i in the usual way on a basis of Gaussian-type orbitals (χ α ) α=1,...,N b . When ±ω + ε i > 0, we discretize δφ ± i on the basis:

δφ + i (r) = N b α=1 a + iα χ α (r) when ω + ε i > 0 δψ ± i (r) = N b α=1 b ± iα χ α (r) otherwise
This sets up a linear system in the coefficients (a, b). Note that from (2) it follows that δφ + i is as localized as V tot ψ + i . Therefore, the decay of the total meanfield potential determines the localization of δφ + i , and therefore the effectiveness of the numerical method. When using semilocal density functionals (with exponentially decaying exchange-correlation potentials), this is determined by the electrostatic potential. We can therefore rank systems by decreasing order of localization: atoms (exponentially decaying potential), nonpolar molecules (potential decaying as 1/r 3 ), polar molecules (1/r 2 ), charged systems (1/r). When using hybrid functionals incorporing Hartree-Fock exchange, the effective potential seen by ionized electrons behaves as 1/r [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF], and we expect our method to have difficulties.

Solution of the linear system

We project the Sternheimer equations ( 1) and (2) on the basis to obtain

(S -K + i )a + i -g i [a, b] = h i when ω + ε i > 0 ((±ω + ε i + iη)S -H)b ± i -g i [a, b] = h i otherwise, (3)
where

S αβ = χ α |χ β H αβ = χ α | -1 2 ∆ + V tot |χ β K + iαβ = χ α |V tot G 0 (+ω + iη + ε i )|χ β g iα [a, b] = χ α |f HXC δρ[a, b]|ψ i h iα = χ α |δV P |ψ i
The computation of S, H and h i are standard. For K, the matrix elements are not analytic. However, the values (G 0 (±ω + iη + ε i )χ β )(r) can be computed analytically for all r (see Appendix). Therefore, we introduce an integration grid with points r l and weights w l , and approximate K αβ as

K αβ ≈ l w l χ α (r l )V tot (r l )(G 0 (±ω + iη + ε i )χ β )(r l ).
The values of V tot (r l ) arising from the nuclei and exchange-correlation terms are computed exactly, as is usual in DFT. Note that since the delocalizing operator G 0 is then multiplied by localized quantities, the grid only needs to be sufficient to integrated localized functions. In practice, we found that a coarse exchangecorrelation grid was often adequate.

The computation of g, assuming that all the (a, b) are known, is more conveniently reformulated as

g iα = δρ[a, b]|f HXC ψ i χ α .
The values of f HXC ψ i χ α on the grid are precomputed, using the same technique as for the computation of V tot . Then, δρ[a, b](r) = 2

N/2 j=1 ψ j (r)(δψ + j (r) + δψ - j (r)) is formed on the grid, by using the values of G 0 (±ω + iη + ε i )χ α on the grid and the (a, b) coefficients.

The linear system (3) is solved with the GMRES iterative solver, preconditioned by S -K + i (for the a block) and (±ω + ε i )S -H (for the b block).

The additional computational cost compared to standard iterative TDDFT computations is summarized in Table 1. Note that these steps are cubic scaling (and, for low-lying excitations where the number of ionized electrons is of the order of unity, quadratically scaling) except for the steps involving f HXC . The scaling can be improved using techniques such as the resolution of the identity. However, in our tests, the step involving the computation of the values of G 0 χ α on the grid dominated the overall computational time, and therefore we did not optimize the other steps.

Operation Cost Precomputations of f HXC ψ i χ α N g N occ N b + N 2 b N g N occ Computation of G 0 (±ω + iη + ε i )χ α on the grid N ω N ionized N g N b Matrix-vector products with K N ω N ionized N iter N g N b
Table 1: Dominant scaling of the main operations compared to standard iterative TDDFT. N b is the number of basis functions, N g the number of grid points, N occ the number of occupied orbitals, N ionized the number of ionized orbitals. Typically, N ionized N occ N b N g . Additionally, N ω is the number of frequencies desired, and N iter is the number of iterations of the iterative solver (typically, ≤ 10).

Results

We implemented the method in the Julia programming language [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF], interfacing with the PySCF package [START_REF] Sun | Recent developments in the pyscf program package[END_REF] to perform the initial setup and DFT run. The integrals described in the Appendix were implemented in the GaIn Fortran library. The code is freely available at https://github.com/antoine-levitt/ PhotoionizationGTO.jl and https://gitlab.maisondelasimulation.fr/beDeft/GaIn.

All results presented below use the LDA exchange-correlation functional, with no spin polarization. The use of the LDA functional is inadequate to obtain accurate photoionization spectra for the systems studied here, but the emphasis in this paper is on the methodology rather than on the particular results.

Unless explicitly mentioned, we used a very coarse exchange-correlation grid (PySCF setting 1), which yielded reasonably accurate results at minimal computational cost. We used the Dunning augmented basis sets [START_REF] Thom | Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen[END_REF][START_REF] Kendall | Electron affinities of the first-row atoms revisited. systematic basis sets and wave functions[END_REF], as provided by the basis set exchange [START_REF] Benjamin P Pritchard | New basis set exchange: An open, up-to-date resource for the molecular sciences community[END_REF]. These basis sets are designed for converging post-Hartree-Fock methods rather than TDDFT properties, and are very suboptimal here. We simply use them to demonstrate that acceptable convergence can be obtained in basis sets not specifically designed for that purpose.

Atoms

Atoms have an exponentially decaying total potential; therefore, we expect the δφ to be exponentially localized, making it an ideal case for our method. On atoms, we are able to compare the results to a reference (black line) computed using the atom-specific method of [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF]. We focus in this section on He and Be atoms.

Helium is the simplest system, with only one occupied orbital (1s). Accordingly, its TDLDA photoionization spectrum has only one threshold, plotted as a dashed vertical line. We see that results are already The result of the standard Casida method with damping (η = 0.5) in the d-aug-cc-pvqz basis set is shown for comparison. Reference data (black line) from the method of [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF].

qualitatively consistent even at extremely coarse basis sets: note for instance that the basis set aug-cc-pvdz only has two basis functions for the p channel relevant here. To check systematic convergence, we used a basis of even-tempered gaussians with N gaussian exponents logarithmically spaced from 0.01 to 10 (these values are taken somewhat arbitrarilly and not optimized).

To appreciate how inadequate the Casida method is to compute photoionization spectra on this system, note that, in the largest basis set used here, d-aug-cc-pvqz and in the frequency range displayed here, there are only two relevant (1s → 2p) excitations, at ω = 0.64 and ω = 0.91. This is clearly not sufficient to reconstruct a full spectrum. The same goes for other approaches such as complex absorbing potentials; these methods, even with optimal parameters, will only move these two poles in the complex plane, which is not sufficient to reproduce the full structure of the spectrum.

To illustrate, we plot Figure 2 for ω = 1 the orbital variation δψ + , as well as the localized δφ + defined by δψ + = G 0 (ω + i0 + )δφ + . It is clearly not possible to represent the delocalized variation δψ + by localized orbitals, but δφ + is relatively short-range and therefore can reasonably be expanded on a gaussian basis set.

Beryllium has two occupied orbitals (1s and 2s), yielding two different thresholds. Here, small standard basis sets are inaccurate, showing a displaced peak after the 2s ionization and unphysical oscillations after the 1s ionization. Using even-tempered basis set (using 10 gaussians with exponents logarithmically spaced between 0.01 and 10) yields an almost perfect ionization spectrum.

The ionization threshold of the 2s orbital lies below the 1s→2p excitation energy, which turns into a resonance. This resonance, extremely hard to capture with standard damping methods, is present and relatively accurate even with very coarse basis sets. Note that, since we are able to compute analytic continuations of the matrix elements of the free Green function as z crosses the positive real axis, we could compute this resonance directly [START_REF] Duchemin | Efficient extraction of resonant states in systems with defects[END_REF], but we do not pursue this direction in this paper.

Nonpolar molecules

We next try our method on nonpolar molecules, on which the asymptotic decay of the total potential is dictated by the quadrupole moment (decaying thus as 1/r 3 ). The test was conducted for H 2 and CH 4 , demonstrating a rapid convergence with respect to the basis set used: in both cases the spectrum is already almost converged at the aug-ccpvdz level. In particular, we find that moving from aug to d-aug basis sets turns out to be much more efficient, convergence-wise, than increasing the zeta label of the basis. This confirms again that the relevant components lies in the asymptotic part of the wave-function, supported by the delocalized atomic orbitals.

Polar molecule

We test the method on the strongly polarized LiH molecule. The strong dipole moment makes it hard for the method to converge, especially after the second ionization. Spurious oscillations after the peaks appear, which are consistent with what has been observed using approximate boundary conditions in radial methods [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF]. However, we remark that these difficulties do not hinder the capture of both ionization thresholds.

Charged systems

We use Li + as a benchmark in this case. As any other charged systems, ions have a long-range total potential decaying in 1/r. In such a case, the asymptotic form of the δψ is modified from a plane wave to a Coulomb wave [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF]Appendix]. Therefore, the functions δφ become much more delocalized, and our method is inadequate.

Accordingly, there are strong oscillations after the ionization threshold, which disappear very slowly when increasing the basis set size; this is consistent with Figure 5 of [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF]. For this example, the coarse integration grids we used before were not sufficient to resolve the fine oscillations, and we used a Gauss-Chebyshev grid with 100 points in the radial direction produced by PySCF.

Conclusion

We have presented a new method to compute photoionization spectra for TDDFT with semilocal exchangecorrelation functionals in gaussian basis sets, and tested it on atoms and small molecules. The method appears to be very efficient on atoms and nonpolar molecules; it struggles on polar molecules and charged systems.

We note the following possible improvements to our methodology

• For simplicity, we used standard Dunning basis sets, with exchange-correlation integration grids. Both of these were designed for a different problem (converging ground states orbitals and exchangecorrelation integrals) than the one addressed here, and it is likely that computational efficiency can be significantly improved by tailoring these to our problem rather than using off-the-shelf technology.

• The current methodology scales formally as the fourth power of the number of electrons, with the possibility of cubic scaling at the cost of a higher prefactor. It should be possible to reduce this scaling even further using techniques like the resolution of the identity.

• Although we tested our scheme on semilocal density functionals, it is in theory straightforwardly adaptable to hybrid functionals. However, since the Sternheimer equation is effectively long-range in these cases [START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF], we expect to face the same difficulties as we did in the ionic case. An efficient treatment of the Coulomb potential is, as far as we know, an open problem in the non-radial case. A partial solution could be to truncate artificially the potential; this would yield better basis set convergence and eliminate the oscillations, at the price of a regularization of the problem (particularly important near ionization thresholds).

• We only explored linear response TDDFT in the frequency domain; it would be interesting to generalize this methodology to the non-perturbative TDDFT equations in the time domain. Recent progress has been made in this direction in [START_REF] Kaye | Eliminating artificial boundary conditions in time-dependent density functional theory using fourier contour deformation[END_REF], using a similar decoupling between the kinetic and potential operators.

Finally, we note that we have focused here on atoms and small molecules. Experience suggests that it might be harder to converge photoionization spectra of small molecules than of large systems, because the relatively fine features computed here average out over a molecule, and because the basis sets used for larger molecules cover a larger region of space. Exploring this further is an interesting topic for future research.

Acknowledgments

We are grateful to Karno Schwinn and Julien Toulouse for the reference data on atoms and ions.

Appendix: matrix elements of the Helmholtz kernel in Gaussian basis sets

Our goal is the computation of integrals of the form I = g 1 , G 0 (ω + i0 + )g 2 = -1 2π R 6 g 1 (r)g 2 (r ) e -λ|r-r | |r -r | .

when g 1 and g 2 are gaussian-type orbitals, for complex values of λ. The case of interest in the present paper is λ = -i √ 2ω (as well as possible analytic continuations of this). This also includes pointwise values of G 0 (ω + i0 + )g 2 by taking for g 1 the limit of a zero-width gaussian. Explicit formulas have been obtained in the context of range-separated hybrids with Yukawa potentials [START_REF] Ten-No | Initiation of explicitly correlated slater-type geminal theory[END_REF][START_REF] Ten-No | New implementation of second-order møller-plesset perturbation theory with an analytic slater-type geminal[END_REF][START_REF] Akinaga | Range-separation by the yukawa potential in long-range corrected density functional theory with gaussian-type basis functions[END_REF]; we recall here the method of computation.

As is standard, it suffices to compute the integral for the primitive radial gaussians g 1 (r) = e -α|r-R1| 2 , g 2 (r) = e -β|r-R2| 2 . Expressions for integrals involving higher angular momenta can be obtained by differentiating the integrals with respect to the centers R 1 and R 2 .

We use the integral representation:

e -λx x = 2 π ∞ 0 dt e -1 2 (x 2 t 2 + λ 2 t 2 )
and relatively straightforward but tedious changes of variables to obtain The computation of derivatives requires special care for numerical stability. We refer to [START_REF] Ten-No | New implementation of second-order møller-plesset perturbation theory with an analytic slater-type geminal[END_REF] for a discussion of the issues and robust methods. We have implemented a similar method to that of [START_REF] Ten-No | New implementation of second-order møller-plesset perturbation theory with an analytic slater-type geminal[END_REF] and generalized it to complex values of λ. Since we have found the method of [START_REF] Molin | Multi-precision computation of the complex error function[END_REF][START_REF] Al | Computation of the complex error function using modified trapezoidal rules[END_REF] to be very fast for the computation of error functions, we found that some of the complexity caused by the avoidance of these computations in [START_REF] Ten-No | New implementation of second-order møller-plesset perturbation theory with an analytic slater-type geminal[END_REF] were not needed, and we use a slightly simplified implementation, using an upward recurrence relation (s1 in [START_REF] Ten-No | New implementation of second-order møller-plesset perturbation theory with an analytic slater-type geminal[END_REF]) for large values of ar, and a Taylor expansion in ar (s3 in [START_REF] Ten-No | New implementation of second-order møller-plesset perturbation theory with an analytic slater-type geminal[END_REF]) for small values.

I = - 1

Figure 1 :

 1 Figure 1: Photoionization of Helium, with standard basis sets (left) and even-tempered basis sets (right).The result of the standard Casida method with damping (η = 0.5) in the d-aug-cc-pvqz basis set is shown for comparison. Reference data (black line) from the method of[START_REF] Schwinn | Photoionization and core resonances from range-separated density-functional theory: General formalism and example of the beryllium atom[END_REF].

Figure 2 :

 2 Figure 2: Orbital variations δψ + and δφ + (see text) for Helium at ω = 1.

Figure 3 :

 3 Figure 3: Photoionization of Beryllium (left), with zoom on the 1s→2p resonance (right).

  Figure 4: H2

Figure 5

 5 Figure 5: CH4

  Figure 6: LiH

Figure 7 :

 7 Figure 7: Li +

  a = αβ/(α + β) the combined Gaussian exponent, r = |R t -R s | and where erfcx denotes the scaled complementary error function: erfcx(z) = e z 2 erfc(z).