Arnaud Wilhelm

Samuel Rivallant
email: samuel.rivallant@isae-supaero.fr

Jean-François Ferrero

Joseph Morlier

Sandwich shield subjected to bird impact: use of surrogate models for influencing parameters analysis and shield behaviour understanding

Keywords: Sandwich shield, Impact design, Bird strike, Parametric Study, Gaussian Process

) are varied and their effects on the shield behaviour and on the target protection are assessed. First, we try to establish an engineer's visualization by varying parameters 2 by 2 using three 5-levels fullfactorial designs of experiments (DOE). These three 2D DOE enable us to visualize precisely the different effects of each parameter. Then a full sensitivity analysis (6D) is performed using a Latin Hypercube sampling, to assess the possible interactions between parameters.

Surrogate models are constructed using the Gaussian Process framework to follow the variation of the outputs in the 6D design space. These surrogate models are finally studied using two statistical methods: the Sobol' method and the Morris method. The methodology developed in this study enables to improve the understanding of the behaviour of a shield under a soft body impact, as a first step towards a shield design tool.

Highlights

• A finite element parametric study is conducted on a sandwich shield under bird impact

• It appears that six parameters have a great influence on the shield behaviour

• The influence of these six design parameters on the sandwich behaviour and the target protection is studied

• The three sandwich core parameters are the most influent

Introduction

During its flight, one of the major threats an aircraft can encounter is the collision with a bird. Such collisions are known to occur frequently, and a 2008 study by the European Aviation Safety Agency [START_REF] Dennis | Bird strike damage & windshield bird strike final report[END_REF] estimates the occurrence at 186 per million flying hours. The possible damages of such an impact can be very diverse, due to the wide range of possible impacting scenarios, from multiple small flock birds (weighting approximatively 50 g) to heavy migratory birds (up to 4 kg). Moreover, the birds can impact all the forward facing structures, i.e. the nose, the windshield, the wings leading edge, the empennage, the engines, etc… To ensure the protection of passengers and crew on a commercial aircraft, the aviation authorities [2] [START_REF]Bird Strike Damage[END_REF] require that the plane should be able to continue its flight and land safely after a 1.82 kg (4 lb) impact at operational speed at sea level (typically around 175 m/s). In the case of an impact on the aircraft nose, the main danger is the failure of the pressurised bulkhead, electrical systems and equipment situated behind the radome instruments. Thus, to meet the certification requirements, a shield is placed in front of the bulkhead, behind the instruments (cf. Fig. 1). Due to their high specific stiffness and the good absorption properties of their core [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF], sandwich structures are the ideal candidates for such shields, and are used for this application by most aircraft manufacturers. In this work, we focus on the study of the behaviour of such a sandwich shield under the certification impact (1.82 kg at 175 m/s).

The first studies about bird strike, in the late 1970s, were focused on understanding the behaviour of a bird impacting a plate at such speed. Using extensive testing, Barber, Taylor and Wilbeck [START_REF] Barber | Characterization of bird impacts on a rigid plate: Part 1[END_REF][6] [START_REF] Wilbeck | Impact behavior of low strength projectiles[END_REF] showed that the bird behave similarly to a fluid during impact: the force transmitted to the target starts with a peak and is followed by a long plateau corresponding to a steady flow. They also showed that, during experiments, the bird can be À soumettre dans International Journal of Solids and Structures 3/27 substituted by an impactor made of gelatine with 10% porosity, with good improvements on repeatability. Consequently, these substitutes are used in tests in most latter studies.

Due to the great number of possible impacting scenarios (size and speed of bird, impacted structure, etc…) the literature on bird strike is quite large. In a 2011 review, Heimbs [START_REF] Heimbs | Computational methods for bird strike simulations: A review[END_REF] lists more than 190 numerical studies on bird strike, and identifies for each the impacted structure (windshield, leading edge, etc…), the impact case (mass of bird, speed, geometry), and the bird modelling strategy used. Three different numerical strategies are commonly used to model a bird: Lagrangian, Arbitrary Lagrangian Eulerian (ALE) or Smoothed Particles Hydrodynamics (SPH). These methods and their differences are presented in detail in [START_REF] Heimbs | Computational methods for bird strike simulations: A review[END_REF].

Concerning the impactor shape and material laws, a study by Airoldi and Cacchione [START_REF] Airoldi | Modelling of impact forces and pressures in Lagrangian bird strike analyses[END_REF] showed that the simulations were closer to test data when the material laws used represented water with 10% porosity and when the impactor was modelled as a cylinder with hemispherical caps and a length to diameter ratio of 1.6. More complex shapes and material laws have also been studied in [START_REF] Hedayati | A new bird model and the effect of bird geometry in impacts from various orientations[END_REF].

Only a few of these studies use a flat sandwich structure as target and focus on its behaviour.

We present here these 5 papers:

In 2006, Hanssen et al. studied a two layers sandwich panel [START_REF] Hanssen | A numerical model for bird strike of aluminium foam-based sandwich panels[END_REF]. Using simulation with the ALE strategy, there were able to represent the observed experimental behaviour, including failure at the clamping bolts. In a second part of their work, they used their model to minimize the core thickness of a simple sandwich with 0.8 mm thick aluminium skins and a core made of 150 kg/m 3 aluminium foam. The lighter sandwich able to stop the bird without front skin tearing had a 150 mm thick core.

In 2012, Hohe et al. reproduced this load case [START_REF] Hohe | Numerical and experimental design of graded cellular sandwich cores for multifunctional aerospace applications[END_REF] and showed that, using a graded core of three layers with increasing density, the front skin strain was spread on a greater surface, with the maximum strain being lower. They concluded that the use of a graded core could improve the resistance of the sandwich shield.

In 2012 also, Liu et al. studied a sandwich with aluminium honeycomb core and a two layers sandwich with the same core material and the same total height [START_REF] Liu | Dynamic Response of Bird Strike on Aluminum Honeycomb-Based Sandwich Panels[END_REF]. They created a finite element model to simulate the load case and were able to represent correctly both impacts.

In 2015, Hedayati et al., in a numerical study, assessed the effect of the position of the middle skin in a two layers sandwich, keeping the total height constant [START_REF] Hedayati | Effect of Using an Inner Plate between Two Faces of a Sandwich Structure in Resistance to Bird-Strike Impact[END_REF]. They showed that the minimal backward deflection is obtained for the shield with balanced cores (i.e. middle plate in the centre). In a second part, they showed that for cores with different densities, the best result was obtained with the lighter core facing the impactor.

The same year, Liu et al. studied the influence of skins and core thickness at constant mass on a simple sandwich with an aluminium foam core [START_REF] Liu | Dynamic response of bird strike on aluminium foam-based sandwich panels[END_REF]. They showed that both the backward deflection and the energy absorption of the panel increase when the core height increases and the skins thickness decreases. Then, they studied the effect of the middle plate position on a two layers sandwich with a constant total height. They showed that the minimal backward deflection is obtained with a first layer height equal to zero, i.e. a simple sandwich with a double front skin. This latter conclusion seems to be in contradiction with the conclusion of [START_REF] Hedayati | Effect of Using an Inner Plate between Two Faces of a Sandwich Structure in Resistance to Bird-Strike Impact[END_REF]. This difference can be explained by the fact that those two studies use slightly different boundary conditions (respectively riveted or clamped) and by the fact that in [START_REF] Hedayati | Effect of Using an Inner Plate between Two Faces of a Sandwich Structure in Resistance to Bird-Strike Impact[END_REF], no bonding is modelled between the skins and the cores.

All these results suggest that the behaviour of the shield is influenced not only by the core and skins properties, but also by the structural coupling between the two. Thus, when designing a new shield, it seems important not to study the skins and core properties separately, but to study all the shield design properties together.

According to this conclusion, in this work the skins and core properties are all considered as design parameters of the shield, along with the different heights and thicknesses. We then study the influence of these parameters on the behaviour of the shield but also on the target protection. This paper is organised as follows:

In Section 2, the case study is defined and the finite element model used is presented. Then, the different outputs which will be followed throughout the study are presented and, according to a previous study [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF], the most influential parameters are determined.

In Section 3, these parameters are studied 2 by 2 using full factorial designs of experiments (DOE). The effect of each parameter and its interactions with the other parameters are studied and physical interpretations are proposed for the observed phenomenon.

In Section 4, these 2D studies are expanded in full 6D by adding 100 new simulation points.

Surrogate models are then created to follow the outputs variations in the 6D design space using one of the main reference framework in the machine learning domain. Statistical methods are then used to analyse the effects of each parameter and the interactions with the others. The results obtained are confronted to the observations made on the 2D studies of the previous section. The radome shields used today in the industry have a geometry which is strongly dependant on all the other surrounding systems, and thus they are different from one aircraft to the other, but these shields are usually flat sandwiches, with an area of about 1 square meter, directly supported by the protected bulkhead. These sandwiches are typically made of two aluminium skins and a core made of metallic foam or honeycomb. The skin thickness range from 1.5 to 5 mm and the core height is usually about 100 mm.

In this work we will study the shield under the impact of a 1.82 kg bird at 175 m/s, representing the certification case. This impact represents an initial kinetic energy of 27.9 kJ.

Finite element model

In order to limit the number of geometrical and boundary conditions parameters, the limit conditions and the shield geometry have been simplified compared to a real engineering case. Moreover, to be able to simulate a great number of different shield designs at an acceptable cost, the finite element model is kept as simple as possible, and the material laws used are kept generic.

Geometry

The geometry chosen is an 800 mm square sandwich, supported by a rigid frame with a 400 mm square aperture in its centre. The sandwich has two skins 3 mm thick and a core with a height of 100 mm. Due to symmetry conditions, only a fourth of the shield is modelled using ABAQUS/EXPLICIT (cf. Fig. 2). The support is modelled with rigid shell elements R3D4. The sandwich shield is modelled as one part, with no possible debonding between the skins and the core. Five layers of reduced integration brick elements (C3D8R) are used for the back skin, ten layers of fully integrated elements (C3D8) are used for the core and five layers of reduced integration elements (C3D8R) are used for the front skin. During the parametric study, the core and skins thicknesses will change but the number of elements used through thickness is kept constant.

The in-plane element size is 10 mm by 10 mm, for a total number of elements for the sandwich of 32000. The impacting bird is modelled as a hemispherical ended cylinder with a radius of 55 mm and a length of 220 mm. 3490 Lagrangian reduced integration elements (C3D8R) are used, with enhanced hourglass control. The contact between the bird and the front skin, and between the back skin and the support, is modelled with ABAQUS general contact without friction. The total simulated time is 6 ms, enough to simulate the whole impact and the shield rebound.

Material models

In order to model the bird material behaviour, a tabulated equation of state representing the behaviour of water with porosity is used. The elements are deleted when their true strain become greater than 500%, in order to limit the decrease of the simulation time step due to highly distorted elements while allowing to represent the flow of the bird. The bird density is 955 kg/m 3 , which gives a 1.82 kg bird.

The shield design parameters being varied during the parametric study, we give here the values for the central point of the design space. This reference case represents a shield with aluminium skins and an aluminium honeycomb core.

The skin material is represented using an isotropic material law. A Young modulus of 72 GPa and a Poisson coefficient of 0.33 are used, with a density of 2800 kg/m 3 . The plastic behaviour is modelled by ABAQUS Johnson-Cook isotropic hardening (equation ( 1)).

𝜎 " = A + B. ε ) * (1)
where A = 300 MPa, B = 400 MPa and n = 0.27. This model is used because it is able to represent the plastic behaviour of many different metals with only three parameters [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF]. No skin failure is modelled, but the maximal skin in-plane strain will be studied to ensure that this hypothesis stay valid.

To simulate the core behaviour, a generic material law is implemented using a VUMAT user routine. An elasto-perfectly plastic law is used for all directions, followed by a densification 7/27 for the out-of-plane compression. Such a strain-stress curve have been shown to be typical for the cellular materials used as sandwich cores [START_REF] Gibson | Cellular Solids: Structure and Properties[END_REF]. Uncoupled shearing and compressive behaviour is assumed. The material law parameters are adapted from [START_REF] Kolopp | Modeling impact on aluminium sandwich including velocity effects in honeycomb core[END_REF]:

-In-plane modulus: 𝐸 , = 𝐸 -= 1 MPa.

-In-plane plateaus: 𝜎 " , = 𝜎 " -= 0.9 MPa.

-In-plane shear modulus: 𝐺 ,-= 5 MPa.

-In-plane shear plateau: 𝜏 " ,-= 1.2 MPa.

-Out-of-plane shear modulus: 𝐺 ,-= 𝐺 -1 = 200 MPa.

-Out-of-plane shear plateaus: 𝜏 " ,-= 𝜏 " -1 = 1.2 MPa.

-Out-of-plane modulus: 𝐸 1 2 = 1000 MPa.

For out-of-plane compression, a plateau followed by densification behaviour has been modelled using equation (2).

𝜎 1 = 𝜎 " 2 + 𝐶. (𝑒 6 6 7 -1) (2)
with 𝜎 " 2 = 0.9 MPa the core crushing plateau, and C = 0.01 MPa and 𝜀 < = 0.37 defining the densification phase. The core density is set to 50 kg/m 3 .

Outputs studied

With this generic finite element model, it is possible to simulate the impact of a bird on different shields, and to study the effect of the shield definition on the target protection. The simulation of one design case takes approximatively 2 hours on 4 CPU, using the CALMIP (Calculateur Midi-Pyrénées) supercomputer. This relatively small computing time allows simulating numerous design points during the parametric study. Fig. 3 presents the evolution of the reference shield shape during impact. In such a parametric study, it is obviously impossible to analyse all the differences between all the finite element simulations, due to the important number of design case tested. It is then necessary to define simple outputs that can be followed during the parametric study and that define the shield behaviour.

Here, we decided to follow four protection criterions, which could be used as design criterion in a real engineering case, and representing the protection capacity of the shield:

-The maximum total force applied to the support F >?)) . The total force is computed at each time step as the sum of all the contact forces between the back skin and the support. This force is then filtered using a moving mean with a window of 0.15 ms and F >?)) is defined as the maximum along time. For the reference shield, F >?)) is given in Fig. 4, with the total force applied by the bird on the front skin.

-The maximum pressure applied to the support P A?)) . This pressure is computed using the contact forces between the back skin of the shield and the support. The contact forces are filtered using a moving average with a window of 0.15 ms along time and a Gaussian filter on a 3*3 nodes window along space. P A?)) is then defined as the maximum along time and space of the pressure obtained. Cf. Fig. 5 for the evolution of P A?)) in the reference case.

-The maximum backward deflection of the back skin δ CA . This criterion is often used in the industry since an important backward deflection can endanger critical systems situated behind the pressurised bulkhead.

-The maximum in-plan strain of the front skin ε DA . Even if no skin rupture has been modelled in this work, this output allows to differentiate two shield designs by knowing how close they are of a potential front skin tearing.

Using these four criterions, it is possible to assess the capacity of a shield to protect the target, but these outputs give no information about the shield deformation. The goal of this parametric study being to understand the behaviour of the shield and its influence on the target protection, it is necessary to have outputs describing how the shield is deformed.

In a previous work [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF], we described a behaviour analysis tool allowing to decompose the deformation of a shield into three modes: Indentation, Bending and Crushing. The main principle of this tool is to project at each time step the vector representing the shield deformation on a basis defined a priori and describing the three deformation modes. It has been shown that this decomposition describes efficiently the deformation of any shield with only a small residue (less than 10% in norm) [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF].

9/27

Using this tool, it is then possible to extract the maximum norm of each deformation mode along time. These three behaviour criterions will give information about how the shield behaves and is deformed during impact without the need to analyse in detail each simulation.

Choice of parameters to study

The finite element model described here presents an important number of parameters. Thus, it is necessary to choose which parameters to take into account before conducting the parametric study. In [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF], we conducted on this finite element model a screening analysis to identify the most influential parameters on the shield behaviour, using the deformation-based behaviour analysis tool presented above. This screening has been extended to identify influences on the shield behaviour (3 outputs) and the four protection criterions (4 additional outputs).

15 parameters or group of parameters have been studied independently. Along with its reference value, each parameter has been assigned a minimal and maximal value. 31 cases have then been simulated: the reference point and the minimal and maximal case for each parameter (One factor at a time DOE). For each of these simulations we measured the -The most influential parameter is clearly the core crushing plateau σ ) F . Any change in the value of this parameter induces important change in the shield behaviour.

-The second most influential parameter is the core out-of-plane shearing plateau τ ) F .

-Then the parameters rank as follows: the core height H F , the front skin yield strength A DA , the support aperture size L J) , the front skin thickness t DA , the elastic modulus of the front skin, the out-of-plane elastic modulus of the core, and then the other parameters.

-It appears that the in-plane properties of the core have very little influence on the shield behaviour, and that the back skin design has little influence compared to the front skin design.

-The effects of the parameters on the shield behaviour can be strongly non-linear, and sometimes even non-monotonic.

According to these conclusions, this parametric study will be focused on the 6 most influential parameters: The core out-of-plane crushing plateau 𝜎 " 2 , the core out-of-plane shearing plateau 𝜏 " 2 , the core height 𝐻 2 , the support aperture size 𝐿 N" , the front skin thickness 𝑡 PQ and the front skin yield strength 𝐴 PQ . Moreover, it is necessary to choose a DOE allowing us to study the possible nonlinearities and interactions, as they were not studied in the screening study.

2D parametric studies

In order to be able to analyse precisely the effects of each parameter, and to be able to visualize the results more easily, the six identified parameters are first studied 2 by 2: the two core out-of-plane plastic plateaus 𝜎 " 2 and 𝜏 " 2 , the two front skin parameters 𝐴 PQ and 𝑡 PQ , and the geometric parameters 𝐻 2 and 𝐿 N" . To be able to assess the nonlinearity and the possible interactions, 5 levels by parameters are chosen. When varying the inputs pairwise, remaining parameters are kept at the baseline values (presented in Section 2.1). The minimal and maximal parameters values used are presented in Table 1 and are the same as in [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF],

representing materials typically seen in the industry for such applications.

À soumettre dans International Journal of Solids and Structures 11/27 The designs of experiments used for these three 2D series are represented in Fig. 6. The reference case is indicated with a full red dot. For the 'core' DOE, the choice was made to simulate more designs toward the lower bound of both 𝜎 " 2 and 𝜏 " 2 because the screening study showed that a decrease of those parameters has more effect on the shield behaviour than an increase.

Effect of the core design

For each of the 25 simulations of this DOE, the seven criterions (3 behaviour criterions and 4 protection criterions) are measured. These results are given on Fig. 7 as a surface reconstructed from the simulation points (black) using linearly interpolated scheme. The red dot indicates the reference case and the value of the outputs at this point is shown by a red line on the colour scales. We can see on this figure that the effects of the core design are strongly non-linear and that an interaction exists between the two design parameters. Two main effects can be observed.

First, there is a strong effect of 𝜎 " 2 , which appears only when 𝜎 " 2 is smaller than 0.8 MPa.

This effect can be observed on Fig. 7 for Indentation, Crushing, 𝐹 TU"" and 𝜀 PQ . For a decreasing value of 𝜎 " 2 , the shield shows increasing values of Indentation and Crushing modes, associated with reducing values of 𝐹 TU"" and 𝜀 PQ . Thus, this strong crushing behaviour is associated with a smaller force transmitted to the support and a smaller strain of the front skin. This effect can be explained by the fact that the force transmitted to the support is linked with the force transmitted between the two skins by the core. With a low crushing plateau of the core, this force is necessarily smaller. Regarding 𝜀 PQ , its decrease may be explained by the fact that a strong crushing behaviour induces an indentation which is less located at the centre of the shield.

So, for a low core out-of-plane crushing plateau, the shield has a deformation dominated by indentation and crushing, with comparatively small bending, whereas for a higher core out-ofplane crushing plateau, the shield shows mainly indentation.

This change of behaviour can be seen clearly on The second effect which can be observed on Fig. 7 is the increase of Bending, 𝑃 TU"" , 𝛿 XQ and 𝜀 PQ for a low core out-of-plane shearing plateau associated with a strong core crushing plateau. The effect of 𝜏 " 2 can be explained by the fact that the core out-of-plane shearing plateau have a direct influence on the sandwich bending resistance. Thus, a decrease of 𝜏 " 2 will lead to an increase of bending, and this effect is enhanced for a high core crushing plateau, which increase the sandwich crushing resistance. This increase of bending leads to an increase of the back skin backward deflection 𝛿 XQ , which can be seen clearly on Fig. 8 where the shield c) (with low shearing plateau and a high crushing plateau) shows an important backward deflection. Similarly, the shield c) shows an indentation which is more located at the centre of the shield, which can explain the increase of the front skin maximum strain 𝜀 PQ .

Another consequence of the increase of bending is the increase of the maximum pressure transmitted to the support 𝑃 TU"" . This consequence is explained by the fact that the bending of the shield localises the contact between the back skin and the support. Thus 𝑃 TU"" increases without any increase of 𝐹 TU"" .

In summary, the core design strongly influences the behaviour of the shield. A low core crushing plateau induces a behaviour dominated by the core crushing. On the contrary, a low core shearing plateau associated with a strong core crushing plateau leads to a shield behaviour with strong bending. For the reference case, the shield is deformed mainly by indentation.

À soumettre dans International Journal of Solids and Structures 14/27

Effect of the front skin design

The results for the 25 simulations studying the front skin design are presented on Fig. 9. First, we can see that nearly all the outputs vary less than on Fig. 8, which shows that, using these ranges of variation, the core design has more influence than the front skin design. This conclusion is coherent with the conclusions of [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF], where the core plateau parameters where identified as the most influent parameters. Only the front skin maximum strain 𝜀 PQ is more influenced by the front skin design. It is clear from Fig. 9 that the front skin thickness 𝑡 PQ and its yield strain 𝐴 PQ play very similar roles. In fact, a shield with a thick skin and a low yield strain seems to behave similarly to a shield with a thin skin and a high yield strain. The only outlier design is the shield with 𝑡 PQ = 2 𝑚𝑚 and 𝐴 PQ = 500 𝑀𝑃𝑎, especially for the Crushing output. When observing this finite element simulation in details, it appears that this is due to the buckling of the front skin, under the compression stresses due to the sandwich bending. This front skin buckling does not appear on any other case simulated in this study. It may come from a thin skin submitted to high stresses, due to a high yield strain.

All the other designs show that a weak front skin (thin and with a low yield strain) increases the indentation and bending of the shield, while a strong front skin reduces it, with nearly no effect on the crushing behaviour. In fact, the front skin strength has an influence on the importance of the shield deformation, as can be seen on Fig. 10, where the three shields shown are the two extreme cases and the reference. On Fig. 9, it appears that the increase of the shield deformation due to a weaker front skin is associated with an increase of 𝑃 TU"" (which is explained as in section 3.1 by an increase of bending localizing the contact between the back skin and the support), but also with a decrease of 𝐹 TU"" . This effect may come from the fact that, due to the important deformation of the shield, the impacting bird is slowed down more progressively by a weak front skin than by a strong one. Indeed, the maximum force transmitted to the support occurs at 𝑡 = 1.86 𝑚𝑠 for the shield a) in Fig. 10, at 𝑡 = 1.62 𝑚𝑠 for the shield b), and at 𝑡 = 1.56 𝑚𝑠 for the shield c).

The increase of the shield deformation is also associated with an increase of the back skin backward deflection (clearly seen on Fig. 10, shield a)), and an increase of the front skin maximum strain. This effect is due to the strong indentation of the sandwich, localizing the front skin strain at its centre. We can notice that the strong indentation also limits the capacity of the bird to flow radially. This phenomenon makes any front skin rupture dramatic since the loss of the front skin rigidity will induce a strong indentation, which will increase the front skin loading by limiting the bird radial flow.

Effect of the geometrical parameters

In the third 2D design of experiments, the influences of the core height 𝐻 2 and the support aperture size 𝐿 N" are studied. The results of this DOE are given on Fig. 11. First, it appears that a decrease of 𝐻 2 leads to a decrease of indentation and an increase of bending. The increase of bending is explained simply by the bending rigidity of a sandwich, which is directly influenced by its core height. This increase of bending is here associated with an increase of the backward deflection and the maximum pressure transmitted to the support, which is coherent with the phenomenon already identified.

Secondly, the increase of 𝐿 N" seems to have the same effect as a decrease of 𝐻 2 on Bending and 𝛿 XQ . This is explained by the fact that this parameter represents the size of the unsupported area of the sandwich. The maximum pressure transmitted 𝑃 TU"" seems to oscillate with 𝐿 N" . This phenomenon doesn't have any physical justification and seems to be caused by the relative position of the support and back skin elements in the simulation, which change when the parameter 𝐿 N" change. Similarly, no physical explanation has been found to justify the evolution of 𝐹 TU"" .

Thus, it seems that the core height and the support aperture size have similar influences on the shield behaviour, with a stronger effect for low values of 𝐻 2 and high values of 𝐿 N" . To illustrate this effect, the deformed shapes of three shields are represented in Fig. 12. The shield a) has a high core and a small aperture, the shield b) is the reference, and the shield c) has a thin core and a large aperture.

Conclusion on the 2D studies

With these three "two factors at a time" DOE, it was possible to explore in detail the influence of each of the six parameters chosen. The use of 5 levels per parameters allowed us to observe the nonlinearity of the different effects, and studying the parameters 2 by 2 showed that important interactions can exists between the different parameters.

These 2D studies confirmed that the most globally influent parameters are the core crushing and shearing plateaus, but they also showed that depending on the output studied, other parameters can be more influent. Different shield behaviours have been identified for different shield designs, and physical explanations for these different behaviours have been proposed.

The links between the behaviour criterions (Indentation, Bending and Crushing), and the target protection criterions (𝑃 TU"" , 𝐹 TU"" , 𝛿 XQ , 𝜀 PQ ) have been identified and explained when possible.

But all these conclusions have been obtained on a strongly reduced design space. In the full 6D space, only three 2D planes have been studied. Moreover, the results showed that interactions can exist between the different parameters. To study these interactions, and to assert the results obtained with these 2D studies, it is necessary to conduct a full 6D study, in which all the parameters are varied together.

Expanding the parametric study in 6D

A 6D study is far more complex to conduct than a 2D study. First, it is impossible to use a simple full factorial DOE, since the number of simulations needed increase exponentially (5 levels for six parameters lead to 5 b = 15625 simulations). Second, the visualisation of the À soumettre dans International Journal of Solids and Structures 18/27 results is not possible is full 6D and can only be done using 2D slices. Third, it is complicated to analyse the results since it is not possible to visualise them globally. It is then necessary to use statistical methods to measure the influence of the different parameters.

To handle these issues, the usual approach is to construct and train a surrogate model. This model is a mathematical function and its goal is to approximate the output studied everywhere in the design space, and not only on the tested points. This approach is widely used in the field of engineering design, and described in details in the revue [START_REF] Queipo | Surrogate-based analysis and optimization[END_REF] and the book [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]. The main steps of this approach are as follows:

-First, the choice of the design points to simulate in order to train the surrogate model. The model will interpolate the output between these known points, and thus the choice of this design of experiment is crucial.

-Second, the construction, training and validation of the surrogate model. It is necessary to choose the form of the mathematical functions which will be used to construct the surrogate model, as it will constrain the precision of the model. The validation consists to check if the model is able to approximate efficiently the output. If not, it is possible either to change the form of the surrogate model or to return to the first step to add more known points.

-Third, the analysis of the model obtained, in order to measure the effects of each parameter and their interactions, and to conduct physical interpretations. Global methods of analysis will be used to validate (or not) the observations made on the 2D studies in Section 3.

Choice of the 6D design of experiments

Usually, the choice of the DOE depends strongly on the form of the surrogate model to train, which is chosen depending on the a priori available knowledge of the output. Here, the 2D studies showed that the outputs behaviour can be strongly non-linear and with interactions. It is then necessary to have a surrogate model able to adapt to a wide range of behaviour.

Moreover, in this work seven different outputs are studied and thus seven different surrogate models will be constructed. As each simulation allows measuring the seven outputs, all seven surrogate models will be trained on the same set of known points. For these reasons, the DOE used needs to be as general as possible, to adapt to the different situations.

A lot of different methods exists to create a design of experiments [START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][21] [START_REF] Simpson | Sampling strategies for computer experiments: design and analysis[END_REF]. We choose here to use a Latin Hypercube Sampling (LHS) strategy, because it allows to choose a priori the number of design points. In a LHS design of experiments with 𝑁 Q points, each parameter has 𝑁 Q levels and takes each level once and only once. This type of DOE is of particular the rows of an identity matrix), but this approach can produce designs with very different uniformity. In order to obtain a good LHS, the original one can be optimized to increase the uniformity. For example, in [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] a strategy using the maximin [START_REF] Johnson | Minimax and maximin distance designs[END_REF] criterion is presented.

In this work, 73 sampling points have already been simulated during the 2D studies presented in Section 3. In order to choose 100 new design points to be simulated, we adapted the strategy presented in [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF] to optimize the maximin criterion of the complete 173 points DOE.

The 100 new points chosen are then simulated and the seven outputs are measured for each case. The 173 known points can then be used to create the surrogate models.

Construction and validation of the surrogate models

As seven models have to be constructed, it would be possible to choose different frameworks for each, but for practical reasons we choose to create the seven surrogate models using the same framework. A lot of different framework of surrogate models exist [START_REF] Queipo | Surrogate-based analysis and optimization[END_REF][20] [START_REF] Simpson | Sampling strategies for computer experiments: design and analysis[END_REF], and we choose here to create our surrogate models using the Gaussian Process (GP) framework, described in detail in the book [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. All the computations for the creation and analysis of these surrogate models are made on Matlab using the GPML toolbox available with the book.

GP are used because they can reproduce many different behaviours. Moreover, their statistical framework allows us to estimate the approximation error at each point of the design space, which can be very useful to estimate the precision of the surrogate models.

Thus, for each of the seven output, a surrogate model is constructed using a Matérn 𝜈 = 3/ 2 kernel with Automated Relevance Determination, and a constant mean [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. The nine hyper-parameters (the mean, the signal variance, the noise variance and the six length-scales) are optimised using the log marginal likelihood maximisation function provided in the GPML toolbox [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. In order to ensure a good convergence of the hyper-parameters optimisation, the data are first normalized using Matlab zscore function, to set their means to zero and their variances to one in each direction.

The obtained surrogate models are then tested using the Leave-One-Out (LOO) method. This method consists of training the surrogate model on all the data points minus one, and then of testing it on the last data point. By repeating this process on each data point, a mean prediction error can be computed. This mean error is presented for each output on Fig. 13, in percentage of the variation range of the output. We can see on Fig. 13 that the mean LOO error is lower than 5% for all the outputs studied.

This shows that the surrogate models constructed are able to represent globally the evolution of the seven outputs. Thus, these models will be used to analyse the influence of each parameter.

Analysis of the full 6D design space

Methods used

In order to analyse a six dimensions space, we choose to use two complementary statistical methods: the Morris method [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] and the Sobol' method [START_REF] Sobol′ | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF].

-The Morris method consists in measuring the elementary effect of each parameter around an initial point, using a one factor at a time design of experiments. By using multiple initial points, evenly distributed in space, it is then possible to compute for each parameter k a mean effect d i and its variance S i . Then, the greater d i is, the more the parameter k is influent, and the greater S i is, the more this influence is non-linear or with interactions with other parameters.

-The Sobol' method allows to compute the parameters sensitivity indices and their total sensitivity indices. The sensitivity indice represents the fraction of the total variance of the output contributed by the parameter individually, while the total sensitivity indice represents the influence of the parameter and all its interaction with other parameters. To estimate these indices, we use the Matlab toolbox GSAT [START_REF] Cannavó | Sensitivity analysis for volcanic source modeling quality assessment and model selection[END_REF].

Thus, the Morris and Sobol' methods both estimate the influence of each parameter, but in slightly different ways: the Morris method allows measuring the linear effect of a parameter and the variance of this effect, without being able to distinguish between non-linearity and interaction with other parameters, while the Sobol' method measure the total influence of a parameter isolated and with interactions, without measuring the non-linearity of this À soumettre dans International Journal of Solids and Structures 21/27 influence. In this work, we used both approach in order to obtain a maximum of knowledge about the output behaviour.

For both these approaches, it is necessary to use the unit hyper-cube as design space. Thus, the parameter space is normalized appropriately for all the surrogate models before any computation.

Results and interpretations

The results given by the Morris method are presented in Fig. 14. For each of the seven surrogate models, we used 100 000 initial points pseudo-randomly distributed to estimate the mean effects 𝑑 l and their variance 𝑆 l . In order to compare the influence of each parameter, they are placed in the 𝑑 , 𝑆 space for each output in Fig. 14, where the dotted diagonal represents the line 𝑑 = 𝑆. To visualize the sign of the mean effects, the parameters are represented by a • when the mean effect is positive and by a * when it is negative. A positive/negative effect means that the output increases/decreases when the parameters increases. First, we can see on Fig. 14 that the parameters are mainly placed near or above the 𝑑 = 𝑆 line, which shows important non-linearity or interactions. On the other hand, Fig. 15 shows that the interactions between parameters are not so important. Thus, we can conclude that the effects of the different parameters are strongly non-linear. This conclusion is coherent with the observations in Section 3.

By observing these two figures, it is clear that the more influential parameters are the core properties and then the core height. This confirms the parameter ranking obtained in [START_REF] Wilhelm | Study of the deformation of a sandwich shield subjected to bird impact: A behaviour analysis tool using vector decomposition[END_REF] and the outputs ranges of variations observed in Section 3. This is particularly true for two behaviour outputs (Indentation and Crushing), where the core out-of-plane crushing plateau 𝜎 " 2 seems to be the only influent parameter. This can be seen more clearly on Fig. 16, where the 173 simulations are represented in the (Indentation, 𝜎 " 2 ) and the (Crushing, 𝜎 " 2 ) planes.

There is a clear correlation between the output value and the 𝜎 " 2 parameter value. The Morris method also enables us to confirm the signs of the effects of the different parameters. For example, the effect of the core height 𝐻 2 is negative on Bending and backward deflection 𝛿 XQ , while positive on front skin maximum strain 𝜀 PQ , which is coherent with the local observations of Section 3.3 (Fig. 12). These signs also show that the two core material parameters 𝜎 " 2 and 𝜏 " 2 have always opposed effects, which confirms that the two 'extreme' designs of core material are either a core with high crushing plateau and low shearing plateau or a core with low crushing plateau and high shearing plateau.

By observing Fig. 15 and Fig. 14, we can also see that the two front skin parameters (its thickness 𝑡 PQ and yield stress 𝐴 PQ ) have effects which are quite similar, and with same sign, on all the outputs. Once again, it confirms the local observations of Section 3.2.

Conclusions on the full 6D study

In this section, the full 6D design space was studied using a set of 100 simulations spread evenly using a Latin Hyper-square DOE. Using these simulations and the 73 simulations already done in Section 3, one surrogate model was constructed for each output studied.

After validation, it was possible to use these models to analyse the design space using two different methods: the Morris method and the Sobol' method. We saw that these two methods are complementary and that using the two in parallel allowed us to measure both the nonlinearity of the parameters effects and the interactions between parameters.

Using the results obtained, we saw that, on the studied design space, the parameters have effects which are often strongly non-linear, but with only small interactions. Then, the results were compared to the local observations of Section 3. Both the ranking of parameters with respect of their influence, and the signs of the parameters effects were coherent with local observations.

Thus this 6D study allowed us to confirm globally the different behaviours observed in Section 3. But it is important to note that, due to the strong non-linearity of the parameters effects, these behaviours can be false locally.

To tackle the problem of bird impact on sandwich shield (numerous parameters, timedependant response…), DOE and machine learning methods were used in this study with finite elements simulations. After an initial screening to reduce the number of parameters to take into account, Gaussian Process were used to create surrogate models in a 6D domain, using smart sampling. This method enables to analyse the effects of the six most influential parameters and their interactions on the shield behaviour (deformed shape and protection criterions) within a reasonable computation time.

The main conclusions of this study are the following:

-First, the two most influential design parameters are clearly the core out-of-plane crushing plateau and the core shearing plateau. This means that during the design of a shield for soft impact, great care should be taken in the choice of the core material. We also saw that these two parameters have effects that are opposed. Thus, a shield with higher crushing plateau and lower shearing plateau will transmit more stresses (and more localized) to the back skin and the support, while reducing the area of core crushed and increasing backward deflection. Yet, when changing the core material designers usually stay in the same category of cellular material (ex. aluminium honeycombs), which mean changing both plateaus parameters in the same direction (either increasing or decreasing).

-The two skin parameters (thickness and material yield stress) have very similar effects. A stronger skin (i.e. higher thickness or higher yield stress) will tend to reduce the shield indentation and flexion, thus spreading more evenly the core crushing and reducing backward deflexion and front skin maximum strain. On the other hand, these two parameters have very different effects on the shield mass. For a lightweight shield, it is then more advantageous to increase the front skin yield strain while reducing its thickness.

This conclusion shows that a soft impact problem is very different from a hard impact problem, where a softer front skin is usually beneficial [START_REF] Kolopp | Experimental study of sandwich structures as armour against medium-velocity impacts[END_REF]. This difference may come from the fact that, for a soft impact case, the strain gradients are smaller. This implies that the deformations inside the shield are more global and it appears that all the parts of the shield are working together at the same time, and not one after the other (as can be the case in hard impact). Thus the sandwich shield has to be studied as a whole and not as a succession of different layers.

-The core height has a strong effect on the shield bending and backward deflection, which is coherent with usual sandwich behaviour. Thus, a higher core is usually beneficial for target protection.

À soumettre dans International Journal of Solids and Structures

25/27

For this work, we deliberately choose a quite simple finite element model in order to simulate numerous possible designs, and the conclusions have to be taken with precautions. In particular, the core material model used here is decoupled, which means that the real interactions between the two core plateaus parameters are probably more important than calculated. Another limitation of this study is the fact that we did not take into account the possible front skin rupture. In the 173 configurations simulated, none showed front skin strains greater than 15%, which shows a posteriori that this hypothesis was valid.

Nonetheless, we have to keep in mind that the front skin rupture would drastically change the shield behaviour, as losing the front skin rigidity would mean losing the "sandwich" behaviour.

At last, thanks to this new understanding, general rules have been proposed to orient the design of better shields. Moreover, the surrogate models created here are a first step towards a design tool to help engineers, as they can be used in an optimisation loop in order to find an optimal shield with only a few more simulations.

Of course, such a numerical study requires validation, and an experimental campaign is currently being conducted to validate the main conclusions.

Fig. 1 .

 1 Fig. 1. Shield in the nose of an A320 aircraft

À 27 2

 27 soumettre dans International Journal of Solids and Structures 5/Case study description

Fig. 2 .

 2 Fig. 2. Finite element model of the sandwich shield

Fig. 3 .

 3 Fig. 3. Reference case: evolution of deformed shape during impact

Fig. 4 .Fig. 5 .

 45 Fig. 4. Reference case: contact forces (bird on front skin and back skin on support)

F

  bird / front skin F back skin / support Time (ms) Force (kN) À soumettre dans International Journal of Solids and Structures 10/27difference of shield behaviour with respect to the reference case, in order to rank the 15 parameters (mean rank on the 7 outputs). The main conclusions of this screening study are as follows:

Fig. 6 .

 6 Fig. 6. Designs of experiments used for the three 2D studies

Fig. 7 .

 7 Fig. 7. Effects of the core design on the seven outputs

Fig. 8 :

 8 the shield in a) (with a low core crushing plateau) shows a strong Crushing behaviour whereas the shield in b) (the reference case) shows mainly Indentation. On this figure, only the outline of the core in the XZ plane is represented (with the second half reconstructed by symmetry for easier visualization). The À soumettre dans International Journal of Solids and Structures 13/27 undeformed profiles of the shields are represented by the dashed line and the rigid support is shown in black.

Fig. 8 .

 8 Fig. 8. Behaviour of the shield during impact for three different core designs (core outline in xz plane)

Fig. 9 .

 9 Fig. 9. Effects of the front skin design on the seven outputs

Fig. 10 .

 10 Fig. 10. Behaviour of the shield during impact for three different front skin designs (core outline in xz plane)

Fig. 11 .

 11 Fig. 11. Effects of the geometrical parameters on the seven outputs

Fig. 12 .

 12 Fig. 12. Behaviour of the shield during impact for three different geometrical designs (core outline in xz plane)

Fig. 13 .

 13 Fig. 13. Mean LOO error in % of the variation range for each output studied

Fig. 14 . 27 Fig. 15 .

 142715 Fig. 14. Analysis of the seven surrogate models using the Morris method

Fig. 16 . 27 Fig. 16

 162716 Fig. 16. Effect of 𝜎 " 2 on Indentation and Crushing

Table 1 :

 1 Minimal and maximal values of the six parameters studied

	Lower limit	Reference	Upper limit	Units

  À soumettre dans International Journal of Solids and Structures 19/27 interest for computer experiments since no point is tested twice, and this property hold true if any parameter appears irrelevant. Constructing a LHS is simple (by random permutation of

Funding

The authors want to thank BPIFrance and the Région Occitanie for their financial support through the FUI project SAMBA (Shock Absorber Material for Bird-shield Application). The authors gratefully acknowledge the members of this project: Stelia, Airbus Group, Cedrem, Esteve, Ateca, Nimitech, Hutchinson and I2M. This work has benefited from access to the HPC resources of CALMIP (Calcul en Midi-Pyrénées).


Xiao Lin

German Bridoux

Adrien Gourgout

Gabriel Seyfarth

Steffen Krämer

Marc Nardone

Benoît Fauqué

Kamran Behnia

Critical Doping for the Onset of a Two-Band Superconducting Ground State in SrTiO 3δ

à 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.

Superconductiviy is induced in insulating SrTiO 3 by introducing n-type charge carriers through chemical doping [1] and survives over three orders of magnitude of carrier concentration. The transition temperature, T c , peaks to 0.45 K around a carrier density of n H ∼ 10 20 cm -3 [2]. A superconducting dome has also been detected in the metallic interfaces of SrTiO 3 [3] when the carrier density is modulated by a gate voltage bias [4]. In unconventional superconductors, such as high-T c cuprates, superconducting domes are often attributed to the proximity of a magnetic order or a Mott insulator. The recent discovery of superconducting dome in gate-tuned MoS 2 [5] in absence of a competing order, however, highlights the limits of our current understanding of the interplay between carrier concentration and superconductivity and motivates a fresh reexamination of superconducting domes. In the specific case of SrTiO 3 , superconductivity occurs in the vicinity of an aborted ferroelectric order [6] and survives deep inside the dilute metallic regime when the Fermi temperature becomes more than one order of magnitude lower than the Debye temperature [7]. This is a second puzzle in addition to the one raised by the drop in T c on the overdoped side. These two questions, raised at the opposite limits of the superconducting dome, remain unsettled.

According to band calculations [9][10][11], doping SrTiO 3 with n-type carriers can fill three bands one after the other. Once the critical threshold for the occupation of a band is attained, a new Fermi surface concentric with the previous one emerges. Previous studies of quantum oscillations in bulk doped SrTiO 3 [7,[12][13][14] have detected both multiple-frequency [7,13,14] and singlefrequency [7,14] oscillations at different doping levels, but did not determine these critical doping levels. Moreover, according to tunneling experiments, doped SrTiO 3 beyond a carrier density of 10 19 cm -3 is a multi-gap superconductor [8]. The interplay between multi-band occupation in the normal state and multi-gap superconductivity has been a subject of recent theoretical attention [20].

We present a systematic study of quantum oscillations and superconducting transition as a function of carrier concentration, n H extended down to 10 17 cm -3 , two orders of magnitude below the range of the tunneling study [8], and find three new results. First of all, the two critical dopings [11], n c1 and n c2 are identified and the magnitude of the cyclotron mass and the Fermi energy of each band are determined. Second, we find that n c1 , the threshold of occupation of the second band, separates two distinct regimes of superconductivity. Below n c1 , the superconductor is single-band with a large T c /T F . Above n c1 , it becomes two-band with a T c , which fails to keep the same pace with increasing T F . Finally, we find that the lowest band presents a deviation from parabolic dispersion near n c1 . We conclude that the attractive interaction between electrons, remarkably strong when the chemical potential is near the bottom of this band, significantly weakens when it shifts upward and other bands are occupied. This feature, combined to the low cutoff frequency imposed by the small Fermi energy, drastically limits the choice of the bosonic mode exchanged by Cooper pairs. The top panel in Fig. 1 shows the variation of the frequency of the quantum oscillations with doping in eighteen different samples of n-doped SrTiO 3 [See the supplement for an extensive discussion of the samples]. The samples studied were either obtained by annealing stoichiometric SrTiO 3 in vacuum [7,15] or commercial niobium-doped samples. In most cases, what was measured was the Shubnikov-de Haas(SdH) effect (quantum oscillations of the magnetoresistance). In seven samples, the Nernst effect was also measured and found to display giant oscillations with a frequency identical to the SdH effect as reported previously [7]. As seen in the figure, the main frequency smoothly evolves with doping. Moreover, at two critical doping levels new frequencies emerge. Panel b shows the evolution of the superconducting transition temperature of the same samples. The most un- derdoped samples (n H < 4 × 10 17 cm -3 ) did not show a superconducting transition down to 60 mK. Nevertheless, they presented a sharp Fermi surface, indicating that superconductivity is preceded by (or concomitant with) the intersection of the chemical potential and the bottom of the conduction band. The most striking feature of the figure is a clear change in the slope of T c (n) in the vicinity of n c1 . This is the first new result of this study.

The bottom panel in Fig. 1 compares our data with the early work reported by Schooleyet al. [2]. For each sample, the error bar represents the width of resistive transition, i.e. the interval between two temperatures FIG. 2: Low-temperature resistivity of reduced (a) and Nbdoped (b) single crystals. Note the continuous evolution of resistive transitions in reduced samples. As the doping increases, the transition first shifts to higher temperatures and then remains more or less constant. corresponding to 0.1 and 0.9 drops in resistivity (See our data in Fig. 2). While the two sets of data match roughly with each other, the structure in the dilute limit is far more clear in our data. The panel shows also the interface dome, restricted to a narrow window far above the concentration range scrutinized by this work. The evolution of resistive superconducting transitions with carrier concentration across n c1 is visible in Fig. 2.

The SdH data leading to an unambiguous identification of n c1 is presented in Fig. 3. The oscillating component of magnetoresistance, obtained after the subtraction of a smooth background (obtained ny a polynomial fit), is plotted as a function of the inverse of magnetic field times the main frequency. The sample with the lowest carrier concentration shows a single set of oscillations with split peaks. We attribute this splitting to the presence of tetragonal domains in our crystals. Recently, Allen et al. [14], by performing an angle-dependent Shubnikov-de Haas study, found that the Fermi surface of the lower band is larger along k z than along either k x or k y . In our study, the magnetic field is applied along the cubic [100] axis. Below 105 K, in the tetragonal phase, our samples can host three domains corresponding to the three possible orientations for the tetragonal z-axis. As illustrated in Fig. 3b, there would be two possible magnitudes for the cross section of the Fermi surface, when the magnetic field is along the nominal [100] axis. In this case, the main frequency corresponds to the smaller cross section, which is two times more likely to occur than the larger one. As discussed in detail in the supplement, this interpretation is in very good agreement with the known topology of the Fermi surface.

As seen in Fig. 3a, at the concentration of 1.2 10 18 cm -3 beating emerges and indicates the presence of a second set of oscillations with a significantly longer pe- FIG. 3: Left: A new frequency emerges in the Shubnikov-de Haas data in the vicinity of nc1. When the carrier density is below nc1, quantum oscillations display a single periodicity with split peaks and each peak is more prominent than its preceding neighbor(panel a). With increasing carrier density, a second frequency emerges on top of the first one (panels b to d). Data is plotted as function of the inverse of magnetic field normalized by the principal frequency, F1, in order to simplify comparison. The Fermi surface of the lower band according to ref. [14] is shown in panel e and its projection along the three high-symmetry axes in panels f to g. In multi-domain crystals, the two possible Fermi cross-sections give rise to split peaks in quantum oscillations.

riod. In samples with n H > n c1 = 1.2 × 10 18 cm -3 , the presence of this second frequency becomes clear and then it rapidly grows with increasing doping. The frequencies were quantified both by a FFT analysis of our data [See the supplement] and direct identification of the two periods in the data.

Analysis of quantum oscillations using the Lifshitz-Kosevitch formalism yields the magnitude of the cyclotron mass. The results of the analysis are shown in the upper panel of Fig. 4. We find that below n c1 , the effective cyclotron mass in the lower band, m 1 keeps a steady magnitude of 1.5 to 1.8m e , as previously reported [7]. At higher doping levels, m 1 becomes as large as 4m e . In this doping range, the middle band hosts significantly lighter carriers (m 2 ∼ 1.3m e ).

Comparing the total density of carriers distributed in the three bands and the carrier concentration obtained by measuring the Hall number, one can check the consistency of the emerging overall picture. Assuming three spherical Fermi surfaces leads to Fig. 4b. As seen in the figure, up to n c2 the two concentrations remain close to each other. The ratio of the two remains lower than unity, which is an expected consequence of the contrast between the assumed sphericity of the Fermi surface and its moderate anisotropy according to both experiment [13,14] and theory [9,11]. As discussed in the supplement, the measured anisotropy of the Fermi surface below n c1 [14] gives a satisfactory account of n SdH n H in our data. At the carrier density of 1.6 10 20 cm -3 , three bands are filled, but experiment detects only two frequencies (See Fig. 1a). It is straightforward to identify both the missing band and the reason of its absence. Only in the two upper bands the mobility is high enough for detectable quantum oscillations in a magnetic field of 17 T. The total density of carriers concentrated in these two visible bands is well below the total carrier concentration (See Fig. 4b). Both these features point to the lower band (which contains most of carriers, but those with less mobility) as the one which becomes invisible at this doping level. Now, we can address a central issue in the description of the band structure in this system. According to the picture first sketched by Mattheiss [9], the t 2g orbitals of titanium atoms give rise to a threefold degenerate band [10,11]. The degeneracy is lifted by a combination of spin-orbit coupling and tetragonal distortion. Recently, van der Marel and co-workers [11] pointed out that the energy gap between the bands and the peculiarity of their dispersion crucially depends on the choice of parameters. Our experiment finds that the effective mass in the lower band begins to increase as soon as the average radius of the Fermi surface exceeds 3 × 10 8 m -1 . As illustrated in panels d and e of Fig. 4, this feature is in good agreement with the choice of parameters in ref. [11], which leads to a non-parabolic dispersion for the lower band in absence of interaction. On the other hand the amplitude of the cyclotron mass and the size of the gaps differs significantly between non-interacting theory and experiment. The experimentally determined mass at the bottom of the band (1.5m e ) is twice as large as the theoretical expectation [11]. As a consequence of mass renormalization, commonly attributed to electronphonon interaction [19], the gap between the lower and middle bands (the Fermi energy of the lower band at n = n c1 ) is 1.8 meV, significantly lower than the expected value of 4.3 meV in the non-interacting band picture [11].

The main outcome of this study is to identify a critical doping, n c1 , below which a single band is occupied and T c rapidly rises with doping. In this regime, the T c /T F ratio is as large as 0.01, comparable to many unconventional superconductors [7]. This process is interrupted when carrier concentration exceeds n c1 . This highlights the specificity of multi-band superconductivity in SrTiO 3 , the subject of a recent theoretical work [20]. Motivated by the experimental reports on variation of T c near n c2 [8], Fernandes and co-workers considered the variation of superconducting transition temperature across a critical doping for different combinations of positive interband (λ ij ) and intraband (λ ii ) superconducting coupling parameters and found that filling an upper band (labeled 2) in addition to a lower band (labeled 1) would enhance T c . The larger the ratios of λ 12 /λ 11 and/or λ 22 /λ 11 , the more drastic is the expected enhancement. This contrasts with our experimental observation of a drop in T c (n) near n c1 .

This discrepancy calls for a theoretical reexamination extended to negative (i.e. repulsive) λ 22 and λ 12 . More importantly, in this analysis, λ ii has been taken to be a constant V i (representing overall attractive interaction between electrons) times N i (0), the density of states at the Fermi energy at band i. This assumption can only generate a monotonously increasing critical temperature, since the density of states unavoidably increases with car-rier concentration. The existence of domes implies that this approximation does not hold and points to a nontrivial evolution of V 1 as the outer Fermi surface expands.

Below n c1 , the normal state of superconducting SrTiO 3 , is a metal with a Fermi energy as low as 1 meV and thus the bosons exchanged between pairing electrons cannot be typical phonons with a Debye energy more than one order of magnitude larger [7]. One candidate to replace them at small wave-vectors is a plasmonic mode suggested by Takada [21]. Another one, invoked by Appel [22], is the soft mode observed below the 105K structural transition [23]. Our results frame the challenge faced by such propositions in a multi-gap context. The attractive interaction should couple electronic states more efficiently below n c1 , i.e. when the Fermi-surface is single component and the dispersion still parabolic.

In spite of its notorious inability to predict new superconductors [24], the BCS theory has been successfully employed to give an account of unexpectedly discovered superconductors such as MgB 2 [25]. In the case of SrTiO 3 , the challenge to theory is framed with exceptional sharpness and few ad hoc parameters. The combination of a simple Fermi surface topology, a very narrow energy window and a wide superconducting dome with fine details provides an exceptional opportunity to test the expected link between the critical temperature, the electronic density of states and the effective electronelectron interaction.

Table S1 lists the properties of the samples used in this study. The table indicates Hall carrier concentration (n H ), Hall mobility (µ H ), the superconducting critical temperature (T c , defined as the temperature at which the resistivity attains half of its normal-state value), oscillation frequencies (F i ) and cyclotron mass (m * i ) for each sample. Oscillation frequencies and cyclotron mass were extracted from the quantum oscillations of magnetoresistance or the Nernst coefficients.

Three different types of electron-doped SrTiO 3 single crystals were utilized in this study. The first category were oxygen deficient SrTiO 3-δ obtained by annealing stoichiometric SrTiO 3 samples in a temperature range of 700-1000 0 C in a vacuum of 10 -6 mbar for 1 to 2 hours [15? ]. The concentration of oxygen vacancies was controlled by finely tuning the annealing temperature. A second category of samples were commerciallybought Nb-doped SrTiO 3 single crystals. Four samples with well-defined Nb content (SrTi 1-x Nb x O 3 ; with x=0.001;0.002;0.01;0.2) were studied. In these four samples, the carrier concentration determined from the Hall coefficient matched the expected carrier concentration assuming that substitution of a titanium atom by a niobium atom introduces a single mobile electron.

A third category of the samples were made by reducing a Nb-doped sample with a nominal Nb concentration of x= 0.0004. Before reduction, the carrier concentration of this nominally SrTi 0.9996 Nb 0.0004 O 3 single crystal was orders of magnitude lower than the expected concentration of 8 × 10 18 cm -3 . Preparing reduced samples cut from this anomalously doped sample produced unexpected samples, which showed high mobility and displayed quantum oscillations of expected frequency, but their transition temperature was anomalously low.

In our analysis, we have focused only on the first two types of samples in which the origin of mobile electrons introduced by chemical substitution is clearly identified. Carriers are introduced either by removing an oxygen atom, or substituting a titanium atom by a niobium atom.

Appendix B: Shubnikov -de Haas oscillations Fig. S1 compares the mobility in our samples with those reported by Spinelli et al. [15]. To observe quantum oscillations electrons are required to be mobile enough to avoid a collision during a cyclotron period (implying µ > B -1 ) and also to avoid localization (k F l > 1, implying µ > e -1 k -2 F ). As seen in the figure, all samples studied here have mobilities high enough to meet these two requirements.

Fig. S2 shows the magnetoresistance (ρ) and its oscillation component(dρ) as a function of B -1 at T=0.55 K for different samples with different carrier concentrations. As seen in the figure, the amplitude of both the total magneto-resistance and its oscillating component decreases as the carrier concentration increases. The oscillation structure is simple and shows single periodicity in the upper panel. With the increase in doping, the oscillations become more complex.

Fig. S3 shows the normalized amplitudes of FFT in different samples. As seen clearly, there is a single frequency below n = 1.65 10 18 cm -3 , implying that only the lowest band is occupied. A second frequency of 8T emerges at n = 1.65 10 18 cm -3 , implying that the middle band starts to be filled. The frequencies become larger by electron doping, indicating an expanding Fermi surface.

1 E 1 7 1 E 1 8 1 E 1 9 1 E 2 0 1 E 2 1
Appendix C: Fermi surface distortion and the split peaks in multi-domain samples

Even at very low doping, near the bottom of the lower band, the Fermi surface of n-doped SrTiO 3 is not a perfect sphere [13]. Working on single-domain La-doped crystals of SrTiO 3 , Allen and co-workers quantified its departure from sphericity [14]. They found that the Fermi surface is squeezed along the c-axis such that the area of its projection on the (k x , k y ) plane is larger than its projection on the (k x , k z ) or the (k y , k z ) planes.

In a multi-domain crystal, the tetragonal distortion can occur along any of the three principal axes. Since the projection of the Fermi surface on the (k x , k z ) and (k y , k z ) planes are identical, with equal populations of the three domains, the main detected frequency would correspond to the projection of the Fermi surface along this direction. In this case, an additional frequency is expected corresponding to the projection of the Fermi surface on the (k x , k y ) plane. This gives a natural explanation for the presence of split peaks in our data. Fig. S4a compares the evolution of main frequencies according to our data and those reported by Allen and co-workers [14]. As seen in the figure, our main frequency is close to the one they detect when the field is along the smaller cross section of the Fermi surface. This confirms the interpretation proposed here.

A second quantitative confirmation comes from the magnitude of the tetragonal distortion determined by Allen and co-workers [14], who directly measured the frequency of oscillations for a field perpendicular to the large face of the Fermi surface (i.e. along k z ) and perpendicular to it. The ratio of these frequencies was found to be 1.5 (at n H = 3.6 × 10 17 cm -3 ) and 1.7 (at n H = 1.2×10 18 cm -3 ). This quantifies the tetragonal distortion of the Fermi surface. If we assume that the Fermi surface is a sphere, our measured frequency of quantum oscillations yields a carrier concentration close to, but systematically lower than, the Hall carrier density. The two set of numbers are given in Table S2. Approximating the Fermi surface to be a squeezed ellipsoid, the magnitude of tetragonal distortion, i.e. the ratio of the two perpendicular projections is the square of the two carrier densities. As seen in Fig. S4b, our data yields a magnitude for tetragonal distortion comparable to the value measured in ref. [14].

  FIG. 1: a) Detected frequencies of quantum oscillations as a function of carrier concentration. At two critical doping levels, designated as nc1 and nc2, a new frequency emerges. At each critical doping a new band starts to be occupied. b) Superconduting resistive transition temperature (on a logarithmic axis) as a function of carrier concentration. Solid squares represent reduced samples (SrTiO 3-δ ) and open squares Nbdoped samples (SrTi1-xNbxO3 with x=0.02; 0.01;0.002 and 0.001). Error bars represent the width of transition. c) Tc (on a linear axis) as a function of carrier per formula unit. Our data is compared with those reported by Schooley and coworkers[2]. The red shaded region shows the rough contours of superconductivity in STO/LAO interface[4, 16, 17].

FIG. 4 2 F

 42 FIG. 4: a) The effective mass in each band as a function of carrier concentration. The masses reported by an ARPES study on a bulk sample at n = 1.4 × 10 20 cm -3 [18] is also shown for comparison. b) Carrier concentration and its distribution among bands according to the detected frequency and assuming spherical Fermi surfaces. The total carrier density falls just below the Hall carrier density. Above 1.6 × 10 20 cm -3 , the lower band becomes invisible because of its reduced mobility. c) Lifshitz-Kosevitch analysis of the oscillations revealing two distinct masses for the lower and middle bands. Each color corresponds to a different oscillation. d) The Fermi energy, defined as 2 k 2 F 2mc in the two lower bands plotted as a function of corresponding kF . The middle band is shifted upward by the energy gap between the two bands. e) Theoretical energy dispersion near the bottom of the bands according to ref.[11].

  FIG.S1: Mobility of Nb-doped and reduced SrTiO3 samples as a function of carrier concentration. The horizontal black lines represent the mobility required to observe quantum oscillations in a given magnetic field. The blue lines represent the mobility threshold to avoid a metal-insulator transition.

3 FIG

 3 FIG. S2:The raw magneto-resistance (ρ) and residual resistance (dρ) as the function of the inverse field (1/B) at T = 550 mK and B // [0,0,1]. dρ is attained from ρ by subtracting the non-oscillating high temperature resistance curve.

  FIG. S3: The normalized Fast Fourier transform (FFT) of the oscillating component at T = 550mK for different samples. The dotted lines show the evolution of frequency of each band as the carrier concentration increases.

TABLE S1 :

 S1 n-doped SrTiO3 single crystals studied in this work.

	Sample	nH	µH	Tc	F1	F2	F3	m * 1	m * 2	m * 3
		cm -3	cm 2 .V -1 .s -1	mK	T	T	T	me	me	me
	SrTiO 3-δ	1.58 10 17	10395	a	8.7	-	-	1.5(0.1)	-	-
	SrTiO 3-δ	5.5 10 17	8941	86	18.2	-	-	1.83(0.07)	-	-
	SrTiO 3-δ	5.7 10 17	N. M.	N. M.	19.2	-	-	N. D.	-	-
	SrTiO 3-δ	6.79 10 17	13670	120.5	21.3	-	-	1.74(0.1)	-	-
	SrTiO 3-δ	7.67 10 17	7409	N. M.	23.5	-	-	N. D.	-	-
	SrTiO 3-δ	1.06 10 18	9949	175	26.5	-	-	1.7(0.3)	-	-
	SrTiO 3-δ	1.2 10 18	9385	190	26.8	-	-	1.8(0.3)	-	-
	SrTiO 3-δ	1.65 10 18	6366	206	31.4	7.1	-	2.3(0.4)	N. D.	-
	SrTiO 3-δ	1.93 10 18	8060	198	38.4	13.6	-	2.4(0.25)	1.7(0.25)	-
	SrTiO 3-δ	2.25 10 18	5868	202	N. M. N. M.	-	N. M.	N. M.	-
	SrTiO 3-δ	2.88 10 18	6078	194	49.8	20.4	-	3.48(0.15) 1.64(0.09)	-
	SrTiO 3-δ	3.88 10 18	5066	153	55.5	23.4	-	3.3(0.2)	1.65(0.15)	-
	SrTiO 3-δ	4.11 10 18	5421	178.5	63	28	-	4(0.5)	N. D.	-
	SrTi0.999Nb0.001O3	1.6 10 19	6290	242	138	69	-	4.2(0.2)	1.65(0.25)	-
	SrTi0.998Nb0.002O3	3.18 10 19	3121	338.6	208	96	7	3.5(0.13)	1.33(0.04)	N. D.
	SrTi0.99Nb0.01O3	1.6 10 20	870	429	N. D.	204	116	N. D.	1.57(0.08) 1.55(0.08)
	SrTi0.98Nb0.02O3	3.47 10 20	406	219	N. D. N. D. 243	N. D.	N. D.	N. D.
	SrTi0.9996Nb0.0004O 3-δ [b] 1.97 10 17	10600	N. M.	10	-	-	N. D.	-	-
	SrTi0.9996Nb0.0004O 3-δ [b] 1.06 10 18	18533	a	28	-	-	1.9 (0.2)	-	-
	SrTi0.9996Nb0.0004O 3-δ [b] 1.46 10 18	14444	a	30.7	-	-	1.9 (0.2)	-	-
	SrTi0.9996Nb0.0004O 3-δ [b]	2.5 10 18	14229	130	N. M. N. M.	-	N. M.	N. M.	-

a. Superconductivity was not detected down to 60 mK b. Nominal content (See text). c. N. D. refers to quantities which could not be determined. d. N. M. refers to quantities which were not measured. e. Irrelevant quantities are specified by -.

TABLE S2 :

 S2 Carrier concentration and the frequency detected below nc1.

	Sample SrTiO 3-δ	n1 = nH (cm -3 ) 1.58 10 17	F1(T ) 8.7	n2 = n SdH (Sphere)(cm -3 ) 1.45 10 17	( n 1 n 2 ) 2 1.19
	SrTiO 3-δ	5.5 10 17	18.2	4.38 10 17	1.58
	SrTiO 3-δ	5.7 10 17	19.2	4.74 10 17	1.45
	SrTiO 3-δ	6.79 10 17	21.3	5.54 10 17	1.50
	SrTiO 3-δ	7.67 10 17	23.5	6.42 10 17	1.42
	SrTiO 3-δ	1.06 10 18	26.5	7.69 10 17	1.90

We thank A. Balatsky, R. Fernandes, D. van der Marel, I. Mazin and J. -M. Triscone for discussions. This work is supported by Agence Nationale de la Recherche as part of QUANTHERM and SUPERFIELD projects and by Eu-roMagNET II under the EU contract number 228043. K. B. acknowledges the hospitality of Aspen Physics center.