M Gerin

N Marcelino

N Biver

E Roueff

L H Coudert
email: coudert@lisa.univ-paris12.fr

M Elkeurti

D C Lis

D Bockelée-Morvan

Detection of 15 NH 2 D in dense cores: A new tool for measuring the 14 N/ 15 N ratio in the cold ISM. ⋆

Keywords: ISM clouds -molecules -individual object (Barnard-1b, L1689N, L134N(S), L1544, NGC1333-IRAS4A)radio lines: ISM

Context. Ammonia is one of the best tracers of cold dense cores. It is also a minor constituent of interstellar ices and, as such, one of the important nitrogen reservoirs in the protosolar nebula, together with the gas phase nitrogen, in the form of N2 and N. An important diagnostic of the various nitrogen sources and reservoirs of nitrogen in the Solar System is the 14 N/ 15 N isotopic ratio. While good data exist for the Solar System, corresponding measurements in the interstellar medium are scarce and of low quality. Aims. Following the successful detection of the singly, doubly, and triply deuterated isotopologues of ammonia, we have searched for 15 NH2D in dense cores, as a new tool for investigating the 14 N/ 15 N ratio in dense molecular gas. Methods. With the IRAM-30m telescope, we have obtained deep integrations of the ortho 15 NH2D (11,1 -10,1) line at 86.4 GHz, simultaneously with the corresponding ortho NH2D line at 85.9 GHz. Results. The ortho 15 NH2D (11,0 -10,1) is detected in Barnard-1b, NGC1333-DCO + , and L1689N, while we obtained upper limits towards LDN1544 and NGC1333-IRAS4A, and a tentative detection towards L134N(S). The para line at 109 GHz remains undetected at the rms noise level achieved. The 14 N/ 15 N abundance ratio in 15 NH2D ranges between 350 and 850, similar to the protosolar value of ∼ 424, and likely higher than the terrestrial ratio of ∼ 270.

Introduction

Nitrogen chemistry is particularly interesting for understanding the connection between the ISM and the formation of the solar nebula, because it is thought that the primitive atmospheres were nitrogen rich, as Titan remains today. Furthermore, the isotopic 15 N/ 14 N ratio has been measured in a variety of Solar System bodies, from the giant planets to the rocky planets, comets, and meteorites. The observed differences in nitrogen fractionation are used to understand how these bodies formed within the protosolar nebula. The combination of nitrogen and hydrogen (D/H) isotopic ratios has been demonstrated to be a very effective way of understanding how the ice mantles were enriched in deuterium and nitrogen. Aléon and Robert (2004) have concluded that a fast condensation Send offprint requests to: e-mail: maryvonne.gerin@ens.fr ⋆ Based on observations obtained with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). of the organic matter, enriched in 15 N and deuterium, is needed in order to keep a significant fractionation in the solid material of the primitive Solar System. They also have evaluated the exothermicity of the fractionation reactions for nitrogen to be 43 ± 10 K. The D fractionation has not been inherited from the native prestellar core, but most likely occurred in the protosolar nebula (author?) (Remusat et al. 2006, Gourier et al. 2008), yet the same physical and chemical processes are thought to operate in the prestellar cores and in the coldest regions of circumstellar disks. In the ISM, recent observations show that, contrary to CO, nitrogen does not deplete from the gas phase in dense cores, except when the density rises significantly above 10 6 cm -3 . Nitrogen species can therefore be very significantly deuterated, with D/H fractionation of several tenths for N 2 D + (Daniel et al. 2007;Pagani et al. 2007) and NH 2 D (Crapsi et al. 2007). Multiply deuterated ammonia in particular can be very abundant (Gerin et al. 2006;Lis et al. 2002a;Lis et al. 2006;Roueff et al. 2005). Nitrogen molecules will therefore be significant molecular reservoirs of deuterium. It is interesting to study whether they could also be enriched in 15 N, and whether signatures from an enrichment at an early evolutionary stage can be identified in primitive matter.

High 15 N enhancements are measured both in HCN and CN cometary gases (Bockelée-Morvan et al. 2008;Schulz et al. 2008), and in primitive carbonaceous meteorites. High 15 N enhancements may have been present in the ammonia ices of the natal presolar cloud according to the fractionation mechanism proposed by Rodgers & Charnley (2008a;2008b;2004) and Charnley & Rodgers (2002). Nitrogen fractionation is not expected to be as efficient as deuterium fractionation in dense cores, yet significant departures from the elemental 14 N/ 15 N ratio may occur in some molecules. As first shown by Terzievia & Herbst (2000), and developed by Charnley & Rodgers (2002) and Rodgers & Charnley (2008a), nitrogen fractionation in the gas phase may operate through ion-molecule reactions involving atomic or ionized nitrogen . Rodgers and Charnley (2008b) have subsequently studied the possible role of neutral-neutral reactions involving 15 N and CN. Little observational interstellar data are available. We have therefore started a survey of the main nitrogen-bearing interstellar NH2D 11,123456789101234561234567891012345species in 5 dense cores and a class 0 source (Table 1). This paper reports the detection of o-15 NH 2 D as the first result of this survey.

Observations

The microwave and far infrared spectra of 15 NH 2 D and 15 NHD 2 have been recently investigated by Elkeurti et al. (2008) and used to produce the corresponding line lists as supplementary data1 , while accurate line lists and partition functions for the 14 N isotopologues of the NH 3 family can be found in Coudert & Roueff (2006). These species are also independently listed in the Cologne Database for Molecular Spectroscopy (CDMS, Müller et al. 2001Müller et al. , 2005)), with small differences in the line frequencies due to different handling of the hamiltonians.

We have chosen to search for the 1 1,1 -1 0,1 line of ortho 15 NH 2 D since the corresponding NH 2 D line is very strong and both the sky transmission and telescope performances are excellent at 86 GHz. The frequency shift introduced by the 15 N substitution is small enough that the two isotopologues can be observed with the same receiver tuning. The line frequencies (Elkeurti et al. 2008), Einstein A coefficients, upper energy levels and critical densities are listed in Table 2. We have used the theoretical estimates of the critical densities from the reduced mass ratio scaling of Machin & Roueff (2006) for the NH 2 D-He values at 10 K, the temperature appropriate for the cold cores we have observed. However these values are probably too large when molecular hydrogen is involved, as found in recent calculations of the NH 3 -H 2 system by Valiron et al. (private communication).

The observations have been performed with the IRAM-30m telescope, during three observing sessions in December 2007, March 2008 and September 2008. We used the A100 and B100 receivers in parallel, tuned to 86.2 GHz in order to detect o-NH 2 D and o-15 NH 2 D with the same detector setting. The weather conditions were average, with 5 -10 mm of water vapor (PWV) . The NH 2 D and 15 NH 2 D lines were observed simultaneously, with the J=1-0 lines of H 15 NC and H 13 CN at 86.055 GHz and 86.338 GHz. We used the VESPA correlator, tuned to a spectral resolution of 40 kHz, and spectral bandpass of 40 MHz for each line. The data were taken using the wobbling secondary reflector, with a beam separation of 240". Telescope pointing was checked on nearby planets and bright radio quasars and was found accurate to ∼ 3 ′′ . Due to rather poor weather conditions during the September run (high PWV and cloudy sky), the pointing accuracy was degraded to ∼ 5 ′′ . Additional observations of the p-15 NH 2 D line at 109.3 GHz were obtained in March 2008. We only searched for this line towards Barnard-1b and detected no signal down to a rms noise level of 18 mK with 0.2 km s -1 velocity resolution. For L134N(S), we combined the data with observations performed in April 2005, as part of the dark cloud line survey project (Marcelino et al. 2009). The weather conditions were excellent (1 -2 mm PWV) and the observations performed in the frequency switching mode.

The data processing was done with the GILDAS 2 software (e.g. [START_REF] Pety | EdP-Sciences[END_REF]. We used the dec08b version of this software, which allows to correct for a minor bug in the frequency calibration during the observations. The IRAM-30m data are presented in main beam temperatures T mb , using the forward and main beam efficiencies F eff and B eff appropriate for 86 GHz, F eff =0.95 and B eff =0.78. The uncertainty in flux calibration is ∼ 10%, as checked by the variation of the intensity of the strong o-NH 2 D and H 13 CN lines in the spectrum. Linear baselines were subtracted.

Because the nuclear spin of 15 N is 1/2, the 15 NH 2 D lines are split into fewer hyperfine components than NH 2 D which makes their detection more favorable. The hyperfine structure of 15 NH 2 D is driven by the quadrupole moment of the deuterium nucleus, which is much smaller than the corresponding value of 14 N. We have checked, by using the nuclear quadrupole constants provided in Garvey et al. (author?) (1976), that the resulting hyperfine splitting is less than 50 kHz. We can thus safely assume that the spectrum reduces to a single component . As shown in Figure 1, the 15 NH 2 D line is clearly detected towards Barnard-1b, and L1689N, while we obtained upper limits towards LDN1544 and NGC1333-IRAS4A and tentative detections towards NGC1333-DCO + and L134N(S). The ratio of peak antenna temperatures of the NH 2 D and 15 NH 2 D lines is 50 -100, and the velocity agreement is excellent. Using the JPL and CDMS spectroscopy data bases, we have checked that no line of known interstellar molecules are expected within ±300 kHz from the 15 NH 2 D line frequency. The identification of the detected feature is therefore secure.

The line parameters were estimated by fitting Gaussian profiles to the detected o-15 NH 2 D lines. For o-NH 2 D we used the HFS routine implemented in CLASS, which allows to take into account the hyperfine components selfconsistently. The opacity of the ortho NH 2 D line is moderate in all sources, with a total opacity for all lines ranging from ∼ 1 to ∼ 5 (Table 3).

Results

NH 2 D and 15 NH 2 D

Results of the fits and derived molecular column densities are listed in Table 3. As we are mostly interested in the ratio of column densities, we computed them under the simple assumption of a single excitation temperature. We used the excitation temperature derived from the NH 2 D fit for both isotopic species. The o-NH 2 D column densities are in good agreement with previously published results for the same sources (Roueff et al. 2005). The [NH 2 D]/[ 15 NH 2 D] abundance ratio range from 360 to 810, with the largest value for 2 See http://www.iram.fr/IRAMFR/GILDAS L1689N. This last source is an interaction region between a molecular outflow and a dense core, and as such may have peculiar properties (Lis et al. 2002b). Given the error bars, the measured [NH 2 D]/[ 15 NH 2 D] ratio is comparable to the 14 N/ 15 N protosolar ratio, as measured in Jupiter (450; Fouchet et al. 2004) and in osbornite-bearing calcium-aluminium-rich inclusions from meteorites (424; Meibom et al. 2007), and likely larger than the terrestrial abundance ratio (270). Although the uncertainty on the [NH 2 D]/[ 15 NH 2 D] ratio remains large, the cold prestellar cores L1689N and LDN1544 seem to have a larger ratio than Barnard-1b and NGC1333-DCO + .

15 N fractionation

Nitrogen fractionation involves two main mechanisms in the gas phase: isotopic dependent photodissociation of molecular N 2 , principally at work in the atmosphere of Titan (Liang et al. 2007) and possible ion-molecule fractionation reactions occurring at low temperatures in cold dense cores as first measured by Adams & Smith (1981) and calculated by Terzieva and Herbst (2000). In this latter case, involved endothermicities values range between a few K up to 36 K for exchange reactions involving 15 N, 15 N + , and 15 NN. Selective photodissociation of N 2 and 14 N 15 N takes place at wavelengths between 80 and 100 nm, a range where cold dense cores are completely opaque. Then, this mechanism does not work in the present context. Charley & Rodgers (2002) and Rodgers & Charnley (2008a) have investigated the nitrogen fractionation in their time dependent, coupled gas/solid chemical models. They conclude that 15 N-rich ammonia and deuterated ammonia can be frozen onto the ice mantles, provided all nitrogen is not converted into N 2 . The gas phase becomes enriched at early times, before the complete freezing of the gas phase molecules.

Additional fractionation reactions may be introduced such as those involving 15 N + with CN and NH 3 and some neutral-neutral reactions between 15 N and CN (Rodgers & Charnley 2008b). However, none of these reactions has been studied in the laboratory and these schemes remain highly hypothetical. We have developed a gas phase chemical code, including ion-molecule fractionation reactions for carbon and nitrogen (Langer 1992;[START_REF] Langer | [END_REF]Terzieva & Herbst 2000) as well as a complete deuterium chemistry (Roueff et al. 2005). We have explicitly introduced D and 13 C on the one hand and D and 15 N on the other hand for NH n , HCN and HNC molecules, in order to directly compare the model results with the observations. The chemical network involves 302 chemical species and 5270 reactions. The maximum number of carbon atoms in a molecule has been limited to 3. We have introduced the additional reaction channels arising from the inclusion of isotopic species. We have also preserved functional groups in dissociative recombination reactions such as : HCND + + e → HCN + D

(1)

HCND + + e → DNC + H (2)
Note that the branching ratios of the dissociative recombination of N 2 H + have been measured again by Molek et al. (2007) with the result that the channel towards N 2 occurs with a probability of at least 90% Fig. 2. Prediction of the gas phase abundances relative to H 2 (left) and 14 N/ 15 N abundance ratio (right) for the main nitrogen species. The models assumes a constant temperature of 10K, and increasing depletions with the gas density, to mimic freezing out. The elemental abundance ratio 14 N/ 15 N is set to 400.

A calculation is shown in Figure 2 for typical dense core parameters, and assuming a 14 N/ 15 N abundance ratio of 400, and an ionization rate of ζ = 2 × 10 -17 s -1 . The model predicts that the 15 N enrichment of ammonia is moderate in the gas phase, while a stronger enrichment is predicted for N 2 H + , and depletion for HCN and CN. Recent models by Rodgers and Charnley (2008a) obtain similar results for the gas phase abundances, the 15 N enrichment of ammonia being more efficient in the solid phase.

Conclusions

We report the detection of heavy deuterated ammonia, 15 NH 2 D, in three cold dense cores. The abundance ratio [NH 2 D]/[ 15 NH 2 D] is compatible with the 14 N/ 15 N protosolar value, and seems larger than the terrestrial value despite the remaining measurement uncertainties. While further observations are needed to improve the accuracy and test our chemical models, ammonia and deuterated ammonia seem to be good probes of the 14 N/ 15 N ratio. Deuterated ammonia is particularly interesting as it probes the coldest and densest regions of prestellar cores which are the reservoirs for the future formation of young stars and their associated protoplanetary disks.

Fig. 1 .

 1 Fig.1. Spectra of the 1 1,1 -1 0,1 lines of o-NH 2 D (grey/red line) and o-15 NH 2 D(black line) . The vertical scale is T mb in K, the horizontal scale is V LSR in km s -1 . The o-NH 2 D spectra have been multiplied by 0.02 except for the L1689N spectra that has been scaled by 0.0125.

Table 1 .

 1 Source list

	Source	RA	Dec	V LSR	n(H 2 ) a
		(2000)	(2000)	( km s -1 ) (cm -3 )
	Barnard 1b	03:33:20.9	31:07:34	6.8	3 × 10
	NGC1333-IRAS4A	03:29:10.5	31:13:31	7.2	2 × 10
	NGC1333-DCO +	03:29:12.3	31:13:25	7.2	1 × 10
	L1544	05:04:16.6	25:10:48	7.4	2 × 10
	L134N(S)	15:54:08.6	-02:52:10	2.4	2 × 10
	L1689N	16:32:29.5	-24:28:53	3.8	2 × 10

a From

Caselli et al. (2008)

Table 2 .

 2 Einstein coefficients, upper level energies and critical densities for the range of temperatures considered in this work

	Molecule	Transition	Frequency	Aij	Eup	ncrit
			(GHz)	(s -1 )	(K) (cm -3 )
	o-NH2D	11,1 -10,1	85926.2780	7.82e-6 20.68 4.2 10 6
	o-15

Table 3 .

 3 Line intensities and molecular column densities computed at LTE with the same Tex for o-NH2D and o-15 NH2D. Tex is derived from the HFS fit of the o-NH2D profile.

				o-NH2D				o-15 NH2D
	Source	T mb ± σ a K	δV kms -1	τ	Tex K	N b 10 14 cm -2	T mb ± σ a mK	I mKkms -1 kms -1 10 11 cm -2 δV N b	[	[NH2D] 15 NH2D]
	Barnard1b N1333-IRAS4A 1.0 ± 0.018 1.38 1.39 ± 0.10 5.0 ± 0.5 2.7 ± 0.6 2.5 ± 0.047 0.79 5.24 ± 0.14 6.0 ± 0.5 4.7 ± 0.5	42 ± 9 ±10	30 ± 4 < 30	0.67 ...	10 ± 2.7 < 10	470 +170 -100 > 270
	N1333-DCO + LDN1544	1.3 ± 0.015 1.15 1.71 ± 0.05 5.3 ± 0.5 2.4 ± 0.4 2.3 ± 0.016 0.47 7.05 ± 0.05 5.3 ± 0.5 4.1 ± 0.5	26 ± 8 ±7	14 ± 3 < 10	0.52 ...	6.7 ± 2.5 < 5.2	360 +260 -110 > 700
	L134N(S) L1689N	2.2 ± 0.033 0.42 4.75 ± 0.10 5.5 ± 0.5 2.4 ± 0.4 5.3 ± 0.030 0.53 6.98 ± 0.02 8.5 ± 0.5 3.4 ± 0.5	24 ± 7 65 ± 17	10 ± 2 26 ± 6	0.40 0.37	4.5 ± 2 4.2 ± 1.5	530 +570 -180 810 +600 -250
	a σ is the rms computed for the original spectral resolution of 40 kHz = 0.136 km s -1 .

b

Available at http://library.osu.edu/sites/msa/suppmat/ v251.i1-2.pp90-101/mmc1.txt

Acknowledgements. We thank the IRAM director for assigning additional time for this program, which helped us to confirm the 15 NH 2 D detection and the 30m staff for their support during the observations. We thank the referee, E. Bergin, for his insightful comments. We acknowledge financial support from the CNRS interdisciplinary program "Origines des Planètes et de la Vie", and from the INSU/CNRS program PCMI. NM is supported by Spanish MICINN through grants AYA2006-14876, by DGU of the Madrid community government under IV-PRICIT project S-0505/ESP-0237 (ASTROCAM), and by Molecular Universe FP6 MCRTN." DCL is supported by the NSF, grant AST-0540882 to the Caltech Submillimeter Observatory.


João Pedro Norenberg
email: jp.norenberg@unesp.br

Americo Cunha
email: americo.cunha@uerj.br

Samuel Da Silva
email: samuel.silva13@unesp.br

Paulo Sergio Varoto
email: varoto@sc.usp.br

Global sensitivity analysis of asymmetric energy harvesters

Keywords: energy harvesting, nonlinear dynamics, sensitivity analysis, Sobol' indices, polynomial chaos expansion

Parametric variability is inevitable in actual energy harvesters. It can significantly affect crucial aspects of the system performance, especially in harvesting systems that present geometric parameters, material properties, or excitation conditions that are susceptible to small perturbations. This work aims to develop an investigation to identify the most critical parameters in the dynamic behavior of asymmetric bistable energy harvesters with nonlinear piezoelectric coupling, considering the variability of their physical and excitation properties. For this purpose, a global sensitivity analysis based on orthogonal variance decomposition, employing Sobol indices, is performed to quantify the effect of the harvester parameters on the variance of the recovered power. This technique quantifies the variance concerning each parameter individually and collectively regarding the total variation of the model. The results indicate that the frequency and amplitude of excitation, asymmetric terms and electrical proprieties of the piezoelectric coupling are the most critical parameters that affect the mean power harvested. It is also

Introduction

Novel technological applications that use small-scale autonomous devices continue to emerge every year. In this context, one of the most significant challenges in developing a self-sufficient machine, i.e., an electromechanical system capable of generating electrical power to sustain its operation or provide additional energy supply. In this perspective, energy harvesting has become a promising research area that comprises dispersed energy in the environment (e.g., solar, wind, heat, vibration, etc.) to convert into electricity using a suitable conversion process.

Among commonly available energy sources present in the environment, the kinetic energy from environmental noise and structural vibration signals, usually wasted, has been major explored for energy harvesting due to a good compromise between energy availability and ease of use. The literature has several recent works with powerful applications involving this idea, such as the generation of bio-cellular energy [START_REF] Catacuzzeno | Energy harvesting from a bio cell[END_REF], medical implants [START_REF] Karami | Powering pacemakers from heartbeat vibrations using linear and nonlinear energy harvesters[END_REF], sensor power [START_REF] Lee | Vortex-induced vibration wind energy harvesting by piezoelectric mems device information[END_REF][START_REF] Wang | A novel humidity resisting and wind direction adapting flag-type triboelectric nanogenerator for wind energy harvesting and speed sensing[END_REF][START_REF] Yi | A highly shape-adaptive, stretchable design based on conductive liquid for energy harvesting and self-powered biomechanical monitoring[END_REF], and even osmotic energy [START_REF] Xin | High-performance silkbased hybrid membranes employed for osmotic energy conversion[END_REF]. The conversion techniques used in vibration energy harvesting are electrostatic [START_REF] Mitcheson | Mems electrostatic micropower generator for low frequency operation[END_REF], electromagnetism [START_REF] Arnold | Review of microscale magnetic power generation[END_REF], and piezoelectricity [START_REF] Erturk | An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations[END_REF]. Piezoelectricity, in particular, has received more attention because of its simplicity and the high-density energy of the resulting electrical signals when the piezoelectric material is used in the sensing mode of operation, which is typically the case in piezoelectric energy harvesting systems.

In a broad sense, energy harvesting devices are modeled as deterministic systems, and, consequently, all corresponding model parameters are free of uncertainty. However, even the most accurate manufacturing processes can present some degree of geometrical variability in large-scale production. Additionally, material property variability can also occur. These, along with other environmental factors, are sources of uncertainty and can affect the overall performance of a given energy harvesting system. Thus, disregarding these sources of uncertainties represents a significant limitation of the deterministic approach, leading to low fidelity or even wrong predictions. Nowadays, all engineering fields that deal with robust and accurate predictions require suitable methods to characterize and quantify these uncertainties, and the general area from applied mathematics and statistics that studies the effect of uncertainties is usually known as Uncertainty Quantification (UQ) [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. This relevant research relies on stochastic techniques to determine the statistical output distribution given the corresponding system's input distribution. UQ is essential in energy harvesting since it can handle low-scale and accurate models. Recent reviews of UQ in the context of energy harvesting are found in [START_REF] Franco | Parameter uncertainties in the design and optimization of cantilever piezoelectric energy harvesters[END_REF][START_REF] Li | Uncertainty analysis of bistable vibration energy harvesters based on the improved interval extension[END_REF][START_REF] Mann | Uncertainty in performance for linear and nonlinear energy harvesting strategies[END_REF].

Sensitivity analysis (SA) is a tool that has been recently receiving increasing attention in the analysis of many engineering systems since it offers a comprehensive evaluation of how the system properties can affect the system's response [START_REF] Abbiati | A global sensitivity analysis framework for hybrid simulation[END_REF]. It allows for identifying the most critical parameters that affect the energy harvesting dynamic response. According to [START_REF] Saltelli | Sensitivity analysis[END_REF], SA aims to determine which factors require a more precise estimation and to identify the weak links in the evaluation chain (those that propagate the most significant variance in the output). Besides assisting UQ techniques in formulating a more straightforward probabilistic model for the system of interest, this analysis can also bring further insight into the system's behavior.

SA techniques can be classified as local and global [START_REF] Cacuci | Sensitivity and uncertainty analysis: theory[END_REF]. The first is based on the approximation of the partial derivative to assess how the variation in one parameter affects the system's response while keeping other parameters fixed according to their nominal values. The second approach provides a comprehensive analysis, evaluating the effect of simultaneously varying the parame-ters within the multidimensional input space. A helpful survey on the state-of-the-art sensitivity analysis techniques is available at [START_REF] Saltelli | Sensitivity analysis[END_REF].

Applications of SA on energy harvesting can be seen, for instance, in Aloui, Larbi, and Chouchane [START_REF] Aloui | Global sensitivity analysis of piezoelectric energy harvesters[END_REF][START_REF] Aloui | Uncertainty quantification and global sensitivity analysis of piezoelectric energy harvesting using macro fiber composites[END_REF], where the authors applied a global sensitivity analysis (GSA) on a linear bimorph piezoelectric energy harvester subjected to a base excitation with a mounted load resistance in series. Ruiz and Meruane [START_REF] Ruiz | Uncertainties propagation and global sensitivity analysis of the frequency response function of piezoelectric energy harvesters[END_REF] conducted GSA based on a variance method concerning the behavior of the frequency response function of the linear energy harvester under unimorph and bimorph piezoelectric configurations. These studies concluded that the effects of physical properties of the piezoelectric materials are fundamental for the understanding of variabilities in the models, and they draw attention to the need to investigate the propagation of uncertainties in piezoelectric energy harvesters. However, to the author's knowledge, there is a gap in the literature for sensitivity analysis in nonlinear energy harvesting systems that could provide in quantitative terms which parameters most affect the variability of the system's response.

Linear models have limited bandwidth around their fundamental resonance frequency. This limitation can be overcome by exploiting nonlinearities, which allow energy scavenges in broadband frequency. Cottone et al. [START_REF] Cottone | Nonlinear energy harvesting[END_REF] and Erturk et al. [START_REF] Erturk | A piezomagnetoelastic structure for broadband vibration energy harvesting[END_REF] were the first works to implement nonlinear bistable energy harvesting systems. They demonstrated a substantial energy increase compared to the linear system in a broadband scenario. This energy harvester has become a classic system, being studied by several authors [START_REF] Cao | Nonlinear time-varying potential bistable energy harvesting from human motion[END_REF][START_REF] Cunha | Enhancing the performance of a bistable energy harvesting device via the cross-entropy method[END_REF][START_REF] Daqaq | On the efficacy of charging a battery using a chaotic energy harvester[END_REF][START_REF] Daqaq | On the role of nonlinearities in vibratory energy harvesting: A critical review and discussion[END_REF][START_REF] Erturk | Broadband piezoelectric power generation on high-energy orbits of the bistable duffing oscillator with electromechanical coupling[END_REF][START_REF] Huang | Nonlinear analysis of multistable energy harvesters for enhanced energy harvesting[END_REF][START_REF] Lopes | Nonlinear characterization of a bistable energy harvester dynamical system[END_REF][START_REF] Yang | Investigation of coupled leverbistable nonlinear energy harvesters for enhancement of inter-well dynamic response[END_REF] due to its rich dynamics and complexity. In addition to these studies, nonlinear energy harvesting systems can be exploited via SA, which can be seen as an exploratory tool capable of revealing individual and joint parametric contributions to the system response. Since nonlinear systems present high sensitivity to input parameter variations and may have complex dynamics with regular and chaotic behaviors and multiple solutions, SA may reveal nontrivial input-output relationships that can improve the energy recovery process.

This work is concerned with the global sensitivity analysis of bistable energy harvesting systems. Several sources of uncertainties and nonlinearities are considered to explore the physical system behavior, such as multiple excitation conditions and geometric asymmetries. Those asymmetries can arise from the manufacturing process, the buckling eccentricity, magnetic forces, and heterogeneity of the materials used. Therefore, accounting for these effects leads to more reliable models for assessing the system's performance in the electromechanical conversion process. Hence, a prob-abilistic model of uncertainties employing a suitable joint probability distribution for the model parameters is conceived. The goal here is to classify which selected parameters most affect the harvested power. The model allows a clear understanding of the physical phenomena involved in the bistable energy harvesting system dynamics. Mainly through the interactions between the corresponding variables, which is an essential step towards reaching robustness of the model employed, serving as a training tool to determine the essential features closely linked to the harvesting output response data.

The remaining parts of this manuscript are organized as follows. Section 2 presents the energy harvesting system of interest. The global sensitivity analysis of this nonlinear system is formulated in Section 3, and the Polynomial Chaos Expansion approach used to compute the Sobol indices is presented in Section 4. In Section 5, the numerical experiments conducted to analyze the sensitivity of a bistable energy harvester are presented and discussed. Section 5.5 summarizes the main results obtained here and the uncertainties propagation of the systems. Finally, in Section 6, important conclusions are presented.

Bistable energy harvesting system

An illustration of the energy harvesting system of interest in this work can be seen in Fig. 1, which consists of a rigid base excited by a harmonic force, a vertical fixedfree beam made of ferromagnetic material, with two permanent magnets located near the bottom part of the base. In the high deformation region of the beam's upper part, a pair of piezoelectric layers is coupled to a resistive circuit responsible for converting the kinetic energy into an electrical potential, which is dissipated in the resistor. The harvester has a bias angle, leading to oscillations around a certain angle φ. Parameter ϑ can be adjusted for different configurations in order to produce asymmetric potential energy functions.

Disregarding the asymmetric terms, the system is reduced to the classic piezo-magneto-elastic beam energy harvester proposed by Erturk, Hoffmann, and Inman [START_REF] Erturk | A piezomagnetoelastic structure for broadband vibration energy harvesting[END_REF]. Besides, this classic model uses a linear relationship between mechanical and electrical domains. Then, classic bistable energy harvester electromechanical equations are stated as

ẍ + 2 ξ ẋ - 1 2 x (1 -x 2 ) -χ v = f cos (Ω t) , (1) v
+ λ v + κ ẋ = 0 , (2)
x(0) = x 0 , ẋ(0) = ẋ0 , v(0) = v 0 , (3)
where x is the beam tip displacement; ξ is the damping ratio; χ is the piezoelectric coupling term in the mechanical equation (depending on the piezoceramic properties and the structure mass); v is the voltage; λ is a reciprocal time constant (λ ∝ 1/R l C p , where R l is the load resistance and C p is the equivalent capacitance of the piezoceramic layers); κ is the piezoelectric coupling term in electrical equation (depending on the piezoceramic properties and the equivalent piezoelectric capacitance); f is the external excitation amplitude; Ω is the external excitation frequency. The initial conditions are x 0 , ẋ0 and v 0 , which, respectively, represent the initial values of the beam edge position, velocity, and voltage through the resistor. All these quantities are dimensionless.

The main quantity of interest (QoI) associated with the dynamical system under analysis is the mean output power given by

P = 1 T t0+T t0 λ v(t) 2 dt , (4)
which is the temporal average of the instantaneous power λ v(t) 2 over a given time interval of length T . There are sources of nonlinearities that are unavoidable when building an energy harvester, for instance, the symmetrical position of the magnets. A quadratic nonlinearity is introduced to the nonlinear restoring force to characterize the asymmetric potentials [START_REF] Halvorsen | Fundamental issues in nonlinear wideband-vibration energy harvesting[END_REF][START_REF] He | Influence of potential function asymmetries on the performance of nonlinear energy harvesters under white noise[END_REF]. Then, the restoring force is described in linear, quadratic, and cubic terms.

Another source of asymmetry often encountered in actual energy harvesters is the sloping plane where the system is attached, known as bias angle φ. In many energy harvesting applications, the system is subjected to rotational movements that cause dynamic asymmetry. Because of the action of gravity, this system suffers changes in its inertial configurations. The resulting bias angle is used in [START_REF] Wang | Nonlinear dynamics and performance enhancement of asymmetric potential bistable energy harvesters[END_REF] to enhance the harvester's performance so that it is added to the dynamic model as an external force from the self-weight of the harvester.

The nonlinear effect between the material strain and the piezoelectric coefficient can be seen in piezoceramic material according to [START_REF] Crawley | Detailed models of piezoceramic actuation of beams[END_REF]. In high deformation conditions, e.g., resonance regions, ignoring the nonlinear effects may cause an underestimation of harvested power [START_REF] Dutoit | Experimental verification of models for microfabricated piezoelectric vibration energy harvesters[END_REF]. Therefore, models and studies considering a nonlinear electromechanical coupling for power harvesting are explored. Triplett and Quinn [START_REF] Triplett | The effect of nonlinear piezoelectric coupling on vibration-based energy harvesting[END_REF] modeled the piezoelectric coupling by a first-order nonlinear relationship between material deformation and piezoelectric coefficient. In contrast, [START_REF] Leadenham | Unified nonlinear electroelastic dynamics of a bimorph piezoelectric cantilever for energy harvesting, sensing, and actuation[END_REF][START_REF] Stanton | Nonlinear piezoelectricity in electroelastic energy harvesters: Modeling and experimental identification[END_REF] also included higher-order terms on the piezoelectric constitutive equation.

In this way, to obtain more accurate estimates of the energy harvested by the system, it is also appropriate to consider a nonlinear electromechanical coupling. According to Triplett and Quinn [START_REF] Triplett | The effect of nonlinear piezoelectric coupling on vibration-based energy harvesting[END_REF] the nonlinear coupling is determined through a dependent function on the linear coupling term θ L , nonlinear coupling term θ N L , and the deformation of the piezoelectric material (s). It follows that, the dimensional piezoelectric coefficient, obtained by [START_REF] Triplett | The effect of nonlinear piezoelectric coupling on vibration-based energy harvesting[END_REF], which is represented as

Θ(s) = θ L (1 + θ N L |s|) . (5)
The following initial value problem describes the governing lumped-parameter equation of motion for the proposed bistable energy harvesting model taking into account asymmetries and nonlinear electromechanical coupling

ẍ + 2 ξ ẋ - 1 2 x (1 + 2δx -x 2 ) -(1 + β |x|)χ v = f cos (Ω t) + p sin φ , (6)
v + λ v + (1 + β |x|)κ ẋ = 0 , (7)
x(0) = x 0 , ẋ(0) = ẋ0 , v(0) = v 0 , (8)
where δ is a coefficient of the quadratic nonlinearity, p is the equivalent dimensionless constant of gravity of ferromagnetic beam, and β is the dimensionless nonlinear coupling term. Figure 2 depicts the displacement time-series for the classical bistable oscillators, with nonlinear electromechanical coupling, asymmetric potential, and bias angle under different amplitudes of excitation. These bistable oscillators have three equilibrium points: two stable (central region of each magnet) and one unstable (between the magnets). Thus, it presents three distinct dynamic behaviors, which mainly depend on the excitation amplitude. The lower energy orbits, i.e., oscillations only at one stable equilibrium point, are known for single-well motion. On the other side, the high-energy orbits, where the system oscillates between the equilibrium points (inter-well motion/ snap-through behavior), can be characterized as periodic and chaotic behaviors. In energy harvesting, vibrations at high-energy orbits can improve the system's performance. This figure also displays that piezoelectric nonlinearities, asymmetric geometry, and bias angle can change the dynamic behavior compared to the classical case, mainly in the unstable region with a low excitation amplitude. In addition, the high sensitivity of this type of system is observed. An animation of these systems is available in the Supplementary Material and also in the STONE-HENGE code repository [START_REF] Norenberg | STONE-HENGE -Suite for Nonlinear Analysis of Energy Harvesting Systems[END_REF], which is shown in detail each behavior for these systems over time. The playlist at [START_REF] Norenberg | STONE-HENGE -Suite for Nonlinear Analysis of Energy Harvesting Systems[END_REF] also presents the linear and nonlinear harvesters' dynamic animations.

Global sensitivity analysis

This section presents the global sensitivity analysis framework to study the bistable energy harvester dynamics. Relevant insights are provided into the system in question that allows model reduction schemes; for instance, it builds a simplified model by fixing relatively unimportant parameters to their nominal values and varying only the most influential parameters [START_REF] Nagel | Principal component analysis and sparse polynomial chaos expansions for global sensitivity analysis and model calibration: Application to urban drainage simulation[END_REF]. The global sensitivity analysis method adopted is a variance decomposition that writes the dynamical system response as a sum of each input variable's individual and combined contributions. The Sobol indices, established in [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] quantifies the contribution of each parameter concerning the total variance of the model. It has recently been used in many works, with high impact research, as in [START_REF] Alemazkoor | Hurricane-induced power outage risk under climate change is primarily driven by the uncertainty in projections of future hurricane frequency[END_REF][START_REF] Lund | Global sensitivity analysis for the design of nonlinear identification experiments[END_REF]. One of its main advantages is dealing with nonlinear and non-parameterized models and providing a quantitative and qualitative classification.

To facilitate mathematical development, consider that the QoI from the dynamical system of interest can be represented by the following functional relationship

Y = M(X) , X = {X  , X  , . . . , X k }, ( 9
)
where X is an input vector with k independent parameters, which is modified by the mathematical operator M, to produce the scalar output (the QoI) Y . Using the Hoeffding-Sobol decomposition [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF][START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF], the QoI can be written as a decomposition into sum-

Y = M  + k i= M i (X i ) + i<j M ij (X i , X j ) + • • • + M ...k (X  . . . X k ) , (10)
which has finitely many sums with finite variance, and according to Hommaa and Saltelli [START_REF] Homma | Importance measures in the global sensitivity analysis of nonlinear models[END_REF] it has an orthogonal decomposition property in terms of conditional expectations

Dx M u (x u ) M v (x v ) dx = 0, ∀ u = v (11)
where u ⊂ {1, . . . , k} and x u def = {x i , . . . , x is }. In order to simplify the notation, the input variables are assumed to follow uniform distributions with a normalized support [0, 1], so that the support of the random input vector is D X = [0, 1] k . Note that it is always possible to normalize a finite support random variable. Thus, all the sums in the expansion can be recursively computed. The first term M  is a constant equal to the expected value of Y , i.e.,

M  = D X M (X) dX . (12)
The M i (X i ) and M ij (X i , X j ) terms are, respectively, the conditional expected value of M given the parameters i and combinations of the parameters ij (i = j), written as

M i (X i ) =   • • •   M (X) dX ∼i -M  , (13)
M ij (X i , X j ) =   • • •   M (X) dX ∼ij -M  - M i (X i ) -M j (X j ) , ( 14
)
where the notation ∼ indicates which variables are excluded from the integration. This allows the meaningfully partition the system response variance and determine the Hoeffding-Sobol decomposition in terms of variance [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. Considering independent input variables, we can quantify the contribution of the individual X i and any combination X ij , . . . , X i...k of variables to the total variance V ar[Y ]. By construction

k i=1 S i + i<j S ij + • • • + S ...k = 1 , (15)
where the indices S are called as Sobol indices and the sum of all individual and joint indices is equal to the unity.

The first-order Sobol indices

S i = V ar[M i (X i )] V ar[M(X)] , (16)
quantify the additive effect of each input separately concerning the total variance, and the second-order Sobol indices S ij

S ij = V ar[M ij (X i , X j )] V ar[M(X)] , (17)
compute the joint-effects of two inputs. Higher-order Sobol indices consider the interaction effects of several parameters and are written in the same way.

Polynomial chaos expansion

Generally, Sobol indices are calculated using the Monte Carlo (MC) method, but this can be computationally expensive due to the slow convergence rate of this sampling method [START_REF] Cunha | Uncertainty quantification through the monte carlo method in a cloud computing setting[END_REF][START_REF] Kroese | Handbook of Monte Carlo Methods[END_REF]. An alternative way to obtain Sobol indices is to employ the Polynomial Chaos Expansion (PCE) surrogate model, which is an approximate model for the original system with a short processing time, reducing the computational cost and keeping calculation accuracy [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF][START_REF] Palar | Global sensitivity analysis via multi-fidelity polynomial chaos expansion[END_REF][START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. PCE was introduced in the engineering community by R. Ghanem, and P. Spanos [START_REF] Ghanem | Stochastic finite elements -a spectral approach[END_REF] in the context of stochastic finite element analysis. The PCE is used for UQ in various domains, e.g., in solid mechanics, thermal and fluid sciences, etc. It can be a non-intrusive method representing uncertain quantities as an expansion, including the decomposition of deterministic coefficients and orthogonal polynomial bases for each random variable probability distribution [START_REF] Xiu | Numerical Methods for Stochastic Computations: A spectral method approach[END_REF]. PCE offers an accurate and efficient way to include nonlinear effects in stochastic analysis and can be seen as an ideal mathematical way to construct and obtain a model response surface in the form of a high-dimensional polynomial in model parameters uncertain [START_REF] Oladyshkin | Data-driven uncertainty quantification using the arbitrary polynomial chaos expansion[END_REF][START_REF] Sepahvand | Uncertainty quantification in stochastic systems using polynomial chaos expansion[END_REF].

Assuming that Y = M(X) is a finite variance random variable, defined in terms of the composition of the random vector X and the operator M, one can write the following PCE

Y ≈ α∈A y α ψ α (X) , (18)
where, ψ α are multivariate orthonormal polynomials concerning the joint probability density function (PDF) f X of X, y α are unknown deterministic coefficients, and the truncation set A ⊂ N M is selected from all possible multi-indices of multivariate polynomials. The unknown coefficients can be obtained from a non-intrusive least-squares regression technique, where samples from the system response are obtained from the full dynamics model. Statistical moments and sensitivity analysis can be computed efficiently via post-processing of the PCE coefficients. Due to orthonormality property of the PCE basis, and the fact that ψ 0 ≡ 1 and E [ψ α (X)] = 0 ∀ α = 0, the mean value and the variance of the system response can estimated as

E [Y ] ≈ y 0 , (19)
and

E (Y -E [Y ]) 2 ≈ α =0 α∈A y 2 α . (20)
Similarly, the Sobol indices of any order can also be determined analytically, e.g., the first-order indices are derived as follows

S i ≈ α =0 α∈Ai y 2 α α =0 α∈A y 2 α , (21)
where the PCE coefficients for the numerator are restricted to the subset A i , according to the parameter i, to compute the conditional variance. The right side of the equation is an approximation due to the truncation scheme used to compute the PCE. Similarly, the second-order Sobol indices are stated as

S ij ≈ α =0 α∈Aij y 2 α α =0 α∈A y 2 α . (22)
Figure 3 outlines a summary of the methodology employed to assess the sensitivity of the system response regarding variations in the parameters.

Results and discussion

Numerical experiments reported in this work adopt the following physical parameters: ξ = 0.01, χ = 0.05, λ = 0.05, κ = 0.5, Ω = 0.8 and f is varied. The initial condition is defined by (x 0 , ẋ0 , v 0 ) = (1, 0, 0). The dynamics is integrated by Runge Kutta scheme (4th order) over the time interval 0 ≤ t ≤ 2000 with relative tolerance of 10 -6 , and absolute tolerance of 10 -9 . The mean output power is computed over the last 50% of these time-series.

Global sensitivity analysis assumes that all the system parameters are independent and uniformly distributed over given intervals, defined by a coefficient of variation (c.v.) of ± 20% around the nominal values. We approached different models to verify specific behaviors for each condition: The classical model without nonlinear electromechanical coupling and asymmetries is evaluated. In the following, the nonlinear electromechanical coupling is included. The asymmetry bistable model with linear piezoelectric is studied. Finally, a proposed model with nonlinear coupling and asymmetries is approached.

Classical bistable energy harvester

The classical bistable oscillator is sensitive to excitation conditions as shown in [START_REF] Lopes | Nonlinear characterization of a bistable energy harvester dynamical system[END_REF]. So, the system's behavior through bifurcation diagrams with nominal parameters is assessed. The goal is to identify regions of the amplitude and frequency of excitation that imply distinct nominal dynamic behavior, being an essential step in the sensitivity analysis. Figure 4a shows the bifurcation diagrams of the classical bistable energy harvester. In the bifurcation diagram, sweeping the excitation amplitude to five frequency values established within the variance range is performed. For excitation frequencies of 0.64 and 0.72, the system presents a periodic behavior for the entire amplitude range. It is possible to identify low-energy orbit regions for lower amplitude values and high-energy orbits for larger values through a phase portrait and Poincare section in Fig. 4b. Alternatively, excitation frequencies of 0.8, 0.88, and 0.96 can lead to chaotic behavior and inter-well motion.

Two different dynamic behaviors from the bifurcation diagram are chosen to validate the metamodel, testing distinct scenarios. The first-order Sobol indices are calculated using the Monte Carlo method. The surrogate PCE model is constructed and tuned using Monte Carlo simulations as a reference. It is worth highlighting that the PCE model requires verification from a reference to guarantee the reliability of convergence and accuracy, as shown in the schematic representation of the methodology, Fig. 3. Finally, the accurately PCE model is used to explore several scenarios of excitation ampli-tude with low-cost processing. To illustrate this process, is selected nominal values f = 0.147 and Ω = 0.8, as a periodic behavior and f = 0.083 and Ω = 0.8, as a chaotic one.

The first-order Sobol indices for the periodic case are shown in Fig. 5a using MC and PCE methods. These indices reveal that 40% of the system sensitivity is due to the piezoelectric coupling acting on the electrical equation κ, then 30% for the frequency Ω, 10% for amplitude f , and 10% for the reciprocal time constant λ. The influence of the piezoelectric coupling on the mechanical system χ and the damping ratio ξ did not cause variability in power harvested.

For a second case, dealing with nominal chaotic behavior, the first-order Sobol indices are shown in Fig. 5c based on MC and PCE methods. For this scenario, the frequency and amplitude are powerful for energyharvested variability, supplying about 40% and 20% of the total variance, respectively. The piezoelectric coupling on the electrical equation has 5%; finally, the reciprocal time constant, piezoelectric coupling on the mechanical equation, and damping have negligible values.

To compare the employed methods, in Tab. 1 shows the Sobol Indices performance between MC and PCE for the periodic and chaotic case. Note that PCE presents a shorter time processing than MC. The mean absolute errors (MAE) for Sobol Indices between the methods are low for both cases. Besides, the Supplementary Material presents the PCE approach by PDF compared to the complete model response, showing good performance. Then, PCE is accurate and can substitute precisely our model, becoming workable to explore several behaviors and high orders of Sobol indices. Figure 5b presents the second-order indices for the first case, where they are calculated only by the validated PCE method. This result shows that only the combined effect of the excitation parameters is relevant, meaning about 10% remaining of the total variance of the system. For higher orders, their evaluations are un-necessary since they have accounted for (almost) the distribution of 100% of the total variance. The secondorder Sobol indices for the second case are shown in Fig. 5d and demonstrate that only the joint effect of excitation parameters Ωf is relevant, about 20%. The increase of the Ω and f effects occur since chaotic behavior is more sensitive to excitation conditions, as proven in first-order Sobol indices.

A wide range of excitation amplitudes using Sobol indices based PCE is explored. Sobol indices are calculated for each excitation amplitude. The indices are presented by a stacked bar plot in Fig. 5e. All the firstorder indices and the relevant second-order indices are inserted. The Ω, f , and κ effects are the parameters that most cause variability in harvested power because they present larger indices for all ranges. In regions with a wide range of stability, where the system has enough energy to oscillate in inter-well motion -high excitation amplitude -the sensitivity to κ is the most relevant. However, the effects of Ω and f increase in the unstable condition, where the system can change the dynamic behavior for a lower amplitude of excitation (as seen in Fig. 4). Finally, Fig. 5e also shows that the higher-order terms have a negligible effect, because the first-and second-order indices account for almost all variance: the bars are reaching close to 1 (see Eq. 15).

The effects of frequency and amplitude of excitation are significant because their variability can change dynamic behaviors, which is seen in the bifurcation diagram under the c.v. = ± 20% at high frequency. The dynamic behavior changes directly affect the power harvested, e.g., high-energy orbits produce much more electricity than chaotic or low-energy orbits. Although the frequency and amplitude of excitation are mainly responsible for behavior changes, the piezoelectric coupling presents greater sensitivity for high amplitude. This is due to the small probability of behavior variation, so the piezoelectric effect predominates. Therefore, whatever scenario adjacent to the orbit changes, f , and Ω are the most influential parameters; otherwise, κ is more relevant.

Classical bistable energy harvester with nonlinear electromechanical coupling

The nonlinear electromechanical coupling is added in the classical bistable energy harvester. First, a parametric investigation of nonlinear electromechanical coupling term β and amplitude excitation range through the mean power performance is performed. Fig. 6 shows a contour map for the mean power over a wide range for the excitation amplitude and the nonlinear electromechanical coupling. The system can undergo an increase in the harvested power as β increases. However, there is a threshold value for the nonlinear term, damaging system performance from that point on. The high values of the nonlinear coupling term harm the power harvested for low excitation amplitude. Besides, in instability regions (chaotic behavior), the system has more significant mean power harvested changes, noticing a discontinuity color region between f = 0.07 to 0.105 (chaotic area). This situation exemplifies how the nonlinear electromechanical coupling requires careful investigation, leading to a more complex analysis.

Figure 7 shows the bifurcation diagram sweeping up and down the amplitude of excitation for different nonlinear coupling values (β = 0.5, 1.0, 1.5, 2.0, 2.5 and 3.0). For low values of the nonlinear coupling, the system maintains almost the same behavior as the absolute amplitude. The nonlinear coupling only decreases the chaotic interval, shifting the chaos onset to higher excitation amplitudes. However, the regions with chaotic intervals vanish as β increases, and the dynamic behavior acquires a large range of low-energy orbits, as seen in Fig. 6 with a large area of low energy as the nonlinear coupling increases (dark blue region). This effect occurs because the nonlinear coupling acts as damping in the system, as explained in [START_REF] Triplett | The effect of nonlinear piezoelectric coupling on vibration-based energy harvesting[END_REF]. The nonlinear piezoelectric configuration harvests more energy in high-energy Fig. 6: Contour map of the mean output power over a wide range of excitation amplitude and nonlinear electromechanical coupling.

orbits, but it demands more external energy for the mechanical system to perform a snap-through behavior.

In Fig. 8, the leading Sobol indices are shown for different β values for several amplitudes excitation scenarios after obtaining a PCE surrogate in the same way as done in the classical case. More detail on PCE accuracy can be seen in the Supplementary Material. The sensitivity to the nonlinear piezoelectric term has a minor influence. However, the higher β values (β ≥ 2) change the individual and combined effect of other parameters, i.e., β indirectly affects the system's sensitivity. The effect of the excitation frequency is dominant throughout the amplitude spectrum. The excitation amplitude has a slight drop in influence, but its contribution is consistently significant in the chaotic regions for low values of β. The electromechanical coupling in the electrical equation also loses some influence, mainly because the nonlinear term drives the system more sensitive to excitation conditions. Higher β values require more external energy to perform interwell motions and achieve a stable region. Also, there is a significant increase in the combined influence of Ω β and ξΩ, previously insignificant. It concludes that when the piezoelectric nonlinearity is high, the system may change its dynamics, affecting its sensitivity and indirectly increasing and decreasing the individual and combined influence of other parameters. These results emphasize that information about nonlinear coupling is crucial. High values for the piezoelectric nonlinearity may change the dynamics, affecting its sensitivity and indirectly interfering with other parameters' individual and joint influence. On the other hand, it is unnecessary to worry about its uncertainties for low nonlinearity. Therefore, we recommend considering the uncertainties for high nonlinearities.

Asymmetric bistable energy harvester with linear electromechanical coupling

This section treats the classical bistable energy harvesting system considering asymmetries. The asymmetric potential, as previously discussed, is defined by a quadratic term, but its coefficient depends on the different variations of each design. The asymmetric term is considered through a random variable described by a uniform distribution ranging from -0.15 to 0.15. Figure 9 shows the potential energy of the system for different δ values. The δ changes the equilibrium point of the model, causing asymmetry in the potential. Left asymmetry occurs for positive values and right asymmetry for negative values. Then, the sensitivity of the system considering the energy potential varying in this uncertainty range is analyzed.

First, the Sobol analysis is performed for the potential asymmetric configuration of the model, disregarding the bias angle, i.e., φ = 0. Figure 10 presents the Sobol indices for the entire spectrum of amplitude. The quadratic nonlinearity term does not cause changes in the system's sensitivity. The influence of the excitation frequency at low amplitudes and the piezoelectric term at high amplitude influences are dominant. A noticeable portion of the variance comes from the secondorder index related to the combined effect of δ and the excitation frequency. However, only the frequency has a significant first-order contribution.

Next, the bias angle is considered in the model. The bias angle φ is also described as a random variable with uniform distribution going from -15 • to 15 • . Positive angles show that the system undergoes a counterclockwise rotation, and the reverse occurs for negative angles. Figure 11 shows the main Sobol indices. The bias angle significantly changes the system's sensitivity, especially for low amplitude cases, when the model exhibits unstable dynamic behavior. At low amplitudes, the model is more sensitive to variations in the bias angle φ, excitation frequency f , and the joint effect (higher-order) from combinations of the parameters with the highest first-order indices. The frequency effect drops because of the bias angle higher effects, which become dominants. At high amplitude regions, the system increases its sensitivity to the piezoelectric coupling κ, but the asymmetry term δ keeps being irrelevant and can be its variability disregarded. Therefore, besides the excitation parameters and electrical piezoelectric properties, the bias angle also becomes an important parameter that can cause variability in the model.

Asymmetric bistable energy harvester with nonlinear electromechanical coupling

Finally, the proposed energy harvesting system considering asymmetry, bias angle, and nonlinear coupling is investigated. The range values of asymmetric coefficient and the bias angle are adopted as previously done. The nominal value of the nonlinear coupling is set equal to 1.0. Figure 12 nonlinear coupling does not demonstrate first or secondorder influences for this scenario. Variability from the asymmetry coefficient, damping, and mechanical coupling terms do not affect the harvested power, so they might be deterministic terms. As well as in the previous cases, for low excitation amplitude, the system is sensitive to parameters that change its behavior, driven by the excitation frequency, amplitude, and bias angle.

For the cases of the high amplitude of excitation, the first-order effect of the piezoelectric coupling is significant. Besides, the excitation frequency becomes the pa- rameter with the greatest influence, progressively overcoming the influence of the piezoelectric coupling. This behavior occurs because of the nonlinear piezoelectric coupling, which demands more external energy for the system to achieve higher orbits, as already discussed. Thus, the frequency variations become very important in changing energy orbits, confirmed by the increase in its Sobol index. Therefore, this complete nonlinear system is still sensitive to parameters that modify the system's orbit. The analysis is complex since the piezoelectric nonlinearity can quantitatively change the other parameters' sensitivity. However, its uncertain effects are only significant for high amounts.

Summary of the findings

Table 2 summarizes the significant results obtained in this study. The results showed that the parameters with high first-order indices control the higher-order sensitivity terms. These results suggest that parameters with low sensitivity indices concerning the harvested power can be treated as fixed variables, reducing the dimension of the probabilistic model -these being the damping ξ, piezoelectric coupling in the mechanical equation χ, and asymmetric coefficient δ. In addition, identifying the parameters with low sensitivity indices allows for manufacturing constraints to be relaxed.

Parameters with medium and high sensitivity have a significant influence. Given the conditions, they should be treated as stochastic variables: the piezoelectric coupling in the electrical equation κ, external conditions Ω and f , bias angle φ, nonlinear piezoelectric coupling β, and reciprocal time constant λ. Also, these parameters become relevant for optimization studies, such that minor changes in their values result in significant changes in the recovered power, which is essential for a proper preliminary optimization design. The innovative finding in this study is that the electrical propriety of piezoelectric coupling (piezoelectric coupling on the electrical equation) can be a starting point to an enhancement process.

Table 2: Summary for the system's sensitivity to each parameter for the bistable energy harvesting models and their conditions to influence. ing to the table, it is noteworthy that uncertainties were only for the sensitivities parameters. Although remarkably nonlinear coupling generates a higher power, as already discussed, the uncertainty range is also higher than the case with nonlinear coupling, highlighting the need to verify the nonlinear term.

Concluding remarks

This work presented a global sensitivity analysis of a bistable energy harvesting system to identify the input parameters that lead to the most variability in the harvested power. Global sensitivity analysis based on Sobol indices was employed to calculate the parameters' individual and combined effect, exploring various dynamic behavior scenarios. We studied a classical bistable energy harvester, then added a nonlinear electromechanical coupling and, finally, a proposed model considering asymmetric coefficient at potential energy, slopping of the system's base, and the nonlinear coupling. All these considerations were done to study a sophisticated model, providing interesting insights for the proposal of a new bistable energy harvesting device.

The results here reaffirm the need to quantify and propagate uncertainties in bistable energy harvesting systems, presenting high sensitivity and a nontrivial interpretation by complex dynamic behavior. Furthermore, the framework proposes a powerful tool to de-velop a robust design, forecast, and optimization in bistable energy harvesting systems.

Fig. 1 :

 1 Fig. 1: Schematic illustration of the Piezo-Magneto-Elastic beam energy harvester.

Fig. 2 :

 2 Fig. 2: The displacement time-series for different amplitude of excitation. The time-series in (a) is classical bistable energy harvester (β = δ = φ = 0), in (b) classical bistable energy harvester with nonlinear coupling (β = 1.0 and δ = φ = 0), in (c) asymmetric bistable energy harvester (δ = 0.15, φ = -10 • and β = 0), and in (d) asymmetric bistable energy harvester with nonlinear coupling (δ = 0.15, φ = -10 • and β = 1.0).

Fig. 3 :

 3 Fig. 3: Methodology schematic for the sensitivity analysis of the energy harvester about parameter variations.

Fig. 4 :

 4 Fig. 4: Numerical simulations results for the bistable energy harvester with linear electromechanical coupling: (a) bifurcation diagrams of voltage as a function of the excitation amplitude several values of the excitation frequency. The sweeping up the amplitude excitations are presented with cool colors while sweeping down appears in hot colors. (b) Phase portrait and Poincare section for different excitation amplitudes and Ω = 0.8.

Fig. 5 :

 5 Fig. 5: Sobol' indices associated with the mean output power for the bistable energy harvester with linear coupling. First for nominal regular steady-state dynamics (f = 0.147 and Ω = 0.8): in (a) the first-order based on MC and PCE methods and in (b) the second-order based on PCE method. Second for nominal irregular steady-state dynamics (f = 0.083 and Ω = 0.8): in (c) the first-order based on MC and PCE methods and in (d) the secondorder based on PCE method. Finally, (e) presents the main Sobol' indices exploring several scenarios with different excitation amplitudes based on the PCE method.

Fig.

  Fig. Bifurcation diagrams of voltage energy harvester with nonlinear electromechanical coupling as a function of the excitation amplitude, for six values of nonlinear coupling term: β ∈ {0.5, 1.0, 1.5, 2.0, 2.5, 3.0} and Ω = 0.8. The sweeping up the amplitude excitation is presented with cool colors, while sweeping down appears in hot colors.

Fig. 8 :

 8 Fig. 8: Sobol indices based on PCE method associated with the mean output power of the bistable energy harvester for several scenarios of excitation amplitude with nonlinear coupling, taken into account: (a) β = 0.5, (b) β = 1.0, (c) β = 1.5, (d) β = 2.0, (e) β = 2.5, (f ) β = 3.0.

Fig. 9 :

 9 Fig.9: Illustration of the potential energy function, of symmetric and asymmetric configurations, for δ equals 0.15 and -0.15. Furthermore, an uncertainty range of δ is considered.

Fig. 10 :

 10 Fig. 10: Sobol indices based on PCE method associated with the mean output power of the bistable energy harvester with asymmetric potential.

Fig. 11 :

 11 Fig. 11: Sobol indices based on PCE method associated with the mean output power of the bistable energy harvester with asymmetric potential and bias angle.

Fig. 12 :

 12 Fig.12: Sobol indices based on PCE method associated with the mean output power of the bistable energy harvester with asymmetric potential, bias angle, and nonlinear piezoelectric coupling.

Figure 13 Fig. 13 :

 1313 Figure13presents the uncertainty propagation of mean power in the excitation amplitude spectrum for the studied cases at different confidence levels. Accord-

Table 1 :

 1 Comparison of the computational cost of the Sobol indices via MC and PCE methods.

	case	method	sample size	CPU time* (seconds)	MAE
	periodic	MC PCE	10 000 1 000	≈ 86 400 ≈ 2 700	1.59%
	chaotic	MC PCE	20 000 2 000	≈ 162 000 ≈ 6 000	2.57%
	*Intel i7-9750H 2.60GHz 8GB 2666GHz DDR4

Acknowledgements

The authors gratefully acknowledge, for the financial support given to this research, the following Brazilian agencies: Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) -Finance Code 001; São Paulo Research Foundation (FAPESP), grant number 19/19684-3; Brazilian National Council for Scientific and Technological Development (CNPq) grant number 306526/2019-0 ; Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) grants 210.167/2019, 211.037/2019 and 201.294/2021.

Code availability

The simulations reported in this paper used the computational code STONEHENGE -Suite for Nonlinear Analysis of Energy Harvesting Systems [START_REF] Norenberg | STONE-HENGE -suite for nonlinear analysis of energy harvesting systems[END_REF]. This code is available for free on GitHub [START_REF] Norenberg | STONE-HENGE -Suite for Nonlinear Analysis of Energy Harvesting Systems[END_REF].

Compliance with ethical standards

Conflict of interest

The authors declare they have no conflict of interest.