Martin Schmelzer
email: m.schmelzer@tudelft.nl

Richard P Dwight
email: r.p.dwight@tudelft.nl

Wouter N Edeling
email: wedeling@stanford.edu

Paola Cinnella
email: paola.cinnella@ensam.eu

Estimation of Model Error Using Bayesian Model-Scenario Averaging with Maximum a Posterori-Estimates

Keywords: CFD, UQ, RANS, Turbulence Modelling, Bayesian Calibration, Bayesian Scenario-Model Averaging

HAL is

Introduction

The understanding of turbulence is one of the key challenges in classical mechanics. A turbulent flow has a three-dimensional, time-dependent and random velocity field [START_REF] Pope | Turbulent Flows[END_REF], which is composed of a wide range of scales varying from the level close to molecular dissipation of heat to the scales determined by the boundary conditions of the flow domain. Despite great developments in the field of high-fidelity turbulence modelling, such as improvements of the large eddy simulation (LES) approach and hybrid-methods combining Reynolds-averaged Navier-Stokes (RANS) and LES, RANS 'continues to be the standard approach used to predict a wide range of flows for very complex configurations across virtually all aerospace product categories' [START_REF] Slotnick | CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences[END_REF]. However, using the less-computationally demanding RANS approach comes at the price of uncertainty due to approximate physical modelling such as turbulence modelling.

In recent years research has focussed on two types of epistemic uncertainties in relation to turbulence modelling, namely parameter uncertainty and model-form uncertainty. Parameter uncertainty arises due to the fact that the closure coefficients of the models are determined via calibration against simple flow configurations using experimental data or scale-resolving simulations. A Bayesian calibration incorporating both the error of the experimental data and model inadequacy revealed how strongly the posterior distributions of the closure coefficients vary even for a simple flow scenario of flow over a flat plate, with different pressure gradients ranging from favourable to strongly adverse [START_REF] Edeling | Bayesian estimates of parameter variability in the k-epsilon turbulence model[END_REF]. This observation makes the generalisation of the coefficient unjustifiable. Similarly, a Bayesian analysis dealing with the predictive performance of commonly used turbulence models, i.e. Launder-Sharma kε, Wilcox kω, Spalart-Allmaras, Baldwin-Lomax and Stress-ω, showed that no superior model could be identified either for the given flow configurations [START_REF] Edeling | Bayesian estimates of parameter variability in the k-epsilon turbulence model[END_REF].

Bayesian Model-Scenario Averaging (BMSA) uses a set of different closure models to compute both an ensemble prediction as well as a-posteriori measures of uncertainty due to the choice of closure model. In addition, the variability of the closure coefficients over different calibration scenarios is included to inject uncertainty caused by applying a set of coefficients to a predictive case for which they were not specifically calibrated [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF]. The full BMSA approach requires the propagation of various posterior distributions through a CFD code, an expensive endeavour when the underlying flow case is computationally expensive. The costs could be mitigated by using surrogate models per turbulence model in order to propagate the distributions more efficiently [START_REF] Edeling | Quantification of modelling uncertainties in turbulent flow simulations[END_REF]. However, also the construction of the surrogate can become expensive for models with many closure coefficients.

In this work we use a major simplification of the full BMSA approach in order to make the technique accessible for industrial flow cases. The reduction of the costs is based on using maximum a posteriori (MAP) estimates of the posterior distributions of the closure coefficients, which means that only a single set of closure coefficients per scenario and per model needs to be propagated through the code [START_REF] Edeling | Pragmatic estimation of model error of rans turbulence models for industrial relevant flow cases[END_REF]. Furthermore, the BMSA method was developed using data of flat plate boundary layer flows, for which the Bayesian inference of the posterior closure coefficient was conducted with a cheap boundary-layer code. The method has recently been applied to airfoil and wing cases as examples in order to assess the predictive capabilities of the BMSA framework for wall-bounded flows in external aerodynamics [START_REF] Edeling | Pragmatic estimation of model error of rans turbulence models for industrial relevant flow cases[END_REF], which represents a natural next step in terms of flow complexity. In this work we apply the method to other flow cases such as flow in a turbulent pipe at Re = 44, 000 and flow over periodic hills at Re H = 5600 in order to assess the robustness of the method also for flow configurations outside of this natural domain. Finally, we also show results for one of the industrial challenges of the UMRIDA project: a generic Falcon Jet (IC-03).

Bayesian predictive methodology 2.1 Calibration

The BMSA framework is based on the Bayesian calibration of the coefficients of closure models. Treating calibration as a stochastic problem we obtain posterior probability distributions for the coefficients, which serve as our uncertain estimates of the coefficients under the measurement error of the reference data and the modelling error. The data used for the calibration consists of boundary-layer data-sets from the 1968 AFOSR-IFP-Stanford conference proceedings [START_REF] Coles | Computation of turbulent boundary layers[END_REF], which are both highly resolved and subject to low measurement noise. The set contains a collection of wide range of favourable and adverse pressure gradients.

A scenario, denoted S, is a particular flow set-up, including boundary-conditions, material parameters, and all other physical properties needed to define the flow, with corresponding experimental data z. The CFD code m CFD (S; M, θ θ θ ) takes as arguments the scenario S ∈ S = {S 1 , . . . , S K }, a turbulence model M ∈ M = {M 1 , . . . , M I } and its closure-coefficients θ θ θ . Given an operator H z (•) that maps the state to the measured quantities z we can define a statistical model to relate z and θ θ θ :

z = δ [H z • m CFD (S; M, θ θ θ )] + ε, (1)
where ε, δ are random-variables (RVs) representing measurement noise and multiplicative model error respectively. The error of the measurements is modelled as zero-mean additive Gaussian noise and the RV δ is a representation for the model error, which following the approach of Cheung et. al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] specifies δ as a correlated Gaussian process, see also [START_REF] Edeling | Bayesian estimates of parameter variability in the k-epsilon turbulence model[END_REF]. The Gaussian choice for both ε and δ yields a Gaussian likelihood function, i.e. the pdf describing the probability of observing the data given a realization of θ θ θ . Finally, an application of Bayes theorem [START_REF] Gelman | Bayesian data analysis[END_REF]

Here, p(θ θ θ |M, S) is the prior distribution which represents our knowledge of θ θ θ prior to training M on z. Since the denominator in (2) does not depend on θ θ θ it is omitted from consideration.

To obtain samples from the posterior distribution p(θ θ θ |z, M, S), we employ the Markov-Chain Monte Carlo method [START_REF] Hastings | Monte carlo sampling methods using markov chains and their applications[END_REF]. To reach convergence of the Markov-chain, we observered that roughly 40,000 code samples were required [START_REF] Edeling | Bayesian estimates of parameter variability in the k-epsilon turbulence model[END_REF]. Ordinarly this would constitute an excessive strain on available computational resources in a CFD context. However, as our experimental data consisted of boundary-layer quantities, we were in a position to use a fast boundary layer code. As such, no real computational bottleneck exists during the calibration phase.

Prediction

Let ∆ be a quantity of interest (QoI) in a particular scenario, which may be a scalar, vector, or functional quantity derived from the flow-state, which can be expressed as

∆ H ∆ • m CFD ( S; M, θ θ θ ), (3)
in which S represents a flow scenario outside of the set of scenarios used within the calibration phase. We wish to obtain a stochastic estimate of ∆ conditional on a set of models M and a set of training scenarios S for the predictive scenario S / ∈ S . The BMSA methodology offers the evaluation of the posterior predictive distribution (ppd) for ∆ in case S conditional on all training data:

p(∆ | S, z) = I ∑ i=1 K ∑ k=1 Θ Θ Θ i i i p(∆ | S, M i , θ θ θ ) p(θ θ θ | S k , M i , z k ) P(M i | S k , z k ) P(S k ) dθ θ θ . (4)
The first term on the right-hand side inside the integral represents the probabilistic equivalent of the simulation results for the QoI given the flow scenario S, a turbulence model M i and closure coefficients θ θ θ . The second and third terms are the posterior probability density distribution of the closure coefficients and the posterior model probabilities respectively, which are the output of the previously conducted calibration procedure summarized in Section 2.1 and detailed in Edeling et al. [START_REF] Edeling | Bayesian estimates of parameter variability in the k-epsilon turbulence model[END_REF][START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF]. The last term represents the scenario probabilities. The solution of the multi-dimensional integral over the closure coefficients is the expensive part of the method, because it requires as many code calls as the entries in the Markov-chain, in this case 40, 000.

Therefore to obtain a practical estimate of p(∆ | S, z) we propose to approximate the marginal posterior probability distributions p(θ θ θ | S k , M i , z k ) with Dirac-δ functions at their maximum a posteriori (MAP) values4

θ θ θ MAP i,k := argmax θ θ θ ∈Θ Θ Θ i i i p(θ θ θ | S k , M i , z k ) (5)
so that

p(θ θ θ | S k , M i , z k ) δ θ θ θ -θ θ θ MAP i,k . (6)
The effect of this approxmation is to neglect the effect of within-model within-scenario variance on the ppd. Thus the ppd variance will be reduced, but still include the effect of multiple models and scenarios. Note that if perfectly plentiful data were available in the training scenarios (and the models were able to fit the data exactly for some values of the closure coefficients), then p(θ θ θ | S k , M i , z k ) would be δ -functions. So one way to think of this approximation is as neglecting the effect of imperfect training. Substituting ( 6) into (4) leads to an approximation of the posterior predictive distribution

p(∆ | z) p(∆ | z) p(∆ | S, z) = I ∑ i=1 K ∑ k=1 P(M i | S k , z k ) P(S k ) Θ Θ Θ i i i p(∆ | S, M i , θ θ θ ) δ θ θ θ -θ θ θ MAP i,k dθ θ θ = I ∑ i=1 K ∑ k=1 P(M i | S k , z k ) P(S k ) p(∆ | S, M i , θ θ θ MAP i,k ) (a) = I ∑ i=1 K ∑ k=1 P(M i | S k , z k ) P(S k ) δ ∆ -m CFD S; M i , θ θ θ MAP i,k . (7)
Equality (a) follows from the fact that the prediction of m CFD is deterministic for deterministic θ θ θ . The approximate ppd is therefore a weighted-sum of I × K δ -functions, one at each prediction of m CFD ( S) for each model, and each scenario's MAP-estimate of θ θ θ . The cost of evaluating the ppd is

I × K runs of m CFD ( S).
The first moment of p(∆ | S, z) can be derived directly from [START_REF] Coles | Computation of turbulent boundary layers[END_REF]. The expectation is

E[∆ | S, z] = ∆ • p(∆ | S, z) d∆ (a) = I ∑ i=1 K ∑ k=1 P(M i | S k , z k ) P(S k ) ∆ • δ ∆ -m CFD S; M i , θ θ θ MAP i,k d∆ (b) = I ∑ i=1 K ∑ k=1 P(M i | S k , z k ) P(S k ) m CFD S; M i , θ θ θ MAP i,k , (8)
where (a) follows from ( 7) and (b) is the integral-identity encoding the statement that the mean of a deterministic quantity is the quantity itself. Unlike the posterior model probabilities P(M i | S k , z k ), which are informed using the reference data during the training phase [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF], the scenario probabilities P(S k ) are of a predictive nature. If reference data for the scenario S would be available, P(S k ) could be optimized accordingly. However, in a predictive setting this data isn't always available. Therefore, P(S k ) is defined based on model agreement per scenario S k : If the models show a high level of agreement regarding the value of ∆ under a specific scenario S k , this scenario receives a higher weight compared to other scenarios, in which the models rather disagree. For that we exploit the principle, that if S is similar to S k the models are expected to give accordingly similar predictions, because θ θ θ MAP i,k has been calibrated under the same data z k . As introduced in [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF] this principle is modelled by

P (S k ) := ξ -p k ∑ K j=1 ξ -p j , ξ k = I ∑ i=1 m CFD S; M i , θ θ θ MAP i,k -E[∆ | S, z k ] 2 , (9)
in which p serves as a tuning parameter scaling the weighting procedure, i.e. for p = 0 all scenarios are equally weighted but for p → ∞ a single scenario is preferred. Given all terms in [START_REF] Coles | Computation of turbulent boundary layers[END_REF] we now want to deduce uncertainty estimates from the ppd. If the full range of P(M i | S k , z k ) P(S k ) is used, the uncertainty bound is largest and determined by the extreme predictions of the entire ensemble. However, a more reasonable approach is to draw samples from the ppd and use percentiles as min/max levels to obtain a confidence interval of the ppd leading potentially to tight uncertainty bounds. This is done in the following.

3 Incompressible pipe flow at Re = 44, 000

We apply the BMSA method to turbulent flow in a straight pipe at Re = 44, 000 with validation data from a DNS simulation [START_REF] Wu | A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow[END_REF]. This flow case is, together with turbulent flow over a flat plate and channel flow, one of the three canonical wall-bounded flow types [START_REF] Pope | Turbulent Flows[END_REF]. The main difference between flow in a pipe or channel and over a flat plate, lies in the development of the boundary layer. Assuming an uniform inflow into a circular straight pipe a boundary layer at the wall develops and its thickness increases downstream similarly to the flat plate case. Forced by the geometry of the pipe the boundary layers from all sides will eventually merge in the middle of the pipe. The distance from the entrance of the pipe to the point, where the boundary layer fills the entire diameter, is known as the entry length. Beyond the entry length the flow is homogeneous in stream-wise direction. In contrast the boundary layer thickness of a flat plate flow with zero pressure gradient is not bounded. However, for all these wall-bounded flow types the physics close to the wall can be assumed to be similar. Therefore, with the application of BMSA based on flat plate flow to pipe flow is a proof-of-concept in order to assess the predictive capabilities of BMSA for wall-bounded flows in general.

Simulation setup

The validation of the BMSA method is based on a comparison with averaged velocity profiles from a DNS simulation [START_REF] Wu | A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow[END_REF]. Therefore, the steady-state RANS simulation mimics the time-averaged DNS simulation by using appropriate boundary conditions given in Table 1. In the DNS simulation the mass flow rate is kept constant via a timevarying adjustment of the pressure gradient in the stream-wise direction [START_REF] Wu | A direct numerical simulation study on the mean velocity characteristics in turbulent pipe flow[END_REF]. For the RANS simulation the mass flow rate is kept constant via a fixed inlet velocity and the pressure gradient is set to zero at the inlet. At the outlet both the velocity gradients are set to zero and pressure is set to zero. In this way the pressure drop is computed according to the fixed mass flow. For the forward-simulations m CFD M i , θ θ θ MAP i,k the results of simulations with nominal values of the coefficients m CFD (M i , θ θ θ o i,k ) are used as an initial condition in order to reduce iteration counts. 5Table 1. The Inlet conditions of velocity U, pressure P/ρ and turbulent quantities. [START_REF] Schwarze | CFD Modellierung. Grundlagen und Anwendungen bei Strömungsprozessen[END_REF] U

[m/s] P/ρ [m 2 /s 2 ] k ε ω ν Inlet (1, 0, 0) T n • ∇P = 0 0.00375 0.00835 25.0516 0.00027 Outlet n • ∇U = 0 (0, 0, 0) T n • ∇k = 0 n • ∇ε = 0 n • ∇ω = 0 n • ∇ ν = 0
Wall (0, 0, 0) T (0, 0, 0) T (0, 0, 0) T (0, 0, 0) T 8367238.00803 (0, 0, 0) T

The geometry of a straight pipe is determined by its diameter and length. While the diameter D is set in order to achieve the target Reynolds number, the length of the pipe needs to be chosen sufficiently longer than the entry length in order to overcome all effects at the inlet, such as the usage of different closure coefficients or the uniform inflow velocity profile. The entry length is estimated by an empirical relation based on the Reynolds number [START_REF] Venkateshan | Mechanical Measurements[END_REF] L turb = 1.359

• D • Re 1/4 D ⇒ L turb ∼ 20D. ( 10
)
For all simulations a conservative value L = 200D was chosen. Due to the symmetry of the case a wedge-shaped mesh with symmetric boundary conditions orthogonal to the stream-wise direction was chosen [START_REF]OpenFOAM: OpenFOAM User Guide[END_REF]. Different meshes were initially studied, with a total number of 2900, 5900 and 9900 cells respectively. The differences for the velocity profiles between the meshes were negligible, but in order to keep the error small the finest mesh was chosen for all forward simulations. In order to study the effect of the turbulence modelling and not additional modelling due to wall-functions the thickness of the smallest cells close to the wall was defined according to y + ≤ 1.0 and Low-Re turbulence models were used. Applying the BMSA framework to the set of simulations leads to the expectation of the radial velocity E(U/U bulk |z) as shown in Figures 1(a ing parameters of the smart-weighting method p = [0, 1, 2]. For the confidence interval 10000 samples were drawn from the posterior predictive mass function given in eq. 7 for each x/R position. The 5th and 95th percentile of the population was used as the lower and upper limit respectively, containing 90% of the samples. For p = 0, i.e. uniform weighting of the scenarios, the impact of outliers on the confidence interval is large leading to an over-prediction of the modelling error. Outliers are caused by sensitivities of the models with respect to certain coefficients leading to predictions very different compared to the ensemble mean. The effect is already reduced for p = 1, for which the confidence interval is covering a tight bound around the expectation. The shape matches with the shape of the majority of the ensemble of forward simulations. For increasing p the confidence interval becomes tighter, since the smart-weighting consecutively increases the weight of S k = 1300 as shown in Figure 2(b) and decreases the influence of the others, reducing the overall variance. Interestingly, the MAP estimates for scenario S k = 1300 are acquired based on a flow over a flat plate with a moderate favourable pressure gradient and a boundary layer close to an equilibrium state [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF], which is very similar to the conditions for the flow in a pipe. A comparison of the DNS data and the expectation for different p via the L 2 -norm

| | U/U bulk | S,z -U o | | 2 (a) L 2 error
||E[U/U bulk | S, z] -U o || 2 (11)
given in Figure 2(a), where the high-resolution DNS data is interpolated by a cubic spline to enable a point-by-point comparison, shows a minimal error for 0 < p < 2, but also an increase of the error for higher values of p. Note that due to the use of MAP estimates, the variance is missing the within-model within-scenario contribution. Therefore, a reduction of the variance with smart-weighting should be done with care in order to not reduce the already underestimated variance even further. Thus, the purpose of the smart-weighting for the BMSA-with-MAP-approach should be to exclude outliers, which can be achieved with relatively low p values. In such a way an optimum of the width of the confidence interval and the correctness of the expectation can be achieved.

BMSA-prediction for pressure drop

Another relevant quantity of interest for the turbulent pipe flow is the pressure drop ∆ P over the pipe length. Table 2 gives the values for the pressure drop ∆ P/ρ normalised by the fluid density ρ and the corresponding Darcy friction coefficient f Darcy

f Darcy = 2 ∆ P D ρ L U 2 (12)
per model. The empirical reference value for the pressure drop is calculated from the Colebrook-White equation, which is an approximation of the empirical Moody-diagram,

1 √ f = -2 log 10 2.51 Re √ f ( 13
)
with zero roughness at Re = 44, 000. For several tuning parameters of the smartweighting methods, i.e. p ∈ [0, 8], the BMSA expectation of the pressure drop E[∆ P | S, z] is shown in Figure 3 as well as the empirical reference data calculated based on the Colebrook-White equation. Up to p = 4 the confidence interval still covers a range including also the empirical reference value, which suggests that the method successfully gives a reasonable uncertainty range. For larger tuning parameters the confidence interval reduces drastically and puts all weight on one single scenario, so that the confidence interval is not visible any more. In line with the results for the velocity profile (see Section 3.2), also for the pressure drop p=2 is a reasonable choice. The expectation consistently underestimates the true value and is therefore worse than the baseline predictions, see Table 2. However, since the BMSA method wasn't trained on pressure drop data this is not unexpected and the main gain of BMSA are the confidence bounds instead of replacing the baseline prediction with the point-estimate of the expectation.

Incompressible flow over periodic hills at Re H = 5600

The periodic hill test case is based on the channel flow case, but with a lower surface modified by a series of periodically arranged hills, and deals with flow separation on the curved surface of the hills and reattachment on the flat plate region between the hills [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF]. Being computationally relatively cheap but still challenging in terms of flow physics it has been used in several workshops as a benchmark test case [START_REF] Manceau | Report on the 10th joint ERCOFTAC (SIG-15)/ IAHR / QNET-CFD Workshop on Refined Turbulence Modelling[END_REF]. We use DNS data of the mean flow field from Breuer et al. [START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] for Re H = 5600, defined using the streamwise bulk velocity U b between the hills crest and the upper surface and the hill height H. The test case is known to be especially challenging for linear eddy-viscosity models, which aren't able to predict the mean effect of the unsteady fluctuation of the separation and reattachment points correctly [START_REF] Jakirlic | Extended excerpt related to the test case: Flow over a periodical arrangement of 2D hills[END_REF]. Thus, the application of BMSA to this challenging flow configuration is a test of the limits of the framework utilizing linear eddy-viscosity models being calibrated for flat plate boundary layer flows.

Simulation setup

In order to mimic the periodicity of the hill-geometry, periodic boundary conditions were applied at the inlet and outlet and no-slip conditions at the walls. A volume forcing is applied to each cell, which maintains a bulk velocity of U b = 1.0 between the hill's crest and the upper surface. The functional form of the lower surface is defined according to the ERCOFTAC test case description 6 . In order to study the effect of the turbulence modelling and not additional modelling due to wall-functions the thickness of the smallest cells close to the wall were defined according to y + ≤ 1.0 and Low-Re turbulence models were used: Spalart-Allmaras, Launder-Sharma kε and Wilcox (2006) kω (for details of these models see [START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF]). Furthermore, a mesh convergence study was conducted for each model using baseline coefficients and the meshes in table 3 were used for the application of the BMSA method. Similar to the pipe flow case in Section 3 the results of the forward-simulations

m CFD M i , θ θ θ MAP i,k
with nominal values of the coefficients m CFD (M i , θ θ θ o i,k ) are used as an initial condition in order to reduce iteration counts. For S 1 4 the simulation using Wilcox (2006) kω didn't converge so that this scenario was excluded from the set leading to 13 used scenarios in total.

BMSA prediction for velocity field

The baseline simulations for the three models show the expected behaviour as reported in the literature for a larger Re-number [START_REF] Jakirlic | Extended excerpt related to the test case: Flow over a periodical arrangement of 2D hills[END_REF]: while both Spalart-Allmaras and the Wilcox (2006) kω over-predict the size of the recirculation zone, characterised by a zero velocity component close to the wall, the Launder-Sharma kε under-predicts this flow feature.

The BMSA approach using MAP estimates of the posterior probability distributions outputs both an expectation of the velocity E[∆ | S, z], for which ∆ = U x , and a confidence interval. For the latter two types were chosen: one obtained by sampling from the posterior predictive distribution and using the range from the 1.0 to the 99.0 percentile, i.e. 98% confidence interval, and another one by using the min/max range of the ppd. We have chosen both a larger confidence interval and the entire range of the distribution in order to take the expected low performance of the linear eddy-viscosity models for this test case into account. Especially the latter can be seen as an approach trying to envelop the true process by all models and for all scenarios. The smart-weighting technique to obtain the scenario probabilities as described in Section 2.2 was applied along the y/H-direction for each streamwise location for the 98% confidence interval. In that manner the predicitive similarity for each scenario per model is evaluated locally for the streamwise direction x.

The BMSA expectation E[∆ | S, z], as shown in Figure 4, over-predicts the recirculation zone, gives similar velocity profiles for x = 3.0 and 4.0 close to the lower surface, but is completely off for other locations. Especially, for x = 0.0 the expectation doesn't capture the local maximum of the velocity close to y/H = 1.0, i.e. at the hills crest, and over-predicts the velocity within the channel for every other location. Throughout the different locations the expectation shows large differences compared to the DNS for the upper part of the velocity profile.

The entire range of the posterior predictive distribution (ppd) and the 98% confidence interval show large differences over the entire domain as shown in Figure 4. The ppd captures the DNS data for x ≥ 3, but doesn't capture the local velocity maxima in the area of the free shear layer in the leeward region of x = 0 to 2.0. Interestingly for x ≤ 3.0 when the range of the ppd shrinks locally in y/H-direction also the DNS is still inside. However, this pattern is not the same for every y/H-position, e.g. at x = 2.0 the expectation and the DNS match for 0.5 ≤ y/H ≤ 1, but the error is high. For the relatively large confidence value 98% the intervals shrinks drastically, so that the DNS data is only captured in a small lower band y/H ≤ 1.0 for 3 ≤ x ≤ 7 and for the part of the velocity profile in the middle of the channel.

Interestingly, in the reattachment region between the two hills, where the case is similar to a flat plate, the BMSA approach in the present form with a confidence of 98% is able to capture the DNS data. However, the interval for the upper surface is always negligible but the DNS data doesn't match with the expectation.

The large effect of the min/max profiles, which determine the min/max range of the ppd, to capture the DNS data proves the aforementioned fact that the used linear eddy-viscosity models suffer from restrictions which inhibit the reproduction of the true flow state for this test case. However, posterior probability distributions of the closure coefficients of the models inferred on other flow scenarios more equal to the one here might lead to a different picture.

Industrial Test-case: Generic Dassault Falcon

The final test-case we consider is the generic Falcon jet -a business jet including tail, engine nacelles and winglets -at transonic conditions (Industrial challenge IC-03 of UMRIDA). The geometry is far more complex than any considered so far, with a resulting spectrum of primary and secondary flows. We provide this case as a demonstration of the applicability of our approach to problems of relevance in the aerospace industry -unfortunately without a reference solution. LES is impractically expensive here, and experimental data is not available. We attempt to justify the results of BMSA based on the expected response of the simulation to the closure modelling.

The CFD code used is Petrov-Galerkin finite-element RANS solver AETHER used within Dassault Aviation. As a preliminary step it was verified that the closure models implemented in AETHER were identical to the models for which closure coefficients were calibrated. The implementations of Spalart-Allmaras and Wilcox kω were found to be sufficiently similar to the calibrated models, and in addition are regarded by Dassault Aviation as suitable for this test-case. Of the 14 sets of closure coefficients computed with each model, only 5 cases were able to be successfully converged with kω, whereas all S-A cases converged without issue. BMSA requires a minimum of 2 models per scenario, and as such only 5 scenarios could be considered. Thus the spectrum of model results is substantially more limited than for both previous test-cases.

The resulting uncertainty in the pressure distribution at a cut on the wing (with p = 0) can be seen in Figures 5(a further. The very low variance can be attributed to two effects: (a) the real lack of sensitivity of the pressure distribution to the turbulent boundary layer in an essentially attached flow, for which an inviscid solution is already satisfactory, and (b) the limited range of models and coefficients used in the study. Indeed by eliminating exactly those coefficients that cause convergence problems (for kω), we are introducing additional bias into the sampling, likely biasing against extreme results. The only alternative -of including unconverged solutions -is even less attractive however. More insight can be gained by looking at force coefficients: Figure 6 shows total, pressure and viscous drag coefficients, all evaluated by surface integration, with mean and 95% confidence intervals, as a function of p. While pressure-and viscous-drag have approximately the same magnitude, the viscous part completely dominates the uncertainty -with C D,p varying at most 3 counts, and C D,v varying between 15 and 25 counts, depending on p. Thus the uncertainty in total drag is driven entirely by C D,v . Given the observed lack of uncertainty in C p -profiles, and the high sensitivity of friction coefficients to closure modelling, this is not surprising -though it does indicate that varying the closure model does not appear have a significant effect on separation behaviour in this case. Examining the relationship of uncertainty against p reveals which quantities are dominated by model-differences (where uncertainty is approximately constant with p), and for which quantities scenario differences are significant. Here C D,p belongs to the former class, and C D,v to the latter.

In summary, the limited number of turbulence models and range of coefficients lead to limited uncertainty in this case, for quantities of interest related to the pressure, which is essentially governed by the inviscid flow behaviour. Significant uncertainty is observed for viscous quantities such as friction drag. Resolving this in practical applications is the subject of ongoing work. 95% conf.

(c) viscous Fig. 6. Total, pressure and viscous drag coefficients for the full Falcon.

Conclusion

The BMSA approach based on MAP estimates of the posterior probability distributions of the closure coefficient has been applied to three different test cases.

For the flow in a straight pipe at Re = 44, 000 the results for the velocity profile and the pressure drop are in line with the validation data for this case, which shows in general that the approach can successfully be applied to other test cases outside of the set used for calibration.

The application of BMSA to the periodic hill test case at Re H = 5600 proves that this test case is challenging for the linear eddy-viscosity assumption and also shows the limits of the BMSA framework in the present form, provided that all models in the chosen set employ the Boussinesq hypothesis. However, the resulting ppd is able to capture most of the validation data successfully. A more suitable test case will be evaluated for BMSA in the future, such as the flow over a backward-facing step, for which the separation is forced by the geometry and not result of the simulation itself. Furthermore, the question of how to incorporate model-form error within a stochastic framework for UQ purposes is the topic of ongoing research.

Finally, the Falcon Jet test case (IC-03) revealed the open problem of how to deal with many non-converged solutions within the BMSA methodology and the successive artificial reduction of uncertainty by excluding these simulations from the set of simulations, which is also a topic for further research.

  yields the expression for the posterior distribution of θ θ θ p (θ θ θ |z, M, S) = p (z|θ θ θ , M, S) p (θ θ θ |M, S) p (z|M, S) ∝ p (z|θ θ θ , M, S) p (θ θ θ |M, S) .

2 Fig. 1 .

 21 Fig. 1. BMSA prediction using different smart weighting tuning parameter p.

  ) to 1(c) for different tun-

Fig. 2 .

 2 Fig. 2.

Fig. 3 .

 3 Fig. 3. BMSA prediction of pressure drop including confidence interval for several tuning parameter p of smart-weighting method.

2 Fig. 4 .

 24 Fig.[START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF]. BMSA for streamwise velocity component U x at several x-locations, E[∆ | S, z] (red), full posterior predictive distribution (ppd) range (shaded grey) and confidence interval of 98% (shaded blue) using p = 2. DNS data of Breuer et al.[START_REF] Breuer | Flow over periodic hills -Numerical and experimental study in a wide range of Reynolds numbers[END_REF] (black).

  Zoom at shock location.

Fig. 5 .

 5 Fig. 5. Pressure distribution at 30% span of the Falcon wing, y = 3848 mm.

  norm of difference between DNS data and BMSA expectation against p-value.

	0.00 0.02 0.04 0.06 0.08 0.10 0.16 0.14 0.12	1400 1300 2700 6300 1100 2100 2500 2400 3300 0141 1200 4400 2133 2134 p =0 p =1 p =2
	(b) Scenario-weights for several p-values. De-
	scription of scenarios in [3].

Table 2 .

 2 Pressure

	Model	∆ P/ρ [ m 2 s 2 ] f Darcy
	Wilcox(2006) k -ω	2.1427	0.02144
	Launder-Sharma k -ε 2.0478	0.02049
	Spalart-Allmaras	2.2074	0.02208
	Empirical value	2.1499	0.021499

drop over pipe length for Launder-Sharma kε, Wilcox (2006) kω and Spalart-Allmaras using nominal closure coefficient

[START_REF] Edeling | Predictive rans simulations via bayesian modelscenario averaging[END_REF]

. The corresponding Darcy friction coefficient f Darcy is calculated using the Darcy-Weisbach equation. The empirical value for f Darcy is based on the Colebrook equation for Re = 44, 000 and the corresponding ∆ P/ρ is calculated using the Darcy-Weisbach equation.

Table 3 .

 3 Mesh size per closure model for flow of periodic hills at Re H = 5600.

		mesh (n x × n y )
	Spalart-Allmaras	100 × 110
	Launder-Sharma k -ε 150 × 140
	Wilcox (2006) k -ω	100 × 110

The MAP estimates are available online at[START_REF] Edeling | MCMC traces of posterior distributions of closure coefficients for 6 turbulence models[END_REF]

The template-cases for each turbulence model are available on GitHub: https:// github.com/shmlzr/UQOpenFOAM

Underlying flow regime 3-30, 2D Periodic Hill Flow: http://qnet-ercoftac. cfms.org.uk

Acknowledgement

We thank Dassault Aviation and especially Gilbert Roge for the collaboration on the test case of the generic falcon jet.


Clément Larquet

Anh-Minh Nguyen

Estelle Glais

Lorenzo Paulatto

Capucine Sassoye

Mohamed Selmane

Pierre Lecante

Clement Maheu

Christophe Geantet

Luis Cardenas

Clément Maheu

Corinne Chanéac

Andrea Gauzzi

Clément Sanchez

Sophie Carenco
email: sophie.carenco@sorbonne-universite.fr

Bandgap Engineering from Cation Balance: the Case of Lanthanide Oxysulfide Nanoparticles

Context and objectives

The crystallographic structure of lanthanide and actinide oxysulfide RE 2 O 2 S (RE: rare earth) was solved for the first time during the 1940's. [START_REF] Pitha | The Preparation of Lanthanum Oxysulfide and Its Properties as a Base Material for Phosphors Stimulated by Infrared 1[END_REF][START_REF] Zachariasen | Crystal Chemical Studies of the 5f-Series of Elements. VII. The Crystal Structure of Ce2O2S, La2O2S and Pu2O2[END_REF] From these pioneering works, numerous studies were published on the preparation of these phases and their optical properties. In particular, doped Ln 2 O 2 S (Ln = La, Gd, Y) turned out to be excellent phosphors with good efficiencies for scintillators, [START_REF] Rossner | Phosphors for X-Ray Detectors in Computed Tomography[END_REF][START_REF] Michail | Light Emission Efficiency and Imaging Performance of Gd2O2S:Eu Powder Scintillator under x-Ray Radiography Conditions[END_REF] upconversion materials [START_REF] Auzel | Upconversion and Anti-Stokes Processes with f and d Ions in Solids[END_REF][START_REF] Hyppänen | Up-Conversion Luminescence Properties of Y2O2S:Yb3+,Er3+ Nanophosphors[END_REF][START_REF] Ajithkumar | Multimodal Bioimaging Using a Rare Earth Doped Gd2O2S:Yb/Er Phosphor with Upconversion Luminescence and Magnetic Resonance Properties[END_REF][START_REF] Li | La0.97RE0.01Yb0.02)2O2S Nanophosphors Converted from Layered Hydroxyl Sulfate and Investigation of Upconversion Photoluminescence (RE=Ho, Er)[END_REF] and bioimaging capabilities. [START_REF] Ajithkumar | Multimodal Bioimaging Using a Rare Earth Doped Gd2O2S:Yb/Er Phosphor with Upconversion Luminescence and Magnetic Resonance Properties[END_REF][START_REF] Rosticher | Gadolinium Oxysulfide Nanoprobes with Both Persistent Luminescent and Magnetic Properties for Multimodal Imaging[END_REF][START_REF] Lei | Controlled Synthesis of Lanthanide-Doped Gd2O2S Nanocrystals with Novel Excitation-Dependent Multicolor Emissions[END_REF] These crystalline phases are stable in air.

With the development of nanotechnologies during the past 20 years, rare-earth based materials, and especially Ln 2 O 2 S oxysulfide nanoparticles, emerged due to a growing interest for nanophosphors with good chemical and thermal stability, lower toxicity and better luminescence quantum efficiency than oxides. Among others, Ding and coworkers developed a versatile colloidal synthesis in organic solvents for Ln 2 O 2 S nanoparticles (Ln = Y, La, Pr, Nd, Sm-Tb). [START_REF] Ding | Sodium Doping Controlled Synthesis of Monodisperse Lanthanide Oxysulfide Ultrathin Nanoplates Guided by Density Functional Calculations[END_REF] The properties of doped Ln 2 O 2 S oxysulfides were widely explored, but only a few studies focused on the properties of these without doping. Ce 2 O 2 S is a peculiar case in the Ln 2 O 2 S series because Ce III can easily oxidize into Ce IV . This unique reactivity, related to the [Xe] 5d 0 4f 1 electronic configuration of Ce III , led to the description of mixed-valence structure of cerium oxysulfides, namely Ce 2 O 2.5 S, [START_REF] Sourisseau | Structure and Bondings in Cerium Oxysulfide Compounds I -Electronic, Infrared and Resonance Raman Spectra of Ce2.0O2.5S[END_REF][START_REF] Sourisseau | Structure and Bonding in Cerium Oxysulfide Compounds. II-Comparative Lattice Dynamics Calculations on Ce2O2S and Ce2.0O2.5S[END_REF] Ce 4 O 4 S 3 , [START_REF] Dugué | Etude Structurale Des Oxysulfures de Cérium(III) et Cérium(IV). I. Structure Cristalline de l'Oxysulfure de Cérium Ce4O4S3[END_REF] and Ce 6 O 6 S 4 . [START_REF] Dugué | Etude Structurale Des Oxysulfures de Cérium(III) et Cérium(IV). II. Structure Cristalline de l'Oxysulfure de Cérium Ce6O6S4[END_REF] In air and at room temperature, bulk Ce 2 O 2 S slowly oxidizes. [START_REF] Flahaut | Les Oxysulfures Me2O2S Des Éléments Du Groupe Des Terres Rares[END_REF][START_REF] Quezel | Proprietes Magnetiques Des Oxysulfures de Terres Rares[END_REF] We showed in a previous work that Ce 2 O 2 S nanoparticles can be obtained, provided that synthesis and storage are performed under inert atmosphere. 15 nm wide and 2 nm thick nanoplates were formed by colloidal synthesis in organic solvent. Because of their high specific surface, brief exposition to air or water swiftly oxidize both cerium and sulfur, triggering the degradation of the nanoparticles. [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] Bimetallic lanthanide oxysulfides (Gd , Ce) 2 O 2 S nanoparticles were also synthesized and their stability in air was found to decrease when the Ce:Gd ratio increase: up to 40 % cerium vs. gadolinium, the nanoparticles were structurally stable, although significant surface oxidation was observed on both sulfur and cerium. [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] Beyond its lack of stability in air, Ce 2 O 2 S is also peculiar because of its electronic and optical properties. Unlike the other lanthanide compounds, Ce 2 O 2 S is not whitish but deep brown. [START_REF] Sourisseau | Structure and Bondings in Cerium Oxysulfide Compounds I -Electronic, Infrared and Resonance Raman Spectra of Ce2.0O2.5S[END_REF][START_REF] Flahaut | Les Oxysulfures Me2O2S Des Éléments Du Groupe Des Terres Rares[END_REF] This absorption of the visible light suggests that Ce 2 O 2 S is a semiconductor whose bandgap is much lower than the others (< 2.0 eV for Ce 2 O 2 S, > 3.0 eV for the other lanthanide oxysulfides).

In this work, the optical and electronic properties of the bimetallic Gd 2(1-y) Ce 2y O 2 S nanoparticles are discussed. First, UV-diffuse reflectance spectra demonstrate the presence of two regimes in the Gd 2(1-y) Ce 2y O 2 S series: (i) a -doping‖ regime for low y values for which two distinct absorption phenomena are observed and (ii) a -solid solution‖ regime for y ≥ 5 % with a large absorption band which shifts with cerium concentration. In the latter case, the nanoparticles bandgap is progressively evolving with their composition. Second, an accurate description of the nanocrystal obtained with a typical colloidal synthesis is presented and accounts for the lack of sulfur in the final compound. Then, DFT calculations on Gd 2 O 2 S and Ce 2 O 2 S show that the reduced thickness of the nanoparticles explain their direct bandgap (while the bulk compound has an indirect bandgap). [START_REF] Mikami | Electronic Structure of Rare-Earth Sesquioxides and Oxysulfides[END_REF] This reconciles the colors of the powders with the bandgaps extracted from absorption spectra in UV-visible range. To the best of our knowledge, this work presents the first example of wide bandgap variation that is easily controlled from the cations ratio in the solid solution of an inorganic phase. It occurs without significant structure evolution.

Results and discussion

Optical properties -Ln/Ce substitution and bandgap evolution

Absorption of Gd 2(1-y) Ce 2y O 2 S nanoparticles

In a previous work, we reported the synthesis of Ce 2 O 2 S nanoparticles that are unstable when exposed to air. However, we showed that it is possible to obtain a large range of cerium-containing Gd 2(1-y) Ce 2y O 2 S nanoparticles that are stable in air. [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] We also showed that the ligands around Gd 2 O 2 S nanoparticles are oleates in bridging and chelating coordination modes. [START_REF] Larquet | Thermal Stability of Oleate-Stabilized Gd 2 O 2 S Nanoplates in Inert and Oxidizing Atmospheres[END_REF] Here, cerium loading is defined as y = n Ce / (n Ce + n Gd ) where n refers to the initial metallic precursor quantity in mol. The Gd/Ce ratios of the precursors introduced at the beginning of the reaction and measured by energy dispersive X-ray spectroscopy (EDX) on the isolated product were consistent (Table S1 and Figure S0), meaning that the ratio in the final powder is similar to the introduced one. The nanoparticles were slightly larger when they were richer in cerium (Table S2) but they presented a similar thickness of three lattice units (Figure S3). More Ce III in the structure in place of Gd III strongly affects light absorption of the Gd 2(1-y) Ce 2y O 2 S nanoparticles powders, from white (y = 0 %) to brown (y = 70 -90 %). The Kubelka-Munk function (F(R) = (1 -R) 2 / 2R where R is the reflectance, see Figure S1) is linearly related to the extinction coefficient of the material and thus to its absorption properties. [START_REF] Kubelka | Ein Beitrag Zur Optik Der Farbanstriche[END_REF] The F(R) spectra of the Gd 2(1-y) Ce 2y O 2 S nanoparticles are plotted in Figure 1A.

The shoulders at ca 300 nm are due to organic ligands. Note that the range of wavelengths where light is absorbed increases with Ce fraction. In agreement with the change of colors observed in the powders (Figure S1), Gd 2 O 2 S absorbs only in the UV region while Ce-rich nanoparticles absorb a large fraction of the visible light (brown powders).

Although this evolution is gradual, two regimes can be distinguished (Figure 1). In the doping regime (0 % < y ≤ 2 %), two inflexion points are observed with a smooth absorption increase before 500 nm and a sharp increase around 270 nm, close to the absorption step of Gd 2 O 2 S. In the solid solution regime (y ≥ 5 %), only the smooth increase of absorption is observed and this feature shifts towards higher wavelengths.

A B

Bandgap of Gd 2(1-y) Ce 2y O 2 S nanoparticles

Thanks to the work by Kubelka et al. [START_REF] Kubelka | Ein Beitrag Zur Optik Der Farbanstriche[END_REF] and later analyzed by Tauc et al., [START_REF] Tauc | Optical Properties and Electronic Structure of Amorphous Germanium[END_REF] it is possible to extract the value of the bandgap from the reflectance spectra using the F(R) function (Figure 1B). This model is parameterized by a power n which depends on the type of bandgap. Then the extrapolation of the linear portion of the (F(R)hν) n vs. hν curves enables to estimate the bandgap.

In a first report, Mikami et al. suggested that the bandgap of Ln 2 O 2 S compounds is indirect. [START_REF] Mikami | Electronic Structure of Rare-Earth Sesquioxides and Oxysulfides[END_REF][START_REF] Mikami | First-Principles Band-Structure Calculation of Yttrium Oxysulfide[END_REF] In this case, n = ½. However, in recent reports on Ln 2 O 2 S compounds, it was argued that the bandgap is direct, which corresponds to n = 2. [START_REF] Som | Spectral and Trapping Parameters of Eu 3+ in Gd 2 O 2 S Nanophosphor[END_REF][START_REF] Bagheri | Synthesis and Characterization of Physical Properties of Gd2O2S:Pr3+ Semi-Nanoflower Phosphor[END_REF] Here, the absorption of the powders are in agreement with the latter picture, as the assumption n = ½ yielded unphysical values of 1.2 -4.1 eV (color mismatch) (see Figure S8). Thus, hereafter we assume a picture of direct bandgap for the Gd 2(1-y) Ce 2y O 2 S nanoparticles.

The analysis of the F(R) spectra of Gd 2(1-y) Ce 2y O 2 S nanoparticles is presented on Figure 1B.

For y ≥ 5% (solid solution regime), the linear zone on the (F(R)hν) 2 curve gives a bandgap below 3 eV. For the two compounds in the range 0 % < y < 5 % (doping regime), two values were extracted from the curves: (i) one related to cerium insertion (low energy region 2.9 -3.1 eV) and (ii) the other close to the bandgap of Gd 2 O 2 S nanoparticles (high energy region 4.3 -4.7 eV). For these two compounds, the aforementioned does not provide a unique value. However, the strong effect of cerium on the optical properties is evident from our data.

Figure 1B inset shows the bandgaps extracted from the absorption spectra. The dark black line in the bottom indicates a linear trend for y > 2 % (solid solution regime), although the reliability of the measurements obtained for y = 80 % and y = 90 % is questionable because of the limited stability of Ce-rich nanoparticles in air. [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] The shift of bandgap between Gd 2 O 2 S and Gd 0.6 Ce 1.4 O 2 S (y = 70 %) is remarkably large (2.6 eV), considering that the crystal structure is preserved. As shown in a previous work, the only structural change between the two samples consists in a 1.5 % increase of lattice parameter a as expected from Vegard's law. [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] In the solid solution regime, the bandgap of the nanocrystals is tunable between 3.0 eV and 2.1 eV by adjusting the nanoparticles composition. To some extent, in the high energy region, the bandgap is also tunable in the doping regime between 4.7 eV and 4.4 eV (inset of Figure 1B, green points). In the doping regime, a relatively large energy gap separates the two energy regions (> 1 eV) corresponding to the two distinct absorption phenomena. We will see in the final sections that replacing Gd III by Ce III disrupts the electronic properties of Gd 2 O 2 S nanoparticles in a large extent. We can suspect that the mixed valence of cerium [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] can be at least partially responsible for the strong absorption in the visible range but so far, we have no experimental evidence to support this hypothesis.

Comparison with other class of compounds

Bandgap engineering enabling large changes of absorption spectra has been discussed for many materials. In general, these changes typically require significant structural changes. In composite semi-conductors, the gap energy variation is related to the interatomic distance.

That is why bandgap shifts can be observed as a consequence of anionic or cationic substitution. The more the structure is modified (variation of the interatomic distances), the higher the expected shift. For example, the anionic substitution was investigated in hybrid lead-halide perovskites, extensively studied for their promising properties for solar cells applications. [START_REF] Jeon | Compositional Engineering of Perovskite Materials for High-Performance Solar Cells[END_REF] The substitution of bromine by iodine in (CHA) 2 PbBr 4(1-x) I 4x (CHA = cyclohexylammonium) is accompanied by a bandgap decrease of ~ 0.7 eV between (CHA) 2 PbBr 4 and (CHA) 2 PbI 4 but with significant structural modifications. [START_REF] Ye | Bandgap Engineering of Lead-Halide Perovskite-Type Ferroelectrics[END_REF] Solid solutions of CdS x Se 1-x are also known to exhibit large absorption shifts. [START_REF] Zhou | Tailoring Band Structure of Ternary CdS x Se 1-x Quantum Dots for Highly Efficient Sensitized Solar Cells[END_REF] Copper zinc tin sulfide (CZST) [START_REF] Malerba | CZTS Stoichiometry Effects on the Band Gap Energy[END_REF] and Copper indium gallium selenide (CIGS) [START_REF] Tinoco | Phase Diagram and Optical Energy Gaps for CuInyGa1-ySe2 Alloys[END_REF] thin films also exhibit a bandgap that is related to the cations stoichiometry and may vary in a broad range. Similar effect was observed on gallium indium oxide nanoparticles. [START_REF] Farvid | Colloidal Gallium Indium Oxide Nanocrystals: A Multifunctional Light-Emitting Phosphor Broadly Tunable by Alloy Composition[END_REF] In the case of oxysulfides, Sm 2 Ti 2 O 7 and Sm 2 Ti 2 O 5 S 2 show a bandgap variation of ΔE g ≈ 1.6 eV but no solid solution can be formed. [START_REF] Ishikawa | Oxysulfides Sm2Ti2S2O5 as Stable Photocatalysts for Water Oxidation and Reduction under Visible-Light Irradiation[END_REF] The substitution of two cations also significantly changes the bandgap: 3.8 eV for ZnS, 2.5 eV for CdS and 2.1 eV for α-HgS. [START_REF] Strehlow | Compilation of Energy Band Gaps in Elemental and Binary Compound Semiconductors and Insulators[END_REF] In addition to strong variations of lattice parameter in isostructural solid solutions, the substitution of zinc by cadmium in Zn 1-x Cd x S drives a cubic blende ZnS to hexagonal CdS structural transition. [START_REF] Li | Zn 1-XCdxS Solid Solutions with Controlled Bandgap and Enhanced Visible-Light Photocatalytic H2-Production Activity[END_REF] In conclusion, the present system is rare because, within well-defined and small nanoplates, the cationic substitution is responsible for an important absorption shift with neither structural nor morphological modifications. In view of this promising result, we investigated whether Ce-substitution constitutes an exception or other substitutions would give similar results.

Comparison with Pr and Eu cation balance

We have chosen Ce 2 O 2 S among all Ln 2 O 2 S compounds because it is the only low bandgap semiconductor. The insertion of other lanthanides is expected to have less effect on the absorption. Also, the question arises whether the insertion of Ce in other Ln 2 O 2 S compounds may lead to absorption shifts similar to those observed in Gd 2(1-y) Ce 2y O 2 S nanoparticles.

Other bimetallic nanoparticles were synthesized to investigate this effect. We selected Pr and Eu, the former because of its cation size (1.20 Å) close to Ce (1.21 Å), the latter because of its two oxidation states (here, +II and +III). Using Ce(acac) obtained. The F(R) spectra of these bimetallic nanoparticles are presented on Figure 2.

For Gd 1.8 Ln 0.2 O 2 S compounds (Ln = Ce, Pr, Eu), the most significant effect on the absorption is obtained with cerium (Figure 2A). For Gd 1.8 Pr 0.2 O 2 S, two absorption edges are observed at 280 nm close to Gd 2 O 2 S edge and 320 nm closer to the absorption edge of Pr 2 O 2 S (around 350 nm, Figure 2B andC). In the case of europium, the edge is shifted to lower wavelengths (below 250 nm) and a broad peak on the curve grows around 340 nm (Figure 2B). Cerium appeared to be unique among the lanthanide investigated. To explain the bandgap nature and its huge sensitivity of the optical properties to the Ce concentration, we underwent an in-depth structural and modeling study.

Gd 2(1-y) Ce 2y O 2 S nanoplates: a sulfur-deficient structure

Nanoparticles morphology

The Gd 2(1-y) Ce 2y O 2 S nanoparticles crystallize in the P-3m1 crystal group. The structure is hexagonal and alternates layers of [Ln 2 O 2 ] 2+ and [S] 2-. [START_REF] Zachariasen | Crystal Chemical Studies of the 5 f -Series of Elements. X. Sulfides and Oxysulfides[END_REF] In this structure, Gd and Ce form a solid solution. [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] In the present case, as in other authors works, [START_REF] Lei | Controlled Synthesis of Lanthanide-Doped Gd2O2S Nanocrystals with Novel Excitation-Dependent Multicolor Emissions[END_REF][START_REF] Ding | Sodium Doping Controlled Synthesis of Monodisperse Lanthanide Oxysulfide Ultrathin Nanoplates Guided by Density Functional Calculations[END_REF] the formed nanoparticles possess a hexagonal anisotropic shape with exacerbated {001} facets.

Transmission Electron Microscopy (TEM) micrographs show nanoparticles either observed lying on this {001} facet (Figure 3A, Figure S3), or exposing their {100} facet (Figure 3B, Figure S3). In both cases, the nanoparticles are packed in fairly regular arrangements that are a consequence of specific interactions between the surface ligands. Gd 2(1-y) Ce 2y O 2 S nanoparticles are 10 to 15 nm wide when y > 5 % and 7 to 10 nm wide otherwise. This is in agreement with the calculated Scherrer's crystal domains obtained from the XRD patterns (see Figure S4 and Table S2), meaning that the nanoparticles are single-crystals. (omitted on the scheme). The number of atoms depicted is in agreement with the nanoparticles average size from TEM micrographs (Figure S3).

Sulfur-defective nanoparticles

First, energy dispersive X-ray spectroscopy (EDS) performed on Gd 2(1-y) Ce 2y O 2 S nanoparticles showed a significant lack of sulfur (Table S1), as reported by us and other groups performing a similar synthesis. particles and completely absent from the {001} facets. In this particular case, the nanoparticle terminates with [Gd 2 O 2 ] 2+ layers on {001} facets and surface gadolinium cations on these layers can be coordinated by oleate ligands. The theoretical S:Gd ratio is not anymore 0.50, as in the bulk compound, but only 0.37. Removing all surface sulfur (not only on {001} facets but also on {100} facets) leads to a lower S:Gd ratio of 0.30 (Gd 2 O 2 S 0.6 nanoparticles, Table S3). The corresponding nanoparticle is represented on Figure 3C.

To balance the charge caused by the lack of sulfur, the presence of organic ligands at the surface of the nanoparticles is necessary. These ligands come from the colloidal synthesis, which requires oleylamine and oleic acid. The presence of oleate moieties in the powder was confirmed by Fourier-transform infrared spectroscopy (FTIR, Figure S5), which showed the typical bands of carboxylate C=O double bonds at 1500 cm -1 and 1383 cm -1 .

Thermogravimetric analysis coupled with mass spectrometry (TGA-MS) indicated around 30 wt% of organic compounds in synthesized Gd 2 O 2 S and confirmed the absence of oleylamine (data not shown). We thus propose that only oleate groups cap the [Ln 2 O 2 ] 2+ layers on each side of the particles, id est on both the {001} and the {100} facets.

Although these models are only ideal cases, they all concord to support that the nanoplates are sulfur-defective and are terminated by [Ln 2 O 2 ] 2+ layers capped with oleates. In order to confirm this, we performed a structural analysis sensitive to defects: pair distribution function analysis. This technique was also used to discuss the localization of sodium cations, which are generally believed to substitute Gd III cations in Ln 2 O 2 S structures. [START_REF] Ding | Sodium Doping Controlled Synthesis of Monodisperse Lanthanide Oxysulfide Ultrathin Nanoplates Guided by Density Functional Calculations[END_REF][START_REF] Zhang | Experimental and Theoretical Studies on the Controlled Synthesis of Alkali-Metal-Doped Rare-Earth Oxysulfide Nanocrystals[END_REF] Pair distribution function analysis X-ray diffraction pattern of Gd 2 O 2 S nanoparticles was measured on a large angular range with a low-background home-made sample-holder. The pair distribution function (PDF) was calculated from total X-Ray scattering as described in the experimental section. It was compared with simulated patterns. As seen on Figure 4, the calculated and experimental PDF are mostly similar (the relative intensity of the second peak at r = 2.86 Å will be discussed in the following paragraph), confirming the validity of the published structure Gd 2 O 2 S as a model for the nanoparticles.

The positions of the experimental PDF peaks match with the positions of the calculated ones, confirming that atomic distances within the particle are similar to these of the bulk structure. This shows that the Na + and SO 4 2-ions detected by EDX, XPS and FTIR [START_REF] Larquet | Synthesis of Ce2O2S and Gd2(1-y)Ce2yO2S Nanoparticles and Reactivity from in Situ X-Ray Absorption Spectroscopy and X-Ray Photoelectron Spectroscopy[END_REF] are not inserted inside the core cells. Despite their close ionic radius (Gd III (VII) = 1,14 Å, Ce III (VII) = 1,21 Å, Na I (VII) = 1,26 Å), [START_REF] Shannon | Revised Effective Ionic Radii and Systematic Studies of Interatomic Distances in Halides and Chalcogenides[END_REF] sodium ions are not in the final crystalline structure.

The 2.86 Å peak, representing the Gd-S pair, appears experimentally with a smaller intensity than the calculated one for Gd 2 O 2 S. To account for this difference, a series of models with decreasing S -II occupancy was simulated (Figure 4B). The S/(Gd 2 O 2 ) ratio mostly affects the relative intensity of the Gd-S peak, with best match around half S atoms from ideal structure.

This value is in the same range than those expected from EDS measurements on several aggregates of particles.

One could expect that such low statistical occupancy for S -II (0.5 occupancy for each S -II ) would lead to clear cell parameter contraction. This is clearly not the case here: whenever allowed to refine, cell parameters variations were balanced and negligible. As discussed previously, the anisotropic shape and size of the particle already lead to a sulfur deficient stoichiometry of Gd 2 O 2 S 0.6 . The experimental 0.5 S:Gd ratio would then correspond to an 85 % occupancy of the available sulfur position in the particle. Again, the absence of cell contraction means that the main core of the particle keeps a Gd [START_REF] Mikami | Electronic Structure of Rare-Earth Sesquioxides and Oxysulfides[END_REF][START_REF] Mikami | First-Principles Band-Structure Calculation of Yttrium Oxysulfide[END_REF] However, we suggest that the calculations accuracy would be significantly improved by including in valence the 4f electrons of the lanthanide metals that are of great importance for their electronic properties.

Furthermore, to be as close as possible to the nanoplates represented on Figure 4C (S:Gd ratio), we also performed DFT on a specific -slab‖ geometry. Here, capping agents that stabilize the surfaces of the nanoplate were not considered because we supposed that only the crystal (core of the nanoplate) was responsible for the optical properties. Yet, the nature of the ligands may slightly influence the band structure.

The calculations were only focused on three layers of the studied Gd 2 O 2 S or Ce 2 O 2 S oxysulfide, that can be stabilized by hydrogen atoms instead of sulfur atoms when needed (sulfur-free {001} surfaces). This model structure mimics the absence of sulfur in the external layers (replaced by oleate ligands in reality) and avoids dangling electronic states that make the system metallic (see in SI the details on the DFT+U calculations). Ln 2 O 2 S nanoparticles. Further analysis using photon-in photon-out spectroscopy such as Resonant inelastic X-ray scattering (RIXS) may further strengthen this argument.

Bandgap nature of (Gd,Ce) 2 O 2 S nanoparticles in relation with the substitution ratio

The cationic substitution between Gd and Ce drastically modifies the electronic properties of the oxysulfide nanoparticles. We performed additional spectroscopic analyses to determine the electronic changes caused by Gd/Ce substitution. Moreover, the band shapes indicate that the Ce III 5d emitting levels are the same for y = 1 -20 %. However, it is hard to distinguish the transitions 5d→ 2 F 5/2 (Ce III ground state) from 5d→ 2 F 7/2 (spin-orbit excited term), so we arbitrarily chose the first case to illustrate our system in Figure 6C. [START_REF] Fu | Compositional Dependence of Scintillation Yield of Glasses with High Gd2O3 Concentrations[END_REF][START_REF] Wang | High Light Yield Ce3+-Doped Dense Scintillating Glasses[END_REF] In addition, ultraviolet photoelectron spectroscopy (UPS) was performed on the samples with y = 0 and 30 %, in order to see how the Ce substitution influence the position of the valence band. To conduct UPS on powders, the dropping method reported lately is used. [START_REF] Maheu | UPS and UV Spectroscopies Combined to Position Energy Levels of TiO 2 Anatase and Rutile Nanopowders[END_REF] The latter requires a solvent which may induce an interfacial dipole between the powder and the ITO substrate and so, a global shift of the spectra. As detailed on Figure S6, we aligned the samples at their Fermi levels. It allows measuring relative position of the valence maximum.

Between y = 0 % and y = 30 %, a significant upward shift of 0.5 eV of the valence band is observed (Figure 6B). This is consistent with the reduced bandgap of the material. The bandgap is however decreased by a larger value (2.2 eV) compared with these of Gd 2 O 2 S. We thus propose that the bandgap engineering caused by Gd/Ce substitution is majorly due to the variation of the conduction band.

A careful background subtraction allows commenting the experimental density of states determined with UPS. Figure 6B inset evidenced that the valence band states are modified by Gd/Ce substitution. For y = 30 %, the intensities of the states between 8 and 10 eV increased while the intensities of the states between 3 and 5 eV decreased. Those states are assigned to the Kohn-Sham bands between -1.5 and -3.0 eV (relatively to the valence band maximum) plotted on Figure 5 and reproduced in the inset. On Figure 6B, Ce states responsible for the valence band shift are observed between 1.5 and 3.5 eV (black arrow). Accordingly, the band structure of GdCeO 2 S phase (y = 50 %) confirms the increased density of these states (Figure S7). Altogether, we describe the Gd 2(1-y) Ce 2y O 2 S nanoparticles as a multi-regime system (Figure 6C). Gd 2 O 2 S nanoparticles exhibit a wide gap semiconductor behavior. With low fractions of cerium (in the -doping regime‖, y ≤ 10 %), the behavior is dual: cerium plays the role of (i) a classical doping ion (constant PL properties) with a small deviation from the high bandgap of Gd 2 O 2 S and (ii) a nucleus of Ce 2 O 2 S with its own absorption properties (low bandgap). The latter is the origin of the absorption properties for higher cerium concentrations (-bimetallic regime‖, y > 10 %) for which the electronic properties of the nanoparticles correspond to a solid solution (Gd,Ce) 2 O 2 S.

Conclusion

(Gd,Ce) 2 O 2 S nanoparticles exhibit colors from white to brown over the whole composition range. Bandgap engineering is thus possible through cation exchange and values varying from 4.7 eV in Gd 2 O 2 S to 2.1 eV in Gd 0.6 Ce 1.4 O 2 S can be obtained without major structural modifications. Surprisingly, the bandgap of the nanoparticles is direct, as validated by density functional theory. This is a result of the limited thickness of the nanocrystals. The influence of the nanoplates width is yet to be explored. Bandgap engineering on (Gd,Ce) 2 O 2 S nanoparticles may suggest applications in solar cells. The fine control of bandgap over a wide range, solely triggered by the cation ratio, is rare in the field of bandgap engineering, especially when considering nanoparticles. The simple idea of using suitable chemical substitutions producing minor structural changes should inspire the design of other bimetallic nanoscaled compounds, which may be promising especially in the field of visible light photocatalysis. (32.5 mmol) were placed in a 100 mL three-neck flask at room temperature. The brown solution was heated to 120 °C under vacuum for 20 min to remove water and other impurities with low boiling points. The mixture was then heated to 310 °C and stirred at this temperature for 30 min under purified N 2 . The transparent solution gradually became turbid starting from 280 °C. Then the mixture was left to cool to room temperature under purified N 2 . The nanoparticles were isolated using ethanol and washed at least three times using a THF/ethanol (1/5 in volume) mixture to remove the remaining reagents and organic matter. 40 to 90 mg of dried Gd 2(1-y) Ce 2y O 2 S particles was thus obtained depending on the initial cerium loading.

Experimental section

Synthesis of Gd

It should be noted that a significant decrease in the reaction yield (measured after several washings of the nanoparticles in air) is observed with increasing cerium concentration in the nanoparticles: starting from 0.5 mmol of lanthanide precursor, a yield of around 100 % vs. Gd (95 mg of powder) was obtained for y = 0 % (Gd 2 O 2 S) while a yield of 48 % (45 mg) was obtained for y = 70 % (Gd 0.6 Ce 1.4 O 2 S). Pr 1.8 Ce 0.2 O 2 S and Eu 1.8 Ce 0.2 O 2 S were synthesized with the same protocol and are stable in air.

X-ray Powder Diffraction (XRD)

The different X-ray diffraction patterns of dry powders were measured on a Bruker D8 ADVANCE diffractometer using the Cu doublet Kα 1 /Kα 2 radiation at 1.5406 Å. Typical diffractograms were collected with steps of 0.05 ° and a scanning speed of 5 s/point with a 1D Lynxeye detector (Bruker). The backgrounds of the patterns were subtracted using the EVA software. When low noise Si monocristals sample holders were used, the angular position 2θ was corrected by adjusting the sample height (correction around 0.5 to 1 mm).

Transmission Electron Microscopy (TEM) and High Resolution TEM (HRTEM)

A drop of a diluted solution of Gd 2(1-y) Ce 2y O 2 S dispersed in pure THF was allowed to dry on an amorphous carbon coated copper grid. TEM images were collected with a TWIN 120 (TECNAI SPIRIT) instrument at 120 kV and with a Gatan Orius 1000 4k x 4k camera.

Energy Dispersive X-ray Spectroscopy (EDS)

A small amount of powder was deposited on carbon adhesive tape on a scanning electron microscope sample holder. EDX analyses were performed on a SEM HITACHI S-3400N instrument at 10 kV. Titanium was chosen as the reference and analyses were performed on at least three different zones on the sample.

Fourier Transform Infrared Spectroscopy (FTIR)

Infrared spectra were collected on a Spectrum 400 (PERKINELMER) spectrometer. The dry sample (1-3 mg) was deposited on an attenuated total reflectance (ATR) crystal.

Transmittance was measured between 4000 cm -1 and 550 cm -1 with steps of 0.5 cm -1 .

UV diffuse reflectance

Dry powders were loaded in the circular sample holder to make a uniform layer of solid.

Powders were used, instead of solutions, first because the stability of the colloidal solutions is limited, second, in order to be able to perform other characterizations as well (XRD, UPS).

The UV-visible diffuse reflectance spectra were measured using an integration sphere of a Cary-WinUV 5000 spectrophotometer (AGILENT) between 250 nm and 800 nm at 1 nm.s -1 and corrected with a sample of BaSO 4 as reference. The F(R) spectra were deducted using the formula ( ) ( ) , where R is the reflectance. Depending on the bandgap (direct or indirect), linear zones could be extrapolated to 0 on the abscissa axis for the Tauc curve (F(R)hν) n = f(hν) where n = 2 for a direct bandgap and n = ½ for an indirect bandgap.

Pair distribution function analysis (PDF)

The reduced atomic PDF, G(r), gives, from a structural model, the probability of finding a pair of atoms at a distance of r as follows:

( ) [ ( ) ]
where ρ(r) is the total atomic pair-density and ρ 0 is the average atomic number density. As defined, peaks in the G(r) function point to specific distances separating a pair of atoms.

G(r) can also be extracted from an experimental XRD pattern. Indeed, G(r) is the sine Fourier transform of the total reduced total scattering function

F(Q) = Q[S(Q) -1]
, S(Q) being the total scattering structure function:

( ) ( ) ∫ [ ( ) ] ( )
where Q is the magnitude of the scattering vector (Q = 4π sin(θ) / λ), 2θ being the angle between the incoming and outgoing X-ray, λ is the X-ray wavelength, and S(Q) is the total scattering structure function which is the corrected and normalized experimental intensity.

The experimental data were recorded from a standard laboratory D8 ADVANCE diffractometer (Bruker) equipped with a Mo anode (doublet K α1 /K α2 radiation at 0.7093 Å), in the reflection mode, using a homemade sample holder designed to minimize its own signal.

The large amount (100 mg) of analyzed powder allows processing the signal with no further treatment such as sample holder background subtraction. Three patterns were recorded with a 1D Lynxeye detector (Bruker) with increasing step time (step size = 0.05°, 2θ range; step time: 8-50 °; 25 s, 50-90 °; 50 s, 90-150 °; 100 s), converted in count per second and merged to obtain the full XRD pattern. The XRD data were then processed using the program PDFgetX3, from Q = 1.8 to 15 Å -1 to obtain G(r) exp . [START_REF] Juhás | PDFgetX3 : A Rapid and Highly Automatable Program for Processing Powder Diffraction Data into Total Scattering Pair Distribution Functions[END_REF] On the other hand, the PDF data were modeled from the Gd 2 O 2 S structure (PDF 04-008-6593)

and compared with the experimental one using PDFGUI. [START_REF] Farrow | PDFfit2 and PDFgui: Computer Programs for Studying Nanostructure in Crystals[END_REF] Calculations in real space were performed for r from 1 to 20 Å, way below the width of the hexagonal crystalline nanoplates, to reduce the effect of anisotropy. A mean size particle factor of 25 Å was introduced for best attenuation in the high r region. [START_REF] Farrow | Relationship between the Atomic Pair Distribution Function and Small-Angle Scattering: Implications for Modeling of Nanoparticles[END_REF][START_REF] Gilbert | Finite Size Effects on the Real-Space Pair Distribution Function of Nanoparticles[END_REF] Real space calculated PDF are here compared without any further structure refinement (neither cell parameters nor atomic positionsthermal agitation were only adjusted to 0.01),

given the complexity of correctly evaluating PDF from finite anisotropic nanoparticles.

DFT calculations

In the ab-initio simulation of rare-earth oxysulfides we have one main difficulty to tackle: how to treat the very localized, but not completely chemically inert, f orbitals of the rare earth.

The 4f orbital contains 1 electron in the conventional isolated-atom structure of cerium, 7 electrons in gadolinium. The 4f orbitals are localized inside a sphere of around 1 Å, centered on the atomic nucleus, but it is only partially filled and its binding energy is comparable with that of the chemically more active 5d and 6s orbitals. Even if these orbitals participate little in chemical bonds, their inclusion in the description of the valence structure is important to have an accurate band structure, especially in the lower energy valence bands.

The very strong self-interaction error, makes a simple DFT calculation impossible, to compensate the error, we have used the DFT+U technique. [START_REF] Anisimov | Theory and Mott Insulators: Hubbard U Instead of Stoner I[END_REF][START_REF] Cococcioni | Linear Response Approach to the Calculation of the Effective Interaction Parameters in the LDA+U Method[END_REF] All calculations were done with the Quantum-ESPRESSO code [START_REF] Giannozzi | QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulations of Materials[END_REF][START_REF] Giannozzi | Advanced Capabilities for Materials Modelling with Quantum ESPRESSO[END_REF] using the PBE generalized gradient approximation. [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] We used ultrasoft pseudopotentials with valence-bound 4f orbitals from the pslibrary. [START_REF] Dal Corso | Pseudopotentials Periodic Table: From H to Pu[END_REF] These pseudopotentials trade an excellent transferability with an extremely high plane-wave cutoff of 200 Ry. Details on the method and results, including the relaxed geometry parameters, are in supplementary information.

Photoluminescence

Excitation and emission spectra were measured by Cary Eclipse Fluorescence Spectrophotometer. Powder samples were stuck on glass slide and excited at 205 nm for emission spectra. For excitation spectra the emission wavelength was fixed at 486 nm. The excitation and emission spectra were respectively normalized by their intensity at 205 nm and 486 nm.

Ultraviolet Photoelectron Spectroscopy

The methodology developed in Maheu et al. was applied to the oxysulfides powders. [START_REF] Maheu | UPS and UV Spectroscopies Combined to Position Energy Levels of TiO 2 Anatase and Rutile Nanopowders[END_REF] They were ultrasonically dispersed in heptane (4 g/L) and several drops were deposited on a conducting ITO substrate (R S = 30 -60 Ω/sq, purchased from Merck). Solvent was removed in the ultra-high vacuum chamber of the apparatus during 48 h. UPS analyses were conducted in an integrated ultrahigh vacuum system, connected to an A Axis Ultra DLD spectrometer (Kratos Analytical). UPS spectra were obtained with the first resonance line of a He lamp (hν = 21.2 eV). To separate the secondary edges of sample and analyzer, a negative potential about 9.1 eV was applied to the sample. This bias was optimized adjusting the work function of a clean Ag foil. The cut-off energy is determined using the first derivative of the UP spectrum and spectra are aligned at the same cut-off, as argued in the supplementary information. At low binding energy, the contribution of the secondary electrons is strong, background subtraction is therefore needed. The background developed by Li et al. is used and I signal is plotted (I signal = I exp -I background ). [START_REF] Li | Inelastic Electron Background Function for Ultraviolet Photoelectron Spectra[END_REF] Spectra were normalized to compare significant variations of the density of states. We determined the valence band maximum (E VB ) with a two times derivation of the spectra. The inflexion point is the E VB .

The intrinsic resolution of the apparatus is 0.2 eV, the error on the E VB is therefore: ΔE VB = |Δ UPS + ε| with ε the error made with the graphical method (ε ≈ 0.2 eV). Note that similar results were obtained using the intersection point between the zero-intensity constant background and a linear extrapolation of the valence band decay. However, as discussed in Maheu et al., [START_REF] Maheu | UPS and UV Spectroscopies Combined to Position Energy Levels of TiO 2 Anatase and Rutile Nanopowders[END_REF] the second derivative method evidence the first electrons that are photogenerated from the top of the valence band.

Table of Content Entry

Figure 1 :

 1 Figure 1: Optical properties of (Gd,Ce) 2 O 2 S nanoparticles. (A) F(R) spectra obtained from

Figure 2 :

 2 Figure 2: Optical properties of (Ln 1 ,Ln 2 ) 2 O 2 S nanoparticles. (A) Sample colors and bandgap

Figure 3 :

 3 Figure 3: Morphology of the nanoparticles. TEM micrographs of frontal (A) and lateral (B)

  [START_REF] Lei | Controlled Synthesis of Lanthanide-Doped Gd2O2S Nanocrystals with Novel Excitation-Dependent Multicolor Emissions[END_REF][START_REF] Ding | Sodium Doping Controlled Synthesis of Monodisperse Lanthanide Oxysulfide Ultrathin Nanoplates Guided by Density Functional Calculations[END_REF] Only 40 to 50 % (S:Gd = 0.2 -0.25) of the expected sulfur is found in the powders of Gd 2(1-y) Ce 2y O 2 S nanocrystals. However, we discussed in a previous study the surface oxidation of sulfur into sulfite and sulfate groups for this series of bimetallic Gd-Ce nanoparticles when they are exposed to air. We can thus imagine an extreme case in which sulfur is only exposed through the {100} facets of the

Figure 4 :

 4 Figure 4: Structure of the nanoparticles. (A) Pair distribution function G(r) obtained from

Figure 5 :

 5 Figure 5: Band structures of Ln 2 O 2 S (Ln = Gd, Ce) in the bulk (top) or for the slab geometry

Figure

  Figure 6A displays the photoluminescence (PL) excitation and emission spectra of

Figure 6 :

 6 Figure 6: Bandgap engineering of the nanoparticles. (A) Photoluminescence excitation and

  3 •xH 2 O, Pr(acac) 3 •xH 2 O, Eu(acetate) 3 •xH 2 O and Gd(acac) 3 •xH 2 O as lanthanide precursors, nanoparticles of Pr 2 O 2 S (light green powder), Eu 2 O 2 S (white powder), Gd 1.8 Pr 0.2 O 2 S (white powder), Gd 1.8 Eu 0.2 O 2 S (white powder), Pr 1.8 Ce 0.2 O 2 S (orange powder), and Eu 1.8 Ce 0.2 O 2 S (yellow powder) were

  Inserting cerium in Pr 2 O 2 S and Eu 2 O 2 S leads to similar effects on light absorption than Gd 2 O 2 S (Ln 1.8 Ce 0.2 O 2 S compounds, Figure2B), even if the bandgap shift caused by Ce III in gadolinium oxysulfide remains by far the largest. The nanoparticles absorb higher wavelengths as the bandgap decreases. Moreover, the absorption peaks specific to transitions between 4f energy levels of the monometallic compounds are no longer present in the bimetallic compounds whose absorption curves are smooth.These results demonstrate that the effect of cerium on the bandgap is not specific to gadolinium compounds. Instead, the incorporation of cerium allows tuning the bandgap of at least two other Ln 2 O 2 S compounds (Ln = Pr or Eu).

	A	Sample color and bandgap (eV) of A 1.8 B 0.2 O 2 S nanoparticles
	A	B	Ce	Pr	Eu	Gd
		Gd	2.8	4.3	[a]	4.7
		Eu	2.5	3.2	3.2
		Pr	2.6	3.6
		Ce	1.8
	B
	C

[a] No measure possible in this absorption range

Investigation of the bandgap nature using density functional calculations DFT

  Moreover, the PDF analysis tends to demonstrate a Ln 2 O 2 S crystalline structure without any defect. It is in accordance with the previous results obtained with our accurate estimations of the crystalline domain sizes of Gd 2(1-y) Ce 2y O 2 S nanoparticles for which we chose to use the Scherrer's formula that supposes perfectly ordered crystals (TableS2). calculations were performed in order to understand the unexpected direct nature of the bandgap. After a series of works,Mikami et al. concluded in 2006 from DFT that the bandgaps of Ln 2 O 2 S (and Ce 2 O 2 S in particular) were indirect.

	2.3
	2 O 2 S structure. The structural
	equivalence of the anisotropic Gd 2 O 2 S ~0.5 nanoparticles with Ln 2 O 2 S cells allows to seek from
	density functional theory (DFT) calculations an explanation of the bandgap and the optical
	properties of (Gd,Ce) 2 O 2 S nanoparticles.

2(1-y) Ce 2y O 2 S nanoparticles

  All products were used as received without further purification. In a typical synthesis of GdCeO 2 S, Gd(acac) 3 •xH 2 O (0.25 mmol), Ce(acac) 3 •xH 2 O (0.25 mmol), S 8 (0.032 mmol), Na(oleate) (0.50 mmol), OM (17 mmol), OA (2.5 mmol) and ODE

	Bimetallic oxysulfides were prepared via a solvothermal reaction including a mixture of
	organic solvents and metallic complexes. Oleylamine (OM; technical grade, 70 %), oleic acid
	(OA; technical grade, 90 %), sulfur (S 8 ; ≥ 99.5 %) and sodium oleate (Na(oleate); ≥ 99 %)
	were purchased from Sigma-Aldrich. 1-Octadecene (ODE; technical grade, 90 %) was
	purchased from Acros Organics. Gadolinium acetylacetonate hydrate (Gd(acac) 3 •xH 2 O;
	99.9 %) and cerium acetylacetonate hydrate (Ce(acac) 3 •xH 2 O; 99.9 %) were purchased from
	Strem Chemicals.

Acknowledgements

This work was supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11-IDEX-0004-02, and more specifically within the framework of the Cluster of Excellence MATISSE led by Sorbonne Université. We acknowledge funding for the PDF detector of our instrument from both the Institut de Matériaux de Paris Centre and the ANR project OxySUN (ANR-16-CE09-0005). The doctoral Program IPV of Sorbonne Université is acknowledged for funding of the PhD thesis of AMN. Sorbonne Université, CNRS, and the Collège de France are acknowledged for financial support. This work was granted HPC resources under the GENCI allocation A0030907320.

Authors contributions

Supporting information

EDS of selected powders, picture of the nanoparticles powder, UV-vis reflectance spectra, complementary data treatment assuming an indirect bandgap in the Kubelka-Munk model, TEM of selected samples, powder XRD patterns, average size of the crystallites, FTIR spectra, models of the nanoparticles, details about the PDF and the DFT+U, complementary information from UPS, calculated band structure of GdCeO 2 S slab. This material is available free of charge via the internet at http://pubs.acs.org.

Competing interests

The authors declare no competing interests.