Joachim Demuynck
email: joachim.demuynck@ugent.be

Michel De Paepe

Henk Huisseune

Roger Sierens

Jeroen Vancoillie

Sebastian Verhelst

Investigation of the influence of engine settings on the heat flux in a hydrogen-and methane-fuelled spark ignition engine

Keywords: hydrogen, methane, internal combustion engine, experimental, heat transfer, thermopile

Hydrogen-fuelled internal combustion engines are a possible solution to make transportation more ecological. Only emissions of oxides of nitrogen (NO x ) occur at high loads, being a constraint for power and efficiency optimization. A thermodynamic model of the engine cycle enables a cheap and fast optimization of engine settings. It needs to accurately predict the heat transfer in the engine because the NO x emissions are influenced by the maximum gas temperature. However, The existing engine heat transfer models in the literature are developed for fossil fuels and they have been cited to be inaccurate for hydrogen. We have measured the heat transfer inside a spark ignited engine with a thermopile to investigate the heat transfer process of hydrogen and to find the differences with a fossil fuel. This paper describes the effects of the compression ratio, ignition timing and mixture richness on the heat transfer process, comparing hydrogen with methane. A convection coefficient is used to separate the effect of the temperature difference between the gas and the wall from the influence of the gas movement and combustion. The paper shows that the convection coefficient gives more insight in the heat transfer process in a combustion engine despite the assumptions involved in its definition. The comparison between hydrogen and methane demonstrates, in contrast to what is believed, that the heat loss of hydrogen can be lower.

Climate change and energy supply both are issues high on the political agenda these days. A solution has to be found to move away from the fossil fuel based energy supply of today. The transportation sector in particular is very dependent on fossil fuels. Therefore, researchers have tried to find alternative fuels for several years.

The hydrogen-fuelled combustion engine is one of the possible alternatives which is investigated at the Transportation Technology group at Ghent University. These engines are an attractive alternative for the current drive trains because the combustion properties of hydrogen enable a high efficiency for all engine loads by using several operational strategies [START_REF] Verhelst | Hydrogen-fueled internal combustion engines[END_REF]. Moreover, hydrogen engines have near-zero noxious and zero greenhouse gas emissions at the tailpipe. The initial research at Ghent University was focused on the experimental optimization of engine operational strategies for maximum power and efficiency, with ultra low emissions of oxides of nitrogen (NO x ) [START_REF] Sierens | Influence of the injection parameters on the efficiency and power output of a hydrogen fueled engine[END_REF][START_REF] Verhelst | Increasing the power output of hydrogen internal combustion engines by means of supercharging and exhaust gas recirculation[END_REF][START_REF] Verhelst | Efficiency comparison between hydrogen and gasoline, on a bi-fuel hydrogen/gasoline engine[END_REF]. The focus shifted to numerical research with the development of a thermodynamic model of the engine cycle, called the GUEST-code (Ghent University Engine Simulation Tool) [START_REF] Verhelst | A study of the combustion in hydrogen-fuelled internal combustion engines[END_REF][START_REF] Verhelst | A quasi-dimensional model for the power cycle of a hydrogen-fuelled ICE[END_REF].

Such a thermodynamic model of the engine cycle enables a cheap and fast optimization of engine settings for operation on hydrogen. Several submodels are necessary to solve the conservation equations of energy and mass: a combustion, a turbulence and a heat transfer model among other things. The last one is important to accurately simulate the NO x emissions which are influenced by the maximum gas temperature. These emissions can occur in hydrogen internal combustion engines at high loads, being an important constraint for power and efficiency optimization.

A lot of research has been performed on the heat transfer in engines, but clear results remain scarce. Consequently, several research groups are still looking into it, e.g. [START_REF] Puzinauskas | Ignition and boost effects on largebore engine in-cylinder heat transfer[END_REF][START_REF] Rakopoulos | Experimental and theoretical study of the short term response temperature transients in the cylinder walls of a diesel engine at various operating conditions[END_REF][START_REF] Galindo | Description of a heat transfer model suitable to calculate transient processes of turbocharged diesel engines with one-dimensional gas-dynamic codes[END_REF][START_REF] Sanli | The influence of engine speed and load on the heat transfer between gases and in-cylinder walls at fired and motored conditions of an idi diesel engine[END_REF][START_REF] Lounici | Investigation on heat transfer evaluation for a more efficient two-zone combustion model in the case of natural gas SI engines[END_REF]. Experimental investigations were mainly done between 1970 and 1990 and the

M A N U S C R I P T A C C E P T E D

ACCEPTED MANUSCRIPT focus has shifted towards numerical work. The results of CFD (Computational Fluid Dynamics) codes are validated with the available measurements in the literature [START_REF] Rakopoulos | Critical evaluation of current heat transfer models used in CFD in-cylinder engine simulations and establishment of[END_REF] or compared with simulation results of semi-empirical models [START_REF] Mohammadi | Estimation of instantaneous local heat transfer coefficient in spark-ignition engines[END_REF]. The agreement with the measurements is poor, therefore, efforts focus on the development of new wall functions [START_REF] Rakopoulos | Critical evaluation of current heat transfer models used in CFD in-cylinder engine simulations and establishment of[END_REF][START_REF] Han | A temperature wall function formulation for variable-density turbulent flows with application to engine convective heat transfer modeling[END_REF][START_REF] Jung | Development of an in-cylinder heat transfer model with compressibility effects on turbulent prandtl number, eddy viscosity ratio and kinematic viscosity variation[END_REF]. Quasi-dimensional models for the engine cycle are still the standard, using semi-empirical submodels for the heat transfer. The latter were mainly developed before 1990 and are reported to be inaccurate for new combustion types like HCCI (Homogeneous Charge Compression Ignition) [START_REF] Chang | New heat transfer correlation for an HCCI engine derived from measurements of instantaneous surface heat flux[END_REF][START_REF] Garcia | A new heat release rate (HRR) law for homogeneous charge compression ignition (HCCI) combustion mode[END_REF][START_REF] Soyhan | Evaluation of heat transfer correlations for HCCI engine modeling[END_REF] and alternative fuels like hydrogen [START_REF] Shudo | Applicability of heat transfer equations to hydrogen combustion[END_REF][START_REF] Wei | A study on transient heat transfer coefficient of in-cylinder gas in the hydrogen fueled engine[END_REF][START_REF] Demuynck | On the applicability of empirical heat transfer models for hydrogen combustion engines[END_REF]. Consequently, experimental research is still needed to understand the heat transfer process inside modern combustion engines.

We are investigating experimentally the differences in the heat transfer of hydrogen and fossil fuels (in this case methane) in order to improve the existing heat transfer models. Several methods to measure the heat flux in an internal combustion engine exist, but new types are still being developed as well [START_REF] Marr | A fast response thermocouple for internal combustion engine surface temperature measurements[END_REF]. Therefore, one of the purposes of our research is to compare several sensors based on the measurement accuracy and construction feasibility as described in [START_REF] Demuynck | Experimental research on the heat transfer inside a hydrogen combustion engine: Evaluation and construction of measurement methods[END_REF]. One of the selected sensors has been used to measure heat transfer rates and wall temperatures in a spark ignition engine [START_REF] Demuynck | Local heat flux measurements in a hydrogen and methane spark ignition engine with a thermopile sensor[END_REF]. A convection coefficient is now used to gain more insight into those measurements. This paper first describes an investigation of motored measurements. Then, it presents the effects of the compression ratio, ignition timing and mixture richness on fired measurements. The paper ultimately compares measurements of hydrogen and methane with the same indicated power output.

Experimental method

Equipment

An overview of the test setup is shown in Figs. 1 and2. The engine used in this research is a four-stroke single-cylinder spark ignited gas engine based on a CFR (Cooperative Fuel Research) engine operated at a constant speed of 600 rpm. It is equipped with PFI (port fuel injection) and has a variable compression ratio. The details of the engine are given in Table 1. Fuel injection and ignition timing are controlled by a MoTeC M4Pro electronic control unit.

The heat flux measurements were carried out with a Vatell HFM-7 sensor which consists of a thermopile (heat flux signal, HFS) and an RTD (substrate temperature signal, RTS). The sensor has a response time of 17 µs. The Vatell AMP-6 amplifier was used as a current source for the RTD and as an amplifier for both output signals. The sensor is calibrated by Vatell (see below) and it has been used before to measure the heat flux in the intake manifold [START_REF] Tillock | Measurement and modeling of thermal flows in an aircooled engine[END_REF][START_REF] Bauer | Averaged and time-resolved heat transfer of steady and pulsating entry flow in intake manifold of a spark-ignition engine[END_REF][START_REF] Wimmer | Heat transfer to the combustion chamber and port walls of IC engines -measurement and prediction[END_REF] and the cylinder [START_REF] Wimmer | Heat transfer to the combustion chamber and port walls of IC engines -measurement and prediction[END_REF] of an internal combustion engine.

As the test engine is easily accessible, the heat flux sensor was successively installed in three different positions under fired operation (P2, P3, P4 as shown in Fig. 1). These openings are at the same height in the cylinder wall and are equally distributed around the circumference of the cylinder. The spark plug was placed in position P1. The heat flux sensor could be mounted in P1 as well in the case of motored operation, because of the absence of the spark plug. The measurement position is mentioned in the caption of all the figures in the following sections.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

In-cylinder pressures were measured with a water-cooled Kistler 701A piezoelectric pressure sensor (in P2 or P4) and inlet pressure with a Kistler 4075A20 piezoresistive pressure sensor. This inlet pressure was used to reference in-cylinder pressure. A 12 bit data acquisition card was used to sample both the heat flux and pressure signals.

It is triggered by a crank angle encoder every 0.5 • CA, resulting in a sampling rate of 7.2 kHz. Next, Gas flows were measured with Bronkhorst Hi-Tec F-201AC (fuel) and F-106BZ (air) flow sensors. Finally, type K thermocouples were used to measure inlet and exhaust gas temperatures.

Data reduction

Most of the literature concerning engine heat transfer calculates a convection coefficient (h) out of the measured heat flux (q) and the difference between the gas temperature (T g ) and the wall temperature (T w ) with equation 1. This equation is actually only valid for steady convection, so the heat transfer process in the engine has to be assumed to be quasi-steady. This means that the instantaneous heat flux is proportional to the instantaneous temperature difference between the gas and the wall. This concept is not entirely justifiable because the heat transfer in an engine is transient and there is a time lag between the temperature difference and the heat flux due to the heat capacity of the gas [START_REF]Heat transfer in the cylinders of reciprocating internal combustion engines[END_REF]. However, an alternative approach has not been proposed so far, so we have calculated a time-resolved convection coefficient with equation 1, evaluating the insight it provides.

h = q (T g -T w ) (1)
The idea is to separate the temperature difference between the gas and the wall from other influencing factors on the heat transfer process. The convection coefficient represents the influence of the gas velocity, the turbulence, the propagating flame front and the gas properties (Prandtl number, viscosity, conductivity and heat capacity). It will be shown further on that the convection coefficient gives more information than the heat flux trace alone.

The combustion gases are assumed to behave like ideal gases. Therefore, the bulk gas temperature is calculated with the following equation of state:

T g = p • V/m • R.
The in-cylinder pressure (p) is measured and the volume (V) can be calculated out of the crank position. The mass (m) can only be determined during the closed part of the combustion cycle, being the sum of the measured incoming mass (air and fuel) and the residuals and neglecting blow-by effects because of the large number of piston rings [START_REF] Verhelst | A study of the combustion in hydrogen-fuelled internal combustion engines[END_REF]. No incoming mass goes directly to the exhaust manifold because the test engine does not have a valve overlap. The residual mass is therefore determined with the equation of state at EVC (exhaust valve closure) using the measured cylinder pressure and assuming that the in-cylinder temperature is equal to the measured exhaust temperature. The specific gas constant (R) at IVC (intake valve closure) can be calculated out of the mass average of the specific gas constants of the air, the fuel and the residual gases. This value is used until the beginning of the combustion. At the end of the combustion R is equal to that of the combustion products. In between, the specific gas constant is calculated with a linear interpolation. The moment where the combustion begins and ends is determined with an apparent rate of heat release analysis.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

An indication of total cycle heat losses (Q h ) is calculated, assuming that the measured heat flux trace at the measurement position occurs evenly over the entire cylinder wall. The measured heat flux is multiplied with the total available in-cylinder surface for each sampling point and all these values are summed up to get the total cycle heat loss.

Error analysis

To judge the quality of the measurement results, a thorough error analysis has been carried out following the methods described in Taylor [START_REF] Taylor | An introduction to error analysis: the study of uncertainties in physical measurements[END_REF]. The error analysis is based on the accuracy of the measurement equipment which is given in Table 2. The pressure measurement circuit is calibrated with a dead-weight tester in the laboratory. The given uncertainty of a pressure sensor is the standard deviation of the calibration which is repeated 50 times.

The HFM sensor is not included in Table 2 because its accuracy depends on several coefficients which are used to calculate the heat flux and wall temperature out of the measured HFS and RTS signal (see HFM manual [START_REF] Vatell | Heat flux microsensor manual[END_REF]). These coefficients and their uncertainty are determined during the calibration of the sensor at the factory of Vatell. The resulting maximum relative error on the heat flux is 3 % during the compression and expansion stroke and 7 % during the intake and exhaust stroke. The maximum relative error on the wall temperature is 5 %.

The uncertainty on the compression ratio ( ) is negligible. The absolute error on the crank angle and ignition timing is 0.5 • CA. A change of 0.5 • CA causes a maximum change of 1 % in the calculated cylinder volume, so this is taken as the relative error on the cylinder volume. The maximum relative error on the gas flow rates is 4 % for air and hydrogen and 9 % for methane. The resulting error on the air-to-fuel equivalence ratio (λ) is 5 % for hydrogen and 10 % for methane.

The indicated work output (W i ) is calculated with an accuracy of 2.5 %. The uncertainty on the heat loss through the cylinder walls (Q h ) is estimated to be 7.6 %. The heat input of the fuel (Q f ) has a relative error of 9 % for methane and 4 % for hydrogen. The relative error on the indicated efficiency (η i = W i /Q f ) is 9.3 % for methane and 4.7 % for hydrogen. The maximum relative error on the heat loss ratio (φ h = Q h /Q f ) is 11.8 % for methane and 8.6 % for hydrogen.

A change of 10 % in λ causes a change of 0.3 % in the calculated specific gas constant (R) of a methane-air mixture. Changing λ with 4 % changes the R of for a hydrogen-air mixture with 0.6 %. These errors on R are negligible. The maximum relative error on the calculated gas temperature is 6 %. The convection coefficient is calculated with an uncertainty of 11 %. These uncertainties are calculated assuming that the cylinder temperature at EVC is equal to the exhaust temperature as explained earlier. A sensitivity analysis has been performed to check the influence of the choice of this temperature. Varying the cylinder temperature at EVC with 20 % changes the residual gas fraction with 10 %, but the calculated bulk gas temperature and convection coefficient only change with 1 %. This influence is acceptable for the calculations in this paper.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Results and discussion

Measurements under motored conditions with a variation in the compression ratio and in-cylinder mass are discussed first. Then, fired measurements on hydrogen are presented. The compression ratio ( ), ignition timing (IGN) and air-to-fuel equivalence ratio (λ) were varied in that case. Ultimately, a comparison is made between hydrogen and methane measurements with the same indicated power output.

The measurements on hydrogen were performed with wide-open-throttle (WOT), varying the power by changing the air-to-fuel equivalence ratio. It is not possible to control the power output in a sufficient range by changing λ in the case of methane, so a throttle in the intake manifold was used. Ignition was at MBT-timing (minimum spark advance for Maximum Brake Torque), unless mentioned otherwise. The ignition timing is always indicated in the graphs.

This paper presents the heat flux traces of the ensemble cycle. The heat flux of the ensemble cycle at a certain crank angle is the average of all the heat fluxes of an entire measurement set at that crank angle. A measurement set always consisted of 35 cycles.

There are three elements in equation 1: the measured heat flux, the calculated convection coefficient and the temperature difference between the calculated gas and the measured wall temperature. These three elements will be plotted in the same graph throughout this paper. The heat flux will be plotted at the top, the temperature difference in the middle and the convection coefficient at the bottom.

Motored operation

The heat flux is first investigated in the engine under motored operation. In this case only air flows into the engine and there is no combustion. This allows to gain insight into the effect of the gas flow and turbulence on the heat transfer process, without the interference of the effect of a propagating flame front. Figure 3 shows that there are only small differences between the heat flux or convection coefficient at the different measurement positions, indicating that bulk processes mainly influence the heat flux under motored operation.

Variation of the compression ratio

The heat flux, the temperature difference and the convection coefficient for a variation in the compression ratio, with constant intake temperature and pressure, are plotted in Fig. 4. The heat flux and wall temperature are measured at P1. Overall, the traces peak just before TDC (top dead centre), shortly before the cylinder pressure peaks. The peak in the heat flux trace (at the top of the graph) rises if the compression ratio increases from 6 to 15, as expected. The maximum in the gas temperature increases from 375 • C to 475 • C, mainly because of the rise in the in-cylinder pressure. The wall temperature only increases by 12 • C. As a result, the peak temperature difference between the gas and the wall increases with the compression ratio. The increase in the heat flux cannot only be caused by the increasing ∆T , because the convection coefficient is significantly higher as well. Consequently, this confirms that gas motion and turbulence have an effect on the heat transfer.

Figure 4 shows that the trend is the other way around near the end of the expansion period, the heat flux decreases with an increasing compression ratio. This can partially

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

be explained by a decrease in ∆T which is caused by two effects. First, the gas temperature is lower (after 415 • CA) because the cylinder pressure is almost the same and the volume is smaller. Second, the wall temperature is almost constant during the cycle and thus still higher. Other effects must have an influence on the heat flux as well, because there is a larger difference in the traces of the convection coefficient.

The trace of h near the end of the expansion period demonstrates that the assumption of a quasi-steady heat transfer process is not justifiable, because the convection coefficients become negative.

Variation of the throttle position

A throttle is used in the intake manifold to investigate the effect of a reduced incylinder mass on the heat transfer. The throttle position can be varied between 0 • and 90 • . The effect of three throttle position on the cylinder pressure ( = 10) is plotted in Fig. 5. Varying its position until 65 • does not have an influence, because the very low engine speed results in a small air flow rate. The corresponding q, ∆T and h at P1 are shown in Figure 6. The heat flux increases if the throttle opening is varied between 0 • and 70 • . This can partially be explained by an increase in the temperature difference between the gas and the wall. The turbulence caused by the throttle must have an influence as well and this effect is visible in the increased convection coefficient. The increase in the gas temperature for a decrease in the pressure is unexpected. However, the in-cylinder mass decreases as well and the ratio of p over m in the equation of state is higher as a consequence.

Further throttling reduces the trapped mass more which causes the cylinder pressure to drop significantly. Consequently, the gas temperature and the heat flux decrease. The convection coefficient decreases as well.

Fired operation 3.2.1. General observations

A typical trace of the wall temperature during an engine cycle under fired operation is plotted in Fig. 7 (at P2 with λ = 1 and = 8). The wall temperature under motored operation is plotted as well to show the effect of combustion. Under fired operation, the wall temperature decreases somewhat during the intake stroke because of the incoming mixture of air and fuel which cools down the wall. The wall temperature then increases by 5 to 10 • C (depending on the engine operation point) during the combustion period and decreases again towards the end of the expansion stroke and during the exhaust stroke. In contrast, under motored condition the temperature almost remains constant over the entire engine cycle. The differences between the wall temperature with a change in the engine operation are the peak-to-peak value and the steady state level. The wall temperature will not be plotted any more in this paper. All the graphs will show the temperature difference (∆T ) between the gas and the wall (T g -T w , see equation 1) because this drives the heat transfer process.

The measured heat flux, the calculated convection coefficient and the temperature difference for a reference measurement are plotted in Fig. 8 to discuss the general observations. The measurement plotted is taken at P2, λ is equal to 1 and the compression ratio is 8.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Two vertical dotted lines are drawn on the graph. Line 1 corresponds with the maximum of the convection coefficient and line 2 coincides with the maximum of the temperature difference. The intersection of these two lines with the trace of the heat flux is of interest. Line 1 corresponds with the end of the initial steep rise in the heat flux trace. Line 2 coincides with the maximum in the heat flux trace. The maximum heat flux thus occurs at the moment that the temperature difference peaks. However, the initial increase in the heat transfer is mainly caused by the propagating flame front that passes over the measurement position. This effect is visible in the trace of the convection coefficient.

Although the definition of a convection coefficient is based on the assumption of a quasi-steady process, which is not the case as shown with the motored measurements, it gives more insight into the heat transfer process than the measured heat flux trace alone. The observations discussed in this section are valid for all the measurements plotted further on, where line 1 and 2 are inserted as well.

In contrast to the motored measurements, there are significant differences in the heat flux between the different measurement positions. Both the instant where the heat flux starts to increase and the peak value differ, caused by the difference in the arrival time of the propagating flame front, see [START_REF] Demuynck | Local heat flux measurements in a hydrogen and methane spark ignition engine with a thermopile sensor[END_REF]. For a variation in the measurement position, the trace of the convection coefficient does not provide extra insight, because the trace of the temperature difference is the same for all the positions.

Variation of the compression ratio

The compression ratio is varied from = 6 to = 10. The traces of q, ∆T and h are plotted in Fig. 9 (at P2 and λ = 1.5). The cylinder pressure of the three measurements is plotted in Fig. 10. The measurements have the same ignition timings (at TDC), so that only the compression ratio is varied.

There is an overall increase in the heat flux for a higher compression ratio, except near the end of the expansion stroke. According to Heywood [START_REF] Sanli | The influence of engine speed and load on the heat transfer between gases and in-cylinder walls at fired and motored conditions of an idi diesel engine[END_REF] several operational properties change with an increasing : the gas pressure and peak burned gas temperature increase, gas motion increases, combustion is faster and the gas temperature late in the expansion stroke decreases. The first four properties increase the heat transfer, the last one decreases it. Fig. 9 shows that the peak in the difference between the gas temperature and the wall temperature does not increase with an increasing compression ratio for these measurements, which is also the case for the gas temperature itself. This is an unexpected result. The higher peak pressure for a higher compression ratio is in the calculation of T g , with the equation of state, countered by the smaller in-cylinder volume, resulting in the same peak gas temperature. The cylinder volume is not only larger at the moment of maximum pressure because of the lower compression ratio, but also because of the fact that the pressure peaks later in the combustion cycle. However, the gas temperature is higher during the beginning of the combustion process, resulting in an increase in the heat flux after TDC. The temperature difference between the gas and the wall for a compression ratio of 6 eventually rises to the same maximum like for the higher compression ratios, boosting the heat flux around a crank angle of 380 • CA to a level which is higher than expected.

The difference in the maximum of the convection coefficient going from = 6 to = 8 is not significant. In contrast, there is a significant increase if is changed

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

between 8 and 10. Consequently, the increase in the heat flux trace from = 6 to = 8 is mainly caused by an increase in the temperature difference around 367 • CA (from 473 to 688 • C). The difference between the gas and the wall temperature does not increase that strongly any more when the compression ratio is further increased (from 8 to 10). Thus, the increase in the heat flux trace for a higher compression ratio in that region is mainly caused by the faster combustion process and the intensified gas motion. This is reflected in an increase in the peak of the convection coefficient (from 1420 to 1805 W/m 2 K).

Variation of the ignition timing

The ignition timing is varied around the MBT-timing to check its influence on the heat transfer, it is advanced and retarded with 12 • CA from MBT. The measured and calculated traces are plotted in Fig. 11 (at P2, = 8 and λ = 1.5). There is a significant increase in the peak of the measured heat flux trace with an advanced ignition timing (from 129 to 158 W/cm 2 ). However, it cannot be concluded which effect causes this increase, because the change in the peak of the temperature difference (from 1455 to 1520 • C) and convection coefficient (from 1335 to 1446 W/m 2 K) are not significant if one considers the estimated errors of 6 and 11 %. Consequently, the shift in time of the heat flux is the most important influence of the changing ignition timing.

Variation of the mixture richness

The air-to-fuel equivalence ratio (λ) is varied from 1 to 2. The measured and calculated traces are plotted in Fig. 12 (at P2 and = 8). The ignition is kept constant again, so there is only a variation in λ. The heat flux increases with a decreasing λ (increasing mixture richness).

Both the convection coefficient and the temperature difference between the gas and the wall increase. There is a larger time difference between the maximum of the convection coefficient and the temperature difference for a lower λ. This causes a more stretched heat flux trace. Consequently, the maximum of the heat flux trace does not occur at the moment of maximum ∆T anymore for λ = 2. The heat flux does stay at a high level until the temperature difference reaches its maximum. A clear time shift in all three plotted traces can be noted if the mixture gets leaner. This visualizes the effect of the mixture richness on the combustion rate. As the mixture gets leaner, the combustion rate slows down, being enforced by the constant ignition timing.

A comparison between methane and hydrogen

A comparison is made between the heat transfer of hydrogen and methane measurements with the same indicated engine power output. The measured heat flux, the calculated convection coefficient and the calculated temperature difference between the gas and the wall temperature for two comparable measurements of hydrogen and methane are plotted in Fig. 13. Extra information about the measurements is given in Table 3. The measurements are taken at P2, the compression ratio being equal to 9. The hydrogen measurements have an air-to-fuel equivalence ratio of 1 and 2. The measurements of methane are always stoichiometric and the power output is equalized with the hydrogen measurements by throttling. Ignition is at MBT-timing.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

The heat flux of the stoichiometric hydrogen measurement (Fig. 13(a)) is clearly higher than the heat flux of the corresponding methane measurement. However, the lean hydrogen measurement (Fig. 13(b)), has a heat flux which is lower than that of methane. In the case of methane, the heat flux does not decrease a lot when the incylinder mass is reduced by throttling. The total cycle heat loss of the measurement with the lowest power output is only 12 % smaller than the one with the highest power output, although the power output is decreased with 26 %. This is also reflected in the increase of φ h from 27 to 30 %. The difference between the two hydrogen measurements is much bigger. The total cycle heat loss decreases with 64 % and φ h with 41 %. The indicated efficiency of the hydrogen measurements increases from 23 to 29 % as a result.

The trace of the temperature difference between the gas and the wall is about the same for the two methane measurements. The reduced in-cylinder mass and pressure cancel each other out in the equation of state. This is the reason why the heat flux does not decrease a lot if the power output is reduced by throttling. The lean hydrogen measurement has a lower peak in the temperature difference and therefore a lower heat flux. Another reason for the decreased heat flux is the slower combustion in the case of the lean hydrogen-air mixture, resulting in a lower convection coefficient.

In the literature it is stated that the heat transfer of hydrogen is higher compared to a fossil fuel (e.g. gasoline [START_REF] Wei | A study on transient heat transfer coefficient of in-cylinder gas in the hydrogen fueled engine[END_REF] or methane [START_REF] Shudo | Thermal efficiency analysis in a hydrogen premixed combustion engine[END_REF]). The measurements discussed above show that this is not always the case, the heat loss of hydrogen combustion can be lower if a lean mixture is used.

The results above demonstrate the advantage of hydrogen. A high engine efficiency can be attained at low engine loads because of three reasons. First, there is no need to use a throttle because of the wide flammability limits, resulting in low pumping losses. Second, the heat losses are low because the peak gas temperature of the lean combustion is low. Third, the thermodynamic properties of the working fluid (κ) change, resulting in an increased theoretical efficiency [START_REF] Verhelst | Hydrogen-fueled internal combustion engines[END_REF].

Conclusions

The heat flux and wall temperatures have been measured in a hydrogen-fuelled spark ignited combustion engine under motored and fired conditions. The effect of a variation in the compression ratio, ignition timing and mixture richness has been investigated. A convection coefficient has been calculated in an attempt to separate the effect of the temperature difference between the gas and the wall on the heat transfer from other effects. The heat flux of hydrogen has been compared to a fossil-fuel, in this case methane, to investigate the differences in the heat transfer process. Measurements with the same indicated power output have been compared. This paper has demonstrated that the convection coefficient gives more insight in the heat transfer process, despite the assumptions involved in its definition. However, the motored measurements showed that the convection coefficient can have negative values which indicate that its definition should be adjusted. The main conclusions of the paper are as follows:

• The convection coefficient represents the influence of the combustion process and in-cylinder movement on the heat transfer process. The maximum in the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
convection coefficient trace corresponded with the end of the initial rising slope of the measured heat flux for all the measurements. The largest part of the initial increase in the heat flux is caused by the combustion process and the gas movement. The heat flux starts to rise at the instant that the flame front passes over the measurement position. The maximum in the heat flux trace occurred together with the maximum in the temperature difference between the in-cylinder gases and the cylinder wall. This was not the case if the instant between the flame passage and the maximum in the temperature difference was spread too much in time. The maximum in the heat flux trace then corresponded with the maximum in the convection coefficient. The heat flux did remain at a high level until the temperature difference reached its maximum.

• The heat flux increased for an increasing compression ratio. This increase is caused by an increase in the temperature difference between the gas and the wall between a compression ratio of 6 and 8. The increase between a compression ratio of 8 and 10 is mainly caused by an increase in the gas movement and a faster combustion process.

• An advanced ignition timing caused an increase in the heat flux. It is not clear if this is caused by an increase in the temperature difference or the convection coefficient, because these differences were not significant. Consequently, the main influence of the changing ignition timing is that it shifts the traces in time.

• The heat flux and convection coefficient increased with an increasing mixture richness. The comparable measurements between hydrogen and methane have demonstrated that the mixture richness has a large influence on the heat transfer process in a combustion engine. The heat flux decreased significantly if the mixture richness of the hydrogen-air mixture was reduced to operate the engine at a low load. Consequently, the total cycle heat loss decreased, resulting in a higher engine efficiency. Reducing the cylinder mass to lower the engine power output did not have the same effect on the heat flux.

  Abbreviations• CA degree crank angle

Figure 1 :Figure 2 :Figure 3 :Figure 4 :igure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :

 12345678910 Figure 1: Cross-section of the CFR engine, P1: spark plug, P1-P4: sensor positions, IV: intake valve, EV: exhaust valve

Figure 11 :

 11 Figure 11: The most important change in the heat flux caused by a varying ignition timing is the shifting in time (P2, = 8, λ = 1.5). The differences in the traces of ∆T and h are not significant if one considers the measurement uncertainties.

Figure 13 :

 13 Figure 13: The heat flux of hydrogen and methane measurements with the same indicated power output are compared (P2, = 9, IGN=MBT).

Table 1 :

 1 Geometrical properties and valve timing of the CFR engine

	Bore	82.55 mm
	Stroke	114.2 mm
	Connecting rod length	254 mm
	Swept volume	611.7 cm 3
	IVO	17 • CA ATDC
	IVC	26 • CA ABDC
	EVO	32 • CA BBDC
	EVC	6 • CA ATDC

Table 2 :

 2 The accuracy of the measurement equipment

	Variable	Device	Accuracy
	Heat flux gain	Vatell AMP-6	±3.6%
	Wall temperature gain	Vatell AMP-6	±1.5%
	In-cylinder pressure	Kistler 701A	±1%
	Intake manifold pressure Kistler 4075A20	±2.5%
	Air flow rate	Bronkhorst F-106BZ ±1%FS
	Hydrogen flow rate	Bronkhorst F-201AC ±1%FS
	Engine speed	ASTM tachometer	±6 rpm
	Atmospheric temperature		±0.5 • C
	Atmospheric pressure		±50 Pa

Table 3 :

 3 Comparison between comparable measurements of hydrogen and methane plotted in Fig.13

	)

fuel λ throttle Q f uel (J) W i (J) η i (%) Q h (J) φ h (%

Acknowledgements

The authors of this paper like to acknowledge the suggestions and technical assistance of Koen Chielens and Patrick De Pue. The research is carried out in the framework of a Ph.D. which is funded by a grant (SB-81139) of the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWT-Vlaanderen). The experimental equipment is funded by a Research Grant (1.5.147.10N) of the Research Foundation -Flanders (FWO). These financial supports are gratefully acknowledged. J. Vancoillie acknowledges the Research Foundation -Flanders (FWO) for the Ph.D. grant 09/ASP/030.


Nicolas Sutton Charani

Sébastien Destercke
email: sebastien.destercke@hds.utc.fr

Nicolas Sutton-Charani
email: nicolas.sutton-charani@hds.utc.fr

Thierry Denoeux
email: thierry.denoeux@hds.utc.fr

Classification trees based on belief functions

published or not. The documents may come

Introduction

Decision trees [START_REF] Breiman | Classification And Regression Trees[END_REF] (classification trees for categorical labels and regression trees for numerical ones) are popular classifiers, due to their simplicity, efficiency and readability. The construction of usual decision trees relies on classical probability theory. However, classical methods are not always fully adequate to deal with some problems. Among these problems are [START_REF] Abellan | Upper entropy of credal sets. Applications to credal classification[END_REF] the fact that all kinds of uncertainties (either in input or output) cannot be modeled faithfully by classical probabilities and (2) the fact that frequencies of occurrence are only sensible to proportions in sample and not to its size.

Outside the fact that the relation between input and outputs may be non-deterministic, a classifier may have to deal with three different possible levels of uncertainty: in inputs, in outputs, and uncertainty due to the fact that the learned classifier is an estimation of the ideal one, due to a limited amount of knowledge or data. In this work, we mainly address the third kind, where the estimation quality translates into imprecision of belief functions.

Belief function theory [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF] offers a convenient framework to deal with all these problems. For instance, Elouedi et al. [START_REF] Elouedi | Belief decision trees: theoretical foundations[END_REF] propose different ways to adapt decision trees in the TBM framework to deal with uncertain outputs during the tree construction. In this work, we extend another approach also using belief functions proposed by Denoeux and Skarstein Bjanger [START_REF] Denoeux | Induction of decision trees from partially classified data using belief functions[END_REF] that can cope with uncertain outputs and imprecision arising from limited sample size. In this sense, this approach is closer to some imprecise probabilistic approaches [START_REF] Abellan | Upper entropy of credal sets. Applications to credal classification[END_REF] that naturally integrate sample size information in their construction.

As the Denoeux and Skarstein Bjanger method only concerns 2-class problems, we extend this methodology to any number of classes. For multi-class problems, we propose three ways of doing such an extension:

• combining belief functions provided by sets of 2-class trees [START_REF] Quost | Pairwise Classifier Combination using Belief Functions[END_REF] • direcltly building multinomial belief functions using the Imprecise Dirichlet Model (IDM) [START_REF] Utkin | Extensions of belief functions and possibility distributions by using the imprecise dirichlet model[END_REF] • direcltly building multinomial belief functions using Denoeux's proposal [START_REF] Denoeux | Constructing belief functions from sample data using multinomial confidence regions[END_REF] Section 2 presents the needed background about decision trees and Denoeux and Skarstein Bjanger method. Section 3 then extends this methodology to the multiclass case. Finally, in Section 4 we compare new classifiers with classical CART and discuss the effects of method parameter on experiment results basis'.

Background

In this section, we present the necessary background.

Decision trees

Let (X,Y ) be a random vector where X = (X 1 , ..., X J ) ∈ X = X 1 × . . . × X J represents the features (continuous or discrete) and Y ∈ Y = {Y 1 , . . . ,Y K } the class to predict. From a sample E = {(X 1 ,Y 1 ), . . . , (X n ,Y n )}, decision tree methods build iteratively a model of (X,Y ) by building a partition of X . Here, we consider binary trees (i.e. CART-like models), where each split provide two children.

The method works as follow: from a root node containing the whole learning sample, the optimal split (among all the variables and their values) in term of information gain is searched. The information gain IG corresponding to splitting on variable X k with value α is computed as follows:

IG(k, α) = i(t 0 ) -p L i(t 1 ) -p R i(t 2 ) (1)
where i(t) is an impurity measure of a node t, t 0 the root node, t 1 and t 2 its child nodes, p L is the proportion of the samples in t 0 verifying the condition X k < α (i.e., p L = n L/n where n is the sample size in t 0 and n L the samples such that X k < α). p R = 1p L is the sample proportion not verifying it. The selected splitting value (k, α) is then the one maximizing IG (computed as a gain in purity).

The method is then applied recursively on each child nodes until no possible information gain superior to a preestablished threshold can be made. In this case, the node becomes a leaf predicting the most frequent class of the leaf sample.

The information gain (or impurity measure) is calculated using Gini-index for CART algorithm or Shanon entropy for C4.5's one (Quinlan [11]). Both of these functions measure the homogeneity in term of classes. They both use the frequencies of the different classes in the node samples, however these frequencies do not depend on the sample size (provided class proportions remain the same). In contrast, the method of Denoeux and Skarstein Bjanger and its impurity measure, that we recall now, do change with the sample size.

Denoeux and Skarstein Bjanger method for 2-class datasets

This method use the same principle as CART, but differs in the computation of information gain: first, they use mass functions instead of simple frequencies and second, they use an impurity measure mixing non-specifity (imprecision) and conflict (variability).

To build the mass functions, they use Dempster's approach to Bernouilli (DaBt) trials that induces the following mass function:

   m DaBt ({Y 1 }) = n 1 n+1 m DaBt ({Y 2 }) = n 2 n+1 m DaBt (Y ) = 1 n+1 , (2)
where n is the number of samples and n 1 , n 2 are the number of samples whose class is Y 1 ,Y 2 , respectively. They then propose to use the following impurity measure [START_REF] Klir | Uncertainty and information: foundations of generalized information theory[END_REF], applied to m DaBt :

U λ (m) = (1 -λ )N(m) + λ D(m) (3)
where

N(m) = ∑ A⊆Y m(A) log 2 |A| is the non-specificity and D(m) = -∑ A⊆Y m(A) log 2 BetP(A)
the variability. The two parts are weighted by λ ∈ [0, 1]. Note that as the size n of the sample increases, m(Y ) (the imprecision) decreases. When using U λ as impurity measure i(t), the information gain (1) can be negative. This gives a natural stopping criterion when building the tree, that is no split is done if all possible information gains are negative. Usually, λ value can be fixed by an optimization procedure (see Section 4). Table 1 shows results obtained with CART-classification trees and with classification trees based on Denoeux and Skarstein Bjanger method. The stopping criteria is the following: keep splitting while IG > β for classical CART-trees (IG > 0 for the one based on U λ ) and while the children nodes of the split contains a minimum of 10 samples. Classical CART are optimized on the threeshold β and U λ -based on the parameter λ , using a 10-fold cross-validation procedure on training samples. Results show that the methods achieve comparable accuracies. Dempster's approach to Bernouilli's trial is not adapted to cases where the output has more than 2 classes. Therefore, we propose 3 ways to handle such cases: break up the classification problem containing K classes (K ≥ 3) into C 2 k 2-classes problems using Quost combination of binary classifiers [START_REF] Quost | Pairwise Classifier Combination using Belief Functions[END_REF], use the IDM approach or the Denoeux multinomial model.

Multi-classes cases

Combinations of binary classifiers

In [START_REF] Quost | Pairwise Classifier Combination using Belief Functions[END_REF], Quost presents a way to handle multi-class classification problems by combining classifiers built on sub-samples containing only 2 classes. That is, he proposes to learn (from the corresponding sub-sample) a conditional belief function for each pair {Y i ,Y j }, 1 ≤ i < j ≤ K of classes and to combine them into a global belief function over Y using an optimisation procedure.

Here, we propose to combine this method with decision trees issued from Denoeux and Skarstein Bjanger method, using the latter as base classifier to learn conditional belief functions (those belief functions are assigned by the DaBt applied to leaves class proportions). Decision trees are well adapted to this kind of combination, since they are simple classifiers. However, note that λ optimization becomes an issue, as K(K-1) /2 (i.e., a quadratic number) of classifiers have to be learned at each optimization step.

IDM

The IDM was introduced in the "imprecise probability" framework by Walley [START_REF] Walley | Inferences from multinomial data: Learning about a bag of marbles[END_REF]. Note that although belief functions can be interpreted as imprecise probabilities, this is far from being their only possible interpretation (and it is the one adopted here). However the IDM turn out to yield a belief function as output, hence it can be used in our framework. The IDM imprecision is controlled by an hyper-parameter s ∈ R + . From a random sample Y 1 , ...,Y n , Walley showed that the lower predictive probability distribution on Y is P(Y k |N, s) = n k/n+s where n k is the number of times Y k has been observed. The corresponding mass function is such that:

m IDM (Y j ) = n j /(n + s) j = 1, ..., K m IDM (Y ) = s/(n + s) (4)
Note that we find back Equation ( 2) for K=2 and s=1. Using m IDM , U λ can be applied to measure the impurity in a node and multi-class trees can thus be created. Analytical form of U λ applied to m IDM can be derived as:

U λ (m IDM ) = (1 -λ )s n + s log 2 (K) - λ n + s K ∑ k=1 n k log 2 [ Kn k + S K(n + s) ] (5)
However even if this model is simple, it's not easy to interpret it within the belief function frameworks. Also, it should be noted that the IDM imprecision only depends on the number of samples n, and not on their distribution over Y . This is not the case for the multinomial construction of Denoeux that offers a tractable and interesting alternative.

Denoeux multinomial model

Denoeux [START_REF] Denoeux | Constructing belief functions from sample data using multinomial confidence regions[END_REF] proposes to use Goodmans confidence intervals [START_REF] Goodman | On simultaneous confidence intervals for multinomial proportions[END_REF] to build a multinomial belief function. The first step is to build probability intervals [START_REF] De Campos | Probability intervals: a tool for uncertain reasoning[END_REF] (probability lower and upper bounds over singletons) and then to transform them into belief functions.

Let (X 1 ,Y 1 ), ..., (X n ,Y n ) be an iid sample where Y k ∈ Y = {Y 1 , ...,Y K }, those probability intervals [P - k , P + k ] are given, for Y k (k=1,...,n), as:

P - k = q + 2n k - √ ∆ k 2(n + q) P + k = q + 2n k + √ ∆ k 2(n + q) ( 6
)
where q is the quantile of order 1α of the chi-square distribution with one degree of freedom, and where ∆ k = q(q + 4n k (n-n k ) /n). As shown in [START_REF] Denoeux | Constructing belief functions from sample data using multinomial confidence regions[END_REF], the lower confidence measure (i.e., P -

(A) = max( ∑ Y k ∈A P - k , 1 -∑ Y k / ∈A P - k )
) built using these regions in the case where K = 2 or 3 are belief functions, hence induce corresponding mass functions.

Note that the built belief functions follow the Hacking principle (see [START_REF] Denoeux | Constructing belief functions from sample data using multinomial confidence regions[END_REF] for details), but that the solution for K = 2 is this time not equivalent to Eq. ( 2).

In the case K > 3, the Möbius inverse of P -takes some negative values so P -is not a belief function. Different methods involving linear programing are proposed in [START_REF] Denoeux | Constructing belief functions from sample data using multinomial confidence regions[END_REF] to approximate it into a belief function. Also, in the special case where the classes are ordinal, Denoeux proposes an algorithm restricted to a certain set of focal elements. A valid predictive bba is obtained. These belief functions can then be used with U λ to build multi-class trees.

Experiments

We start by comparing the classifiers performances, and then discuss the effect of the λ parameter in U λ .

Comparison between classifiers

We compare the three proposed extensions with the classical probabilistic CART algorithm. Table 2 shows 3 multi-class UCI datasets characteristics. Table 3 presents experimental results on the previous datasets comparing the accuracy of 4 types of classifier: The tree growing strategy is the following: keep splitting while

•
• IG > β for CART and IG > 0 for the tree based on U λ • the children nodes sample size is greater than 10 CART trees and the ones using U λ are optimized by a 10-fold cross-validation procedure: for CART we optimize the threshold β and for trees based on U λ we optimize λ . None of the trees are post-pruned, as we are only interested in accuracies of each models, and not in their simplicity (defining a proper pruning strategy for U λ based decision trees remains the matter of further research). As we can see the different classifiers' accuracies are competitive although the time computations are longer with the ones using belief functions probably because of the number of focal elements wich can be higher than with probabilities (where they are restricted only to singletons)

Discussion about λ

Figures 1 study the impact of λ in terms of tree complexity (using the classical number of leafs criterion) and in terms of accuracy on the UCI dataset "Breast Tissue". They show that this complexity increases with λ , confirming that 1λ can be interpreted as the importance given to the lack of samples in a node (i.e., to non-specifity N(m)) and to the propensity of IG to be negative. This suggests that optimization (here, a 10-fold cross-validation procedure) should also integrate tree complexity as a criterion. The best accuracy level is here obtained for small lambdas which give small trees. Accuracy seems to be little influenced by lambda value, suggesting that lambda value should be kept low.

Conclusion

In this paper, we have extended Denoeux and Skarstein Bjanger method to the multiclass case, proposing three ways to do so. The IDM is not really based on belief function and may result in too simple belief functions; Denoeux multinomial model is more elaborated, fits better with a belief function approach, but requires heavy computational efforts (asking to solve a linear program for each split test); 2-class decomposition is efficient, but makes the interpretability of results possibly harder (and, in any case, longer), as it builds a quadratic number of decision trees.

We have shown that the presented methods have prediction power comparable to classical methods. However the present work is only a starting point with many perspectives: one of the major interest of using belief functions is the ability to handle uncertain data in inputs or outputs, a feature we shall integrate to the present methods in future works (using, for example, extensions of EM-algorithm to learn trees [START_REF] Ciampi | Growing a tree classifier with imprecise data[END_REF][START_REF] Denoeux | Maximum likelihood estimation from uncertain data in the belief function framework[END_REF]). Another interesting extension would be to adapt this model to continuous outputs and to regression problems. Additional experiments exploring the method behavior should also be done.

  standard CART trees based on Gini index (CART) • trees based on U λ with m IDM (IDM) • combination of 2-classes trees based on U λ (combination) • trees based on U λ with m Multinomial (multi)

Fig. 1

 1 Fig. 1 number of nodes/λ (top) and error rate/number of nodes (bottom)

Table 1

 1 Accuracies of trees depending of the used impurity measure

	Data set	Number of features	standard CART error rate	trees based on U λ error rate
	Blood transfusion	4	23.5%	24.2%
	Statlog heart	13	28%	25.7%
	Tic-tac	9	21.5%	11.5%
	Breast-cancer	10	5.9%	4.7%
	Pima	8	27.3%	25.1%
	Haberman	3	26.6%	26%

Table 3

 3 Accuracies of trees depending of the masses assignement model

	Data set	Number of features	Number of classes	learning sets size	test sets size
	Iris	4	3	113	37
	Balance scale	4	3	469	156
	Wine	13	3	134	44
	Table 2 UCI data sets used in experiments


Luciano Fusina
email: fusina.luciano@unibo.it

Giandomenico Nivellini

Silvia Spezzano

Prof Luciano Fusina

The ν ν ν ν 1 and ν ν ν ν 3 band system of 15 NH 3

Keywords: 15NH3, high-resolution infrared spectra, stretching vibrations, rovibrational perturbations, spectroscopic parameters 15 NH 3, high resolution infrared spectra, stretching vibrations, rovibrational analysis, rovibrational perturbations, spectroscopic parameters

published or not. The documents may come

Introduction

Ammonia is one of the most extensively studied molecules since it collects the interest of experimental, theoretical and astrophysical scientists. It is the prototype molecule for the study of inversion doubling in C 3v molecules. The interpretation of its spectroscopic behaviour still poses a difficult challenge to both experimental and theoretical approaches.

Ammonia has been detected in the interstellar medium, in remote stars [1], and in the atmospheres of the giant planets of the solar system and of their moons, i.e. the 14 N/ 15 N isotopic ratio in Jupiter's atmosphere from observation of NH 3 transitions has been recently reported [2]. Most of the spectroscopic investigations, based on spectra recorded at high resolution and with different spectroscopic techniques, have concerned 14 NH 3 and comparatively a few high resolution studies have been performed on its isotopologues, 15 NH 3 , 14 ND 3 and 15 ND 3 . In the last ten years there has been a revival of interest in the computation of the energies of vibration and ro-vibration levels up to high excitation, fostered by the improvement of the computational methodologies [3][4][5][6][7][8][9][10][11]. Accurate reproduction of the observed transitions supports the reliability of predictions for unobserved ones and can help the interpretation of data collected by the recent or programmed astronomical missions [START_REF] Wotten | Proc. SPIE[END_REF][START_REF] Becklin | Proc. SPIE[END_REF][START_REF] Pilbratt | Proc. SPIE[END_REF]. Experimental studies have been devoted to the enlargement of the spectroscopic data-bases of 14 NH 3 and 15 NH 3 [START_REF] Urban | [END_REF][16][17][18][19][20][21][22][23]. The need for spectroscopic data of ammonia isotopologues has been evidenced in the latest vibrational energy calculation [11]. In fact, these data are essential to test the ability of the calculated potential energy surfaces to predict, with the same precision achieved for 14 NH 3 , the energy pattern of the other isotopologues. Up to now the spectroscopic investigations have been mainly devoted to the ground state or to the low energy bending states of ammonia, since the stretching levels are strongly perturbed by nearby v 2 = 4, v 4 = 2, or v 2 = 2, v 4 = 1 levels [24].

Precise values for the ground state parameters of 15 NH 3 have been determined by the simultaneous analysis of inversion and inversion-rotation transitions measured in the microwave [25] and far infrared [26][27][28] regions, and of combination differences (GSCD) from forbidden transitions to the v 4 = 1 and v 2 = 2 states [29]. The v 2 = 1, 2, 3 (A' and A" symmetry), v 4 = 1 and v 2 = v 4 = 1 (E' and E" symmetry) bending states have been characterised by means of Fourier transform spectra [29][30][31], laser Stark [32,33], diode laser [34] and two-photon [35,36] spectra of the ν 2 band, and of inversion spectra in the v 4 = 1 state [37], measured in the microwave region. Far-infrared laser emissions have also been observed by optically pumping 15 NH 3 samples with a tunable CO 2 -laser [38,39]. Contrary, the v 1 = 1 (A' and A" symmetry) and v 3 = 1 (E' and E" symmetry) stretching states are less characterized. In fact, only the Q branch of ν 1 recorded with the coherent anti-Stokes Raman technique has been measured and assigned [40]. At higher energies, the v 1 = v 2 = 1 [41], v 1 = v 3 = 1 [19], v 1 = 1, v 4 = 2 [20], and v 1 = 3, v 4 = 2 [22] states have been reported.

Aim of the present study is to improve the understanding of the fundamental stretching vibrations in 15 NH 3 , thus contributing to determine the complete anharmonic potential surface for ammonia. We report for the first time the analysis of the high resolution infrared spectra recorded in the region of the ν 1 and ν 3 band system of 15 NH 3 . All the transitions between levels with the same or the opposite inversion symmetry assigned in the spectra will be analysed simultaneously, taking into account all the symmetry-allowed interactions between vibration-inversion-rotation levels in the theoretical model [42 and refs. therein].

The paper is structured as follows. The experimental details are given in section 2; the description of the spectra and the procedure adopted for the assignments are given in section

Experimental details

The spectra of 15 NH 3 were recorded in the region between 3050 and 3650 cm -1 on a BOMEM DA3.002 Fourier transform spectrometer capable of a resolution of 0.004 cm -1 . The achieved unapodized resolution was about 0.010 cm -1 . The Doppler width of 15 NH 3 at room temperature is 0.0095 cm -1 at 3300 cm -1 . Spectra were recorded using a 0.18 m long glass cell with NaCl windows, with sample pressure of 2.0 mbar, or using a multipass cell with an optical path length of 5.0 m and 8.0 m, at pressures of 2.0 and 10.0 mbar, respectively. An isotopically pure (99%) sample of 15 NH 3 was purchased from Stohler Isotope Chemicals.

Absorption lines of residual H 2 O in the spectrometer were used for calibration [43]. Both the wavenumber precision and accuracy of the spectra are estimated to be about 0.5×10 -3 cm -1 .

Description of the spectrum and assignments

The ν 1 parallel band of 15 NH 3 consists of a crowded Q-branch region, which extends for about 6 cm -1 around 3331 cm -1 , due to the proximity of the origins of s ← a and a ← s bands, at 3332.4 and 3334.1 cm -1 , respectively. The P and R branches appear as groups of lines separated by about 20 cm -1 (≈2B) with resolved, blue-degraded, K structures for both a ← s and s ← a components up to J' = 5 or 7, respectively. Transitions with the same J" and K" values are ca. 1.7 cm -1 apart, a value which corresponds to the sum of the inversion splittings in the ground and in the vibration excited states, owing to the a ← s and s ← a selection rule. Transitions are easily identified up to P (7) and R (6). For higher J'' the K branch structures seem to be perturbed, particularly at high K. In fact, with increasing K the lines are red-degraded at low K, form a branch head and reverse their trend. In addition, the spacing between the lines and their intensity become irregular, see Fig. 1. As a consequence, transitions for K > 6 have been identified only with the assistance of precise GSCD calculated with the parameters in [28].

[insert figure 1 about here] The spectrum of the ν 3 perpendicular band is characterised by strong p P K (J) branches which are red-degraded and by r R K (J) branches which are blue-degraded, for increasing J.

Intense p Q K (J) and r Q K (J) transitions are also observed towards higher and lower wavenumbers with respect to the band centres at about 3435 cm -1 . Lines are red-degraded for increasing J in both r Q and p Q sub-branches. Transitions belonging to the ν 3 band are easily identified since they always appear as close doublets, corresponding to the two components, s ← s and a ← a, of the inversion splitting. The separation of the s and a doublet components is ca. 0.40 cm -1 and, due to the selection rule of the transitions, corresponds to the difference of the inversion splitting in the ground and v 3 = 1 states. In fact, the energy differences between the a and s levels for J = K = 0 in the ground and in v 3 = 1 are 0.758 cm -1 and 0.332 cm -1 , respectively. Fig. 2 illustrates portions of the spectrum where p Q K (J) and r Q K (J) transitions are present. For the p Q 4 (J) and r Q 2 (J) lines the separation between s and a components decreases slowly with increasing J, with the s component at higher wavenumbers, because of the larger inversion splitting in the ground state compared to that in the excited state. A more irregular pattern was observed for the s and a separation within the doublets of the p Q 3 (J) sub-branch, see Fig. 3. The separation is alternatively large and small for odd and even J, up to J = 6, but for J = 7 is larger than expected, the s component being in all cases at higher wavenumbers. For J = 8, corresponding to the last observed p Q 3 (J) doublet, the separation becomes anomalously small and the a component appears at higher wavenumbers. The p P 3 (J) and p R 3 (J) sub-branches exhibit anomalies analogous to those observed in the p Q 3 (J) sub-branch.

[insert figures 2 and 3 about here] These anomalies could not be rationalised in terms of vibration-rotation-inversion interactions between the v 1 = 1 and v 3 = 1 levels only, and should be attributed also to interactions with other nearby vibration excited states. A similar behaviour was observed in the spectra of the stretching fundamentals for 14 NH 3 , and was in part rationalized taking into account the interactions with the v 4 = 2, l = 0, ±2 states [24]. Transitions involving these excited states were in fact observed in the spectra. Other interactions with v 2 = 2, v 4 = 1 (s and a) and v 2 = 4 (s) were suggested, but only a limited number of transitions to the latter state were identified. The overtones and combination states present in the region between 3000 and 3600 cm -1 for 15 NH 3 are listed in Table 1.

Lines which have not been assigned either to ν 1 or ν 3 transitions are scattered all over the investigated spectral region, the most intense being located in the lower wavenumber range. The feature observed between 3210 and 3230 cm -1 , in good agreement with the theoretical predictions for 2ν 4 0 [11], was assigned to the Q branch transitions. A series of doublets separated by about 1 cm -1 , from about 3238 cm -1 towards higher wavenumbers, was observed. They were assigned to p Q 1 (J) transitions and represented the starting point for the analysis of the whole 2ν Extensive use of GSCD, calculated from the parameters in [28], has been made to check the correctness of the assignments, or to predict the wavenumbers of weak P, Q, or R lines from the corresponding assigned R, Q, or P partners. An important aid for the initial assignments was the relative intensity of the lines, which is governed by the nuclear spin-statistical weights, and by the absence of K'' = 0 lines for even J'' (s) or odd J'' (a) in the various sub-branches. Finally, toward the end of the analysis, a large number of perturbation allowed transitions were identified in all bands, once their positions could be calculated with reliable precision from improved values of the parameters. However, only 256 lines were contained in the data set, i.e. those pertaining to series of at least 3 perturbation allowed transitions. In total, 1606 transitions have been assigned, 243 of which to the ν 1 band with J' max /K' max values up to 9/8, 666 to the ν 3 band, with J' max /K' max up to 13/13, 222 to the 2ν 4 0 band, with J' max /K' max up to 9/6, and 425 to the 2ν 4 ±2 band, with J' max /K' max up to 12/11. The maximum values of J' and K' of the transitions were due essentially to the decreasing line intensity for high quantum numbers, and to the effect of the interactions.

Discussion

Details of the analysis

The anomalies observed in the spectral features of 15 NH 3 are similar to those reported in a previous investigation of the same spectral range for 14 NH 3 [24]. That study evidenced the necessity of analysing simultaneously the experimental data adopting a theoretical model All the assigned transitions were analysed on the basis of a model similar to that adopted for the same system of bands in 14 NH 3 [24] using an upgraded version of the program described in details in Ref. [44]. For each value of J the program constructs two Hamiltonian matrices of dimensions 4×(2J+1): one is set up for levels of A 1 ', A 2 ' and A 1 ", A 2 " symmetry, and the other corresponds to one of the two (degenerate) blocks of E' and E" symmetry. The two blocks are numerically diagonalised to obtain the energy levels. The diagonal matrix elements, containing the usual contributions up to the 6 th power in the angular momentum operators, were given by (i)

E v (J,K,l) = (i) 0 v E + (i) B v [J(J + 1) -k 2 ] + (i) C v k 2 -2 (i) (Cζ) v kl -(i) D vJ [J(J + 1)] 2 -(i) D vJK [J(J + 1)]k 2 -(i) D vK k 4 + (i) η vJ [J(J + 1)]kl + (i) η vK k 3 l + (i) H vJ [J(J + 1)] 3 + (i) H vJK [J(J + 1)] 2 k 2 + (i) H vKJ [J(J + 1)]k 4 + (i) H vK k 6 (1)
where (i) refers to the parity of the level with respect to the inversion, s or a ; k = ±K; v = 1, 3 or 4 ; l = l 3 or l 4 and all terms containing l vanish for the v 1 = 1 and v 4 = 2 0 levels. The off-diagonal matrix elements are reported in the Appendix, where H mn represents a group of terms in the Hamiltonian containing m vibrational operators q r and/or p r , and n rotational operators J α [45]. We remind that, except the Fermi constants, the subscript of a constant represents the ∆k and ∆l selection rules for the levels connected by the interaction term in the Hamiltonian. Moreover, a constant without the left superscript (i) corresponds to an s ↔ a Our model included interaction terms only up to the fourth-order of approximation, but the ∆k = ±3 and ∆k = ±6 interaction terms containing the coefficients α, β and η 3 , were neglected since no A 1 /A 2 splitting for K = 3n was observed.

A weighted least-squares analysis was performed simultaneously on all the assigned transitions to derive the parameters of the excited states. In the least-squares fits a unitary statistical weight has been attributed to each transition that corresponds to an isolated line, while the weight has been reduced to 1/n if n transitions are assigned to the same wavenumber. Finally, all the transition wavenumbers that differed from their corresponding calculated values by more than 0.1 cm -1 were excluded from the data set in the final cycle of the fitting procedure.

The parameters obtained from a polynomial fit reported in Table II of [28] were used to derive the energy of the ground state levels necessary to calculate the wavenumbers of the transitions.

Spectroscopic parameters

The analysis was accomplished in an iterative process, which consists in improving the values of the excited state parameters by means of the assigned lines, predicting the wavenumbers of higher J and K unassigned transitions to be searched in the spectra, enlarging the data set, and so on. Once all the transitions that could be unambiguously assigned were C , η J , and η K were refined. The other parameters, for instance the distortion constants D and H, were constrained to their ground state values or to zero. In the following steps, higher order parameters were included in the model according to some selective criteria. For example, both the J and K dependences of a lower order constant already determined, i.e. J tt q and K tt q , were added simultaneously and, if present, both the s and a constants. The significant fourth order interaction parameters were determined together with the distortion constants D while the sextic distortion constants H were kept constrained. At the end of each fit we checked that the newly, as well as the previously determined parameters, were statistically determined, that they improved the standard deviation of the fit and that anomalous correlations between the parameters were not present. The parameters not fulfilling these requirements were constrained either to zero or to their corresponding ground state values in the following fits. Empirical damping factors were introduced to approach the convergence region in cases of severe correlations. They were set equal to zero in the final cycles of the refinement procedure.

Table 2 lists the parameters derived from the analysis of 1514 transitions contained in the final fit, 20 of v 1 =1, 26 of v 3 = 1, 20 of v 4 = 2 0 , 25 of v 4 = 2 ±2 and 9 interaction constants, including the H's which were not refined. All the parameters in Eq. ( 1) and in the Appendix which are not reported in Table 2 were set equal to zero and constrained in the refinement procedure. The number of transitions rejected in the final fit was equal to 92, ca. 6% of the assigned transitions. The chosen rejection limit, 0.1 cm -1 , is much larger than the estimated experimental precision, 0.0005 cm -1 , to discard only poorly measured blended lines or transitions perturbed by interactions with levels other than those considered in our model. The exclusion of lines corresponding to high J or K values in a given sub-branch, which could be Parameters in Table 2 can be compared to those for 14 NH 3 (see Table 7 in [24]). The unperturbed vibration term values follow the expected decreasing trend due to the 15 N isotopic substitution. The K dependent parameters C and Cζ have nearly identical values.

Contrary, with few exceptions, the distortion constants D differ substantially. Different choices of refined interaction parameters prevent a meaningful comparison. However, it is worth noting that the v 1 = 1/ v 4 = 2 0 Fermi interaction constants of the two isotopologues differ by only about 3%. Finally, the vibration energies of the analysed states, including the values for v 1 = 1, 3333.22 (s) and 3334.16 (a) cm -1 , and v 4 = 2 0 , 3210.43 (s) and 3212.12 (a)

cm -1 , perturbed by the Fermi interaction, are in excellent agreement with those calculated from the potential energy surface, see Table 1, thus supporting its reliability. It should be pointed out that while the differences between calculated and experimental term values are about 0.02 cm -1 for the E states, equal to the standard deviation in Ref. [11], they are larger for the A states affected by the Fermi resonance, 0.04-0.09 cm -1 . Since the values calculated from the potential energy surface are expected to have the same precision, the differences for 3 seem to support the last suggestion. Our results confirm the difficulties encountered in the analysis of the stretching region for ammonia and reaffirm the need for further experimental work to identify and characterise all the vibrational states present in the studied spectral range.

A table with a list of the transition wavenumbers, assignments and residuals from the simultaneous fit associated with all the analysed bands is available from the Bologna authors. 15 NH 3 illustrating the K progression of a few q R K (J) and q P K (J) sub-branches. The J sub-branches at lower wavenumbers are s ← a transitions, those at higher wavenumbers a ← s transitions. The K assignment is given. Lines identified by an asterisk are r R 3 (3) transitions of ν 3 . Note also that the K = 0 lines are observed depending on even/odd and a/s symmetry. Experimental conditions: pressure, 2.0 mbar; path length, 5.0 m.

Figure 2. Details of the ν 3 fundamental of 15 NH 3 illustrating the K progression of a few p Q K (J) and r Q K (J) sub-branches. The J assignment and the a/s symmetry are given.

Experimental conditions: pressure, 2.0 mbar; path length, 5.0 m. 2 )

a a J K s s K l J H H hc l K K q q J J q K J F ν ν ± + ± ± = ± + + + ± ( ) 4

  theoretical model used in the analyses and the results of the fits are described in section 4 which concludes the manuscript.

  into account all the symmetry-allowed interactions between and within the vibrationally excited states [42 and ref. therein]. Fermi-type resonances are effective between v 1 = 1 and v 4 = 2 0 and between v 3 = 1 and v 4 = 2 ±2 , Coriolis interactions between v 1 = 1 and v 3 = 1, v 1 = 1 and v 4 = 2 ±2 , v 3 = 1 and v 4 = 2 0 , l-type resonances within v 3 = 1 and v 4 = 2 ±2 and k-type resonances within all vibration states.

  an (s) or (a) superscript represents interaction between levels of the same parity, i.e. (s) q 22 , constant for both the <(a); v 1 = 1; v 3 = 0 0 ; J k | H 21 /hc | (s); v 1 = 0; v 3 = 1 ±1 ; J k ± 1> and <(s); v 1 = 1; v 3 = 0 0 ; J k | H 21 /hc | (a); v 1 = 0; v 3 = 1 ±1 ; J k ± 1> matrix elements.

  systematic investigation to derive the set of parameters that best reproduces the experimental data was carried out. All the constants up to the third order of approximation, E 0 , B, C, Cζ,

3 = 1 ,

 31 the refinement of further constants, was carefully avoided. The standard deviation of the fit, 0.034 cm -1 , is about 70 times larger than the estimated experimental precision, 0.0005 cm -1 . This result clearly suggests that the other bands present in the investigated region should be analysed and the interaction model adopted in the analysis should be enlarged to include the v 2 = 2 / v 4 = 1 and v 2 = 4 excited bending states.The values of the refined s and a rotation parameters are very similar in all the vibration states and close to the corresponding ones in the ground state. However, all the parameters should be considered as effective, including the contributions due to interactions with neighbouring levels not considered in our model. With the exception of the values for v large differences are observed between the centrifugal distortion constants D for the s and a states, and between corresponding parameters of the ground and excited states.

  could be ascribed to several factors. These could be the high correlation between the Fermi constants and the deperturbed term values and the impossibility to determine the J and K dependences affecting the Fermi constant value, leading to a worse reproduction by the model of the transitions assigned to the A states with respect to the E states. The statistics for the fitted line positions of the separate bands listed in Table

Figure 1 .

 1 Figure captions

Figure 3 .

 3 Figure 3. Details of the ν 3 fundamental of 15 NH 3 illustrating the K progression of the p Q 3 (J)

  In the high wavenumber region lines without a recognizable pattern are present: they probably belong to the 4ν 2 s band but could not be assigned with certainty.

	their identification was easier.
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	8

4 ±2 band. Fortunately, both the parallel 2ν 4 0 and perpendicular 2ν 4 ±2 bands are located in a spectral region where only lines of ν 1 already assigned are present, so

Table 1

 1 Vibrational states in the region between 3000 and 3600 cm -1 for 15 NH 3

		Molecular Physics	Page 18 of 24
	Vibrational state	Calculated / cm -1 [11]	This work
	v 2 = 1, v 4 = 1 ±1 (s) (E)	3179.54
	v 2 = 1, v 4 = 1 ±1 (a) (E)	3487.77
	v 4 = 2 0 (s) (A)	3210.51	3210.43
	v 4 = 2 0 (a) (A)	3212.08	3212.12
	v 4 = 2 ±2 (s) (E)	3233.94	3233.92533(631)
	v 4 = 2 ±2 (a) (E)	3235.32	3235.33849(707)
	v 1 = 1 (s) (A)	3333.27	3333.22
	v 1 = 1 (a) (A)	3334.25	3334.16
	v 3 = 1 (s) (E)	3435.12	3435.14291(547)
	v 3 = 1 (a) (E)	3435.45	3435.47468(542)
	v 2 = 4 (s) (A)	3438.70

URL: http://mc.manuscriptcentral.com/tandf/tmph

Table 2 .

 2 Spectroscopic parameters (in cm -1 ) for the ν 1 /ν 3 /2ν 4 system of15 NH 3

	a

a The uncertainties as 1σ given in parentheses refer to the last significant digits. Parameters without uncertainty have been constrained to the corresponding ground state values.

Table 3

 3 Statistics for fitted line positions a

					Molecular Physics	Page 20 of 24
	Band		Number of fitted / assigned lines	rms (cm -1 )
	ν 1			s ← a,s	275 / 293 135 / 148	0.037 0.030
				a ← s,a	140 / 145	0.042
	2ν	0 4		s ← a,s	213 / 222 125 / 127	0.042 0.045
				a ← s,a	88 / 95	0.037
	ν 3			s ← a,s	642 / 666 321 / 326	0.029 0.027
				a ← s,a	321 / 340	0.030
	2ν	4 ±	2	s ← a,s	384 / 425 212 / 228	0.036 0.040
				a ← s,a	172 / 197	0.030
	Global fit		1514 / 1606	0.034
	a The two inversion parities of the lower state indicate symmetry allowed and perturbation allowed
	transitions, respectively

URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgements

L.F. and G.N. acknowledge the financial support from the Università di Bologna and from the Ministero della Ricerca e dell'Università, PRIN Trasferimenti di energia e di carica: dalle collisioni ai processi dissipativi.