C Antoine

E Bruand

M Guitreau

J-L Devidal

Understanding preservation of primary signatures in apatite by comparing matrix and zircon-hosted crystals from the Eoarchean Acasta Gneiss Complex (Canada)

Keywords: Acasta, apatite, zircon, Archean, primary composition

A novel way to investigate the petrogenesis of ancient poly-metamorphosed terranes is to study zircon-hosted mineral inclusions, which are sensitive to melt evolution such as apatite. Recent contributions on such inclusions in unmetamorphosed granitoids can provide valuable petrogenetic information and, in turn, represent a way to circumvent effects of metamorphism. Yet, the impact of metamorphism on apatite inclusion has never been studied in detail. To address the issue of chemical and isotopic preservation of primary signals in apatite crystals both in the matrix and armored within zircons, we have studied apatite crystals from four 3.6-4.0 Ga TTG granitoids from the Acasta Gneiss Complex (Canada). Our results demonstrate that U-Th-Pb isotope systematics in matrix apatite crystals were reset at 1.8-1.7 Ga (Wopmay orogen) whereas primary REE signatures were preserved in many crystals. In contrast, zircon-hosted apatite inclusions all preserved primary REE signatures despite variable ages between 1.7 and 4.0 Ga. We interpret reset ages to be a consequence of metamorphism that managed to affect U-Th-Pb systematics because of advanced radiation damage accumulation in host-zircon lattices. Only the most pristine zircon crystal has an apatite inclusion with a concordant age consistent with the magmatic age of the zircon (4.0 Ga). In addition, our results show that apatite crystals from TTG have distinct REE composition from post-Archean granitoids apatites, that is preserved even in some apatites with reset ages. This capacity to retain primary information and discriminate granitoid types makes apatite a very valuable tool for reconstructing the nature and evolution of ancient crustal rocks through the use of detrital minerals.

Introduction

Rare-Earth bearing minerals such as apatite (Ca5(PO4)3(OH, Cl, F)), have proven to be extremely valuable in retrieving the nature of their host magmas and the crystallization history of granitoids [START_REF] Chu | Apatite Composition: Tracing Petrogenetic Processes in Transhimalayan Granitoids[END_REF][START_REF] Jennings | Characterization of magma from inclusions in zircon: Apatite and biotite work well, feldspar less so[END_REF][START_REF] Bruand | An apatite for progress: Inclusions in zircon and titanite constrain petrogenesis and provenance[END_REF]. Hence, they can be very useful complements to zircon in crustal evolution studies (e.g., [START_REF] Iizuka | Improvements of precision and accuracy in in situ Hf isotope microanalysis of zircon using the laser ablation-MC-ICPMS technique[END_REF][START_REF] Belousova | The growth of the continental crust: constraints from zircon Hf-isotope data[END_REF][START_REF] Dhuime | A change in the geodynamics of continental growth 3 billion years ago[END_REF]. Apatite is ubiquitous in different rock types, and contains measurable amounts of a large variety of trace elements (e.g., LREE-MREE, Sr, Pb, Mn, halogens; e.g., [START_REF] Prowatke | Trace element partitioning between apatite and silicate melts[END_REF][START_REF] Gregory | Tracing the evolution of calc-alkaline magmas: in-situ Sm-Nd isotope studies of accessory minerals in the Bergell Igneous Complex, Italy[END_REF][START_REF] Miles | Apatite: A new redox proxy for silicic magmas?[END_REF]. The concentration of trace elements in apatite is sensitive to melt evolution (in-situ crystal fractionation, mixing, SiO2 and ASI contents; e.g., [START_REF] Sha | Apatite chemical composition, determined by electron microprobe and laser-ablation inductively coupled plasma mass spectrometry, as a probe into granite petrogenesis[END_REF][START_REF] Belousova | Trace-element signatures of apatites in granitoids from the Mt Isa Inlier, northwestern Queensland[END_REF][START_REF] Prowatke | Trace element partitioning between apatite and silicate melts[END_REF][START_REF] Zirner | Rare earth elements in apatite as a monitor of magmatic and metasomatic processes: The Ilímaussaq complex, South Greenland[END_REF]. Furthermore, apatite is also datable via U-Pb geochronology although prone to resetting during metamorphism, which makes it problematic to recover meaningful information about the magmatic history of old crustal segments. One way to overcome this limitation is to study apatite inclusions armored within zircon and, thus, combine the potential of apatite in retrieving magmatic information with the resilience of zircon to metamorphic processes. Previous workers have shown that zircon-hosted apatite inclusions can indeed retain important information about its host-rock [START_REF] Jennings | Characterization of magma from inclusions in zircon: Apatite and biotite work well, feldspar less so[END_REF][START_REF] Bruand | An apatite for progress: Inclusions in zircon and titanite constrain petrogenesis and provenance[END_REF][START_REF] Emo | Evidence for evolved Hadean crust from Sr isotopes in apatite within Eoarchean zircon from the Acasta Gneiss Complex[END_REF]. For instance, [START_REF] Bruand | An apatite for progress: Inclusions in zircon and titanite constrain petrogenesis and provenance[END_REF] have studied an unmetamorphosed suite of granitoids and showed that zirconhosted inclusions shared the same trace element signatures as cogenetic apatite crystals within the matrix. They have also shown that Sr/Sm ratio in apatite can discriminate felsic magmas from mafic ones and, as previously suggested by [START_REF] Belousova | Trace-element signatures of apatites in granitoids from the Mt Isa Inlier, northwestern Queensland[END_REF], that Sr content in apatite inclusions can be used as a proxy for the degree of differentiation of their parental magma (evaluation of the SiO2 and Sr contents).

Very recently, [START_REF] Emo | Evidence for evolved Hadean crust from Sr isotopes in apatite within Eoarchean zircon from the Acasta Gneiss Complex[END_REF] demonstrated that zircon-hosted and matrix apatite crystals in two 3. 6-3.7 Ga orthogneisses from the Acasta Gneiss Complex (AGC, Canada) exhibit a general discrepancy in their Sr isotope signatures. The authors interpreted the Sr isotopic signature in apatite inclusions as the first direct estimate of initial 87 Sr/ 86 Sr for these ancient rocks whereas Sr isotope signatures in matrix apatite were reset during later metamorphic event. This interpretation is in line with recent findings of [START_REF] Fisher | Disturbances in the Sm-Nd isotope system of the Acasta Gneiss Complex-Implications for the Nd isotope record of the early Earth[END_REF] who demonstrated that Sm-Nd isotope systematics in matrix apatite from Acasta gneisses were reset during the Wopmay orogen at about 1.85 Ga. However, it is presently unclear to which extent zircon-hosted apatite inclusions can resist metamorphic resetting, especially in ancient zircon crystals that experienced significant radiation damage such as those in the AGC (e.g. [START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF]. Consequently, the aim of this paper is to study the preservation potential of trace element and U-Th-Pb isotope signatures in zircon-hosted inclusions and, hence, better constrain the extent to which such approach can allow retrieval of information about the formation of Earth's first crusts. For this purpose, we carefully and systematically measured trace element contents and U-Th-Pb isotope systematics in zirconhosted apatite as well as in matrix apatite crystals in well-characterized 3. 6-3.96 Ga samples from the Acasta Gneiss Complex [START_REF] Guitreau | Hafnium isotope evidence from Archean granitic rocks for deep-mantle origin of continental crust[END_REF][START_REF] Guitreau | Lu-Hf isotope systematics of the hadean-Eoarchean Acasta gneiss complex (northwest territories, Canada)[END_REF][START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF][START_REF] Mojzsis | Component geochronology in the polyphase ca. 3920 Ma Acasta Gneiss[END_REF].

Geological background and sample description

The Acasta Gneiss Complex (AGC) is located on the western edge of the Slave craton (Northwest Territories, Canada; Fig. 1; e.g. [START_REF] Bowring | 3.96 Ga gneisses from the Slave province, Northwest Territories, Canada[END_REF]. It is bordered to the west by the Wopmay orogen (1900( -1700 Ma; Ma;[START_REF] Bleeker | The Slave Craton: geology and metallogenic evolution. Mineral deposits of Canada: a synthesis of major deposit-types, district metallogeny, the evolution of geological provinces, and exploration methods[END_REF]. The key outcrop locality of the AGC covers an area of ~30 km 2 although its true extent is presently unknown (e.g. [START_REF] Iizuka | Geology and zircon geochronology of the Acasta Gneiss Complex, northwestern Canada: new constraints on its tectonothermal history[END_REF]. This latter is divided into a western and an eastern domain by a major northeast-trending shear-zone with the western part being highly-strained and the eastern part being less deformed [START_REF] Iizuka | Geology and zircon geochronology of the Acasta Gneiss Complex, northwestern Canada: new constraints on its tectonothermal history[END_REF]Reimink et al., 2016aReimink et al., , 2016b[START_REF] Reimink | Petrogenesis and tectonics of the Acasta Gneiss Complex derived from integrated petrology and 142Nd and 182W extinct nuclidegeochemistry[END_REF]. The AGC is mostly made of orthogneisses of dioritic to granitic affinity but some mafic-ultramafic gneisses are also present. [START_REF] Bowring | 3.96 Ga gneisses from the Slave province, Northwest Territories, Canada[END_REF][START_REF] Iizuka | Geology and zircon geochronology of the Acasta Gneiss Complex, northwestern Canada: new constraints on its tectonothermal history[END_REF][START_REF] Guitreau | Lu-Hf isotope systematics of the hadean-Eoarchean Acasta gneiss complex (northwest territories, Canada)[END_REF][START_REF] Mojzsis | Component geochronology in the polyphase ca. 3920 Ma Acasta Gneiss[END_REF][START_REF] Koshida | Petrology and geochemistry of mafic rocks in the Acasta Gneiss Complex: Implications for the oldest mafic rocks and their origin[END_REF]Reimink et al., 2016b).

The earliest magmatic events that occurred in the AGC had been dated between 3962 ± 3 [START_REF] Bowring | 3.96 Ga gneisses from the Slave province, Northwest Territories, Canada[END_REF] and 4031 ± 3 Ma [START_REF] Bowring | Priscoan (4.00-4.03 Ga) orthogneisses from northwestern Canada[END_REF], and recently revised to 4019. 6 ± 1.8 Ma by (Reimink et al., 2016b) The AGC, therefore, hosts Earth's oldest rocks but also exhibits evidence for a crustal pre-history as argued by the presence of a 4200 Ma inherited zircon core [START_REF] Iizuka | 4.2 Ga zircon xenocryst in an Acasta gneiss from northwestern Canada: Evidence for early continental crust[END_REF] and by signatures of both short-and long-lived isotope systematics (e.g., [START_REF] Bowring | The Earth's early evolution[END_REF][START_REF] Iizuka | Improvements of precision and accuracy in in situ Hf isotope microanalysis of zircon using the laser ablation-MC-ICPMS technique[END_REF][START_REF] Iizuka | Reworking of Hadean crust in the Acasta gneisses, northwestern Canada: evidence from in-situ Lu-Hf isotope analysis of zircon[END_REF][START_REF] Guitreau | Lu-Hf isotope systematics of the hadean-Eoarchean Acasta gneiss complex (northwest territories, Canada)[END_REF][START_REF] Mojzsis | Component geochronology in the polyphase ca. 3920 Ma Acasta Gneiss[END_REF][START_REF] Roth | Combined147, 146Sm-143,142 Nd constraints on the longevity and residence time of early terrestrial crust[END_REF][START_REF] Bauer | Coupled zircon Lu-Hf and U-Pb isotopic analyses of the oldest terrestrial crust, the >4.03 Ga Acasta Gneiss Complex[END_REF][START_REF] Emo | Evidence for evolved Hadean crust from Sr isotopes in apatite within Eoarchean zircon from the Acasta Gneiss Complex[END_REF][START_REF] Reimink | Petrogenesis and tectonics of the Acasta Gneiss Complex derived from integrated petrology and 142Nd and 182W extinct nuclidegeochemistry[END_REF];

The AGC evolution was marked by several geological events that started out with the formation of the Idiwhaa unit at 4020 Ma (Reimink et al., 2016b) and ended around 1700-1800 Ma with a pervasive metamorphic event related to the Wopmay orogen [START_REF] Hoffman | Birthdate for the Coronation paleocean: age of initial rifting in Wopmay orogen, Canada[END_REF][START_REF] St-Onge | Wopmay orogen revisited: Phase equilibria modeling, detrital zircon geochronology, and U-Pb monazite dating of a regional Buchan-type metamorphic sequence[END_REF]Table A.1). The igneous activity of the AGC essentially consists in several granitic pulses that occurred between 4020 and 2900 Ma, with a major event at 3600 Ma that involved both juvenile addition and crustal anatexis (Table A.1). The metamorphic history of the AGC is less constrained and mostly derived from U-Pb dating of zircon, apatite and titanite (e.g., [START_REF] Hodges | 40Ar39Ar thermochronology of isotopically zoned micas: Insights from the southwestern USA proterozoic orogen[END_REF][START_REF] Bleeker | The Acasta gneisses: an imperfect sample of Earth's oldest crust[END_REF][START_REF] Sano | Palaeoproterozoic thermal events recorded in the∼ 4.0 Ga Acasta gneiss, Canada: Evidence from SHRIMP U-Pb dating of apatite and zircon[END_REF][START_REF] Iizuka | Geology and zircon geochronology of the Acasta Gneiss Complex, northwestern Canada: new constraints on its tectonothermal history[END_REF][START_REF] Mojzsis | Component geochronology in the polyphase ca. 3920 Ma Acasta Gneiss[END_REF][START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF][START_REF] Fisher | Disturbances in the Sm-Nd isotope system of the Acasta Gneiss Complex-Implications for the Nd isotope record of the early Earth[END_REF][START_REF] Kirkland | Modelling U-Pb discordance in the Acasta Gneiss: Implications for fluid-rock interaction in Earth's oldest dated crust[END_REF]. However, with the available dataset, it clearly appears that during the AGC geological history, most metamorphic events are closely linked to the magmatic activity (e.g., [START_REF] Bleeker | The Acasta gneisses: an imperfect sample of Earth's oldest crust[END_REF]. Major identified metamorphic events occurred at 3750-3700 Ma, 3400-3300 Ma, 3200-2900 Ma, and the latest important metamorphic event identified is linked to the Wopmay orogeny (1800-1700 Ma) with Wopmay associated P-T estimates ranging from 3.3-9 kbar and 575-745 °C [START_REF] Sano | Palaeoproterozoic thermal events recorded in the∼ 4.0 Ga Acasta gneiss, Canada: Evidence from SHRIMP U-Pb dating of apatite and zircon[END_REF][START_REF] St-Onge | Wopmay orogen revisited: Phase equilibria modeling, detrital zircon geochronology, and U-Pb monazite dating of a regional Buchan-type metamorphic sequence[END_REF]. Some more discrete events also likely occurred afterwards (e.g. [START_REF] Hodges | 40Ar39Ar thermochronology of isotopically zoned micas: Insights from the southwestern USA proterozoic orogen[END_REF][START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF].

The samples studied in this contribution (AG09009, AG09014, AG09015, and AG09016) are orthogneisses located in the eastern domain of the AGC (Fig. 1) and have been already described and characterized in [START_REF] Guitreau | Hafnium isotope evidence from Archean granitic rocks for deep-mantle origin of continental crust[END_REF] and [START_REF] Mojzsis | Component geochronology in the polyphase ca. 3920 Ma Acasta Gneiss[END_REF]. Samples AG09014, AG09015, and AG09016 have major and trace element compositions that are similar to typical TTG series (e.g., low K2O/ Na2O, high La/Yb, and strong depletion in heavy REE; e.g. [START_REF] Moyen | Forty years of TTG research[END_REF] and they fall within the field of tonalite and trondhjemite in a normative An-Ab-Or diagram [START_REF] Barker | Trondhjemite: definition, environment and hypotheses of origin, Developments in petrology[END_REF]. Sample AG09009 is more potassic (K2O= 3.64 wt %; K2O/ Na2O ~ 1) despite high La/Yb, and strong depletion in HREE, and therefore falls within the field of granite in a normative An-Ab-Or plot. Samples AG09009, AG09014, and AG09015 were dated by zircon U-Pb geochronology at 3600 Ma [START_REF] Guitreau | Hafnium isotope evidence from Archean granitic rocks for deep-mantle origin of continental crust[END_REF], whereas AG09016 was dated at 3947 Ma [START_REF] Guitreau | Hafnium isotope evidence from Archean granitic rocks for deep-mantle origin of continental crust[END_REF][START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF].

Analytical Methods

Sample preparation

Samples were crushed with a jaw crusher and sieved to recover the fraction < 400 µm. This fraction was subsequently run on a Wilfley table and heavy mineral further sorted using heavy liquids.

Recovered zircon crystals were handpicked, mounted in epoxy resin, and polished to their approximate center. Zircons and apatites in mounts and thick sections were imaged for backscattered electrons (BSE) and cathodoluminescence (CL) using a scanning electron microscope (SEM) JEOL JSM-5910 LV at Laboratoire Magmas et Volcans (LMV; Clermont-Ferrand, France). The accelerating voltage was set to 15 kV during image acquisition.

Electron probe microanalysis (EPMA)

Apatite crystals occurring as zircon-hosted inclusions and matrix grains were analyzed for major and minor elements with an EPMA CAMECA SX100 at LMV. Analyses were made with an accelerating voltage of 20kV and a 10 nA beam current. The crystals used for analysis were PET for Ca, Ce and S, TAP for Si, Na and Y, LPET for P, Cl, La and Sr, and PC0 for F and O. Apatite analysis was done with a defocused beam of 5 μm to minimize potential halogen migration [START_REF] Pyle | Electron microprobe analysis of REE in apatite, monazite and xenotime: protocols and pitfalls[END_REF]. Durango apatite reference material was used as a secondary standard for quality control during analysis [START_REF] Marks | The volatile inventory (F, Cl, Br, S, C) of magmatic apatite: An integrated analytical approach[END_REF].

Trace element analysis by LA-ICP-MS

Trace element analysis of zircon-hosted apatite inclusions was done by laser ablation-inductively coupled plasma mass spectrometry (LA-ICP-MS) using the Resonetics Resolution M-50E 193nm excimer laser coupled to the Thermo Element XR at LMV. The spot sizes used ranged between 4 and 12 μm (depending on the phases analyzed) with a repetition rate of 1 Hz and a fluence of 2.88 J/cm². We purposely set a low frequency to avoid drilling through apatite too quickly. Each run lasted 2 minutes with 20 seconds of background acquisition. Each sequence was bracketed by two analyses of GSE-1G standards (9 μm spot size) at the beginning and at the end of the sequence, and one GSC-1G in the middle of the sequence [START_REF] Jochum | Chemical Characterisation of the USGS Reference Glasses GSA-1G, GSC-1G, GSD-1G, GSE-1G, BCR-2G, BHVO-2G and BIR-1G Using EPMA, ID-TIMS[END_REF]. Hf content was monitored during analysis of apatite inclusions in order to identify potential concurrent sampling of the host zircon. All analyses with Hf content above the background values were discarded.

Matrix apatites were analyzed by LA-ICP-MS using the same laser as for inclusions but coupled to an Agilent 7500 ICP-MS with a spot size of 40 μm, a frequency of 4 Hz, and a fluence of 4.5 J/cm². The same amount of time as for inclusion apatites was used for data acquisition and each run started and finished with two measurements of NIST 610 standard [START_REF] Pearce | A compilation of new and published major and trace element data for NIST SRM 610 and NIST SRM 612 glass reference materials[END_REF].

For all runs (matrix and inclusions), Durango apatite [START_REF] Marks | The volatile inventory (F, Cl, Br, S, C) of magmatic apatite: An integrated analytical approach[END_REF] was regularly analyzed throughout the analytical session to monitor data quality at the same spot size as the analyzed unknowns to account for downhole fractionation (Table A.2). Internal normalization for LA-ICP-MS data was done using 43 Ca and data processing was conducted using the GLITTER software [START_REF] Griffin | GLITTER: data reduction software for laser ablation ICP-MS[END_REF].

U-Th-Pb geochronology by LA-ICP-MS

Apatite crystals were dated by LA-ICP-MS at LMV using a Resonetics Resolution M-50E laserablation system coupled to a Thermo Element XR ICP-MS. During U-Th-Pb measurements, the laser frequency was set to 2 Hz, the fluence to 2.5 J.cm -2 and the spot size between 40 and 9 μm so as to adapt to the size of the crystals targeted for analysis. Analyses consisted in ~20s background counting followed by 60s of signal acquisition in peak-jumping mode (5 ms dwell time and 20 μs settling time). Bahia and Xuxa apatite were analyzed for quality check [START_REF] Rosa | Geochronology (U-Pb/Pb-Pb) and isotopic signature (Rb-Sr/Sm-Nd) of the Paleoproterozoic Guanambi batholith, southwest Bahia State (NE Brazil)[END_REF][START_REF] Schuch | Caracterização de apatitas da Província Borborema, NE do Brasil, como materiais de referência para geocronologia U-Pb via[END_REF]. The apatite standard MAD [START_REF] Thomson | Routine low-damage apatite U-Pb dating using laser ablation-multicollector-ICPMS[END_REF] was used for external normalization. Unknowns were analyzed as batches of 6 and were bracketed by analyses of MAD apatite. Data were processed using in-house spreadsheets based on the t0 intercept method [START_REF] Chew | U-Pb and Th-Pb dating of apatite by LA-ICPMS[END_REF][START_REF] Cochrane | High temperature (> 350 C) thermochronology and mechanisms of Pb loss in apatite[END_REF]. Because we do not know the true 207 Pb/ 206 Pb of MAD apatite, which is a combination of radiogenic Pb and common Pb, and that it is not reported in [START_REF] Thomson | Routine low-damage apatite U-Pb dating using laser ablation-multicollector-ICPMS[END_REF], one cannot determine the mass bias for 207 Pb/ 206 Pb during LA-ICP-MS analyses. Consequently, we derived the Pb isotope mass bias by using the 208 Pb/ 206 Pb expected from the "true" 232 Th/ 238 U and the age of MAD apatite. Similarly, we estimated the amount of common Pb within analyzed apatite crystals using the difference between the mass-bias corrected 208 Pb/ 206 Pb and that expected from the mass-bias corrected 232 Th/ 238 U in a way akin to that of [START_REF] Chew | U-Pb and Th-Pb dating of apatite by LA-ICPMS[END_REF]. Common Pb corrections were further done using the Pb evolution model of [START_REF] Stacey | Approximation of terrestrial lead isotope evolution by a twostage model[END_REF] and a convergence method similar to that of [START_REF] Thomson | Routine low-damage apatite U-Pb dating using laser ablation-multicollector-ICPMS[END_REF] except that we monitored 206 Pb/ 238 U and 208 Pb/ 232 Th for convergence. U-Th-Pb age calculations were done using a present-day 238 U/ 235 U of 137.818 [START_REF] Hiess | 238 U/ 235 U Systematics in Terrestrial Uranium-Bearing Minerals[END_REF] and decay constants for 238 U, 235 U and 232 Th of 1.55125.10 -10 , 9.8485.10 -10 , and 4.9475.10 -11 , respectively [START_REF] Jaffey | Precision measurement of half-lives and specific activities of U 235 and U 238[END_REF][START_REF] Le Roux | Half-life of 232Th[END_REF]. The Isoplot software [START_REF] Ludwig | User's manual for Isoplot 3.6: A geochronological toolkit for Microsoft Excel[END_REF]) was used to build Tera-Wasserburg plots.

The 2056 ± 3 Ma Bahia [START_REF] Rosa | Geochronology (U-Pb/Pb-Pb) and isotopic signature (Rb-Sr/Sm-Nd) of the Paleoproterozoic Guanambi batholith, southwest Bahia State (NE Brazil)[END_REF] and 573.0 ± 0.8 Ma Xuxa [START_REF] Schuch | Caracterização de apatitas da Província Borborema, NE do Brasil, como materiais de referência para geocronologia U-Pb via[END_REF] apatite standards were analyzed using 40, 20, 15, and 9 μm spot sizes. The results for apatite standards are detailed in Table A.3 andFig 2048 ± 31 (2σ; MSWD = 2.7; n = 30), 2086 ± 34 (2σ; MSWD = 1.5; n = 13), and2127 ± 40 Ma (2σ; MSWD = 2.7; n = 13), respectively, for Bahia, and 569.4 ± 2.9 (2σ; MSWD = 1.1; n = 20), 574.9 ± 4.4 (2σ; MSWD = 0.39; n = 17), and 575.1 ± 3.2 Ma (2σ; MSWD = 0.97; n = 17), respectively, for Xuxa. Therefore, our measurements are in good agreement with previously reported ages.

As illustrated in Fig. A.1, common Pb correction results in migration of data points towards the Concordia curve but not complete converge to a concordant age, especially for Bahia apatite data. This can be explained by our correction method, which is well suited for minerals with low Th/U but becomes less precise and less accurate as Th/U increases as a result of radiogenic 208 Pb/ 206 Pb getting closer to that of common Pb (Ireland and Williams, 2003). Consequently, the influence of common Pb becomes less detectable in high Th/U mineral. As an illustration, a 10% contamination of a 2056 Ma Bahia apatite (Th/U = 15), a 573 Ma Xuxa apatite (Th/U = 30), and a 1700 Ma Acasta apatite (Th/U = 0.5) decreases determined 208 Pb/ 206 Pb ratios by 5 and 7 % in Bahia and Xuxa apatites, respectively, while increasing the 208 Pb/ 206 Pb by almost 150% in the Acasta apatite. This observation can explain why, later in the manuscript, the common Pb corrections work much better for Acasta apatites, which commonly have low Th/U, as most data points become concordant once corrected for common Pb.

Apatite crystals from Acasta gneisses were analyzed in three different contexts: (1) as separate grains within an epoxy mounts (40, 20, and 9 μm spots on the same apatite crystals), (2) as grains within polished sections (matrix apatite), and (3) as inclusions within zircon crystals mounted in epoxy. All results are available in Table A.4.

Petrography

Sample AG09009

Sample AG09009 is mainly made of quartz, alkali feldspar, plagioclase, biotite, muscovite, chlorite and accessory phases such as zircon, apatite, monazite, allanite, thorite, ilmenite and rutile (Fig. 2A). Locally, chlorite partly replaces micas and can occur as small veins. Felspars are slightly sericitized. In this sample, apatite is the main accessory phase and most commonly occurs as euhedral crystals (<150 microns in length) either directly in contact with biotite or as inclusions in quartz and feldspar. Apatite can be zoned from light gray (core) to dark gray (rim) in cathodoluminescence (Fig. 2D) and can contain zircon inclusions. Apatite occasionally appears as secondary mineral associated with allanite and monazite. In these rare cases, a monazite grain is surrounded by a corona of apatite, itself surrounded by a corona of fibrous allanite (Fig. 2B).

Numerous small inclusions of monazite are persistent inside the apatite corona.

Biotite, alkali feldspar, plagioclase, chlorite, apatite, chromite and chalcopyrite have been identified as inclusions in zircon (Fig. 2C) and are usually < 20 microns. Most zircon crystals are metamict (damaged by radioactivity), which leads to the formation of fractures that can be filled by several minerals (e.g. biotite, chlorite, plagioclase) thereby forming polymineralic inclusions.

Most apatite inclusions are monomineralic and euhedral except for one case where apatite is surrounded by chlorite (Fig. 2C).

Sample AG09014

Sample AG09014 is made of quartz, biotite, plagioclase and accessory phases. These latter are zircon, apatite, pyrite, monazite and rutile crystals. In this sample, apatite is the main accessory phase (Fig. 2I-J), it is euhedral, small (< 100 microns) and can contain zircon inclusions. Apatite can be found either as inclusion in quartz and plagioclase or in contact with biotite. Apatite is usually unzoned (Fig. 2I-J). In this sample, inclusions in zircon are monomineralic and small (< 10 microns). They are mostly apatite (Table 1), quartz, chlorite and monazite crystals.

Sample AG09015

In sample AG09015, the main rock forming minerals are quartz, biotite, plagioclase and oxides (rutile, ilmenite; Table 1). Accessory phases are zircon, apatite, sulfides (pyrite, chalcopyrite), titanite and allanite. As for samples AG09014, AG09015 is only slightly altered. Apatite is the main accessory phase; it is euhedral and < 200 microns. It is either in direct contact with biotite or as inclusion within quartz and feldspar. Most crystals are zoned with a light gray core and a dark gray rim in CL (Fig. 2H). They can also be found as clusters where zonation is oscillatory with a dark gray core (Fig. 2G).

In this sample, inclusions found in zircon are apatite, biotite, plagioclase and chlorite.

These inclusions are small (< 5 microns) and monomineralic.

Sample AG09016

Sample AG09016 is made of quartz, plagioclase, biotite and muscovite. Accessory phases are zircon, apatite, rutile and titanite. This sample shows limited alteration. Apatite is ubiquitous, euhedral and < 250 microns in size and occurs in contact with biotite or in inclusions in quartz and plagioclase. Apatite crystals have two types of zonation: (i) a light gray core with a dark rim or (ii) an oscillatory zoning with a dark core (Fig. 2F). Some apatite grains are occasionally fractured with biotite filling these cracks.

In this sample, inclusions in zircon can be either mono-or polymineralic with grain sizes < 20 microns. The identified inclusions are apatite (Fig. 2E), chlorite, biotite, plagioclase, quartz and thorite. For this sample, apatite inclusions are always monomineralic.

Results

Apatite chemistry

Apatite crystals from the studied samples are all fluorapatites (Table A.5). Based on EPMA data, apatite chemistries define two groups: a first group (samples AG09009 and AG09016) with high Y2O3 (0.04 to 0.6 wt %) and low Cl (Cl <0.045 wt %) and a second group (AG09014 and AG0915) with lower Y2O3 (<0.2 wt %) and higher Cl (0.048-0.2 wt%; Table A .5).

Contrary to samples AG09009 and AG09016, only matrix apatite has been analyzed in AG09014 and AG09015 samples because inclusions were too small (< 10 µm in size). 2) whereas group 2 patterns are sub-parallel to group 1 but with lower REE content and no Eu anomaly. Apatite inclusions have the same REE patterns as group 1 with some grains having a higher REE concentration (Fig. 3A,4A). Overall, apatites from AG09009 have the highest concentrations in La (~166 ppm) and Yb (~353 ppm, Table A.4) among all samples analyzed in this study (Fig. 4A). These highest values are systematically associated with analyses done within crystal cores whereas lower concentrations are linked to analyses done within rims (Fig. 4B-C). On the other hand, group 2 apatite compositions have only been found in the outer-part of grains that had core compositions belonging to group 1 apatites (Fig. 2D) and have the lowest REE concentrations (Fig. 4A).

Sample AG09014. Chondrite-normalized REE patterns for apatite from sample AG09014

show similar LREE depletion compare to HREE. In this sample, two groups of apatite can be distinguished. The first group (Group 1) is similar to group 1 identified in other samples from this study (Fig. 3C). The second group is depleted in all REE compared to group 1 and has a distinctive impoverishment in LREE.

Sample AG09015. Contrary to previous samples, apatite from AG09015 displays multiple and strongly distinguishable compositions. Five groups can be defined based on chondritenormalized patterns (Fig. 3D). Group 1 REE pattern is comparable to group 1 from the other samples with a depletion in LREE compared to HREE (Fig. 4A) and a negative Eu anomaly (0.4< Eu/Eu*<0.6). All other groups have a strong depletion in LREE, which distinguishes them from group 1. In Fig. 4A, it can be seen that the other groups (2-5) have low La and Yb contents and are undistinguishable from each other (gray dots in Fig. 4D). These groups have highly depleted LREE chondrite-normalized patterns (LaN<100) compared to HREE (YbN<1000). The best chemical criterion that allows these groups to be discriminated is their respective Eu anomaly (Fig. 3D).

Group 2 has no Eu anomaly (0.9<Eu/Eu*<1.3), Group 3 a slightly negative Eu anomaly (Eu/Eu*≈0.5) and Group 4 a positive Eu anomaly (Eu/Eu*>1.2). Group 5 apatite has slightly higher LREE chondrite-normalized values compared to group 4, and a small negative Eu anomaly (Eu/Eu*≈0.75). Group 2, 3 and 5 are always neoformed (no zonation, low REE content) or recrystallized grains, while group 4 can be found as new generation grains or as an overgrowth surrounding group 2 apatites. Sample AG09016. Much as for sample AG09009, chondrite-normalized REE patterns of matrix apatites from AG09016 define two groups. Group 1 is depleted in LREE compared to HREE (mean (La/Yb)N = 0.14) with a consistent negative Eu anomaly (~0.3; Fig. 3B,Table A.4.).

Group 2 defines a sub-parallel pattern to group 1, with lower REE content, a stronger depletion in LREE compared to HREE (c.a. (La/Yb)N=0.03) and small Eu anomalies (0.6 <Eu/Eu*<0.7).

Chondrite-normalized REE patterns of zircon-hosted apatite inclusions have strong similarities with group 1 apatite, although the LREE content in the former is more variable (Fig. 3B).

Apatite U-Th

Zircon-hosted apatite

Over the large quantity of zircons mounted in epoxy and examined within thin sections (~200), only a few apatite inclusions larger than 10 μm were found. One apatite was analyzed in sample AG09009 whereas eight were analyzed in sample AG09016. Most zircon crystals in which we found apatite inclusions have non-pristine internal textures, except for zircon AG09016-88 which exhibits fine-oscillatory zoning. Ages for 6 of the 9 host-zircons are reported in Figure 6. The 3 remaining zircon crystals were not dated because of their disturbed internal zoning which would have resulted in meaningless discordant ages. Apatite ages mentioned in this section correspond to data corrected for common Pb unless specified otherwise. Overall, U-Pb isotope systematics in apatite inclusions exhibit contrasting behavior compared to matrix apatite. The former essentially distribute along the Concordia curve between 1700 and 4000 Ma whereas the latter forms a robust trend between 1700-1800 Ma and common Pb (Fig. 7). In detail, zircon-hosted apatite inclusions display three distribution behaviors: (i) on or close to the Concordia curve, (ii) on the matrix apatite trend, and (iii) inverse discordance.

The oldest inclusion (AG09016-88) has a concordant age of 4003 ± 89 Ma ( 207 Pb/ 206 Pb age corrected for common Pb; 2 SE; 0.7% discordance; Fig. 7) which is consistent with data obtained on the pristine-looking host-zircon (3967 ± 49 My, [START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF]. A second analysis in the rim of the same inclusion shows a reversely discordant age, though consistent with the first one (4111 ± 90 Ma). Two other data points from apatites AG09016-100 and 108 plot on the left of the Concordia curve at 207 Pb/ 206 Pb ages corresponding to ~4000 Ma. Apatites AG09016-24 and 37 show concordant ages between 3200 and 4000 Ma. Apatite AG09016-15 is the biggest inclusion and is not fully enclosed within the host-zircon (Fig. 6). The zircon exhibits a dark patchy zoning indicative of advanced radiation damage accumulation. In this apatite inclusion, 3 analyses from core to rim were done. Obtained U-Th-Pb ages decrease progressively from core (~3200 My) to rim (~1800 My old; Fig. 7). Apatite AG09016-62 displays two ages consistent with the oldest domains of apatite AG09016-15 (i.e., 3200 Ma and 2400 Ma). The single data point for AG09009 (08) corresponds to an age of 1867 ± 143 Ma ( 207 Pb/ 206 Pb age; 2 SE) and comply well with the major trend defined by matrix apatites, much as inclusion AG09016-11 (Fig. 7). It is important to note that inclusion AG09009-08 is surrounded by chlorite and chlorite-filled cracks, suggesting later fluid circulation within the zircon crystal (Fig. 6). Similarly, apatite AG09016-11 is traversed by a major crack that runs through the zircon.

Discussion

Preservation of magmatic chemical signature

In all studied TTG samples, apatite exhibits at least two chemical groups of which Group 1 has the same characteristics across the entire range of samples. Group 1 apatites are characterized by chondrite-normalized REE patterns depleted in LREE compare to HREE and a negative Eu anomaly (Fig. 3). In a Nd vs Eu/Eu* diagram, group 1 apatites (both matrix and inclusions) have systematically lower Eu/Eu* values and higher Nd content than that of other apatite groups (Fig. 8). Texturally, this group corresponds to core composition when a core to rim systematic is present (Fig. 4B,C). Moreover, all zircon-hosted inclusions have group 1 chemical characteristics. Finally, recent data on apatite from adakite (TTG-like signature; Bruand et al., 2020 and TTGs from two other cratons (Kaapvaal, South Africa and Karelia, Finland) show comparable chemistry to group 1 apatite. All these results indicate that group 1 apatite composition represents primary (magmatic) apatite REE signature.

Matrix apatite

Our new U-Th-Pb ages on matrix apatite indicate that the vast majority of analyzed crystals exhibit ages that are consistent with the latest pervasive metamorphic event recorded by Acasta gneisses, i.e., the 1800-1700 Ma Wopmay orogen (Fig. 5). This is in line with apatite ages obtained by [START_REF] Sano | Palaeoproterozoic thermal events recorded in the∼ 4.0 Ga Acasta gneiss, Canada: Evidence from SHRIMP U-Pb dating of apatite and zircon[END_REF] and [START_REF] Fisher | Disturbances in the Sm-Nd isotope system of the Acasta Gneiss Complex-Implications for the Nd isotope record of the early Earth[END_REF] and indicate that metamorphic conditions associated with the Wopmay orogen exceeded the closure temperature of apatite (> 350-600°C e.g. [START_REF] Cherniak | Lead diffusion in apatite and zircon using ion implantation and Rutherford Backscattering techniques[END_REF][START_REF] Chamberlain | Apatite-feldspar U-Pb thermochronometer: a reliable, mid-range (∼450°C), diffusion-controlled system[END_REF][START_REF] Cochrane | High temperature (> 350 C) thermochronology and mechanisms of Pb loss in apatite[END_REF]. We can, therefore, assume that any previously recorded metamorphic events had been overprinted by that related to the Wopmay orogen (Table A.1). Rare concordant matrix apatite ages at 3197 ± 82 and 2721 ± 63 Ma can point to possible local preservation of older ages, in line with titanite geochronology [START_REF] Davidek | Conventional U-Pb geochronology of the Acasta gneisses using single crystal fragmentation techniques[END_REF][START_REF] Fisher | Disturbances in the Sm-Nd isotope system of the Acasta Gneiss Complex-Implications for the Nd isotope record of the early Earth[END_REF]; Fig. 5), but overall magmatic ages are not preserved in matrix apatites despite the local preservation of primary REE signature in the same apatites (group 1). The dichotomy between the primary chemical signature of group 1 matrix apatite and the secondary record of its U-Th-Pb data is consistent with diffusion work done on REE by [START_REF] Cherniak | Rare earth element diffusion in apatite[END_REF] who showed that REE in apatite can diffuse relatively slowly even in the case of high temperature metamorphic events. This study demonstrates that two different types of substitution are responsible for diffusion in apatite; simple substitution (REE 3+ <->REE 3+ ) and coupled substitutions (i.e. REE 3+ + Si 4+ <-> Ca 2+ + P 5+ ). In diffusion experiments, [START_REF] Cherniak | Rare earth element diffusion in apatite[END_REF] showed that Sr, Pb and O are more readily reset in apatite but can keep their primary signature in the core of the crystal (10 μm zoning) for metamorphic episodes below 750°C and for a duration of 1.05 . 10 4 years. On the other hand, apatite grains (100μm in size) keep their primary REE signature for metamorphic episodes of less than 10 5 years in simple substitutions and 10 7 years in coupled substitutions for the same temperature.

Apatite inclusions

Zircon-hosted apatite inclusions within AG09016, all exhibit primary REE signatures and, contrary to matrix apatites, essentially distribute along the Concordia curve with ages ranging between ~1700 and 4000 Ma (Fig. 7; Table A .6). This demonstrates that apatite inclusions have not been completely overprinted by the Wopmay orogen but, on the contrary, kept the memory of older metamorphic events and even original igneous U-Pb ages in some cases. This is precisely the case for apatite AG09016-88 which displays a concordant age of 4003 ± 89 Ma ( 207 Pb/ 206 Pb age, 0.7% discordant) that is very consistent with that obtained on the host-zircon (3967 ± 49 Ma; [START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF]. The large range of ages exhibited by apatite inclusions is consistent with the fact that the AGC recorded several tectono-thermal events from 3750 to 1700 Ma, although it is difficult to associate all ages to specific events (Table 2). Yet, three data points correspond to an age of ~3200 Ma consistent with that of titanite and garnet within Acasta samples [START_REF] Maneiro | Development of a detrital garnet geochronometer and the search for Earth's oldest garnet[END_REF][START_REF] Fisher | Disturbances in the Sm-Nd isotope system of the Acasta Gneiss Complex-Implications for the Nd isotope record of the early Earth[END_REF]. We interpret the general distribution to represent punctual thermally activated resetting during metamorphic events (e.g., [START_REF] Cochrane | High temperature (> 350 C) thermochronology and mechanisms of Pb loss in apatite[END_REF]. Three other data points are consistent with the matrix apatite trend, which suggest that the U-Pb system in these crystals were reset in a similar way as matrix apatite, despite encapsulation within zircon. Therefore, apatite inclusions have been either reset within zircon without any influence from the external media, or in connection with this latter via cracks and porous zircon domains.

Reversely discordant data can be seen as either analytical artifacts or evidence for Pb contamination, which would significantly lower the measured U/Pb ratio. An analytical artifact would mean odd fractionation between U and Pb, which is inconsistent with the relatively good agreement between U, Th, and Pb concentrations determined during trace element analyses and U-Th-Pb isotope measurements. In addition, such extreme displacement from the Concordia curve would mean an instrumental fractionation much too large compared to the rest of our U/Pb data.

Therefore, we think that the U-Pb isotope signatures of these inclusions are real and come from Pb contamination of apatite inclusions. Common Pb cannot be responsible for this contamination because the 208 Pb/ 206 Pb of reversely discordant apatite inclusions is too low. Figure 9 shows that matrix apatite Pb isotope composition is controlled to a large extent by common Pb contamination because data points form a positive correlation towards the common Pb end-member. In contrast, apatite inclusions have a horizontal distribution and, hence, no visible influence from common Pb, except for a few apatites that exhibited U-Pb data distribution consistent with the matrix apatite trend in Figure 7. Therefore, we conclude that common Pb contamination cannot account for the reversely discordant apatite data. Instead, we propose that reversely discordant apatite inclusions were contaminated by radiogenic Pb coming from the host zircon during a thermal event such as the Wopmay orogen. Radiogenic Pb migration has already been described in Archean zircons, notably from the Napier Complex, Antarctica [START_REF] Williams | Unsupported radiogenic Pb in zircon: a cause of anomalously high Pb-Pb, U-Pb and Th-Pb ages[END_REF][START_REF] Kusiak | Mobilization of radiogenic Pb in zircon revealed by ion imaging: implications for early Earth geochronology[END_REF][START_REF] Lyon | Pb nanospheres in ancient zircon yield model ages for zircon formation and Pb mobilization[END_REF]. The evolution of the Pb isotope composition of a zircon with a Th/U of 0.5 is presented in Figure 9 and it is clearly visible that contamination by such radiogenic Pb would result in horizontal spread of data in a 208 Pb/ 206 Pb versus 207 Pb/ 206 Pb plot.

The zircons that host apatite inclusions analyzed in this study are often not pristine-looking (Fig. 6). Instead, they show advanced signs of radiation damage, a common feature in Archean zircon, have variably discordant ages, and some of them even display metamorphic ages (See BSE images in Fig. 6). Because zircons underwent disturbances in their U-Pb isotope systematics by thermal events [START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF], one can easily expect that apatite inclusions might also be affected. This would be consistent with our geochronological results. An important conclusion from our apatite dating is that compared to matrix apatite, no common Pb seems to influence the data distribution in AG09016 (Table A .6. and Fig. 9). This means that host-zircons satisfactorily isolated apatite inclusions from external fluids possibly bringing common Pb, which accounts for both the resetting of U-Pb ages, the existence of reversely discordant U-Pb data and the fact that Sr isotope data for apatite inclusions from [START_REF] Emo | Evidence for evolved Hadean crust from Sr isotopes in apatite within Eoarchean zircon from the Acasta Gneiss Complex[END_REF] are less radiogenic than those in the matrix despite zircons showing advanced radiation damage in CL images. If external fluids were to penetrate within zircon, apatite Sr isotope composition should be influenced by the fluids, which are more radiogenic. In summary, our results highlight that U-Th-Pb ages can be compromised in apatite inclusions hosted in zircons with advanced radiation damage but this is not the case for REE concentrations.

Secondary apatite signatures

The secondary groups of apatites are chemically distinct from the primary group described above.

All of them show a strong decrease in REE compared to group 1 apatites (or magmatic apatites), with HREE depleted to a lesser extent than LREE. All apatites from groups 2 and more (except AG09014), display no or more rarely a positive Eu anomaly (Fig. 3D). These groups texturally appear as new or recrystallized grains or are rimming group 1 cores. Their growth are related to secondary processes such as metamorphism and/or metasomatism which in the light of current literature on those phases are sometimes difficult to distinguish [START_REF] Bau | Rare-earth element mobility during hydrothermal and metamorphic fluid-rock interaction and the significance of the oxidation state of europium[END_REF][START_REF] Henrichs | The trace element and U-Pb systematics of metamorphic apatite[END_REF].

In the studied samples, we identified different late apatite generations:

(i) In one sample (AG09009), a core (group 1) to rim (group 2) systematic has been identified (Fig. 4B). The depletion in LREE and the absence of Eu anomaly characterizing the rim is interpreted as metamorphic (Fig. 3A). This type of overgrowth is also observed in metamorphosed post-Archean apatites and has similar low REE+Y, U+Th rims around igneous apatite cores [START_REF] Henrichs | The trace element and U-Pb systematics of metamorphic apatite[END_REF].

(ii) For the three other samples, second generation of apatite crystals (Fig. 3 group 2 to 5 depending on samples have been observed. In samples AG09015 and AG09016, these groups are interpreted as metamorphic and/or metasomatic. On the other hand, group 2 from sample AG09014 is the only second generation displaying a strong negative Eu anomaly and rather magmatic characteristics. This could be related to a late magmatic interaction with a different magma or the consequence of the crystallization of another LREE-mineral during the crystallization of the granitoid [START_REF] Hoskin | Identifying accessory mineral saturation during differentiation in granitoid magmas: an integrated approach[END_REF][START_REF] Bruand | Accessory mineral chemistry of high Ba-Sr granites from northern Scotland: constraints on petrogenesis and records of whole-rock signature[END_REF]; Table 2). Here it probably points to a lower grade of metamorphism (than for AG09015) because of its higher content in Cl and lower content in Y (Table A.5.).

(iii) Apatites derived from primary monazite. These secondary apatites have only been identified in sample AG09009 (Fig. 2B). In this case, primary monazite is destabilized in a metamorphic corona of apatite and fibrous allanite. This type of secondary textures has been described occasionally in natural samples [START_REF] Gasser | The behaviour of monazite from greenschist facies phyllites to anatectic gneisses: An example from the Chugach Metamorphic Complex, southern Alaska[END_REF] in metapelite under amphibolite facies; [START_REF] Finger | Replacement of primary monazite by apatite-allanite-epidote coronas in an amphibolite facies granite gneiss from the eastern Alps[END_REF] metamorphosed I and S-type granites under amphibolite facies) and in experimental studies [START_REF] Budzyń | Experimental constraints on the relative stabilities of the two systems monazite-(Ce)-allanite-(Ce)-fluorapatite and xenotime-(Y)-(Y, HREE)-rich epidote-(Y, HREE)-rich fluorapatite, in high Ca and Na-Ca environments under PT conditions of 200-1000 MPa and 450-750[END_REF], and is characteristic of amphibolite facies metamorphic conditions.

All secondary apatite compositions have U-Th-Pb ages corresponding to the Wopmay orogen 1800-1700 Ma and have nucleated or were reset during this event (Fig. 5). Again, this is consistent with young ages obtained by [START_REF] Sano | Palaeoproterozoic thermal events recorded in the∼ 4.0 Ga Acasta gneiss, Canada: Evidence from SHRIMP U-Pb dating of apatite and zircon[END_REF] and [START_REF] Fisher | Disturbances in the Sm-Nd isotope system of the Acasta Gneiss Complex-Implications for the Nd isotope record of the early Earth[END_REF].

Apatite chemistry and identification of host-granitoid

Archean TTGs are known for their specific geochemical signature, notably their high (La/Yb)N and low HREE concentrations (Fig. 10C; e.g., [START_REF] Moyen | Forty years of TTG research[END_REF]. They can be easily discriminated from post-Archean granitoids in (La/Yb)N vs YbN diagrams. In this study, we found that the same diagram discriminates very well TTG apatite from Post-Archean granitoid apatite (Fig. 10A) alike in Bruand et al. 2020. However, the chemical characteristics define an opposite trend when compared to the whole rock data. Indeed, apatites from TTGs have systematically low (La/Yb)N ratios, high HREE contents and notably high Y concentrations, whereas apatites from post-Archean granitoids commonly show the opposite. Interestingly, apatite chemistry from peraluminous post-Archean granitoids are marginally overlapping with the apatite from TTG (Fig 10a-d). Their (La/Yb)N are nonetheless systematically higher than apatite TTG from Archean cratons (Fig. 10D). Yet, a very important conclusion of our study is that we observed distinct generations of apatite that we attribute to magmatic, late-stage, metamorphic and/or metasomatic processes. The late events affecting the AGC have an impact on the LREE chemistry of apatite whereas the HREE signature is similar. The consequence of this, is that the (La/Yb)N ratio of secondary apatite in TTG will be even lower than the magmatic one (Fig. 10B). However, no matter which generation analyzed apatites belong to, they all show the peculiar chemistry of TTG apatite in a (La/Yb)N vs YbN diagram, hence, allowing identification of their source rock (Fig. 10). This signature is not specific to Acasta samples as shown by Bruand et al. (2020). Apatites analyzed in TTG from Kaapvaal and Karelia cratons (South-Africa and Finland respectively) have the same chemical characteristics.

Sr/Y vs Y diagram is another well-known discriminating diagram for TTG whole-rock. However, this diagram only moderately discriminates the TTG from other granitoids when using apatite chemistry. Indeed, Sr and Y contents in apatite have been shown to be strongly influenced by SiO2, Sr and the ASI of the melt [START_REF] Bea | Behavior of accessory phases and redistribution of Zr, REE, Y, Th, and U during metamorphism and partial melting of metapelites in the lower crust: an example from the Kinzigite Formation of Ivrea-Verbano, NW Italy[END_REF][START_REF] Sha | Apatite chemical composition, determined by electron microprobe and laser-ablation inductively coupled plasma mass spectrometry, as a probe into granite petrogenesis[END_REF]Belousova et al., 2002;[START_REF] Jennings | Characterization of magma from inclusions in zircon: Apatite and biotite work well, feldspar less so[END_REF][START_REF] Bruand | Accessory mineral chemistry of high Ba-Sr granites from northern Scotland: constraints on petrogenesis and records of whole-rock signature[END_REF]. Consequently, in a Sr/Y vs Y diagram, apatite from peraluminous and TTG granitoid are undistinguishably defining the same field.REE partitioning in apatite are intimitaly thlinked to the host rock ASI, substitutions and co-crystallizing phases. To explain the dichotomy between the LREE depleted and enriched HREE characters of TTG apatite and the LREE-enriched and HREE-depleted characters of TTG (whole rock; Fig. 10 A-C), we need to invoke the presence of other LREE-bearing phases. Recently, Bruand et al. 2020 have highlighted that the presence of monazite is not rare in TTG granitoids. In the studied samples, only few occurrences of monazite interpreted as primary have been identified. The occurrences of monazite in one of the samples suggest that the presence of monazite could be partially responsible for the specific chemistry of TTG apatite. In post Archean granitoids cocrystallizing monazite would mean a reduced host magma with ASI >1 and relatively low Ca content [START_REF] Hsieh | Petrogenesis of the Nanling Mountains granites from South China: Constraints from systematic apatite geochemistry and whole-rock geochemical and Sr-Nd isotope compositions[END_REF].Trace elements allow making the distinction between mafic I-type granitoids and I/S-type granitoids [START_REF] Sha | Apatite chemical composition, determined by electron microprobe and laser-ablation inductively coupled plasma mass spectrometry, as a probe into granite petrogenesis[END_REF] but the different types of granites can only be discriminated by using Sr and Th. The concentrations in both elements are rather different from the ones from post-Archean apatite, which shows that the composition of apatite in TTG are most likely due to different melt composition and probably fO2 during crystallization of the apatite in the Archean granitoid. However, the competition between REE-bearing phases in TTG are not well known and the absence of partition coefficient for those REE-bearing phases in such granitoid preclude further interpretation of the data but highlight that more work is needed on this subject.

Implications

The results of the present study demonstrate that matrix apatites in polymetamorphosed Acasta gneisses can preserve primary REE signatures, whereas U-Th-Pb isotope systematics have been fully reset. On the contrary, zircon-hosted apatite inclusions can both preserve primary REE signatures and U-Th-Pb isotope systematics provided that the host-zircon experienced limited radiation damage and remained isolated from the external medium. Interestingly, apatite can be used as a discriminating tool for different granitoid types (This study; [START_REF] Jennings | Characterization of magma from inclusions in zircon: Apatite and biotite work well, feldspar less so[END_REF][START_REF] Bruand | Accessory mineral chemistry of high Ba-Sr granites from northern Scotland: constraints on petrogenesis and records of whole-rock signature[END_REF][START_REF] Bruand | An apatite for progress: Inclusions in zircon and titanite constrain petrogenesis and provenance[END_REF]Bruand et al., , 2020)). Therefore, our results strongly suggest that the geochemistry and geochronology of detrital apatite and/or detrital zircon-hosted apatite inclusions can be used as a new tool to unravel the origin, nature and evolution of complex multistage Archean cratons as well as the type of granitoid Hadean detrital zircons crystallized from. However, one needs to be aware of the limitations of this approach that essentially pertain to metamorphic resetting of isotope systematics. This limits the use of detrital apatite but also suggests that when one focuses on zircon-hosted inclusions, the host-zircon should be well characterized for, at least, internal textures and U-Pb ages so as to understand the pristineness of the information encapsulated in the apatite inclusion. This is because all inclusions within zircon are not systematically primary but could grow within old radiation-damaged detrital zircons (e.g. [START_REF] Rasmussen | Metamorphic replacement of mineral inclusions in detrital zircon from Jack Hills, Australia: Implications for the Hadean Earth[END_REF]. Consequently, analyzed inclusions should be monomineralic to avoid perturbation from later events and should be contained within a zircon with limited radiation-damage.

Conclusions

The aim of this study was to decipher whether or not, matrix apatites and apatite inclusions in zircons from Archean TTGs from the Acasta Gneiss Complex (Canada) preserved primary compositions. Here, we were able to demonstrate that primary REE composition of apatite in the matrix can be preserved for large crystals (in this study mostly over 40 µm), whereas their U-Th-Pb isotope systematics were fully reset by HT metamorphism around 1.8-1.7 Ga (Wopmay Orogen). Our study also shows that apatites from TTGs have geochemical signatures (La/Yb)N vsYbN) that are very distinct from those of apatite from other granitoid type, regardless of their primary or secondary origin. In addition, our results show that Archean apatites trapped within zircon that have withstood multiple HT metamorphic events were effectively protected and kept their primary isotopic and geochemical signature, provided that zircon lattice underwent moderate radiation damage. Altogether, these new results have important implications for studies of Archean and Hadean crustal evolution because this tool can be applied to detrital apatites, knowing the limitations of such an approach, but more importantly to detrital zircon-hosted apatites. In the latter case, it is important to target the most pristine-looking zircon crystals, to characterize them for, at least, U-Pb ages, and look for monomineralic apatite inclusions in priority.

Acknowlegdements

We are grateful to Claire Fonquernie for mineral separation and to Jean-Marc Hénot for technical help during SEM sessions. The authors acknowledge financial support from LabEx ClerVolc (ANR-10-LABX-0006), Région Auvergne (Auvergne Fellowship program), and the French Agence Nationale de la Recherche (ANR) through the funded project Zircontinents (PI: M.

Guitreau; ANR-17-CE31-0021). This is a Laboratory of Excellence ClerVolc contribution number 404.. We thank Steve Mozjsis and The John Templeton Foundation-FfAME Origins program in the support of CRiO. Finally, we deeply thank the editor W. Behr, T. Iizuka and an anonymous reviewer for their input. containing apatite inclusions analyzed in this study. Also reported are locations and labels of U-Th-Pb isotope measurements together with zircon ages, when available [START_REF] Guitreau | Hafnium isotope evidence from Archean granitic rocks for deep-mantle origin of continental crust[END_REF][START_REF] Guitreau | Crystallization and Disturbance Histories of Single Zircon Crystals From Hadean-Eoarchean Acasta Gneisses Examined by LA-ICP-MS U-Pb Traverses[END_REF]. Pb plot for apatite data from this study. Also reported are the evolution of the terrestrial lead, as given in [START_REF] Stacey | Approximation of terrestrial lead isotope evolution by a twostage model[END_REF], and the Pb isotope evolution of a 4 Ga zircon that has a present-day Th/U of 0.5. Matrix and zircon-hosted apatites exhibit a clear dichotomy with the former defining a positive trend towards common Pb whereas the latter distribute horizontally, except for a few data points that have U-Pb ages consistent with those of matrix apatite (see Fig. 6). Note that apatite 207 Pb/ 206 Pb ratios vary between values that correspond to the age of the Wopmay orogen and that of the oldest Acasta gneisses, except for three data points that we surmise incorporated zircon-derived radiogenic Pb during metamorphism related to Wopmay orogen (see text for details). and post-Archean granitoids. Literature data from [START_REF] Hoskin | Identifying accessory mineral saturation during differentiation in granitoid magmas: an integrated approach[END_REF][START_REF] Belousova | Trace-element signatures of apatites in granitoids from the Mt Isa Inlier, northwestern Queensland[END_REF][START_REF] Dempster | Magmatic zoning in apatite: a monitor of porosity and permeability change in granites[END_REF][START_REF] Chu | Apatite Composition: Tracing Petrogenetic Processes in Transhimalayan Granitoids[END_REF][START_REF] Jennings | Characterization of magma from inclusions in zircon: Apatite and biotite work well, feldspar less so[END_REF][START_REF] Bruand | Accessory mineral chemistry of high Ba-Sr granites from northern Scotland: constraints on petrogenesis and records of whole-rock signature[END_REF]Bruand et al., , 2020)). B) Details on the distinction between magmatic and metamorphic apatites from this study. C) (La/Yb)N vs YbN discriminating diagram between TTG and post-archean rocks modified from [START_REF] Hansen | Crustal melting in the lower parts of island arcs: an example from the Bremanger Granitoid Complex, west Norwegian Caledonides[END_REF]. D) Log(La/Yb)N vs YbN discriminating diagram between TTG and post-archean granitoids.

Figure Captions

  . A.1. With all spot sizes combined, we obtained U-Pb Concordia intercept and weighted average 206 Pb * / 238 U and 208 Pb * / 232 Th ages of

  Sample AG09009. Chondrite-normalized REE patterns of apatites from the matrix define two groups: group 1 shows depleted LREE compared to HREE and a negative Eu anomaly (Eu/Eu*= EuN/√[(SmN).(Gd)N]< 0.5; Fig. 3A, Table

  Most of matrix apatite crystals in AG09009 define a major trend that intersects the Concordia curve in a Tera-Wasserburg plot (Fig.5A) at 1797 ± 27 Ma (MSWD = 3.6; n = 27). Once corrected for common Pb (red ellipses in Fig.5A), most data points plot onto the Concordia curve with a weighted average 206 Pb * / 238 U age of 1794 ± 24 Ma (MSWD = 1.4). One data point lingers outside the major trend and has a Concordia age of 3197 ± 82 Ma (MSWD = 0.09) and 3184 ± 120 Ma (MSWD = 0.11; Fig.5A) before and after common Pb correction, respectively. Data for AG09014 apatite crystals form a robust trend with a lower intercept at 1700 ± 15 Ma (MSWD = 1.6; n = 12; Fig.5B). Once corrected for common Pb, most data points plot on the Concordia curve and exhibit a weighted average 206 Pb * / 238 U age of 1710 ± 20 Ma (MSWD = 0.92).With the exception of one data point, all analyses from AG09015, show a trend that gives a lower-intercept age of 1713 ± 21 Ma (MSWD = 4.1; n = 35; Fig.5C). Once corrected for common Pb, all data points plot on or very near the Concordia curve and give a weighted average 206 Pb * / 238 U age of 1723 ± 20 Ma (MSWD = 1.9). One data point plots on the Concordia and gives an age of 2721 ± 63 Ma.All data points for apatites from AG09016 plot along a linear trend towards a lowerintercept of 1746 ± 46 Ma (MSWD = 3.7; n = 37; Fig.5D). Once corrected for common Pb, most data points are concordant and have a weighted average 206 Pb * / 238 U age of 1721 ± 23 Ma (MSWD = 2.2).

Figure 1 :

 1 Figure 1: Map of the AGC with the location of the samples for this study (red dot). Map modified

Figure 2 :

 2 Figure 2: BSE and CL images of the studied samples A) Representative picture of AG09009

Figure 3 :

 3 Figure 3: Chondrite-normalized REE patterns (normalization factor :[START_REF] Sun | Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes[END_REF]

Figure 4 :

 4 Figure 4: A) La vs Yb diagram for all samples with first and secondary groups of matrix apatites.

Figure 5 :

 5 Figure 5: Tera-Wasserburg diagrams presenting U-Pb results for matrix apatites from AG09009

Figure 6 :

 6 Figure 6: Cathodoluminescence (CL) and back-scattered-electron (BSE) images for zircons

Figure 7 :

 7 Figure 7: Tera-Wasserburg diagram for apatite inclusions in AG09009 and AG09016 corrected for

Figure 8 :

 8 Figure 8: Discriminating Nd vs Eu/Eu* diagram for all samples, matrix and inclusions.

Figure 9 :

 9 Figure 9: 208 Pb/ 206 Pb versus 207 Pb/ 206 Pb plot for apatite data from this study. Also reported are the

Figure 10 :

 10 Figure 10: A) (La/Yb)N vs YbN discriminating diagram for apatites from Archean granitoids (TTG)

Table Captions Table 1 :

 Captions1 Recapitulative table of the sample characteristics.

Table 2 :

 2 Representative LA-ICP-MS analysis for all apatite samples from the matrix and in inclusions within zircons.


Agathe Reille
email: agathe.reille@ensam.eu

Nicolas Hascoet
email: nicolas.hascoet@ensam.eu

Chady Ghnatios
email: cghnatios@ndu.edu.lb

Amine Ammar
email: amine.ammar@ensam.eu

Elias Cueto
email: ecueto@unizar.es

Jean Louis Duval

Francisco Chinesta
email: francisco.chinesta@ensam.eu

Roland Keunings
email: roland.keunings@uclouvain.be

Incremental dynamic mode decomposition: A reduced-model learner operating at the low-data limit

Keywords:

The present work aims at proposing a new methodology for learning reduced models from a small amount of data. It is based on the fact that discrete models, or their transfer function counterparts, have a low rank and then they can be expressed very efficiently using few terms of a tensor decomposition. An efficient procedure is proposed as well as a way for extending it to nonlinear settings while keeping limited the impact of data noise. The proposed methodology is then validated by considering a nonlinear elastic problem and constructing the model relating tractions and displacements at the observation points.

Introduction

In physics in general and in engineering in particular, when addressing a generic problem, the first step consists in selecting the best suited model for the description of the evolution of the system under consideration. Nowadays, a variety of models exists that are well established and validated, covering most of the domains of physics, e.g., solid and fluid mechanics, electromagnetism, heat transfer... These models are based on the combination of two types of equations of very different nature. From one side, the so-called balance equations, whose validity in the framework of classical physics is out of discussion. The second type of equations are the so-called constitutive models. These relate primal variables and dual variables, e.g., stress and strain, temperature and heat flux. Their expression and validity depend on the nature of the considered system and the applied loading-in the most general sense.

Constitutive equations are often phenomenological, and in general involve some parameters assumed to be intrinsic to the system under consideration.

For this reason, and even in the context of classical engineering, before solving a given problem, one must decide on the most appropriate model. There obviously exists a risk on the pertinence of the chosen model for addressing the phenomena under study. By performing some engineered experiments in order to calibrate the model, that is, to extract the value of the parameters that it involves, these risks can be alleviated.

The calibrated model is expressed generally as a system of coupled nonlinear partial differential equations and must be solved for given loading and boundary conditions.

Very often, the solution of these problems is not tractable by an analytical procedure. Numerical methods were introduced in the twentieth century for that purpose, and are nowadays widely employed, in industry as well as in academia. Thus, approximate procedures were proposed, most of them computationally manipulable. The first option consists in assuming the solution expressible as a combination of a reduced number of functions with a physical or mathematical foundation. The Ritz method follows this rationale. Thus, if these functions are noted by G i (x), i = 1, . . . , G, the solution approximation reads

u(x, t) ≈ G i=1 α i (t)G i (x) (1)
where the time-dependent coefficients α i (t) are calculated by invoking a projection method, like the Galerkin method, for instance.

However, the knowledge of these functions becomes a tricky issue in many cases. When the choice of an appropriate, physically sound basis is unavailable, the best option is considering a general purpose approximation basis, for example consisting of polynomials (Lagrange, Chebyshev, Legendre, Fourier...): P i (x), i = 1, . . . , P, from which the approximation now reads u(x, t) ≈

P i=1 β i (t)P i (x) (2)
and where, for ensuring accuracy, a sufficient number of terms must be considered, implying a quite large number of terms in the finite sum [START_REF] Loève | Probability Theory[END_REF].

In order to better adapt to complex domains, global functions P i (x) are often replaced by compactly supported functions N i (x), leading to the so-called Finite Element Method, FEM. If we consider N nodes, the approximation reads

u(x, t) ≈ N i=1 U i (t)N i (x) (3)
where now U i (t) represents the searched solution at node x i and at time t. This formulation becomes general enough at the price of computing many-in some cases too many-unknowns, the nodal values U i (t), ∀i.

In the sequel, for the sake of notational simplicity, the dependence on time is not made explicit, and it is assumed that the discrete system that results from introducing approximation (3) into the weak form of the problem (for the sake of simplicity assumed to be linear) leads to the linear system KU = F (4)

Here, matrix K represents the discrete form of the model, whereas vectors U and F contain the nodal unknowns and loads, respectively. Their respective sizes are N × N, N × 1 and N × 1.

Reduced modelling

In many applications of practical interest, despite the richness of Eq. ( 4) that is able to consider any choice of the loading in a vector space of dimension N, usual loadings contain in practice many correlations. This results in the solution U being defined in a subspace of dimension n, with n N.

Proper Orthogonal Decomposition proceeds by embedding the solution onto that subspace of dimension n

[1] [2] [3] [4].
Nonlinear manifold learning constitutes its nonlinear counterpart. For that purpose, from a set of snapshots of the solution U 1 , . . . , U S , a reduced basis R i (x), i = 1, ..., n, with n N is extracted. The solution is finally expressed by

u(x, t) ≈ n i=1 γ i (t)R i (x) (5)
whose matrix form reads

U = Bγ (6)
The columns of matrix B contain the reduced functions at node locations. Thus, the component B ij = R j (x i ) with i = 1, . . . , N and j = 1, ..., n. Now, by injecting ( 6) into (4) and premultiplying by the transpose of B (something equivalent to a Galerkin projection in the reduced basis), it results

(B KB)γ = B F (7)
whose respective sizes are n × n, n × 1 and n × 1.

In the nonlinear case, different procedures have been proposed for alleviating the construction of K. Among them, one can find the Discrete Empirical Interpolation Method [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF] or the Hyper-Reduction [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF], to cite but a few.

In order to avoid the necessity of performing an "a priori" learning, a simultaneous searching of space and time functions was proposed in the context of the so-called Proper Generalized Decomposition, PGD [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF] [8] [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations[END_REF] [10] [START_REF] Cueto | Proper Generalized Decompositions: An Introduction to Computer Implementation with Matlab[END_REF]. It was then extended for addressing space separation [START_REF] Bognet | Separated representations of 3d elastic solutions in shell geometries[END_REF] or parametric solutions [START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF] [14] [START_REF] Díez | Generalized parametric solutions in Stokes flow[END_REF] [START_REF] Chinesta | On the use of proper generalized decompositions for solving the multidimensional chemical master equation[END_REF] in some of our numerous previous works. Thus, the space-time separated representation now reads

u(x, t) ≈ D i=1 T i (t)X i (x) (8)
where both groups of functions, T i and X i , are obtained by injecting the approximation (8) into the problem weak form and then using a rank-one enrichment for incrementally constructing the separated representation. At each enrichment step of that greedy algorithm, an alternated direction strategy is considered for addressing the nonlinearity arising from the product of both unknown functions T i and X i .

Despite its generality, the computed solution (8) depends on the loading term. In the case of POD-based techniques, the reduced algebraic system can easily be inverted, leading to the reduced transfer function (B KB) -1 . Note that the reduced model matrix (B KB) does not depend explicitly on the loading. Thus, as soon as a new loading F comes into play, it is to be projected on the reduced basis B T F ≡ f. Then, the reduced solution is obtained from γ = (B KB) -1 f. To recover the solution in the original finite element basis, one simply has to perform the projection U = Bγ . PGD, on the contrary, does not compute such a reduced model (or its transfer function counterpart), but instead it proceeds by either (i) calculating the solution for each loading, or (ii) expressing the loading in an reduced basis (e.g., POD) according to

f (x, t) = F i=1 μ i F i (x, t) (9)
that allows transforming the initial problem into a parametric one. Within the PGD rationale, the solution is searched in the parametric form u(x, t, μ), with μ = (μ 1 , ..., μ F ).

A criticism that is usually attributed to the PGD technique, if compared to POD-based procedures, is precisely the necessity of the load to be representable in terms of the reduced basis F 1 , . . . , F F . However, it is important to note that even if such a constraint is less explicit when using the POD rationale, it also applies implicitly [START_REF] Gonzalez | Real-time direct integration of reduced solid dynamics equations[END_REF]. Thus, it is important to note that the POD rationale is based on the existence of a reduced basis B, and that this basis was constructed from a particular choice of snapshots U 1 , . . . , U S , associated with a particular loading F 1 , . . . , F S . Thus, if a quite different loading comes into play-one whose solution cannot be accurately approximated by the reduced basis R i involved in B-the transfer function can be applied, but nothing guarantees the accuracy of the resulting solution. Finally, the applicability of Model Order Reduction techniques remains subordinate to the fulfillment of the conditions assumed during their construction.

All the just discussed techniques belong to the vast family of Model Order Reduction techniques. In fact, models were not really reduced, only the solution procedure.

Learning models

We argued in some of our former works that, in some circumstances, models are too uncertain to represent the physical system [START_REF] Ibañez | Data-driven non-linear elasticity: constitutive manifold construction and problem discretization[END_REF]. In that case, a data-driven procedure becomes an appealing choice [19] [20] [21] [START_REF] Kevrekidis | Equation-free modeling[END_REF]. For the sake of simplicity, we assume the problem to be linear and that N couples (U i , F i )-assumed independent, that is, spanning the whole vector space of dimension N-are available. We also assume that nothing is known on the model whose discrete expression consists of the matrix K.

Thus, with all the couples fulfilling the algebraic relationship KU i = F i , ∀i [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] one could construct the vector K (vector expression of the matrix K), the extended loading vector F that includes all the loads, i.e. F = (F 1 , . . . , F N ) and finally matrix Ũ from the different U i organized in an adequate manner so as to enable Eqs. [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] to be expressed as Ũ K = F [START_REF] Cueto | Proper Generalized Decompositions: An Introduction to Computer Implementation with Matlab[END_REF] They can be solved to obtain the model K, or its matrix counterpart K.

There are many ways of performing such an identification. In the linear or the nonlinear cases as, for example, deeplearning (based on neural networks, NN), dynamic mode decomposition, DMD... to cite but a few.

It is important to note that the complexity of standard solutions involving known models K, scales with the number of unknowns N, however, when trying to identify the model the complexity scales with N 2 (number of component of K), justifying that a lot of data is mandatory for extracting the hidden model, situation that becomes even more stringent in the nonlinear case that requires identifying a model K| U around any possible value of U.

This apparent complexity justifies the fact of associating the term big-data to machine learning. However, in engineering, the available data is in some circumstances very scarce: collecting data could be technologically challenging, sometimes simply impossible, and in all cases, generally expensive.

Thus, the big-data is being, or should be, transformed in a smarter counterpart, with smart-data paradigm rationalizing the amount of needed data, while driving its acquisition (collection): which data, where and when?

The present work tries to exploit the fact previously discussed, that despite the richness of loadings and responses (being the model the application that relates both), in general both are living in low-dimensional manifolds, and therefore the model relating both is expected to reflect this fact. This consideration was exploited when proposing the so-called hyper-reduction techniques, however they were applied on pre-assumed models, whereas in what follows we apply it while assuming the model unknown.

After this introduction, the next section proposes a procedure for extracting the reduced model from a small amount of data, then the procedure will be illustrated from the analysis of a quite simple nonlinear problem, before finishing by addressing some general conclusions.

Methods

For the sake of simplicity, we start by assuming the discrete linear problem KU = F [START_REF] Bognet | Separated representations of 3d elastic solutions in shell geometries[END_REF] where F and U represent the input and output vectors respectively. In a general case, they could have different nature and represent the values of their respective fields at the observation points. The respective sizes are N × 1 and N × 1.

As in the case of reduced-order modeling, we assume that inputs and outputs are (to a certain degree of approximation) living in a sub-space of dimension n, much smaller than N. Thus, the rank of K is expected to be also n, even if a priori, it was ready to operate in a larger space of dimension N.

The question is therefore how to extract the reduced model? Among the numerous possibilities we explored, we comment here on two of them.

Rank-n constructor

Here, we consider a set of S input-output couples (F i , U i ), and assume the model to be expressible from its low-rank form

K LR K ≈ K LR = n j=1 C j ⊗ R j = n j=1 C j R j ( 13
)
where ⊗ denotes the tensor product, and C j and R j are the so-called column and row vectors. This expression is somehow similar to the separated representation used in the PGD, the SVD (Singular Value Decomposition) or the CUR decomposition [23] [24].

Remark 1.

If the model is expected to be symmetric, one could enforce that symmetry in the solution representation, i.e.

K ≈ K LR = n j=1 (C j R j + R j C j ) ( 14
)
Let us define the functional E(K LR ) according to

E(K LR ) = S i=1 F i -K LR U i p ( 15
)
The choice of many different norms could be envisaged here. For exploiting sparsity, for instance, the choice should be the

L1-norm.

In what follows, we consider the standard L2-norm, i.e. we take p = 2. By considering matrices F and U containing in their columns vectors F i and U i respectively, the previous expression can be rewritten using the Frobenius norm (related to p = 2) as [START_REF] Chinesta | On the use of proper generalized decompositions for solving the multidimensional chemical master equation[END_REF] whose minimization results in K LR (UU T ) = FU [START_REF] Gonzalez | Real-time direct integration of reduced solid dynamics equations[END_REF] or, equivalently,

E(K LR ) = F -K LR U F
K LR = (FU )(UU T ) -1 (18)
This proves that K LR and, more particularly, its column and row vectors, correspond to the SVD (or rank-n truncated SVD) decomposition of (FU )(UU T ) -1 . Remark 2. Within the standard PGD rationale, one could compute the separated representation of K LR by computing progressively the separated representation using the standard greedy algorithm applied on Eq. ( 17) from K * : K LR (UU T ) = K * : (FU ) [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF] with ":" the twice-contracted tensor product, and K LR approximated by [START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF]. The test matrix K * is obtained from:

K * = n j=1 C * R j ( 20
)
or its symmetrized counterpart.

This procedure is specially suitable when N becomes large, thus making difficult the calculation of the SVD decomposition of (FU )(UU T ) -1 .

Remark 3. This procedure ensures that when addressing a linear problem, N linearly independent loadings F j , j = 1, ..., N, ensure the construction of a rank-N model K. If the model is computed by incorporating progressively the available datafrom one loading leading to model K 1 (a single data model), to N leading to K N -we can define the model enrichment j K = K j -K j-1 , j = 2, ..., N.

Note that adding more data, F j , j > N, to the model does not make it evolve anymore, as expected. In other words: K j = K N , j > N. In the nonlinear case, where a model is a priori expected to exist around each U, the model continues to evolve when j > N.

With the low-rank model thus calculated, as soon as a new datum arrives, F, the solution is evaluated from

U = K LR -1 F (21)
Remark 4. In section 2.3 we propose an alternative methodology based on the calculation of transfer functions that avoids model inversion.

Progressive greedy construction

In this case we proceed progressively. We consider the first available datum, the pair (F 1 , U 1 ). Thus, the first, one-rank, reduced model reads [START_REF] Ch | A rational strategy for the resolution of parametrized problems in the {PGD} framework[END_REF] and one possible solution, ensuring symmetry, consists of R 1 = F 1 and

K 1 = (C 1 R 1 ) (22) ensuring (C 1 R 1 )U 1 = F 1
C 1 = F 1 ϕ 1
, with ϕ i = F 1 U 1 . Many other solutions exist, since there are N available data for computing 2N unknowns (the components of the row and column vectors). Of course, the problem can be written as a minimization problem using an adequate norm.

Suppose now that a second datum arrives, (F 2 , U 2 ), from which we can also compute its associated rank-one approximation, and so on, for any new datum (F i , U i ). This leads to

K i = (C i R i ) (24)
ensuring at its turn

(C i R i )U i = F i ( 25
)
where

C i = F i ϕ i (with ϕ i = F i U i ) and R i = F i .
For any other U, the model could be interpolated from the just defined rank-one models, K i , i = 1, ..., S, according to

K| U ≈ S i=1 K i I i (U) (26)
with I i (U) the interpolation functions operating in the space of the data U. This constructor is not appropriate when addressing linear behaviors. If we assume known two local rank-one behaviors K 1 and K 2 ensuring the fulfillment of relations K 1 U 1 = F 1 and K 2 U 2 = F 2 , it follows that if for example F = 0.5(F 1 + F 2 ), by considering K = 0.5(K 1 + K 2 ), the resulting output does not satisfy linearity, i.e. U = 0.5(U 1 + U 2 ).

However, in the nonlinear case, by defining the secant behavior at the middle point associated with F = 0.5(F 1 + F 2 ) as K = 0.5(K 2 -K 1 ), we will have K = K 1 + K = 0.5(K 1 + K 2 ), fact that allows viewing the progressive greedy construction and its associated interpolation as a linearization procedure.

General remarks

-All the previous analysis was based on the calculation of the reduced model K LR . However, this procedure needs its inversion K LR -1 for evaluating the solution U = K LR -1 F. This issue could be easily circumvented as follows.

The discrete model KU = F can be rewritten as U = K -1 F and by introducing the discrete transfer function T (that represents the inverse of K), one could learn directly, by using all the previous rationale, the reduced expression of the transfer function, i.e. T LR .

The advantage of learning the reduced discrete transfer function is that, as soon as a datum F is available-and under the assumption that it is living in the same subspace that served for constructing the reduced model-the evaluation becomes straightforward: a simple matrix-vector product. It can even be more simplified, since matrix T LR is expressed in a separate format, U = T LR F, with

T ≈ T LR = n j=1 Ĉ j R j ( 27
)
-Nothing really changes when considering nonlinear models, except the fact that the models, K or T and, in turn, the vectors they involve, must be assigned to the neighborhood of the data U or F. In practice, data can be classified, creating a number clusters where the models are constructed [START_REF] Amsallem | An interpolation method for adapting reduced-order models and application to aeroelasticity[END_REF].

Thus, each time a datum F arrives, the cluster to which it belongs-say, κ-is first identified. Then, the low-rank discrete transfer functions of that cluster, T LR κ , is chosen and the solution evaluated according to

U = T LR κ F (28)
-Again within the PGD rationale, all the discussion above can be extended to parametric settings. The description and results constitute a work in progress that will be reported in ongoing publications. Another work in progress concerns the obtention of error bounds, needed for certifying the constructed models and their predictions. There is a vast corps of literature on verification and validation, but they where associated with model-based simulations. Here, for the best identified model, the existing error estimation procedures can be applied. However, the remaining question concerns the estimation of the error associated with the model itself, whose quality depends on the collected data, i.e. validation, more than verification. This issue also applies in model-based simulations: how to be sure that the considered models are the appropriate ones? Despite the apparent difficulty, since the model is defined with the only support of a discrete basis, some estimators could be easily defined, and will be reported in future works. -All the previous developments could be accomplished by first extracting the reduced bases associated with input and output vectors (action and reaction) and then learning the associated reduced models within the framework of a progressive POD. However, the procedure here considered seems more physically sound. A fully reduced counterpart of the procedures here proposed and described constitutes a work in progress, whose results will be reported in ongoing publications. -The present procedure becomes quite close to the so-called dynamic model decomposition as soon as the model is assumed to be expressed in a tensor form and extended to nonlinear settings from an appropriate clustering. -The main issue when considering data-driven models is the noise impact on the predictions. There are many possibilities, most of them based on the use of filters. Here, we considered, as proved later, a simple procedure. The data (F i , U i ) is obtained with a resolution N, on which noise applies. A quite efficient filter consists in simply considering the model, i.e. its row and column vectors, described with a coarser resolution M < N. In fact, inputs and outputs could exhibit fast fluctuations, but in many applications, in which model order reduction applies, these fast fluctuations are almost due to noise, and for that reason, the model could be expected being described with a coarser resolution.

Results

A nonlinear example

In order to prove the potential interest of the methodology here described, we consider a nonlinear elastic model with Poisson coefficient and Young's modulus given by respectively by ν = 0.3 and E = 10 5 (1 + 1000 ε I I ), with ε I I the second invariant of the strain tensor ε. Units are MPa for stresses, N/mm for applied tractions and mm for lengths.

The mechanical problem is defined in the two-dimensional domain depicted in Fig. 1, that also shows the finite element mesh employed for generating the synthetic pseudo-experimental data. Displacements are prevented on its bottom boundary. Left and right boundaries are traction-free and a distributed tension applies on its upper boundary, given by t = (0, a 0 + a 1 x).

Fig. 2 depicts, for a 0 = 100 and a 1 = 0, the displacement field and the Young modulus distribution, while Fig. 3 depicts the displacement along the upper boundary.

From the reference elastic model we explore the parameter space (a 0 , a 1 ) by computing realizations generated from (a 0 + r 0 , a 1 + r 1 ), with r 0 and r 1 two random numbers uniformly distributed in the interval [-1, 1] around points (a 0 = Thus, the nonlinear elastic model was employed for generating the synthetic data. In the present case, it consists of traction and displacement on the upper boundary of the domain. For example, for the loading cluster around (a 0 = 100, a 1 = 0), the input (loading) and output (vertical displacement) data are shown in Fig. 4.

From all the available data, the discrete, low-rank transfer function T LR is calculated at each loading cluster. Fig. 5 depicts the column and row vectors involved in the rank-two model associated with the cluster (a 0 = 100, a 1 = 0). The obtained rank, two, corresponds to the expected one, since the loading is defined in a manifold of dimension 2, consisting of parameters a 0 and a 1 .

When evaluating the model at any point in the parametric domain (a 0 , a 1 ), the model T LR is interpolated from the four computed models

T LR = 4 κ=1 T LR κ I κ (a 0 , a 1 ) (29)
where I κ (a 0 , a 1 ) denotes the bilinear interpolation functions in the parametric space (a 0 , a 1 ) associated with cluster κ.

Influence of noisy data

In the present case the loading was randomly perturbed, as well the computed displacement field (using the unperturbed loading). As mentioned in Section 2, a simple filtering procedure consists in truncating the approximation to a number of M < N terms. Our experience indicates that this simple procedure calculates first those modes with a lower frequency content, thus filtering noise in practice. Fig. 7 compares column and row vectors involved in the discrete low-rank transfer function (whose dimensionality increased due to the noise but was truncated to rank 3) in absence of any noise filtering (i.e. M = N) and when using M N. A performing filter capability is appreciated, with a beneficial effect on the accuracy of the predictions.

Equivalent neural network

It is easy to transform the proposed procedure into an equivalent artificial neural-network (ANN). Fig. 8 compares a standard ANN with one emulating the proposed reduced model learning, combined with the model interpolation for inferring responses from applied inputs. The proposed NN tries to emulate a linear combination of locally linear learned behaviors. As discussed in Section 2.2 the linear combination of behaviors can be viewed as a linearization of a nonlinear behavior.

In the proposed NN four clusters were considered (even if for the sake of clarity the figure only depicts two) representing the four locally linear behaviors (two parameters -a 0 and a 1 -with to classes each), while the output layer ensures a bilinear interpolation of them.

The standard ANN consists of three neuron layers, the input and output involving 66 and the intermediate with 100. The intermediate layer for the physics-informed neural network contained 100 neurons per cluster.

Both NN were trained by 300 input/output vectors generated t = (0, a 0 + a 1 x), with a 0 and a 1 two uniformly distributed random number numbers in the intervals [0, 100] and [-10, 10] respectively. Fig. 9 compares the displacement predictions obtained by the trained neural networks, the standard one (ANN) and the physics-informed one for a 0 = 100 and a 1 = 0, with the reference solution, from which the superiority of the physicsinformed NN can be stressed.

Conclusions

The present work succeeded at proposing a new methodology for learning reduced models from a small amount of data by assuming a tensor decomposition of the discrete reduced model or its transfer function counterpart. It was extended to nonlinear settings and proved that the associated physics-informed NN exhibits a higher accuracy than a standard one.

Our present works in progress concern the validation and verification as well as its extension to parametric settings.

Fig. 1 .

 1 Fig. 1. Domain equipped with a finite-element mesh.

Fig. 2 .

 2 Fig. 2. Displacement field (left) and Young modulus distribution (right).

Fig. 3 .

 3 Fig. 3. Displacement on the upper boundary.

Fig. 4 .

 4 Fig. 4. Top: input (applied tension). Bottom: output (vertical displacement).

Fig. 6

 6 Fig. 6 compares the reference solution and the one obtained by interpolating models at point (a 0 = 50, a 1 = 5), T LR , and then computing the response from U = T LR F.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Column and row vectors defining the discrete low-rank transfer function T LR .

Fig. 7 .

 7 Fig. 7. Comparing columns and row vectors involved in the discrete low-rank discrete transfer function with M = N (top) and M N enabling filtering

Fig. 8 .

 8 Fig. 8. Standard one-layer NN (top) and the physics-informed one related to the low-rank reduced model (bottom).

Fig. 9 .

 9 Fig. 9. Predictions obtained from the standard one-layer ANN and the physics-informed one (Low-Rank Reduced Model) related to the low-rank reduced model for a 0 = 100 and a 1 = 0. Two clusters are depicted in the so-called hidden layer, even if the NN contains four.

Acknowledgements

The work of E.C. has been partially supported by the Spanish Ministry of Economy and Competitiveness through Grant number DPI2017-85139-C2-1-R and by the Regional Government of Aragon and the European Social Fund, research group T88. The other authors acknowledge the ANR (Agence nationale de la recherche, France) through its grant AAPG2018 DataBEST.