R J Koval

Joshua M Pearce

A S Ferlauto

R W Collins

C R Wronski

Evolution of the Mobility Gap with Thickness in Hydrogen-Diluted Intrinsic Si:H Materials in the Phase Transition Region and Its Effect on p-i-n Solar Cell Characteristics

published or not. The documents may come

INTRODUCTION

Over the past decade, the deposition of a-Si:H i-layers from hydrogen diluted silane has been established as a means of improving the performance and stability of a-Si:H based p-i-n and n-i-p solar cells [START_REF] Bennett | Conf. Record of the 23 rd IEEE PVSC[END_REF][START_REF] Yang | [END_REF]. It is found that both the materials and the solar cells prepared with H 2 -dilution of SiH 4 exhibit improved optoelectronic properties and degradation kinetics distinctly different from their counterparts prepared without H 2 -dilution. They exhibit less degradation under 1 sun illumination and reach a degraded steady state in significantly shorter times [START_REF] Lee | Conf. Record of 25 th IEEE PVSC[END_REF]. The use of such cells as the wide bandgap component of multi-junction devices imposes requirements on their bandgaps and thicknesses brought about by the required matching of the currents generated in the component cells. Therefore there is great interest in using hydrogen dilution to vary the gap of Si:H intrinsic layers while at the same time optimizing the stabilized performance by adjusting their thicknesses.

The recently established evolution of the microstructural and optical properties with thickness for the "hydrogen diluted" materials [START_REF] Koh | [END_REF], coupled with the strong effect that the substrate material imposes on such properties for the overdeposited layers [5,6] greatly complicates the optimization of stable a-Si:H based solar cells. The complexity is further exacerbated in cases of cell structures optimized using a twostep i-layer process in which the p/i interface region has a higher hydrogen dilution ratio R (R≡[H 2 ]/[SiH 4 ]) than that of the bulk i-layer [START_REF] Koh | [END_REF][START_REF] Koval | Conf. Record of the 28 th IEEE PVSC[END_REF]. Furthermore, the inhomogeneity and anisotropic nature of these evolving materials give rise to difficulties in correlating thin film material properties with device characteristics due to the inherent differences between the transverse carrier transport in the cell configuration and the in-plane transport in the standard coplanar configuration [7]. The real-time in-situ studies of Si:H film growth which established an evolutionary phase diagram [8], together with the characterization and numerical modeling of the various cell structures described in this work, provides a unique method of circumventing these problems and allows reliable correlations to be made between the intrinsic material properties and corresponding cell characteristics.

EXPERIMENTAL DETAILS

The specular SnO 2 /p-a-SiC:H:B/i-Si:H/n-µc-Si:H:P/Cr solar cells were deposited in a multichamber rf-PECVD system at 200°C under previously described conditions [9]. To minimize the possibility of shunt currents in the dark J-V characteristics at low forward bias, small device areas were defined by evaporating 2 mm 2 Cr electrodes on the p-i-n structure and subsequently reactive ion etching the exposed n-µc-Si:H to isolate the cells. The dark J-V characteristics of the cells were measured over a temperature range from 25-125°C. The evolution of the microstructure in the various Si:H layers was characterized using a rotating compensator multichannel ellipsometer [10]. The cell characteristics were modeled using AMPS (Analysis of Microelectronic and Photonic Structures) simulation code [11]. The gap state distribution used for the a-Si:H materials includes charged defects and is similar to that used in the analysis of results in the detailed studies of a-Si:H films [12]. The amorphous-to-(mixed-phasemicrocrystalline) [a→(a+µc)] transitions were modeled as abrupt changes in the mobility gap, with the band discontinuities split equally between the valence and conduction bands [13], and midgap defect densities of 10 16 cm -3 for the µc-Si:H [14].

RESULTS AND DISCUSSION

Effects on performance and stability of the thickness-dependent a→(a+µc) phase transitions in the i-layers of Si:H based solar cells have been previously reported [9]. The initial fill factors (FF) of cells prepared with i-layers from 1100 Å to 1 µm, a thickness range over which such transitions occur, are shown in Fig. 1 for cells having 0≤R≤20 but otherwise identical deposition conditions. At a thickness of 1100 Å, the initial FF of all cells are virtually the same irrespective of R. For R=0, where no phase transition occurs over the entire thickness range of the i-layer, the FF decreases as expected in a manner closely matching that of the numerical simulation for a uniform homogeneous i-layer. However, for the cells prepared with dilution, there are clearly discernible thickness-dependent deviations from the R=0 FFs that can be attributed to the phase transitions in the bulk layers as indicated by the phase diagram developed for Si:H film deposition on amorphous film surfaces [8]. It is important to note here that the phase transition in the Si:H materials is not abrupt, but passes through an amorphous + microcrystalline (a+µc) mixed-phase region, in which the microcrystallites embedded in the amorphous matrix increase in in-plane scale with increasing thickness at a rate that depends on R.

The mobility gap (E µ ) of the bulk i-layers and their transition regions plays a key role in transport and collection of photogenerated carriers. In p-i-n cells whose characteristics are dominated by the bulk, the dark current characteristics are very sensitive to the i-layer E µ . This is illustrated in Fig. 2a with the forward dark J-V characteristics of p-i-n cells having 4000 Å thick R=0 and R=10 i-layers (shown as symbols). The J-V characteristics exhibit the classical exponential diode behavior with a reverse saturation current density (J o ) that can be obtained from the intercept of the linear region of this semilog plot. The R=0 and R=10 materials have the same densities of midgap states (D o ) and E µ values of 1.80 eV and 1.87 eV, respectively, as determined from internal photoemission measurements [15]. The sensitivity of the J-V characteristics to E µ is reflected in the values of J o , which differ by approximately one order of magnitude, and the excellent agreement with simulation results (shown as solid lines) obtained using the corresponding values of E µ . It is found that reliable E µ values can be obtained from the dependence of J o on temperature in p-i-n cells where the dark currents are dominated by recombination-generation (R-G) in the bulk i-layer [16,17,18]. In this case, the dependence of J o on the i-layer mobility gap is:

        - = kT E A T J o 2 exp 3 µ (1)
where A is a temperature-independent pre-factor and E µ is the i-layer mobility gap. The dependences of J o /T 3 on temperature for the cells in Fig. 2a are shown in Fig. 2b, where the well-defined slopes yield values of E µ that are in excellent agreement with the values quoted above.

The correlations among carrier recombination kinetics in the bulk i-layer, J o , and E µ established for p-i-n cells with homogeneous i-layers have been extended to include the more complicated case of cells with evolving i-layers. This effort was carried out by studying a series of p-i-n cells each with a total i-layer thickness of 4000 Å. An R=40 layer was incorporated at the p/i interface of each cell, leading to a two-step i-layer with R=40 layer thicknesses ranging from 100 Å to 400 Å. The cells exhibit optimum performance when the R=40 interface layer thickness is 200 Å, whereas with increasing thickness above this optimum both FF and V oc decrease [START_REF] Koval | Conf. Record of the 28 th IEEE PVSC[END_REF]. Insights into the electronic properties of the transition material from amorphous to mixed-phase Si:H and then to single -phase µc-Si:H that give rise to deleterious effects on cell performance can be obtained from the dark J-V characteristics of the cells as shown in Fig. 3a. The increase in carrier R-G in the vicinity of the phase transition region can be directly inferred from the systematic increase of the dark currents in the exponential region as the R=40 Si:H layer thickness increases from 100 to 400 Å and traverses the phase boundary. The increase by two orders of magnitude in the room temperature value of J o clearly indicates a large increase in R-G currents that can be attributed primarily to changes in the mobility gap. This is evident from Fig. 3b wherein plots of J o /T 3 versus (2kT) -1 are shown for the results in Fig. 3a. Because the value of 1.88 eV obtained from the cell with the 200 Å R=40 interface layer matches that for the bulk homogeneous i-layer (see Fig. 2b), any decreases in activation energies from this value with increasing interface layer thickness can be used to obtain the changes in the effective mobility gap that occur across the phase transition. The decreases to 1.66 eV for the 300 Å interface layer and to 1.22 eV for the 400 Å layer indicate that over a thickness of 200 Å there is a transition from the amorphous phase to a material having a mobility very close to that of bulk microcrystalline Si:H. However, the fact that the value of J o increases by only two orders of magnitude across the phase transition is a clear indication that the microcrystalline phase does not extend throughout the i-layer but that the subsequent growth of the R=10 material proceeds in the amorphous phase as has been substantiated by RTSE studies [8].

To obtain additional insights into the transitions as well as their effect on cell characteristics, studies were carried out on a cell in which a 400 Å R=40 p/i interface (similar to that described above) was followed by an R=20 bulk i-layer region (referred to here as cell C). The characteristics of this cell are compared with those of the optimized structure described above having a 200 Å R=40 p/i interface and an R=10 bulk region (referred to as cell A). The comparison also includes a cell deposited with a one-step R=20 i-layer (referred to as cell B). While cell A has a relatively homogeneous amorphous ilayer, ellipsometry studies show that the a→(a+µc) transition in the R=20 layer deposited on an amorphous film occurs after an accumulated film thickness of approximately 1000 Å. Simila r studies have also demonstrated that the growth of the same R=20 layer on an R=40 layer that has traversed the phase boundary proceeds in the microcrystalline phase so that the i-layer is fully microcrystalline with the exception of the initial 400 Å. In fact, the first 200 Å of the i-layer of cell C is amorphous whereas the next 200 Å is mixed-phase (a+µc). The experimental J-V results at room temperature are shown in Figs. 4a and4b under 1 sun illumination and in the dark, respectively, for cells A, B, and C. It can be seen in Fig. 4a that the FF and V oc are the highest for cell A since the thickness and R values in the two step ilayer are chosen so that protocrystallinity [i.e., Si:H deposition just on the amorphous side of the a→(a+µc) phase boundary] is maintained. The cell performance becomes progressively worse for cells B and C as the position of the phase boundary in the i-layer approaches the p-contact (indicated by the direction of the arrows).

This drop in cell performance is attributed to narrow gap material in the i-layer of cells B and C whose presence is reflected in the dark J-V characteristics as shown in Fig. 4b. In p-i-n structures most of the carrier recombination occurs in a narrow region, approximately halfway between the p-and ncontacts at low forward bias. At higher forward bias, this region extends toward both contacts. Fig. 5 shows such a bias-dependence and spatial non-uniformity of the carrier recombination rates in a p-i-n cell with a homogeneous amorphous i-layer. Inflection points in the dark J-V characteristics observed in Fig. 4b can thus be explained by the large contribution of the narrow gap material to the recombination currents at low forward bias. The eventual merging of the currents in cells B and C with that of cell A at higher forward bias in Fig. 4b is due to a shift of the recombination from the bulk i-layer to the vicinity of the contact layer regions that are similar for all three cells. The forward bias voltage at which the dark current of cell C merges with that of cell A is larger than the corresponding bias for cell B. This behavior occurs because the phase boundary in cell C is closer to the p-layer so that the carrier recombination in the narrow gap material dominates over a larger voltage range. The simulation result shown in Fig. 4b (solid line) was carried out for cell C assuming a 400 Å wide gap Si:H layer with the remainder of the ilayer (3600 Å) as a narrow gap µc-Si:H material having a mobility gap of 1.15 eV. The latter value is obtained experimentally from plot of J o /T 3 vs. (2kT) -1 for cell C. The use of an abrupt change in E µ is justified by the data shown in Fig. 3b as well as by RTSE. Both results suggest that the transition from the amorphous phase to a fully coalesced microcrystalline phase in the R=40 film occurs over a thickness of only a few hundred angstroms. The large decreases of the FF in cells B and C are a direct consequence of the phase transitions. The higher photo-carrier densities generated in the µc-Si:H portion of the i-layer lead to a collapse of the electric field in that region in order to maintain current continuity across the entire i-layer. As a consequence, most of the voltage drop then occurs across the amorphous region. The weaker average field across the i-layer results in the anomalous inflection point in the light J-V characteristics as seen in Fig. 4a, with the effect becoming more pronounced as the phase transition approaches the p-contact. Such characteristics could only be modeled for i-layers consisting of two distinctly different layersa widegap amorphous layer and a narrow gap microcrystalline one. Such an approach was also successful in simulating the results of earlier studies of intrinsic µc-Si:H growth on p-a-SiC:H [START_REF] Siebke | [END_REF]. The results of the simulation for cell C is shown as a solid line in Fig. 4a, where the less than perfect fit is probably due to the necessary oversimplification of the highly anisotropic microcrystalline material when a onedimensional simulator such as AMPS is being used.

CONCLUSIONS

Correlating a-Si:H material properties with cell characteristics has been an ongoing challenge. The evolution of the growth of the Si:H films prepared with H 2 -dilution from an amorphous phase to the highly anisotropic microcrystalline phase makes the task even more difficult. Hence cell optimization requires a systematic approach such as that described here. In this unique approach, p-i-n cells were selfconsistently modeled assuming a change in the i-layer mobility gap at the thickness determined by RTSE, whereby the values of the mobility gaps are obtained from measurements on the actual p-i-n structures. As a result, meaningful correlations can be made between the materials and solar cell characteristic s since the uncertainties associated with the initial nucleation and subsequent growth of Si:H films on different substrate materials are eliminated. Furthermore, the simulations corroborate the thickness position of the phase boundary as determined by RTSE. The insights provided here into the transition from protocrystalline Si:H to mixed-phase Si:H and then to single -phase µc-Si:H are of critical importance to the engineering of cell structures with high performance and improved stability.

Figure 1 .

 1 Figure 1. Dependence of initial fill factor of p-i-n cells on the intrinsic layer thickness for R=0, 10, and 20. For comparison, simulation results for an R=0 structure having a homogeneous i-layer are also included.

Figure 2 .

 2 Figure 2. (a) Experimental dark J-V characteristics (symbols) of R=0 and R=10 p-i-n cells having 4000 Å thick i-layers and equivalent contact layers. AMPS simulation results are shown as solid lines. (b) Thermal activation plots of J o /T 3 for the cells shown in (a). The symbols are experimental results and the solid lines are exponential fits used to extract E µ .

Figure 3 .

 3 Figure 3. (a) Experimental dark J-V characteristics of 4000Å p-i-n cells having R=10 bulk i-layers and different thickness R=40 p/i interface layers. (b) Thermal activation plots of J o /T 3 for the cells shown in (a). The symbols are experimental results and the solid lines are exponential fits used to extract E µ .

Figure 5 .

 5 Figure 5. Recombination profiles in a typical a-Si:H p-i-n solar cell in the dark having a homogeneous i-layer. Different values of forward bias voltage are shown (V=0.1, 0.2, 0.4, 0.8, 1.0 V) with the direction of the arrow indicating increasing forward bias voltage.

Figure 4 .

 4 Figure 4. Experimental data (symbols) and simulations (solid lines) for the J-V characteristics of cells A, B, and C measured at room temperature (a) under 1 sun illumination and (b) in the dark. Direction of the arrows indicate decreasing distance between the a→(a+µc) transition and p-contact.


J M Pearce

R J Koval

R W Collins

C R Wronski

Quantitative Correlation of High Quality a-Si:H p-i-n Solar Cell Characteristics with Properties of the Bulk and p/i Interface Region

Studies have been carried out on high quality hydrogenated amorphous silicon (a-Si:H) p-i-n solar cells with protocrystalline i-layers to establish the nature of p/i interfaces and to quantify their contributions to various solar cell characteristics. The p-a-SiC:H:B/i-a-Si:H/n-µc-Si:H:P cell structures used had the a-Si:H i-layers deposited from hydrogen diluted silane with R≡[H2]/[SiH4]=10. The high quality p/i interface regions obtained with R=10, indicated by the high and stable open circuit voltage (Voc) values, were further improved by increasing R in the 200Å of a-Si:H adjacent to the a-SiC:H layer. From the systematic improvement and ability to obtain p/i interface regions with outstanding quality, it was possible to track their contributions to cell characteristics relative to those from the bulk. Results of dark current voltage (JD-V) and short circuit current-open circuit voltage (Jsc-Voc) characteristics are presented which clearly demonstrate that even high quality interface regions in p-in cells can mask some contributions of protocrystalline bulk layers. Results are also presented and discussed on how the relative contributions of bulk and p/i interfaces can be isolated and quantified so they can be used as inputs for reliable analysis of solar cell characteristics.

INTRODUCTION

Open circuit voltage (Voc) is one of the central factors that determine the overall performance of amorphous silicon solar cells. There has been significant progress in increasing the open circuit voltages in both p-i-n and n-i-p a-Si:H solar cells, however, the physical phenomena that underlie their magnitude and stability are not fully understood. The open circuit voltage can also be limited by intrinsic layer recombination. Thus, there have been a large number of attempts to relate the properties of the intrinsic materials in a-Si:H p-i-n and n-i-p solar cells to the magnitude and stability of their Voc [1,2]. Unfortunately, the success of these attempts have been incomplete due predominantly to the inability to quantify the individual contributions of the bulk i-layer and those from the p/i interface region to the cell characteristics. This is partially due to the fact that in such studies both the p/i interface regions and the bulk i layer were changed simultaneously. Such, absence of reliably quantified information about the individual contributions of the bulk and interface seriously limits the ability of systematically (rather than empirically) improving the Voc. In the work reported here the contributions of the same protocrystalline i-layers [3,4] with different p/i interface regions to the dark J-V and Jsc-Voc characteristics were quantified. Because of the high quality of the intrinsic material, the effect of good p/i interface regions can be identified on both the JD-V and Voc characteristics. This allows the contributions to the Voc of different p/i interface regions as well as that of the bulk to be quantified for cell structures that exhibit both high Vocs and high Fill Factors (0.72 for 4000 Å i layers).

Even though large improvements have been made in the quality of both the i-layers and the p/i interfaces, a better understanding of how they affect cell characteristics is still required before reliable analyses can be carried out of cell performance and stability. This study establishes and quantifies the contributions of p/i interface regions to various cell characteristics in high quality p-i-n a-Si:H solar cells with protocrystalline ilayers.

EXPERIMENTAL

Baseline, state-of-the-art solar cells were fabricated with intrinsic protocrystalline Si:H by RF plasma enhanced chemical vapor deposition (PECVD) with SiH4 diluted with hydrogen with a dilution ration R≡[H2]/[SiH4] = 10 under conditions previously described [4]. The cells were fabricated on specular SnO2 consisting of 250Å a-SiC:H:B(p)/ 4000Å a-Si:H (i)/350Å μc-Si:H:P(n)/ 1000Å Cr. The high quality of the protocrystalline i-layer, deposited at 200˚C has been clearly established from detailed studies of cells with different thicknesses [5] and corresponding thin film materials [6]. In these p-i-n cells 200 Å p/i interface regions were systematically changed by increasing their R from 10 to 40 [START_REF] Lee | Conf. Record of 27 th IEEE PVSC[END_REF]. To minimize possible contributions of shunts to the dark J-V's at low forward bias these cells had small areas (0.02 cm 2 ) isolated by reactive ion etching of the n μc-Si:H layers. The JD-V and Jsc-Voc characteristics were measured at 25˚C with the latter being obtained with illuminations between ~10

-7 and 1 sun from an Oriel solar simulator. Results were obtained on cells from the annealed state (4 hours at 170˚C) to the 1 sun degraded steady state (DSS) obtained after 100 hours of illumination [START_REF] Lee | Conf. Record of the 25 th IEEE PVSC[END_REF].

RESULTS

The p/i interfaces in the baseline structures were of high quality, as reflected by the high values of fill factor (FF) and Voc and the liner dependence of Jsc on intensity over the entire range of illumination. However, it was found that the Voc could be further improved by incorporating a thin (100-200Å) material obtained with higher H2 dilutions in the interface between the p and i regions [START_REF] Lee | Proc. 2 nd World Conf. on Photovoltaic SolarEnergy Conversion[END_REF]. Systematic improvements were obtained in the p/i interface regions, which resulted in even higher stabilized Voc's for higher diluted interfaces. For a 200Å thick interface R=40 material resulted in the highest Voc as seen in Table 1. If the dilution in the interface layer is th IEEE Photovoltaic Specialists Conference, Alaska 2000. increased further the material undergoes the amorphous to microcrystalline phase transition [3] and the Voc of the cell declines. From the systematic improvements in the interface regions it was possible to track their contributions to cell characteristics relative to those from the bulk. It is well known that increases in Voc of p-i-n or n-i-p solar cells can be obtained with bulk i layers having wider bandgaps. Unfortunately, the wider bandgap materials show less optical absorption of the solar spectrum and thus cause a reduction in the short circuit current density (Jsc), so they are unsuitable for use in the bulk intrinsic layer of high efficiency single junction cells. The absence of limitations due to the built in potential for cells studied here where established with results of Voc at low temperatures and at intensities greater than one sun [START_REF] Lee | Conf. Record of 27 th IEEE PVSC[END_REF].

p/i interface R=[H2]/[

The effect of the improved interfaces on the annealed state dark J-V characteristic is illustrated in Figure 1, where results for the baseline and optimized (R=40 p/i interface) are shown for the annealed state. It can be seen in Figure 1, that introducing the optimized p/i. interface region has a profound effect on the dark J-V characteristics, where the exponential currents for voltages between 0.4V and 1V are reduced up to an order of magnitude.

Since these currents are associated with the generation and recombination of carriers this decrease can be attributed to a significant reduction of recombination in the p/i interface region, which originally partially obscures the contributions of the bulk. It is important to note here that the extent to which p/i interfaces can dominate the dark J-V characteristics has not been fully taken into account in the numerous modeling results of annealed p-i-n characteristics. In such analysis arbitrary assumptions are made about the possible contributions of the p/i interface relative to the bulk in order to ascertain their properties. Such importance of the p/i interface in determining the dark J-V characteristics is acutely evident in these high quality protocrystalline i-layers and cannot be ignored particularly in the annealed state. It is expected that as the quality of the bulk declines the contributions of the p/i interface to the dark J-V characteristic will be reduced. This is indeed the case and is illustrated in Figure 2, where dark J-V characteristics for the cells in Figure 1 are shown in the degraded steady state (DSS) obtained with 1 sun illumination at 25°C. Now, the dark J-V's are nearly identical since the generation-recombination currents are in this case dominated by the light induced defects in the bulk even for the baseline cell.

It is well known that Voc, and in particular its stability, are very sensitive to the properties of the p layers and the p/i interface regions. Nevertheless, it still depends on the properties of the bulk intrinsic layers and in particular their bandgap. We find that it is possible to separate the contributions for the bulk and the p/i interface region from Jsc-Voc characteristics obtained over a wide range of illumination intensities. This is illustrated in Figure 3, where the Jsc-Voc characteristics for the baseline and optimized cells in the annealed state are shown for illumination levels between 1 and 10 -7 suns. The solid line is a guide to the eye indicating the intensities over which the Jsc-Voc characteristics for the two cell overlap. It can be seen in Figure 3 that, even though there is a significant difference between the Voc's for the cells with and without the R=40 interface region at 1 sun illumination, at intensities less than about 10 -4 suns the Voc's are identical. This indicates that in this regime the Voc in both cells is determined by the properties of the same protocrystalline bulk layer, whereas, at the higher intensities the contributions of the different p/i interface regions can be identified. It is interesting to note here that both the Jsc-Voc and the JD-V (Figure 1) characteristics overlap up to ~0.4V. Such a contribution of the bulk is further confirmed by the results obtained after the two cells have reached a degraded steady state (DSS) under one sun illumination. It is expected that as the quality of the bulk declines the region over which the contribution of the p/i interface can bee seen should decrease, and in the case of poor quality material could disappear completely. The light induced degradation to DSS reduces the quality of the bulk by introducing about 10 17 cm -3 defects, where these defects increase the generation-recombination currents about by a about a factor of ten.

The annealed state and DSS for the baseline R=10 cell are shown in Figure 4. Although there is an absence of any major change in the Voc with 1 sun illumination, at lower illumination intensities there is a marked difference between the two characteristics. Due to the light induced degradation of the i-layer the region over which the bulk dominates is extended to intensities of 10 -3 suns for this baseline cell. The same phenomenon is seen in Figure 5 for the cell with an R=40 p/i interface region. In this case, however, there is a small drop in the 1 sun Voc at 1 sun illumination. The light induced degradation in the i layer extends the bulk domination to 1 sun for this optimized cell. It should be noted that at higher intensities (~50 suns) the annealed and DSS characteristics for this cell also merge together, similar to the baseline case, and there is no decrease in Voc.

DISCUSSION

The results presented and discussed above show how the bulk and p/i interface contributions to dark J-V and Jsc-Voc characteristics can be identified and quantified. The extended exponential dependence of the currents on voltage seen in both the dark J-V and Jsc-Voc characteristics can be expressed by a classical diode equation as:

J= JO[exp(qV/mkT)-1]

where Jo is an effective reverse saturation current density and m is an operational diode quality factor over that region. Since the values of m are determined by the voltage dependence of recombination of injected carriers such operational diode quality factors have several components. This include the effect of i layer thickness, and the density and distributions of recombination centers and in addition to this bulk component there is one due to recombination at the interface. A seen in Figure 3, in the th IEEE Photovoltaic Specialists Conference, Alaska 2000. annealed state at low intensities both the baseline cell and the one with the optimized p/i interface cell have m = 1.6. In the case of the baseline cell, however, a value of m = 1.3 occurs at Vocs greater than 0.4 V, where in the case of the R=40 p/i interface cell this does not occur until ~0.7V. The decrease in m value can be attributed to recombination that is dominated by the interfaces for these cells. It can be seen in Figure 5 that for the R=40 p/i interface cell in the degraded steady state, a single value of m now extends to 1 sun illumination indicating that there is a single contribution to the recombination processes. This implies that the recombination due to the higher density of defects, introduced by light, now dominates over that at the p/i interface all the way to Voc for 1 sun illumination. Separating and quantifying such m values in jsc-Voc characteristics offers important information about the nature of the bulk and p/i interfaces. However, it should be noted that neither the p/i interface nor the bulk values of m are unique for all high quality a-Si:H. Hence, it is found that in high performance solar cells the recombination and Voc need not be limited by the p/i interface. The results presented here also clearly indicate that estimates for the maximum 1 sun Voc and the built in potentials cannot be estimated given results on cells whose Jsc-Voc characteristics are dominated by generation-recombination at the p/i interface [2]. Detailed results of the role of carrier recombination on these characteristics as well as their effect on limiting the maximum Voc that can be obtained with given bulk layers will be discussed elsewhere.

The separation of bulk and p/i interface contributions to the characteristics discussed here offer meaningful inputs into self consistent analysis of solar cell characteristics which does not rely on arbitrary and ad hoc assumptions about the nature of the p/i interface. In addition, the results illustrate that Voc's at 1 sun can be stable when they are limited by the p/i interface regions and not the bulk i layers and that any valid analysis must be based on results over a wide range of intensities. It is also imperative in any analysis to ensure that the modeling be carried out first on cells whose characteristics are indeed dominated by the bulk and only then can the contributions of various p/i interfaces subsequently be realistically analyzed. Such analysis has been undertaken on the characteristics, which have been shown to be bulk-dominated and self consistency is being established based on "operational" gap state parameters obtained from Schottky barrier [START_REF] Lu | [END_REF] and thin film studies [11].

Fig. 1 .

 1 Fig.1. The annealed state dark J-V characteristics of for an R=10 cell and R=10 cell with an R=40 p/i interface region.

Fig. 2 .

 2 Fig.2. The 1 sun degraded steady state dark J-V characteristics of an R=10 cell and an R=10 cell with an R=40 p/i interface region.

Fig. 3 .

 3 Fig.3. The Voc-Jsc (light intensity) characteristics of an R=10 cell and an R=10 cell with an R=40 p/i interface region in the annealed state.

Fig. 4 .

 4 Fig. 4. The Voc-Jsc (light intensity) characteristics of an R=10 cell in the annealed and degraded steady states

Fig. 5 .

 5 Fig.5. The Voc-Jsc (light intensity) characteristics of the cell with an R=40 p/i interface region in the annealed and 1 sun degraded steady state.

th IEEE Photovoltaic Specialists Conference, Alaska 2000.

ACKNOWLEDGEMENTS

We would like to acknowledge the support by the National Renewable Energy Laboratory (NREL) under subcontracts HXAF-8-17619-22 and the New Energy Development Organization (NEDO) International Joint Research Grant. REFERNCES


E Sarrouy
email: emmanuelle.sarrouy@centrale-marseille.fr

Phase driven study for stochastic linear multi-dofs dynamic response

Keywords: Structural dynamics, Uncertainty, Polynomial Chaos, Frequency response function

This work addresses the computation of dynamic responses of stochastic linear systems using polynomial chaos expansion. As is now well known, polynomial chaos does not offer an accurate representation of dynamic response around resonances when the responses are evaluated for several frequency values. A new parametrization of the frequency response function is then proposed: instead of considering the frequency as the main parameter, a "total phase" parameter is defined and used to define the dynamical system to be solved. It is shown via two applications that this approach offers very accurate results when conjugated to polynomial chaos with low degree.

Introduction

Evaluation of the frequency response function (frf) is one of the main steps of a study in the frame of structural dynamics. In the case of deterministic linear systems, this presents no major difficulty except for the computational cost when processing large systems. However, when uncertainties are introduced the robust Monte Carlo Simulation (MCS) method leads to prohibitive computational costs. Hence, several works propose to overcome this problem by using different approaches: building surrogate mode-based models [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF], Kriging [START_REF] Wang | Frequency response function-based model updating using Kriging model[END_REF] or polynomial chaos for example. In this latter case, multiple studies showed that for stochastic systems, the classical generalized Polynomial Chaos (gPC, [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF]) fails to offer an accurate representation of the stochastic dynamic response in the vicinity of resonances, that is in frequency ranges where accuracy is especially required. This phenomenon is emphasized over a single degree of freedom (dof) in [START_REF] Pagnacco | Polynomial chaos for modeling multimodal dynamical systems -Investigations on a single degree of freedom system[END_REF] and in the frame of rotor dynamics in [START_REF] Sinou | Influence of Polynomial Chaos expansion order on an uncertain asymmetric rotor system response[END_REF]. It is analysed in [START_REF] Jacquelin | Polynomial Chaos Expansion and Steady-State Response of a Class of Random Dynamical Systems[END_REF] as "PC-resonances'" that is resonances due to Polynomial Chaos Expansion (PCE) which are located close to the dynamical system resonances.

To overcome this problem, Jacquelin et al. proposed to use an accelerating scheme to get faster convergence of PCE for the first two moments [START_REF] Jacquelin | Polynomial chaos expansion in structural dynamics: accelerating the convergence of the first two statistical moment sequences[END_REF]; Multi-Element generalized Polynomial Chaos method [START_REF] Wan | Multi-element generalized polynomial chaos for arbitrary probability measures[END_REF] can also be applied successfully as tested in [START_REF] Pagnacco | Polynomial chaos for modeling multimodal dynamical systems -Investigations on a single degree of freedom system[END_REF]. Another possibility is to use another parametrization of the frf as proposed in [START_REF] Sarrouy | A constant phase approach for the frequency response of stochastic linear oscillators[END_REF] for a single dof system: instead of computing the dynamic stochastic responses for different frequencies, the stochastic responses are evaluated for different response phase values; this leads to a simpler dependency of the response with respect to the random variables and hence a more accurate representation using low degree polynomials. This work is a generalization of this phase parametrization approach to linear systems with multiple dofs. It is close to a recent work by Yaghoubi et al. who propose a piecewise linear transform of frequency intervals between several selected frequencies (resonances and anti-resonance) [START_REF] Yaghoubi | Sparse polynomial chaos expansions of frequency response functions using stochastic frequency transformation[END_REF]. This work proposes an analytic function instead of a two steps process (frequencies selection and piecewise linear transformations). Considering multiple dofs, a "total phase" parameter is defined and the stochastic dynamical problem formulation is redefined accordingly. The idea is that for a same value of the scalar "total phase" different realizations will have close responses which will be accurately represented using low degree polynomials whereas for a same frequency value, the responses from one realization to another can be very different and hence require higher degree polynomials.

The theoretical steps are detailed in Sec. 2. Sec. 3 is devoted to computational considerations and defines an approximate expression of the "total phase" which can save time for large systems. Finally, the classical approach and the proposed one are compared via two numerical applications in Sec. 4.

Problem statement and theoretical points

Let's consider a general linear damped multi-dofs system in its classical formulation:

M ü + C u + Ku = f (t) (1)
This is a second order linear ordinary differential equation with constant coefficients; let us denote its size n. It is assumed throughout the paper that the mass and stiffness matrices M and K are symmetric definite positive and that the damping matrix C is such that it reflects a sub-critical damping. When working in the field of structural dynamics, one usually wants to establish the frf of system [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF]. The excitation f (t) is assumed to be periodic with a single frequency component:

f (t) = f c cos(ωt) + f s sin(ωt), (f c , f s ) ∈ (R n ) 2 (2)
The stationary response u is then sought in a similar form:

u(t) = a cos(ωt) + b sin(ωt), (a, b) ∈ (R n ) 2 (3)
Components a and b can be evaluated using a complex reformulation of equation [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF] as described in Appendix A:

K(ω)û = f with K(ω) = K -ω 2 M + ωC (4)
where

 2 = -1, û = a -b ∈ C n and f = f c -f s ∈ C n .
To build a frf, û is evaluated for a given range of circular frequency ω values; to be more specific, it is evaluated for N discrete values ω (j) , 1 ≤ j ≤ N over the range of interest. It is assumed here that the system is damped so that the dynamical flexibility matrix Ŝ(ω) = K(ω) -1 always exists (there are no undamped modes). Hence, for each circular frequency value ω there is one and only one stationary solution expressed in its complex description û satisfying (1). Let's now consider uncertainties that are modelled in the probabilistic frame: random variables or random fields are used to render observed or assumed variations of the structure; the excitation vector f is considered free of any uncertainty. In the case of random fields, it is assumed here that they are approximated by a collection of uncorrelated random variables via a Karhunen-Loève decomposition [START_REF] Loève | Probability Theory II[END_REF][START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] for example. Hence uncertainties are eventually rendered via a set of ñ random variables, denoted ξ i , 1 ≤ i ≤ ñ. This leads to random matrices and system response which are denoted as follows:

M ü + C u + K ũ = f (t) (5)
To establish the frf in the stochastic context, the deterministic approach is usually expanded: for several values of ω within the frequency range of interest, ũ is evaluated, most frequently in its complex description ũ as exposed above.

The system which is solved is then:

K(ω) ũ = f with K(ω) = K -ω 2 M + ω C (6)
A way to express ũ is to expand it into a truncated polynomial series: this is what gPC does [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF][START_REF] Lucor | Adaptive Generalized Polynomial Chaos for Nonlinear Random Oscillators[END_REF][START_REF] Ernst | On the convergence of generalized polynomial chaos expansions[END_REF]. However, it was demonstrated that this representation is non accurate in the vicinity of resonances which are the frequency ranges of interest [START_REF] Pagnacco | Polynomial chaos for modeling multimodal dynamical systems -Investigations on a single degree of freedom system[END_REF][START_REF] Sinou | Influence of Polynomial Chaos expansion order on an uncertain asymmetric rotor system response[END_REF].

The idea developed in this work to overcome such a problem is to read Eq. ( 4) as n relationships between n + 1 variables ûk , 1 ≤ k ≤ n (components of û) and ω. The classical way to process this dynamical system is formalized as the addition of an (n + 1)-th equation C(û, ω) = 0 where C(û, ω) = ω -ω (j) with ω (j) a given circular frequency. We propose to modify this latest equation so as to draw smoother surfaces considering ξ i variations in the random case, and hence getting a better convergence rate of the gPC for this new problem. First, this new constraint equation is defined based on the modal properties of system [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF]. Then, its combination with gPC is exposed. A last subsection explains how to get classical -for constant frequencies -results from results of the proposed approach.

Quadratic eigenvalue problem and total phase definition

Let's consider the quadratic eigenvalue problem (QEP) associated with Eq. ( 1):

(λ 2 k M + λ k C + K)ϕ k = 0 (7)
where λ k and ϕ k , 1 ≤ k ≤ 2n, denote the eigenvalues and eigenvectors respectively.

When there is no damping (C = 0), the spectral theorem provides the existence of a basis of real eigenvectors that simultaneously diagonalize M and K matrices. For classically damped systems, that is systems with a damping matrix verifying the Caughey condition [15]

CM -1 K = KM -1 C (8)
the set of real eigenvectors obtained for the undamped system also diagonalize C. The eigenvalue problem of non-classically damped systems is not as simple as the two cases mentioned above [START_REF] Tisseur | The Quadratic Eigenvalue Problem[END_REF][START_REF] Veselić | Damped Oscillations of Linear Systems: A Mathematical Introduction[END_REF]. In this latter case, there does not always exist a basis (even a complex one) which simultaneously diagonalize mass, damping and stiffness matrices [START_REF] Ma | The decoupling of damped linear systems in oscillatory free vibration[END_REF].

It is assumed in the rest of the paper that the QEP (7) has 2n distinct eigenvalues which are complex (not purely real nor purely imaginary): each mode is damped, none is critically or overdamped. In the case of complex eigenvalues, as matrices are real, if (λ k , ϕ k ) is a solution, then (λ k , ϕ k ) is also a solution: eigenvalues come in complex pairs. Each eigenvalue λ k can be split into real and imaginary contributions:

λ k = -α k + ω k (9)
with α k > 0 and  2 = -1. One chooses to index λ k = -α k + ω k with positive and increasing ω k contribution from 1 to n and their conjugate from n + 1 to 2n. Hence, ω k > 0 and

λ k+n = λ k , 1 ≤ k ≤ n.
This said, one can consider the characteristic polynomial q(λ) associated to the linear operator in [START_REF] Jacquelin | Polynomial chaos expansion in structural dynamics: accelerating the convergence of the first two statistical moment sequences[END_REF] [16, p. 250]:

q(λ) = det(λ 2 M + λC + K) = det(M ) n k=1 (λ -λ k )(λ -λ k ) (10)
Finally, the total phase is defined as follows:

φ tot (ω) := -∠(det M ) - n k=1 ∠((ω -λ k )(ω -λ k )) (11)
where ∠(z) denotes the angle of the complex number z lying in ] -π, π]. In others words, φ tot (ω) is the "unwrapped" version of ∠(q(ω)). The next few lines demonstrate that the total phase defined in ( 11) is one to one and onto from

R + onto ] -nπ, 0].
First, as M is assumed to be definite and positive, det M > 0 and ∠ det M = 0.

Then, the range of each partial factor

φ tot,k (ω) := -∠((ω -λ k )(ω -λ k )) (12)
can be easily established by rewriting it in the following way

φ tot,k (ω) = -∠(|λ k | 2 -ω 2 + (2α k )ω) ( 13
)
This function is well known in the frame of structural dynamics; it is strictly decreasing from [0, +∞[ onto ] -π, 0]. It follows that the sum φ tot of these n partial factors φ tot,k is a strictly decreasing function from [0, +∞[ onto ]-nπ, 0].

This implies that to each total phase value φ (j)

tot ∈] -nπ, 0] one and only one ω value can be associated and subsequently one and only one û value by virtue of Eq. ( 4). Hence, the φ tot variable can be used to parametrize the dynamical problem instead of a classical circular frequency parametrization: varying it continuously from 0 to -nπ let us compute the frf for the first n modes.

Instead of studying the frf by solving the following system

K(ω)û = f ω = ω (j) (14)
for given ω (j) values in [0, +∞[, one proposes here to consider

K(ω)û = f φ tot (ω) = φ (j) tot (15)
for given φ

(j)
tot values in ] -nπ, 0]. This new formulation implies to solve a real nonlinear constraint equation instead of a trivial equation to define ω value. Moreover, the current definition of the total phase requires the evaluation of all the eigenvalues λ k which can be a costly and inaccurate procedure for large systems: Sec. 3 proposes a preconditioning method and an approximation of φ tot (ω) for low frequencies which requires only the first eigenvalues evaluation. Finally, it is important to notice that, when uncertainties are introduced, the classical formulation (14) let the circular frequency be a non random variable while the constraint equation in system (15) depends on the system random matrices via eigenvalues λk and hence, ω becomes a random variable:

φ tot (ω) = -∠(det M ) - n k=1 ∠((ω -λk )(ω -λk )) (16)
The systems considered in the random case are then:

K(ω) ũ = f ω = ω (j) (17)
for given ω (j) values in [0, +∞[, and [START_REF] Ma | The decoupling of damped linear systems in oscillatory free vibration[END_REF] for given φ

K(ω) ũ = f φtot (ω) = φ (j) tot
(j) tot values in ] -nπ, 0].

Combination with Polynomial Chaos

First, the main points of a Polynomial Chaos expansion are briefly reminded. Then, the combination with the problem considered -Eq. [START_REF] Veselić | Damped Oscillations of Linear Systems: A Mathematical Introduction[END_REF] or Eq. ( 18) -is exposed leading to two different global algorithms.

Polynomial Chaos: main points and notations

Only the guidelines of the Polynomial Chaos expansion (PCE) are recalled here. A dimension-one stochastic space is considered: only one random variable ξ is used to introduce randomness in the system. The reader is referred to [START_REF] Xiu | Numerical Methods for Stochastic Computations -A Spectral Method Approach[END_REF] and references therein for a complete presentation of PCE.

Let (Θ, A, P) denote a probability space with Θ the event space, A the σalgebra on Θ and P a probability measure. The probability density function (pdf) associated to a random variable U is denoted p U . U expected value is defined as:

µ U = E[U ] = Θ U (θ) dP(θ) = R u p U (u) du (19)
and its standard deviation is denoted σ U :

σ U = E[(U -µ U ) 2 ] ( 20
)
With these notations, the random variable ξ is a function from Θ onto R:

ξ : Θ -→ R θ -→ ξ(θ). (21)
Considering a second-order random process X, the Polynomial Chaos expansion proposes to express it as a polynomial series X p using a set of N p orthogonal polynomials in the variable ξ and denoted ψ k :

X(θ) = X p (ξ(θ)) := Np-1 n=0 x n ψ k (ξ(θ)) (22)
where the order N p is theoretically infinite for general situations.

The deterministic coefficients x n can be evaluated using either an intrusive method or a non-intrusive one. The intrusive method follows a Galerkin approach: Eq. ( 22) is introduced in the equations governing X and theses equations are projected onto the set of orthogonal polynomials ψ k . The non-intrusive method uses the orthogonality of the polynomials with respect to a inner product denoted < •, • >:

x n = < X, ψ k > < ψ k , ψ k > ( 23
)
where the numerator is usually evaluated using a quadrature rule.

The intrusive method provides a set of n × N p coupled algebraic equations and often requires a special implementation while the non-intrusive approach determines the set of coefficients x n one after the other in an independent manner.

Already existing codes can then be used to evaluate X realizations needed for the quadrature.

The choice of the polynomial basis is somehow arbitrary; however some bases can be considered as optimal to describe common distributions [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF]. In the present case and throughout the paper, it will be assumed that the random variable ξ follows a uniform distribution over [-1, 1]:

p ξ (x) = 1 2 if -1 ≤ x ≤ 1 0 otherwise ( 24
)
This makes the Legendre polynomial basis the most natural choice. The first 6 polynomials (with degree less or equal to 5) are:

ψ 0 (x) = 1 ψ 3 (x) = 1 2 (5x 3 -3x) ψ 1 (x) = x ψ 4 (x) = 1 8 (35x 4 -30x 2 + 3) ψ 2 (x) = 1 2 (3x 2 -1) ψ 5 (x) = 1 8 (63x 5 -70x 3 + 15x) (25)
This set of polynomials is orthogonal with respect to the following inner product

< f, g >= 1 2 1 -1 f (x)g(x)dx (26)
and hence,

< X, ψ k >= E[X ψ k ] (27)
In the Illustrations section (Sec. 4), PCE coefficients x n will be evaluated following a non-intrusive approach by using a Gauss-Legendre quadrature rule to compute Eq. ( 23) numerator.

Once PCE coefficients x n are evaluated, there are two ways to post-process the result. First, the mean and variance can be directly computed (provided

ψ 0 (x) = 1): E[X p ] = x 0 < ψ 0 , ψ 0 > and E (X p -E[X p ]) 2 = Np-1 n=1 x 2 n < ψ k , ψ k > (28)
Second, the cumulative distribution function (cdf) and the pdf can be evaluated based on MC simulations. The difference with the usual processing is that X realizations are computed using its PCE X p (i.e. Eq. ( 22)) rather than solving the direct problem. This saves a lot of computational time and resource when the samples are large.

Global algorithms

In this paper, a non-intrusive method is used to evaluate PCE coefficients:

a Gauss-Legendre quadrature with N GL nodes is used to evaluate the inner product Eq. ( 26) involved in Eq. ( 23). Depending on the method appliedconstant frequency or constant total phase -the quantities whose PCE has to be evaluated differ.

In the classical case (Eq. ( 17)), for N frequency values ω (j) of interest, the PCE of ũ has to be evaluated: the linear system [START_REF] Veselić | Damped Oscillations of Linear Systems: A Mathematical Introduction[END_REF] has to be solved N GL times for a dynamic stiffness matrix evaluated at the quadrature nodes ξ

(i) GL , 1 ≤ i ≤ N GL . Algorithm 1 depicts the global procedure. û(j,i) ω denotes the realization of ũ for ω (j) and ξ (i)
GL while ũ(j) ω denotes the random vector solution of Eq. ( 17) for ω (j) .

Algorithm 1: Global procedure considering a constant frequency Data: System random matrices, Gauss-Legendre quadrature nodes ξ

(i) GL and ω (j) discrete values of interest. 1 foreach ω (j) , 1 ≤ j ≤ N do 2 foreach ξ (i) GL , 1 ≤ i ≤ NGL do 3 Compute û(j,i)
ω solution of Eq. ( 17). The procedure using a constant total phase is depicted in Algorithm 2. In this case, given values are those of total phase: φ (j) tot , 1 ≤ j ≤ N . For each of these values, one has to compute the PCE of both ũ and ω, solutions of Eq. [START_REF] Ma | The decoupling of damped linear systems in oscillatory free vibration[END_REF]. As the computation of φtot implies the computation of { M , C, K} eigenvalues, it is obviously more interesting to switch foreach instructions of the previous algorithm. û(j,i)

φ (resp. ω (j,i) φ
) denotes the realization of ũ (resp.

ωφ ) for φ (j) tot and ξ

(i)
GL while ũ(j) φ (resp. ω(j) φ ) denotes the random vector (resp.

random variable) solution of Eq. ( 18) for φ (j) tot .

Algorithm 2: Global procedure considering a constant total phase Data: System random matrices, Gauss-Legendre quadrature nodes ξ (k) and φ

(j) tot discrete values of interest. 1 foreach ξ (i) GL , 1 ≤ i ≤ NGL do 2 Compute eigenvalues λ k solutions of det(λ 2 k M + λ k C + K) = 0. 3 foreach φ (j) tot , 1 ≤ j ≤ N do 4 Compute ω (j,i) φ and û(j,i)
φ solution of Eq. ( 18).

5 end 6 end 7 foreach φ (j) tot , 1 ≤ j ≤ N do 8 Compute ũ(j)
φ and ω(j) φ PCE using quadrature rule (23). 9 end

Going back to a constant frequency description from a constant total phase study

A last point which will be illustrated in Sec. 4 is the possibility to rebuilt a constant frequency view of the stochastic response when doing the primary calculations using a constant total phase approach described in Eq. [START_REF] Ma | The decoupling of damped linear systems in oscillatory free vibration[END_REF]. Basically, its relies on a linear interpolation of the results obtained for discrete total phase values φ (j) tot in order to evaluate ũ(ω t , ξ (i) ) for each ξ (i) of the sample where ω t denotes the target frequency. First, one looks for the total phase discrete values which embrace ω t for the given ξ (i) value:

find j 0 ∈ {1, . . . , N } such that ω (j0,i) φ ≤ ω t ≤ ω (j0+1,i) φ (29) then, θ ∈ [0; 1] is defined such that ω t = (1 -θ)ω (j0,i) φ + θω (j0+1,i) φ (30)
and finally ũ is evaluated using the following interpolation:

ũ(ω t , ξ (i) ) = (1 -θ)û (j0,i) φ + θ û(j0+1,i) φ ( 31
)
Algorithm 3 states the different steps to get

ũ(i) t = ũ(ω t , ξ (i) ) values.
Algorithm 3: Evaluating constant frequency stochastic response from a constant total phase study

Data: ũ(j) φ and ω(j) φ PCE for φ (j) tot discrete values, 1 ≤ j ≤ N ; ξ (i) , sample of random variable values, 1 ≤ i ≤ Ns; ωt frequency of interest. 1 foreach ξ (i) , 1 ≤ i ≤ Ns do 2 Compute ω(j,i) φ = ω(j) φ (ξ (i) ), 1 ≤ j ≤ N using ω(j) φ PCE. 3 Find j0 such that ω(j 0 ,i) φ ≤ ωt ≤ ω(j 0 +1,i) φ . 4 Compute θ = (ωt - ω(j 0 ,i) φ
)/(ω

(j 0 +1,i) φ - ω(j 0 ,i) φ ). 5 Compute ũ(j 0 ,i) φ = ũ(j 0 ) φ (ξ (i) ) using ũ(j) φ PCE. 6 Compute ũ(j 0 +1,i) φ = ũ(j 0 +1) φ (ξ (i) ) using ũ(j+1) φ PCE. 7 Compute ũ(i) t = ũ(ωt, ξ (i) ) using Eq. (31). 8 end 9 Process ũ(i) t , 1 ≤ i ≤ Ns sample.

Computational considerations

First, as complex numbers angles are defined modulo 2π, it is obvious that the total phase should not be computed using the angle of the product of the n factors (ω -λ k )(ω -λ k ) but as the sum of the angle of each factor, that is the sum of each φ tot,k defined in Eq. ( 12). This said, this section addresses first the choice of a linearisation used to process the QEP and then the expression of an approximation which requires only the smallest eigenvalues. The last subsection is devoted to the efficiency and accuracy of this approximation.

Choice of a linearisation

Most methods used to compute the eigenvalues solution of Eq. ( 10) rely on a linearisation of the QEP [START_REF] Jacquelin | Polynomial chaos expansion in structural dynamics: accelerating the convergence of the first two statistical moment sequences[END_REF] [START_REF] Tisseur | The Quadratic Eigenvalue Problem[END_REF]Sec. 5]. Different linearisations are available in the literature [START_REF] Duncan | Matrices applied to the motions of damped systems[END_REF], [START_REF] Tisseur | The Quadratic Eigenvalue Problem[END_REF]Sec. 3.4], [START_REF] Saad | Numerical methods for large eigenvalue problems[END_REF]Sec. 9.3]. The usual formulation is as follows:

U 0 0 M B ẋ + 0 -U K C A x = 0 f ( 32
)
where U has to be invertible so as to ensure that

x k+n = dx k /dt, 1 ≤ k ≤ n.
The choice of U depends on the properties of system matrices [START_REF] Afolabi | Linearization of the quadratic eigenvalue problem[END_REF]. The one linearisation used here is the energy phase space linearisation proposed by Veselić [17, Sec. 3.2]: let L K and L M be such that in what follows the Cholesky decomposition is used to compute L • matrices. Then, use the following change of variables

K = L K L T K and M = L M L T M ( 33
u = L -T K x 1 , u = L -T M x 2 (34)
to rewrite the dynamical equation ( 1):

M L -T M ẋ2 + CL -T M x 2 + KL -T K x 1 = f (t) (35a) L -T M x 2 = L -T K ẋ1 (35b)
After pre-multiplying Eq. (35a) by L -1 M and a few transformations, one gets an equivalent linear system with size 2n:

ẋ = 0 L T K L -T M -L -1 M L K -L -1 M CL -T M Ae x + 0 L -1 M f , x = x 1 x 2 (36)
Linear problems (1), ( 32) and (36) have the same eigenvalues [17, Sec. 14] but not the same numerical properties as emphasized in Fig. 1. A simple clampedfree bending beam finite element model with a varying number of dofs is generated; modal damping is introduced once mass and stiffness matrices are assembled. The characteristics of the beam are the same than the clamped beam example studied in Sec. 4.2 except for limit conditions which are a perfect clamping on one end and with the other end being perfectly free; hence this model is fully deterministic. For each system size, linearisation (32) with U = -K and (36) are compared using first, the reciprocal condition number of B -1 A (resp. A e ) evaluated via Matlab function rcond, then the CPU times for computing all the eigenvalues using Matlab command eig(-B\A) (resp. eig(Ae))

or to compute the first 12 smallest eigenvalues -that is the first 6 λ k and their conjugates λ k -using eigs(-B\A,'sm',12) (resp. eigs(Ae,'sm',12)). CPU times include A and B matrices creation (resp. A e as well as L M and L K evaluation); displayed CPU times are the means over 5 identical computations for each case.

First, the reciprocal condition number is higher when using the energy phase space linearisation. Second, energy phase space linearisation leads to smaller CPU times when one evaluates the eigenvalues, especially when only the first smallest eigenvalues are to be extracted which is interesting considering the total phase approximation detailed below.

Similar results are obtained when using U = M .

Total phase approximation

In order to avoid the computation of all the eigenvalues λ k , one can first, use a truncation to eigenvalues matching the range of interest, that is, low frequencies: only N λ eigenvalues are kept, using the following criterion:

|ω k | ≤ 2π f max r (37)
where f max denotes the maximum frequency of interest and r a security coefficient greater than 1. Doing so, one also decreases the number of factors φ tot,k (ω) involved in the summation [START_REF] Loève | Probability Theory II[END_REF] which leads to a faster evaluation and an easier resolution of the constraint equation. The second step is, as when dealing with model reduction, to use a "static compensation" to take higher order modes influence into account. Considering small circular frequency values ω compared to eigen circular frequencies ω k for k > N λ , one can write, for each partial factor

φ tot,k (ω) = - 2α k |λ k | 2 ω + o(ω 2 ) (38)
Hence, an approximation of φ tot using a linear correction can be written:

φ tot (ω) ≈ - N λ k=1 ∠((ω -λ k )(ω -λ k )) + n k=N λ +1 - 2α k |λ k | 2 ω (39)
This approximation still requires the evaluation of higher order eigenvalues λ k .

To avoid their computation one uses one of the following equalities which are demonstrated in Appendix B, depending on the linearisation used (Eq. (32) or Eq. ( 36)):

n k=1 - 2α k |λ k | 2 = -tr(K -1 C) = -tr(L -1 K CL -T K ) (40)
This finally leads to an approximation of the total phase which requires the evaluation of K -1 C or L -1 K CL -T K and the first N λ eigenvalues λ k only; this final expression is denoted φ N λ :

φ N λ (ω) = - N λ k=1 ∠((ω -λ k )(ω -λ k ))+ -tr(L -1 K CL -T K ) - N λ k=1 - 2α k |λ k | 2 ω (41)
The previous expression is written in the case of an energy phase space linearisation: this is the only one used in the rest of the paper.

Efficiency and accuracy results

The numerical efficiency and accuracy of the several ways to evaluate the total phase are tested using the same clamped-free beam model as in Sec. 3.1. Uncertainty is added in a very simple way:

K(ξ) = K(1 + 0.05ξ) (42)
First, the complexity of several algorithms is examined. On the one hand, all the eigenvalues are evaluated using either linearisation (32) with U = -K or the energy phase space linearisation (36). Matlab functions used are eig(A,-B) and eig(Ae) respectively. On the other hand, following the definition of the approximated total phase (41) with N λ = 6, only the first 12 smallest eigenvalues λ k and λ k , 1 ≤ k ≤ N λ = 6 are evaluated using eigs(Ae,12,'sm') together with the matrix L -1 K CL -T K whose trace is required. These evaluations are realized for n ∈ {10, 20, . . . , 100, 200, . . . , 1000} and five ξ values each time, ξ ∈ {-1, -0.5, 0, 0.5, 1}. Results are depicted in Fig. 2. Dashed lines match the fitted linear interpolation in log-log scales for each case. First, as expected from results of Sec. 3.1, using an energy phase space linearisation is always cheaper than using a classical linearisation. Second, direct evaluation of all eigenvalues is cheaper for small systems (here, with less than n = 100 dofs), while for large systems, computing only some of them becomes a lot cheaper as the difference between the two slopes (3.44 in the first case and 2.01 in the second one) points out even though

L -1
K CL -T K has to be evaluated. The accuracy of the approximation is now tested: considering the clampedfree beam system for n = 300 and the nominal system only (ξ = 0), the total phase is evaluated at 10 000 points equally spaced over [0, 1.1 × ω 6 ] using either initial definition [START_REF] Loève | Probability Theory II[END_REF] or approximation (41) with N λ = 6. Both curves are displayed in Fig. 3 along with the relative error generated by the approximation (right-side ordinate axis, dashed dotted line). It appears that the relative error remains very small (less than 0.05%) over the whole frequency range. This example hence seems to indicate that using approximated total phase expression for low frequencies is very reasonable considering the accuracy criterion.

A final test considers the time used for solving the constraint equation in system [START_REF] Caughey | Classical Normal Modes in Damped Linear Dynamic Systems[END_REF]: when the total phase is evaluated using the approximated formulation, less factors are involved and one can expect that it is easier to solve. Figure 2: Total phase: complexity with respect to system size n for several computational methods

! 1 ! 2 ! 3 ! 4 ! 5 ! 6 -7:
-6:

-5:

-4:

-3:

-2:

-: To test this hypothesis, the total phase range [-6π, 0] is divided in 6 000 points linearly spaced φ (j) tot . The constraint equation is then solved for each point in decreasing order using fzero Matlab function and the previous solution value as the initial guess for the current resolution (first guess to initialize the loop when j = 1 is ω = 0). CPU time is collected once this resolution loop is applied 10 times successively. When using the initial definition of total phase [START_REF] Loève | Probability Theory II[END_REF], the whole procedure costs 55.7 s; when using the approximated expression (41), it lasts only 37.8 s, that is 32% less.

0 0 0.5 |φ tot -φ N λ | |φ tot | × 100 φ N λ φtot Relative error % ω φtot, φ N λ
Comparison of the cost of global procedures including Polynomial Chaos expansions when considering either Eq. ( 17) or Eq. ( 18) and using initial total phase [START_REF] Loève | Probability Theory II[END_REF] or approximated total phase (41) expressions will be addressed through the numerical examples developed in Sec. 4.

Illustrations

The proposed method is applied to two systems. The first one is a very simple 2 mass-spring damped system which will illustrate the ability of the proposed method to return the expected accurate results. The second one is a beam with imperfect clamping conditions at both ends; it will illustrate the the accuracy of the truncated expression for the total phase and the relative cost of the methods.

For both systems, the usual approach, using a constant frequency and the proposed approach, using a constant total phase will be compared. For the sake of brevity the two approaches will be further referred to as CFA for Constant Frequency Approach and CTPA for Constant Total Phase Approach.

For both examples and both approaches, only one random variable is introduced and degree 5 Polynomial Chaos with Legendre polynomials expressed in Eq. ( 25) is used. Coefficients are evaluated using a Gauss-Legendre quadrature rule with 11 nodes for numerator and denominator expressed in Eq. ( 23). When considering CTPA, an energy phase space linearisation is used to solve the QEP.

Stochastic dimension is kept small on purpose so as not to mix the "curse of dimensionality" related to polynomial chaos when multiple variables are introduced in the model and the point of the paper which is a new parametrization of the dynamical system response in order to get more accurate PCE results.

2-dofs system

The first system used to illustrate and compare both approaches is a simple 2-dofs system depicted in Fig. 4 whose matrices are:

M = m 0 0 m , C = c 1 + c 2 -c 2 -c 2 c 2 , K = k 1 + k 2 -k 2 -k 2 k 2 (43)
Only the second mass undergoes an external load: f (t) = {0, 1} T cos(ωt). Numerical values for nominal system are m 1 = 1 kg, In the random case, the first spring stiffness depends on a random variable ξ:

c 1 = 1 kg.s -1 , c 2 = m k1 c 1 k 2 m c 2 u 1 u 2 f 2 (t) = cos(ωt)
k1 = k 1 (1 + ξ∆k 1 ) (44)
with ∆k 1 = 0.75; ξ is assumed to have a uniform distribution over [-1, 1].

A complex formulation is used to solve the linear dynamical problem: unknowns are the complex variables û1 and û2 .

A similar case with a Gaussian random variable instead of a uniform one is illustrated in Appendix D.

MCS results

First, Monte Carlo simulations are run in order to visualize the differences between the two approaches. The dynamical system is solved for 100 001 values of ξ equally spaced over [-1, 1] and for 201 values of ω = 2πf equally spaced over [0, 5] × 2π to illustrate the CFA and for 201 values of φ tot over [-1.9318 × π, 0] which are not equally spaced (see Appendix C) to illustrate the CTPA. These simulations will serve as references to evaluate the accuracy of each approach. Fig. 5 shows the complex response of the second dof (û 2 ) for different values of ξ: amplitude |û 2 | and angle ∠û 2 (that is the phase between u 2 and the excitation) are plotted versus the excitation frequency f and the total phase φ tot is superimposed to the angle in the lower plots. This shows first that contrarily to the angle, the total phase is continuously decreasing when the frequency increases and that important total phase slopes happen around resonances and with values centred around -π/2 and -3π/2. Second, two frequencies are shown on the plots (1.05 Hz and 3.4 Hz). Considering different values of ξ, it is obvious that for a given frequency, very different responses happen in terms of amplitude, angle, and total phase. On the contrary, considering the response for a given total phase value (-0.5π or -1.5π for instance) narrows the variability of the responses. This latter observation is even more obvious when looking at Fig. 6. Amplitude and angle values are rendered via a grey scale and visualized in either a (f, ξ) plane in column (a) or in (φ tot , ξ) plane in column (b). For both views, the paths of a constant frequency study and a constant total phase study are displayed: constant total phase paths do not cross very different grey tones that is, tend to describe similar states of the system in terms of amplitude and angle, whereas constant frequency paths embrace large ranges of grey tones which implies larger variations to be rendered via PCE. Similar results could be plotted for the first dof û1 .

PCE results

Only the amplitudes |û k | (not angles) are fitted using polynomial chaos. Both algorithms (Algorithm 1 for CFA and Algorithm 2 for CTPA) are applied using the same samples in ω or φ tot as for MCS. For CTPA, total phase formula used is Eq. ( 11) without truncation as the system is small and all its modes are within the frequency range of interest. Fig. 7, 8 and 9 compare results obtained via both methods with MCS results to evaluate their global accuracy. The upper row in Fig. 7 shows the values of |û 2 | evaluated from PCE for each of the 201 ω (j) values (column (a)) or φ (j) tot values (column (b)). The lower row shows the relative error between PCE and MCS. The error scales for both methods are very different: more than 100% error arise with the CFA while errors are less than 5% in the case of CTPA. This shows that CTPA is accurate even for low PCE orders. Fig. 8 provides more general indicators, namely error on mean and standard deviation for each method. It shows that CFA errors in mean and variance increase especially around resonances while the errors stay contained in the CTPA case. As ω becomes a random variable in the CTPA case, the same error plots are provided in Fig. 9 to show that the nonlinear constraint equation of system ( 18) is solved with a great accuracy too. resulting pdfs do not match the expected ones at all. On the contrary, in case of a CTPA, |û k | polynomial representation match the real dependencies very accurately leading to very accurate pdfs as depicted in Fig. 11.

Finally, Algorithm 3 is applied to evaluate |û k | distributions for eigenfrequencies from CTPA results. Fig. 12 shows than it leads to much more accurate results than when using CFA directly (Fig. 10).

The case when a Gaussian random variable is used illustrated in Appendix D demonstrates that the accuracy does not depend on the law that the random variable follows: the same qualitative results are found.

Clamped beam

The second application is a bending beam with imperfect clamping conditions at both ends; this can arise for example when assemblies with bolted joints are studied. The system is depicted in Fig. 13. It consists in a beam with length L = 1 m, width b = 1 cm and height h = 2 cm made of steel (Young modulus E = 2 • 10 11 Pa, density ρ = 7800 kg.m -3 ). A finite element model with 100 nodes and 2 dofs per node is computed leading to mass and stiffness matrices, M and K b respectively. One translational and one rotational spring (with respective stiffness k t and k r ) link each extremity to the ground; these stiffness are assumed to introduce uncertainty modelled as follows: kt = 10 7 (1+0.5ξ) N.m -1 and kr = 10 5 (1 + 0.8ξ) N.m.rad -1 where ξ follows a uniform law over [-1, 1]. This leads to random translational and rotational stiffness matrices Kt (ξ) and Kr (ξ) which are added to K b to get the system random stiffness matrix:

K(ξ) = K b + Kt (ξ) + Kr (ξ) (45)
To ensure sub-critical damping for all the modes, modal damping is applied. First, undamped modes are evaluated using mass matrix M and nominal system stiffness matrix K(ξ = 0). This provides 200 circular eigenfrequencies ω k and eigenvectors ϕ k . Then the diagonal damping matrix in the modal basis is created C = diag(2η k ω k ) with damping rates η k logarithmically spaced between η 1 = 0.02 and η 200 = 0.2. Finally, the damping matrix in the physical basis C is computed. An ascending excitation is applied on the 149-th dof (located approximately at 0.75 m from the left end): f (t) = 100 cos(ωt). The frequency range of interest [0, 600] Hz embraces the first 3 modes of this 200 dofs model. Fig. 14 provides an overview of the system response variation over the frequency range of interest and shows that each of the three modes varies. Two methods are used for total phase evaluation: full as defined by Eq. ( 11) and truncated with N λ = 5 and a "static compensation" as defined by Eq. (41).

This will demonstrate the accuracy of the truncated formula and provide elements for cost comparison of all the methods. As in the previous application only amplitudes |û k | (not angles) are approximated using polynomial chaos. Two dofs will be monitored: the displacement u 149 where the force is applied and the displacement u 101 of the node located approximately at 0.51 m of the left end of the beam.

Fig. 15 provides an overview of the system response variations by displaying dof 149 amplitude for various ξ values when evaluated using direct computation. Constant frequency and constant total phase lines are drawn for the first three resonances (thick lines) and the two anti-resonances (thin lines). This already shows that constant frequency lines cross responses with larger variations than constant total phase lines.

Accuracy

Fig. 16 displays relative error on mean and variance over the frequency range of interest in the CFA case (graphics (a) and (b)) and over the total phase range of interest in the CTPA case (graphics (c) and (d)). The grey patches match areas around resonances depicted in Fig. 17. These areas are defined as follows: for each ξ value in {-1, 0, +1}, the frf is computed over the whole frequency range and the maximum value over all translational displacements |û 2k+1 | is stored for each frequency. Then, each resonance is located and frequencies or total phase matching neighbouring points with amplitude greater than half this resonance amplitude are used to define the areas.

Another way to evaluate the accuracy of the results is to use the Wasserstein distance [START_REF] Villani | Optimal transport. Old and new[END_REF]Chap. 6]. This distance provides a means to compare the quality of the approximation considering simultaneously the convergence on the mean, the variance and the distribution law. The Wasserstein distance (Wd) of order 2 between two probability measures P 1 and P 2 having a second order moment is defined as [START_REF] Villani | Optimal transport. Old and new[END_REF]: where L(X) denotes the law of X. It can be evaluated by:

Wd 2 (P 1 , P 2 ) = inf E (X 1 -X 2 ) 2 1/2 : L(X 1 ) = P 1 , L(X 2 ) = P 2 , (46)
Wd 2 (P 1 , P 2 ) = 1 0 (F -1 X1 (t) -F -1 X2 (t)) 2 dt 1/2 , (47)
where F X denotes the cumulative distribution function associated to X. The distance between the CFA or CTPA approximations and the reference sample U ref can then be evaluated using Eq. (47). To render this distance dimensionless with respect to the order of magnitude of U [START_REF] Del Barrio | Tests of goodness of fit based on the L 2 -Wasserstein distance[END_REF], the Wasserstein distance of order 2 is divided by the standard deviation evaluated using the reference sample; this quantity is denoted Wd:

Wd(P U , P U ref ) = 1 0 (F -1 U (t) -F -1 U ref (t)) 2 dt σ 2 U ref 1/2 . ( 48
)
These distances are plotted on Fig. 18 which shows that truncated and full CTPA methods provide the same accuracy and a greater accuracy than CFA.

To get other illustrations of the Wasserstein distance in a mechanical frame, the reader is referred to [START_REF] Chouvion | Development of error criteria for adaptive multielement polynomial chaos approaches[END_REF][START_REF] Pagnacco | Pitfalls in the frequency response represented onto polynomial chaos for random SDOF mechanical systems[END_REF].

As in the previous application, the CFA obtained by CFA while the right column displays results obtained by CTPA with N λ = 5 and using the constant frequency rebuild process from Sec. 2.3. This shows that the CTPA method which is designed to provide a greater accuracy around resonances is able to provide accurate results around anti-resonances too.

Numerical cost

No proper study with multiple runs was conducted. CPU times provided are hence given for guidance only. 601 ω (j) values over [0, 600 × 2π] and 601 φ (j) tot values between [-3.0298π, 0] were used. As a Gauss-Legendre quadrature rule with 11 nodes is used to evaluated PCE coefficients, each problem [START_REF] Veselić | Damped Oscillations of Linear Systems: A Mathematical Introduction[END_REF] or [START_REF] Ma | The decoupling of damped linear systems in oscillatory free vibration[END_REF] has to be solved for 11 different ξ (i) GL values. In the case of CTPA the constraint equation requires a numerical resolution and, for each ξ (i) GL value, the 2 × N λ = 10 smallest eigenvalues (if truncated) or all of them (if not truncated) have to be evaluated. MCS were run using 10001 ξ values in the sample. Time for MCS is about 4 h for CFA, and 5 h 30 min for CTPA (+37.5%).

Most of the additional time is due to the total phase function creation which implies computing L K and L M matrices for energy phase space linearisation and evaluating eigenvalues. Resolution of the nonlinear constraint equation using Matlab fzero function only accounts for about 3% of this additional time. Time for PCE is about 15.5 s for CFA, 22.3 s for CTPA without truncation (+43.9%) and 21.5 s for CTPA truncated with N λ = 5 (+38.7%). The gain of evaluating only some of the eigenvalues is not important as expected considering Fig. 2; however, it would increase for larger systems. If CTPA implies higher costs than CFA, it must be kept in mind that it returns accurate results and is more than 600 times cheaper than a MCS. The time to rebuild the three constant frequency plots in Fig. 20 from CTPA result is about 0.9 s.

Conclusions

A new parametrization of linear dynamical problems is proposed for systems with sub-critical damping. It uses a total phase parameter which involves the eigenvalues of the damped problem. It is demonstrated that one and only one dynamical solution (couple frequency/displacements) is associated to any value of this parameter. Stochastic dynamic responses are hence computed using given values for this total phase parameter rather than using given values for the frequency.

Two examples show that Polynomial Chaos Expansion with low degree then provides accurate results even in the vicinity of resonances that is where it usually looses accuracy when considering given frequency values. These illustrations show that the proposed parametrization gathers similar responses contrarily to 29 classical constant frequency parametrization which embraces very different dynamical responses. Even though the current examples involve only one random variable, they prove that low degree PCE provide accurate results because the response variations are simpler to describe when considering a given total phase value rather than a constant frequency value. Hence, in the case when multiple random variables are involved, the "curse of dimensionality" the PCE method suffers from will be limited as fewer monomial will have to be considered.

The limits of the method are not explored here. Further work should test it on systems for which modes switch for different realizations which can happen when considering rotating machineries for example. One can guess that in this case, the accuracy will not be as good as for the systems tested in this paper and that some improvements should be developed. The case when some modes are over-damped, should be examined too. If high frequency over-damped modes should not be a problem for low frequency studies, the case when low frequency modes do not respect the subcritical damping rule could be more difficult to handle. resonances and a spacing such that one gets an approximate linear spacing in frequencies in between these ranges. Let's assume that one wants to study the response of a structure over its first K modes included in [0, ω max ] range for circular frequency. The equivalent range in total phase is [φ tot max , 0] where φ tot max is evaluated from Eq. ( 11) or Eq. ( 41) and ω max . First, one has to define ranges around resonances: resonance for mode k occurs when the k-th partial total phase factor φ tot,k is close to -π/2. Let's define the limits around -π/2 for which most of the resonance peak is embraced: φ - lim = -π/10 and φ + lim = -9π/10. ω - k and ω + k such that φ tot,k (ω ± k ) = φ ± lim can be computed from these values, using Eq. ( 12). This is illustrated in Fig. C.23(a). Corresponding total phase limit values φ ± tot k can then be easily computed using either Eq. [START_REF] Loève | Probability Theory II[END_REF] or Eq. (41). This is illustrated via with ∆k 2 = 0.15 and ξ ∼ N (0, 1). ξ realizations leading to negative stiffness are put aside. The computations are the same than in Sec. 4.1 except that a Hermite polynomial basis is used for PCE (still with degree 5).

This configuration leads to an almost constant first mode but a second mode with large variations as depicted in Fig. D.26. This figure shows again that while constant frequency cuts account for very different responses for the second mode, constant total phase cuts gather very similar responses. Hence, the response approximation by PCE is more accurate when considering CTPA than CFA. This is illustrated in Fig. D.27 which shows the results obtained considering a constant frequency cut around the second mode. While the direct degree 5 polynomial approximation of û2 in the case of CFA presents high errors leading to erroneous pdfs (column (a)), the degree 5 approximations in the CTPA case have good enough accuracies which provide a mean to rebuild a very accurate constant frequency response (column (b)) using the steps explained in Sec. 2.3.

4 end 5

 45 Computeũ(j) ω PCE using quadrature rule[START_REF] Villani | Optimal transport. Old and new[END_REF]. 6 end

Figure 1 :

 1 Figure 1: Reciprocal condition number (a) and CPU time to solve the full eigenvalue problem using eig (b) or to compute only the first 2×6 smallest eigenvalues using eigs(•,'sm',12) (c) depending on the system size n and the linearisation. Solid (blue) line: classical linearisation; Dashed (red) line: energy phase space linearisation; Dashed dotted (black) line: ratio between energy phase space and classical linearisation.

  eigs(Ae,12,'sm'), L -1 K CL -T K slope: 2.50 eig(Ae) slope: 3.44 eig(A,-B) n CPU time (s)

Figure 3 :

 3 Figure 3: Total phase: accuracy of approximation φ N λ (N λ = 6)

Figure 4 : 2 -

 42 Figure 4: 2-dofs system with random stiffness k1

Figure 5

 5 Figure 5: 2-dofs, MCS -System second dof response (û 2 ) for different ξ values. Vertical grid lines denote the two frequencies used for further constant frequency study while the horizontal grid lines in the phases plots denotes the two total phase values used for further constant total phase study.

Fig. 10 compare

 10 MCS and PCE results for CFA at the eigenfrequencies of the nominal system: |û k | dependencies to ξ are very poorly represented and the

Figure 6 : 2 -Figure 7

 627 Figure 6: 2-dofs, MCS -System second dof response (û 2 ) versus (a) f and ξ, (b) φtot and ξ. Dashed (red) lines: constant frequency paths; Solid (green) lines: constant total phase paths.

Figure 8 :

 8 Figure 8: 2-dofs, PCE, Accuracy comparison -Comparison of MCS and PCE results using global indicators over the full frequency range when using (a,b) CFA and (c,d) CTPA. (a,c) mean and (b,d) standard deviation of |û k |; thick lines depict MCS results while thin lines denote the absolute difference between MCS and PCE results.

Figure 9 : 2 -Figure 10 : 2 -

 92102 Figure 9: 2-dofs, PCE, Accuracy comparison -Comparison of MCS and PCE results using global indicators over the full frequency range when using CTPA. (a) mean and (b) standard deviation of ω; thick lines depict MCS results (values are to be read on left scale) while thin lines denote the absolute difference between MCS and PCE results (values are to be read on right scales).

Figure 11 : 2 -

 112 Figure 11: 2-dofs, CTPA -Comparison of MCS and PCE results for constant total phases. Upper plots: |û k | versus ξ; lower plots: pdf of |û k |.

Figure 12 : 2 -

 122 Figure 12: 2-dofs, Constant frequency results rebuilt from CTPA results -Comparison of MCS and PCE results for constant frequencies. Upper plots: |û k | versus ξ; lower plots: pdf of |û k |.

Figure 13 :Figure 14 :

 1314 Figure 13: Clamped beam with random stiffness kt and kr.

Figure 15 :

 15 Figure 15: Beam, MCS -System 149-th dof response amplitude (û 149 ) versus (a) f and ξ, (b) φtot and ξ. Dashed (red) lines: constant frequency paths; Solid (green) lines: constant total phase paths.

  Relative error on µ| ûk | (%)

Figure 16 :

 16 Figure 16: Beam, PCE, Accuracy comparison -Comparison of MCS and PCE results using global indicators over the full frequency range when using (a,b) CFA, (c,d) CTPA truncated with N λ = 5. Relative error in percent for (a,c) mean and (b,d) standard deviation of |û k |. Solid line is relative to |û 101 |; dashed dotted line is relative to |û 149 |.

Figure 17 :Figure 18 :

 1718 Figure 17: Beam, PCE, Accuracy comparison -Patches: Areas with amplitudes greater than half the resonance amplitude for (a) frequency and (b) total phase. Thin lines: realizations for ξ ∈ {-1, 0, +1}. Thick lines: portion of curves with amplitudes greater than half the resonance amplitude.

Figure 19 :

 19 Figure 19: Beam, CFA -Comparison of MCS and PCE results for constant frequencies. Upper plots: |û k | versus ξ; lower plots: pdf of |û k |. |û 101 | is not plotted in column (b) as its amplitude range is very low compared to |û 149 |.

Figure 20 :

 20 Figure 20: Beam, Constant frequency results rebuilt from CTPA truncated with N λ = 5 -Comparison of MCS and PCE results for constant frequencies. Upper plots: |û k | versus ξ; lower plots: pdf of |û k |. |û 101 | is not plotted in column (b) as its amplitude range is very low compared to |û 149 |.

Figure 21 :

 21 Figure 21: Beam -Comparison of MCS and PCE results around the second anti-resonance (f = 430 Hz). Upper plots: |û k | versus ξ; lower plots: pdf of |û k |. (a) Results from CFA; (b) Results rebuilt from CTPA truncated with N λ = 5.

Figure C. 22 : 2 -

 222 Figure C.22: 2-dofs dynamical system response evaluated at 17 points linearly spaced (a) in frequency and (b) in total phase.

Figure C. 23 :ξ

 23 Figure C.23: Sampling zones definition

  Figure D.27: 2-dofs Gaussian -Comparison of MCS and PCE results around the second mode (f = 3.4 Hz). Upper plots: |û k | versus ξ; lower plots: pdf of |û k |. (a) Results from CFA; (b) Results rebuilt from CTPA.

Appendix A. Evaluating the response using a complex rewriting of the dynamical equation Injecting excitation and amplitude decompositions [START_REF] Wang | Frequency response function-based model updating using Kriging model[END_REF] and (3) into the dynamical equation [START_REF] Pichler | A mode-based meta-model for the frequency response functions of uncertain structural systems[END_REF] and balancing cosine and sine terms, one gets 2 real equations with size n:

To compute a and b, one can either solve this linear system with size 2n or consider the following composition in the complex space C n , (A.1a) -(A.1b):

By analogy with the static case (ω = 0), K(ω) is often called the dynamic stiffness matrix and its inverse, denoted Ŝ(ω) the dynamical flexibility matrix.

Appendix B. Demonstration of equation (40)

Case of a classical linearisation (32). The first step is to express the inverse of matrix A defined id Eq. (32). Since K is supposed to be symmetric definite positive, it is invertible and one can write

whatever U matrix is, provided it is invertible. The next step is to note that since λ k are the eigenvalues of the linear eigenvalue problem associated to Eq. (32):

Finally the detailed expression of A -1 B,

Case of an energy phase space linearisation (36). In this second case, only one matrix A e is involved. Its inverse is:

The result is straighforward:

Appendix C. Phase sampling

The method to create samples of total phase values for which system response should be evaluated is not obvious. Indeed, using equally spaced values would not return a relevant mesh of the dynamic response as illustrated via The good point of an equally spaced total phase sample is that it focuses on the resonance peaks; the bad point is that there is not enough points in between these ranges. On the contrary, an equally spaced frequency sample does not focuses on resonance peaks and can somehow miss them but describes ranges outside resonances more properly. A mix is then proposed to build an efficient total phase sampling: a linear spacing of total phase values is used around

• Z a : zone after the last resonance. ω ∈ [ω + K , ω max ] and

In zones Z k m , total phase points are equally spaced in total phase; in zones Z b , Z k i and Z a , total phase sampling is such that points are equally spaced in frequency (this is done by computing images of linearly spaced samples in ω using either Eq. [START_REF] Loève | Probability Theory II[END_REF] or Eq. ( 41)).

The last thing to determine is the ratio of points in each zone. These ratios are up to the user and may vary from one study to another. The ratios used throughout the paper are 25% for Z k m zones, and the other 75% are distributed such that Z b and Z a are described with the same number of points and each Z k i zone has twice the number of points of Z b zone. The resulting percentage of points in each zone for K = 2 and K = 3 modes are provided in Tab. C.1. To illustrate the interest of such a phase sampling, dynamic responses of the two systems used for applications in Sec. 4 are computed using either linearly spaced values in frequency (as is usually done) and total phase samples computed as explained above, using the same number of points in both samples. The error in the rebuilt responses is evaluated for one dof by the integral of the absolute difference between the resulting curve and a reference curve obtained via a very fine sampling in frequency (area between both curves). Results are displayed on As expected, considering a same number of points N , the resonance peaks are described more accurately using the proposed total phase sampling method while keeping enough points in between these areas. Moreover, the error with a reference curve is very much smaller in the general case which indicates that it may be interesting to use such a sampling to compute frf with a minimal number of points even in the deterministic case.

Zone

Appendix D. 2-dof system with a Gaussian random variable

The 2-dof system illustrated in Sec. 4.1 is considered but with a deterministic k 1 value and a random k 2 variable which follows a normal law: