Pierre Bordet

Florian Kergourlay

Ariane Pinto

Nils Blanc

Ariane Pinto Ab

Pauline Martinetto

Applying multivariate analysis to X-ray diffraction computed tomography: the study of medieval applied brocades

teaching and research institutions in France or abroad, or from public or private research centers.

ARTICLE

Please do not adjust margins Please do not adjust margins a. Univ. Grenoble Alpes, CNRS, Institut Néel, 38000 Grenoble, France.

Introduction

For almost a century, X-ray powder diffraction has been the technique of choice for the structural investigation of solid matter, starting from simple phase identification to the ab initio structure solution of complex organic molecules, quantification of multiphase mixtures, or investigation of microstructures. These developments have been closely linked to the improved performances of x-ray sources, instruments and analytical methods, being in the laboratory or at dedicated synchrotron facilities. More recently, with the progresses of powerful 3rd generation synchrotron sources offering micron size monochromatic beams of very high intensities, measurements allowing 3D computed tomographic reconstructions from x-ray diffraction data have become possible 1,2,3 . In such XRD-CT experiments, a small pencil x-ray beam is scanned through the sample and at each scan position, XRD patterns are recorded while the sample is rotated over an axis perpendicular to the beam. Reconstruction algorithms allow building 3D graphical representations of the crystalline and amorphous phase distribution in the sample, based on the selection of a characteristic diffraction peak or diffuse scattering signal. The quality of the reconstructions mainly depends on the size of the rotation angular step, and the spatial resolution is determined by the beam size and hence the translation scan step. Allowing to obtain a 3D map of the crystalline as well as amorphous phases in a sample, XRD-CT is complementary to other x-ray based tomographic techniques such as absorptionbased x-ray tomography (X-CT) and fluorescence x-ray tomography (XRF-CT). The former informs on the variations of electron density in the sample, while the latter provides the distribution of chemical elements. Since the appearance of the XRD-CT technique, a number of efficient software has been developed to perform the reconstruction and examination of the measurements 4,5 . The data are generally organized as a 3D matrix (or data cube), named as (T,R,D) hereafter in the text where 2 dimensions represent the coordinates of a scan points in translation (T) and rotation (R) and the 3rd dimension represents the XRD pattern (D) as function of the scattering parameter (e.g. the 2θ angle or Q=4.π.sinθ/λ, λ being the radiation wavelength) measured at each (T, R) point. The data cube can be represented graphically as 2D R-T sinograms for a particular slice in 2θ of the XRD pattern, showing the distribution of the crystalline phases having characteristic Bragg peaks in the selected slice. The same can be also visualized in real space after tomographic reconstruction using well-known algorithms. This procedure allows identifying and locating the different crystalline phases in the sample. Since nowadays 2D diffraction images are generally recorded, a preliminary step consists in azimuthally integrating the data to yield 1D diffraction patterns which can then be inserted in the data This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins cube. In addition, the XRD patterns corresponding to an arbitrary volume of the sample can be obtained by summing the patterns from the voxels chosen in the real space reconstruction. By separating the contributions from the different phases, this so-called reverse analysis provides a way to identify the minor ones with better accuracy than would be feasible on the pattern from the whole sample, and to dramatically improve the detection limit of the technique. Cultural Heritage Sciences are developing a growing interest for the technical skills shown by ancient artists and craftsmen in the fabrication of artworks and artefacts. In the need to identify recipes and material sources and to understand their degradation over time for restoration purpose, the contribution of x-ray based techniques has become more and more important. With the recent advent of mobile XRF and XRD instruments, non-invasive measurements can now be performed in situ, avoiding the need to move the artwork from its location. 2D mapping of artworks with flat surfaces (paintings, frescoes, manuscripts...) now allows observing the distribution of pigments, binders, etc. at micron size resolutions 6,7 . However, in many cases the detailed investigations of the sample design, stratigraphy or degradation damages require obtaining a three-dimensional picture from microsamples carefully selected from the investigated artwork. X-ray based tomographic methods are then very well suited to thoroughly determine the phase composition and stratigraphy of such microsamples in a nondestructive way, leaving them untouched for complementary analyses. They allow, for example, to visualize the full 3D arrangement of pictorial layers, without having to destructively prepare the microsample 8,9 .

With the development of new generation synchrotron machines, the availability to submicron size beams allowing to rapidly perform such kind of experiments at higher resolutions is quite an enthusiastic perspective for the investigation of complex materials such as those from the Cultural Heritage field. However, with the foreseen new beamlines' performances, scanning a 1 mm size sample with 1 µm translation and 360 rotation steps (a quite typical case) yields in a few minutes 360 000 XRD patterns. This overwhelming amount of data has then to be processed in an efficient way in order to extract in a realistic timespan a usable, if possible quantitative information. Furthermore, despite the obvious interest of XRD-CT for Cultural Heritage Science, several technical difficulties still exist which can limit its use. Due to the relatively large gauge size, the voxels may contain more than one phase. It must be noted here that decreasing the gauge size for XRD-CT has a limit, since a sufficient number of grains must be contained in the gauge volume for powder averaging. Micron-size layers are not rare, for example in gilding layers. Hence, some quantitative (e.g. Rietveld) refinement of the XRD patterns calculated for each voxel is necessary to reach a precise, at least semi-quantitative description of the variation of sample composition and structure. Even when a quantitative knowledge of phase proportion is not required, such fitting may be necessary to identify minor phases. If sequential Rietveld refinement would be possible in theory, in practice it seems quite unrealistic given the very large amount of data, the random variation of phase contents (phases can appear and disappear abruptly from one voxel to the next) and the possibly large number of phases. Furthermore, such heavy data processing with no guarantee of rapid success requires a long time and cannot be considered for quality check of the data during experiments. Bearing these problems in mind, we have investigated the use of multivariate analysis techniques for the processing of XRD-CT data collected on microsamples taken from ancient artworks. Here we compare the use of various algorithms, like Principal Component Analysis (PCA), supervised and nonsupervised Multivariate Curve Resolution -Alternating Least-Squares (MCR-ALS) and Non-Negative Matrix Factorization (NMF). Using the latter, we show that a typical data cube can be automatically decomposed in a small number of crystallographically significant components within few minutes. These components which may contain more than one phase, can be submitted to quantitative analysis using Rietveld refinement, leading to a complete description of the 3D phase distribution in the sample. Multivariate analysis techniques have been used recently in Materials Science in a wealth of problems describing the evolution of a chemical system as a function of an external parameter, such as time, temperature, pressure, etc., which involves processing large experimental data sets coming from various types of measurements. For these techniques to be relevant, it must be possible to find a small number of invariable components so that each individual measurement representing a state of the system can be described as a linear combination of the invariable components. Only the concentrations (or weights) of the components vary with the external parameter. The problem is thus to determine from the measurements the nature of the components and their weights for each value of the external parameter. This can be expressed mathematically by the relation: D = CS+E where D is a n.m matrix representing the measurements, n being the number of data sets (spectra, diffractograms, etc.) and m the number of points in each of them (as a function of 2θ, Q, etc.), C is a n.d concentration matrix containing the concentrations of each or the d components for each of the n measurements, S is a d.m matrix containing the d components of dimension m. E is the error matrix which represents the difference between the experimental and computed data and must be minimized by the algorithm. Among many other applications, one of the most widespread technique in Chemometrics, PCA, has been used for more than 20 years for the processing of X-ray absorption data (XAS) 10 . For PCA, the components are considered as orthogonal and do not necessarily have a direct physical meaning. The optimal number of components to use for the decomposition (d) can be estimated from the norm of E(d) and by visual inspection of the components. For higher rank components, the signal becomes essentially noise and does not carry relevant information. This procedure is indeed a way to reduce the statistical noise from a dataset by reconstructing it using a Please do not adjust margins Please do not adjust margins limited number of components. The selected number of components can be then used to launch a MCR-ALS minimization [11][12][13][14] . Like PCA, the MCR-ALS algorithm maximizes the explained variance in the data, but physical or chemical constraints can be applied to C and/or S matrices, rather than the orthogonality constraint used in PCA. In the ALS process, C is least-squares refined with fixed D and S, and then S is refined with fixed D and the previously refined C. The refinements are alternately carried out until convergence is reached. The main interest of the process lies in the possibility to input constraints such as e.g. positivity to S and/or C, which at the end will lead to physically/chemically sensible results

15

. A set of components defining a starting S matrix must be provided to launch the process, for which purpose, several kinds of algorithms have been devised 12 . One interesting feature of MCR-ALS is the possibility to introduce predefined components in the starting S matrix, for example diffraction patterns from isolated phases measured in the same condition. In this so-called "supervised" case, the components are not refined and only their weights are computed. However, although they have been largely used for various types of spectral analyses, these techniques have been only scarcely applied to diffraction data. Burley . Opposite to PCA or MCR, the number of components must be supplied a priori to NMF but no starting values for the components are required. Therefore, no supervised process is possible. Due to the non-negativity constraints on the S and C matrices inherent to the algorithm, the components contained in S are expected to bear a physical meaning, i.e. they will consist in the scattering signal from one or a few phases from the mixture. We report herein the application of multivariate analysis to XRD-CT data obtained from microsamples taken from Cultural Heritage artefacts. Such samples can be composed of many different phases of different natures, from mineral to organic, and with various crystallinity states, from fully amorphous to large crystal grains. They can be intimately mixed at the micron scale and the knowledge of their respective location, layer's thicknesses, etc., is important for the identification of the recipe used by the ancient craftsmen. Therefore, Cultural Heritage samples constitute an excellent playground for developing and testing these techniques.

Here, this experimental strategy has been applied to investigate the local specificities of manufacturing so-called applied brocade decors in the Savoy duchy. Applied brocade is a tin relief medieval polychromy technique developed during the first half of the 15th century in the Low Countries which then spread into all western Europe, which allowed to mimic the rendering of the sumptuous medieval clothes embroidered with gold or/and silver yarns. Its basic component is generally a sheet of tin foil pressed into a mold previously incised with the decorative pattern. The tin sheet is then filled with some pasty material, unmolded and applied on the surface to be decorated. It is then gilded with a gold/silver foil and sometimes highlighted with a colored glaze. These complex artefacts have been submitted to several kinds of degradations over time, such as extensive over-paintings, oxidation of metallic layers, etc. As part of an ongoing interdisciplinary study, our work aims at a detailed description of applied brocade in 15th and 16th century Savoy Duchy 22-25 .

Experimental

Sample

For this study, we have focused our attention on a wooden Pièta statue coming from the Montrottier castle located in Lovagny, France (Figure S1), onto which applied brocade decorations had been identified by visual inspection. In the course of its restoration at the ARC-Nucleart laboratory (CEA, Grenoble), the statue was submitted to detailed in situ inspection using a portable XRD/XRF instrument, confirming at specific places the presence of tin-containing phases and gold, characteristic of the presence of applied brocades 25 . In order to gain insight into the composition, stratigraphy, manufacturing and degradation of these artefacts, several microsamples were collected for further investigations using synchrotron XRD/XRF-CT, at various locations were the presence of brocades had been confirmed using the portable x-ray instrument. These microsamples are chips of millimetric size in their larger dimension and a few hundred microns thick. They were sampled so as to preserve the entire stratigraphy from the preparation layers to the brocade system itself, including possible overpaintings. The selected microsample locations on the statue are shown in Figure S1, and a typical one is shown in Figure 1. Please do not adjust margins Please do not adjust margins using a Zeiss Ultra+ microscope. Microanalysis (EDX) measurements were performed using a Silicon Drift (SDD) detector. In order to identify the organic materials, Fourier Transform Infrared (FTIR) measurements were also performed.

XRD/XRF-CT measurements

XRD/XRF-CT measurements were performed at the BM02-D2AM French CRG beamline 26 of the European Synchrotron Radiation Facility (ESRF, Grenoble, France). The microsamples were mounted on a Mitegen micromount covered with a thin layer of wax avoiding any invasive intervention and allowing the possibility to perform further complementary analyses, fixed on a goniometer head and aligned at the center of the D2AM 6-circle diffractometer. For XRD measurements, a 2D ImXPAD D5 pixel detector was attached to the diffractometer 2θ arm at a sample-detector distance of 215 mm and a 2θ angle of 3.5°. The geometry and beam wavelength were calibrated from the XRD measurement of a NIST-660c LaB 6 capillary sample using the PyFAI library 27 . The XRF detector (KETEK SDD AXAS-M1) was located at 100 mm from the sample in reflection geometry at an angle of 40° from the direct beam. The energy scale was calibrated using the PyMCA software 28 from the positions of a few well identified fluorescence peaks. The energy of the monochromatic beam was tuned at 20 keV (or 0.61992 Å) with a flux of about 10 10 photons/s. This energy has been chosen in order to allow good detection efficacy for a large number of elements using XRF, and to access a wide range of reciprocal space for XRD. The beam was focused by a set of Kirkpatrick-Baez mirrors down to a probe size of 25 x 40 µm 2 offering a gauge volume containing a sufficiently large number of crystallites for powder diffraction averaging and allowing to perform in a reasonable time frame a thorough study of the microsamples. With these conditions, each measurement point required a 1s exposure time. In order to improve the resolution of the image, and since we used the smallest available beam size on the instrument, we applied the oversampling method 3 and decreased the translation step size down from 40 µm to 25 µm, leading to a final voxel size of 25 x 25 x 25 µm 3 in the tomographic reconstructions. In order to reduce the time spent for sample rotation, a continuous rotation scan was used with a 2° rotation step. Thus, the tomographic data collection of a microsample slice of 1 mm maximum length, using e.g. 45 translation points each 25 µm (a few points are added on each side out of the sample), required approximately 2.5 hours. Five microsamples from the Montrottier Pièta were investigated, in which several slices were measured to obtain a more precise 3D description, their locations being chosen by visual inspection. For example, the positions of the slices for the microsample S2018-157 are shown in Figure 1. In the following, we will concentrate on the slice #2 of sample 2018-157.

Pre-treatment of the XRD/XRF-CT data

After collection, the data were treated to account for various effects and organized as data cubes (i.e. 3D matrices) for their subsequent investigation using multivariate analysis. All treatments were made using scripts written in Python 3, and the initial and treated data were stored in HDF5 format. For the XRF data, the recorded spectra (4096 points each) were stored as data cubes without modification. The XRD data required first the transformation from 2D to 1D diffraction patterns which was carried out by radial integration using PyFAI. The 1D patterns contain 3000 points extending from 2θ =0° to 33 ° with a step size of ≈0.011°. During this process, it was realized that many diffraction images showed the presence of strong Bragg spots coming from relatively large isolated grains. Although these grains (identified as mostly gypsum grains) represent a minor volume of the sample, the strong Bragg peaks they produce can greatly reduce the quality of the powder patterns and tomographic reconstructions. Therefore, their contribution was filtered out using the function 'separate' from the PyFAI library. The 1D XRD patterns were then stored as (T, R, D) data cubes. For our microsample S2018-157 which required a 2.05 mm translation scan, a data cube contained 82 x 180 = 14760 pattern of 3000 points each, i.e. 44 280 000 points. The XRD data cubes were submitted to further pre-treatments including:

-limiting the 2θ range; -automatic background subtraction. In such tomographic experiment, the background level and shape can fluctuate quite strongly during the rotation of complex shaped samples, due to variations in x-ray absorption, shadowing of diffuse xray scattering, etc. This can greatly complicate the task of identifying a limited number of components for the multivariate analysis algorithms. Therefore, it is mandatory to remove background beforehand. After several attempts, we adapted a 1D rolling ball algorithm applied sequentially to all patterns, after Savitzky-Golay smoothing. This proved to be efficient and robust, especially to take account of abrupt background changes between adjacent points in the sinogram;

-zeroing the data outside the sample area, which improves the quality of the tomographic reconstruction;

-applying a center-of-mass correction to account for off centering of the sample rotation axis with respect to the x-ray beam, leading to a dissymmetry of the sinogram and artefacts in the reconstructed image. The sinograms were then transformed to reconstructed images and displayed using the inverse radon transform iradon implemented in the scikit-image python library 29 (Figure 2). At this stage, it is possible to examine the sinograms or corresponding reconstructions to look for the location of specific features as phases (for XRD) or chemical elements (XRF).

Please do not adjust margins Please do not adjust margins

Element / phase identification and Rietveld refinements of summed data

The XRF spectrum summed over all voxels is shown in Figure S2 and was used to identify the chemical elements present in the sample. The spectrum is strongly dominated by heavy elements such as lead and mercury. Gold seems to be present as a shoulder on the left side of the mercury peaks. Iron and calcium are the most prominent lighter elements, with minor ones such as copper, nickel and chromium probably coming from contamination by the experimental setup rather than from the sample. Finally, sulfur, tin and traces of chlorine are detectable by their Lα lines at 3.44 keV (Sn) and Kα lines at 2.31 keV (S) and 2.62 keV (Cl), respectively. The peak at 2.98 keV may be attributed to ArKα from air scattering, though the Lα lines of silver may contribute at the same position. These data being the sum of a layered sample rotated in the beam during experiment, a tentative XRF quantification of the elements would be rather inaccurate. Moreover, the elements detected can participate to several phases located at different places and with different purposes. In order to get further insight into the sample organization, it is necessary to make use of the diffraction data.

The XRD data summed over all voxels of the sinogram were used for a first identification of the crystalline phases with the EVA software [START_REF] Bruker-Axs | software DIFFRAC[END_REF] and the PDF4+ database 31 (Figure S3). For complex samples such as those investigated here, which are made of several metallic and painted layers deposited on a substrate, the number of phases can be quite large. It may be difficult to identify all the phases, since Bragg peaks from minor constituents can be hidden among the strong peaks from major ones (see Figure S3). Here, seven phases could be identified: cinnabar (HgS), gypsum (CaSO 4 ,2H 2 O), the two tin oxides cassiterite (SnO 2 ) and romarchite (SnO), two lead carbonates: cerussite (PbCO 3 ) and hydrocerussite (Pb 3 (CO 2 ) 2 ,(OH) 2 ) and beeswax identified as n-paraffin. The presence of gold, which is expected from visual inspection of the sample and from the XRF data, can be only assessed by a small shoulder on the right side of a peak from cassiterite at 2θ = 15.1°. To ascertain the phase identification, a Rietveld refinement was carried out using the FullProf software 32 allowing to better identify the minor phases (Figure 3). Thanks to quantitative analysis, it also allows obtaining an estimate of the mass fraction of the crystalline phases in the sample. Finally, diffractograms from the individual phases can be extracted for subsequent use in the supervised MCR-ALS procedure. For these refinements, a Pseudo-Voigt profile function was used and the background was determined by linear interpolation between selected points. Only the scale factor, an overall atomic displacement parameter, the cell and profile parameters were refined for each phase, the atomic structure parameters being kept fixed. To compare with the results from multivariate analysis (see below) three additional phases were also included in the refinement, although they were not clearly identified in the search-match analysis: minium (Pb 3 O 4 ), goethite (FeOOH) and chlorargyrite (AgCl) resulting in a total of 11 phases. The latter two phases are also compatible with the observation of iron, silver and chlorine in the XRF data. In decreasing order of weight fraction: beeswax (41%), gypsum (23%), cassiterite (15%), romarchite (8%), cinnabar (8%), goethite (2%), hydrocerussite, cerussite and chlorargyrite (1%), minium and gold (<1%). In the case of beeswax, for which no crystal structure has been published, we recorded a XRPD pattern on a Bruker D8 Endeavor diffractometer, with CuKα radiation in Bragg-Brentano geometry, for a natural beeswax sample (purchased from C.T.S., Italy) mixed with corundum powder. This allowed us to check that the XRPD pattern of natural beeswax is very close to those observed in our samples as well as to that of n-paraffin in the PDF4+ database 33 . The corundum powder was used to calibrate the beeswax pattern Bragg peak intensities for the determination of its weight fraction by Rietveld analysis of the tomographic data.

3D analysis of the XRD-CT data using PCA/MCR-ALS and NMF

The data cube was used for PCA decomposition using the scikit-learn python library 34 . Figure 4a shows that most of the data variance is explained by a small number of components. The plot of the first two components (figure 4b&c) also demonstrates that this type of decomposition is not suitable for phase identification, with strong negative signal in component 2, since the orthogonality constraint used in PCA is not adapted to this kind of physical data. Nevertheless, this can be used to estimate the number d of components to be used in a further analysis using e.g. MCR-ALS. MCR-ALS was applied using the pyMCR python package 35 and non-negative least-squares regressor.

Non-negativity constraints were applied to both the C and S matrices (i.e. concentrations and components). This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins To provide a starting estimate of the S matrix, we used the first d components found by PCA and change all negative values to a small positive one to respect the non-negativity constraint. Different numbers of components were tested. In principle, d cannot exceed the number of phases in the sample, i.e. close to 9 for the present sample if we consider the 11 phases used in the Rietveld refinement of the sum data. However, given the large size of the voxels, the presence of components consisting in a mixture of phases expected to be adjacent in the sample (e.g. the two tin oxides) is highly probable, which will reduce the number of components to be found. Figure S4 shows that the data can be successfully described using either 6 or 9 components. The individual patterns for the components can be seen in Figures S5 andS6, as well as their concentration maps for the sinogram and real space reconstructions. For the 9 components case, components 1 and 4 consist almost exclusively in gypsum, while components 0 and 8 contain mainly cinnabar, components 6 and 7 are mixtures of cassiterite and romarchite, etc. As suggested above, the process is not able to decompose the data into individual phases. Figure S7 shows a part of the XRD patterns for components 1 and 4. The patterns are very similar and can be indexed as gypsum, but the peaks are obviously shifted between the two patterns. This might be attributed to an apparent change of position of the diffraction center occurring during sample rotation, due to the very anisotropic sample shape, which modifies the distance to the detector. Since the data can be equally well described with only 6 components a choice of a larger number is not justified. In any case, the location of individual phases in the sample requires an additional step consisting in the identification and quantification (e.g. by Rietveld refinement) of the different phases in each component pattern. However, this task is now much easier than the previous one carried on the sum data, since the number of phases in each single pattern is much smaller and less peak overlapping occurs. This also allows identifying additional phases which could have been overlooked previously. In principle, the MCR-ALS method with these positivity constraints used here should give equivalent results as NMF. In order to check this, we have used the NMF python implementation in the scikit-learn library 29 . Here only the number of components has to be provided. As shown in figure S8 for the case of 6 components, the set of patterns obtained after NMF analysis and their location in the reconstructed slice are similar to those obtained by MCR-ALS, although the order and scaling are different.

This provides an excellent test of the procedure reproducibility. It is also quite interesting to compare the execution times for both methods: using a standard laptop computer under Windows10® (i7 at 2.6 GHz, 16 Go of RAM) the MCR-ALS process took 863 s and the NMF 9.3 s. Weather such a large difference is due to the algorithm or to the programming or implementation efficiency is beyond the scope of this paper, but obviously, the NMF algorithm used here is able to decompose a large data cube in a quite short time, well compatible with the requirements for assessing data quality during a synchrotron experiment. The physico-chemical significance of the obtained component patterns found by NMF with 6 components has now to be checked by using them for phase identification and quantification. The phase matching search results can be seen in figure S9, while the Rietveld plot are shown in Figure 5. In figure S9, the main phases found in the sum pattern analysis can be clearly identified in the component patterns. Moreover, minor phases which cannot be directly observed on the sum pattern can also be detected here. For example, on the pattern of component 5, gold can be clearly seen with 3 peaks at 2θ=15.11°, 17.5° and 24.8° and chlorargyrite with 3 peaks at 2θ=11.1°, 12.8° and 18.2°. A small peak at 2θ=5.7° can be tentatively attributed to minium. In components 0, 2 and 4, two peaks at 2θ=8.5° and 14.5° suggest the presence of goethite. Rietveld refinements of the 6 component patterns were carried out to validate their presence and determine the mass fractions of the phases in each component, which are reported in Table 1. From this table, the phase concentration matrix C p representing the concentration of each phase in each voxel of the sample tomographic reconstruction can be obtained by the C p =C.

t P matrix multiplication, where C is the concentration matrix of the components found by NMF and t P is the transpose of the P matrix containing the phase concentrations of each components. The P matrix, corresponding to Table 1, is of dimensions p.d, p being the number of phases and d the number of components. However, before proceeding to this calculation, the components have to be properly scaled with respect to each other, since no physical scaling is provided by the algorithm.

Please do not adjust margins

Please do not adjust margins In order to scale the components, we derived a method similar to the Reference Intensity Ratio (RIR) method which uses the ratio of the strongest peak of a given phase a to the strongest peak of a standard c (usually corundum) [START_REF] Madsen | Powder Diffraction: Theory and Practice[END_REF] . However, using the strongest peak height here proved to be quite inaccurate because of large variations between the peak widths of different phases, and we rather used the scale factors obtained from the Rietveld refinement instead. For each component, the scaling is calculated using the phase a with the largest weight fraction. The ratio between the scale factors of phase a and corundum c for an equimass mixture can be easily computed as:

where ATZ and V are the molecular mass and the volume of the unit cell as computed by the Fullprof program, respectively. The scale factor to be applied to each component i can then be obtained as:

where

Discussion

Supervised MCR-ALS analysis can provide a good check for the validity of the phase concentrations extracted using (unsupervised) NMF. For this analysis, the patterns of the 11 individual phases introduced in the S matrix were simulated using the results of the previous Rietveld refinement of the sum pattern, with a 0 intensity background. MCR-ALS was applied using the pyMCR python package 35 and non-negative least-squares regressor. Non-negativity constraints were applied to the C matrix, while the S matrix was kept fixed. One can see in Figure 7 that the global agreement on the sum pattern is similar to the one obtained by Rietveld refinement. Moreover, this process directly yields for each voxel of the sample the concentration of the phases which have been introduced in the calculation. A small remaining background is This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins still present at low angles which cannot be accounted for by the component set for which the background is fixed at 0. The observed background is the sum of the remaining signal after background subtraction of the 7380 individual patterns, and will be hard to completely eliminate. Since it varies in an erratic way from pattern to pattern, it cannot be accounted for by introducing an additional background component in the calculation. Also, the peak profile from the Rietveld refinement of the sum data may have difficulty to represent small variations of profile or scattering center during the tomographic data collection. Figure S10 shows the distribution of phase concentrations to be compared with those of Figure 6 found by NMF. The distributions of the phases with the strongest concentrations (gypsum, cinnabar, beeswax, romarchite, cassiterite, cerussite and hydrocerussite) are quasi identical for the two methods. Gold, goethite and chlorargyrite are only partly similar, and the distribution of minium for the supervised analysis is somewhat erratic. Given the small concentrations of these phases and the only approximate agreement of the decomposition using the supervised method as seen from Figure 7, it is no surprise that these phases could not be modelized accurately. On the other hand, the very good agreement for the major phases determined independently by the two methods constitute a satisfactory validation test.

Table 2 shows the mass fractions of the phases as determined by Rietveld refinement of the sum data and by multivariate analysis using NMF decomposition (6 components) and supervised MCR-ALS. The overall agreement is quite satisfactory given the large differences between the methods. The NMF values show a somewhat poorer agreement with the Rietveld ones, especially due to an overestimation of beeswax, a phase that mainly appears as 2 Bragg peaks in the diffraction patterns and may be difficult to estimate accurately. As phase quantification through NMF requires an additional step (Rietveld fitting of the components) as compared to supervised MCR-ALS, the obtained precision maybe decreased.

Although it can look appealing, using the supervised method seems to be prone to several drawbacks. It first requires to correctly identify all the phases present in the sample, through search-match analysis of the sum pattern and to prepare the individual patterns for all these phases. This necessitates to perform the Rietveld refinement of the sum pattern, since the line profiles depend on the experimental resolution function and on the crystalline quality of each phase. This work can be rather long and the description of minor phase profiles cannot be very accurate. Furthermore, as we have seen above, the diffraction center of a given phase can vary during the tomographic experiment, depending on the phase distribution in the sample. The peak shifts due to these variations will yield a broadening and distortion of the peak profiles of the sum pattern. Finally, the presence of a remaining background in the sum pattern may lead to inaccuracies in the localization of minor phases. This method should thus be restricted to relatively simple and small samples where such drawbacks would be limited. The phase distribution in the sample chip can now be examined on the basis of the NMF multivariate analysis results. From bottom to top of the chip (the bottom is on the right side in Figure 8a) the first feature consists in a preparation layer made of gypsum, followed by a layer of cinnabar mixed with some gypsum and goethite. One does not observe a strong signal in the center of the sample, which may be due to the presence of a non-diffracting/amorphous phase, mixed locally with cinnabar, lead white and goethite. More to the top, one observes a strong layer of beeswax, mixed with a somewhat disturbed distribution of tin oxides. Gold mixed with silver chloride and maybe minium are located as two strong spots on top of the tin oxide layer. Finally, on top of the gold spot, a mixture of cinnabar and lead white is observed. This stratigraphy can be interpreted in the following way: a preparation layer made of gypsum was covered by a vermillion paint, on which the brocade system was applied. It consisted in two filling layers, one made of an amorphous phase mixed with various pigments, and a second one made of beeswax on which the tin foil and gilding were attached. The tin layer is now fully degraded into a mixture of the two tin oxides. The gilding was made of gold and silver, possibly zwischgold, consisting in two thin sheets of gold and silver hammered together 37 . Finally, the brocade is covered by an overpainting layer containing lead white and cinnabar. Despite their limited spatial resolution, these observations (such as the nature of the preparation layer, the presence and composition of two brocade filling layers and the use of zwischgold) already constitute important markers which can be used to identify fabrication techniques specific to local craftsmen or schools.

In order to compare these results with those which can be obtained by more classic investigation techniques, we sacrificed the sample 2018-157 which was cut at the level of layer 2, for SEM-EDX and FTIR measurements. The results of optical microscopy and SEM are shown in Figure 8, and the results of the analyses are summarized in Table 3. The Please do not adjust margins Please do not adjust margins observations made with these techniques are in very good agreement with the results of the XRD-CT investigation. The organic phases present in the two brocade filling layers could be identified by FTIR as beeswax (as found by XRD-CT) and an organic layer loaded with lead soaps. For x-ray diffraction, lead soap compounds are known to yield strong Bragg peaks only at d spacing larger than 40 Å 38 , leading to too small angles to be detected in our geometrical configuration. This is also seen as a dip in phase concentrations observed on the profile of Figure 8f showing the major phase concentrations obtained by NMF. On this profile taken from a virtual cut at the middle of the sample, the succession of phases and the layers thicknesses correspond well to what is observed on the optical microscopy image. This demonstrates that the NMF method for analysis of XRD-CT data is able to perform quantitative phases reconstruction in very complex systems without the need to destroy the sample. 3. Right FTIR spectra measured at points 3 (e) and 4 (d), compared with a measurement of pure beeswax (in red). In f, a profile taken at the middle of the sample showing the major phase concentrations obtained by NMF (the scale is in pixel, 1 pixel=24 µm). The concentrations are scaled to 1.

Table 3: results from the optical microscopy, SEM and FTIR investigations of sample S2018-157 layer 2.

Conclusion

The results shown above demonstrate the utility of using multivariate analysis methods, and particularly NMF, for the treatment of large data sets from XRF/XRD-CT experiments. Especially concerning XRD, the method allows to obtain in very reasonable computing time and in an automatic process the decomposition of the data cube into a small set of physically meaningful components representing diffraction patterns which can be used for search-match identification of phases and even Rietveld refinements and quantitative analysis. Depending on the sample complexity and experimental resolution (i.e. beam size vs characteristic size of sample inhomogeneity) the extracted components will consist of single or mixed phases. Using the quantitative analysis results for the components, the complete set of phase concentrations can be retrieved. Applied here to XRD-CT, the same method can be used to any kind of XRD data collection leading to large data sets, such as 2D mapping of large areas which is now often performed on artworks such as paintings or frescoes, etc. and can also be applied to XRF data.

17 .

 17 Rodriguez et al. used MCR-ALS for an operando XRD study of battery discharge 18 , and Taris et al. for time-resolved XRD data of cementing reactions 19 . Long et al. used a similar technique: Non-Negative Matrix factorization (NMF) for the identification of crystalline phase mixtures from thin film by XRD 20 . More recently, the same was successfully applied to combinatorial phase diagram investigations 21

Figure 1 :

 1 Figure 1: Reference microsample, Montrottier Pièta, Saint Mary Magdalene's dress, "continuous brocade", S2018-157. Left: front view showing the presence of gold traces. Center: Back view showing the preparation layer. The nylon Mitegen micromount is visible at the bottom. Right: side views revealing the sample stratigraphyIn addition, SEM-EDX measurements were also performed on some microsamples. The sample was cut at the level of a chosen layer, imbedded in resin and polished for Scanning electron microscopy (SEM) examination. SEM was carried out

Figure 2 :

 2 Figure 2: XRD pattern (top), sinogram (bottom left) and tomographic reconstruction (bottom right) for the sum of all XRD-CT data after data reduction. The approximate shape of the sample is shown as a grey line.

Figure 3 :

 3 Figure 3: Rietveld refinement of the summed XRD-CT data including 11 phases. Rp=12.4, Rwp=16.0, Rexp=9.78, Chi2=2.69.

Figure 4 :

 4 Figure 4: Results of the PCA decomposition: a) cumulative sum of explained variance vs component number; b & c) first and second components.

Figure 5 :

 5 Figure 5: Rietveld refinements of the 6 component patterns extracted by NMF. Component # 0 to 5 from left to right and top to bottom. Chi2=19.5, 6.0, 12.4, 11.3, 11.4 and 14.6, respectively.

  Wa is the mass fraction of phase a in the component i. The intensity values from the pattern of each component i is divided by the corresponding Si value, while the component concentrations from matrix C are multiplied by the same values, so that the calculated C.S matrix containing the reconstructed patterns for all pixels remains unchanged. The correctly scaled C matrix can then be used to compute the Cp matrix containing the phase concentrations in all voxels as shown above. The distribution of the 11 phases in the layer are shown in Figure 6. The composition of any voxel or its evolution along line cuts can be obtained to facilitate the interpretation of the results.

Figure 6 :

 6 Figure 6: Tomographic reconstruction of the concentrations of the 11 phases quantified by Rietveld analysis of the 6 component patterns obtained by NMF. From left to right and top to bottom: gypsum, cinnabar, beeswax, romarchite, cassiterite, cerussite, hydrocerussite, goethite, minium, gold, chlorargyrite. The approximate sample shape is shown as a grey line. The color hue is scaled to maximum concentration for each plot.

Figure 7 :

 7 Figure 7: Comparison of the experimental sum data with the reconstructed sum pattern after supervised MCR-ALS decomposition using 11 components from Rietveld refinement. Blue: experimental pattern, orange: calculated, green: difference (shifted by -10% of the maximum).

Figure 8 :

 8 Figure 8: Observation with optical microscopy (a) and SEM (b & c) of the stratigraphy of sample 2018-157 at the level of layer 2. The number refers to analysis points detailed in Table3. Right FTIR spectra measured at points 3 (e) and 4 (d), compared with a measurement of pure beeswax (in red). In f, a profile taken at the middle of the sample showing the major phase concentrations obtained by NMF (the scale is in pixel, 1 pixel=24 µm). The concentrations are scaled to 1.

Table 1 :

 1 Mass fractions in % of the different phases for each of the 6 components found by NMF after Rietveld refinement.

Table 2 :

 2 Mass fractions in % of the different phases determined by Rietveld refinement of the sum data and by multivariate analysis using NMF decomposition (6 components) and supervised MCR-ALS.

	Phases	Rietveld (sum data)	NMF (6 components)	MCR-ALS supervised
	Gypsum	23.1(1)	18.0	25.6
	Cinnabar	7.8(3)	1.3	8.3
	Beeswax	40.9(4)	57.5	40.0
	Goethite	1.9(1)	1.8	3.4
	Cerussite	0.9(7)	0.5	1.1
	Minium	<0.01	0.1	0.3
	Hydrocerussite 1.3(1)	1.6	1.5
	Chlorargyrite 1.0(1)	0.5	0.8
	Gold	0.2(2)	0.2	0.3
	Romarchite	8.3(4)	7.6	8.8
	Cassiterite	14.6(7)	10.9	9.9

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 7

Acknowledgements

The authors wish to thank the team at ARC-Nucléart (Florence Lelong and Fédéric Fabre) for the choice and preparation of the micro samples. They are grateful to Sophie Omère, Heritage Curator at the Direction Régionale des Affaires Culturelles Auvergne-Rhône Alpes and Julien Coppier, member of the Académie Florimontane (Annecy) and owner of the castle of Montrottier (Lovagny, Haute-Savoie, France), for authorizing access to the sculpture. We thank Pierre-Olivier Autran for his help in the center of mass corrections. This work was financially supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02, Cross Disciplinary Program Patrimalp). The microfocusing set-up was funded by the French National Research Agency (ANR) under "Investissements d'avenir" program with the grant number: ANR-11-EQPX-0010.

Author Contributions

The manuscript was written through contributions of all authors.

Conflicts of interest

There are no conflicts to declare.


Antonio Alguacil
email: antonio.alguacil-cabrerizo@isae-supaero.fr

Michael Bauerheim

Marc C Jacob

Stéphane Moreau

A generic deep learning framework for propagation and scattering of acoustic waves in quiescent flows

A deep learning surrogate for the direct numerical prediction of two-dimensional acoustic waves propagation and scattering with obstacles is developed through an auto-regressive spatiotemporal convolutional neural network. A single database of high-fidelity lattice Boltzmann temporal simulations is employed in the training of the network, achieving accurate predictions for long simulation times for a variety of test cases, representative of bounded and unbounded configurations. The capacity of the network to extrapolate outside the manifold of examples seen during the training phase is demonstrated by the obtaining of accurate acoustic predictions for relevant applications, such as the scattering of acoustic waves on an airfoil trailing edge, an engine nacelle or in-duct propagation. The method is tested for two types of input normalizations, coupled with an a-posteriori correction which improves the acoustic energy conservation of the predictions. The use of an adaptive local normalization along with the physics-based energy conservation results in an error reduction for all the studied cases.

I. Introduction

T propagation of acoustic waves in complex media constitutes a challenging task in the context of aerodynamically generated noise. The presence of complex boundaries, mean flow or background medium inhomogeneities affects the trajectory of propagating waves, through known physical mechanisms such as scattering, dispersion, reflection or absorption. Taking into account such phenomena is required to design quiet aeronautical devices, however creating such pre-design tools with a low cost is still an open challenge. To do so, two types of approaches exist for predicting noise propagation [START_REF] Lele | A second golden age of aeroacoustics?[END_REF]: on the one hand, direct numerical computational aeroacoustics (CAA) calculate simultaneously both acoustic sources and the subsequent wave propagation. Direct computations yield highly accurate results but are highly demanding in terms of computational costs. On the other hand, hybrid methods constitute a more affordable way to perform aeroacoustic predictions. These methods rely on the separate computations of the source region and the propagation into the far field. The propagation step can be achieved either through semi-analytical means (acoustic analogies) or fully numerical approaches (Linearized Euler Equations, Acoustic Perturbation Equations).

The coupling between the source and the propagation regions can become cumbersome, especially when considering the complex media effects on propagation. For example, acoustic analogies for non-homogeneous media, e.g. Lilley's analogy for sheared mean flows [START_REF] Lilley | On the noise from jets. Noise Mechanism[END_REF]), require an additional effort in the computation of the associated Green's function, because of the presence of hydrodynamic modes and critical layers. The same observation applies for complex boundary conditions as in the Ffowcs Williams and Hawkings analogy [START_REF] Williams | Sound generation by turbulence and surfaces in arbitrary motion[END_REF], where new integrals appear when accounting for the presence of moving surfaces. Thus, the more mean flow or boundaries effects are taken into account by the acoustic analogy, the more complex becomes the calculation of the associated Green's function or its subsequent convolution with sources. In the most complex cases, relevant for industrial application, the Green's function can no longer be fully analytically calculated (e.g. annular duct with swirling mean flow analogy [START_REF] Posson | The acoustic analogy in an annular duct with swirling mean flow[END_REF]), and numerical methods must be employed instead. Therefore, in the general cases of complex mean flows or boundary conditions, propagation methods must rely on expensive numerical methods, using for instance high-order numerical schemes to ensure low dissipation and dispersion properties, which prevents these methods from being applied efficiently at a pre-design stage.

A recent work by the authors [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF][START_REF] Alguacil | Predicting the Propagation of Acoustic Waves using Deep Convolutional Neural Networks[END_REF] developed an alternative method to replicate, at lower cost, a high-fidelity numerical code for the propagation of acoustic waves using a surrogate model: instead of resolving the discretized propagation equations, a neural network was trained with examples of high-fidelity direct Lattice-Boltzmann (LBM) simulations, resulting in a data-driven convolutional neural network for acoustic wave propagation. It was found that such an approach managed to reproduce wave propagation in simple academic cases of 2D quiescent closed domains with reflecting walls. These promising results suggest that a similar approach might be applied in open-domain problems, in the presence of complex boundary conditions (such as ducts or scattering objects). In particular, the method has demonstrated good generalization capabilities (e.g. being able to predict wave propagation of unseen source distribution) while diminishing classical numerical constraints (e.g. the CFL number) allowing faster simulations without loss of accuracy.

Yet, the previous framework [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF] implicitly encoded the boundary conditions of the problem, and was thus not applicable to more complex geometries and representative configurations.

Several studies have dealt with the problem of encoding boundary condition information into neural network surrogates for modelling partial differential equations (PDEs). Some make use of built-in tools of convolutional networks, such as periodic padding to encode simple boundary conditions such as periodic conditions [START_REF] Geneva | Modeling the Dynamics of PDE Systems with Physics-Constrained Deep Auto-Regressive Networks[END_REF]. Other works dealing with more complex geometries employ binary inputs [START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF] or signed distance functions [START_REF] Rüttgers | Prediction of Acoustic Fields Using a Lattice-Boltzmann Method and Deep Learning[END_REF] to encode the geometry information. However, such works only tackle steady-state problems, such as RANS predictions [START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF], or Sound Pressure Level (SPL) estimation [START_REF] Rüttgers | Prediction of Acoustic Fields Using a Lattice-Boltzmann Method and Deep Learning[END_REF].

Closer to the problem of interest, works such as Lino et al. [START_REF] Lino | Simulating Surface Wave Dynamics with Convolutional Networks[END_REF] predict the full spatio-temporal propagation and scattering of seismic waves with a convolutional network. They employ a binary field to encode geometry information and use an auto-regressive strategy to perform the temporal prediction. However, few quantitative results are shown, in particular regarding the ability of the network to perform accurately outside the training dataset. Furthermore, it seems that they train the network on different datasets, depending on the different types boundary conditions. A general database capable of training the neural network for different types of boundaries conditions would be beneficial to such kind of methods, as it would avoid the cost of retraining the network each time. A typical example of such a pathology can be found in Physics-Informed Neural Networks (PINNs) [START_REF] Raissi | Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations[END_REF], which need to be retrained for each new set of initial or boundary conditions. Some works, such as Moseley et al. [START_REF] Moseley | Solving the wave equation with physics-informed deep learning[END_REF] employ PINNs to solve the wave equation for complex background media [START_REF] Moseley | Solving the wave equation with physics-informed deep learning[END_REF], and parameterize the sources input position in order to avoid retraining the model for each new set of initial conditions. However, the problem of using a single neural network for various types of boundary conditions, such as reflective and non-reflective acoustic conditions, is still an open question.

Consequently, the present work proposes a unified generic framework based on an deep auto-regressive convolutional neural network to tackle both bounded and unbounded acoustic problems in quiescent flows, with potential scattering by inner objects. The work focuses on two mains aspects: (i) the ability of the trained neural network to perform in various physics regimes dependent on the choice of boundary conditions, some of which differ significantly from the ones present in the chosen database, (ii) the practical choices for pre-and post-processing the data when predicting in the auto-regressive loop. In particular, two types of data normalization are compared, along with a modified Energy Preserving Correction, which was introduced in [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves using Deep Convolutional Neural Networks[END_REF]. The latter intends to improve the conservation of acoustic energy over the temporal neural network predictions, and thus improving significantly its long-term predictions.

First the methodology regarding the choice of neural network architecture and its training are presented in Section II, Then, section III describes the unique dataset employed for training the neural network and its validation. In section IV, several test cases representative of unbounded problems are evaluated using the present framework, and compared with the LBM reference. In particular, Sections IV.C and IV.D explore two additional cases representative of typical industrial applications, namely a NACA airfoil and an engine nacelle. Finally, duct configurations, as well as a completely closed domain, are also studied in Section V.

II. Methodology

A. Modeling dynamical systems with a learned surrogate

The objective of this work is to efficiently approximate the complete space-time evolution of a dynamical system with a surrogate, model trained on high-fidelity data. Here, the system of interest is the propagation and scattering of acoustic waves with hard-reflecting obstacles in quiescent flows, which can be described by the following Cauchy problem for the acoustic density wave equation:

2 ( , ) 2 + 2 0 Δ ( , ) = ( , ), ∈ D, ∈ [0, ] B ( ( , )) = 0, ∈ D ( , 0) = I 0 ( ) (1)
where is defined in the domain D with boundaries D. B is the boundary operator that enforces the boundary conditions (here, B = 0 expresses the absence of sources along the boundaries). Typically for hard-reflecting walls it simply reads ∇ ( • )| D = 0. 0 is the speed of sound, a source term and I 0 is the initial condition expressed as a volume distribution. For an arbitrary boundary condition operator B, Eq. ( 1) is usually resolved numerically, by discretizing the state variable in space and time and subsequently integrating the time evolution from to + Δ using an iterative solver, where Δ is the time-step between to iterations.

The goal of the present study is to recast this numerical integration as an optimization problem, in order to learn the space-time evolution of the acoustic density. Discretizing the state variable in space and time, let ∈ R 1 × 2 be the solution of Eq. ( 1) at time-step on a uniformly spaced discretization of D defined by a grid with nodes in each direction.

The learned PDE operator , defined by its trainable parameters , must fulfill the following time-invariant equation,

+1 = ( +1 , ) (2)
where +1 = { , -1 , ..., -} is the input of the surrogate model, corresponding to the previous time-states. Note that following [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF], the present study will use 4 input snapshots ( = 3). The operator is obtained by solving the following optimization problem:

arg min L ( ( +1 , ), +1 ) (3)
where L is a metric of the distance between the approximated prediction from the learned operator ˆ +1 = ( +1 , ) and the actual solution +1 . Here +1 is obtained through a high-fidelity numerical simulation database as described in section III, but could be obtained by other means: analytical solutions (if existing), experimental observation or a blending of the previous. Such a space-time operator can be learned through a highly non-linear neural network regressor. For high-dimensional databases, such as the ones encountered in Computational Aeroacoustics or Computational Fluid Dynamics, the use of convolutional networks is an efficient solution to simultaneously learn the spatial and time integration, as shown in earlier works [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves using Deep Convolutional Neural Networks[END_REF][START_REF] Geneva | Modeling the Dynamics of PDE Systems with Physics-Constrained Deep Auto-Regressive Networks[END_REF][START_REF] Lino | Simulating Surface Wave Dynamics with Convolutional Networks[END_REF]. The main advantage with respect to traditional numerical schemes is two-fold: in comparison with explicit schemes, such as Euler or Runge-Kutta integrators, the learned surrogate is not constrained by the classical time-stepping stability constraints [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF][START_REF] Liu | Hierarchical Deep Learning of Multiscale Differential Equation Time-Steppers[END_REF]. Furthermore, it does not require any matrix inversion as implicit integrators do.

B. Autoregressive spatio-temporal prediction of wave propagation

Once trained, the resulting neural network is employed in an auto-regressive strategy to predict the complete spatio-temporal evolution of acoustic fields. Although the network is only trained to predict one single time-step ahead from the input data, this auto-regressive strategy is used as a time-integrator to propagate signals up to an arbitrary time horizon.

Formally, the acoustic density at time = Δ can be calculated as

= • • ... • --1 times ( +1 ) (4)
where +1 = { , , ..., 0 } contains the initial condition. In practice, = 3 is chosen, so that the last predicted frame ( + 1 in Fig. 1) is added to a new series of inputs (composed of -2, -1, and + 1 in this example).

C. Neural network: U-Net convolutional network

A U-Net convolutional neural network [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] is employed as shown in Fig. 1. Such an architecture is capable of treating large structured image-like inputs with few trainable parameters by sliding small local trainable convolutional filters around the input (typically 3 × 3 or 5 × 5 filters). In comparison with classical fully-connected networks, convolutional networks naturally encode neighbourhood information between pixels which is key when resolving Partial Differential Equations. This can be related to the stencils used to calculate the spatial gradients in a discrete grid.

Furthermore, convolutions can be combined with down-sampling operations to increase the receptive field [START_REF] Luo | Understanding the Effective Receptive Field in Deep Convolutional Neural Networks[END_REF], in order to capture long-range spatial information. The U-net architecture combines convolutions with pairs of downand up-sampling operations to treat the input in a multi-scale way. In this case, max pooling operation along with bi-linear up-sampling is employed. The multi-scale treatment separately processes the different scales of the problem, specialising the different filters to focus either on local or global flow features. Unlike for the Multi-Scale used in previous works [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF], it was found that the U-Net network converged quicker than the Multi-Scale one. This can be attributed to the presence of skip connections that facilitate the backward flow of gradients during the optimization of the network. Moreover, U-Net have better inference time performance compared with the Multi-Scale net of [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF], which further accelerate the acoustic predictions.

The input to the neural network corresponds to the previously defined +1 vector, concatenated with an additional obstacle grid indicating the presence of the reflecting obstacles with a Boolean encoding (ones assigned to obstacles pixels and zeros to fluid ones). Thus, the position of the boundary conditions are explicitly given as inputs to the neural network, as employed in other works using convolutional neural networks in fluid-mechanics related problems [START_REF] Thuerey | Deep Learning Methods for Reynolds-Averaged Navier-Stokes Simulations of Airfoil Flows[END_REF][START_REF] Ajuria Illarramendi | Towards an hybrid computational strategy based on Deep Learning for incompressible flows[END_REF].

The time-step Δ at which the prediction → + Δ is performed is fixed by the temporal spacing between input and output, thus in order to change its value, a new training is required.

Appendix A details the parametric functions used for each one of the U-Net scales, as well as the number of trainable parameters.

III. Dataset D1: Propagation of Gaussian pulses in unbounded domain with hard wall obstacles

As in [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves using Deep Convolutional Neural Networks[END_REF], a LBM code is used to generate the data for the optimization of the network on the training dataset. Such a dataset, employed to feed inputs into the neural network and compare its outputs with supervised references, is decribed in the following section. A particular attention is paid to the validation of the LBM code on the wave scattering by objects to guarantee the accuracy of the training database. The dataset (indicated as D1) consists of 100 2D-simulations of time-propagating acoustic waves on square domains with non-reflecting boundary conditions and random obstacles with hard reflecting walls (cylinders and rectangles) placed inside the calculation domain, as seen in Fig. 2. This dataset aims at studying the capacity of the neural network to predict scattering of acoustic waves by such hard-wall reflecting obstacles. The physical dimension of the computational domain is chosen as (see Fig. 3a), each direction being discretized with = 200. Similar to the dataset used in [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves using Deep Convolutional Neural Networks[END_REF], the wave sources are composed of a random number of Gaussian pulses which are used as initial conditions for the density field, namely:

( , 0) = 0 + =1 exp - log 2 2 || -|| (5)
where 0 is the mean flow density, is sampled from a uniform distribution U [START_REF] Lele | A second golden age of aeroacoustics?[END_REF][START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF], is the pulse amplitude, here fixed to = 10 -3 , is the pulse half-width, fixed to /Δ = 12 and is the pulse center, such that = U [0.2 , 0.8 ] 2 . Both the number of cylinders and rectangles are drawn from a random uniform discrete distribution U [0, 3]. The rectangle sizes (height and width) vary from 0.01 to 0.3 , while cylinder radii vary from 0.08 to 0.15 . Boundary conditions at the obstacle walls are modelled as bounce-back conditions. This type of boundary is not generally fitted for inclined boundaries, such as the cylinder wall. However, in absence of mean flow, if the wavelength of acoustic waves is sufficiently large with respect to the cell size Δ , then the bounce-back "staircase" will be compact and have no effect on the accuracy of numerical results. This is demonstrated in section III.B, where the LBM method is validated for wave scattering by a cylinder. Non-reflecting boundary conditions are modelled using a perfectly matched layer (PML) that prescribes a damping term in the governing equations in order to attenuate outgoing acoustic modes [START_REF] Brogi | Hermite regularization of the lattice Boltzmann method for open source computational aeroacoustics[END_REF]. The output fields from the LBM simulations are cropped to only feed the physical domain to the neural network.

Results are reported using a normalized time, defined as

= 0 / (6)
A grid spacing of Δ = /( -1) = 0.5025 is employed. Notice that in LBM the value of the time-step Δ is proportional to the grid spacing, here Δ = / 0 Δ . This acts as a Courant-Friedrichs-Lewy (CFL) condition, depending on the sound speed in lattice units , which is determined by the lattice choice, here equal to ( ) = 1/ √ 3 0.57 in the 2 9 lattice [START_REF] Marié | Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF].

The time corresponds to the propagation time of an acoustic wave from one boundary to the other. Each training simulation is stopped at = 0.5, when the Gaussian pulse has reached the absorbing layer and density fields are recorded at time-steps that are multiples of Δ = 0.0087 (i.e. each 3 LBM iterations). The latter time-step is the one used by the CNN, in order to demonstrate the CNN capability to replicate high-fidelity codes while mitigating their CFL constraint. In practice, the computed LBM fields are packed into groups of 4+1 frames (inputs + target) for the CNN training, as detailed in section II.B.

Note that Convolutional network are translation-invariant by construction [START_REF] Lecun | Backpropagation Applied to Handwritten Zip Code Recognition[END_REF]. However, they are not invariant by rotation. Thus, data augmentation is performed during the training of the network: input-target pairs are randomly flipped in angles of 0, /2, and 3 /2 radians. The LBM code is validated with respect to the scattering of a Gaussian pulse by a cylinder. A schematic is shown in Fig. 3b. A rigid cylinder of radius is located at ( , ) = (0, 0) inside a domain of size 30 × 30 , plus an additional 2 of sponge zones at each boundary. A density Gaussian pulse of half-width / = 0.4 is located at ( , ) = (8 , 0) and a probe is located at ( , ) = (0, 10 ). A recursive and regularized BGK collision model (rrBGK-LBM) [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: recursivity and regularization[END_REF] is employed to maintain code stability with low numerical dissipation. Two cases are compared: a first one with the pulse half-width discretized by 12 lattice points (901 lattice points in the domain excluding sponge zones), and a second one with 18 points per half-width (1351 lattice points in the domain excluding sponge zones).

B. Validation of the dataset

Results are compared to the analytical solution derived by Tam [START_REF] Tam | Computational aeroacoustics: A wave number approach[END_REF]. For a pulse of half-width and pulse parameter = (log 2)/ 2 located at ( , ), the total unknown potential can be divided in two components related respectively to the incident (subscript ) and reflected waves (subscript ):

( , , ) = ( , , ) + ( , , ) (7)
The total pressure field is found as follows:

( , , ) = - = -( + ) = ∫ ∞ 0 ( ( , , ) + ( , , ) - (8)
where ( , , ) is given by

( , , ) = 1 2 - 2 (4 ) 0 ( ) (9)
and = ( -

) 2 + ( -) 2 .
represents the Bessel function of first kind of order k. ( , , ) is given, in polar coordinates ( , ):

( , , ) = ∞ =0 ( ) (1) ( ) cos( ) (10)
where

( ) = 2 - 2 (4 (1) ( 0 ) ∫ 0 1 ( 0 ) 0 -cos( ) -sin( ) 0 cos( ) (11)
with 0 = 1 and = 2 for ≠ 0. 0 is defined as:

0 = 2 0 -2 0 cos -2 0 sin + 2 + 2 (12)
and (1) is the derivative of the Hankel function of the first kind of order . Figure 4 shows the time evolution of the density at the probe location for two different resolutions of the LBM lattice. A resolution of 12 point per half-width is sufficient to capture all the physics of the problem accurately. The first amplitude peak represents the direct propagation of the initial Gaussian pulse, while the other two peaks at 0 / = 8 and 10 arise from the pulse scattered by the cylinder wall. Therefore, the pulses in the dataset are resolved with /Δ = 12 points per initial half-width.

C. Data normalization and energy-correction

Neural networks benefit from scaling the inputs so that the mean and the standard deviation of their inputs are close to 0 and 1 respectively [START_REF] Goodfellow | Deep Learning[END_REF], in order to improve the stability of the neural network training. However, scaling all the inputs with their statistical moments generates a loss of physical information for the neural network, most importantly in time-dependent problems, as the network may no longer capture the relative changes of amplitude between different time-steps, accounting for phenomena such as diffusion or in the present case, energy decay (e.g. the acoustic signals leave the computational domain after some time). Thus, two types of input normalizations are proposed in this work, in order to evaluate the performance of both strategies. Furthermore, an a-posteriori correction based on the conservation of acoustic energy is also presented to decrease the error accumulation over time of the learned model.

Global normalization :

The first normalization, used in a U-Net network dubbed as GlobUnet, used a global scaling to normalized all inputs by 1/ where is the initial amplitude of one Gaussian pulse, such that:

˜ = -0 (13)
This scaling intends to preserve the scale of the original dataset, which may be important to account for time-decaying phenomena.

Local normalization :

The second normalization, denominated LocUnet, employs an input-to-input scaling. In particular, each input is normalized by ( -where is the standard deviation of first frame for each input:

˜ = -0 ( -) (14)
Because four frames ( = 3) are used at each input, the network is presented with information about amplitude decay between all the different snapshots, although more locally than for the global normalization. The disadvantage of this method is that some of the network inputs at late simulations time (e.g. when the acoustic signal has mostly left the computational domain) may have a small amplitude compared with the initial time-steps, with greater amplitude. Thus, the network may struggle when processing such kind of low-signal inputs, and focus only on high-energy ones. In order to study such claims, both normalizations are employed for all test cases of Sections IV and V. Energy preserving correction: In [START_REF] Alguacil | Predicting the propagation of acoustic waves using deep convolutional neural networks[END_REF], authors found the the network predictions could be improved by using an a posteriori acoustic energy-preserving correction, tailored to the boundary conditions of that study, namely computational domains closed by reflecting walls. The acoustic intensity fluxes being nil, the acoustic energy is preserved over time in the inviscid case with no mean flow, thus a spatially uniform correction yielded:

= 0 -ˆ (15)
where

ˆ = ∫ Ω ˆ ∫ Ω 0 = ∫ Ω 0 ∫ Ω . (16)
and ˆ ( ) is the prediction of neural network at time . In the case of radiative boundaries (e.g. numerically treated as non-reflecting boundary conditions), the previous assumptions no longer holds, as the acoustic intensity fluxes may no longer be nil, as shown by the acoustic energy equation in conservation form:

∫ Ω + ∫ Ω I • n = 0 ( 17
)
where represents the acoustic energy density and I the acoustic intensity:

= 2 2 0 2 0 + 0 2 2 I = v. (18)
Here, a modified version of the energy correction for open domain cases is discussed. Here / 0 , and thus the correction yields:

= ˆ -ˆ +1 (19)
where

ˆ = ∫ Ω ˆ ∫ Ω . ( 20
)
This correction only remains valid for small time steps or low frequency acoustic sources. In total, four strategies are tested by training the presented U-Net CNN on the dataset 1, and using the two normalization strategies described above as well as the energy preserving correction. The four models are summarised in table 1

D. Training the neural network

The supervised training of the neural network is performed over the dataset composed of pairs of input-target samples ( +1 , +1 )) ≤ ≤ .

Equation ( 3) is solved by defining the loss function L used to train the neural network, which is computed at the output of , and reads:

L = 1 =1 || ˆ +1 , +1 || 2 2 ( 21
)
where the summation is done over the total number of data samples of prediction-reference pairs. The loss minimizes the L2-norm mean-square error (MSE) of the prediction ˆ +1 with respect to the reference field +1 .

In practice, the optimization of the neural network is performed using a variant of the classical mini-batch stochastic gradient descent algorithm, the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], with an initial learning rate of 4.10 -4 and a learning rate schedule which decreases its value by 15% each time the network validation loss reaches a plateau for more than 10 optimizer passes on the full dataset (called epochs). The Pytorch open-sourced framework [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] is employed to perform the training operations. The algorithm is stopped when the loss converges to a plateau for 100 consecutive epochs, which is achieved after a total of 1500 epochs. The total training time in a NVIDIA Tesla V100 GPU takes 48 wall-clock hours, with a batch size of 32.

IV. Density-field prediction in unbounded test cases

The autoregressive methodology is tested in a series of unbounded cases to demonstrate the ability of the framework to predict acoustic propagation for several configurations with a quiescent mean flow. The objective is to demonstrate that a single neural network, trained on a dataset such as the one described in section III, can perform accurately on all the previous test scenarios, without the need to build a tailored database for each case. The interest of such approach is to avoid the repetition of the costly training process for each new configuration, and to demonstrate that the employed neural network is able to perform well on new cases not seen during the training process.

A. Free-field propagation

First, the network is tested on the free-field propagation case, for a domain of size × without any obstacles and with non-reflecting boundary conditions (Fig. 5a). The initial condition corresponds to a Gaussian pulse of half-width / = 0.06, located in the center of the domain. The neural network is provided with the first four frames ( = 3) of the simulation, issued from a LBM simulation with the same boundary and initial conditions, and the auto-regressive method unrolls subsequently the complete prediction, until the non-dimensional time = 2 is reached.

Snapshots at several times are shown in Fig. 6, for the four studied networks and compared with the reference LBM method. Results show that all neural network perform well during the first phase of the propagation, before the pulse reaches the domain outer boundaries. Once the pulse interacts with the non-reflective boundary, some networks (GlobUnet and LocUnet without correction) seem to struggle keeping the symmetry of the pulse while LocUnet with correction manages to keep this symmetry much better, even when interacting with corners.

More quantitative results are shown in Fig. 7, where a slice at constant / = 0.5 is extracted. The energy preserving correction seems to benefit the LocUnet method (green curves) whereas it hampers the GlobUnet one (red curves). This can be attributed to the lack of uniform drift in the GlobUnet case, thus the hypothesis for correction by some constant value no longer holds (see section ??). However, in the LocUnet case, such a correction seems to improve clearly the accuracy of the results. This can been for example in Fig. 8, where the time evolution of two global metrics is shown, namely the evolution of the Mean Squared Error L 2 (relative to the initial pulse amplitude ) and the evolution of the acoustic pressure energy over time. The use of the energy preserving correction for the local-normalized U-Net manages to mitigate the error accumulation over time, typical of auto-regressive methods. However it does the contrary for the globally scaled U-Net, suggesting that no error drift is occurs for such model. In fact, the uncorrected version of GlobUnet behaves very similarly in terms of L2 error evolution than the corrected version of the LocUnet method. Furthermore, it can be observed from the energy plot that for the uncorrected version, both networks show similar trends in the constant energy regime ( < 0.5) corresponding to times before the pulse reaches the domain boundaries. However, the sudden change of regime towards an energy-decaying one after = 0.5 seems to be better captured by the global normalization. The addition of the EPC makes the local normalization the best performing method, both in terms of error and energy fit.

This first test highlights one of the key challenges of such data-driven methods: the change from energy-preserving to energy-decaying state with very low residual levels in the domain and decreasing signal-to-noise levels, even though the physics behavior remains identical (i.e. the propagation of waves in free-field).

B. Scattering by a cylinder

Next, as in the validation case of the dataset in section III.B, the scattering of a Gaussian pulse by a cylinder is studied. The computational domain is sketched in Fig. 5b / = 0.06, located at a distance / = 0.25 from the cylinder of radius / = 0.075, located in the center of the domain, while the four boundaries are left as non-reflecting ones.

Results for such a test-case are shown in Fig. 9. Density fields of both the reference LBM simulation and the neural network predictions are compared, along with the local error normalized by the ∞ norm of the error, in order to show where the error is concentrated. The challenge for the neural network is to correctly predict the scattering of the pulse by the cylinder, and also to correctly model the non-reflecting boundary conditions, as already discussed for the free-field propagation case.

Slices of density fields at two positions / = 0.5 and / = 0.25 are shown in Fig. 10, as well as evolution of two global quantities in Fig. 11, namely the mean-square error between the LBM target and the network predictions, and the temporal evolution of the total fluxes integrated over the four domain outer boundaries over time, namely:

( , ) = ∬ Ω (22)
In order to calculate such fluxes, the acoustic velocity has been inferred through the linearized conservation of the Euler equation, which yields:

( ) = ∫ 0 ∇ ( ) (23)
It can be seen that for all four tested networks, a similar behaviour is observed as in the previous test: the error grows initially next to the cylinder wall, while at later times the regions of high error are concentrated near the boundaries. The LocUnet+C (i.e. normalized by the standard deviation and with energy preserving correction) achieves a lower error than the rest of methods. Here the GlobUnet encounters more difficulties, typically for times after = 0.5. In particular, it can be seen that the network does not manage to correctly predict the absorption by the rightmost domain boundary at / = 1 of the wave reflected by the cylinder. The error is locally higher in that region near the boundary and the EPC is unable to correct such a behavior (e.g. for = 0.76). On the contrary, the energy-corrected LocUnet manages to closely follow the signal oscillations even for such long time horizons. The evolution of the fluxes in time confirms this observation: the GlobUnet strategy only manages to capture the high amplitude fluxes at early simulation times, but fails for the lower amplitude peaks at 0.6 and 0.8. This suggests that the global scaling forces the network to focus on high-amplitude signals, while lower amplitude ones have been neglected during the training process. Thus, it demonstrates the superiority of the local normalization, which can adapt to varying signal amplitudes. In combination with the energy correction, this local normalization is able to capture accurately the energy-decaying dynamics.

C. Scattering with NACA0012 airfoil

The third test case corresponds to the scattering of a Gaussian pulse wave with a NACA0012 airfoil. The computational domain is sketched in Fig. 5c. The airfoil leading edge is located at coordinate ( , ) = (0.3 , 0.5 ) while the chord and thickness are 0.4 and 0.08 respectively. The initial acoustic source is a Gaussian pulse located at ( , ) = (0.5 , 0.8 ), above the airfoil trailing edge. This is the first test case where the network encounters a geometry which was not present in the database. The purpose is to study the ability of the network to predict the strong scattering typical of sharp edges, as in the airfoil trailing edge. Since the corrected GlobUnet and the uncorrected LocUnet have been found ineffective on the two previous cases, only the two relevant networks (GlobUnet and LocUnet+C) are shown in the following.

Results in Figures 12 and13 show that for low times of propagation ( < 0.4) the error concentrates at the trailing edge region: as seen for example in the slices for = 0.28 and = 0.35 both GlobUnet and LocUnet+C do not fit closely the density levels at the trailing edge ( / = 0.7), indicating the difficulty of the network to handle sharp edges, which it was not trained for. The leading edge is on the contrary well predicted. This offset in density creates some dispersion phenomena visible in the error fields of Fig. 12: the back-scattered wave from the trailing edge is not in phase with the reference simulation. Interestingly, the prediction error near the sharp trailing edge has a pattern similar to a dipolar source term. Note that the GlobUnet strategies exhibits artifacts close to the trailing edge, even at early times (e.g. = 0.10). A better understanding of this localized error and possible corrections are left for future work.

D. Scattering by a 2D Nacelle

In order to demonstrate the capabilities of the presented framework in an industrial-relevant case, a simplified duct-nacelle configuration is tested, as shown in Fig. 5d. Each side of the nacelle is modeled as a symmetric NACA0012 A Gaussian pulse source is initially located inside the duct at ( , ) = (0.5 , 0.5 ), and the neural network is unrolled to propagate the wave through the duct and into the far-field.

This case contains most of the physics studied in this work: wave reflection, duct acoustics and scattering of waves at both ends of the duct. Results shown in Fig. 14 reveal the complex patterns that appear in such a case: axial-traveling waves propagate both upstream and downstream of the duct, while transverse waves remains trapped inside the duct, forcing the appearance of duct-like modes. Both GlobUnet and LocUnet+C manage to predict accurately the density patterns appearing in such a scenario. A slice along the axis of the nacelle ( / = 0.5) (Fig. 15) shows the overall fitting of the LBM reference, even after the trapped modes become established after few iterations ( > 0.24).

However, as was seen in section IV.C, the scattering by both trailing edges creates a local error accumulation, which creates a phase shift in the subsequent propagation. Again, early artifacts appear for GlobUnet at the sharp trailing edges (e.g. = 0.1), which then propagate and pollute the rest of the time predictions. This can be seen in the plot of density. After the in-duct modes haven been established ( > 0.5), this phase shift becomes clear.

Furthermore, the evolution of acoustic energy and fluxes in the domain, shown in Fig. 16, evidences the presence of this in-duct modes through the periodic oscillation of energy and fluxes. Both neural networks show signs of phase shifts as they do not follow the oscillatory patterns after > 0.5. However, LocUnetD1+C manages to follow the global trend of oscillating energy, while the GlobUnet diverges faster from the reference simulation.

Thus, this test highlights that the studied method is able to accurately predict complex acoustic propagation cases which differ significantly from the training database. It also shows two of the main difficulties encountered by the surrogate model: the scattering with sharp edges and the difficulty to predict in-duct propagation, because of the presence of trapped modes. The scattering seems to be the cause of the in-duct propagation discrepancy as it generates dispersion errors.

V. Acoustic prediction in bounded domain

A. In-duct propagation

In order to complete the observations made in the previous case regarding the difficulty of neural networks to tackle the prediction of in-duct wave propagation, a fifth test case is presented, where the propagation of a Gaussian pulse inside a two-dimensional duct of constant diameter is considered. It allows the study of the long-time propagation by neural network of trapped modes, without taking into account difficulties due to sharp edges Three configurations are studied as shown in Fig. 17. The size of the domain is kept constant at × while 3 diameter ℎ of the duct are investigated: ℎ/ = 0.4, 0.6 and 0.8. Figure 18 compares the results for both GlobUnet and LocUnet+C which were identified as the best performing methods in the previous tests. The initial bouncing on the duct walls is seen to be well predicted ( = 0.38). However, once the wavefronts reach simultaneously the four edges of the simulation, where the non-reflecting condition coincides with the duct wall, the error accumulates quickly in that region. This behavior is clearly seen in the GlobUnet case, which creates some spurious back-reflected waves that pollute the entire domain. While this phenomenon is less present in the LocUnet+C case, it it still visible. This mismatch creates a dispersion error that can be observed in the energy evolution for the three studied cases (Fig. 19). In the three cases, the LocUnet+C method demonstrates its superiority as it manages to track the energy evolution much closer than the globally-scaled network. However, the phase dispersion is also visible in both cases, which is an effect of the error created by the edges of the simulation. These can be attributed to the fact that such a type of junction between non-reflecting and wall boundary conditions is never seen during training. Thus it is hard for the network to infer the behavior in such points. Here, it is hypothesised that even if the duct walls are of infinite length in the axial direction, the network has always learned on finite-length obstacles. Thus, it may be possible that it sees the duct-end as an open-end, thus creating some reflected waves because of the wave scattering at a duct opening. In order to mitigate this behavior, the dataset could be augmented with some examples of infinite length ducts.

Furthermore, the analysis of the energy change over time in Fig. 19 reveals that the GlobUnet is not robust with respect to an increase of the duct aspect-ratio. For values ℎ/ > 0.6, such a normalization does not seem capable of following the changes of acoustic energy over time. On the other hand, the LocUnet+C method fits much closer the reference data, although the dispersion phenomena is clearly visible in such results.

B. Closed box propagation

The last test case corresponds to the same type of boundary condition studied in [START_REF] Alguacil | Predicting the Propagation of Acoustic Waves using Deep Convolutional Neural Networks[END_REF], namely a hard reflecting wall on the four boundaries of the computational domain. The initial condition corresponds to a Gausian pulse centered at ( , ) = (0.5 , 0.5 ). In practice, the neural network is given a geometry mask with one continuous layer of "obstacle" pixels around the computational domain. Results in Fig. 20 show the evolution of the neural network predictions for both tested networks (GlobUnet and LocUnet+C). In both cases, both networks manage to reproduce the highly symmetric patterns until times around = 0.6. For later times, a slight phase shift, originating again from the corners, unbalances the symmetry in both cases. As the energy content of this problem remains constant, it is a hard problem for networks that have only seen decaying problems in the training database. The hypothesis is made that the network may need further temporal context, i.e. having access to longer histories of previous snapshots at the network input, or training on a database of longer simulations, in order to accurately capture the long temporal dependencies of the problem. The latter is tested by training GlobUnet on a new dataset D2, identical to D1 except that the LBM simulations are stopped at later times ( 2 training = 1.44 instead of 1 training = 0.87 ). Results in Fig. 20 show that this new training is able to improve significantly the results for this particular test case. However, all the previous tests were performed with this longer dataset but no improvement was found for the other cases. The explanantion for such a phenomena remains an open question.

VI. Conclusion

In this work, a method to treat a variety of acoustic propagation cases on quiescent media with a deep learning surrogate model is presented. The neural networks are trained on a single database of high-fidelity LBM simulations which contains examples of acoustic wave propagation, reflection and scattering with obstacles, as well as free-field propagation. Multiple test cases demonstrate that a convolutional neural network trained on such a database, with some specific normalization (namely local normalization by the input standard deviation) manages to reproduce closely the results of the LBM reference. Accurate results are obtained even on boundary conditions and scattering configurations not seen during the learning phase. The combination of a locally-normalized training with an a posteriori correction, which imposes the conservation of acoustic energy, significantly improves the results and reduces the long-term error accumulation of the autoregressive method. The different examples also highlight some of the challenges encountered by the network, that are also typical of traditional numerical solvers: the accurate treatment of non-reflective boundary conditions, the scattering of acoustic waves by sharp edges and the prediction of duct modes in confined configurations. It also demonstrates that these highly efficient non-linear networks may be capable of extrapolating to changes in the underlying statistical data distribution (e.g. energy-decaying versus energy-conserving flows).

While neural network techniques for fluid dynamics and aeroacoustics are in their early phases of development, studying their capability in more complex cases remains a crucial step for their applicability to real world problems. The insights gained in this work suggests that besides the studied problem of treating the temporal integration accuracy [START_REF] Geneva | Modeling the Dynamics of PDE Systems with Physics-Constrained Deep Auto-Regressive Networks[END_REF][START_REF] Azencot | Forecasting sequential data using consistent koopman autoencoders[END_REF], care should be taken when modeling problems with complex boundary conditions. Thus, techniques from the CFD community could be employed to improve the capabilities of such data driven methods (e.g. non-reflective conditions etc). Figure 21 shows the details of operations. The notation is follows the following convention • # represents the size of the 2D convolutional kernel. For instance, 3 is a 3 × 3 kernel, with 9 learnable parameters. • # denotes the stride of the convolution operator.

A. U-Net convolutional network

• # is the number of output channels after each convolution.

• The non-linear activation function: either "ReLU" for Rectifying Linear Unit, or the idendity function ("Id").

Fig. 1

 1 Fig. 1 Schematics of U-Net multi-scale CNN with three scales of convolutions at full, half and quarter resolution. Input is composed of 4 consecutive frames = -3, = -2, ..., = of acoustic density and of a Boolean mask representing the position of the obstacles. The output is the next frame at = + 1 of resolution NxN voxels. Each rectangular block is proportional to the number of convolutional filters per layer. Skip connections are employed to facilitate the training of the network.

Fig. 2

 2 Fig. 2 Example of one simulation from the dataset at times (from left to right) = 0.03, 0.17, 0.35 and 0.52. Initial conditions are composed of Gaussian pulses of acoustic density located at random positions. Each initial condition ( = 0) is stepped for 300 LBM iterations in total ( = 0.87). Acoustic density fields are saved every 3 LBM iterations.

Fig. 3

 3 Fig. 3 Schematic of the dataset configuration. (a) Sketch of boundary and initial conditions, made of Gaussian pulses of half-width propagating in a domain of size × with arbitrary cylindrical and rectangular obstacles located inside the physical domain. A sponge layer of width is added to absorb the outgoing acoustic waves. (b) Configuration used for the validation case: a cylinder of radius located at the domain center scatters a Gaussian pulse. Measurements are performed at the probe location.

Fig. 4

 4 Fig. 4 Fluctuating density at probe location ( , ) = (0, 10 ) for dimensionless time = 0 / .

Fig. 5

 5 Fig. 5 Schematic view of test setups, for cases of (a) free-field propagation, (b) scattering with cylinder, (c) scattering with NACA0012 airfoil and (d) proapgation inside nacelle and scattering into far field. Gray hashed zones represent non-reflecting boundary conditions, while black zones represent reflecting waves.

Fig. 6 Fig. 7 Fig. 8

 678 Fig. 6 Acoustic density prediction for test case 1, for the propagation of a Gaussian pulse in free-field.

Fig. 10 Fig. 11

 1011 Fig. 10 Slice of density fields at different times (top) at / = 0.5 and (bottom) / = 0.25 for networks: ( ) LBM reference, ( ) GlobUnet, ( ) GlobUnet+C, ( ) LocUnet and, ( ) LocUnet+C for test case 2.

Fig. 13 Fig. 14

 1314 Fig. 13 Slices of density fields at / = 0.5 for test case 3, for different times and networks ( ) LBM reference, ( ) GlobUnet, and ( ) LocUnet+C.

Fig. 15 Fig. 16

 1516 Fig. 15 Slices of density fields at / = 0.5 for test case 4, for different times and networks ( ) LBM reference, ( ) GlobUnet, and ( ) LocUnet+C.

Fig. 17

 17 Fig. 17 Schematic view of test setups, for cases in-duct wave propagation for different diameter-to-domain length ratios ℎ/ = (a) 0.4, (b) 0.6, (c) 0.8. . Gray hashed zones represent non-reflecting boundary conditions, while black zones represent reflecting waves.

Fig. 18 Fig. 19 Fig. 20

 181920 Fig. 18 Acoustic density prediction for test case 5, of in-duct wave propagation with ℎ/ = 0.6.

Fig. 21

 21 Fig. 21 Details of the employed U-Net [14] convolutional architecture.

Table 1 Summary of the different strategies for testing the neural network

 1 .

	Name	Dataset Normalization EPC
	GlobUnetD1	D1	No
	GlobUnetD1+C	D1	Yes
	LocUnetD1	D1	No
	LocUnetD1+C	D1	Yes

VII. Acknowledgments

The authors wish to thank Wagner Gonçalves Pinto, from ISAE-Supaero for fruitful discussion on the topic. This work was partly supported by the French "Programme d'Investissements d'avenir" ANR-17-EURE-0005. We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) through the Discovery Grant program. M. Bauerheim is also supportted by the French Direction générale de l'armement (DGA) through the AID POLA3 project. The simulations were carried out using the CALMIP supercomputing facilities (project P20035), France and Compute Canada clusters.