Eric Laloy
email: elaloy@sckcen.be

Niklas Linde

Cyprien Ruffino

Romain Hérault

Gilles Gasso

Diederik Jacques

Gradient-based deterministic inversion of geophysical data with generative adversarial networks: Is it feasible?

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

1 Introduction and Scope [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] recently proposed to use generative adversarial networks (GANs), a game changer data generation algorithm (e.g., [START_REF] Goodfellow | Generative adversarial networks[END_REF][START_REF] Goodfellow | Deep learning[END_REF], to define a low-dimensional parameterization encoding complex geologic prior models, thereby allowing efficient and accurate geostatistical inversion with Markov chain Monte Carlo (MCMC) methods [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF]. GANs have permitted impressive advancements for a wide range of applications such as image and texture synthesis, image-to-image translation and super-resolution [START_REF] Creswell | Generative adversarial networks: An overview[END_REF]. In the near future, we expect to witness a dramatic increase in development and use of GAN-inspired algorithms for geostatistical simulation (e.g., [START_REF] Mosser | Reconstruction of three-dimensional porous media using generative adversarial neural networks[END_REF][START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] and inversion (e.g., [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF][START_REF] Mosser | Stochastic seismic waveform inversion using generative adversarial networks as a geological prior[END_REF][START_REF] Richardson | Generative adversarial networks for model order reduction in seismic full-waveform inversion[END_REF]. A key element of the GAN approach is that the dimensions of the learned low-dimensional or "latent" space are independent from each other and have a known probability distribution, typically either a uniform or a standard Gaussian distribution. [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] have shown that inversions based on such parameterizations work well for global probabilistic inference of complex binary 2D and 3D prior subsurface models. However, exploring the GANderived latent space with state-of-the-art MCMC sampling [START_REF] Vrugt | Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF][START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance computing[END_REF] still necessitates tens of thousands (or more) forward evaluations [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF]. Such a computational expense can be prohibitive when using computationallydemanding forward solvers encountered in the geosciences. However, the latent space might also lend itself to a conventional deterministic gradient-based local search, which is often much more computationally frugal than global and probabilistic inversion.

In the geosciences, gradient-based deterministic inversion methods can be roughly divided into (1) methods that make use of the so-called sensitivity or Jacobian matrix, J, that is, the matrix of derivatives of the forward model outputs, F (m), with respect to the model parameters, m (e.g., physical properties of the model grid blocks): J i,j = δF (m) i δm j , and (2) methods that only require the gradient vector of the misfit (that is, objective or loss) function, L, with respect to m, ∇L with ∇L j = δL δm j , such as the steepest descent method. Although ∇L is often cheaper to calculate than J, it is generally understood that methods relying on J are more effective and robust than those based on ∇L because J can be used to (implicitely) approximate the Hessian matrix, H with H i,j = δJ i,j δm j , which contains information about the local curvature of the misfit function [START_REF] Aster | Parameter estimation and inverse problems[END_REF].

Complex spatial prior models are usually represented by a so-called training image (TI). A TI is a large gridded 2D or 3D unconditional representation of the expected target spatial field that can be either continuous or categorical (e.g., geologic facies image) and is typically used to guide geostatistical simulation by multiple-point statistics (MPS) algorithms (e.g, [START_REF] Strebelle | Conditional simulation of complex geological structures using multiple point statistics[END_REF][START_REF] Mariethoz | The Direct Sampling method to perform multiple-point geostatistical simulations[END_REF]. When using a GAN (see, e.g., [START_REF] Goodfellow | Deep learning[END_REF], for details about the GAN architecture and training) to encode the prescribed TI (or, alternatively, MPS realizations from it), a GAN-based model realization is produced by feeding the so-called generator, G, with a latent vector, z: m = G (z), and the inversion is performed within the latent space p (z). For gradientbased inversion, the required Jacobian is then the matrix of derivatives of F [G (z)] with respect to z: J z with J z i,j = δF [G (z)] i δz j . In principle, the derivatives of G (z) with respect to z, V with V i,j = δG (z) i δz j , can be directly computed by autodifferentiation [START_REF] Al-Rfou | A Python framework for fast computation of mathematical expressions[END_REF][START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems[END_REF][START_REF] Paszke | Automatic differentiation in PyTorch[END_REF]. When J is available, for instance through the use of an adjoint model, J z can be obtained as J z = JV. In our experience, calculating each element of V using the reverse-mode autodifferentiation engine that equips current deep learning (DL) libraries such as TensorFlow [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Distributed Systems[END_REF] and PyTorch [START_REF] Paszke | Automatic differentiation in PyTorch[END_REF] can be slow, especially for large m and when run on a CPU. For instance, for the PyTorch implementation of the spatial GAN [START_REF] Jetchev | Texture synthesis with spatial generative adversarial networks[END_REF] used herein it takes about 5 minutes to construct the V array on a last generation Intel ® i7 CPU and about 1 minute using a NVIDIA Quadro M2000M GPU, when z is 15-dimensional and G (z) is of size 129 × 65. More importantly, generating high-quality categorical geologic structures with a GAN may require postprocessing of the realizations through filtering and/or thresholding (e.g., [START_REF] Mosser | Reconstruction of three-dimensional porous media using generative adversarial neural networks[END_REF][START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF]. If not differentiable, this postprocessing, pp (•), decouples the V array computed by the DL code from the actual realization, pp [G (z)].

Another solution that alleviates the need of having an adjoint model and is not perturbed by postprocessing operations consists of estimating J z directly using a finitedifference approximation. This incurs a relatively low computational cost since it requires 2N z + 1 forward runs to evaluate J z with a central difference scheme, with N z the length of the z vector, which by construction is low-dimensional (say two orders of magnitude less than m).

Lastly, when using the gradient vector, ∇L with respect to z, ∇L (z), instead of the Jacobian, J z , ∇L (z) and the trained G (z) can be integrated within a single DL computational graph [START_REF] Richardson | Generative adversarial networks for model order reduction in seismic full-waveform inversion[END_REF][START_REF] Mosser | Stochastic seismic waveform inversion using generative adversarial networks as a geological prior[END_REF]. This is possible because G (z) is fully differentiable end-to-end and DL libraries are equipped with advanced gradient descent algorithms such as Adam [START_REF] Kingma | Adam: a Method for Stochastic Optimization[END_REF]. If F [G (z)] is fully differentiable then the DL library can compute ∇L (z) internally by autodifferentiation. This is for instance the case when F (•) is linear and the continuous G (z) realizations do not need postprocessing such as thresholding. As shown by [START_REF] Richardson | Generative adversarial networks for model order reduction in seismic full-waveform inversion[END_REF], when pp [G (z)] or F (•) is not differentiable ∇L (m) can however be externally computed and integrated within the DL graph. In this work, we used an adjoint model giving ∇L (m) and backpropagated ∇L (z) = ∇L (m) V where the ∇L (m) V product can be efficiently computed within a single backward pass without explicitely calculating the V array (see our companion code).

No matter whether ∇L (z) or J z is used or regardless of how they are calculated, an important open question is to what extent the nonlinearity in the G (z) transform adversely affects the convergence of gradient-based deterministic inversion performed in the latent space p (z). Working with p (z) rather the original model space p (m) conveniently ensures that (1) any generated model honors the prior TI and (2) the parameter dimensionality is reduced by orders of magnitude compared to traditional deterministic inversion (e.g., [START_REF] De Groot-Hedlin | Occam's inversion to generate smooth, twodimensional models from magnetotelluric data[END_REF]. Nevertheless, the non-linear relationship between G (z) and z adds significant nonlinearities to the inverse problem, in addition to those of the forward model, F (m). In this work, we consider a best case scenario and use a linear tomography problem based on ground penetrating radar (GPR) data to show that even if F (m) is linear with respect to m, the nonlinearity in G (z) may frequently affect the performance of gradient-based deterministic inversions. In contrast, probabilistic global optimization using the DREAM (ZS) code [START_REF] Vrugt | Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling[END_REF][START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance computing[END_REF]) is found to work well for the considered problem, although at a rather high computational cost. The overall motivation of this work is to raise awareness that although model parameterizations based on GAN algorithms can be very powerful in representing prior information and reducing model dimensionality, the resulting inverse problems may become too non-linear to enable adequate convergence of deterministic gradient-based inversions even when the forward solver itself is only weakly-nonlinear or even linear.

The remainder of this paper is organized as follows. Section 2 summarizes related work before section 3 briefly describes our used GAN and considered inversion algorithms, and details the considered inverse problem. This is followed by section 4 which presents our inversion results. In section 5, we discuss our main findings and outline possible future developments. Finally, section 6 provides a short conclusion.

Related work

To the best of our knowledge, [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] were the first to introduce and demonstrate the idea of using the latent space learned by a GAN to perform probabilistic inversion of hydrologic or geophysical data. [START_REF] Richardson | Generative adversarial networks for model order reduction in seismic full-waveform inversion[END_REF] recently proposed to embed ∇L (m) and a generator, G (z), trained beforehand within a single (fully differentiable) DL computational graph to deterministically solve a 2D seismic full-waveform inversion problem by making use of the stochastic gradient descent algorithms implemented in DL libraries. [START_REF] Mosser | Stochastic seismic waveform inversion using generative adversarial networks as a geological prior[END_REF] recently used a trained generator for solving a 2D seismic inversion problem within a probabilistic framework. They combined the estimation of ∇L (m) obtained from an adjoint model together with G (z) and the DL-based calculation of V within the same (fully differentiable) computational graph to perform Metropolis-adjusted Langevin MCMC sampling (e.g., [START_REF] Roberts | Optimal scaling of discrete approximations to langevin diffusions[END_REF]. Outside of the geosciences, [START_REF] Creswell | Inverting the generator of a generative adversarial network[END_REF] and [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF] inverted full images produced by a trained generator using usual benchmarks in computer vision and studied to what extent latent vectors can be recovered. In addition, [START_REF] Bora | Compressed sensing using generative models[END_REF] showed that gradient descent within the latent space of a GAN works overall well to recover an image from a set of measurements obtained by applying a linear operator to the true image. Conceptually, the [START_REF] Bora | Compressed sensing using generative models[END_REF] study has a lot in commons with ours. However, Bora and coworkers consider a totally different training set than us: the CelebA dataset which consists of face images of celebrities. In section 5, we compare our findings to those of [START_REF] Bora | Compressed sensing using generative models[END_REF].

Methods

Generative Adversarial Networks

We use a PyTorch implementation of the spatial GAN (SGAN, Jetchev et al., 2016) used by [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] to generate realizations from the 2D channelized aquifer training image (TI) depicted in Figure 1a. Our method is not limited to such channel images, but this example is selected to be similar to benchmark images that are typically used to test MPS algorithms. For brevity, we refer the reader to [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF] for a description of the basic network architecture and training procedure, and only provide below the main principles behind the used SGAN.

The building blocks of a GAN are the generator and discriminator. The generator, G (z), is fed with a low-dimensional latent vector, z, and produces a model realization, G (z) = m, whose spatial statstics match those found in the TI provided that the GAN training performed beforehand (see below) was successful. The z vector is trained to obey either a standard normal distribution, z ∝ N (0, I), or a truncated uniform distribution, z ∝ U (-1, 1), and the z variables are independent of each other. For image generation, the discriminator component only serves for training G (z) as detailed below.

Due to its purely convolutional nature (see, e.g., [START_REF] Goodfellow | Deep learning[END_REF], for more details on convolutions and convolutional layers in deep neural networks), for an ergodic TI our SGAN can be trained at relatively low computational cost using a small realization domain, before being used to generate arbitrarily large realizations. The latent space of our SGAN has a spatial structure, with z being reshaped into an m×n×o×q array Z for the 3D case. As detailed in [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF], the m, n and o dimensions are directly related to the three spatial dimensions while q, which we set to 1 in this work, is an extra dimension that can encode additional information about the data representation (see geostatistical simulation case study 2 in [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF]. For a square (w × h) or cubic (w × h × l) generation domain, the relationship between z x = m = n = o and m x = w = h = l is given by m

x = [z x -1] 2 dp + 1, (1)
where dp is the number of stacked convolutional layers in G (z). This allows for a rather strong dimensionality reduction. For example, when z x = 5 and dp = 5, we have m x = 129. We refer to [START_REF] Jetchev | Texture synthesis with spatial generative adversarial networks[END_REF] for more details on the (2D) SGAN architecture.

In a GAN, the generator and discriminator are trained (or "learned") simultaneously with opposing goals. The discriminator, D (m), is fed with samples from the "real" training set, which from now on will be referred to as m true with distribution p data (m), and "fake" samples (i.e., realizations) created by the generator: m = G (z). In our case, the real samples m true are a set of patches randomly cut from the TI. The discriminator tries to distinguish between m true and m by computing, for each received sample, the probability that it belongs to p data (m true ). Conversely, the generator, G (z), aims at fooling D (•) into labeling m as a sample from p data (m true ) [START_REF] Goodfellow | Generative adversarial networks[END_REF]. This translates into the following minimization-maximization problem min

G(•) max D(•) E mtrue∼p data (mtrue) [log (D (m true ))] + E z∼pz(z) [log (1 -D (G (z)))] , (2)
Compared to the SGAN code used by [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF], we replaced batch normalization by instance normalization and added two transposed dilated convolutional layers to the generator, before its output (see section 8). Upon selection of the appropriate training epoch (see below), this removes the need to post-process the realizations by median filtering to eliminate small impurities, thereby leaving us with thresholding as the unique postprocessing operation. The realizations produced by our trained G (z) are continuous on the [0, 1] range and they are converted into binary images by thresholding at the 0.5 level.

Training was performed on a GPU Tesla K40 for 50 epochs with 64 mini-batches of 100 training samples per epoch. Epoch 36 was deemed to produce the best realizations. Figures 1b-c show two (randomly chosen) realizations generated by the trained SGAN before thresholding, tr [G (z)]. These 513 × 513 realization are obtained by sampling a 289-dimensional z vector from p (z) ∝ U (-1, 1) (see [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF], for details and a performance comparison with a popular MPS algorithm). Even if these images are continuous, they appear almost categorical in this representation, which also suggests that the subsequent thresholding is a relatively mild operation.

Gradient-Based Deterministic Inversion

A common representation of the forward problem is

d = F (m) + e, (3)
L (m) = Φ D + λΦ M , (4)
Φ D = 1 2 [F (m) -d] T C -1 e [F (m) -d] , (5)
Φ M = 1 2 [m -m prior ] T C -1 m [m -m prior ] , (6)
where m prior and C m signify the a priori model and its covariance matrix, respectively, and where λ weights the influence of honoring prior Gaussian information about m on L (m). In case of a linear relationship between d and F (m), F (m) reduces to the product Fm. In practice, it is common to replace C m by a regularization operator, L. This leads to the joint minimization of the misfit, Φ D , and

Φ M = ||L (m -m prior ) || 2 2 L (m) = 1 2 [F (m) -d] T C -1 e [F (m) -d] + λ [m -m prior ] T L T L [m -m prior ] . (7)
The popular Gauss-Newton (GN) method for minimizing equation ( 7) iteratively updates m until either the target data misfit or a maximum number of iterations has been reached using the following update mechanism

m k+1 = m k + ∆m = J T k C -1 e J T k + λL T L -1 J T k C -1 e dk + m prior , (8)
with dk = d -F (m k ) + J k (m -m prior ) , (9)
where k denotes the iteration number. In the context of our GAN-based dimensionality reduction, we learn a p (z) = U (-1, 1) for the latent vector z. Limited testing with learning a standard normal distribution, N (0, 1), is found to provide good results but the generated realizations are nevertheless of slightly lower quality than when using U (-1, 1). Fortunately, when the components of a random vector, x are independent, they can be converted into standard normal variables using the iso-probabilistic transform: labeling as τ x (x i ) (Ψ (ξ i )) the cumulative density function (CDF) of x i (a standard normal variable ξ i ), the direct and inverse transform are given by

ξ i = Ψ -1 [τ x (x i )] and x i = τ -1 x [Ψ (ξ i )] . (10)
To perform the inversion within the GAN latent space, we therefore place a standard normal uncorrelated prior (so-called damping regularization) on a vector z SN of the same size as z and achieve GN updates of z SN while converting each proposed z SN into the corresponding z ∝ U (-1, 1) with equation ( 10) before creating the resulting G (z) realization. Under the standard normal prior on z SN , the GN update simplifies to

z SN,k +1 = J T z,k C -1 e J T z,k + λI -1 J T z,k C -1 e d -F G τ -1 z [Ψ (z SN,k )] + J z,k z SN,k (11) where J z,k = J z SN k
is the Jacobian matrix of z SN,k , and I is the identity matrix. In this work, we simply set λ = 1 as adjusting λ dynamically did not show any advantages over fixing λ = 1.

When only the gradient of the objective or loss function with respect to the model parameters, ∇L (m), is used, a steepest descent type of update can be performed

m k+1 = m k -α∇L (m) k , ( 12
)
where α is a step length parameter that is generally adapted dynamically.

In this study, we used the Adam algorithm [START_REF] Kingma | Adam: a Method for Stochastic Optimization[END_REF] implemented in Pytorch to update z k using the p (z) = U (-1, 1) prior directly. Adam is an advanced gradient-based algorithm that only makes use of the loss function gradient (see Kingma and Ba, 2015, for details). We used the following algorithmic parameters, a learning rate of 0.01 and the default values of 0.9 and 0.999 for the β 1 and β 2 momentum parameters. Since our considered forward operator, F (m), is linear (see section 3.4), it implies that for continuous m = G (z) realizations it is straightforward to have the DL library internally computing ∇L (z) (at least for common loss functions L (•)). Yet this is no longer possible if the G (z) realizations require thresholding. To circumvent this issue, we implemented an adjoint model giving ∇L (m) within the Pytorch computational graph [START_REF] Richardson | Generative adversarial networks for model order reduction in seismic full-waveform inversion[END_REF]. Here ∇L (m) is backpropagated through the generator network to update, using Adam, the latent vector, z, that produced the model, m = G (z). As explained earlier, this is possible because the generator is fully differentiable end-to-end and ∇L (z) = ∇L (m) V (see also our companion code for more details).

A simple sum of squared errors (SSR, see equation 14) was chosen for the loss function, L (m). Consequently, specific steps are needed to deal with the fact that the estimates in z can move outside the support of the U (-1, 1) prior. Three options to deal with this situation have been investigated by [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF] in the context of direct inversion of a full (face) image: (i) allow z values to leave the [-1, 1] interval, (ii) replace components that are too large with the maximum allowed value and components that are too small with the minimum allowed value (i.e., standard clipping), (iii) reassign the exceeding z components uniformly at random in the [-1, 1] range (i.e., stochastic clipping). In agreement with the findings by [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF], limited testing revealed that stochastic clipping produces slightly better results for our considered case studies. Therefore, stochastic clipping is adopted throughout.

Global Optimization

For comparison purposes, we used the DREAM (ZS) algorithm [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance computing[END_REF][START_REF] Vrugt | Markov chain Monte Carlo simulation using the DREAM software package: theory, concepts, and MATLAB implementation[END_REF] as global optimizer. This is a MCMC sampler designed to sample the posterior density function of the parameters. Yet it is herein used to solely find an appropriate solution rather than sampling the full posterior distribution. The DREAM (ZS) scheme runs parallel interacting Markov chains and is thus population-based, and makes use of differential evolution principles to propose candidate points. It has been proven efficient in many hydrologic and geophysical applications (see, e.g., Vrugt, 2016, for references). Full algorithmic details can be found in [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance computing[END_REF] and [START_REF] Vrugt | Markov chain Monte Carlo simulation using the DREAM software package: theory, concepts, and MATLAB implementation[END_REF]. Using a population-based search method requires finding the right trade-off between the size of the population (or "exploration") and the number of generations (or "exploitation") for a given number of forward model runs. After some testing, we selected 8 Markov chains, meaning a 8-member population, and performed as many MCMC iterations, that is, population generations, as necessary to fit the data to the assigned noise level.

Inverse Problem

We consider a synthetic 2D straight-ray cross-hole GPR tomography exemple. Given that the physics is assumed to be fully linear, any convergence problems arising in solving the inverse problem can be attributed to the nonlinear relationship between z and G (z). We use two binary images with two facies of homogeneous GPR velocity (0.06 m ns -1 and 0.08 m ns -1 ) as true models (see Figures 2a, 3a, 4a, and5a). Both are produced using m = G (z) with z randomly sampled from U (-1, 1). Cross-hole GPR imaging uses a transmitter antenna to emit a high-frequency electromagnetic wave at a location in one borehole and a receiver antenna to record the arriving energy at a location in another borehole (e.g., [START_REF] Annan | GPR methods for hydrogeological studies[END_REF]. The considered measurement data are first-arrival traveltimes for several transmitter and receiver locations. These data contain information about the GPR velocity distribution between the boreholes. The GPR velocity primarily depends on dielectric permittivity, which is strongly influenced by volumetric water content and, consequently, porosity in saturated media. Our setup consists of two vertical boreholes that are located 6.0 m apart. Sources (left) and receivers (right) are located between 0.5 and 12.0 m depth with 0.5 m spacing (Figures 2a, 3a, 4a, and5a), leading to a total dataset of d N = 576 traveltimes. Under the linear physics assumption, synthetic traveltime data, d, are simulated as

d = Am + e, (13)
where A = J contains the path length in each model cell, e represents independent random draws from a zero-mean homoscedastic Gaussian distribution with a typical standard deviation σ e = 1 ns. When a simple SSR loss function (without regularization) is used

L (m) = ||d -Am|| 2 ≡ L (z) = ||d -AG (z) || 2 , ( 14
)
then the gradient vector of L (m) with respect to m becomes

∇L (m) = -2A T [d -Am] T ≡ ∇L (G (z)) = -2A T [d -AG (z)] T . ( 15
)
As explained above, this allows us to use Pytorch to backpropagate ∇L (m) (equation ( 15)) through the GAN generator and optimize z such that equation ( 14) is minimized within a single computational graph using gradient descent. From now on, this will be referred to a strategy 1. In addition, estimating J z SN k at each iteration k by a centered 2point finite difference scheme together with using a GN search as described by equations 7-11 will be called strategy 2. Here J z SN is derived using a perturbation factor of 0.1. This perturbation value was found to be the most appropriate after testing with several values in the [0.01, 0.5] range. A perturbation of 0.1 for a N (0, 1) prior might seem large. Yet one must keep in mind that the necessary thresholding of the produced realizations, G (z), can cause one or more sensitivies in J z SN to be zero or near-zero if the used perturbation is too small. Such excessively small J z SN i,j values will induce instabilities in the GN search described by equation ( 11).

Lastly, note that besides true models I and II we use two different realizations of e in equation ( 13) which leads to two different "true" datasets for each true model. Realization I induces a true root mean square error (RMSE) between "true" (corrupted) and uncorrupted measurements of 1.006 ns. For realization II, the corresponding RMSE is 1.001 ns. Repeating the inversions with different noise realizations is useful to check the robustness of the inversion results against the exact noise values. Overall, this leads to 4 inverse case studies: 2 true models × 2 noise realizations. To ease the comparisons between the different inversion runs, we use the weighted root mean square error (WRMSE) which we define as the ratio of the achieved RMSE by the inversion to the true RMSE. A WRMSE > 1 ns thus basically means that the data are underfitted while WRMSE = 1 denotes an appropriate data misfit.

Results

Inverting the Full Image

Before considering application of the forward operator, we tested whether gradient descent within the latent space of our trained GAN is able to recover the true (nonthresholded) G (z) model and corresponding z vector from the inversion of G (z). In other words, is it possible to recover the values in z given exhaustive knowledge of the true G (z)? This task is similar to the one considered by [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF], but in a geophysical context rather than for face images. The considered two true models in [0, 1] (I and II) generated by our trained GAN are depicted in Figures 2a and3a. No data corruption was used and 50 inversion trials using 10,000 Adam iterations were performed with randomly chosen starting points. Table 1 andFigures 23summarize the results. In Table 1, a RMSE of 0.02 between the true and reconstructed images is selected to define a successfull reconstruction. This choice is based on the observation that, visually, the retrieved models are either very similar to the true images with a RMSE < 0.02 or largely differ from the true images with a RMSE > 0.2 -0.3. Furthermore, similarly to [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF], we report the percentage of inversion results for which the mean squared error (MSE) between the true and inferred z values is < 0.01. Overall, 37% of the inversion runs recover the true models (Table 1). Furthermore, each of these successful runs corresponds to a successful recovery in the latent space for true model I (Table 1 andFigure 2) but less for model II (Table 1 andFigure 3). As illustrated by Figures 3b and3e, this is because a slight linear bias in the z values can nevertheless produce good image reconstructions. Our results somewhat differ from those obtained by [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF] who report 100% recovery in the z-space to arbitrary precision. There are two major differences between our work and that of [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF]. First, they consider face images which are much more clustered compared to our channelized model: the eyes, nose, ears, etc. tend to occupy the same portion of the image and they are located at globally similar distances from each other. Second, [START_REF] Lipton | Precise recovery of latent vectors from generative adversarial networks[END_REF] performed 100,000 gradient descent iterations while we performed 10 times less iterations.

Geophysical Inversion

We now turn our attention to the actual geophysical inversion, in which a subset of linear averages are observed under noise. For each strategy (1 and 2) and each combination of true model (I or II) and noise realization (I and II), the inversion was repeated 100 times using 100 different (randomly chosen) starting models. For strategy 1, the maximum total number of iterations (i.e., proposal steps) was set to 10,000. Here each Adam iteration consumes one forward solve. Since the resulting computational expense may be overly large for realistic nonlinear geophysical inverse problems, inversion results are also looked at after a total of 100 and 1000 iterations. Due to the used centered 2-point finite-difference approximation of J z SN (see above), for strategy 2 one iteration incurs a cost of 31 forward solves. Therefore, the maximum number of iterations for strategy is set to 5000 which translates into 155,000 forward model runs. In addition, inversion performance is also considered after 100 (3100 forward solves) and 1000 (31,000 forward solves) iterations. Moreover, to verify the impact of thresholding on the quality of the inferred models, each case study is also run without thresholding.

Strategy 1 without tresholding relies on the direct computation of ∇L (z) by autodifferentiation with no adjoint model involved and is therefore basically the same as the approach by [START_REF] Bora | Compressed sensing using generative models[END_REF]. For the continuous case, we consider the continuous versions of the true models shown in Figures 2a and3a and the same noise realizations as for the binary case. To convert a generated continuous model m = G (z) in [0, 1] into a continuous velocity field m VEL (m ns -1 ) in [0.06, 0.08], the following relationship is used:

m VEL = 0.06 + 0.02 (1 -m)
The results of the 100 repetitions of the 8 different inversions are summarized in Table 2 for the continuous case and in Table 3 for the binary case. Furthermore, Figures 4567depict for each strategy applied to the binary case, true models I and II, and noise realization II, the true model (subplot (a)), two best-fitting models found among the 100 repetitions (subplots (b) and (c)), the best-fitting model across the allowed iterations of a randomly chosen trial for which the corresponding best WRMSE is > 1.2 ns (subplot (d)), the sampled WRMSE trajectory corresponding to the best-fitting model shown in subplot (b) (subplot (e)) and the sampled WRMSE trajectory associated with the model displayed in subplot (d) (subplot (f)).

Table 2: Continuous case: number of inversion runs among the 100 performed ones that achieved a WRMSE ≤ 1.2, ≤ 1.1 and ≤ 1, respectively, for each inversion strategy (1 and 2), total number of iterations, N iter , combination of true model (I and II) and noise realization (I and II). The WRMSE is defined for each inversion run as the ratio of the best RMSE (over the iterations of a given run) to the RMSE of the true data (1.006 ns for noise realization I and 1.001 ns for noise realization II). It is observed for both strategies 1 and 2, that the quality of the best-fitting solution heavily depends on the starting model and the allowed number of iterations. Indeed for the smallest budget of 100 iterations, no satisfying solutions are found no matter the inversion approach. We find that strategy 1 (Figures 4-5 and Tables 2-3) shows a much higher success rate than strategy 2 (Figures 6-7 and Tables 23) for the largest number of iterations allowed to each method. Hence strategy 1 retrieves solutions with a WRMSE ≤ 1.01 in 4% to 22% of the trials for the continuous case and in 7% to 19% of the trials for the binary case. On the contrary, strategy 2 never does so (Tables 23).

Strategy

Strikingly, there appears to be very little negative impact of thresholding (compare Table 2 to Table 3), except for strategy 2 and true model I (compare the bottom row of Table 2 with those of Table 3). This indicates that plugging the externally computed ∇L (m) into the DL graph to permit gradient descent works well for the considered examples.

When the maximum number of iterations is reduced to 1000, the proportion of solutions associated with a WRMSE ≤ 1.1 becomes small but is overall two times larger for strategy 1 (2% on average over the 4 datasets) than for strategy 2 (1% on average over the 4 datasets). Since a model with a WRMSE as large as 1.2 may still visually ressemble the true model well (see Figure 6c), it is instructive to also consider a WRMSE ≤ 1.2 as threshold. Here strategy 2 becomes more attractive, but for true model I only (Tables 2-3). We have tried to improve strategy 2 by (i) damping the model update at the early iteration before increasingly adding more contrast to the proposed models and (ii) decreasingly smooth the proposed models after each iteration as the search progresses. Nevertheless, neither of those schemes proved successful.

There is a systematic difference in peformance between the cases with true models I and II, with true model I always inducing a larger success rate no matter the inversion strategy and computational effort (Tables 2 and3). This is due to the fact that model I has higher prior probability or, to use deep learning terminology, is closer to the range of the generator [START_REF] Bora | Compressed sensing using generative models[END_REF] than model II. This can be informally seen by comparing models I (Figure 2a) and II (Figure 3a) with the original TI (Figure 1a): the patterns displayed by model I are occuring more frequently in the TI than the eye-type of channel branching contained in model II. It is seen that strategy 1 offers greater robustness than strategy 2 against variations in prior model probability, as shown by the smaller differences in success rates between the model I and II cases for strategy 1 (Tables 2 and3).

With respect to sampling trajectories of the inversions, for strategy 1 most of the data misift reduction in a productive run (i.e., a run with a favorable starting model) is achieved within the first 3000-5000 iterations (Figures 4e and5e). The main reduction in data misfit for strategy 2 can occur both rather quickly (< 500 iterations) or much later (Figures 6e and7e), while this strategy hardly recovers models with WRMSE ≤ 1.1. In addition, the sampling behavior of both algorithms can be unstable (Figures 4e, 5e, 6e and 7e). In accordance with those findings, it is also noted for strategy 1 that if WRMSE values larger than 2-3 are still sampled after some 2000 forward simulations, then the considered run will be unable to obtain a high-quality solution (Figures 4f and5f).

It is worth noting that using strategy 2 with a trained GAN for which the latent space obeys a standard normal distribution, N (0, 1), thereby removing the additional nonlinearity caused by the Ψ (z SN,k ) transform, leads to globally similar results as those presented in Tables 2 and3 for strategy 2.

Even for our best performing strategy 1 and considered maximum number of iterations of 10,000, the overall success (WRMSE ≤ 1.01) rate remains relatively low (Tables 2 and3). As posited earlier in this paper, we argue that the reason for this is the high degree of nonlinearity associated with the relationship between z and G (z). More insights into the effect of the nonlinear G (z) transform on the considered simple inverse problem are provided by Figure 8. The latter displays a 2D slice in the WRMSE land-scape corresponding to the true reference solution obtained by the combination of the true model I (see Figure 2a, 4a or 6a) and noise realization II. To construct Figure 8, the z 3 to z 15 components of the true z vector that induces model I are kept fixed while the z 1 and z 2 dimensions are gradually varied between -1 and 1. While the WRMSE response surface appears smooth when considered at a coarse resolution (Figure 8a), a zoom into the basin of attraction (Figure 8b) reveals that the problem is in fact multimodal with many local minima. Also, this WRMSE landscape contains many small spikes and some large flat areas. Lastly, we find that performing a global optimization with DREAM (ZS) parameterized as described in section 3.3 provides a solution with a WRMSE of 1.000 for all of the 4 scenarios (true models I and II combined with noise realizations I and II). Nevertheless, this is at the cost of a serial total of 200,000 to 800,000 forward runs. Among these amounts, some 50,000 serial forward calls are consumed to reach a WRMSE of about 1.1 while the rest is spent by moving from 1.1 to ≤ 1. This points out that traveling across the relatively rough misfit terrain caused by the nonlinearities in G (z) may require global optimization or search methods even in the best case of a linear forward problem. Although not done herein, note also that population-based methods such as DREAM (ZS) are straightforward to evolve in parallel using one CPU per population member, which can save a large amount of time [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF].

Discussion

This paper demonstrates that even if the forward model is linear, the high degree of nonlinearity in the relationship between latent vector, z, and the generator of a GAN, G (z), can cause gradient-based deterministic inversion performed in the GAN latent space p (z) to fail. This is here illustrated for a classical binary channelized aquifer TI.

Less structured and/or continuous TIs might be more amenable to gradient-based deterministic inversion within p (z) since the G (z) transform might then be less nonlinear. The poor convergence to the target misfit is in stark contrast to the results obtained when inverting the full image without noise (Table 2, Figures 23). This suggests that limited sampling of spatial averages under noise further complicates the recovery of the true z values.

Our findings might appear to contradict those of [START_REF] Bora | Compressed sensing using generative models[END_REF] who investigated whether gradient descent within the latent space of a trained GAN can recover an image from a set of measurements obtained by applying a linear operator to the true image. In cases where this true image is itself created by running the trained GAN for a given z vector (similarly as in our study), [START_REF] Bora | Compressed sensing using generative models[END_REF] systematically find small reconstruction errors. Nonetheless, Bora and coworkers consider face images (from the Celeba dataset) which are thus very different than our channelized aquifer training image. Face images are comparatively more clustered as the same structures (nose, ears, eyes,...) are more or less always located in the same areas of the image. Yet the channels in our training set can be anywhere (see [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF], for the full 2500 × 2500 training image from which the training set is sampled). This larger degree of freedom in our images (for this aspect) might make the inversion problem harder to solve. In addition, [START_REF] Bora | Compressed sensing using generative models[END_REF] do not give much details about their exact experimental setup which makes an in-depth comparison with our results difficult.

It might seem surprising that, if a sufficient number of iterations is allowed, strategy 1, which is based on the objective function gradient, ∇L (m), outperforms a GN search that uses a finite-difference approximation of the full Jacobian matrix in the p (z) space, J z . In the training of deep neural networks, stochastic gradient descent has proven highly successful. [START_REF] Goodfellow | Deep learning[END_REF] argue that such first-order methods may more easily escape saddle points, while second-order methods tend to get attracted and stuck in them. The situation might be similar here where we try to estimate z for an already trained network, namely that the less "stable" gradient-based methods can more easily move away from saddle-points and possibly local minima. Furthermore, strategy 1 relies on a rather powerful stochastic gradient descent algorithm (Adam, [START_REF] Kingma | Adam: a Method for Stochastic Optimization[END_REF] and uses a perfectly accurate estimation of ∇L (m) since in the considered case, ∇L (m) is available analytically. On the contrary, when approximated by finite differencing J z is inevitably fraught with errors.

Note that for the gradient descent used with strategy 1, no optimization of the Adam's hyperparameters, that is, the learning rate, and β 1 and β 2 momentum parameters (see [START_REF] Kingma | Adam: a Method for Stochastic Optimization[END_REF], was done. Only the learning rates of 0.1 and 0.001 were found to lead to a worse performance than using our choice of 0.01. Other Adam's hyperparameter settings or other descent algorithms such as SGD, RMSPROP (see, e.g., [START_REF] Goodfellow | Deep learning[END_REF] or Adabound [START_REF] Luo | Adaptive gradient methods with dynamic bound of learning rate[END_REF] might enabled improved inversion quality for the considered case studies.

Lastly, notice that some recent work suggests that projecting the learned latent space onto a so-called "Riemannian" manifold may lead to a reduced non-linearity of G (z) with respect to the projected z [START_REF] Shao | The Riemannian geometry of deep generative models[END_REF][START_REF] Arvanitidis | Latent space oddity: on the curvature of deep generative models[END_REF][START_REF] Chen | Metrics for deep generative models[END_REF]. In the future, it might be worthwhile to attempt inversion within such a manifold.

Concluding remarks

Performing inversion within the low-dimensional, latent space of a trained GAN, p (z), rather than within the original model space, p (m), offers important advantages in that (1) any generated model proposal honors the prior training image (TI) and (2) the parameter space is reduced by orders of magnitude compared to the original model space. Global probabilistic inversion within p (z) has been recently shown to work well [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF]. However, such probabilistic inversion still incurs a large computational cost while the learned p (z) could also be possibly used within computationally less demanding gradient-based deterministic inversions. In this work, we show that owing to the highly nonlinear relationship between the GAN generator, G (z), and latent vector, z, gradient-based deterministic inversion may fail even though the physics of the forward problem is linear. For a channelized aquifer binary TI and a synthetic linear GPR tomography problem involving 576 measurements with low noise, we find that when allowing for a total of 10,000 iterations, about 13% of the trials by the considered gradient descent algorithm locate a solution that has the required data misfit, in average over the considered test cases. When restricting the maximum allowed number of iterations to 1000, the percentage of success becomes 0%, though approximately 2 % of the inferred models still closely ressemble their associated true model and induce a nearly correct data misfit. Furthermore, the tested Gauss-Newton inversion approach that approximates the Jacobian matrix to create updates proved to be unable to recover a solution with the appropriate data misfit. Overall, deterministic inversion performance is found to depend on the inversion approach, starting model, true reference model, number of iterations and noise realization. In contrast, costly global optimization with DREAM (ZS) always finds an appropriate solution.

Computer Code Availability

The GAN and associated inversion codes used in this study are available at https: //github.com/elaloy/gan_for_gradient_based_inv.

Used GAN architecture

Compared to the generator architecture detailed in [START_REF] Laloy | Training-image based geostatistical inversion using a spatial generative adversarial neural network[END_REF], we have the following differences.

• Batch normalization was replaced by instance normalization • The nonlinearity in the 5 th layer of the generator is a ReLU instead of an hyperbolic tangent.

• A first transposed dilated convolutional layer takes the output of the 5 th layer as input. This 6 th layer is parameterized with 64 output channels, a kernel size of 5 pixels, a stride of 1 pixel, a padding of 6 pixels, an output padding of 0 pixels and dilation coefficient of 3. The activation function in this 6 th layer is a ReLU and instance normalization is applied after the nonlinearity.

• A second transposed dilated convolutional layer takes the output of the 6 th layer as input. This 7 th layer is parameterized with 1 output channel, a kernel size of 5 pixels, a stride of 1 pixel, a padding of 10 pixels, an output padding of 0 pixels and dilation coefficient of 5. The nonlinearity in this 7 th layer is a hyperbolic tangent.

Figure 1 :

 1 Figure 1: (a) Fraction of size 513 × 513 of the used 2500 × 2500 binary TI and (b) -(c) randomly chosen 513 × 513 realizations derived by our SGAN. Each continuous realization is generated by sampling 289 random numbers from a uniform distribution, U (-1, 1), without any postprocessing.

  where d = (d 1 , . . . , d N ) ∈ R N , N ≥ 1 is the number of measurement data, F (m) denotes a deterministic forward model with parameters or model m and the noise term e lumps all sources of errors. If the probability distribution of e can be assumed to be zero-mean Gaussian with covariance matrix C e , and m is assigned a multiGaussian prior distribution, then the optimal model, m, minimizes the following objective or loss function, L (m)

Figure 2 :

 2 Figure 2: Inversion results for the full image inversion using the true model I: (a) true model, (b) and (c): retrieved models with successful image recovery, (d) retrieved model with unsuccessful image recovery, (e -g) scatter plot of the recovered z vector associated with the model shown in subplot (b -d) and the true z vector, respectively.

Figure 3 :

 3 Figure 3: Inversion results for the full image inversion using the true model II: (a) true model, (b) and (c): retrieved models with successful image recovery, (d) retrieved model with unsuccessful image recovery, (e -g) scatter plot of the recovered z vector associated with the model shown in subplot (b -d) and the true z vector, respectively.

Figure 4 :

 4 Figure 4: Binary case: inversion results for strategy 1, truel model I and noise realization II: (a) true model, (b) and (c) two randomly chosen best-fitting models found over 100 repetitions, (d) best-fitting derived model across the allowed iterations of a randomly chosen repetition for which the corresponding best WRMSE is > 1.2 ns, (e) sampled WRMSE trajectory for the best model found among the 100 repetitions (displayed in subplot(b)), and (f) sampled WRMSE trajectory associated with the model depicted in subplot (d). The red triangles and the green squares in subfigures (a-d) represent the GPR source and receiver positions, respectively.

Figure 5 :

 5 Figure 5: Binary case: inversion results for strategy 1, truel model II and noise realization II: (a) true model, (b) and (c) two randomly chosen best-fitting models found over 100 repetitions, (d) best-fitting derived model across the allowed iterations of a randomly chosen repetition for which the corresponding best WRMSE is > 1.2 ns, (e) sampled WRMSE trajectory for the best model found among the 100 repetitions (displayed in subplot(b)), and (f) sampled WRMSE trajectory associated with the model depicted in subplot (d). The red triangles and the green squares in subfigures (a-d) represent the GPR source and receiver positions, respectively.

Figure 6 :

 6 Figure 6: Binary case: inversion results for strategy 2, truel model I and noise realization II: (a) true model, (b) and (c) two randomly chosen best-fitting models found over 100 repetitions, (d) best-fitting derived model across the allowed iterations of a randomly chosen repetition for which the corresponding best WRMSE is > 1.2 ns, (e) sampled WRMSE trajectory for the best model found among the 100 repetitions (displayed in subplot(b)), and (f) sampled WRMSE trajectory associated with the model depicted in subplot (d). The red triangles and the green squares in subfigures (a-d) represent the GPR source and receiver positions, respectively.

Figure 7 :

 7 Figure 7: Binary case: inversion results for strategy 2, truel model II and noise realization II: (a) true model, (b) and (c) two randomly chosen best-fitting models found over 100 repetitions, (d) best-fitting derived model across the allowed iterations of a randomly chosen repetition for which the corresponding best WRMSE is > 1.2 ns, (e) sampled WRMSE trajectory for the best model found among the 100 repetitions (displayed in subplot(b)), and (f) sampled WRMSE trajectory associated with the model depicted in subplot (d). The red triangles and the green squares in subfigures (a-d) represent the GPR source and receiver positions, respectively.

Figure 8 :

 8 Figure 8: 2D slice in the WRMSE landscape corresponding to the true reference solution obtained by the combination of the true model I (see Figure 2a, 4a or 6a) and noise realization II. The (a) subplot shows a coarse scale resolution over the full WRMSE range while the (b) subplot presents a zoom into the main basin of attraction. The red cross in the (b) subplot denotes the true values of z 1 and z 2 in this case.

Table 1 :

 1 Percentage (%) of inversion runs that achieved a RMSE ≤ 0.02 in recontructing the full 125 × 60 image (G (z)) and a MSE ≤ 0.01 in recovering the true 15-dimensional z vector, for each true model (I and II).

	True model RMSE < 0.02 MSE < 0.01
	1	34	32
	2	40	8

Table 3 :

 3 Binary case: number of inversion runs among the 100 performed ones for the continuous case that achieved a WRMSE ≤ 1.2, ≤ 1.1 and ≤ 1, respectively, for each inversion strategy (1 and 2), total number of iterations, N iter , combination of true model (I and II) and noise realization (I and II). The WRMSE is defined for each inversion run as the ratio of the best RMSE (over the iterations of a given run) to the RMSE of the true data (1.006 ns for noise realization I and 1.001 ns for noise realization II).

	Strategy	N iter	True model Noise WRMSE ≤ 1.2 WRMSE ≤ 1.1 WRMSE ≤ 1.01
	1	100	I	I	0	0	0
	1	100	I	II	1	0	0
	1	100	II	I	0	0	0
	1	100	II	II	0	0	0
	1	1000	I	I	5	3	0
	1	1000	I	II	4	3	0
	1	1000	II	I	3	2	0
	1	1000	II	II	1	0	0
	1	10,000	I	I	21	19	13
	1	10,000	I	II	24	23	19
	1	10,000	II	I	10	9	7
	1	10,000	II	II	24	17	12
	2	100	I	I	22	0	0
	2	100	I	II	8	3	0
	2	100	II	I	0	0	0
	2	100	II	II	0	0	0
	2	1000	I	I	100	0	0
	2	1000	I	II	57	4	0
	2	1000	II	I	3	0	0
	2	1000	II	II	2	0	0
	2	5000	I	I	100	0	0
	2	5000	I	II	97	5	0
	2	5000	II	I	6	0	0
	2	5000	II	II	9	1	0


Alicja Domaracka

Grégoire Danger

Multiple beam irradiation platform MIRRPLA: origin and evolution of organic matter in the solar system

"Where do we come from?" is one of the important fascinating open questions of science and philosophy. How did life emerge? What is the origin of organic matter in the universe? Could life also emerge on other worlds than our Earth? Complex organic molecules have indeed been observed in space (comets, meteorites, molecular clouds). In 2023, France launched the Programmes et équipements prioritaires de recherche (PEPR) Origins (life and universe), through the program leader the National Centre for Scientific Research, to address these questions [START_REF]PEPR Origins[END_REF]. Understanding the origin of primitive organic matter during the formation and evolution of the solar system is fundamental because the contribution of extraterrestrial organic matter via asteroids and comets is one of the possible sources of organic matter available to primitive Earth [START_REF] Caselli | Our astrochemical heritage[END_REF].

In astrophysical environments, icy layers are omnipresent (e.g., micrograins in molecular clouds) and are constantly exposed to multiple radiations (ultraviolet [UV], X-rays, electrons, and ions from solar winds, cosmic rays, and ions trapped in magnetospheres) [START_REF] Caro | Prebiotic chemistry in icy grain mantles in space. An experimental and observational approach[END_REF]. Following radiolysis, complex organic molecules can be synthesized [START_REF] Caselli | Our astrochemical heritage[END_REF][START_REF] Caro | Prebiotic chemistry in icy grain mantles in space. An experimental and observational approach[END_REF]. These different types of radiation act simultaneously, but existing laboratory astrophysics experiments are generally performed using a single radiation source and synergistic effects are not studied. It is therefore necessary to carry out simultaneous irradiation with several beams (photons, ions, electrons) to simulate a realistic space radiation field. Within the PEPR Origins, a unique multibeam irradiation platform (UV photons, keV electrons and keV-GeV ions delivered by the Grand Accélérateur National d'Ions Lourds (GANIL, Caen France) equipped with an infrared spectrometer and a cold head to prepare layers of ices made of mixtures of small molecules will be constructed. Complex organic molecules formed during the irradiation and/or heating of ice layers will be detected using gas phase chromatography and in-situ mass spectrometry. This new versatile instrument, the MIRRPLA platform (Multiplebeam IRRadiation PLAtform to investigate the origin and evolution of organic matter of the Solar System), will be open to the various scientific communities (astrophysics, radiobiology, environmental and materials sciences) as well as to industry via the GANIL-Centre de Recherche sur les Ions, les Materiaux et la Photonique (CIMAP) -Centre for Interdisciplinary Research with Heavy Ions (CIRIL) user platform [START_REF] Rangama | CIRIL: Interdisciplinary Research at GANIL[END_REF]. Figure 1 shows the involved partners and organizations. Figure 2 the unique principle of MIRRPLA.

The uniqueness of the proposed setup consists in the combination of different types of ionizing radiation (photons, electrons, ions), including the possibility of simultaneous irradiation of samples: for example, GANIL beam plus photons or UV electrons to simulate a complex irradiation field, as in space (Figure 2). The modifications induced by the energy treatment will be monitored by Fourier transform infrared absorption spectroscopy. The species emitted to the gaseous phase via a sputtering or desorption mechanism will be monitored by quadrupole mass spectrometry. In addition, a high-resolution mass spectrometer coupled with gas chromatography (CG-Orbitrap) will be connected to the irradiation chamber to study the volatile organic compounds synthesized during the irradiation of ice astrophysical analogs.

In addition to studies of synergistic effects with realistic space radiation fields on icy surfaces, MIRRPLA will also provide organic residues (analogs of primary matter from space) to other PEPR Origins partners for benchmarking chemical analysis techniques or comparison to returned space samples. Beyond the questions related to the origin of complex (prebiotic) organic species in the universe, MIRRPLA also will have an impact in other fields of societal importance. To give some examples, health issues such as fundamental processes in mixed irradiation of biomolecules and deoxyribonucleic acid fragments to develop new cancer treatment protocols can be addressed. Concerning environmental science, polycyclic aromatic hydrocarbons are produced in combustion and are key intermediates in the formation and growth of soot particles. In addition, carbon nanoparticles can be covered with water ice in the upper part of the atmosphere and be exposed to different types of ionizing radiation. For materials sciences, this unique facility makes it possible to study (1) the aging of materials under mixed irradiation environments (nuclear fuel cells, reactor materials, ceramics/glass for the storage of radioactive waste, polymers), and (2) simultaneous multiradiation effects on materials used in spacecraft and space missions (hull, electronics).

The kickoff meetings of MIRRPLA and Origins were held on 26 June 2023 at CIMAP-GANIL, Caen, and on 18-21 September 2023, at Collège de France, Paris, respectively.

Figure 1 .

 1 Figure 1. Partners involved in MIRRPLA: CIMAP (Centre de Recherche sur les Ions, les Materiaux et la Photonique, Caen, France) and laboratory PIIM (Physique des Interactions Ioniques et Moléculaires, Marseille, France).

Figure 2 .

 2 Figure 2. Schematic representaion of the unique principle of MIRRPLA.

Acknowledgment: This work benefits from state aid managed by the French National Research Agency under France 2030 investment plan (reference ANR-22-EXOR-0012).


David Babonneau

Sophie Camélio

Lionel Simonot

Frédéric Pailloux

Philippe Guérin

Bruno Lamongie

Olivier Lyon

S Camelio

Tunable plasmonic dichroism of Au nanoparticles

Keywords: 61.46.-w -Structure of nanoscale materials 78.67.-n -Optical properties of low-dimensional, mesoscopic, and nanoscale materials and structures 81.16.-c -Methods of micro-and nanofabrication and processing

come

Introduction. -Planar arrays of metallic nanostructures such as nanoparticles or nanowires periodically aligned on dielectric surfaces are ideal for fundamental studies of coupling and ordering phenomena as well as for applied research on plasmonics [1], photovoltaics [2], and magnetism [3]. Especially, noble metal nanostructures have received considerable attention lately due to a range of possible applications as dichroic filters, molecular sensors, and subwavelength photonic waveguides exploiting the selective photon absorption and local field enhancement arising at wavelengths near the so-called surface plasmon resonance (SPR) condition [4]. Despite the progress in nanofabrication techniques, the realization of nanostructures with well-controled composition, size, shape, and distribution remains a significant challenge, especially for the large-scale production of regular arrays on substrates that may be applied as dichroic filters [5][6][7]. Lithographic techniques could potentially offer a solution [8,9], but are not suitable for most practical purposes owing to the low output and high production costs. Recent advances in nanostructure preparation using colloidal chemistry show that self-alignment of metallic nanorods may be realized, providing an easy and inexpensive way for tailoring the plasmonic properties of metamaterials [10][11][12][13]. Another effective and low-cost method to fabricate large-area ar-rays of Au and Ag nanostructures consists in combining substrate prepatterning by defocused ion-etching with subsequent Volmer-Weber growth by physical vapor deposition at glancing angles [14][15][16][17]. It has been proved that these materials possess promising optical properties due to their morphological anisotropy reflected in a sharp polarization-dependent excitation of their localized surface plasmons. The ability to tune the optical dichroism of Au [15] and Ag [17] nanowires by varying the amount of deposited metal (t metal ) has been proved recently. However, there is still no report on the influence of t metal on the growth mechanisms as well as the resulting optical properties in the low metal-coverage regime (i.e., when disconnected nanoparticles are formed). In this article, we demonstrate that the optical dichroism of Au nanoparticle arrays grown on rippled alumina thin films can be easily tuned either in the visible (t Au ≤ 1 nm) or both in the visible and the near infrared regime (t Au ≥ 2 nm). The relationship between such optical properties and the structure (morphology and organization) of the nanoparticles is discussed.

Experimental. -Rippled surfaces were produced on 90 nm-thick amorphous Al 2 O 3 thin films deposited by ionbeam sputtering onto fused-silica substrates and carbon- coated copper grids. Ion-etching was performed at room temperature with Xe + ions at a sputtering angle of θ = 55 • with respect to the surface normal, with an ion energy E = 1000 eV, and a fluence of about 3.6 × 10 16 ions cm -2 (ion flux ∼ 2 × 10 14 ions cm -2 s -1 ), leading to ripples that are oriented perpendicular to the ion beam [16]. The rippled thin films were used as templates with the aim of creating self-organized chains of Au nanoparticles. Deposition of Au was accomplished at a temperature of 200 • C and a rate of 0.01 nm s -1 by DC magnetron sputtering under a glancing incidence of 5 • from the surface, with the atomic flux being oriented opposite to the Xe + beam. The effective Au thickness was varied between t Au = 0.5 nm and t Au = 3.0 nm in order to avoid the formation of Au nanowires. The rippled thin films before and after Au deposition were examined by grazing incidence small-angle x-ray scattering (GISAXS) on the SWING beamline at SOLEIL synchrotron with different orientations of the incident x-ray beam with respect to the ripples. The energy of the incident x-ray beam was 11.857 keV, the angle of incidence was α i = 0.3 • , and the scattered photons were collected as a function of the in-plane 2θ f and out-ofplane α f exit angles with a two-dimensional CCD detector [Fig. 1(a)]. The planar arrays of Au nanoparticles were also characterized by high-angle annular dark-field scanning transmission electronic microscopy (HAADF-STEM) with a JEOL 2200FS microscope using an acceleration voltage of 200 kV, a probe size of 0.7 nm, and an inner collection angle of 50 mrad. The optical properties were investigated by means of spectroscopic transmission of polarized light in the range between 400 and 950 nm by using a QE65000 Ocean Optics TM spectrophotometer.

Results. -Figure 1(b) shows the GISAXS pattern of an as-etched Al 2 O 3 thin film recorded with the incident xray beam along the x-direction perpendicular to the Xe + beam. Two vertical streaks on either side of the z-axis are readily observed, indicating a periodic surface corrugation (namely ripples) with a wave vector along the y-direction parallel to the Xe + beam. The position and width of the streaks are related to the period L y and correlation length ξ y of the surface corrugation, while the vertical extend of the streaks is related to the ripple height h [16]. The quantitative analysis of the GISAXS intensity was performed in the framework of the distorted wave-Born approximation with the FitGISAXS program [18] assuming a rippled surface with a symmetric sawtooth profile as shown in Fig. 1(c) where tan γ = 2h/L y . The full GISAXS map was fitted within the local monodisperse approximation by using the one-dimensional paracrystal model [19] to retrieve L y and ξ y , and by assuming a log-normal distribution of L y with a mean period ⟨L⟩ = L y exp[ 3 2 (ln σ) 2 ], and a dispersion ln σ. , the in-plane distribution of the Au nanoparticles becomes strongly anisotropic with short-range order in the x-direction and long-range order in the y-direction. Moreover, the period of the chains of nanoparticles corresponds to that of the ripples, L y , suggesting that the periodicity is imposed by the prepatterned Al 2 O 3 template [14][15][16][17]. The in-plane size and shape distributions of the Au nanoparticles were analyzed from the plane-view HAADF-STEM images by assuming hemiellipsoids with in-plane projected sizes D x and D y along the xand y-directions, respectively. The mean values ⟨D x ⟩ and ⟨D y ⟩ together with the corresponding standard deviations σ x and σ y are reported in Tab. 1. It is worth noting that not only the effective in-plane size D = √ ⟨D x ⟩⟨D y ⟩, but also the shape anisotropy of the nanoparticles (⟨D x ⟩/⟨D y ⟩) increase with t Au . Actually, for t Au = 3.0 nm, agglomeration and coarsening proceed leading to the formation of elongated nanoparticles predominantly aligned along the direction of the ripples (i.e., perpendicular to the Au atomic flux). Furthermore, it is worth noting that whereas the size dispersion in the y-direction (σ y /⟨D y ⟩) tends to reach a saturation value of ∼ 0.6 with increasing t Au , the size distribution in the x-direction (σ x /⟨D x ⟩) dramatically broadens. As seen in Fig. 2(d)-2(f), the structural anisotropy of such planar arrays of Au nanoparticles results in a dependence of their optical transmittance on the polarization of the incident light. For y-polarization (transverse excitation), the transmittance T y shows a localized SPR whose spectral position changes from 540 nm to 610 nm when t Au increases from 0.5 nm [Fig. 2(d)] to 3.0 nm [Fig. 2(f)]. This behavior is characteristic of disconnected noble-metal islands supported on dielectric films and the spectral shift can be ascribed to a decrease of the effective out-of-plane ratio of the nanoparticles defined by H/D [20]. In addition, a progressive decrease of the transmittance minimum is readily observed, which is due to the increase of the nanoparticle volume and of the surface coverage. For xpolarization (longitudinal excitation), the transmittance T x still shows a localized minimum when t Au = 0.5 nm and t Au = 1.0 nm. However, rotating the polarization from transverse to longitudinal results in a slight broadening of the SPR band together with a red-shift by approximately 30 nm, which cause an optical dichroism in the visible range only (Fig. 3). For t Au = 3.0 nm, although the nanoparticles are still disconnected, the longitudinal transmittance T x displays a spectral behavior characteristic of infinite metal nanowires, with a strong decrease toward the near infrared region and no evidence of localized SPR [21,22]. The resulting optical dichroism of the nanocomposite films is evidenced in Fig. 3, which exhibits the difference in transmission ∆T = T y -T x divided by the transmittance T y measured with transverse excitation.

It appears that at high coverage (t Au ≥ 2 nm), depending on the polarization, an efficient harvesting of light over a broadband spectrum [23] might be achieved either in the visible (y-polarization) or in the near infrared regime (x-polarization).

To go further in the structural characterization of the Au nanoparticles, Figs. 4(a) and 4(b) show the GISAXS patterns of the nanocomposite film with t Au = 0.5 nm recorded with the incident x-ray beam along the directions parallel (x-direction) and perpendicular (y-direction) to the ripples. In contrast to the HAADF-STEM observations [Fig. 2(a)], the GISAXS measurements reveal a strong structural anisotropy. While the GISAXS pattern obtained with the x-ray beam along the y-direction exhibits two broad lobes expected for a planar distribution of nanoparticles with short-range order, the GISAXS pattern obtained with the x-ray beam along the x-direction displays intense and sharp streaks located at a position in excellent agreement with the corresponding GISAXS pattern before Au deposition [Fig. 1(b)]. However, it can be seen that the intensity distribution becomes strongly asymmetric after Au deposition, thus demonstrating that the scattering of the Au nanoparticles dominates the total GISAXS intensity as expected from the high contrast Table 1: Structural parameters retrieved from HAADF-STEM and GISAXS analyses as a function of the Au effective thickness tAu. In-plane mean projected sizes ⟨Dx⟩ and ⟨Dy⟩ and corresponding standard deviations σx and σy (from HAADF-STEM). Out-of-plane sizes of the nanoparticles H, interparticle distances Lx, and corresponding correlation length ratios ξy/ξx (from GISAXS).

t Au (nm) ⟨D x ⟩ (nm) σ x (nm) σ x /⟨D x ⟩ ⟨D y ⟩ (nm) σ y (nm) σ y /⟨D y ⟩ H (nm) L x (nm) ξ y /ξ x 0.5 3.8( 5 , with the scattered intensity concentrating toward the origin of the reciprocal space owing to the growth of the nanoparticles. The GISAXS patterns were analyzed assuming chains of hemiellipsoidal nanoparticles with in-plane dimensions ⟨D x ⟩ and ⟨D y ⟩, as determined from HAADF-STEM (Tab. 1). Also, it has been considered that the distance between chains is L y and the correlation length is ξ y as extracted by GISAXS from as-etched Al 2 O 3 thin films. First, the GISAXS maps recorded with the incident x-ray beam along the y-direction were fitted to determine the out-of-plane dimension H of the nanoparticles as well as the interparticle distance L x and the corresponding correlation length ξ x . The results reported in Tab. 1 confirm that the effective out-of-plane ratio of the nanoparticles H/D decreases with t Au , while the interparticle distance L x increases progressively. Furthermore, in agreement with the HAADF-STEM analyses, it appears that the ξ x values are much smaller than the corresponding ξ y values. Also, it can be seen that the proximity of the nanoparticles in the x-direction is strongly affected, with the dimensionless parameter ⟨D x ⟩/L x increasing from 0.69 to 0.97. In contrast, ⟨D y ⟩/L y which increases from 0.15 to 0.51 is always smaller than 2/3, the value above which interparticle coupling effects become meaningful [24]. It is thus tempting to ascribe the red-shift of the longitudinal SPR to the increase of both the in-plane shape anisotropy and of the interparticle coupling effects, and its broadening to the strong increase of the dispersions of particle size and separation in the x-direction. Finally, the parameters obtained previously were used to calculate the GISAXS intensity with the incident x-ray beam along the x-direction. Neglecting the ripple scattering and assuming the Au nanoparticles to be in the ripple valleys leads to a symmetric GISAXS pattern [Fig. 5(a)], which is in disagreement with the experimental results. However, by considering nanoparticles tilted from the surface normal by a tilt angle of γ, the asymmetry of the experimental GISAXS map is well reproduced [Fig. 5(b)]. These results indicate that, although the planar distribution is almost isotropic at low coverage rates, a number of Au nanoparticles nucleate and then grow preferentially on the slopes of the rippled surface illuminated by the Au atomic flux during glancing-angle deposition, thus leading to a replication of the long-range order between ripples in the y-direction at higher coverage rates.

Conclusion. -

In conclusion, we have presented evidence for the production of arrays of aligned Au nanoparticles using glancing-angle deposition on self-organized rippled Al 2 O 3 templates. HAADF-STEM and GISAXS experiments combined with quantitative analysis provide complementary information on the influence of the shadowing effects on the growth mechanisms. We showed that the (in-plane and out-of-plane) shape anisotropy of the nanoparticles and their proximity to each other can be modified by varying the amount of deposited metal, which allows a dichroic and tunable optical response due to the anisotropic excitation of localized plasmons. For low Au coverage (t Au ≤ 1 nm), almost isotropic nucleation and growth of the nanoparticles lead to the achievement of a weak optical dichroism in the visible range. For high Au coverage (t Au ≥ 2 nm) and longitudinal excitation, these nanostructured materials exhibit a broadband spectrum typical of infinite nanowire arrays, although the Au nanoparticles are electrically disconnected. Accordingly, it is possible to design plasmonic nanostructures capable of an efficient harvesting of light either in the visible or in the near infrared regime depending on the polarization of the incident light. The method offers many advantages (applicable to a wide range of materials, possibility to tune the periodicity and the morphology of the aligned nanoparticles by changing the prepatterning and deposition conditions, fast output and low production costs, etc.) that make it attractive for numerous plasmonic applications. * * *

Fig. 1 :

 1 Fig. 1: (a) Schematic drawing of the GISAXS geometry. (b) Experimental GISAXS pattern of a rippled Al2O3 thin film with the incident x-ray beam along the x-direction. (c) Schematic drawing of the model surface profile. (d) Simulated GISAXS pattern with the parameters obtained from the fit of the experimental map.

  Figures 2(a)-2(c) show the plane-view HAADF-STEM images and corresponding autocorrelation functions of similar rippled Al 2 O 3 thin films after glancing-angle Au deposition, with effective thicknesses of t Au = 0.5 nm, 1.0 nm, and 3.0 nm, respectively. Our observations indicate that despite shadowing effects, the initial stages of Au growth [Fig.2(a)] is characterized by an almost isotropic nucleation owing to the high surface diffusivity of Au on Al 2 O 3 at 200 • C. In contrast, at higher coverage rates [Figs. 2(b) and 2(c)], the in-plane distribution of the Au nanoparticles becomes strongly anisotropic with short-range order in the x-direction and long-range order in the y-direction. Moreover, the period of the chains of nanoparticles corresponds to that of the ripples, L y , suggesting that the periodicity is imposed by the prepatterned Al 2 O 3 template[14][15][16][17]. The in-plane size and shape distributions of the Au nanoparticles were analyzed from the plane-view HAADF-STEM images by assuming hemiellipsoids with in-plane projected sizes D x and D y along the xand y-directions, respectively. The mean values ⟨D x ⟩ and ⟨D y ⟩ together with the corresponding standard deviations σ x and σ y are reported in Tab. 1. It is worth noting that not only the effective in-plane size D = √ ⟨D x ⟩⟨D y ⟩, but also the shape anisotropy of the nanoparticles (⟨D x ⟩/⟨D y ⟩) increase with t Au . Actually, for t Au = 3.0 nm, agglomeration and coarsening proceed leading to the formation of elongated nanoparticles predominantly aligned along the direction of the ripples (i.e., perpendicular to the Au atomic flux). Furthermore, it is worth noting that whereas the size dispersion in the y-direction (σ y /⟨D y ⟩) tends to reach a saturation value of ∼ 0.6 with increasing t Au , the size distribution in the x-direction (σ x /⟨D x ⟩) dramatically broadens.

Fig. 2 :

 2 Fig. 2: Plane-view HAADF-STEM images of rippled Al2O3 thin films after glancing-angle Au deposition. The effective Au thicknesses are (a) tAu = 0.5 nm, (b) tAu = 1.0 nm, and (c) tAu = 3.0 nm. The directions of the Xe + beam and Au flux are indicated by the arrows, and the autocorrelation functions of the images are shown in the insets. (d-f) Corresponding transmission measurements at normal incidence with the polarization of the incoming light being oriented parallel (xpolarization) and perpendicular (y-polarization) to the ripples.

Fig. 3 :

 3 Fig. 3: Difference in transverse and longitudinal transmittances ∆T = Ty -Tx divided by the transverse transmittance Ty with increasing effective thicknesses from tAu = 0.5 nm to tAu = 3.0 nm.

Fig. 4 :

 4 Fig.4: Experimental GISAXS patterns of rippled Al2O3 thin films after glancing-angle Au deposition. (a) tAu = 0.5 nm and x-rays along the x-direction, (b) tAu = 0.5 nm and x-rays along the y-direction, (c) tAu = 1.0 nm and x-rays along the xdirection, (d) tAu = 1.0 nm and x-rays along the y-direction, (e) tAu = 3.0 nm and x-rays along the x-direction, (f) tAu = 3.0 nm and x-rays along the y-direction.

Fig. 5 :

 5 Fig. 5: Simulated GISAXS patterns (tAu = 1.0 nm and x-rays along the x-direction) by assuming that the Au nanoparticles (a) in the ripple valleys and (b) tilted from the surface normal by a tilt angle of γ.

We would like to acknowledge assistance from SOLEIL staff during the GISAXS experiments at the SWING beamline.